added creation of thickness mask for Isocortex

Change-Id: I866e1b26b3ca5946ecbeb240381cf9db2078a3f4

src/synthesis_workflow/circuit.py

@@ -1,14 +1,17 @@
 """Functions for slicing mvd3 circuit files to place specific cells only."""
+import logging
 import pandas as pd
 from tqdm import tqdm
 from voxcell import CellCollection
+from voxcell.nexus.voxelbrain import LocalAtlas
 import numpy as np
 
 from atlas_analysis.planes.planes import create_planes as _create_planes
 from atlas_analysis.planes.planes import create_centerline as _create_centerline
 from atlas_analysis.planes.planes import _smoothing
 from brainbuilder.app.cells import _place as place
 
+L = logging.getLogger(__name__)
 LEFT = 0
 RIGHT = 1
 
@@ -49,11 +52,53 @@ def halve_atlas(annotated_volume, axis=0, side=LEFT):
     return annotated_volume
 
 
+def create_atlas_thickness_mask(atlas_dir):
+    """Create a mask on an atlas to to select voxels with small enough thicknesses values.
+
+    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),
+        )
+
+    L.info(
+        "%s percentage of bad voxels in total",
+        np.round(100 * too_thick[isocortex_mask].mean(), 3),
+    )
+
+    return brain_regions.with_data(
+        np.logical_and(~too_thick, isocortex_mask).astype(np.uint8)
+    )
+
+
 def build_circuit(
     cell_composition_path,
     mtype_taxonomy_path,
     atlas_path,
     density_factor=0.01,
+    mask=None,
     seed=None,
 ):
     """Builds a new circuit by calling ``brainbuilder.app.cells._place``.
@@ -72,7 +117,7 @@ def build_circuit(
         mini_frequencies_path=None,
         atlas_cache=None,
         region=None,
-        mask_dset=None,
+        mask_dset=mask,
         density_factor=density_factor,
         soma_placement="basic",
         atlas_properties=None,
@@ -185,8 +230,8 @@ def create_planes(
 
     elif plane_type == "aligned":
         _n_points = 10
-        bbox = layer_annotation.bbox
         centerline = np.zeros([_n_points, 3])
+        bbox = layer_annotation.bbox
         centerline[:, centerline_axis] = np.linspace(
             bbox[0, centerline_axis],
             bbox[1, centerline_axis],

src/synthesis_workflow/tasks/circuit.py

@@ -14,6 +14,7 @@
 from synthesis_workflow.circuit import create_planes
 from synthesis_workflow.circuit import halve_atlas
 from synthesis_workflow.circuit import slice_circuit
+from synthesis_workflow.circuit import create_atlas_thickness_mask
 from synthesis_workflow.tasks.config import AtlasLocalTarget
 from synthesis_workflow.tasks.config import CircuitConfig
 from synthesis_workflow.tasks.config import CircuitLocalTarget
@@ -158,6 +159,9 @@ class BuildCircuit(WorkflowTask):
         left_op=luigi.parameter.operator.lt,
         description="The density of positions generated from the atlas",
     )
+    mask_path = luigi.Parameter(
+        default=None, description="path to save thickness mask (NCx only)"
+    )
     seed = luigi.IntParameter(default=None, description="pseudo-random generator seed")
 
     def requires(self):
@@ -168,12 +172,20 @@ def run(self):
         """"""
         cell_composition_path = self.input().ppath / self.cell_composition_path
         mtype_taxonomy_path = self.input().ppath / self.mtype_taxonomy_path
+
+        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 = Path(self.mask_path).stem
+
         cells = build_circuit(
             cell_composition_path,
             mtype_taxonomy_path,
             CircuitConfig().atlas_path,
             self.density_factor,
-            self.seed,
+            mask=thickness_mask_path,
+            seed=self.seed,
         )
         ensure_dir(self.output().path)
         cells.save(self.output().path)

src/synthesis_workflow/tasks/synthesis.py

@@ -284,13 +284,22 @@ class BuildAxonMorphologies(WorkflowTask):
     placement_seed = luigi.IntParameter(default=0)
     nb_jobs = luigi.IntParameter(default=20)
 
+    def get_neuron_db_path(self, db_name):
+        """Helper function to fix neuronDB vs neurondb in file names."""
+        return (
+            Path(self.axon_cells_path) / Path(self.neurondb_basename).parent / db_name
+        ).with_suffix(".xml")
+
     def requires(self):
         """"""
         tasks = {"circuit": SliceCircuit()}
+
+        neurondb_path = self.get_neuron_db_path("neurondb")
+        if not neurondb_path.exists():
+            neurondb_path = self.get_neuron_db_path("neuronDB")
+
         tasks["axon_cells"] = BuildAxonMorphsDF(
-            neurondb_path=(
-                Path(self.axon_cells_path) / self.neurondb_basename
-            ).with_suffix(".xml"),
+            neurondb_path=neurondb_path,
             morphology_dirs={"clone_path": self.axon_cells_path},
         )
         return tasks

src/version.py

@@ -1,2 +1,2 @@
 """Package version."""
-VERSION = "0.0.6"
+VERSION = "0.0.7.dev0"