From 422b327402a5d6ab8d467da5b4efd9e7b3ecb7f1 Mon Sep 17 00:00:00 2001 From: Axel Garcia Date: Thu, 8 Oct 2026 13:24:36 +0200 Subject: [PATCH] BUG: Use 64-bit index/offset types in CUDA filters to avoid overflow --- include/rtkCudaAverageOutOfROIImageFilter.hcu | 4 +- include/rtkCudaBackProjectionImageFilter.hcu | 20 ++--- include/rtkCudaBackProjectionImageFilter.hxx | 4 +- include/rtkCudaConstantVolumeSeriesSource.hcu | 4 +- include/rtkCudaConstantVolumeSource.hcu | 4 +- .../rtkCudaCyclicDeformationImageFilter.hcu | 16 ++-- .../rtkCudaFDKBackProjectionImageFilter.hcu | 12 +-- .../rtkCudaForwardProjectionImageFilter.hcu | 4 +- .../rtkCudaForwardProjectionImageFilter.hxx | 18 ++-- include/rtkCudaForwardWarpImageFilter.hcu | 24 ++--- include/rtkCudaInterpolateImageFilter.hcu | 8 +- include/rtkCudaSplatImageFilter.hcu | 8 +- include/rtkCudaUtilities.hcu | 17 +++- .../rtkCudaWarpBackProjectionImageFilter.hcu | 28 +++--- ...tkCudaWarpForwardProjectionImageFilter.hcu | 32 +++---- include/rtkCudaWarpImageFilter.hcu | 20 ++--- src/rtkCudaAverageOutOfROIImageFilter.cu | 21 +++-- src/rtkCudaAverageOutOfROIImageFilter.cxx | 11 ++- src/rtkCudaBackProjectionImageFilter.cu | 42 ++++----- src/rtkCudaConjugateGradientImageFilter.cu | 56 +++++++++++- src/rtkCudaConstantVolumeSeriesSource.cu | 17 ++-- src/rtkCudaConstantVolumeSeriesSource.cxx | 11 ++- src/rtkCudaConstantVolumeSource.cu | 15 ++-- src/rtkCudaConstantVolumeSource.cxx | 10 +-- src/rtkCudaCropImageFilter.cu | 15 ++-- src/rtkCudaCyclicDeformationImageFilter.cu | 16 ++-- src/rtkCudaCyclicDeformationImageFilter.cxx | 20 ++--- src/rtkCudaDisplacedDetectorImageFilter.cu | 12 +-- src/rtkCudaFDKBackProjectionImageFilter.cu | 30 +++---- src/rtkCudaFDKBackProjectionImageFilter.cxx | 4 +- src/rtkCudaFDKWeightProjectionFilter.cu | 31 ++++--- ...udaFFTProjectionsConvolutionImageFilter.cu | 44 +++++----- src/rtkCudaFirstOrderKernels.cu | 30 ++++--- src/rtkCudaForwardProjectionImageFilter.cu | 56 ++++++------ src/rtkCudaForwardWarpImageFilter.cu | 88 +++++++++---------- src/rtkCudaForwardWarpImageFilter.cxx | 6 +- src/rtkCudaInterpolateImageFilter.cu | 10 ++- src/rtkCudaInterpolateImageFilter.cxx | 10 +-- src/rtkCudaLagCorrectionImageFilter.cu | 26 +++--- ...CudaLastDimensionTVDenoisingImageFilter.cu | 5 +- src/rtkCudaParkerShortScanImageFilter.cu | 27 +++--- ...CudaPolynomialGainCorrectionImageFilter.cu | 34 ++++--- ...rtkCudaRayCastBackProjectionImageFilter.cu | 48 ++++++---- ...tkCudaRayCastBackProjectionImageFilter.cxx | 10 +-- src/rtkCudaSplatImageFilter.cu | 10 ++- src/rtkCudaSplatImageFilter.cxx | 10 +-- ...aTotalVariationDenoisingBPDQImageFilter.cu | 38 ++++---- src/rtkCudaUtilities.cu | 4 +- src/rtkCudaWarpBackProjectionImageFilter.cu | 64 +++++++------- src/rtkCudaWarpBackProjectionImageFilter.cxx | 6 +- ...rtkCudaWarpForwardProjectionImageFilter.cu | 69 +++++++-------- ...tkCudaWarpForwardProjectionImageFilter.cxx | 14 +-- src/rtkCudaWarpImageFilter.cu | 28 +++--- src/rtkCudaWarpImageFilter.cxx | 6 +- ...rtkCudaWeidingerForwardModelImageFilter.cu | 12 +-- 55 files changed, 632 insertions(+), 557 deletions(-) diff --git a/include/rtkCudaAverageOutOfROIImageFilter.hcu b/include/rtkCudaAverageOutOfROIImageFilter.hcu index d4d770044..9c413ffd8 100644 --- a/include/rtkCudaAverageOutOfROIImageFilter.hcu +++ b/include/rtkCudaAverageOutOfROIImageFilter.hcu @@ -19,9 +19,9 @@ #ifndef rtkCudaAverageOutOfROIImageFilter_hcu #define rtkCudaAverageOutOfROIImageFilter_hcu -#include +#include "rtkCudaUtilities.hcu" void -CUDA_average_out_of_ROI(int size[4], float * input, float * output, float * roi); +CUDA_average_out_of_ROI(itk::SizeValueType size[4], float * input, float * output, float * roi); #endif diff --git a/include/rtkCudaBackProjectionImageFilter.hcu b/include/rtkCudaBackProjectionImageFilter.hcu index c7f0499fc..1f03f09fe 100644 --- a/include/rtkCudaBackProjectionImageFilter.hcu +++ b/include/rtkCudaBackProjectionImageFilter.hcu @@ -22,15 +22,15 @@ #include "RTKExport.h" void RTK_EXPORT -CUDA_back_project(int projSize[3], - int volSize[3], - float * matrices, - float * volIndexToProjPPs, - float * projPPToProjIndex, - float * dev_vol_in, - float * dev_vol_out, - float * dev_img, - double radiusCylindricalDetector, - unsigned int vectorLength); +CUDA_back_project(itk::SizeValueType projSize[3], + itk::SizeValueType volSize[3], + float * matrices, + float * volIndexToProjPPs, + float * projPPToProjIndex, + float * dev_vol_in, + float * dev_vol_out, + float * dev_img, + double radiusCylindricalDetector, + unsigned int vectorLength); #endif diff --git a/include/rtkCudaBackProjectionImageFilter.hxx b/include/rtkCudaBackProjectionImageFilter.hxx index c7193d406..653f1864d 100644 --- a/include/rtkCudaBackProjectionImageFilter.hxx +++ b/include/rtkCudaBackProjectionImageFilter.hxx @@ -67,11 +67,11 @@ CudaBackProjectionImageFilter::GPUGenerateData() } // Cuda convenient format for dimensions - int projectionSize[3]; + itk::SizeValueType projectionSize[3]; projectionSize[0] = this->GetInput(1)->GetBufferedRegion().GetSize()[0]; projectionSize[1] = this->GetInput(1)->GetBufferedRegion().GetSize()[1]; - int volumeSize[3]; + itk::SizeValueType volumeSize[3]; volumeSize[0] = this->GetOutput()->GetBufferedRegion().GetSize()[0]; volumeSize[1] = this->GetOutput()->GetBufferedRegion().GetSize()[1]; volumeSize[2] = this->GetOutput()->GetBufferedRegion().GetSize()[2]; diff --git a/include/rtkCudaConstantVolumeSeriesSource.hcu b/include/rtkCudaConstantVolumeSeriesSource.hcu index 1206a506b..1a9ccc5d9 100644 --- a/include/rtkCudaConstantVolumeSeriesSource.hcu +++ b/include/rtkCudaConstantVolumeSeriesSource.hcu @@ -19,9 +19,9 @@ #ifndef rtkCudaConstantVolumeSeriesSource_hcu #define rtkCudaConstantVolumeSeriesSource_hcu -#include +#include "rtkCudaUtilities.hcu" void -CUDA_generate_constant_volume_series(int size[4], float * dev_out, float constantValue); +CUDA_generate_constant_volume_series(itk::SizeValueType size[4], float * dev_out, float constantValue); #endif diff --git a/include/rtkCudaConstantVolumeSource.hcu b/include/rtkCudaConstantVolumeSource.hcu index 327a20a2f..1ebd08c2b 100644 --- a/include/rtkCudaConstantVolumeSource.hcu +++ b/include/rtkCudaConstantVolumeSource.hcu @@ -19,9 +19,9 @@ #ifndef rtkCudaConstantVolumeSource_hcu #define rtkCudaConstantVolumeSource_hcu -#include +#include "rtkCudaUtilities.hcu" void -CUDA_generate_constant_volume(int size[3], float * dev_out, float constantValue); +CUDA_generate_constant_volume(itk::SizeValueType size[3], float * dev_out, float constantValue); #endif diff --git a/include/rtkCudaCyclicDeformationImageFilter.hcu b/include/rtkCudaCyclicDeformationImageFilter.hcu index c5ac08700..9216be0ac 100644 --- a/include/rtkCudaCyclicDeformationImageFilter.hcu +++ b/include/rtkCudaCyclicDeformationImageFilter.hcu @@ -19,13 +19,15 @@ #ifndef rtkCudaCyclicDeformationImageFilter_hcu #define rtkCudaCyclicDeformationImageFilter_hcu +#include "rtkCudaUtilities.hcu" + void -CUDA_linear_interpolate_along_fourth_dimension(unsigned int inputSize[4], - float * input, - float * output, - unsigned int frameInf, - unsigned int frameSup, - double weightInf, - double weightSup); +CUDA_linear_interpolate_along_fourth_dimension(itk::SizeValueType inputSize[4], + float * input, + float * output, + unsigned int frameInf, + unsigned int frameSup, + double weightInf, + double weightSup); #endif diff --git a/include/rtkCudaFDKBackProjectionImageFilter.hcu b/include/rtkCudaFDKBackProjectionImageFilter.hcu index 262425d4f..2772b6450 100644 --- a/include/rtkCudaFDKBackProjectionImageFilter.hcu +++ b/include/rtkCudaFDKBackProjectionImageFilter.hcu @@ -20,11 +20,11 @@ #define rtkCudaFDKBackProjectionImageFilter_hcu void -CUDA_reconstruct_conebeam(int img_dim[3], - int vol_dim[3], - float * matrices, - float * dev_vol_in, - float * dev_vol_out, - float * dev_img); +CUDA_reconstruct_conebeam(itk::SizeValueType img_dim[3], + itk::SizeValueType vol_dim[3], + float * matrices, + float * dev_vol_in, + float * dev_vol_out, + float * dev_img); #endif diff --git a/include/rtkCudaForwardProjectionImageFilter.hcu b/include/rtkCudaForwardProjectionImageFilter.hcu index 52177f817..a40360af1 100644 --- a/include/rtkCudaForwardProjectionImageFilter.hcu +++ b/include/rtkCudaForwardProjectionImageFilter.hcu @@ -22,8 +22,8 @@ #include "RTKExport.h" void RTK_EXPORT -CUDA_forward_project(int projSize[3], - int volSize[3], +CUDA_forward_project(itk::SizeValueType projSize[3], + itk::SizeValueType volSize[3], float * translatedProjectionIndexTransformMatrices, float * translatedVolumeTransformMatrices, float * dev_proj_in, diff --git a/include/rtkCudaForwardProjectionImageFilter.hxx b/include/rtkCudaForwardProjectionImageFilter.hxx index e01a50d20..5989e1ea4 100644 --- a/include/rtkCudaForwardProjectionImageFilter.hxx +++ b/include/rtkCudaForwardProjectionImageFilter.hxx @@ -49,11 +49,11 @@ CudaForwardProjectionImageFilter::GPUGenerateData() const typename Superclass::GeometryType * geometry = this->GetGeometry(); const unsigned int Dimension = TInputImage::ImageDimension; - const unsigned int iFirstProj = this->GetInput(0)->GetRequestedRegion().GetIndex(Dimension - 1); - const unsigned int nProj = this->GetInput(0)->GetRequestedRegion().GetSize(Dimension - 1); - const unsigned int nPixelsPerProj = this->GetOutput()->GetBufferedRegion().GetSize(0) * - this->GetOutput()->GetBufferedRegion().GetSize(1) * - itk::NumericTraits::GetLength(); + const unsigned int iFirstProj = this->GetInput(0)->GetRequestedRegion().GetIndex(Dimension - 1); + const unsigned int nProj = this->GetInput(0)->GetRequestedRegion().GetSize(Dimension - 1); + const itk::SizeValueType nPixelsPerProj = this->GetOutput()->GetBufferedRegion().GetSize(0) * + this->GetOutput()->GetBufferedRegion().GetSize(1) * + itk::NumericTraits::GetLength(); itk::Vector source_position; @@ -83,11 +83,11 @@ CudaForwardProjectionImageFilter::GPUGenerateData() } // Cuda convenient format for dimensions - int projectionSize[3]; + itk::SizeValueType projectionSize[3]; projectionSize[0] = this->GetOutput()->GetBufferedRegion().GetSize()[0]; projectionSize[1] = this->GetOutput()->GetBufferedRegion().GetSize()[1]; - int volumeSize[3]; + itk::SizeValueType volumeSize[3]; volumeSize[0] = this->GetInput(1)->GetBufferedRegion().GetSize()[0]; volumeSize[1] = this->GetInput(1)->GetBufferedRegion().GetSize()[1]; volumeSize[2] = this->GetInput(1)->GetBufferedRegion().GetSize()[2]; @@ -166,8 +166,8 @@ CudaForwardProjectionImageFilter::GPUGenerateData() source_positions[(iProj - iFirstProj) * 3 + d] = source_position[d]; // Ignore the 4th component } - int projectionOffset = 0; - const unsigned int vectorLength = itk::PixelTraits::Dimension; + itk::OffsetValueType projectionOffset = 0; + const unsigned int vectorLength = itk::PixelTraits::Dimension; for (unsigned int i = 0; i < nProj; i += SLAB_SIZE) { diff --git a/include/rtkCudaForwardWarpImageFilter.hcu b/include/rtkCudaForwardWarpImageFilter.hcu index 2f61d9a2c..5e277ac59 100644 --- a/include/rtkCudaForwardWarpImageFilter.hcu +++ b/include/rtkCudaForwardWarpImageFilter.hcu @@ -20,17 +20,17 @@ #define rtkCudaForwardWarpImageFilter_hcu void -CUDA_ForwardWarp(int input_vol_dim[3], - int input_dvf_dim[3], - int output_vol_dim[3], - float IndexInputToPPInputMatrix[12], - float IndexInputToIndexDVFMatrix[12], - float PPOutputToIndexOutputMatrix[12], - float * dev_input_vol, - float * dev_input_xdvf, - float * dev_input_ydvf, - float * dev_input_zdvf, - float * dev_output_vol, - bool isLinear); +CUDA_ForwardWarp(itk::SizeValueType input_vol_dim[3], + itk::SizeValueType input_dvf_dim[3], + itk::SizeValueType output_vol_dim[3], + float IndexInputToPPInputMatrix[12], + float IndexInputToIndexDVFMatrix[12], + float PPOutputToIndexOutputMatrix[12], + float * dev_input_vol, + float * dev_input_xdvf, + float * dev_input_ydvf, + float * dev_input_zdvf, + float * dev_output_vol, + bool isLinear); #endif diff --git a/include/rtkCudaInterpolateImageFilter.hcu b/include/rtkCudaInterpolateImageFilter.hcu index 5c0ba262d..98001a113 100644 --- a/include/rtkCudaInterpolateImageFilter.hcu +++ b/include/rtkCudaInterpolateImageFilter.hcu @@ -19,9 +19,13 @@ #ifndef rtkCudaInterpolateImageFilter_hcu #define rtkCudaInterpolateImageFilter_hcu -#include +#include "rtkCudaUtilities.hcu" void -CUDA_interpolation(const int4 & inputSize, float * input, float * output, int projectionNumber, float ** weights); +CUDA_interpolation(const itk::SizeValueType inputSize[4], + float * input, + float * output, + int projectionNumber, + float ** weights); #endif diff --git a/include/rtkCudaSplatImageFilter.hcu b/include/rtkCudaSplatImageFilter.hcu index dd85716ea..2ea3d5ef3 100644 --- a/include/rtkCudaSplatImageFilter.hcu +++ b/include/rtkCudaSplatImageFilter.hcu @@ -19,9 +19,13 @@ #ifndef rtkCudaSplatImageFilter_hcu #define rtkCudaSplatImageFilter_hcu -#include +#include "rtkCudaUtilities.hcu" void -CUDA_splat(const int4 & outputSize, float * input, float * output, int projectionNumber, float ** weights); +CUDA_splat(const itk::SizeValueType outputSize[4], + float * input, + float * output, + int projectionNumber, + float ** weights); #endif diff --git a/include/rtkCudaUtilities.hcu b/include/rtkCudaUtilities.hcu index 4dc5ec4d8..5a2ca23ca 100644 --- a/include/rtkCudaUtilities.hcu +++ b/include/rtkCudaUtilities.hcu @@ -23,10 +23,19 @@ #include #include #define ITK_STATIC +#include #include #undef ITK_STATIC #include +// Keep CUDA vector types as unconditional 64-bit (unsigned for sizes, signed for +// indices/offsets). Conditional mapping based on ITK_USE_64BITS_IDS like itk::SizeValueType can flip +// between 32 and 64-bit and reintroduce integer overflows on large data sizes. +using SizeValueType3 = ulonglong3; +using SizeValueType4 = ulonglong4; +using IndexValueType3 = longlong3; +using IndexValueType4 = longlong4; + #define CUDA_CHECK_ERROR \ { \ cudaError_t err = cudaGetLastError(); \ @@ -112,9 +121,9 @@ operator+=(float3 & a, float3 b) a.z += b.z; } inline __host__ __device__ int -iDivUp(int a, int b) +iDivUp(itk::SizeValueType a, itk::SizeValueType b) { - return (a % b != 0) ? (a / b + 1) : (a / b); + return static_cast((a % b != 0) ? (a / b + 1) : (a / b)); } inline __host__ __device__ float dot_vector(float3 u, float3 v) @@ -123,7 +132,7 @@ dot_vector(float3 u, float3 v) } __host__ void -prepareScalarTextureObject(int size[3], +prepareScalarTextureObject(itk::SizeValueType size[3], float * dev_ptr, cudaArray *& threeDArray, cudaTextureObject_t & tex, @@ -132,7 +141,7 @@ prepareScalarTextureObject(int size[3], const cudaTextureAddressMode texAddressMode = cudaAddressModeBorder); __host__ void -prepareVectorTextureObject(int size[3], +prepareVectorTextureObject(itk::SizeValueType size[3], const float * dev_ptr, std::vector & componentArrays, const unsigned int nComponents, diff --git a/include/rtkCudaWarpBackProjectionImageFilter.hcu b/include/rtkCudaWarpBackProjectionImageFilter.hcu index 325efe84a..6d28929d3 100644 --- a/include/rtkCudaWarpBackProjectionImageFilter.hcu +++ b/include/rtkCudaWarpBackProjectionImageFilter.hcu @@ -20,19 +20,19 @@ #define rtkCudaWarpBackProjectionImageFilter_hcu void -CUDA_warp_back_project(int projSize[3], - int volSize[3], - int dvf_size[3], - float * matrices, - float * volIndexToProjPPs, - float * projPPToProjIndex, - float * dev_vol_in, - float * dev_vol_out, - float * dev_proj, - float * dev_input_dvf, - float IndexInputToIndexDVFMatrix[12], - float PPInputToIndexInputMatrix[12], - float IndexInputToPPInputMatrix[12], - double radiusCylindricalDetector); +CUDA_warp_back_project(itk::SizeValueType projSize[3], + itk::SizeValueType volSize[3], + itk::SizeValueType dvf_size[3], + float * matrices, + float * volIndexToProjPPs, + float * projPPToProjIndex, + float * dev_vol_in, + float * dev_vol_out, + float * dev_proj, + float * dev_input_dvf, + float IndexInputToIndexDVFMatrix[12], + float PPInputToIndexInputMatrix[12], + float IndexInputToPPInputMatrix[12], + double radiusCylindricalDetector); #endif diff --git a/include/rtkCudaWarpForwardProjectionImageFilter.hcu b/include/rtkCudaWarpForwardProjectionImageFilter.hcu index 9745c70b8..a1f8e79c2 100644 --- a/include/rtkCudaWarpForwardProjectionImageFilter.hcu +++ b/include/rtkCudaWarpForwardProjectionImageFilter.hcu @@ -20,20 +20,20 @@ #define rtkCudaWarpForwardProjectionImageFilter_hcu void -CUDA_warp_forward_project(int projSize[3], - int volSize[3], - int dvfSize[3], - float * matrices, - float * dev_proj_in, - float * dev_proj_out, - float * dev_vol, - float t_step, - float * source_positions, - float box_min[3], - float box_max[3], - float spacing[3], - float * dev_input_dvf, - float IndexInputToIndexDVFMatrix[12], - float PPInputToIndexInputMatrix[12], - float IndexInputToPPInputMatrix[12]); +CUDA_warp_forward_project(itk::SizeValueType projSize[3], + itk::SizeValueType volSize[3], + itk::SizeValueType dvfSize[3], + float * matrices, + float * dev_proj_in, + float * dev_proj_out, + float * dev_vol, + float t_step, + float * source_positions, + float box_min[3], + float box_max[3], + float spacing[3], + float * dev_input_dvf, + float IndexInputToIndexDVFMatrix[12], + float PPInputToIndexInputMatrix[12], + float IndexInputToPPInputMatrix[12]); #endif diff --git a/include/rtkCudaWarpImageFilter.hcu b/include/rtkCudaWarpImageFilter.hcu index 0823e83b2..f95ed32fc 100644 --- a/include/rtkCudaWarpImageFilter.hcu +++ b/include/rtkCudaWarpImageFilter.hcu @@ -20,15 +20,15 @@ #define rtkCudaWarpImageFilter_hcu void -CUDA_warp(int input_vol_dim[3], - int input_dvf_dim[3], - int output_vol_dim[3], - float IndexOutputToPPOutputMatrix[12], - float IndexOutputToIndexDVFMatrix[12], - float PPInputToIndexInputMatrix[12], - float * dev_input_vol, - float * dev_output_vol, - float * dev_DVF, - bool isLinear); +CUDA_warp(itk::SizeValueType input_vol_dim[3], + itk::SizeValueType input_dvf_dim[3], + itk::SizeValueType output_vol_dim[3], + float IndexOutputToPPOutputMatrix[12], + float IndexOutputToIndexDVFMatrix[12], + float PPInputToIndexInputMatrix[12], + float * dev_input_vol, + float * dev_output_vol, + float * dev_DVF, + bool isLinear); #endif diff --git a/src/rtkCudaAverageOutOfROIImageFilter.cu b/src/rtkCudaAverageOutOfROIImageFilter.cu index ee2706541..e3ad22f65 100644 --- a/src/rtkCudaAverageOutOfROIImageFilter.cu +++ b/src/rtkCudaAverageOutOfROIImageFilter.cu @@ -28,7 +28,7 @@ // TEXTURES AND CONSTANTS // -__constant__ int4 c_Size; +__constant__ SizeValueType4 c_Size; //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ // K E R N E L S -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ @@ -37,18 +37,18 @@ __constant__ int4 c_Size; __global__ void -average_along_dim_4(float * in, float * out, float * roi, unsigned int strideInFloats) +average_along_dim_4(float * in, float * out, float * roi, itk::SizeValueType strideInFloats) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int k = blockIdx.z * blockDim.z + threadIdx.z; + itk::IndexValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::IndexValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::IndexValueType k = blockIdx.z * blockDim.z + threadIdx.z; if (i >= c_Size.x || j >= c_Size.y || k >= c_Size.z) return; // Compute the index of the initial voxel - long int id = (k * c_Size.y + j) * c_Size.x + i; - long int strided_id = id; // strided_id will run along the 4th dimension + itk::OffsetValueType id = (k * c_Size.y + j) * c_Size.x + i; + itk::OffsetValueType strided_id = id; // strided_id will run along the 4th dimension // Compute the average along last dimension float avg = 0; @@ -77,14 +77,13 @@ average_along_dim_4(float * in, float * out, float * roi, unsigned int strideInF //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ void -CUDA_average_out_of_ROI(int size[4], float * input, float * output, float * roi) +CUDA_average_out_of_ROI(itk::SizeValueType size[4], float * input, float * output, float * roi) { - int4 dev_Size = make_int4(size[0], size[1], size[2], size[3]); - cudaMemcpyToSymbol(c_Size, &dev_Size, sizeof(int4)); + cudaMemcpyToSymbol(c_Size, size, sizeof(SizeValueType4)); // Compute the stride (in floats, not Bytes) to jump from one voxel // to the next one along 4th dimension - unsigned int strideInFloats = size[0] * size[1] * size[2]; + itk::SizeValueType strideInFloats = size[0] * size[1] * size[2]; // Thread Block Dimensions dim3 dimBlock = dim3(8, 8, 8); diff --git a/src/rtkCudaAverageOutOfROIImageFilter.cxx b/src/rtkCudaAverageOutOfROIImageFilter.cxx index 7eca6f93c..1656d0e65 100644 --- a/src/rtkCudaAverageOutOfROIImageFilter.cxx +++ b/src/rtkCudaAverageOutOfROIImageFilter.cxx @@ -29,12 +29,11 @@ CudaAverageOutOfROIImageFilter ::CudaAverageOutOfROIImageFilter() = default; void CudaAverageOutOfROIImageFilter ::GPUGenerateData() { - int size[4]; - - for (int i = 0; i < 4; i++) - { - size[i] = this->GetOutput()->GetBufferedRegion().GetSize()[i]; - } + itk::SizeValueType size[4]; + size[0] = this->GetOutput()->GetBufferedRegion().GetSize()[0]; + size[1] = this->GetOutput()->GetBufferedRegion().GetSize()[1]; + size[2] = this->GetOutput()->GetBufferedRegion().GetSize()[2]; + size[3] = this->GetOutput()->GetBufferedRegion().GetSize()[3]; float * pin = static_cast(this->GetInput()->GetCudaDataManager()->GetGPUBufferPointer()); float * pout = static_cast(this->GetOutput()->GetCudaDataManager()->GetGPUBufferPointer()); diff --git a/src/rtkCudaBackProjectionImageFilter.cu b/src/rtkCudaBackProjectionImageFilter.cu index 1d4c24e55..91a61ba95 100644 --- a/src/rtkCudaBackProjectionImageFilter.cu +++ b/src/rtkCudaBackProjectionImageFilter.cu @@ -40,11 +40,11 @@ #include // Constant memory -__constant__ float c_matrices[SLAB_SIZE * 12]; // Can process stacks of at most SLAB_SIZE projections -__constant__ float c_volIndexToProjPP[SLAB_SIZE * 12]; -__constant__ float c_projPPToProjIndex[9]; -__constant__ int3 c_projSize; -__constant__ int3 c_volSize; +__constant__ float c_matrices[SLAB_SIZE * 12]; // Can process stacks of at most SLAB_SIZE projections +__constant__ float c_volIndexToProjPP[SLAB_SIZE * 12]; +__constant__ float c_projPPToProjIndex[9]; +__constant__ SizeValueType3 c_projSize; +__constant__ SizeValueType3 c_volSize; //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ // K E R N E L S -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ @@ -55,9 +55,9 @@ template __global__ void kernel_backProject(float * dev_vol_in, float * dev_vol_out, float radius, cudaTextureObject_t * dev_tex_proj) { - itk::SizeValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - itk::SizeValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; - itk::SizeValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; + itk::IndexValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::IndexValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; + itk::IndexValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; if (i >= c_volSize.x || j >= c_volSize.y || k >= c_volSize.z) { @@ -65,7 +65,7 @@ kernel_backProject(float * dev_vol_in, float * dev_vol_out, float radius, cudaTe } // Index row major into the volume - itk::SizeValueType vol_idx = i + (j + k * c_volSize.y) * (c_volSize.x); + itk::OffsetValueType vol_idx = i + (j + k * c_volSize.y) * (c_volSize.x); float3 ip, pp; float voxel_data[VVectorLength]; @@ -123,20 +123,20 @@ kernel_backProject(float * dev_vol_in, float * dev_vol_out, float radius, cudaTe /////////////////////////////////////////////////////////////////////////// // FUNCTION: CUDA_back_project ///////////////////////////// void -CUDA_back_project(int projSize[3], - int volSize[3], - float * matrices, - float * volIndexToProjPPs, - float * projPPToProjIndex, - float * dev_vol_in, - float * dev_vol_out, - float * dev_proj, - double radiusCylindricalDetector, - unsigned int vectorLength) +CUDA_back_project(itk::SizeValueType projSize[3], + itk::SizeValueType volSize[3], + float * matrices, + float * volIndexToProjPPs, + float * projPPToProjIndex, + float * dev_vol_in, + float * dev_vol_out, + float * dev_proj, + double radiusCylindricalDetector, + unsigned int vectorLength) { // Copy the size of inputs into constant memory - cudaMemcpyToSymbol(c_projSize, projSize, sizeof(int3)); - cudaMemcpyToSymbol(c_volSize, volSize, sizeof(int3)); + cudaMemcpyToSymbol(c_projSize, projSize, sizeof(SizeValueType3)); + cudaMemcpyToSymbol(c_volSize, volSize, sizeof(SizeValueType3)); // Copy the projection matrices into constant memory cudaMemcpyToSymbol(c_matrices, &(matrices[0]), 12 * sizeof(float) * projSize[2]); diff --git a/src/rtkCudaConjugateGradientImageFilter.cu b/src/rtkCudaConjugateGradientImageFilter.cu index 42b53f68f..7978fb4c1 100644 --- a/src/rtkCudaConjugateGradientImageFilter.cu +++ b/src/rtkCudaConjugateGradientImageFilter.cu @@ -50,7 +50,7 @@ CUDA_subtract(size_t numberOfElements, float * out, float * toBeSubtracted) cublasCreate(&handle); const float alpha = -1.0; -#if CUDA_VERSION < 12000 +#if CUBLAS_VER_MAJOR < 12 cublasSaxpy(handle, (int)numberOfElements, &alpha, toBeSubtracted, 1, out, 1); #else cublasSaxpy_64(handle, numberOfElements, &alpha, toBeSubtracted, 1, out, 1); @@ -69,7 +69,11 @@ CUDA_subtract(size_t numberOfElements, double * out, double * toBeSubtracted) cublasCreate(&handle); const double alpha = -1.0; +#if CUBLAS_VER_MAJOR < 12 cublasDaxpy(handle, (int)numberOfElements, &alpha, toBeSubtracted, 1, out, 1); +#else + cublasDaxpy_64(handle, numberOfElements, &alpha, toBeSubtracted, 1, out, 1); +#endif // Destroy Cublas context cublasDestroy(handle); @@ -87,24 +91,44 @@ CUDA_conjugate_gradient(size_t numberOfElements, float * Xk, float * Rk, float * // Compute Rk_square = sum(Rk(:).^2) by cublas float Rk_square = 0; +#if CUBLAS_VER_MAJOR < 12 cublasSdot(handle, (int)numberOfElements, Rk, 1, Rk, 1, &Rk_square); +#else + cublasSdot_64(handle, numberOfElements, Rk, 1, Rk, 1, &Rk_square); +#endif // Compute alpha_k = Rk_square / sum(Pk(:) .* APk(:)) float Pk_APk = 0; +#if CUBLAS_VER_MAJOR < 12 cublasSdot(handle, (int)numberOfElements, Pk, 1, APk, 1, &Pk_APk); +#else + cublasSdot_64(handle, numberOfElements, Pk, 1, APk, 1, &Pk_APk); +#endif const float alpha_k = Rk_square / (Pk_APk + eps); const float minus_alpha_k = -alpha_k; // Compute Xk+1 = Xk + alpha_k * Pk +#if CUBLAS_VER_MAJOR < 12 cublasSaxpy(handle, (int)numberOfElements, &alpha_k, Pk, 1, Xk, 1); +#else + cublasSaxpy_64(handle, numberOfElements, &alpha_k, Pk, 1, Xk, 1); +#endif // Compute Rk+1 = Rk - alpha_k * APk +#if CUBLAS_VER_MAJOR < 12 cublasSaxpy(handle, (int)numberOfElements, &minus_alpha_k, APk, 1, Rk, 1); +#else + cublasSaxpy_64(handle, numberOfElements, &minus_alpha_k, APk, 1, Rk, 1); +#endif // Compute beta_k = sum(Rk+1(:).^2) / Rk_square float Rkplusone_square = 0; +#if CUBLAS_VER_MAJOR < 12 cublasSdot(handle, (int)numberOfElements, Rk, 1, Rk, 1, &Rkplusone_square); +#else + cublasSdot_64(handle, numberOfElements, Rk, 1, Rk, 1, &Rkplusone_square); +#endif float beta_k = Rkplusone_square / (Rk_square + eps); float one = 1.0; @@ -112,8 +136,13 @@ CUDA_conjugate_gradient(size_t numberOfElements, float * Xk, float * Rk, float * // Compute Pk+1 = Rk+1 + beta_k * Pk // This requires two cublas functions, // since axpy would store the result in the wrong array +#if CUBLAS_VER_MAJOR < 12 cublasSscal(handle, (int)numberOfElements, &beta_k, Pk, 1); cublasSaxpy(handle, (int)numberOfElements, &one, Rk, 1, Pk, 1); +#else + cublasSscal_64(handle, numberOfElements, &beta_k, Pk, 1); + cublasSaxpy_64(handle, numberOfElements, &one, Rk, 1, Pk, 1); +#endif // Destroy Cublas context cublasDestroy(handle); @@ -131,24 +160,44 @@ CUDA_conjugate_gradient(size_t numberOfElements, double * Xk, double * Rk, doubl // Compute Rk_square = sum(Rk(:).^2) by cublas double Rk_square = 0; +#if CUBLAS_VER_MAJOR < 12 cublasDdot(handle, (int)numberOfElements, Rk, 1, Rk, 1, &Rk_square); +#else + cublasDdot_64(handle, numberOfElements, Rk, 1, Rk, 1, &Rk_square); +#endif // Compute alpha_k = Rk_square / sum(Pk(:) .* APk(:)) double Pk_APk = 0; +#if CUBLAS_VER_MAJOR < 12 cublasDdot(handle, (int)numberOfElements, Pk, 1, APk, 1, &Pk_APk); +#else + cublasDdot_64(handle, numberOfElements, Pk, 1, APk, 1, &Pk_APk); +#endif const double alpha_k = Rk_square / (Pk_APk + eps); const double minus_alpha_k = -alpha_k; // Compute Xk+1 = Xk + alpha_k * Pk +#if CUBLAS_VER_MAJOR < 12 cublasDaxpy(handle, (int)numberOfElements, &alpha_k, Pk, 1, Xk, 1); +#else + cublasDaxpy_64(handle, numberOfElements, &alpha_k, Pk, 1, Xk, 1); +#endif // Compute Rk+1 = Rk - alpha_k * APk +#if CUBLAS_VER_MAJOR < 12 cublasDaxpy(handle, (int)numberOfElements, &minus_alpha_k, APk, 1, Rk, 1); +#else + cublasDaxpy_64(handle, numberOfElements, &minus_alpha_k, APk, 1, Rk, 1); +#endif // Compute beta_k = sum(Rk+1(:).^2) / Rk_square double Rkplusone_square = 0; +#if CUBLAS_VER_MAJOR < 12 cublasDdot(handle, (int)numberOfElements, Rk, 1, Rk, 1, &Rkplusone_square); +#else + cublasDdot_64(handle, numberOfElements, Rk, 1, Rk, 1, &Rkplusone_square); +#endif double beta_k = Rkplusone_square / (Rk_square + eps); double one = 1.0; @@ -156,8 +205,13 @@ CUDA_conjugate_gradient(size_t numberOfElements, double * Xk, double * Rk, doubl // Compute Pk+1 = Rk+1 + beta_k * Pk // This requires two cublas functions, // since axpy would store the result in the wrong array +#if CUBLAS_VER_MAJOR < 12 cublasDscal(handle, (int)numberOfElements, &beta_k, Pk, 1); cublasDaxpy(handle, (int)numberOfElements, &one, Rk, 1, Pk, 1); +#else + cublasDscal_64(handle, numberOfElements, &beta_k, Pk, 1); + cublasDaxpy_64(handle, numberOfElements, &one, Rk, 1, Pk, 1); +#endif // Destroy Cublas context cublasDestroy(handle); diff --git a/src/rtkCudaConstantVolumeSeriesSource.cu b/src/rtkCudaConstantVolumeSeriesSource.cu index a420cb973..190ff9c98 100644 --- a/src/rtkCudaConstantVolumeSeriesSource.cu +++ b/src/rtkCudaConstantVolumeSeriesSource.cu @@ -27,7 +27,7 @@ // TEXTURES AND CONSTANTS // -__constant__ int4 c_Size; +__constant__ SizeValueType4 c_Size; //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ // K E R N E L S -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ @@ -38,14 +38,14 @@ __constant__ int4 c_Size; __global__ void set_volume_series_to_constant(float * out, float value) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int k = blockIdx.z * blockDim.z + threadIdx.z; + itk::IndexValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::IndexValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::IndexValueType k = blockIdx.z * blockDim.z + threadIdx.z; - if (i >= c_Size.x || j >= c_Size.y || k >= c_Size.z * c_Size.w) + if (i >= c_Size.x || j >= c_Size.y || k >= c_Size.z) return; - long int id = (k * c_Size.y + j) * c_Size.x + i; + itk::OffsetValueType id = (k * c_Size.y + j) * c_Size.x + i; out[id] = value; } @@ -57,10 +57,9 @@ set_volume_series_to_constant(float * out, float value) //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ void -CUDA_generate_constant_volume_series(int size[4], float * dev_out, float constantValue) +CUDA_generate_constant_volume_series(itk::SizeValueType size[4], float * dev_out, float constantValue) { - int4 dev_Size = make_int4(size[0], size[1], size[2], size[3]); - cudaMemcpyToSymbol(c_Size, &dev_Size, sizeof(int4)); + cudaMemcpyToSymbol(c_Size, size, sizeof(SizeValueType4)); // NOTE : memset sets every BYTE of the memory to the value given in // argument (here, 0). With 0, it's fine, as all number formats seem diff --git a/src/rtkCudaConstantVolumeSeriesSource.cxx b/src/rtkCudaConstantVolumeSeriesSource.cxx index 32f51196c..c74e5e37c 100644 --- a/src/rtkCudaConstantVolumeSeriesSource.cxx +++ b/src/rtkCudaConstantVolumeSeriesSource.cxx @@ -29,12 +29,11 @@ CudaConstantVolumeSeriesSource ::CudaConstantVolumeSeriesSource() = default; void CudaConstantVolumeSeriesSource ::GPUGenerateData() { - int outputSize[4]; - - for (int i = 0; i < 4; i++) - { - outputSize[i] = this->GetOutput()->GetRequestedRegion().GetSize()[i]; - } + itk::SizeValueType outputSize[4]; + outputSize[0] = this->GetOutput()->GetRequestedRegion().GetSize()[0]; + outputSize[1] = this->GetOutput()->GetRequestedRegion().GetSize()[1]; + outputSize[2] = this->GetOutput()->GetRequestedRegion().GetSize()[2]; + outputSize[3] = this->GetOutput()->GetRequestedRegion().GetSize()[3]; float * pout = static_cast(this->GetOutput()->GetCudaDataManager()->GetGPUBufferPointer()); diff --git a/src/rtkCudaConstantVolumeSource.cu b/src/rtkCudaConstantVolumeSource.cu index c0a6ddbae..2f7659201 100644 --- a/src/rtkCudaConstantVolumeSource.cu +++ b/src/rtkCudaConstantVolumeSource.cu @@ -27,7 +27,7 @@ // TEXTURES AND CONSTANTS // -__constant__ int3 c_Size; +__constant__ SizeValueType3 c_Size; //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ // K E R N E L S -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ @@ -38,14 +38,14 @@ __constant__ int3 c_Size; __global__ void set_volume_to_constant(float * out, float value) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int k = blockIdx.z * blockDim.z + threadIdx.z; + itk::IndexValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::IndexValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::IndexValueType k = blockIdx.z * blockDim.z + threadIdx.z; if (i >= c_Size.x || j >= c_Size.y || k >= c_Size.z) return; - long int id = (k * c_Size.y + j) * c_Size.x + i; + itk::OffsetValueType id = (k * c_Size.y + j) * c_Size.x + i; out[id] = value; } @@ -57,10 +57,9 @@ set_volume_to_constant(float * out, float value) //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ void -CUDA_generate_constant_volume(int size[3], float * dev_out, float constantValue) +CUDA_generate_constant_volume(itk::SizeValueType size[3], float * dev_out, float constantValue) { - int3 dev_Size = make_int3(size[0], size[1], size[2]); - cudaMemcpyToSymbol(c_Size, &dev_Size, sizeof(int3)); + cudaMemcpyToSymbol(c_Size, size, sizeof(SizeValueType3)); // NOTE : memset sets every BYTE of the memory to the value given in // argument (here, 0). With 0, it's fine, as all number formats seem diff --git a/src/rtkCudaConstantVolumeSource.cxx b/src/rtkCudaConstantVolumeSource.cxx index 4717d0d07..501c5efb9 100644 --- a/src/rtkCudaConstantVolumeSource.cxx +++ b/src/rtkCudaConstantVolumeSource.cxx @@ -29,12 +29,10 @@ CudaConstantVolumeSource ::CudaConstantVolumeSource() = default; void CudaConstantVolumeSource ::GPUGenerateData() { - int outputSize[3]; - - for (int i = 0; i < 3; i++) - { - outputSize[i] = this->GetOutput()->GetRequestedRegion().GetSize()[i]; - } + itk::SizeValueType outputSize[3]; + outputSize[0] = this->GetOutput()->GetRequestedRegion().GetSize()[0]; + outputSize[1] = this->GetOutput()->GetRequestedRegion().GetSize()[1]; + outputSize[2] = this->GetOutput()->GetRequestedRegion().GetSize()[2]; float * pout = static_cast(this->GetOutput()->GetCudaDataManager()->GetGPUBufferPointer()); diff --git a/src/rtkCudaCropImageFilter.cu b/src/rtkCudaCropImageFilter.cu index 75995c469..a899dcd38 100644 --- a/src/rtkCudaCropImageFilter.cu +++ b/src/rtkCudaCropImageFilter.cu @@ -33,17 +33,18 @@ crop_kernel(float * input, const uint3 inputDim, const unsigned int Blocks_Y) { - unsigned int blockIdx_z = blockIdx.y / Blocks_Y; - unsigned int blockIdx_y = blockIdx.y - __umul24(blockIdx_z, Blocks_Y); - unsigned int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - unsigned int j = __umul24(blockIdx_y, blockDim.y) + threadIdx.y; - unsigned int k = __umul24(blockIdx_z, blockDim.z) + threadIdx.z; + unsigned int blockIdx_z = blockIdx.y / Blocks_Y; + unsigned int blockIdx_y = blockIdx.y - __umul24(blockIdx_z, Blocks_Y); + itk::IndexValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::IndexValueType j = __umul24(blockIdx_y, blockDim.y) + threadIdx.y; + itk::IndexValueType k = __umul24(blockIdx_z, blockDim.z) + threadIdx.z; if (i >= cropDim.x || j >= cropDim.y || k >= cropDim.z) return; - unsigned long int out_idx = i + j * cropDim.x + k * cropDim.y * cropDim.x; - unsigned long int in_idx = (cropIdx.x + i) + (cropIdx.y + j) * inputDim.x + (cropIdx.z + k) * inputDim.y * inputDim.x; + itk::OffsetValueType out_idx = i + j * cropDim.x + k * cropDim.y * cropDim.x; + itk::OffsetValueType in_idx = + cropIdx.x + i + (cropIdx.y + j) * inputDim.x + (cropIdx.z + k) * inputDim.y * inputDim.x; output[out_idx] = input[in_idx]; } diff --git a/src/rtkCudaCyclicDeformationImageFilter.cu b/src/rtkCudaCyclicDeformationImageFilter.cu index 53bcce3e7..0d458a189 100644 --- a/src/rtkCudaCyclicDeformationImageFilter.cu +++ b/src/rtkCudaCyclicDeformationImageFilter.cu @@ -29,21 +29,19 @@ // TEXTURES AND CONSTANTS // -__constant__ int4 c_inputSize; - //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ // K E R N E L S -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_( S T A R T )_ //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ void -CUDA_linear_interpolate_along_fourth_dimension(unsigned int inputSize[4], - float * input, - float * output, - unsigned int frameInf, - unsigned int frameSup, - double weightInf, - double weightSup) +CUDA_linear_interpolate_along_fourth_dimension(itk::SizeValueType inputSize[4], + float * input, + float * output, + unsigned int frameInf, + unsigned int frameSup, + double weightInf, + double weightSup) { cublasHandle_t handle; cublasCreate(&handle); diff --git a/src/rtkCudaCyclicDeformationImageFilter.cxx b/src/rtkCudaCyclicDeformationImageFilter.cxx index 5234d60e4..827d998dd 100644 --- a/src/rtkCudaCyclicDeformationImageFilter.cxx +++ b/src/rtkCudaCyclicDeformationImageFilter.cxx @@ -34,17 +34,17 @@ CudaCyclicDeformationImageFilter ::GPUGenerateData() this->Superclass::BeforeThreadedGenerateData(); // Prepare the data to perform the linear interpolation on GPU - unsigned int inputSize[4]; - - for (unsigned int i = 0; i < 4; i++) - inputSize[i] = this->GetInput()->GetBufferedRegion().GetSize()[i]; - for (unsigned int i = 0; i < 3; i++) + itk::SizeValueType inputSize[4]; + inputSize[0] = this->GetInput()->GetBufferedRegion().GetSize()[0]; + inputSize[1] = this->GetInput()->GetBufferedRegion().GetSize()[1]; + inputSize[2] = this->GetInput()->GetBufferedRegion().GetSize()[2]; + inputSize[3] = this->GetInput()->GetBufferedRegion().GetSize()[3]; + if ((this->GetOutput()->GetRequestedRegion().GetSize()[0] != inputSize[0]) || + (this->GetOutput()->GetRequestedRegion().GetSize()[1] != inputSize[1]) || + (this->GetOutput()->GetRequestedRegion().GetSize()[2] != inputSize[2])) { - if (this->GetOutput()->GetRequestedRegion().GetSize()[i] != inputSize[i]) - { - itkExceptionMacro("In rtk::CudaCyclicDeformationImageFilter: the output's requested region must have the same " - "size as the input's buffered region on the first 3 dimensions"); - } + itkExceptionMacro("In rtk::CudaCyclicDeformationImageFilter: the output's requested region must have the same " + "size as the input's buffered region on the first 3 dimensions"); } float * pin = static_cast(this->GetInput()->GetCudaDataManager()->GetGPUBufferPointer()); diff --git a/src/rtkCudaDisplacedDetectorImageFilter.cu b/src/rtkCudaDisplacedDetectorImageFilter.cu index b888487b3..f2fa95eb0 100644 --- a/src/rtkCudaDisplacedDetectorImageFilter.cu +++ b/src/rtkCudaDisplacedDetectorImageFilter.cu @@ -60,11 +60,12 @@ kernel_displaced_weight(int3 proj_idx_in, ) { // compute thread index - int3 tIdx; + IndexValueType3 tIdx; tIdx.x = blockIdx.x * blockDim.x + threadIdx.x; tIdx.y = blockIdx.y * blockDim.y + threadIdx.y; tIdx.z = blockIdx.z * blockDim.z + threadIdx.z; - long int tIdx_comp = tIdx.x + tIdx.y * proj_size_out.x + tIdx.z * proj_size_out_buf.x * proj_size_out_buf.y; + itk::OffsetValueType tIdx_comp = + tIdx.x + tIdx.y * proj_size_out.x + tIdx.z * proj_size_out_buf.x * proj_size_out_buf.y; // check if outside of projection grid if (tIdx.x >= proj_size_out.x || tIdx.y >= proj_size_out.y || tIdx.z >= proj_size_out.z) @@ -72,9 +73,10 @@ kernel_displaced_weight(int3 proj_idx_in, // compute projection index from thread index int3 pIdx = make_int3(tIdx.x + proj_idx_out.x, tIdx.y + proj_idx_out.y, tIdx.z + proj_idx_out.z); - // combined proj. index -> use thread index in z because accessing memory only with this index - long int pIdx_comp = (pIdx.x - proj_idx_in.x) + (pIdx.y - proj_idx_in.y) * proj_size_in_buf.x + - (pIdx.z - proj_idx_in.z) * proj_size_in_buf.x * proj_size_in_buf.y; + // combined proj. index + itk::OffsetValueType pIdx_comp = (tIdx.x - proj_idx_in.x + proj_idx_out.x) + + (tIdx.y - proj_idx_in.y + proj_idx_out.y) * proj_size_in_buf.x + + (tIdx.z - proj_idx_in.z + proj_idx_out.z) * proj_size_in_buf.x * proj_size_in_buf.y; // check if outside overlapping region if (pIdx.x < proj_idx_in.x || pIdx.x >= (proj_idx_in.x + proj_size_in.x) || pIdx.y < proj_idx_in.y || diff --git a/src/rtkCudaFDKBackProjectionImageFilter.cu b/src/rtkCudaFDKBackProjectionImageFilter.cu index 50d8532d3..1b51b223c 100644 --- a/src/rtkCudaFDKBackProjectionImageFilter.cu +++ b/src/rtkCudaFDKBackProjectionImageFilter.cu @@ -40,9 +40,9 @@ #include // Constant memory -__constant__ float c_matrices[SLAB_SIZE * 12]; // Can process stacks of at most SLAB_SIZE projections -__constant__ int3 c_projSize; -__constant__ int3 c_vol_size; +__constant__ float c_matrices[SLAB_SIZE * 12]; // Can process stacks of at most SLAB_SIZE projections +__constant__ SizeValueType3 c_projSize; +__constant__ SizeValueType3 c_vol_size; //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ // K E R N E L S -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ @@ -52,9 +52,9 @@ __constant__ int3 c_vol_size; __global__ void kernel_fdk_3Dgrid(float * dev_vol_in, float * dev_vol_out, cudaTextureObject_t tex_proj) { - itk::SizeValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - itk::SizeValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; - itk::SizeValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; + itk::IndexValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::IndexValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; + itk::IndexValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; if (i >= c_vol_size.x || j >= c_vol_size.y || k >= c_vol_size.z) { @@ -62,7 +62,7 @@ kernel_fdk_3Dgrid(float * dev_vol_in, float * dev_vol_out, cudaTextureObject_t t } // Index row major into the volume - itk::SizeValueType vol_idx = i + (j + k * c_vol_size.y) * (c_vol_size.x); + itk::OffsetValueType vol_idx = i + (j + k * c_vol_size.y) * (c_vol_size.x); float3 ip; float voxel_data = 0; @@ -93,16 +93,16 @@ kernel_fdk_3Dgrid(float * dev_vol_in, float * dev_vol_out, cudaTextureObject_t t /////////////////////////////////////////////////////////////////////////// // FUNCTION: CUDA_back_project ///////////////////////////// void -CUDA_reconstruct_conebeam(int proj_size[3], - int vol_size[3], - float * matrices, - float * dev_vol_in, - float * dev_vol_out, - float * dev_proj) +CUDA_reconstruct_conebeam(itk::SizeValueType proj_size[3], + itk::SizeValueType vol_size[3], + float * matrices, + float * dev_vol_in, + float * dev_vol_out, + float * dev_proj) { // Copy the size of inputs into constant memory - cudaMemcpyToSymbol(c_projSize, proj_size, sizeof(int3)); - cudaMemcpyToSymbol(c_vol_size, vol_size, sizeof(int3)); + cudaMemcpyToSymbol(c_projSize, proj_size, sizeof(SizeValueType3)); + cudaMemcpyToSymbol(c_vol_size, vol_size, sizeof(SizeValueType3)); // Copy the projection matrices into constant memory cudaMemcpyToSymbol(c_matrices, &(matrices[0]), 12 * sizeof(float) * proj_size[2]); diff --git a/src/rtkCudaFDKBackProjectionImageFilter.cxx b/src/rtkCudaFDKBackProjectionImageFilter.cxx index 7a65ac0c9..062fc95bc 100644 --- a/src/rtkCudaFDKBackProjectionImageFilter.cxx +++ b/src/rtkCudaFDKBackProjectionImageFilter.cxx @@ -59,11 +59,11 @@ CudaFDKBackProjectionImageFilter ::GPUGenerateData() } // Cuda convenient format for dimensions - int projectionSize[3]; + itk::SizeValueType projectionSize[3]; projectionSize[0] = this->GetInput(1)->GetBufferedRegion().GetSize()[0]; projectionSize[1] = this->GetInput(1)->GetBufferedRegion().GetSize()[1]; - int volumeSize[3]; + itk::SizeValueType volumeSize[3]; volumeSize[0] = this->GetOutput()->GetBufferedRegion().GetSize()[0]; volumeSize[1] = this->GetOutput()->GetBufferedRegion().GetSize()[1]; volumeSize[2] = this->GetOutput()->GetBufferedRegion().GetSize()[2]; diff --git a/src/rtkCudaFDKWeightProjectionFilter.cu b/src/rtkCudaFDKWeightProjectionFilter.cu index 5f274980f..0f9404cd3 100644 --- a/src/rtkCudaFDKWeightProjectionFilter.cu +++ b/src/rtkCudaFDKWeightProjectionFilter.cu @@ -40,38 +40,37 @@ kernel_weight_projection(int2 proj_idx, cudaTextureObject_t tex_geom // geometry texture object ) { - // compute projection index (== thread index) - int3 pIdx; - pIdx.x = blockIdx.x * blockDim.x + threadIdx.x; - pIdx.y = blockIdx.y * blockDim.y + threadIdx.y; - pIdx.z = blockIdx.z * blockDim.z + threadIdx.z; - long int pIdx_comp_in = pIdx.x + (pIdx.y + pIdx.z * proj_size_buf_in.y) * (proj_size_buf_in.x); - long int pIdx_comp_out = pIdx.x + (pIdx.y + pIdx.z * proj_size_buf_out.y) * (proj_size_buf_out.x); + // compute projection index (== thread index), 64-bit to avoid 32-bit overflow of combined indices + itk::IndexValueType pIdx_x = blockIdx.x * blockDim.x + threadIdx.x; + itk::IndexValueType pIdx_y = blockIdx.y * blockDim.y + threadIdx.y; + itk::IndexValueType pIdx_z = blockIdx.z * blockDim.z + threadIdx.z; + itk::OffsetValueType pIdx_comp_in = pIdx_x + (pIdx_y + pIdx_z * proj_size_buf_in.y) * proj_size_buf_in.x; + itk::OffsetValueType pIdx_comp_out = pIdx_x + (pIdx_y + pIdx_z * proj_size_buf_out.y) * proj_size_buf_out.x; // check if outside of projection grid - if (pIdx.x >= proj_size.x || pIdx.y >= proj_size.y || pIdx.z >= proj_size.z) + if (pIdx_x >= proj_size.x || pIdx_y >= proj_size.y || pIdx_z >= proj_size.z) return; - const float sdd = tex1Dfetch(tex_geom, pIdx.z * 7); - const float sid = tex1Dfetch(tex_geom, pIdx.z * 7 + 1); - const float wFac = tex1Dfetch(tex_geom, pIdx.z * 7 + 5); + const float sdd = tex1Dfetch(tex_geom, pIdx_z * 7); + const float sid = tex1Dfetch(tex_geom, pIdx_z * 7 + 1); + const float wFac = tex1Dfetch(tex_geom, pIdx_z * 7 + 5); if (sdd == 0) // parallel { dev_proj_out[pIdx_comp_out] = dev_proj_in[pIdx_comp_in] * wFac; } else // divergent { - const float pOffX = tex1Dfetch(tex_geom, pIdx.z * 7 + 2); - const float pOffY = tex1Dfetch(tex_geom, pIdx.z * 7 + 3); - const float sOffY = tex1Dfetch(tex_geom, pIdx.z * 7 + 4); - const float tAngle = tex1Dfetch(tex_geom, pIdx.z * 7 + 6); + const float pOffX = tex1Dfetch(tex_geom, pIdx_z * 7 + 2); + const float pOffY = tex1Dfetch(tex_geom, pIdx_z * 7 + 3); + const float sOffY = tex1Dfetch(tex_geom, pIdx_z * 7 + 4); + const float tAngle = tex1Dfetch(tex_geom, pIdx_z * 7 + 6); const float sina = sin(tAngle); const float cosa = cos(tAngle); const float tana = tan(tAngle); // compute projection point from index float2 pPoint = - TransformIndexToPhysicalPoint(make_int2(pIdx.x + proj_idx.x, pIdx.y + proj_idx.y), proj_orig, proj_row, proj_col); + TransformIndexToPhysicalPoint(make_int2(pIdx_x + proj_idx.x, pIdx_y + proj_idx.y), proj_orig, proj_row, proj_col); pPoint.x = pPoint.x + pOffX + tana * (sdd - sid); pPoint.y = pPoint.y + pOffY - sOffY; diff --git a/src/rtkCudaFFTProjectionsConvolutionImageFilter.cu b/src/rtkCudaFFTProjectionsConvolutionImageFilter.cu index da33d935c..425605ae2 100644 --- a/src/rtkCudaFFTProjectionsConvolutionImageFilter.cu +++ b/src/rtkCudaFFTProjectionsConvolutionImageFilter.cu @@ -29,16 +29,16 @@ __global__ void multiply_kernel(cufftComplex * projFFT, int3 fftDimension, cufftComplex * kernelFFT, unsigned int Blocks_Y) { - unsigned int blockIdx_z = blockIdx.y / Blocks_Y; - unsigned int blockIdx_y = blockIdx.y - __umul24(blockIdx_z, Blocks_Y); - unsigned int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - unsigned int j = __umul24(blockIdx_y, blockDim.y) + threadIdx.y; - unsigned int k = __umul24(blockIdx_z, blockDim.z) + threadIdx.z; + unsigned int blockIdx_z = blockIdx.y / Blocks_Y; + unsigned int blockIdx_y = blockIdx.y - __umul24(blockIdx_z, Blocks_Y); + itk::IndexValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::IndexValueType j = __umul24(blockIdx_y, blockDim.y) + threadIdx.y; + itk::IndexValueType k = __umul24(blockIdx_z, blockDim.z) + threadIdx.z; if (i >= fftDimension.x || j >= fftDimension.y || k >= fftDimension.z) return; - long int proj_idx = i + (j + k * fftDimension.y) * fftDimension.x; + itk::OffsetValueType proj_idx = i + (j + k * fftDimension.y) * fftDimension.x; cufftComplex result; result.x = projFFT[proj_idx].x * kernelFFT[i].x - projFFT[proj_idx].y * kernelFFT[i].y; @@ -49,17 +49,17 @@ multiply_kernel(cufftComplex * projFFT, int3 fftDimension, cufftComplex * kernel __global__ void multiply_kernel2D(cufftComplex * projFFT, int3 fftDimension, cufftComplex * kernelFFT, unsigned int Blocks_Y) { - unsigned int blockIdx_z = blockIdx.y / Blocks_Y; - unsigned int blockIdx_y = blockIdx.y - __umul24(blockIdx_z, Blocks_Y); - unsigned int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - unsigned int j = __umul24(blockIdx_y, blockDim.y) + threadIdx.y; - unsigned int k = __umul24(blockIdx_z, blockDim.z) + threadIdx.z; + unsigned int blockIdx_z = blockIdx.y / Blocks_Y; + unsigned int blockIdx_y = blockIdx.y - __umul24(blockIdx_z, Blocks_Y); + itk::IndexValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::IndexValueType j = __umul24(blockIdx_y, blockDim.y) + threadIdx.y; + itk::IndexValueType k = __umul24(blockIdx_z, blockDim.z) + threadIdx.z; if (i >= fftDimension.x || j >= fftDimension.y || k >= fftDimension.z) return; - long int kernel_idx = i + j * fftDimension.x; - long int proj_idx = kernel_idx + k * fftDimension.y * fftDimension.x; + itk::OffsetValueType kernel_idx = i + j * fftDimension.x; + itk::OffsetValueType proj_idx = kernel_idx + k * fftDimension.y * fftDimension.x; cufftComplex result; result.x = projFFT[proj_idx].x * kernelFFT[kernel_idx].x - projFFT[proj_idx].y * kernelFFT[kernel_idx].y; @@ -134,16 +134,16 @@ padding_kernel(float * input, float * truncationWeights, size_t sizeWeights) { - unsigned int blockIdx_z = blockIdx.y / Blocks_Y; - unsigned int blockIdx_y = blockIdx.y - __umul24(blockIdx_z, Blocks_Y); - int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - int j = __umul24(blockIdx_y, blockDim.y) + threadIdx.y; - int k = __umul24(blockIdx_z, blockDim.z) + threadIdx.z; + unsigned int blockIdx_z = blockIdx.y / Blocks_Y; + unsigned int blockIdx_y = blockIdx.y - __umul24(blockIdx_z, Blocks_Y); + itk::IndexValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::IndexValueType j = __umul24(blockIdx_y, blockDim.y) + threadIdx.y; + itk::IndexValueType k = __umul24(blockIdx_z, blockDim.z) + threadIdx.z; if (i >= paddingDim.x || j >= paddingDim.y || k >= paddingDim.z) return; - unsigned long int out_idx = i + (j + k * paddingDim.y) * paddingDim.x; + itk::OffsetValueType out_idx = i + (j + k * paddingDim.y) * paddingDim.x; i -= paddingIdx.x; j -= paddingIdx.y; k -= paddingIdx.z; @@ -157,14 +157,14 @@ padding_kernel(float * input, // left mirroring (equation 3a in [Ohnesorge et al, Med Phys, 2000]) else if (i < 0 && -i < sizeWeights) { - int begRow = (j + k * inputDim.y) * inputDim.x; + itk::OffsetValueType begRow = (j + k * inputDim.y) * inputDim.x; output[out_idx] = (2 * input[begRow + 1] - input[-i + begRow]) * truncationWeights[-i]; } // right mirroring (equation 3b in [Ohnesorge et al, Med Phys, 2000]) else if ((i >= inputDim.x) && (i - inputDim.x + 1) < sizeWeights) { - unsigned int borderDist = i - inputDim.x + 1; - int endRow = inputDim.x - 1 + (j + k * inputDim.y) * inputDim.x; + itk::OffsetValueType borderDist = i - inputDim.x + 1; + itk::OffsetValueType endRow = inputDim.x - 1 + (j + k * inputDim.y) * inputDim.x; output[out_idx] = (2 * input[endRow] - input[endRow - borderDist]) * truncationWeights[borderDist]; } // zero padding diff --git a/src/rtkCudaFirstOrderKernels.cu b/src/rtkCudaFirstOrderKernels.cu index 3d3a788a1..add1c8cc0 100644 --- a/src/rtkCudaFirstOrderKernels.cu +++ b/src/rtkCudaFirstOrderKernels.cu @@ -18,20 +18,22 @@ #include "rtkCudaFirstOrderKernels.hcu" +#include "rtkCudaUtilities.hcu" + __global__ void divergence_kernel(float * grad_x, float * grad_y, float * grad_z, float * out, int3 c_Size, float3 c_Spacing) { - int i = blockIdx.x * blockDim.x + threadIdx.x; - int j = blockIdx.y * blockDim.y + threadIdx.y; - int k = blockIdx.z * blockDim.z + threadIdx.z; + itk::IndexValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::IndexValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::IndexValueType k = blockIdx.z * blockDim.z + threadIdx.z; if (i >= c_Size.x || j >= c_Size.y || k >= c_Size.z) return; - long int id = (k * c_Size.y + j) * c_Size.x + i; - long int id_x = (k * c_Size.y + j) * c_Size.x + i - 1; - long int id_y = (k * c_Size.y + j - 1) * c_Size.x + i; - long int id_z = ((k - 1) * c_Size.y + j) * c_Size.x + i; + itk::OffsetValueType id = (k * c_Size.y + j) * c_Size.x + i; + itk::OffsetValueType id_x = (k * c_Size.y + j) * c_Size.x + i - 1; + itk::OffsetValueType id_y = (k * c_Size.y + j - 1) * c_Size.x + i; + itk::OffsetValueType id_z = ((k - 1) * c_Size.y + j) * c_Size.x + i; float3 A; float3 B; @@ -69,17 +71,17 @@ divergence_kernel(float * grad_x, float * grad_y, float * grad_z, float * out, i __global__ void gradient_kernel(float * in, float * grad_x, float * grad_y, float * grad_z, int3 c_Size, float3 c_Spacing) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int k = blockIdx.z * blockDim.z + threadIdx.z; + itk::IndexValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::IndexValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::IndexValueType k = blockIdx.z * blockDim.z + threadIdx.z; if (i >= c_Size.x || j >= c_Size.y || k >= c_Size.z) return; - long int id = (k * c_Size.y + j) * c_Size.x + i; - long int id_x = (k * c_Size.y + j) * c_Size.x + i + 1; - long int id_y = (k * c_Size.y + j + 1) * c_Size.x + i; - long int id_z = ((k + 1) * c_Size.y + j) * c_Size.x + i; + itk::OffsetValueType id = (k * c_Size.y + j) * c_Size.x + i; + itk::OffsetValueType id_x = (k * c_Size.y + j) * c_Size.x + i + 1; + itk::OffsetValueType id_y = (k * c_Size.y + j + 1) * c_Size.x + i; + itk::OffsetValueType id_z = ((k + 1) * c_Size.y + j) * c_Size.x + i; if (i == (c_Size.x - 1)) grad_x[id] = 0; diff --git a/src/rtkCudaForwardProjectionImageFilter.cu b/src/rtkCudaForwardProjectionImageFilter.cu index 6494bc461..59cbf2ff7 100644 --- a/src/rtkCudaForwardProjectionImageFilter.cu +++ b/src/rtkCudaForwardProjectionImageFilter.cu @@ -39,13 +39,13 @@ #include // CONSTANTS // -__constant__ int3 c_projSize; -__constant__ float3 c_boxMin; -__constant__ float3 c_boxMax; -__constant__ float3 c_spacing; -__constant__ int3 c_volSize; -__constant__ float c_tStep; -__constant__ float c_radius; +__constant__ SizeValueType3 c_projSize; +__constant__ float3 c_boxMin; +__constant__ float3 c_boxMax; +__constant__ float3 c_spacing; +__constant__ SizeValueType3 c_volSize; +__constant__ float c_tStep; +__constant__ float c_radius; __constant__ float c_translatedProjectionIndexTransformMatrices[SLAB_SIZE * 12]; // Can process stacks of at most SLAB_SIZE projections __constant__ float @@ -63,9 +63,9 @@ template __global__ void kernel_forwardProject(float * dev_proj_in, float * dev_proj_out, cudaTextureObject_t * dev_tex_vol) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int numThread = j * c_projSize.x + i; + itk::SizeValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::SizeValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::SizeValueType numThread = j * c_projSize.x + i; if (i >= c_projSize.x || j >= c_projSize.y) return; @@ -75,7 +75,7 @@ kernel_forwardProject(float * dev_proj_in, float * dev_proj_out, cudaTextureObje float3 pixelPos; float tnear, tfar; - for (unsigned int proj = 0; proj < c_projSize.z; proj++) + for (itk::SizeValueType 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]); @@ -96,7 +96,7 @@ kernel_forwardProject(float * dev_proj_in, float * dev_proj_out, cudaTextureObje ray.d = pixelPos - ray.o; - int projOffset = numThread + proj * c_projSize.x * c_projSize.y; + itk::SizeValueType projOffset = numThread + proj * c_projSize.x * c_projSize.y; // Detect intersection with box if (!intersectBox(ray, &tnear, &tfar, c_boxMin, c_boxMax) || tnear >= 1.0f || tfar <= 0.f || tfar == tnear) @@ -160,27 +160,27 @@ kernel_forwardProject(float * dev_proj_in, float * dev_proj_out, cudaTextureObje /////////////////////////////////////////////////////////////////////////// // FUNCTION: CUDA_forward_project() ////////////////////////////////// void -CUDA_forward_project(int projSize[3], - int volSize[3], - float * translatedProjectionIndexTransformMatrices, - float * translatedVolumeTransformMatrices, - float * dev_proj_in, - float * dev_proj_out, - float * dev_vol, - float t_step, - float * source_positions, - float radiusCylindricalDetector, - float box_min[3], - float box_max[3], - float spacing[3], - unsigned int vectorLength) +CUDA_forward_project(itk::SizeValueType projSize[3], + itk::SizeValueType volSize[3], + float * translatedProjectionIndexTransformMatrices, + float * translatedVolumeTransformMatrices, + float * dev_proj_in, + float * dev_proj_out, + float * dev_vol, + float t_step, + float * source_positions, + float radiusCylindricalDetector, + float box_min[3], + float box_max[3], + float spacing[3], + unsigned int vectorLength) { // Constant memory - cudaMemcpyToSymbol(c_projSize, projSize, sizeof(int3)); + cudaMemcpyToSymbol(c_projSize, projSize, sizeof(SizeValueType3)); cudaMemcpyToSymbol(c_boxMin, box_min, sizeof(float3)); cudaMemcpyToSymbol(c_boxMax, box_max, sizeof(float3)); cudaMemcpyToSymbol(c_spacing, spacing, sizeof(float3)); - cudaMemcpyToSymbol(c_volSize, volSize, sizeof(int3)); + cudaMemcpyToSymbol(c_volSize, volSize, sizeof(SizeValueType3)); cudaMemcpyToSymbol(c_tStep, &t_step, sizeof(float)); cudaMemcpyToSymbol(c_radius, &radiusCylindricalDetector, sizeof(float)); diff --git a/src/rtkCudaForwardWarpImageFilter.cu b/src/rtkCudaForwardWarpImageFilter.cu index 17e7030f8..8e5c78ec5 100644 --- a/src/rtkCudaForwardWarpImageFilter.cu +++ b/src/rtkCudaForwardWarpImageFilter.cu @@ -53,19 +53,18 @@ __constant__ float c_PPOutputToIndexOutputMatrix[12]; __global__ void fillHoles_3Dgrid(float * dev_vol_out, float * dev_accumulate_weights, int3 out_dim) { - int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - int j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; - int k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; + itk::OffsetValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::OffsetValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; + itk::OffsetValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; if (i >= out_dim.x || j >= out_dim.y || k >= out_dim.z) { return; } - // Index row major into the volume - long int out_idx = i + (j + k * out_dim.y) * (out_dim.x); - long int current_idx; - int radius = 3; + itk::OffsetValueType out_idx = i + (j + k * out_dim.y) * out_dim.x; + itk::OffsetValueType current_idx; + int radius = 3; float eps = 1e-6; @@ -87,7 +86,8 @@ fillHoles_3Dgrid(float * dev_vol_out, float * dev_accumulate_weights, int3 out_d if ((i + delta_i >= 0) && (i + delta_i < out_dim.x) && (j + delta_j >= 0) && (j + delta_j < out_dim.y) && (k + delta_k >= 0) && (k + delta_k < out_dim.z)) { - current_idx = i + delta_i + (j + delta_j + (k + delta_k) * out_dim.y) * (out_dim.x); + itk::OffsetValueType ni = i + delta_i, nj = j + delta_j, nk = k + delta_k; + current_idx = ni + (nj + nk * out_dim.y) * out_dim.x; sum += dev_vol_out[current_idx] * dev_accumulate_weights[current_idx]; sum_weights += dev_accumulate_weights[current_idx]; } @@ -104,9 +104,9 @@ fillHoles_3Dgrid(float * dev_vol_out, float * dev_accumulate_weights, int3 out_d __global__ void normalize_3Dgrid(float * dev_vol_out, float * dev_accumulate_weights, int3 out_dim) { - unsigned int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - unsigned int j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; - unsigned int k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; + itk::SizeValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::SizeValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; + itk::SizeValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; if (i >= out_dim.x || j >= out_dim.y || k >= out_dim.z) { @@ -114,7 +114,7 @@ normalize_3Dgrid(float * dev_vol_out, float * dev_accumulate_weights, int3 out_d } // Index row major into the volume - long int out_idx = i + (j + k * out_dim.y) * (out_dim.x); + itk::SizeValueType out_idx = i + (j + k * out_dim.y) * (out_dim.x); float eps = 1e-6; @@ -134,9 +134,9 @@ linearSplat_3Dgrid(float * dev_vol_in, cudaTextureObject_t tex_ydvf, cudaTextureObject_t tex_zdvf) { - unsigned int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - unsigned int j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; - unsigned int k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; + itk::SizeValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::SizeValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; + itk::SizeValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; if (i >= in_dim.x || j >= in_dim.y || k >= in_dim.z) { @@ -144,8 +144,8 @@ linearSplat_3Dgrid(float * dev_vol_in, } // Index row major into the volume - float3 idx = make_float3(i, j, k); - long int in_idx = i + (j + k * in_dim.y) * (in_dim.x); + float3 idx = make_float3(i, j, k); + itk::SizeValueType in_idx = i + (j + k * in_dim.y) * (in_dim.x); // Matrix multiply to get the index in the DVF texture of the current point in the output volume float3 IndexInDVF = matrix_multiply(idx, c_IndexInputToIndexDVFMatrix); @@ -188,23 +188,17 @@ linearSplat_3Dgrid(float * dev_vol_in, float weight110 = Distance.x * Distance.y * (1 - Distance.z); float weight111 = Distance.x * Distance.y * Distance.z; - // Compute indices in the volume - long int out_idx000 = (BaseIndexInOutput.x + 0) + (BaseIndexInOutput.y + 0) * out_dim.x + - (BaseIndexInOutput.z + 0) * out_dim.x * out_dim.y; - long int out_idx001 = (BaseIndexInOutput.x + 0) + (BaseIndexInOutput.y + 0) * out_dim.x + - (BaseIndexInOutput.z + 1) * out_dim.x * out_dim.y; - long int out_idx010 = (BaseIndexInOutput.x + 0) + (BaseIndexInOutput.y + 1) * out_dim.x + - (BaseIndexInOutput.z + 0) * out_dim.x * out_dim.y; - long int out_idx011 = (BaseIndexInOutput.x + 0) + (BaseIndexInOutput.y + 1) * out_dim.x + - (BaseIndexInOutput.z + 1) * out_dim.x * out_dim.y; - long int out_idx100 = (BaseIndexInOutput.x + 1) + (BaseIndexInOutput.y + 0) * out_dim.x + - (BaseIndexInOutput.z + 0) * out_dim.x * out_dim.y; - long int out_idx101 = (BaseIndexInOutput.x + 1) + (BaseIndexInOutput.y + 0) * out_dim.x + - (BaseIndexInOutput.z + 1) * out_dim.x * out_dim.y; - long int out_idx110 = (BaseIndexInOutput.x + 1) + (BaseIndexInOutput.y + 1) * out_dim.x + - (BaseIndexInOutput.z + 0) * out_dim.x * out_dim.y; - long int out_idx111 = (BaseIndexInOutput.x + 1) + (BaseIndexInOutput.y + 1) * out_dim.x + - (BaseIndexInOutput.z + 1) * out_dim.x * out_dim.y; + // Compute indices in the volume (64-bit to avoid 32-bit overflow of the products) + itk::OffsetValueType bx0 = BaseIndexInOutput.x, by0 = BaseIndexInOutput.y, bz0 = BaseIndexInOutput.z; + itk::OffsetValueType bx1 = bx0 + 1, by1 = by0 + 1, bz1 = bz0 + 1; + itk::OffsetValueType out_idx000 = bx0 + by0 * out_dim.x + bz0 * out_dim.x * out_dim.y; + itk::OffsetValueType out_idx001 = bx0 + by0 * out_dim.x + bz1 * out_dim.x * out_dim.y; + itk::OffsetValueType out_idx010 = bx0 + by1 * out_dim.x + bz0 * out_dim.x * out_dim.y; + itk::OffsetValueType out_idx011 = bx0 + by1 * out_dim.x + bz1 * out_dim.x * out_dim.y; + itk::OffsetValueType out_idx100 = bx1 + by0 * out_dim.x + bz0 * out_dim.x * out_dim.y; + itk::OffsetValueType out_idx101 = bx1 + by0 * out_dim.x + bz1 * out_dim.x * out_dim.y; + itk::OffsetValueType out_idx110 = bx1 + by1 * out_dim.x + bz0 * out_dim.x * out_dim.y; + itk::OffsetValueType out_idx111 = bx1 + by1 * out_dim.x + bz1 * out_dim.x * out_dim.y; // Determine whether they are indeed in the volume bool isInVolume_out_idx000 = (BaseIndexInOutput.x + 0 >= 0) && (BaseIndexInOutput.x + 0 < out_dim.x) && @@ -350,18 +344,18 @@ nearestNeighborSplat_3Dgrid(float * dev_vol_in, /////////////////////////////////////////////////////////////////////////// // FUNCTION: CUDA_ForwardWarp ///////////////////////////// void -CUDA_ForwardWarp(int input_vol_dim[3], - int input_dvf_dim[3], - int output_vol_dim[3], - float IndexInputToPPInputMatrix[12], - float IndexInputToIndexDVFMatrix[12], - float PPOutputToIndexOutputMatrix[12], - float * dev_input_vol, - float * dev_input_xdvf, - float * dev_input_ydvf, - float * dev_input_zdvf, - float * dev_output_vol, - bool isLinear) +CUDA_ForwardWarp(itk::SizeValueType input_vol_dim[3], + itk::SizeValueType input_dvf_dim[3], + itk::SizeValueType output_vol_dim[3], + float IndexInputToPPInputMatrix[12], + float IndexInputToIndexDVFMatrix[12], + float PPOutputToIndexOutputMatrix[12], + float * dev_input_vol, + float * dev_input_xdvf, + float * dev_input_ydvf, + float * dev_input_zdvf, + float * dev_output_vol, + bool isLinear) { // Prepare channel description for arrays static cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(); @@ -386,7 +380,7 @@ CUDA_ForwardWarp(int input_vol_dim[3], /////////////////////////////////////// /// Initialize the output - size_t memorySizeOutput = sizeof(float) * output_vol_dim[0] * output_vol_dim[1] * output_vol_dim[2]; + itk::SizeValueType memorySizeOutput = sizeof(float) * output_vol_dim[0] * output_vol_dim[1] * output_vol_dim[2]; cudaMemset((void *)dev_output_vol, 0, memorySizeOutput); ////////////////////////////////////// diff --git a/src/rtkCudaForwardWarpImageFilter.cxx b/src/rtkCudaForwardWarpImageFilter.cxx index 9b680c383..ae40d1831 100644 --- a/src/rtkCudaForwardWarpImageFilter.cxx +++ b/src/rtkCudaForwardWarpImageFilter.cxx @@ -49,17 +49,17 @@ CudaForwardWarpImageFilter ::GPUGenerateData() } // Cuda convenient format for dimensions - int inputVolumeSize[3]; + itk::SizeValueType inputVolumeSize[3]; inputVolumeSize[0] = this->GetInput(0)->GetBufferedRegion().GetSize()[0]; inputVolumeSize[1] = this->GetInput(0)->GetBufferedRegion().GetSize()[1]; inputVolumeSize[2] = this->GetInput(0)->GetBufferedRegion().GetSize()[2]; - int inputDVFSize[3]; + itk::SizeValueType inputDVFSize[3]; inputDVFSize[0] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[0]; inputDVFSize[1] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[1]; inputDVFSize[2] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[2]; - int outputVolumeSize[3]; + itk::SizeValueType outputVolumeSize[3]; outputVolumeSize[0] = this->GetOutput()->GetBufferedRegion().GetSize()[0]; outputVolumeSize[1] = this->GetOutput()->GetBufferedRegion().GetSize()[1]; outputVolumeSize[2] = this->GetOutput()->GetBufferedRegion().GetSize()[2]; diff --git a/src/rtkCudaInterpolateImageFilter.cu b/src/rtkCudaInterpolateImageFilter.cu index 9760f0119..203250cd2 100644 --- a/src/rtkCudaInterpolateImageFilter.cu +++ b/src/rtkCudaInterpolateImageFilter.cu @@ -28,19 +28,23 @@ #include void -CUDA_interpolation(const int4 & inputSize, float * input, float * output, int projectionNumber, float ** weights) +CUDA_interpolation(const itk::SizeValueType inputSize[4], + float * input, + float * output, + int projectionNumber, + float ** weights) { cublasHandle_t handle; cublasCreate(&handle); // CUDA device pointers - size_t nVoxelsOutput = inputSize.x * inputSize.y * inputSize.z; + size_t nVoxelsOutput = inputSize[0] * inputSize[1] * inputSize[2]; size_t memorySizeOutput = nVoxelsOutput * sizeof(float); // Reset output volume cudaMemset((void *)output, 0, memorySizeOutput); - for (int phase = 0; phase < inputSize.w; phase++) + for (itk::SizeValueType phase = 0; phase < inputSize[3]; phase++) { float weight = weights[phase][projectionNumber]; if (weight != 0) diff --git a/src/rtkCudaInterpolateImageFilter.cxx b/src/rtkCudaInterpolateImageFilter.cxx index f9a69e42a..06f8be1ea 100644 --- a/src/rtkCudaInterpolateImageFilter.cxx +++ b/src/rtkCudaInterpolateImageFilter.cxx @@ -29,11 +29,11 @@ CudaInterpolateImageFilter ::CudaInterpolateImageFilter() = default; void CudaInterpolateImageFilter ::GPUGenerateData() { - int4 inputSize; - inputSize.x = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[0]; - inputSize.y = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[1]; - inputSize.z = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[2]; - inputSize.w = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[3]; + itk::SizeValueType inputSize[4]; + inputSize[0] = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[0]; + inputSize[1] = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[1]; + inputSize[2] = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[2]; + inputSize[3] = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[3]; float * pvolseries = static_cast(this->GetInputVolumeSeries()->GetCudaDataManager()->GetGPUBufferPointer()); float * pvol = static_cast(this->GetOutput()->GetCudaDataManager()->GetGPUBufferPointer()); diff --git a/src/rtkCudaLagCorrectionImageFilter.cu b/src/rtkCudaLagCorrectionImageFilter.cu index a3ae51a15..c114cc419 100644 --- a/src/rtkCudaLagCorrectionImageFilter.cu +++ b/src/rtkCudaLagCorrectionImageFilter.cu @@ -37,25 +37,23 @@ kernel_lag_correction(int3 proj_idx_in, { constexpr int modelOrder = 4; - // compute thread index - int3 tIdx; - tIdx.x = blockIdx.x * blockDim.x + threadIdx.x; - tIdx.y = blockIdx.y * blockDim.y + threadIdx.y; - tIdx.z = blockIdx.z * blockDim.z + threadIdx.z; - long int tIdx_comp = tIdx.x + tIdx.y * proj_size_out.x + tIdx.z * proj_size_out_buf.x * proj_size_out_buf.y; + // compute thread index (64-bit to avoid 32-bit overflow of combined indices) + itk::SizeValueType tIdx_x = blockIdx.x * blockDim.x + threadIdx.x; + itk::SizeValueType tIdx_y = blockIdx.y * blockDim.y + threadIdx.y; + itk::SizeValueType tIdx_z = blockIdx.z * blockDim.z + threadIdx.z; + itk::SizeValueType tIdx_comp = tIdx_x + tIdx_y * proj_size_out.x + tIdx_z * proj_size_out_buf.x * proj_size_out_buf.y; // check if outside of projection grid - if (tIdx.x >= proj_size_out.x || tIdx.y >= proj_size_out.y || tIdx.z >= proj_size_out.z) + if (tIdx_x >= proj_size_out.x || tIdx_y >= proj_size_out.y || tIdx_z >= proj_size_out.z) return; - // compute projection index from thread index - int3 pIdx = make_int3(tIdx.x + proj_idx_out.x, tIdx.y + proj_idx_out.y, tIdx.z + proj_idx_out.z); - // combined proj. index -> use thread index in z because accessing memory only with this index - long int pIdx_comp = (pIdx.x - proj_idx_in.x) + (pIdx.y - proj_idx_in.y) * proj_size_in_buf.x + - (pIdx.z - proj_idx_in.z) * proj_size_in_buf.x * proj_size_in_buf.y; + // combined projection index (arithmetic in 64-bit to avoid overflow of the products) + itk::SizeValueType pIdx_comp = (tIdx_x - proj_idx_in.x + proj_idx_out.x) + + (tIdx_y - proj_idx_in.y + proj_idx_out.y) * proj_size_in_buf.x + + (tIdx_z - proj_idx_in.z + proj_idx_out.z) * proj_size_in_buf.x * proj_size_in_buf.y; - long int sIdx_comp = tIdx.x + tIdx.y * proj_size_out.x; - unsigned idx_s = sIdx_comp * modelOrder; + itk::SizeValueType sIdx_comp = tIdx_x + tIdx_y * proj_size_out.x; + itk::SizeValueType idx_s = sIdx_comp * modelOrder; float yk = static_cast(dev_proj_in[pIdx_comp]); float xk = yk; diff --git a/src/rtkCudaLastDimensionTVDenoisingImageFilter.cu b/src/rtkCudaLastDimensionTVDenoisingImageFilter.cu index 5a30b0cef..d89483d4b 100644 --- a/src/rtkCudaLastDimensionTVDenoisingImageFilter.cu +++ b/src/rtkCudaLastDimensionTVDenoisingImageFilter.cu @@ -45,8 +45,9 @@ denoise_oneD_TV_kernel(float * in, float * out, float beta, float gamma, int nit float * gradient = &shared[3 * c_Size.w]; // Each thread reads one element into the shared buffer - long int gindex = ((threadIdx.x * c_Size.z + blockIdx.z) * c_Size.y + blockIdx.y) * c_Size.x + blockIdx.x; - int lindex = threadIdx.x; + itk::SizeValueType tx = threadIdx.x; + itk::SizeValueType gindex = ((tx * c_Size.z + blockIdx.z) * c_Size.y + blockIdx.y) * c_Size.x + blockIdx.x; + int lindex = threadIdx.x; input[lindex] = in[gindex]; __syncthreads(); diff --git a/src/rtkCudaParkerShortScanImageFilter.cu b/src/rtkCudaParkerShortScanImageFilter.cu index be4656371..c3af49606 100644 --- a/src/rtkCudaParkerShortScanImageFilter.cu +++ b/src/rtkCudaParkerShortScanImageFilter.cu @@ -57,26 +57,25 @@ kernel_parker_weight(int2 proj_idx, cudaTextureObject_t tex_geom // geometry texture object ) { - // compute projection index (== thread index) - int3 pIdx; - pIdx.x = blockIdx.x * blockDim.x + threadIdx.x; - pIdx.y = blockIdx.y * blockDim.y + threadIdx.y; - pIdx.z = blockIdx.z * blockDim.z + threadIdx.z; - long int pIdx_comp_in = pIdx.x + (pIdx.y + pIdx.z * proj_size_buf_in.y) * (proj_size_buf_in.x); - long int pIdx_comp_out = pIdx.x + (pIdx.y + pIdx.z * proj_size_buf_out.y) * (proj_size_buf_out.x); + // compute projection index (== thread index), 64-bit to avoid 32-bit overflow of combined indices + itk::SizeValueType pIdx_x = blockIdx.x * blockDim.x + threadIdx.x; + itk::SizeValueType pIdx_y = blockIdx.y * blockDim.y + threadIdx.y; + itk::SizeValueType pIdx_z = blockIdx.z * blockDim.z + threadIdx.z; + itk::SizeValueType pIdx_comp_in = pIdx_x + (pIdx_y + pIdx_z * proj_size_buf_in.y) * proj_size_buf_in.x; + itk::SizeValueType pIdx_comp_out = pIdx_x + (pIdx_y + pIdx_z * proj_size_buf_out.y) * proj_size_buf_out.x; // check if outside of projection grid - if (pIdx.x >= proj_size.x || pIdx.y >= proj_size.y || pIdx.z >= proj_size.z) + if (pIdx_x >= proj_size.x || pIdx_y >= proj_size.y || pIdx_z >= proj_size.z) return; - float sdd = tex1Dfetch(tex_geom, pIdx.z * 5 + 0); - float sx = tex1Dfetch(tex_geom, pIdx.z * 5 + 1); - float px = tex1Dfetch(tex_geom, pIdx.z * 5 + 2); - float sid = tex1Dfetch(tex_geom, pIdx.z * 5 + 3); + float sdd = tex1Dfetch(tex_geom, pIdx_z * 5 + 0); + float sx = tex1Dfetch(tex_geom, pIdx_z * 5 + 1); + float px = tex1Dfetch(tex_geom, pIdx_z * 5 + 2); + float sid = tex1Dfetch(tex_geom, pIdx_z * 5 + 3); // convert actual index to point float pPoint = - TransformIndexToPhysicalPoint(make_int2(pIdx.x + proj_idx.x, pIdx.y + proj_idx.y), proj_orig, proj_row, proj_col); + TransformIndexToPhysicalPoint(make_int2(pIdx_x + proj_idx.x, pIdx_y + proj_idx.y), proj_orig, proj_row, proj_col); // alpha projection angle float hyp = sqrtf(sid * sid + sx * sx); // to untilted situation @@ -85,7 +84,7 @@ kernel_parker_weight(int2 proj_idx, float alpha = atan(-1 * l * invsid); // beta projection angle: Parker's article assumes that the scan starts at 0 - float beta = tex1Dfetch(tex_geom, pIdx.z * 5 + 4); + float beta = tex1Dfetch(tex_geom, pIdx_z * 5 + 4); beta -= firstAngle; if (beta < 0) beta += (2.f * CUDART_PI_F); diff --git a/src/rtkCudaPolynomialGainCorrectionImageFilter.cu b/src/rtkCudaPolynomialGainCorrectionImageFilter.cu index 2955a3056..53b473d51 100644 --- a/src/rtkCudaPolynomialGainCorrectionImageFilter.cu +++ b/src/rtkCudaPolynomialGainCorrectionImageFilter.cu @@ -37,39 +37,37 @@ kernel_gain_correction(int3 proj_idx_in, float * dev_gain_in, float * powerlut) { - // compute thread index - int3 tIdx; - tIdx.x = blockIdx.x * blockDim.x + threadIdx.x; - tIdx.y = blockIdx.y * blockDim.y + threadIdx.y; - tIdx.z = blockIdx.z * blockDim.z + threadIdx.z; - long int tIdx_comp = tIdx.x + tIdx.y * proj_size_out.x + tIdx.z * proj_size_out_buf.x * proj_size_out_buf.y; + // compute thread index (64-bit to avoid 32-bit overflow of combined indices) + itk::SizeValueType tIdx_x = blockIdx.x * blockDim.x + threadIdx.x; + itk::SizeValueType tIdx_y = blockIdx.y * blockDim.y + threadIdx.y; + itk::SizeValueType tIdx_z = blockIdx.z * blockDim.z + threadIdx.z; + itk::SizeValueType tIdx_comp = tIdx_x + tIdx_y * proj_size_out.x + tIdx_z * proj_size_out_buf.x * proj_size_out_buf.y; // check if outside of projection grid - if (tIdx.x >= proj_size_out.x || tIdx.y >= proj_size_out.y || tIdx.z >= proj_size_out.z) + if (tIdx_x >= proj_size_out.x || tIdx_y >= proj_size_out.y || tIdx_z >= proj_size_out.z) return; - // compute projection index from thread index - int3 pIdx = make_int3(tIdx.x + proj_idx_out.x, tIdx.y + proj_idx_out.y, tIdx.z + proj_idx_out.z); - // combined proj. index -> use thread index in z because accessing memory only with this index - long int pIdx_comp = (pIdx.x - proj_idx_in.x) + (pIdx.y - proj_idx_in.y) * proj_size_in_buf.x + - (pIdx.z - proj_idx_in.z) * proj_size_in_buf.x * proj_size_in_buf.y; + // combined projection index (arithmetic in 64-bit to avoid overflow of the products) + itk::SizeValueType pIdx_comp = (tIdx_x - proj_idx_in.x + proj_idx_out.x) + + (tIdx_y - proj_idx_in.y + proj_idx_out.y) * proj_size_in_buf.x + + (tIdx_z - proj_idx_in.z + proj_idx_out.z) * proj_size_in_buf.x * proj_size_in_buf.y; int modelOrder = static_cast(cst_coef[0]); - long int sIdx_comp = tIdx.x + tIdx.y * proj_size_out.x; // in-slice index + itk::SizeValueType sIdx_comp = tIdx_x + tIdx_y * proj_size_out.x; // in-slice index // Correct for dark field unsigned short xk = 0; if (dev_proj_in[pIdx_comp] > dev_dark_in[sIdx_comp]) xk = dev_proj_in[pIdx_comp] - dev_dark_in[sIdx_comp]; - float yk = 0.f; - int lutidx = xk * modelOrder; // index to powerlut - int projsize = proj_size_in.x * proj_size_in.y; + float yk = 0.f; + itk::SizeValueType lutidx = static_cast(xk) * modelOrder; // index to powerlut + itk::SizeValueType projsize = static_cast(proj_size_in.x) * proj_size_in.y; for (int n = 0; n < modelOrder; n++) { - int gainidx = n * projsize + sIdx_comp; - float gainM = dev_gain_in[gainidx]; + itk::SizeValueType gainidx = n * projsize + sIdx_comp; + float gainM = dev_gain_in[gainidx]; yk += gainM * powerlut[lutidx + n]; } diff --git a/src/rtkCudaRayCastBackProjectionImageFilter.cu b/src/rtkCudaRayCastBackProjectionImageFilter.cu index 8898230ec..00ebb534b 100644 --- a/src/rtkCudaRayCastBackProjectionImageFilter.cu +++ b/src/rtkCudaRayCastBackProjectionImageFilter.cu @@ -55,10 +55,10 @@ __constant__ float c_sourcePos[SLAB_SIZE * 3]; // Can process stacks of __global__ void kernel_ray_cast_back_project(float * dev_vol_out, float * dev_proj) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int numThread = j * c_projSize.x + i; - unsigned int proj = blockIdx.z * blockDim.z + threadIdx.z; + itk::SizeValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::SizeValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::SizeValueType numThread = j * c_projSize.x + i; + itk::SizeValueType proj = blockIdx.z * blockDim.z + threadIdx.z; if (i >= c_projSize.x || j >= c_projSize.y || proj >= c_projSize.z) return; @@ -100,10 +100,10 @@ kernel_ray_cast_back_project(float * dev_vol_out, float * dev_proj) float mmStep = vStep * dirLengthInMM; // Skip rays intersecting less half a step - float toSplat; - long int indices[8]; - float weights[8]; - int3 floor_pos; + float toSplat; + itk::SizeValueType indices[8]; + float weights[8]; + int3 floor_pos; // First position in the box float halfVStep = 0.5f * vStep; @@ -144,18 +144,27 @@ kernel_ray_cast_back_project(float * dev_vol_out, float * dev_proj) pos_high.y = min(floor_pos.y + 1, c_volSize.y - 1); pos_high.z = min(floor_pos.z + 1, c_volSize.z - 1); - // Compute indices in the volume - indices[0] = pos_low.x + pos_low.y * c_volSize.x + pos_low.z * c_volSize.x * c_volSize.y; // index000 - indices[1] = pos_low.x + pos_low.y * c_volSize.x + pos_high.z * c_volSize.x * c_volSize.y; // index001 - indices[2] = pos_low.x + pos_high.y * c_volSize.x + pos_low.z * c_volSize.x * c_volSize.y; // index010 - indices[3] = pos_low.x + pos_high.y * c_volSize.x + pos_high.z * c_volSize.x * c_volSize.y; // index011 - indices[4] = pos_high.x + pos_low.y * c_volSize.x + pos_low.z * c_volSize.x * c_volSize.y; // index100 - indices[5] = pos_high.x + pos_low.y * c_volSize.x + pos_high.z * c_volSize.x * c_volSize.y; // index101 - indices[6] = pos_high.x + pos_high.y * c_volSize.x + pos_low.z * c_volSize.x * c_volSize.y; // index110 - indices[7] = pos_high.x + pos_high.y * c_volSize.x + pos_high.z * c_volSize.x * c_volSize.y; // index111 + // Compute indices in the volume (using 64-bit arithmetic to avoid overflow on large volumes) + const itk::SizeValueType volSize_X = c_volSize.x; + const itk::SizeValueType volSize_XY = volSize_X * c_volSize.y; + const itk::SizeValueType pxl = pos_low.x; + const itk::SizeValueType pxh = pos_high.x; + const itk::SizeValueType pyl = pos_low.y; + const itk::SizeValueType pyh = pos_high.y; + const itk::SizeValueType pzl = pos_low.z; + const itk::SizeValueType pzh = pos_high.z; + indices[0] = pxl + pyl * volSize_X + pzl * volSize_XY; // index000 + indices[1] = pxl + pyl * volSize_X + pzh * volSize_XY; // index001 + indices[2] = pxl + pyh * volSize_X + pzl * volSize_XY; // index010 + indices[3] = pxl + pyh * volSize_X + pzh * volSize_XY; // index011 + indices[4] = pxh + pyl * volSize_X + pzl * volSize_XY; // index100 + indices[5] = pxh + pyl * volSize_X + pzh * volSize_XY; // index101 + indices[6] = pxh + pyh * volSize_X + pzl * volSize_XY; // index110 + indices[7] = pxh + pyh * volSize_X + pzh * volSize_XY; // index111 // Compute the value to be splatted - toSplat = dev_proj[numThread + proj * c_projSize.x * c_projSize.y] * mmStep; + itk::SizeValueType projOffset = numThread + proj * c_projSize.x * c_projSize.y; + toSplat = dev_proj[projOffset] * mmStep; atomicAdd(&dev_vol_out[indices[0]], toSplat * weights[0]); atomicAdd(&dev_vol_out[indices[1]], toSplat * weights[1]); atomicAdd(&dev_vol_out[indices[2]], toSplat * weights[2]); @@ -170,7 +179,8 @@ kernel_ray_cast_back_project(float * dev_vol_out, float * dev_proj) } // Last position - toSplat = dev_proj[numThread + proj * c_projSize.x * c_projSize.y] * (tfar - t + halfVStep) * dirLengthInMM; + itk::SizeValueType projOffsetLast = numThread + proj * c_projSize.x * c_projSize.y; + toSplat = dev_proj[projOffsetLast] * (tfar - t + halfVStep) * dirLengthInMM; atomicAdd(&dev_vol_out[indices[0]], toSplat * weights[0]); atomicAdd(&dev_vol_out[indices[1]], toSplat * weights[1]); atomicAdd(&dev_vol_out[indices[2]], toSplat * weights[2]); diff --git a/src/rtkCudaRayCastBackProjectionImageFilter.cxx b/src/rtkCudaRayCastBackProjectionImageFilter.cxx index c2c638bb0..efbb1c762 100644 --- a/src/rtkCudaRayCastBackProjectionImageFilter.cxx +++ b/src/rtkCudaRayCastBackProjectionImageFilter.cxx @@ -42,10 +42,10 @@ CudaRayCastBackProjectionImageFilter ::GPUGenerateData() { itkGenericExceptionMacro(<< "Error, ThreeDCircularProjectionGeometry expected"); } - constexpr unsigned int Dimension = 3; - const unsigned int iFirstProj = this->GetInput(1)->GetRequestedRegion().GetIndex(Dimension - 1); - const unsigned int nProj = this->GetInput(1)->GetRequestedRegion().GetSize(Dimension - 1); - const unsigned int nPixelsPerProj = + constexpr unsigned int Dimension = 3; + const unsigned int iFirstProj = this->GetInput(1)->GetRequestedRegion().GetIndex(Dimension - 1); + const unsigned int nProj = this->GetInput(1)->GetRequestedRegion().GetSize(Dimension - 1); + const itk::SizeValueType nPixelsPerProj = this->GetInput(1)->GetBufferedRegion().GetSize(0) * this->GetInput(1)->GetBufferedRegion().GetSize(1); itk::Vector source_position; @@ -155,7 +155,7 @@ CudaRayCastBackProjectionImageFilter ::GPUGenerateData() source_positions[(iProj - iFirstProj) * 3 + d] = source_position[d]; // Ignore the 4th component } - int projectionOffset = 0; + itk::OffsetValueType projectionOffset = 0; for (unsigned int i = 0; i < nProj; i += SLAB_SIZE) { // If nProj is not a multiple of SLAB_SIZE, the last slab will contain less than SLAB_SIZE projections diff --git a/src/rtkCudaSplatImageFilter.cu b/src/rtkCudaSplatImageFilter.cu index 726aa3e43..d3710042b 100644 --- a/src/rtkCudaSplatImageFilter.cu +++ b/src/rtkCudaSplatImageFilter.cu @@ -28,14 +28,18 @@ #include void -CUDA_splat(const int4 & outputSize, float * input, float * output, int projectionNumber, float ** weights) +CUDA_splat(const itk::SizeValueType outputSize[4], + float * input, + float * output, + int projectionNumber, + float ** weights) { cublasHandle_t handle; cublasCreate(&handle); - size_t numel = outputSize.x * outputSize.y * outputSize.z; + size_t numel = outputSize[0] * outputSize[1] * outputSize[2]; - for (int phase = 0; phase < outputSize.w; phase++) + for (itk::SizeValueType phase = 0; phase < outputSize[3]; phase++) { float weight = weights[phase][projectionNumber]; if (weight != 0) diff --git a/src/rtkCudaSplatImageFilter.cxx b/src/rtkCudaSplatImageFilter.cxx index 14d198df8..bff171f3a 100644 --- a/src/rtkCudaSplatImageFilter.cxx +++ b/src/rtkCudaSplatImageFilter.cxx @@ -29,11 +29,11 @@ CudaSplatImageFilter ::CudaSplatImageFilter() = default; void CudaSplatImageFilter ::GPUGenerateData() { - int4 outputSize; - outputSize.x = this->GetOutput()->GetLargestPossibleRegion().GetSize()[0]; - outputSize.y = this->GetOutput()->GetLargestPossibleRegion().GetSize()[1]; - outputSize.z = this->GetOutput()->GetLargestPossibleRegion().GetSize()[2]; - outputSize.w = this->GetOutput()->GetLargestPossibleRegion().GetSize()[3]; + itk::SizeValueType outputSize[4]; + outputSize[0] = this->GetOutput()->GetLargestPossibleRegion().GetSize()[0]; + outputSize[1] = this->GetOutput()->GetLargestPossibleRegion().GetSize()[1]; + outputSize[2] = this->GetOutput()->GetLargestPossibleRegion().GetSize()[2]; + outputSize[3] = this->GetOutput()->GetLargestPossibleRegion().GetSize()[3]; float * pvolseries = static_cast(this->GetOutput()->GetCudaDataManager()->GetGPUBufferPointer()); float * pvol = static_cast(this->GetInputVolume()->GetCudaDataManager()->GetGPUBufferPointer()); diff --git a/src/rtkCudaTotalVariationDenoisingBPDQImageFilter.cu b/src/rtkCudaTotalVariationDenoisingBPDQImageFilter.cu index db32815a4..e76ec13a7 100644 --- a/src/rtkCudaTotalVariationDenoisingBPDQImageFilter.cu +++ b/src/rtkCudaTotalVariationDenoisingBPDQImageFilter.cu @@ -38,14 +38,14 @@ __constant__ float3 c_Spacing; __global__ void magnitude_threshold_kernel(float * grad_x, float * grad_y, float * grad_z, float gamma) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int k = blockIdx.z * blockDim.z + threadIdx.z; + itk::SizeValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::SizeValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::SizeValueType k = blockIdx.z * blockDim.z + threadIdx.z; if (i >= c_Size.x || j >= c_Size.y || k >= c_Size.z) return; - long int id = (k * c_Size.y + j) * c_Size.x + i; + itk::SizeValueType id = (k * c_Size.y + j) * c_Size.x + i; float norm = sqrt(grad_x[id] * grad_x[id] + grad_y[id] * grad_y[id] + grad_z[id] * grad_z[id]); if (norm > gamma) @@ -60,17 +60,17 @@ magnitude_threshold_kernel(float * grad_x, float * grad_y, float * grad_z, float __global__ void gradient_and_subtract_kernel(float * in, float * grad_x, float * grad_y, float * grad_z) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int k = blockIdx.z * blockDim.z + threadIdx.z; + itk::SizeValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::SizeValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::SizeValueType k = blockIdx.z * blockDim.z + threadIdx.z; if (i >= c_Size.x || j >= c_Size.y || k >= c_Size.z) return; - long int id = (k * c_Size.y + j) * c_Size.x + i; - long int id_x = (k * c_Size.y + j) * c_Size.x + i + 1; - long int id_y = (k * c_Size.y + j + 1) * c_Size.x + i; - long int id_z = ((k + 1) * c_Size.y + j) * c_Size.x + i; + itk::SizeValueType id = (k * c_Size.y + j) * c_Size.x + i; + itk::SizeValueType id_x = (k * c_Size.y + j) * c_Size.x + i + 1; + itk::SizeValueType id_y = (k * c_Size.y + j + 1) * c_Size.x + i; + itk::SizeValueType id_z = ((k + 1) * c_Size.y + j) * c_Size.x + i; if (i != (c_Size.x - 1)) grad_x[id] -= ((in[id_x] - in[id]) / c_Spacing.x); @@ -83,14 +83,14 @@ gradient_and_subtract_kernel(float * in, float * grad_x, float * grad_y, float * __global__ void multiply_by_beta_kernel(float * input, float * output, float beta) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int k = blockIdx.z * blockDim.z + threadIdx.z; + itk::SizeValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::SizeValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::SizeValueType k = blockIdx.z * blockDim.z + threadIdx.z; if (i >= c_Size.x || j >= c_Size.y || k >= c_Size.z) return; - long int id = (k * c_Size.y + j) * c_Size.x + i; + itk::SizeValueType id = (k * c_Size.y + j) * c_Size.x + i; output[id] = input[id] * beta; } @@ -98,14 +98,14 @@ multiply_by_beta_kernel(float * input, float * output, float beta) __global__ void subtract_kernel(float * in1, float * in2, float * out) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int k = blockIdx.z * blockDim.z + threadIdx.z; + itk::SizeValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::SizeValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::SizeValueType k = blockIdx.z * blockDim.z + threadIdx.z; if (i >= c_Size.x || j >= c_Size.y || k >= c_Size.z) return; - long int id = (k * c_Size.y + j) * c_Size.x + i; + itk::SizeValueType id = (k * c_Size.y + j) * c_Size.x + i; out[id] = in1[id] - in2[id]; } diff --git a/src/rtkCudaUtilities.cu b/src/rtkCudaUtilities.cu index 753d94cd5..780d61881 100644 --- a/src/rtkCudaUtilities.cu +++ b/src/rtkCudaUtilities.cu @@ -83,7 +83,7 @@ GetFreeGPUGlobalMemory(int device) } __host__ void -prepareScalarTextureObject(int size[3], +prepareScalarTextureObject(itk::SizeValueType size[3], float * dev_ptr, cudaArray *& threeDArray, cudaTextureObject_t & tex, @@ -136,7 +136,7 @@ prepareScalarTextureObject(int size[3], } __host__ void -prepareVectorTextureObject(int size[3], +prepareVectorTextureObject(itk::SizeValueType size[3], const float * dev_ptr, std::vector & componentArrays, const unsigned int nComponents, diff --git a/src/rtkCudaWarpBackProjectionImageFilter.cu b/src/rtkCudaWarpBackProjectionImageFilter.cu index 778cd7c46..a73f4b629 100644 --- a/src/rtkCudaWarpBackProjectionImageFilter.cu +++ b/src/rtkCudaWarpBackProjectionImageFilter.cu @@ -41,14 +41,14 @@ #include // CONSTANTS ////////////////////////////////////////////////////////////// -__constant__ float c_matrices[SLAB_SIZE * 12]; // Can process stacks of at most SLAB_SIZE projections -__constant__ float c_volIndexToProjPP[SLAB_SIZE * 12]; -__constant__ float c_projPPToProjIndex[9]; -__constant__ int3 c_projSize; -__constant__ int3 c_volSize; -__constant__ float c_IndexInputToIndexDVFMatrix[12]; -__constant__ float c_PPInputToIndexInputMatrix[12]; -__constant__ float c_IndexInputToPPInputMatrix[12]; +__constant__ float c_matrices[SLAB_SIZE * 12]; // Can process stacks of at most SLAB_SIZE projections +__constant__ float c_volIndexToProjPP[SLAB_SIZE * 12]; +__constant__ float c_projPPToProjIndex[9]; +__constant__ SizeValueType3 c_projSize; +__constant__ SizeValueType3 c_volSize; +__constant__ float c_IndexInputToIndexDVFMatrix[12]; +__constant__ float c_PPInputToIndexInputMatrix[12]; +__constant__ float c_IndexInputToPPInputMatrix[12]; //////////////////////////////////////////////////////////////////////////// //_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_ @@ -64,9 +64,9 @@ kernel_warp_back_project_3Dgrid(float * dev_vol_in, cudaTextureObject_t tex_zdvf, cudaTextureObject_t tex_proj) { - unsigned int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - unsigned int j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; - unsigned int k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; + itk::SizeValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::SizeValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; + itk::SizeValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; if (i >= c_volSize.x || j >= c_volSize.y || k >= c_volSize.z) { @@ -74,7 +74,7 @@ kernel_warp_back_project_3Dgrid(float * dev_vol_in, } // Index row major into the volume - long int vol_idx = i + (j + k * c_volSize.y) * (c_volSize.x); + itk::SizeValueType vol_idx = i + (j + k * c_volSize.y) * c_volSize.x; float3 IndexInDVF, Displacement, PP, IndexInInput, ip; float voxel_data = 0; @@ -120,9 +120,9 @@ kernel_warp_back_project_3Dgrid_cylindrical_detector(float * dev_vol cudaTextureObject_t tex_zdvf, cudaTextureObject_t tex_proj) { - unsigned int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - unsigned int j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; - unsigned int k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; + itk::SizeValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::SizeValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; + itk::SizeValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; if (i >= c_volSize.x || j >= c_volSize.y || k >= c_volSize.z) { @@ -130,7 +130,7 @@ kernel_warp_back_project_3Dgrid_cylindrical_detector(float * dev_vol } // Index row major into the volume - long int vol_idx = i + (j + k * c_volSize.y) * (c_volSize.x); + itk::SizeValueType vol_idx = i + (j + k * c_volSize.y) * c_volSize.x; float3 IndexInDVF, Displacement, PP, IndexInInput, ip, pp; float voxel_data = 0; @@ -183,24 +183,24 @@ kernel_warp_back_project_3Dgrid_cylindrical_detector(float * dev_vol /////////////////////////////////////////////////////////////////////////// // FUNCTION: CUDA_back_project ///////////////////////////// void -CUDA_warp_back_project(int projSize[3], - int volSize[3], - int dvf_size[3], - float * matrices, - float * volIndexToProjPPs, - float * projPPToProjIndex, - float * dev_vol_in, - float * dev_vol_out, - float * dev_proj, - float * dev_input_dvf, - float IndexInputToIndexDVFMatrix[12], - float PPInputToIndexInputMatrix[12], - float IndexInputToPPInputMatrix[12], - double radiusCylindricalDetector) +CUDA_warp_back_project(itk::SizeValueType projSize[3], + itk::SizeValueType volSize[3], + itk::SizeValueType dvf_size[3], + float * matrices, + float * volIndexToProjPPs, + float * projPPToProjIndex, + float * dev_vol_in, + float * dev_vol_out, + float * dev_proj, + float * dev_input_dvf, + float IndexInputToIndexDVFMatrix[12], + float PPInputToIndexInputMatrix[12], + float IndexInputToPPInputMatrix[12], + double radiusCylindricalDetector) { // Copy the size of inputs into constant memory - cudaMemcpyToSymbol(c_projSize, projSize, sizeof(int3)); - cudaMemcpyToSymbol(c_volSize, volSize, sizeof(int3)); + cudaMemcpyToSymbol(c_projSize, projSize, sizeof(SizeValueType3)); + cudaMemcpyToSymbol(c_volSize, volSize, sizeof(SizeValueType3)); // Copy the projection matrices into constant memory cudaMemcpyToSymbol(c_matrices, &(matrices[0]), 12 * sizeof(float) * projSize[2]); diff --git a/src/rtkCudaWarpBackProjectionImageFilter.cxx b/src/rtkCudaWarpBackProjectionImageFilter.cxx index 5478de8ff..889e4bb73 100644 --- a/src/rtkCudaWarpBackProjectionImageFilter.cxx +++ b/src/rtkCudaWarpBackProjectionImageFilter.cxx @@ -146,16 +146,16 @@ CudaWarpBackProjectionImageFilter ::GPUGenerateData() } // Cuda convenient format for dimensions - int projectionSize[3]; + itk::SizeValueType projectionSize[3]; projectionSize[0] = this->GetInputProjectionStack()->GetBufferedRegion().GetSize()[0]; projectionSize[1] = this->GetInputProjectionStack()->GetBufferedRegion().GetSize()[1]; - int volumeSize[3]; + itk::SizeValueType volumeSize[3]; volumeSize[0] = this->GetOutput()->GetBufferedRegion().GetSize()[0]; volumeSize[1] = this->GetOutput()->GetBufferedRegion().GetSize()[1]; volumeSize[2] = this->GetOutput()->GetBufferedRegion().GetSize()[2]; - int inputDVFSize[3]; + itk::SizeValueType inputDVFSize[3]; inputDVFSize[0] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[0]; inputDVFSize[1] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[1]; inputDVFSize[2] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[2]; diff --git a/src/rtkCudaWarpForwardProjectionImageFilter.cu b/src/rtkCudaWarpForwardProjectionImageFilter.cu index 4b73e135f..1f6b6f86d 100644 --- a/src/rtkCudaWarpForwardProjectionImageFilter.cu +++ b/src/rtkCudaWarpForwardProjectionImageFilter.cu @@ -39,14 +39,14 @@ #include // CONSTANTS -__constant__ int3 c_projSize; -__constant__ float3 c_boxMin; -__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__ SizeValueType3 c_projSize; +__constant__ float3 c_boxMin; +__constant__ float3 c_boxMax; +__constant__ float3 c_spacing; +__constant__ SizeValueType3 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_IndexInputToPPInputMatrix[12]; __constant__ float c_IndexInputToIndexDVFMatrix[12]; @@ -66,9 +66,9 @@ kernel_warped_forwardProject(float * dev_proj_in, cudaTextureObject_t tex_zdvf, cudaTextureObject_t tex_vol) { - unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; - unsigned int j = blockIdx.y * blockDim.y + threadIdx.y; - unsigned int numThread = j * c_projSize.x + i; + itk::SizeValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::SizeValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::SizeValueType numThread = j * c_projSize.x + i; if (i >= c_projSize.x || j >= c_projSize.y) return; @@ -78,7 +78,7 @@ kernel_warped_forwardProject(float * dev_proj_in, float3 pixelPos; float tnear, tfar; - for (unsigned int proj = 0; proj < c_projSize.z; proj++) + for (itk::SizeValueType 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]); @@ -91,8 +91,8 @@ kernel_warped_forwardProject(float * dev_proj_in, // Detect intersection with box if (!intersectBox(ray, &tnear, &tfar, c_boxMin, c_boxMax) || tfar < 0.f) { - dev_proj_out[numThread + proj * c_projSize.x * c_projSize.y] = - dev_proj_in[numThread + proj * c_projSize.x * c_projSize.y]; + itk::OffsetValueType projOffset = numThread + proj * c_projSize.x * c_projSize.y; + dev_proj_out[projOffset] = dev_proj_in[projOffset]; } else { @@ -139,9 +139,8 @@ kernel_warped_forwardProject(float * dev_proj_in, sum += sample; pos += step; } - dev_proj_out[numThread + proj * c_projSize.x * c_projSize.y] = - dev_proj_in[numThread + proj * c_projSize.x * c_projSize.y] + - (sum + (tfar - t + halfVStep) / vStep * sample) * c_tStep; + itk::OffsetValueType projOffset = numThread + proj * c_projSize.x * c_projSize.y; + dev_proj_out[projOffset] = dev_proj_in[projOffset] + (sum + (tfar - t + halfVStep) / vStep * sample) * c_tStep; } } } @@ -154,29 +153,29 @@ kernel_warped_forwardProject(float * dev_proj_in, /////////////////////////////////////////////////////////////////////////// // FUNCTION: CUDA_forward_project() ////////////////////////////////// void -CUDA_warp_forward_project(int projSize[3], - int volSize[3], - int dvfSize[3], - float * matrices, - float * dev_proj_in, - float * dev_proj_out, - float * dev_vol, - float t_step, - float * source_positions, - float box_min[3], - float box_max[3], - float spacing[3], - float * dev_input_dvf, - float IndexInputToIndexDVFMatrix[12], - float PPInputToIndexInputMatrix[12], - float IndexInputToPPInputMatrix[12]) +CUDA_warp_forward_project(itk::SizeValueType projSize[3], + itk::SizeValueType volSize[3], + itk::SizeValueType dvfSize[3], + float * matrices, + float * dev_proj_in, + float * dev_proj_out, + float * dev_vol, + float t_step, + float * source_positions, + float box_min[3], + float box_max[3], + float spacing[3], + float * dev_input_dvf, + float IndexInputToIndexDVFMatrix[12], + float PPInputToIndexInputMatrix[12], + float IndexInputToPPInputMatrix[12]) { // constant memory - cudaMemcpyToSymbol(c_projSize, projSize, sizeof(int3)); + cudaMemcpyToSymbol(c_projSize, projSize, sizeof(SizeValueType3)); cudaMemcpyToSymbol(c_boxMin, box_min, sizeof(float3)); cudaMemcpyToSymbol(c_boxMax, box_max, sizeof(float3)); cudaMemcpyToSymbol(c_spacing, spacing, sizeof(float3)); - cudaMemcpyToSymbol(c_volSize, volSize, sizeof(int3)); + cudaMemcpyToSymbol(c_volSize, volSize, sizeof(SizeValueType3)); cudaMemcpyToSymbol(c_tStep, &t_step, sizeof(float)); // Copy the source position matrix into a float3 in constant memory diff --git a/src/rtkCudaWarpForwardProjectionImageFilter.cxx b/src/rtkCudaWarpForwardProjectionImageFilter.cxx index 1756ec964..b0dceb103 100644 --- a/src/rtkCudaWarpForwardProjectionImageFilter.cxx +++ b/src/rtkCudaWarpForwardProjectionImageFilter.cxx @@ -145,9 +145,9 @@ CudaWarpForwardProjectionImageFilter ::GPUGenerateData() const Superclass::GeometryType * geometry = this->GetGeometry(); 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 = + const unsigned int iFirstProj = this->GetInputProjectionStack()->GetRequestedRegion().GetIndex(Dimension - 1); + const unsigned int nProj = this->GetInputProjectionStack()->GetRequestedRegion().GetSize(Dimension - 1); + const itk::SizeValueType nPixelsPerProj = this->GetOutput()->GetBufferedRegion().GetSize(0) * this->GetOutput()->GetBufferedRegion().GetSize(1); itk::Vector source_position; @@ -172,16 +172,16 @@ CudaWarpForwardProjectionImageFilter ::GPUGenerateData() } // Cuda convenient format for dimensions - int projectionSize[3]; + itk::SizeValueType projectionSize[3]; projectionSize[0] = this->GetOutput()->GetBufferedRegion().GetSize()[0]; projectionSize[1] = this->GetOutput()->GetBufferedRegion().GetSize()[1]; - int volumeSize[3]; + itk::SizeValueType volumeSize[3]; volumeSize[0] = this->GetInputVolume()->GetBufferedRegion().GetSize()[0]; volumeSize[1] = this->GetInputVolume()->GetBufferedRegion().GetSize()[1]; volumeSize[2] = this->GetInputVolume()->GetBufferedRegion().GetSize()[2]; - int inputDVFSize[3]; + itk::SizeValueType inputDVFSize[3]; inputDVFSize[0] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[0]; inputDVFSize[1] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[1]; inputDVFSize[2] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[2]; @@ -254,7 +254,7 @@ CudaWarpForwardProjectionImageFilter ::GPUGenerateData() source_positions[(iProj - iFirstProj) * 3 + d] = source_position[d]; // Ignore the 4th component } - int projectionOffset = 0; + itk::OffsetValueType projectionOffset = 0; for (unsigned int i = 0; i < nProj; i += SLAB_SIZE) { // If nProj is not a multiple of SLAB_SIZE, the last slab will contain less than SLAB_SIZE projections diff --git a/src/rtkCudaWarpImageFilter.cu b/src/rtkCudaWarpImageFilter.cu index 0038ec72d..ebefea8b4 100644 --- a/src/rtkCudaWarpImageFilter.cu +++ b/src/rtkCudaWarpImageFilter.cu @@ -58,9 +58,9 @@ kernel_3Dgrid(float * dev_vol_out, cudaTextureObject_t tex_zdvf, cudaTextureObject_t tex_input_vol) { - unsigned int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - unsigned int j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; - unsigned int k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; + itk::SizeValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::SizeValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; + itk::SizeValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; if (i >= vol_dim.x || j >= vol_dim.y || k >= vol_dim.z) { @@ -68,7 +68,7 @@ kernel_3Dgrid(float * dev_vol_out, } // Index row major into the volume - long int vol_idx = i + (j + k * vol_dim.y) * (vol_dim.x); + itk::SizeValueType vol_idx = i + (j + k * vol_dim.y) * vol_dim.x; // Matrix multiply to get the index in the DVF texture of the current point in the output volume float3 IndexInDVF = matrix_multiply(make_float3(i, j, k), c_IndexOutputToIndexDVFMatrix); @@ -102,16 +102,16 @@ kernel_3Dgrid(float * dev_vol_out, /////////////////////////////////////////////////////////////////////////// // FUNCTION: CUDA_warp ///////////////////////////// void -CUDA_warp(int input_vol_dim[3], - int input_dvf_dim[3], - int output_vol_dim[3], - float IndexOutputToPPOutputMatrix[12], - float IndexOutputToIndexDVFMatrix[12], - float PPInputToIndexInputMatrix[12], - float * dev_input_vol, - float * dev_output_vol, - float * dev_DVF, - bool isLinear) +CUDA_warp(itk::SizeValueType input_vol_dim[3], + itk::SizeValueType input_dvf_dim[3], + itk::SizeValueType output_vol_dim[3], + float IndexOutputToPPOutputMatrix[12], + float IndexOutputToIndexDVFMatrix[12], + float PPInputToIndexInputMatrix[12], + float * dev_input_vol, + float * dev_output_vol, + float * dev_DVF, + bool isLinear) { // Prepare DVF textures std::vector DVFComponentArrays; diff --git a/src/rtkCudaWarpImageFilter.cxx b/src/rtkCudaWarpImageFilter.cxx index d68137baa..649bcc29d 100644 --- a/src/rtkCudaWarpImageFilter.cxx +++ b/src/rtkCudaWarpImageFilter.cxx @@ -49,17 +49,17 @@ CudaWarpImageFilter ::GPUGenerateData() } // Cuda convenient format for dimensions - int inputVolumeSize[3]; + itk::SizeValueType inputVolumeSize[3]; inputVolumeSize[0] = this->GetInput(0)->GetBufferedRegion().GetSize()[0]; inputVolumeSize[1] = this->GetInput(0)->GetBufferedRegion().GetSize()[1]; inputVolumeSize[2] = this->GetInput(0)->GetBufferedRegion().GetSize()[2]; - int inputDVFSize[3]; + itk::SizeValueType inputDVFSize[3]; inputDVFSize[0] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[0]; inputDVFSize[1] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[1]; inputDVFSize[2] = this->GetDisplacementField()->GetBufferedRegion().GetSize()[2]; - int outputVolumeSize[3]; + itk::SizeValueType outputVolumeSize[3]; outputVolumeSize[0] = this->GetOutput()->GetBufferedRegion().GetSize()[0]; outputVolumeSize[1] = this->GetOutput()->GetBufferedRegion().GetSize()[1]; outputVolumeSize[2] = this->GetOutput()->GetBufferedRegion().GetSize()[2]; diff --git a/src/rtkCudaWeidingerForwardModelImageFilter.cu b/src/rtkCudaWeidingerForwardModelImageFilter.cu index 6fee318d2..a76f18d51 100644 --- a/src/rtkCudaWeidingerForwardModelImageFilter.cu +++ b/src/rtkCudaWeidingerForwardModelImageFilter.cu @@ -64,9 +64,9 @@ kernel_forward_model(float * pMatProj, unsigned int nProjSpectrum, int nIdxProj) { - unsigned int i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; - unsigned int j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; - unsigned int k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; + itk::SizeValueType i = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; + itk::SizeValueType j = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; + itk::SizeValueType k = __umul24(blockIdx.z, blockDim.z) + threadIdx.z; if (i >= c_projSize.x || j >= c_projSize.y || k >= c_projSize.z) { @@ -74,9 +74,9 @@ kernel_forward_model(float * pMatProj, } // Index row major in the projection - long int first_proj_idx = - i + (j + (nIdxProj + k) % nProjSpectrum * c_projSize.y) * c_projSize.x; // To determine the efficient spectrum - long int proj_idx = i + (j + k * c_projSize.y) * (c_projSize.x); // For all the rest + itk::OffsetValueType first_proj_idx = + i + (j + ((nIdxProj + k) % nProjSpectrum) * c_projSize.y) * c_projSize.x; // To determine the efficient spectrum + itk::SizeValueType proj_idx = i + (j + k * c_projSize.y) * c_projSize.x; // For all the rest // Compute the efficient spectrum at the current pixel float efficientSpectrum[VBins * VEnergies];