Merge branch 'minnititino-master'

added BHC and zero_clipping_value functions

neutompy/preproc/preproc.py

@@ -16,7 +16,7 @@
 import sys
 from time import sleep
 
-__author__  = "Davide Micieli"
+__author__  = "Davide Micieli, Triestino Minniti"
 __all__     = ['draw_ROI',
 			   'normalize_proj',
 			   'log_transform',
@@ -25,9 +25,12 @@
 			   'remove_outliers',
 			   'remove_outliers_stack',
 			   'remove_stripe',
-			   'remove_stripe_stack'
+			   'remove_stripe_stack',
+			   'simple_BHC',
+			   'zero_clipping_value'
 			  ]
 
+
 def draw_ROI(img, title, ratio=0.85):
 	"""
 	This function allows to select interactively a rectangular region of interest
@@ -1243,3 +1246,71 @@ def remove_stripe_stack(arr, level, wname='db5', sigma=1.5, axis=1, out=None):
 	out_arr = out_arr.swapaxes(0,axis)
 
 	return out_arr
+
+
+def simple_BHC(norm, a0=0., a1=0., a2=0.02, a3=0., out=None):
+    """
+	This function performs the simple beam hardening correction (BHC) corresponding
+	to a polynomial correction.
+	The user can choose 4 parameters which define a 5-th order polynomial.
+
+	s = -ln(I/I 0 )
+    s'= s + a0*s^2 + a1*s^3 + a2*s^4 + a3*s^5
+
+	Parameters
+	----------
+	norm :  3d array
+			Three-dimensional stack of the normalized projections
+
+	a0 :    float, optional
+	        The first term of the polynomial
+
+	a1 :    float, optional
+	        The second term of the polynomial
+
+	a2 :    float, optional
+	        The third term of the polynomial
+
+	a3 :    float, optional
+	        The fourth term of the polynomial
+
+	Returns
+	-------
+	out : 3d array
+	      Beam hardening correction (BHC) of the input array.
+    """
+
+    a0 = np.float32(a0)
+    a1 = np.float32(a1)
+    a2 = np.float32(a2)
+    a3 = np.float32(a3)
+    out = ne.evaluate('norm + a0*norm**2 + a1*norm**3 + a2*norm**4 + a3*norm**5', out=out)
+
+    return out
+
+
+def zero_clipping_value(norm, cl=0.01, out=None):
+    """
+	This function clips the values in an array.
+	Values below the threshold value cl are replaced with cl.
+	This function prevents values close to (or below) zero to introduce new artefacts.
+
+	s' = MAX ( s, cl)
+
+	Parameters
+	----------
+	norm :  3d array
+			Three-dimensional stack of the normalized projections
+
+	cl :    float, optional
+	        Clipping value
+
+	Returns
+	-------
+	out : 3d array
+	      Clipping correction of the input array.
+    """
+
+    cl  = np.float32(cl)
+    out = ne.evaluate( 'where(norm<cl, cl, norm)', out=out)
+    return out