mikedh / trimesh

Python library for loading and using triangular meshes.
https://trimesh.org
MIT License
3.02k stars 583 forks source link

intersection fails for specific test case #2308

Open mahmouddinar opened 3 weeks ago

mahmouddinar commented 3 weeks ago

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

mahmouddinar commented 3 weeks ago

Another thing I found was that the collet1 test part returns False when mesh.is_volume but also has a value for mesh.volume np.float64(2.074216113966994), and it is much less than the bounding box volume How can it has volume but not being a volume? Thank you.