Feat: Add boundary thickness mask

src/synthesis_workflow/circuit.py

@@ -1,55 +1,47 @@
 """Functions for slicing mvd3 circuit files to place specific cells only."""
 import logging
+from pathlib import Path
 
 import numpy as np
 import pandas as pd
 import yaml
 from brainbuilder.app.cells import _place as place
+from neurocollage.mesh_helper import MeshHelper
 from neurocollage.planes import get_cells_between_planes
 from neurocollage.planes import slice_n_cells
 from neurocollage.planes import slice_per_mtype
 from region_grower.atlas_helper import AtlasHelper
+from trimesh.voxel import VoxelGrid
 from voxcell import CellCollection
 from voxcell.nexus.voxelbrain import LocalAtlas
 
 L = logging.getLogger(__name__)
 
 
-def create_atlas_thickness_mask(atlas_dir):
-    """Create a mask on an atlas to to select voxels with small enough thicknesses values.
+def create_boundary_mask(atlas_dir, region, boundary_thickness=10):
+    """Create a mask on an atlas to to select voxels away from boundary with some voxel distance."""
+    atlas = LocalAtlas(atlas_dir["atlas"])
+    mesh = MeshHelper(atlas_dir, region)
 
-    WARNING: this function can only be used for Isocortex atlas for now.
-    """
-    atlas = LocalAtlas(atlas_dir)
-    brain_regions = atlas.load_data("brain_regions", memcache=True)
-    tolerance = 2.0 * brain_regions.voxel_dimensions[0]
-
-    # Layer thicknesses from J. Defilipe 2017 (unpublished), see Section 5.1.1.4
-    # of the release report "Neocortex Tissue Reconstruction",
-    # https://github.com/BlueBrain/ncx_release_report.git
-    max_thicknesses = [210.639, 190.2134, 450.6398, 242.554, 670.2, 893.62]
-
-    isocortex_mask = atlas.get_region_mask("Isocortex", memcache=True).raw
-
-    too_thick = np.full(isocortex_mask.shape, False, dtype=np.bool)
-    for i, max_thickness in enumerate(max_thicknesses, 1):
-        ph = atlas.load_data(f"[PH]{i}", memcache=True)
-        with np.errstate(invalid="ignore"):
-            invalid_thickness = (ph.raw[..., 1] - ph.raw[..., 0]) > (max_thickness + tolerance)
-        too_thick = np.logical_or(too_thick, invalid_thickness)
-
-        L.info(
-            "Layer %s with %s percentage of bad voxels",
-            i,
-            np.round(100 * invalid_thickness[isocortex_mask].mean(), 3),
-        )
+    region_mask = atlas.get_region_mask(region, memcache=True)
 
-    L.info(
-        "%s percentage of bad voxels in total",
-        np.round(100 * too_thick[isocortex_mask].mean(), 3),
-    )
+    boundary_mesh = mesh.get_boundary_mesh()
+    data = np.zeros_like(region_mask.raw, dtype=int)
+    data_vg = VoxelGrid(data)
+    # pylint: disable=no-member
+    data[tuple(data_vg.points_to_indices(boundary_mesh.triangles_center).T)] = 1
+    for _ in range(boundary_thickness):
+        data[:, :, :-1] += data[:, :, 1:]
+        data[:, :, 1:] += data[:, :, :-1]
+        data[:, :-1] += data[:, 1:]
+        data[:, 1:] += data[:, :-1]
+        data[:-1] += data[1:]
+        data[1:] += data[:-1]
 
-    return brain_regions.with_data(np.logical_and(~too_thick, isocortex_mask).astype(np.uint8))
+    mask = 1 - np.clip(data, 0, 1)
+    region_mask = atlas.get_region_mask(region, memcache=True)
+    mask[~region_mask.raw] = 0
+    return region_mask.with_data(mask.astype(np.uint8))
 
 
 def get_regions_from_composition(cell_composition_path):
