Hi,
It might be caused by the line
double projectionsOrigin[3] = { 0.0 };
With this setting, unless you specify an offset in x and y of your
detector in the geometry, the source-isocenter ray will point to the
lower-left corner of your projections, so you may end up in a setting
where none of the rays go through your volume. To get centered
projections, you should set their origin to (Size - 1) * Spacing / 2.
Obviously only the first two dimensions matter, the third is unused.
Have a look at the metadata of the projections generated by wiki
examples, e.g. http://wiki.openrtk.org/index.php/RTK/Scripts/FDK : their
origin is not [0,0,0].
Regards,
Cyril
On 22/06/2018 16:32, Slartibartfast Fjords wrote:
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++;
}
_______________________________________________
Rtk-users mailing list
Rtk-users@public.kitware.com
https://public.kitware.com/mailman/listinfo/rtk-users
_______________________________________________
Rtk-users mailing list
Rtk-users@public.kitware.com
https://public.kitware.com/mailman/listinfo/rtk-users
