From 304167fe8deac184b01cb8d201adc4443883f9a0 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Tue, 20 Aug 2019 16:26:44 -0500 Subject: [PATCH 01/35] STY: Use standard logging instead of print messages --- source/tomopy/recon/algorithm.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/source/tomopy/recon/algorithm.py b/source/tomopy/recon/algorithm.py index 115fdaaa1..eefb997be 100644 --- a/source/tomopy/recon/algorithm.py +++ b/source/tomopy/recon/algorithm.py @@ -151,7 +151,7 @@ def recon( Total Variation reconstruction technique :cite:`Chambolle:11`. 'grad' - Gradient descent method. + Gradient descent method. 'tikh' Tikhonov regularization with identity Tikhonov matrix. @@ -376,11 +376,13 @@ def _dist_recon(tomo, center, recon, algorithm, args, kwargs, ncore, nchunk): # check if ncore is limited by env variable pythreads = os.environ.get("TOMOPY_PYTHON_THREADS") if pythreads is not None and ncore > int(pythreads): - print("Warning! 'TOMOPY_PYTHON_THREADS' has been set to '{0}', which is less than" - " specified ncore={1}. Limiting ncore to {0}...".format(pythreads, ncore)) + logger.warning("'TOMOPY_PYTHON_THREADS' has been set to '{0}'," + " which is less than specified ncore={1}." + " Limiting ncore to {0}...".format(pythreads, ncore)) ncore = int(pythreads) - print("Reconstructing {} slice groups with {} master threads...".format(len(slcs), ncore)) + logger.info("Reconstructing {} slice groups with {} master threads..." + "".format(len(slcs), ncore)) # this is used internally to prevent oversubscription os.environ["TOMOPY_PYTHON_THREADS"] = "{}".format(ncore) From eb339177f39eec27d647e7674f6fbf8924317496 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Tue, 20 Aug 2019 16:36:50 -0500 Subject: [PATCH 02/35] STY: NotImplementedError instead of a logger.warning(...) Python has a specific exception for unimplemented functions. https://docs.python.org/3/library/exceptions.html#NotImplementedError --- source/tomopy/sim/project.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/source/tomopy/sim/project.py b/source/tomopy/sim/project.py index 85566dde9..d0d205d3e 100644 --- a/source/tomopy/sim/project.py +++ b/source/tomopy/sim/project.py @@ -295,7 +295,7 @@ def project2( center = get_center(shape, center) extern.c_project2(objx, objy, center, tomo, theta) - + # NOTE: returns sinogram order with emmission=True if not emission: # convert data to be transmission type @@ -310,7 +310,7 @@ def project2( def project3( - objx, objy, objz, theta, center=None, + objx, objy, objz, theta, center=None, emission=True, pad=True, sinogram_order=False, axis=0, ncore=None, nchunk=None): """ @@ -363,7 +363,7 @@ def project3( center = get_center(shape, center) extern.c_project3(objx, objy, objz, center, tomo, theta, axis) - + # NOTE: returns sinogram order with emmission=True if not emission: # convert data to be transmission type @@ -373,7 +373,7 @@ def project3( tomo = np.swapaxes(tomo, 0, 1) # doesn't copy data # copy data to sharedmem tomo = dtype.as_sharedmem(tomo, copy=True) - + return tomo @@ -409,7 +409,7 @@ def fan_to_para(tomo, dist, geom): ndarray Transformed 3D tomographic data. """ - logger.warning('Not implemented.') + raise NotImplementedError('Not implemented.') def para_to_fan(tomo, dist, geom): @@ -436,7 +436,7 @@ def para_to_fan(tomo, dist, geom): ndarray Transformed 3D tomographic data. """ - logger.warning('Not implemented.') + raise NotImplementedError('Not implemented.') def add_focal_spot_blur(tomo, spotsize): @@ -454,7 +454,7 @@ def add_focal_spot_blur(tomo, spotsize): spotsize : float Focal spot size of circular x-ray source. """ - logger.warning('Not implemented.') + raise NotImplementedError('Not implemented.') def _get_magnification(r1, r2): @@ -499,4 +499,4 @@ def _get_otf(dx, dy, px, py, spotsize): array 2D OTF function. """ - logger.warning('Not implemented.') + raise NotImplementedError('Not implemented.') From d56a998438554a0fc1d7d544d7f0997c6e8196d3 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Sun, 19 Sep 2021 21:02:55 -0500 Subject: [PATCH 03/35] TST: Add atol param to ordered subset tests --- test/test_tomopy/test_recon/test_algorithm.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/test_tomopy/test_recon/test_algorithm.py b/test/test_tomopy/test_recon/test_algorithm.py index dba9ad40e..21da52d69 100644 --- a/test/test_tomopy/test_recon/test_algorithm.py +++ b/test/test_tomopy/test_recon/test_algorithm.py @@ -100,6 +100,7 @@ def test_bart(self): ), read_file('bart.npy'), rtol=1e-2, + atol=1e-5, ) def test_fbp(self): @@ -211,6 +212,7 @@ def test_osem(self): ), read_file('osem.npy'), rtol=1e-2, + atol=1e-5, ) def test_ospml_hybrid(self): @@ -228,6 +230,7 @@ def test_ospml_hybrid(self): ), read_file('ospml_hybrid.npy'), rtol=1e-2, + atol=1e-5, ) def test_ospml_quad(self): @@ -245,6 +248,7 @@ def test_ospml_quad(self): ), read_file('ospml_quad.npy'), rtol=1e-2, + atol=1e-5, ) def test_pml_hybrid(self): From a8f24575a917feed99551bbf3e1b0a3839afabc0 Mon Sep 17 00:00:00 2001 From: Doga Gursoy Date: Wed, 29 Sep 2021 11:47:35 -0500 Subject: [PATCH 04/35] Convert upsample_factor from positional to key argument (#563) This is to overcome the error: TypeError: phase_cross_correlation() takes 2 positional arguments but 3 were given --- source/tomopy/prep/alignment.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/tomopy/prep/alignment.py b/source/tomopy/prep/alignment.py index cc39c67a3..edd0f0ab9 100644 --- a/source/tomopy/prep/alignment.py +++ b/source/tomopy/prep/alignment.py @@ -190,7 +190,7 @@ def align_seq( # Register current projection in sub-pixel precision shift, error, diffphase = phase_cross_correlation( - _prj[m], _sim[m], upsample_factor) + _prj[m], _sim[m], upsample_factor=upsample_factor) err[m] = np.sqrt(shift[0]*shift[0] + shift[1]*shift[1]) sx[m] += shift[0] sy[m] += shift[1] @@ -329,7 +329,7 @@ def align_joint( # Register current projection in sub-pixel precision shift, error, diffphase = phase_cross_correlation( - _prj[m], _sim[m], upsample_factor) + _prj[m], _sim[m], upsample_factor=upsample_factor) err[m] = np.sqrt(shift[0]*shift[0] + shift[1]*shift[1]) sx[m] += shift[0] sy[m] += shift[1] From eb14a80417a74c9999762b0c5d1b3edde6937402 Mon Sep 17 00:00:00 2001 From: "Jonathan R. Madsen" Date: Tue, 12 Oct 2021 20:17:57 -0500 Subject: [PATCH 05/35] Fixed mismatched function prototype vs. implementation --- source/include/gridrec/gridrec.h | 2 +- source/include/recon/recon.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/source/include/gridrec/gridrec.h b/source/include/gridrec/gridrec.h index b589e4bf1..c84c2a13c 100644 --- a/source/include/gridrec/gridrec.h +++ b/source/include/gridrec/gridrec.h @@ -57,7 +57,7 @@ extern "C" DLL void gridrec(const float* data, int dy, int dt, int dx, const float* center, const float* theta, float* recon, int ngridx, int ngridy, - const char fname[16], const float* filter_par); + const char* fname, const float* filter_par); #ifdef __cplusplus } diff --git a/source/include/recon/recon.h b/source/include/recon/recon.h index 68acd4bad..178e93478 100644 --- a/source/include/recon/recon.h +++ b/source/include/recon/recon.h @@ -60,7 +60,7 @@ void DLL void DLL fbp(const float* data, int dy, int dt, int dx, const float* center, const float* theta, - float* recon, int ngridx, int ngridy, const char name[16], const float* filter_par); + float* recon, int ngridx, int ngridy, const char* name, const float* filter_par); void DLL grad(const float* data, int dy, int dt, int dx, const float* center, const float* theta, From a11840905ba4db232bdddfb501da3789b3bc8ade Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Mon, 1 Nov 2021 13:36:07 -0500 Subject: [PATCH 06/35] Adding in median filter for non-finite values in 3D array --- source/tomopy/misc/corr.py | 129 +++++++++++++++++++++++++++++++++++++ 1 file changed, 129 insertions(+) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index bbb13f980..9930395a0 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -63,6 +63,7 @@ import numexpr as ne import concurrent.futures as cf from scipy.signal import medfilt2d +from tqdm import tqdm logger = logging.getLogger(__name__) @@ -267,6 +268,134 @@ def median_filter_cuda(arr, size=3, axis=0): return out +def find_nonfinite_values_prj(projection, slc_idx): + """Allows the user to easily obtain the Z, x, y coordinates of nonfinite + values for a given 2D projection in a 3D array. + + Parameters + ---------- + projection : 2D nd.array + The 2D projection for a given 3D data array. + + slc_idx : int + The projection index inside the main 3D array. + + Returns + ------- + The indices for the current projection used, the x coordinates, and y + coordinates of all non-finite values. + """ + + # Determining where the nonfinite values are + boolean_of_nonfinite = ~np.isfinite(projection) + + # Obtaining there x, y locations + idx_of_nonfinite_prj = np.nonzero(boolean_of_nonfinite) + + # Creating an identical shaped array for the z position in the 3D array the prj comes from + idx_of_nonfinite_data = np.asarray([slc_idx] * len(idx_of_nonfinite_prj[0])) + + return idx_of_nonfinite_data, idx_of_nonfinite_prj[0], idx_of_nonfinite_prj[1] + + +def determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): + """Determine the proper kernel bounds for a given x, y index, and kernel size. + + Parameters + ---------- + x_idx : int + The x index for a given nonfinite value. + + y_idx : int + The y index for a given nonfinite value. + + shape_x : int + The max shape of the 2D projection in the x dimension. + + shape_y : int + The max shape of the 2D projection in the y dimension. + + Returns + ------- + The integer kernel bounds to surround the nonfinite value with. + """ + + # Determining x kernel + x_idx_lower = x_idx - kernel + x_idx_higher = x_idx + kernel + 1 + + if x_idx_lower < 0: + x_idx_lower = 0 + if x_idx_higher > shape_x-1: + x_idx_higher = shape_x-1 + + # Determining y kernel + y_idx_lower = y_idx - kernel + y_idx_higher = y_idx + kernel + 1 + + if y_idx_lower < 0: + y_idx_lower = 0 + if y_idx_higher > shape_y-1: + y_idx_higher = shape_y-1 + + # Forcing values to be integers + integer_values_pre = [x_idx_lower, x_idx_higher, y_idx_lower, y_idx_higher] + integer_values = [int(i) for i in integer_values_pre] + + return integer_values + +def median_filter_nonfinite(data, kernel=1): + """Apply a selective median filter with size kernel to all nonfinite values in the 3D data array. + + Parameters + ---------- + data : 3D nd.array + The 3D array of data with nonfinite values in it. + + kernel : int + The size of the kernel to be used for a local median filter. + + Returns + ------- + The corrected 3D array with all nonfinite values removed based upon the local + median value defined by the kernel size. + """ + + # Creating store arrays for the prj, x, y indicies for bad values + nonfinite_idx_store = [[],[],[]] + + # Iterating throug each projection to save on RAM + for idx, projection in tqdm(enumerate(data), desc='Finding Bad Values', total=len(data)): + nonfinite_idx_store = np.concatenate((nonfinite_idx_store, + find_nonfinite_values(projection, int(idx))), axis=1) + + shape = np.shape(data) + + # Iterating through each bad value and replace it with finite median + for indices in tqdm(zip(nonfinite_idx_store.transpose()), desc='Removing Bad Values'): + + # Turning index values to integers + prj_idx, x_idx, y_idx = [int(ii) for ii in indices[0]] + + # Determining the lower and upper bounds for kernel + x_lower, x_higher, y_lower, y_higher = determine_nonfinite_kernel_idxs(x_idx, + y_idx, + kernel, + shape[1], + shape[2]) + + # Extracting kernel data and fining finite median + kernel_cropped_data = data[prj_idx, x_lower:x_higher, y_lower:y_higher] + kernel_cropped_data = np.asarray(kernel_cropped_data) + median_corrected_data = np.median(kernel_cropped_data[np.isfinite(kernel_cropped_data)]) + + # Replacing bad data with finite median + data[prj_idx, x_idx, y_idx] = median_corrected_data + + return data + + + def sobel_filter(arr, axis=0, ncore=None): """ Apply Sobel filter to 3D array along specified axis. From cbe645e1be1b6c3313bb0d33607e8111f7f8b9a5 Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Mon, 1 Nov 2021 13:40:08 -0500 Subject: [PATCH 07/35] Formatting for PEP8 --- source/tomopy/misc/corr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 9930395a0..ab9e02c2c 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -344,6 +344,7 @@ def determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): return integer_values + def median_filter_nonfinite(data, kernel=1): """Apply a selective median filter with size kernel to all nonfinite values in the 3D data array. @@ -395,7 +396,6 @@ def median_filter_nonfinite(data, kernel=1): return data - def sobel_filter(arr, axis=0, ncore=None): """ Apply Sobel filter to 3D array along specified axis. From bd53fba1df1b1157e34ebcb142febee8a47731ae Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Mon, 1 Nov 2021 13:48:15 -0500 Subject: [PATCH 08/35] Updating PEP8v2 --- source/tomopy/misc/corr.py | 61 ++++++++++++++++++++------------------ 1 file changed, 32 insertions(+), 29 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index ab9e02c2c..f9f6a25ef 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -271,55 +271,56 @@ def median_filter_cuda(arr, size=3, axis=0): def find_nonfinite_values_prj(projection, slc_idx): """Allows the user to easily obtain the Z, x, y coordinates of nonfinite values for a given 2D projection in a 3D array. - + Parameters ---------- projection : 2D nd.array The 2D projection for a given 3D data array. - + slc_idx : int The projection index inside the main 3D array. - + Returns ------- The indices for the current projection used, the x coordinates, and y coordinates of all non-finite values. """ - + # Determining where the nonfinite values are boolean_of_nonfinite = ~np.isfinite(projection) - + # Obtaining there x, y locations idx_of_nonfinite_prj = np.nonzero(boolean_of_nonfinite) - + # Creating an identical shaped array for the z position in the 3D array the prj comes from - idx_of_nonfinite_data = np.asarray([slc_idx] * len(idx_of_nonfinite_prj[0])) - + idx_of_nonfinite_data = np.asarray( + [slc_idx] * len(idx_of_nonfinite_prj[0])) + return idx_of_nonfinite_data, idx_of_nonfinite_prj[0], idx_of_nonfinite_prj[1] def determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): """Determine the proper kernel bounds for a given x, y index, and kernel size. - + Parameters ---------- x_idx : int The x index for a given nonfinite value. - + y_idx : int The y index for a given nonfinite value. - + shape_x : int The max shape of the 2D projection in the x dimension. - + shape_y : int The max shape of the 2D projection in the y dimension. - + Returns ------- The integer kernel bounds to surround the nonfinite value with. """ - + # Determining x kernel x_idx_lower = x_idx - kernel x_idx_higher = x_idx + kernel + 1 @@ -337,25 +338,26 @@ def determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): y_idx_lower = 0 if y_idx_higher > shape_y-1: y_idx_higher = shape_y-1 - + # Forcing values to be integers integer_values_pre = [x_idx_lower, x_idx_higher, y_idx_lower, y_idx_higher] integer_values = [int(i) for i in integer_values_pre] - + return integer_values def median_filter_nonfinite(data, kernel=1): - """Apply a selective median filter with size kernel to all nonfinite values in the 3D data array. - + """Apply a selective median filter with size kernel to all nonfinite + values in the 3D data array. + Parameters ---------- data : 3D nd.array The 3D array of data with nonfinite values in it. - + kernel : int The size of the kernel to be used for a local median filter. - + Returns ------- The corrected 3D array with all nonfinite values removed based upon the local @@ -363,7 +365,7 @@ def median_filter_nonfinite(data, kernel=1): """ # Creating store arrays for the prj, x, y indicies for bad values - nonfinite_idx_store = [[],[],[]] + nonfinite_idx_store = [[], [], []] # Iterating throug each projection to save on RAM for idx, projection in tqdm(enumerate(data), desc='Finding Bad Values', total=len(data)): @@ -374,21 +376,22 @@ def median_filter_nonfinite(data, kernel=1): # Iterating through each bad value and replace it with finite median for indices in tqdm(zip(nonfinite_idx_store.transpose()), desc='Removing Bad Values'): - + # Turning index values to integers prj_idx, x_idx, y_idx = [int(ii) for ii in indices[0]] - + # Determining the lower and upper bounds for kernel x_lower, x_higher, y_lower, y_higher = determine_nonfinite_kernel_idxs(x_idx, - y_idx, - kernel, - shape[1], - shape[2]) - + y_idx, + kernel, + shape[1], + shape[2]) + # Extracting kernel data and fining finite median kernel_cropped_data = data[prj_idx, x_lower:x_higher, y_lower:y_higher] kernel_cropped_data = np.asarray(kernel_cropped_data) - median_corrected_data = np.median(kernel_cropped_data[np.isfinite(kernel_cropped_data)]) + median_corrected_data = np.median( + kernel_cropped_data[np.isfinite(kernel_cropped_data)]) # Replacing bad data with finite median data[prj_idx, x_idx, y_idx] = median_corrected_data From 7e8de9322bcb672bcba7f52e1c12f09ba9300f40 Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Mon, 1 Nov 2021 17:43:51 -0500 Subject: [PATCH 09/35] Update source/tomopy/misc/corr.py Adding _ to the find_nonfinite_values_prj function since it is not intended to be used by the User Co-authored-by: Daniel Ching --- source/tomopy/misc/corr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index f9f6a25ef..bd77822ce 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -268,7 +268,7 @@ def median_filter_cuda(arr, size=3, axis=0): return out -def find_nonfinite_values_prj(projection, slc_idx): +def _find_nonfinite_values_prj(projection, slc_idx): """Allows the user to easily obtain the Z, x, y coordinates of nonfinite values for a given 2D projection in a 3D array. From 79be52b09af2f6711dd528b3348d3ffccc9506aa Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Mon, 1 Nov 2021 17:44:20 -0500 Subject: [PATCH 10/35] Update source/tomopy/misc/corr.py Adding _ to the determine_nonfinite_kernel_idxs function since it is not intended to be used by the User Co-authored-by: Daniel Ching --- source/tomopy/misc/corr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index bd77822ce..d3b066187 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -299,7 +299,7 @@ def _find_nonfinite_values_prj(projection, slc_idx): return idx_of_nonfinite_data, idx_of_nonfinite_prj[0], idx_of_nonfinite_prj[1] -def determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): +def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): """Determine the proper kernel bounds for a given x, y index, and kernel size. Parameters From 040be55f8aeb43155a9621378b4f71504a27d60b Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Mon, 1 Nov 2021 18:01:15 -0500 Subject: [PATCH 11/35] Updating median_filter_nonfinite to not require tqdm --- source/tomopy/misc/corr.py | 74 ++++++++++++++++++++++++-------------- 1 file changed, 48 insertions(+), 26 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index d3b066187..656a639cc 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -63,7 +63,6 @@ import numexpr as ne import concurrent.futures as cf from scipy.signal import medfilt2d -from tqdm import tqdm logger = logging.getLogger(__name__) @@ -268,22 +267,27 @@ def median_filter_cuda(arr, size=3, axis=0): return out +import numpy as np +from tqdm import tqdm + + def _find_nonfinite_values_prj(projection, slc_idx): - """Allows the user to easily obtain the Z, x, y coordinates of nonfinite - values for a given 2D projection in a 3D array. + """ + Allows the user to easily obtain the Z, x, y coordinates of nonfinite + values for a given 2D projection in a 3D array. Parameters ---------- projection : 2D nd.array The 2D projection for a given 3D data array. - slc_idx : int The projection index inside the main 3D array. Returns ------- - The indices for the current projection used, the x coordinates, and y - coordinates of all non-finite values. + ndarray + The indices for the current projection used, the x coordinates, and y + coordinates of all non-finite values. """ # Determining where the nonfinite values are @@ -300,25 +304,24 @@ def _find_nonfinite_values_prj(projection, slc_idx): def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): - """Determine the proper kernel bounds for a given x, y index, and kernel size. + """ + Determine the proper kernel bounds for a given x, y index, and kernel size. Parameters ---------- x_idx : int The x index for a given nonfinite value. - y_idx : int The y index for a given nonfinite value. - shape_x : int The max shape of the 2D projection in the x dimension. - shape_y : int The max shape of the 2D projection in the y dimension. Returns ------- - The integer kernel bounds to surround the nonfinite value with. + ndarray + The integer kernel bounds to surround the nonfinite value with. """ # Determining x kernel @@ -346,56 +349,75 @@ def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): return integer_values -def median_filter_nonfinite(data, kernel=1): - """Apply a selective median filter with size kernel to all nonfinite - values in the 3D data array. +def median_filter_nonfinite(data, kernel=1, pbar=None): + """ + Apply a selective 2D median filter, slice by slice, for a 3D data array + with filter_size=kernel to all nonfinite values in the 3D data array. Parameters ---------- data : 3D nd.array The 3D array of data with nonfinite values in it. - kernel : int The size of the kernel to be used for a local median filter. + pbar : tqdm instance + A tqdm instance allowing the User to display a useful progressbar. Returns ------- - The corrected 3D array with all nonfinite values removed based upon the local - median value defined by the kernel size. + ndarray + The corrected 3D array with all nonfinite values removed based upon the local + median value defined by the kernel size. """ # Creating store arrays for the prj, x, y indicies for bad values nonfinite_idx_store = [[], [], []] + # Allowing the User to use a tqdm progressbar if desired + if pbar is not None: + pbar.total = total = len(data) + pbar.set_description('Finding Bad Values') + # Iterating throug each projection to save on RAM - for idx, projection in tqdm(enumerate(data), desc='Finding Bad Values', total=len(data)): + for idx, projection in enumerate(data): nonfinite_idx_store = np.concatenate((nonfinite_idx_store, - find_nonfinite_values(projection, int(idx))), axis=1) + _find_nonfinite_values_prj(projection, + int(idx))), axis=1) + if pbar is not None: + pbar.update(1) shape = np.shape(data) + # Allowing the User to use a tqdm progressbar if desired + if pbar is not None: + pbar.reset() + pbar.set_description('Removing Bad Values') + pbar.total = len(nonfinite_idx_store.transpose()) + # Iterating through each bad value and replace it with finite median - for indices in tqdm(zip(nonfinite_idx_store.transpose()), desc='Removing Bad Values'): + for indices in zip(nonfinite_idx_store.transpose()): # Turning index values to integers prj_idx, x_idx, y_idx = [int(ii) for ii in indices[0]] # Determining the lower and upper bounds for kernel - x_lower, x_higher, y_lower, y_higher = determine_nonfinite_kernel_idxs(x_idx, - y_idx, - kernel, - shape[1], - shape[2]) + x_lower, x_higher, y_lower, y_higher = _determine_nonfinite_kernel_idxs(x_idx, + y_idx, + kernel, + shape[1], + shape[2]) # Extracting kernel data and fining finite median kernel_cropped_data = data[prj_idx, x_lower:x_higher, y_lower:y_higher] - kernel_cropped_data = np.asarray(kernel_cropped_data) median_corrected_data = np.median( kernel_cropped_data[np.isfinite(kernel_cropped_data)]) # Replacing bad data with finite median data[prj_idx, x_idx, y_idx] = median_corrected_data + if pbar is not None: + pbar.update(1) + return data From e5ee919bc29efc78b24e3bc09455ed3fe7577967 Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Mon, 1 Nov 2021 18:08:08 -0500 Subject: [PATCH 12/35] Accidentally left tqdm inside imports - removing with commit --- source/tomopy/misc/corr.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 656a639cc..aaa496b85 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -267,10 +267,6 @@ def median_filter_cuda(arr, size=3, axis=0): return out -import numpy as np -from tqdm import tqdm - - def _find_nonfinite_values_prj(projection, slc_idx): """ Allows the user to easily obtain the Z, x, y coordinates of nonfinite From 832a4b2d535773172474136cd23884e34d59806e Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Mon, 1 Nov 2021 18:14:19 -0500 Subject: [PATCH 13/35] Adding function to __all__ and adding authors name --- source/tomopy/misc/corr.py | 1 + 1 file changed, 1 insertion(+) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index aaa496b85..b9fbf1dde 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -76,6 +76,7 @@ 'gaussian_filter', 'median_filter', 'median_filter_cuda', + 'median_filter_nonfinite', 'sobel_filter', 'remove_nan', 'remove_neg', From d7d3e3ef6ced3da49ac6c25bf98cf9355ea5d35d Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Mon, 1 Nov 2021 19:18:07 -0500 Subject: [PATCH 14/35] Adding a CallbackClass to median_filter_nonfinite and adding testing for function --- source/tomopy/misc/corr.py | 48 +++++++++++++++---------- test/test_tomopy/test_misc/test_corr.py | 19 ++++++++-- 2 files changed, 47 insertions(+), 20 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index b9fbf1dde..420ba82ea 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -67,7 +67,7 @@ logger = logging.getLogger(__name__) -__author__ = "Doga Gursoy" +__author__ = "Doga Gursoy, William Judge" __credits__ = "Mark Rivers, Xianghui Xiao" __copyright__ = "Copyright (c) 2015, UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' @@ -346,7 +346,7 @@ def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): return integer_values -def median_filter_nonfinite(data, kernel=1, pbar=None): +def median_filter_nonfinite(data, kernel=1, callback=None): """ Apply a selective 2D median filter, slice by slice, for a 3D data array with filter_size=kernel to all nonfinite values in the 3D data array. @@ -357,8 +357,10 @@ def median_filter_nonfinite(data, kernel=1, pbar=None): The 3D array of data with nonfinite values in it. kernel : int The size of the kernel to be used for a local median filter. - pbar : tqdm instance - A tqdm instance allowing the User to display a useful progressbar. + callback : func + A callback function to be run for each loop iteration. + Use self.iterator_len and self.iterator_desc to obtain the + length and description for each loop Returns ------- @@ -367,29 +369,39 @@ def median_filter_nonfinite(data, kernel=1, pbar=None): median value defined by the kernel size. """ + # Defining a callback class for a User defined progressbar + class CallbackClass(): + + def __init__(self, iterator_len, iterator_desc): + self.iterator_len = iterator_len + self.iterator_desc = iterator_desc + # Creating store arrays for the prj, x, y indicies for bad values nonfinite_idx_store = [[], [], []] - # Allowing the User to use a tqdm progressbar if desired - if pbar is not None: - pbar.total = total = len(data) - pbar.set_description('Finding Bad Values') + # Allowing the User to use a progressbar if desired + if callback is not None: + data_len_loop1 = len(data) + callback_loop1 = CallbackClass( + iterator_len=data_len_loop1, iterator_desc='Finding Bad Values') + callback_loop1.callback = callback # Iterating throug each projection to save on RAM for idx, projection in enumerate(data): nonfinite_idx_store = np.concatenate((nonfinite_idx_store, _find_nonfinite_values_prj(projection, - int(idx))), axis=1) - if pbar is not None: - pbar.update(1) + int(idx))), axis=1) + if callback is not None: + callback_loop1.callback(callback_loop1) shape = np.shape(data) - # Allowing the User to use a tqdm progressbar if desired - if pbar is not None: - pbar.reset() - pbar.set_description('Removing Bad Values') - pbar.total = len(nonfinite_idx_store.transpose()) + # Allowing the User to use a progressbar if desired + if callback is not None: + data_len_loop2 = len(nonfinite_idx_store.transpose()) + len(data) + callback_loop2 = CallbackClass(iterator_len=data_len_loop2, + iterator_desc='Removing Bad Values') + callback_loop2.callback = callback # Iterating through each bad value and replace it with finite median for indices in zip(nonfinite_idx_store.transpose()): @@ -412,8 +424,8 @@ def median_filter_nonfinite(data, kernel=1, pbar=None): # Replacing bad data with finite median data[prj_idx, x_idx, y_idx] = median_corrected_data - if pbar is not None: - pbar.update(1) + if callback is not None: + callback_loop2.callback(callback_loop2) return data diff --git a/test/test_tomopy/test_misc/test_corr.py b/test/test_tomopy/test_misc/test_corr.py index 53f2df45f..72172e1ed 100644 --- a/test/test_tomopy/test_misc/test_corr.py +++ b/test/test_tomopy/test_misc/test_corr.py @@ -49,12 +49,12 @@ unicode_literals) import unittest -from tomopy.misc.corr import gaussian_filter, median_filter, remove_neg, remove_nan, remove_outlier, circ_mask +from tomopy.misc.corr import gaussian_filter, median_filter, median_filter_nonfinite, remove_neg, remove_nan, remove_outlier, circ_mask from ..util import read_file, loop_dim import numpy as np from numpy.testing import assert_allclose -__author__ = "Doga Gursoy" +__author__ = "Doga Gursoy, William Judge" __copyright__ = "Copyright (c) 2015, UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' @@ -66,6 +66,21 @@ def test_gaussian_filter(self): def test_median_filter(self): loop_dim(median_filter, read_file('cube.npy')) + def test_median_filter_nonfinite(): + data_org = np.ones(shape=(100, 100, 100)) + + for i in range(50): + x = np.random.randint(0, 100) + y = np.random.randint(0, 100) + z = np.random.randint(0, 100) + data_org[z, x, y] = np.inf + + data_post_corr = median_filter_nonfinite( + data_org.copy(), kernel=1, callback=None) + + assert np.max(data_org) == np.inf + assert np.max(data_post_corr) == 1.0 + def test_remove_neg(self): assert_allclose( remove_neg( From 13578a863f7acc07d62fb2fb234a1f8b8b7a8dbe Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Mon, 1 Nov 2021 19:20:17 -0500 Subject: [PATCH 15/35] Fixing test formatting for median_filter_nonfinite() --- test/test_tomopy/test_misc/test_corr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_tomopy/test_misc/test_corr.py b/test/test_tomopy/test_misc/test_corr.py index 72172e1ed..b74bf84e9 100644 --- a/test/test_tomopy/test_misc/test_corr.py +++ b/test/test_tomopy/test_misc/test_corr.py @@ -66,7 +66,7 @@ def test_gaussian_filter(self): def test_median_filter(self): loop_dim(median_filter, read_file('cube.npy')) - def test_median_filter_nonfinite(): + def test_median_filter_nonfinite(self): data_org = np.ones(shape=(100, 100, 100)) for i in range(50): From 6f79b9e261c78bc4fd529ecd728b6b7ea0e79130 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Tue, 2 Nov 2021 12:08:07 -0500 Subject: [PATCH 16/35] DOC: use consistent type descriptions --- source/tomopy/misc/corr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 420ba82ea..5857ba647 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -353,7 +353,7 @@ def median_filter_nonfinite(data, kernel=1, callback=None): Parameters ---------- - data : 3D nd.array + data : ndarray The 3D array of data with nonfinite values in it. kernel : int The size of the kernel to be used for a local median filter. From b8698337a58ebca447c91808e161f96f15fedbaa Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Tue, 2 Nov 2021 13:13:13 -0500 Subject: [PATCH 17/35] Use callback funtion instead of CallbackClass for median_filter_nonfinite --- source/tomopy/misc/corr.py | 77 +++++++++++++++----------------------- 1 file changed, 31 insertions(+), 46 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 420ba82ea..cb43bf0e7 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -275,7 +275,7 @@ def _find_nonfinite_values_prj(projection, slc_idx): Parameters ---------- - projection : 2D nd.array + projection : nd.array The 2D projection for a given 3D data array. slc_idx : int The projection index inside the main 3D array. @@ -353,7 +353,7 @@ def median_filter_nonfinite(data, kernel=1, callback=None): Parameters ---------- - data : 3D nd.array + data : nd.array The 3D array of data with nonfinite values in it. kernel : int The size of the kernel to be used for a local median filter. @@ -369,63 +369,48 @@ def median_filter_nonfinite(data, kernel=1, callback=None): median value defined by the kernel size. """ - # Defining a callback class for a User defined progressbar - class CallbackClass(): - - def __init__(self, iterator_len, iterator_desc): - self.iterator_len = iterator_len - self.iterator_desc = iterator_desc + # Defining a callback function if None is provided + if callback is None: + def callback(total, description, unit): + pass # Creating store arrays for the prj, x, y indicies for bad values nonfinite_idx_store = [[], [], []] - # Allowing the User to use a progressbar if desired - if callback is not None: - data_len_loop1 = len(data) - callback_loop1 = CallbackClass( - iterator_len=data_len_loop1, iterator_desc='Finding Bad Values') - callback_loop1.callback = callback + # Obtaining information for user defined progressbar + shape = np.shape(data) + pbar_total = len(data) + pbar_desc = 'Nonfinite Median Filter' + pbar_unit = ' Prjs' # Iterating throug each projection to save on RAM for idx, projection in enumerate(data): - nonfinite_idx_store = np.concatenate((nonfinite_idx_store, - _find_nonfinite_values_prj(projection, - int(idx))), axis=1) - if callback is not None: - callback_loop1.callback(callback_loop1) - - shape = np.shape(data) - - # Allowing the User to use a progressbar if desired - if callback is not None: - data_len_loop2 = len(nonfinite_idx_store.transpose()) + len(data) - callback_loop2 = CallbackClass(iterator_len=data_len_loop2, - iterator_desc='Removing Bad Values') - callback_loop2.callback = callback + nonfinite_idx = _find_nonfinite_values_prj(projection, int(idx)) + nonfinite_idx = np.asarray(nonfinite_idx) - # Iterating through each bad value and replace it with finite median - for indices in zip(nonfinite_idx_store.transpose()): + # Iterating through each bad value and replace it with finite median + for indices in zip(nonfinite_idx.transpose()): - # Turning index values to integers - prj_idx, x_idx, y_idx = [int(ii) for ii in indices[0]] + # Turning index values to integers + prj_idx, x_idx, y_idx = [int(ii) for ii in indices[0]] - # Determining the lower and upper bounds for kernel - x_lower, x_higher, y_lower, y_higher = _determine_nonfinite_kernel_idxs(x_idx, - y_idx, - kernel, - shape[1], - shape[2]) + # Determining the lower and upper bounds for kernel + x_lower, x_higher, y_lower, y_higher = _determine_nonfinite_kernel_idxs(x_idx, + y_idx, + kernel, + shape[1], + shape[2]) - # Extracting kernel data and fining finite median - kernel_cropped_data = data[prj_idx, x_lower:x_higher, y_lower:y_higher] - median_corrected_data = np.median( - kernel_cropped_data[np.isfinite(kernel_cropped_data)]) + # Extracting kernel data and fining finite median + kernel_cropped_data = data[prj_idx, + x_lower:x_higher, y_lower:y_higher] + median_corrected_data = np.median( + kernel_cropped_data[np.isfinite(kernel_cropped_data)]) - # Replacing bad data with finite median - data[prj_idx, x_idx, y_idx] = median_corrected_data + # Replacing bad data with finite median + data[prj_idx, x_idx, y_idx] = median_corrected_data - if callback is not None: - callback_loop2.callback(callback_loop2) + callback(pbar_total, pbar_desc, pbar_unit) return data From aa3a4c00f992690ee70e9ed722737007d3777f71 Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Tue, 2 Nov 2021 13:19:57 -0500 Subject: [PATCH 18/35] Correcting the median_filter_nonfinite function description --- source/tomopy/misc/corr.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index c41cbc4f4..1779da8fc 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -359,8 +359,10 @@ def median_filter_nonfinite(data, kernel=1, callback=None): The size of the kernel to be used for a local median filter. callback : func A callback function to be run for each loop iteration. - Use self.iterator_len and self.iterator_desc to obtain the - length and description for each loop + Define total, description, and unit variables inside the function. + Where total is the total length of the progress bar, + description is the description of the progress bar, + and unit is the iteration unit of the progress bar. Returns ------- From 70cf4e86895cbb23e0b49a2cd3c8b2e3b197ea99 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Tue, 2 Nov 2021 15:47:24 -0500 Subject: [PATCH 19/35] TST: Test all sorts of nonfinite values --- test/test_tomopy/test_misc/test_corr.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/test_tomopy/test_misc/test_corr.py b/test/test_tomopy/test_misc/test_corr.py index b74bf84e9..3a788b7bc 100644 --- a/test/test_tomopy/test_misc/test_corr.py +++ b/test/test_tomopy/test_misc/test_corr.py @@ -74,6 +74,14 @@ def test_median_filter_nonfinite(self): y = np.random.randint(0, 100) z = np.random.randint(0, 100) data_org[z, x, y] = np.inf + x = np.random.randint(0, 100) + y = np.random.randint(0, 100) + z = np.random.randint(0, 100) + data_org[z, x, y] = -np.inf + x = np.random.randint(0, 100) + y = np.random.randint(0, 100) + z = np.random.randint(0, 100) + data_org[z, x, y] = np.nan data_post_corr = median_filter_nonfinite( data_org.copy(), kernel=1, callback=None) From afa195d6d8159d818e8fe8cda49e672f385a0935 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Tue, 2 Nov 2021 15:48:23 -0500 Subject: [PATCH 20/35] DOC: Refactor docstrings for median filter nonfinite --- source/tomopy/misc/corr.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 1779da8fc..2177091ca 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -348,8 +348,10 @@ def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): def median_filter_nonfinite(data, kernel=1, callback=None): """ - Apply a selective 2D median filter, slice by slice, for a 3D data array - with filter_size=kernel to all nonfinite values in the 3D data array. + Remove nonfinite values from a 3D array using an in-place 2D median filter. + + The 2D selective median filter is applied along the last two axes of + the array. Parameters ---------- @@ -357,12 +359,11 @@ def median_filter_nonfinite(data, kernel=1, callback=None): The 3D array of data with nonfinite values in it. kernel : int The size of the kernel to be used for a local median filter. - callback : func - A callback function to be run for each loop iteration. - Define total, description, and unit variables inside the function. - Where total is the total length of the progress bar, - description is the description of the progress bar, - and unit is the iteration unit of the progress bar. + callback : func(total, description, unit) + A function called after every internal loop iteration. + total is number of loop iterations. + description is 'Nonfinite median filter'. + unit is ' prjs'. Returns ------- From f33ebac6e20b666f9253787e3ff9d13b61be328f Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Tue, 2 Nov 2021 16:09:01 -0500 Subject: [PATCH 21/35] Updating median_filter_nonfinite test --- test/test_tomopy/test_misc/test_corr.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test/test_tomopy/test_misc/test_corr.py b/test/test_tomopy/test_misc/test_corr.py index 3a788b7bc..4f7ee7a7b 100644 --- a/test/test_tomopy/test_misc/test_corr.py +++ b/test/test_tomopy/test_misc/test_corr.py @@ -66,7 +66,7 @@ def test_gaussian_filter(self): def test_median_filter(self): loop_dim(median_filter, read_file('cube.npy')) - def test_median_filter_nonfinite(self): + def test_median_filter_nonfinite(): data_org = np.ones(shape=(100, 100, 100)) for i in range(50): @@ -85,8 +85,10 @@ def test_median_filter_nonfinite(self): data_post_corr = median_filter_nonfinite( data_org.copy(), kernel=1, callback=None) - - assert np.max(data_org) == np.inf + + assert np.nanmax(data_org) == np.inf + assert np.nanmin(data_org) == -np.inf + assert np.isnan(np.max(data_org)) assert np.max(data_post_corr) == 1.0 def test_remove_neg(self): From 2383519bc5164fba18dd8f4a171a04a5240fb3ad Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Tue, 2 Nov 2021 16:12:05 -0500 Subject: [PATCH 22/35] REF: Parameterize filter size consistently with other functions It is convention to use filter size as the width of the filter not the radius. --- source/tomopy/misc/corr.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 2177091ca..1d251f8d3 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -346,10 +346,10 @@ def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): return integer_values -def median_filter_nonfinite(data, kernel=1, callback=None): +def median_filter_nonfinite(data, size=3, callback=None): """ Remove nonfinite values from a 3D array using an in-place 2D median filter. - + The 2D selective median filter is applied along the last two axes of the array. @@ -357,8 +357,8 @@ def median_filter_nonfinite(data, kernel=1, callback=None): ---------- data : ndarray The 3D array of data with nonfinite values in it. - kernel : int - The size of the kernel to be used for a local median filter. + size : int, optional + The size of the filter. callback : func(total, description, unit) A function called after every internal loop iteration. total is number of loop iterations. @@ -400,7 +400,7 @@ def callback(total, description, unit): # Determining the lower and upper bounds for kernel x_lower, x_higher, y_lower, y_higher = _determine_nonfinite_kernel_idxs(x_idx, y_idx, - kernel, + size//2, shape[1], shape[2]) From 2e8b430c4d02b47fb5095678320565d14ff6d25c Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Tue, 2 Nov 2021 16:14:06 -0500 Subject: [PATCH 23/35] TST: Test that all values are correct instead of just one --- test/test_tomopy/test_misc/test_corr.py | 43 ++++++++++++++----------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/test/test_tomopy/test_misc/test_corr.py b/test/test_tomopy/test_misc/test_corr.py index 3a788b7bc..cd2f2bf02 100644 --- a/test/test_tomopy/test_misc/test_corr.py +++ b/test/test_tomopy/test_misc/test_corr.py @@ -45,14 +45,21 @@ # POSSIBILITY OF SUCH DAMAGE. # # ######################################################################### -from __future__ import (absolute_import, division, print_function, - unicode_literals) - import unittest -from tomopy.misc.corr import gaussian_filter, median_filter, median_filter_nonfinite, remove_neg, remove_nan, remove_outlier, circ_mask -from ..util import read_file, loop_dim + import numpy as np from numpy.testing import assert_allclose +from tomopy.misc.corr import ( + gaussian_filter, + median_filter, + median_filter_nonfinite, + remove_neg, + remove_nan, + remove_outlier, + circ_mask, +) + +from ..util import read_file, loop_dim __author__ = "Doga Gursoy, William Judge" __copyright__ = "Copyright (c) 2015, UChicago Argonne, LLC." @@ -67,8 +74,8 @@ def test_median_filter(self): loop_dim(median_filter, read_file('cube.npy')) def test_median_filter_nonfinite(self): + # Add some random non-finite values to an array of all ones data_org = np.ones(shape=(100, 100, 100)) - for i in range(50): x = np.random.randint(0, 100) y = np.random.randint(0, 100) @@ -84,27 +91,25 @@ def test_median_filter_nonfinite(self): data_org[z, x, y] = np.nan data_post_corr = median_filter_nonfinite( - data_org.copy(), kernel=1, callback=None) - - assert np.max(data_org) == np.inf - assert np.max(data_post_corr) == 1.0 + data_org.copy(), + size=5, + callback=None, + ) + # All the post filtering values should be 1 because all of the finite + # values are 1. + assert np.all(data_post_corr == 1.0) def test_remove_neg(self): - assert_allclose( - remove_neg( - [-2, -1, 0, 1, 2]), [0, 0, 0, 1, 2]) + assert_allclose(remove_neg([-2, -1, 0, 1, 2]), [0, 0, 0, 1, 2]) def test_remove_nan(self): - assert_allclose( - remove_nan( - [np.nan, 1.5, 2, np.nan, 1]), [0, 1.5, 2, 0, 1]) + assert_allclose(remove_nan([np.nan, 1.5, 2, np.nan, 1]), + [0, 1.5, 2, 0, 1]) def test_remove_outlier(self): proj = read_file('proj.npy') proj[8][4][6] = 20 - assert_allclose( - remove_outlier(proj, dif=10), - read_file('proj.npy')) + assert_allclose(remove_outlier(proj, dif=10), read_file('proj.npy')) def test_circ_mask(self): loop_dim(circ_mask, read_file('obj.npy')) From 74aa99f48d6aacec2f412cea08d09453133d38c4 Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Tue, 2 Nov 2021 16:23:38 -0500 Subject: [PATCH 24/35] Adding self argument to test --- test/test_tomopy/test_misc/test_corr.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_tomopy/test_misc/test_corr.py b/test/test_tomopy/test_misc/test_corr.py index 4f7ee7a7b..036e36774 100644 --- a/test/test_tomopy/test_misc/test_corr.py +++ b/test/test_tomopy/test_misc/test_corr.py @@ -66,7 +66,7 @@ def test_gaussian_filter(self): def test_median_filter(self): loop_dim(median_filter, read_file('cube.npy')) - def test_median_filter_nonfinite(): + def test_median_filter_nonfinite(self): data_org = np.ones(shape=(100, 100, 100)) for i in range(50): @@ -85,7 +85,7 @@ def test_median_filter_nonfinite(): data_post_corr = median_filter_nonfinite( data_org.copy(), kernel=1, callback=None) - + assert np.nanmax(data_org) == np.inf assert np.nanmin(data_org) == -np.inf assert np.isnan(np.max(data_org)) From af94d9d9f8c76f555c2f34205b073abbec087e76 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Tue, 2 Nov 2021 16:42:40 -0500 Subject: [PATCH 25/35] REF: Remove unnecessary variable declarations --- source/tomopy/misc/corr.py | 81 +++++++------------------------------- 1 file changed, 15 insertions(+), 66 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 1d251f8d3..abd74dff9 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -268,38 +268,6 @@ def median_filter_cuda(arr, size=3, axis=0): return out -def _find_nonfinite_values_prj(projection, slc_idx): - """ - Allows the user to easily obtain the Z, x, y coordinates of nonfinite - values for a given 2D projection in a 3D array. - - Parameters - ---------- - projection : ndarray - The 2D projection for a given 3D data array. - slc_idx : int - The projection index inside the main 3D array. - - Returns - ------- - ndarray - The indices for the current projection used, the x coordinates, and y - coordinates of all non-finite values. - """ - - # Determining where the nonfinite values are - boolean_of_nonfinite = ~np.isfinite(projection) - - # Obtaining there x, y locations - idx_of_nonfinite_prj = np.nonzero(boolean_of_nonfinite) - - # Creating an identical shaped array for the z position in the 3D array the prj comes from - idx_of_nonfinite_data = np.asarray( - [slc_idx] * len(idx_of_nonfinite_prj[0])) - - return idx_of_nonfinite_data, idx_of_nonfinite_prj[0], idx_of_nonfinite_prj[1] - - def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): """ Determine the proper kernel bounds for a given x, y index, and kernel size. @@ -317,7 +285,7 @@ def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): Returns ------- - ndarray + x_lower, x_higher, y_lower, y_higher : int The integer kernel bounds to surround the nonfinite value with. """ @@ -339,11 +307,7 @@ def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): if y_idx_higher > shape_y-1: y_idx_higher = shape_y-1 - # Forcing values to be integers - integer_values_pre = [x_idx_lower, x_idx_higher, y_idx_lower, y_idx_higher] - integer_values = [int(i) for i in integer_values_pre] - - return integer_values + return x_idx_lower, x_idx_higher, y_idx_lower, y_idx_higher def median_filter_nonfinite(data, size=3, callback=None): @@ -371,49 +335,34 @@ def median_filter_nonfinite(data, size=3, callback=None): The corrected 3D array with all nonfinite values removed based upon the local median value defined by the kernel size. """ - # Defining a callback function if None is provided if callback is None: def callback(total, description, unit): pass - # Creating store arrays for the prj, x, y indicies for bad values - nonfinite_idx_store = [[], [], []] - - # Obtaining information for user defined progressbar - shape = np.shape(data) - pbar_total = len(data) - pbar_desc = 'Nonfinite Median Filter' - pbar_unit = ' Prjs' - # Iterating throug each projection to save on RAM - for idx, projection in enumerate(data): - nonfinite_idx = _find_nonfinite_values_prj(projection, int(idx)) - nonfinite_idx = np.asarray(nonfinite_idx) + for projection in data: + nonfinite_idx = np.nonzero(~np.isfinite(projection)) # Iterating through each bad value and replace it with finite median - for indices in zip(nonfinite_idx.transpose()): - - # Turning index values to integers - prj_idx, x_idx, y_idx = [int(ii) for ii in indices[0]] + for x_idx, y_idx in zip(*nonfinite_idx): # Determining the lower and upper bounds for kernel x_lower, x_higher, y_lower, y_higher = _determine_nonfinite_kernel_idxs(x_idx, y_idx, size//2, - shape[1], - shape[2]) + data.shape[1], + data.shape[2]) # Extracting kernel data and fining finite median - kernel_cropped_data = data[prj_idx, - x_lower:x_higher, y_lower:y_higher] + kernel_cropped_data = projection[x_lower:x_higher, y_lower:y_higher] median_corrected_data = np.median( kernel_cropped_data[np.isfinite(kernel_cropped_data)]) # Replacing bad data with finite median - data[prj_idx, x_idx, y_idx] = median_corrected_data + projection[x_idx, y_idx] = median_corrected_data - callback(pbar_total, pbar_desc, pbar_unit) + callback(data.shape[0], 'Nonfinite median filter', ' prjs') return data @@ -830,8 +779,8 @@ def enhance_projs_aps_1id(imgstack, median_ks=5, ncore=None): """ Enhance the projection images with weak contrast collected at APS 1ID - This filter uses a median fileter (will be switched to enhanced recursive - median fileter, ERMF, in the future) for denoising, and a histogram + This filter uses a median fileter (will be switched to enhanced recursive + median fileter, ERMF, in the future) for denoising, and a histogram equalization for dynamic range adjustment to bring out the details. Parameters @@ -865,7 +814,7 @@ def enhance_projs_aps_1id(imgstack, median_ks=5, ncore=None): def _enhance_img(img, median_ks, normalized=True): """ - Enhance the projection image from aps 1ID to counter its weak contrast + Enhance the projection image from aps 1ID to counter its weak contrast nature Parameters @@ -878,7 +827,7 @@ def _enhance_img(img, median_ks, normalized=True): specify whether the enhanced image is normalized between 0 and 1, default is True - Returns + Returns ------- ndarray enhanced projection image @@ -901,7 +850,7 @@ def _calc_histequal_wgt(img): Returns ------- ndarray - histogram euqalization weights (0-1) in the same shape as original + histogram euqalization weights (0-1) in the same shape as original image """ return (np.sort(img.flatten()).searchsorted(img) + 1)/np.prod(img.shape) From 7140bd7c07a1da6db4debb87d3a9ca8bba0c637d Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Tue, 2 Nov 2021 16:45:38 -0500 Subject: [PATCH 26/35] BUG: it is valid for a range end to equal the shape --- source/tomopy/misc/corr.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index abd74dff9..e44f50d70 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -295,8 +295,8 @@ def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): if x_idx_lower < 0: x_idx_lower = 0 - if x_idx_higher > shape_x-1: - x_idx_higher = shape_x-1 + if x_idx_higher > shape_x: + x_idx_higher = shape_x # Determining y kernel y_idx_lower = y_idx - kernel @@ -304,8 +304,8 @@ def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): if y_idx_lower < 0: y_idx_lower = 0 - if y_idx_higher > shape_y-1: - y_idx_higher = shape_y-1 + if y_idx_higher > shape_y: + y_idx_higher = shape_y return x_idx_lower, x_idx_higher, y_idx_lower, y_idx_higher From 102d60eaaa7fe5844fc64e20ae640d7ae5bd5961 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Tue, 2 Nov 2021 16:49:38 -0500 Subject: [PATCH 27/35] STY: Format whitespace --- source/tomopy/misc/corr.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index e44f50d70..996c6e772 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -337,6 +337,7 @@ def median_filter_nonfinite(data, size=3, callback=None): """ # Defining a callback function if None is provided if callback is None: + def callback(total, description, unit): pass @@ -348,14 +349,22 @@ def callback(total, description, unit): for x_idx, y_idx in zip(*nonfinite_idx): # Determining the lower and upper bounds for kernel - x_lower, x_higher, y_lower, y_higher = _determine_nonfinite_kernel_idxs(x_idx, - y_idx, - size//2, - data.shape[1], - data.shape[2]) + ( + x_lower, + x_higher, + y_lower, + y_higher, + ) = _determine_nonfinite_kernel_idxs( + x_idx, + y_idx, + size // 2, + data.shape[1], + data.shape[2], + ) # Extracting kernel data and fining finite median - kernel_cropped_data = projection[x_lower:x_higher, y_lower:y_higher] + kernel_cropped_data = projection[x_lower:x_higher, + y_lower:y_higher] median_corrected_data = np.median( kernel_cropped_data[np.isfinite(kernel_cropped_data)]) From 05a887c0e2b4ec307410860c2df4fb5a7a154c02 Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Thu, 4 Nov 2021 11:53:54 -0500 Subject: [PATCH 28/35] Replacing _determine_nonfinite_kernel_idxs with min()/max() --- source/tomopy/misc/corr.py | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 996c6e772..3f1eaaf06 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -349,18 +349,10 @@ def callback(total, description, unit): for x_idx, y_idx in zip(*nonfinite_idx): # Determining the lower and upper bounds for kernel - ( - x_lower, - x_higher, - y_lower, - y_higher, - ) = _determine_nonfinite_kernel_idxs( - x_idx, - y_idx, - size // 2, - data.shape[1], - data.shape[2], - ) + x_lower = max(0, x_idx - (size // 2)) + x_higher = min(data.shape[1], x_idx + (size // 2) + 1) + y_lower = max(0, y_idx - (size // 2)) + y_higher = min(data.shape[2], y_idx + (size // 2) + 1) # Extracting kernel data and fining finite median kernel_cropped_data = projection[x_lower:x_higher, From 2deaae488a6e367dc7944db03081100324acadaf Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Thu, 4 Nov 2021 13:43:05 -0500 Subject: [PATCH 29/35] Adding kernel ValueError to median_filter_nonfinite() + coorresponding tests + removing unused functions --- source/tomopy/misc/corr.py | 56 ++++++------------------- test/test_tomopy/test_misc/test_corr.py | 39 ++++++++++------- 2 files changed, 38 insertions(+), 57 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 3f1eaaf06..39765fda4 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -268,49 +268,7 @@ def median_filter_cuda(arr, size=3, axis=0): return out -def _determine_nonfinite_kernel_idxs(x_idx, y_idx, kernel, shape_x, shape_y): - """ - Determine the proper kernel bounds for a given x, y index, and kernel size. - - Parameters - ---------- - x_idx : int - The x index for a given nonfinite value. - y_idx : int - The y index for a given nonfinite value. - shape_x : int - The max shape of the 2D projection in the x dimension. - shape_y : int - The max shape of the 2D projection in the y dimension. - - Returns - ------- - x_lower, x_higher, y_lower, y_higher : int - The integer kernel bounds to surround the nonfinite value with. - """ - - # Determining x kernel - x_idx_lower = x_idx - kernel - x_idx_higher = x_idx + kernel + 1 - - if x_idx_lower < 0: - x_idx_lower = 0 - if x_idx_higher > shape_x: - x_idx_higher = shape_x - - # Determining y kernel - y_idx_lower = y_idx - kernel - y_idx_higher = y_idx + kernel + 1 - - if y_idx_lower < 0: - y_idx_lower = 0 - if y_idx_higher > shape_y: - y_idx_higher = shape_y - - return x_idx_lower, x_idx_higher, y_idx_lower, y_idx_higher - - -def median_filter_nonfinite(data, size=3, callback=None): +def median_filter_nonfinite(data, size=3, callback=None, tshoot=False): """ Remove nonfinite values from a 3D array using an in-place 2D median filter. @@ -334,6 +292,13 @@ def median_filter_nonfinite(data, size=3, callback=None): ndarray The corrected 3D array with all nonfinite values removed based upon the local median value defined by the kernel size. + + Raises + ------ + ValueError + If the filter comes across a kernel only containing non-finite values a ValueError + is raised for the user to increase their kernel size to avoid replacing the current + pixel with 0. """ # Defining a callback function if None is provided if callback is None: @@ -357,6 +322,11 @@ def callback(total, description, unit): # Extracting kernel data and fining finite median kernel_cropped_data = projection[x_lower:x_higher, y_lower:y_higher] + + if len(kernel_cropped_data[np.isfinite(kernel_cropped_data)]) == 0: + raise ValueError("Found kernel containing only non-finite values.\ + Please increase kernel size") + median_corrected_data = np.median( kernel_cropped_data[np.isfinite(kernel_cropped_data)]) diff --git a/test/test_tomopy/test_misc/test_corr.py b/test/test_tomopy/test_misc/test_corr.py index cd2f2bf02..4fb6b4296 100644 --- a/test/test_tomopy/test_misc/test_corr.py +++ b/test/test_tomopy/test_misc/test_corr.py @@ -74,31 +74,42 @@ def test_median_filter(self): loop_dim(median_filter, read_file('cube.npy')) def test_median_filter_nonfinite(self): - # Add some random non-finite values to an array of all ones + + # Set a standard random value to make reproducible + np.random.seed(1) + data_org = np.ones(shape=(100, 100, 100)) - for i in range(50): - x = np.random.randint(0, 100) - y = np.random.randint(0, 100) - z = np.random.randint(0, 100) - data_org[z, x, y] = np.inf - x = np.random.randint(0, 100) - y = np.random.randint(0, 100) - z = np.random.randint(0, 100) - data_org[z, x, y] = -np.inf - x = np.random.randint(0, 100) - y = np.random.randint(0, 100) - z = np.random.randint(0, 100) - data_org[z, x, y] = np.nan + + # Add some random non-finite values to an array of all ones + for non_finite in [np.nan, np.inf, -np.inf]: + for i in range(50): + x = np.random.randint(0, 100) + y = np.random.randint(0, 100) + z = np.random.randint(0, 100) + data_org[z, x, y] = non_finite data_post_corr = median_filter_nonfinite( data_org.copy(), size=5, callback=None, ) + # All the post filtering values should be 1 because all of the finite # values are 1. assert np.all(data_post_corr == 1.0) + # Making sure filter raises ValueError when function finds a filter + # filled with non-finite values. + for non_finite in [np.nan, np.inf, -np.inf]: + data_org = np.empty((10, 10)) + data_org[:] = non_finite + self.assertRaises(ValueError, + median_filter_nonfinite, + test_array.copy(), + size=3, + callback=None, + ) + def test_remove_neg(self): assert_allclose(remove_neg([-2, -1, 0, 1, 2]), [0, 0, 0, 1, 2]) From 4d2ce3bddc8698f77ae229c4e426f363321aa851 Mon Sep 17 00:00:00 2001 From: WilliamJudge94 <37912237+WilliamJudge94@users.noreply.github.com> Date: Thu, 4 Nov 2021 14:05:23 -0500 Subject: [PATCH 30/35] BUG: Fixing variable name bug --- test/test_tomopy/test_misc/test_corr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_tomopy/test_misc/test_corr.py b/test/test_tomopy/test_misc/test_corr.py index 4fb6b4296..5b86c76bb 100644 --- a/test/test_tomopy/test_misc/test_corr.py +++ b/test/test_tomopy/test_misc/test_corr.py @@ -105,7 +105,7 @@ def test_median_filter_nonfinite(self): data_org[:] = non_finite self.assertRaises(ValueError, median_filter_nonfinite, - test_array.copy(), + data_org[:].copy(), size=3, callback=None, ) From 5c020e27795182fb33e4da8e77b6ea32e2b20644 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Wed, 10 Nov 2021 18:32:25 -0600 Subject: [PATCH 31/35] BUG: Prevent filter sweep in median_nonfinite Modifies tests so that they force the correct ValueError. Ensure that ValueError is always raised when there is a group of nonfinites larger than or equal to the filter size. Without making a copy of the current slice, filters would use already filtered values in the filter. --- source/tomopy/misc/corr.py | 148 ++++++++++++++---------- test/test_tomopy/test_misc/test_corr.py | 16 +-- 2 files changed, 93 insertions(+), 71 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 39765fda4..a41493fb4 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -45,7 +45,6 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # # POSSIBILITY OF SUCH DAMAGE. # # ######################################################################### - """ Module for data correction and masking functions. """ @@ -66,26 +65,26 @@ logger = logging.getLogger(__name__) - __author__ = "Doga Gursoy, William Judge" __credits__ = "Mark Rivers, Xianghui Xiao" __copyright__ = "Copyright (c) 2015, UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' -__all__ = ['adjust_range', - 'circ_mask', - 'gaussian_filter', - 'median_filter', - 'median_filter_cuda', - 'median_filter_nonfinite', - 'sobel_filter', - 'remove_nan', - 'remove_neg', - 'remove_outlier', - 'remove_outlier1d', - 'remove_outlier_cuda', - 'remove_ring', - 'enhance_projs_aps_1id', - ] +__all__ = [ + 'adjust_range', + 'circ_mask', + 'gaussian_filter', + 'median_filter', + 'median_filter_cuda', + 'median_filter_nonfinite', + 'sobel_filter', + 'remove_nan', + 'remove_neg', + 'remove_outlier', + 'remove_outlier1d', + 'remove_outlier_cuda', + 'remove_ring', + 'enhance_projs_aps_1id', +] def adjust_range(arr, dmin=None, dmax=None): @@ -152,11 +151,14 @@ def gaussian_filter(arr, sigma=3, order=0, axis=0, ncore=None): ncore = mproc.mp.cpu_count() with cf.ThreadPoolExecutor(ncore) as e: - slc = [slice(None)]*arr.ndim + slc = [slice(None)] * arr.ndim for i in range(arr.shape[axis]): slc[axis] = i - e.submit(filters.gaussian_filter, arr[tuple(slc)], sigma, - order=order, output=out[tuple(slc)]) + e.submit(filters.gaussian_filter, + arr[tuple(slc)], + sigma, + order=order, + output=out[tuple(slc)]) return out @@ -187,10 +189,12 @@ def median_filter(arr, size=3, axis=0, ncore=None): ncore = mproc.mp.cpu_count() with cf.ThreadPoolExecutor(ncore) as e: - slc = [slice(None)]*arr.ndim + slc = [slice(None)] * arr.ndim for i in range(arr.shape[axis]): slc[axis] = i - e.submit(filters.median_filter, arr[tuple(slc)], size=(size, size), + e.submit(filters.median_filter, + arr[tuple(slc)], + size=(size, size), output=out[tuple(slc)]) return out @@ -228,25 +232,25 @@ def median_filter_cuda(arr, size=3, axis=0): winAllow = range(2, 16) - if(axis != 0): + if (axis != 0): arr = np.swapaxes(arr, 0, axis) if size in winAllow: - loffset = int(size/2) - roffset = int((size-1)/2) + loffset = int(size / 2) + roffset = int((size - 1) / 2) prjsize = arr.shape[0] imsizex = arr.shape[2] imsizey = arr.shape[1] filter = tomocuda.mFilter(imsizex, imsizey, prjsize, size) - out = np.zeros( - shape=(prjsize, imsizey, imsizex), dtype=np.float32) + out = np.zeros(shape=(prjsize, imsizey, imsizex), dtype=np.float32) for step in range(prjsize): # im_noisecu = arr[:][step][:].astype(np.float32) im_noisecu = arr[step].astype(np.float32) im_noisecu = np.lib.pad(im_noisecu, ((loffset, roffset), - (loffset, roffset)), 'symmetric') + (loffset, roffset)), + 'symmetric') im_noisecu = im_noisecu.flatten() filter.setCuImage(im_noisecu) @@ -255,7 +259,7 @@ def median_filter_cuda(arr, size=3, axis=0): results = results.reshape(imsizey, imsizex) out[step] = results - if(axis != 0): + if (axis != 0): out = np.swapaxes(out, 0, axis) else: warnings.warn("Window size not support, using cpu median filter") @@ -309,6 +313,7 @@ def callback(total, description, unit): # Iterating throug each projection to save on RAM for projection in data: nonfinite_idx = np.nonzero(~np.isfinite(projection)) + projection_copy = projection.copy() # Iterating through each bad value and replace it with finite median for x_idx, y_idx in zip(*nonfinite_idx): @@ -320,11 +325,12 @@ def callback(total, description, unit): y_higher = min(data.shape[2], y_idx + (size // 2) + 1) # Extracting kernel data and fining finite median - kernel_cropped_data = projection[x_lower:x_higher, - y_lower:y_higher] + kernel_cropped_data = projection_copy[x_lower:x_higher, + y_lower:y_higher] if len(kernel_cropped_data[np.isfinite(kernel_cropped_data)]) == 0: - raise ValueError("Found kernel containing only non-finite values.\ + raise ValueError( + "Found kernel containing only non-finite values.\ Please increase kernel size") median_corrected_data = np.median( @@ -363,7 +369,7 @@ def sobel_filter(arr, axis=0, ncore=None): ncore = mproc.mp.cpu_count() with cf.ThreadPoolExecutor(ncore) as e: - slc = [slice(None)]*arr.ndim + slc = [slice(None)] * arr.ndim for i in range(arr.shape[axis]): slc[axis] = i e.submit(filters.sobel, arr[slc], output=out[slc]) @@ -455,14 +461,16 @@ def remove_outlier(arr, dif, size=3, axis=0, ncore=None, out=None): ncore, chnk_slices = mproc.get_ncore_slices(arr.shape[axis], ncore=ncore) - filt_size = [size]*arr.ndim + filt_size = [size] * arr.ndim filt_size[axis] = 1 with cf.ThreadPoolExecutor(ncore) as e: - slc = [slice(None)]*arr.ndim + slc = [slice(None)] * arr.ndim for i in range(ncore): slc[axis] = chnk_slices[i] - e.submit(filters.median_filter, arr[tuple(slc)], size=filt_size, + e.submit(filters.median_filter, + arr[tuple(slc)], + size=filt_size, output=tmp[tuple(slc)]) arr = dtype.as_float32(arr) @@ -510,17 +518,20 @@ def remove_outlier1d(arr, dif, size=3, axis=0, ncore=None, out=None): other_axes = [i for i in range(arr.ndim) if i != axis] largest = np.argmax([arr.shape[i] for i in other_axes]) lar_axis = other_axes[largest] - ncore, chnk_slices = mproc.get_ncore_slices( - arr.shape[lar_axis], ncore=ncore) - filt_size = [1]*arr.ndim + ncore, chnk_slices = mproc.get_ncore_slices(arr.shape[lar_axis], + ncore=ncore) + filt_size = [1] * arr.ndim filt_size[axis] = size with cf.ThreadPoolExecutor(ncore) as e: - slc = [slice(None)]*arr.ndim + slc = [slice(None)] * arr.ndim for i in range(ncore): slc[lar_axis] = chnk_slices[i] - e.submit(filters.median_filter, arr[slc], size=filt_size, - output=tmp[slc], mode='mirror') + e.submit(filters.median_filter, + arr[slc], + size=filt_size, + output=tmp[slc], + mode='mirror') with mproc.set_numexpr_threads(ncore): out = ne.evaluate('where(arr-tmp>=dif,tmp,arr)', out=out) @@ -568,13 +579,13 @@ def remove_outlier_cuda(arr, dif, size=3, axis=0): winAllow = range(2, 16) - if(axis != 0): + if (axis != 0): arr = np.swapaxes(arr, 0, axis) if size in winAllow: prjsize = arr.shape[0] - loffset = int(size/2) - roffset = int((size-1)/2) + loffset = int(size / 2) + roffset = int((size - 1) / 2) imsizex = arr.shape[2] imsizey = arr.shape[1] @@ -584,7 +595,8 @@ def remove_outlier_cuda(arr, dif, size=3, axis=0): for step in range(prjsize): im_noisecu = arr[step].astype(np.float32) im_noisecu = np.lib.pad(im_noisecu, ((loffset, roffset), - (loffset, roffset)), 'symmetric') + (loffset, roffset)), + 'symmetric') im_noisecu = im_noisecu.flatten() filter.setCuImage(im_noisecu) @@ -593,7 +605,7 @@ def remove_outlier_cuda(arr, dif, size=3, axis=0): results = results.reshape(imsizey, imsizex) out[step] = results - if(axis != 0): + if (axis != 0): out = np.swapaxes(out, 0, axis) else: warnings.warn("Window size not support, using cpu outlier removal") @@ -606,9 +618,18 @@ def remove_outlier_cuda(arr, dif, size=3, axis=0): return out -def remove_ring(rec, center_x=None, center_y=None, thresh=300.0, - thresh_max=300.0, thresh_min=-100.0, theta_min=30, - rwidth=30, int_mode='WRAP', ncore=None, nchunk=None, out=None): +def remove_ring(rec, + center_x=None, + center_y=None, + thresh=300.0, + thresh_max=300.0, + thresh_min=-100.0, + theta_min=30, + rwidth=30, + int_mode='WRAP', + ncore=None, + nchunk=None, + out=None): """ Remove ring artifacts from images in the reconstructed domain. Descriptions of parameters need to be more clear for sure. @@ -659,9 +680,9 @@ def remove_ring(rec, center_x=None, center_y=None, thresh=300.0, dz, dy, dx = rec.shape if center_x is None: - center_x = (dx - 1.0)/2.0 + center_x = (dx - 1.0) / 2.0 if center_y is None: - center_y = (dy - 1.0)/2.0 + center_y = (dy - 1.0) / 2.0 if int_mode.lower() == 'wrap': int_mode = 0 @@ -673,14 +694,14 @@ def remove_ring(rec, center_x=None, center_y=None, thresh=300.0, if not 0 <= theta_min < 180: raise ValueError("theta_min should be in the range [0 - 180)") - args = (center_x, center_y, dx, dy, dz, thresh_max, thresh_min, - thresh, theta_min, rwidth, int_mode) + args = (center_x, center_y, dx, dy, dz, thresh_max, thresh_min, thresh, + theta_min, rwidth, int_mode) axis_size = rec.shape[0] ncore, nchunk = mproc.get_ncore_nchunk(axis_size, ncore, nchunk) with cf.ThreadPoolExecutor(ncore) as e: for offset in range(0, axis_size, nchunk): - slc = np.s_[offset:offset+nchunk] + slc = np.s_[offset:offset + nchunk] e.submit(extern.c_remove_ring, out[slc], *args) return out @@ -768,17 +789,18 @@ def enhance_projs_aps_1id(imgstack, median_ks=5, ncore=None): ndarray 3D enhanced image stacks. """ - ncore = mproc.mp.cpu_count()-1 if ncore is None else ncore + ncore = mproc.mp.cpu_count() - 1 if ncore is None else ncore # need to use multiprocessing to speed up the process tmp = [] with cf.ProcessPoolExecutor(ncore) as e: for n_img in range(imgstack.shape[0]): - tmp.append(e.submit(_enhance_img, - imgstack[n_img, :, :], - median_ks, - ) - ) + tmp.append( + e.submit( + _enhance_img, + imgstack[n_img, :, :], + median_ks, + )) return np.stack([me.result() for me in tmp], axis=0) @@ -806,7 +828,7 @@ def _enhance_img(img, median_ks, normalized=True): wgt = _calc_histequal_wgt(img) img = medfilt2d(img, kernel_size=median_ks).astype(np.float64) img = ne.evaluate('(img**2)*wgt', out=img) - return img/img.max() if normalized else img + return img / img.max() if normalized else img def _calc_histequal_wgt(img): @@ -824,4 +846,4 @@ def _calc_histequal_wgt(img): histogram euqalization weights (0-1) in the same shape as original image """ - return (np.sort(img.flatten()).searchsorted(img) + 1)/np.prod(img.shape) + return (np.sort(img.flatten()).searchsorted(img) + 1) / np.prod(img.shape) diff --git a/test/test_tomopy/test_misc/test_corr.py b/test/test_tomopy/test_misc/test_corr.py index 5b86c76bb..fdbb7599b 100644 --- a/test/test_tomopy/test_misc/test_corr.py +++ b/test/test_tomopy/test_misc/test_corr.py @@ -101,14 +101,14 @@ def test_median_filter_nonfinite(self): # Making sure filter raises ValueError when function finds a filter # filled with non-finite values. for non_finite in [np.nan, np.inf, -np.inf]: - data_org = np.empty((10, 10)) - data_org[:] = non_finite - self.assertRaises(ValueError, - median_filter_nonfinite, - data_org[:].copy(), - size=3, - callback=None, - ) + data_org = np.empty((1, 3, 3)) + data_org[:, -2:, -2:] = non_finite + with self.assertRaises(ValueError): + result = median_filter_nonfinite( + data_org.copy(), + size=3, + callback=None, + ) def test_remove_neg(self): assert_allclose(remove_neg([-2, -1, 0, 1, 2]), [0, 0, 0, 1, 2]) From e5f82f7d4c6e4dcb6d608f168fd9f56ffaae5c67 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Wed, 10 Nov 2021 18:35:06 -0600 Subject: [PATCH 32/35] STY: Wrap lines at 80 chars --- source/tomopy/misc/corr.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index a41493fb4..e0008b065 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -294,15 +294,14 @@ def median_filter_nonfinite(data, size=3, callback=None, tshoot=False): Returns ------- ndarray - The corrected 3D array with all nonfinite values removed based upon the local - median value defined by the kernel size. + The corrected 3D array with all nonfinite values removed based upon the + local median value defined by the kernel size. Raises ------ ValueError - If the filter comes across a kernel only containing non-finite values a ValueError - is raised for the user to increase their kernel size to avoid replacing the current - pixel with 0. + If the filter comes across a kernel only containing non-finite values a + ValueError is raised for the user to increase their kernel size. """ # Defining a callback function if None is provided if callback is None: From 4f6644afc2ba653eb1a8171aa11153ec37852da5 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Wed, 10 Nov 2021 18:35:28 -0600 Subject: [PATCH 33/35] STY: Remove unused parameters --- source/tomopy/misc/corr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index e0008b065..847202141 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -272,7 +272,7 @@ def median_filter_cuda(arr, size=3, axis=0): return out -def median_filter_nonfinite(data, size=3, callback=None, tshoot=False): +def median_filter_nonfinite(data, size=3, callback=None): """ Remove nonfinite values from a 3D array using an in-place 2D median filter. From b6422656c6b1fdd6bf9da544a74e46d6117740f0 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Thu, 11 Nov 2021 16:56:20 -0600 Subject: [PATCH 34/35] DOC: Update docs and match new API with existing filter funcs --- doc/source/api/tomopy.misc.corr.rst | 3 ++- source/tomopy/misc/corr.py | 26 +++++++++++++------------- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/doc/source/api/tomopy.misc.corr.rst b/doc/source/api/tomopy.misc.corr.rst index 682c7fd56..967847275 100644 --- a/doc/source/api/tomopy.misc.corr.rst +++ b/doc/source/api/tomopy.misc.corr.rst @@ -9,12 +9,13 @@ .. rubric:: **Functions:** .. autosummary:: - + adjust_range circ_mask gaussian_filter median_filter median_filter_cuda + median_filter_nonfinite sobel_filter remove_nan remove_neg diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 847202141..f5d616025 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -272,7 +272,7 @@ def median_filter_cuda(arr, size=3, axis=0): return out -def median_filter_nonfinite(data, size=3, callback=None): +def median_filter_nonfinite(arr, size=3, callback=None): """ Remove nonfinite values from a 3D array using an in-place 2D median filter. @@ -281,8 +281,8 @@ def median_filter_nonfinite(data, size=3, callback=None): Parameters ---------- - data : ndarray - The 3D array of data with nonfinite values in it. + arr : ndarray + The 3D array with nonfinite values in it. size : int, optional The size of the filter. callback : func(total, description, unit) @@ -310,7 +310,7 @@ def callback(total, description, unit): pass # Iterating throug each projection to save on RAM - for projection in data: + for projection in arr: nonfinite_idx = np.nonzero(~np.isfinite(projection)) projection_copy = projection.copy() @@ -319,28 +319,28 @@ def callback(total, description, unit): # Determining the lower and upper bounds for kernel x_lower = max(0, x_idx - (size // 2)) - x_higher = min(data.shape[1], x_idx + (size // 2) + 1) + x_higher = min(arr.shape[1], x_idx + (size // 2) + 1) y_lower = max(0, y_idx - (size // 2)) - y_higher = min(data.shape[2], y_idx + (size // 2) + 1) + y_higher = min(arr.shape[2], y_idx + (size // 2) + 1) # Extracting kernel data and fining finite median - kernel_cropped_data = projection_copy[x_lower:x_higher, + kernel_cropped_arr = projection_copy[x_lower:x_higher, y_lower:y_higher] - if len(kernel_cropped_data[np.isfinite(kernel_cropped_data)]) == 0: + if len(kernel_cropped_arr[np.isfinite(kernel_cropped_arr)]) == 0: raise ValueError( "Found kernel containing only non-finite values.\ Please increase kernel size") - median_corrected_data = np.median( - kernel_cropped_data[np.isfinite(kernel_cropped_data)]) + median_corrected_arr = np.median( + kernel_cropped_arr[np.isfinite(kernel_cropped_arr)]) # Replacing bad data with finite median - projection[x_idx, y_idx] = median_corrected_data + projection[x_idx, y_idx] = median_corrected_arr - callback(data.shape[0], 'Nonfinite median filter', ' prjs') + callback(arr.shape[0], 'Nonfinite median filter', ' prjs') - return data + return arr def sobel_filter(arr, axis=0, ncore=None): From 301ee367d18ca6d18f2b9b18e2c531c33d4739e4 Mon Sep 17 00:00:00 2001 From: Daniel Ching Date: Thu, 11 Nov 2021 17:11:58 -0600 Subject: [PATCH 35/35] DOC: versionadded tags to documentation --- source/tomopy/misc/corr.py | 3 +++ source/tomopy/prep/alignment.py | 22 +++++++++++++--------- source/tomopy/prep/stripe.py | 16 +++++++++------- 3 files changed, 25 insertions(+), 16 deletions(-) diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index f5d616025..1c3eaf2c8 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -279,6 +279,8 @@ def median_filter_nonfinite(arr, size=3, callback=None): The 2D selective median filter is applied along the last two axes of the array. + .. versionadded:: 1.11 + Parameters ---------- arr : ndarray @@ -302,6 +304,7 @@ def median_filter_nonfinite(arr, size=3, callback=None): ValueError If the filter comes across a kernel only containing non-finite values a ValueError is raised for the user to increase their kernel size. + """ # Defining a callback function if None is provided if callback is None: diff --git a/source/tomopy/prep/alignment.py b/source/tomopy/prep/alignment.py index edd0f0ab9..215c7f4c2 100644 --- a/source/tomopy/prep/alignment.py +++ b/source/tomopy/prep/alignment.py @@ -82,7 +82,7 @@ 'remove_slits_aps_1id', 'distortion_correction_proj', 'distortion_correction_sino', - 'load_distortion_coefs', + 'load_distortion_coefs', ] @@ -915,6 +915,8 @@ def distortion_correction_proj(tomo, xcenter, ycenter, list_fact, Apply distortion correction to projections using the polynomial model. Coefficients are calculated using Vounwarp package :cite:`Vo:15`. + .. versionadded:: 1.7 + Parameters ---------- tomo : ndarray @@ -922,7 +924,7 @@ def distortion_correction_proj(tomo, xcenter, ycenter, list_fact, xcenter : float Center of distortion in x-direction. From the left of the image. ycenter : float - Center of distortion in y-direction. From the top of the image. + Center of distortion in y-direction. From the top of the image. list_fact : list of floats Polynomial coefficients of the backward model. ncore : int, optional @@ -948,17 +950,17 @@ def distortion_correction_proj(tomo, xcenter, ycenter, list_fact, def _unwarp_image_backward(mat, xcenter, ycenter, list_fact): """ Unwarp an image using the polynomial model. - + Parameters ---------- mat : 2D array. - xcenter : float + xcenter : float Center of distortion in x-direction. From the left of the image. ycenter : float Center of distortion in y-direction. From the top of the image. - list_fact : list of floats + list_fact : list of floats Polynomial coefficients of the backward model. - + Returns ------- 2D array @@ -982,7 +984,7 @@ def _unwarp_image_backward(mat, xcenter, ycenter, list_fact): def _distortion_correction_proj(tomo, xcenter, ycenter, list_fact): for m in np.arange(tomo.shape[0]): proj = tomo[m, :, :] - proj = _unwarp_image_backward(proj, xcenter, ycenter, list_fact) + proj = _unwarp_image_backward(proj, xcenter, ycenter, list_fact) tomo[m, :, :] = proj @@ -992,6 +994,8 @@ def distortion_correction_sino(tomo, ind, xcenter, ycenter, list_fact): the polynomial model. Coefficients are calculated using Vounwarp package :cite:`Vo:15`. + .. versionadded:: 1.7 + Parameters ---------- tomo : ndarray @@ -1001,7 +1005,7 @@ def distortion_correction_sino(tomo, ind, xcenter, ycenter, list_fact): xcenter : float Center of distortion in x-direction. From the left of the image. ycenter : float - Center of distortion in y-direction. From the top of the image. + Center of distortion in y-direction. From the top of the image. list_fact : list of floats Polynomial coefficients of the backward model. @@ -1021,7 +1025,7 @@ def distortion_correction_sino(tomo, ind, xcenter, ycenter, list_fact): yd_list = np.clip(ycenter + flist * yu, 0, height - 1) yd_min = np.int16(np.floor(np.amin(yd_list))) yd_max = np.int16(np.ceil(np.amax(yd_list))) + 1 - yd_list = yd_list - yd_min + yd_list = yd_list - yd_min sino = np.zeros((depth, width), dtype=np.float32) indices = yd_list, xd_list for i in np.arange(depth): diff --git a/source/tomopy/prep/stripe.py b/source/tomopy/prep/stripe.py index 26303c0f1..cd80d47b4 100644 --- a/source/tomopy/prep/stripe.py +++ b/source/tomopy/prep/stripe.py @@ -436,7 +436,7 @@ def remove_stripe_based_filtering( tomo : ndarray 3D tomographic data. sigma : float - Sigma of the Gaussian window which is used to separate + Sigma of the Gaussian window which is used to separate the low-pass and high-pass components of the intensity profiles of each column. Recommended values: 3->10. size : int @@ -550,7 +550,7 @@ def _create_2d_window(nrow, ncol, sigma, pad): ---------- nrow : int Height of the window. - ncol: int + ncol: int Width of the window. sigma: tuple of 2 floats Sigmas of the window. @@ -590,7 +590,7 @@ def _2d_filter(mat, win2d, matsign, pad): mat : 2D array of floats nrow : int Height of the window. - ncol: int + ncol: int Width of the window. sigma: tuple of 2 floats Sigmas of the window. @@ -598,7 +598,7 @@ def _2d_filter(mat, win2d, matsign, pad): Padding. Returns - ------- + ------- Filtered image. """ matpad = np.pad(mat, ((0, 0), (pad, pad)), mode='edge') @@ -649,7 +649,7 @@ def remove_large_stripe(tomo, snr=3, size=51, drop_ratio=0.1, norm=True, 3D tomographic data. snr : float Ratio used to locate of large stripes. - Greater is less sensitive. + Greater is less sensitive. size : int Window size of the median filter. drop_ratio : float, optional @@ -754,7 +754,7 @@ def remove_dead_stripe(tomo, snr=3, size=51, norm=True, 3D tomographic data. snr : float Ratio used to detect locations of large stripes. - Greater is less sensitive. + Greater is less sensitive. size : int Window size of the median filter. norm : bool, optional @@ -826,7 +826,7 @@ def remove_all_stripe(tomo, snr=3, la_size=61, sm_size=21, dim=1, 3D tomographic data. snr : float Ratio used to locate large stripes. - Greater is less sensitive. + Greater is less sensitive. la_size : int Window size of the median filter to remove large stripes. sm_size : int @@ -868,6 +868,8 @@ def remove_stripe_based_interpolation(tomo, snr=3, size=31, drop_ratio=0.1, Remove most types of stripe artifacts from sinograms based on interpolation. Derived from algorithm 4, 5, and 6 in :cite:`Vo:18`. + .. versionadded:: 1.9 + Parameters ---------- tomo : ndarray