function GUI_get_corr_points(morphfile)

% Author: Farshad Moradi, MD PhD
% Place: California Institute of Technology
% Date: Sep 08, 2008 
%
% A graphic user interface program to identify key points on two images 
% which will be used for a morph. 
% Points can be added by left click. Left click near an existing point
% marks it for moving. Left click at the new location moves the point to
% the new location. Right click cancels the move. Right click near a point
% that is not selected removes it (from both images).
% The corresponding points have the same color. 
% press "r" to assign random colors to each pair of points
% press "s" to save morph file
% press "m" to preview 50% morph
% press Esc to quit (a save dialog will ask if you would like to save)
% press "q" to abort.

global img1 img2 points1 points2
if ~exist('morphfile','var')
    [filename,pathname] = uiputfile('*.*','Choose A morph file');
    if isequal(filename,0) || isequal(pathname,0)
    return;
    end
    morphfile = fullfile(pathname, filename);
end
try
    load(morphfile);
    try
        dummy = double(img1);
        dummy = double(img2);
        dummy = points1+points2;
    catch
        error('Invalid morph file');
    end
catch
    [filename,pathname] = uigetfile('*.*','Choose Picture 1');
    img1= imread(fullfile(pathname, filename));
    [filename,pathname] = uigetfile('*.*','Choose Picture 2');
    img2= imread(fullfile(pathname, filename));
    points1 = [];
    points2 = [];
end
quit = 0;
try
    load('GUI_get_corr_points_config.cmat', '-MAT');
    figure(figno);
    set(figno, 'Position', fig_pos);
catch
    figno = 2;
    disp('could not load config file. Using defaults');
end
while ~quit
    [quit, save_status] = edit_morph(figno);
    if save_status==2 % ask
        yesno = questdlg('Save?', mfilename);
        if strcmpi(yesno, 'yes')
            save_status = 1;
        elseif strcmpi(yesno, 'no')
            quit = 1;
        end
    end
    if save_status == 1
        save(morphfile, 'img1', 'img2', 'points1', 'points2');
    else
        disp('Data not saved');
    end
end
fig_pos = get(figno, 'Position');
save GUI_get_corr_points_config.cmat figno fig_pos -MAT
close(figno);
    

function [quit, save_status] =  edit_morph(figno)

global img1 img2 points1 points2

% add on corner points
[y,x,z]=size(img1);
corners=[1,1;1,y;x,1;x,y]; porig1 = cat(1,points1,corners);
tri1 = delaunay(porig1(:,1), porig1(:,2), {'QJ'});

[y,x,z]=size(img2);
corners=[1,1;1,y;x,1;x,y]; porig2 = cat(1,points2,corners);
tri2 = delaunay(porig2(:,1), porig2(:,2), {'QJ'});

figure(figno); clf;
h1 = subplot(1,2,1);
imshow(img1);
hold on;

h2 = subplot(1,2,2);
imshow(img2);
hold on;

quit = 0;
clickthreshold = 5;
clrs = 'bgkmcy';


btn = 'r';

set(figno, 'pointer', 'crosshair');
subplot(1,2,1);

save_status = 2; % ask
pt = [];
newfig_pos = [];

