Fix assignment of compartment grid points (#192)

* set surface distances in Gradient.py

* Change grid point size to be relative to the grid spacing

* sort grid points by stored values

* ignore nan values while scaling distances

* set `scale_to_next_surface` as an option

* do not calculate scaled surface distances if not needed

* only reassign select grid points

* * remove nans and sort grid values
* user relative size for grid points in simularium

* reduce chunk size for mesh point calculations

* update peroxisome packing config

* formatting

cellpack/autopack/Compartment.py

@@ -1077,7 +1077,9 @@ def prepare_buildgrid_box(self, env):
             len(env.grid.masterGridPositions),
         )
 
-    def set_surface_distances(self, env, master_grid_positions):
+    def set_surface_distances(
+        self, env, master_grid_positions, calc_distance_between_surfaces=False
+    ):
         surface_mask = numpy.equal(self.number, env.grid.compartment_ids)
         surface_ids = numpy.nonzero(surface_mask)
         surface_positions = master_grid_positions[surface_ids]
@@ -1094,7 +1096,11 @@ def set_surface_distances(self, env, master_grid_positions):
         all_surface_distances = numpy.full(master_grid_positions.shape[0], numpy.nan)
         all_surface_distances[grid_pt_indexes] = surface_distances
 
-        if self.parent is not None and parent_id != 0:
+        if (
+            calc_distance_between_surfaces
+            and self.parent is not None
+            and parent_id != 0
+        ):
             grid_pts_between_surfaces = numpy.equal(
                 env.grid.compartment_ids, -parent_id
             )
@@ -1319,7 +1325,10 @@ def BuildGrid_trimesh(
         self.log.info(f"GOT POINTS IN SPHERE {len(points_in_encap_sphere)}")
 
         point_compartment_ids = compartment_ids[points_in_encap_sphere]
-        point_ids_to_assign = points_in_encap_sphere[point_compartment_ids == 0]
+        # largest compartments need to be created first for this to work
+        point_ids_to_assign = points_in_encap_sphere[
+            numpy.abs(point_compartment_ids) < number
+        ]
         point_positions = numpy.float16(master_grid_positions[point_ids_to_assign])
 
         # check surface points

cellpack/autopack/Environment.py

@@ -1194,7 +1194,6 @@ def BuildCompartmentsGrids(self):
             )  # return inside and surface point
             aInteriorGrids.append(points_inside_compartments)
             aSurfaceGrids.append(points_on_compartment_surfaces)
-
         self.grid.aInteriorGrids = aInteriorGrids
         self.grid.aSurfaceGrids = aSurfaceGrids
         self.log.info("I'm out of the loop and have build my grid with inside points")
@@ -1310,7 +1309,9 @@ def buildGrid(self, rebuild=True):
                         gradient.mode_settings["object"], "surface_distances"
                     ):
                         gradient.mode_settings["object"].set_surface_distances(
-                            self, self.grid.masterGridPositions
+                            self,
+                            self.grid.masterGridPositions,
+                            gradient.mode_settings.get("scale_to_next_surface", False),
                         )
                 self.gradients[g].build_weight_map(
                     boundingBox, self.grid.masterGridPositions

cellpack/autopack/Gradient.py

@@ -117,8 +117,8 @@ def get_normalized_values(self, values):
         """
         Scale values between 0 and 1
         """
-        max_value = max(values)
-        min_value = min(values)
+        max_value = numpy.nanmax(values)
+        min_value = numpy.nanmin(values)
         return (values - min_value) / (max_value - min_value)
 
     def pickPoint(self, listPts):
@@ -175,10 +175,7 @@ def build_surface_distance_weight_map(self):
                 "object"
             ].scaled_distance_to_next_surface
         else:
-            self.distances = (
-                self.mode_settings["object"].surface_distances
-                / self.mode_settings["object"].max_distance
-            )
+            self.distances = self.mode_settings["object"].surface_distances
         self.set_weights_by_mode()
 
     def build_directional_weight_map(self, bb, master_grid_positions):

cellpack/autopack/MeshStore.py

@@ -4,7 +4,7 @@
 import trimesh
 from cellpack import autopack
 
-CHUNK_SIZE = 100000
+CHUNK_SIZE = 50000
 
 
 class MeshStore:
@@ -261,6 +261,7 @@ def contains_point(self, geomname, point):
 
     def contains_points_mesh(self, geomname, points):
         mesh = self.get_object(geomname)
+        inside = numpy.full(len(points), False)
         if mesh is not None:
             if len(points) <= CHUNK_SIZE:
                 return mesh.contains(points)  # TODO: check for memory leak
