function [img2,sfimg_safe,img1,allTRI,upts,uompts,allTRIoriginals,triangulation_used_for_warp]=warp_to_mean_shape...
    (startpts,ompts,img,indices,critical_triangle_expansion,edges,msiz,allTRI);
%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
if nargin<6
    edges=[];
end
if nargin<7
    %msiz=[];
    msiz = [fliplr(round(max(ompts) - min(ompts))), size(img,3)];

end
if nargin<8
    allTRI=[];
end
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
if ~isempty(edges)
    [values,ii,jj]=intersect(edges(:,1),indices);
    remaining_edges=edges(ii,:);
    [values,ii,jj]=intersect(remaining_edges(:,2),indices);
    remaining_edges=remaining_edges(ii,:);
end
if nargin<4
    indices = 1:size(startpts,1);
end
temp_ind=indices; %1:length(pts);



% need to get the vertices associated with landmarks and used in the
% triangulation before all the extra vertices are added
%[sfimg_ignore, allTRIoriginals] = triwarp(img, startpts, ompts); % warp to mean shape
allTRIoriginals = delaunay(ompts(:,1), ompts(:,2));%Removed ,{'QJ','QJ','QJ'});
    %fprintf(1,'delaunay options Removed ,{''QJ'',''QJ'',''QJ''}\n')



% ignore this sfimg
allTRI=[];
%allTRIoriginals=[];
upts=startpts; % whilst working with edges, maintain all indexes (temp_ind,:);
uompts=ompts;  % (temp_ind,:);
%[sfimg, TRI] = triwarp(img, upts, uompts); % warp to mean shape



% now need to increase the number of points by adding new ones along the
% remaining_edges
% add points to each existing line segment
NextraPts=1;
interp_pts = [];
interp_uompts=[];



threshold=0.05;
image_area=msiz(1)*msiz(2);
updated_pts=[upts;interp_pts]; %upts;
updated_uompts=[uompts;interp_uompts]; %uompts;

all_upts=updated_pts;
all_uompts=updated_uompts;

upts=all_upts;
uompts=all_uompts;

all_upts=upts;
all_uompts=uompts;


TRI = delaunay(all_uompts(:,1), all_uompts(:,2));%,{'QJ','QJ','QJ'});
    %fprintf(1,'delaunay options Removed ,{''QJ'',''QJ'',''QJ''}\n')




triangulation_used_for_warp=TRI;
%size(TRI)
TRI_mean = TRI';
TRI_mean = TRI_mean(:);
% use the triangulation from the mean
[sfimg] = triwarp(img, all_upts, all_uompts,TRI_mean); % warp to mean shape
%[sfimg, w, a] = TPS_im_warp(img, all_uompts, all_upts);
if any(size(sfimg)~=msiz)
    % so we have to put sfimg into the right sized image
    im=uint8(zeros(msiz));
    im(1:size(sfimg,1),1:size(sfimg,2),:)=sfimg(1:end,1:end,:);
    sfimg=im;
end
sfimg_safe=sfimg;




upts=all_upts;
uompts=all_uompts;
[rows,cols,planes]=size(sfimg);
rimg2=double(reshape(sfimg,rows*cols,planes));
[rows1,cols1,planes1]=size(img);
rimg1=double(reshape(img,rows1*cols1,planes1));


%----------------
scale_ratio=(size(sfimg,1)*size(sfimg,2))/(size(img,1)*size(img,2));
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)));
    %if area>0.001
        if area>10^10
        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
end
%----------------

img2=reshape(rimg2,[rows,cols,planes]);
img1=reshape(rimg1,[rows1,cols1,planes1]);