while ~quit && save_status == 2
    if btn=='r'
        % refresh
        hp1 = [];
        hp2 = [];
        for i=1:size(points1,1)
            clr = clrs(ceil(rand(1)*length(clrs)));
            hp1(end+1) = plot(h1, points1(i,1), points1(i,2), ['.' clr], 'Markersize', 10);
            hp2(end+1) = plot(h2, points2(i,1), points2(i,2), ['.' clr], 'Markersize', 10);
        end
    end
    figure(figno);
    keydown = waitforbuttonpress;
    if gcf~=figno
        continue;
    end
    g = gca;
    if keydown
        btn = get(gcf, 'CurrentCharacter');
        quit = (btn==27);
    else
        pt = get(g, 'CurrentPoint');
        btn = get(gcf, 'SelectionType');
        if iscell(btn)
           disp(btn);
           btn = btn{1};
        end
        pt = pt(1,1:2);
        if min(pt)<0
            msgbox('Mouse pointer outside figures!');
            continue;
        end
    end
    
    if quit
        break
    elseif strcmpi(btn,'q')
        save_status = 0;
        quit = 1;
    elseif strcmpi(btn,'s')
        save_status = 1;
    elseif strcmpi(btn,'r')
        % refresh
        delete(hp1);
        delete(hp2);
    elseif strcmpi(btn,'m')
        try
            newfig = figure;
            dummy = create_morphs(double(img1), double(img2), points1, points2, 0.5);
            dummy = double(dummy{1});
            f=find(dummy<0); dummy(f)=0;
            figure(newfig);
            if ~isempty(newfig_pos)
                set(newfig, 'Position', newfig_pos);
            end
            image(dummy/max(dummy(:)));
            pause;
        catch
            disp('error morphing/displaying morph');
        end
        try
            newfig_pos = get(newfig, 'Position');
            close(newfig);
        catch
            disp('error closing morph display');
        end
    elseif strcmpi(btn,'f')
        %load fiducial points
        [filename,pathname] = uigetfile('*.*','Choose A morph file');
        try
            fid_data = load( fullfile(pathname, filename) );
            if g==h1
                plot(g, fid_data.points1(:,1), fid_data.points1(:,2), 'o', 'Markersize', 8, 'Color', [0.3 0.6 0.6]);
            elseif g==h2
                plot(g, fid_data.points2(:,1), fid_data.points2(:,2), 'o', 'Markersize', 8, 'Color', [0.3 0.6 0.6]);
            end
        catch
            disp('Error loading fiducial data');
        end
    elseif strcmpi(btn,'z') && ~isempty(pt) %zoom
        ax = axis(h1);
        % if in zoom mode then zoom out
        if ax(2)-ax(1)<size(img1,2)
            axis(h1, [0 size(img1,2) 0 size(img1,1)]+0.5);
            axis(h2, [0 size(img2,2) 0 size(img2,1)]+0.5);
        else
            % zoom in
            if g==h1
                tri = tri1(tsearch(porig1(:,1), porig1(:,2), tri1, pt(1), pt(2)),:);
                Xinv = inv([porig1(tri,1), porig1(tri,2), ones(3,1)]');
                M = [porig2(tri,1), porig2(tri,2), ones(3,1)]' * Xinv;
                p1 = pt;
                q1 = [pt,1]*M';
            else
                tri = tri2(tsearch(porig2(:,1), porig2(:,2), tri2, pt(1), pt(2)), :);
                Xinv = inv([porig2(tri,1), porig2(tri,2), ones(3,1)]');
                M = [porig1(tri,1), porig1(tri,2), ones(3,1)]' * Xinv;
                q1 = pt;
                p1 = [pt,1]*M'; 
            end
            [h w zz] = size(img1);
            p1 = min(p1(1:2)+[w,h]/4, [w h]);
            p0 = max(p1-[w,h]/2, [0 0]);
            axis(h1, [p0(1), p1(1), p0(2), p1(2)]);
            [h w zz] = size(img2);
            q1 = min(q1(1:2)+[w,h]/4, [w h]);
            q0 = max(q1-[w,h]/2, [0 0]);
            axis(h2, [q0(1), q1(1), q0(2), q1(2)]);
        end
    elseif (g==h1 || g==h2) && ~isempty(pt) 
        if isempty(points1)
            d = [];
        elseif g==h1
            d=sqrt((points1(:,1)-pt(1)).^2+(points1(:,2)-pt(2)).^2);
        else
            d=sqrt((points2(:,1)-pt(1)).^2+(points2(:,2)-pt(2)).^2);
        end
        [d,i] = min(d);
        clr = clrs(ceil(rand(1)*length(clrs)));
        if d<clickthreshold
            % move/delete
            if g==h1
                dummy = [...
                    plot(h1, points1(i,1), points1(i,2), 'xr', 'Markersize', 10), ...
                    plot(h2, points2(i,1), points2(i,2), 'dr', 'Markersize', 10)];
            else
                dummy = [...
                    plot(h1, points1(i,1), points1(i,2), 'dr', 'Markersize', 10), ...
                    plot(h2, points2(i,1), points2(i,2), 'xr', 'Markersize', 10)];
            end
            if strcmpi(btn,'alt')
                % delete
                points1(i,:) = NaN;
                points2(i,:) = NaN;
                delete(dummy);
            else
                % move
                keydown = waitforbuttonpress;
                if keydown
                    btn = get(gcf, 'CurrentCharacter');
                    quit = (btn==27);
                else
                    pt = get(g, 'CurrentPoint');
                    pt = pt(1,1:2);
                    btn = get(gcf, 'SelectionType');
                end
                delete(dummy);
                if ~quit && gca==g && ~strcmpi(btn,'alt')
                    dummy = get(hp1(i), 'Color');
                    if g==h1
                        delete(hp1(i));
                        points1(i,:) = pt;
                        hp1(i) = plot(h1, pt(1), pt(2), '.', 'Color', dummy, 'Markersize', 10);
                    else
                        delete(hp2(i));
                        points2(i,:) = pt;
                        hp2(i) = plot(h2, pt(1), pt(2), '.', 'Color', dummy, 'Markersize', 10);
                    end
                end
            end
        else
            % add
            if g==h1
                tri = tri1(tsearch(porig1(:,1), porig1(:,2), tri1, pt(1), pt(2)),:);
                Xinv = inv([porig1(tri,1), porig1(tri,2), ones(3,1)]');
                M = [porig2(tri,1), porig2(tri,2), ones(3,1)]' * Xinv;
                dummy = M*[pt';1];
                points1(end+1,:) = pt;
                points2(end+1,:) = dummy(1:2);
            else
                tri = tri2(tsearch(porig2(:,1), porig2(:,2), tri2, pt(1), pt(2)), :);
                Xinv = inv([porig2(tri,1), porig2(tri,2), ones(3,1)]');
                M = [porig1(tri,1), porig1(tri,2), ones(3,1)]' * Xinv;
                dummy = M*[pt';1];
                points2(end+1,:) = pt;
                points1(end+1,:) = dummy(1:2);
            end
            hp1(end+1) = plot(h1, points1(end,1), points1(end,2), ['.' clr], 'Markersize', 10);
            hp2(end+1) = plot(h2, points2(end,1), points2(end,2), ['.' clr], 'Markersize', 10);
        end
    end
    if ~isempty(points1)
        f = find(isnan(points1(:,1)));
        delete(hp1(f));
        delete(hp2(f));
        f = find(~isnan(points1(:,1)));
        points1=points1(f,:);
        points2=points2(f,:);
        hp1 = hp1(f);
        hp2 = hp2(f);
    end
end
set(figno, 'pointer', 'arrow');