From 24842bfb8b0d019ed955f5749e3a6c873257dc2b Mon Sep 17 00:00:00 2001
From: Mathis Rasmussen <mathis@oncology.au.dk>
Date: Fri, 20 Jan 2023 11:45:41 +0100
Subject: [PATCH 01/10] added smoothing factor, which is just skipping every x
 coord. Also added adjustments for high res structure sets

---
 rt_utils/image_helper.py | 11 ++++++++---
 rt_utils/rtstruct.py     |  7 +++++++
 rt_utils/utils.py        |  6 +++++-
 3 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/rt_utils/image_helper.py b/rt_utils/image_helper.py
index b030987..9431ed0 100644
--- a/rt_utils/image_helper.py
+++ b/rt_utils/image_helper.py
@@ -60,7 +60,10 @@ def get_contours_coords(roi_data: ROIData, series_data):
             mask_slice = create_pin_hole_mask(mask_slice, roi_data.approximate_contours)
 
         # Get contours from mask
-        contours, _ = find_mask_contours(mask_slice, roi_data.approximate_contours)
+        contours, _ = find_mask_contours(mask_slice,
+                                         roi_data.approximate_contours,
+                                         smoothing_factor=roi_data.smoothing_factor,
+                                         scaling_factor=roi_data.scaling_factor)
         validate_contours(contours)
 
         # Format for DICOM
@@ -82,7 +85,7 @@ def get_contours_coords(roi_data: ROIData, series_data):
     return series_contours
 
 
-def find_mask_contours(mask: np.ndarray, approximate_contours: bool):
+def find_mask_contours(mask: np.ndarray, approximate_contours: bool, smoothing_factor=1, scaling_factor=1):
     approximation_method = (
         cv.CHAIN_APPROX_SIMPLE if approximate_contours else cv.CHAIN_APPROX_NONE
     )
@@ -93,8 +96,10 @@ def find_mask_contours(mask: np.ndarray, approximate_contours: bool):
     contours = list(
         contours
     )  # Open-CV updated contours to be a tuple so we convert it back into a list here
+
+    assert smoothing_factor > 0
     for i, contour in enumerate(contours):
-        contours[i] = [[pos[0][0], pos[0][1]] for pos in contour]
+        contours[i] = [[(contour[i][0][0]/scaling_factor), (contour[i][0][1]/scaling_factor)] for i in range(0, len(contour), smoothing_factor)]
     hierarchy = hierarchy[0]  # Format extra array out of data
 
     return contours, hierarchy
diff --git a/rt_utils/rtstruct.py b/rt_utils/rtstruct.py
index dfe82be..ee4583a 100644
--- a/rt_utils/rtstruct.py
+++ b/rt_utils/rtstruct.py
@@ -35,6 +35,7 @@ def add_roi(
         use_pin_hole: bool = False,
         approximate_contours: bool = True,
         roi_generation_algorithm: Union[str, int] = 0,
+        smoothing_factor: int = 1,
     ):
         """
         Add a ROI to the rtstruct given a 3D binary mask for the ROI's at each slice
@@ -43,6 +44,10 @@ def add_roi(
         If approximate_contours is set to False, no approximation will be done when generating contour data, leading to much larger amount of contour data
         """
 
+        ## If upscaled coords are given, they should be adjusted accordingly
+        rows = self.series_data[0][0x00280010].value
+        scaling_factor = int(mask.shape[0] / rows)
+
         # TODO test if name already exists
         self.validate_mask(mask)
         roi_number = len(self.ds.StructureSetROISequence) + 1
