From 8780abcb253432a008eb55d598782c2ae7397f65 Mon Sep 17 00:00:00 2001 From: Austin Soplata Date: Tue, 21 May 2024 14:43:47 +0200 Subject: [PATCH] Update thalamus (#11) * Overhaul the thalamus into steps, and add docs Load thalamus meshes from CLI input, clean CLI * Improve Blender instructions * Change thalamus region list regex This updates the regular expression used for thalamus placement-hints to be in a different format that has been tested successfully, excludes habenular and peripeduncular subregions, and to be valid for the hierarchy/annotation used at its appropriate step in the Atlas pipeline. For information on which regions were chosen and this list was created, see the internal BBP Confluence page located at "Circuits > Mouse Thalamus > Atlas-based Whole-thalamus subregion selection". This regex has been built from the region list of the desired and present thalamus regions as of the "final" version of the hierarchy and annotation built by the Atlas pipeline, which is the output of the rule `split_barrel_ccfv3_l23split`. This change is meant to go in tandem with https://github.com/BlueBrain/atlas-direction-vectors/pull/27 . * Update layer names to be Atlas-Pipeline-compatible * Fix formatting errors * Attempt to update tests for new thal workflow * Replace part of test anno with region in metadata * fix test + format * format * Fix final linting issues This does a lot of small things for passing the linting. For mypy, I had to add additional ignores for the Trimesh returned types, since the ignore on the module as a whole wasn't preventing mypy from expecting `load_mesh` to return a Geometry object, which is a grandparent of Trimesh objects. I don't know if Trimesh changed their API, I couldn't figure it out from the docs, and I don't know why mypy was raising this now. In all the cases I could test or see, a proper "Trimesh" object was returned instead of the more generic Geometry. I don't think we need to worry about this. For the pylint disable W0613 (unused-argument), I needed some polymorphism for the thalamus case, but I wasn't sure how to handle that alongside the linters' type-checking. I think this is the simplest solution and is harmless. Everything else is minor. * Make Alexis changes to CLI * Apply MG code review changes --------- Co-authored-by: Austin E. Soplata Co-authored-by: arnaudon --- .../app/metadata/thalamus_metadata.json | 7 +- atlas_placement_hints/app/placement_hints.py | 87 ++++- .../compute_placement_hints.py | 5 + .../distances/create_watertight_mesh.py | 6 +- atlas_placement_hints/layered_atlas.py | 344 ++++++++++++++++-- tests/app/test_placement_hints.py | 41 ++- tests/placement_hints/mocking_tools.py | 12 +- tests/placement_hints/test_layered_atlas.py | 2 +- 8 files changed, 440 insertions(+), 64 deletions(-) diff --git a/atlas_placement_hints/app/metadata/thalamus_metadata.json b/atlas_placement_hints/app/metadata/thalamus_metadata.json index ed0522c..5633e36 100644 --- a/atlas_placement_hints/app/metadata/thalamus_metadata.json +++ b/atlas_placement_hints/app/metadata/thalamus_metadata.json @@ -7,9 +7,10 @@ }, "layers": { - "names": ["RT", "THnotRT"], - "queries": ["RT", "@^\\bAD|AM|AMd|AMv|ATN|AV|CL|CM|DORpm|DORsm|EPI|Eth|GENd|GENv|IAD|IAM|IGL|ILM|IMD|IntG|LAT|LD|LGd|LGd-co|LGd-ip|LGd-sh|LGv|LGvl|LGvm|LH|LP|MD|MDc|MDl|MDm|MED|MG|MGd|MGm|MGv|MH|MTN|PCN|PF|PIL|PIN|PO|POL|PP|PR|PT|PVT|PoT|RE|REth|RH|SGN|SMT|SPA|SPF|SPFm|SPFp|SubG|TH|VAL|VENT|VM|VP|VPL|VPLpc|VPM|VPMpc|Xi\\b$"], + "names": ["layer_1", "layer_2"], + "queries": ["RT", +"@^(?:AD|AMd|AMv|AV|CL|CM|Eth|IAD|IAM|IGL|IMD|IntG|LD|LGd-co|LGd-ip|LGd-sh|LGv_O|LP|MD_O|MGd|MGm|MGv|PCN|PF|PIL|PO|POL|PR|PT|PVT|PoT|RE|RH|SGN|SMT|SPA|SPFm|SPFp|SubG|TH_O|VAL|VM|VPL|VPLpc|VPM|VPMpc|Xi)$"], "attribute": "acronym", "with_descendants": false } -} \ No newline at end of file +} diff --git a/atlas_placement_hints/app/placement_hints.py b/atlas_placement_hints/app/placement_hints.py index 9a195bf..0d4af04 100644 --- a/atlas_placement_hints/app/placement_hints.py +++ b/atlas_placement_hints/app/placement_hints.py @@ -74,6 +74,7 @@ def _placement_hints( # pylint: disable=too-many-locals max_thicknesses: Optional[List[float]] = None, flip_direction_vectors: bool = False, has_hemispheres: bool = False, + thalamus_meshes_dir: str = "", ) -> None: """ Compute the placement hints for a laminar region of the mouse brain. @@ -93,6 +94,10 @@ def _placement_hints( # pylint: disable=too-many-locals has_hemispheres: (optional) If True, split the volume into halves along the z-axis and handle each of theses 'hemispheres' separately. Otherwise the whole volume is handled. Defaults to True. + thalamus_meshes_dir: (optional) Path of the directory to load thalamus meshes + from. Currently only used for thalamus. Required if you are producing thalamus + placement-hints. Defaults to None. + """ direction_vectors = voxcell.VoxelData.load_nrrd(direction_vectors_path) assert_meta_properties([direction_vectors, atlas.region]) @@ -106,6 +111,7 @@ def _placement_hints( # pylint: disable=too-many-locals max_thicknesses, flip_direction_vectors=flip_direction_vectors, has_hemispheres=has_hemispheres, + thalamus_meshes_dir=thalamus_meshes_dir, ) if not Path(output_dir).exists(): os.makedirs(output_dir) @@ -305,41 +311,84 @@ def isocortex( required=True, help="path of the directory to write. It will be created if it doesn't exist.", ) +@click.option( + "--thalamus-meshes-dir", + required=True, + help="""Path of the directory to use for either saving or loading thalamus + meshes. It will be created if it doesn't exist.""", +) +@click.option( + "--load-cut-thalamus-meshes", + required=False, + help="""(Optional) Flag to load your custom thalamus meshes, and then use + them to calculate placement-hints.""", + default=False, + is_flag=True, +) @log_args(L) +# pylint: disable=too-many-arguments def thalamus( - verbose, annotation_path, hierarchy_path, metadata_path, direction_vectors_path, output_dir + verbose, + annotation_path, + hierarchy_path, + metadata_path, + direction_vectors_path, + output_dir, + thalamus_meshes_dir, + load_cut_thalamus_meshes, ): """Generate and save the placement hints of the mouse thalamus. - Placement hints are saved under the names sepecified in `app/metadata/thalamus_metadata.json`. - Default to: + First, call this without passing '--load-cut-thalamus-meshes' to + create your region meshes, but not your placement-hints. Then, hand-cut + your meshes according to the documentation in + 'atlas_placement_hints/layered_atlas.py::ThalamusAtlas'. Finally, call + this again while passing '--load-cut-thalamus-meshes'. + + Placement hints are saved under the names specified in + `app/metadata/thalamus_metadata.json`. These default to: \b - `[PH]y.nrrd` - - `[PH]Rt.nrrd`, `[PH]VPL.nrrd` + - `[PH]layer_1.nrrd` (This is the "top-most" layer, equivalent to "RT". + This previously used the filename `[PH]RT.nrrd`) + - `[PH]layer_2.nrrd` (This is the "deepest" layer, equivalent to the + thalamus except the habenular, peripeduncular, and reticular regions. + This previously used the filename `[PH]THnotRT.nrrd`) - A report together with an nrrd volume on problematic distance computations are generated - in `output_dir` under the names: + A report together with an nrrd volume on problematic distance computations + are generated in `output_dir` under the names: \b - `distance_report.json` - - `_problematic_voxel_mask.nrrd` (mask of the voxels for which the computed - placement hints cannot be trusted). is the region name specified in - thalamus_metadata.json. Defaults to "Thalamus". - - The annotation file can contain the thalamus or a superset. - For the algorithm to work properly, some space should separate the boundary - of the thalamus from the boundary of its enclosing array. + - `_problematic_voxel_mask.nrrd` (mask of the voxels for which + the computed placement hints cannot be trusted). is the + region name specified in thalamus_metadata.json. Defaults to "Thalamus". + + The annotation file can contain the thalamus or a superset. For the + algorithm to work properly, some space should separate the boundary of the + thalamus from the boundary of its enclosing array. + + For instructions on all steps necessary to generate the thalamus' placement + hints, see + 'atlas-placement-hints/atlas_placement_hints/layered_atlas.py::ThalamusAtlas' + and its methods for details. """ set_verbose(L, verbose) atlas = _create_layered_atlas(annotation_path, hierarchy_path, metadata_path) - _placement_hints( - atlas, - direction_vectors_path, - output_dir, - has_hemispheres=True, - ) + + if load_cut_thalamus_meshes: + _placement_hints( + atlas, + direction_vectors_path, + output_dir, + thalamus_meshes_dir=thalamus_meshes_dir, + has_hemispheres=True, + ) + else: + Path(thalamus_meshes_dir).mkdir(parents=True, exist_ok=True) + atlas.create_uncut_thalamus_meshes(thalamus_meshes_dir) @app.command() diff --git a/atlas_placement_hints/compute_placement_hints.py b/atlas_placement_hints/compute_placement_hints.py index 2a2bed4..08c261a 100644 --- a/atlas_placement_hints/compute_placement_hints.py +++ b/atlas_placement_hints/compute_placement_hints.py @@ -27,6 +27,7 @@ def compute_placement_hints( max_thicknesses: Optional[List[float]] = None, flip_direction_vectors: bool = False, has_hemispheres: bool = True, + thalamus_meshes_dir: str = "", ) -> Tuple[DistanceInfo, Dict]: """ Compute the placement hints for a laminar region of the mouse brain. @@ -43,6 +44,9 @@ def compute_placement_hints( has_hemispheres: (optional) If True, split the volume into halves along the z-axis and handle each of theses 'hemispheres' separately. Otherwise the whole volume is handled. Defaults to True. + thalamus_meshes_dir: (optional) Path of the directory to load thalamus meshes + from. Currently only used for thalamus. Required if you are producing thalamus + placement-hints. Defaults to None. Returns: Tuple with the following items. @@ -70,6 +74,7 @@ def compute_placement_hints( direction_vectors, flip_direction_vectors=flip_direction_vectors, has_hemispheres=has_hemispheres, + thalamus_meshes_dir=thalamus_meshes_dir, ) distances_to_meshes = distances_info["distances_to_layer_boundaries"] diff --git a/atlas_placement_hints/distances/create_watertight_mesh.py b/atlas_placement_hints/distances/create_watertight_mesh.py index 24bce5c..5125ab6 100644 --- a/atlas_placement_hints/distances/create_watertight_mesh.py +++ b/atlas_placement_hints/distances/create_watertight_mesh.py @@ -210,7 +210,7 @@ def create_watertight_trimesh( optimized_mesh = trimesh.load_mesh(output_filepath_opt) if optimization_info: - log_mesh_optimization_info(optimized_mesh, unoptimized_mesh) + log_mesh_optimization_info(optimized_mesh, unoptimized_mesh) # type: ignore - optimized_mesh.fix_normals() - return optimized_mesh + optimized_mesh.fix_normals() # type: ignore + return optimized_mesh # type: ignore diff --git a/atlas_placement_hints/layered_atlas.py b/atlas_placement_hints/layered_atlas.py index f0a7137..c9c9826 100644 --- a/atlas_placement_hints/layered_atlas.py +++ b/atlas_placement_hints/layered_atlas.py @@ -4,20 +4,24 @@ This module is used for the computation of placement hints in the mouse isocortex and in the mouse Hippocampus CA1 region. """ + from __future__ import annotations import logging import os +import sys from abc import ABC, abstractmethod from enum import IntEnum -from typing import TYPE_CHECKING, Dict, List, Union +from typing import Dict, List, Union import numpy as np +import trimesh # type: ignore from atlas_commons.typing import BoolArray, FloatArray, NDArray from atlas_commons.utils import create_layered_volume, query_region_mask, split_into_halves from cached_property import cached_property # type: ignore from cgal_pybind import estimate_thicknesses from tqdm import tqdm # type: ignore +from voxcell import RegionMap, VoxelData # type: ignore from atlas_placement_hints.distances.create_watertight_mesh import create_watertight_trimesh from atlas_placement_hints.distances.distances_to_meshes import ( @@ -31,10 +35,6 @@ get_convex_hull_boundary, ) -if TYPE_CHECKING: # pragma: no cover - import trimesh # type: ignore - from voxcell import RegionMap, VoxelData # type: ignore - logging.basicConfig(level=logging.INFO) L = logging.getLogger(__name__) @@ -213,16 +213,23 @@ def create_layer_meshes(self, layered_volume: NDArray[np.integer]) -> List["trim return meshes - def _compute_dists_and_obtuse_angles(self, volume, direction_vectors): + def _compute_dists_and_obtuse_angles( + self, volume, direction_vectors, hemisphere=LEFT, thalamus_meshes_dir: str = "" + ): + # pylint: disable=unused-argument layer_meshes = self.create_layer_meshes(volume) # pylint: disable=fixme # TODO: compute max_smooth_error and use it as the value of rollback_distance # in the call of distances_from_voxels_to_meshes_wrt_dir() return distances_from_voxels_to_meshes_wrt_dir(volume, layer_meshes, direction_vectors) - def _dists_and_obtuse_angles(self, direction_vectors, has_hemispheres=False): + def _dists_and_obtuse_angles( + self, direction_vectors, has_hemispheres=False, thalamus_meshes_dir: str = "" + ): if not has_hemispheres: - return self._compute_dists_and_obtuse_angles(self.volume, direction_vectors) + return self._compute_dists_and_obtuse_angles( + self.volume, direction_vectors, thalamus_meshes_dir=thalamus_meshes_dir + ) # Processing each hemisphere individually hemisphere_distances = [] hemisphere_volumes = split_into_halves(self.volume) @@ -237,7 +244,10 @@ def _dists_and_obtuse_angles(self, direction_vectors, has_hemispheres=False): self.metadata["region"]["name"], ) dists_to_layer_meshes, obtuse = self._compute_dists_and_obtuse_angles( - hemisphere_volumes[hemisphere], direction_vectors + hemisphere_volumes[hemisphere], + direction_vectors, + thalamus_meshes_dir=thalamus_meshes_dir, + hemisphere=hemisphere, ) hemisphere_distances.append(dists_to_layer_meshes) hemisphere_obtuse_angles.append(obtuse) @@ -257,6 +267,7 @@ def compute_distances_to_layer_boundaries( direction_vectors: FloatArray, has_hemispheres: bool = True, flip_direction_vectors: bool = False, + thalamus_meshes_dir: str = "", ) -> Dict[str, Union[FloatArray, BoolArray]]: """ Compute distances from voxels to layers boundaries wrt to direction vectors. @@ -273,6 +284,9 @@ def compute_distances_to_layer_boundaries( flip_direction_vectors: True if the direction vectors should be reverted, False otherwise. This flag needs to be set to True depending on the algorithm used to generated orientation.nrrd. + thalamus_meshes_dir: (optional) Path of the directory to load thalamus meshes + from. Currently only used for thalamus. Required if you are producing thalamus + placement-hints. Defaults to None. Returns: distances_info: dict with the following entries. @@ -287,7 +301,7 @@ def compute_distances_to_layer_boundaries( direction_vectors = -direction_vectors distances_to_layer_meshes, obtuse_angles = self._dists_and_obtuse_angles( - direction_vectors, has_hemispheres + direction_vectors, has_hemispheres, thalamus_meshes_dir=thalamus_meshes_dir ) L.info("Fixing disordered distances ...") fix_disordered_distances(distances_to_layer_meshes) @@ -303,33 +317,303 @@ class ThalamusAtlas(MeshBasedLayeredAtlas): """ Class holding the data of a two-layer atlas for the mouse thalamus. - The second layer of the thalamus, that is, the complement of the reticular nucleus, - cannot be defined via a simple regular expression because the thalamus (id = 549, non-leaf) - has voxels with labels 549 in both AIBS CCFv2 and CCFv3 mouse brain models. + The second layer of the thalamus, that is, the complement of the reticular + nucleus, cannot be defined via a simple regular expression because the + thalamus (id = 549, non-leaf) has voxels with labels 549 in both AIBS CCFv2 + and CCFv3 mouse brain models. + + (AES, <2023-06-28 Wed>: I'm not sure the above + comment applies anymore due to our recent, in-progress creation of + leaf-only annotations, but nonetheless the thalamus is a special case for + making its placement-hints.) + + If you're reading this, you're probably making new thalamus meshes because + the annotation has changed or some other reason. Generating placement-hints + for the thalamus has been changed, and now requires a manual step. You + should do the following steps: + + 1. Pass the argument '--thalamus-meshes-dir /your/folder/here' to the + top-level CLI command 'atlas-placement-hints thalamus'. This will create + the meshes, but NOT the placement-hints, and the program will exit. + + 2. MANUALLY cut the reticular meshes into the 'top' and 'bottom' halves + (aka the inner and outer halves if looking outwards from the center of the + thalamus) using software like 'Blender'. Do this for each hemisphere. The + RT region 'top' and 'bottom' halves are curvy and complex enough that it + was decided that manual cutting was the most effective / efficient way to + get a good mesh layer for the thalamus placement-hints, since the previous + computational way included too many holes due to the curviness. To do the + cutting in Blender, you can follow these instructions: + + A. Import each reticular mesh ('File > Import > Stl (.stl)'). + + B. Click the dropdown menu in the upper left corner that says 'Object + Mode' and change it to 'Edit Mode'. + + C. In the central window where the meshes are displayed, (you may have + to zoom out) select either the 'inner' (bottom) or 'outer' (top) half + of that hemisphere's reticular mesh. + + D. Click 'Mesh > Separate > Selection' (or press hotkey P then click + Selection). + + E. In the 'Scene Collection' window to the top right, which lists all + the meshes, you should now see a new entry that is named similarly to + 'reticular_nucleus_mesh_right_hemisphere_original.001', depending on + your input file. Click on its entry in the 'Scene Collection' window to + select it. + + F. Select 'File > Export > Stl (.stl)' to export this new mesh into its + own file, making sure to click the box that says 'Selection Only' in + the export prompt. Name the file appropriate to the hemisphere and + top/bottom selection that you've just done. If you're unsure which mesh + you have selected, you can click the 'Eye' symbols in the 'Scene + Collection' window to toggle which meshes are shown in the main view. + + G. Since you have 'separated' your mesh into two, the original mesh + object now consists of the remainder of the mesh which you did not + select in the previous steps. In other words, if you previously + selected, 'separated', and exported the 'bottom' of + 'reticular_nucleus_mesh_right_hemisphere_original', the entry in 'Scene + Collection' for the mesh + 'reticular_nucleus_mesh_right_hemisphere_original' now only consists of + the 'top' part of the mesh. Click the entry for it in 'Scene + Collection', then 'File > Export' like above, making sure to select + 'Selection Only' and name it appropriately. + + H. Repeat the process for the other hemisphere. + + Note that this 'separation' is a different operation than 'splitting'. + There are many tutorials on Youtube for how to do this. In the outgoing + filename, change 'original' to 'handcut'. In the end, you should end up + with four new files: + 'reticular_nucleus_mesh_left_hemisphere_bottom_handcut.stl', + 'reticular_nucleus_mesh_left_hemisphere_top_handcut.stl', + 'reticular_nucleus_mesh_right_hemisphere_bottom_handcut.stl', and + 'reticular_nucleus_mesh_right_hemisphere_top_handcut.stl'. After you have + made the 4 new files, I recommend you view them individually (using + software like Paraview) to double check that you named them correctly, etc. + + 3. Re-run the top-level command 'atlas-placement-hints thalamus' but this + time with both the argument '--thalamus-meshes-dir /your/folder/here' and + the flag '--load-cut-thalamus-meshes'. This will NOT create the meshes, but + WILL create the placement-hints using your newly handcut meshes! Note that + this must be done using a compute node with at least as much RAM as that of + 'memory_bound' in + 'atlas-placement-hints/atlas_placement_hints/distances/ + utils.py:memory_efficient_intersection()', + otherwise this WILL fail silently! With 'memory_bound' set to 300, if your + meshes are not very optimized (e.g. using ultraliser's default 1 + '--optimization-iterations'), and have a similar number of faces to RT mask + surface voxels, this should take approximately 1.5 hours for the whole + thing. If your meshes are more optimized, this can significantly speed up + the process. + """ - def create_layer_meshes(self, layered_volume: NDArray[np.integer]) -> List["trimesh.Trimesh"]: + def create_uncut_thalamus_meshes(self, thalamus_meshes_dir: str): """ Create meshes representing the upper boundary of each layer of the thalamus atlas. + + Because the lower boundary of the thalamus is too irregular to obtain + meaningful ray-mesh interesections, we consider instead its convex + hull, which provides us with a smooth approximation. See pictures and + discussion of https://bbpteam.epfl.ch/project/issues/browse/NSETM-1433 + + AES <2023-06-28 Wed>: This now creates + meshes that are expected to be manually cut using Blender/etc., as + described in the 'ThalamusAtlas' class docstring. Also, note that this + creates the 'upper boundary of each layer', but ALSO creates the + bottom-most layer as well (see the docstring for + 'ThalamusAtlas.load_layer_meshes()' for details). """ - # Because the lower boundary of the thalamus is too irregular to obtain meaningful - # ray-mesh interesections, we consider instead its convex hull, which provides us with a - # smooth approximation. See pictures and discussion of - # https://bbpteam.epfl.ch/project/issues/browse/NSETM-1433 - thalamus_convex_hull_boundary = get_convex_hull_boundary(layered_volume) - hull_mask = detailed_mesh_mask(thalamus_convex_hull_boundary, layered_volume.shape) - reticular_nucleus_mesh = create_watertight_trimesh(layered_volume == 1) - reticular_nucleus_mesh_top = clip_mesh(reticular_nucleus_mesh, hull_mask) - reticular_nucleus_mesh_bottom = clip_mesh(reticular_nucleus_mesh, hull_mask, remainder=True) - reticular_nucleus_mesh_bottom.invert() - overall_bottom = clip_mesh( - thalamus_convex_hull_boundary, - ~detailed_mesh_mask(reticular_nucleus_mesh, layered_volume.shape), + + L.info( + """Currently set to CREATE meshes, but NOT calculate thalamus + placement-hints. See + 'atlas-placement-hints/atlas_placement_hints/layered_atlas.py:ThalamusAtlas' + for details.""" ) - overall_bottom.fix_normals() - overall_bottom.invert() + hemisphere_volumes = split_into_halves(self.volume) + for hemisphere in [LEFT, RIGHT]: + if hemisphere == LEFT: + hemisphere_string = "left" + elif hemisphere == RIGHT: + hemisphere_string = "right" + L.info( + "Creating uncut meshes for the %s hemisphere (hemisphere %d) of the %s region ...", + hemisphere_string, + hemisphere, + self.metadata["region"]["name"], + ) + + # Because the lower boundary of the thalamus is too irregular to obtain meaningful + # ray-mesh interesections, we consider instead its convex hull, which provides us with a + # smooth approximation. See pictures and discussion of + # https://bbpteam.epfl.ch/project/issues/browse/NSETM-1433 + thalamus_convex_hull_boundary = get_convex_hull_boundary(hemisphere_volumes[hemisphere]) + hull_mask = detailed_mesh_mask( + thalamus_convex_hull_boundary, hemisphere_volumes[hemisphere].shape + ) + reticular_nucleus_mesh = create_watertight_trimesh(hemisphere_volumes[hemisphere] == 1) + reticular_nucleus_mesh_top = clip_mesh(reticular_nucleus_mesh, hull_mask) + reticular_nucleus_mesh_bottom = clip_mesh( + reticular_nucleus_mesh, hull_mask, remainder=True + ) + reticular_nucleus_mesh_bottom.invert() + overall_bottom = clip_mesh( + thalamus_convex_hull_boundary, + ~detailed_mesh_mask(reticular_nucleus_mesh, hemisphere_volumes[hemisphere].shape), + ) + overall_bottom.fix_normals() + overall_bottom.invert() + + thalamus_convex_hull_boundary.export( + os.path.join( + thalamus_meshes_dir, + f"thalamus_convex_hull_boundary_{hemisphere_string}_hemisphere_original.stl", + ) + ) + reticular_nucleus_mesh.export( + os.path.join( + thalamus_meshes_dir, + f"reticular_nucleus_mesh_{hemisphere_string}_hemisphere_original.stl", + ) + ) + reticular_nucleus_mesh_bottom.export( + os.path.join( + thalamus_meshes_dir, + f"reticular_nucleus_mesh_{hemisphere_string}_hemisphere_bottom_original.stl", + ) + ) + reticular_nucleus_mesh_top.export( + os.path.join( + thalamus_meshes_dir, + f"reticular_nucleus_mesh_{hemisphere_string}_hemisphere_top_original.stl", + ) + ) + overall_bottom.export( + os.path.join( + thalamus_meshes_dir, + f"overall_bottom_{hemisphere_string}_hemisphere_original.stl", + ) + ) - return [reticular_nucleus_mesh_top, reticular_nucleus_mesh_bottom, overall_bottom] + L.info("Finished creating and saving meshes for both hemispheres. Exiting.") + sys.exit() + + # pylint: disable=W0613 + def load_layer_meshes( + self, layered_volume: NDArray[np.integer], thalamus_meshes_dir: str, hemisphere=LEFT + ) -> List["trimesh.Trimesh"]: + """ + Load thalamus meshes meshes for bottom-most layer and upper boundary of other layers. + + Long-winded explanation by AES: + + Inside 'layered_volume', voxels marked 0 belong to 'outside' the + relevant volume, and can be thought of as an unused 'Layer 0' that + corresponds to the 'uppermost' space that is not considered a 'real' + layer (like the 'pia' in cortex). + + Inside 'layered_volume', voxels marked 1 belong to the reticular + nucleus (RT), and can be thought of as 'Layer 1' or the 'uppermost' + layer we actually care about, similar to cortex L1. In this function, + we seek to load a hand-cut 'upper boundary mesh' for this layer, + corresponding to the outermost half (loosely speaking) of a mesh of the + RT. This upper boundary mesh for Layer 1 (RT) is the first object in + the list that is returned by the function. + + Inside 'layered_volume', voxels marked 2 belong to non-RT thalamus, and + can be thought of as 'Layer 2' or the next-deeper layer, similar to + cortex L2. In this function, we seek to load a hand-cut 'upper boundary + mesh' for this layer, corresponding to the innermost half (loosely + speaking) of a mesh of the RT (since we're assuming there's no gap + between RT and nonRT, in general). This upper boundary mesh for Layer 2 + (nonRT) is the second object in the list that is returned by this + function. + + The final object returned by this function is the 'lower boundary mesh' + for the entire volume that we care about, aka the 'lowest boundary + mesh'. (Using the previous upper boundaries and this single lower + boundary is enough information to compute the distance from anywhere to + the nearest upper and lower boundaries of each layer.) Confusingly, and + unlike cortex, our lowest boundary is essentially the same as the + boundary of Layer 0: the convex hull mesh of the entire thalamus. This + is because all the direction vectors in nonRT (our Layer 2) point to RT + (our Layer 1), and the origins of those direction vectors are + distributed along almost a full half-sphere if looking from RT, at much + wider angles than a more laminar structure like the cortex. + + Currently, since the lowest boundary mesh is actually the entire convex + hull, that means there are some portions of the mesh where the 'bottom' + layer boundary actually sits 'above' the uppermost boundary of the + topmost layer! However, this shouldn't cause any problems. + """ + if thalamus_meshes_dir == "": + L.info( + """\n --> ERROR: You must pass a directory containing the + appropriate meshes to '--thalamus-meshes-dir'. See + 'atlas-placement-hints/atlas_placement_hints/layered_atlas.py:ThalamusAtlas' + for details. Exiting.\n""" + ) + sys.exit() + + if hemisphere == LEFT: + hemisphere_string = "left" + elif hemisphere == RIGHT: + hemisphere_string = "right" + + L.info( + """Loading hand-sliced thalamic reticular nucleus top and bottom + meshes of %s hemisphere from '%s'""", + hemisphere_string, + thalamus_meshes_dir, + ) + + reticular_nucleus_mesh_top = trimesh.load_mesh( + os.path.join( + thalamus_meshes_dir, + f"reticular_nucleus_mesh_{hemisphere_string}_hemisphere_top_handcut.stl", + ) + ) + + reticular_nucleus_mesh_bottom = trimesh.load_mesh( + os.path.join( + thalamus_meshes_dir, + f"reticular_nucleus_mesh_{hemisphere_string}_hemisphere_bottom_handcut.stl", + ) + ) + # We have to invert the normals of the bottom mesh faces after hand-cutting, since by + # default they face "outward", not "inward" like we want (so they can align with the + # direction vectors) + reticular_nucleus_mesh_bottom.invert() # type: ignore # type: ignore + + overall_bottom = trimesh.load_mesh( + os.path.join( + thalamus_meshes_dir, + f"overall_bottom_{hemisphere_string}_hemisphere_original.stl", + ) + ) + + return [ + reticular_nucleus_mesh_top, # type: ignore # type: ignore + reticular_nucleus_mesh_bottom, # type: ignore # type: ignore + overall_bottom, # type: ignore # type: ignore + ] + + def _compute_dists_and_obtuse_angles( + self, volume, direction_vectors, hemisphere=LEFT, thalamus_meshes_dir: str = "" + ): + # Note that this is LOADING meshes, not creating them! + L.info( + """Currently set to LOAD meshes and calculate thalamus + placement-hints. See + 'atlas-placement-hints/atlas_placement_hints/layered_atlas.py:ThalamusAtlas' + for details.""" + ) + layer_meshes = self.load_layer_meshes(volume, thalamus_meshes_dir, hemisphere) + return distances_from_voxels_to_meshes_wrt_dir(volume, layer_meshes, direction_vectors) class VoxelBasedLayeredAtlas(AbstractLayeredAtlas): diff --git a/tests/app/test_placement_hints.py b/tests/app/test_placement_hints.py index 6561dd6..be85213 100644 --- a/tests/app/test_placement_hints.py +++ b/tests/app/test_placement_hints.py @@ -1,5 +1,6 @@ """test app/placement_hints""" import json +from shutil import copyfile import numpy as np import numpy.testing as npt @@ -36,10 +37,48 @@ def test_thalamus(): "direction_vectors.nrrd", "--output-dir", "placement_hints", + "--thalamus-meshes-dir", + "thalamus_meshes", + "--create-uncut-thalamus-meshes-flag", ], ) assert result.exit_code == 0, str(result.output) + copyfile( + "thalamus_meshes/reticular_nucleus_mesh_left_hemisphere_bottom_original.stl", + "thalamus_meshes/reticular_nucleus_mesh_left_hemisphere_bottom_handcut.stl", + ) + copyfile( + "thalamus_meshes/reticular_nucleus_mesh_left_hemisphere_top_original.stl", + "thalamus_meshes/reticular_nucleus_mesh_left_hemisphere_top_handcut.stl", + ) + copyfile( + "thalamus_meshes/reticular_nucleus_mesh_right_hemisphere_bottom_original.stl", + "thalamus_meshes/reticular_nucleus_mesh_right_hemisphere_bottom_handcut.stl", + ) + copyfile( + "thalamus_meshes/reticular_nucleus_mesh_right_hemisphere_top_original.stl", + "thalamus_meshes/reticular_nucleus_mesh_right_hemisphere_top_handcut.stl", + ) + + result2 = runner.invoke( + tested.thalamus, + [ + "--annotation-path", + "annotation.nrrd", + "--hierarchy-path", + "hierarchy.json", + "--direction-vectors-path", + "direction_vectors.nrrd", + "--output-dir", + "placement_hints", + "--thalamus-meshes-dir", + "thalamus_meshes", + "--load-cut-thalamus-meshes-flag", + ], + ) + assert result2.exit_code == 0, str(result2.output) + # The values selected below as upper bounds are surprisingly large, which can be explained # as follows. Due to the shape and the size of the simplified brain region under test, # voxels close to the boundary of the volume are problematic (rays issued from them @@ -107,7 +146,7 @@ def test_thalamus(): ph_y = VoxelData.load_nrrd("placement_hints/[PH]y.nrrd") npt.assert_array_equal(ph_y.raw.shape, thalamus_mock.annotation.raw.shape) - ph_th_no_rt = VoxelData.load_nrrd("placement_hints/[PH]THnotRT.nrrd") + ph_th_no_rt = VoxelData.load_nrrd("placement_hints/[PH]layer_2.nrrd") npt.assert_array_equal(ph_th_no_rt.raw.shape, thalamus_mock.annotation.raw.shape + (2,)) diff --git a/tests/placement_hints/mocking_tools.py b/tests/placement_hints/mocking_tools.py index d5647f5..dd8e355 100644 --- a/tests/placement_hints/mocking_tools.py +++ b/tests/placement_hints/mocking_tools.py @@ -1,6 +1,7 @@ """ Mocking tools for the unit tests related to placement hints. """ + from typing import Tuple import numpy as np @@ -209,7 +210,8 @@ def __init__(self, padding: int, shape: Tuple[int, int, int], layer_thickness_ra raw = np.zeros(shape, dtype=int) reticular_nucleus_thickness = int(layer_thickness_ratio * shape[0]) raw[:reticular_nucleus_thickness, ...] = 262 # Region id of the reticular nucleus (RT) - raw[reticular_nucleus_thickness:, ...] = 549 # Region id of the thalamus (TH) + # Region id of the Anterodorsal nucleus (AD) of the thalamus + raw[reticular_nucleus_thickness:, ...] = 64 self.volume = shape[0] * shape[1] * shape[2] # Number of voxels with positive labels @@ -224,14 +226,10 @@ def __init__(self, padding: int, shape: Tuple[int, int, int], layer_thickness_ra "acronym": "TH", "children": [ { - "id": 856, - "acronym": "DORpm", + "id": 64, + "acronym": "AD", "children": [{"acronym": "RT", "id": 262}], }, - { - "id": 864, - "acronym": "DORsm", - }, ], }, ], diff --git a/tests/placement_hints/test_layered_atlas.py b/tests/placement_hints/test_layered_atlas.py index 5ab19c6..c4e7c7e 100644 --- a/tests/placement_hints/test_layered_atlas.py +++ b/tests/placement_hints/test_layered_atlas.py @@ -165,7 +165,7 @@ def setUpClass(cls): def test_create_layered_volume(self): expected = self.thalamus_mock.annotation.raw.copy() - expected[expected == 549] = 2 + expected[expected == 64] = 2 expected[expected == 262] = 1 npt.assert_array_equal(self.thalamus_atlas.volume, expected)