Unable to work with Raw Projection Images in ITK-RTK

Hi, I was trying to do a comparative study with TIGRE and RTK for which I have stored Shepp-Logan data from TIGRE in .raw format. With this data and its corresponding parameters, I am trying to reconstruct it in RTK. While debugging Abort error occurs. I am not sure if it is because I have not provided the geometry parameters properly. The code I have written is shown below:
#include <rtkConstantImageSource.h>
#include <rtkThreeDCircularProjectionGeometryXMLFileWriter.h>
#include <rtkCudaFDKConeBeamReconstructionFilter.h>
#include <rtkFieldOfViewImageFilter.h>
#include <itkImageFileReader.h>
#include <itkRawImageIO.h>
#include <itkTimeProbe.h>
#include <rtkThreeDCircularProjectionGeometryXMLFileReader.h>

int main(int argc, char **argv)
{
    if (argc < 3)
    {
        std::cout << "Usage: FirstReconstruction <outputimage> <outputgeometry>" << std::endl;
        return EXIT_FAILURE;
    }

    // Create a time probe
    itk::TimeProbe clock;

    // Defines the image type
    using ImageType = itk::CudaImage<float, 3>;

    // Defines the RTK geometry object
    using GeometryType = rtk::ThreeDCircularProjectionGeometry;
    GeometryType::Pointer geometry = GeometryType::New();
    unsigned int numberOfProjections = 100;
    /*double firstAngle = 0;
    double angularArc = 360;
    unsigned int sid = 600;  // source to isocenter distance
    unsigned int sdd = 1200;*/ // source to detector distance
    for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)
    {
        // double angle = firstAngle + noProj * angularArc / numberOfProjections;
        geometry->AddProjection(1000, 1536, noProj * 360. / numberOfProjections);
    }

    // Write the geometry to disk
    rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
    xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
    xmlWriter->SetFilename(argv[2]);
    xmlWriter->SetObject(geometry);
    xmlWriter->WriteFile();

    // Read the projection data from a .raw file
    using ProjectionImageType = itk::CudaImage<float, 2>;
    using ReaderType = itk::ImageFileReader<ProjectionImageType>;
    ReaderType::Pointer reader = ReaderType::New();
    itk::RawImageIO<float>::Pointer rawImageIO = itk::RawImageIO<float>::New();
    reader->SetImageIO(rawImageIO);
    std::string rawFileName = "D:\\Navami\\TIGREOutput\\SheppLoganTrial_projections.raw";  // Replace with the actual file path
    reader->SetFileName(rawFileName);
    //    try
    //    {
    reader->Update();
    //    }
    //    catch (itk::ExceptionObject &ex)
    //    {
    //        std::cerr << "Error reading the .raw file: " << ex << std::endl;
    //        return EXIT_FAILURE;
    //    }
    ProjectionImageType::Pointer projectionData = reader->GetOutput();

    // Create a 3D volume to hold the projection data
    using VolumeType = itk::CudaImage<float, 3>;
    VolumeType::Pointer volumeData = VolumeType::New();
    VolumeType::RegionType region;
    VolumeType::SizeType size;
    VolumeType::IndexType index;

    size[0] = projectionData->GetLargestPossibleRegion().GetSize()[0];
    size[1] = projectionData->GetLargestPossibleRegion().GetSize()[1];
    size[2] = numberOfProjections;

    index.Fill(0);
    region.SetSize(size);
    region.SetIndex(index);

    volumeData->SetRegions(region);
    volumeData->Allocate();

    // Copy the projection data into the 3D volume
    using VolumeIteratorType = itk::ImageRegionIterator<VolumeType>;
    VolumeIteratorType volumeIt(volumeData, volumeData->GetLargestPossibleRegion());

    while (!volumeIt.IsAtEnd())
    {
        // Use a 3D index instead of a 2D index
        VolumeType::IndexType index3D;
        index3D[0] = volumeIt.GetIndex()[0];
        index3D[1] = volumeIt.GetIndex()[1];
        index3D[2] = volumeIt.GetIndex()[2];

        // Use the raw pointer to access the pixel value
        volumeIt.Set(projectionData->GetBufferPointer()[volumeData->ComputeOffset(index3D)]);
        ++volumeIt;
    }

    // Create a stack of empty projection images
    using ConstantImageSourceType = rtk::ConstantImageSource<ImageType>;
    ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New();
    ConstantImageSourceType::PointType origin;
    ConstantImageSourceType::SpacingType spacing;
    ConstantImageSourceType::SizeType sizeOutput;

    origin[0] = 0;
    origin[1] = 0;
    origin[2] = 0;

    sizeOutput[0] = 256;
    sizeOutput[1] = 256;
    sizeOutput[2] = 100;

    spacing.Fill(0.8);

    constantImageSource->SetOrigin(origin);
    constantImageSource->SetSpacing(spacing);
    constantImageSource->SetSize(sizeOutput);
    constantImageSource->SetConstant(0.);

    // Create reconstructed image
    ConstantImageSourceType::Pointer constantImageSource2 = ConstantImageSourceType::New();
    origin.Fill(0);
    spacing.Fill(2.);
    sizeOutput[0] = 128;
    sizeOutput[1] = 128;
    sizeOutput[2] = 128;
    constantImageSource2->SetOrigin(origin);
    constantImageSource2->SetSpacing(spacing);
    constantImageSource2->SetSize(sizeOutput);
    constantImageSource2->SetConstant(0.);

    // FDK reconstruction
    std::cout << "Reconstructing..." << std::endl;
    // Start the clock
    clock.Start();
    using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter;
    FDKGPUType::Pointer feldkamp = FDKGPUType::New();
    feldkamp->SetInput(0, constantImageSource2->GetOutput());
    feldkamp->SetInput(1, volumeData);  // Use the projection data read from the .raw file
    feldkamp->SetGeometry(geometry);
    feldkamp->GetRampFilter()->SetTruncationCorrection(0.);
    feldkamp->GetRampFilter()->SetHannCutFrequency(0.0);
    feldkamp->Update();
    // Stop the clock after the reconstruction
    clock.Stop();

    // Field-of-view masking
    using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>;
    FOVFilterType::Pointer fieldofview = FOVFilterType::New();
    fieldofview->SetInput(0, feldkamp->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(argv[1]);
    writer->SetInput(fieldofview->GetOutput());
    writer->Update();

    // Output the time taken
    std::cout << "Time taken for reconstruction: " << clock.GetTotal() << " seconds." << std::endl;
    std::cout << "Done!" << std::endl;
    return EXIT_SUCCESS;
}
Please help me to solve this problem. Thanks in advance.
Regards, Navami.
InsightSoftwareConsortium / ITK

Unable to work with Raw Projection Images in ITK-RTK #4447
