[Rtk-users] Ring artifact even with homogenous valued projection

[Benjamin W. Maloney]

Hi Simon,

Thanks. This is almost certainly what the issue is
It's not obvious to me how to do this in C++

I tried:
feldkamp->GetRampFilter()->SetZeroPadFactors(1);
and
feldkamp->GetRampFilter()->SetZeroPadFactors(2);
but neither seemed to make a difference.

I didn't see any other pad functions though

Thanks!
Ben
________________________________
From: Simon Rit <simon.rit at creatis.insa-lyon.fr>
Sent: Friday, April 2, 2021 4:51 PM
To: Benjamin W. Maloney <Benjamin.W.Maloney.TH at dartmouth.edu>
Cc: rtk-users at public.kitware.com <rtk-users at public.kitware.com>
Subject: Re: [Rtk-users] Ring artifact even with homogenous valued projection

Hi,
If you pass homogeneous non-zero images, you expect such a ring artefact. The reason is that FDK assumes non-truncated projections, i.e., the projections should be 0 at the border. If not, then the ramp filtering will make a transition from the value to 0.
Probably, your projections are non-zero on the borders either. This can be because they have not been calibrated to line integrals by taking the log of the flat field divided by the measure, because they have not been well calibrated or because they are truncated (detector too small for the object). You can then mitigate the effect by padding with data, see the --pad option in rtkfdk.
Simon

On Fri, Apr 2, 2021 at 7:35 PM Benjamin W. Maloney <Benjamin.W.Maloney.TH at dartmouth.edu<mailto:Benjamin.W.Maloney.TH at dartmouth.edu>> wrote:
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
Rtk-users at public.kitware.com<mailto:Rtk-users at public.kitware.com>
https://public.kitware.com/mailman/listinfo/rtk-users<https://nam12.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpublic.kitware.com%2Fmailman%2Flistinfo%2Frtk-users&data=04%7C01%7CBenjamin.W.Maloney.TH%40dartmouth.edu%7C4d50d70e01ba4520ce5808d8f6192adf%7C995b093648d640e5a31ebf689ec9446f%7C0%7C0%7C637529935259373895%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=pq9c3k0Rh60WeSzEvtwqmsOnAL3MFmUGv1yZvz560lo%3D&reserved=0>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://public.kitware.com/pipermail/rtk-users/attachments/20210406/96943848/attachment.htm>