Missed intersections (embree-related)

tomas16 commented 7 months ago

Summary

I encountered a case where rays seem to pass through solid geometry. 2 of the rays are at fairly glancing angles, but 3 others aren't. In all 5 cases, the rays seem to hit the edge between 2 triangles.

System configuration

System information:

OS: 14.2 CPU: b'Intel(R) Core(TM) i9-9880H CPU @ 2.30GHz' Python: 3.11.8 | packaged by conda-forge | (main, Feb 16 2024, 20:51:20) [Clang 16.0.6 ] LLVM: 16.0.6

Dr.Jit: 0.4.4 Mitsuba: 3.5.0 Is custom build? False Compiled with: AppleClang 14.0.0.14000029 Variants: scalar_rgb scalar_spectral llvm_ad_rgb

Description

I've extracted the relevant geometry from the original mesh and created a standalone example. In the example, I compare to:

trimesh.ray.ray_pyembree.RayMeshIntersector: trimesh uses Embree version 2. (I refer to this as EmbreeIntersector)
trimesh.ray.ray_triangle.RayMeshIntersector: a python fallback, similar to mitsuba without Embree. (I refer to this as PyIntersector)

Coloring in the plots:

red = mitsuba
blue = trimesh EmbreeIntersector
green = trimesh PyIntersector

Screenshot 2024-04-25 at 4 45 37 PM

In above screenshot, the left and middle geometry are parts of a sphere. According to mitsuba, the rays hit on the inside. Comparing to trimesh, the reason is they first pass through a different triangle at a shallow angle.

On the right, there are 3 rays passing through a first surface and then hitting a surface behind it. trimesh correctly finds the intersections with the first surface.

Note that even the EmbreeIntersector in trimesh gets this right. But trimesh uses Embree 2, vs mitsuba's Embree 3.

Some close-ups:

Screenshot 2024-04-25 at 3 34 43 PM

Screenshot 2024-04-25 at 3 35 59 PM

I tried to narrow down the root cause, so I compiled custom versions of mitsuba:

no embree: only the ray on the left still intersects the wrong triangle. This seems to be the most challenging case. Since the ray hits an edge, in some cases the prim_index will be different from trimesh output, essentially the other triangle that shares the edge. (which is fine)
no embree and double precision: everything matches trimesh PyIntersector results.

Another thing I observed is that in all these cases, the rays and the edges of the triangles they're intersecting form a plane. And moreover, these are exactly the XZ and YZ planes...

To investigate this further, I added a rotation around the Z axis of both the rays and mesh. Results were mixed:

1 or 2 deg rotation: 1 of the 3 rays on the right still misses the first intersection. The other 4 cases are correct.
3 deg rotation: the ray on the left still hits the inside of the sphere. The rest are correct.
5 deg rotation: everything is correct.

Steps to reproduce

import drjit as dr
import mitsuba as mi
import numpy as np
import trimesh
import vedo
from scipy.spatial.transform import Rotation
from trimesh.ray.ray_pyembree import RayMeshIntersector as EmbreeIntersector
from trimesh.ray.ray_triangle import RayMeshIntersector as PyIntersector
from typing import cast

mi.set_variant("llvm_ad_rgb")

def load_rays():
    direction_bytes = (b"u\x16d\xa4o{\xe8\xbeu\x16d\xbf\xbb\x9a\xbd>\x00\x00\x00\x00\xa4\xccm\xbf"
                       b"\x03\x11\xbb>\x00\x00\x00\x00&Mn\xbf\x96]\xb9>\x00\x00\x00\x00'\xa2n\xbf"
                       b"\xb39\xe9\xbe\x00\x00\x00\x00\xe1\xe5c\xbf")
    direction = np.frombuffer(direction_bytes, dtype=np.float32).reshape((-1, 3))
    return np.broadcast_arrays((0, 0, 20), direction)

