動作環境

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)