intersection fails for specific test case

[mahmouddinar]

Hello and thank you for the work you have done with Trimesh.

I am working on a function creating voxels by finding an intersection of cuboids with a mesh. I copy the function below and attach the failing test case. Earlier I had an issue with intersections not being found when the precision of vertices was too high. Setting "process=False" in creating the trimesh objects (cuboids) solved that issue, but I wonder if it creates another one. The test part is not too small (extents: [2.34324074, 4.51999998, 2.34349799]). But it has skinny sections. Had to attach as zip; there is another case that works without this issue.
test.zip

Below is the code:

`import trimesh
import numpy as np
import pyvista as pv
import math
import pickle
from IPython.display import Image, display
import os
import subprocess
from collections import deque
import itertools
from decimal import Decimal, getcontext
import json
import plotly.graph_objects as go
import binvox_rw
def geometric_progression_sum(p, k):
if p == 1:
return k
return (1 - p**k) / (1 - p)
def exp_voxel_vf(mesh, progression, resolution, stock_mesh=None, for_channels=False):
segments = resolution//2
if stock_mesh is not None:
bounding_box_center = stock_mesh.bounding_box.centroid
bb_x, bb_y, bb_z = stock_mesh.bounding_box.extents
else:
bounding_box_center = mesh.bounding_box.centroid
bb_x, bb_y, bb_z = mesh.bounding_box.extents

ss_x, ss_y, ss_z = bb_x / geometric_progression_sum(progression, segments)/2, \
                bb_y / geometric_progression_sum(progression, segments)/2, \
                bb_z / geometric_progression_sum(progression, segments)/2

seq_x = [sum(ss_x * progression ** (segments - 1 - j) for j in range(i+1)) for i in range(segments)]
seq_y = [sum(ss_y * progression ** (segments - 1 - j) for j in range(i+1)) for i in range(segments)]
seq_z = [sum(ss_z * progression ** (segments - 1 - j) for j in range(i+1)) for i in range(segments)]
x_neg_to_pos = np.concatenate([np.flip(bounding_box_center[0] - seq_x), np.array([bounding_box_center[0]]), seq_x + bounding_box_center[0]])
y_neg_to_pos = np.concatenate([np.flip(bounding_box_center[1] - seq_y), np.array([bounding_box_center[1]]), seq_y + bounding_box_center[1]])
z_neg_to_pos = np.concatenate([np.flip(bounding_box_center[2] - seq_z), np.array([bounding_box_center[2]]), seq_z + bounding_box_center[2]])

faces = np.array([
    [0, 4, 5, 1],  # front
    [2, 3, 7, 6],  # back
    [0, 2, 6, 4],  # bottom
    [1, 5, 7, 3],  # top
    [0, 1, 3, 2],  # left
    [4, 6, 7, 5],  # right
])
voxel_vf_grid = np.zeros((resolution, resolution, resolution))
cuboids = []

not_volume_count = 0
for i in range(resolution):
    for j in range(resolution):
        for k in range(resolution):
            x_levels = [x_neg_to_pos[i], x_neg_to_pos[i+1]]
            y_levels = [y_neg_to_pos[j], y_neg_to_pos[j+1]]
            z_levels = [z_neg_to_pos[k], z_neg_to_pos[k+1]]
           
            vertices =  np.array(list(itertools.product(x_levels, y_levels, z_levels)))
           
            cuboid = trimesh.Trimesh(vertices=vertices, faces=faces, process=False)
            cuboids.append(cuboid)
           
            try:
                intersection = mesh.intersection(cuboid)
                # difference = cuboid.difference(intersection)
                if not intersection.is_volume or not cuboid.is_volume:
                    # raise ValueError(f"Intersection at index {i,j,k} is not a volume!")
                    # err_cuboids.append(([i, j, k], "Intersection is not a volume"))
                    voxel_vf_grid[i, j, k] = 0
                    not_volume_count += 1
                    continue
                voxel_vf_grid[i, j, k] = intersection.volume/cuboid.volume
           
            except Exception as e:
             
                print(f"{mesh.metadata['file_name'].split('.')[0]}_Error at index {i,j,k}: {e}")

if for_channels:
    return voxel_vf_grid, cuboids
relative_dir = os.path.relpath(os.path.dirname(mesh.metadata['file_path']), start=os.getcwd())
mesh_base = os.path.splitext(os.path.basename(mesh.metadata['file_path']))[0]


np.save(relative_dir+"/"+mesh_base+f"_vg_{progression}_{resolution}", voxel_vf_grid)

with open(relative_dir+"/"+mesh_base+f"_cb_{progression}_{resolution}", 'wb') as f:
    pickle.dump(cuboids, f)
print(f"Count of invalid cuboids!': {not_volume_count}")
return voxel_vf_grid, cuboids

mesh = trimesh.load_mesh('test/collet1.stl')
exp_voxel_vf(mesh, 2, 8) `