[Rtk-users] Back-projection filter outputs only zeros

[Slartibartfast Fjords]

Hi all,

I've been trying for a while to use RTK in my own software to perform
reconstructions. I've looked at the applications and examples to understand
how to set it up but I can't get it working.

I have my own data object which is stored in a single array. I used ITK's
Import Image Filter to import the data in RTK. Right now I'm trying to get
everything working with my own Shepp-Logan projections. So my pipeline
looks like this for now:
- Create the geometry
- Create a Constant Image Source
- Import my projection data
- Set up the backprojection
- Take the output of the back-projection and export it to my own data object

The problem is that all that I get as the output is zeros everywhere, on
the whole volume. I've tried skipping the back-projection to see if it was
my data that was the problem, but then it works fine and I just get back my
projections as the output. I've tried a few different back-projection
methods, but it's always giving me zeros as the output.

Has anybody experienced this before? What did I miss in my setup? I'll put
a bit of my code below. My guess is that I didn't properly set the input
image with the projections.

Thanks a lot for your time!

Adam

Code:

// Geometry
rtk::ThreeDCircularProjectionGeometry::Pointer geometry;
geometry = rtk::ThreeDCircularProjectionGeometry::New();
for (uint32_t p = 0; p < numberOfProjections; p++)
{
double angle = beginAngle + (double)p * (endAngle - beginAngle) /
(double)numberOfProjections;
geometry->AddProjection(sid, sdd, angle, isox, isoy, outAngle, inAngle,
sourceOffsetX, sourceOffsetY);
}
TRY_AND_EXIT_ON_ITK_EXCEPTION(geometry->Update());

// Constant Image Source
ConstantImageSourceType::SizeType reconstructionSize;
reconstructionSize[0] = outputChannelSizeVector.getX();
reconstructionSize[1] = outputChannelSizeVector.getY();
reconstructionSize[2] = outputChannelSizeVector.getZ();

ConstantImageSourceType::PointType origin;
origin[0] = outputOriginVector.getX();
origin[1] = outputOriginVector.getY();
origin[2] = outputOriginVector.getZ();

ConstantImageSourceType::SpacingType reconstructionSpacing;
reconstructionSpacing[0] = outputSpacingVector.getX();
reconstructionSpacing[1] = outputSpacingVector.getY();
reconstructionSpacing[2] = outputSpacingVector.getZ();

ConstantImageSourceType::Pointer constantImageSource =
ConstantImageSourceType::New();
constantImageSource->SetOrigin(origin);
constantImageSource->SetSpacing(reconstructionSpacing);
constantImageSource->SetSize(reconstructionSize);
constantImageSource->SetConstant(0.0);

// ITK Import Image Filter
ImportFilterType::SizeType projectionsSize;
projectionsSize[0] = projectionsXSize;
projectionsSize[1] = projectionsYSize;
projectionsSize[2] = projectionsZSize;

ImportFilterType::SpacingType projectionsSpacing;
projectionsSpacing[0] = projectionsXSpacing;
projectionsSpacing[1] = projectionsYSpacing;
projectionsSpacing[2] = projectionsZSpacing; // not that it matters much...

const bool importImageFilterWillOwnTheBuffer = false;

double projectionsOrigin[3] = { 0.0 };

ImportFilterType::Pointer importFilter = ImportFilterType::New();
ImportFilterType::IndexType importFilterStartIndex;
importFilterStartIndex.Fill(0);

ImportFilterType::RegionType region;
region.SetIndex(importFilterStartIndex);
region.SetSize(projectionsSize);

importFilter->SetRegion(region);
importFilter->SetOrigin(projectionsOrigin);
importFilter->SetSpacing(projectionsSpacing);

PixelType *projectionsDataBuffer = (PixelType *)projectionsChannelFloat.
getTRawDataChunk(0);
importFilter->SetImportPointer(projectionsDataBuffer, projectionsSize[0]*
projectionsSize[1] * projectionsSize[2], importImageFilterWillOwnTheBuffer);

TRY_AND_EXIT_ON_ITK_EXCEPTION(importFilter->Update());

// Backrprojection setup
rtk::BackProjectionImageFilter<OutputImageType, OutputImageType>::Pointer
bp;
bp = rtk::FDKBackProjectionImageFilter<OutputImageType,
OutputImageType>::New();

bp->SetInput(0, constantImageSource->GetOutput());
bp->SetInput(1, importFilter->GetOutput());
bp->SetGeometry(geometry);
TRY_AND_EXIT_ON_ITK_EXCEPTION(bp->Update());

OutputImageType::PixelContainer::Element * iterator;
iterator = bp->GetOutput()->GetBufferPointer();

PixelType *outputBuffer = (PixelType *)outputChannel.getTRawDataChunk(0);
uint64_t numberOfPixels = xSize * ySize * zSize;
for (uint64_t i = 0; i < numberOfPixels; i++)
{
outputBuffer [i] = *iterator;
iterator++;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://public.kitware.com/pipermail/rtk-users/attachments/20180622/505eea34/attachment.html>