# DVFs interpolation by BSpline

**URL:** https://discourse.itk.org/t/dvfs-interpolation-by-bspline/1325
**Category:** Beginner Questions
**Created:** [October 3, 2018, 6:26pm UTC](https://discourse.itk.org/t/dvfs-interpolation-by-bspline/1325 "2018-10-03T18:26:40Z")
**Posts on this page:** 8
**Page:** 1

<div class="post-metadata">

### Author: ![kjwjj627](https://discourse.itk.org/letter_avatar_proxy/v4/letter/k/bcef8e/32.png) [@kjwjj627](https://discourse.itk.org/u/kjwjj627)
#### Post date: [October 3, 2018, 6:26pm UTC](https://discourse.itk.org/t/dvfs-interpolation-by-bspline/1325/1 "2018-10-03T18:26:40Z")

</div>

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;  
};

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [October 3, 2018, 6:29pm UTC](https://discourse.itk.org/t/dvfs-interpolation-by-bspline/1325/2 "2018-10-03T18:29:34Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![kjwjj627](https://discourse.itk.org/letter_avatar_proxy/v4/letter/k/bcef8e/32.png) [@kjwjj627](https://discourse.itk.org/u/kjwjj627)
#### Post date: [October 3, 2018, 7:21pm UTC](https://discourse.itk.org/t/dvfs-interpolation-by-bspline/1325/3 "2018-10-03T19:21:49Z")

</div>

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

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [October 3, 2018, 7:34pm UTC](https://discourse.itk.org/t/dvfs-interpolation-by-bspline/1325/4 "2018-10-03T19:34:33Z")

</div>

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](https://files.fm/). Then paste a link to it here.

---

<div class="post-metadata">

### Author: ![kjwjj627](https://discourse.itk.org/letter_avatar_proxy/v4/letter/k/bcef8e/32.png) [@kjwjj627](https://discourse.itk.org/u/kjwjj627)
#### Post date: [October 4, 2018, 2:03pm UTC](https://discourse.itk.org/t/dvfs-interpolation-by-bspline/1325/6 "2018-10-04T14:03:30Z")

</div>

> **[File upload, sharing and cloud backup online service.](https://files.fm/u/635ax46c)**
>
> Visit this link to see the content

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [October 4, 2018, 6:50pm UTC](https://discourse.itk.org/t/dvfs-interpolation-by-bspline/1325/7 "2018-10-04T18:50:03Z")

</div>

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](https://discourse.itk.org/uploads/default/original/1X/904cf6fc3a35103c29b878dbafbddc0879032508.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](https://itk.org/ITKExamples/src/Filtering/ImageGrid/ResampleAVectorImage/Documentation.html) the input images to a lower resolution.

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

---

<div class="post-metadata">

### Author: ![kjwjj627](https://discourse.itk.org/letter_avatar_proxy/v4/letter/k/bcef8e/32.png) [@kjwjj627](https://discourse.itk.org/u/kjwjj627)
#### Post date: [October 5, 2018, 2:10pm UTC](https://discourse.itk.org/t/dvfs-interpolation-by-bspline/1325/8 "2018-10-05T14:10:37Z")

</div>

Hi, Zukic

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

 ![error](https://discourse.itk.org/uploads/default/original/1X/c16a0631c6af5dcd2903ff4b6ea356d3e7383b9a.jpeg)

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [October 8, 2018, 3:59pm UTC](https://discourse.itk.org/t/dvfs-interpolation-by-bspline/1325/9 "2018-10-08T15:59:31Z")

</div>

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.