diff --git a/.github/workflows/build_docs.yml b/.github/workflows/build_docs.yml index 857afe5..c8bc78d 100644 --- a/.github/workflows/build_docs.yml +++ b/.github/workflows/build_docs.yml @@ -7,17 +7,19 @@ on: workflow_dispatch: env: - DEB_PYTHON_INSTALL_LAYOUT: deb_system + HDF5_MPI: "ON" + HDF5_DIR: "/usr/local/" DISPLAY: ":99.0" - CI: 1 + DEB_PYTHON_INSTALL_LAYOUT: deb_system + LIBGL_ALWAYS_SOFTWARE: 1 jobs: build: runs-on: ubuntu-22.04 + container: ghcr.io/fenics/dolfinx/lab:stable env: PUBLISH_DIR: ./_build/html - DISPLAY: ":99.0" PYVISTA_TRAME_SERVER_PROXY_PREFIX: "/proxy/" PYVISTA_TRAME_SERVER_PROXY_ENABLED: "True" PYVISTA_OFF_SCREEN: false @@ -28,12 +30,7 @@ jobs: uses: actions/checkout@v4 - name: Install dependencies for pyvista - run: sudo apt-get update && sudo apt-get install -y libgl1-mesa-dev xvfb - - - name: Setup python - uses: actions/setup-python@v5 - with: - python-version: 3.12 + run: apt-get update && apt-get install -y libxrender1 xvfb - name: Install dependencies run: python3 -m pip install ".[docs]" diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index f3a1464..396b2ac 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -14,27 +14,44 @@ jobs: fail-fast: false matrix: python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"] - os: [ubuntu-latest, windows-latest, macos-latest] + os: [ubuntu-latest, macos-latest] steps: - uses: actions/checkout@v4 - - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v5 + + - name: Cache conda + uses: actions/cache@v3 + env: + # Increase this value to reset cache if environment.yml has not changed + CACHE_NUMBER: 0 + with: + path: ~/conda_pkgs_dir + key: + ${{ runner.os }}-conda-${{ env.CACHE_NUMBER }}-py-${{ matrix.python-version }}${{ + hashFiles('environment.yml') }} + + - uses: conda-incubator/setup-miniconda@v3 with: python-version: ${{ matrix.python-version }} - allow-prereleases: true + auto-update-conda: false + channels: conda-forge + activate-environment: mri2mesh + environment-file: environment.yml - name: Install mri2mesh + shell: bash -el {0} run: | python -m pip install -e ".[test]" - name: Test with pytest + shell: bash -el {0} run: | python -m pytest --cov=mri2mesh --cov-report html --cov-report xml --cov-report term-missing -v - name: Coverage report if: matrix.python-version == '3.10' && matrix.os == 'ubuntu-latest' + shell: bash -el {0} run: | python3 -m coverage report | sed 's/^/ /' >> $GITHUB_STEP_SUMMARY python3 -m coverage json diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index e963a88..98a7b35 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -29,6 +29,7 @@ repos: - id: mypy files: ^src/|^tests/ args: ["--config-file", "pyproject.toml"] + additional_dependencies: ['types-toml', 'types-PyYAML'] - repo: https://github.com/kynan/nbstripout diff --git a/_toc.yml b/_toc.yml index bc7a9db..de2f8a3 100644 --- a/_toc.yml +++ b/_toc.yml @@ -5,6 +5,7 @@ parts: - caption: Examples chapters: - file: "examples/synthseg.ipynb" + - file: "examples/cli.ipynb" - caption: Community chapters: - file: "CONTRIBUTING" diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..4463aac --- /dev/null +++ b/environment.yml @@ -0,0 +1,20 @@ +name: mri2mesh +channels: + - conda-forge +dependencies: + - numpy + - fenics-dolfinx + - ipython + - matplotlib + - pyvista + - scikit-image + - meshio + - h5py + - seaborn + - nibabel + - imageio-ffmpeg + - napari + - jupyter + - pip + - pip: + - wildmeshing diff --git a/examples/cli.ipynb b/examples/cli.ipynb new file mode 100644 index 0000000..29222d9 --- /dev/null +++ b/examples/cli.ipynb @@ -0,0 +1,356 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Command line interface\n", + "\n", + "The package comes with a command line interface called `mri2mesh`.\n", + "\n", + "*Note: In this documentation we will start the commands with `!` but this is only to make it run. In the terminal you should ommit the `!`*\n", + "\n", + "You can list the different commands with `--help`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "!mri2mesh --help" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Listing segmentation labels\n", + "\n", + "One simeple thing you can do is to list the segmentation labels for different segmentation software. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "!mri2mesh labels --help" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "For example for [synthseg](https://github.com/BBillot/SynthSeg) you can do" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "!mri2mesh labels synthseg" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "## Generating surfaces\n", + "\n", + "To see all options for generating surfaces you can again pass the help command" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "!mri2mesh surface --help" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "Let ut create an indealized brain and put the surfaces in a new folder" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "!mri2mesh surface idealized --help" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "!rm -rf idealized-brain\n", + "!mri2mesh surface idealized -o idealized-brain" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "Let us see which files that was generated" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "!ls idealized-brain" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "We can for example take a look at the parameters the was used" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "!cat idealized-brain/surface_parameters.json" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "## Generating meshes\n", + "\n", + "We will now show how to create a mesh from the surfaces" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "!mri2mesh mesh --help" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "We will start by creating a mesh with the `create` command" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "!mri2mesh mesh create --help" + ] + }, + { + "cell_type": "markdown", + "id": "19", + "metadata": {}, + "source": [ + "We see that this command takes a config file as input. We can generate a template for this config file using the `template` command" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20", + "metadata": {}, + "outputs": [], + "source": [ + "!mri2mesh mesh template ." + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "Let us have a look at this file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "!cat config.json" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "Most of these key-value pairs are paramerters to the meshing algorithm. The `csg_tree` describes how to combine the surfaces into a mesh and also create proper subdomains. In our case we would like to take the union of all surfaces. We can do this by first taking the union of `LV.ply` and `V34.ply` into a right component and `skull.ply` and `parenchyma_incl_ventr.ply` into a left component, and then finally take the union of the left and right component. In the template however, we need to update the paths because the surfaces are located in the folder `idealized-brain`. Let us also change the output directory to the same folder." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "import json\n", + "\n", + "config_fname = Path(\"config.json\")\n", + "config = json.loads(config_fname.read_text())\n", + "surface_folder = \"idealized-brain\"\n", + "config[\"outdir\"] = surface_folder\n", + "csg_tree = {\n", + " \"operation\": \"union\",\n", + " \"right\": {\n", + " \"operation\": \"union\",\n", + " \"left\": f\"{surface_folder}/LV.ply\",\n", + " \"right\": f\"{surface_folder}/V34.ply\"\n", + " },\n", + " \"left\": {\n", + " \"operation\": \"union\",\n", + " \"left\": f\"{surface_folder}/skull.ply\",\n", + " \"right\": f\"{surface_folder}/parenchyma_incl_ventr.ply\"\n", + " }\n", + "}\n", + "config[\"csg_tree\"] = csg_tree\n", + "config_fname.write_text(json.dumps(config))" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": {}, + "source": [ + "We can now create the mesh by passing in the config file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26", + "metadata": {}, + "outputs": [], + "source": [ + "!mri2mesh mesh create config.json" + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": {}, + "source": [ + "Let us have a look at the mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28", + "metadata": {}, + "outputs": [], + "source": [ + "from mpi4py import MPI\n", + "import dolfinx\n", + "\n", + "with dolfinx.io.XDMFFile(MPI.COMM_WORLD, f\"{surface_folder}/mesh.xdmf\", \"r\") as xdmf:\n", + " mesh = xdmf.read_mesh(name=\"mesh\")\n", + " cell_tags = xdmf.read_meshtags(mesh, name=\"cell_tags\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29", + "metadata": {}, + "outputs": [], + "source": [ + "import pyvista as pv \n", + "\n", + "pv.start_xvfb()\n", + "vtk_mesh = dolfinx.plot.vtk_mesh(mesh, cell_tags.dim, cell_tags.indices)\n", + "bgrid = pv.UnstructuredGrid(*vtk_mesh)\n", + "bgrid.cell_data[\"Cell tags\"] = cell_tags.values\n", + "bgrid.set_active_scalars(\"Cell tags\")\n", + "p = pv.Plotter(window_size=[800, 800])\n", + "p.add_mesh_clip_plane(bgrid)\n", + "if not pv.OFF_SCREEN:\n", + " p.show()\n", + "else:\n", + " figure = p.screenshot(\"idealized_brain.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/synthseg.ipynb b/examples/synthseg.ipynb index c886d60..6fd3591 100644 --- a/examples/synthseg.ipynb +++ b/examples/synthseg.ipynb @@ -22,10 +22,7 @@ "from pathlib import Path\n", "import mri2mesh\n", "import numpy as np\n", - "import pyvista as pv\n", - "\n", - "pv.start_xvfb()\n", - "# pv.set_jupyter_backend('trame')" + "import pyvista as pv" ] }, { @@ -53,7 +50,10 @@ "metadata": {}, "outputs": [], "source": [ + "pv.start_xvfb()\n", "# Path to the Nifty file\n", + "outdir = Path(\"results-synthseg\")\n", + "outdir.mkdir(exist_ok=True)\n", "path = Path(\"201_t13d_synthseg.nii.gz\")\n", "mri2mesh.viz.volume_clip.main(path)" ] @@ -298,6 +298,27 @@ "cell_type": "markdown", "id": "28", "metadata": {}, + "source": [ + "Let us also smooth the surface and save it as a `.ply` file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29", + "metadata": {}, + "outputs": [], + "source": [ + "par_surf = par_surf.smooth_taubin(n_iter=20, pass_band=0.05)\n", + "par_surf = par_surf.clip_closed_surface(normal=(0, 0, 1), origin=(0, 0, 1))\n", + "par_surf.compute_normals(inplace=True, flip_normals=False)\n", + "pv.save_meshio(outdir / \"par.ply\", par_surf)" + ] + }, + { + "cell_type": "markdown", + "id": "30", + "metadata": {}, "source": [ "## Extracting surfaces of the lateral ventricles \n", "\n", @@ -307,7 +328,7 @@ { "cell_type": "code", "execution_count": null, - "id": "29", + "id": "31", "metadata": {}, "outputs": [], "source": [ @@ -317,7 +338,7 @@ }, { "cell_type": "markdown", - "id": "30", + "id": "32", "metadata": {}, "source": [ "We can also plot the isosurface with pyvista" @@ -326,7 +347,7 @@ { "cell_type": "code", "execution_count": null, - "id": "31", + "id": "33", "metadata": {}, "outputs": [], "source": [ @@ -339,7 +360,7 @@ }, { "cell_type": "markdown", - "id": "32", + "id": "34", "metadata": {}, "source": [ "We can now generate a surface of this mask using `pyvista`" @@ -348,7 +369,7 @@ { "cell_type": "code", "execution_count": null, - "id": "33", + "id": "35", "metadata": {}, "outputs": [], "source": [ @@ -357,7 +378,7 @@ }, { "cell_type": "markdown", - "id": "34", + "id": "36", "metadata": {}, "source": [ "Let us plot the surface with `pyvista`" @@ -366,7 +387,7 @@ { "cell_type": "code", "execution_count": null, - "id": "35", + "id": "37", "metadata": {}, "outputs": [], "source": [ @@ -378,7 +399,7 @@ }, { "cell_type": "markdown", - "id": "36", + "id": "38", "metadata": {}, "source": [ "We see that the surface is not very smooth, but we can use the `smooth_taubin` method to smooth it\n" @@ -387,7 +408,7 @@ { "cell_type": "code", "execution_count": null, - "id": "37", + "id": "39", "metadata": {}, "outputs": [], "source": [ @@ -396,14 +417,6 @@ "plotter.add_mesh(surf_lateral_ventricles)\n", "plotter.show()" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "38", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { @@ -427,7 +440,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.4" + "version": "3.12.9" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index 44d2dd3..abce634 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,19 +5,20 @@ build-backend = "setuptools.build_meta" [project] name = "mri2mesh" version = "0.1.0" -dependencies = ["pyvista", "numpy", "matplotlib", "nibabel", "scikit-image", "scipy", ] +dependencies = ["pyvista", "numpy", "matplotlib", "nibabel", "scikit-image", "scipy", "meshio"] classifiers = [ "License :: OSI Approved :: MIT License", "Programming Language :: Python :: 3", ] description = "Tool for converting labeled MRI data to a mesh" authors = [ - {name = "Marius Causemann", email = "mariusca@simula.no"}, {name = "Henrik Finsberg", email = "henriknf@simula.no"}, + {name = "Marius Causemann", email = "mariusca@simula.no"}, ] -keywords = ["mri", "fem", "brain", "fluid", "mechanics", "meshing"] -urls = {Homepage = "https://github.com/kent-and/mri2mesh"} +license = {text = "MIT"} +keywords = ["mri", "fem", "brain", "meshing"] +urls = {Homepage = "https://github.com/scientificcomputing/mri2mesh"} [project.readme] @@ -29,8 +30,9 @@ content-type = "text/markdown" mri2mesh = "mri2mesh.cli:main" [project.optional-dependencies] -test = ["pytest", "pytest-cov"] -docs = ["pyvista[jupyter]", "jupyter-book"] +mesh = ["wildmeshing", "h5py"] +test = ["pytest", "pytest-cov", "mri2mesh[mesh]"] +docs = ["pyvista[jupyter]", "jupyter-book", "mri2mesh[mesh]"] @@ -122,7 +124,7 @@ tag = true sign_tags = false tag_name = "v{new_version}" tag_message = "Bump version: {current_version} → {new_version}" -current_version = "1.1.3" +current_version = "0.1.0" [[tool.bumpversion.files]] diff --git a/src/mri2mesh/__init__.py b/src/mri2mesh/__init__.py index 641960b..49ae005 100644 --- a/src/mri2mesh/__init__.py +++ b/src/mri2mesh/__init__.py @@ -1,9 +1,19 @@ +from importlib.metadata import metadata + from . import viz from . import morphology from . import cli from . import vtk_utils from . import surface +from . import mesh from . import reader from .reader import read -__all__ = ["viz", "morphology", "cli", "vtk_utils", "surface", "reader", "read"] +__all__ = ["viz", "morphology", "cli", "vtk_utils", "surface", "reader", "read", "mesh"] + +meta = metadata("mri2mesh") +__version__ = meta["Version"] +__author__ = meta["Author-email"] +__license__ = meta["License"] +__email__ = meta["Author-email"] +__program_name__ = meta["Name"] diff --git a/src/mri2mesh/cli.py b/src/mri2mesh/cli.py index 1bed5b6..3de48e8 100644 --- a/src/mri2mesh/cli.py +++ b/src/mri2mesh/cli.py @@ -1,7 +1,8 @@ import logging import argparse +from typing import Sequence, Optional -from . import viz, surface +from . import viz, surface, mesh, segmentation_labels def setup_parser(): @@ -27,6 +28,16 @@ def setup_parser(): # Surface generation parser surface_parser = subparsers.add_parser("surface", help="Generate surfaces") surface.add_surface_parser(surface_parser) + + # Mesh generation parser + mesh_parser = subparsers.add_parser("mesh", help="Generate meshes") + mesh.add_mesh_parser(mesh_parser) + + label_parser = subparsers.add_parser("labels", help="List segmentation labels") + label_parser.add_argument( + "name", type=str, choices=["synthseg", "neuroquant"], help="Name of the labels to display" + ) + return parser @@ -35,8 +46,8 @@ def _disable_loggers(): logging.getLogger(libname).setLevel(logging.WARNING) -def dispatch(parser: argparse.ArgumentParser) -> int: - args = vars(parser.parse_args()) +def dispatch(parser: argparse.ArgumentParser, argv: Optional[Sequence[str]] = None) -> int: + args = vars(parser.parse_args(argv)) logging.basicConfig(level=logging.DEBUG if args.pop("verbose") else logging.INFO) _disable_loggers() @@ -54,6 +65,10 @@ def dispatch(parser: argparse.ArgumentParser) -> int: viz.dispatch(args.pop("viz-command"), args) elif command == "surface": surface.dispatch(args.pop("surface-command"), args) + elif command == "mesh": + mesh.dispatch(args.pop("mesh-command"), args) + elif command == "labels": + segmentation_labels.dispatch(args.pop("name")) else: logger.error(f"Unknown command {command}") parser.print_help() @@ -64,6 +79,6 @@ def dispatch(parser: argparse.ArgumentParser) -> int: return 0 -def main() -> int: +def main(argv: Optional[Sequence[str]] = None) -> int: parser = setup_parser() - return dispatch(parser) + return dispatch(parser, argv) diff --git a/src/mri2mesh/mesh/__init__.py b/src/mri2mesh/mesh/__init__.py new file mode 100644 index 0000000..9650ae2 --- /dev/null +++ b/src/mri2mesh/mesh/__init__.py @@ -0,0 +1,59 @@ +import typing +import logging +import argparse +from pathlib import Path + +from . import basic, idealized_brain +from .basic import create_mesh + +__all__ = ["basic", "idealized_brain", "create_mesh"] + +logger = logging.getLogger(__name__) + + +def add_mesh_parser(parser: argparse.ArgumentParser) -> None: + subparsers = parser.add_subparsers(dest="mesh-command") + + template_parser = subparsers.add_parser( + "template", help="Create a sample configuration file called" + ) + template_parser.add_argument("outdir", type=Path, help="Output directory") + template_parser.add_argument( + "--name", type=str, default="config.json", help="Name of the configuration file" + ) + + create_parser = subparsers.add_parser("create", help="Create a mesh") + create_parser.add_argument("filename", type=Path, help="Path to the configuration file") + + convert_parser = subparsers.add_parser("convert", help="Convert mesh to dolfinx") + convert_parser.add_argument("mesh_dir", type=Path, help="Directory containing mesh files") + + idealized_parser = subparsers.add_parser( + "idealized", + help="Generate idealized surface", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + idealized_brain.add_arguments(idealized_parser) + + +def dispatch(command, args: dict[str, typing.Any]) -> int: + if command == "template": + basic.generate_sameple_config(**args) + + elif command == "create": + mesh_dir = basic.create_mesh_from_config(**args) + try: + basic.convert_mesh_dolfinx(mesh_dir=mesh_dir) + except ImportError: + logger.debug("dolfinx not installed, skipping conversion to dolfinx") + + elif command == "convert": + basic.convert_mesh_dolfinx(**args) + + elif command == "idealized": + idealized_brain.main(**args) + + else: + raise ValueError(f"Unknown command {command}") + + return 0 diff --git a/src/mri2mesh/mesh/basic.py b/src/mri2mesh/mesh/basic.py new file mode 100644 index 0000000..d996e82 --- /dev/null +++ b/src/mri2mesh/mesh/basic.py @@ -0,0 +1,238 @@ +import typing +import json +import logging +import pprint +from pathlib import Path +import numpy as np + +logger = logging.getLogger(__name__) + + +def sample_config(): + tree = { + "operation": "union", + "right": { + "operation": "union", + "left": "surfaces/LV.ply", + "right": "surfaces/V34.ply", + }, + "left": { + "operation": "union", + "left": "surfaces/skull.ply", + "right": "surfaces/parenchyma_incl_ventr.ply", + }, + } + + config = { + "outdir": "mesh", + "epsilon": 0.0009, + "edge_length_r": 0.015, + "coarsen": True, + "stop_quality": 8, + "max_its": 30, + "loglevel": 10, + "disable_filtering": False, + "use_floodfill": False, + "smooth_open_boundary": False, + "manifold_surface": False, + "csg_tree": tree, + } + return config + + +def default_config() -> dict[str, typing.Any]: + config = sample_config() + config["csg_tree"] = {} + return config + + +def generate_sameple_config(outdir: Path, name: str = "config.json"): + logger.info("Generating sample configuration file") + config = sample_config() + path = outdir / name + if path.suffix == ".json": + path.write_text(json.dumps(config, indent=2)) + elif path.suffix == ".toml": + import toml + + path.write_text(toml.dumps(config)) + elif path.suffix == ".yaml": + import yaml + + path.write_text(yaml.dump(config)) + else: + raise ValueError(f"Unknown file extension {path.suffix}") + logger.info("Configuration file written to %s", path) + + +def read_config(path: Path) -> dict[str, typing.Any]: + logger.info("Reading configuration file %s", path) + if path.suffix == ".json": + return json.loads(path.read_text()) + elif path.suffix == ".toml": + import toml + + return toml.loads(path.read_text()) + elif path.suffix == ".yaml": + import yaml + + return yaml.load(path.read_text(), Loader=yaml.FullLoader) + else: + raise ValueError(f"Unknown file extension {path.suffix}") + + +class CSGTree(typing.TypedDict): + operation: str + left: typing.Union[str, "CSGTree"] + right: typing.Union[str, "CSGTree"] + + +def create_mesh_from_config(filename: Path) -> Path: + if not filename.exists(): + raise FileNotFoundError(f"File {filename} not found") + config = default_config() + config.update(read_config(filename)) + create_mesh(**config) + return Path(config["outdir"]) + + +def create_mesh( + outdir: typing.Union[Path, str], + csg_tree: CSGTree, + epsilon: float = 0.0009, + edge_length_r: float = 0.015, + skip_simplify: bool = False, + coarsen: bool = True, + stop_quality: int = 8, + max_its: int = 30, + loglevel: int = 10, + disable_filtering: bool = False, + use_floodfill: bool = False, + smooth_open_boundary: bool = False, + manifold_surface: bool = False, +): + logger.info("Creating mesh, with the following configuration") + params = { + "outdir": str(outdir), + "csg_tree": csg_tree, + "epsilon": epsilon, + "edge_length_r": edge_length_r, + "skip_simplify": skip_simplify, + "coarsen": coarsen, + "stop_quality": stop_quality, + "max_its": max_its, + "loglevel": loglevel, + "disable_filtering": disable_filtering, + "use_floodfill": use_floodfill, + "smooth_open_boundary": smooth_open_boundary, + "manifold_surface": manifold_surface, + } + logger.info(pprint.pformat(params)) + + import wildmeshing as wm + import meshio + import pyvista as pv + + outdir = Path(outdir) + outdir.mkdir(parents=True, exist_ok=True) + (outdir / "mesh_params.json").write_text(json.dumps(params, indent=2)) + + tetra = wm.Tetrahedralizer( + epsilon=epsilon, + edge_length_r=edge_length_r, + coarsen=coarsen, + stop_quality=stop_quality, + max_its=max_its, + skip_simplify=skip_simplify, + ) + tetra.set_log_level(loglevel) + + tetra.load_csg_tree(json.dumps(csg_tree)) + tetra.tetrahedralize() + point_array, cell_array, marker = tetra.get_tet_mesh( + all_mesh=disable_filtering, + smooth_open_boundary=smooth_open_boundary, + floodfill=use_floodfill, + manifold_surface=manifold_surface, + correct_surface_orientation=True, + ) + coords, connections = tetra.get_tracked_surfaces() + + tetra_mesh = meshio.Mesh( + point_array, [("tetra", cell_array)], cell_data={"cell_tags": [marker.ravel()]} + ) + tetra_mesh_pv = pv.from_meshio(tetra_mesh).clean() + pv.save_meshio(outdir / "tetra_mesh.xdmf", tetra_mesh_pv) + + for i, coord in enumerate(coords): + np.save(outdir / f"coords_{i}.npy", coord) + + for i, conn in enumerate(connections): + np.save(outdir / f"connections_{i}.npy", conn) + + +def convert_mesh_dolfinx(mesh_dir: Path): + logger.info("Converting mesh to dolfinx in %s", mesh_dir) + from mpi4py import MPI + from scipy.spatial.distance import cdist + import dolfinx + + threshold = 1.0 + fdim = 2 + + coords = [] + for path in sorted(mesh_dir.glob("coords_*.npy"), key=lambda x: int(x.stem.split("_")[-1])): + coords.append(np.load(path)) + logger.debug(f"Found {len(coords)} coordinates") + + connections = [] + for path in sorted( + mesh_dir.glob("connections_*.npy"), key=lambda x: int(x.stem.split("_")[-1]) + ): + connections.append(np.load(path)) + logger.debug(f"Found {len(connections)} connections") + + assert len(connections) == len(coords) + + logger.debug("Loading mesh") + comm = MPI.COMM_WORLD + with dolfinx.io.XDMFFile(comm, mesh_dir / "tetra_mesh.xdmf", "r") as xdmf: + mesh = xdmf.read_mesh(name="Grid") + cell_tags = xdmf.read_meshtags(mesh, name="Grid") + + logger.debug("Mesh loaded") + + facets = [] + values = [] + for i, coord in enumerate(coords, start=1): + logger.debug(f"Processing coord {i}") + + def locator(x): + # Find the distance to all coordinates + distances = cdist(x.T, coord) + # And return True is they are close + return np.any(distances < threshold, axis=1) + + f = dolfinx.mesh.locate_entities_boundary(mesh, dim=fdim, marker=locator) + v = np.full(f.shape[0], i, dtype=np.int32) + facets.append(f) + values.append(v) + + logger.debug("Create meshtags") + facet_tags = dolfinx.mesh.meshtags( + mesh, + fdim, + np.hstack(facets), + np.hstack(values), + ) + facet_tags.name = "facet_tags" + cell_tags.name = "cell_tags" + mesh.name = "mesh" + + logger.debug("Save files") + with dolfinx.io.XDMFFile(comm, mesh_dir / "mesh.xdmf", "w") as xdmf: + xdmf.write_mesh(mesh) + xdmf.write_meshtags(facet_tags, mesh.geometry) + xdmf.write_meshtags(cell_tags, mesh.geometry) + + logger.info("Mesh saved to %s", mesh_dir / "mesh.xdmf") diff --git a/src/mri2mesh/mesh/idealized_brain.py b/src/mri2mesh/mesh/idealized_brain.py new file mode 100644 index 0000000..21660f6 --- /dev/null +++ b/src/mri2mesh/mesh/idealized_brain.py @@ -0,0 +1,224 @@ +from pathlib import Path +import logging + +logger = logging.getLogger(__name__) + + +def add_arguments(parser): + parser.add_argument( + "-o", + "--outdir", + type=Path, + required=True, + help="Output directory", + ) + parser.add_argument( + "--r", + type=float, + default=0.1, + help="Radius of the brain", + ) + parser.add_argument( + "--parenchyma-factor", + type=float, + default=0.6, + help="Parenchyma factor", + ) + parser.add_argument( + "--lv-factor", + type=float, + default=0.2, + help="Left ventricle factor", + ) + parser.add_argument( + "--v34-center-factor", + type=float, + default=0.4, + help="Ventricle 3/4 center factor", + ) + parser.add_argument( + "--v34-height-factor", + type=float, + default=0.42, + help="Ventricle 3/4 height factor", + ) + parser.add_argument( + "--v34-radius-factor", + type=float, + default=0.06, + help="Ventricle 3/4 radius factor", + ) + parser.add_argument( + "--skull-x0", + type=float, + default=0.0103891, + help="Skull x0", + ) + parser.add_argument( + "--skull-x1", + type=float, + default=0.155952, + help="Skull x1", + ) + parser.add_argument( + "--skull-y0", + type=float, + default=0.0114499, + help="Skull y0", + ) + parser.add_argument( + "--skull-y1", + type=float, + default=0.173221, + help="Skull y1", + ) + parser.add_argument( + "--skull-z0", + type=float, + default=0.001, + help="Skull z0", + ) + parser.add_argument( + "--skull-z1", + type=float, + default=0.154949, + help="Skull z1", + ) + parser.add_argument( + "--epsilon", + type=float, + default=0.0009, + help="Epsilon", + ) + parser.add_argument( + "--edge-length-r", + type=float, + default=0.015, + help="Edge length", + ) + parser.add_argument( + "--skip-simplify", + action="store_true", + help="Skip simplification", + ) + parser.add_argument( + "--coarsen", + action="store_true", + help="Coarsen", + ) + parser.add_argument( + "--stop-quality", + type=int, + default=8, + help="Stop quality", + ) + parser.add_argument( + "--max-its", + type=int, + default=30, + help="Max iterations", + ) + parser.add_argument( + "--loglevel", + type=int, + default=10, + help="Log level", + ) + parser.add_argument( + "--disable-filtering", + action="store_true", + help="Disable filtering", + ) + parser.add_argument( + "--use-floodfill", + action="store_true", + help="Use floodfill", + ) + parser.add_argument( + "--smooth-open-boundary", + action="store_true", + help="Smooth open boundary", + ) + parser.add_argument( + "--manifold-surface", + action="store_true", + help="Manifold surface", + ) + + +def main( + outdir: Path, + r: float = 0.1, + parenchyma_factor: float = 0.6, + lv_factor: float = 0.2, + v34_center_factor: float = 0.4, + v34_height_factor: float = 0.42, + v34_radius_factor: float = 0.06, + skull_x0: float = 0.0103891, + skull_x1: float = 0.155952, + skull_y0: float = 0.0114499, + skull_y1: float = 0.173221, + skull_z0: float = 0.001, + skull_z1: float = 0.154949, + epsilon: float = 0.0009, + edge_length_r: float = 0.015, + skip_simplify: bool = False, + coarsen: bool = True, + stop_quality: int = 8, + max_its: int = 30, + loglevel: int = 10, + disable_filtering: bool = False, + use_floodfill: bool = False, + smooth_open_boundary: bool = False, + manifold_surface: bool = False, +) -> None: + logger.info("Generating idealized brain surface") + from ..surface.idealized_brain import main as main_surface + + main_surface( + outdir, + r=r, + parenchyma_factor=parenchyma_factor, + lv_factor=lv_factor, + v34_center_factor=v34_center_factor, + v34_height_factor=v34_height_factor, + v34_radius_factor=v34_radius_factor, + skull_x0=skull_x0, + skull_x1=skull_x1, + skull_y0=skull_y0, + skull_y1=skull_y1, + skull_z0=skull_z0, + skull_z1=skull_z1, + ) + from .basic import create_mesh, CSGTree, convert_mesh_dolfinx + + csg_tree: CSGTree = { + "operation": "union", + "right": { + "operation": "union", + "left": f"{outdir}/LV.ply", + "right": f"{outdir}/V34.ply", + }, + "left": { + "operation": "union", + "left": f"{outdir}/skull.ply", + "right": f"{outdir}/parenchyma_incl_ventr.ply", + }, + } + + create_mesh( + outdir, + csg_tree=csg_tree, + epsilon=epsilon, + edge_length_r=edge_length_r, + skip_simplify=skip_simplify, + coarsen=coarsen, + stop_quality=stop_quality, + max_its=max_its, + loglevel=loglevel, + disable_filtering=disable_filtering, + use_floodfill=use_floodfill, + smooth_open_boundary=smooth_open_boundary, + manifold_surface=manifold_surface, + ) + convert_mesh_dolfinx(mesh_dir=outdir) diff --git a/src/mri2mesh/morphology.py b/src/mri2mesh/morphology.py index c1ac793..dabbd32 100644 --- a/src/mri2mesh/morphology.py +++ b/src/mri2mesh/morphology.py @@ -3,7 +3,7 @@ import numpy.typing as npt import skimage.morphology as skim import skimage -import scipy.ndimage as ndi + logger = logging.getLogger(__name__) @@ -34,6 +34,8 @@ def get_closest_point(a: npt.NDArray, b: np.ndarray) -> tuple[np.int64, ...]: >>> assert np.allclose(a_b_index, (5, 5, 5)) """ + import scipy.ndimage as ndi + dist = ndi.distance_transform_edt(np.logical_not(a)) assert isinstance(dist, np.ndarray) dist[np.logical_not(b)] = np.inf diff --git a/src/mri2mesh/reader.py b/src/mri2mesh/reader.py index 291d0ae..ddd7fa8 100644 --- a/src/mri2mesh/reader.py +++ b/src/mri2mesh/reader.py @@ -2,12 +2,10 @@ from pathlib import Path import logging import typing - import numpy as np import numpy.typing as npt from dataclasses import dataclass, field -import nibabel as nib from .segmentation_labels import NEUROQUANT_LABELS, SYNTHSEG_LABELS @@ -54,6 +52,8 @@ def read( label_name: typing.Literal["synthseg", "neuroquant"] = "synthseg", padding: int = 5, ) -> Segmentation: + import nibabel as nib + logger.info(f"Loading segmentation from {input}") seg = nib.load(input) img = np.pad(seg.get_fdata(), padding) diff --git a/src/mri2mesh/segmentation_labels.py b/src/mri2mesh/segmentation_labels.py index 8698c13..8be0d12 100644 --- a/src/mri2mesh/segmentation_labels.py +++ b/src/mri2mesh/segmentation_labels.py @@ -1,3 +1,20 @@ +import logging + +logger = logging.getLogger(__name__) + + +def dispatch(name: str) -> None: + logger.info(f"Listing segmentation labels for {name}") + if name == "synthseg": + for label, values in SYNTHSEG_LABELS.items(): + print(f"{label}: {values}") + elif name == "neuroquant": + for label, values in NEUROQUANT_LABELS.items(): + print(f"{label}: {values}") + else: + raise ValueError(f"Unknown segmentation labels {name}") + + SYNTHSEG_LABELS = { "BACKGROUND": [0], "LEFT_CEREBRAL_WHITE_MATTER": [2], diff --git a/src/mri2mesh/surface/__init__.py b/src/mri2mesh/surface/__init__.py index 3d9a974..b7b09d0 100644 --- a/src/mri2mesh/surface/__init__.py +++ b/src/mri2mesh/surface/__init__.py @@ -2,22 +2,34 @@ import argparse import logging -from . import parenchyma +from . import parenchyma, idealized_brain logger = logging.getLogger(__name__) def add_surface_parser(parser: argparse.ArgumentParser) -> None: - subparsers = parser.add_subparsers(dest="surface-command") + subparsers = parser.add_subparsers( + dest="surface-command", + ) parenchyma_parser = subparsers.add_parser("parenchyma", help="Generate parenchyma surface") parenchyma.add_arguments(parenchyma_parser) + idealized_parser = subparsers.add_parser( + "idealized", + help="Generate idealized surface", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + idealized_brain.add_arguments(idealized_parser) + def dispatch(command, args: dict[str, typing.Any]) -> int: if command == "parenchyma": parenchyma.main(**args) + elif command == "idealized": + idealized_brain.main(**args) + else: raise ValueError(f"Unknown command {command}") diff --git a/src/mri2mesh/surface/idealized_brain.py b/src/mri2mesh/surface/idealized_brain.py new file mode 100644 index 0000000..685238c --- /dev/null +++ b/src/mri2mesh/surface/idealized_brain.py @@ -0,0 +1,154 @@ +from pathlib import Path +import argparse +import logging +import json +import datetime +import numpy as np + +logger = logging.getLogger(__name__) + + +def add_arguments(parser: argparse.ArgumentParser) -> None: + parser.add_argument( + "-o", + "--outdir", + type=Path, + required=True, + help="Output directory", + ) + parser.add_argument( + "--r", + type=float, + default=0.1, + help="Radius of the brain", + ) + parser.add_argument( + "--parenchyma-factor", + type=float, + default=0.6, + help="Parenchyma factor", + ) + parser.add_argument( + "--lv-factor", + type=float, + default=0.2, + help="Left ventricle factor", + ) + parser.add_argument( + "--v34-center-factor", + type=float, + default=0.4, + help="Ventricle 3/4 center factor", + ) + parser.add_argument( + "--v34-height-factor", + type=float, + default=0.42, + help="Ventricle 3/4 height factor", + ) + parser.add_argument( + "--v34-radius-factor", + type=float, + default=0.06, + help="Ventricle 3/4 radius factor", + ) + parser.add_argument( + "--skull-x0", + type=float, + default=0.0103891, + help="Skull x0", + ) + parser.add_argument( + "--skull-x1", + type=float, + default=0.155952, + help="Skull x1", + ) + parser.add_argument( + "--skull-y0", + type=float, + default=0.0114499, + help="Skull y0", + ) + parser.add_argument( + "--skull-y1", + type=float, + default=0.173221, + help="Skull y1", + ) + parser.add_argument( + "--skull-z0", + type=float, + default=0.001, + help="Skull z0", + ) + parser.add_argument( + "--skull-z1", + type=float, + default=0.154949, + help="Skull z1", + ) + + +def main( + outdir: Path, + r: float = 0.1, + parenchyma_factor: float = 0.6, + lv_factor: float = 0.2, + v34_center_factor: float = 0.4, + v34_height_factor: float = 0.42, + v34_radius_factor: float = 0.06, + skull_x0: float = 0.0103891, + skull_x1: float = 0.155952, + skull_y0: float = 0.0114499, + skull_y1: float = 0.173221, + skull_z0: float = 0.001, + skull_z1: float = 0.154949, +) -> None: + import pyvista as pv + + logger.info("Creating idealized brain surface in %s", outdir) + outdir.mkdir(exist_ok=True, parents=True) + + z = np.array([0, 0, 1]) + skull = pv.Box((skull_x0, skull_x1, skull_y0, skull_y1, skull_z0, skull_z1)) + c = skull.center_of_mass() + par = pv.Sphere(r * parenchyma_factor, center=c) + LV = pv.Sphere(r * lv_factor, center=c) + V34 = pv.Cylinder( + c - v34_center_factor * r * z, + direction=z, + height=r * v34_height_factor, + radius=v34_radius_factor * r, + ) + ventricles = pv.merge([LV, V34]) + + for s, n in zip( + [V34, LV, par, skull, ventricles], + ["V34", "LV", "parenchyma_incl_ventr", "skull", "ventricles"], + ): + pv.save_meshio(outdir / f"{n}.ply", s) + + from .. import __version__ + + (outdir / "surface_parameters.json").write_text( + json.dumps( + { + "r": r, + "parenchyma_factor": parenchyma_factor, + "lv_factor": lv_factor, + "v34_center_factor": v34_center_factor, + "v34_height_factor": v34_height_factor, + "v34_radius_factor": v34_radius_factor, + "skull_x0": skull_x0, + "skull_x1": skull_x1, + "skull_y0": skull_y0, + "skull_y1": skull_y1, + "skull_z0": skull_z0, + "skull_z1": skull_z1, + "version": __version__, + "timestamp": datetime.datetime.now().isoformat(), + }, + indent=2, + ) + ) diff --git a/src/mri2mesh/viz/mpl_slice.py b/src/mri2mesh/viz/mpl_slice.py index 9121605..c148a6a 100644 --- a/src/mri2mesh/viz/mpl_slice.py +++ b/src/mri2mesh/viz/mpl_slice.py @@ -2,8 +2,6 @@ import argparse import typing from pathlib import Path -import nibabel as nib -import matplotlib.pyplot as plt import numpy as np @@ -58,6 +56,8 @@ def main( input: Path, output: Path, tag: str, cmap: str, add_colorbar: bool, nx: int, ny: int ) -> int: output.mkdir(parents=True, exist_ok=True) + import nibabel as nib + img = nib.load(input).get_fdata() plot_slices( img, @@ -82,6 +82,8 @@ def plot_slice( ny: int, slice: str, ): + import matplotlib.pyplot as plt + minval = np.min(img) maxval = np.max(img) fig, ax = plt.subplots(nx, ny, figsize=(20, 20)) diff --git a/tests/test_cli.py b/tests/test_cli.py new file mode 100644 index 0000000..054a41f --- /dev/null +++ b/tests/test_cli.py @@ -0,0 +1,37 @@ +import pytest + +from mri2mesh import cli + + +@pytest.mark.parametrize("name", ["synthseg", "neuroquant"]) +def test_labels_cli(name, capsys): + cli.main(["labels", name]) + captured = capsys.readouterr() + # Both labels have Background as the first label + assert captured.out.startswith("BACKGROUND: [0]") + + +def test_mesh_idealized_cli(tmp_path): + cli.main(["mesh", "idealized", "-o", str(tmp_path)]) + for name in [ + "connections_3.npy", + "coords_0.npy", + "coords_1.npy", + "connections_2.npy", + "connections_0.npy", + "coords_3.npy", + "coords_2.npy", + "connections_1.npy", + "mesh.h5", + "ventricles.ply", + "skull.ply", + "tetra_mesh.xdmf", + "mesh.xdmf", + "LV.ply", + "surface_parameters.json", + "tetra_mesh.h5", + "parenchyma_incl_ventr.ply", + "V34.ply", + "mesh_params.json", + ]: + assert (tmp_path / name).exists()