動作環境

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
[]