@@ -56,6 +61,8 @@ def add_roi(
             use_pin_hole,
             approximate_contours,
             roi_generation_algorithm,
+            smoothing_factor,
+            scaling_factor,
         )
 
         self.ds.ROIContourSequence.append(
diff --git a/rt_utils/utils.py b/rt_utils/utils.py
index f04089e..626bdd1 100644
--- a/rt_utils/utils.py
+++ b/rt_utils/utils.py
@@ -41,7 +41,6 @@ class SOPClassUID:
 @dataclass
 class ROIData:
     """Data class to easily pass ROI data to helper methods."""
-
     mask: str
     color: Union[str, List[int]]
     number: int
@@ -51,6 +50,11 @@ class ROIData:
     use_pin_hole: bool = False
     approximate_contours: bool = True
     roi_generation_algorithm: Union[str, int] = 0
+    smoothing_factor: int = 1
+    scaling_factor: int = 1
+
+
+
 
     def __post_init__(self):
         self.validate_color()

From 22fe96699b592647d62a7bfab00a1ad46e55bdaa Mon Sep 17 00:00:00 2001
From: Mathis Rasmussen <mathis@oncology.au.dk>
Date: Mon, 12 Jun 2023 21:24:09 +0200
Subject: [PATCH 02/10] Add smoothing pipeline

---
 rt_utils/image_helper.py | 26 ++++++++++--------
 rt_utils/rtstruct.py     | 20 +++++++++-----
 rt_utils/smoothing.py    | 58 ++++++++++++++++++++++++++++++++++++++++
 rt_utils/utils.py        |  7 +++--
 4 files changed, 90 insertions(+), 21 deletions(-)
 create mode 100644 rt_utils/smoothing.py

diff --git a/rt_utils/image_helper.py b/rt_utils/image_helper.py
index 9431ed0..636835c 100644
--- a/rt_utils/image_helper.py
+++ b/rt_utils/image_helper.py
@@ -1,7 +1,8 @@
 import os
-from typing import List
+from typing import List, Union
 from enum import IntEnum
 
+import SimpleITK
 import cv2 as cv
 import numpy as np
 from pydicom import dcmread
@@ -62,8 +63,9 @@ def get_contours_coords(roi_data: ROIData, series_data):
         # Get contours from mask
         contours, _ = find_mask_contours(mask_slice,
                                          roi_data.approximate_contours,
-                                         smoothing_factor=roi_data.smoothing_factor,
                                          scaling_factor=roi_data.scaling_factor)
+        if not contours:
+            continue
         validate_contours(contours)
 
         # Format for DICOM
@@ -85,7 +87,8 @@ def get_contours_coords(roi_data: ROIData, series_data):
     return series_contours
 
 
-def find_mask_contours(mask: np.ndarray, approximate_contours: bool, smoothing_factor=1, scaling_factor=1):
+
+def find_mask_contours(mask: np.ndarray, approximate_contours: bool, scaling_factor: int):
     approximation_method = (
         cv.CHAIN_APPROX_SIMPLE if approximate_contours else cv.CHAIN_APPROX_NONE
     )
@@ -97,10 +100,10 @@ def find_mask_contours(mask: np.ndarray, approximate_contours: bool, smoothing_f
         contours
     )  # Open-CV updated contours to be a tuple so we convert it back into a list here
 
-    assert smoothing_factor > 0
     for i, contour in enumerate(contours):
-        contours[i] = [[(contour[i][0][0]/scaling_factor), (contour[i][0][1]/scaling_factor)] for i in range(0, len(contour), smoothing_factor)]
-    hierarchy = hierarchy[0]  # Format extra array out of data
+        contours[i] = [[(contour[i][0][0] / scaling_factor), (contour[i][0][1] / scaling_factor)] for i in
+                       range(0, len(contour))]
+    # hierarchy = hierarchy[0]  # Format extra array out of data
 
     return contours, hierarchy
 
@@ -131,7 +134,7 @@ def create_pin_hole_mask(mask: np.ndarray, approximate_contours: bool):
 
 
 def draw_line_upwards_from_point(
-    mask: np.ndarray, start, fill_value: int
+        mask: np.ndarray, start, fill_value: int
 ) -> np.ndarray:
     line_width = 2
     end = (start[0], start[1] - 1)
@@ -201,7 +204,7 @@ def get_patient_to_pixel_transformation_matrix(series_data):
 
 
 def apply_transformation_to_3d_points(
-    points: np.ndarray, transformation_matrix: np.ndarray
+        points: np.ndarray, transformation_matrix: np.ndarray
 ):
     """
     * Augment each point with a '1' as the fourth coordinate to allow translation
@@ -224,7 +227,7 @@ def get_slice_directions(series_slice: Dataset):
     slice_direction = np.cross(row_direction, column_direction)
 
     if not np.allclose(
-        np.dot(row_direction, column_direction), 0.0, atol=1e-3
+            np.dot(row_direction, column_direction), 0.0, atol=1e-3
     ) or not np.allclose(np.linalg.norm(slice_direction), 1.0, atol=1e-3):
         raise Exception("Invalid Image Orientation (Patient) attribute")
 
@@ -268,7 +271,7 @@ def get_slice_contour_data(series_slice: Dataset, contour_sequence: Sequence):
 
 
 def get_slice_mask_from_slice_contour_data(
-    series_slice: Dataset, slice_contour_data, transformation_matrix: np.ndarray
+        series_slice: Dataset, slice_contour_data, transformation_matrix: np.ndarray
 ):
     # Go through all contours in a slice, create polygons in correct space and with a correct format 
     # and append to polygons array (appropriate for fillPoly) 
@@ -280,9 +283,10 @@ def get_slice_mask_from_slice_contour_data(
         polygon = np.array(polygon).squeeze()
         polygons.append(polygon)
     slice_mask = create_empty_slice_mask(series_slice).astype(np.uint8)
-    cv.fillPoly(img=slice_mask, pts = polygons, color = 1)
+    cv.fillPoly(img=slice_mask, pts=polygons, color=1)
     return slice_mask
 
+
 def create_empty_series_mask(series_data):
     ref_dicom_image = series_data[0]
     mask_dims = (
diff --git a/rt_utils/rtstruct.py b/rt_utils/rtstruct.py
index ee4583a..abef4d9 100644
--- a/rt_utils/rtstruct.py
+++ b/rt_utils/rtstruct.py
@@ -1,11 +1,11 @@
-from typing import List, Union
+from typing import List, Union, Dict
 
 import numpy as np
 from pydicom.dataset import FileDataset
 
 from rt_utils.utils import ROIData
-from . import ds_helper, image_helper
-
+from . import ds_helper, image_helper, smoothing
+from typing import Tuple
 
 class RTStruct:
     """
@@ -35,7 +35,9 @@ def add_roi(
         use_pin_hole: bool = False,
         approximate_contours: bool = True,
         roi_generation_algorithm: Union[str, int] = 0,
-        smoothing_factor: int = 1,
+        apply_smoothing: Union[str, None] = None, # strings can be "2d" or "3d"
+        smoothing_parameters: Dict = smoothing.default_smoothing_parameters
+
     ):
         """
         Add a ROI to the rtstruct given a 3D binary mask for the ROI's at each slice
@@ -43,6 +45,13 @@ def add_roi(
         If use_pin_hole is set to true, will cut a pinhole through ROI's with holes in them so that they are represented with one contour
         If approximate_contours is set to False, no approximation will be done when generating contour data, leading to much larger amount of contour data
         """
+        if apply_smoothing:
+            if apply_smoothing == "3d":
+                mask = smoothing.pipeline_3d(mask, **smoothing_parameters)
+            elif apply_smoothing == "2d":
+                mask = smoothing.pipeline_2d(mask, **smoothing_parameters)
+            else:
+                print("Invalid input. Use '2d' or '3d'. Otherwise leave as None")
 
         ## If upscaled coords are given, they should be adjusted accordingly
         rows = self.series_data[0][0x00280010].value
@@ -61,8 +70,7 @@ def add_roi(
             use_pin_hole,
             approximate_contours,
             roi_generation_algorithm,
-            smoothing_factor,
-            scaling_factor,
+            scaling_factor
         )
 
         self.ds.ROIContourSequence.append(
diff --git a/rt_utils/smoothing.py b/rt_utils/smoothing.py
new file mode 100644
index 0000000..6a20f6a
--- /dev/null
+++ b/rt_utils/smoothing.py
@@ -0,0 +1,58 @@
+import SimpleITK as sitk
+import numpy as np
+from scipy import ndimage, signal
+from typing import Union, Tuple
+default_smoothing_parameters = {
+    "scaling_factor": (3, 3, 1),
+    "sigma": 2,
+    "threshold": 0.4,
+    "kernel_size": (3, 3, 1),
+    "iterations": 2
+}
+def kron_upscale(mask: np.ndarray, scaling_factor: Tuple[int, ...]):
+    return np.kron(mask, np.ones(scaling_factor))
+
+
+def gaussian_blur(mask: np.ndarray, sigma: float):
+    return ndimage.gaussian_filter(mask, sigma=sigma)
+
+
+def binary_threshold(mask: np.ndarray, threshold: float):
+    return mask > threshold
+
+def median_filter(mask: np.ndarray, kernel_size: Union[int, Tuple[int, ...]]):
+    return ndimage.median_filter(mask.astype(float), size=kernel_size, mode="nearest")
+
+def pipeline_3d(mask: np.ndarray,
+                iterations: int,
+                scaling_factor: int,
+                sigma: float,
+                threshold: float,
+                kernel_size: Union[int, Tuple[int, ...]]):
+    scaling_factor = (scaling_factor, scaling_factor, 1)
+    for i in range(iterations):
+        mask = kron_upscale(mask=mask, scaling_factor=scaling_factor)
+        mask = gaussian_blur(mask=mask, sigma=sigma)
+        mask = binary_threshold(mask=mask, threshold=threshold)
+        mask = median_filter(mask=mask, kernel_size=kernel_size)
+        mask = mask.astype(bool)
+    return mask
+
+def pipeline_2d(mask: np.ndarray,
+                iterations: int,
+                scaling_factor: int,
+                sigma: float,
+                threshold: float,
+                kernel_size: Union[int, Tuple[int, ...]]):
+    scaling_factor = (scaling_factor, scaling_factor, 1)
+    for i in range(iterations):
+        mask = kron_upscale(mask=mask, scaling_factor=scaling_factor)
+        for z in range(mask.shape[2]):
+            slice = mask[:, : , z]
+            slice = gaussian_blur(mask=slice, sigma=sigma)
+            slice = binary_threshold(mask=slice, threshold=threshold)
+            slice = median_filter(mask=slice, kernel_size=kernel_size)
+            mask[:, :, z] = slice
+        mask = mask.astype(bool)
+    return mask
+
diff --git a/rt_utils/utils.py b/rt_utils/utils.py
index 626bdd1..3afb5c6 100644
--- a/rt_utils/utils.py
+++ b/rt_utils/utils.py
@@ -1,4 +1,4 @@
-from typing import List, Union
+from typing import List, Union, Tuple
 from random import randrange
 from pydicom.uid import PYDICOM_IMPLEMENTATION_UID
 from dataclasses import dataclass
@@ -50,10 +50,9 @@ class ROIData:
     use_pin_hole: bool = False
     approximate_contours: bool = True
     roi_generation_algorithm: Union[str, int] = 0
-    smoothing_factor: int = 1
     scaling_factor: int = 1
-
-
+    smooth_radius: Union[int, None] = None
+    smooth_scale: Union[int, None] = None
 
 
     def __post_init__(self):

From 5db67646e1d8a59e2e570debd00b9e26ac2a9a86 Mon Sep 17 00:00:00 2001
From: Mathis Rasmussen <mathis@oncology.au.dk>
Date: Wed, 14 Jun 2023 10:31:57 +0200
Subject: [PATCH 03/10] Add memory efficient smoothing

---
 rt_utils/rtstruct.py  |  11 +--
 rt_utils/smoothing.py | 158 +++++++++++++++++++++++++++++++-----------
 2 files changed, 119 insertions(+), 50 deletions(-)

diff --git a/rt_utils/rtstruct.py b/rt_utils/rtstruct.py
index abef4d9..6b0582a 100644
--- a/rt_utils/rtstruct.py
+++ b/rt_utils/rtstruct.py
@@ -36,8 +36,7 @@ def add_roi(
         approximate_contours: bool = True,
         roi_generation_algorithm: Union[str, int] = 0,
         apply_smoothing: Union[str, None] = None, # strings can be "2d" or "3d"
-        smoothing_parameters: Dict = smoothing.default_smoothing_parameters
-
+        smoothing_parameters: Union[Dict, None] = None
     ):
         """
         Add a ROI to the rtstruct given a 3D binary mask for the ROI's at each slice
@@ -46,12 +45,8 @@ def add_roi(
         If approximate_contours is set to False, no approximation will be done when generating contour data, leading to much larger amount of contour data
         """
         if apply_smoothing:
-            if apply_smoothing == "3d":
-                mask = smoothing.pipeline_3d(mask, **smoothing_parameters)
-            elif apply_smoothing == "2d":
-                mask = smoothing.pipeline_2d(mask, **smoothing_parameters)
-            else:
-                print("Invalid input. Use '2d' or '3d'. Otherwise leave as None")
+            mask = smoothing.pipeline(mask=mask, apply_smoothing=apply_smoothing,
+                                      smoothing_parameters=smoothing_parameters)
 
         ## If upscaled coords are given, they should be adjusted accordingly
         rows = self.series_data[0][0x00280010].value
diff --git a/rt_utils/smoothing.py b/rt_utils/smoothing.py
index 6a20f6a..fae1aac 100644
--- a/rt_utils/smoothing.py
+++ b/rt_utils/smoothing.py
@@ -1,58 +1,132 @@
 import SimpleITK as sitk
 import numpy as np
 from scipy import ndimage, signal
-from typing import Union, Tuple
+from typing import Union, Tuple, Dict
+import logging
+
 default_smoothing_parameters = {
-    "scaling_factor": (3, 3, 1),
-    "sigma": 2,
-    "threshold": 0.4,
-    "kernel_size": (3, 3, 1),
-    "iterations": 2
+    "iterations": 3,
+    "np_kron": {"scaling_factor": 2},
+    "ndimage_gaussian_filter": {"sigma": 2,
+                                "radius": 3},
+    "threshold": {"threshold": 0.4},
+    "ndimage_median_filter": {"size": 3,
+                              "mode": "nearest"}
 }
-def kron_upscale(mask: np.ndarray, scaling_factor: Tuple[int, ...]):
+
+def kron_upscale(mask: np.ndarray, **kwargs):
+    scaling_factor = (kwargs["scaling_factor"], kwargs["scaling_factor"], 1)
     return np.kron(mask, np.ones(scaling_factor))
 
+def gaussian_blur(mask: np.ndarray, **kwargs):
+    return ndimage.gaussian_filter(mask, **kwargs)
 
-def gaussian_blur(mask: np.ndarray, sigma: float):
-    return ndimage.gaussian_filter(mask, sigma=sigma)
+def binary_threshold(mask: np.ndarray, **kwargs):
+    return mask > kwargs["threshold"]
 
+def median_filter(mask: np.ndarray, **kwargs):
+    return ndimage.median_filter(mask.astype(float), **kwargs)
 
-def binary_threshold(mask: np.ndarray, threshold: float):
-    return mask > threshold
+def crop_mask(mask: np.ndarray):
+    three_d = (len(mask.shape) == 3)
+    if three_d:
+        x, y, z = np.nonzero(mask)
+        x_max = (x.max() + 1) if x.max() < mask.shape[0] else x.max()
+        y_max = (y.max() + 1) if y.max() < mask.shape[1] else y.max()
+        z_max = (z.max() + 1) if z.max() < mask.shape[2] else z.max()
+        bbox = np.array([x.min(), x_max, y.min(), y_max, z.min(), z_max])
 
-def median_filter(mask: np.ndarray, kernel_size: Union[int, Tuple[int, ...]]):
-    return ndimage.median_filter(mask.astype(float), size=kernel_size, mode="nearest")
+        return mask[bbox[0]: bbox[1],
+                    bbox[2]: bbox[3],
+                    bbox[4]: bbox[5]], bbox
+    else:
+        x, y = np.nonzero(mask)
+        x_max = (x.max() + 1) if x.max() < mask.shape[0] else x.max()
+        y_max = (y.max() + 1) if y.max() < mask.shape[1] else y.max()
+        bbox = [x.min(), x_max, y.min(), y_max]
 
-def pipeline_3d(mask: np.ndarray,
-                iterations: int,
-                scaling_factor: int,
-                sigma: float,
-                threshold: float,
-                kernel_size: Union[int, Tuple[int, ...]]):
-    scaling_factor = (scaling_factor, scaling_factor, 1)
-    for i in range(iterations):
-        mask = kron_upscale(mask=mask, scaling_factor=scaling_factor)
-        mask = gaussian_blur(mask=mask, sigma=sigma)
-        mask = binary_threshold(mask=mask, threshold=threshold)
-        mask = median_filter(mask=mask, kernel_size=kernel_size)
-        mask = mask.astype(bool)
-    return mask
+        return mask[bbox[0]: bbox[1],
+                    bbox[2]: bbox[3]], bbox
+
+def restore_mask_dimensions(cropped_mask: np.ndarray, new_shape, bbox):
+    new_mask = np.zeros(new_shape)
+
+    new_mask[bbox[0]: bbox[1], bbox[2]: bbox[3], bbox[4]: bbox[5]] = cropped_mask
+    return new_mask.astype(bool)
+
+def iteration_2d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, ndimage_median_filter):
+    cropped_mask = kron_upscale(mask=cropped_mask, **np_kron)
+
+    for z_idx in range(cropped_mask.shape[2]):
+        slice = cropped_mask[:, :, z_idx]
+        slice = gaussian_blur(mask=slice, **ndimage_gaussian_filter)
+        slice = binary_threshold(mask=slice, **threshold)
+        slice = median_filter(mask=slice, **ndimage_median_filter)
+
+        cropped_mask[:, :, z_idx] = slice
+    return cropped_mask
+
+def iteration_3d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, ndimage_median_filter):
+    cropped_mask = kron_upscale(mask=mask, **np_kron)
+    cropped_mask = gaussian_blur(mask=cropped_mask, **ndimage_gaussian_filter)
+    cropped_mask = binary_threshold(mask=cropped_mask, **threshold)
+    cropped_mask = median_filter(mask=cropped_mask, **ndimage_median_filter)
 
