動作環境
GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04.4 LTS desktop amd64
TensorFlow v1.7.0
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)
scipy v0.19.1
geopandas v0.3.0
MATLAB R2017b (Home Edition)
ADDA v.1.3b6
gnustep-gui-runtime v0.24.0-3.1
PyMieScatt v1.7.0
check_intersection_3D_180512.py v0.1: geometry (3D) > rayとtriangleの交差判定 > check_intersection_3D_180512.py > v0.1 > C: intersect3D_RayTriangle()からの移植, v0.2: 戻り値の変更
処理概要
- rayとtriangleの交差判定(3D)を行う
対象データ
対象データ > 座標関係
に記載の３ケースを使用。
geometry_intersection_3D_180513.py v0.3
- テスト関数 Test_reachingIntersection()追加
geometry_intersection_3D_180513.py
import numpy as np
import sys
'''
v0.3 May 18, 2018
- add Test_reachingIntersection()
v0.2 May 18, 2018
- add check_reachingIntersection()
v0.1 May 18, 2018
- separated from [check_intersection_3D_180512.py] v0.3
'''
# return values
RET_DEGENERATE = -1
RET_DISJOINT = 0
RET_INTESECT_UNIQ_POINT = 1
RET_IN_SAME_PLANE = 2
RET_NOT_REACHING_PLANE = 3
def intersect3D_RayTriangle(ray, tri):
# ported from
# http://geomalgorithms.com/a06-_intersect-2.html#intersect3D_RayTriangle()
retNULL = np.array([]) # intersect point of ray and plane
# get triangle edge vectors and plane normal
vecU = tri[1] - tri[0]
vecV = tri[2] - tri[0]
vecN = np.cross(vecU, vecV)
if (np.count_nonzero(vecN) == 0): # triangle is degenerate
return RET_DEGENERATE, retNULL # do not deal with thise case
vecRayDir = ray[1] - ray[0] # ray direction vector
w0 = ray[0] - tri[0]
dota = -np.dot(vecN, w0)
dotb = np.dot(vecN, vecRayDir)
if (abs(dotb) < sys.float_info.epsilon):
# ray is parallel to triangle plane
if (abs(dota) < sys.float_info.epsilon): # ray lies in triangle plane
return RET_IN_SAME_PLANE, retNULL
else:
return RET_DISJOINT, retNULL # ray disjoint from plane
# get intersect point of ray with triangle plane
wrk = dota / dotb
if (wrk < 0.0): # ray goes away from triangle
return RET_DISJOINT, retNULL
# for a segment, also test if (wrk > 1.0) => no intersect
retI = ray[0] + wrk * vecRayDir # intersect point of ray and plane
# is retI inside triangle?
uu = np.dot(vecU, vecU)
uv = np.dot(vecU, vecV)
vv = np.dot(vecV, vecV)
vecW = retI - tri[0]
wu = np.dot(vecW, vecU)
wv = np.dot(vecW, vecV)
valD = uv * uv - uu * vv
# get and test parametric coords
valS = (uv * wv - vv * wu) / valD
if (valS < 0.0 or valS > 1.0): # retI is oustide tri
return RET_DISJOINT, retNULL
valT = (uv * wu - uu * wv) / valD
if (valT < 0.0 or (valS + valT) > 1.0): # retI is outside tri
return RET_DISJOINT, retNULL
return RET_IN_SAME_PLANE, [ray[0], retI] # retI is in T
def calc_norm_of_ray(aray):
wrk = aray[1] - aray[0]
return np.linalg.norm(wrk)
def check_intersection(ray, tri):
retNULL = np.array([]) # intersect point of ray and plane
res, retI = intersect3D_RayTriangle(ray, tri)
if len(retI) == 0:
return res, retI
nrm1 = calc_norm_of_ray(ray) # original ray
nrm2 = calc_norm_of_ray(retI) # ray reaching to the plane
if nrm1 < nrm2: # not reaching the plane
return RET_NOT_REACHING_PLANE, retNULL
return res, retI
def check_reachingIntersection(ray, tri):
res, retI = check_intersection(ray, tri)
return len(retI) > 0
def Test_reachingIntersection():
# Triangular vertices
TV1 = np.array([-42, -8, 60])
TV2 = np.array([43, 41, 60])
TV3 = np.array([21, -54, 60])
tri = [TV1, TV2, TV3]
# ray with intersection (reaching to the triangle)
ray_with_intsct = [np.array([0, 0, 0]), np.array([8.0666, -7.7, 66])]
# ray with intersection (not reaching to the triangle)
ray_without_intsct1 = [np.array([0, 0, 0]), np.array([3.666, -3.5, 30])]
# ray without intersection
ray_without_intsct2 = [np.array([0, 0, 0]), np.array([-5, -5, -20])]
print('========')
for elem in ray_with_intsct, ray_without_intsct1, ray_without_intsct2:
res, retI = check_intersection(elem, tri)
isIntsct = check_reachingIntersection(elem, tri)
print(isIntsct, res, retI)
check_intersection_3D_180512.py v0.5
check_intersection_3D_180512.py
import numpy as np
import sys
import geometry_intersection_3D_180513 as GI3
'''
v0.5 May 13, 2018
- use check_reachingIntersection() of [GI3]
v0.4 May 13, 2018
- move functions to [geometry_intersection_3D_180513.py]
v0.3 May 13, 2018
- add [RET_NOT_REACHING_PLANE]
- add check_intersection()
- rename check_intersection() to intersect3D_RayTriangle()
- add calc_norm_of_ray()
v0.2 May 12, 2018
- return [ray[0], retI] instead of [retI]
v0.1 May 12, 2018
- check_intersection()
- add get_ray()
- add calc_ray_end_point()
'''
def calc_ray_end_point(tv1, tv2, tv3, expand):
# ray directed to the triangular plane's center
# expand > 1.0 : intersecting ray
# expand < 1.0 : non-intersecting ray
return expand * (tv1 + tv2 + tv3) / 3.0
def get_ray(rst, tv1, tv2, tv3, expand):
red = calc_ray_end_point(tv1, tv2, tv3, expand)
return (RST, red)
# Triangular vertices
TV1 = np.array([-42, -8, 60])
TV2 = np.array([43, 41, 60])
TV3 = np.array([21, -54, 60])
tri = [TV1, TV2, TV3]
# Ray starting point
RST = np.array([0, 0, 0])
# rays (intersecting | non-intersecting)
ray_with_intsct = get_ray(RST, TV1, TV2, TV3, expand=1.1)
ray_without_intsct1 = get_ray(RST, TV1, TV2, TV3, expand=0.5)
ray_without_intsct2 = (RST, np.array([-5, -5, -20]))
print(tri)
print(ray_with_intsct)
print(ray_without_intsct1)
print('===')
for elem in ray_with_intsct, ray_without_intsct1, ray_without_intsct2:
res, retI = GI3.check_intersection(elem, tri)
isIntsct = GI3.check_reachingIntersection(elem, tri)
print(isIntsct)
print(res)
print(retI)
print('---')
# same as above
GI3.Test_reachingIntersection()
実行例
run
$ python3 check_intersection_3D_180512.py
[array([-42, -8, 60]), array([43, 41, 60]), array([ 21, -54, 60])]
(array([0, 0, 0]), array([ 8.06666667, -7.7 , 66. ]))
(array([0, 0, 0]), array([ 3.66666667, -3.5 , 30. ]))
===
True
2
[array([0, 0, 0]), array([ 7.33333333, -7. , 60. ])]
---
False
3
[]
---
False
0
[]
---
========
True 2 [array([0, 0, 0]), array([ 7.33327273, -7. , 60. ])]
False 3 []
False 0 []
- ray_with_intsct
- >> 交差と判定 (2 = RET_IN_SAME_PLANE)
- ray_without_intsct1
- >> 非交差と判定 (3 = RET_NOT_REACHING_PLANE)
- ray_without_intsct2
- >> 非交差と判定 (0 = RET_DISJOINT)