Fitting smooth closed spline to scattered 2D points

Hi everyone,

I have a set of points scattered in a 2D plane. The points are scattered in a somwhat circular/elipptical way, but I would like to fit some arbitrarily (but smoothed) shaped spline. The points are distributed like:

To enable closing of a spline, the last and first data points are the same.

As said I would like to fit a closed smoothed spline to this curve. I've tried the cscvn-option giving me the below, but I would like something a bit more smoothened out (cscvn forces the spline to go exactly through each and every point):

Reading up in the help, csaps seem to enable smoothing (to be fair, I'm not really sure how), but if trying this on my data I get a smoothed non-closed curve like below:

The last option I've looked at is creating a new parameter to which I fit both my (x,y)-data (as suggested in an earlier thread: http://www.mathworks.com/matlabcentral/newsreader/view_thread/33868 ), but this actually gives me an even more jagged curve:

However I might be applying the parametrization in an incorrect way...

So, to my question. Does anyone have a good way of creating a closed smoothed spline to my data? Any suggestion or help would be very much appreciated. (P.S: I'm trying to create an automatized code, so I'm preferably looking for executable code than GUI in the curve fitting toolbox (but I understand that they are somewhat connected...)).

Thanks for the help!

/David

Answer 1

[Image Analyst]

If you don't want or require the curve to go through every point exactly, you can use a Savitzky-Golay filter on the parameterized data.

% Now smooth with a Savitzky-Golay sliding polynomial filter
windowWidth = 35
polynomialOrder = 2
smoothX = sgolayfilt(x, polynomialOrder, windowWidth);
smoothY = sgolayfilt(y, polynomialOrder, windowWidth);

[David]

Hi,

thanks for the suggestion. Applying such a filter on my parameterized (cscvn-generated) curve I end up with a smoother spline:

However I seem to always get a non-continuity at the end even though the first and last entry of x and y are given as identical (I've tested for higher order polynomial filters but still get this discontinuity).

Is there a clever way of getting a closed smoothed spline, maybe using the Savitzky-Golay filter?

Thanks again for the help.

/David

[Image Analyst]

Did you try making the first and last point the same? Or maybe even replicate more than that, like half a window width's worth of points? It would make it easier for people to help you if you gave the x and y arrays.

[David]

Of course: the x and y data in the above case is:

x = [   70.7137
 72.5285
 73.5333
 75.2102
 76.1013
 74.0646
 70.6159
 65.8768
 62.0974
 57.4835
 52.0009
 44.6378
 44.0573
 40.1012
 42.1993
 43.7036
 43.4872
 46.6293
 50.1801
 54.1760
 58.7491
 63.6772
 70.7137];

y = [   68.4054
 72.5215
 76.7557
 80.9434
 85.8935
 90.3746
 93.9230
 95.2436
 97.9917
 97.1133
100.3399
 92.8572
 87.8427
 84.5356
 79.8470
 76.0016
 71.3394
 68.1848
 65.6478
 63.6865
 61.0153
 62.0705
 68.4054];

In which the first and last data point is the same (x(1) = x(end)). As I'm not that familiar with the filter I'm perhaps not really getting the window width. What does it represent? If understanding I'm hopping to adapth the replication the way you suggest...

By the way, the smoothed and (non-closed) curve that I'm now getting is obtained simply by:

    data(:,1) = x;
    data(:,2) = y;
    [spline_points, t] = fnplt(cscvn(data),'r',2);
    windowWidth = 45;
    polynomialOrder = 3;
    spline_points(1,:) = sgolayfilt(spline_points(1,:), polynomialOrder, windowWidth);
    spline_points(2,:) = sgolayfilt(spline_points(2,:), polynomialOrder, windowWidth);

If you have any conrete example on how to treat the data, that'd be great!

[Image Analyst]

The window width is the number of points to be used when fitting a polynomial. For example with 45 points, and order 3, it will fir a cubic function to the points 22 before, the point itself, and the 22 points after. Then it will replace the point with the value of the cubic at that location. For your plot with so few points, you might try a window width of 3 or 5. The larger the window, the smoother the results is and the smaller the window, the more it hugs the training points. The larger the order of the polynomial to use, the closer it will hug the training points (of course), and the lower the order, the smoother the fitted shape will be.

[Image Analyst]

Another option to consider: http://www.mathworks.com/matlabcentral/answers/103916-how-can-i-smooth-a-signal-while-preserving-peaks-in-matlab-r2013a

[David]

Hmmmm, I see. Playing around with the WindowWidth and PolynomialOrder I can find some optimal fitting for this specific dataset. But since I have multiple sets of similar shape, I'm still looking at finding some automatized optimization.

So I was more thinking of perhaps enforcing continuity in the derivative, but that probably belongs under a new question topic.

[Image Analyst]

David you can use activecontour() - do you know what active contours (also know as snakes or balloons or alpha shapes, or at least related to them) are? Attached is an example of active contour. You can use poly2mask() to turn your points into a solid mask, then use activecontour to create a smooth boundary. This uses an energy minimization method that perhaps you'll feel more comfortable with than polynomial regression (which is basically what sgolayfilt() is).

[David]

Hi again, thanks for coming with new suggestions.

No I'm not that familiar with active contours. I did some tries but seem not to find the right way of getting a smoother boundary. According to the help of activecontour(), an example command would be:

bw = activecontour(I, mask, 200, 'edge','SmoothFactor',1.5);

but when I try to perform the same (see code below) I get an error prompt:

if true
BW = poly2mask(xy(1,:),xy(2,:),120,120);
mask = zeros(size(BW));
mask(10:end-10,10:end-10) = 1;
bw = activecontour(BW,mask,'SmoothFactor',1.5);

Error message:

Error using activecontour>parse_inputs (line 287)
Too many input arguments were passed to activecontour.

So it seems like 'SmoothFactor' cannot be passed to the activecontour. Has this been changed in a new Matlab version (I'm running R2014a)? Or am I running the commands incorrectly?

Lastly: if I manage to get a smoothed contour map (i.e. binary image), can I then retrieve the edge points (i.e. going back to some spline/polygon shape) again?

Again, thanks for all the help.

Answer 2

[Praful Agrawal]

the original task can be achieved by converting the piecewise polynomial fit to bspline fit using fn2fm function:

X - set of points

pp = cscvn(X); bsp = fn2fm(pp, 'B-');