-def pipeline_2d(mask: np.ndarray,
-                iterations: int,
-                scaling_factor: int,
-                sigma: float,
-                threshold: float,
-                kernel_size: Union[int, Tuple[int, ...]]):
-    scaling_factor = (scaling_factor, scaling_factor, 1)
+    return cropped_mask
+
+def pipeline(mask: np.ndarray,
+             apply_smoothing: str,
+             smoothing_parameters: Union[Dict, None]):
+
+    if not smoothing_parameters:
+        smoothing_parameters = default_smoothing_parameters
+
+    iterations = smoothing_parameters["iterations"]
+    np_kron = smoothing_parameters["np_kron"]
+    ndimage_gaussian_filter = smoothing_parameters["ndimage_gaussian_filter"]
+    threshold = smoothing_parameters["threshold"]
+    ndimage_median_filter = smoothing_parameters["ndimage_median_filter"]
+
+    logging.info(f"Original mask shape {mask.shape}")
+    logging.info(f"Cropping mask to non-zero")
+    cropped_mask, bbox = crop_mask(mask)
+    final_shape, final_bbox = get_final_mask_shape_and_bbox(mask=mask,
+                                                            scaling_factor=np_kron["scaling_factor"],
+                                                            iterations=iterations,
+                                                            bbox=bbox)
+    logging.info(f"Final scaling with factor: {np_kron['scaling_factor']} in {iterations} iterations")
     for i in range(iterations):
