Merge pull request #2524 from SCIInstitute/amorris/2523-group-pvalues

Fix multi-domain group p-value display cross-domain bleed

Author: akenmorris

Libs/Analyze/Shape.cpp

@@ -819,8 +819,24 @@ void Shape::set_point_features(std::string feature, Eigen::VectorXd values) {
   auto group = get_meshes(DisplayMode::Reconstructed);
 
   if (group.valid()) {
-    for (auto mesh : group.meshes()) {
-      mesh->interpolate_scalars_to_mesh(feature, get_local_correspondence_points(), values);
+    auto local_particles = particles_.get_local_particles();
+    int num_domains = local_particles.size();
+
+    if (num_domains > 1) {
+      // Split particles and values per domain so that each mesh only interpolates
+      // from its own domain's particles, preventing cross-domain bleed in the KD-tree
+      // lookup (e.g. pelvis p-values bleeding into femur mesh near the hip joint).
+      int value_offset = 0;
+      for (int d = 0; d < num_domains && d < group.meshes().size(); d++) {
+        int num_domain_particles = local_particles[d].size() / 3;
+        Eigen::VectorXd domain_values = values.segment(value_offset, num_domain_particles);
+        group.meshes()[d]->interpolate_scalars_to_mesh(feature, local_particles[d], domain_values);
+        value_offset += num_domain_particles;
+      }
+    } else {
+      for (auto mesh : group.meshes()) {
+        mesh->interpolate_scalars_to_mesh(feature, get_local_correspondence_points(), values);
+      }
     }
   }
 }

Python/shapeworks/shapeworks/stats.py

@@ -16,7 +16,7 @@
 '''
 
 
-def compute_pvalues_for_group_difference(data, number_of_particles, permutations=100):
+def compute_pvalues_for_group_difference(data, number_of_particles, permutations=100, seed=42):
     group_id = data["group_ids"].unique()
     group_0_filenames = list(data[data["group_ids"] == group_id[0]]["filename"])
     group_1_filenames = list(data[data["group_ids"] == group_id[1]]["filename"])
@@ -26,22 +26,24 @@ def compute_pvalues_for_group_difference(data, number_of_particles, permutations
     group_1_data = sw.ParticleSystem(group_1_filenames).Particles()
     group_1_data = np.reshape(group_1_data, (number_of_particles, 3, -1))
 
-    return compute_pvalues_for_group_difference_data(group_0_data, group_1_data, permutations)
+    return compute_pvalues_for_group_difference_data(group_0_data, group_1_data, permutations, seed)
 
 
-def compute_pvalues_for_group_difference_studio(group_0_data, group_1_data, permutations=100):
+def compute_pvalues_for_group_difference_studio(group_0_data, group_1_data, permutations=100, seed=42):
     number_of_particles = int(group_0_data.shape[0] / 3)
     group_0 = np.reshape(group_0_data, (number_of_particles, 3, -1))
     group_1 = np.reshape(group_1_data, (number_of_particles, 3, -1))
-    return compute_pvalues_for_group_difference_data(group_0, group_1, permutations)
+    return compute_pvalues_for_group_difference_data(group_0, group_1, permutations, seed)
 
 
-def compute_pvalues_for_group_difference_data(group_0_data, group_1_data, permutations=100):
+def compute_pvalues_for_group_difference_data(group_0_data, group_1_data, permutations=100, seed=42):
     number_of_particles = group_0_data.shape[0]
     group_0_size = group_0_data.shape[-1]
     group_1_size = group_1_data.shape[-1]
     subset_size = min(group_0_size, group_1_size)
 
+    rng = np.random.RandomState(seed)
+
     pvalues_matrix = np.zeros((number_of_particles, permutations))
     idx = 0
     for p in range(permutations):
@@ -51,8 +53,8 @@ def compute_pvalues_for_group_difference_data(group_0_data, group_1_data, permut
             sw_message("Aborted")
             return
 
-        group_0_index = np.random.choice(group_0_size, subset_size, replace=False)
-        group_1_index = np.random.choice(group_1_size, subset_size, replace=False)
+        group_0_index = rng.choice(group_0_size, subset_size, replace=False)
+        group_1_index = rng.choice(group_1_size, subset_size, replace=False)
 
         group_0_subset = group_0_data[:, :, group_0_index]
         group_1_subset = group_1_data[:, :, group_1_index]