New features robert (#70)

* add remove_to_label_for_cc_filter

* add remove_to

* add itk global_cord support

* add referecne for loading POI (POI.load)

* enable loading via Path

* add set_above_3_point_plane

* make an internal function to make vert-seg

* do not corp when empty instead of failing

* fix issue with empty images and cropout

* bug fixes

* ruff

* remove deprecated function

* x

---------

Co-authored-by: ga84mun <[email protected]>

Author: robert-graf

TPTBox/core/nii_wrapper.py

@@ -747,18 +747,18 @@ def pad_to(self,target_shape:list[int]|tuple[int,int,int] | Self, mode:MODES="co
             s = s.apply_crop(tuple(crop),inplace=inplace)
         return s.apply_pad(padding,inplace=inplace,mode=mode)
 
-    def apply_pad(self,padd:Sequence[tuple[int|None,int]],mode:MODES="constant",inplace = False):
+    def apply_pad(self,padd:Sequence[tuple[int|None,int]],mode:MODES="constant",inplace = False,verbose:logging=True):
         #TODO add other modes
         #TODO add testcases and options for modes
         transform = np.eye(4, dtype=int)
         for i, (before,_) in enumerate(padd):
             #transform[i, i] = pad_slice.step if pad_slice.step is not None else 1
             transform[i, 3] = -before  if before is not None else 0
         affine = self.affine.dot(transform)
-        print(mode)
         args = {}
         if mode == "constant":
             args["constant_values"]=self.get_c_val()
+        log.print(f"Padd {padd}; {mode=}, {args}",verbose=verbose)
         arr = np.pad(self.get_array(),padd,mode=mode,**args) # type: ignore
 
         nii:_unpacked_nii = (arr,affine,self.header)
@@ -1260,21 +1260,24 @@ def get_segmentation_connected_components(self, labels: int |list[int], connecti
             cc = {i: self.set_array(k) for i,k in cc.items()}
         return cc, cc_n
 
-    def get_connected_components(self, labels: int |list[int]=1, connectivity: int = 3, verbose: bool=False,inplace=False) -> Self:
-        arr = self.get_seg_array()
-        cc, _ = np_connected_components(arr, connectivity=connectivity, label_ref=labels, verbose=verbose)
-        out = None
-
-        for i,k in cc.items():
-            if out is None:
-                out = k
-            else:
-                out += i*k
-        if out is None:
-            return self if inplace else self.copy()
+    def get_connected_components(self, labels: int |list[int]=1, connectivity: int = 3, verbose: bool=False,inplace=False) -> Self:  # noqa: ARG002
+        out = np_get_largest_k_connected_components(self.get_seg_array(), label_ref=labels, connectivity=connectivity, return_original_labels=False)
         return self.set_array(out,inplace=inplace)
 
-    def filter_connected_components(self, labels: int |list[int]|None,min_volume:int=0,max_volume:int|None=None, max_count_component = None, connectivity: int = 3,keep_label=False, inplace=False):
+        #arr = self.get_seg_array()
+        #cc, _ = np_connected_components(arr, connectivity=connectivity, label_ref=labels, verbose=verbose)
+        #out = None
+
+        #for i,k in cc.items():
+        #    if out is None:
+        #        out = k
+        #    else:
+        #        out += i*k
+        #if out is None:
+        #    return self if inplace else self.copy()
+        #return self.set_array(out,inplace=inplace)
+
+    def filter_connected_components(self, labels: int |list[int]|None,min_volume:int=0,max_volume:int|None=None, max_count_component = None, connectivity: int = 3,removed_to_label=0,keep_label=False, inplace=False):
         """
         Filter connected components in a segmentation array based on specified volume constraints.
 
@@ -1289,7 +1292,7 @@ def filter_connected_components(self, labels: int |list[int]|None,min_volume:int
         Returns:
         None
         """
-        nii = self.get_largest_k_segmentation_connected_components(max_count_component,labels,connectivity=connectivity,return_original_labels=keep_label,min_volume=min_volume,max_volume=max_volume)
+        nii = self.get_largest_k_segmentation_connected_components(max_count_component,labels,connectivity=connectivity,return_original_labels=keep_label,min_volume=min_volume,max_volume=max_volume,removed_to_label=removed_to_label)
         if keep_label and labels is not None:
             if isinstance(labels,int):
                 labels = [labels]
@@ -1319,7 +1322,7 @@ def get_segmentation_connected_components_center_of_mass(self, label: int, conne
         return np_get_connected_components_center_of_mass(arr, label=label, connectivity=connectivity, sort_by_axis=sort_by_axis)
 
 
-    def get_largest_k_segmentation_connected_components(self, k: int | None, labels: int | list[int] | None = None, connectivity: int = 1, return_original_labels: bool = True,inplace=False,min_volume:int=0,max_volume:int|None=None):
+    def get_largest_k_segmentation_connected_components(self, k: int | None, labels: int | list[int] | None = None, connectivity: int = 1, return_original_labels: bool = True,inplace=False,min_volume:int=0,max_volume:int|None=None,removed_to_label=0):
         """Finds the largest k connected components in a given array (does NOT work with zero as label!)
 
         Args:
@@ -1329,7 +1332,7 @@ def get_largest_k_segmentation_connected_components(self, k: int | None, labels:
             return_original_labels (bool): If set to False, will label the components from 1 to k. Defaults to True
         """
         msk_i_data = self.get_seg_array()
-        out = np_get_largest_k_connected_components(msk_i_data, k=k, label_ref=labels, connectivity=connectivity, return_original_labels=return_original_labels,min_volume=min_volume,max_volume=max_volume)
+        out = np_get_largest_k_connected_components(msk_i_data, k=k, label_ref=labels, connectivity=connectivity, return_original_labels=return_original_labels,min_volume=min_volume,max_volume=max_volume,removed_to_label=removed_to_label)
         return self.set_array(out,inplace=inplace)
 
     def compute_surface_mask(self, connectivity: int, dilated_surface: bool = False):

TPTBox/core/np_utils.py

@@ -53,7 +53,7 @@ def np_extract_label(
         arr = arr.copy()
 
     if use_crop:
-        crop = np_bbox_binary(arr, px_dist=1)
+        crop = np_bbox_binary(arr, px_dist=1, raise_error=False)
         arrc = arr[crop]
     else:
         arrc = arr
@@ -82,10 +82,13 @@ def np_extract_label(
 
 def cc3dstatistics(arr: UINTARRAY, use_crop: bool = True) -> dict:
     assert np.issubdtype(arr.dtype, np.unsignedinteger), f"cc3dstatistics expects uint type, got {arr.dtype}"
-    if use_crop:
-        crop = np_bbox_binary(arr)
-        arrc = arr[crop]
-        return _cc3dstats(arrc)
+    try:
+        if use_crop:
+            crop = np_bbox_binary(arr, raise_error=False, px_dist=2)
+            arrc = arr[crop]
+            return _cc3dstats(arrc)
+    except ValueError as e:
+        print(e)
     return _cc3dstats(arr)
 
 
@@ -303,10 +306,7 @@ def np_dilate_msk(
         # try:
         arr_bin = arr.copy()
         arr_bin[np.isin(arr_bin, labels, invert=True)] = 0
-        try:
-            crop = np_bbox_binary(arr_bin, px_dist=1 + n_pixel)
-        except AssertionError:
-            crop = tuple([slice(None)] * arr.ndim)
+        crop = np_bbox_binary(arr_bin, px_dist=1 + n_pixel, raise_error=False)
         arrc = arr[crop]
     else:
         arrc = arr
@@ -379,10 +379,7 @@ def np_erode_msk(
         # try:
         arr_bin = arr.copy()
         arr_bin[np.isin(arr_bin, labels, invert=True)] = 0
-        try:
-            crop = np_bbox_binary(arr_bin, px_dist=1 + n_pixel)
-        except AssertionError:
-            crop = tuple([slice(None)] * arr.ndim)
+        crop = np_bbox_binary(arr_bin, px_dist=1 + n_pixel, raise_error=False)
         arrc = arr[crop]
     else:
         arrc = arr
@@ -461,9 +458,9 @@ def np_calc_crop_around_centerpoint(
     n_dim = len(poi)
     if isinstance(pad_to_size, int):
         pad_to_size = np.ones(n_dim) * pad_to_size
-    assert (
-        n_dim == len(arr.shape) == len(cutout_size) == len(pad_to_size)
-    ), f"dimension mismatch, got dim {n_dim}, poi {poi}, arr shape {arr.shape}, cutout {cutout_size}, pad_to_size {pad_to_size}"
+    assert n_dim == len(arr.shape) == len(cutout_size) == len(pad_to_size), (
+        f"dimension mismatch, got dim {n_dim}, poi {poi}, arr shape {arr.shape}, cutout {cutout_size}, pad_to_size {pad_to_size}"
+    )
 
     poi = tuple(int(i) for i in poi)
     shape = arr.shape
@@ -491,7 +488,7 @@ def np_calc_crop_around_centerpoint(
     )
 
 
-def np_bbox_binary(img: np.ndarray, px_dist: int | Sequence[int] | np.ndarray = 0) -> tuple[slice, ...]:
+def np_bbox_binary(img: np.ndarray, px_dist: int | Sequence[int] | np.ndarray = 0, raise_error=True) -> tuple[slice, ...]:
     """calculates a bounding box in n dimensions given a image (factor ~2 times faster than compute_crop)
 
     Args:
@@ -502,7 +499,11 @@ def np_bbox_binary(img: np.ndarray, px_dist: int | Sequence[int] | np.ndarray =
         list of boundary coordinates as slices tuple
     """
     assert img is not None, "bbox_nd: received None as image"
-    assert np_count_nonzero(img) > 0, "bbox_nd: img is empty, cannot calculate a bbox"
+    if np_count_nonzero(img) == 0:
+        if raise_error:
+            assert AssertionError("bbox_nd: img is empty, cannot calculate a bbox")
+        return tuple([slice(None)] * img.ndim)
+
     n = img.ndim
     shp = img.shape
     if isinstance(px_dist, int):
@@ -704,6 +705,8 @@ def np_get_largest_k_connected_components(
     return_original_labels: bool = True,
     min_volume: float = 0,
     max_volume: float | None = None,
+    removed_to_label=0,
+    _return_unsorted=False,
 ) -> UINTARRAY:
     """finds the largest k connected components in a given array (does NOT work with zero as label!)
 
@@ -731,6 +734,8 @@ def np_get_largest_k_connected_components(
     arr2[np.isin(arr, labels, invert=True)] = 0  # type:ignore
 
     labels_out, n = connected_components(arr2, connectivity=connectivity, return_N=True)
+    if _return_unsorted:
+        return labels_out
     if k is None:
         k = n
     k = min(k, n)  # if k > N, will return all N but still sorted
@@ -750,7 +755,8 @@ def np_get_largest_k_connected_components(
         if k == i:
             break
         i += 1
-
+    if removed_to_label != 0:
+        arr[np.logical_and(labels_out != 0, arr == 0)] = removed_to_label
     if return_original_labels:
         arr *= cc_out > 0  # to get original labels
         return arr
@@ -843,10 +849,7 @@ def np_translate_arr(arr: np.ndarray, translation_vector: tuple[int, int] | tupl
 
 
 def np_fill_holes(
-    arr: np.ndarray,
-    label_ref: LABEL_REFERENCE = None,
-    slice_wise_dim: int | None = None,
-    use_crop: bool = True,
+    arr: np.ndarray, label_ref: LABEL_REFERENCE = None, slice_wise_dim: int | None = None, use_crop: bool = True, pbar=False
 ) -> np.ndarray:
     """Fills holes in segmentations
 
@@ -863,16 +866,19 @@ def np_fill_holes(
     labels: Sequence[int] = _to_labels(arr, label_ref)
 
     if use_crop:
-        gcrop = np_bbox_binary(arr, px_dist=1)
+        gcrop = np_bbox_binary(arr, px_dist=1, raise_error=False)
         arrc = arr[gcrop]
     else:
         arrc = arr
+    if pbar:
+        from tqdm import tqdm
 
+        labels = tqdm(labels, desc="fill_holes")  # type: ignore
     for l in labels:  # type:ignore
         arr_l = arrc.copy()
         arr_l = np_extract_label(arr_l, l)
         if use_crop:
-            crop = np_bbox_binary(arr_l, px_dist=1)
+            crop = np_bbox_binary(arr_l, px_dist=1, raise_error=False)
             arr_lc = arr_l[crop]
         else:
             arr_lc = arr_l

TPTBox/core/poi.py

@@ -182,6 +182,10 @@ def zoom(self, value):
         else:
             self._zoom = tuple(round(float(v), ROUNDING_LVL) for v in value)  # type: ignore
 
+    @spacing.setter
+    def spacing(self, value):
+        self.zoom = value
+
     def clone(self, **qargs):
         return self.copy(**qargs)
 
@@ -560,12 +564,7 @@ def rescale(self, voxel_spacing: ZOOMS = (1, 1, 1), decimals=ROUNDING_LVL, verbo
         return self.copy(centroids=points, zoom=voxel_spacing, shape=shp)
 
     def rescale_(self, voxel_spacing: ZOOMS = (1, 1, 1), decimals=3, verbose: logging = False) -> Self:
-        return self.rescale(
-            voxel_spacing=voxel_spacing,
-            decimals=decimals,
-            verbose=verbose,
-            inplace=True,
-        )
+        return self.rescale(voxel_spacing=voxel_spacing, decimals=decimals, verbose=verbose, inplace=True)
 
     def to_global(self):
         """Converts the Centroids object to a global POI_Global object.
@@ -666,8 +665,8 @@ def make_point_cloud_nii(self, affine=None, s=8, sphere=False):
             affine = self.affine
         arr = np.zeros(self.shape_int)
         arr2 = np.zeros(self.shape_int)
-        s1 = min(s // 2, 1)
-        s2 = min(s - s1, 1)
+        s1 = max(s // 2, 1)
+        s2 = max(s - s1, 1)
         from math import ceil, floor
 
         if sphere:
@@ -723,7 +722,7 @@ def filter_points_inside_shape(self, inplace=False) -> Self:
         return self.copy(filtered_centroids)
 
     @classmethod
-    def load(cls, poi: POI_Reference):
+    def load(cls, poi: POI_Reference, reference: Has_Grid | None = None):
         """Load a Centroids object from various input sources.
 
         This method provides a convenient way to load a Centroids object from different sources,
@@ -761,7 +760,18 @@ def load(cls, poi: POI_Reference):
             >>> existing_poi = POI(...)
             >>> loaded_poi = POI.load(existing_poi)
         """
-        return load_poi(poi)
+        poi_obj = load_poi(poi)
+        if reference is not None:
+            if poi_obj.spacing is None:
+                poi_obj.spacing = reference.spacing
+            if poi_obj.rotation is None:
+                poi_obj.rotation = reference.rotation
+            if poi_obj.shape is None:
+                poi_obj.shape = reference.shape
+            if poi_obj.origin is None:
+                poi_obj.origin = reference.origin
+            reference.assert_affine(poi_obj)
+        return poi_obj
 
     def assert_affine(
         self,
@@ -1468,12 +1478,7 @@ def _is_not_yet_computed(ids_in_arr: Sequence[int], extend_to: POI | None, subre
 
 
 def calc_centroids(
-    msk: Image_Reference,
-    decimals=3,
-    first_stage=-1,
-    second_stage: int | Location = 50,
-    extend_to: POI | None = None,
-    inplace: bool = False,
+    msk: Image_Reference, decimals=3, first_stage=-1, second_stage: int | Location = 50, extend_to: POI | None = None, inplace: bool = False
 ) -> POI:
     """
     Calculates the centroid coordinates of each region in the given mask image.
@@ -1531,10 +1536,7 @@ def calc_centroids(
 ######## Utility #######
 
 
-def calc_poi_average(
-    pois: list[POI],
-    keep_points_not_present_in_all_pois: bool = False,
-) -> POI:
+def calc_poi_average(pois: list[POI], keep_points_not_present_in_all_pois: bool = False) -> POI:
     """Calculates average of POI across list of POIs and removes all points that are not fully present in all given POIs
 
     Args:

TPTBox/core/poi_fun/poi_global.py

@@ -17,7 +17,17 @@ class POI_Global(Abstract_POI):
     Inherits from the `Abstract_POI` class and contains methods for converting the POI to different coordinate systems.
     """
 
-    def __init__(self, input_poi: poi.POI | POI_Descriptor):
+    def __init__(self, input_poi: poi.POI | POI_Descriptor | dict[str, dict[str, tuple[float, ...]]], itk_coords: bool = False):
+        self.itk_coords = itk_coords
+
+        if isinstance(input_poi, dict):
+            global_points = POI_Descriptor()
+            self.info = {}
+            self.format = None
+            for k1, d1 in input_poi.items():
+                for k2, v in d1.items():
+                    global_points[k1:k2] = v
+
         if isinstance(input_poi, POI_Descriptor):
             global_points = input_poi
             self.info = {}
@@ -78,7 +88,7 @@ def to_other_poi(self, ref: poi.POI_Reference) -> poi.POI:
         """
         return self.to_other(poi.POI.load(ref))
 
-    def to_other(self, msk: Has_Grid) -> poi.POI:
+    def to_other(self, msk: Has_Grid, verbose=False) -> poi.POI:
         """
         Convert the POI to another coordinate system.
 
@@ -90,7 +100,12 @@ def to_other(self, msk: Has_Grid) -> poi.POI:
         """
         out = poi.POI_Descriptor(definition=self._get_centroids().definition)
         for k1, k2, v in self.items():
+            if self.itk_coords:
+                assert len(v) == 3, "n-d vec not implemented for n != 3"
+                v = (-v[0], -v[1], v[2])  # noqa: PLW2901
             v_out = msk.global_to_local(v)
+            if verbose:
+                print(v, "-->", v_out)
             out[k1, k2] = tuple(v_out)
 
         return poi.POI(centroids=out, **msk._extract_affine(), info=self.info, format=self.format)
@@ -101,4 +116,5 @@ def copy(self, centroids: POI_Descriptor | None = None) -> Self:
         p = POI_Global(centroids)
         p.format = self.format
         p.info = deepcopy(self.info)
+        p.itk_coords = self.itk_coords
         return p  # type: ignore