geometry (3D) > rayとtriangleの交差判定 > check_intersection_3D_180512.py > v0.3: triangleへ達しているか判定追加 | v0.4:処理の[geometry_intersection_3D_180513.py] v0.1への分離

```txt:動作環境
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
```

v0.1: [geometry (3D) > rayとtriangleの交差判定 > check_intersection_3D_180512.py > v0.1 > C: intersect3D_RayTriangle()からの移植, v0.2: 戻り値の変更](https://qiita.com/7of9/items/83fc9517c01c040f06e4)

### 処理概要

- rayとtriangleの交差判定を3Dで行う


### 対象データ

[対象データ > 座標関係](https://qiita.com/7of9/items/83fc9517c01c040f06e4#%E5%AF%BE%E8%B1%A1%E3%83%87%E3%83%BC%E3%82%BF--%E5%BA%A7%E6%A8%99%E9%96%A2%E4%BF%82)
に記載の３つのデータ。

### code v0.3

- rayがtriangleに達しているかの判定を加えた
- [RET_NOT_REACHING_PLANE]を追加

```check_intersection_3D_180512.py
import numpy as np
import sys

'''
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()
'''

# return values
RET_DEGENERATE = -1
RET_DISJOINT = 0
RET_INTESECT_UNIQ_POINT = 1
RET_IN_SAME_PLANE = 2
RET_NOT_REACHING_PLANE = 3


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


# 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)

```

```txt: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.        ])]
3
[]
0
[]
```

- ray_with_intsct
    - >> 交差と判定 (2 = RET_IN_SAME_PLANE)
- ray_without_intsct1
    - >> 非交差と判定 (3 = RET_NOT_REACHING_PLANE)
- ray_without_intsct2
    - >> 非交差と判定 (0 = RET_DISJOINT)
    

### code v0.4

- 処理を[geometry_intersection_3D_180513.py]へ分離した

```geometry_intersection_3D_180513.py
import numpy as np
import sys

'''
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

```


```check_intersection_3D_180512.py
import numpy as np
import sys
import geometry_intersection_3D_180513 as 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('---')

res, retI = GI3.check_intersection(ray_with_intsct, tri)
print(res)
print(retI)
res, retI = GI3.check_intersection(ray_without_intsct1, tri)
print(res)
print(retI)
res, retI = GI3.check_intersection(ray_without_intsct2, tri)
print(res)
print(retI)

```