Hi all,
I'm running into an issue where I'm getting a large ring artifact at the edge
of my reconstructions. It is in the same place regardless of the input images
(including homogeneous images), so it is something I'm doing wrong on the
reconstruction side rather than a defective pixel on the detector.
I'm happy to share my input images and output reconstructions if that helps
Relevant code:
using GeometryType = rtk::ThreeDCircularProjectionGeometry;
GeometryType::Pointer geometry = GeometryType::New();
unsigned int numberOfProjections = (steps_taken+1);
float firstAngle = 0.;
float angularArc = 360.132;
const double sid = 306.6; // source to isocenter distance, 306.6mm, 0.3 pixel
size, image binned by 4 in x and y
const double sdd = 561.3; // 561.3 source to detector distance
const double offsetX_d = 0.;
const double offsetY_d = 0.;
const double outOfPlaneAngle = -2.55;
const double inPlaneAngle = 0.0;
const double offsetX_s = 0.;
const double offsetY_s = 0.;
for (unsigned int noProj = 0; noProj < numberOfProjections; noProj++)
{
double angle = 205 + ((double)noProj * ((double)angularArc /
((double)steps_taken + 1))); //Adjust the first number (23 before +) to rotate
output of reconstruction
geometry->AddProjection(sid, sdd, angle, offsetX_d, offsetY_d, outOfPlaneAngle,
inPlaneAngle, offsetX_s, offsetY_s);
}
// Write the geometry to disk
rtk::ThreeDCircularProjectionGeometryXMLFileWriter::Pointer xmlWriter;
xmlWriter = rtk::ThreeDCircularProjectionGeometryXMLFileWriter::New();
xmlWriter->SetFilename("..\\..\\..\\data\\geometry_newSystem_full.xml");
xmlWriter->SetObject(geometry);
xmlWriter->WriteFile();
// Defines the image type
using ImageType = itk::CudaImage< float, 3 >;
// Projection Images
using ConstantImageSourceType = rtk::ConstantImageSource< ImageType >;
ConstantImageSourceType::PointType origin;
ConstantImageSourceType::SpacingType spacing;
ConstantImageSourceType::SizeType sizeOutput;
int size_x = 3072 / binn;
int size_y = 864 / binn;
int size_z = numberOfProjections;
//Binned x4
origin[0] = -0.5 * (size_x - 1) * 0.3/4.*binn;
origin[1] = -0.5 * (size_y - 1) * 0.3/4.*binn;
origin[2] = 0.;
sizeOutput[0] = size_x;
sizeOutput[1] = size_y;
sizeOutput[2] = numberOfProjections;
spacing.Fill(0.3/4.*binn);
// Load in my data
typedef float PixelType;
const unsigned char Dimension = 3;
typedef itk::CudaImage< PixelType, Dimension > ImageType;
ImageType::Pointer image;
typedef itk::ImageFileReader< ImageType > ReaderType;
ReaderType::Pointer reader;
reader = ReaderType::New();
reader->SetFileName(projection_file_r); // bin x4
reader->Update();
reader->GetOutput()->SetOrigin(origin);
reader->GetOutput()->SetSpacing(spacing);
// Create reconstructed image
ConstantImageSourceType::Pointer constantImageSource2 =
ConstantImageSourceType::New();
// x4
sizeOutput.Fill(size_x);
sizeOutput[1] = size_y;
origin.Fill(-0.5 * (size_x - 1) * 0.3/4.*binn*sid / sdd);
origin[1] = (-0.5 * (size_y - 1) * 0.3/4.*binn*sid / sdd);
spacing.Fill(0.3/4.*binn*sid / sdd);
constantImageSource2->SetOrigin(origin);
constantImageSource2->SetSpacing(spacing);
constantImageSource2->SetSize(sizeOutput);
constantImageSource2->SetConstant(0.);
// FDK reconstruction
//GPU
using FDKGPUType = rtk::CudaFDKConeBeamReconstructionFilter;
FDKGPUType::Pointer feldkamp = FDKGPUType::New();
feldkamp->SetInput(0, constantImageSource2->GetOutput());
feldkamp->SetInput(1, reader->GetOutput());
feldkamp->SetGeometry(geometry);
using WriterType = itk::ImageFileWriter< ImageType >;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(reconstruction_file_w);
writer->SetInput(feldkamp->GetOutput());
writer->Update();
_______________________________________________
Rtk-users mailing list
[email protected]
https://public.kitware.com/mailman/listinfo/rtk-users
