[Rtk-users] Back Projection
Adarsh S Sunil
adarshto3 at gmail.com
Tue Aug 16 13:22:59 UTC 2022
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/20220816/4e060deb/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/20220816/4e060deb/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/20220816/4e060deb/attachment-0003.png>
More information about the Rtk-users
mailing list