diff --git a/include/rtkCudaForwardProjectionImageFilter.hcu b/include/rtkCudaForwardProjectionImageFilter.hcu index 52177f817..86eb29895 100644 --- a/include/rtkCudaForwardProjectionImageFilter.hcu +++ b/include/rtkCudaForwardProjectionImageFilter.hcu @@ -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], diff --git a/include/rtkCudaForwardProjectionImageFilter.hxx b/include/rtkCudaForwardProjectionImageFilter.hxx index 6220572f6..262729781 100644 --- a/include/rtkCudaForwardProjectionImageFilter.hxx +++ b/include/rtkCudaForwardProjectionImageFilter.hxx @@ -41,12 +41,6 @@ template void CudaForwardProjectionImageFilter::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); @@ -55,8 +49,6 @@ CudaForwardProjectionImageFilter::GPUGenerateData() this->GetOutput()->GetBufferedRegion().GetSize(1) * itk::NumericTraits::GetLength(); - itk::Vector 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, @@ -116,11 +108,17 @@ CudaForwardProjectionImageFilter::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++) @@ -158,12 +156,25 @@ CudaForwardProjectionImageFilter::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 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(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; @@ -185,6 +196,8 @@ CudaForwardProjectionImageFilter::GPUGenerateData() pvol, m_StepSize, (&(source_positions[3 * i])), + (&(source_to_pixels[3 * i])), + isParallel, radiusCylindricalDetector, boxMin, boxMax, @@ -195,6 +208,7 @@ CudaForwardProjectionImageFilter::GPUGenerateData() delete[] translatedProjectionIndexTransformMatrices; delete[] translatedVolumeTransformMatrices; delete[] source_positions; + delete[] source_to_pixels; } } // end namespace rtk diff --git a/include/rtkCudaRayCastBackProjectionImageFilter.hcu b/include/rtkCudaRayCastBackProjectionImageFilter.hcu index 8d6dbdf5f..94748d1b6 100644 --- a/include/rtkCudaRayCastBackProjectionImageFilter.hcu +++ b/include/rtkCudaRayCastBackProjectionImageFilter.hcu @@ -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 diff --git a/include/rtkCudaWarpForwardProjectionImageFilter.hcu b/include/rtkCudaWarpForwardProjectionImageFilter.hcu index 9745c70b8..5409489cb 100644 --- a/include/rtkCudaWarpForwardProjectionImageFilter.hcu +++ b/include/rtkCudaWarpForwardProjectionImageFilter.hcu @@ -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], diff --git a/src/rtkCudaForwardProjectionImageFilter.cu b/src/rtkCudaForwardProjectionImageFilter.cu index f2ea32c1b..555735fef 100644 --- a/src/rtkCudaForwardProjectionImageFilter.cu +++ b/src/rtkCudaForwardProjectionImageFilter.cu @@ -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; //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ @@ -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; @@ -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], @@ -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, diff --git a/src/rtkCudaRayCastBackProjectionImageFilter.cu b/src/rtkCudaRayCastBackProjectionImageFilter.cu index 34a28aa18..a993f956f 100644 --- a/src/rtkCudaRayCastBackProjectionImageFilter.cu +++ b/src/rtkCudaRayCastBackProjectionImageFilter.cu @@ -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) @@ -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 @@ -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)); @@ -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, diff --git a/src/rtkCudaRayCastBackProjectionImageFilter.cxx b/src/rtkCudaRayCastBackProjectionImageFilter.cxx index 77c35625f..1233498c5 100644 --- a/src/rtkCudaRayCastBackProjectionImageFilter.cxx +++ b/src/rtkCudaRayCastBackProjectionImageFilter.cxx @@ -31,12 +31,6 @@ namespace rtk void CudaRayCastBackProjectionImageFilter ::GPUGenerateData() { - if (!this->GetGeometry()->GetSourceToDetectorDistances().empty() && - this->GetGeometry()->GetSourceToDetectorDistances()[0] == 0) - { - itkGenericExceptionMacro(<< "Parallel geometry is not handled by CUDA forward projector."); - } - const auto * geometry = dynamic_cast(this->GetGeometry()); if (!geometry) { @@ -48,8 +42,6 @@ CudaRayCastBackProjectionImageFilter ::GPUGenerateData() const unsigned int nPixelsPerProj = this->GetInput(1)->GetBufferedRegion().GetSize(0) * this->GetInput(1)->GetBufferedRegion().GetSize(1); - itk::Vector 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, @@ -105,11 +97,17 @@ CudaRayCastBackProjectionImageFilter ::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++) @@ -147,12 +145,25 @@ CudaRayCastBackProjectionImageFilter ::GPUGenerateData() static_cast(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 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(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; @@ -170,7 +181,9 @@ CudaRayCastBackProjectionImageFilter ::GPUGenerateData() pout, pproj + nPixelsPerProj * projectionOffset, m_StepSize, - reinterpret_cast(&(source_positions[3 * i])), + &(source_positions[3 * i]), + &(source_to_pixels[3 * i]), + isParallel, radiusCylindricalDetector, boxMin, boxMax, @@ -183,6 +196,7 @@ CudaRayCastBackProjectionImageFilter ::GPUGenerateData() delete[] translatedProjectionIndexTransformMatrices; delete[] translatedVolumeTransformMatrices; delete[] source_positions; + delete[] source_to_pixels; } } // end namespace rtk diff --git a/src/rtkCudaWarpForwardProjectionImageFilter.cu b/src/rtkCudaWarpForwardProjectionImageFilter.cu index 4b73e135f..e0e113833 100644 --- a/src/rtkCudaWarpForwardProjectionImageFilter.cu +++ b/src/rtkCudaWarpForwardProjectionImageFilter.cu @@ -45,8 +45,10 @@ __constant__ float3 c_boxMax; __constant__ float3 c_spacing; __constant__ int3 c_volSize; __constant__ float c_tStep; -__constant__ float c_matrices[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_matrices[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; __constant__ float c_IndexInputToPPInputMatrix[12]; __constant__ float c_IndexInputToIndexDVFMatrix[12]; @@ -80,12 +82,20 @@ kernel_warped_forwardProject(float * dev_proj_in, 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]); - pixelPos = matrix_multiply(make_float3(i, j, 0), &(c_matrices[12 * proj])); - ray.d = pixelPos - ray.o; + if (c_isParallel) + { + // Start outside the volume and follow the shared parallel ray direction. + ray.d = make_float3(c_sourceToPixel[3 * proj], c_sourceToPixel[3 * proj + 1], c_sourceToPixel[3 * proj + 2]); + ray.o = pixelPos - ray.d; + } + else + { + // Setting ray origin + ray.o = make_float3(c_sourcePos[3 * proj], c_sourcePos[3 * proj + 1], c_sourcePos[3 * proj + 2]); + ray.d = pixelPos - ray.o; + } ray.d = ray.d / sqrtf(dot(ray.d, ray.d)); // Detect intersection with box @@ -163,6 +173,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], @@ -178,9 +190,12 @@ CUDA_warp_forward_project(int projSize[3], cudaMemcpyToSymbol(c_spacing, spacing, sizeof(float3)); cudaMemcpyToSymbol(c_volSize, volSize, sizeof(int3)); cudaMemcpyToSymbol(c_tStep, &t_step, 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_matrices, &(matrices[0]), 12 * sizeof(float) * projSize[2]); diff --git a/src/rtkCudaWarpForwardProjectionImageFilter.cxx b/src/rtkCudaWarpForwardProjectionImageFilter.cxx index 1756ec964..a66e71d5e 100644 --- a/src/rtkCudaWarpForwardProjectionImageFilter.cxx +++ b/src/rtkCudaWarpForwardProjectionImageFilter.cxx @@ -127,12 +127,6 @@ CudaWarpForwardProjectionImageFilter ::GenerateInputRequestedRegion() void CudaWarpForwardProjectionImageFilter ::GPUGenerateData() { - if (!this->GetGeometry()->GetSourceToDetectorDistances().empty() && - this->GetGeometry()->GetSourceToDetectorDistances()[0] == 0) - { - itkGenericExceptionMacro(<< "Parallel geometry is not handled by CUDA forward projector."); - } - itk::Matrix matrixIdxInputVol; itk::Matrix indexInputToPPInputMatrix; itk::Matrix indexInputToIndexDVFMatrix; @@ -144,14 +138,17 @@ CudaWarpForwardProjectionImageFilter ::GPUGenerateData() } const Superclass::GeometryType * geometry = this->GetGeometry(); - const unsigned int Dimension = InputImageType::ImageDimension; + bool isParallel = geometry->GetSourceToDetectorDistances()[0] == 0.; + if (isParallel && geometry->GetRadiusCylindricalDetector() != 0.) + { + itkGenericExceptionMacro(<< "Parallel geometry assumes a flat panel detector."); + } + const unsigned int Dimension = InputImageType::ImageDimension; const unsigned int iFirstProj = this->GetInputProjectionStack()->GetRequestedRegion().GetIndex(Dimension - 1); const unsigned int nProj = this->GetInputProjectionStack()->GetRequestedRegion().GetSize(Dimension - 1); const unsigned int nPixelsPerProj = this->GetOutput()->GetBufferedRegion().GetSize(0) * this->GetOutput()->GetBufferedRegion().GetSize(1); - itk::Vector 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, @@ -232,8 +229,9 @@ CudaWarpForwardProjectionImageFilter ::GPUGenerateData() } // Compute matrices to transform projection index to volume index, one per projection - auto * matrices = new float[12 * nProj]; - auto * source_positions = new float[4 * nProj]; + auto * matrices = new float[12 * nProj]{}; + auto * source_positions = new float[3 * nProj]{}; + auto * source_to_pixels = new float[3 * nProj]{}; // Go over each projection for (unsigned int iProj = iFirstProj; iProj < iFirstProj + nProj; iProj++) @@ -246,12 +244,25 @@ CudaWarpForwardProjectionImageFilter ::GPUGenerateData() for (int k = 0; k < 4; k++) matrices[(j + 3 * (iProj - iFirstProj)) * 4 + k] = static_cast(d_matrix[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 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(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; @@ -270,6 +281,8 @@ CudaWarpForwardProjectionImageFilter ::GPUGenerateData() pvol, m_StepSize, &(source_positions[3 * i]), + &(source_to_pixels[3 * i]), + isParallel, boxMin, boxMax, spacing, @@ -281,6 +294,7 @@ CudaWarpForwardProjectionImageFilter ::GPUGenerateData() delete[] matrices; delete[] source_positions; + delete[] source_to_pixels; } } // end namespace rtk diff --git a/test/rtkcudaraycastadjointoperatorstest.cxx b/test/rtkcudaraycastadjointoperatorstest.cxx index 110f30349..e63314889 100644 --- a/test/rtkcudaraycastadjointoperatorstest.cxx +++ b/test/rtkcudaraycastadjointoperatorstest.cxx @@ -136,5 +136,25 @@ rtkcudaraycastadjointoperatorstest(int, char *[]) randomVolumeSource->GetOutput(), bp->GetOutput(), randomProjectionsSource->GetOutput(), fw->GetOutput()); std::cout << "\n\nTest PASSED! " << std::endl; + // Start over with parallel geometry + geometry = rtk::ThreeDCircularProjectionGeometry::New(); + for (unsigned int noProj = 0; noProj < NumberOfProjectionImages; noProj++) + geometry->AddProjection(600., 0., noProj * 360. / NumberOfProjectionImages); + + std::cout << "\n\n****** CUDA ray cast Forward projector, parallel geometry ******" << std::endl; + + fw->SetGeometry(geometry); + TRY_AND_EXIT_ON_ITK_EXCEPTION(fw->Update()); + + std::cout << "\n\n****** CUDA ray cast Back projector, parallel geometry ******" << std::endl; + + bp->SetGeometry(geometry.GetPointer()); + + TRY_AND_EXIT_ON_ITK_EXCEPTION(bp->Update()); + + CheckScalarProducts( + randomVolumeSource->GetOutput(), bp->GetOutput(), randomProjectionsSource->GetOutput(), fw->GetOutput()); + std::cout << "\n\nTest PASSED! " << std::endl; + return EXIT_SUCCESS; } diff --git a/test/rtkforwardprojectiontest.cxx b/test/rtkforwardprojectiontest.cxx index 4d56c7b13..fbf9f8efc 100644 --- a/test/rtkforwardprojectiontest.cxx +++ b/test/rtkforwardprojectiontest.cxx @@ -206,5 +206,21 @@ rtkforwardprojectiontest(int, char *[]) CheckImageQuality(stream->GetOutput(), slp->GetOutput(), 1.28, 44, 255.0); std::cout << "\n\nTest PASSED! " << std::endl; + std::cout << "\n\n****** Case 6: Shepp-Logan, parallel geometry, flat panel detector ******" << std::endl; + + // Parallel geometry: SDD = 0 + geometry = GeometryType::New(); + for (unsigned int i = 0; i < NumberOfProjectionImages; i++) + geometry->AddProjection(1000., 0., i * 8.); + + slp->SetGeometry(geometry); + slp->Update(); + + jfp->SetGeometry(geometry); + stream->Update(); + + CheckImageQuality(stream->GetOutput(), slp->GetOutput(), 1.28, 40, 255.0); + std::cout << "\n\nTest PASSED! " << std::endl; + return EXIT_SUCCESS; }