function [interpolated_pts,interpolated_uompts]=interpolate_pts_in_triangle_centres(upts,uompts,image_area,threshold)
%%%%%%%%%
%%%%%%%%%start
% Andy H. Instead of interpolating along edges, interpolate center of
% triangles.
interp_uompts = [];
interp_pts = [];
TRI = delaunay(upts(:,1), upts(:,2),{'QJ','QJ','QJ'});
for t = 1:size(TRI,1)
    x1 = upts(TRI(t,1),1);
    x2 = upts(TRI(t,2),1);
    x3 = upts(TRI(t,3),1);
    y1 = upts(TRI(t,1),2);
    y2 = upts(TRI(t,2),2);
    y3 = upts(TRI(t,3),2);
    v1 = [x2-x1; y2-y1; 0];
    v2 = [x3-x1; y3-y1; 0];
    area = .5*(norm(cross(v1,v2))); % area of triangle
  %  if area>threshold*image_area
        base_orig_pts = [upts(TRI(t,1:3),1), upts(TRI(t,1:3),2)];
        base_inter_pts = [(base_orig_pts(1,:)+ base_orig_pts(2,:))/2; (base_orig_pts(2,:)+ base_orig_pts(3,:))/2; (base_orig_pts(3,:)+ base_orig_pts(1,:))/2];
        base_pts = [base_orig_pts; base_inter_pts];
        target_orig_pts = [uompts(TRI(t,1:3),1), uompts(TRI(t,1:3),2)];
        target_inter_pts = [(target_orig_pts(1,:)+ target_orig_pts(2,:))/2; (target_orig_pts(2,:)+ target_orig_pts(3,:))/2; (target_orig_pts(3,:)+ target_orig_pts(1,:))/2];
        target_pts = [target_orig_pts; target_inter_pts];
        T = cp2tform(base_pts, target_pts,'linear conformal');
        ip = (base_pts(1,:) + base_pts(2,:) + base_pts(3,:))/3;
        [x, y] = tformfwd(T, ip(1), ip(2));
        interp_pts = [interp_pts, ip'];
        interp_uompts = [interp_uompts, [x;y]];
   % end
end
interpolated_pts=interp_pts';
interpolated_uompts=interp_uompts';