Trouble getting the right alignment using python

Hello!, I’m trying to use shapeworks via the python API, but I’m having trouble with getting the transform / alignment right. My issue is similar to this topic, but the advice from there has not solved it.

I’ve set up a simple script to demonstrate my issue. My use case is roughly similar to the hip_multiple_domain example, but when I run it with some test meshes, the alignment looks wrong in shapeworks studio (see the image) and the model output is garbage.

Interestingly, when I go back to the groom menu in shapeworks studio and run with just the alignment box ticked, everything works.

Here is my code:

import shapeworks as sw
from glob import glob
from pathlib import Path
import numpy as np
import subprocess
import pyvista as pv
import os


def main():
    data_dirs = [Path("data/001"), Path("data/002")]
    for data_dir in data_dirs:
        if not os.path.exists(data_dir):
            os.makedirs(data_dir)

    # Create meshes
    # ------------------------------------------------------------------------------------------------------------------
    mesh_1_left = pv.Sphere(radius=1.5, center=[-2, 0, 0]).scale([1, 1, 1], inplace=False)
    mesh_1_right = pv.Sphere(radius=1.5, center=[2, 0, 0]).scale([1, 1, 1], inplace=False)

    mesh_2_left = pv.Sphere(radius=1.5, center=[-4, -1, 0]).scale([1, 1, 1.5], inplace=False)
    mesh_2_right = pv.Sphere(radius=1.5, center=[1, -1, 0]).scale([1, 1, 2], inplace=False)

    mesh_1_left.save(f"{data_dirs[0]}/left.vtk")
    mesh_1_right.save(f"{data_dirs[0]}/right.vtk")
    mesh_2_left.save(f"{data_dirs[1]}/left.vtk")
    mesh_2_right.save(f"{data_dirs[1]}/right.vtk")

    # Find reference
    # ------------------------------------------------------------------------------------------------------------------
    all_mesh_files = []
    for data_dir in data_dirs:
        mesh_files = glob(str(data_dir / "*"))
        all_mesh_files.append(mesh_files)

    domains_per_shape = [len(files) for files in all_mesh_files][0]  # 2 Files in each data directory

    # Find the combined reference mesh
    all_meshes = [sw.Mesh(file) for file in np.array(all_mesh_files).ravel()]
    ref_index, combined_meshes = sw.find_reference_mesh_index(all_meshes, domains_per_shape)
    combined_reference_mesh = combined_meshes[ref_index]

    # Find the individual domain meshes at the reference mesh index:
    # find_reference_mesh_index destroys all_meshes during the combining process. So we have to create them again
    all_meshes = [sw.Mesh(file) for file in np.array(all_mesh_files).ravel()]
    domain_reference_meshes = []
    for i in range(domains_per_shape):
        domain_reference_mesh = all_meshes[ref_index * domains_per_shape + i]
        domain_reference_meshes.append(domain_reference_mesh)

    # Compute transforms
    # ------------------------------------------------------------------------------------------------------------------
    all_combined_transforms = []
    all_domain_transforms = []
    for data_dir in data_dirs:
        mesh_files = glob(str(data_dir / "*"))
        meshes = [sw.Mesh(file) for file in mesh_files]

        combined_mesh = meshes[0].copy()  # Copy so we don't destroy meshes[0] for later!
        for mesh in meshes[1:]:
            combined_mesh += mesh

        combined_transform = combined_mesh.createTransform(combined_reference_mesh, sw.Mesh.AlignmentType.Rigid,
                                                           10).flatten()

        domain_transforms = []
        for mesh, ref_mesh in zip(meshes, domain_reference_meshes):
            domain_transform = mesh.createTransform(mesh, sw.Mesh.AlignmentType.Rigid, 10).flatten()
            domain_transforms.append(domain_transform)

        all_combined_transforms.append(combined_transform)
        all_domain_transforms.append(domain_transforms)

    # Run optimize
    # ------------------------------------------------------------------------------------------------------------------
    params = {
        "checkpointing_interval": 200,
        "keep_checkpoints": 0,
        "iterations_per_split": 2500,
        "optimization_iterations": 200,
        "starting_regularization": 1000,
        "ending_regularization": 10,
        "relative_weighting": 4,
        "initial_relative_weighting": 0.03,
        "save_init_splits": 0,
        "verbosity": 0,
        "use_normals": 1,
        "normals_strength": 5.0,
        "procrustes": 0,
        "procrustes_scaling": 1,
        "procrustes_rotation_translation": 1,
        "number_of_particles": [32, 32]
    }

    subjects = []
    for i, data_dir in enumerate(data_dirs):
        subject = sw.Subject()
        mesh_files = glob(str(data_dir / "*"))
        subject.set_number_of_domains(len(mesh_files))

        combined_transform = all_combined_transforms[i]
        domain_transforms = all_domain_transforms[i]
        transforms = [*domain_transforms, combined_transform]

        subject.set_groomed_transforms(transforms)
        subject.set_groomed_filenames(mesh_files)
        subject.set_original_filenames(mesh_files)
        subject.set_display_name(str(Path(data_dir).stem))
        subjects.append(subject)

    project = sw.Project()
    project.set_subjects(subjects)
    parameters = sw.Parameters()
    for key, value in params.items():
        parameters.set(key, sw.Variant(value))

    project.set_parameters("optimize", parameters)

    spreadsheet_file = "sphere_test.swproj"
    project.save(str(spreadsheet_file))

    optimize_cmd = ("shapeworks optimize --progress --name " + str(spreadsheet_file)).split()
    subprocess.check_call(optimize_cmd)

    analyze_cmd = ('ShapeWorksStudio ' + spreadsheet_file).split()
    subprocess.check_call(analyze_cmd)


