function ccurve2 = TaubinSmooth( ccurve, LAMBDA, ITER)
%function xy = curveSmooth( ccurve )
%function xy = curveSmooth( ccurve, , , , )
% Smooths the given curve
% Precondition: curveXYPair is already ordered (i.e. clockwise or
% counter-clockwise) from within 'ccurve'; this should be the case
% if you used the ClosedCurve constructor, and didn't modify curveXYpair.
% @param ccurve ClosedCurve object to smooth
% @param LAMBDA parameters to set to control smoothing, as
% described in the article below
% @param ITER ""
% @returns ccurve2 ClosedCurve object smoothed
% Reference:
% Taubin, Gabriel. "A Signal Processing Approach to Fair Surface
% Design." IBM T.J. Watson Research Center.
%
% Taubin also has a website online with all sorts of articles:
% (BOOT INTO LINUX AND COPY FROM BOOKMARKS).
% (earlier, linux was crashing like no other...)
% We're using Taubin's mesh approach, only to a 2D curve...
% at least we know there's a 3d extension as well!
% (c) Moo K. Chung. 2010.
% Modified from Tomas Hoffmann's original code.
xyPair = ccurve.curveXYpair;
NUM_NEIGHBOR_RINGS = 1;
numPairs = length( xyPair(:,1) );
I = eye(numPairs);
W = zeros( numPairs );
for i=1:numPairs
% get the indices of the neighbors
neighbors = [(i-NUM_NEIGHBOR_RINGS):(i-1) (i+1):(i+NUM_NEIGHBOR_RINGS)];
neighbors = neighborFix( neighbors, numPairs );
% get the weights
for j=1:length(neighbors)
W(i,neighbors(j)) = phi(i,neighbors(j),xyPair);
end
W(i,:) = W(i,:) ./ sum(W(i,:));
end
% calculate K
K = I - W;
fk = (I - LAMBDA*K); % Gaussian filter
%fk = (I - LAMBDA*K)*(I - MU*K);
xyPair = fk^ITER * xyPair;
% return a new ccurve
ccurve2 = xyPair;
return; % END FUNCTION
function p = phi( i, j, xyPair )
%function p = phi( i, j )
% returns phi(v_i,v_j)
p = sqrt( (xyPair(i,1)-xyPair(j,1))^2 + (xyPair(i,2)-xyPair(j,2))^2 );
return; % END FUNCTION
function neighbor = neighborFix(neighbor, numPairs)
%function neighbor = neighborFix(neighbor, numPairs)
% fixes neighbors to the range (1,numPairs)
for n=1:length(neighbor)
neighbor(n) = modPlus( neighbor(n), numPairs );
end
return; % END FUNCTION
function n = modPlus(x, m)
%function n = modPlus(x, m)
n = mod(x-1,m) + 1;
return % END FUNCTION