Trouble getting the right alignment using python

1

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.

Untitled
Untitled673×592 62 KB

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

2

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

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

4

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.

5

Thanks, I’ve added an issue to track: