2360 - Fix ellipsoid_multiple_domain_mesh example not working Use Case (#2362)

* Fix #2360 - ellipsoid_multiple_domain_mesh example not working

* Change order of alignment so that global doesn't use the local alignment.

* Fix MLPCA when not using only X,Y,Z

Author: akenmorris

Examples/Python/ellipsoid_multiple_domain_mesh.py

@@ -38,9 +38,7 @@ def Run_Pipeline(args):
     else:
         dataset_name = "ellipsoid_multiple_domain_mesh"
         sw.download_dataset(dataset_name, output_directory)
-        dataset_name = "ellipsoid_joint_rotation"
-        mesh_files = sorted(glob.glob(output_directory +
-                                     dataset_name + "/meshes/*.vtk"))
+        mesh_files = sorted(glob.glob(output_directory + "/meshes/*.vtk"))
 
         if args.use_subsample:
             inputMeshes =[sw.Mesh(filename) for filename in mesh_files]

Libs/Groom/Groom.cpp

@@ -529,58 +529,9 @@ bool Groom::run_alignment() {
 
   // per-domain alignment
   for (size_t domain = 0; domain < num_domains; domain++) {
-    if (abort_) {
-      return false;
-    }
-
     auto params = GroomParameters(project_, project_->get_domain_names()[domain]);
-    if (params.get_alignment_enabled()) {
-      any_alignment = true;
-    }
-
     if (params.get_use_icp()) {
       global_icp = true;
-      std::vector<Mesh> reference_meshes;
-      std::vector<Mesh> meshes;
-      for (size_t i = 0; i < subjects.size(); i++) {
-        if (!subjects[i]->is_excluded()) {
-          Mesh mesh = get_mesh(i, domain, true);
-          // if fixed subjects are present, only add the fixed subjects
-          if (subjects[i]->is_fixed() || !project_->get_fixed_subjects_present()) {
-            reference_meshes.push_back(mesh);
-          }
-          meshes.push_back(mesh);
-        }
-      }
-
-      reference_index = params.get_alignment_reference();
-      subset_size = params.get_alignment_subset_size();
-
-      Mesh reference_mesh = vtkSmartPointer<vtkPolyData>::New();
-      if (reference_index < 0 || reference_index >= subjects.size()) {
-        reference_index = MeshUtils::findReferenceMesh(reference_meshes, subset_size);
-        reference_index = reference_meshes[reference_index].get_id();
-      }
-      reference_mesh = get_mesh(reference_index, domain, true);
-
-      params.set_alignment_reference_chosen(reference_index);
-      params.save_to_project();
-
-      auto transforms = Groom::get_icp_transforms(meshes, reference_mesh);
-      assign_transforms(transforms, domain);
-    } else if (params.get_use_landmarks()) {
-      global_landmarks = true;
-      std::vector<vtkSmartPointer<vtkPoints>> landmarks;
-      for (size_t i = 0; i < subjects.size(); i++) {
-        landmarks.push_back(get_landmarks(i, domain));
-      }
-
-      int reference_index = Groom::find_reference_landmarks(landmarks);
-      params.set_alignment_reference_chosen(reference_index);
-      params.save_to_project();
-
-      auto transforms = Groom::get_landmark_transforms(landmarks, reference_index);
-      assign_transforms(transforms, domain);
     }
   }
 
@@ -652,6 +603,63 @@ bool Groom::run_alignment() {
     }
   }
 
+  // per-domain alignment
+  for (size_t domain = 0; domain < num_domains; domain++) {
+    if (abort_) {
+      return false;
+    }
+
+    auto params = GroomParameters(project_, project_->get_domain_names()[domain]);
+    if (params.get_alignment_enabled()) {
+      any_alignment = true;
+    }
+
+    if (params.get_use_icp()) {
+      global_icp = true;
+      std::vector<Mesh> reference_meshes;
+      std::vector<Mesh> meshes;
+      for (size_t i = 0; i < subjects.size(); i++) {
+        if (!subjects[i]->is_excluded()) {
+          Mesh mesh = get_mesh(i, domain, true);
+          // if fixed subjects are present, only add the fixed subjects
+          if (subjects[i]->is_fixed() || !project_->get_fixed_subjects_present()) {
+            reference_meshes.push_back(mesh);
+          }
+          meshes.push_back(mesh);
+        }
+      }
+
+      reference_index = params.get_alignment_reference();
+      subset_size = params.get_alignment_subset_size();
+
+      Mesh reference_mesh = vtkSmartPointer<vtkPolyData>::New();
+      if (reference_index < 0 || reference_index >= subjects.size()) {
+        reference_index = MeshUtils::findReferenceMesh(reference_meshes, subset_size);
+        reference_index = reference_meshes[reference_index].get_id();
+      }
+      reference_mesh = get_mesh(reference_index, domain, true);
+
+      params.set_alignment_reference_chosen(reference_index);
+      params.save_to_project();
+
+      auto transforms = Groom::get_icp_transforms(meshes, reference_mesh);
+      assign_transforms(transforms, domain);
+    } else if (params.get_use_landmarks()) {
+      global_landmarks = true;
+      std::vector<vtkSmartPointer<vtkPoints>> landmarks;
+      for (size_t i = 0; i < subjects.size(); i++) {
+        landmarks.push_back(get_landmarks(i, domain));
+      }
+
+      int reference_index = Groom::find_reference_landmarks(landmarks);
+      params.set_alignment_reference_chosen(reference_index);
+      params.save_to_project();
+
+      auto transforms = Groom::get_landmark_transforms(landmarks, reference_index);
+      assign_transforms(transforms, domain);
+    }
+  }
+
   return true;
 }
 

Libs/Particles/ParticleShapeStatistics.cpp

@@ -212,12 +212,13 @@ void ParticleShapeStatistics::compute_multi_level_analysis_statistics(std::vecto
     m += (this->num_particles_array_[idx]);
   }
   super_matrix_.resize(m, n);
+
   for (unsigned int i = 0; i < points.size(); i++) {
     unsigned int p = super_matrix_.rows();
     for (unsigned int j = 0; j < p; j++) {
-      super_matrix_(j, i * values_per_particle_) = points[i][j * values_per_particle_];
-      super_matrix_(j, i * values_per_particle_ + 1) = points[i][j * values_per_particle_ + 1];
-      super_matrix_(j, i * values_per_particle_ + 2) = points[i][j * values_per_particle_ + 2];
+      for (int k = 0; k < values_per_particle_; k++) {
+        super_matrix_(j, i * values_per_particle_ + k) = points[i][j * values_per_particle_ + k];
+      }
     }
   }
 
@@ -263,9 +264,9 @@ void ParticleShapeStatistics::compute_multi_level_analysis_statistics(std::vecto
   for (unsigned int i = 0; i < num_samples_; i++) {
     unsigned int p = m;
     for (unsigned int j = 0; j < p; j++) {
-      z_shape_dev_objective(j * values_per_particle_, i) = z_shape_dev_centered(j, i * values_per_particle_);
-      z_shape_dev_objective(j * values_per_particle_ + 1, i) = z_shape_dev_centered(j, i * values_per_particle_ + 1);
-      z_shape_dev_objective(j * values_per_particle_ + 2, i) = z_shape_dev_centered(j, i * values_per_particle_ + 2);
+      for (int k = 0; k < values_per_particle_; k++) {
+        z_shape_dev_objective(j * values_per_particle_ + k, i) = z_shape_dev_centered(j, i * values_per_particle_ + k);
+      }
     }
   }
   points_minus_mean_shape_dev_.resize(M, num_samples_);
@@ -277,9 +278,9 @@ void ParticleShapeStatistics::compute_multi_level_analysis_statistics(std::vecto
   z_rel_pose_objective.resize(domains_per_shape_ * values_per_particle_, num_samples_);
   for (unsigned int k = 0; k < domains_per_shape_; k++) {
     for (unsigned int i = 0; i < num_samples_; i++) {
-      z_rel_pose_objective(k * values_per_particle_, i) = z_rel_pose_centered(k, i * values_per_particle_);
-      z_rel_pose_objective(k * values_per_particle_ + 1, i) = z_rel_pose_centered(k, i * values_per_particle_ + 1);
-      z_rel_pose_objective(k * values_per_particle_ + 2, i) = z_rel_pose_centered(k, i * values_per_particle_ + 2);
+      for (int j = 0; j < values_per_particle_; j++) {
+        z_rel_pose_objective(k * values_per_particle_ + j, i) = z_rel_pose_centered(k, i * values_per_particle_ + j);
+      }
     }
   }
   points_minus_mean_rel_pose_.resize(domains_per_shape_ * values_per_particle_, num_samples_);