Skip to content

BUG: Avoid 32-bit index/offset overflow in CUDA filters for large images #989

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

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

BUG: Avoid 32-bit index/offset overflow in CUDA filters for large images #989

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

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions include/rtkCudaAverageOutOfROIImageFilter.hcu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,9 @@
#ifndef rtkCudaAverageOutOfROIImageFilter_hcu
#define rtkCudaAverageOutOfROIImageFilter_hcu

#include <vector_types.h>
#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
4 changes: 2 additions & 2 deletions include/rtkCudaConstantVolumeSeriesSource.hcu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,9 @@
#ifndef rtkCudaConstantVolumeSeriesSource_hcu
#define rtkCudaConstantVolumeSeriesSource_hcu

#include <vector_types.h>
#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
4 changes: 2 additions & 2 deletions include/rtkCudaConstantVolumeSource.hcu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,9 @@
#ifndef rtkCudaConstantVolumeSource_hcu
#define rtkCudaConstantVolumeSource_hcu

#include <vector_types.h>
#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
16 changes: 9 additions & 7 deletions include/rtkCudaCyclicDeformationImageFilter.hcu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
14 changes: 7 additions & 7 deletions include/rtkCudaForwardProjectionImageFilter.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,11 +49,11 @@ CudaForwardProjectionImageFilter<TInputImage, TOutputImage>::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<typename TInputImage::PixelType>::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<typename TInputImage::PixelType>::GetLength();

itk::Vector<double, 4> source_position;

Expand Down Expand Up @@ -166,8 +166,8 @@ CudaForwardProjectionImageFilter<TInputImage, TOutputImage>::GPUGenerateData()
source_positions[(iProj - iFirstProj) * 3 + d] = source_position[d]; // Ignore the 4th component
}

int projectionOffset = 0;
const unsigned int vectorLength = itk::PixelTraits<typename TInputImage::PixelType>::Dimension;
itk::OffsetValueType projectionOffset = 0;
const unsigned int vectorLength = itk::PixelTraits<typename TInputImage::PixelType>::Dimension;

for (unsigned int i = 0; i < nProj; i += SLAB_SIZE)
{
Expand Down
8 changes: 6 additions & 2 deletions include/rtkCudaInterpolateImageFilter.hcu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,13 @@
#ifndef rtkCudaInterpolateImageFilter_hcu
#define rtkCudaInterpolateImageFilter_hcu

#include <vector_types.h>
#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
4 changes: 2 additions & 2 deletions include/rtkCudaSplatImageFilter.hcu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,9 @@
#ifndef rtkCudaSplatImageFilter_hcu
#define rtkCudaSplatImageFilter_hcu

#include <vector_types.h>
#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
8 changes: 6 additions & 2 deletions include/rtkCudaUtilities.hcu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,14 @@
#include <string>
#include <vector>
#define ITK_STATIC
#include <itkIntTypes.h>
#include <itkMacro.h>
#undef ITK_STATIC
#include <cuda.h>

using SizeValueType3 = ulonglong3;
using SizeValueType4 = ulonglong4;

#define CUDA_CHECK_ERROR \
{ \
cudaError_t err = cudaGetLastError(); \
Expand Down Expand Up @@ -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<int>((a % b != 0) ? (a / b + 1) : (a / b));
}
inline __host__ __device__ float
dot_vector(float3 u, float3 v)
Expand Down
27 changes: 13 additions & 14 deletions src/rtkCudaAverageOutOfROIImageFilter.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@

// TEXTURES AND CONSTANTS //

__constant__ int4 c_Size;
__constant__ SizeValueType4 c_Size;

//_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
// K E R N E L S -_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
Expand All @@ -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;
Expand Down Expand Up @@ -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);

Expand Down
11 changes: 5 additions & 6 deletions src/rtkCudaAverageOutOfROIImageFilter.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<float *>(this->GetInput()->GetCudaDataManager()->GetGPUBufferPointer());
float * pout = static_cast<float *>(this->GetOutput()->GetCudaDataManager()->GetGPUBufferPointer());
Expand Down
56 changes: 55 additions & 1 deletion src/rtkCudaConjugateGradientImageFilter.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

@SimonRit SimonRit Sep 4, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess we have a good reason to question the compatibility with Cuda 11 now.

@SimonRit SimonRit Sep 11, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's keep a cleanup for later and keep it as is

cublasSaxpy(handle, (int)numberOfElements, &alpha, toBeSubtracted, 1, out, 1);
#else
cublasSaxpy_64(handle, numberOfElements, &alpha, toBeSubtracted, 1, out, 1);
Expand All @@ -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);
Expand All @@ -87,33 +91,58 @@ 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;

// 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);
Expand All @@ -131,33 +160,58 @@ 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;

// 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);
Expand Down
Loading
Loading