function ccurve2 = TaubinSmooth( ccurve, LAMBDA, ITER)
%function xy = curveSmooth( ccurve )
%function xy = curveSmooth( ccurve, <LAMBDA>, <KB>, <ITER>, <NUM_NEIGHBOR_RINGS> )
% 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
