function [eg_res] = disp_moverlay2(emesh, m_ref, m_res, varargin)
% disp_moverlay(emesh, m_ref, m_res, gmm_stain, conn)
% displays composites of reference (m_ref) with others (m_res)
% emesh is standard emesh struct
% m_ref and m_res are indices for emesh
% m_ref must have size of 1, m_res is an array
% m_ref can be also logical, then a synthetic mesh is constructed
% m_res contains the results that are displayed overlayed
%       can be either an index to emesh or a struct
%       if it's a struct, idx contains the indices and stain the meshes.
%       The meshes in emesh will be replaced by the meshes in m_res.
% gmm_stain can be either a mixture model for each triangle or the
%       precomputed GMM probabilities for each m_res
% conn is optional, if given will query the gene symbols


   % Initialize variables
   counter = 1;
   stainshow = 'Staining';
   showref = 0;
   gmmthresh = 0.5 * 255;
   conn = [];
   selt = 0;
   eg_res = [];

   
    if ~isempty(varargin)
        
        i = 1;
        while i <= length(varargin),
            switch(varargin{i}),
                case 'db'
                    i = i+1;
                    conn = varargin{i};
                case 'source'
                    i = i+1;
                    mesh_src = varargin{i};
                case 'binary'
                    i = i+1;
                    eg_res = varargin{i};
                case 'symbols'
                    i = i+1;
                    sym_all = varargin{i};
                    if length(sym_all) == 1,
                        sym_name = sym_all;
                        sym_idx = 1:length(sym_all);
                    else
                        sym_name = sym_all{1};
                        sym_idx = sym_all{2};
                    end
                otherwise
                    error('Argument %s not known', varargin{i});
            end
            i = i+1;
        end
    end

    m_resfiles = cell(length(m_res), 1);

    % If m_res is a vector, assume it's an index for emesh entries
    if isvector(m_res)
       if isempty(mesh_src),
           for i=1:length(m_res),
               % Test if we query DB for symbols
               if isempty(conn),
                   m_resfiles{i} = emesh.files{m_res(i)};
               elseif isempty(sym_all)
                   fn = emesh.files{m_res(i)};
                   sym = pg_query(conn,0, pg_insitus('symbol'),myinsitu2fn(fn));
                   m_resfiles{i} = sprintf('%s (%s)', sym{1}, fn);
               end
           end
           resmesh = emesh.stain(:, m_res);
       else
           resmesh = mesh_src(:, m_res);
           sym_idx = sym_idx(m_res);
       end
    else
       resmesh = m_res; 
    end
    
    if ~isempty(sym_all),
        if length(sym_idx) ~= size(resmesh, 2),
            error('Symbol idx has different size');
        end
        for i=1:length(m_res),
            if isempty(m_resfiles{i}),
                m_resfiles{i} = sprintf('%s', sym_name{sym_idx(i)});
            else
                m_resfiles{i} = sprintf('%s (%s)', sym_name{sym_idx(i)}, m_resfiles{i});
            end
        end
    end   
   
   if isempty(eg_res),
       t_conn = fully_connected(emesh);
       eg_res = (1 - tibin(t_conn, resmesh)) .* 255;
   else
       if isvector(m_res)
           eg_res = (1 - eg_res(:, m_res)) .* 255;
       elseif size(eg_res,2) ~= size(resmesh,2)
           error('Binarized instances have different length than results');
       else
           eg_res = (1 - eg_res) .* 255;
       end
   end
   m_selr = 1:length(m_res);
   
   if isempty(m_ref),
       refstain = [];
   else
       if size(m_ref,1) == size(emesh.stain,1),
           if islogical(m_ref),
               refstain = zeros(size(emesh.stain(:,1),1), 1);
               refstain(m_ref == 0) = 255;
           else
               refstain = m_ref;
           end
       else
           refstain = emesh.stain(:, m_ref);
       end
   end
   
   triannot{1} = 'All triangles';
   triannot{2} = 'Triangle %d';
   
   ol_method{1} = 'median';
   ol_method{2} = 'mean';
   ol_method{3} = 'geomean';
   ol_method{4} = 'harmmean';
   ol_method{5} = 'mode';
   ol_method{6} = 'quantile';
   ol_method{7} = 'SVD';
   ol_meth_no = 1;
   
   %  Create and then hide the GUI as it is being constructed.
   % minimum [360,500,600,300] w/o list
   f = figure('Visible','off','Position', [360,500,700,430], ...
       'Units', 'pixels', 'WindowButtonDownFcn', {@mousebutton_Callback});
 
   %  Construct the components
   htext = uicontrol('Style','text','String', m_resfiles{counter},...
          'Position',[500,10,120,15]);
   hlist = uicontrol('Style', 'listbox', 'String', m_resfiles,...
          'Min', 1, 'Max', 1, 'Position',[500, 30, 150, 160],...
          'Callback',{@listbox_Callback});
   htri = uicontrol('Style','text','String', triannot{1},...
          'Position',[500,200,120,15]);
   hsref = uicontrol('Style', 'checkbox', 'String', 'Show reference',...
          'Min', 0, 'Max', 1, 'Value', 0, 'Position',[500, 260, 150,30],...
          'Callback', {@showref_Callback});
      
   hoverl = uicontrol('Style', 'popupmenu', 'String', ol_method,...
          'Value', ol_meth_no, 'Position',[500, 210, 150, 40],...
          'Callback', {@overlay_Callback});
      
   hstain = uibuttongroup('visible','off', 'Units', 'pixels', 'Position',[500 300 150 110]);
   hs_u1 = uicontrol('Style','Radio','String','Staining',...
          'pos',[10 70 100 30],'parent',hstain,'HandleVisibility','off');
   hs_u2 = uicontrol('Style','Radio','String','Binary',...
          'pos',[10 40 100 30],'parent',hstain,'HandleVisibility','off');
   set(hstain,'SelectionChangeFcn',@selstain_Callback);
   set(hstain,'SelectedObject', hs_u1); 
   set(hstain,'Visible','on');      
   
   hp = uipanel('Units','Pixels','Position',[50,10,400,400]); 
   h_aref = axes('Parent', hp, 'Units','Pixels','Position',[20,20,340,160]);
   h_aov = axes('Parent', hp, 'Units','Pixels','Position',[20,220,340,160]);
