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

[Cyril Mory]

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 at public.kitware.com
> https://public.kitware.com/mailman/listinfo/rtk-users

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://public.kitware.com/pipermail/rtk-users/attachments/20180622/acadc2d9/attachment-0001.html>