[Rtk-users] Reconstruction of CBCT using CudaFDKConeBeamReconstrctionFilter

[yx deng]

Hi Simon,
*I got a better reconstruction result with your advice last week,but there
still are two mainly questions.on the one hand, as the picture1 showing,
the question that the shape of the region corresponding to the
reconstructed target changes from a rectangle to a diamond in the slice is
still exist.on the other hand,there are much noise  and artifact in the
result. the projection images used were processed with dark image and flat
image.The picture 2 is the reconstruction slice using astra toolbox with
the same iamges. * I want to get the same effect as Figure 2 by using the
RTK , please help me to point out which code i need to fix  in my code* .*
*.*[image: picture1.png]           [image: picture2.png]
                 picture1
 picture2

The codes as follows:
  int CT_ReconAndDetec(string projPath,string  ReconPath,
    int    beginAng,int    endAng,int    rotateAngle,
    int    projNum, int    DetecX,int    DetecZ,
    int    nx,      int    ny,    int    nz,
    float  pixel_size){
   using ImageType = itk::CudaImage<float, 3>;
  using ProjImageType=itk::Image<unsigned short,3>;
  using GeometryType = rtk::ThreeDCircularProjectionGeometry;
  GeometryType::Pointer geometry = GeometryType::New();
  unsigned int          numberOfProjections = 1201;
  double                firstAngle = 0;
  double                angularArc = 360;

  unsigned int sid = 66;
  unsigned int sdd = 199;
  for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)
  {
    double angle = firstAngle + noProj * angularArc / numberOfProjections;
    double ProjOffsetX = 0.;
    double ProjOffsetY = 0.;
    double OutoffPlanAngle = 0;
    double InPlanAngle = 0;
    double SourceOffSetX = 0.;
    double SourceOffSetY = 0.;
    geometry->AddProjection(sid, sdd, angle, ProjOffsetX, ProjOffsetY,
              OutoffPlanAngle, InPlanAngle, SourceOffSetX, SourceOffSetY);
  }

  using NamesGeneratorType = itk::NumericSeriesFileNames;
  NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
  nameGenerator->SetSeriesFormat(projPath+"/scan_%06d.tif");
  nameGenerator->SetStartIndex(0);
  nameGenerator->SetEndIndex(1200);
  nameGenerator->SetIncrementIndex(1);
  std::vector<std::string> fileNames = nameGenerator->GetFileNames();

 using Reader3DType = rtk::ProjectionsReader<ProjImageType>;
  Reader3DType::Pointer readerT = Reader3DType::New();
  using TIFFIOType = itk::TIFFImageIO;
  TIFFIOType::Pointer tiffIO = TIFFIOType::New();

  readerT->SetImageIO(tiffIO);
  readerT->SetFileNames(fileNames);
  try
  {
    readerT->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << error << std::endl;
    return EXIT_FAILURE;
  }
  using CastFilterType = itk::CastImageFilter<ProjImageType, ImageType>;
  CastFilterType::Pointer castFilter = CastFilterType::New();
  castFilter->SetInput(readerT->GetOutput());
  try
  {
    castFilter->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error casting the image: " << error << std::endl;
    return EXIT_FAILURE;
  }

  ImageType::Pointer inputImage = castFilter->GetOutput();
  ImageType::PointType newOrigin;
  newOrigin[0]= -DetecX / 2.0 * pixel_size;
  newOrigin[1] = -DetecZ / 2.0 * pixel_size;
  newOrigin[2] = 0.;
  inputImage->SetOrigin(newOrigin);

 ImageType::DirectionType newDirection;
 newDirection.SetIdentity()
 inputImage->SetDirection(newDirection);

  ImageType::SpacingType newSpacing;
  newSpacing[0] = pixel_size;
  newSpacing[1] = pixel_size;
  newSpacing[2] = pixel_size;
  inputImage->SetSpacing(newSpacing);

  using ConstantImageSourceType = rtk::ConstantImageSource<ImageType>;
  ConstantImageSourceType::Pointer     constantImageSource =
ConstantImageSourceType::New();
  ConstantImageSourceType::PointType   origin;
  ConstantImageSourceType::SpacingType spacing;
  ConstantImageSourceType::SizeType    sizeOutput;
  sizeOutput[0] = nx;
  sizeOutput[1] = ny;
  sizeOutput[2] = nz;
  spacing.Fill(pixel_size);
  origin[0] = -(sizeOutput[0] / 2.0)*spacing[0];
  origin[1] = -(sizeOutput[1] / 2.0)*spacing[1];
  origin[2] = -(sizeOutput[2] / 2.0)*spacing[2];
  constantImageSource->SetOrigin(origin);
  constantImageSource->SetSpacing(spacing);
  constantImageSource->SetSize(sizeOutput);
  constantImageSource->SetConstant(0.5);

  using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter;
  FDKGPUType::Pointer feldkamp = FDKGPUType::New();
   feldkamp->SetInput(0, constantImageSource->GetOutput());
  feldkamp->SetInput(1, inputImage);
  feldkamp->SetGeometry(geometry);
  feldkamp->GetRampFilter()->SetTruncationCorrection(0.5);
  feldkamp->GetRampFilter()->SetHannCutFrequency(0.86);

  using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>;
  FOVFilterType::Pointer fieldofview = FOVFilterType::New();
  fieldofview->SetInput(0, feldkamp->GetOutput());
  fieldofview->SetProjectionsStack(inputImage); // readerT->GetOutput()
  fieldofview->SetGeometry(geometry);

  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(ReconPath);
  writer->SetInput(fieldofview->GetOutput());
  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << error << std::endl;
  }
   std::cout << "Done!" << std::endl;
  return 0;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20240625/b1e9d239/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: picture1.png
Type: image/png
Size: 170811 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20240625/b1e9d239/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: picture2.png
Type: image/png
Size: 299476 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20240625/b1e9d239/attachment-0003.png>