Issue in ITK-RTK reconstruction with Shepp-Logan phantom

[navami123nair]

Hi,
I was trying to do a comparative study between TIGRE and RTK for which I have extracted the Shepp-Logan Raw projection data from TIGRE and I am trying to reconstruct it in RTK. Initially I tried to read the raw data as such but was not successful. So i generated MetaImage Header file of the raw data and tried to read it instead and was successful. Now my code has run completely and an output image is generated in mhd format but unfortunately when I try to visualize the output using 3D Slicer, the image is completely black, ie, there are no pixel values. So i tried to print the pixel values of both projection and reconstrued images and found that the pixel values of projection images had both zero and non-zero values but reconstructed image had both zero and negative values. So I have done some post processing which resulted in pixel values with zeros and values like 1, 2 or 3. But still no image can be seen. I don't know what to do and how to rectify this issue. I will provide my code below.

#include <iostream>
#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 <QtWidgets>
#include <itkImage.h>
#include <itkMetaImageIO.h>
#include <itkImageIOFactory.h>
#include <itkCastImageFilter.h>
#include <itkIntensityWindowingImageFilter.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<uint16_t, 3>;

    // Defines the RTK geometry object
    using GeometryType = rtk::ThreeDCircularProjectionGeometry;
    GeometryType::Pointer geometry = GeometryType::New();
    unsigned int numberOfProjections = 100;


    for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)
    {
        //double angle = static_cast<double>(noProj) * 360. / static_cast<double>(numberOfProjections);
        const int num_angles = 100;
        const double start_angle = 0.0;
        const double end_angle = 2.0 * M_PI;

        // Create an array of angles
        double angles[num_angles];

        // Calculate evenly spaced angles
        for (int i = 0; i < num_angles; ++i) {
            double fraction = static_cast<double>(i) / (num_angles - 1);
            angles[i] = start_angle + fraction * (end_angle - start_angle);

            geometry->AddProjection(1000, 1536, angles[i]);
        }

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


    //using ProjectionImageType = rtk::CudaFDKBackProjectionImageFilter::ProjectionImageType;
    using ProjectionImageType = itk::CudaImage<uint16_t, 3>;

    // Read the projection data from the .mha file
    using ReaderType = itk::ImageFileReader<ImageType>;
    ReaderType::Pointer reader = ReaderType::New();

    // Create MetaImageIO object and set component type explicitly
    using MetaImageIOType = itk::MetaImageIO;
    MetaImageIOType::Pointer metaImageIO = MetaImageIOType::New();
    metaImageIO->SetComponentType(itk::ImageIOBase::USHORT); // Set the component type to unsigned short (16-bit)
    reader->SetImageIO(metaImageIO);

    reader->SetFileName("D:/Navami/TIGRE_OUT_NEW/SheppLoganTrial_projections.mhd"); // Update the path accordingly

    try {
        reader->Update();

//        //        // Accessing the image
//        ImageType::Pointer image   = reader->GetOutput();


//        // Define an iterator to iterate over the image
//               itk::ImageRegionIterator<ImageType> it(image, image->GetLargestPossibleRegion());

//               // Iterate through the image and print pixel values
//               for (it.GoToBegin(); !it.IsAtEnd(); ++it) {
//                   std::cout << "Pixel value: " << static_cast<int>(it.Get()) << std::endl;
//               }
//        // Print projection image pixel values after reading
//        std::cout << "Projection image pixel values after reading:" << std::endl;
//        itk::ImageRegionIterator<ImageType> itUint16(reader->GetOutput(), reader->GetOutput()->GetLargestPossibleRegion());
//        for (itUint16.GoToBegin(); !itUint16.IsAtEnd(); ++itUint16) {
//            std::cout << "Pixel value: " << itUint16.Get() << std::endl;
//        }

    } catch (itk::ExceptionObject &ex) {
        std::cerr << "Error reading the .mhd file: " << ex << std::endl;
        return EXIT_FAILURE;
    }


////// TO WRITE PROJECTION DATA //////

//    // Writer for projection data
//       using ProjectionWriterType = itk::ImageFileWriter<ImageType>;
//       ProjectionWriterType::Pointer projectionWriter = ProjectionWriterType::New();
//       projectionWriter->SetFileName("D:/Navami/TIGRE_OUT_NEW/projection_data.mha"); // Set the filename for the projection data
//       projectionWriter->SetInput(reader->GetOutput());      // Set the input as the projection data
//       try
//       {
//           std::cout << "Writing projection data..." << std::endl;
//           projectionWriter->Update(); // Write the projection data
//           std::cout << "Projection data written successfully." << std::endl;
//       }
//       catch (itk::ExceptionObject &ex)
//       {
//           std::cerr << "Error writing projection data: " << ex << std::endl;
//           return EXIT_FAILURE;
//       }



    using CastFilterType = itk::CastImageFilter<ImageType, itk::CudaImage<float, 3>>;
    CastFilterType::Pointer castFilter = CastFilterType::New();
    castFilter->SetInput(reader->GetOutput());
    castFilter->Update();

    // Define the type of the output image after casting
    using FloatImageType = itk::CudaImage<float, 3>;
    // Get the output image after casting
    FloatImageType::Pointer castOutputImage = castFilter->GetOutput();

//    // Define an iterator to iterate over the output image after casting
//    using FloatIteratorType = itk::ImageRegionIterator<FloatImageType>;
//    FloatIteratorType it(castOutputImage, castOutputImage->GetLargestPossibleRegion());

//    // Print the pixel values after casting
//    std::cout << "Pixel values after casting to float:" << std::endl;
//    for (it.GoToBegin(); !it.IsAtEnd(); ++it) {
//        std::cout << "Pixel value: " << it.Get() << std::endl;
//    }

    // Create reconstructed image
    using ConstantImageSourceType = rtk::ConstantImageSource<ImageType>;
    ConstantImageSourceType::Pointer constantImageSource2 = ConstantImageSourceType::New();

    ConstantImageSourceType::PointType origin;
    origin.Fill(0);

    ConstantImageSourceType::SpacingType spacing;
    spacing.Fill(1.6);

    ConstantImageSourceType::SizeType sizeOutput;
    sizeOutput[0] = 128;
    sizeOutput[1] = 128;
    sizeOutput[2] = 128;

    constantImageSource2->SetOrigin(origin);
    constantImageSource2->SetSpacing(spacing);
    constantImageSource2->SetSize(sizeOutput);
    constantImageSource2->SetConstant(0.0);
    constantImageSource2->Update();


    // Define the type of the output image after casting for constantImageSource2
    using FloatConstantImageType = itk::CudaImage<float, 3>;

    // Define the casting filter for constantImageSource2
    using ConstantImageCastFilterType = itk::CastImageFilter<ImageType, FloatConstantImageType>;
    ConstantImageCastFilterType::Pointer constantImageCastFilter = ConstantImageCastFilterType::New();
    constantImageCastFilter->SetInput(constantImageSource2->GetOutput());
    constantImageCastFilter->Update(); // Perform the conversion


    // FDK reconstruction
    std::cout << "Reconstructing..." << std::endl;
    // Start the clock
    clock.Start();
    using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter;
    FDKGPUType::Pointer feldkamp = FDKGPUType::New();

     // try {
    //using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter;
    //FDKGPUType::Pointer feldkamp = FDKGPUType::New();
    feldkamp->SetInput(0, constantImageCastFilter->GetOutput());
    feldkamp->SetInput(1, castFilter->GetOutput());
    feldkamp->SetGeometry(geometry);
    feldkamp->GetRampFilter()->SetTruncationCorrection(300.);
    feldkamp->GetRampFilter()->SetHannCutFrequency(1.5);


    // FDK reconstruction
    feldkamp->Update();






    // Normalize the reconstructed image
        using OutputImageType = itk::CudaImage<float, 3>;
        using IntensityWindowingFilterType = itk::IntensityWindowingImageFilter<OutputImageType, ImageType>;
        IntensityWindowingFilterType::Pointer windowingFilter = IntensityWindowingFilterType::New();
        windowingFilter->SetInput(feldkamp->GetOutput());
        windowingFilter->SetWindowMinimum(0);
        windowingFilter->SetWindowMaximum(65535);
        windowingFilter->SetOutputMinimum(0);
        windowingFilter->SetOutputMaximum(65535);
        windowingFilter->Update();

        // Accessing the normalized output image
        ImageType::Pointer normalizedOutputImage = windowingFilter->GetOutput();

        // Define an iterator to iterate over the image and print pixel values
        itk::ImageRegionIterator<ImageType> it(normalizedOutputImage, normalizedOutputImage->GetLargestPossibleRegion());
        for (it.GoToBegin(); !it.IsAtEnd(); ++it) {
            std::cout << "Normalized Pixel value: " << it.Get() << std::endl;
        }

    // Stop the clock after the reconstruction
    clock.Stop();
    // Output the time taken
    std::cout << "Time taken for reconstruction: " << clock.GetTotal() << " seconds." << std::endl;

    // Define the type for the output of the FDK reconstruction
    using FloatImageType = itk::CudaImage<float, 3>;

    // Define the type of the output image after casting for FOV filter
    using UShortImageType = itk::CudaImage<unsigned short, 3>;

    // Define the type for the output of the FDK reconstruction
    //  using FloatImageType = itk::CudaImage<float, 3>;

    // Define the casting filter for the output of the FDK reconstruction
    using FOVCastFilterType = itk::CastImageFilter<FloatImageType, UShortImageType>;
    FOVCastFilterType::Pointer fovCastFilter = FOVCastFilterType::New();
    fovCastFilter->SetInput(feldkamp->GetOutput());
    fovCastFilter->Update(); // Perform the conversion




    // Field-of-view masking
    using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>;
    FOVFilterType::Pointer fieldofview = FOVFilterType::New();
    fieldofview->SetInput(fovCastFilter->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(argv[1]);
    writer->SetInput(fieldofview->GetOutput());
    writer->Update();
    std::cout << "Done!" << std::endl;


    return EXIT_SUCCESS;
}

Please help me with a solution. Thanks in advance.

[SimonRit]

User issues are addressed on the mailing list or discourse.