DVFs interpolation by BSpline


(wonny) #1

Hi, ITK users

I really really need your help, please help me…
I am trying to do BSpline interpolation of DVFs obtained by image
registrations of 10 images.
I’m trying to modify the BSplineScatteredDataPointSetToImageFilter from the
Insight Journal. However,
I am having trouble to find the correct parameter in order to make the
filter work.
I’ll have scattered 3-D points of several DVFs over time and I want to fit
these DVFs by a 4-D Bspline. Afterwards I want to evaluate output DVF to one
time step.

Here you can find the code
Please help me…please.

Thanks in advance.

#include “itkBSplineScatteredDataPointSetToImageFilter.h”
#include “itkPointSet.h”
#include “itkImage.h”
#include “itkVectorImage.h”
#include “itkImageFileReader.h”
#include “itkImageFileWriter.h”

int main (int argc, char* argv[])
{
const unsigned int ImageDimension = 3;

using VectorType = itk::Vector<float, ImageDimension>;
using InputImageType = itk::Image<VectorType, ImageDimension>;
using InputPointType = InputImageType::PointType;
using OutputImageType = itk::Image<VectorType, ImageDimension + 1>;
using PointSetType = itk::PointSet<VectorType, ImageDimension + 1>;
using OutputPointType = PointSetType::PointType;
using ReaderType = itk::ImageFileReader<InputImageType>;
using WriterType = itk::ImageFileWriter<OutputImageType>;

PointSetType::Pointer pointSet = PointSetType::New();
unsigned long pointId = 0;
InputImageType::PointType origin;
InputImageType::SpacingType spacing;
InputImageType::SizeType size;
for (int i = 1; i < argc; ++i)
{
ReaderType::Pointer reader = ReaderType::New();
std::cout << "Reading file " << argv[i] << std::endl;
reader->SetFileName(argv[i]);
try
{
reader->Update();
}
catch (itk::ExceptionObject & err)
{
std::cout << “ExceptionObject caught !” << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
using IteratorType = itk::ImageRegionConstIterator<InputImageType>;
const InputImageType * image = reader->GetOutput();
IteratorType it(image, image->GetBufferedRegion());
it.GoToBegin();
InputPointType inPoint;
OutputPointType outPoint;
int iCount = 0;
while (!it.IsAtEnd())
{
inPoint = image->GetPixel(it.GetIndex());
for (int j = 0; j < ImageDimension; ++j)
outPoint[j] = inPoint[i];
outPoint[ImageDimension] = i - 1;
pointSet->SetPoint(pointId, outPoint);
// Transfer the pixel data to the value associated with the point.
pointSet->SetPointData(pointId, it.Get());
++it;
++pointId;
++iCount;
}
std::cout << "Number of points in " << argv[i] << " = " << iCount <<
std::endl;
std::cout << "Total number of points = " <<
pointSet->GetNumberOfPoints() << std::endl;
origin = image->GetOrigin();
spacing = image->GetSpacing();
size = image->GetLargestPossibleRegion().GetSize();
}

typedef itk::BSplineScatteredDataPointSetToImageFilter < PointSetType,
OutputImageType > SplineFilterType;
SplineFilterType::Pointer splineFilter = SplineFilterType::New();

int splineorder=3; // complexity of the spline

SplineFilterType::ArrayType ncontrol;
ncontrol[0]=splineorder + 1;
SplineFilterType::ArrayType closedim;
closedim[0]= 0;

OutputImageType::PointType parametricDomainOrigin;
OutputImageType::SpacingType parametricDomainSpacing;
OutputImageType::SizeType parametricDomainSize;
for (int i = 0; i < ImageDimension; ++i)
{
parametricDomainOrigin[i] = origin[i];
parametricDomainSpacing[i] = spacing[i];
parametricDomainSize[i] = size[i];
}
parametricDomainOrigin[ImageDimension] = 0;
parametricDomainSize[ImageDimension] = argc - 2;
parametricDomainSpacing[ImageDimension] = 10.0;

splineFilter->SetGenerateOutputImage( true ); // the only reason to turn
this off is if one only wants to use the control point lattice for further
processing
splineFilter->SetInput ( pointSet );
splineFilter->SetSplineOrder ( splineorder );
splineFilter->SetNumberOfControlPoints ( ncontrol );
splineFilter->SetNumberOfLevels( 3 );
splineFilter->SetCloseDimension ( closedim );
splineFilter->SetSize( parametricDomainSize );
splineFilter->SetSpacing( parametricDomainSpacing );
splineFilter->SetOrigin( parametricDomainOrigin );
std::cout << “Before update spline filter” << std::endl;
splineFilter->Update();
std::cout << “After update spline filter” << std::endl;

WriterType::Pointer writer = WriterType::New();
writer->SetInput(splineFilter->GetOutput());
writer->SetFileName(“Output.mhd”);
writer->Update();
std::cout << “After write image filter” << std::endl;

return EXIT_SUCCESS;
};


(Dženan Zukić) #2

To easily reproduce your problem please provide the input file too. If it crashes with different input files, please provide the smallest one. If the input is part of data from the Insight Journal, just point us to it.


(wonny) #3

I have 5 DVFs files but they are total 2.6 GB and I don’t know how I put it here.


(Dženan Zukić) #4

You should compress the files into .7z, .bz2, .xz or other similar archive format using maximum compression ratio option. Then upload it to a file sharing service, one of which is files.fm. Then paste a link to it here.


(wonny) #6

(Dženan Zukić) #7

Your initialization of variables ncontrol and closedim was faulty - only the first element was initialized. Here is a slightly updated code, which crashes in splineFilter->Update(); due to insufficient memory (9 GB private bytes on my machine).

tester.cpp (4.7 KB)

The word scattered in the filter’s name BSplineScatteredDataPointSetToImageFilter implies that it is not really meant to work with densely sampled data. If it also crashes for you, you could try downsampling the input images to a lower resolution.

I don’t know if there is a more appropriate filter for what you want to do.


(wonny) #8

Hi, Zukic

Thank you for your help. I tried your updated code. it came out error, and can you see it?


(Dženan Zukić) #9

I think your input fields contain infinities. In Modules\Filtering\ImageGrid\include\itkBSplineScatteredDataPointSetToImageFilter.hxx, line 555, for n=155058176 the variable point has content {-inf,-inf,-inf,3.0}. But you could have discovered this yourself have you tried to debug your program.