if __name__ == "__main__":
    main()

Can anyone help me understand this

Problem solved! The issue was that all groomed files need to have unique filenames.

Making the following changes to my test code demsonstrates this fix:

    mesh_1_left.save(f"{data_dirs[0]}/left_001.vtk")
    mesh_1_right.save(f"{data_dirs[0]}/right_001.vtk")
    mesh_2_left.save(f"{data_dirs[1]}/left_002.vtk")
    mesh_2_right.save(f"{data_dirs[1]}/right_002.vtk")

Now it runs as expected.

By the way, to get transforms that are exactly the same as the “Center” option in studio, the following code works:

    # Compute transforms
    # ------------------------------------------------------------------------------------------------------------------
    all_combined_transforms = []
    all_domain_transforms = []
    for data_dir in data_dirs:
        mesh_files = glob(str(data_dir / "*"))
        meshes = [sw.Mesh(file) for file in mesh_files]

        combined_mesh = meshes[0].copy()  # Copy so we don't destroy meshes[0] for later!
        for mesh in meshes[1:]:
            combined_mesh += mesh

        com1 = combined_reference_mesh.centerOfMass()
        com2 = combined_mesh.centerOfMass()
        centered_mesh = combined_mesh.copy().translate(
            [-(com2[0] - com1[0]), -(com2[1] - com1[1]), -(com2[2] - com1[2])])

        combined_transform = combined_mesh.createTransform(centered_mesh, sw.Mesh.AlignmentType.Rigid,
                                                           1).flatten()

        domain_transforms = []
        for mesh, ref_mesh in zip(meshes, domain_reference_meshes):
            com = mesh.centerOfMass()
            centered_mesh = mesh.copy().translate([-com[0], -com[1], -com[2]])
            domain_transforms.append(mesh.createTransform(centered_mesh, sw.Mesh.AlignmentType.Rigid, 1).flatten())

        all_combined_transforms.append(combined_transform)
        all_domain_transforms.append(domain_transforms)`

Is this a bug that we need to fix in ShapeWorks or was it just an issue with the script you were using?

A bug in ShapeWorks I believe. In my project I use a dataset structure where files are identified by directory and I would prefer not to need to give the files globally unique names. But the workaround is ok for now.

Thanks, I’ve added an issue to track: