RTKConsortium / RTK

Reconstruction Toolkit
Apache License 2.0
246 stars 145 forks source link

Member function not found in rtk::ThreeDCircularProjectionGeometry #567

Closed navami123nair closed 12 months ago

navami123nair commented 12 months ago

Hi, I was trying to write a code for performing image reconstruction using RTK Toolkit. As a beginner, I was modifying the example code in RTK (First Reconstruction). While writing the code an issue occurred regarding "Member not found in rtk::ThreeDCircularProjectionGeometry". My code is:

// RTK includes
#include <rtkConstantImageSource.h>
#include <rtkThreeDCircularProjectionGeometryXMLFileWriter.h>
//#include <rtkRayEllipsoidIntersectionImageFilter.h>
#include <rtkFDKConeBeamReconstructionFilter.h>
#include <rtkFieldOfViewImageFilter.h>
#include <rtkThreeDCircularProjectionGeometry.h>

// ITK includes
#include <itkImageFileWriter.h>
#include <iostream>
#include <fstream>
#include <vector>

// Function to read raw data from a binary file
std::vector<float> readRawData(const std::string& filename, size_t dataSize)
{
    std::vector<float> rawData;
    rawData.resize(dataSize);

    std::ifstream file(filename, std::ios::binary);

    if (!file.is_open())
    {
        std::cerr << "Error opening file: " << filename << std::endl;
        // Handle the error, return an empty vector, or throw an exception
        return rawData;
    }

    file.read(reinterpret_cast<char*>(rawData.data()), dataSize * sizeof(float));

    if (file.fail())
    {
        std::cerr << "Error reading file: " << filename << std::endl;
        // Handle the error, return an empty vector, or throw an exception
    }

    file.close();

    return rawData;
}

int
main(int argc, char ** argv)
{

  if (argc < 3)
  {
    std::cout << "Usage: FirstReconstruction <outputimage> <outputgeometry>" << std::endl;
    return EXIT_FAILURE;
  }

  // Define geometry parameters
    const unsigned int numberOfProjections = 351;
    const unsigned int detectorRows = 1536;
    const unsigned int detectorColumns = 1536;
    const double detectorPixelSize = 0.139;
    const unsigned int imageRows = 1536;
    const unsigned int imageColumns = 1536;
    const double spacingBetweenSlices = 0.5;
    const unsigned int numberOfSlices = 512;
    const double DSD = 1690.0;
    const double DSO = 921.0;
    const double detectorShiftX = 0.0;
    const double detectorShiftY = 0.0;
    const double detectorPixelPitch = 0.139;

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

    // Defines the RTK geometry object
    using GeometryType = rtk::ThreeDCircularProjectionGeometry;
        GeometryType::Pointer geometry = GeometryType::New();
        geometry->SetRadiusCylindricalDetector(detectorColumns * detectorPixelSize / 2.0);
        geometry->SetSourceToDetectorDistance(DSD);
        geometry->SetSourceToIsocenterDistance(DSO);
        geometry->SetProjectionMatrix(1, 0, 0, 0, 1, 0);

        // Read angles from angles.txt
           std::ifstream anglesFile("D:\\Navami\\CBCT_TIGRE CODE\\Data\\CBCT\\Raw1\\angle.txt");
           if (!anglesFile.is_open())
           {
               std::cerr << "Error opening angles file." << std::endl;
               return EXIT_FAILURE;
           }

           std::vector<double> angles;
           double angle;
           while (anglesFile >> angle)
           {
               angles.push_back(angle);
           }
           anglesFile.close();

           if (angles.size() != numberOfProjections)
           {
               std::cerr << "Error: Number of angles in angles.txt does not match the expected number of projections." << std::endl;
               return EXIT_FAILURE;
           }

      for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)
        {
          double angle = noProj * 360.0 / numberOfProjections;
          geometry->AddProjection(DSD, DSO, angle, detectorShiftX, detectorShiftY, detectorPixelPitch);
        }

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

  // 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] = -127;
//  origin[1] = -127;
//  origin[2] = 0.;

  origin[0] = -0.5 * detectorColumns * detectorPixelSize;
    origin[1] = -0.5 * detectorRows * detectorPixelSize;
    origin[2] = 0.;

//  sizeOutput[0] = 128;
//  sizeOutput[1] = 128;
//  sizeOutput[2] = numberOfProjections;

    sizeOutput[0] = imageColumns;
      sizeOutput[1] = imageRows;
      sizeOutput[2] = numberOfProjections;

  //spacing.Fill(2.);
  spacing.Fill(spacingBetweenSlices);

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

  // Loop through raw data files and reconstruct
      for (unsigned int projIndex = 0; projIndex < numberOfProjections; projIndex++)
      {
          // Read raw data for each projection
          std::string rawDataFilename = "D:\\Navami\\CBCT_TIGRE CODE\\Data\\CBCT\\Raw1\\projD1_%d_" + std::to_string(projIndex) + ".bin";
          size_t projectionSize = detectorRows * detectorColumns;
          std::vector<float> rawData = readRawData(rawDataFilename, projectionSize);

          // Set the raw data as the input for ConstantImageSource
              std::copy(rawData.begin(), rawData.end(), constantImageSource->GetOutput()->GetBufferPointer());
          // Create reconstructed image
          ConstantImageSourceType::Pointer constantImageSource2 = ConstantImageSourceType::New();
          sizeOutput.Fill(imageColumns);
          origin.Fill(-0.5 * imageColumns * spacingBetweenSlices);
          spacing.Fill(spacingBetweenSlices);
          constantImageSource2->SetOrigin(origin);
          constantImageSource2->SetSpacing(spacing);
          constantImageSource2->SetSize(sizeOutput);
          constantImageSource2->SetConstant(0.);

// FDK reconstruction
 std::cout << "Reconstructing..." << std::endl;
 using FDKCPUType = rtk::FDKConeBeamReconstructionFilter<ImageType>;
 FDKCPUType::Pointer feldkamp = FDKCPUType::New();
 feldkamp->SetInput(0, constantImageSource2->GetOutput());
 feldkamp->SetInput(1, constantImageSource->GetOutput());
 feldkamp->SetGeometry(geometry);
 feldkamp->GetRampFilter()->SetTruncationCorrection(0.);
 feldkamp->GetRampFilter()->SetHannCutFrequency(0.0);

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

         // Writer for each slice
         for (unsigned int sliceIndex = 0; sliceIndex < numberOfSlices; sliceIndex++)
         {
             std::cout << "Writing output image for slice " << sliceIndex << "..." << std::endl;
             using WriterType = itk::ImageFileWriter<ImageType>;
             WriterType::Pointer writer = WriterType::New();
             writer->SetFileName("output_image_proj_" + std::to_string(projIndex) + "_slice_" + std::to_string(sliceIndex) + ".mha");

             // Extract the 2D slice using itk::ExtractImageFilter
                 using ExtractFilterType = itk::ExtractImageFilter<ImageType, ImageType>;
                 ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
                 extractFilter->SetDirectionCollapseToIdentity();
                 extractFilter->SetInput(fieldofview->GetOutput());

                             // Define the extraction region
                                 ExtractFilterType::InputImageRegionType extractionRegion;
                                 extractionRegion.SetSize(2, 0); // Set the size of the third dimension to 0
                                 extractionRegion.SetIndex(2, sliceIndex); // Set the index of the third dimension to the current slice

                                 extractFilter->SetExtractionRegion(extractionRegion);
                                 extractFilter->Update();

                 // Set the input of the writer to the extracted slice
                    writer->SetInput(extractFilter->GetOutput());
                    writer->Update();
         }
      }

  std::cout << "Done!" << std::endl;
  return EXIT_SUCCESS;
}

In this code the issue occurs while defining the RTK Geometry Object.

// Defines the RTK geometry object
    using GeometryType = rtk::ThreeDCircularProjectionGeometry;
        GeometryType::Pointer geometry = GeometryType::New();
        geometry->SetRadiusCylindricalDetector(detectorColumns * detectorPixelSize / 2.0);
        geometry->SetSourceToDetectorDistance(DSD);
        geometry->SetSourceToIsocenterDistance(DSO);
        geometry->SetProjectionMatrix(1, 0, 0, 0, 1, 0); 

The error occurring is attached below. Issue in RTK

I humbly request for any suggestions or solutions to solve this issue. Thank You.

SimonRit commented 12 months ago

Please post your questions on the mailing list. These member functions don't exist indeed, see Doxygen documentation.