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 image 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.