Commit 3adf19c9 authored by Pierre Paleo's avatar Pierre Paleo
Browse files

Merge branch 'halftomo_cor_outside' into 'master'

Halftomo: support CoR outside detector region

See merge request !129
parents d0d799ba 7f2131a3
Pipeline #46179 failed with stages
in 5 minutes and 39 seconds
__version__ = "2020.5.18-dev"
__version__ = "2020.5.19-dev"
__nabu_modules__ = [
"app",
"cuda",
......
......@@ -30,8 +30,6 @@ __global__ void halftomo_kernel(
output[y * n_x2 + x] = sinogram[y * n_x + x];
}
// output[:, nx - d : nx] = (1 - weights) * sino[:n_a2, nx - d :]
// + weights * sino[n_a2:, ::-1][:, d : 2 * d]
else if (x < n_x) { // x in [n_x - d, n_x [
......@@ -51,3 +49,49 @@ __global__ void halftomo_kernel(
}
}
// Same as previous kernel, but now the CoR is outside the image support
// i.e rotation_axis_position >= n_x
// Weigting is different !
__global__ void halftomo_kernel_cor_outside(
float* sinogram,
float* output,
float* weights,
int n_angles,
int n_x,
int rotation_axis_position
) {
int x = blockDim.x * blockIdx.x + threadIdx.x;
int y = blockDim.y * blockIdx.y + threadIdx.y;
int n_a2 = n_angles / 2;
int d = rotation_axis_position - n_x;
int n_x2 = 2 * rotation_axis_position;
if ((x >= n_x2) || (y >= n_a2)) return;
// output[:, :nx] = sino[:n_a2, :nx]
if (x < n_x) {
output[y * n_x2 + x] = sinogram[y * n_x + x];
}
// output[:, nx : nx + d] = (1 - weights) * sino[:n_a2, -d :][:, ::-1]
else if (x < n_x + d) {
float w = weights[x - n_x];
// output[y * n_x2 + x] = (1.0f - w) * sinogram[y*n_x + n_x - 1 - (n_x - x)];
output[y * n_x2 + x] = (1.0f - w) * sinogram[y*n_x + 2 * n_x - 1 - x];
}
// output[:, nx + d : nx + 2*d] = weights * sino[n_a2:, ::-1][:, :d] = sino[n_a2:, -d-1:-n_x-1:-1]
else if (x < n_x + 2 * d) {
float w = weights[x - (n_x + d)];
output[y * n_x2 + x] = w * sinogram[(n_a2 + y)*n_x + x - 2*d];
}
// output[:, nx+2*d:] = sino[n_a2:, ::-1][:, 2 * d :] = sino[n_a2:, -2*d-1:-n_x-1:-1]
else {
output[y * n_x2 + x] = sinogram[(n_a2 + y)*n_x + (n_x2 - 1 - x)];
}
}
\ No newline at end of file
......@@ -113,11 +113,6 @@ class SinoProcessing(object):
def _radios_to_sinos_halftomo(self, radios, sinos):
# TODO
if not(np.isscalar(self.rot_center)):
raise NotImplementedError("Half tomo with varying rotation axis position is not implemented yet")
#
n_a, n_z, n_x = radios.shape
n_a2 = n_a // 2
out_dwidth = 2 * self._rot_center_int
......@@ -176,7 +171,11 @@ def convert_halftomo(sino, rotation_axis_position):
"""
assert sino.ndim == 2
assert (sino.shape[0] % 2) == 0
na, nx = sino.shape
if rotation_axis_position >= nx:
return _convert_halftomo_right(sino, rotation_axis_position)
na2 = na // 2
r = rotation_axis_position
d = nx - r
......@@ -194,6 +193,32 @@ def convert_halftomo(sino, rotation_axis_position):
return res
def _convert_halftomo_right(sino, rotation_axis_position):
"""
Converts a sinogram into a sinogram with extended FOV with the "half tomography"
setting, with a CoR outside the image support.
"""
assert sino.ndim == 2
assert (sino.shape[0] % 2) == 0
assert rotation_axis_position > sino.shape[1]
na, nx = sino.shape
na2 = na // 2
r = rotation_axis_position
d = r - nx
res = np.zeros((na2, 2 * r), dtype="f")
sino1 = sino[:na2, :]
sino2 = sino[na2:, ::-1]
w = np.linspace(0, 1, d, endpoint=True) # has to be adapted when the sample actually goes beyond this point
res[:, :nx] = sino1[:, :]
res[:, nx:nx+d] = sino1[:, -d:][:, ::-1] * (1 - w)
res[:, nx+d:nx+2*d] = sino2[:, :d][:, ::-1] * w
res[:, nx+2*d:] = sino2[:, :]
return res
class SinoNormalization(SinoProcessing):
"""
......
......@@ -26,8 +26,12 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing):
def _init_cuda_halftomo(self):
if not(self.halftomo):
return
if self._rot_center_int < self.n_x:
kernel_name = "halftomo_kernel"
else:
kernel_name = "halftomo_kernel_cor_outside"
self.halftomo_kernel = CudaKernel(
"halftomo_kernel",
kernel_name,
get_cuda_srcfile("halftomo.cu"),
signature="PPPiii",
)
......@@ -40,7 +44,7 @@ class CudaSinoProcessing(SinoProcessing, CudaProcessing):
1
)
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.halftomo_weights = np.linspace(0, 1, abs(d), endpoint=True, dtype="f")
self.d_halftomo_weights = garray.to_gpu(self.halftomo_weights)
if self._halftomo_flip:
self.xflip_kernel = CudaKernel(
......
Supports Markdown
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