function create_mask_for_template
% function create_mask_for_template
%
%
% Dr. A. I. Hanna (2007)

[filename, pathname] = uigetfile('*.temp_dat', 'Pick a Template file', ['Templates', filesep]);
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel');
    return
end

template = load([pathname, filesep, filename], '-mat');
template = template.PMTemplate;
[a, template_name] = fileparts(filename);

I = get(template, 'image');
pts = get(template, 'pts');
edges = get(template, 'loops');
if iscell(edges)
    edges = edges{1};
end
figh = figure;
axish = gca;
plotTemplate(I, pts, edges, axish);
use_predefined_interior = 0;
interior_h = [];
boundary_filename = [pathname, filesep, template_name,'.boundary'];
if exist(boundary_filename, 'file')
    boundary_pts_ind = load(boundary_filename, '-mat');
    boundary_index = boundary_pts_ind.boundary_index;
    interior_h = plotBoundaryPoints(axish, pts, boundary_index);
    answer = questdlg('Use found boundary points?', '', ...
        'Yes', 'No', 'Yes');
    switch answer,
        case 'Yes',
            use_predefined_interior = 1;
        case 'No'
            use_predefined_interior= 0;
    end
end
if use_predefined_interior==0
    delete(interior_h);
    boundary_index = select_boundary_points(I, pts, edges, axish);
end
close(figh);

[figh] = plotImageAndMask(I, pts, edges, boundary_index);
answer = questdlg('Would you like to calculate the ordering?', 'Ordering', 'Yes', 'No', 'Yes');
close(figh);
switch answer,
    case 'Yes',
        str = sprintf('%s\n%s', 'Please roughly outline the shape of your organ.', 'This allows us to accurately determine the ordering of the points');
        uiwait(msgbox(str,'Outline','modal'));
        not_done = 1;
        while not_done
            boundary_index = calcPointOrdering(I, pts, edges, boundary_index);
            [figh] = plotImageAndMask(I, pts, edges, boundary_index);
            answer = questdlg('Finished?', 'Finished ', ...
                'Yes', 'No', 'Yes');
            close(figh);
            switch answer,
                case 'Yes',
                    not_done = 0;
                case 'No'
                    not_done= 1;
            end
        end

end

save(boundary_filename,'boundary_index');
return
%%%
%
%%%
function figh = plotImageAndMask(I, pts, edges, boundary_index)
figh = figure;
axish = subplot(1,2,1);
plotTemplate(I, pts, edges, axish);
title('Original Template with Points');

subplot(1,2,2);
mask = roipoly(I, pts(1,boundary_index), pts(2,boundary_index));
imagesc(mask); axis(gca, 'image', 'ij');
colormap gray;
title('Current Mask');


%%%
%
%%%
function boundary_index = select_boundary_points(I, pts, edges, axish)
title(axish, sprintf('%s', 'Please select all the INTERIOR points'));
boundary_index = [];
not_done = 1;
interior_h = [];
while not_done
    delete(interior_h);
    interior_h = [];
    [x, y] = getline(axish, 'closed');
    in = find(inpolygon(pts(1,:),pts(2,:),x,y));
    boundary_index = setdiff(1:size(pts,2), in);
    interior_h = plotBoundaryPoints(axish, pts, boundary_index);
    answer = questdlg('Finished?', 'Finished ', ...
        'Yes', 'No', 'Yes');
    switch answer,
        case 'Yes',
            not_done = 0;
        case 'No'
            not_done= 1;
    end
end
return;
%%%
%
%%%
function interior_h = plotBoundaryPoints(axish, pts, in)
interior_h =  plot(axish, pts(1,in), pts(2,in), 'o', 'markeredgecolor', 'g', 'markerfacecol', 'g');

%%%
%
%%%
function plotTemplate(I, pts, edges, axish)
linecol = [1 1 1];
facecol = 'y';
edgecol = 'r';
hold(axish, 'on');
imagesc(I, 'parent', axish);
axis(axish, 'image', 'ij');
for i=1:size(edges,1)
    plot(axish, [pts(1, edges(i,1)), pts(1, edges(i,2))], [pts(2, edges(i,1)), pts(2, edges(i,2))], 'Color', linecol);
end
plot(axish, pts(1,:), pts(2,:), 'o', 'markeredgecolor', edgecol, 'markerfacecol', facecol);

%%%
%
%%%
function sorted_boundary_index = calcPointOrdering(I, pts, edges, boundary_index)
figh = figure;
axish = gca;
ordering = 1:size(boundary_index);
plotTemplate(I, pts, edges, axish);
[x, y] = getline(axish, 'closed');
close(figh);
[bw] = roipoly(I, x, y);
edge = bwboundaries(bw);
edge = fliplr(edge{1});
boundary_pts = pts(:, boundary_index)';
d = matdistance(edge, boundary_pts);
[val, map] = min(d, [], 1);
[smap, ind] = sort(map);
sorted_boundary_index = boundary_index(ind);
s_boundary_pts = pts(:, sorted_boundary_index)';
return;