function [img2,sfimg_safe,img1,TRI,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 = [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);

upts=startpts; 
uompts=ompts;  

all_upts=upts;
all_uompts=uompts;

%Use the following to form the triangulation on the source image
%TRI = delaunay(all_upts(:,1), all_upts(:,2),{'QJ','QJ','QJ'});
%Use the following to form the triangulation on the destination (mean) image
% this has the advantage that any distortions introduced warping to the
% mean are removed when warping back to form the reconstructed image
TRI = delaunay(uompts(:,1), uompts(:,2),{'QJ','QJ','QJ'});

if remove_flipped_triangles
    flipped_tri_ind = find_flipped_triangles(TRI, upts, uompts);
    TRI(flipped_tri_ind, :) = [];
end

triangulation_used_for_warp=TRI;
%size(TRI)
TRI = TRI';
TRI = TRI(:);
% use the triangulation from the mean
[sfimg] = triwarp(img, upts, uompts,TRI); % warp to mean shape
%[sfimg, w, a] = TPS_im_warp(img, uompts, 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;



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


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