Skip to content

ENH: Support CUDA ray cast projectors with parallel geometry #966

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
axel-grc wants to merge 1 commit into
base: main
Choose a base branch
from
Open

ENH: Support CUDA ray cast projectors with parallel geometry #966

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions include/rtkCudaForwardProjectionImageFilter.hcu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@ CUDA_forward_project(int projSize[3],
float * dev_vol,
float t_step,
float * source_positions,
float * source_to_pixels,
bool is_parallel,
float radiusCylindricalDetector,
float box_min[3],
float box_max[3],
Expand Down
48 changes: 31 additions & 17 deletions include/rtkCudaForwardProjectionImageFilter.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,12 +41,6 @@ template <class TInputImage, class TOutputImage>
void
CudaForwardProjectionImageFilter<TInputImage, TOutputImage>::GPUGenerateData()
{
if (this->GetGeometry()->GetSourceToDetectorDistances().size() &&
this->GetGeometry()->GetSourceToDetectorDistances()[0] == 0)
{
itkGenericExceptionMacro(<< "Parallel geometry is not handled by CUDA forward projector.");
}

const typename Superclass::GeometryType * geometry = this->GetGeometry();
const unsigned int Dimension = TInputImage::ImageDimension;
const unsigned int iFirstProj = this->GetInput(0)->GetRequestedRegion().GetIndex(Dimension - 1);
Expand All @@ -55,8 +49,6 @@ CudaForwardProjectionImageFilter<TInputImage, TOutputImage>::GPUGenerateData()
this->GetOutput()->GetBufferedRegion().GetSize(1) *
itk::NumericTraits<typename TInputImage::PixelType>::GetLength();

itk::Vector<double, 4> source_position;

// Setting BoxMin and BoxMax
// SR: we are using cuda textures where the pixel definition is not center but corner.
// Therefore, we set the box limits from index to index+size instead of, for ITK,
Expand Down Expand Up @@ -116,11 +108,17 @@ CudaForwardProjectionImageFilter<TInputImage, TOutputImage>::GPUGenerateData()
}

// Compute matrices to transform projection index to volume index, one per projection
auto * translatedProjectionIndexTransformMatrices = new float[12 * nProj];
auto * translatedVolumeTransformMatrices = new float[12 * nProj];
auto * source_positions = new float[4 * nProj];
auto * translatedProjectionIndexTransformMatrices = new float[12 * nProj]{};
auto * translatedVolumeTransformMatrices = new float[12 * nProj]{};
auto * source_positions = new float[3 * nProj]{};
auto * source_to_pixels = new float[3 * nProj]{};

float radiusCylindricalDetector = geometry->GetRadiusCylindricalDetector();
bool isParallel = geometry->GetSourceToDetectorDistances()[0] == 0.;
if (isParallel && radiusCylindricalDetector != 0.)
{
itkGenericExceptionMacro(<< "Parallel geometry assumes a flat panel detector.");
}

// Go over each projection
for (unsigned int iProj = iFirstProj; iProj < iFirstProj + nProj; iProj++)
Expand Down Expand Up @@ -158,12 +156,25 @@ CudaForwardProjectionImageFilter<TInputImage, TOutputImage>::GPUGenerateData()
(float)translatedVolumeTransformMatrix[j][k];
}

// Compute source position in volume indices
source_position = volPPToIndex * geometry->GetSourcePosition(iProj);

// Copy it into a single large array
for (unsigned int d = 0; d < 3; d++)
source_positions[(iProj - iFirstProj) * 3 + d] = source_position[d]; // Ignore the 4th component
if (isParallel)
{
// Parallel projections have no finite source; store the common ray direction instead.
itk::Vector<double, 3> sourceToPixel;
for (unsigned int d = 0; d < 3; d++)
sourceToPixel[d] = geometry->GetRotationMatrices()[iProj][2][d];
sourceToPixel *= -2. * geometry->GetSourceToIsocenterDistances()[iProj];
auto sourceToPixelVnl = volPPToIndex.GetVnlMatrix().extract(3, 3) * sourceToPixel.GetVnlVector();

for (unsigned int d = 0; d < 3; d++)
source_to_pixels[(iProj - iFirstProj) * 3 + d] = static_cast<float>(sourceToPixelVnl[d]);
}
else
{
// Cone-beam projections use the finite source position in volume indices.
auto sourcePosition = volPPToIndex * geometry->GetSourcePosition(iProj);
for (unsigned int d = 0; d < 3; d++)
source_positions[(iProj - iFirstProj) * 3 + d] = sourcePosition[d];
}
}

