<html><head></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">Thanks for that. That's weird. Let me try to run your code and I'll let you know what I find.<div><br></div><div><br><div><div>On Dec 15, 2010, at 12:25 PM, Ian Malone wrote:</div><br class="Apple-interchange-newline"><blockquote type="cite">
<div bgcolor="#ffffff" text="#000000">
Nicholas Tustison wrote:
<blockquote cite="mid:61580405-01B7-44B8-A63E-0B1E95E3E975@gmail.com" type="cite">
<pre wrap="">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?
</pre>
</blockquote>
<br>
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:<br>
t=0.5<br>
original [0.825336, 0.52007, 0.219882] magnitude 1<br>
fitted (evaluating spline at t=0.5) [0.928681, 0.585141, 0.247394]
magnitude 1.12519<br>
<br>
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.<br>
<br>
Ian<br>
<br>
<blockquote cite="mid:61580405-01B7-44B8-A63E-0B1E95E3E975@gmail.com" type="cite">
<pre wrap="">
</pre>
<blockquote type="cite">
<pre wrap="">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 <a class="moz-txt-link-abbreviated" href="http://www.kitware.com/">www.kitware.com</a>
Visit other Kitware open-source projects at
<a class="moz-txt-link-freetext" href="http://www.kitware.com/opensource/opensource.html">http://www.kitware.com/opensource/opensource.html</a>
Kitware offers ITK Training Courses, for more information visit:
<a class="moz-txt-link-freetext" href="http://www.kitware.com/products/protraining.html">http://www.kitware.com/products/protraining.html</a>
Please keep messages on-topic and check the ITK FAQ at:
<a class="moz-txt-link-freetext" href="http://www.itk.org/Wiki/ITK_FAQ">http://www.itk.org/Wiki/ITK_FAQ</a>
Follow this link to subscribe/unsubscribe:
<a class="moz-txt-link-freetext" href="http://www.itk.org/mailman/listinfo/insight-users">http://www.itk.org/mailman/listinfo/insight-users</a>
</pre>
</blockquote>
<pre wrap=""><!---->
</pre>
</blockquote>
<br>
</div>
<span><itk-curvebsplinetest2.png></span></blockquote></div><br></div></body></html>