LoginSignup
0
0

More than 5 years have passed since last update.

@7of9

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

Last updated at Posted at 2018-05-13
動作環境
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: 戻り値の変更

処理概要

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

対象データ

対象データ > 座標関係
に記載の３つのデータ。

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

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0