SCIInstitute / ShapeWorks

ShapeWorks
http://sciinstitute.github.io/ShapeWorks/
Other
101 stars 32 forks source link

python function find_reference_mesh_index not working. #2152

Closed acreegan closed 8 months ago

acreegan commented 8 months ago

I'm having trouble getting find_reference_mesh_index to produce the results I expect.

Code

import pyvista as pv
import tempfile
import shapeworks as sw

# Snippet to check that find_reference_mesh works how we expect it to

# Create some pyvista meshes. Center should be the medoid
left = pv.Box().translate([-4, -4, -4], inplace=False)
center = pv.Box().translate([1, 1, 1], inplace=False)
right = pv.Box().translate([3, 3, 3], inplace=False)

# Save to .vtk and read to convert to Shapeworks mesh
sw_meshes = []
for pv_mesh in [left, center, right]:
    with tempfile.TemporaryDirectory() as d:
        pv_mesh.save(f"{d}/temp_mesh.vtk")
        sw_meshes.append(sw.Mesh(f"{d}/temp_mesh.vtk"))

# Demonstrate that the Shapeworks meshes are loaded with the correct center of masses
for i, mesh in enumerate(sw_meshes):
    print(f"Mesh {i} center of mass: {mesh.centerOfMass()}")

# Try find_reference_mesh_index with different ordered mesh lists
ref_index_1 = sw.find_reference_mesh_index(sw_meshes, domains_per_shape=1)
ref_index_2 = sw.find_reference_mesh_index([sw_meshes[1], sw_meshes[2], sw_meshes[0]], domains_per_shape=1)
ref_index_3 = sw.find_reference_mesh_index([sw_meshes[2], sw_meshes[0], sw_meshes[1]], domains_per_shape=1)

# Print what we found. We expect the indices to be 1, 0, 2
for i, ref_index in enumerate([ref_index_1, ref_index_2, ref_index_3]):
    print(f"Ref index {i}: {ref_index}")

Output

Mesh 0 center of mass: [-4. -4. -4.]
Mesh 1 center of mass: [1. 1. 1.]
Mesh 2 center of mass: [3. 3. 3.]
Ref index 0: 0
Ref index 1: 0
Ref index 2: 0

Expected Output

Mesh 0 center of mass: [-4. -4. -4.]
Mesh 1 center of mass: [1. 1. 1.]
Mesh 2 center of mass: [3. 3. 3.]
Ref index 0: 1
Ref index 1: 0
Ref index 2: 2
akenmorris commented 8 months ago

find_reference_mesh_index (and findReferenceMesh, etc) use pairwise ICP registrations to try to find the shape with the minimal surface to surface distance after registration. The goal is to find something of a median shape such that aligning all others to the reference shape provides a good alignment and starting point for correspondence optimization.

It appears that all of your shapes are the same, just translated. The ICP registration will start with centering on the center of mass. In this snippet, the surface to surface distances will all be zero and any shape can be the reference. Likely it will return the first one (index 0) just because it's first in the list of minimum distances.

acreegan commented 8 months ago

Got it, thanks for clearing that up.

Below is a fixed snippet showing a correct interpretation of the function. Instead of translating I have scaled the meshes.

import pyvista as pv
import tempfile
import shapeworks as sw

# Snippet to check that find_reference_mesh works how we expect it to

# Create some pyvista meshes. Center should be the medoid
small = pv.Box().scale([100, 100, 100], inplace=False)
middle = pv.Box().scale([101, 101, 101], inplace=False)
large = pv.Box().scale([102, 102, 102], inplace=False)

# Save to .vtk and read to convert to Shapeworks mesh
sw_meshes = []
for pv_mesh in [small, middle, large]:
    with tempfile.TemporaryDirectory() as d:
        pv_mesh.save(f"{d}/temp_mesh.vtk")
        sw_meshes.append(sw.Mesh(f"{d}/temp_mesh.vtk"))

# Try find_reference_mesh_index with different ordered mesh lists
ref_index_1 = sw.find_reference_mesh_index(sw_meshes, domains_per_shape=1)
ref_index_2 = sw.find_reference_mesh_index([sw_meshes[1], sw_meshes[2], sw_meshes[0]], domains_per_shape=1)
ref_index_3 = sw.find_reference_mesh_index([sw_meshes[2], sw_meshes[0], sw_meshes[1]], domains_per_shape=1)

# Print what we found. We expect the indices to be 1, 0, 2
for i, ref_index in enumerate([ref_index_1, ref_index_2, ref_index_3]):
    print(f"Ref index {i}: {ref_index}")