[Rtk-users] Output Error
Navami S
navamisnair30 at gmail.com
Tue Dec 26 09:53:30 CET 2023
Hi,
I was trying to do image reconstruction with cbct data using RTK
toolkit (CPP). Well I wrote a code and executed it, but unfortunately I am
not getting the output correctly. I am having doubts regarding the geometry
I am creating. The code written and geometry file are attached below.
Please help me with a solution. Thanks in advance.
Regards,
Navami
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20231226/c1fa23bd/attachment-0001.htm>
-------------- next part --------------
#include <QCoreApplication>
#include <QImage>
#include <QLabel>
#include <QVBoxLayout>
#include <QPixmap>
#include <QApplication>
#include <cmath> // for M_PI
#include <QScrollArea>
#include <algorithm>
#include <QDir>
#include <QStringList>
//#define M_PI 3.14159265358979323846
//#include <itkVersion.h>
// RTK includes
#include <rtkThreeDCircularProjectionGeometry.h>
#include <rtkConstantImageSource.h>
#include <rtkThreeDCircularProjectionGeometryXMLFileWriter.h>
#include <rtkFDKConeBeamReconstructionFilter.h>
#include <rtkFieldOfViewImageFilter.h>
// ITK includes
#include <itkImageFileWriter.h>
#include <itkExtractImageFilter.h>
#include <itkImageFileReader.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <stdexcept> // for std::runtime_error
#include <sstream> // For std::istringstream
#include <cstdint> // for uint16_t
#include <itkPasteImageFilter.h>
#include <itkImageRegionIterator.h>
int
main(int argc, char ** argv)
{
//std::cout << "ITK Version: " << itk::Version::GetITKVersion() << std::endl;
if (argc < 3)
{
std::cout << "Usage: FirstReconstruction <outputimage> <outputgeometry>" << std::endl;
return EXIT_FAILURE;
}
//QString outputImagePath = QString(argv[1]);
// QString outputGeometryPath = QString(argv[2]);
QDir outputDir(QFileInfo(argv[1]).absoluteDir());
if (!outputDir.exists() && !outputDir.mkpath(outputDir.path()))
{
std::cerr << "Error creating output directory: " << outputDir.path().toStdString() << std::endl;
return EXIT_FAILURE;
}
QApplication app(argc, argv);
// // Create a widget to hold the images
QWidget widget;
QVBoxLayout layout(&widget); // Declare 'layout' outside the loop
QVBoxLayout mainLayout;
// // Define geometry parameters
const unsigned int numberOfProjections = 351;
const unsigned int numberOfSlices = 512;
const unsigned int detectorRows = 1536;
const unsigned int detectorColumns = 1536;
const double detectorPixelSize = 0.139;
const unsigned int imageRows = 512;
const unsigned int imageColumns = 512;
const double spacingBetweenSlices = 0.5;
const double SourceToDetectorDistance = 1690.0;
const double SourceToIsocenterDistance = 921.0;
const double detectorShiftX = 0.0;
const double detectorShiftY = 0.0;
//const double detectorPixelPitch = 0.139;
// Defines the image type
using ImageType = itk::Image<unsigned short, 3>;
using GeometryType = rtk::ThreeDCircularProjectionGeometry;
GeometryType::Pointer geometry = GeometryType::New();
// // Read angles from angles.txt
std::ifstream anglesFile("D:\\Navami\\CBCT_TIGRE CODE\\Data\\CBCT\\angles\\angle.txt");
if (!anglesFile.is_open())
{
std::cerr << "Error opening angles file." << std::endl;
return EXIT_FAILURE;
}
std::vector<double> angles;
double angle;
// while (anglesFile >> angle)
// {
// //angles.push_back(angle);
// angles.push_back(angle * M_PI / 180.0); // Convert degrees to radians
// }
// anglesFile.close();
// std::vector<double> angles;
// double angle;
while (anglesFile >> angle) {
angles.push_back(angle);
}
anglesFile.close();
if (angles.size() != numberOfProjections)
{
std::cerr << "Error: Number of angles in angles.txt does not match the expected number of projections." << std::endl;
return EXIT_FAILURE;
}
for (unsigned int noProj = 1; noProj <= numberOfProjections; noProj++)
{
double angle = noProj * 360.0 / numberOfProjections;
geometry->AddProjection( SourceToIsocenterDistance, SourceToDetectorDistance, angle, detectorShiftX, detectorShiftY/*, detectorPixelPitch*/);
}
// // Write the geometry to disk
rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
xmlWriter->SetFilename(argv[2]);
xmlWriter->SetObject(geometry);
xmlWriter->WriteFile();
// // Set up FDK reconstruction
using ConstantImageSourceType = rtk::ConstantImageSource<ImageType>;
ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New();
ConstantImageSourceType::PointType origin;
ConstantImageSourceType::SpacingType spacing;
ConstantImageSourceType::SizeType sizeOutput;
sizeOutput[0] = imageColumns;
sizeOutput[1] = imageRows;
sizeOutput[2] = numberOfProjections;
// sizeOutput[2] = numberOfSlices;
// spacing.Fill(spacingBetweenSlices);
spacing[0] = 0.5; // Adjust the spacing based on your requirements
spacing[1] = 0.5;
spacing[2] = 0.5;
// To place the origin at the exact center of the image
origin[0] = -0.5 * (imageColumns - 1) * spacing[0];
origin[1] = -0.5 * (imageRows - 1) * spacing[1];
origin[2] = 0.0;
std::cout << "Before setting in constantImageSource:" << std::endl;
std::cout << " Origin: " << origin << std::endl;
std::cout << " Spacing: " << spacing << std::endl;
std::cout << " Size: " << sizeOutput << std::endl;
constantImageSource->SetOrigin(origin);
constantImageSource->SetSpacing(spacing);
constantImageSource->SetSize(sizeOutput);
constantImageSource->SetConstant(0.);
constantImageSource-> Update();
// Get output spacing from constantImageSource
ConstantImageSourceType::SpacingType outputSpacingArray = constantImageSource->GetOutput()->GetSpacing();
std::cout << "After setting in constantImageSource:" << std::endl;
std::cout << " Size: " << constantImageSource->GetOutput()->GetLargestPossibleRegion().GetSize() << std::endl;
//std::cout << " Spacing: " << constantImageSource->GetOutput()->GetSpacing() << std::endl;
std::cout << " Origin: " << constantImageSource->GetOutput()->GetOrigin() << std::endl;
// Print the output spacing values
std::cout << "Output Spacing Array: " << outputSpacingArray << std::endl;
// Create a 3D volume to store all slices
using VolumeType = itk::Image<unsigned short, 3>;
VolumeType::Pointer volume = VolumeType::New();
// Set the size and spacing of the volume based on the slices
VolumeType::SizeType volumeSize;
VolumeType::SpacingType volumeSpacing;
volumeSize[0] = imageColumns;
volumeSize[1] = imageRows;
volumeSize[2] = numberOfSlices;
volumeSpacing[0] = spacing[0];
volumeSpacing[1] = spacing[1];
volumeSpacing[2] = spacingBetweenSlices;
volume->SetRegions(volumeSize);
volume->SetSpacing(volumeSpacing);
volume->Allocate();
// Loop through raw data files and reconstruct
for (unsigned int projIndex = 1; projIndex <= numberOfProjections; projIndex++)
{
// // Read raw data for each projection
std::string rawDataFilename = "D:\\Navami\\CBCT_TIGRE CODE\\Data\\CBCT\\Raw1\\projD1_" + std::to_string(projIndex) + ".raw";
size_t projectionSize = detectorRows * detectorColumns;
// // std::cout <<rawDataFilename<<std::endl;
// std::vector<uint16_t> rawData = readRawData(rawDataFilename, projectionSize);
// Open the file
std::ifstream file(rawDataFilename, std::ios::binary);
if (!file.is_open())
{
throw std::runtime_error("Error opening file: " + rawDataFilename);
}
// Get the file size
file.seekg(0, std::ios::end);
size_t fileSize = file.tellg();
// Calculate the expected size
size_t expectedSize = projectionSize * sizeof(uint16_t);
std::cout << "File size: " << fileSize << ", Expected size: " << expectedSize << std::endl;
if (fileSize != expectedSize)
{
file.close(); // Close the file before throwing an exception
throw std::runtime_error("Invalid file size for: " + rawDataFilename);
}
// Allocate vector based on file size
std::vector<uint16_t> rawData(projectionSize);
// Read data into the vector
file.seekg(0, std::ios::beg);
file.read(reinterpret_cast<char*>(rawData.data()), expectedSize);
// Close the file
file.close();
// Set the raw data as input for FDK reconstruction
constantImageSource->SetConstant(rawData[0]); // Assuming rawData is a 1D vector
// Create reconstructed image
ConstantImageSourceType::Pointer constantImageSource2 = ConstantImageSourceType::New();
sizeOutput.Fill(imageColumns);
origin.Fill(-0.5 * (imageColumns - 1) * spacingBetweenSlices);
spacing.Fill(spacingBetweenSlices);
std::cout << "Before setting in constantImageSource2:" << std::endl;
std::cout << " Origin: " << origin << std::endl;
std::cout << " Spacing: " << spacing << std::endl;
std::cout << " Size: " << sizeOutput << std::endl;
constantImageSource2->SetOrigin(origin);
constantImageSource2->SetSpacing(spacing);
constantImageSource2->SetSize(sizeOutput);
constantImageSource2->SetConstant(0.);
constantImageSource2-> Update();
// // FDK reconstruction
std::cout << "Reconstructing..." << std::endl;
using FDKCPUType = rtk::FDKConeBeamReconstructionFilter<ImageType>;
FDKCPUType::Pointer feldkamp = FDKCPUType::New();
feldkamp->SetInput(0, constantImageSource2->GetOutput());
feldkamp->SetInput(1, constantImageSource->GetOutput());
feldkamp->SetGeometry(geometry);
feldkamp->GetRampFilter()->SetTruncationCorrection(0.);
feldkamp->GetRampFilter()->SetHannCutFrequency(0.0);
feldkamp->Update();
/
// // Field-of-view masking
using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>;
FOVFilterType::Pointer fieldofview = FOVFilterType::New();
fieldofview->SetInput(0, feldkamp->GetOutput());
fieldofview->SetProjectionsStack(constantImageSource2->GetOutput());
fieldofview->SetGeometry(geometry);
fieldofview->Update();
//feldkamp->GetRampFilter()->SetTruncationCorrection(0.);
// feldkamp->GetRampFilter()->SetHannCutFrequency(0.0);
// // Writer for each slice
for (unsigned int sliceIndex = 1; sliceIndex <= numberOfSlices; sliceIndex++)
{
std::cout << "Writing output image for slice " << sliceIndex << "..." << std::endl;
// Extract the 2D slice using itk::ExtractImageFilter
using ExtractFilterType = itk::ExtractImageFilter<ImageType, ImageType>;
ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
extractFilter->SetDirectionCollapseToIdentity();
extractFilter->SetInput(fieldofview->GetOutput());
// Define the extraction region
ExtractFilterType::InputImageRegionType extractionRegion;
std::cout << "Image Dimensionality: " << fieldofview->GetOutput()->GetImageDimension() << std::endl;
extractionRegion.SetSize(0, imageColumns); // Set the size for the first dimension
extractionRegion.SetSize(1, imageRows); // Set the size for the second dimension
extractionRegion.SetSize(2, 1); // Set the size for the third dimension (slice index)
extractionRegion.SetIndex(0, 0); // Set the index of the first dimension to the current slice
extractionRegion.SetIndex(1, 0); // Set the index of the second dimension to 0
extractionRegion.SetIndex(2, sliceIndex); // Set the index of the third dimension to the current slice
std::cout << "Before SetExtractionRegion" << std::endl;
// Add print statements for debugging
std::cout << "Extraction Region Size: " << extractionRegion.GetSize() << std::endl;
std::cout << "Extraction Region Index: " << extractionRegion.GetIndex() << std::endl;
extractFilter->SetExtractionRegion(extractionRegion);
//extractFilter->Update();
extractFilter->UpdateLargestPossibleRegion();
// Copy the slice to the 3D volume
itk::ImageRegionIterator<ImageType> sliceIt(extractFilter->GetOutput(), extractFilter->GetOutput()->GetLargestPossibleRegion());
itk::ImageRegionIterator<VolumeType> volumeIt(volume, volume->GetLargestPossibleRegion());
while (!sliceIt.IsAtEnd() && !volumeIt.IsAtEnd())
{
volumeIt.Set(sliceIt.Get());
++sliceIt;
++volumeIt;
}
/*//std::string outputFileName = QString(argv[1]).append("/slice").append(QString::number(sliceIndex)).append(".mha").toStdString();
using WriterType = itk::ImageFileWriter<ImageType>;
WriterType::Pointer writer = WriterType::New();
// writer->SetFileName(outputFileName);
writer->SetFileName(argv[1]);
writer->SetInput(extractFilter->GetOutput());
writer->Update()*/;
// Save the 3D volume to a single file
std::string outputFileName = outputDir.absolutePath().toStdString() + "/output_volume.mha";
using VolumeWriterType = itk::ImageFileWriter<VolumeType>;
VolumeWriterType::Pointer volumeWriter = VolumeWriterType::New();
volumeWriter->SetFileName(outputFileName);
volumeWriter->SetInput(volume);
volumeWriter->Update();
// extractFilter->SetExtractionRegion(extractionRegion);
//std::cout << "After SetExtractionRegion" << std::endl;
}
// // // Display the QImage in a QLabel
// QImage sliceImage(imageColumns, imageRows, QImage::Format_Grayscale8);
// //std::cout << "Done!" << std::endl;
// // int y; // Declaration of the variable y
// for (int y = 0; y < sliceImage.height(); ++y)
// // std::cout << "Done!" << std::endl;
// {
// for (int x = 0; x < sliceImage.width(); ++x)
// {
// uint16_t pixelValue = fieldofview->GetOutput()->GetPixel({ x, y, 0 });
// uint8_t pixelValue8Bit = static_cast<uint8_t>((pixelValue * 255) / std::numeric_limits<uint16_t>::max());
// sliceImage.setPixel(x, y, qRgb(pixelValue8Bit, pixelValue8Bit, pixelValue8Bit));
// }
// }
// QLabel* sliceLabel = new QLabel;
// sliceLabel->setPixmap(QPixmap::fromImage(sliceImage));
// layout.addWidget(sliceLabel);
// scrollLayout->addWidget(sliceLabel);
// mainLayout.addWidget(sliceLabel);
// // Save the image to a file
// QString outputFileName = outputDir.absolutePath() + "/slice" + QString::number(projIndex) + ".png";
// sliceImage.save(outputFileName);
// }
// // // Set the layout of the scroll widget
// scrollWidget->setLayout(scrollLayout);
// // Set the widget to the scroll area
// scrollArea.setWidget(scrollWidget);
// // Set up the main widget layout
// QVBoxLayout mainLayout;
// mainLayout.addWidget(&scrollArea);
// // Set the main widget layout
// widget.setLayout(&mainLayout);
// widget.show();
// // Create a scroll area to hold the images
// QScrollArea scrollArea;
// QWidget* scrollWidget = new QWidget;
// QVBoxLayout* scrollLayout = new QVBoxLayout(scrollWidget); // Add this line
// scrollLayout->addLayout(&mainLayout); // Add this line
// scrollWidget->setLayout(scrollLayout); // Add this line
// scrollArea.setWidget(scrollWidget);
// scrollArea.show();
// return app.exec();
}
//std::cout << "Done!!!" << std::endl;
return EXIT_SUCCESS;
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Geometry_CBCT.xml
Type: text/xml
Size: 130678 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20231226/c1fa23bd/attachment-0001.xml>
More information about the Rtk-users
mailing list