SimonRit
commented
4 months ago This is a user issue, prefer the mailing list for this. To start with, you can remove the rei filter (no need to draw ellipsoid projections in your input projections). And specify the template type of ImageSeriesReader, I think you want ImageType.
I am using the "3D Cone Beam Tomography dataset of a walnut" from the Open Xray Tomographic Datasets with the following parameters: Source to Image Distance (SID): 210.66mm Source to Detector Distance (SDD): 553.74mm Number of images: 360 Scan starting angle: 0 degrees Scan ending angle: 180 degrees Image dimensions: 2240 x 2368 pixels but reconstructed slices are returning as zeros and NaN.
Detail codes information as follows: int main(int argc, char * argv) { if (argc < 3) { std::cout << "Usage: FirstReconstruction " << std::endl;
return EXIT_FAILURE;
}
// Defines the image type
using ImageType = itk::CudaImage<float, 3>;
// Defines the RTK geometry object
using GeometryType = rtk::ThreeDCircularProjectionGeometry;
GeometryType::Pointer geometry = GeometryType::New();
unsigned int numberOfProjections = 360;
double firstAngle = 0;
double angularArc = 180;
unsigned int sid = 210.66;//70;// source to isocenter distance
unsigned int sdd = 553.74; //425;; // source to detector distance
for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)
{
double angle = firstAngle + noProj angularArc / numberOfProjections;
double ProjOffsetX = 0.;
double ProjOffsetY = 0.;
double OutoffPlanAngle = 0.;
double InPlanAngle = 0.;
double SourceOffSetX = 0.;
double SourceOffSetY = 0.;
geometry->AddProjection(
sid, sdd, angle, ProjOffsetX, ProjOffsetY, OutoffPlanAngle, InPlanAngle, SourceOffSetX, SourceOffSetY);
}
// Write the geometry to disk
rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
xmlWriter->SetFilename(argv[2]);
xmlWriter->SetObject(geometry);
xmlWriter->WriteFile();
// Create a stack of empty projection images
using ConstantImageSourceType = rtk::ConstantImageSource;
ConstantImageSourceType::Pointer constantImageSource = ConstantImageSourceType::New();
ConstantImageSourceType::PointType origin;
ConstantImageSourceType::SpacingType spacing;
ConstantImageSourceType::SizeType sizeOutput;
//加载TIF图像
using Reader3DType = itk::ImageSeriesReader;
//using Reader3DType = itk::ImageFileReader;
Reader3DType::Pointer readerT = Reader3DType::New();
using NamesGeneratorType = itk::NumericSeriesFileNames;
NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
nameGenerator->SetSeriesFormat("D:/RTK/Build/test/rtkFDKCudaTest.dir/walnut_projections/20201111walnut%04d.tif");
nameGenerator->SetStartIndex(1);
nameGenerator->SetEndIndex(360);
nameGenerator->SetIncrementIndex(1);
using TIFFIOType=itk::TIFFImageIO;
TIFFIOType::Pointer tiffIO = TIFFIOType::New();
std::vector fileNames = nameGenerator->GetFileNames();
readerT->SetImageIO(tiffIO);
readerT->SetFileNames(fileNames); // fileNames
try
{
readerT->Update();
}
catch (itk::ExceptionObject & error)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << error << std::endl;
return EXIT_FAILURE;
}
// Create projections of an ellipse
using REIType = rtk::RayEllipsoidIntersectionImageFilter<ImageType, ImageType>;
REIType::Pointer rei = REIType::New();
REIType::VectorType semiprincipalaxis, center;
center[0] = -1120;
center[1] = -1184;
center[2] = -128.;
rei->SetDensity(1);
rei->SetAngle(0.);
rei->SetCenter(center);
rei->SetGeometry(geometry);
rei->SetInput(readerT->GetOutput());
// Create reconstructed image
ConstantImageSourceType::Pointer constantImageSource2 = ConstantImageSourceType::New();
sizeOutput[0] = 2048;
sizeOutput[1] = 2048;
sizeOutput[2] = 512;
//origin.Fill(-1024.);
origin[0] = -1023.5;
origin[1] = -1023.5;
origin[2] = -256;
spacing.Fill(1.);
constantImageSource2->SetOrigin(origin);
constantImageSource2->SetSpacing(spacing);
constantImageSource2->SetSize(sizeOutput);
constantImageSource2->SetConstant(1.);
// FDK reconstruction
std::cout << "Reconstructing..." << std::endl;
using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter;
FDKGPUType::Pointer feldkamp = FDKGPUType::New();
feldkamp->SetInput(0, constantImageSource2->GetOutput());
feldkamp->SetInput(1, rei->GetOutput());
feldkamp->SetGeometry(geometry);
feldkamp->GetRampFilter()->SetTruncationCorrection(0.);
feldkamp->GetRampFilter()->SetHannCutFrequency(0.0);
// Field-of-view masking using FOVFilterType = rtk::FieldOfViewImageFilter<ImageType, ImageType>; FOVFilterType::Pointer fieldofview = FOVFilterType::New(); fieldofview->SetInput(0, feldkamp->GetOutput()); fieldofview->SetProjectionsStack(rei->GetOutput()); // readerT->GetOutput() fieldofview->SetGeometry(geometry); // Writer std::cout << "Writing output image..." << std::endl; using WriterType = itk::ImageFileWriter;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(argv[1]);
writer->SetInput(fieldofview->GetOutput()); //
try
{
writer->Update();
}
catch (itk::ExceptionObject & error)
{
std::cerr << error << std::endl;
}
std::cout << "Done!" << std::endl;
return EXIT_SUCCESS;
}