function mrphs = create_morphs(img1, img2, points1, points2, strengths, distort)

% function mrphs = create_morphs(img1, img2, points1, points2, strengths)
%
% creats morphs of img1 and img2 assuming that the coordinates specified in
% points1 and poits2 correspond to each other.
% morph strenghts are between zero and one.
%
% Author: Farshad Moradi, MD PhD
% Place: California Institute of Technology
% Date: Sep 03, 2008 

if ~exist('distort', 'var')
   distort = 0;
end

[h,w, dummy]=size(img1); %adds on corner points 
points1=[points1; 1,1; 1,h; w,1; w,h];

[h,w, nchannels]=size(img2); %adds on corner points
points2=[points2; 1,1; 1,h; w,1; w,h];

[X,Y] = meshgrid(1:w, 1:h);
I = cell(nchannels,1);
for channel = 1:nchannels
    I{channel} = sub2ind([h,w,nchannels], Y(:), X(:), ones(h*w,1)*channel);
end

if nargout<1
    figure(1); clf;
    image(img1); colormap(gray(256));
end

mrphs = cell(length(strengths), 1);

for sx=1:length(strengths)
    if nargout<1
        tic;
    end
    s = strengths(sx);
    % interim points
    p = points1+s*(points2-points1);
    % find optimal triangulation of corresponding points
    tri = delaunay(p(:,1), p(:,2), {'QJ'});
    n = size(tri,1);

    
    % find the triangle corresponding to each point in the resulting image
    T = tsearch(p(:,1), p(:,2), tri, X(:), Y(:));
    
    temp = ones(h,w,nchannels);
    % src and dst correspond to img1 and img2, respectively
    src = cell(n,1);
    dst = cell(n,1);
    mix = cell(nchannels,n);

    for i=1:n
        % find transformations
        Xinv = inv([p(tri(i,:),1), p(tri(i,:),2), ones(3,1)]');
        M1 = [points1(tri(i,:),1), points1(tri(i,:),2), ones(3,1)]' * Xinv;
        M2 = [points2(tri(i,:),1), points2(tri(i,:),2), ones(3,1)]' * Xinv;

        f=find(T==i);
        dummy = [X(f), Y(f), ones(length(f),1)]';
        src{i} = M1*dummy;
        dst{i} = M2*dummy;
        for c=1:nchannels
            mix{c,i} = I{c}(f);
        end
    end
    
    src = cat(2, src{:});
    dst = cat(2, dst{:});
    for c=1:nchannels
        % do the interpolation
        if distort
            temp(cat(1,mix{c,:})) = interp2( img1(:,:,c), src(1,:), src(2,:) , '*cubic', 1);
        else
            temp(cat(1,mix{c,:})) = ...
                (1-s)*interp2( img1(:,:,c), src(1,:), src(2,:), '*cubic', 1 )+ ...
                s*interp2(img2(:,:,c), dst(1,:), dst(2,:), '*cubic', 1 );
        end
    end
    if nargout<1
        figure(1);
        image(temp);
        hold on;
        triplot(tri, p(:,1), p(:,2), 'LineStyle', ':', 'color', [0.5 0.7 0.7]);
        toc;
        pause;
    end
    fprintf('%d out of %d\n', sx, length(strengths));
    mrphs{sx} = temp;
end