function [Vsmooth, numBoundaryEdges, A] = smoothMesh(V, F, k) # V: vertices, F: face data, k: num iterations

Vsmooth = V;
V_s = size(V,2);
F_s = size(F,2);
numBoundaryEdges = 0;

# compute adjacency matrix
A = [];

for i=1:F_s
    A(F(1,i), F(2, i)) = 1;
    A(F(2,i), F(3, i)) = 1;
    A(F(3,i), F(1, i)) = 1;
end

# print number of boundary edges
boundaryEdges = find((A.-A')>0);
numBoundaryEdges = length(boundaryEdges)
# for k iterations: compute every vertex as the mean of its neighbors


for i=1:k
    Vsmooth_temp = Vsmooth;
    for j=1:V_s
        neighbours_x = find(A(:,j)>0);
        neighbours_y = find(A(j,:)>0);

        neighbours_index = unique([neighbours_x', neighbours_y]);

        sum_neighbors = 0;
        for l=1:size(neighbours_index, 2);
            sum_neighbors += Vsmooth(:,neighbours_index(l));
        end
        sum_neighbors /= size(neighbours_index,2);

        Vsmooth_temp(:,j) = sum_neighbors;
    end
    Vsmooth = Vsmooth_temp;
end

endfunction