Hi,
Please avoid multiple posts, one is enough (ref
https://discourse.itk.org/t/issue-in-itk-rtk-reconstruction-with-shepp-logan-phantom/6464
).
Simon
On Sat, Feb 24, 2024 at 9:17 AM Navami S <navamisnai...@gmail.com> wrote:
> Hi,
> I was trying to do a comparative study between TIGRE and RTK for which I
> have extracted the Shepp-Logan Raw projection data from TIGRE and I am
> trying to reconstruct it in RTK. Initially I tried to read the raw data as
> such but was not successful. So I generated MetaImage Header file of the
> raw data and tried to read it instead and was successful. Now my code has
> run completely and an output image is generated in mhd format but
> unfortunately when I try to visualize the output using 3D Slicer, the image
> is completely black, ie, there are no pixel values. So I tried to print the
> pixel values of both projection and reconstrued images and found that the
> pixel values of projection images had both zero and non-zero values but
> reconstructed images had both zero and negative values. So I have done some
> post processing which resulted in pixel values with zeros and values like
> 1, 2 or 3. But still no image can be seen. I don't know what to do and how
> to rectify this issue. I will provide my code below.
>
> #include <iostream>
> #include <itkCudaImage.h>
> #include <rtkConstantImageSource.h>
> #include <rtkThreeDCircularProjectionGeometryXMLFileWriter.h>
> #include <rtkCudaFDKConeBeamReconstructionFilter.h>
> #include <rtkFieldOfViewImageFilter.h>
> #include <itkImageFileReader.h>
> #include <itkTimeProbe.h>
> #include <rtkThreeDCircularProjectionGeometryXMLFileReader.h>
> #include <QtWidgets>
> #include <itkImage.h>
> #include <itkMetaImageIO.h>
> #include <itkImageIOFactory.h>
> #include <itkCastImageFilter.h>
> #include <itkIntensityWindowingImageFilter.h>
>
>
>
> int main(int argc, char **argv)
> {
> if (argc < 3)
> {
> std::cout << "Usage: FirstReconstruction <outputimage>
> <outputgeometry>" << std::endl;
> return EXIT_FAILURE;
> }
>
>
>
> // Create a time probe
> itk::TimeProbe clock;
>
> // Defines the image type
> using ImageType = itk::CudaImage<uint16_t, 3>;
>
> // Defines the RTK geometry object
> using GeometryType = rtk::ThreeDCircularProjectionGeometry;
> GeometryType::Pointer geometry = GeometryType::New();
> unsigned int numberOfProjections = 100;
>
>
> for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)
> {
> //double angle = static_cast<double>(noProj) * 360. /
> static_cast<double>(numberOfProjections);
> const int num_angles = 100;
> const double start_angle = 0.0;
> const double end_angle = 2.0 * M_PI;
>
> // Create an array of angles
> double angles[num_angles];
>
> // Calculate evenly spaced angles
> for (int i = 0; i < num_angles; ++i) {
> double fraction = static_cast<double>(i) / (num_angles - 1);
> angles[i] = start_angle + fraction * (end_angle - start_angle);
>
> geometry->AddProjection(1000, 1536, angles[i]);
> }
>
> }
> // Write the geometry to disk
> rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
> xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
> xmlWriter->SetFilename(argv[2]);
> xmlWriter->SetObject(geometry);
> xmlWriter->WriteFile();
>
>
> //using ProjectionImageType =
> rtk::CudaFDKBackProjectionImageFilter::ProjectionImageType;
> using ProjectionImageType = itk::CudaImage<uint16_t, 3>;
>
> // Read the projection data from the .mha file
> using ReaderType = itk::ImageFileReader<ImageType>;
> ReaderType::Pointer reader = ReaderType::New();
>
> // Create MetaImageIO object and set component type explicitly
> using MetaImageIOType = itk::MetaImageIO;
> MetaImageIOType::Pointer metaImageIO = MetaImageIOType::New();
> metaImageIO->SetComponentType(itk::ImageIOBase::USHORT); // Set the
> component type to unsigned short (16-bit)
> reader->SetImageIO(metaImageIO);
>
>
> reader->SetFileName("D:/Navami/TIGRE_OUT_NEW/SheppLoganTrial_projections.mhd");
> // Update the path accordingly
>
> try {
> reader->Update();
>
> // // // Accessing the image
> // ImageType::Pointer image = reader->GetOutput();
>
>
> // // Define an iterator to iterate over the image
> // itk::ImageRegionIterator<ImageType> it(image,
> image->GetLargestPossibleRegion());
>
> // // Iterate through the image and print pixel values
> // for (it.GoToBegin(); !it.IsAtEnd(); ++it) {
> // std::cout << "Pixel value: " <<
> static_cast<int>(it.Get()) << std::endl;
> // }
> // // Print projection image pixel values after reading
> // std::cout << "Projection image pixel values after reading:" <<
> std::endl;
> // itk::ImageRegionIterator<ImageType> itUint16(reader->GetOutput(),
> reader->GetOutput()->GetLargestPossibleRegion());
> // for (itUint16.GoToBegin(); !itUint16.IsAtEnd(); ++itUint16) {
> // std::cout << "Pixel value: " << itUint16.Get() << std::endl;
> // }
>
> } catch (itk::ExceptionObject &ex) {
> std::cerr << "Error reading the .mhd file: " << ex << std::endl;
> return EXIT_FAILURE;
> }
>
>
> ////// TO WRITE PROJECTION DATA //////
>
> // // Writer for projection data
> // using ProjectionWriterType = itk::ImageFileWriter<ImageType>;
> // ProjectionWriterType::Pointer projectionWriter =
> ProjectionWriterType::New();
> //
> projectionWriter->SetFileName("D:/Navami/TIGRE_OUT_NEW/projection_data.mha");
> // Set the filename for the projection data
> // projectionWriter->SetInput(reader->GetOutput()); // Set the
> input as the projection data
> // try
> // {
> // std::cout << "Writing projection data..." << std::endl;
> // projectionWriter->Update(); // Write the projection data
> // std::cout << "Projection data written successfully." <<
> std::endl;
> // }
> // catch (itk::ExceptionObject &ex)
> // {
> // std::cerr << "Error writing projection data: " << ex <<
> std::endl;
> // return EXIT_FAILURE;
> // }
>
>
>
> using CastFilterType = itk::CastImageFilter<ImageType,
> itk::CudaImage<float, 3>>;
> CastFilterType::Pointer castFilter = CastFilterType::New();
> castFilter->SetInput(reader->GetOutput());
> castFilter->Update();
>
> // Define the type of the output image after casting
> using FloatImageType = itk::CudaImage<float, 3>;
> // Get the output image after casting
> FloatImageType::Pointer castOutputImage = castFilter->GetOutput();
>
> // // Define an iterator to iterate over the output image after casting
> // using FloatIteratorType = itk::ImageRegionIterator<FloatImageType>;
> // FloatIteratorType it(castOutputImage,
> castOutputImage->GetLargestPossibleRegion());
>
> // // Print the pixel values after casting
> // std::cout << "Pixel values after casting to float:" << std::endl;
> // for (it.GoToBegin(); !it.IsAtEnd(); ++it) {
> // std::cout << "Pixel value: " << it.Get() << std::endl;
> // }
>
> // Create reconstructed image
> using ConstantImageSourceType = rtk::ConstantImageSource<ImageType>;
> ConstantImageSourceType::Pointer constantImageSource2 =
> ConstantImageSourceType::New();
>
> ConstantImageSourceType::PointType origin;
> origin.Fill(0);
>
> ConstantImageSourceType::SpacingType spacing;
> spacing.Fill(1.6);
>
> ConstantImageSourceType::SizeType sizeOutput;
> sizeOutput[0] = 128;
> sizeOutput[1] = 128;
> sizeOutput[2] = 128;
>
> constantImageSource2->SetOrigin(origin);
> constantImageSource2->SetSpacing(spacing);
> constantImageSource2->SetSize(sizeOutput);
> constantImageSource2->SetConstant(0.0);
> constantImageSource2->Update();
>
>
> // Define the type of the output image after casting for
> constantImageSource2
> using FloatConstantImageType = itk::CudaImage<float, 3>;
>
> // Define the casting filter for constantImageSource2
> using ConstantImageCastFilterType = itk::CastImageFilter<ImageType,
> FloatConstantImageType>;
> ConstantImageCastFilterType::Pointer constantImageCastFilter =
> ConstantImageCastFilterType::New();
> constantImageCastFilter->SetInput(constantImageSource2->GetOutput());
> constantImageCastFilter->Update(); // Perform the conversion
>
>
> // FDK reconstruction
> std::cout << "Reconstructing..." << std::endl;
> // Start the clock
> clock.Start();
> using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter;
> FDKGPUType::Pointer feldkamp = FDKGPUType::New();
>
> // try {
> //using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter;
> //FDKGPUType::Pointer feldkamp = FDKGPUType::New();
> feldkamp->SetInput(0, constantImageCastFilter->GetOutput());
> feldkamp->SetInput(1, castFilter->GetOutput());
> feldkamp->SetGeometry(geometry);
> feldkamp->GetRampFilter()->SetTruncationCorrection(300.);
> feldkamp->GetRampFilter()->SetHannCutFrequency(1.5);
>
>
> // FDK reconstruction
> feldkamp->Update();
>
>
>
>
>
>
> // Normalize the reconstructed image
> using OutputImageType = itk::CudaImage<float, 3>;
> using IntensityWindowingFilterType =
> itk::IntensityWindowingImageFilter<OutputImageType, ImageType>;
> IntensityWindowingFilterType::Pointer windowingFilter =
> IntensityWindowingFilterType::New();
> windowingFilter->SetInput(feldkamp->GetOutput());
> windowingFilter->SetWindowMinimum(0);
> windowingFilter->SetWindowMaximum(65535);
> windowingFilter->SetOutputMinimum(0);
> windowingFilter->SetOutputMaximum(65535);
> windowingFilter->Update();
>
> // Accessing the normalized output image
> ImageType::Pointer normalizedOutputImage =
> windowingFilter->GetOutput();
>
> // Define an iterator to iterate over the image and print pixel values
> itk::ImageRegionIterator<ImageType> it(normalizedOutputImage,
> normalizedOutputImage->GetLargestPossibleRegion());
> for (it.GoToBegin(); !it.IsAtEnd(); ++it) {
> std::cout << "Normalized Pixel value: " << it.Get() << std::endl;
> }
>
> // Stop the clock after the reconstruction
> clock.Stop();
> // Output the time taken
> std::cout << "Time taken for reconstruction: " << clock.GetTotal() << "
> seconds." << std::endl;
>
> // Define the type for the output of the FDK reconstruction
> using FloatImageType = itk::CudaImage<float, 3>;
>
> // Define the type of the output image after casting for FOV filter
> using UShortImageType = itk::CudaImage<unsigned short, 3>;
>
> // Define the type for the output of the FDK reconstruction
> // using FloatImageType = itk::CudaImage<float, 3>;
>
> // Define the casting filter for the output of the FDK reconstruction
> using FOVCastFilterType = itk::CastImageFilter<FloatImageType,
> UShortImageType>;
> FOVCastFilterType::Pointer fovCastFilter = FOVCastFilterType::New();
> fovCastFilter->SetInput(feldkamp->GetOutput());
> fovCastFilter->Update(); // Perform the conversion
>
>
>
>
> // Field-of-view masking
> using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>;
> FOVFilterType::Pointer fieldofview = FOVFilterType::New();
> fieldofview->SetInput(fovCastFilter->GetOutput());
> fieldofview->SetProjectionsStack(reader->GetOutput());
> fieldofview->SetGeometry(geometry);
> fieldofview->Update();
>
> // Writer
> std::cout << "Writing output image..." << std::endl;
> using WriterType = itk::ImageFileWriter<ImageType>;
> WriterType::Pointer writer = WriterType::New();
> writer->SetFileName(argv[1]);
> writer->SetInput(fieldofview->GetOutput());
> writer->Update();
> std::cout << "Done!" << std::endl;
>
>
> return EXIT_SUCCESS;
> }
>
> Please help me with a solution. Thanks in advance.
> _______________________________________________
> Rtk-users mailing list
> rtk-us...@openrtk.org
> https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users
>
_______________________________________________
Rtk-users mailing list
rtk-us...@openrtk.org
https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users
