Issue with GPU Memory Not Releasing After Reconstruction in RTK FDK

[SS99aaNN]

Dear Developer,

For certain reasons, my program requires performing reconstruction in a loop. However, it seems that after each reconstruction, the GPU memory is not being released. I am wondering if this is a bug or if it could be related to the version of CUDA I am using.

I am currently using ITK-5.4 with the RTK module enabled, and my CUDA version is 12.4.
The image below shows the memory usage when I loop through FDK reconstruction with dimensions 896x896x720, where all values in the projection are set to 1.

Thank you for your assistance with this matter.

Best regards,
sean

[SS99aaNN]

here i prepared a simplified example to reproduce my issue.

//#include "RTKRecon.h"
// RTK includes
#include "rtkConstantImageSource.h"
#include "rtkFieldOfViewImageFilter.h"
#include "rtkProgressCommands.h"
#include <rtkImportImageFilter.h>
// CUDA
#include "rtkCudaFDKConeBeamReconstructionFilter.h"



float* rtkCudaFDK_simplified(float* projection_input)
{
	// 开始重建
	// Read the RTK geometry object
	// Write the geometry to dis
	using OutputPixelType = float;
	constexpr unsigned int Dimension = 3;
	using CPUOutputImageType = itk::Image<OutputPixelType, Dimension>;
	using OutputImageType = itk::CudaImage<OutputPixelType, Dimension>;

	using ConstantImageSourceType = rtk::ConstantImageSource<OutputImageType>;
	using GeometryType = rtk::ThreeDCircularProjectionGeometry;
	using ImportFilterType = rtk::ImportImageFilter<OutputImageType>;

	using FDKCUDAType = rtk::CudaFDKConeBeamReconstructionFilter;
	using PercentageProgressCommandType = rtk::PercentageProgressCommand<FDKCUDAType>;
	using CopyImageType = itk::Image<OutputPixelType, Dimension>;
	using CastFilterType = itk::CastImageFilter<OutputImageType, CopyImageType>;

	//setting up the geometry
	GeometryType::Pointer geometry = GeometryType::New();
	for (int noProj = 0; noProj < 720; noProj++)
		geometry->AddProjection(
			240.219,
			359.98,
			0 + noProj * 360. / 720,
			0.0628873,
			-0.0217019,
			0.,
			-0.0605461,
			0.,
			0.);

	// read the input projection for reconstruction
	ConstantImageSourceType::PointType   input_projection_origin;
	ConstantImageSourceType::SpacingType input_projection_spacing;
	ConstantImageSourceType::SizeType    input_projection_sizeOutput;

	input_projection_origin[0] = -89.5;
	input_projection_origin[1] = -89.5;
	input_projection_origin[2] = -71.9;

	input_projection_spacing[0] = 0.2;
	input_projection_spacing[1] = 0.2;
	input_projection_spacing[2] = 0.2;

	input_projection_sizeOutput[0] = 896;
	input_projection_sizeOutput[1] = 896;
	input_projection_sizeOutput[2] = 720;

	long long proj_img_len = (long long)896 * 896 * 720;
	ImportFilterType::Pointer projection_import_filter = ImportFilterType::New();

	OutputImageType::IndexType start;
	start.Fill(0);
	OutputImageType::RegionType region;
	region.SetSize(input_projection_sizeOutput);
	region.SetIndex(start);
	projection_import_filter->SetRegion(region);

	projection_import_filter->SetOrigin(input_projection_origin);
	projection_import_filter->SetSpacing(input_projection_spacing);
	projection_import_filter->SetImportPointer(projection_input, proj_img_len, false);


	// create a volume of reconstruction image
	ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New();
	ConstantImageSourceType::PointType   output_volume_origin;
	ConstantImageSourceType::SpacingType output_volume_spacing;
	ConstantImageSourceType::SizeType    output_volume_sizeOutput;

	output_volume_origin[0] = -51.15;
	output_volume_origin[1] = -51.15;
	output_volume_origin[2] = -51.15;

	output_volume_spacing[0] = 0.1;
	output_volume_spacing[1] = 0.1;
	output_volume_spacing[2] = 0.1;

	output_volume_sizeOutput[0] = 1024;
	output_volume_sizeOutput[1] = 1024;
	output_volume_sizeOutput[2] = 1024;

	constantImageSource->SetOrigin(output_volume_origin);
	constantImageSource->SetSpacing(output_volume_spacing);
	constantImageSource->SetSize(output_volume_sizeOutput);
	constantImageSource->SetConstant(0.);

	// FDK reconstruction
	FDKCUDAType::Pointer feldkamp = FDKCUDAType::New();

	feldkamp->SetInput(0, constantImageSource->GetOutput());
	feldkamp->SetInput(1, projection_import_filter->GetOutput());
	feldkamp->SetGeometry(geometry);
	feldkamp->GetRampFilter()->SetTruncationCorrection(0); // default = 0
	feldkamp->GetRampFilter()->SetHannCutFrequency(0); // default = 0
	feldkamp->GetRampFilter()->SetHannCutFrequencyY(0); // default = 0

	CastFilterType::Pointer castFilter = CastFilterType::New();
	castFilter->SetInput(feldkamp->GetOutput());
	castFilter->Update();


	long long recon_img_len = (long long)1024 * 1024 * 1024;
	float* output_test = (new float[recon_img_len]);
	std::copy(castFilter->GetOutput()->GetBufferPointer(), castFilter->GetOutput()->GetBufferPointer() + castFilter->GetOutput()->GetBufferedRegion().GetNumberOfPixels(), output_test);
	return(output_test);
}


void main()
{
	long long proj_img_len = (long long)896 * 896 * 720;
	long long recon_img_len = (long long)1024 * 1024 * 1024;
	float *proj = new float[proj_img_len];
	for (long long i = 0; i < proj_img_len; i++) { proj[i] = 1.f; }
	float *test = new float[proj_img_len];

	for (int i = 0; i < 5; i++) {
		size_t free, total;
		cudaMemGetInfo(&free, &total);
		std::cout << "Free memory: " << free / (1024 * 1024) << " MB" << std::endl;
		std::cout << "Total memory: " << total / (1024 * 1024) << " MB" << std::endl;
		std::cout << "recon: " << i << std::endl;
		test = rtkCudaFDK_simplified(proj);
		std::cout << std::endl;
	}
}

[SimonRit]

Thanks for the report. I think I found the issue, see #614.

[SimonRit]

Fixed by #614.