Code covered by the BSD License

# Smoothing 2D Contours Using Local Regression Lines

### Tolga Birdal (view profile)

17 Mar 2011 (Updated )

This submission contains both the implementation and the test function for 2D contour smoothing.

test_smooth_contours()
```% by Tolga Birdal
%{
This file contains both the implementation and the test function for
contour smoothing. So to test, just press F5.
The contour smoothing is done by projecting all the contour points
onto the local regression line. The radius defines the number of points
on each side of the contour, which will contribute to the computation
of the local regression line. The higher the number of points, the
smoother the curve.

Because of the linear nature of fitting, when too much smoothing is
desired,  the algorithm loses important features such as corners, and
confuses in such critical regions.

Actually, the algorithm is fully parallelizable. If a parallel
implemtation is desired the main loop in smooth_contours function can
easily be parallelized.

You can check further details in the comments.
%}
function test_smooth_contours()

close all;

% to binary
L=logical(I);

% get the start point for bwtraceboundary
[y x]=find(L==1);

% get the contour
contour = bwtraceboundary(L, [47,187], 'E', 8);

% smooth the contour
[Ys Xs]=smooth_contours(contour(:,1), contour(:,2), 21);

out=zeros(size(I));
ind=sub2ind(size(I),round(Ys), round(Xs));
out(ind)=1;
figure,imshow(out);

% show both
figure, subplot(1,2,1), scatter(x,y);
subplot(1,2,2), scatter(Xs, Ys);

end

% The actual computation

Xs=zeros(length(X),1);
Ys=zeros(length(X),1);

% copy out-of-bound points as they are

% obtain the bounding box
maxX=max(max(X));
minX=min(min(X));
maxY=max(max(Y));
minY=min(min(Y));

% smooth now
xLocal=X(ind);
yLocal=Y(ind);

% local regression line
%p=polyfit(xLocal,yLocal,1);
[a b c] = wols(xLocal,yLocal,gausswin(length(xLocal),5));
p(1)=-a/b;
p(2)=-c/b;

% project point on local regression line
[x2, y2]=project_point_on_line(p(1), p(2), X(i), Y(i));

% check erronous smoothing
% points should stay inside the bounding box
if (x2>=minX && y2>minY && x2<=maxX && y2<=maxY)
Xs(i)=x2;
Ys(i)=y2;
else
Xs(i)=X(i);
Ys(i)=Y(i);
end
end

end

% Projects the point (x1, y1) onto the line defined as y=m1*x+b1
function [x2, y2]=project_point_on_line(m1, b1, x1, y1)

m2=-1./m1;
b2=-m2*x1+y1;
x2=(b2-b1)./(m1-m2);
y2=m2.*x2+b2;

end

function [a b c] = wols(x,y,w)
% Weighted orthogonal least squares fit of line a*x+b*y+c=0 to a set of 2D points with coordiantes given by x and y and weights w
n = sum(w);
meanx = sum(w.*x)/n;
meany = sum(w.*y)/n;
x = x - meanx;
y = y - meany;
y2x2 = sum(w.*(y.^2 - x.^2));
xy = sum(w.*x.*y);
alpha = 0.5 * acot(0.5 * y2x2 / xy) + pi/2*(y2x2 > 0);
%if y2x2 > 0, alpha = alpha + pi/2; end
a = sin(alpha);
b = cos(alpha);
c = -(a*meanx + b*meany);
end```