@@ -86,15 +78,14 @@ def build_circuit(
     """
     if seed is not None:
         np.random.seed(seed)
-
     return place(
         composition_path=cell_composition_path,
         mtype_taxonomy_path=mtype_taxonomy_path,
         atlas_url=atlas_path,
         mini_frequencies_path=None,
         atlas_cache=None,
         region=region or get_regions_from_composition(cell_composition_path),
-        mask_dset=mask,
+        mask_dset=Path(mask).resolve().with_suffix("") if mask else None,
         density_factor=density_factor,
         soma_placement="basic",
         atlas_properties=None,

src/synthesis_workflow/tasks/circuit.py

@@ -7,7 +7,6 @@
 import luigi
 import pandas as pd
 import yaml
-from luigi.parameter import OptionalPathParameter
 from luigi.parameter import PathParameter
 from luigi_tools.parameter import RatioParameter
 from luigi_tools.task import ParamRef
@@ -19,7 +18,7 @@
 
 from synthesis_workflow.circuit import build_circuit
 from synthesis_workflow.circuit import circuit_slicer
-from synthesis_workflow.circuit import create_atlas_thickness_mask
+from synthesis_workflow.circuit import create_boundary_mask
 from synthesis_workflow.circuit import slice_circuit
 from synthesis_workflow.tasks.config import AtlasLocalTarget
 from synthesis_workflow.tasks.config import CircuitConfig
@@ -152,6 +151,37 @@ def output(self):
         return AtlasLocalTarget(f"{self.atlas_planes_path}{suffix}.npz")
 
 
+class CreateBoundaryMask(WorkflowTask):
+    """Create mask at the boundary of the region to prevent placing somata too close to boundary."""
+
+    boundary_thickness = luigi.IntParameter(
+        default=0,
+        description=":int: Thickness to create a mask to prevent placing cells near boundary"
+        ", in units of voxel size.",
+    )
+    mask_path = luigi.Parameter(default="boundary_mask.nrrd")
+
+    def requires(self):
+        """Required input tasks."""
+        return GetSynthesisInputs()
+
+    def run(self):
+        """Actual process of the task."""
+        thickness_mask = create_boundary_mask(
+            {
+                "atlas": CircuitConfig().atlas_path,
+                "structure": self.input().pathlib_path / CircuitConfig().region_structure_path,
+            },
+            CircuitConfig().region,
+            self.boundary_thickness,
+        )
+        thickness_mask.save_nrrd(self.output().path)
+
+    def output(self):
+        """Outputs of the task."""
+        return CircuitLocalTarget(self.mask_path)
+
+
 class BuildCircuit(WorkflowTask):
     """Generate cell positions and me-types from atlas, compositions and taxonomy."""
 
@@ -160,9 +190,6 @@ class BuildCircuit(WorkflowTask):
         left_op=luigi.parameter.operator.lt,
         description=":float: The density of positions generated from the atlas.",
     )
-    mask_path = OptionalPathParameter(
-        default=None, description=":str: Path to save thickness mask (NCx only)."
-    )
     seed = luigi.OptionalIntParameter(default=42, description=":int: Pseudo-random generator seed.")
     mtype_taxonomy_file = luigi.Parameter(
         default="mtype_taxonomy.tsv",
@@ -171,19 +198,17 @@ class BuildCircuit(WorkflowTask):
 
     def requires(self):
         """Required input tasks."""
-        return {"synthesis": GetSynthesisInputs(), "composition": GetCellComposition()}
+        tasks = {"synthesis": GetSynthesisInputs(), "composition": GetCellComposition()}
+        if CreateBoundaryMask().boundary_thickness > 0:
+            tasks["boundary"] = CreateBoundaryMask()
+
+        return tasks
 
     def run(self):
         """Actual process of the task."""
         # pylint: disable=no-member
         mtype_taxonomy_path = self.input()["synthesis"].pathlib_path / self.mtype_taxonomy_file
 
-        thickness_mask_path = None
-        if self.mask_path is not None:
-            thickness_mask = create_atlas_thickness_mask(CircuitConfig().atlas_path)
-            thickness_mask.save_nrrd(self.mask_path)
-            thickness_mask_path = self.mask_path.stem
-
         # create uniform densities if they don't exist
         mtypes = SynthesisConfig().mtypes
         if not mtypes:
@@ -196,7 +221,6 @@ def run(self):
                 if data["traits"]["mtype"] == mtype:
                     name = "".join(list(data["density"])[1:-1])
                     nrrd_path = Path(CircuitConfig().atlas_path) / f"{name}.nrrd"
-
             if not nrrd_path.exists():
                 annotation = get_layer_annotation(
                     {
@@ -219,7 +243,9 @@ def run(self):
             mtype_taxonomy_path,
             CircuitConfig().atlas_path,
             self.density_factor,
-            mask=thickness_mask_path,
+            mask=self.input()["boundary"].path
+            if CreateBoundaryMask().boundary_thickness > 0
+            else None,
             seed=self.seed,
             region=CircuitConfig().region,
         )

tox.ini

@@ -82,7 +82,7 @@ commands =
     pylint -j {env:PYLINT_NPROCS:1} {[base]files}
 
 [testenv:format]
-basepython = python3.8
+basepython = python3
 skip_install = true
 deps =
     codespell