Hi,
I am trying to reconstruct a CBCT raw data using RTK. I was able to execute the code and generate output, but unfortunately I am facing issues with the output generated. I don’t know if the issue is due to the geometry or any other parameters. I will write my code below and will also attach my xml file and output here. I humbly for request any advices and solutions. Thanks in advance.
#include <iostream>
#include <fstream>
#include <itkCudaImage.h>
#include <rtkConstantImageSource.h>
#include <rtkThreeDCircularProjectionGeometryXMLFileWriter.h>
#include <rtkCudaFDKConeBeamReconstructionFilter.h>
#include <rtkFieldOfViewImageFilter.h>
#include <itkImageFileReader.h>
#include <itkTimeProbe.h>
#include <rtkThreeDCircularProjectionGeometryXMLFileReader.h>
#include <itkCastImageFilter.h>
#include <itkMetaImageIO.h>
#include <cmath>
#include <cstdint>
#define M_PI 3.14159265358979323846
//// Function to convert degrees to radians
//double degreesToRadians(double degrees) {
// return degrees * M_PI / 180.0;
//}
int main(int argc, char **argv)
{
// if (argc < 4)
// {
// std::cout << "Usage: CbctReconstruction <inputmhdfile> <outputimage> <outputgeometry>" << std::endl;
// return EXIT_FAILURE;
// }
// Create a time probe
itk::TimeProbe clock;
// Define the image type
using ImageType = itk::CudaImage<uint16_t, 3>;
// Define the RTK geometry object
using GeometryType = rtk::ThreeDCircularProjectionGeometry;
GeometryType::Pointer geometry = GeometryType::New();
unsigned int numberOfProjections = 351; // Adjusted number of projections
// System geometry parameters
const double DSD = 1800.0; // Distance Source Detector (mm)
const double DSO = 921.0; // Distance Source Origin (mm)
// // const double detectorPixelPitch = 0.139; // Detector pixel pitch (mm)
// double firstAngle = 0;
// double angularArc = 351;
// for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)
// {
// double angle = firstAngle + noProj * angularArc / numberOfProjections;
// // double angleRadians = degreesToRadians(angle); // Convert angle to radians
// geometry->AddProjection(DSO, DSD, angle);
// geometry->AddProjection(DSO, DSD, angle);
// }
// Assuming the file contains angles separated by newlines
std::ifstream angleFile("D:/Navami/TIGRE_OUT_NEW/Data/CBCT/angles/angle.txt");
if (!angleFile.is_open()) {
std::cerr << "Error opening angle file." << std::endl;
return EXIT_FAILURE;
}
////////-----TO READ ANGLES IN RADIANS--------////////
// Create the geometry
for (unsigned int noProj = 1; noProj <= numberOfProjections; ++noProj)
{
double angle;
angleFile >> angle;
// double angleRadians = degreesToRadians(angle); // Convert angle to radians
geometry->AddProjection(DSO, DSD, angle);
// std::cout << "Angle " << noProj << ": " << angle << std::endl;
}
// Read the input image from MetaImage file
using ReaderType = itk::ImageFileReader<ImageType>;
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName("D:/Navami/TIGRE_OUT_NEW/CBCT_Projections/combined_projections_meta.mhd");
// Get the MetaDataDictionary from the ImageIO
itk::MetaDataDictionary &metadata = reader->GetMetaDataDictionary();
// Set detector row, column, and pixel pitch
itk::EncapsulateMetaData<std::string>(metadata, "detector_row", "1536");
itk::EncapsulateMetaData<std::string>(metadata, "detector_column", "1536");
itk::EncapsulateMetaData<std::string>(metadata, "detector_pixel_pitch", "0.139");
try {
reader->Update();
} catch (itk::ExceptionObject &excp) {
std::cerr << "Exception thrown while reading the image file " << argv[1] << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
// Define the image type
using InputImageType = itk::CudaImage<uint16_t, 3>;
using OutputImageType = itk::CudaImage<float, 3>;
// Cast the image type from uint16_t to float
using CastFilterType = itk::CastImageFilter<InputImageType, OutputImageType>;
CastFilterType::Pointer caster = CastFilterType::New();
caster->SetInput(reader->GetOutput());
caster->Update();
// Write the geometry to disk
rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
xmlWriter->SetFilename("D:/Navami/TIGRE_OUT_NEW/CBCT_Out_Geo/CBCT_GeoOut_13.03.24/2.xml");
xmlWriter->SetObject(geometry);
xmlWriter->WriteFile();
// Create a constant image source for reconstruction
using ConstantImageSourceType = rtk::ConstantImageSource<ImageType>;
ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New();
// Set image origin, spacing, and size based on provided parameters
ImageType::PointType origin;
origin.Fill(0.5 * (512 - 1) * -0.5);
constantImageSource->SetOrigin(origin);
ImageType::SpacingType spacing;
spacing.Fill(0.5); // Adjusted spacing between slices
constantImageSource->SetSpacing(spacing);
ImageType::SizeType sizeOutput;
sizeOutput[0] = 512; // Image Columns
sizeOutput[1] = 512; // Image Rows
sizeOutput[2] = 512; // Slices
constantImageSource->SetSize(sizeOutput);
// Set the constant value for the image source
constantImageSource->SetConstant(0.0);
constantImageSource->Update();
// Perform casting from unsigned 16-bit integer to float
using CastFilterType = itk::CastImageFilter<InputImageType, OutputImageType>;
CastFilterType::Pointer castconstantImageSource = CastFilterType::New();
castconstantImageSource->SetInput(constantImageSource->GetOutput());
castconstantImageSource->Update();
// Perform the reconstruction
std::cout << "Reconstructing..." << std::endl;
clock.Start();
// FDK reconstruction
// std::cout << "Reconstructing..." << std::endl;
using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter;
FDKGPUType::Pointer feldkamp = FDKGPUType::New();
feldkamp->SetInput(0, castconstantImageSource->GetOutput());
feldkamp->SetInput(1, caster->GetOutput());
feldkamp->SetGeometry(geometry);
feldkamp->GetRampFilter()->SetTruncationCorrection(0.);
feldkamp->GetRampFilter()->SetHannCutFrequency(0.0);
clock.Stop();
std::cout << "Time taken for reconstruction: " << clock.GetTotal() << " seconds." << std::endl;
using OutputImageTypeUnsignedShort = itk::CudaImage<uint16_t, 3>;
// Cast the output of the FDK reconstruction filter to unsigned short
using CastFilterTypeUnsignedShort = itk::CastImageFilter<OutputImageType, OutputImageTypeUnsignedShort>;
CastFilterTypeUnsignedShort::Pointer casterUnsignedShort = CastFilterTypeUnsignedShort::New();
casterUnsignedShort->SetInput(feldkamp->GetOutput());
casterUnsignedShort->Update();
// Field-of-view masking
using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>;
FOVFilterType::Pointer fieldofview = FOVFilterType::New();
fieldofview->SetInput(0, casterUnsignedShort->GetOutput());
fieldofview->SetProjectionsStack(reader->GetOutput());
fieldofview->SetGeometry(geometry);
fieldofview->Update();
// Writer
std::cout << "Writing output image..." << std::endl;
using WriterType = itk::ImageFileWriter<ImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName("D:/Navami/TIGRE_OUT_NEW/CBCT_Out_Geo/CBCT_GeoOut_13.03.24/4.mhd");
writer->SetInput(casterUnsignedShort->GetOutput());
writer->Update();
std::cout << "Reconstruction done." << std::endl;
return EXIT_SUCCESS;
}
2.mhd (323 Bytes)
2.xml (124.0 KB)
combined_projections_meta.mhd (278 Bytes)
