Search code examples
matlabinterpolationcurve-fittingnormals

Matlab generate smooth curve between scatter points


I need to generate a curve between scatter points then identify the unit normal of the curve at each point. Here is an example of a point cloud

figure
x = [1 2 0 0 1 1 2 3 4 2];
y = [4 6 9 1 1 2 4 9 2 3];
scatter(x,y)
hold on
line(x,y)
xlim([0 4])
ylim([0 10])

enter image description here

NOTE: the 2 points along the y-axis are connected

Instead of a line between the points, I'd like to create a smooth curve. I'm not sure how to do this when points in x and y repeat. An attempt using spline fails. After I know the curve, I need to find the unit normals at each point. How do I go about this?

EDIT: Basically I want to do what is show here for polyfit in the matlab docs. Assuming that x was unique in my case, this wouldn't be an issue. I could identify the polynomial and then, I believe, determine the unit normals from the polynomial function evaluated at that point. But in my case, the x and y data repeat so a straight forward application doesn't work.


Solution

  • One way to get a smooth path is to treat this as a parametric function and interpolate x and y separately.

    x = [1 2 0 0 1 1 2 3 4 2];
    y = [4 6 9 1 1 2 4 9 2 3];
    t = 1:numel(x);
    
    tq = 1:0.1:t(end);
    xq = interp1(t,x,tq,'v5cubic');
    yq = interp1(t,y,tq,'v5cubic');
    
    plot(x,y,' ob',xq,yq,'-r');
    

    To estimate the normals you can take the average normal of the two line segments around the sample points. This code is a bit ugly but it gets the job done.

    n = zeros(2,numel(x));
    for tidx = 1:numel(t)
        tt = t(tidx);
        idx1 = find(tq <= tt,1,'last');
        idx0 = idx1 - 1;
        idx2 = idx1 + 1;
        if idx0 > 0
            n1 = [yq(idx1) - yq(idx0); xq(idx0) - xq(idx1)];
            n(:,tidx) = n(:,tidx) + n1/norm(n1);
        end
        if idx2 <= numel(tq)
            n2 = [yq(idx2) - yq(idx1); xq(idx1) - xq(idx2)];
            n(:,tidx) = n(:,tidx) + n2/norm(n2);
        end
        n(:,tidx) = n(:,tidx) / norm(n(:,tidx));
    end
    
    plot(x,y,' ob',xq,yq,'-r',[x.' x.'+n(1,:).'].', [y.' y.'+n(2,:).'].',' -k');
    axis equal;
    

    enter image description here

    If you use pchip instead of v5cubic for the interpolation method then you get more symmetry around the sample points. However, it appears that any sharp turns (90 degrees or greater) are not smoothed.

    enter image description here