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

Merge branch 'polar_backproj' into 'master'

Polar FBP

See merge request !114
parents 51d30ba4 e8a14908
Pipeline #47051 passed with stages
in 5 minutes and 55 seconds
__version__ = "2021.1.2-dev" __version__ = "2021.1.3-dev"
__nabu_modules__ = [ __nabu_modules__ = [
"app", "app",
"cuda", "cuda",
......
#ifndef SHARED_SIZE
#define SHARED_SIZE 256
#endif
texture<float, 2, cudaReadModeElementType> tex_projections;
__global__ void backproj_polar(
float* d_slice,
int num_projs,
int num_bins,
float axis_position,
int n_x,
int n_y,
int offset_x,
int offset_y,
float* d_cos,
float* d_msin,
float scale_factor
)
{
int i_r = offset_x + blockDim.x * blockIdx.x + threadIdx.x;
int i_theta = offset_y + blockDim.y * blockIdx.y + threadIdx.y;
float r = i_r - axis_position;
float x = r * d_cos[i_theta];
float y = - r * d_msin[i_theta];
/*volatile*/ __shared__ float s_cos[SHARED_SIZE];
/*volatile*/ __shared__ float s_msin[SHARED_SIZE];
int next_fetch = 0;
int tid = threadIdx.y * blockDim.x + threadIdx.x;
float costheta, msintheta;
float h1;
float sum1 = 0.0f;
for (int proj = 0; proj < num_projs; proj++) {
if (proj == next_fetch) {
// Fetch SHARED_SIZE values to shared memory
__syncthreads();
if (next_fetch + tid < num_projs) {
s_cos[tid] = d_cos[next_fetch + tid];
s_msin[tid] = d_msin[next_fetch + tid];
}
next_fetch += SHARED_SIZE;
__syncthreads();
}
costheta = s_cos[proj - (next_fetch - SHARED_SIZE)];
msintheta = s_msin[proj - (next_fetch - SHARED_SIZE)];
float c1 = fmaf(costheta, x, axis_position); // cos(theta)*xr + axis_pos
float s1 = fmaf(msintheta, y, 0.0f); // -sin(theta)*yr
h1 = c1 + s1;
if (h1 >= 0 && h1 < num_bins) sum1 += tex2D(tex_projections, h1 + 0.5f, proj + 0.5f);
}
// useful only if n_x < blocksize_x or n_y < blocksize_y
if (i_r >= n_x) return;
if (i_theta >= n_y) return;
d_slice[i_theta*(n_x) + i_r] = sum1 * scale_factor;
}
...@@ -14,6 +14,10 @@ class Backprojector(CudaProcessing): ...@@ -14,6 +14,10 @@ class Backprojector(CudaProcessing):
""" """
Cuda Backprojector. Cuda Backprojector.
""" """
cuda_fname = "backproj.cu"
cuda_kernel_name = "backproj"
def __init__( def __init__(
self, self,
sino_shape, sino_shape,
...@@ -237,8 +241,8 @@ class Backprojector(CudaProcessing): ...@@ -237,8 +241,8 @@ class Backprojector(CudaProcessing):
shared_size += int(np.prod(block)) shared_size += int(np.prod(block))
shared_size *= 4 # sizeof(float32) shared_size *= 4 # sizeof(float32)
self._kernel_options = { self._kernel_options = {
"file_name": get_cuda_srcfile("backproj.cu"), "file_name": get_cuda_srcfile(self.cuda_fname),
"kernel_name": "backproj", "kernel_name": self.cuda_kernel_name,
"kernel_signature": self._get_kernel_signature(), "kernel_signature": self._get_kernel_signature(),
"texture_name": tex_name, "texture_name": tex_name,
"sourcemodule_options": sourcemodule_options, "sourcemodule_options": sourcemodule_options,
...@@ -324,3 +328,40 @@ class Backprojector(CudaProcessing): ...@@ -324,3 +328,40 @@ class Backprojector(CudaProcessing):
fbp = filtered_backprojection # shorthand fbp = filtered_backprojection # shorthand
class PolarBackprojector(Backprojector):
"""
Cuda Backprojector with output in polar coordinates.
"""
cuda_fname = "backproj_polar.cu"
cuda_kernel_name = "backproj_polar"
# patch parent method: force slice_shape to (n_angles, n_x)
def _set_angles(self, angles, n_angles):
Backprojector._set_angles(self, angles, n_angles)
self.slice_shape = (self.n_angles, self.n_x)
# patch parent method:
def _set_slice_roi(self, slice_roi):
if slice_roi is not None:
raise ValueError("slice_roi is not supported with this class")
Backprojector._set_slice_roi(self, slice_roi)
# patch parent method: don't do the 4X compute-workload optimization for this kernel
def _get_kernel_options(self):
Backprojector._get_kernel_options(self)
block = self._kernel_options["block"]
self._kernel_options["grid"] = (
updiv(self.n_x, block[0]),
updiv(self.n_y, block[1])
)
# patch parent method: update kernel args
def _compile_kernels(self):
n_y = self.n_y
self.n_y = self.n_angles
Backprojector._compile_kernels(self)
self.n_y = n_y
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