smoothsurf.m

function p = smoothsurf(node, mask, conn, iter, useralpha, usermethod, userbeta)
%
% p=smoothsurf(node,mask,conn,iter,useralpha,usermethod,userbeta)
%
% smoothing a surface mesh
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2007/11/21
%
% input:
%    node:  node coordinates of a surface mesh
%    mask:  flag whether a node is movable: 0 movable, 1 non-movable
%           if mask=[], it assumes all nodes are movable
%    conn:  input, a cell structure of length size(node), conn{n}
%           contains a list of all neighboring node ID for node n,
%           this can be computed from meshconn function
%    iter:  smoothing iteration number
%    useralpha: scaler, smoothing parameter, v(k+1)=(1-alpha)*v(k)+alpha*mean(neighbors)
%    usermethod: smoothing method, including 'laplacian','laplacianhc' and 'lowpass'
%    userbeta: scaler, smoothing parameter, for 'laplacianhc'
%
% output:
%    p: output, the smoothed node coordinates
%
% recommendations
%    Based on [Bade2006], 'Lowpass' method outperforms 'Laplacian-HC' in volume
%    preserving and both are significantly better than the standard Laplacian method
%
%    [Bade2006]  R. Bade, H. Haase, B. Preim, "Comparison of Fundamental Mesh
%                Smoothing Algorithms for Medical Surface Models,"
%                Simulation and Visualization, pp. 289-304, 2006.
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

p = node;
if (isempty(mask))
    nn = size(node, 1);
    idx = 1:nn;
else
    idx = find(mask == 0)';
    nn = length(idx);
end
alpha = 0.5;
method = 'laplacian';
beta = 0.5;
if (nargin > 4)
    alpha = useralpha;
    if (nargin > 5)
        method = usermethod;
        if (nargin > 6)
            beta = userbeta;
        end
    end
end
ibeta = 1 - beta;
ialpha = 1 - alpha;

for i = 1:nn
    if (length(conn{idx(i)}) == 0)
        idx(i) = 0;
    end
end
idx = idx(idx > 0);
nn = length(idx);

if (strcmp(method, 'laplacian'))
    for j = 1:iter
        for i = 1:nn
            p(idx(i), :) = ialpha * p(idx(i), :) + alpha * mean(node(conn{idx(i)}, :));
        end
        node = p;
    end
elseif (strcmp(method, 'laplacianhc'))
    for j = 1:iter
        q = p;
        for i = 1:nn
            p(idx(i), :) = mean(q(conn{idx(i)}, :));
        end
        b = p - (alpha * node + ialpha * q);
        for i = 1:nn
            p(idx(i), :) = p(idx(i), :) - (beta * b(i, :) + ibeta * mean(b(conn{idx(i)}, :)));
        end
    end
elseif (strcmp(method, 'lowpass'))
    beta = -1.02 * alpha;
    ibeta = 1 - beta;
    for j = 1:iter
        for i = 1:nn
            p(idx(i), :) = ialpha * node(idx(i), :) + alpha * mean(node(conn{idx(i)}, :));
        end
        node = p;
        for i = 1:nn
            p(idx(i), :) = ibeta * node(idx(i), :) + beta * mean(node(conn{idx(i)}, :));
        end
        node = p;
    end
end