diff --git a/src/synthesis_workflow/circuit.py b/src/synthesis_workflow/circuit.py index fd15652..3ce65bd 100644 --- a/src/synthesis_workflow/circuit.py +++ b/src/synthesis_workflow/circuit.py @@ -1,7 +1,9 @@ """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 @@ -9,6 +11,7 @@ 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], diff --git a/src/synthesis_workflow/tasks/circuit.py b/src/synthesis_workflow/tasks/circuit.py index 100d5b2..335f56c 100644 --- a/src/synthesis_workflow/tasks/circuit.py +++ b/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) diff --git a/src/synthesis_workflow/tasks/synthesis.py b/src/synthesis_workflow/tasks/synthesis.py index f32c3d7..9a1d331 100644 --- a/src/synthesis_workflow/tasks/synthesis.py +++ b/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 diff --git a/src/version.py b/src/version.py index 01978b3..189f240 100644 --- a/src/version.py +++ b/src/version.py @@ -1,2 +1,2 @@ """Package version.""" -VERSION = "0.0.6" +VERSION = "0.0.7.dev0"