{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Visualize the Moeller-Trumbore Algorithm {#moeller_ray_trace_example}\n\nThis example demonstrates the Moeller-Trumbore intersection algorithm\nusing pyvista.\n\nFor additional details, please reference the following:\n\n-   [Moeller-Trumbore intersection\n    algorithm](https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm)\n-   [Fast Minimum Storage Ray Triangle\n    Intersectio](https://cadxfem.org/inf/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf)\n\nFirst, define the ray triangle intersection method.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n\nimport pyvista as pv\n\n\ndef ray_triangle_intersection(ray_start, ray_vec, triangle):\n    \"\"\"Moeller-Trumbore intersection algorithm.\n\n    Parameters\n    ----------\n    ray_start : np.ndarray\n        Length three numpy array representing start of point.\n\n    ray_vec : np.ndarray\n        Direction of the ray.\n\n    triangle : np.ndarray\n        ``3 x 3`` numpy array containing the three vertices of a\n        triangle.\n\n    Returns\n    -------\n    bool\n        ``True`` when there is an intersection.\n\n    tuple\n        Length three tuple containing the distance ``t``, and the\n        intersection in unit triangle ``u``, ``v`` coordinates.  When\n        there is no intersection, these values will be:\n        ``[np.nan, np.nan, np.nan]``\n\n    \"\"\"\n    # define a null intersection\n    null_inter = np.array([np.nan, np.nan, np.nan])\n\n    # break down triangle into the individual points\n    v1, v2, v3 = triangle\n    eps = 0.000001\n\n    # compute edges\n    edge1 = v2 - v1\n    edge2 = v3 - v1\n    pvec = np.cross(ray_vec, edge2)\n    det = edge1.dot(pvec)\n\n    if abs(det) < eps:  # no intersection\n        return False, null_inter\n    inv_det = 1.0 / det\n    tvec = ray_start - v1\n    u = tvec.dot(pvec) * inv_det\n\n    if u < 0.0 or u > 1.0:  # if not intersection\n        return False, null_inter\n\n    qvec = np.cross(tvec, edge1)\n    v = ray_vec.dot(qvec) * inv_det\n    if v < 0.0 or u + v > 1.0:  # if not intersection\n        return False, null_inter\n\n    t = edge2.dot(qvec) * inv_det\n    if t < eps:\n        return False, null_inter\n\n    return True, np.array([t, u, v])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Create a basic triangle within pyvista\npoints = np.array([[0, 0, 0], [0, 1, 0], [1, 0, 0]])\nfaces = np.array([3, 0, 1, 2])\ntri = pv.PolyData(points, faces)\n\n# cast a ray above pointed downwards\nstart = np.array([0.3, 0.25, 1])\ndirection = np.array([0, 0, -1])\n\n# compute if the intersection exists\ninter, tuv = ray_triangle_intersection(start, direction, points)\nt, u, v = tuv\n\nprint('Intersected', inter)\nprint('t:', t)\nprint('u:', u)\nprint('v:', v)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the problem setup and the intersection\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "if inter:\n    # reconstruct intersection point in barycentric coordinates.  See\n    # https://en.wikipedia.org/wiki/Barycentric_coordinate_system\n    a, b, c = (1 - u - v), u, v\n    point = tri.points[0] * a + tri.points[1] * b + tri.points[2] * c\n\n    pl = pv.Plotter()\n    pl.add_text(f'Intersected at ({point[0]:.3}, {point[0]:.3}, {point[0]:.3})', font_size=26)\n    pl.add_mesh(tri)\n    _ = pl.add_arrows(\n        np.array([start]),\n        np.array([direction]),\n        show_scalar_bar=False,\n        color='r',\n        style='wireframe',\n    )\n    pl.add_points(np.array([point]), point_size=20, render_points_as_spheres=True, color='b')\n    pl.add_point_labels(tri, [f'a = {1 - u - v:.3}', f'b = {u:.3}', f'c = {v:.3}'], font_size=40)\n    pl.show_bounds()\n    pl.camera_position = 'xy'\n    pl.show()\n\nelse:  # no intersection\n    pl = pv.Plotter()\n    pl.add_text('No intersection')\n    _ = pl.add_arrows(\n        np.array([start]),\n        np.array([direction]),\n        show_scalar_bar=False,\n        color='r',\n        style='wireframe',\n    )\n    pl.add_mesh(tri)\n\n    pl.show_bounds()\n    pl.camera_position = 'xy'\n\n    pl.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}