Adding median3d and dezinger3d filters into tomopy.misc.corr

source/tomopy/misc/corr.py

@@ -348,25 +348,24 @@ def callback(total, description, unit):
     return arr
 
 
-def median_filter3d(arr, kernel_half_size=1, ncore=None):
+def median_filter3d(arr, size=3, ncore=None):
     """
     Apply 3D median filter to 3D array.
 
     Parameters
     ----------
     arr : ndarray
-        Input 3D array.
-    kernel_half_size : int, optional
-        The half size of the filter's kernel, i.e. 1 results in the full kernel
-        size of 3 x 3 x 3.
+        Input 3D array either float32 or uint16 data type.
+    size : int, optional
+        The size of the filter's kernel.
     ncore : int, optional
         Number of cores that will be assigned to jobs. All cores will be used
         if unspecified.
 
     Returns
     -------
     ndarray
-        Median filtered 3D array.
+        Median filtered 3D array either float32 or uint16 data type.
     Raises
     ------
     ValueError
@@ -375,7 +374,7 @@ def median_filter3d(arr, kernel_half_size=1, ncore=None):
     """
     input_type = arr.dtype
     if (input_type != 'float32') and (input_type != 'uint16'):
-        arr = dtype.as_float32(arr)  # silently convert to float32 data type
+        arr = dtype.as_float32(arr)  # silent convertion to float32 data type
     out = np.empty_like(arr)
     dif = 0.0  # set to 0 to avoid selective filtering
 
@@ -387,7 +386,16 @@ def median_filter3d(arr, kernel_half_size=1, ncore=None):
     if ncore is None:
         ncore = mproc.mp.cpu_count()
 
-    #deal with different data types
+    # convert the full kernel size (odd int) to a half size as the C function requires
+    if size < 3:
+        size = 3 # check if the kernel size is not too small
+    if (size % 2) == 0:
+        # dealing with even integers
+        kernel_half_size = (int)(0.5*size)
+    else:
+        kernel_half_size = (int)(0.5*(size-1))
+
+    # deal with different data types
     if (input_type == 'float32'):
         extern.c_median_filt3d_float32(arr, out, kernel_half_size, dif, ncore,
                                        dx, dy, dz)
@@ -397,29 +405,27 @@ def median_filter3d(arr, kernel_half_size=1, ncore=None):
     return out
 
 
-def remove_outlier3d(arr, kernel_half_size=1, dif=0.1, ncore=None):
+def remove_outlier3d(arr, dif=0.1, size=3, ncore=None):
     """
-    Also a so-called dezinger. Selectively applies 3D median filter to 3D array
-    to remove outliers specifically.
+    Selectively applies 3D median filter to a 3D array to remove outliers specifically. Also called a dezinger.
 
     Parameters
     ----------
     arr : ndarray
-        Input 3D array.
-    kernel_half_size : int, optional
-        The half size of the filter's kernel, i.e. 1 results in the full kernel
-        size of 3 x 3 x 3.
+        Input 3D array either float32 or uint16 data type.
     dif : float
         Expected difference value between outlier value and the median value of
         the array.
+    size : int, optional
+        The size of the filter's kernel.        
     ncore : int, optional
         Number of cores that will be assigned to jobs. All cores will be used
         if unspecified.
 
     Returns
     -------
     ndarray
-        Dezingered 3D array.
+        Dezingered 3D array either float32 or uint16 data type.
     Raises
     ------
     ValueError
@@ -428,7 +434,7 @@ def remove_outlier3d(arr, kernel_half_size=1, dif=0.1, ncore=None):
     """
     input_type = arr.dtype
     if (input_type != 'float32') and (input_type != 'uint16'):
-        arr = dtype.as_float32(arr)  # silently convert to float32 data type
+        arr = dtype.as_float32(arr)  # silent convertion to float32 data type
     out = np.empty_like(arr)
 
     if np.ndim(arr) == 3:
@@ -439,7 +445,16 @@ def remove_outlier3d(arr, kernel_half_size=1, dif=0.1, ncore=None):
     if ncore is None:
         ncore = mproc.mp.cpu_count()
 
-    #deal with different data types
+    # convert the full kernel size (odd int) to a half size as the C function requires
+    if size < 3:
+        size = 3 # check if the kernel size is not too small
+    if (size % 2) == 0:
+        # dealing with even integers
+        kernel_half_size = (int)(0.5*size)
+    else:
+        kernel_half_size = (int)(0.5*(size-1))
+
+    # deal with different data types
     if (input_type == 'float32'):
         extern.c_median_filt3d_float32(arr, out, kernel_half_size, dif, ncore,
                                        dx, dy, dz)

test/test_tomopy/test_misc/test_corr.py

@@ -48,7 +48,7 @@
 import unittest
 
 import numpy as np
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_equal
 from tomopy.misc.corr import (
     gaussian_filter,
     median_filter,
@@ -113,32 +113,17 @@ def test_median_filter_nonfinite(self):
                 )
 
     def test_median_filter3d(self):
-        # generate a small noncubic volume and apply random noise
-        A = np.float32(np.zeros((30,40,20)))
-        A[:] = 0.1
-        A[10:20,20:30,:] = 0.5
-        np.random.seed(1)
-        A_noise = A + np.random.normal(loc=0,
-                          scale=0.03 * A,
-                          size=np.shape(A))
-        A_filtered = median_filter3d(np.float32(A_noise))
-        max_filtered = 0.513833
-        assert_allclose(np.max(A_filtered), max_filtered, rtol=1e-06)
+        A = np.arange(4*5*6).reshape(4,5,6)
+        A_median = median_filter3d(np.float32(A))
+        assert_equal(read_file('median_filter3d.npy'), A_median)
 
     def test_remove_outlier3d(self):
-        # generate a small noncubic volume and apply random noise with an outlier
-        A = np.float32(np.zeros((30,40,20)))
-        A[:] = 0.1
-        A[10:20,20:30,:] = 0.5        
-        np.random.seed(1)
-        A_noise = A + np.random.normal(loc=0,
-                          scale=0.03 * A,
-                          size=np.shape(A))
-        A[15,25,0] = 1.5 # placing an outlier
-        A_filtered = remove_outlier3d(np.float32(A_noise), kernel_half_size=1, dif = 0.5)
-        max_filtered = 0.56252176
-        assert_allclose(np.max(A_filtered), max_filtered, rtol=1e-06)        
-
+        A = np.arange(4*5*6).reshape(4,5,6)
+        A[2,2,2] = 1000.0 # introduce an outlier
+        A_dezinged = remove_outlier3d(np.float32(A), dif = 500, size=3)
+        A[2,2,2] = 75 # substituted value by dezinger
+        assert_equal(np.float32(A), A_dezinged)
+
     def test_remove_neg(self):
         assert_allclose(remove_neg([-2, -1, 0, 1, 2]), [0, 0, 0, 1, 2])