[Rtk-users] Back Projection
Adarsh S Sunil
adarshto3 at gmail.com
Tue Aug 30 06:42:49 UTC 2022
Hi,
Thanks for the reply. I tried based on the suggestion and I am able to
centre the Forward Projection.
Now the forward projection looks like this:
[image: forward.PNG]
Backward Projection looks like this:
[image: backward.PNG]
I tried to use the AddImageFilter in different ways like:
using AddImageFilterType = itk::AddImageFilter<ImageType, ImageType,
ImageType>;
AddImageFilterType::Pointer addImageFilter = AddImageFilterType::New();
addImageFilter->SetInput1(changeFilter->GetOutput());
addImageFilter->SetInput2(ForwardProj->GetOutput());
addImageFilter->Update();
using AddFilterType = itk::AddImageFilter<ImageType, ImageType,
ImageType>;
AddFilterType::Pointer addFilter = AddFilterType::New();
addFilter->SetInput1(outputImage);
addFilter->SetInput2(changeFilter->GetOutput());
addFilter->Update();
But it is not working. When I use it, the program flow ends at
AddImageFilter.
What do I need to do to get the correct Backward projection? Please help.
Thanks & Regards
Adarsh S S
On Thu, 25 Aug 2022 at 02:20, Simon Rit <simon.rit at creatis.insa-lyon.fr>
wrote:
> Hi,
> To make air at 0, you need to add the opposite of the air value (1000 for
> HU) to your image, e.g. with itk::AddImageFilter.
> To shift your image, you can change its origin with
> itk::ChangeInformationImageFilter.
> Simon
>
> On Mon, Aug 22, 2022 at 5:26 PM Adarsh S Sunil <adarshto3 at gmail.com>
> wrote:
>
>> Hi Simon,
>>
>> Thanks for the support. I tried based on your suggestions.
>>
>>
>>
>> First I tried to align my volume to center by adjusting values in SetOrigin()
>> and geometry->AddProjection(sid, sdd, 0, 0, 0, angle, 0, 0 , 0).
>>
>> But it only shifted the projection even I tried to change offset for
>> source.
>>
>> Also I couldn’t find anything that shows how to make the air at 0.
>>
>> Could you please suggest a reference which shows how to do them both??
>>
>>
>>
>>
>>
>> Thanks and Regards
>>
>> Adarsh S S
>>
>> On Thu, Aug 18, 2022, 3:29 AM Simon Rit <simon.rit at creatis.insa-lyon.fr>
>> wrote:
>>
>>> I ran it and the output is not completely black. But you probably need
>>> to modify your projections to make the air at 0 (add -1024). You also may
>>> want to center your volume, the forward projections look weird currently.
>>> Simon
>>>
>>> On Wed, Aug 17, 2022 at 8:41 AM Adarsh S Sunil <adarshto3 at gmail.com>
>>> wrote:
>>>
>>>> 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/20220830/ad2083f3/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/20220830/ad2083f3/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/20220830/ad2083f3/attachment-0003.png>
More information about the Rtk-users
mailing list