[ITK] [ITK-users] DVFs interpolation by BSpline

[wonnybro]

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, unfortunately it crashes.

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



--
Sent from: http://itk-users.7.n7.nabble.com/
The ITK community is transitioning from this mailing list to discourse.itk.org. Please join us there!
________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
https://itk.org/mailman/listinfo/insight-users