[Rtk-users] Back Projection
Adarsh S Sunil
adarshto3 at gmail.com
Fri Aug 12 12:36:41 UTC 2022
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/20220812/8395bad3/attachment-0001.htm>
More information about the Rtk-users
mailing list