int projectionOffset = 0;
Expand All @@ -185,6 +196,8 @@ CudaForwardProjectionImageFilter<TInputImage, TOutputImage>::GPUGenerateData()
pvol,
m_StepSize,
(&(source_positions[3 * i])),
(&(source_to_pixels[3 * i])),
isParallel,
radiusCylindricalDetector,
boxMin,
boxMax,
Expand All @@ -195,6 +208,7 @@ CudaForwardProjectionImageFilter<TInputImage, TOutputImage>::GPUGenerateData()
delete[] translatedProjectionIndexTransformMatrices;
delete[] translatedVolumeTransformMatrices;
delete[] source_positions;
delete[] source_to_pixels;
}

} // end namespace rtk
Expand Down
28 changes: 15 additions & 13 deletions include/rtkCudaRayCastBackProjectionImageFilter.hcu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,18 +20,20 @@
#define rtkCudaRayCastBackProjectionImageFilter_hcu

void
CUDA_ray_cast_back_project(int projSize[3],
int volSize[3],
float * translatedProjectionIndexTransformMatrices,
float * translatedVolumeTransformMatrices,
float * dev_vol_in,
float * dev_vol_out,
float * dev_proj,
float t_step,
double * source_positions,
float radiusCylindricalDetector,
float box_min[3],
float box_max[3],
float spacing[3]);
CUDA_ray_cast_back_project(int projSize[3],
int volSize[3],
float * translatedProjectionIndexTransformMatrices,
float * translatedVolumeTransformMatrices,
float * dev_vol_in,
float * dev_vol_out,
float * dev_proj,
float t_step,
float * source_positions,
float * source_to_pixels,
bool is_parallel,
float radiusCylindricalDetector,
float box_min[3],
float box_max[3],
float spacing[3]);

