[Rtk-users] Back Projection
Simon Rit
simon.rit at creatis.insa-lyon.fr
Wed Aug 24 20:50:12 UTC 2022
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/20220824/4ff76798/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/20220824/4ff76798/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/20220824/4ff76798/attachment-0003.png>
More information about the Rtk-users
mailing list