function CompareWarpedCloneTool(varargin)
% function CompareWarpedCloneTool(varargin)
%
% A tool that takes a list of binary images and a corresponding list of rgb
% images and then lets the user click on sectors in any image. As the user
% clicks on a sector, the best fitted ellipse for that sector will be
% shown. If the user then clicks on another sector, the corresponding
% ellipse for that sector is displayed and the growth tensor that
% corresponds to the two ellipses is then calculated and shown in another
% figure.
%
% Inputs:
%
%  'bimgfilenames' - a cell array containing the full pathnames of the
%                    binary images (usually tif no compression).
%
%  'imgfilenames' - a cell array containing the full pathnames of the
%                   rgb images.
%
% Dr. Andrew I. Hanna, CMP & JIC, 2008.

error(nargchk(0,inf,nargin));
if checkSATDir(pwd) == 0
    fprintf('You are not in a valid SectorAnalysisToolbox Project directory\n');
    return;
end
data = parseargs(varargin);

% Open, move, and get the handles from the figure file.
ad.figMain = openfig(mfilename, 'reuse');
movegui(ad.figMain, 'center');
set(ad.figMain, 'Visible', 'on');
ad.handles = guihandles(ad.figMain);
set(ad.figMain, 'Color', get(ad.handles.panelA, 'BackgroundColor'));
set(ad.figMain, 'ResizeFcn', {@doResize});

ad = setupcallbacks(ad);
ad = initadstructure(ad, data);

ad.data.uitableh = uitable('Parent', ad.handles.tablepanel);

% Remove the save icon from the toolbar
%customtoolbar(ad.figMain);
setappdata(0,'CompareWarpedCloneToolData',ad);

updateSectorInfoList();
changeSectorA(ad.handles.sectorsA);
changeSectorB(ad.handles.sectorsB);
showMeanShape();
updateSectorA();
updateSectorB();
try
    uiwait(ad.figMain);
catch
    if ishandle(ad.figMain)
        delete(ad.figMain)
    end
end
value = 0;
return;
%%
function doResize(aa, evd)
updateSectorInfoList();
return

%%
function ad = setupcallbacks(ad)
% Set all the callback functions
set(ad.handles.sectorsA, 'callback', {@changeSectorA});
set(ad.handles.sectorsB, 'callback', {@changeSectorB});
set(ad.handles.compbtn, 'callback', {@compareSectors});
set(ad.handles.addbtn, 'callback', {@addCurrGTToList});
set(ad.handles.clearallbtn, 'callback', {@clearAll});
set(ad.handles.savebtn, 'callback', {@saveInfo});
set(ad.handles.loadbtn, 'callback', {@loadInfo});
set(ad.handles.updateMeshBtn, 'callback', {@updateMesh});

return
%%
function updateMesh(btn, evd)
ad = getappdata(0,'CompareWarpedCloneToolData');
ad.data.dx = str2double(get(ad.handles.dxtxt, 'String'));
ad.data.dy = str2double(get(ad.handles.dytxt, 'String'));
setappdata(0,'CompareWarpedCloneToolData',ad);
showMeanShape();
updateSectorA();
updateSectorB();
return
%%
function ad = initadstructure(ad, data)

ad.pathnames.StagedDirString =  'Staged';
ad.data.sectorA = initdata();
ad.data.sectorB = initdata();
ad.data.dx = 4;
ad.data.dy = 4;
ad.data.basegrid.x = [];
ad.data.basegrid.y = [];
ad.data.meanshape = satio_load_mean_shape;
ad.data.meanshape  = cat(1, ad.data.meanshape , ad.data.meanshape (1,:));
ad.data.D.gt = zeros(3,3);
ad.data.D.ellipses = cell(2,1);

ad.data.sector_info_cell = {};
ad.handles.gtaxes = -1;
ad.handles.sectorah = -1;
ad.handles.sectorbh = -1;

ad.data.ellipvec = {};
ad.data.linwidth = 2;
ad.data.sector_files = data.imgfilenames;
rgbfilenames = cell(1, length(ad.data.sector_files));
for i=1:length(ad.data.sector_files)
   [junk, rgbfilenames{i}] = fileparts(ad.data.sector_files{i});
end
set(ad.handles.sectorsA, 'String', rgbfilenames);
set(ad.handles.sectorsB, 'String', rgbfilenames);
return
%%
function data = initdata()
data.warped_sectors = [];
data.shape = [];
data.filename = [];
return
%%
function s = loadShape(imfilename)
s = [];
sfilename = fileparts(imfilename);
sfilename = [sfilename, filesep, 'edge_xy.mat'];
if ~exist(sfilename, 'file')
    return