#endif
2 changes: 2 additions & 0 deletions include/rtkCudaWarpForwardProjectionImageFilter.hcu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@ CUDA_warp_forward_project(int projSize[3],
float * dev_vol,
float t_step,
float * source_positions,
float * source_to_pixels,
bool is_parallel,
float box_min[3],
float box_max[3],
float spacing[3],
Expand Down
41 changes: 29 additions & 12 deletions src/rtkCudaForwardProjectionImageFilter.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ __constant__ float
__constant__ float
c_translatedVolumeTransformMatrices[SLAB_SIZE * 12]; // Can process stacks of at most SLAB_SIZE projections
__constant__ float c_sourcePos[SLAB_SIZE * 3]; // Can process stacks of at most SLAB_SIZE projections
__constant__ float c_sourceToPixel[SLAB_SIZE * 3]; // Can process stacks of at most SLAB_SIZE projections
__constant__ int c_isParallel;


//_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
Expand All @@ -77,24 +79,34 @@ kernel_forwardProject(float * dev_proj_in, float * dev_proj_out, cudaTextureObje

for (unsigned int proj = 0; proj < c_projSize.z; proj++)
{
// Setting ray origin
ray.o = make_float3(c_sourcePos[3 * proj], c_sourcePos[3 * proj + 1], c_sourcePos[3 * proj + 2]);

if (c_radius == 0)
if (c_isParallel)
{
// Start outside the volume and follow the shared parallel ray direction.
pixelPos = matrix_multiply(make_float3(i, j, 0), &(c_translatedProjectionIndexTransformMatrices[12 * proj]));
ray.d = make_float3(c_sourceToPixel[3 * proj], c_sourceToPixel[3 * proj + 1], c_sourceToPixel[3 * proj + 2]);
ray.o = pixelPos - ray.d;
}
else
{
float3 posProj;
posProj = matrix_multiply(make_float3(i, j, 0), &(c_translatedProjectionIndexTransformMatrices[12 * proj]));
double a = posProj.x / c_radius;
posProj.x = sin(a) * c_radius;
posProj.z += (1. - cos(a)) * c_radius;
pixelPos = matrix_multiply(posProj, &(c_translatedVolumeTransformMatrices[12 * proj]));
}
// Setting ray origin
ray.o = make_float3(c_sourcePos[3 * proj], c_sourcePos[3 * proj + 1], c_sourcePos[3 * proj + 2]);

if (c_radius == 0)
{
pixelPos = matrix_multiply(make_float3(i, j, 0), &(c_translatedProjectionIndexTransformMatrices[12 * proj]));
}
else
{
float3 posProj;
posProj = matrix_multiply(make_float3(i, j, 0), &(c_translatedProjectionIndexTransformMatrices[12 * proj]));
double a = posProj.x / c_radius;
posProj.x = sin(a) * c_radius;
posProj.z += (1. - cos(a)) * c_radius;
pixelPos = matrix_multiply(posProj, &(c_translatedVolumeTransformMatrices[12 * proj]));
}

ray.d = pixelPos - ray.o;
ray.d = pixelPos - ray.o;
}
ray.d = ray.d / sqrtf(dot(ray.d, ray.d));

int projOffset = numThread + proj * c_projSize.x * c_projSize.y;
Expand Down Expand Up @@ -170,6 +182,8 @@ CUDA_forward_project(int projSize[3],
float * dev_vol,
float t_step,
float * source_positions,
float * source_to_pixels,
bool is_parallel,
float radiusCylindricalDetector,
float box_min[3],
float box_max[3],
Expand All @@ -184,12 +198,15 @@ CUDA_forward_project(int projSize[3],
cudaMemcpyToSymbol(c_volSize, volSize, sizeof(int3));
cudaMemcpyToSymbol(c_tStep, &t_step, sizeof(float));
cudaMemcpyToSymbol(c_radius, &radiusCylindricalDetector, sizeof(float));
int isParallel = is_parallel ? 1 : 0;
cudaMemcpyToSymbol(c_isParallel, &isParallel, sizeof(int));

dim3 dimBlock = dim3(16, 16, 1);
dim3 dimGrid = dim3(iDivUp(projSize[0], dimBlock.x), iDivUp(projSize[1], dimBlock.x));

// Copy the source position matrix into a float3 in constant memory
cudaMemcpyToSymbol(c_sourcePos, &(source_positions[0]), 3 * sizeof(float) * projSize[2]);
cudaMemcpyToSymbol(c_sourceToPixel, &(source_to_pixels[0]), 3 * sizeof(float) * projSize[2]);

// Copy the projection matrices into constant memory
cudaMemcpyToSymbol(c_translatedProjectionIndexTransformMatrices,
Expand Down
67 changes: 42 additions & 25 deletions src/rtkCudaRayCastBackProjectionImageFilter.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ __constant__ float
__constant__ float
c_translatedVolumeTransformMatrices[SLAB_SIZE * 12]; // Can process stacks of at most SLAB_SIZE projections
__constant__ float c_sourcePos[SLAB_SIZE * 3]; // Can process stacks of at most SLAB_SIZE projections
__constant__ float c_sourceToPixel[SLAB_SIZE * 3]; // Can process stacks of at most SLAB_SIZE projections
__constant__ int c_isParallel;

__global__ void
kernel_ray_cast_back_project(float * dev_vol_out, float * dev_proj)
Expand All @@ -68,24 +70,34 @@ kernel_ray_cast_back_project(float * dev_vol_out, float * dev_proj)
float3 pixelPos;
float tnear, tfar;

// Setting ray origin
ray.o = make_float3(c_sourcePos[3 * proj], c_sourcePos[3 * proj + 1], c_sourcePos[3 * proj + 2]);

if (c_radius == 0)
if (c_isParallel)
{
// Start outside the volume and follow the shared parallel ray direction.
pixelPos = matrix_multiply(make_float3(i, j, 0), &(c_translatedProjectionIndexTransformMatrices[12 * proj]));
ray.d = make_float3(c_sourceToPixel[3 * proj], c_sourceToPixel[3 * proj + 1], c_sourceToPixel[3 * proj + 2]);
ray.o = pixelPos - ray.d;
}
else
{
float3 posProj;
posProj = matrix_multiply(make_float3(i, j, 0), &(c_translatedProjectionIndexTransformMatrices[12 * proj]));
double a = posProj.x / c_radius;
posProj.x = sin(a) * c_radius;
posProj.z += (1. - cos(a)) * c_radius;
pixelPos = matrix_multiply(posProj, &(c_translatedVolumeTransformMatrices[12 * proj]));
}
// Setting ray origin
ray.o = make_float3(c_sourcePos[3 * proj], c_sourcePos[3 * proj + 1], c_sourcePos[3 * proj + 2]);

if (c_radius == 0)
{
pixelPos = matrix_multiply(make_float3(i, j, 0), &(c_translatedProjectionIndexTransformMatrices[12 * proj]));
}
else
{
float3 posProj;
posProj = matrix_multiply(make_float3(i, j, 0), &(c_translatedProjectionIndexTransformMatrices[12 * proj]));
double a = posProj.x / c_radius;
posProj.x = sin(a) * c_radius;
posProj.z += (1. - cos(a)) * c_radius;
pixelPos = matrix_multiply(posProj, &(c_translatedVolumeTransformMatrices[12 * proj]));
}

ray.d = pixelPos - ray.o;
ray.d = pixelPos - ray.o;
}
ray.d = ray.d / sqrtf(dot(ray.d, ray.d));

// Detect intersection with box
Expand Down Expand Up @@ -184,19 +196,21 @@ kernel_ray_cast_back_project(float * dev_vol_out, float * dev_proj)
}

void
CUDA_ray_cast_back_project(int projSize[3],
int volSize[3],
float * translatedProjectionIndexTransformMatrices,
float * translatedVolumeTransformMatrices,
float * dev_vol_in,
float * dev_vol_out,
float * dev_proj,
float t_step,
double * source_positions,
float radiusCylindricalDetector,
float box_min[3],
float box_max[3],
float spacing[3])
CUDA_ray_cast_back_project(int projSize[3],
int volSize[3],
float * translatedProjectionIndexTransformMatrices,
float * translatedVolumeTransformMatrices,
float * dev_vol_in,
float * dev_vol_out,
float * dev_proj,
float t_step,
float * source_positions,
float * source_to_pixels,
bool is_parallel,
float radiusCylindricalDetector,
float box_min[3],
float box_max[3],
float spacing[3])
{
// Constant memory
cudaMemcpyToSymbol(c_projSize, projSize, sizeof(int3));
Expand All @@ -206,9 +220,12 @@ CUDA_ray_cast_back_project(int projSize[3],
cudaMemcpyToSymbol(c_volSize, volSize, sizeof(int3));
cudaMemcpyToSymbol(c_tStep, &t_step, sizeof(float));
cudaMemcpyToSymbol(c_radius, &radiusCylindricalDetector, sizeof(float));
int isParallel = is_parallel ? 1 : 0;
cudaMemcpyToSymbol(c_isParallel, &isParallel, sizeof(int));

// Copy the source position matrix into a float3 in constant memory
cudaMemcpyToSymbol(c_sourcePos, &(source_positions[0]), 3 * sizeof(float) * projSize[2]);
cudaMemcpyToSymbol(c_sourceToPixel, &(source_to_pixels[0]), 3 * sizeof(float) * projSize[2]);

// Copy the projection matrices into constant memory
cudaMemcpyToSymbol(c_translatedProjectionIndexTransformMatrices,
Expand Down
Loading
Loading