diff --git a/httomolibgpu/misc/morph.py b/httomolibgpu/misc/morph.py index f34d0d4a..e58d32d6 100644 --- a/httomolibgpu/misc/morph.py +++ b/httomolibgpu/misc/morph.py @@ -45,6 +45,7 @@ __all__ = [ "sino_360_to_180", "data_resampler", + "average_projection_frames", ] @@ -305,3 +306,65 @@ def data_resampler( if expanded: scaled_data = cp.squeeze(scaled_data, axis=axis) return scaled_data + + +def average_projection_frames( + data: cp.ndarray, + projection_averaging_factor: int = 2, +) -> cp.ndarray: + """ + This method averages/downsamples by averaging data along the angular direction based on the provided projection_averaging_factor. + + Parameters + ---------- + data : cp.ndarray + 3d cupy array given as (angles, detY, detX). + projection_averaging_factor : int + to average every (defined by the provided factor) consecutive projections into one effective projection. + + Raises + ---------- + ValueError: When data is not 3D + + Returns + ------- + cp.ndarray: 3D cupy array with averaged projection data + """ + + ### Data and parameters checks ### + methods_name = "average_projection_frames" + __check_if_data_3D_array(data, methods_name) + __check_if_data_correct_type( + data, accepted_type=["float32", "uint16"], methods_name=methods_name + ) + __check_if_positive_nonzero( + projection_averaging_factor, + "projection_averaging_factor", + True, + True, + methods_name, + ) + + ################################### + if projection_averaging_factor > 1: + k = projection_averaging_factor + n_proj = data.shape[0] + + n_full = n_proj // k + remainder = n_proj % k + + n_out = n_full + (remainder > 0) + + averaged = cp.empty((n_out, *data.shape[1:]), dtype=data.dtype) + + if n_full: + averaged[:n_full] = ( + data[: n_full * k].reshape(n_full, k, *data.shape[1:]).mean(axis=1) + ) + + if remainder: + averaged[-1] = data[n_full * k :].mean(axis=0) + + return averaged + else: + return data diff --git a/tests/test_misc/test_morph.py b/tests/test_misc/test_morph.py index d7b015a1..065cab06 100644 --- a/tests/test_misc/test_morph.py +++ b/tests/test_misc/test_morph.py @@ -4,7 +4,11 @@ from cupy.cuda import nvtx import pytest from numpy.testing import assert_allclose -from httomolibgpu.misc.morph import sino_360_to_180, data_resampler +from httomolibgpu.misc.morph import ( + sino_360_to_180, + data_resampler, + average_projection_frames, +) @pytest.mark.parametrize( @@ -34,9 +38,9 @@ def test_sino_360_to_180_wrong_dims(ensure_clean_memory, shape): @pytest.mark.parametrize("axis", [0, 1, 2]) def test_data_resampler(data, axis, ensure_clean_memory): newshape = [60, 80] - scaled_data = data_resampler( - data, newshape=newshape, axis=axis, interpolation="linear" - ).get() + scaled_data = cp.asnumpy( + data_resampler(data, newshape=newshape, axis=axis, interpolation="linear") + ) assert scaled_data.ndim == 3 if axis == 0: @@ -52,6 +56,58 @@ def test_data_resampler(data, axis, ensure_clean_memory): assert scaled_data.flags.c_contiguous +@pytest.mark.parametrize("projection_averaging_factor", [1, 2, 3, 7]) +def test_average_projection_frames_testdata( + data, projection_averaging_factor, ensure_clean_memory +): + averaged_data = cp.asnumpy( + average_projection_frames( + data, projection_averaging_factor=projection_averaging_factor + ) + ) + + assert averaged_data.ndim == 3 + + if projection_averaging_factor == 1: + assert averaged_data.shape == (180, 128, 160) + if projection_averaging_factor == 2: + assert averaged_data.shape == (90, 128, 160) + assert_allclose(np.max(averaged_data), 1089) + if projection_averaging_factor == 3: + assert averaged_data.shape == (60, 128, 160) + assert_allclose(np.max(averaged_data), 1070) + if projection_averaging_factor == 7: + assert averaged_data.shape == (26, 128, 160) + assert_allclose(np.max(averaged_data), 1044) + assert averaged_data.dtype == np.uint16 + assert averaged_data.flags.c_contiguous + + +@pytest.mark.parametrize("projection_averaging_factor", [2, 3]) +def test_average_projection_frames_randdata( + projection_averaging_factor, ensure_clean_memory +): + data_host = ( + np.random.random_sample(size=(1801, 5, 2560)).astype(np.float32) * 2.0 + 0.001 + ) + data = cp.asarray(data_host, dtype=np.float32) + + averaged_data = cp.asnumpy( + average_projection_frames( + data, projection_averaging_factor=projection_averaging_factor + ) + ) + + assert averaged_data.ndim == 3 + + if projection_averaging_factor == 2: + assert averaged_data.shape == (901, 5, 2560) + if projection_averaging_factor == 3: + assert averaged_data.shape == (601, 5, 2560) + assert averaged_data.dtype == np.float32 + assert averaged_data.flags.c_contiguous + + @pytest.mark.parametrize("rotation", ["left", "right"]) @pytest.mark.perf def test_sino_360_to_180_performance(ensure_clean_memory, rotation):