-        mask = kron_upscale(mask=mask, scaling_factor=scaling_factor)
-        for z in range(mask.shape[2]):
-            slice = mask[:, : , z]
-            slice = gaussian_blur(mask=slice, sigma=sigma)
-            slice = binary_threshold(mask=slice, threshold=threshold)
-            slice = median_filter(mask=slice, kernel_size=kernel_size)
-            mask[:, :, z] = slice
-        mask = mask.astype(bool)
+        logging.info(f"Iteration {i} out of {iterations}")
+        logging.info(f"Applying filters")
+
+        if apply_smoothing == "2d":
+            mask = iteration_2d(mask,
+                                np_kron=np_kron,
+                                ndimage_gaussian_filter=ndimage_gaussian_filter,
+                                threshold=threshold,
+                                ndimage_median_filter=ndimage_median_filter)
+        elif apply_smoothing == "3d":
+            mask = iteration_3d(mask,
+                                np_kron=np_kron,
+                                ndimage_gaussian_filter=ndimage_gaussian_filter,
+                                threshold=threshold,
+                                ndimage_median_filter=ndimage_median_filter)
+        else:
+            raise Exception("Wrong dimension parameter. Use '2d' or '3d'.")
+
+    # Restore dimensions
+    logging.info("Restoring original mask shape")
+    mask = restore_mask_dimensions(cropped_mask, final_shape, final_bbox)
     return mask
 
+def get_final_mask_shape_and_bbox(mask, bbox, scaling_factor, iterations):
+    final_scaling_factor = pow(scaling_factor, iterations)
+
+    final_shape = np.array(mask.shape)
+    final_shape[:2] *= final_scaling_factor
+    bbox[:4] *= final_scaling_factor
+    logging.info("Final shape: ", final_shape)
+    logging.info("Final bbox: ", bbox)
+    return final_shape, bbox
+
+

From fa379b90871e9c9f72ea998649c212353af93b74 Mon Sep 17 00:00:00 2001
From: Mathis Rasmussen <mathis@oncology.au.dk>
Date: Wed, 14 Jun 2023 10:47:01 +0200
Subject: [PATCH 04/10] Small bug fixed in 2d pipeline

---
 rt_utils/smoothing.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/rt_utils/smoothing.py b/rt_utils/smoothing.py
index fae1aac..a22c503 100644
--- a/rt_utils/smoothing.py
+++ b/rt_utils/smoothing.py
@@ -55,7 +55,7 @@ def restore_mask_dimensions(cropped_mask: np.ndarray, new_shape, bbox):
     return new_mask.astype(bool)
 
 def iteration_2d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, ndimage_median_filter):
-    cropped_mask = kron_upscale(mask=cropped_mask, **np_kron)
+    cropped_mask = kron_upscale(mask=mask, **np_kron)
 
     for z_idx in range(cropped_mask.shape[2]):
         slice = cropped_mask[:, :, z_idx]

From 22f4ba2cde11f7f07a516a93a03f8be2b3fbaaca Mon Sep 17 00:00:00 2001
From: Mathis Rasmussen <mathis@oncology.au.dk>
Date: Wed, 14 Jun 2023 13:57:48 +0200
Subject: [PATCH 05/10] Fix

---
 rt_utils/image_helper.py |  3 +-
 rt_utils/smoothing.py    | 81 +++++++++++++++++++++-------------------
 2 files changed, 45 insertions(+), 39 deletions(-)

diff --git a/rt_utils/image_helper.py b/rt_utils/image_helper.py
index 636835c..63c5b43 100644
--- a/rt_utils/image_helper.py
+++ b/rt_utils/image_helper.py
@@ -103,7 +103,8 @@ def find_mask_contours(mask: np.ndarray, approximate_contours: bool, scaling_fac
     for i, contour in enumerate(contours):
         contours[i] = [[(contour[i][0][0] / scaling_factor), (contour[i][0][1] / scaling_factor)] for i in
                        range(0, len(contour))]
-    # hierarchy = hierarchy[0]  # Format extra array out of data
+    
+    hierarchy = hierarchy[0]  # Format extra array out of data
 
     return contours, hierarchy
 
diff --git a/rt_utils/smoothing.py b/rt_utils/smoothing.py
index a22c503..6e95041 100644
--- a/rt_utils/smoothing.py
+++ b/rt_utils/smoothing.py
@@ -1,17 +1,17 @@
 import SimpleITK as sitk
 import numpy as np
 from scipy import ndimage, signal
-from typing import Union, Tuple, Dict
+from typing import List, Union, Tuple, Dict
 import logging
 
 default_smoothing_parameters = {
     "iterations": 3,
+    "crop_margins": [10, 10, 1],
     "np_kron": {"scaling_factor": 2},
     "ndimage_gaussian_filter": {"sigma": 2,
                                 "radius": 3},
     "threshold": {"threshold": 0.4},
-    "ndimage_median_filter": {"size": 3,
-                              "mode": "nearest"}
+    "ndimage_median_filter": {"size": 3}
 }
 
 def kron_upscale(mask: np.ndarray, **kwargs):
@@ -27,26 +27,29 @@ def binary_threshold(mask: np.ndarray, **kwargs):
 def median_filter(mask: np.ndarray, **kwargs):
     return ndimage.median_filter(mask.astype(float), **kwargs)
 
-def crop_mask(mask: np.ndarray):
-    three_d = (len(mask.shape) == 3)
-    if three_d:
-        x, y, z = np.nonzero(mask)
-        x_max = (x.max() + 1) if x.max() < mask.shape[0] else x.max()
-        y_max = (y.max() + 1) if y.max() < mask.shape[1] else y.max()
-        z_max = (z.max() + 1) if z.max() < mask.shape[2] else z.max()
-        bbox = np.array([x.min(), x_max, y.min(), y_max, z.min(), z_max])
-
-        return mask[bbox[0]: bbox[1],
-                    bbox[2]: bbox[3],
-                    bbox[4]: bbox[5]], bbox
-    else:
-        x, y = np.nonzero(mask)
-        x_max = (x.max() + 1) if x.max() < mask.shape[0] else x.max()
-        y_max = (y.max() + 1) if y.max() < mask.shape[1] else y.max()
-        bbox = [x.min(), x_max, y.min(), y_max]
-
-        return mask[bbox[0]: bbox[1],
-                    bbox[2]: bbox[3]], bbox
+def get_new_margin(column, margin, column_length):
+    new_min = column.min() - margin
+    if new_min < 0:
+        new_min = 0
+    
+    new_max = column.max() + margin
+    if new_max > column_length:
+        new_max = column_length
+
+    return new_min, new_max
+    
+def crop_mask(mask: np.ndarray, crop_margins: List[int]):  
+    x, y, z = np.nonzero(mask)
+
+    x_min, x_max = get_new_margin(x, crop_margins[0], mask.shape[0])
+    y_min, y_max = get_new_margin(y, crop_margins[1], mask.shape[1])
+    z_min, z_max = get_new_margin(z, crop_margins[2], mask.shape[2])
+    
+    bbox = np.array([x_min, x_max, y_min, y_max, z_min, z_max])
+
+    return mask[bbox[0]: bbox[1],
+                bbox[2]: bbox[3],
+                bbox[4]: bbox[5]], bbox
 
 def restore_mask_dimensions(cropped_mask: np.ndarray, new_shape, bbox):
     new_mask = np.zeros(new_shape)
@@ -56,14 +59,15 @@ def restore_mask_dimensions(cropped_mask: np.ndarray, new_shape, bbox):
 
 def iteration_2d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, ndimage_median_filter):
     cropped_mask = kron_upscale(mask=mask, **np_kron)
