[Insight-users] itkBSplineScatteredDataImageFilter exception thrown

[Nicholas Tustison]

For B-splines you shouldn't expect to see the control points on the sphere and depending on how noisy you're data is, I don't know that you should expect it on the sphere either.  Are these scattered points actually on the surface?  

Also, I'm having difficulty understanding the offset issue.  Can you send me a picture?



> Thanks, I've tried making that fix and also just copying the .txx and .h out of that git tree-ish, the crash still occurs with 1000 points and 0.001 spacing, but I assume that's simply because that range doesn't cover [0,1] so shouldn't be expected to work.
> 
> What's currently confusing me is that I'm trying to fit curves (1D parameter, R^3 values) on the surface of a radius 1 sphere at the moment and the control points and spline positions I'm getting out are not on the sphere at level 1, higher levels improve things, but using the same number of control points at level 1 does not. I assume this is simply related to the fitting and I'm not too concerned about it. The really puzzling thing is this curve at constant theta (code attached):
> pointval[0] = cos(theta);
> pointval[1] = sin(theta)*cos(phi);
> pointval[2] = sin(theta)*sin(phi);
> 
> While pointval[0] is constant (0.825336 in my example), when I fit at level 1 with 20 control points I get a constant offset along the middle of the curve (with a value of 0.928681, changes slightly if I add more input points). This leaves me wondering if I'm doing something wrong.
> 
> imalone
> 
> 
> 
> #include "itkBSplineScatteredDataPointSetToImageFilter.h"
> #include "itkPointSet.h"
> #include <math.h>
> 
> const unsigned VectorDim=3;
> 
> typedef itk::Vector < double , VectorDim > VectorType ;
> typedef itk::Image < VectorType , 1 > SplineCurveImageType ;
> typedef itk::PointSet < VectorType , 1 > SplineCurvePointSetType ;
> typedef itk::BSplineScatteredDataPointSetToImageFilter
>    < SplineCurvePointSetType , SplineCurveImageType > SplineFilterType ;
> 
> 
> VectorType getpointattt(double tt){
>  VectorType pointval;
>  double phistart=0.2;
>  double phirange=0.4;
>  double phi=phistart+tt*phirange;
>  double theta=0.6;
>  pointval[0] = cos(theta);
>  pointval[1] = sin(theta)*cos(phi);
>  pointval[2] = sin(theta)*sin(phi);
>  return pointval;
> }
> 
> int main (void) {
>  SplineFilterType::Pointer spline = SplineFilterType::New();
>  SplineCurvePointSetType::Pointer datapoints = SplineCurvePointSetType::New();
> 
>  SplineCurveImageType::SpacingType spacing ;
>  SplineCurveImageType::SizeType size ;
>  SplineCurveImageType::PointType origin ;
> 
>  unsigned numberofinputpoints = 100;
> 
>  spacing.Fill(1.0/(1000));
>  size.Fill ( 1000+1 );
>  origin.Fill ( 0.0 );
>  double pointspacing = 1.0/(numberofinputpoints-1);
> 
> 
>  double tstep=1.0/(numberofinputpoints-1);
>  for (unsigned ii=0 ; ii<numberofinputpoints ; ii++ ) {
>    VectorType midpointval;
> 
>    double tt=tstep*ii;
>    midpointval=getpointattt(tt);
> 
>    double mag=0;
>    for ( unsigned element=0 ; element<3 ; element++) {
>      mag += midpointval[element]*midpointval[element];
>    }
>    mag=sqrt(mag);
> 
>    SplineCurvePointSetType::PointType thispointparam;
>    thispointparam[0] = tt;
>    datapoints->SetPoint(ii, thispointparam);
>    datapoints->SetPointData(ii,midpointval);
>  }
> 
>  SplineFilterType::ArrayType ncontrol ;
>  ncontrol[0]=20;
>  SplineFilterType::ArrayType closedim;
>  closedim[0]= 0;
>  int splineorder=2;
>  spline->SetGenerateOutputImage(false);
>  spline->SetOrigin ( origin );
>  spline->SetSize ( size );
>  spline->SetSpacing ( spacing );
>  spline->SetInput ( datapoints );
>  spline->SetSplineOrder ( splineorder );
>  spline->SetNumberOfControlPoints ( ncontrol );
>  spline->SetNumberOfLevels(1);
>  spline->SetCloseDimension ( closedim );
>  try {
>    spline->Update();
>  }
>  catch (itk::ExceptionObject &err) {
>    std::cerr << "Failure in spline fitting" << std::endl;
>    std::cerr << spline << std::endl;
>    std::cerr << err << std::endl;
>    return 0 ;
>  }
> 
>  for (unsigned ii=0 ; ii<=20; ii++) {
>    std::cout<<"Point "<< ii << std::endl;
> 
>    double tt=ii/20.0;
>    VectorType loc = getpointattt(tt);
>    VectorType splineloc;
>    SplineCurvePointSetType::PointType ttpoint;
>    ttpoint[0] = tt;
> 
>    spline->Evaluate(ttpoint, splineloc);
>    double splinemag=0, mag=0;
>    for (unsigned el=0; el<3; el++ ) {
>      splinemag+=splineloc[el]*splineloc[el];
>      mag+=loc[el]*loc[el];
>    }
>    splinemag=sqrt(splinemag);
>    mag=sqrt(mag);
>    std::cout << tt << " "
> 	      << loc
> 	      << " mag " << mag << " "
> 	      << splineloc
> 	      << " mag " << splinemag
> 	      << std::endl;
>  }
> 
>  std::cout << "Control points " << spline->GetCurrentNumberOfControlPoints()
> 	    << std::endl;
> 
>  int showcontrolpoints=1;
>  if (showcontrolpoints) {
>    SplineCurveImageType::Pointer phiLattice = spline->GetPhiLattice();
>    SplineCurveImageType::RegionType::SizeType latticesize = phiLattice->GetLargestPossibleRegion().GetSize();
>    for ( unsigned ii=0; ii< latticesize[0] ; ii++) {
>      VectorType controlloc;
>      SplineCurveImageType::IndexType index;
>      index[0]=ii;
>      controlloc=phiLattice->GetPixel(index);
>      std::cout << controlloc << std::endl;
>    }
>  }
> 
>  return 0;
> };
> _____________________________________
> 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.html
> 
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
> 
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users