Merge branch 'development_robert' of github.com:Hendrik-code/TPTBox into development_robert

Author: robert-graf

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/nii_wrapper.py

@@ -880,7 +880,7 @@ def reorient_same_as(self, img_as: Nifti1Image | Self, verbose:logging=False, in
         return self.reorient(axcodes_to=axcodes_to, verbose=verbose, inplace=inplace)
     def reorient_same_as_(self, img_as: Nifti1Image | Self, verbose:logging=False) -> Self:
         return self.reorient_same_as(img_as=img_as,verbose=verbose,inplace=True)
-    def rescale(self, voxel_spacing:float|tuple[float,...]=(1, 1, 1), c_val:float|None=None, verbose:logging=False, inplace=False,mode:MODES='nearest',order: int |None = None,align_corners:bool=False):
+    def rescale(self, voxel_spacing:float|tuple[float,...]=(1, 1, 1), c_val:float|None=None, verbose:logging=False, inplace=False,mode:MODES='nearest',order: int |None = None,align_corners:bool=False,atol=0.001):
         """
         Rescales the NIfTI image to a new voxel spacing.
 
@@ -895,6 +895,7 @@ def rescale(self, voxel_spacing:float|tuple[float,...]=(1, 1, 1), c_val:float|No
             mode (str, optional): One of the supported modes by scipy.ndimage.interpolation (e.g., "constant", "nearest",
                 "reflect", "wrap"). See the documentation for more details. Defaults to "constant".
             align_corners (bool|default): If True or not set and seg==True. Aline corners for scaling. This prevents segmentation mask to shift in a direction.
+            atol: absolute tolerance for skipping if already close in voxel_spacing
         Returns:
             NII: A new NII object with the resampled image data.
         """
@@ -916,7 +917,7 @@ def rescale(self, voxel_spacing:float|tuple[float,...]=(1, 1, 1), c_val:float|No
             order = 0 if self.seg else 3
         #print(aff.shape,shp,zms,voxel_spacing)
         voxel_spacing = tuple([v if v != -1 else z for v,z in zip_strict(voxel_spacing,zms)])
-        if voxel_spacing == self.zoom:
+        if np.isclose(voxel_spacing, self.zoom,atol=atol).all():
             log.print(f"Image already resampled to voxel size {self.zoom}",verbose=verbose)
             return self.copy() if inplace else self
 

TPTBox/images/poi_preview.png

@@ -121,9 +121,9 @@ def make_snapshot3D(
 
 
 def make_snapshot3D_parallel(
-    imgs: list[Path],
+    imgs: list[Path | str],
     output_paths: list[Image_Reference],
-    orientation: VIEW | list[VIEW] = "A",
+    view: VIEW | list[VIEW] = "A",
     ids_list: list[Sequence[int]] | None = None,
     smoothing=20,
     resolution: float = 2,
@@ -141,7 +141,7 @@ def make_snapshot3D_parallel(
                 kwds={
                     "output_path": out_path,
                     "img": img,
-                    "view": orientation,
+                    "view": view,
                     "ids_list": ids_list,
                     "smoothing": smoothing,
                     "resolution": resolution,

TPTBox/images/snp2D_example.png

@@ -79,14 +79,18 @@ def run_inference_on_file(
         nnunet_path = next(next(iter(model_path.glob(f"*{idx:03}*"))).glob("*__nnUNet*ResEnc*"))
     except StopIteration:
         nnunet_path = next(next(iter(model_path.glob(f"*{idx:03}*"))).glob("*__nnUNetPlans*"))
-    folds = [int(f.name.split("fold_")[-1]) for f in nnunet_path.glob("fold*")]
+    folds = sorted([int(f.name.split("fold_")[-1]) for f in nnunet_path.glob("fold*")])
     if max_folds is not None:
         folds = max_folds if isinstance(max_folds, list) else folds[:max_folds]
 
     # if idx in _unets:
     #    nnunet = _unets[idx]
     # else:
     print("load model", nnunet_path.name, "; folds", folds) if verbose else None
+    with open(Path(nnunet_path, "plans.json")) as f:
+        plans_info = json.load(f)
+    with open(Path(nnunet_path, "dataset.json")) as f:
+        ds_info = json.load(f)
     nnunet = load_inf_model(
         nnunet_path,
         allow_non_final=True,
@@ -101,10 +105,6 @@ def run_inference_on_file(
     )
 
     #    _unets[idx] = nnunet
-    with open(Path(nnunet_path, "plans.json")) as f:
-        plans_info = json.load(f)
-    with open(Path(nnunet_path, "dataset.json")) as f:
-        ds_info = json.load(f)
     if "orientation" in ds_info:
         orientation = ds_info["orientation"]
     zoom = None