diff --git a/include/rtkCudaAverageOutOfROIImageFilter.hcu b/include/rtkCudaAverageOutOfROIImageFilter.hcu index d4d770044..497f007ca 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(SizeValueType4 size, float * input, float * output, float * roi); #endif diff --git a/include/rtkCudaConstantVolumeSeriesSource.hcu b/include/rtkCudaConstantVolumeSeriesSource.hcu index 1206a506b..4951b2bf2 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(SizeValueType4 size, float * dev_out, float constantValue); #endif diff --git a/include/rtkCudaConstantVolumeSource.hcu b/include/rtkCudaConstantVolumeSource.hcu index 327a20a2f..bc04ac3e8 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(SizeValueType3 size, float * dev_out, float constantValue); #endif diff --git a/include/rtkCudaCyclicDeformationImageFilter.hcu b/include/rtkCudaCyclicDeformationImageFilter.hcu index c5ac08700..72c0a708f 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(SizeValueType4 inputSize, + float * input, + float * output, + unsigned int frameInf, + unsigned int frameSup, + double weightInf, + double weightSup); #endif diff --git a/include/rtkCudaForwardProjectionImageFilter.hxx b/include/rtkCudaForwardProjectionImageFilter.hxx index e01a50d20..f24cc025d 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; @@ -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/rtkCudaInterpolateImageFilter.hcu b/include/rtkCudaInterpolateImageFilter.hcu index 5c0ba262d..edd8cb453 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 SizeValueType4 & inputSize, + float * input, + float * output, + int projectionNumber, + float ** weights); #endif diff --git a/include/rtkCudaSplatImageFilter.hcu b/include/rtkCudaSplatImageFilter.hcu index dd85716ea..740c788ca 100644 --- a/include/rtkCudaSplatImageFilter.hcu +++ b/include/rtkCudaSplatImageFilter.hcu @@ -19,9 +19,9 @@ #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 SizeValueType4 & outputSize, float * input, float * output, int projectionNumber, float ** weights); #endif diff --git a/include/rtkCudaUtilities.hcu b/include/rtkCudaUtilities.hcu index 4dc5ec4d8..90d18ea92 100644 --- a/include/rtkCudaUtilities.hcu +++ b/include/rtkCudaUtilities.hcu @@ -23,10 +23,14 @@ #include #include #define ITK_STATIC +#include #include #undef ITK_STATIC #include +using SizeValueType3 = ulonglong3; +using SizeValueType4 = ulonglong4; + #define CUDA_CHECK_ERROR \ { \ cudaError_t err = cudaGetLastError(); \ @@ -112,9 +116,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) diff --git a/src/rtkCudaAverageOutOfROIImageFilter.cu b/src/rtkCudaAverageOutOfROIImageFilter.cu index ee2706541..c3fd77a2b 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::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; // 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::SizeValueType id = (k * c_Size.y + j) * c_Size.x + i; + itk::SizeValueType strided_id = id; // strided_id will run along the 4th dimension // Compute the average along last dimension float avg = 0; @@ -77,21 +77,20 @@ 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(SizeValueType4 size, 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.x * size.y * size.z; // Thread Block Dimensions dim3 dimBlock = dim3(8, 8, 8); - int blocksInX = iDivUp(size[0], dimBlock.x); - int blocksInY = iDivUp(size[1], dimBlock.y); - int blocksInZ = iDivUp(size[2], dimBlock.z); + int blocksInX = iDivUp(size.x, dimBlock.x); + int blocksInY = iDivUp(size.y, dimBlock.y); + int blocksInZ = iDivUp(size.z, dimBlock.z); dim3 dimGrid = dim3(blocksInX, blocksInY, blocksInZ); diff --git a/src/rtkCudaAverageOutOfROIImageFilter.cxx b/src/rtkCudaAverageOutOfROIImageFilter.cxx index 7eca6f93c..13e2d2d08 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]; - } + SizeValueType4 size; + size.x = this->GetOutput()->GetBufferedRegion().GetSize()[0]; + size.y = this->GetOutput()->GetBufferedRegion().GetSize()[1]; + size.z = this->GetOutput()->GetBufferedRegion().GetSize()[2]; + size.w = 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/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..d41fbd41e 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::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 * c_Size.w) 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] = 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(SizeValueType4 size, 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 @@ -74,7 +73,7 @@ CUDA_generate_constant_volume_series(int size[4], float * dev_out, float constan // run a kernel to replace the zeros with constantValue. // Reset output volume - size_t memorySizeOutput = size[0] * size[1] * size[2] * size[3] * sizeof(float); + size_t memorySizeOutput = size.x * size.y * size.z * size.w * sizeof(float); cudaMemset((void *)dev_out, 0, memorySizeOutput); if (!(constantValue == 0)) @@ -82,9 +81,9 @@ CUDA_generate_constant_volume_series(int size[4], float * dev_out, float constan // Thread Block Dimensions dim3 dimBlock = dim3(4, 4, 16); - int blocksInX = iDivUp(size[0], dimBlock.x); - int blocksInY = iDivUp(size[1], dimBlock.y); - int blocksInZ = iDivUp(size[2] * size[3], dimBlock.z); + int blocksInX = iDivUp(size.x, dimBlock.x); + int blocksInY = iDivUp(size.y, dimBlock.y); + int blocksInZ = iDivUp(size.z * size.w, dimBlock.z); dim3 dimGrid = dim3(blocksInX, blocksInY, blocksInZ); diff --git a/src/rtkCudaConstantVolumeSeriesSource.cxx b/src/rtkCudaConstantVolumeSeriesSource.cxx index 32f51196c..e005c7988 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]; - } + SizeValueType4 outputSize; + outputSize.x = this->GetOutput()->GetRequestedRegion().GetSize()[0]; + outputSize.y = this->GetOutput()->GetRequestedRegion().GetSize()[1]; + outputSize.z = this->GetOutput()->GetRequestedRegion().GetSize()[2]; + outputSize.w = 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..67da41de2 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::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] = 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(SizeValueType3 size, 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 @@ -74,7 +73,7 @@ CUDA_generate_constant_volume(int size[3], float * dev_out, float constantValue) // run a kernel to replace the zeros with constantValue. // Reset output volume - size_t memorySizeOutput = size[0] * size[1] * size[2] * sizeof(float); + size_t memorySizeOutput = size.x * size.y * size.z * sizeof(float); cudaMemset((void *)dev_out, 0, memorySizeOutput); if (!(constantValue == 0)) @@ -82,9 +81,9 @@ CUDA_generate_constant_volume(int size[3], float * dev_out, float constantValue) // Thread Block Dimensions dim3 dimBlock = dim3(16, 4, 4); - int blocksInX = iDivUp(size[0], dimBlock.x); - int blocksInY = iDivUp(size[1], dimBlock.y); - int blocksInZ = iDivUp(size[2], dimBlock.z); + int blocksInX = iDivUp(size.x, dimBlock.x); + int blocksInY = iDivUp(size.y, dimBlock.y); + int blocksInZ = iDivUp(size.z, dimBlock.z); dim3 dimGrid = dim3(blocksInX, blocksInY, blocksInZ); diff --git a/src/rtkCudaConstantVolumeSource.cxx b/src/rtkCudaConstantVolumeSource.cxx index 4717d0d07..96604c9a7 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]; - } + SizeValueType3 outputSize; + outputSize.x = this->GetOutput()->GetRequestedRegion().GetSize()[0]; + outputSize.y = this->GetOutput()->GetRequestedRegion().GetSize()[1]; + outputSize.z = 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..1c74edaf8 100644 --- a/src/rtkCudaCropImageFilter.cu +++ b/src/rtkCudaCropImageFilter.cu @@ -33,17 +33,17 @@ 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::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 >= 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::SizeValueType out_idx = i + j * cropDim.x + k * cropDim.y * cropDim.x; + itk::SizeValueType 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..bcbc8e741 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(SizeValueType4 inputSize, + float * input, + float * output, + unsigned int frameInf, + unsigned int frameSup, + double weightInf, + double weightSup) { cublasHandle_t handle; cublasCreate(&handle); @@ -51,7 +49,7 @@ CUDA_linear_interpolate_along_fourth_dimension(unsigned int inputSize[4], float wInf = (float)weightInf; float wSup = (float)weightSup; - size_t numel = inputSize[0] * inputSize[1] * inputSize[2] * 3; + size_t numel = inputSize.x * inputSize.y * inputSize.z * 3; cudaMemset((void *)output, 0, numel * sizeof(float)); diff --git a/src/rtkCudaCyclicDeformationImageFilter.cxx b/src/rtkCudaCyclicDeformationImageFilter.cxx index 5234d60e4..ef5e71b59 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++) + SizeValueType4 inputSize; + inputSize.x = this->GetInput()->GetBufferedRegion().GetSize()[0]; + inputSize.y = this->GetInput()->GetBufferedRegion().GetSize()[1]; + inputSize.z = this->GetInput()->GetBufferedRegion().GetSize()[2]; + inputSize.w = this->GetInput()->GetBufferedRegion().GetSize()[3]; + if ((this->GetOutput()->GetRequestedRegion().GetSize()[0] != inputSize.x) || + (this->GetOutput()->GetRequestedRegion().GetSize()[1] != inputSize.y) || + (this->GetOutput()->GetRequestedRegion().GetSize()[2] != inputSize.z)) { - 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..e4aa9cb65 100644 --- a/src/rtkCudaDisplacedDetectorImageFilter.cu +++ b/src/rtkCudaDisplacedDetectorImageFilter.cu @@ -60,11 +60,11 @@ kernel_displaced_weight(int3 proj_idx_in, ) { // compute thread index - int3 tIdx; + SizeValueType3 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::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) @@ -72,9 +72,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::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; // 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/rtkCudaFDKWeightProjectionFilter.cu b/src/rtkCudaFDKWeightProjectionFilter.cu index 5f274980f..bcbcb702f 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::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; - 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..bdeee4cbe 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::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 >= fftDimension.x || j >= fftDimension.y || k >= fftDimension.z) return; - long int proj_idx = i + (j + k * fftDimension.y) * fftDimension.x; + itk::SizeValueType 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::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 >= 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::SizeValueType kernel_idx = i + j * fftDimension.x; + itk::SizeValueType 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::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 >= 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..a66ec5e0d 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::OffsetValueType i = blockIdx.x * blockDim.x + threadIdx.x; + itk::OffsetValueType j = blockIdx.y * blockDim.y + threadIdx.y; + itk::OffsetValueType 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::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] = 0; diff --git a/src/rtkCudaForwardProjectionImageFilter.cu b/src/rtkCudaForwardProjectionImageFilter.cu index 6494bc461..c729d5d19 100644 --- a/src/rtkCudaForwardProjectionImageFilter.cu +++ b/src/rtkCudaForwardProjectionImageFilter.cu @@ -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::OffsetValueType 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) diff --git a/src/rtkCudaForwardWarpImageFilter.cu b/src/rtkCudaForwardWarpImageFilter.cu index 17e7030f8..59466809e 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) && @@ -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/rtkCudaInterpolateImageFilter.cu b/src/rtkCudaInterpolateImageFilter.cu index 9760f0119..df83e1728 100644 --- a/src/rtkCudaInterpolateImageFilter.cu +++ b/src/rtkCudaInterpolateImageFilter.cu @@ -28,7 +28,11 @@ #include void -CUDA_interpolation(const int4 & inputSize, float * input, float * output, int projectionNumber, float ** weights) +CUDA_interpolation(const SizeValueType4 & inputSize, + float * input, + float * output, + int projectionNumber, + float ** weights) { cublasHandle_t handle; cublasCreate(&handle); @@ -40,7 +44,7 @@ CUDA_interpolation(const int4 & inputSize, float * input, float * output, int pr // Reset output volume cudaMemset((void *)output, 0, memorySizeOutput); - for (int phase = 0; phase < inputSize.w; phase++) + for (itk::SizeValueType phase = 0; phase < inputSize.w; phase++) { float weight = weights[phase][projectionNumber]; if (weight != 0) diff --git a/src/rtkCudaInterpolateImageFilter.cxx b/src/rtkCudaInterpolateImageFilter.cxx index f9a69e42a..75b5d86d8 100644 --- a/src/rtkCudaInterpolateImageFilter.cxx +++ b/src/rtkCudaInterpolateImageFilter.cxx @@ -29,7 +29,7 @@ CudaInterpolateImageFilter ::CudaInterpolateImageFilter() = default; void CudaInterpolateImageFilter ::GPUGenerateData() { - int4 inputSize; + SizeValueType4 inputSize; inputSize.x = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[0]; inputSize.y = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[1]; inputSize.z = this->GetInputVolumeSeries()->GetBufferedRegion().GetSize()[2]; 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..2a7a5a490 100644 --- a/src/rtkCudaSplatImageFilter.cu +++ b/src/rtkCudaSplatImageFilter.cu @@ -28,14 +28,14 @@ #include void -CUDA_splat(const int4 & outputSize, float * input, float * output, int projectionNumber, float ** weights) +CUDA_splat(const SizeValueType4 & outputSize, float * input, float * output, int projectionNumber, float ** weights) { cublasHandle_t handle; cublasCreate(&handle); size_t numel = outputSize.x * outputSize.y * outputSize.z; - for (int phase = 0; phase < outputSize.w; phase++) + for (itk::SizeValueType phase = 0; phase < outputSize.w; phase++) { float weight = weights[phase][projectionNumber]; if (weight != 0) diff --git a/src/rtkCudaSplatImageFilter.cxx b/src/rtkCudaSplatImageFilter.cxx index 14d198df8..4b821258f 100644 --- a/src/rtkCudaSplatImageFilter.cxx +++ b/src/rtkCudaSplatImageFilter.cxx @@ -29,7 +29,7 @@ CudaSplatImageFilter ::CudaSplatImageFilter() = default; void CudaSplatImageFilter ::GPUGenerateData() { - int4 outputSize; + SizeValueType4 outputSize; outputSize.x = this->GetOutput()->GetLargestPossibleRegion().GetSize()[0]; outputSize.y = this->GetOutput()->GetLargestPossibleRegion().GetSize()[1]; outputSize.z = this->GetOutput()->GetLargestPossibleRegion().GetSize()[2]; 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..d479797a3 100644 --- a/src/rtkCudaUtilities.cu +++ b/src/rtkCudaUtilities.cu @@ -166,7 +166,7 @@ prepareVectorTextureObject(int size[3], // Allocate an intermediate memory space to extract the components of the input volume float * singleComponent; - size_t numel = size[0] * size[1] * size[2]; + size_t numel = static_cast(size[0]) * size[1] * size[2]; if (nComponents > 1) { cudaMalloc(&singleComponent, numel * sizeof(float)); diff --git a/src/rtkCudaWarpBackProjectionImageFilter.cu b/src/rtkCudaWarpBackProjectionImageFilter.cu index 778cd7c46..2ca616a54 100644 --- a/src/rtkCudaWarpBackProjectionImageFilter.cu +++ b/src/rtkCudaWarpBackProjectionImageFilter.cu @@ -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; diff --git a/src/rtkCudaWarpForwardProjectionImageFilter.cu b/src/rtkCudaWarpForwardProjectionImageFilter.cu index 4b73e135f..f6965f09a 100644 --- a/src/rtkCudaWarpForwardProjectionImageFilter.cu +++ b/src/rtkCudaWarpForwardProjectionImageFilter.cu @@ -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; } } } diff --git a/src/rtkCudaWarpForwardProjectionImageFilter.cxx b/src/rtkCudaWarpForwardProjectionImageFilter.cxx index 1756ec964..aca1284a8 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; @@ -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..965120d3e 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); 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];