Commit e5eb936f authored by Pierre Paleo's avatar Pierre Paleo
Browse files

Halftomo: support CoR on the left

parent 86dba6ba
......@@ -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):
......
......@@ -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;
}
......@@ -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
......
......@@ -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
......@@ -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"
......@@ -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)
......
......@@ -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),
......
......@@ -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)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment