diff --git a/nabu/app/fullfield.py b/nabu/app/fullfield.py index 957c77342d731c86fae13cd72d13ba839a932bc4..95bfbd51d3c5d1aa0c15071c3e8dee416e66199a 100644 --- a/nabu/app/fullfield.py +++ b/nabu/app/fullfield.py @@ -459,6 +459,7 @@ class FullFieldPipeline: self._rec_roi = (x_s, x_e, y_s, y_e) self._allocate_recs(y_e - y_s, x_e - x_s) + @use_options("reconstruction", "reconstruction") def _init_reconstruction(self): options = self.processing_options["reconstruction"] @@ -468,6 +469,8 @@ class FullFieldPipeline: if options["enable_halftomo"]: rot_center = options["rotation_axis_position_halftomo"] + if self.sino_builder._halftomo_flip: + rot_center = self.sino_builder.rot_center else: rot_center = options["rotation_axis_position"] if options.get("cor_estimated_auto", False): diff --git a/nabu/cuda/src/ElementOp.cu b/nabu/cuda/src/ElementOp.cu index 22bd2da195ef1ba4bfc28618d6ac3a89cac76e0f..909a09df8cf16bcc4dec67a92d6ed61df1737356 100644 --- a/nabu/cuda/src/ElementOp.cu +++ b/nabu/cuda/src/ElementOp.cu @@ -86,3 +86,20 @@ __global__ void nlog(float* array, int Nx, int Ny, int Nz, float clip_min, float #endif array[pos] = -logf(val); } + + + +// Reverse elements of a 2D array along "x", i.e: +// arr = arr[:, ::-1] +// launched with grid (Nx/2, Ny) +__global__ void reverse2D_x(float* array, int Nx, int Ny) { + uint x = blockDim.x * blockIdx.x + threadIdx.x; + uint y = blockDim.y * blockIdx.y + threadIdx.y; + if ((x >= Nx/2) || (y >= Ny)) return; + uint pos = y*Nx + x; + uint pos2 = y*Nx + (Nx - 1 - x); + float tmp = array[pos]; + array[pos] = array[pos2]; + array[pos2] = tmp; +} + diff --git a/nabu/preproc/sinogram.py b/nabu/preproc/sinogram.py index d1ff14ac3b64269edae5c9562b76da9ad3dfbae5..a9a6242abd86d09b39e9be55726d4a1447d5817e 100644 --- a/nabu/preproc/sinogram.py +++ b/nabu/preproc/sinogram.py @@ -71,7 +71,15 @@ class SinoProcessing(object): ) self.rot_center = rot_center self._rot_center_int = int(round(self.rot_center)) - + # If CoR is on the left: "flip" the logic + self._halftomo_flip = False + rc = self._rot_center_int + if rc < (self.n_x - 1)/2: + rc = self.n_x - 1 - rc + self._rot_center_int = rc + self.rot_center = self.n_x - self.rot_center + self._halftomo_flip = True + # def _set_halftomo(self, halftomo): self.halftomo = halftomo @@ -123,7 +131,12 @@ class SinoProcessing(object): else: sinos = np.zeros(self.sinos_halftomo_shape, dtype=np.float32) for i in range(n_z): - sinos[i][:] = convert_halftomo(radios[:, i, :], self._rot_center_int) + radio = radios[:, i, :] + if self._halftomo_flip: + radio = radio[:, ::-1] + sinos[i][:] = convert_halftomo(radio, self._rot_center_int) + if self._halftomo_flip: + sinos[i][:] = sinos[i][:, ::-1] return sinos diff --git a/nabu/preproc/sinogram_cuda.py b/nabu/preproc/sinogram_cuda.py index aca5caed78e7a58a99bcd6312b8f767c24c08590..309a0c0557d0cb7787d970fbf4e0cc684448514f 100644 --- a/nabu/preproc/sinogram_cuda.py +++ b/nabu/preproc/sinogram_cuda.py @@ -42,6 +42,20 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing): d = self.n_x - rc # will have to be adapted for varying axis pos self.halftomo_weights = np.linspace(0, 1, d, endpoint=True, dtype="f") self.d_halftomo_weights = garray.to_gpu(self.halftomo_weights) + if self._halftomo_flip: + self.xflip_kernel = CudaKernel( + "reverse2D_x", + get_cuda_srcfile("ElementOp.cu"), + signature="Pii" + ) + blk = (32, 32, 1) + self._xflip_blksize = blk + self._xflip_gridsize_1 = ( + updiv(self.n_x, blk[0]), + updiv(self.n_angles, blk[1]), + 1 + ) + self._xflip_gridsize_2 = self._halftomo_gridsize # Overwrite parent method @@ -80,15 +94,28 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing): else: sinos = garray.zeros(self.sinos_halftomo_shape, dtype=np.float32) + # need to use a contiguous 2D, array otherwise kernel does not work + d_radio = radios[:, 0, :].copy() for i in range(n_z): + d_radio[:] = radios[:, i, :].copy() + if self._halftomo_flip: + self.xflip_kernel( + d_radio, n_x, n_a, + grid=self._xflip_gridsize_1, block=self._xflip_blksize + ) self.halftomo_kernel( - radios[:, i, :].copy(), # need to copy this 2D, array otherwise kernel does not work + d_radio, sinos[i], self.d_halftomo_weights, n_a, n_x, rc, grid=self._halftomo_gridsize, block=self._halftomo_blksize ) + if self._halftomo_flip: + self.xflip_kernel( + sinos[i], 2*rc, n_a2, + grid=self._xflip_gridsize_2, block=self._xflip_blksize + ) return sinos @@ -125,7 +152,3 @@ class CudaSinoNormalization(SinoNormalization, CudaProcessing): sino, *self._cuda_kernel_args, **self._cuda_kernel_kwargs ) return sino - - - - diff --git a/nabu/preproc/tests/test_halftomo.py b/nabu/preproc/tests/test_halftomo.py index 627706fa0d9a8f0f462e7ffdcab6496d37dd6acf..2a739ed522e89a7851db81d392120a83565f3f11 100644 --- a/nabu/preproc/tests/test_halftomo.py +++ b/nabu/preproc/tests/test_halftomo.py @@ -5,12 +5,13 @@ import h5py from nabu.testutils import compare_arrays, utilstest from nabu.preproc.sinogram import SinoProcessing, convert_halftomo from nabu.cuda.utils import __has_pycuda__ + if __has_pycuda__: import pycuda.gpuarray as garray from nabu.preproc.sinogram_cuda import CudaSinoProcessing -@pytest.fixture(scope='class') +@pytest.fixture(scope="class") def bootstrap(request): cls = request.cls radio, sino_ref, cor = get_data_h5("halftomo.h5") @@ -26,37 +27,87 @@ def bootstrap(request): def get_data_h5(*dataset_path): dataset_relpath = os.path.join(*dataset_path) dataset_path = utilstest.getfile(dataset_relpath) - with h5py.File(dataset_path, 'r') as hf: + with h5py.File(dataset_path, "r") as hf: radio = hf["entry/radio/results/data"][()] sino = hf["entry/sino/results/data"][()] cor = hf["entry/sino/configuration/configuration/rotation_axis_position"][()] return radio, sino, cor - -@pytest.mark.usefixtures('bootstrap') +@pytest.mark.usefixtures("bootstrap") class TestHalftomo: - def test_halftomo(self): sino_processing = SinoProcessing( - radios_shape=self.radios.shape, - rot_center=self.rot_center, - halftomo=True + radios_shape=self.radios.shape, rot_center=self.rot_center, halftomo=True ) sinos_halftomo = sino_processing.radios_to_sinos(self.radios) - _, err = compare_arrays(sinos_halftomo[0], self.sino_ref, self.tol, return_residual=True) - assert err < self.tol, "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" + _, err = compare_arrays( + sinos_halftomo[0], self.sino_ref, self.tol, return_residual=True + ) + assert ( + err < self.tol + ), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" - @pytest.mark.skipif(not(__has_pycuda__), reason="Need pycuda for this test") + @pytest.mark.skipif(not (__has_pycuda__), reason="Need pycuda for this test") def test_cuda_halftomo(self): sino_processing = CudaSinoProcessing( - radios_shape=self.radios.shape, - rot_center=self.rot_center, - halftomo=True + radios_shape=self.radios.shape, rot_center=self.rot_center, halftomo=True ) d_radios = garray.to_gpu(self.radios) d_sinos = garray.zeros(sino_processing.sinos_halftomo_shape, "f") sino_processing.radios_to_sinos(d_radios, output=d_sinos) sino_halftomo = d_sinos.get()[0] - _, err = compare_arrays(sino_halftomo, self.sino_ref, self.tol, return_residual=True) - assert err < self.tol, "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" + _, err = compare_arrays( + sino_halftomo, self.sino_ref, self.tol, return_residual=True + ) + assert ( + err < self.tol + ), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" + + @staticmethod + def _flip_array(arr): + if arr.ndim == 2: + return np.fliplr(arr) + res = np.zeros_like(arr) + for i in range(arr.shape[0]): + res[i] = np.fliplr(arr[i]) + return res + + def test_halftomo_left(self): + na, nz, nx = self.radios.shape + left_cor = nx - 1 - self.rot_center + radios = self._flip_array(self.radios) + sino_processing = SinoProcessing( + radios_shape=radios.shape, rot_center=left_cor, halftomo=True + ) + sinos_halftomo = sino_processing.radios_to_sinos(radios) + _, err = compare_arrays( + sinos_halftomo[0], + self._flip_array(self.sino_ref), + self.tol, + return_residual=True, + ) + assert ( + err < self.tol + ), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" + + @pytest.mark.skipif(not (__has_pycuda__), reason="Need pycuda for this test") + def test_cuda_halftomo_left(self): + na, nz, nx = self.radios.shape + left_cor = nx - 1 - self.rot_center + radios = self._flip_array(self.radios) + sino_processing = CudaSinoProcessing( + radios_shape=radios.shape, rot_center=left_cor, halftomo=True + ) + d_radios = garray.to_gpu(radios) + d_sinos = garray.zeros(sino_processing.sinos_halftomo_shape, "f") + sino_processing.radios_to_sinos(d_radios, output=d_sinos) + sino_halftomo = d_sinos.get()[0] + _, err = compare_arrays( + sino_halftomo, self._flip_array(self.sino_ref), + self.tol, return_residual=True + ) + assert ( + err < self.tol + ), "Something wrong with SinoProcessing.radios_to_sino, halftomo=True" + diff --git a/nabu/resources/computations.py b/nabu/resources/computations.py index ec7b600ba4e7a2f0a7975c55d8336fcfee65e9cb..79a7ff68ed57db65de527f46c48f6e597506c37b 100644 --- a/nabu/resources/computations.py +++ b/nabu/resources/computations.py @@ -69,7 +69,9 @@ def estimate_required_memory(process_config, chunk_size=None, radios_and_slices= if "rotation_axis_position_halftomo" in rec_config: # Slice has a different shape in half acquisition. # Not sure if start_x, ..., end_y are supported - Nx_rec = 2 * rec_config["rotation_axis_position_halftomo"] + rc = rec_config["rotation_axis_position_halftomo"] + Nx_rec = 2 * rc + Nx_rec = max(2 * rc, 2 * (Nx - rc)) Ny_rec = Nx_rec else: Nx_rec = (rec_config["end_x"] - rec_config["start_x"] + 1) diff --git a/nabu/resources/dataset_validator.py b/nabu/resources/dataset_validator.py index 22dd5d31d542dfd28d61da35a3737b200dfcae2f..4cc2bf9b7abf8f78963fe644f16b49a712b7ed3c 100644 --- a/nabu/resources/dataset_validator.py +++ b/nabu/resources/dataset_validator.py @@ -50,6 +50,16 @@ class NabuValidator(object): return res + def _get_nx_ny(self): + if self.nabu_config["reconstruction"]["enable_halftomo"]: + cor = int(round(self.dataset_infos.axis_position)) + nx = max(2*cor, 2 * (self.dataset_infos.radio_dims[0] - 1 - cor)) + else: + nx, nz = self.dataset_infos.radio_dims + ny = nx + return nx, ny + + def convert_negative_indices(self): """ Convert any negative index to the corresponding positive index. @@ -59,8 +69,7 @@ class NabuValidator(object): if self.nabu_config["reconstruction"]["enable_halftomo"]: if self.dataset_infos.axis_position is None: raise ValueError("Cannot use rotation axis position in the middle of the detector when half tomo is enabled") - cor = int(round(self.dataset_infos.axis_position)) - ny = nx = 2*cor + nx, ny = self._get_nx_ny() what = ( ("reconstruction", "start_x", nx), ("reconstruction", "end_x", nx), @@ -199,8 +208,7 @@ class NabuValidator(object): nx, nz = self.dataset_infos.radio_dims rec_params = self.nabu_config["reconstruction"] if rec_params["enable_halftomo"]: - cor = int(round(self.dataset_infos.axis_position)) - ny = nx = 2*cor + ny, nx = self._get_nx_ny() what = ( ("start_x", "end_x", nx), ("start_y", "end_y", nx), diff --git a/nabu/resources/processconfig.py b/nabu/resources/processconfig.py index 40309d1619843e3c4c5f7c54eb1a62f320469151..a8645b58997b4c028d2fe46f3b620ef04e912454 100644 --- a/nabu/resources/processconfig.py +++ b/nabu/resources/processconfig.py @@ -220,7 +220,7 @@ class ProcessConfig: if rec_options["enable_halftomo"]: rec_options["angles"] = rec_options["angles"][:rec_options["angles"].size//2] cor_i = int(round(rec_options["rotation_axis_position"])) - # New key + # New keys rec_options["rotation_axis_position_halftomo"] = (2*cor_i-1)/2. # New key rec_options["cor_estimated_auto"] = isinstance(nabu_config["reconstruction"]["rotation_axis_position"], str)