Added Gaussian blurring to smooth structure boundary

cerr/dataclasses/structure.py

@@ -698,6 +698,33 @@ def importStructureMask(mask3M, assocScanNum, structName, planC, structNum=None)
     for slc in range(dim[2]):
         if not np.any(paddedMask3M[:,:,slc]):
             continue
+        # if upsamplingRate > 1:
+        #     img = sitk.GetImageFromArray(paddedMask3M[:,:,slc].astype(float))
+        #     resample = sitk.ResampleImageFilter()
+        #     # resample.SetOutputDirection(output_direction)
+        #     # resample.SetOutputOrigin(output_origin)
+        #     # resample.SetOutputSpacing(output_spacing)
+        #     origSiz = paddedMask3M[:,:,slc].shape
+        #     outputSize = [siz * upsamplingRate for siz in origSiz]
+        #     outputSpacing = [siz/upsamplingRate for siz in img.GetSpacing()]
+        #     resample.SetReferenceImage(img)
+        #     resample.SetSize(outputSize)
+        #     resample.SetOutputSpacing(outputSpacing)
+        #     resample.SetOutputOrigin(img.GetOrigin())
+        #     resample.SetDefaultPixelValue(0)
+        #     resample.SetInterpolator(getattr(sitk,upsamplingMethod))
+        #     sitkResampImg = resample.Execute(img)
+        #     upsampled2M = sitk.GetArrayFromImage(sitkResampImg)
+        #     gaussian = sitk.SmoothingRecursiveGaussianImageFilter()
+        #     gaussian.SetSigma(outputSpacing[0])
+        #     blurImage = gaussian.Execute(sitkResampImg)
+        #     blurMask2M = sitk.GetArrayFromImage(blurImage)
+        #     contours = measure.find_contours(blurMask2M, 0.5)
+        #     origRowV = np.linspace(0,origSiz[0])
+        #     origColV = np.linspace(0,origSiz[1])
+        #     upsampRowV = np.linspace(0,outputSize[0])
+        #     upsampColV = np.linspace(0,outputSize[1])
+        # else:
         contours = measure.find_contours(paddedMask3M[:,:,slc] == 1, 0.5)
         if scn.flipSliceOrderFlag(planC.scan[assocScanNum]):
             slcNum = numSlcs - slc - 1
@@ -716,6 +743,10 @@ def importStructureMask(mask3M, assocScanNum, structName, planC, structNum=None)
         for iContour, contour in enumerate(contours):
             contObj = Contour()
             segment = np.empty((contour.shape[0],3))
+            # if upsamplingRate > 1:
+            #     colV = np.interp(contour[:, 1] , upsampColV, origColV, 0, 0) - 1
+            #     rowV = np.interp(contour[:, 0] , upsampRowV, origRowV, 0, 0) - 1
+            # else:
             colV = contour[:, 1] - 1 # - 1 to account for padding
             rowV = contour[:, 0] - 1 # - 1 to account for padding
             slcV = np.ones_like(rowV) * slcNum
@@ -1048,6 +1079,53 @@ def getLargestConnComps(structNum, numConnComponents, planC=None, saveFlag=None,
 
     return maskOut3M, planC
 
+def getGaussianBlurredMask(structNum, sigmaVoxel, planC, saveFlag=False,\
+              replaceFlag=None, procSructName=None):
+    """
+    Function for Gaussian blurring of binary masks
+
+    Args:
+        structNum (int): int for index of structure in planC.
+        sigmaVoxel (float): sigma used in Gaussian in units of number of voxels.
+         corresponding to each image dimension.
+        planC: pyCERR plan container object.
+        saveFlag: [optional, default=False] bool flag for saving processed mask to planC.
+        replaceFlag: [optional, default=False] bool flag for replacing input mask with
+                    processed mask in planC.
+        procSructName: [optional, default=None] string for output structure name.
+                      Original structure name is used if None.
+    Returns:
+        filledMask3M: np.ndarray(dtype=bool) for filled mask.
+        planC: pyCERR plan container object.
+    """
+
+    # Get binary mask of structure
+    mask3M = rs.getStrMask(structNum,planC)
+
+    # Get mask resolution
+    assocScanNum = scn.getScanNumFromUID(planC.structure[structNum].assocScanUID,\
+                                              planC)
+    #inputResV = planC.scan[assocScanNum].getScanSpacing() * sigmaVoxel
+
+    blurredMask3M = maskUtils.gaussianBlurring(mask3M,sigmaVoxel)
+
+    # Save to planC
+    if saveFlag:
+        if procSructName is None:
+            procSructName = planC.structure[structNum].structureName
+
+        assocScanNum = scn.getScanNumFromUID(planC.structure[structNum].assocScanUID,\
+                                                  planC)
+        newStructNum = None
+        if replaceFlag:
+            # Delete structNum
+            #del planC.structure[structNum]
+            newStructNum = structNum
+        planC = importStructureMask(blurredMask3M >= 0.5, assocScanNum, procSructName, planC, newStructNum)
+
+
+    return blurredMask3M, planC
+
 
 def getLabelMap(strNumV, planC, labelDict=None, dim=3):
     """

cerr/utils/mask.py

@@ -9,6 +9,7 @@
 from skimage import exposure, filters, morphology, transform
 from skimage.morphology import square, octagon
 import cerr.utils.statistics as statUtil
+import SimpleITK as sitk
 
 
 def getDown2Mask(inM, sample):
@@ -146,11 +147,39 @@ def morphologicalClosing(binaryMask, structuringElement):
         structuringElement (np.array): Flat morphological structuring element.
 
     Returns:
-        closedMask (np.ndarray(dtype=bool)) Closed mask using input structuring element.
+        numpy.ndarray(dtype=bool): Closed mask using input structuring element.
     """
+
     closedMask = ndimage.binary_closing(binaryMask, structure=structuringElement)
     return closedMask
 
+
+def gaussianBlurring(binaryMask, sigmaVox):
+    """
+    Function for Gaussian blurring of input binary mask
+
+    Args:
+        binaryMask (numpy.array): Binary mask to blur.
+        sigmaVox (float): Sigma for Gaussian in units of voxels.
+
+    Returns:
+        numpy.ndarray(dtype=bool): Blurred mask using Gaussian blur with input sigma.
+    """
+
+    gaussian = sitk.SmoothingRecursiveGaussianImageFilter()
+    gaussian.SetSigma(sigmaVox)
+    dim = binaryMask.shape
+    blurredMask3M = np.empty_like(binaryMask, dtype=float)
+    for slc in range(dim[2]):
+        if not np.any(binaryMask[:,:,slc]):
+            blurredMask3M[:,:,slc] = binaryMask[:,:,slc]
+            continue
+        img = sitk.GetImageFromArray(binaryMask[:,:,slc].astype(float))
+        blurImage = gaussian.Execute(img)
+        blurredMask3M[:,:,slc] = sitk.GetArrayFromImage(blurImage)
+    return blurredMask3M
+
+
 def computeBoundingBox(binaryMaskM, is2DFlag=False, maskFlag=0):
     """
     Function for finding extents of bounding box given a binary mask