adding badPixelCorrection to preprocessing algorithms

demo_leapctype/d14_outlierCorrection.py

@@ -7,7 +7,7 @@
 from leap_preprocessing_algorithms import *
 
 '''
-The script demonstrates two different methods to handle outliers (zingers) in your projection data
+The script demonstrates two different methods to handle outliers (zingers) and bad pixels in your projection data
 '''
 
 # Specify the number of detector columns which is used below
@@ -48,24 +48,44 @@
 I_0 = 50000.0
 g[:] = -np.log(np.random.poisson(I_0*np.exp(-g))/I_0)
 
-# Make some of the detector pixels have zero value
-ind = np.abs(np.random.normal(0,1,g.shape)) > 3.0
-g[ind] = 0.0
-
 # Choose which method you'd like to test
-whichMethod = 1
-#whichMethod = 2
+#whichMethod = 1
+whichMethod = 2
+#whichMethod = 3
 
 if whichMethod == 1:
+
+    # Make some of the detector pixels have zero value
+    ind = np.abs(np.random.normal(0,1,g.shape)) > 3.0
+    g[ind] = 0.0
+
     # Perform median filter based outlier correction
     # This is a very fast method which is good for isolated bad pixels
     outlierCorrection(leapct,g)
 
+    # Reconstruct the data
+    f = leapct.allocateVolume()
+    leapct.FBP(g,f)
+elif whichMethod == 2:
 
+    # Make some detector pixels dead for all projections
+    ind = np.abs(np.random.normal(0,1,(g.shape[1], g.shape[2]))) > 3.0
+    g[:,ind] = 0.0
+
+    # Correct for bad pixels
+    badPixelMap = np.zeros((g.shape[1], g.shape[2]), dtype=np.float32)
+    badPixelMap[ind] = 1.0
+    badPixelCorrection(leapct, g, badPixelMap)
+
     # Reconstruct the data
     f = leapct.allocateVolume()
     leapct.FBP(g,f)
 else:
+
+    # Make some of the detector pixels have zero value
+    ind = np.abs(np.random.normal(0,1,g.shape)) > 3.0
+    g[ind] = 0.0
+
     # Perform RWLS reconstruction where we set the weights to be zero where there are bad pixels
     # This method is computationally expensive, but is good when you have large regions of bad pixels
     # or want to perform so-called Metal Artifact Reduction (MAR)

docs/source/preprocessing_algorithms.rst

@@ -2,6 +2,7 @@ Preprocessing Algorithms
 ========================
 
 .. autofunction:: leap_preprocessing_algorithms.makeAttenuationRadiographs
+.. autofunction:: leap_preprocessing_algorithms.badPixelCorrection
 .. autofunction:: leap_preprocessing_algorithms.outlierCorrection
 .. autofunction:: leap_preprocessing_algorithms.outlierCorrection_highEnergy
 .. autofunction:: leap_preprocessing_algorithms.detectorDeblur_FourierDeconv

src/leap_preprocessing_algorithms.py

@@ -87,10 +87,39 @@ def makeAttenuationRadiographs(leapct, g, air_scan=None, dark_scan=None, ROI=Non
 
     return True
 
+def badPixelCorrection(leapct, g, badPixelMap, windowSize=3, isAttenuationData=True):
+    r"""Removes bad pixels from CT projections
+    
+    LEAP CT geometry parameters must be set prior to running this function 
+    and can be applied to any CT geometry type.
+    This algorithm processes each projection independently
+    and removes bad pixels specified by the user using a median filter
+    
+    Args:
+        leapct (tomographicModels object): This is just needed to access LEAP algorithms
+        g (contiguous float32 numpy array or torch tensor): attenuation or transmission projection data
+        badPixelMap (C contiguous float32 numpy array or torch tensor): 2D bad pixel map (numRows x numCols) where a value of 1.0 marks a pixel as bad
+        windowSize (int): the window size; can be 3, 5, or 7
+        isAttenuationData (bool): True if g is attenuation data, False otherwise
+        
+    Returns:
+        True if successful, False otherwise
+    """
+
+    if g is None or badPixelMap is None:
+        return False
+
+    # This algorithm processes each transmission
+    if isAttenuationData:
+        leapct.expNeg(g)
+    leapct.badPixelCorrection(g, badPixelMap, windowSize)
+    if isAttenuationData:
+        leapct.negLog(g)
+    return True
+
 def outlierCorrection(leapct, g, threshold=0.03, windowSize=3, isAttenuationData=True):
     r"""Removes outliers (zingers) from CT projections
     
-    Assumes the input data is in attenuation space.
     No LEAP parameters need to be set for this function to work 
     and can be applied to any CT geometry type.
     This algorithm processes each projection independently
@@ -121,7 +150,6 @@ def outlierCorrection(leapct, g, threshold=0.03, windowSize=3, isAttenuationData
 def outlierCorrection_highEnergy(leapct, g, isAttenuationData=True):
     """Removes outliers (zingers) from CT projections
     
-    Assumes the input data is in attenuation space.
     No LEAP parameters need to be set for this function to work 
     and can be applied to any CT geometry type.
     This algorithm processes each projection independently
@@ -152,7 +180,6 @@ def outlierCorrection_highEnergy(leapct, g, isAttenuationData=True):
 def LowSignalCorrection(leapct, g, threshold=0.03, windowSize=3, signalThreshold=0.001, isAttenuationData=True):
     r"""Corrects detector pixels that have very low transmission (photon starvation)
     
-    Assumes the input data is in attenuation space.
     No LEAP parameters need to be set for this function to work 
     and can be applied to any CT geometry type.
     This algorithm processes each projection independently

src/leapctype.py

@@ -4543,10 +4543,12 @@ def badPixelCorrection(self, g, badPixelMap, windowSize=3):
         r"""Bad Pixel Correction
         
         The provided input must be projection data and the CT geometry parameters must be set.
+        This algorithm processes each projection independently
+        and removes bad pixels specified by the user using a median filter.
         
         Args:
             g (C contiguous float32 numpy array or torch tensor): 3D projection data array
-            badPixelMap (C contiguous float32 numpy array or torch tensor): 2D bad pixel map
+            badPixelMap (C contiguous float32 numpy array or torch tensor): 2D bad pixel map (numRows x numCols)
             windowSize (int): the window size; can be 3, 5, or 7
         
         Returns: