diff --git a/MANIFEST.in b/MANIFEST.in index 5ba9d7b1..d7ae6fbc 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -11,3 +11,4 @@ graft brainglobe_atlasapi *.py prune tests prune tutorials prune brainglobe_atlasapi/atlas_generation/atlas_scripts/* +prune brainglobe_atlasapi/atlas_generation/hashes diff --git a/brainglobe_atlasapi/atlas_generation/atlas_scripts/kim_devccf_mouse.py b/brainglobe_atlasapi/atlas_generation/atlas_scripts/kim_devccf_mouse.py new file mode 100644 index 00000000..eae75af6 --- /dev/null +++ b/brainglobe_atlasapi/atlas_generation/atlas_scripts/kim_devccf_mouse.py @@ -0,0 +1,355 @@ +import argparse +import time +from pathlib import Path + +import numpy as np +import pandas as pd +import pooch +from brainglobe_utils.IO.image import load_nii +from rich.progress import track + +from brainglobe_atlasapi import utils +from brainglobe_atlasapi.atlas_generation.mesh_utils import ( + Region, + create_region_mesh, +) +from brainglobe_atlasapi.atlas_generation.wrapup import wrapup_atlas_from_data +from brainglobe_atlasapi.config import DEFAULT_WORKDIR +from brainglobe_atlasapi.structure_tree_util import get_structures_tree + +ATLAS_NAME = "kim_dev_mouse" +SPECIES = "Mus musculus" +ATLAS_LINK = "https://kimlab.io/brain-map/DevCCF/" +CITATION = "Kronman, F.N., Liwang, J.K., Betty, R. et al. 2024, https://doi.org/10.1038/s41467-024-53254-w" +ORIENTATION = ["left", "superior", "posterior"] +ROOT_ID = 15564 +VERSION = 3 +PACKAGER = "Carlo Castoldi " +ATLAS_FILE_URL = "https://doi.org/10.6084/m9.figshare.26377171.v1" + +TIMEPOINTS = ("E11.5", "E13.5", "E15.5", "E18.5", "P04", "P14", "P56") +MODALITIES = ( + "LSFM", # Light Sheet Fluorescence Microscopy + "MRI-adc", # MRI Apparent Diffusion Coefficient + "MRI-dwi", # MRI Difusion Weighted Imaging + "MRI-fa", # MRI Fractional Anisotropy + "MRI-MTR", # MRI Magnetization Transfer Ratio + "MRI-T2", # MRI T2-weighted +) + + +def pooch_init(download_dir_path: Path, timepoints: list[str]) -> pooch.Pooch: + utils.check_internet_connection() + + empty_registry = { + a + ".zip": None + for a in [*timepoints, "DevCCFv1_OntologyStructure.xlsx"] + } + p = pooch.create( + path=download_dir_path, + base_url="doi:10.6084/m9.figshare.26377171.v1/", + registry=empty_registry, + ) + p.load_registry( + Path(__file__).parent.parent / "hashes" / (ATLAS_NAME + ".txt") + ) + return p + + +def fetch_animal(pooch_: pooch.Pooch, age: str, modality: str): + assert age in TIMEPOINTS, f"Unknown age timepoint: '{age}'" + archive = age + ".zip" + if modality == "LSFM": + resolution_um = 20 + elif modality in MODALITIES: + match age: + case "E11.5": + resolution_um = 31.5 + case "E13.5": + resolution_um = 34 + case "E15.5": + resolution_um = 37.5 + case "E18.5": + resolution_um = 40 + case _: + resolution_um = 50 + else: + raise RuntimeError(f"Unknown reference image modality: {modality}") + members = [ + f"{age}/{age.replace('.','-')}_DevCCF_Annotations_{resolution_um}um.nii.gz", + f"{age}/{age.replace('.','-')}_{modality}_{resolution_um}um.nii.gz", + ] + fetched_paths = pooch_.fetch( + archive, + progressbar=True, + processor=pooch.Unzip(extract_dir=".", members=members), + ) + # the file paths returned by pooch.ExtractorProcessor (superclass of Unzip) + # may not respect the order of the given members. + # see: https://github.com/fatiando/pooch/issues/457 + annotations_path = next(p for p in fetched_paths if p.endswith(members[0])) + reference_path = next(p for p in fetched_paths if p.endswith(members[1])) + annotations = load_nii(annotations_path, as_array=True) + reference = load_nii(reference_path, as_array=True) + return annotations, reference, resolution_um + + +def fetch_ontology(pooch_: pooch.Pooch): + devccfv1_path = pooch_.fetch( + "DevCCFv1_OntologyStructure.xlsx", progressbar=True + ) + xl = pd.ExcelFile(devccfv1_path) + # xl.sheet_names # it has two excel sheets + # 'DevCCFv1_Ontology', 'README' + df = xl.parse("DevCCFv1_Ontology", header=1) + df = df[["Acronym", "ID16", "Name", "Structure ID Path16", "R", "G", "B"]] + df.rename( + columns={ + "Acronym": "acronym", + "ID16": "id", + "Name": "name", + "Structure ID Path16": "structure_id_path", + "R": "r", + "G": "g", + "B": "b", + }, + inplace=True, + ) + structures = list(df.to_dict(orient="index").values()) + for structure in structures: + if structure["acronym"] == "mouse": + structure["acronym"] = "root" + structure_path = structure["structure_id_path"] + structure["structure_id_path"] = [ + int(id) for id in structure_path.strip("/").split("/") + ] + structure["rgb_triplet"] = [ + structure["r"], + structure["g"], + structure["b"], + ] + del structure["r"] + del structure["g"] + del structure["b"] + return structures + + +def create_meshes( + output_path: str | Path, + structures, + annotation_volume, + root_id, + decimate_fraction, +): + if not isinstance(output_path, Path): + output_path = Path(output_path) + output_path.mkdir(exist_ok=True) + + tree = get_structures_tree(structures) + + labels = np.unique(annotation_volume).astype(np.uint16) + + for key, node in tree.nodes.items(): + if key in labels: + is_label = True + else: + is_label = False + node.data = Region(is_label) + + # Mesh creation + closing_n_iters = 2 + smooth = True # smooth meshes after creation + start = time.time() + for node in track( + tree.nodes.values(), + total=tree.size(), + # list(tree.nodes.values())[:5], + # total=5, + description="Creating meshes", + ): + output_file = output_path / f"{node.identifier}.obj" + if output_file.exists(): + continue + create_region_mesh( + ( + output_path, + node, + tree, + labels, + annotation_volume, + root_id, + closing_n_iters, + decimate_fraction, + smooth, + ) + ) + + print( + "Finished mesh extraction in: ", + round((time.time() - start) / 60, 2), + " minutes", + ) + return output_path + + +def create_mesh_dict(structures, meshes_dir_path): + meshes_dict = dict() + structures_with_mesh = [] + for s in structures: + # Check if a mesh was created + mesh_path = meshes_dir_path / f'{s["id"]}.obj' + if not mesh_path.exists(): + print(f"No mesh file exists for: {s}, ignoring it") + continue + else: + # Check that the mesh actually exists (i.e. not empty) + if mesh_path.stat().st_size < 512: + print(f"obj file for {s} is too small, ignoring it.") + continue + structures_with_mesh.append(s) + meshes_dict[s["id"]] = mesh_path + print( + f"In the end, {len(structures_with_mesh)} " + "structures with mesh are kept" + ) + return meshes_dict, structures_with_mesh + + +class HideChoicesRawTextHelpFormatter(argparse.RawTextHelpFormatter): + def _get_help_string(self, action): + help_text = action.help + if action.choices: + help_text = help_text.split("(choices")[0].strip() + return help_text + + def _format_action_invocation(self, action): + if action.option_strings: + return ", ".join(action.option_strings) + if action.metavar: + return action.metavar + return super()._format_action_invocation(action) + + +arg_parser = argparse.ArgumentParser( + formatter_class=HideChoicesRawTextHelpFormatter +) +timepoints_help = """the age timepoint at which the atlas will be generated. + Options are: + - E11.5 ⅂ + - E13.5 |- Embryonic day + - E15.5 | + - E18.5 ⅃ + - P04 ⅂ + - P14 |- Postnatal day + - P56 ⅃ + +By default, it generates an atlas for each timepoint. +""" +modalities_help = """the reference image acquisition modality. + Options are: + - LSFM, Light Sheet Fluorescence Microscopy + - MRI-adc MRI Apparent Diffusion Coefficient + - MRI-dwi MRI Difusion Weighted Imaging + - MRI-fa MRI Fractional Anisotropy + - MRI-MTR MRI Magnetization Transfer Ratio + - MRI-T2 MRI T2-weighted + +By default, LSFM +""" +decimate_fraction_help = ( + "fraction of original meshes' vertices to be kept." + " Must be a number > 0 and <= 1.\nBy default, 0.2" +) + + +def decimate_fraction_type(arg): + try: + f = float(arg) + if 0 < f <= 1: + return f + except ValueError: + pass + raise argparse.ArgumentTypeError( + f"invalid value: '{arg}' (must be a number > 0 and <= 1)" + ) + + +arg_parser.add_argument( + "-t", + "--timepoints", + default=TIMEPOINTS, + type=str, + nargs="+", + choices=TIMEPOINTS, + help=timepoints_help, +) +arg_parser.add_argument( + "-m", + "--modality", + default="LSFM", + type=str, + choices=MODALITIES, + help=modalities_help, +) +arg_parser.add_argument( + "-d", + "--decimate-fraction", + default=0.2, + type=decimate_fraction_type, + help=decimate_fraction_help, +) + +if __name__ == "__main__": + params = vars(arg_parser.parse_args()) + timepoints = params["timepoints"] + modality = params["modality"] + decimate_fraction = params["decimate_fraction"] + bg_root_dir = DEFAULT_WORKDIR / ATLAS_NAME + download_dir_path = bg_root_dir / "downloads" + download_dir_path.mkdir(exist_ok=True, parents=True) + pooch_ = pooch_init(download_dir_path, timepoints) + structures = fetch_ontology(pooch_) + + for age in timepoints: + atlas_name = f"{ATLAS_NAME}_{age.replace('.', '-')}_{modality}" + annotation_volume, reference_volume, resolution_um = fetch_animal( + pooch_, age, modality + ) + atlas_dir = bg_root_dir / atlas_name + atlas_dir.mkdir(exist_ok=True) + print(f"Saving atlas data at {atlas_dir}") + # Create meshes: + meshes_dir_path = atlas_dir / "meshes" + create_meshes( + meshes_dir_path, + structures, + annotation_volume, + ROOT_ID, + decimate_fraction, + ) + meshes_dict, structures_with_mesh = create_mesh_dict( + structures, meshes_dir_path + ) + # Wrap up, compress, and remove file: + print("Finalising atlas") + output_filename = wrapup_atlas_from_data( + atlas_name=atlas_name, + atlas_minor_version=VERSION, + citation=CITATION, + atlas_link=ATLAS_LINK, + species=SPECIES, + resolution=(resolution_um,) * 3, + orientation=ORIENTATION, + root_id=ROOT_ID, + reference_stack=reference_volume, + annotation_stack=annotation_volume, + structures_list=structures_with_mesh, + meshes_dict=meshes_dict, + working_dir=atlas_dir, + atlas_packager=PACKAGER, + hemispheres_stack=None, # it is symmetric + cleanup_files=False, + compress=True, + scale_meshes=True, + ) + print("Done. Atlas generated at: ", output_filename) diff --git a/brainglobe_atlasapi/atlas_generation/hashes/kim_dev_mouse.txt b/brainglobe_atlasapi/atlas_generation/hashes/kim_dev_mouse.txt new file mode 100644 index 00000000..d5c89d9d --- /dev/null +++ b/brainglobe_atlasapi/atlas_generation/hashes/kim_dev_mouse.txt @@ -0,0 +1,8 @@ +DevCCFv1_OntologyStructure.xlsx 505d518b1a3b638475782150acf7750357244039585e1cc52aa48c79af005b87 +E11.5.zip 3dcdb118ae84f1b7935e66557e8be8139f9f0a489b34e6b339e971bb7dba9ff1 +E13.5.zip 8a795213be03d251198e88a4b48b94e8249049ec16c6b4e4852677ec1f92d7f6 +E15.5.zip 55fef88e7cd25b46b6d41415a2ad522301b0d4c72da48e1c90cf7a7b15d47cdf +E18.5.zip 7575614200087ce2012b226cc630ac1d826f3f75a89f66b4f4cc190639fd69fb +P04.zip 269812a1dfb41b6b25f817d1f250f1292370b75abd207e293e71193633ac4ba0 +P14.zip 79b6ea12c69fe8aebf643bac038747f42d1fd3b7ce4d9b593ae52b0b542c6c61 +P56.zip ec5cd5da9747a11fb3e266bdbde9767dc90bf1b457f6c030ed9a06f122be886a