-
+    return cropped_mask
     for z_idx in range(cropped_mask.shape[2]):
-        slice = cropped_mask[:, :, z_idx]
+        slice = cropped_mask[:, :, z_idx]            
         slice = gaussian_blur(mask=slice, **ndimage_gaussian_filter)
         slice = binary_threshold(mask=slice, **threshold)
         slice = median_filter(mask=slice, **ndimage_median_filter)
 
         cropped_mask[:, :, z_idx] = slice
+
     return cropped_mask
 
 def iteration_3d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, ndimage_median_filter):
@@ -82,40 +86,41 @@ def pipeline(mask: np.ndarray,
         smoothing_parameters = default_smoothing_parameters
 
     iterations = smoothing_parameters["iterations"]
+    crop_margins = smoothing_parameters["crop_margins"]
+
     np_kron = smoothing_parameters["np_kron"]
     ndimage_gaussian_filter = smoothing_parameters["ndimage_gaussian_filter"]
     threshold = smoothing_parameters["threshold"]
     ndimage_median_filter = smoothing_parameters["ndimage_median_filter"]
 
-    logging.info(f"Original mask shape {mask.shape}")
-    logging.info(f"Cropping mask to non-zero")
-    cropped_mask, bbox = crop_mask(mask)
+    print(f"Original mask shape {mask.shape}")
+    print(f"Cropping mask to non-zero")
+    cropped_mask, bbox = crop_mask(mask, crop_margins=crop_margins)
     final_shape, final_bbox = get_final_mask_shape_and_bbox(mask=mask,
                                                             scaling_factor=np_kron["scaling_factor"],
                                                             iterations=iterations,
                                                             bbox=bbox)
-    logging.info(f"Final scaling with factor: {np_kron['scaling_factor']} in {iterations} iterations")
+    print(f"Final scaling with factor of {np_kron['scaling_factor']} for {iterations} iterations")
     for i in range(iterations):
-        logging.info(f"Iteration {i} out of {iterations}")
-        logging.info(f"Applying filters")
-
+        print(f"Iteration {i+1} out of {iterations}")
+        print(f"Applying filters")
         if apply_smoothing == "2d":
-            mask = iteration_2d(mask,
+            cropped_mask = iteration_2d(cropped_mask,
                                 np_kron=np_kron,
                                 ndimage_gaussian_filter=ndimage_gaussian_filter,
                                 threshold=threshold,
                                 ndimage_median_filter=ndimage_median_filter)
         elif apply_smoothing == "3d":
-            mask = iteration_3d(mask,
+            cropped_mask = iteration_3d(cropped_mask,
                                 np_kron=np_kron,
                                 ndimage_gaussian_filter=ndimage_gaussian_filter,
                                 threshold=threshold,
                                 ndimage_median_filter=ndimage_median_filter)
         else:
             raise Exception("Wrong dimension parameter. Use '2d' or '3d'.")
-
+ 
     # Restore dimensions
-    logging.info("Restoring original mask shape")
+    print("Restoring original mask shape")
     mask = restore_mask_dimensions(cropped_mask, final_shape, final_bbox)
     return mask
 
@@ -125,8 +130,8 @@ def get_final_mask_shape_and_bbox(mask, bbox, scaling_factor, iterations):
     final_shape = np.array(mask.shape)
     final_shape[:2] *= final_scaling_factor
     bbox[:4] *= final_scaling_factor
-    logging.info("Final shape: ", final_shape)
-    logging.info("Final bbox: ", bbox)
+    print("Final shape: ", final_shape)
+    print("Final bbox: ", bbox)
     return final_shape, bbox
 
 

From 4efd7939f5047d617dad24d4d7c7e29f4ec72bd6 Mon Sep 17 00:00:00 2001
From: Mathis Rasmussen <mathis@oncology.au.dk>
Date: Wed, 14 Jun 2023 14:57:15 +0200
Subject: [PATCH 06/10] Default Parameters works really well

---
 rt_utils/smoothing.py | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/rt_utils/smoothing.py b/rt_utils/smoothing.py
index 6e95041..f5dc398 100644
--- a/rt_utils/smoothing.py
+++ b/rt_utils/smoothing.py
@@ -5,9 +5,10 @@
 import logging
 
 default_smoothing_parameters = {
-    "iterations": 3,
-    "crop_margins": [10, 10, 1],
-    "np_kron": {"scaling_factor": 2},
+    "iterations": 2,
+    "crop_margins": [20, 20, 1],
+    "initial_shift": [0, 0, 0],
+    "np_kron": {"scaling_factor": 3},
     "ndimage_gaussian_filter": {"sigma": 2,
                                 "radius": 3},
     "threshold": {"threshold": 0.4},
@@ -38,7 +39,7 @@ def get_new_margin(column, margin, column_length):
 
     return new_min, new_max
     
-def crop_mask(mask: np.ndarray, crop_margins: List[int]):  
+def crop_mask(mask: np.ndarray, crop_margins: np.ndarray):  
     x, y, z = np.nonzero(mask)
 
     x_min, x_max = get_new_margin(x, crop_margins[0], mask.shape[0])
@@ -59,7 +60,7 @@ def restore_mask_dimensions(cropped_mask: np.ndarray, new_shape, bbox):
 
 def iteration_2d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, ndimage_median_filter):
     cropped_mask = kron_upscale(mask=mask, **np_kron)
-    return cropped_mask
+
     for z_idx in range(cropped_mask.shape[2]):
         slice = cropped_mask[:, :, z_idx]            
         slice = gaussian_blur(mask=slice, **ndimage_gaussian_filter)
@@ -86,8 +87,7 @@ def pipeline(mask: np.ndarray,
         smoothing_parameters = default_smoothing_parameters
 
     iterations = smoothing_parameters["iterations"]
-    crop_margins = smoothing_parameters["crop_margins"]
-
+    crop_margins = np.array(smoothing_parameters["crop_margins"])
     np_kron = smoothing_parameters["np_kron"]
     ndimage_gaussian_filter = smoothing_parameters["ndimage_gaussian_filter"]
     threshold = smoothing_parameters["threshold"]
@@ -129,7 +129,9 @@ def get_final_mask_shape_and_bbox(mask, bbox, scaling_factor, iterations):
 
     final_shape = np.array(mask.shape)
     final_shape[:2] *= final_scaling_factor
+
     bbox[:4] *= final_scaling_factor
+    bbox[:4] -= round(final_scaling_factor * 0.5)  # Shift organ 
     print("Final shape: ", final_shape)
     print("Final bbox: ", bbox)
     return final_shape, bbox

From bb833e84ae2ba91d0f13ac2bafd54b2e5b701d86 Mon Sep 17 00:00:00 2001
From: Mathis Rasmussen <mathis@oncology.au.dk>
Date: Thu, 15 Jun 2023 08:56:58 +0200
Subject: [PATCH 07/10] Cut out redundant filter

---
 rt_utils/smoothing.py | 108 ++++++++++++++++++++++++++++--------------
 1 file changed, 72 insertions(+), 36 deletions(-)

diff --git a/rt_utils/smoothing.py b/rt_utils/smoothing.py
index f5dc398..7d042d4 100644
--- a/rt_utils/smoothing.py
+++ b/rt_utils/smoothing.py
@@ -4,31 +4,53 @@
 from typing import List, Union, Tuple, Dict
 import logging
 
+# A set of parameters that is know to work well
 default_smoothing_parameters = {
     "iterations": 2,
     "crop_margins": [20, 20, 1],
-    "initial_shift": [0, 0, 0],
     "np_kron": {"scaling_factor": 3},
     "ndimage_gaussian_filter": {"sigma": 2,
                                 "radius": 3},
     "threshold": {"threshold": 0.4},
-    "ndimage_median_filter": {"size": 3}
 }
+# A set of parameters that is know to work well
+default_smoothing_parameters_2 = {
+    "iterations": 3,
+    "crop_margins": [20, 20, 1],
+    "np_kron": {"scaling_factor": 2},
+    "ndimage_gaussian_filter": {"sigma": 2,
+                                "radius": 5},
+    "threshold": {"threshold": 0.4},
+}
+
 
-def kron_upscale(mask: np.ndarray, **kwargs):
-    scaling_factor = (kwargs["scaling_factor"], kwargs["scaling_factor"], 1)
-    return np.kron(mask, np.ones(scaling_factor))
+def kron_upscale(mask: np.ndarray, params):
+    """
+    This function upscales masks like so
 
-def gaussian_blur(mask: np.ndarray, **kwargs):
-    return ndimage.gaussian_filter(mask, **kwargs)
+    1|2       1|1|2|2
+    3|4  -->  1|1|2|2
+              3|3|4|4
+              3|3|4|4
+    
+    Scaling only in x and y direction          
+    """
+
+    scaling_array = (params["scaling_factor"], params["scaling_factor"], 1)
+    
+    return np.kron(mask, np.ones(scaling_array))
 
-def binary_threshold(mask: np.ndarray, **kwargs):
-    return mask > kwargs["threshold"]
+def gaussian_blur(mask: np.ndarray, params):
+    return ndimage.gaussian_filter(mask, **params)
 
-def median_filter(mask: np.ndarray, **kwargs):
-    return ndimage.median_filter(mask.astype(float), **kwargs)
+def binary_threshold(mask: np.ndarray, params):
+    return mask > params["threshold"]
 
 def get_new_margin(column, margin, column_length):
+    """
+    This functions takes a column (of x, y, or z) coordinates and adds a margin.
+    If margin exceeds mask size, the margin is returned to most extreme possible values
+    """
     new_min = column.min() - margin
     if new_min < 0:
         new_min = 0
@@ -39,7 +61,11 @@ def get_new_margin(column, margin, column_length):
 
     return new_min, new_max
     
-def crop_mask(mask: np.ndarray, crop_margins: np.ndarray):  
+def crop_mask(mask: np.ndarray, crop_margins: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: 
+    """
+    This function crops masks to non-zero pixels padded by crop_margins.
+    Returns (cropped mask, bounding box)
+    """ 
     x, y, z = np.nonzero(mask)
 
     x_min, x_max = get_new_margin(x, crop_margins[0], mask.shape[0])
@@ -53,36 +79,45 @@ def crop_mask(mask: np.ndarray, crop_margins: np.ndarray):
                 bbox[4]: bbox[5]], bbox
 
 def restore_mask_dimensions(cropped_mask: np.ndarray, new_shape, bbox):
+    """
+    This funtion restores mask dimentions to the given shape.
+    """
     new_mask = np.zeros(new_shape)
 
     new_mask[bbox[0]: bbox[1], bbox[2]: bbox[3], bbox[4]: bbox[5]] = cropped_mask
     return new_mask.astype(bool)
 
 def iteration_2d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, ndimage_median_filter):
-    cropped_mask = kron_upscale(mask=mask, **np_kron)
+    """
+    This is the actual set of filters. Applied iterative over z direction
+    """
+    cropped_mask = kron_upscale(mask=mask, params=np_kron)
 
     for z_idx in range(cropped_mask.shape[2]):
         slice = cropped_mask[:, :, z_idx]            
-        slice = gaussian_blur(mask=slice, **ndimage_gaussian_filter)
-        slice = binary_threshold(mask=slice, **threshold)
-        slice = median_filter(mask=slice, **ndimage_median_filter)
+        slice = gaussian_blur(mask=slice, params=ndimage_gaussian_filter)
+        slice = binary_threshold(mask=slice, params=threshold)
 
         cropped_mask[:, :, z_idx] = slice
 
     return cropped_mask
 
 def iteration_3d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, ndimage_median_filter):
-    cropped_mask = kron_upscale(mask=mask, **np_kron)
-    cropped_mask = gaussian_blur(mask=cropped_mask, **ndimage_gaussian_filter)
-    cropped_mask = binary_threshold(mask=cropped_mask, **threshold)
-    cropped_mask = median_filter(mask=cropped_mask, **ndimage_median_filter)
+    """
+    This is the actual filters applied iteratively in 3d.
+    """
+    cropped_mask = kron_upscale(mask=mask, params=np_kron)
+    cropped_mask = gaussian_blur(mask=cropped_mask, params=ndimage_gaussian_filter)
+    cropped_mask = binary_threshold(mask=cropped_mask, params=threshold)
 
     return cropped_mask
 
 def pipeline(mask: np.ndarray,
              apply_smoothing: str,
              smoothing_parameters: Union[Dict, None]):
-
+    """
+    This is the entrypoint for smoothing a mask.
+    """
     if not smoothing_parameters:
         smoothing_parameters = default_smoothing_parameters
 
@@ -91,49 +126,50 @@ def pipeline(mask: np.ndarray,
     np_kron = smoothing_parameters["np_kron"]
     ndimage_gaussian_filter = smoothing_parameters["ndimage_gaussian_filter"]
     threshold = smoothing_parameters["threshold"]
-    ndimage_median_filter = smoothing_parameters["ndimage_median_filter"]
 
-    print(f"Original mask shape {mask.shape}")
-    print(f"Cropping mask to non-zero")
+    logging.info(f"Original mask shape {mask.shape}")
+    logging.info(f"Cropping mask to non-zero")
     cropped_mask, bbox = crop_mask(mask, crop_margins=crop_margins)
     final_shape, final_bbox = get_final_mask_shape_and_bbox(mask=mask,
                                                             scaling_factor=np_kron["scaling_factor"],
                                                             iterations=iterations,
                                                             bbox=bbox)
-    print(f"Final scaling with factor of {np_kron['scaling_factor']} for {iterations} iterations")
+    logging.info(f"Final scaling with factor of {np_kron['scaling_factor']} for {iterations} iterations")
     for i in range(iterations):
-        print(f"Iteration {i+1} out of {iterations}")
-        print(f"Applying filters")
+        logging.info(f"Iteration {i+1} out of {iterations}")
+        logging.info(f"Applying filters")
         if apply_smoothing == "2d":
             cropped_mask = iteration_2d(cropped_mask,
                                 np_kron=np_kron,
                                 ndimage_gaussian_filter=ndimage_gaussian_filter,
-                                threshold=threshold,
-                                ndimage_median_filter=ndimage_median_filter)
+                                threshold=threshold)
         elif apply_smoothing == "3d":
             cropped_mask = iteration_3d(cropped_mask,
                                 np_kron=np_kron,
                                 ndimage_gaussian_filter=ndimage_gaussian_filter,
-                                threshold=threshold,
-                                ndimage_median_filter=ndimage_median_filter)
+                                threshold=threshold)
         else:
             raise Exception("Wrong dimension parameter. Use '2d' or '3d'.")
  
     # Restore dimensions
-    print("Restoring original mask shape")
+    logging.info("Restoring original mask shape")
     mask = restore_mask_dimensions(cropped_mask, final_shape, final_bbox)
     return mask
 
 def get_final_mask_shape_and_bbox(mask, bbox, scaling_factor, iterations):
+    """
+    This function scales image shape and the bounding box which should be used for the final mask
+    """
+
     final_scaling_factor = pow(scaling_factor, iterations)
 
     final_shape = np.array(mask.shape)
     final_shape[:2] *= final_scaling_factor
 
-    bbox[:4] *= final_scaling_factor
-    bbox[:4] -= round(final_scaling_factor * 0.5)  # Shift organ 
-    print("Final shape: ", final_shape)
-    print("Final bbox: ", bbox)
+    bbox[:4] *= final_scaling_factor # Scale bounding box to final shape
+    bbox[:4] -= round(final_scaling_factor * 0.5)  # Shift volumes to account for the shift that occurs as a result of the scaling
+    logging.info("Final shape: ", final_shape)
+    logging.info("Final bbox: ", bbox)
     return final_shape, bbox
 
 

From 12c3ba7caea994fc886f004908a11b9ade92814b Mon Sep 17 00:00:00 2001
From: Mathis Rasmussen <mathis@oncology.au.dk>
Date: Thu, 15 Jun 2023 09:05:26 +0200
Subject: [PATCH 08/10] Add comments

---
 rt_utils/rtstruct.py  | 12 +++++++++---
 rt_utils/smoothing.py |  4 ++--
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/rt_utils/rtstruct.py b/rt_utils/rtstruct.py
index 6b0582a..55c204d 100644
--- a/rt_utils/rtstruct.py
+++ b/rt_utils/rtstruct.py
@@ -35,8 +35,14 @@ def add_roi(
         use_pin_hole: bool = False,
         approximate_contours: bool = True,
         roi_generation_algorithm: Union[str, int] = 0,
-        apply_smoothing: Union[str, None] = None, # strings can be "2d" or "3d"
-        smoothing_parameters: Union[Dict, None] = None
+        apply_smoothing: Union[str, None] = None, # strings can be "2d" or "3d" or something else if a different smoothing function is used
+        smoothing_function = smoothing.pipeline,  # Can be any function/set of functions that takes the following parameters
+            #                                     # smoothing_function(mask=mask, apply_smoothing=apply_smoothing,
+                                                  #                    smoothing_parameters=smoothing_parameters) -> np.ndarray
+                                                  # The returned np.ndarray can be of any integer scalar shape in x and y of the used dicom image.
+                                                  # Note that Z direction should not be scaled. For instance CT_image.shape == (512, 512, 150).
+                                                  # Smoothed returned array can be (1024, 1024, 150) or (5120, 5120, 150), though you RAM will suffer with the latter.
+        smoothing_parameters: Union[Dict, None] = None,
     ):
         """
         Add a ROI to the rtstruct given a 3D binary mask for the ROI's at each slice
@@ -45,7 +51,7 @@ def add_roi(
         If approximate_contours is set to False, no approximation will be done when generating contour data, leading to much larger amount of contour data
         """
         if apply_smoothing:
-            mask = smoothing.pipeline(mask=mask, apply_smoothing=apply_smoothing,
+            mask = smoothing_function(mask=mask, apply_smoothing=apply_smoothing,
                                       smoothing_parameters=smoothing_parameters)
 
         ## If upscaled coords are given, they should be adjusted accordingly
diff --git a/rt_utils/smoothing.py b/rt_utils/smoothing.py
index 7d042d4..8ab6b79 100644
--- a/rt_utils/smoothing.py
+++ b/rt_utils/smoothing.py
@@ -87,7 +87,7 @@ def restore_mask_dimensions(cropped_mask: np.ndarray, new_shape, bbox):
     new_mask[bbox[0]: bbox[1], bbox[2]: bbox[3], bbox[4]: bbox[5]] = cropped_mask
     return new_mask.astype(bool)
 
-def iteration_2d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, ndimage_median_filter):
+def iteration_2d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold):
     """
     This is the actual set of filters. Applied iterative over z direction
     """
