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)