[Insight-users] itkBSplineScatteredDataImageFilter exception thrown

[Nicholas Tustison]

Thanks for that.  That's weird. Let me try to run your code and I'll let you know what I find.


On Dec 15, 2010, at 12:25 PM, Ian Malone wrote:

> 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
>>>     
>> 
>>   
> 
> <itk-curvebsplinetest2.png>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20101215/ed1a4d08/attachment.htm>