@@ -102,7 +102,7 @@ def iteration_2d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold,
 
     return cropped_mask
 
-def iteration_3d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, ndimage_median_filter):
+def iteration_3d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold):
     """
     This is the actual filters applied iteratively in 3d.
     """

From bbc5d4825daf16684eef0c4fa2903a9cbf7f420c Mon Sep 17 00:00:00 2001
From: Mathis Rasmussen <mathis@oncology.au.dk>
Date: Thu, 15 Jun 2023 09:06:36 +0200
Subject: [PATCH 09/10] Implement smoothing as option

---
 rt_utils/image_helper.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/rt_utils/image_helper.py b/rt_utils/image_helper.py
index 63c5b43..87b9ef6 100644
--- a/rt_utils/image_helper.py
+++ b/rt_utils/image_helper.py
@@ -100,6 +100,7 @@ def find_mask_contours(mask: np.ndarray, approximate_contours: bool, scaling_fac
         contours
     )  # Open-CV updated contours to be a tuple so we convert it back into a list here
 
+    # Coordinates are rescaled to image grid by dividing with scaling factor
     for i, contour in enumerate(contours):
         contours[i] = [[(contour[i][0][0] / scaling_factor), (contour[i][0][1] / scaling_factor)] for i in
                        range(0, len(contour))]

