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

[Dženan Zukić]

Hi JongWon,

welcome to ITK! To reproduce your problem we will need the input file too.
If it crashes with different input files, please provide the smallest one.
Also, we have migrated to the forum <https://discourse.itk.org/>, so please
post the updated question there.

Regards,
Dženan

On Wed, Oct 3, 2018 at 11:30 AM wonnybro <kjwjj627 at gmail.com> wrote:

> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://itk.org/pipermail/community/attachments/20181003/bc5face6/attachment.html>
-------------- next part --------------
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