Search code examples
pythonmeshraycastingraytracing

Ray-tracing (or casting) between rays and STL file


I have a 3D shape (a bone) as a triangular mesh in STL format and I'm working with Python language.

I then have a ray with its origin and direction and I need to know:

  • which face of the mesh collides with this ray
  • coordinates of intersection point
  • the angle between my ray and the normal vector of the collided face

I know 'pycaster' ( https://pyscience.wordpress.com/2014/09/21/ray-casting-with-python-and-vtk-intersecting-linesrays-with-surface-meshes/ ) but it only works with Python 2.x and I need to use Python 3.x.

I know 'trimeshgeom' class ( http://cgkit.sourceforge.net/doc2/_sources/trimeshgeom.txt ) from cgkit but PyCharm can't install it as an interpreter (can't understand why).

Does anybody know how I could do what I'm trying to do?

Thanks, Cheers,

Davide


Solution

  • We had a similar task once and ended up implementing the rather simple Möller-Trumbore-Algorithm.

    The code can be stolen from printrun:

    def ray_triangle_intersection(ray_near, ray_dir, v123):
        """
        Möller–Trumbore intersection algorithm in pure python
        Based on http://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
        """
        v1, v2, v3 = v123
        eps = 0.000001
        edge1 = v2 - v1
        edge2 = v3 - v1
        pvec = numpy.cross(ray_dir, edge2)
        det = edge1.dot(pvec)
        if abs(det) < eps:
            return False, None
        inv_det = 1. / det
        tvec = ray_near - v1
        u = tvec.dot(pvec) * inv_det
        if u < 0. or u > 1.:
            return False, None
        qvec = numpy.cross(tvec, edge1)
        v = ray_dir.dot(qvec) * inv_det
        if v < 0. or u + v > 1.:
            return False, None
    
        t = edge2.dot(qvec) * inv_det
        if t < eps:
            return False, None
    
    return True, t
    

    For the angle calculation you can use e.g. angle_between_vectors() from transformations.py