From cafdfe0d8a03d7b829d7bc621ce8fda9d6eee3c4 Mon Sep 17 00:00:00 2001
From: Mathis Rasmussen <mathis@oncology.au.dk>
Date: Thu, 15 Jun 2023 09:06:36 +0200
Subject: [PATCH 10/10] Implement smoothing as option

---
 rt_utils/image_helper.py |  1 +
 rt_utils/smoothing.py    | 55 ++++++++++++++++++++++------------------
 2 files changed, 32 insertions(+), 24 deletions(-)

diff --git a/rt_utils/image_helper.py b/rt_utils/image_helper.py
index 63c5b43..87b9ef6 100644
--- a/rt_utils/image_helper.py
+++ b/rt_utils/image_helper.py
@@ -100,6 +100,7 @@ def find_mask_contours(mask: np.ndarray, approximate_contours: bool, scaling_fac
         contours
     )  # Open-CV updated contours to be a tuple so we convert it back into a list here
 
+    # Coordinates are rescaled to image grid by dividing with scaling factor
     for i, contour in enumerate(contours):
         contours[i] = [[(contour[i][0][0] / scaling_factor), (contour[i][0][1] / scaling_factor)] for i in
                        range(0, len(contour))]
diff --git a/rt_utils/smoothing.py b/rt_utils/smoothing.py
index 8ab6b79..9731cff 100644
--- a/rt_utils/smoothing.py
+++ b/rt_utils/smoothing.py
@@ -5,8 +5,9 @@
 import logging
 
 # A set of parameters that is know to work well
