Bugfinder

TPTBox/__init__.py

@@ -27,6 +27,7 @@
 from TPTBox.core.poi import calc_poi_from_two_segs
 from TPTBox.core.poi import calc_poi_from_two_segs as calc_poi_labeled_buffered
 from TPTBox.core.poi_fun.poi_global import POI_Global
+from TPTBox.core.poi_fun.save_load import load_poi
 from TPTBox.core.vert_constants import ZOOMS, Location, Vertebra_Instance, v_idx2name, v_idx_order, v_name2idx
 
 # Logger

TPTBox/core/bids_files.py

@@ -1147,10 +1147,9 @@ def get_changed_path(  # noqa: C901
                 return out_path
             if not out_path.exists():
                 return out_path
-            if "run" in info:
-                info["run"] += 1
-            else:
-                info["run"] = 2
+            # "run" is a decimal entity and must stay a string: validate_entities() calls
+            # .isdecimal() on it, which an int does not support.
+            info["run"] = str(int(info["run"]) + 1) if "run" in info else "2"
 
     def save_changed_path(
         self,
@@ -1337,7 +1336,7 @@ def open_json(self) -> dict:
         Raises:
             KeyError: If no JSON file is registered in :attr:`file`.
         """
-        with open(self.file["json"]) as f:
+        with open(self.file["json"], encoding="utf-8") as f:
             return json.load(f)
 
     def open_poi(self, nii: TPTBox.Image_Reference | None = None) -> TPTBox.POI:

TPTBox/core/dicom/dicom2nii_utils.py

@@ -284,7 +284,7 @@ def test_name_conflict(json_ob: dict, file: str | Path) -> bool:
         ``False`` otherwise (file does not exist or content matches).
     """
     if Path(file).exists():
-        with open(file) as f:
+        with open(file, encoding="utf-8") as f:
             js = json.load(f)
             if "grid" in js:
                 del js["grid"]
@@ -328,7 +328,7 @@ def load_json(file: str | Path) -> dict:
     Returns:
         Parsed contents of the JSON file.
     """
-    with open(file) as file_handel:
+    with open(file, encoding="utf-8") as file_handel:
         return json.load(file_handel)
 
 

TPTBox/core/nii_wrapper.py

@@ -72,7 +72,6 @@
 if TYPE_CHECKING:
     from stl.mesh import Mesh
     from torch import device
-MODES = Literal["constant", "nearest", "reflect", "wrap"]
 _unpacked_nii = tuple[np.ndarray, AFFINE, nib.nifti1.Nifti1Header]
 _formatwarning = warnings.formatwarning
 
@@ -99,6 +98,33 @@ def formatwarning_tb(*args, **kwargs) -> str:
 _dtype_non_u = {"int8", "int16"}
 
 
+def _smallest_int_dtype(arr: np.ndarray, unsigned: bool) -> type:
+    """Smallest integer dtype that can represent every value in ``arr``.
+
+    Both bounds are considered: picking on ``max()`` alone silently wraps negative
+    values, because the casts here use ``casting="unsafe"``.
+
+    Args:
+        arr: Array whose value range determines the dtype.
+        unsigned: If True, choose an unsigned type (requires ``arr.min() >= 0``).
+
+    Returns:
+        The selected numpy dtype.
+    """
+    mi = arr.min()
+    ma = arr.max()
+    if unsigned:
+        assert mi >= 0, f"an unsigned dtype requires non-negative values, but the minimum is {mi}"
+        for cand, limit in ((np.uint8, 256), (np.uint16, 65536), (np.uint32, 2**32)):
+            if ma < limit:
+                return cand
+        return np.uint64
+    for cand, limit in ((np.int8, 128), (np.int16, 32768), (np.int32, 2**31)):
+        if ma < limit and mi >= -limit:
+            return cand
+    return np.int64
+
+
 def _check_if_nifty_is_lying_about_its_dtype(self: NII):
     """Infers the correct dtype by inspecting the actual value range of the NIfTI dataobj."""
     change_dtype = False
@@ -240,7 +266,7 @@ def __init__(self, nii: Nifti1Image|_unpacked_nii, seg=False,c_val=None, desc:st
         self.set_description(desc)
         if seg:
             self._unpack()
-            if isinstance(self.dtype,np.floating):
+            if np.issubdtype(self.dtype,np.floating):
                 self.set_dtype_("smallest_uint")
 
 
@@ -672,7 +698,7 @@ def set_array(self, arr: np.ndarray | Self, inplace=False, verbose: logging = Fa
             arr = arr.astype(np.uint8)
         if arr.dtype == np.float16:
             arr = arr.astype(np.float32)
-        if self.seg and isinstance(arr, (np.floating, float)):
+        if self.seg and np.issubdtype(arr.dtype, np.floating):
             arr = arr.astype(np.int32)
         #if self.dtype == arr.dtype: #type: ignore
         nii:_unpacked_nii = (arr,self.affine,self.header.copy())
@@ -712,30 +738,20 @@ def set_dtype(self, dtype: type | Literal['smallest_int', 'smallest_uint'] = np.
             The NII with the new dtype (``self`` when ``inplace=True``, a new NII otherwise).
         """
         sel = self if inplace else self.copy()
-        if dtype == "smallest_uint":
+        arr = None  # get_array() copies the whole volume; fetch it at most once
+        if dtype in ("smallest_uint", "smallest_int"):
             arr = self.get_array()
-            if arr.max()<256:
-                dtype = np.uint8
-            elif arr.max()<65536:
-                dtype = np.uint16
-            else:
-                dtype = np.int32
-        elif dtype == "smallest_int":
-            arr = self.get_array()
-            if arr.max()<128:
-                dtype = np.int8
-            elif arr.max()<32768:
-                dtype = np.int16
-            else:
-                dtype = np.int32
+            dtype = _smallest_int_dtype(arr, unsigned=dtype == "smallest_uint")
         if self.__unpacked:
             self._unpack()
             sel._arr = sel._arr.astype(dtype)
             sel.header.set_data_dtype(dtype)
         else:
             sel.nii.set_data_dtype(dtype)
             if sel.nii.get_data_dtype() != self.dtype: #type: ignore
-                sel.nii = Nifti1Image(self.get_array().astype(dtype,casting=casting,order=order),self.affine,self.header)
+                if arr is None:
+                    arr = self.get_array()
+                sel.nii = Nifti1Image(arr.astype(dtype,casting=casting,order=order),self.affine,self.header)
 
         return sel
     def set_dtype_(self, dtype: type | Literal['smallest_uint', 'smallest_int'] = np.float32, order: Literal["C", "F", "A", "K"] = 'K', casting: Literal["no", "equiv", "safe", "same_kind", "unsafe"] = "unsafe") -> Self:
@@ -808,7 +824,7 @@ def reorient(self:Self, axcodes_to: AX_CODES|str|None = ("P", "I", "R"), verbose
             new_img = arr, new_aff,self.header
             log.print("Image reoriented from", nio.ornt2axcodes(ornt_fr), "to", axcodes_to,verbose=verbose)
         else:
-            return self if not inplace else self.copy()
+            return self if inplace else self.copy()
         if inplace:
             self.nii = new_img
             return self
@@ -1118,13 +1134,13 @@ def rescale(self, voxel_spacing:float|tuple[float,...]=(1, 1, 1), c_val:float|No
             NII: A new NII object with the resampled image data.
         """
         if isinstance(voxel_spacing, (int,float)):
-            voxel_spacing =(voxel_spacing for _ in range(min(3,self.affine.shape[0]-1)))
+            voxel_spacing =tuple(voxel_spacing for _ in range(min(3,self.affine.shape[0]-1)))
         n = self.dims
         while  n> len(voxel_spacing):
             voxel_spacing = (*voxel_spacing, -1)
         if all(a in (-1, b) for a,b in zip(voxel_spacing, self.zoom)):
             log.print(f"Image already resampled to voxel size {self.zoom}",verbose=verbose)
-            return self.copy() if inplace else self
+            return self if inplace else self.copy()
 
         c_val = self.get_c_val(c_val)
         # resample to new voxel spacing based on the current x-y-z-orientation
@@ -1137,7 +1153,7 @@ def rescale(self, voxel_spacing:float|tuple[float,...]=(1, 1, 1), c_val:float|No
         voxel_spacing = tuple([v if v != -1 else z for v,z in zip_strict(voxel_spacing,zms)])
         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
+            return self if inplace else self.copy()
 
         # Calculate new shape
         new_shp = tuple(np.rint([shp[i] * zms[i] / voxel_spacing[i] for i in range(len(voxel_spacing))]).astype(int))
@@ -2059,6 +2075,11 @@ def is_segmentation_in_border(self,minimum=0,voxel_tolerance: int = 2,use_mm: bo
         - bool: True if the segmentation is within the defined tolerance of the
         border, False otherwise.
         """
+        # compute_crop(raise_error=False) returns full-extent slices for an empty mask rather
+        # than None, so an explicit emptiness check is needed - otherwise "nothing segmented"
+        # is reported as "touching the border".
+        if self.is_empty:
+            return False
         slices = self.compute_crop(minimum,dist=0,use_mm=use_mm,raise_error=False)
         if slices is None:
             return False
@@ -2361,9 +2382,10 @@ def save(self, file: str | Path, make_parents=True, verbose: logging = True, dty
             return self.save_nrrd(file,verbose=verbose)
 
         arr = self.get_array() if not self.seg else self.get_seg_array()
-        if isinstance(arr,np.floating) and self.seg:
-            self.set_dtype_("smallest_uint")
-            arr = self.get_array() if not self.seg else self.get_seg_array()
+        if self.seg and np.issubdtype(arr.dtype, np.floating):
+            # A segmentation must never be written out as float. Cast the local array:
+            # `save` is a query and must not mutate `self`.
+            arr = arr.astype(_smallest_int_dtype(arr, unsigned=True))
 
         self.header.set_data_dtype(arr.dtype)
         out = Nifti1Image(arr, self.affine,self.header)#,dtype=arr.dtype)
@@ -2601,7 +2623,8 @@ def __getitem__(self, key)-> Any:
         elif isinstance(key,np.ndarray):
             return self.get_array()[key]
         elif isinstance(key,slice):
-            self.__getitem__((key,Ellipsis,Ellipsis))
+            # pad with full slices for the trailing dimensions; Ellipsis is rejected above
+            return self.__getitem__((key, *(slice(None) for _ in range(len(self.shape) - 1))))
         else:
             raise TypeError("Invalid argument type:", type(key))
     def __setitem__(self, key,value):
@@ -2697,19 +2720,19 @@ def extract_label(self,label:int|Enum|Sequence[int]|Sequence[Enum]|None, keep_la
         assert self.seg, "extracting a label only makes sense for a segmentation mask"
         if label is None:
             if keep_label:
-                return self.copy() if inplace else self
+                return self if inplace else self.copy()
             else:
                 return self.clamp(0,1,inplace=inplace)
         seg_arr = self.get_seg_array()
 
+        if isinstance(label,str):
+            label = int(label)  # a str is also a Sequence, so this must come first
         if isinstance(label, Sequence):
             labels:int|list[int] = [idx.value if isinstance(idx,Enum) else idx for idx in label]
             assert 0 not in labels, 'Zero label does not make sense. This is the background'
         else:
             if isinstance(label,Enum):
                 label = label.value
-            if isinstance(label,str):
-                label = int(label)
 
             assert label != 0, 'Zero label does not make sense. This is the background'
             labels = label
@@ -2740,6 +2763,8 @@ def extract_label_(self, label: int | Enum | Sequence[int] | Sequence[Enum], kee
     def remove_labels(self,label:int|Enum|Sequence[int]|Sequence[Enum], inplace=False, verbose:logging=True, removed_to_label=0) -> Self:
         """If this NII is a segmentation you can single out one label."""
         assert label != 0, 'Zero label does not make sens.  This is the background'
+        if isinstance(label,str):
+            label = int(label)  # a str is also a Sequence, so this must come first
         if not isinstance(label,Sequence):
             label = [label] # type: ignore
         flat: list[int] = []

TPTBox/core/nii_wrapper_math.py

@@ -5,8 +5,6 @@
 from typing import TYPE_CHECKING, Union
 
 import numpy as np
-from skimage.metrics import peak_signal_noise_ratio as psnr
-from skimage.metrics import structural_similarity as ssim
 from typing_extensions import Self
 
 from TPTBox.core.np_utils import np_dice
@@ -286,12 +284,19 @@ def normalize(self,min_out = 0, max_out = 1, quantile = 1., clamp_lower:float|No
         """
         arr = self.get_array()
         max_v = np.quantile(arr[arr>0],q=quantile)
-        arr = self.clamp(clamp_lower,max_v,inplace=inplace)
-        arr -= arr.min() - min_out/max_out
-        arr /= arr.max() *max_out
-        assert arr.max() == max_out, f"{arr.max()} == {max_out}"
-        assert arr.min() == min_out
-        return self.set_array(arr.get_array(),inplace)
+        clamped = self.clamp(clamp_lower,max_v,inplace=inplace).get_array()
+        if not np.issubdtype(clamped.dtype, np.floating):
+            clamped = clamped.astype(np.float32)
+        mi = clamped.min()
+        ma = clamped.max()
+        if ma == mi:
+            # Degenerate (constant) image: nothing to spread out, map onto the lower bound.
+            clamped = np.full_like(clamped, min_out)
+        else:
+            clamped = (clamped - mi) / (ma - mi) * (max_out - min_out) + min_out
+            assert np.isclose(clamped.max(), max_out), f"{clamped.max()} == {max_out}"
+            assert np.isclose(clamped.min(), min_out), f"{clamped.min()} == {min_out}"
+        return self.set_array(clamped,inplace)
     def normalize_(self,min_out = 0, max_out = 1, quantile = 1., clamp_lower:float|None=None)->Self:
         """In-place variant of `normalize`."""
         return self.normalize(min_out = min_out, max_out = max_out, quantile = quantile, clamp_lower=clamp_lower,inplace=True)
@@ -413,10 +418,7 @@ def threshold(self,threshold=0.5, inplace=False)->Self:
         Returns:
             Self: Binarised segmentation instance.
         """
-        arr = self.get_array()
-        arr2 = arr.copy()
-        arr[arr2>=threshold] = 1
-        arr[arr2<=threshold] = 0
+        arr = (self.get_array() >= threshold).astype(np.uint8)
         nii = self if inplace else self.copy()
         nii.seg = True
         nii:NII = nii.set_array(arr,inplace,verbose=False)
@@ -450,9 +452,12 @@ def ssim(self, nii:NII_Proxy, min_v = 0)->float:
         Returns:
             float: SSIM score in the range [-1, 1] (1 = identical).
         """
+        # imported here: skimage.metrics pulls in scipy.stats, ~35% of `import TPTBox`
+        from skimage.metrics import structural_similarity as ssim
+
         img_1 = nii.get_array() - min_v
         img_2 = self.get_array() - min_v
-        img_1/= img_1.max()
+        img_1 = img_1/ img_1.max()  # out-of-place: /= fails on integer arrays
         img_1[img_1<=0] = 0
         img_2= img_2/ img_2.max()
         img_2[img_2<=0] = 0
@@ -474,9 +479,12 @@ def psnr(self,nii: NII_Proxy,min_v=0)->float:
         Returns:
             float: PSNR score in dB (higher is better; inf when images are identical).
         """
+        # imported here: skimage.metrics pulls in scipy.stats, ~35% of `import TPTBox`
+        from skimage.metrics import peak_signal_noise_ratio as psnr
+
         img_1 = nii.get_array() - min_v
         img_2 = self.get_array() - min_v
-        img_1/= img_1.max()
+        img_1 = img_1/ img_1.max()  # out-of-place: /= fails on integer arrays
         img_1[img_1<=0] = 0
         img_2= img_2/img_2.max()
         img_2[img_2<=0] = 0