%   align([hforw,hback,htext],'Center','None');
      
   % Initialize the GUI.
   %Change units to normalized so components resize automatically.
   set([f,h_aref, h_aov, hp, htext, hlist, hsref, hoverl, hstain, hs_u1, hs_u2],...
   'Units','normalized');
   %Create a plot in the axes.
   displaypanels({'overlay', 'gene'});
   % Assign the GUI a name to appear in the window title.
   set(f,'Name','Overlayed insitus')
   % Move the GUI to the center of the screen.
   movegui(f,'center')
   % Make the GUI visible.
   set(f,'Visible','on');
 
   
   %  Callbacks for simple_gui. These callbacks automatically have 
   %  access to component handles and initialized data because
   %  they are nested at a lower level.
   
    function listbox_Callback(source, eventdata)
         % Determine the selected data set.
         str = get(source, 'String');
         val = get(source,'Value');
         % Set current data to the selected data set.
        counter = val;
        set(htext, 'String', m_resfiles{m_selr(counter)});
        displaypanels({'gene'});
    end
   
    function showref_Callback(source, eventdata)
         showref = get(source,'Value');
         displaypanels({'overlay'});
    end

    function overlay_Callback(source, eventdata)
         ol_meth_no = get(source,'Value');
         displaypanels({'overlay'});
    end


    function selstain_Callback(source, eventdata)
