[Insight-users] itkBSplineScatteredDataImageFilter exception thrown

[Ian Malone]

Nicholas Tustison wrote:
> 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?
>   

Hi, I've attached a picture which shows the section being input (a short 
curve section on the sphere, blue) and the resulting spline (green), the 
red line is a continuation of the curve on the surface, together with 
the sphere to get an idea of the relative positions. For a point in the 
middle of the curve:
t=0.5
original [0.825336, 0.52007, 0.219882] magnitude 1
fitted (evaluating spline at t=0.5)   [0.928681, 0.585141, 0.247394] 
magnitude 1.12519

This is purely fitting to the algebraic function with 100 input points 
evenly spaced over the blue arc. I take your point about the control 
points not lying on the sphere, but the change in x value is odd when it 
is constant across the input points. Possibly this is due to some kind 
of ringing from the filter background value? SetNumberOfLevels(3) 
results in being very close to the surface.

Ian

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20101215/5b564079/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: itk-curvebsplinetest2.png
Type: image/png
Size: 44761 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20101215/5b564079/attachment.png>