function p_embryo = defmesh2(ti, task, im1, lmn)
% DEFMESH - mesh calculations to landmarks
% 
% p_embryo = defmesh2(ti, task, im1, lmn)
%
%   ti:  TI struct (must contain p,t)
%   im1: RGB or Grayscale image (not required for task 1 & 2)
%   lmn:  predefined landmarks to use
%           1 = 4 corners of ellipses
%           2 = 8 corners of ellipses (default if no argument given)
%           3 = 16 corners of ellipses
%   task: task to do
%       currently:
%           1: Show the numbers of the vertices in the reference mesh
%           2: Show the current landmark points in the reference mesh
%           3: Interactively enter landmarks
%   p_embryo: returns the adjusted coordinates of the mesh
%

% predefined landmark points (compatible with fembryo)
%   3: 4 cornerpoints
dlmp{3} = [ 1, 87, 177, 88 ];
%   1: 8 cornerpoints
dlmp{1} = [ 1, 87, 177, 88, 25, 23, 161, 143 ];
%   2: 16 cornerpoints
dlmp{2} = [ 1, 87, 177, 88, 14, 23, 54, 121, 161, 175, 163, 143, 111, 44, 25, 9 ];
% other landmark points
dlmc{1} = [ 94 ];
dlmc{2} = [ 94 ];
dlmc{3} = [ 94 ];

if nargin < 4
    lmn = 2;
end

lmp = dlmp{lmn};
lmc = dlmc{lmn};
lm = [lmp, lmc];


t = ti.t;
p = ti.p;

[m,n]=size(t);
NPs = size(p,1);

figure; hold on; axis ij;

if task > 2
    % Draw the image before the mesh
    p_scale = [ (p(:,1) + 4) *25+100 (p(:,2)+2)*25+50 ];
    if size(im1,3) == 1
        im2(:,:,1) = im1;
        im2(:,:,2) = im1;
        im2(:,:,3) = im1;
    else
        im2 = im1;
    end;
    
    image(im2)
else
    p_scale = [ (p(:,1) + 4) *150 (p(:,2)+2)*150 ];
end


% draw triangles
trimeshc(p_scale, t);


switch task
    case 1
        % show numbers of triangle corner points
        for i = 1:NPs
            text(p_scale(i,1),p_scale(i,2),num2str(i));
        end;
        % show selected landmark points
        for i = 1:size(lm,2)
            scatter(p_scale(lm(i), 1), p_scale(lm(i), 2), 25, 'r', 'filled');
        end;
    case 2
        % show numbers of triangles
        for i = 1:size(t,1)
            % centroid is the mean of all vertices
            ty = mean(p_scale(t(i,:),1));
            tx = mean(p_scale(t(i,:),2));
            text(ty, tx, num2str(i),'VerticalAlignment','middle',...
                'HorizontalAlignment','center', 'FontSize',9);
        end
        
    case 3
        % Interactively enter the landmarks
        pos = enterLandmarks(p_scale, lm);
        if isempty(pos),
            p_embryo = [];
            return;
        end
        p_embryo = alignmesh(p, t, pos, lm);
        trimeshc(p_embryo, t)
end;


hold off
end

function trimeshc(p_scale, t)
% draw triangles
% homemade simplified version of trimesh

[m,n]=size(t);

for i=1:m
    plot(p_scale(t(i,[1:n,1]),1),p_scale(t(i,[1:n,1]),2),'-');
end;

end


function pos = enterLandmarks(p_scale, lmp)

% find corresponding coordinates of landmark points
fprintf('Please mark the red points to the corresponding location on the embryo:\n');
fprintf('  left mouse button: mark point\n');
fprintf('  right mouse button: cancel last input');
fprintf('  middle mouse button: abort!');

pos = zeros(2, size(lmp,2));
i = 1;
while i <= size(lmp,2)
    fprintf('Mark %d. point (edge number %d)\n', i, lmp(i));
    scatter(p_scale(lmp(i), 1), p_scale(lmp(i), 2), 25, 'r', 'filled');
    [Yp,Xp,button]=ginput(1);
    scatter(p_scale(lmp(i), 1), p_scale(lmp(i), 2), 25, 'g', 'filled');
    if button == 1,
        scatter(Yp, Xp, 25, 'b', 'filled');
        pos(:,i) = [Yp Xp]';
        i = i + 1;
    elseif button == 2,
        % abort
        pos = [];
        break;
    elseif button == 3,
        i = i - 1;
        if i > 0,
            scatter(pos(1,i), pos(2,i), 60, 'r', 'x');
        else
            i = 1;
        end
        continue;
    end
end;
end



function d = distance(a,b)
% DISTANCE - computes Euclidean distance matrix
%
% E = distance(A,B)
%
%    A - (DxM) matrix 
%    B - (DxN) matrix
%
% Returns:
%    E - (MxN) Euclidean distances between vectors in A and B
%
%
% Description : 
%    This fully vectorized (VERY FAST!) m-file computes the 
%    Euclidean distance between two vectors by:
%
%                 ||A-B|| = sqrt ( ||A||^2 + ||B||^2 - 2*A.B )
%    
%    The output is the same as MathWorks' (Neural Network Toolbox) 
%    'dist' funtion (ie, d = dist(A',B), where A is a (DxM) matrix and B 
%    a (DxN) matrix, returns the same as my d = distance(A,B) ), 
%    but this function executes much faster. 
%
% Example : 
%    A = rand(400,100); B = rand(400,200);
%    d = distance(A,B);

% Author   : Roland Bunschoten
%            University of Amsterdam
%            Intelligent Autonomous Systems (IAS) group
%            Kruislaan 403  1098 SJ Amsterdam
%            tel.(+31)20-5257524
%            bunschot@wins.uva.nl
% Last Rev : Oct 29 16:35:48 MET DST 1999
% Tested   : PC Matlab v5.2 and Solaris Matlab v5.3
% Thanx    : Nikos Vlassis

% Copyright notice: You are free to modify, extend and distribute 
%    this code granted that the author of the original code is 
%    mentioned as the original author of the code.

if (nargin ~= 2)
   error('Not enough input arguments');
end

if (size(a,1) ~= size(b,1))
   error('A and B should be of same dimensionality');
end

aa=sum(a.*a,1); bb=sum(b.*b,1); ab=a'*b; 
d = sqrt(abs(repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab));
end