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 <adarsh...@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....@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 <adarsh...@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@public.kitware.com
>>> https://public.kitware.com/cgi-bin/mailman/listinfo/rtk-users
>>>
>>
_______________________________________________
Rtk-users mailing list
Rtk-users@public.kitware.com
https://public.kitware.com/cgi-bin/mailman/listinfo/rtk-users