@@ -273,7 +274,7 @@ def contains_points_mesh(self, geomname, points):
                     else:
                         inside = numpy.append(inside, mesh.contains(chunk))
                 return inside
-        return numpy.full(len(points), False)
+        return inside
 
     def get_smallest_radius(self, geomname, center):
         mesh = self.get_object(geomname)

cellpack/autopack/upy/simularium/simularium_helper.py

@@ -338,10 +338,33 @@ def concatObjectMatrix(self):
     def GetAbsPosUntilRoot(self, obj):
         return [0, 0.0, 0.0]
 
-    def add_grid_data_to_scene(self, incoming_name, positions, values):
+    @staticmethod
+    def remove_nans(positions, values):
+        naninds = np.isnan(values)
+        values = values[~naninds]
+        positions = positions[~naninds]
+        return positions, values
+
+    @staticmethod
+    def sort_values(positions, values):
+        inds = np.argsort(values)
+        values = values[inds]
+        positions = positions[inds]
+        return positions, values
+
+    def add_grid_data_to_scene(self, incoming_name, positions, values, radius=0.5):
+
+        positions, values = self.remove_nans(positions, values)
+        if len(values) == 0:
+            print("no values to display")
+            return
+
+        positions, values = self.sort_values(positions, values)
+
         colormap = matplotlib.cm.Reds(values)
+
         for index, value in enumerate(values):
-            name = f"{incoming_name}#{value}"
+            name = f"{incoming_name}#{value:.3f}"
             self.display_data[name] = DisplayData(
                 name=name,
                 display_type=DISPLAY_TYPE.SPHERE,
@@ -353,7 +376,7 @@ def add_grid_data_to_scene(self, incoming_name, positions, values):
                 name,
                 None,
                 f"{incoming_name}-{index}",
-                0.5,
+                radius,
                 point_pos,
                 np.identity(4),
                 None,
@@ -457,6 +480,7 @@ def init_scene_with_objects(
         grid_point_positions=None,
         grid_point_compartment_ids=None,
         show_sphere_trees=False,
+        grid_pt_radius=0.5,
     ):
         self.time = 0
         instance_number = 0
@@ -533,7 +557,7 @@ def init_scene_with_objects(
                         name,
                         None,
                         f"{name}-{index}",
-                        0.5,
+                        grid_pt_radius,
                         point_pos,
                         np.identity(4),
                         None,

cellpack/autopack/writers/__init__.py

@@ -155,7 +155,11 @@ def save_as_simularium(self, env, all_ingr_as_array, compartments):
         grid_positions = env.grid.masterGridPositions if env.show_grid_spheres else None
         compartment_ids = env.grid.compartment_ids if env.show_grid_spheres else None
         env.helper.init_scene_with_objects(
-            all_ingr_as_array, grid_positions, compartment_ids, env.show_sphere_trees
+            objects=all_ingr_as_array,
+            grid_point_positions=grid_positions,
+            grid_point_compartment_ids=compartment_ids,
+            show_sphere_trees=env.show_sphere_trees,
+            grid_pt_radius=env.grid.gridSpacing / 4,
         )
 
         if compartments is not None:
@@ -167,10 +171,16 @@ def save_as_simularium(self, env, all_ingr_as_array, compartments):
         if grid_positions is not None and len(env.gradients):
             for _, gradient in env.gradients.items():
                 env.helper.add_grid_data_to_scene(
-                    f"{gradient.name}-distances", grid_positions, gradient.distances
+                    f"{gradient.name}-distances",
+                    grid_positions,
+                    gradient.distances,
+                    env.grid.gridSpacing / 4,
                 )
                 env.helper.add_grid_data_to_scene(
-                    f"{gradient.name}-weights", grid_positions, gradient.weight
+                    f"{gradient.name}-weights",
+                    grid_positions,
+                    gradient.weight,
+                    env.grid.gridSpacing / 4,
                 )
         file_name = env.helper.writeToFile(
             env.result_file, env.boundingBox, env.name, env.version

examples/packing-configs/peroxisome_packing_config.json

@@ -4,15 +4,15 @@
     "inner_grid_method": "trimesh",
     "live_packing": false,
     "ordered_packing": false,
-    "number_of_packings": 10,
+    "number_of_packings": 1,
     "out": "out/analyze",
     "overwrite_place_method": true,
     "place_method": "spheresSST",
     "save_analyze_result": true,
     "show_grid_plot": false,
-    "save_plot_figures": false,
+    "save_plot_figures": true,
     "spacing": 2.5,
     "use_periodicity": false,
     "show_sphere_trees": false,
-    "load_from_grid_file": true
+    "load_from_grid_file": false
 }