%        disp(source);
%        disp([eventdata.EventName,'  ',... 
%             get(eventdata.OldValue,'String'),'  ', ...
%             get(eventdata.NewValue,'String')]);
%        disp(get(get(source,'SelectedObject'),'String'));
         % Determine the selected data set.
         stainshow = get(get(source,'SelectedObject'),'String');
         displaypanels({'overlay', 'gene'});
    end

    function mousebutton_Callback(source, eventdata)
        save_units = get(f, 'Units');
        set(f, 'Units', 'pixels');
        %figpos = get(f, 'Position');
        point = get(f,'CurrentPoint');
        set(source, 'Units', save_units);
        %x = point(1);
        %y = point(2);

        save_units = get(h_aov, 'Units');
        set(h_aov, 'Units', 'pixels');
        axespos = get(h_aov, 'Position');
        set(h_aov, 'Units', save_units);
        
        save_units = get(hp, 'Units');
        set(hp, 'Units', 'pixels');
        panelpos = get(hp, 'Position');
        set(hp, 'Units', save_units);

        t = emesh.t;
        p_scale = emesh.p_scale;
        
        arect = [ axespos(1) + panelpos(1), axespos(2) + panelpos(2), axespos(3), axespos(4) ];
        if rectint(arect, [point 1 1]) == 0
            return;
        end
        
        prect = [min(p_scale), max(p_scale) - min(p_scale)];
        x = (point(1) - arect(1)) * prect(3) / arect(3) + prect(1);
        y = (arect(2) + arect(4) - point(2)) * prect(4) / arect(4) + prect(2);
                
        x1 = p_scale(t(:,1),1); y1 = p_scale(t(:,1),2);
        x2 = p_scale(t(:,2),1); y2 = p_scale(t(:,2),2);
        x3 = p_scale(t(:,3),1); y3 = p_scale(t(:,3),2);
        x0 = repmat(x, [length(t),1]); y0 = repmat(y, [length(t),1]);

        b0 =  (x2 - x1) .* (y3 - y1) - (x3 - x1) .* (y2 - y1);
        b123 = [ ((x2 - x0) .* (y3 - y0) - (x3 - x0) .* (y2 - y0)) ./ b0, ...
        ((x3 - x0) .* (y1 - y0) - (x1 - x0) .* (y3 - y0)) ./ b0, ...
        ((x1 - x0) .* (y2 - y0) - (x2 - x0) .* (y1 - y0)) ./ b0 ];

        number = find(min(b123') > 0);
        selt = number;
        if number > 0,
            annot = sprintf(triannot{2}, number);
            set(htri, 'String', annot);
            m_selr = selectimages(number);
        else
            set(htri, 'String', triannot{1});
            m_selr = 1:length(m_res);
        end
        displaypanels({'gene', 'overlay'});
        
    end

    function [selected] = selectimages(number)
        selected = find(eg_res(number,:) == 0);
        if isempty(selected),
            questdlg('No genes with expression. Showing full list.', ' OK ');
            selected = 1:length(m_res);
        end
        
        % update the listbox
        for i=1:length(selected),
            selfiles{i} = m_resfiles{selected(i)};
        end
        counter = 1;
        set(hlist, 'String', selfiles);
        set(hlist, 'Value', m_selr(counter));
        set(htext, 'String', m_resfiles{m_selr(counter)});
    end


    function [overlaymesh] = overlaystaining(mesh)
        if strcmp(ol_method{ol_meth_no}, 'median'),
            overlaymesh = median(mesh, 2);
        elseif strcmp(ol_method{ol_meth_no}, 'mean'),
            overlaymesh = mean(mesh, 2);
        elseif strcmp(ol_method{ol_meth_no}, 'geomean'),
            overlaymesh = geomean(mesh, 2);
        elseif strcmp(ol_method{ol_meth_no}, 'harmmean'),
            overlaymesh = harmmean(mesh, 2);
        elseif strcmp(ol_method{ol_meth_no}, 'mode'),
            overlaymesh = mode(mesh, 2);
        elseif strcmp(ol_method{ol_meth_no}, 'quantile'),
            overlaymesh = quantile(mesh, 0.75, 2);
        elseif strcmp(ol_method{ol_meth_no}, 'SVD'),
            [U, S, V] = svd(mesh);
            curr_stain = U(:,1);
            overlaymesh = adj_mesh(curr_stain, 255, 0);
        end
    end


    function displaypanels(panels)
        
        if strmatch('overlay', panels),
            switch stainshow,
                case 'Staining'
                    overlaymesh = overlaystaining(resmesh(:,m_selr));
                case 'Binary'
                    overlaymesh = overlaystaining(eg_res(:,m_selr));
                    % TODO: overlays
            end
            if showref == 1,
                axes(h_aov); displaymeshExt(overlaymesh, refstain, selt, emesh);
            else
                axes(h_aov); displaymeshExt(overlaymesh, [], selt, emesh);
            end
%            rectangle('Position',[0,0,400,200])
        end
        if strmatch('gene', panels),
            switch stainshow,
                case 'Staining'
                    compmesh = resmesh(:,m_selr(counter));
                case 'Binary'
                    compmesh = eg_res(:, m_selr(counter));
            end
            axes(h_aref); displaymeshExt(compmesh, [], 0, emesh);
%            rectangle('Position',[0,0,400,200])
        end
    end



end
    
    
function displaymeshExt(compmesh, refmesh, highlite, emesh)

    axis ij; axis off;
    hold on;
    
    if isempty(refmesh),
        cmesh = repmat(compmesh ./ 255, [1 3]);
    else
        cmesh = zeros(length(refmesh), 3);
        cmesh(:,1) = (255 - refmesh) ./ 255;
        cmesh(:,2) = (255 - compmesh) ./ 255;
    end
    
    for i = 1:size(emesh.t,1)
        roic = fix(emesh.p_scale(emesh.t(i,:),1)');
        roir = fix(emesh.p_scale(emesh.t(i,:),2)');
        fill(roic,roir,cmesh(i,:), 'EdgeColor', cmesh(i,:));
    end;

    if highlite > 0,
        i = highlite;
        roic = fix(emesh.p_scale(emesh.t(i,:),1)');
        roir = fix(emesh.p_scale(emesh.t(i,:),2)');
        fill(roic,roir,cmesh(i,:), 'EdgeColor', [0 0 1]);
    end
    
end


% find the binarized pattern
function [bin] = tibin(t_conn, stain)

    j_e_mrf = exist('JTImrf2');
    if j_e_mrf == 8,
    	jmrf = JTImrf2(t_conn);
    else
        error('Java MRF not found - check your classpath!\n');
    end

    fprintf('Calculating binarized patterns... ');

    % datapoints
    dp = size(emesh.stain,1);
    dlen = length(selection);
        
    bpat = zeros(dp, dlen);
    for i=1:dlen, 
        jmrf.twoMesh(2, stain(:,i), stain(:,i), 2); 
        bpat(:,i) = jmrf.getLabels(1); 
    end

    bin = zeros(dp, dlen);
    for i=1:dlen, 
        if mean(stain(bpat(:,i) == 1,i)) < mean(stain(bpat(:,i) == 2,i)),
            bin(st4tmesh_bpat(:,i) == 1,i) = 1;
        else
            bin(st4tmesh_bpat(:,i) == 2,i) = 1;
        end
    end

    fprintf('done\n');

end


% from find_restricted
function [t_conn] = fully_connected(emesh)

    for i=1:length(emesh.t)
        tnfr{i} = [];   
        tnf{i} = []; 
        for j=1:size(emesh.t,2)       
            [prow, pcol] = find(emesh.t==emesh.t(i,j));
            tnfr{i} = [tnfr{i}, prow(prow ~= i)'];
        end
        for j=1:length(emesh.tn{i})
           if length(find(tnfr{i} == emesh.tn{i}(j))) > 1,
            tnf{i} = [tnf{i}, emesh.tn{i}(j)];
          end
        end
    end

    t_conn=zeros(311,311);
    for i=1:length(tnf), 
        for j=1:length(tnf{i}), 
            t_conn(i,tnf{i}(j)) = 1; 
            t_conn(i,i) = 1; 
        end;
    end

    t_conn = logical(t_conn);
end

% from displaymesh
function [adj] = adj_mesh(toadj, rmin, rmax)

    if min(toadj) == max(toadj),
        adj = repmat(255, [length(toadj), 1]);
        return;
    end

    to_r = range(toadj);
    to_min = min(toadj);
    adj_r = abs(rmax - rmin);
    
    adj = abs(rmin - ((toadj - to_min) .* (adj_r / to_r)));
    adj(adj > 255) = 255;
end