Merge pull request #88 from Hendrik-code/development_robert

Development robert

Author: robert-graf

.gitignore

@@ -165,4 +165,5 @@ tutorials/*PixelPandemonium/*
 tutorials/dataset-PixelPandemonium/*
 *.html
 _*.py
-dicom_select
+dicom_select
+examples

README.md

@@ -87,23 +87,36 @@ nii.shape #shape
 ```
 ### Stitching
 Python function and script for arbitrary image stitching. [See Details](TPTBox/stitching/)
+
+![Example of a stitching](TPTBox/stitching/stitching.jpg)
 ### Spineps and Points of Interests (POI) integration
+
+![Example of two lumbar vertebrae. The left example is derived from 1 mm isotropic CT, the right from sagittal MRI with a resolution of 3.3 mm in the left–right direction. Top row: Subregion of the vertebra used for analysis. Middle row: Extreme points. Bottom row: Corpus edge and ligamentum flavum points.](TPTBox/images/poi_preview.png)
 For our Spine segmentation pipline follow the installation of [SPINEPS](https://github.com/Hendrik-code/spineps).
+Image Source: Rule-based Key-Point Extraction for MR-Guided Biomechanical Digital Twins of the Spine; 
+
+
 
 SPINEPS will produce two mask: instance and semantic labels. With these we can compute our POIs. There are either center of mass points or surface points with bioloical meaning. See [Validation of a Patient-Specific Musculoskeletal Model for Lumbar Load Estimation Generated by an Automated Pipeline From Whole Body CT](https://pubmed.ncbi.nlm.nih.gov/35898642/)
 ```python
-from TPTBox import NII, POI, Location, calc_poi_from_subreg_vert
+from TPTBox import NII, POI, Location, POI_Global, calc_poi_from_subreg_vert
 from TPTBox.core.vert_constants import v_name2idx
+from TPTBox.segmentation.spineps import run_spineps_single
+
+# This requires that spineps is installed
+output_paths = run_spineps_single(
+    "file-path-of_T2w.nii.gz",
+    model_semantic="t2w",
+    ignore_compatibility_issues=True,
+)
+out_spine = output_paths["out_spine"]
+out_vert = output_paths["out_vert"]
+semantic_nii = NII.load(out_spine, seg=True)
+instance_nii = NII.load(out_vert, seg=True)
 
-p = "/dataset-DATASET/derivatives/A/"
-semantic_nii = NII.load(f"{p}sub-A_sequ-stitched_acq-sag_mod-T2w_seg-spine_msk.nii.gz", seg=True)
-instance_nii = NII.load(f"{p}sub-A_sequ-stitched_acq-sag_mod-T2w_seg-vert_msk.nii.gz", seg=True)
-poi_path = f"{p}sub-A_sequ-stitched_acq-sag_mod-T2w_seg-spine_ctd.json"
-poi = POI.load(poi_path)
 poi = calc_poi_from_subreg_vert(
     instance_nii,
     semantic_nii,
-    # buffer_file=poi_path,
     subreg_id=[
         Location.Vertebra_Full,
         Location.Arcus_Vertebrae,
@@ -157,18 +170,34 @@ poi = calc_poi_from_subreg_vert(
         Location.Vertebra_Direction_Right,
     ],
 )
-# poi.save(poi_path)
 poi = poi.round(2)
 print("Vertebra T4 Vertebra Corpus Center of mass:", poi[v_name2idx["T4"], Location.Vertebra_Corpus])
-# rescale/reorante like nii
+print("The id number of T4 Vertebra_Corpus is ", v_name2idx["T4"], Location.Vertebra_Corpus.value)
+
+# rescale/reorante local poi like nii
 poi_new = poi.reorient(("P", "I", "R")).rescale((1, 1, 1))
+# Local and global POIs can be rescaled to a target spacing with:
 poi_new = poi.resample_from_to(other_nii_or_poi)
 
+# local to global poi
+global_poi = poi.to_global(itk_coords=True)
+# You can save global pois in mrk.json format for import and editing in slicer.
+global_poi.save_mrk("FILE.mrk.json", glyphScale=3.0)
+# Import as a Markup in slicer; To make points editable you must click on the "lock" symbol under Markups - Control Points - Interaction
+
+# Save in our format:
+poi.save(poi_path)
+# Loading local/global Poi
+poi = POI.load(poi_path)
+poi = POI_Global.load(poi_path)
+
+
+
 ```
 
 
 ### Snapshot2D Spine
-
+![Snapshot2D Spine example](TPTBox/images/snp2D_example.png)
 The snapshot function automatically generates sag, cor, axial cuts in the center of a segmentation.
 
 ```python
@@ -188,17 +217,21 @@ create_snapshot(snp_path="snapshot.jpg", frames=[ct_frame, mr_frame])
 
 
 ### Snapshot3D
-
+![Snapshot3D example](TPTBox/images/snp3D_example.jpg)
 Requires additonal python packages: vtk fury xvfbwrapper
 
 ```python
-from TPTBox.mesh3D.snapshot3D import make_snapshot3D
+from TPTBox.mesh3D.snapshot3D import make_snapshot3D, make_snapshot3D_parallel
+
 # all segmentation; view give the rotation of an image
-make_snapshot3D("sub-101000_msk.nii.gz","snapshot3D.png",view=["A", "L", "P", "R"])
+make_snapshot3D("sub-101000_msk.nii.gz", "snapshot3D.png", view=["A", "L", "P", "R"])
 # Select witch segmentation per panel are chosen.
-make_snapshot3D("sub-101000_msk.nii.gz","snapshot3D_v2.png",view=["A"], ids_list=[[1,2],[3]])
+make_snapshot3D("sub-101000_msk.nii.gz", "snapshot3D_v2.png", view=["A"], ids_list=[[1, 2], [3]])
+# we proviede a implementation to process multiple images at the same time.
+make_snapshot3D_parallel(["a.nii.gz", "b.nii.gz"], ["snp_a.png", "snp_b.png"], view=["A"])
 ```
 
+<!---
 ### Logger
 
 ```python
@@ -219,3 +252,4 @@ TBD
 
 > [!IMPORTANT]
 > Importantly
+-->

TPTBox/core/bids_constants.py

@@ -152,6 +152,7 @@
     "logit",
     "localizer",
     "difference",
+    "labels",
 ]
 # https://bids-specification.readthedocs.io/en/stable/appendices/entity-table.html
 formats_relaxed = [*formats, "t2", "t1", "t2c", "t1c", "cta", "mr", "snapshot", "t1dixon", "dwi"]
@@ -175,6 +176,7 @@
 file_types = [
     "nii.gz",
     "json",
+    "mrk.json",
     "png",
     "jpg",
     "tsv",
@@ -190,6 +192,7 @@
     "xlsx",
     "bvec",
     "bval",
+    "html",
 ]
 # Description see: https://bids-specification.readthedocs.io/en/stable/99-appendices/09-entities.html
 

TPTBox/core/bids_files.py

@@ -3,6 +3,7 @@
 
 import json
 import os
+import random
 import sys
 import typing
 from collections.abc import Sequence
@@ -358,10 +359,14 @@ def add_file_2_subject(self, bids: BIDS_FILE | Path, ds=None) -> None:
         ) if self.verbose else None
         self.subjects[subject].add(bids)
 
-    def enumerate_subjects(self, sort=False) -> list[tuple[str, Subject_Container]]:
+    def enumerate_subjects(self, sort=False, shuffle=False) -> list[tuple[str, Subject_Container]]:
         # TODO Enumerate should put out numbers...
         if sort:
             return sorted(self.subjects.items())
+        if shuffle:
+            s = list(self.subjects.items())
+            random.shuffle(s)
+            return s
         return self.subjects.items()  # type: ignore
 
     def iter_subjects(self, sort=False) -> list[tuple[str, Subject_Container]]:
@@ -686,7 +691,7 @@ def get_changed_path(  # noqa: C901
         from_info=False,
         auto_add_run_id=False,
         additional_folder: str | None = None,
-        dataset_path: str | None = None,
+        dataset_path: str | Path | None = None,
         make_parent=False,
         no_sorting_mode: bool = False,
         non_strict_mode: bool = False,

TPTBox/core/nii_poi_abstract.py

@@ -7,6 +7,7 @@
 import nibabel as nib
 import nibabel.orientations as nio
 import numpy as np
+from scipy.spatial.transform import Rotation
 from typing_extensions import Self
 
 from TPTBox.core.np_utils import np_count_nonzero
@@ -113,6 +114,65 @@ def _extract_affine(self: Has_Grid, rm_key=(), **args):
             out.pop(k)
         return out
 
+    def change_affine(
+        self,
+        translation=None,
+        rotation_degrees=None,
+        scaling=None,
+        degrees=True,
+        inplace=False,
+    ):
+        """
+        Apply a transformation (translation, rotation, scaling) to the affine matrix.
+
+        Parameters:
+            translation: (n,) array-like in mm
+            rotation_degrees: (n,) array-like (pitch, yaw, roll) in degrees
+            scaling: (n,) array-like scaling factors along x, y, z
+        """
+        warnings.warn("change_affine is untested", stacklevel=2)
+        n = self.affine.shape[0]
+        transform = np.eye(n)
+
+        # Scaling
+        if scaling is not None:
+            assert len(scaling) == n - 1, f"Scaling must be a {n - 1}-element array-like."
+            S = np.diag([*list(scaling), 1])
+            transform = S @ transform
+
+        # Rotation
+        if rotation_degrees is not None:
+            assert len(rotation_degrees) == n - 1, f"Rotation must be a {n - 1}-element array-like."
+            rot = Rotation.from_euler("xyz", rotation_degrees, degrees=degrees).as_matrix()
+            R_mat = np.eye(n)
+            R_mat[: n - 1, : n - 1] = rot
+            transform = R_mat @ transform
+
+        # Translation
+        if translation is not None:
+            T = np.eye(n)
+            T[: n - 1, n - 1] = translation
+            transform = T @ transform
+        if not inplace:
+            self = self.copy()  # noqa: PLW0642
+        # Update the affine
+        self.affine = transform @ self.affine
+        return self
+
+    def change_affine_(self, translation=None, rotation_degrees=None, scaling=None, degrees=True):
+        return self.change_affine(
+            translation=translation,
+            rotation_degrees=rotation_degrees,
+            scaling=scaling,
+            degrees=degrees,
+            inplace=True,
+        )
+
+    def copy(self) -> Self:
+        raise NotImplementedError(
+            "The copy method must be implemented in the subclass. It should return a new instance of the same type with the same attributes."
+        )
+
     def assert_affine(
         self,
         other: Self | Has_Grid | None = None,
@@ -181,11 +241,11 @@ def assert_affine(
                 found_errors.append(f"rotation mismatch {self.rotation}, {rotation}") if not rotation_match else None
         if zoom is not None and (not ignore_missing_values or self.zoom is not None):
             if self.zoom is None:
-                found_errors.append(f"zoom mismatch {self.zoom}, {zoom}")
+                found_errors.append(f"spacing mismatch {self.zoom}, {zoom}")
             else:
                 zms_diff = (self.zoom[i] - zoom[i] for i in range(3))
                 zms_match = np.all([abs(a) <= error_tolerance for a in zms_diff])
-                found_errors.append(f"zoom mismatch {self.zoom}, {zoom}") if not zms_match else None
+                found_errors.append(f"spacing mismatch {self.zoom}, {zoom}") if not zms_match else None
         if orientation is not None and (not ignore_missing_values or self.affine is not None):
             if self.orientation is None:
                 found_errors.append(f"orientation mismatch {self.orientation}, {orientation}")
@@ -256,7 +316,14 @@ def get_empty_POI(self, points: dict | None = None):
         from TPTBox import POI
 
         p = {} if points is None else points
-        return POI(p, orientation=self.orientation, zoom=self.zoom, shape=self.shape, rotation=self.rotation, origin=self.origin)
+        return POI(
+            p,
+            orientation=self.orientation,
+            zoom=self.zoom,
+            shape=self.shape,
+            rotation=self.rotation,
+            origin=self.origin,
+        )
 
     def make_empty_POI(self, points: dict | None = None):
         from TPTBox import POI
@@ -267,7 +334,15 @@ def make_empty_POI(self, points: dict | None = None):
             args["level_one_info"] = self.level_one_info
             args["level_two_info"] = self.level_two_info
 
-        return POI(p, orientation=self.orientation, zoom=self.zoom, shape=self.shape, rotation=self.rotation, origin=self.origin, **args)
+        return POI(
+            p,
+            orientation=self.orientation,
+            zoom=self.zoom,
+            shape=self.shape,
+            rotation=self.rotation,
+            origin=self.origin,
+            **args,
+        )
 
     def make_empty_nii(self, seg=False, _arr=None):
         from TPTBox import NII
@@ -332,6 +407,30 @@ def to_deepali_grid(self, align_corners: bool = True):
         grid = grid.align_corners_(align_corners)
         return grid
 
+    @classmethod
+    def from_deepali_grid(cls, grid):
+        try:
+            from deepali.core import Grid as dp_Grid
+        except Exception:
+            log.print_error()
+            log.on_fail("run 'pip install hf-deepali' to install deepali")
+            raise
+        grid_: dp_Grid = grid
+        size = grid_.size()
+        spacing = grid_.spacing().cpu().numpy()
+        origin = grid_.origin().cpu().numpy()
+        direction = grid_.direction().cpu().numpy()
+        # Convert to ITK LPS convention
+        origin[:2] *= -1
+        direction[:2] *= -1
+        # Replace small values and -0 by 0
+        epsilon = sys.float_info.epsilon
+        origin[np.abs(origin) < epsilon] = 0
+        direction[np.abs(direction) < epsilon] = 0
+        grid = Grid(shape=size, origin=origin, spacing=spacing, rotation=direction)  # type: ignore
+
+        return grid
+
     def get_num_dims(self):
         return len(self.shape)
 
@@ -342,9 +441,14 @@ def __init__(self, **qargs) -> None:
         for k, v in qargs.items():
             if k == "spacing":
                 k = "zoom"  # noqa: PLW2901
+            if k == "direction":
+                k = "rotation"  # noqa: PLW2901
             if k == "rotation":
                 v = np.array(v)  # noqa: PLW2901
                 if len(v.shape) == 1:
                     s = int(np.sqrt(v.shape[0]))
                     v = v.reshape(s, s)  # noqa: PLW2901
             setattr(self, k, v)
+
+        ort = nio.io_orientation(self.affine)
+        self.orientation = nio.ornt2axcodes(ort)  # type: ignore