def trace_mitsuba(meshfile, ray_origin, ray_direction, zrotation_deg=0):
    rays = mi.Ray3f(o=ray_origin, d=ray_direction)
    mmesh = mi.load_dict({
        'type': 'ply',
        'filename': meshfile,
        'face_normals': False,
        'to_world': mi.ScalarTransform4f.rotate([0, 0, 1], angle=zrotation_deg),
    })
    scene = mi.load_dict({
        'type': 'scene',
        'my_mesh': mmesh
    })
    scene = cast(mi.Scene, scene)   # enables auto-complete
    si = scene.ray_intersect(rays, mi.RayFlags.ShadingFrame, coherent=False)
    assert dr.all(si.is_valid())
    return np.asarray(si.p), np.asarray(si.prim_index)

def trace_trimesh(intersector, ray_origin, ray_direction):
    xyz, _, prim_index = intersector.intersects_location(ray_origins=ray_origin,
                                                         ray_directions=ray_direction,
                                                         multiple_hits=False)
    return xyz, prim_index

# Inputs
zrotation_deg = 0
meshfile = "submesh.ply"

# Load data
ray_origin, ray_direction = load_rays()
tmesh = trimesh.load(meshfile)

# Rotate around Z => this is for checking whether it matters to be aligned with the XZ or YZ
# plane to trigger this issue
rotation = Rotation.from_rotvec([0, 0, zrotation_deg], degrees=True)
transform = np.zeros((4, 4))
transform[:3, :3] = rotation.as_matrix()
tmesh.apply_transform(transform)
ray_direction = rotation.apply(ray_direction)

# Run the ray tracers
mi_xyz, mi_prim_index = trace_mitsuba(meshfile, ray_origin, ray_direction, zrotation_deg)
em_xyz, em_prim_index = trace_trimesh(EmbreeIntersector(tmesh), ray_origin, ray_direction)
py_xyz, py_prim_index = trace_trimesh(PyIntersector(tmesh), ray_origin, ray_direction)

# Plot results
cases = ((em_xyz, em_prim_index, "EmbreeIntersector", (0, 0, 255)),
         (py_xyz, py_prim_index, "PyIntersector", (0, 200, 0)))
vedo.settings.immediate_rendering = False
plot = vedo.Plotter(N=2, size=(2000, 1000))
lines = np.concatenate(np.broadcast_arrays(ray_origin[0, :], mi_xyz[:, None, :]), axis=1)
for ind, (tm_xyz, tm_prim_index, intersector_name, color) in enumerate(cases):
    plot.at(ind)
    vmesh = vedo.Mesh(tmesh).c('grey', alpha=0.5).linewidth(1)
    vmesh.cellcolors[mi_prim_index] = [255, 0, 0, 128]
    vmesh.cellcolors[tm_prim_index] = color + (128,)
    plot += vmesh
    plot += vedo.Points(mi_xyz, c='red', r=6)
    plot += vedo.Points(tm_xyz, c=color, r=6)
    plot += vedo.Lines(lines, lw=2)
    plot += vedo.Text2D("trimesh " + intersector_name, pos="top-center")
    plot.show(axes=True)
plot.interactive()

submesh.ply.zip

njroussel commented 7 months ago

Hi @tomas16

I believe Embree 2, by default, does robust ray intersections whereas Embree turned that into a feature flag that Mitsuba does not enable.

To enable robust ray intersections in Mitsuba, you need to replace this line with: rtcSetSceneFlags(s.accel, RTC_SCENE_FLAG_ROBUST); and then recompile the project.

I haven't tried it myself on your reproducer yet, could you give this a go?

tomas16 commented 7 months ago

Hi @njroussel, thanks. Setting that flag solves the issue.

I've been building custom versions of mitsuba by invoking cmake and ninja, but that doesn't build a fully finished python package. If I build through python -m build, how can I specify which variants to build, or any additional cmake flags, like e.g. MI_ENABLE_EMBREE? So far I've been mixing and matching and manually copying files, but I might want to build a customized version of mitsuba to host on our company-internal pypi server.

merlinND commented 6 months ago

@tomas16 this sounds related to #1147, I suggest continuing the install / package-related discussion over there.

njroussel commented 5 months ago

Hi @tomas16

I'll close this issue, we've resolved the original problem.

I've put this on our internal feature tracker - it would be practical to have this a feature flag one can enable at runtime. If anyone comes across this, we'd welcome a PR for it - it should be fairly simple.

dvicini commented 2 months ago

Addressed in https://github.com/mitsuba-renderer/mitsuba3/pull/1300

mitsuba-renderer / mitsuba3