[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