end
s = load(sfilename);
if ~isfield(s, 'E')
    return;
end
s = s.E;
s = cat(1, s, s(1,:));
return
%%
function s = loadWarpedSectorInfo(imfilename)
s = [];
sfilename = fileparts(imfilename);
sfilename = [sfilename, filesep, 'warped_sector_info.mat'];
if ~exist(sfilename, 'file')
    return
end
s = load(sfilename);
if ~isfield(s, 'warped_sector_info')
    return;
end
s = s.warped_sector_info;
return
%%
function data = loadData(imfilename)
wh = waitbar(0, 'Loading...');
data.warped_sectors = loadWarpedSectorInfo(imfilename);
waitbar(0.5, wh);
data.shape = loadShape(imfilename);
data.filename = imfilename;
delete(wh);
return
%%
function changeSectorA(popup, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
val = get(popup, 'Value');
data = loadData(ad.data.sector_files{val});
ad.data.sectorA = data;
setappdata(0, 'CompareWarpedCloneToolData', ad);
updateSectorA();
return
%%
function changeSectorB(popup, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
val = get(popup, 'Value');
str = get(popup, 'String');
data = loadData(ad.data.sector_files{val});
ad.data.sectorB = data;
setappdata(0, 'CompareWarpedCloneToolData', ad);
updateSectorB();
return
%%
function showMeanShape()
ad = getappdata(0, 'CompareWarpedCloneToolData');
if isempty(ad.data.meanshape)
    return;
end
plotGrid();
S = ad.data.meanshape;
S = cat(1, S, S(1,:));
plot(ad.handles.meanshapeaxes, S(:,2), S(:,1), '-o');

axis(ad.handles.meanshapeaxes, 'image', 'ij');
return
%%
function plotGrid()
ad = getappdata(0, 'CompareWarpedCloneToolData');
if isempty(ad.data.meanshape)
    return;
end
cla(ad.handles.meanshapeaxes);
S = ad.data.meanshape;
xrange = [min(S(:,2)), max(S(:,2))];
yrange = [min(S(:,1)), max(S(:,1))];
xind = linspace(xrange(1), xrange(2), ad.data.dx+1);
yind = linspace(yrange(1), yrange(2), ad.data.dy+1);
[x, y] = meshgrid(xind, yind);

ad.data.basegrid.centroids = calcCentroids(x, y);
ad.data.basegrid.x = x;
ad.data.basegrid.y = y;
hold(ad.handles.meanshapeaxes, 'on');
mesh(x, y, -1*ones(size(x)), 'Parent', ad.handles.meanshapeaxes);
setappdata(0, 'CompareWarpedCloneToolData', ad);
plotCentroids(ad.handles.meanshapeaxes);
return
%%
function plotCentroids(ah)
ad = getappdata(0, 'CompareWarpedCloneToolData');
c = ad.data.basegrid.centroids;
%plot(ah, c(:, 1), c(:, 2), 'rd');
str = cell(size(c,1),1);
for i=1:size(c,1)
    str{i} = num2str(i);
end
text(c(:,1), c(:,2), str, 'Parent', ah, 'Color', 'k');
setappdata(0, 'CompareWarpedCloneToolData', ad);
return
%%
function c = calcCentroids(x, y)
%x = x';
%y = y'
cx = x + (circshift(x,[0 -1])-x)/2; cx = cx(1:end-1,:);  cx = cx(:,1:end-1);
cy = y + (circshift(y,[-1 0])-y)/2; cy = cy(1:end-1,:);  cy = cy(:, 1:end-1);
cx = cx';
cy = cy';
c =[cx(:), cy(:)];
return
%%
function plotSectorAGrid()
ad = getappdata(0, 'CompareWarpedCloneToolData');
mh = plotSectorGrid(ad.handles.sectorAaxes);
setappdata(0, 'CompareWarpedCloneToolData', ad);
return;
%%
function plotSectorBGrid()
ad = getappdata(0, 'CompareWarpedCloneToolData');
mh = plotSectorGrid(ad.handles.sectorBaxes);
setappdata(0, 'CompareWarpedCloneToolData', ad);
return;
%%
function mh = plotSectorGrid(ah)
ad = getappdata(0, 'CompareWarpedCloneToolData');
mh = [];
if isempty(ad.data.basegrid.x) || isempty(ad.data.basegrid.y)
    return;
end
ecol = [0 0 0];
S = ad.data.meanshape;
x = ad.data.basegrid.x;
y = ad.data.basegrid.y;
mh = mesh(x, y, ones(size(x)), 'Parent', ah, 'FaceColor', 'none', 'EdgeColor', ecol, 'HitTest', 'off');
plotCentroids(ah)
return
%%
function updateSectorA(popup, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
cla(ad.handles.sectorAaxes);
if ~isempty(ad.data.sectorA.warped_sectors)
    ad.handles.sectorah = showSectors(ad.data.sectorA.warped_sectors, ad.handles.sectorAaxes);
end
if ~isempty(ad.data.sectorA.shape)
    ad.handles.shapeah = showShape(ad.data.meanshape, ad.handles.sectorAaxes);
end
setappdata(0, 'CompareWarpedCloneToolData', ad);
compareSectors(ad.handles.compbtn);
plotSectorAGrid();
axis(ad.handles.sectorAaxes, 'image', 'ij');
return
%%
function sh = showShape(shape, ah)
sh = plot(ah, shape(:,2), shape(:,1), '-o');
return
%%
function updateSectorB(popup, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
cla(ad.handles.sectorBaxes);


if ~isempty(ad.data.sectorB.warped_sectors)
    ad.handles.sectorah = showSectors(ad.data.sectorB.warped_sectors, ad.handles.sectorBaxes);
end
if ~isempty(ad.data.sectorB.shape)
    ad.handles.shapebh = showShape(ad.data.meanshape, ad.handles.sectorBaxes);
end
setappdata(0, 'CompareWarpedCloneToolData', ad);
compareSectors(ad.handles.compbtn);
plotSectorBGrid();
axis(ad.handles.sectorBaxes, 'image', 'ij');
return
%%
function imh = showSectors(sectors, ah)
ad = getappdata(0, 'CompareWarpedCloneToolData');
imh = [];
hold(ah, 'on');
imh = [];
for i=1:length(sectors)
    gt = sectors(i).growthTensor;
    mu = sectors(i).Centroid;
    imh = cat(1, imh, gtlib_plotGrowthTensor('growth_tensor', gt, 'offset', mu, 'Parent', ah, 'linewidth', ad.data.linwidth));
end
axis(ah, 'image', 'ij');
return
%%
function saveInfo(btn, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
sector_info = ad.data.sector_info_cell;
[filename, pathname] = uiputfile('*.xls', 'Save as...');
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
    return;
end
SUCCESS = xlswrite([pathname, filesep, filename],sector_info);
return
%%
function loadInfo(btn, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
[filename, pathname] = uigetfile('*.xls', 'Load...', [pwd, filesep]);
if isequal(filename,0) || isequal(pathname,0)
    disp('User pressed cancel')
    return;
end
[num, txt, raw] = xlsread([pathname, filesep, filename]);
ad.data.sector_info_cell = raw;
setappdata(0, 'CompareWarpedCloneToolData', ad);
updateSectorInfoList();
return
%%
function clearAll(btn, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
ad.data.sector_info_cell = {};
setappdata(0, 'CompareWarpedCloneToolData', ad);
updateSectorInfoList();
%%
function addCurrGTToList(btn, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
if isempty(ad.data.D)
    return;
end
sector_info = ad.data.sector_info_cell;
deltaSt = 20;
growthrates = growthParameters2GrowthRates(ad.data.D.gt, deltaSt);
[R, C] = size(sector_info);
if R==0
    R = 1;
    sector_info{1, 1} = 'Smaj';
    sector_info{1, 2} = 'Smin';
    sector_info{1, 3} = 'Sarea';
    sector_info{1, 4} = 'Kmaj';
    sector_info{1, 5} = 'Kmin';
    sector_info{1, 6} = 'Karea';
    sector_info{1, 7} = 'Kanis';
    sector_info{1, 8} = 'Rmaj';
    sector_info{1, 9} = 'Rmin';
    sector_info{1, 10} = 'Rarea';
    sector_info{1, 11} = 'Dmaj';
    sector_info{1, 12} = 'Dmin';
    sector_info{1, 13} = 'Darea';
    sector_info{1, 14} = 'theta';
    sector_info{1, 15} = 'From Region - To Region';
    sector_info{1, 16} = 'Grid Size';
    sector_info{1, 17} = 'From Image - To Image';
    
    
    
end
sector_info{R+1, 1} = growthrates.Smaj;
sector_info{R+1, 2} = growthrates.Smin;
sector_info{R+1, 3} = growthrates.Sarea;
sector_info{R+1, 4} = growthrates.Kmaj;
sector_info{R+1, 5} = growthrates.Kmin;
sector_info{R+1, 6} = growthrates.Karea;
sector_info{R+1, 7} = growthrates.Kanis;
sector_info{R+1, 8} = growthrates.Rmaj;
sector_info{R+1, 9} = growthrates.Rmin;
sector_info{R+1, 10} = growthrates.Rarea;
sector_info{R+1, 11} = growthrates.Dmaj;
sector_info{R+1, 12} = growthrates.Dmin;
sector_info{R+1, 13} = growthrates.Darea;
sector_info{R+1, 14} = growthrates.theta;
sector_info{R+1, 15} = [num2str(ad.data.D.ellipses{1}.RegionNumber), '  - ', num2str(ad.data.D.ellipses{2}.RegionNumber)];
sector_info{R+1, 16} = num2str(ad.data.D.ellipses{1}.GridSize);
sector_info{R+1, 17} = [ad.data.D.ellipses{1}.ImageName, '  - ', ad.data.D.ellipses{2}.ImageName];
ad.data.sector_info_cell = sector_info;
setappdata(0, 'CompareWarpedCloneToolData', ad);
updateSectorInfoList();
return
%%
function updateSectorInfoList()
ad = getappdata(0, 'CompareWarpedCloneToolData');
sector_info = ad.data.sector_info_cell;
S = cell(size(sector_info,1),1);
delete(ad.data.uitableh);
set(ad.handles.tablepanel, 'Units', 'pixels');
pos = get(ad.handles.tablepanel, 'Pos');
pos(1:2) = 1;
ad.data.uitableh = uitable('Parent', ad.handles.tablepanel, 'ColumnNames',  {sector_info{1, :}}, 'Data', sector_info,  'Position', pos);
setappdata(0, 'CompareWarpedCloneToolData', ad);
%%
function compareSectors(btn, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
if get(btn, 'Value')==1
    set(ad.handles.sectorAaxes, 'ButtonDownFcn', {@sectorAClicked});
    set(ad.handles.sectorBaxes, 'ButtonDownFcn', {@sectorBClicked});
else
    set(ad.handles.sectorAaxes, 'ButtonDownFcn', '');
    set(ad.handles.sectorBaxes, 'ButtonDownFcn', '');
end
setappdata(0, 'CompareWarpedCloneToolData', ad);
%%
function ind = findClosestEllipse(pt, sectors, rad)
ind = [];

centroids = [sectors.Centroid];
centroids = reshape(centroids, 2, length(centroids)/2)';
N = size(centroids,1);
d = sqrt(sum((centroids - ones(N,1)*pt).^2,2));
[val, ind] = min(d);

if val<rad
    ind = ind(1);
else
    ind = [];
end


return
%%
function R = getRegionNumber(pt)
ad = getappdata(0, 'CompareWarpedCloneToolData');
x = ad.data.basegrid.x;
y = ad.data.basegrid.y;
cols = ad.data.dx;
rows = ad.data.dy;
colnum = sum(x(1,:)<pt(1));
rownum = sum(y(:,1)<pt(2));
R = (rownum-1)*cols + colnum;
return
%%
function sectorAClicked(ah, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
if isempty(ad.data.sectorA.warped_sectors)
    return;
end
if length(ad.data.ellipvec)==2
    if ishandle(ad.data.ellipvec{1}.h)
        delete(ad.data.ellipvec{1}.h);
    end
    if ishandle(ad.data.ellipvec{2}.h)
        delete(ad.data.ellipvec{2}.h);
    end
    ad.data.ellipvec = {};
end
pt = get(ah, 'CurrentPoint');
pt = round(pt(1,1:2));
ind = findClosestEllipse(pt, ad.data.sectorA.warped_sectors, Inf);
if ~isempty(ind)
   gt = ad.data.sectorA.warped_sectors(ind);
   gt.RegionNumber = getRegionNumber(pt);
   gt.GridSize = [ad.data.dx, ad.data.dy];
   [junk, fname] = fileparts(ad.data.sectorA.filename);
   gt.ImageName = fname;
   mu = gt.Centroid;
   col = 'y';
   if isempty(ad.data.ellipvec)
       col = 'r';
   end
   gt.h = plotGT(ah, gt, mu, col);
    if length(ad.data.ellipvec)<3
        ad.data.ellipvec = cat(2, ad.data.ellipvec, gt);
    end
end
setappdata(0, 'CompareWarpedCloneToolData', ad);
if length(ad.data.ellipvec)==2
    compareEllipses();
end
%%
function sectorBClicked(ah, evd)
ad = getappdata(0, 'CompareWarpedCloneToolData');
if isempty(ad.data.sectorB.warped_sectors)
    return;
end
if length(ad.data.ellipvec)==2
    if ishandle(ad.data.ellipvec{1}.h)
        delete(ad.data.ellipvec{1}.h);
    end
    if ishandle(ad.data.ellipvec{2}.h)
        delete(ad.data.ellipvec{2}.h);
    end
    ad.data.ellipvec = {};
end

pt = get(ah, 'CurrentPoint');
pt = round(pt(1,1:2));
ind = findClosestEllipse(pt, ad.data.sectorB.warped_sectors, Inf);
if ~isempty(ind)
    gt = ad.data.sectorB.warped_sectors(ind);
    gt.RegionNumber = getRegionNumber(pt);
    gt.GridSize = [ad.data.dx, ad.data.dy];
    [junk, fname] = fileparts(ad.data.sectorB.filename);
    gt.ImageName = fname;
     col = 'y';
   if isempty(ad.data.ellipvec)
       col = 'r';
   end
    gt.h = plotGT(ah, gt, gt.Centroid, col);
    if length(ad.data.ellipvec)<3
        ad.data.ellipvec = cat(2, ad.data.ellipvec, gt);
    end
end
setappdata(0, 'CompareWarpedCloneToolData', ad);
if length(ad.data.ellipvec)==2
        compareEllipses();
end
%%
function compareEllipses()
ad = getappdata(0, 'CompareWarpedCloneToolData');
ellipses = ad.data.ellipvec;
E1 = ellipses{1}.growthTensor;
E2 = ellipses{2}.growthTensor;

ellipses{1}
ellipses{2}

[smax1, smin1, theta1] = gtlib_growthTensor2Params(E1);
[smax2, smin2, theta2] = gtlib_growthTensor2Params(E2);
ungrow_smax = smax2/smax1;
ungrow_smin = smin2/smin1;
D2 = gtlib_growthParams2Tensor(ungrow_smax, ungrow_smin, theta2);

%D2 = zeros(3,3);
%D2(1:2, 1:2) = E2(1:2,1:2)*inv(E1(1:2,1:2));
% keep it symmetric
%D2 = (D2 + D2')/2;

%D2*E1
%E2
ad.data.D.gt = D2;
ad.data.D.ellipses = ellipses;
%[U, S, V] = svd(D2);
%[a, b, c] = gtlib_growthTensor2Params(D2)



%ad.data.ellipvec = {};
setappdata(0, 'CompareWarpedCloneToolData', ad);
plotDiffGT();
return
%%
function plotDiffGT()
ad = getappdata(0, 'CompareWarpedCloneToolData');
if ~ishandle(ad.handles.gtaxes)
    figh = figure;
    set(figh, 'Position', [0 100 300 300]);
    ad.handles.gtaxes = gca;
end
ah = ad.handles.gtaxes;
gt = ad.data.D.gt;
cla(ah); 
gtlib_plotGrowthTensor('growth_tensor', gt, 'colour', 'g', 'Parent',ah , 'linewidth', 3);
[smax, smin, theta] = gtlib_growthTensor2Params(gt);
axis(ah, 'image', 'ij');
t = sprintf('Gmax: %3.3f, Gmin: %3.3f, Theta: %3.3f', smax, smin, theta);
fprintf('%s\n', t);
title(ah, t, 'FontSize', 8);
setappdata(0, 'CompareWarpedCloneToolData', ad);
return
%%
function h = plotGT(ah, gt, pt, col)
ad = getappdata(0, 'CompareWarpedCloneToolData');
hold(ah, 'on');
if nargin<4
    col = rand(1,3);
    col = [0 1 1];
end
h = gtlib_plotGrowthTensor('growth_tensor', gt.growthTensor, 'offset', pt, 'colour', col, 'linewidth', ad.data.linwidth);
return
%%
function data = parseargs(varargin)
varargin = varargin{1};

if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error(['Arguments to ', mfilename, ' must come param/value in pairs.'])
end
data.imgfilenames = {};
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'imgfilenames'
            data.imgfilenames = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to ', mfilename, '.  Name was ' varargin{i}])
    end
end
if isempty(data.imgfilenames)
    error('Number of images must be >0');
end
return;
