smoothing 2D scattered points

Hi Friends,

Is there anyway to smooth 2D scattered points?

I drew my question in this picture:

on the left we have an scattered path of points, but, rationally one can think of a smoothed representing set of points (black dots) in the right picture.

Thank you very much.

Answer 1

[Image Analyst]

Yes. You can use something like spline or sgolayfilt() to smooth both the x and y coordinates independently. See attached demo. Here is the main part:

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

Attach your x and y data if you can't figure it out from the well commented demo code.

[Alireza Ahani]

I used the lines, the result is as follow:

some parts the right wing of data is what i'm looking for, but the rest is not.

I attched the data.

[Image Analyst]

This is what your data looks like when plotted:

Plus if you put lines in between the dots, you'll see it's not parameterized so that the dots follow the shape. What I would suggest is to find the dots above and below 363 and fit lines through them

s = load('xy.mat')
x = s.X;
y = s.Y;
plot(x, y, 'b.');
grid on;

% Find line through top data.
above363 = y > 363;
x2 = x(above363);
y2 = y(above363);
% Fit a line through them
coefficients2 = polyfit(x2, y2, 1)
% Get the fitted values
x2Fit = unique(x);
y2Fit = polyval(coefficients2, x2Fit);
hold on;
plot(x2Fit, y2Fit, 'r-', 'LineWidth', 2);

% Find line through bottom data.
below363 = y <= 363;
x3 = x(below363);
y3 = y(below363);
% Fit a line through them
coefficients3 = polyfit(x3, y3, 1)
% Get the fitted values
x3Fit = unique(x3);
y3Fit = polyval(coefficients3, x3Fit);
hold on;
plot(x3Fit, y3Fit, 'r-', 'LineWidth', 2);

Is this from an image? If so then you should use bwskel() on the binary image instead.

[Alireza Ahani]

Thank you very much,

But I don't have the pre-knowledge about this place (y==363) or any other turning points.

[Alireza Ahani]

oh, I saw the "bwskel" you mention, that is so much interesting,

Highly likely it is the solution for this problem.

Thank you very much @Image Analyst . God bless.

Answer 2

[Craig Stewart]

One approach would be to do a sequential search on neighboring points, averaging these - demo code below. This could be extended to also limit maximum point-to-point step size

function [xs,ys] = points2line %(x,y,catchmentRad)

