Issue with the RTK output generated using a head phantom


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();

        // 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 {
        } 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();

    // Write the geometry to disk
    rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
    xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();

    // 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);

    ImageType::SpacingType spacing;
    spacing.Fill(0.5); // Adjusted spacing between slices

    ImageType::SizeType sizeOutput;
    sizeOutput[0] = 512; // Image Columns
    sizeOutput[1] = 512; // Image Rows
    sizeOutput[2] = 512; // Slices

    // Set the constant value for the image source

    // Perform casting from unsigned 16-bit integer to float
    using CastFilterType = itk::CastImageFilter<InputImageType, OutputImageType>;
    CastFilterType::Pointer castconstantImageSource = CastFilterType::New();

    // Perform the reconstruction
    std::cout << "Reconstructing..." << std::endl;

    // 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());

    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();

       // Field-of-view masking
       using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>;
       FOVFilterType::Pointer fieldofview = FOVFilterType::New();
       fieldofview->SetInput(0, casterUnsignedShort->GetOutput());

    // Writer
    std::cout << "Writing output image..." << std::endl;
    using WriterType = itk::ImageFileWriter<ImageType>;
    WriterType::Pointer writer = WriterType::New();

    std::cout << "Reconstruction done." << std::endl;

    return EXIT_SUCCESS;

2.mhd (323 Bytes)
2.xml (124.0 KB)
combined_projections_meta.mhd (278 Bytes)

You only provided the mhd (header) file. Maybe you can share some screenshot of the reconstructed image?

I am attaching few images of the reconstructed slices and also 3D image generated of the same.

I have never seen such artifacts. One thing is that you don’t have a full 360° acquisition, only ~210° so you need to use a ParkerShortScanImageFilter as in rtkfdk which implements [Parker, Med Phys, 1982]. Beside, combined_projections_meta.mhd has no Offset so the origin of your projections is [0,0,0] and you need to offset them in the geometry.

Sorry, but I did not understand what you meant by “offset the origin in the geometry”.

Please have a look at the geometry doc. Currently, the origin of your projection is not set so the point (0,0) is located in the first pixel. You need to either change it (typically, to -0.5*(nbpixel-1)*spacing) or you need to set it with ProjectionOffsetX / ProjectionOffsetY.

But even though I have not set the values for ProjectionOffsetX / ProjectionOffsetY as 0, by default its value automatically becomes zero right.

Right, which is a problem if the projection origin is also at 0 (the projection origin is recorded as “Origin” in the mhd file, sorry for the confusion). Then, the projection corner is in front of the center of rotation.

Sorry, I did not understand “the projection corner is in front of the center of rotation” this statement.

Did you try to read the doc? Without any offset set, the center of rotation is (0,0,0) and it projects on the point with coordinate (0,0) in the projections. If the projection origin (in the sense of m_Origin in itk::Image) is set at (0,0), the first pixel of the projection has coordinates (0,0). This is why m_Origin of the projections is generally set such that point (0,0) is in the center of the image., see here.