function [tri, vertices] = interp_tri(tri, vertices, iter)
% function [tri, vertices] = interp_tri(tri, vertices, iter)
% 
% A function to interpolate the triangulation. For each triangle, find the
% center and then add this to the set of vertices.
%
% Input Params: 
%    tri - a triangulation e.g. delaunay
%    vertices - a Mx2 matrix of vertices
%    iter - the number of iterations to perfrom, (default = 1)
%
% Output Params:
%    tri - the new triangulation with extra vertices
%    vertices - the new set of points.
%
%  Note: if you want to use this with another set of points, then you had
%  better warp the interpolation pts to the new frame of reference.
%  
% Dr. A. I. Hanna (2006)
if nargin <2
    return;
end
if size(tri,1)==0
    return;
end
if size(vertices,2) ~=2
    return;
end
if nargin <3
    iter = 1;
end
for i=1:iter
    newd = [];
    for i=1:size(tri,1)
        v1 = [vertices(tri(i,2),1) - vertices(tri(i,1),1),  vertices(tri(i,2),2) - vertices(tri(i,1),2)];
        v2 = [vertices(tri(i,3),1) - vertices(tri(i,2),1),  vertices(tri(i,3),2) - vertices(tri(i,2),2)];
        newd = [newd; sum(vertices(tri(i,:),1))/3, sum(vertices(tri(i,:),2))/3];
    end
    vertices = [vertices; newd];
    tri = delaunay(vertices(:,1), vertices(:,2));
end
return;