-default_smoothing_parameters = {
-    "iterations": 2,
+default_smoothing_parameters_2d = {
+    "scaling_iterations": 2,
+    "filter_iterations": 1,
     "crop_margins": [20, 20, 1],
     "np_kron": {"scaling_factor": 3},
     "ndimage_gaussian_filter": {"sigma": 2,
@@ -14,8 +15,8 @@
     "threshold": {"threshold": 0.4},
 }
 # A set of parameters that is know to work well
-default_smoothing_parameters_2 = {
-    "iterations": 3,
+default_smoothing_parameters_2d_2 = {
+    "scaling_iterations": 3,
     "crop_margins": [20, 20, 1],
     "np_kron": {"scaling_factor": 2},
     "ndimage_gaussian_filter": {"sigma": 2,
@@ -87,28 +88,30 @@ def restore_mask_dimensions(cropped_mask: np.ndarray, new_shape, bbox):
     new_mask[bbox[0]: bbox[1], bbox[2]: bbox[3], bbox[4]: bbox[5]] = cropped_mask
     return new_mask.astype(bool)
 
-def iteration_2d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold):
+def iteration_2d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, filter_iterations):
     """
     This is the actual set of filters. Applied iterative over z direction
     """
     cropped_mask = kron_upscale(mask=mask, params=np_kron)
 
-    for z_idx in range(cropped_mask.shape[2]):
-        slice = cropped_mask[:, :, z_idx]            
-        slice = gaussian_blur(mask=slice, params=ndimage_gaussian_filter)
-        slice = binary_threshold(mask=slice, params=threshold)
+    for filter_iteration in range(filter_iterations):
+        for z_idx in range(cropped_mask.shape[2]):
+            slice = cropped_mask[:, :, z_idx]            
+            slice = gaussian_blur(mask=slice, params=ndimage_gaussian_filter)
+            slice = binary_threshold(mask=slice, params=threshold)
 
-        cropped_mask[:, :, z_idx] = slice
+            cropped_mask[:, :, z_idx] = slice
 
     return cropped_mask
 
-def iteration_3d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold):
+def iteration_3d(mask: np.ndarray, np_kron, ndimage_gaussian_filter, threshold, filter_iterations):
     """
     This is the actual filters applied iteratively in 3d.
     """
-    cropped_mask = kron_upscale(mask=mask, params=np_kron)
-    cropped_mask = gaussian_blur(mask=cropped_mask, params=ndimage_gaussian_filter)
-    cropped_mask = binary_threshold(mask=cropped_mask, params=threshold)
+    for filter_iteration in range(filter_iterations):
+        cropped_mask = kron_upscale(mask=mask, params=np_kron)
+        cropped_mask = gaussian_blur(mask=cropped_mask, params=ndimage_gaussian_filter)
+        cropped_mask = binary_threshold(mask=cropped_mask, params=threshold)
 
     return cropped_mask
 
@@ -119,9 +122,11 @@ def pipeline(mask: np.ndarray,
     This is the entrypoint for smoothing a mask.
     """
     if not smoothing_parameters:
-        smoothing_parameters = default_smoothing_parameters
+        smoothing_parameters = default_smoothing_parameters_2d
+
+    scaling_iterations = smoothing_parameters["scaling_iterations"]
+    filter_iterations = smoothing_parameters["filter_iterations"]
 
-    iterations = smoothing_parameters["iterations"]
     crop_margins = np.array(smoothing_parameters["crop_margins"])
     np_kron = smoothing_parameters["np_kron"]
     ndimage_gaussian_filter = smoothing_parameters["ndimage_gaussian_filter"]
@@ -132,22 +137,24 @@ def pipeline(mask: np.ndarray,
     cropped_mask, bbox = crop_mask(mask, crop_margins=crop_margins)
     final_shape, final_bbox = get_final_mask_shape_and_bbox(mask=mask,
                                                             scaling_factor=np_kron["scaling_factor"],
-                                                            iterations=iterations,
+                                                            scaling_iterations=scaling_iterations,
                                                             bbox=bbox)
-    logging.info(f"Final scaling with factor of {np_kron['scaling_factor']} for {iterations} iterations")
-    for i in range(iterations):
-        logging.info(f"Iteration {i+1} out of {iterations}")
+    logging.info(f"Final scaling with factor of {np_kron['scaling_factor']} for {scaling_iterations} scaling_iterations")
+    for i in range(scaling_iterations):
+        logging.info(f"Iteration {i+1} out of {scaling_iterations}")
         logging.info(f"Applying filters")
         if apply_smoothing == "2d":
             cropped_mask = iteration_2d(cropped_mask,
                                 np_kron=np_kron,
                                 ndimage_gaussian_filter=ndimage_gaussian_filter,
-                                threshold=threshold)
+                                threshold=threshold,
+                                filter_iterations=filter_iterations)
         elif apply_smoothing == "3d":
             cropped_mask = iteration_3d(cropped_mask,
                                 np_kron=np_kron,
                                 ndimage_gaussian_filter=ndimage_gaussian_filter,
-                                threshold=threshold)
+                                threshold=threshold,
+                                filter_iterations=filter_iterations)
         else:
             raise Exception("Wrong dimension parameter. Use '2d' or '3d'.")
  
@@ -156,12 +163,12 @@ def pipeline(mask: np.ndarray,
     mask = restore_mask_dimensions(cropped_mask, final_shape, final_bbox)
     return mask
 
-def get_final_mask_shape_and_bbox(mask, bbox, scaling_factor, iterations):
+def get_final_mask_shape_and_bbox(mask, bbox, scaling_factor, scaling_iterations):
     """
     This function scales image shape and the bounding box which should be used for the final mask
     """
 
-    final_scaling_factor = pow(scaling_factor, iterations)
+    final_scaling_factor = pow(scaling_factor, scaling_iterations)
 
     final_shape = np.array(mask.shape)
     final_shape[:2] *= final_scaling_factor