[Rtk-users] Back Projection
Adarsh S Sunil
adarshto3 at gmail.com
Wed Aug 17 06:40:36 UTC 2022
Hi,
I’m attaching my code here. I’m trying to get the Forward
projection and backward projection from dicom data.
CT input :
https://drive.google.com/drive/folders/1Vc3HfqvfBdtviuVdbBhDQev2ZDEH7NYX?usp=sharing
Usage: Run from output location on cmd as
FirstCudaReconstruction.exe <outputFile.mha> <geoemtry.xml> <CT DICOM
directory location>
Please check it and let me know the problems in it.
Thanks and Regards
Adarsh SS
Thanks & Regards
Adarsh S S
On Tue, 16 Aug 2022 at 19:05, Simon Rit <simon.rit at creatis.insa-lyon.fr>
wrote:
> Hi,
> Can you maybe share a complete script with an input CT? I also don't
> understand what's going on but we never get a full picture of what you're
> doing.
> Simon
>
> On Tue, Aug 16, 2022 at 3:23 PM Adarsh S Sunil <adarshto3 at gmail.com>
> wrote:
>
>> Hi,
>>
>> Thank you for the answer.
>>
>>
>>
>> The following is the output for Forward Projection (ForwardProj->
>> GetOutput()) :
>>
>>
>> [image: image002.png]
>>
>>
>>
>> As you suggested I tried with changeFilter->GetOutput() instead of outputImage1
>> as follows:
>>
>>
>>
>> using BackFilter = rtk::CudaFDKConeBeamReconstructionFilter;
>>
>> auto BackProjFilter = BackFilter::New();
>>
>> // FDK reconstruction
>>
>> std::cout << "Reconstructing..." << std::endl;
>>
>> BackProjFilter->SetInput(0, changeFilter->GetOutput());
>>
>> BackProjFilter->SetInput(1, ForwardProj->GetOutput());
>>
>> BackProjFilter->SetGeometry(geometry);
>>
>> BackProjFilter->Update();
>>
>> BackProjFilter->UpdateOutputInformation();
>>
>>
>>
>> // Writer
>>
>> std::cout << "Writing output image..." << std::endl;
>>
>> using WriterType = itk::ImageFileWriter<ImageType>;
>>
>> WriterType::Pointer writer = WriterType::New();
>>
>> writer->SetFileName(argv[1]);
>>
>> writer->SetInput(BackProjFilter->GetOutput());
>>
>> writer->Update();
>>
>>
>>
>> But I’m still getting the output BackProjFilter->GetOutput() as a black
>> image.
>>
>>
>> [image: image004.png]
>>
>>
>>
>> I checked the value for changeFilter->GetOutput()and it gives the dicom
>> data.
>>
>>
>>
>> I’m using ITK version 5.2.1 and RTK version 2.0. Can you please help me
>> understand why I’m not getting a back projection?
>>
>> Are there any extra options to set for getting the back projection? What
>> do I have to do to get the correct back projection?
>>
>>
>>
>>
>>
>> Thanks and Regards
>>
>> Adarsh S S
>>
>>
>>
>>
>>
>> On Sat, Aug 13, 2022, 10:43 AM Simon Rit <simon.rit at creatis.insa-lyon.fr>
>> wrote:
>>
>>> Your approach is correct and there is no need to do multiple files. The
>>> outputImage1 (for backprojection) should have the same information as the
>>> tomography (changeFilter->GetOutput()), not the same information as the
>>> projections (outputImage). Maybe that's the problem here?
>>> Simon
>>>
>>> On Fri, Aug 12, 2022 at 2:37 PM Adarsh S Sunil <adarshto3 at gmail.com>
>>> wrote:
>>>
>>>> Hi,
>>>>
>>>>
>>>>
>>>> Thanks for the reply. The outputImage passed as input to
>>>> the backprojection does not have the same information as
>>>> reader->GetOutput().
>>>>
>>>>
>>>>
>>>> I’m writing the forwardProjection in to a single file named Forward.mha.
>>>>
>>>>
>>>>
>>>> rtk::ForwardProjectionImageFilter<ImageType, ImageType>::Pointer ForwardProj
>>>> =
>>>>
>>>> rtk::CudaForwardProjectionImageFilter<ImageType, ImageType>::New();
>>>>
>>>> dynamic_cast<rtk::CudaForwardProjectionImageFilter<ImageType,
>>>> ImageType> *>(ForwardProj.GetPointer())->SetStepSize(1);
>>>>
>>>> ForwardProj->SetInput(0, outputImage);
>>>>
>>>> ForwardProj->SetInput(1, changeFilter->GetOutput());
>>>>
>>>> ForwardProj->SetGeometry(geometry);
>>>>
>>>> ForwardProj->Update();
>>>>
>>>>
>>>>
>>>> // Writer
>>>>
>>>> std::cout << "Writing output image..." << std::endl;
>>>>
>>>> using WriterType1 = itk::ImageFileWriter<ImageType>;
>>>>
>>>> WriterType1::Pointer writer1 = WriterType1::New();
>>>>
>>>> writer1->SetFileName("Forward.mha");
>>>>
>>>> writer1->SetInput(ForwardProj->GetOutput());
>>>>
>>>> writer1->Update();
>>>>
>>>>
>>>>
>>>> Also, I’m setting the same regions, spacing, origin etc. to forward
>>>> and backward projections as follows:
>>>>
>>>>
>>>>
>>>> using RegionType = itk::ImageRegion<3>;
>>>>
>>>> ImageType::Pointer inputImage = reader->GetOutput();
>>>>
>>>> RegionType inputRegion =
>>>> inputImage->GetLargestPossibleRegion();
>>>>
>>>> ImageType::Pointer outputImage = ImageType::New();
>>>>
>>>> ImageType::Pointer outputImage1 = ImageType::New();
>>>>
>>>>
>>>>
>>>> using changeImageFilter = itk::ChangeInformationImageFilter<ImageType
>>>> >;
>>>>
>>>> auto changeFilter = changeImageFilter::New();
>>>>
>>>> double rotationX = 90;
>>>>
>>>> itk::Versor<double> rotation;
>>>>
>>>> const double angleInRadians = rotationX * vnl_math::pi /
>>>> 180.0;
>>>>
>>>> rotation.SetRotationAroundY(angleInRadians);
>>>>
>>>> const ImageType::DirectionType direction =
>>>> inputImage->GetDirection();
>>>>
>>>> const ImageType::DirectionType newDirection = direction *
>>>> rotation.GetMatrix();
>>>>
>>>> changeFilter->SetOutputDirection(newDirection);
>>>>
>>>> changeFilter->ChangeAll();
>>>>
>>>> changeFilter->SetInput(inputImage);
>>>>
>>>> changeFilter->UpdateOutputInformation();
>>>>
>>>> changeFilter->Update();
>>>>
>>>>
>>>>
>>>> const unsigned int NumberOfProjectionImages = 180;
>>>>
>>>> ImageType::IndexType start;
>>>>
>>>> start[0] = inputRegion.GetIndex()[0];
>>>>
>>>> start[1] = inputRegion.GetIndex()[1];
>>>>
>>>> start[2] = inputRegion.GetIndex()[2];
>>>>
>>>>
>>>>
>>>> RegionType::SizeType size;
>>>>
>>>>
>>>>
>>>> size[0] = inputRegion.GetSize()[0];
>>>>
>>>> size[1] = inputRegion.GetSize()[1];
>>>>
>>>> size[2] = NumberOfProjectionImages;
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> RegionType region;
>>>>
>>>> region.SetSize(size);
>>>>
>>>> region.SetIndex(start);
>>>>
>>>> ImageType::SpacingType spacing;
>>>>
>>>>
>>>>
>>>> spacing[0] = inputImage->GetSpacing()[0];
>>>>
>>>> spacing[1] = inputImage->GetSpacing()[1];
>>>>
>>>> spacing[2] = inputImage->GetSpacing()[2];
>>>>
>>>>
>>>>
>>>> ImageType::PointType origin;
>>>>
>>>>
>>>>
>>>> origin[0] = inputImage->GetOrigin()[0];
>>>>
>>>> origin[1] = inputImage->GetOrigin()[1];
>>>>
>>>> origin[2] = inputImage->GetOrigin()[2];
>>>>
>>>>
>>>>
>>>> //For Forward Projection
>>>>
>>>> outputImage->SetRegions(region);
>>>>
>>>> outputImage->SetSpacing(spacing);
>>>>
>>>> outputImage->SetOrigin(origin);
>>>>
>>>> outputImage->Allocate();
>>>>
>>>> outputImage->FillBuffer(size[0] * size[1] * size[2]);
>>>>
>>>>
>>>>
>>>> //For Backward Projection
>>>>
>>>> outputImage1->SetRegions(region);
>>>>
>>>> outputImage1->SetSpacing(spacing);
>>>>
>>>> outputImage1->SetOrigin(origin);
>>>>
>>>> outputImage1->Allocate();
>>>>
>>>> outputImage1->FillBuffer(size[0] * size[1] * size[2]);
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> I tried to give the output of forward projection directly. But it did
>>>> not work. I was only getting a blank image.
>>>>
>>>> So I tried a different approach. I read somewhere that using Projection
>>>> reader is better. So, I’m trying using projection readers.
>>>>
>>>> Then I put the Forward Projection reader (ForwardProj->GetOutput()) as
>>>> the input of Backward Projection.
>>>>
>>>>
>>>>
>>>> using FileNamesContainer = std::vector<std::string>;
>>>>
>>>> FileNamesContainer fileNames =
>>>> nameGenerator->GetFileNames(seriesIdentifier);
>>>>
>>>>
>>>>
>>>> using ReaderType1 = rtk::ProjectionsReader<ImageType>;
>>>>
>>>> ReaderType1::Pointer reader1 = ReaderType1::New();
>>>>
>>>> reader1->SetFileNames(fileNames);
>>>>
>>>> reader1->Update();
>>>>
>>>> using BackFilter = rtk::CudaFDKConeBeamReconstructionFilter;
>>>>
>>>> auto BackProjFilter = BackFilter::New();
>>>>
>>>> // FDK reconstruction
>>>>
>>>> std::cout << "Reconstructing..." << std::endl;
>>>>
>>>> BackProjFilter->SetInput(0, outputImage1);
>>>>
>>>> BackProjFilter->SetInput(1, ForwardProj->GetOutput());
>>>>
>>>> BackProjFilter->SetGeometry(geometry);
>>>>
>>>> BackProjFilter->Update();
>>>>
>>>>
>>>>
>>>> // Writer
>>>>
>>>> std::cout << "Writing output image..." << std::endl;
>>>>
>>>> using WriterType = itk::ImageFileWriter<ImageType>;
>>>>
>>>> WriterType::Pointer writer = WriterType::New();
>>>>
>>>> writer->SetFileName("BackImage.mha");
>>>>
>>>> writer->SetInput(outputImage1);
>>>>
>>>> writer->Update();
>>>>
>>>>
>>>>
>>>> Still I’m not getting any valid output.
>>>>
>>>> Is my approach to get backward projection from forward projection
>>>> correct?
>>>>
>>>> Do we have to write forward and backward projections to multiple files?
>>>> Then how to do it?
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> Thanks and regards
>>>>
>>>> Adarsh S S
>>>>
>>>>
>>>>
>>>>
>>>> On Fri, Aug 12, 2022, 2:16 PM Simon Rit <simon.rit at creatis.insa-lyon.fr>
>>>> wrote:
>>>>
>>>>> Hi,
>>>>> It should work. UpdateOutputInformation is not required (it's done by
>>>>> Update already). Does the outputImage passed as input to the backprojection
>>>>> have the same information as reader->GetOutput()?
>>>>> Simon
>>>>>
>>>>> On Fri, Aug 12, 2022 at 8:57 AM Adarsh S Sunil <adarshto3 at gmail.com>
>>>>> wrote:
>>>>>
>>>>>>
>>>>>> Hi all,
>>>>>>
>>>>>> I'm a beginner in rtk. I want to use the forward and backward
>>>>>> projection in CT images. For that I defined the RTK geometry object and
>>>>>> seriesIdContainer.
>>>>>>
>>>>>> Then I read CT data and assign it to a reader and set it as
>>>>>> inputImage.
>>>>>>
>>>>>>
>>>>>> * while (seriesItr != seriesUID.end())*
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> * { seriesIdentifier = seriesItr->c_str(); seriesItr++;
>>>>>> std::cout << "\nReading: "; std::cout << seriesIdentifier << std::endl;
>>>>>> using FileNamesContainer = std::vector<std::string>;
>>>>>> FileNamesContainer fileNames =
>>>>>> nameGenerator->GetFileNames(seriesIdentifier); using ImageIOType =
>>>>>> itk::GDCMImageIO; auto dicomIO = ImageIOType::New();
>>>>>> reader->SetImageIO(dicomIO); reader->SetFileNames(fileNames);
>>>>>> reader->ForceOrthogonalDirectionOff(); // properly read CTs with gantry
>>>>>> tilt TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->Update()) }*
>>>>>>
>>>>>>
>>>>>>
>>>>>> * using RegionType = itk::ImageRegion<3>; ImageType::Pointer
>>>>>> inputImage = reader->GetOutput(); RegionType inputRegion =
>>>>>> inputImage->GetLargestPossibleRegion(); ImageType::Pointer outputImage =
>>>>>> ImageType::New();*
>>>>>>
>>>>>> Then I created Imageregion and imageFilter for outputImage.
>>>>>> Then I set the regions, spacing, origin etc.
>>>>>>
>>>>>> Then I get the Forward Projection like the following:
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> * rtk::ForwardProjectionImageFilter<ImageType, ImageType>::Pointer
>>>>>> ForwardProj = rtk::CudaForwardProjectionImageFilter<ImageType,
>>>>>> ImageType>::New();
>>>>>> dynamic_cast<rtk::CudaForwardProjectionImageFilter<ImageType, ImageType>
>>>>>> *>(ForwardProj.GetPointer()) ->SetStepSize( 1 ); ForwardProj->SetInput( 0,
>>>>>> outputImage ); ForwardProj->SetInput( 1, changeFilter->GetOutput());
>>>>>> ForwardProj->SetGeometry( geometry ); ForwardProj->Update();*
>>>>>>
>>>>>> Then I wrote the ForwardProj to an mha file.
>>>>>>
>>>>>> Then I put the output from Forward Projection as input of Back
>>>>>> Projection.
>>>>>> Now I am trying to get the back projection like the following:
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> * using ReaderType1 = rtk::ProjectionsReader<ImageType>;
>>>>>> ReaderType1::Pointer reader1 = ReaderType1::New(); reader1->SetFileNames(
>>>>>> fileNames ); reader1->Update(); using BackFilter =
>>>>>> rtk::CudaFDKConeBeamReconstructionFilter; auto BackProjFilter =
>>>>>> BackFilter::New(); // FDK reconstruction BackProjFilter->SetInput( 0,
>>>>>> outputImage ); BackProjFilter->SetInput( 1, reader1->GetOutput() );
>>>>>> BackProjFilter->SetGeometry( geometry ); BackProjFilter->Update();
>>>>>> BackProjFilter->UpdateOutputInformation();*
>>>>>>
>>>>>> But I'm not getting the back projection. Are there any additional
>>>>>> options I have to set??
>>>>>> Is there anything I have missed??
>>>>>>
>>>>>> Thanks & Regards
>>>>>> Adarsh S S
>>>>>> _______________________________________________
>>>>>> Rtk-users mailing list
>>>>>> Rtk-users at public.kitware.com
>>>>>> https://public.kitware.com/cgi-bin/mailman/listinfo/rtk-users
>>>>>>
>>>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/rtk-users/attachments/20220817/84d15d9e/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image002.png
Type: image/png
Size: 119970 bytes
Desc: not available
URL: <http://public.kitware.com/pipermail/rtk-users/attachments/20220817/84d15d9e/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image004.png
Type: image/png
Size: 61241 bytes
Desc: not available
URL: <http://public.kitware.com/pipermail/rtk-users/attachments/20220817/84d15d9e/attachment-0003.png>
-------------- next part --------------
// RTK includes
#include <rtkConstantImageSource.h>
#include <rtkThreeDCircularProjectionGeometryXMLFileWriter.h>
#include <rtkRayEllipsoidIntersectionImageFilter.h>
#include <rtkCudaFDKConeBeamReconstructionFilter.h>
#include <rtkFieldOfViewImageFilter.h>
#include "rtkForwardProjectionImageFilter.h"
#include "rtkCudaForwardProjectionImageFilter.h"
#include <itkGDCMSeriesFileNames.h>
#include <itkImageSeriesReader.h>
#include <rtkCudaFDKBackProjectionImageFilter.h>
#include <rtkCudaBackProjectionImageFilter.h>
#include <rtkCudaFDKConeBeamReconstructionFilter.h>
// ITK includes
#include <itkImageFileWriter.h>
#include <itkGDCMImageIO.h>
#include <itkVersor.h>
#include <itkChangeInformationImageFilter.h>
#include <vnl/vnl_math.h>
int
main(int argc, char ** argv)
{
if (argc < 3)
{
std::cout << "Usage: FirstReconstruction <outputimage> <outputgeometry>" << std::endl;
return EXIT_FAILURE;
}
// 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 = 245;
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(sid, sdd, 0, 0, 0, angle, 0);
}
// Write the geometry to disk
rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
xmlWriter->SetFilename(argv[2]);
xmlWriter->SetObject(geometry);
xmlWriter->WriteFile();
using PixelType = signed short;
constexpr unsigned int Dimension = 3;
using NamesGeneratorType = itk::GDCMSeriesFileNames;
auto nameGenerator = NamesGeneratorType::New();
std::string dirName = "."; // current directory by default
if (argc >= 3)
{
dirName = argv[3];
}
nameGenerator->SetUseSeriesDetails(true);
nameGenerator->AddSeriesRestriction("0008|0021");
nameGenerator->SetGlobalWarningDisplay(false);
nameGenerator->SetDirectory(dirName);
using SeriesIdContainer = std::vector<std::string>;
const SeriesIdContainer & seriesUID = nameGenerator->GetSeriesUIDs();
auto seriesItr = seriesUID.begin();
auto seriesEnd = seriesUID.end();
if (seriesItr != seriesEnd)
{
std::cout << "The directory: ";
std::cout << dirName << std::endl;
std::cout << "Contains the following DICOM Series: ";
std::cout << std::endl;
}
else
{
std::cout << "No DICOMs in: " << dirName << std::endl;
return EXIT_SUCCESS;
}
while (seriesItr != seriesEnd)
{
std::cout << seriesItr->c_str() << std::endl;
++seriesItr;
}
using ReaderType = itk::ImageSeriesReader<ImageType>;
auto reader = ReaderType::New();
seriesItr = seriesUID.begin();
std::string seriesIdentifier;
while (seriesItr != seriesUID.end())
{
seriesIdentifier = seriesItr->c_str();
seriesItr++;
std::cout << "\nReading: ";
std::cout << seriesIdentifier << std::endl;
using FileNamesContainer = std::vector<std::string>;
FileNamesContainer fileNames = nameGenerator->GetFileNames(seriesIdentifier);
using ImageIOType = itk::GDCMImageIO;
auto dicomIO = ImageIOType::New();
reader->SetImageIO(dicomIO);
reader->SetFileNames(fileNames);
reader->ForceOrthogonalDirectionOff(); // properly read CTs with gantry tilt
TRY_AND_EXIT_ON_ITK_EXCEPTION(reader->Update())
}
using RegionType = itk::ImageRegion<3>;
ImageType::Pointer inputImage = reader->GetOutput();
RegionType inputRegion = inputImage->GetLargestPossibleRegion();
ImageType::Pointer outputImage = ImageType::New();
ImageType::Pointer outputImage1 = ImageType::New();
using changeImageFilter = itk::ChangeInformationImageFilter<ImageType>;
auto changeFilter = changeImageFilter::New();
double rotationX = 90;
itk::Versor<double> rotation;
const double angleInRadians = rotationX * vnl_math::pi / 180.0;
rotation.SetRotationAroundY(angleInRadians);
const ImageType::DirectionType direction = inputImage->GetDirection();
const ImageType::DirectionType newDirection = direction * rotation.GetMatrix();
changeFilter->SetOutputDirection(newDirection);
changeFilter->ChangeAll();
changeFilter->SetInput(inputImage);
changeFilter->UpdateOutputInformation();
changeFilter->Update();
const unsigned int NumberOfProjectionImages = 245;
ImageType::IndexType start;
start[0] = inputRegion.GetIndex()[0];
start[1] = inputRegion.GetIndex()[1];
start[2] = inputRegion.GetIndex()[2];
RegionType::SizeType size;
size[0] = inputRegion.GetSize()[0];
size[1] = inputRegion.GetSize()[1];
size[2] = NumberOfProjectionImages;
RegionType region;
region.SetSize(size);
region.SetIndex(start);
ImageType::SpacingType spacing;
spacing[0] = inputImage->GetSpacing()[0];
spacing[1] = inputImage->GetSpacing()[1];
spacing[2] = inputImage->GetSpacing()[2];
ImageType::PointType origin;
origin[0] = inputImage->GetOrigin()[0];
origin[1] = inputImage->GetOrigin()[1];
origin[2] = inputImage->GetOrigin()[2];
outputImage->SetRegions(region);
outputImage->SetSpacing(spacing);
outputImage->SetOrigin(origin);
outputImage->Allocate();
outputImage->FillBuffer(size[0] * size[1] * size[2]);
outputImage1->SetRegions(region);
outputImage1->SetSpacing(spacing);
outputImage1->SetOrigin(origin);
outputImage1->Allocate();
outputImage1->FillBuffer(size[0] * size[1] * size[2]);
rtk::ForwardProjectionImageFilter<ImageType, ImageType>::Pointer ForwardProj =
rtk::CudaForwardProjectionImageFilter<ImageType, ImageType>::New();
dynamic_cast<rtk::CudaForwardProjectionImageFilter<ImageType, ImageType> *>(ForwardProj.GetPointer())->SetStepSize(1);
ForwardProj->SetInput(0, outputImage);
ForwardProj->SetInput(1, changeFilter->GetOutput());
ForwardProj->SetGeometry(geometry);
ForwardProj->Update();
// Writer
std::cout << "Writing output image..." << std::endl;
using WriterType1 = itk::ImageFileWriter<ImageType>;
WriterType1::Pointer writer1 = WriterType1::New();
writer1->SetFileName("Forward.mha");
writer1->SetInput(ForwardProj->GetOutput());
writer1->Update();
using FileNamesContainer = std::vector<std::string>;
FileNamesContainer fileNames = nameGenerator->GetFileNames(seriesIdentifier);
using ReaderType1 = rtk::ProjectionsReader<ImageType>;
ReaderType1::Pointer reader1 = ReaderType1::New();
reader1->SetFileNames(fileNames);
reader1->Update();
using BackFilter = rtk::CudaFDKConeBeamReconstructionFilter;
auto BackProjFilter = BackFilter::New();
// FDK reconstruction
std::cout << "Reconstructing..." << std::endl;
BackProjFilter->SetInput(0, changeFilter->GetOutput());
BackProjFilter->SetInput(1, ForwardProj->GetOutput());
BackProjFilter->SetGeometry(geometry);
BackProjFilter->Update();
BackProjFilter->UpdateOutputInformation();
// Writer
std::cout << "Writing output image..." << std::endl;
using WriterType = itk::ImageFileWriter<ImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(argv[1]);
writer->SetInput(BackProjFilter->GetOutput());
writer->Update();
std::cout << "Done!" << std::endl;
return EXIT_SUCCESS;
}
More information about the Rtk-users
mailing list