doplot = true;
x = [-35.2257 -24.2666 -24.2666 -15.9168 -14.3512 -11.22 -7.0451 -7.0451 -6.0014 6.0014 8.6107 12.7856 21.6573 22.701 30.5289 38.3568 46.1848 47.7504 47.7504 46.7066 43.0536 39.4006 39.4006 39.4006 39.4006 38.8787 33.6601 35.2257 35.2257 24.2666 31.0508 34.1819 35.2257 38.3568 36.2694 35.7475 34.1819 32.0945 32.0945 32.6164 33.6601 32.6164 32.6164 32.6164 35.2257 48.2722 51.9252 58.1876 74.8872 74.8872 76.4528 78.0183 79.5839 82.1932 83.237 85.8463 85.8463 86.3681 86.3681 82.1932 77.4965 75.9309 75.9309 63.9281 63.9281 55.5783 55.5783 55.5783 44.0973 44.0973 33.1382 33.1382 30.5289 25.8322 25.3103 21.6573 19.048 19.048 14.8731 11.22 2.8702 2.8702 0.7828 0.7828 0.7828 -4.9577 -6.5233 -7.0451 -10.1763 -14.3512 -14.3512 -18.0042 -24.2666 -24.2666 -28.4415 -28.4415 -47.2285 -47.7504 -53.4908 -58.1876 -58.7094 -62.8843 -64.9718 -65.4937 -64.4499 -64.4499 -64.4499 -64.4499 -64.9718 -65.4937 -63.4062 -63.4062 -59.7532 -59.7532 -61.3188 -60.275 -60.7969 -62.3625 -63.9281 -63.9281 -63.9281 -69.1467 -70.1904 -70.1904 -70.7123 -70.7123 -73.3216 -75.409 -76.4528 -76.4528 -83.7588 -83.7588 -81.6714 -81.6714 -81.6714 -81.1495 -82.1932 -82.1932 -82.1932 -82.1932 -63.9281 -64.9718 -64.4499 -63.9281 -62.8843 -62.8843 -61.8406 -61.8406 34.1819 37.3131 38.3568 33.6601 33.6601 34.7038 33.1382 33.1382 33.1382 33.1382 33.1382 34.7038 34.7038 32.6164 32.6164 35.7475 35.7475 36.2694 38.8787 39.4006 39.9224 39.9224 39.9224 39.9224];
y = [-647.8209 -651.4739 -651.4739 -651.4739 -650.952 -650.4302 -652.5176 -652.5176 -653.5614 -654.0832 -651.9958 -649.9083 -650.952 -651.9958 -653.0395 -651.9958 -651.4739 -650.4302 -649.9083 -645.7334 -639.4711 -630.0776 -627.9901 -623.8152 -623.2934 -622.2496 -606.5938 -599.8096 -599.8096 -581.5444 -575.2821 -570.5853 -567.976 -560.1481 -552.3202 -551.2764 -548.1453 -541.3611 -540.8392 -538.7518 -528.8364 -517.3554 -516.8336 -514.2243 -502.2214 -503.787 -503.787 -502.2214 -496.481 -496.481 -487.0874 -486.0437 -480.8251 -475.0846 -472.9972 -462.0381 -461.5162 -452.6446 -452.1227 -444.8166 -432.292 -429.6827 -429.6827 -415.0705 -415.0705 -406.1989 -406.1989 -406.1989 -396.2835 -396.2835 -392.6305 -392.6305 -395.2398 -396.2835 -396.8054 -399.4147 -402.5458 -402.5458 -405.677 -407.2426 -406.1989 -406.1989 -407.7645 -407.7645 -407.7645 -414.0268 -414.0268 -414.0268 -414.0268 -416.1142 -416.1142 -416.1142 -419.2454 -419.2454 -420.811 -420.811 -426.5515 -426.5515 -432.8138 -436.4669 -436.9887 -442.2073 -453.1664 -454.2102 -457.8632 -457.8632 -462.5599 -462.5599 -464.6474 -468.8223 -477.1721 -477.1721 -482.3907 -482.3907 -498.5684 -499.6121 -514.2243 -518.3992 -523.0959 -524.1396 -528.3145 -539.7955 -541.8829 -542.9267 -548.6671 -550.2327 -555.9732 -557.5388 -563.8011 -564.8449 -572.6728 -574.2384 -583.6319 -591.4598 -595.1128 -599.2877 -605.0282 -608.1594 -609.725 -610.7687 -531.4457 -526.7489 -523.0959 -518.3992 -515.268 -511.0931 -503.787 -502.2214 -537.1862 -546.0578 -550.2327 -559.1044 -560.1481 -561.1918 -566.4104 -573.7165 -573.7165 -573.7165 -582.0663 -586.763 -589.8942 -597.2003 -597.7221 -604.5063 -605.0282 -605.0282 -610.7687 -612.3343 -620.1622 -621.2059 -621.2059 -621.2059];
catchmentRad = 15;

if doplot
    figure
    plot(x,y,'k.','MarkerSize',20)
    hold on
end

% Trace the outline by finding next closest point (in chunks)
xy = x + 1i*y;
ii = 1;
jj = 0;
used = false(size(xy));
while any(~used)
    jj = jj + 1;
    d = abs(xy-xy(ii));
    inrange = d <= catchmentRad & ~used; %   
    XY(jj) = mean(xy(inrange)); %#ok<AGROW>
    used(inrange) = true;
    if doplot
        if jj == 1
            h1 = plot(x(inrange),y(inrange),'g.','MarkerSize',10,'DisplayName','In range');
            h2 = plot(x(ii),y(ii),'r.','MarkerSize',20,'DisplayName','Point');
            h3 = plot(real(XY(jj)),imag(XY(jj)),'b.','MarkerSize',20,'DisplayName','Mean');
        else
            set(h1,'XData',x(inrange),'YData',y(inrange))
            set(h2,'XData',x(ii),'YData',y(ii))
            set(h3,'XData',real(XY(jj)),'YData',imag(XY(jj)))
        end
        pause(0.05)
    end
    
    % find next closest point
    [~,ii] = min(abs((xy+1e6*double(used)-xy(ii))));
end
X = real(XY);
Y = imag(XY);

% Interpolate to smooth
D = cumsum([0, abs(diff(XY))]);
d = linspace(0,max(D),1000);
xs = interp1(D,X,d,'spline');
ys = interp1(D,Y,d,'spline');

if doplot
    delete([h1 h2 h3])
    %plot(X,Y,'g.','MarkerSize',2)
    plot(xs,ys,'b','LineWidth',2)
end

[Alessio Cislaghi]

Dear Craig, your script is very useful. Is the methodology explained in a document? What is catchmentRad?

[Image Analyst]

It's the "catchment radius", or in plain English the maximum distance away that another point is allowed to be and still be considered part of the same line. If it's farther away than that, it's considered too far away to be the next point.

Answer 3

[KSSV]

Read about polyfit.

[Alireza Ahani]

Hi,

an injective function can not be fitted into this type of data (X1--->Y1 && X1--->Y2).

smoothing seems to be the only remedy.