[Rtk-users] Back Projection
Simon Rit
simon.rit at creatis.insa-lyon.fr
Tue Aug 16 13:35:28 UTC 2022
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/20220816/f5f04af3/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/f5f04af3/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/f5f04af3/attachment-0003.png>
More information about the Rtk-users
mailing list