動作環境
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
処理概要
- rayとtriangleの交差判定を3Dで行う
参考
上記の実装(C)をNumpy実装にしてみた。
対象データ > 座標関係
以下の３種類のデータを使ってみる。
各データはTinkercadで可視化しているので、座標関係はそちらで確認ください。
(triangleの頂点:赤球、ray:青球からオレンジ球方向)
- A. 交差する座標 (ray_with_intsct)
- triangle方向を向いているray
- triangleまで達するray
- https://www.tinkercad.com/things/j1XHIjnvZV7
- B. 交差しない座標 (ray_without_intsct1)
- triangle方向を向いているray
- triangleまで達しないray
- https://www.tinkercad.com/things/h8BoX9qKbY8
- C. 交差しない座標 (ray_without_intsct2)
- triangle方向を向いていないray
- https://www.tinkercad.com/things/3o8ogNnN5o8
code v0.1
check_intersection_3D_180512.py
import numpy as np
import sys
'''
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)
def check_intersection(ray, tri):
# ported from
# http://geomalgorithms.com/a06-_intersect-2.html#intersect3D_RayTriangle()
# return values
RET_DEGENERATE = -1
RET_DISJOINT = 0
RET_INTESECT_UNIQ_POINT = 1
RET_IN_SAME_PLANE = 2
retI = 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, retI # 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, retI
else:
return RET_DISJOINT, retI # 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, retI
# 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, retI
valT = (uv * wu - uu * wv) / valD
if (valT < 0.0 or (valS + valT) > 1.0): # retI is outside tri
return RET_DISJOINT, retI
return RET_IN_SAME_PLANE, retI # retI is in T
# 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('---')
res, retI = check_intersection(ray_with_intsct, tri)
print(res)
print(retI)
res, retI = check_intersection(ray_without_intsct1, tri)
print(res)
print(retI)
res, retI = check_intersection(ray_without_intsct2, tri)
print(res)
print(retI)
実行結果
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. ]))
---
2
[ 7.33333333 -7. 60. ]
2
[ 7.33333333 -7. 60. ]
0
[]
- ray_with_intsct
- >> 交差と判定 (2 = RET_IN_SAME_PLANE)
- ray_without_intsct1
- >> 交差と判定 (2)
- ray_without_intsct2
- >> 非交差と判定 (0 = RET_DISJOINT)
TODO
- ray_without_intsct1 の非交差判定
- retIとrayの長さの比較
code v0.2
- 戻り値にrayの始点を追加した
check_intersection_3D_180512.py
import numpy as np
import sys
'''
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)
def check_intersection(ray, tri):
# ported from
# http://geomalgorithms.com/a06-_intersect-2.html#intersect3D_RayTriangle()
# return values
RET_DEGENERATE = -1
RET_DISJOINT = 0
RET_INTESECT_UNIQ_POINT = 1
RET_IN_SAME_PLANE = 2
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
# 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('---')
res, retI = check_intersection(ray_with_intsct, tri)
print(res)
print(retI)
res, retI = check_intersection(ray_without_intsct1, tri)
print(res)
print(retI)
res, retI = check_intersection(ray_without_intsct2, tri)
print(res)
print(retI)
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. ]))
---
2
[array([0, 0, 0]), array([ 7.33333333, -7. , 60. ])]
2
[array([0, 0, 0]), array([ 7.33333333, -7. , 60. ])]
0
[]