function [img2,sfimg_safe,img1,allTRI,upts,uompts,allTRIoriginals]=warp_to_mean_shape(startpts,ompts,img,indices,critical_triangle_expansion,edges);
%function img2=warp_to_mean_shape(pts,ompts,img,indices);
%
%pts, pointmodel from file (one column, x;y;x;y;...)
%ompts, mean shape model points (two columns x,y;x,y;...)
%img, image to be warped rgb
%indices, indices of points to be used for the model
%
%img2, image warped to mean shape (usually much smaller than img)
%      NaN is used to indicate triangles that have expanded to a
%      point that makes the pixels unreliable
%
if nargin<5
    critical_triangle_expansion=3;
end
allTRI=[];
startpts = startpts(:);
startpts = reshape(startpts, 2, length(startpts)/2)';
% keep only those landmarks that are part of the current model
%pts = startpts(indices,:);
% keep only those edges that include the current set of landmarks
[values,ii,jj]=intersect(edges(:,1),indices);
remaining_edges=edges(ii,:);
[values,ii,jj]=intersect(remaining_edges(:,2),indices);
remaining_edges=remaining_edges(ii,:);

% remove approximately duplicate points so as to minimise the
% number of small triangles
%ptstemp=startpts;
removelist=[];
temp_ind=indices; %1:length(pts);
for i=1:length(startpts) % identify any that have points nearby
    dists=((startpts(temp_ind,1)-startpts(i,1)).^2)+((startpts(temp_ind,2)-startpts(i,2)).^2);
    ind=find(dists<4); % i.e. within sqrt(2)
    if length(ind)>1
        sd=setdiff(ind', i); % exclude all but current point
        temp_ind=setdiff(temp_ind,sd); % update list
    end
end
% now have to remove edges associated with those landmarks that have now
% been removed
[values,ii,jj]=intersect(remaining_edges(:,1),temp_ind);
remaining_edges=remaining_edges(ii,:);
[values,ii,jj]=intersect(remaining_edges(:,2),temp_ind);
remaining_edges=remaining_edges(ii,:);

[sfimg, allTRIoriginals] = triwarp(img, startpts, ompts); % warp to mean shape

upts=startpts(temp_ind,:);
uompts=ompts(temp_ind,:);
% now identify the triangulation for the landmarks themselves to allow
% labelling of diagnotic images using allTRIoriginals

% now need to increase the number of points by adding new ones along the
% remaining_edges
% add 10 points to each existing line segment
NextraPts=5;
interp_pts = [];
interp_uompts=[];
for i=1:length(remaining_edges)
    xp1=upts(remaining_edges(i,1),1);
    yp1=upts(remaining_edges(i,1),2);
    xp2=upts(remaining_edges(i,2),1);
    yp2=upts(remaining_edges(i,2),2);
    direc_v = [xp2-xp1; yp2-yp1];
    if norm(direc_v)~=0
        direc_v = (direc_v/norm(direc_v))*(norm(direc_v)/(NextraPts+1));
        for j=1:NextraPts
            jj= [xp1; yp1] + j*direc_v;
            if ~isnan(jj)
                interp_pts = [interp_pts,jj];
            end
        end
    end
    xp1=uompts(remaining_edges(i,1),1);
    yp1=uompts(remaining_edges(i,1),2);
    xp2=uompts(remaining_edges(i,2),1);
    yp2=uompts(remaining_edges(i,2),2);
    direc_v = [xp2-xp1; yp2-yp1];
    if norm(direc_v)~=0
        direc_v = (direc_v/norm(direc_v))*(norm(direc_v)/(NextraPts+1));
        for j=1:NextraPts
            jj= [xp1; yp1] + j*direc_v;
            if ~isnan(jj)
                interp_uompts = [interp_uompts, [xp1; yp1] + j*direc_v];
            end
        end
    end
end
all_upts=[upts;interp_pts'];
all_uompts=[uompts;interp_uompts'];
upts=all_upts;
uompts=all_uompts;
[sfimg, TRI] = triwarp(img, upts, uompts); % warp to mean shape
%[sfimg, TRI] = triwarp(img, all_upts, all_uompts); % warp to mean shape
% now identify those triangles that are larger than expected
% scale_ratio=sqrt(size(sfimg,1)*size(sfimg,2))/sqrt(size(img,1)*size(img,2));
scale_ratio=(size(sfimg,1)*size(sfimg,2))/(size(img,1)*size(img,2));
sfimg_safe=sfimg;
[rows,cols,planes]=size(sfimg);
rimg2=double(reshape(sfimg,rows*cols,planes));
[rows1,cols1,planes1]=size(img);
rimg1=double(reshape(img,rows1*cols1,planes1));
threshold_area_ratio=critical_triangle_expansion*scale_ratio;
% go through all the triangles looking for those that expand
% too much
for i=1:size(TRI,1)
    x1 = upts(TRI(i,1),1);
    x2 = upts(TRI(i,2),1);
    x3 = upts(TRI(i,3),1);
    y1 = upts(TRI(i,1),2);
    y2 = upts(TRI(i,2),2);
    y3 = upts(TRI(i,3),2);
    v1 = [x2-x1; y2-y1; 0];
    v2 = [x3-x1; y3-y1; 0];
    area = .5*(norm(cross(v1,v2)));
    x1 = uompts(TRI(i,1),1);
    x2 = uompts(TRI(i,2),1);
    x3 = uompts(TRI(i,3),1);
    y1 = uompts(TRI(i,1),2);
    y2 = uompts(TRI(i,2),2);
    y3 = uompts(TRI(i,3),2);
    v1 = [x2-x1; y2-y1; 0];
    v2 = [x3-x1; y3-y1; 0];
    marea = .5*(norm(cross(v1,v2)));
    ratio=marea/area;
    if ratio>threshold_area_ratio % then ignore this triangle
        if nargout>2 %nargin>1
            allTRI=[allTRI; TRI(i,:)];
            mxi=upts(TRI(i,:),1); %
            myi=upts(TRI(i,:),2); % (TRI(i,:));
            mm=roipoly(img,mxi,myi); % find enclosed pixels
            marea=sum(mm(:));
            ind=find(mm); % and substitute NaN
            rimg1(ind,:)=repmat(NaN*ones(1,planes),length(ind),1);
        end
        mxi=uompts(TRI(i,:),1); %
        myi=uompts(TRI(i,:),2); % (TRI(i,:));
        mm=roipoly(sfimg,mxi,myi); % find enclosed pixels
        marea=sum(mm(:));
        ind=find(mm); % and substitute NaN
        rimg2(ind,:)=repmat(NaN*ones(1,planes),length(ind),1);
    end
end
img2=reshape(rimg2,[rows,cols,planes]);
img1=reshape(rimg1,[rows1,cols1,planes1]);
