[Rtk-users] Back Projection
Simon Rit
simon.rit at creatis.insa-lyon.fr
Wed Aug 17 21:59:15 UTC 2022
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/20220817/1235ebec/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/1235ebec/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/1235ebec/attachment-0003.png>
More information about the Rtk-users
mailing list