[Rtk-users] Back Projection
Adarsh S Sunil
adarshto3 at gmail.com
Mon Aug 22 15:26:13 UTC 2022
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/20220822/0e2a7282/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/20220822/0e2a7282/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/20220822/0e2a7282/attachment-0003.png>
More information about the Rtk-users
mailing list