function [res] = study_clusters(emesh, conn, src, clus_iassign, clus_idx, varargin)
% [res] = study_clusters(emesh, conn, source, clus_assign, clus_idx, ...)
% Just a collection of tools to study clustering/similar meshes
% e.g. 
% study_clusters(emesh, conn, st4idx_X50, apc_r1, 0, 'display');
%   displays all clusters in apc_r1 derived from st4idx_X50
% study_clusters(emesh, conn, st4idx_X50, apc_r1, 1078, 'display', 'images', 'mrf', st4dmesh_X50);
%   displays meshes of cluster centers, MRF meshes and source images of cluster labeled 1078
% study_clusters(emesh, conn, st4idx_X50, apc_r1, 2463, 'overlay', 'mrf', st4dmesh_X50);
%   displays composite tool for cluster labeled 2463 supplying the MRF
%   meshes in st4dmes_X50
% emesh   emesh struct
% conn    database connection
% src     source of cluster images:
%         src scalar = index in emesh
%         or
%         src is emesh like structure
% clus_assign cluster assignments in src
% clus_idx -1, 0 or index in clus_assign, behaviour depends on other arguments
% 'display'   displays clusters:
%             clus_idx == 0:   displays all meshes
%             clus_idx == -1:  displays median of meshes or mrf
%             clus_idx == index in clus_assign: display contents of cluster clus_idx
%             Optional:
%             'images' displays the source images (requires conn)
%             'mrf'    displays the MRF meshes too
% 'composite' displays composite of cluster members of clus_idx using disp_composite
% 'overlay'   displays overlay of cluster members of clus_idx using disp_overlay
%             Optional 'mrf' argument supplies the MRF meshes
% 


    show_comp = 0;
    show_clus = 0;
    show_overl = 0;
    show_pcs = 0;
    show_res = 0;
    show_img = 0;
    show_svd = 0;
    bicluster_memb = [];
    pcs_iarg = {};
    res = [];
    disp_method = [];
    
    pclus = [];
    
    if clus_idx > 0 && length(find(clus_iassign == clus_idx)) == 0,
        error('Cluster %d not in clus_assign', clus_idx);
    end
    
    
    % check if source is just an index, if yes, create artificial source
    if ~isstruct(src),
        source.idx = src;
        source.stain = emesh.stain(:, src);
        source.dmesh = [];
        source.sym = [];
    else
        source = src;
    end
    

    if ~isempty(varargin)
        i = 1;
        while i <= length(varargin),
            switch(varargin{i}),
                case 'display'
                    show_clus = 1;
                case 'images'
                    show_img = 1;
                case 'composite'
                    show_comp = 1;
                case 'overlay'
                    show_overl = 1;
                case 'nogui'
                    show_res = 1;
                case 'svd'
                    show_svd = 1;
                case 'cpoly_d'
                    show_pcs = 2;
                case 'cpoly'
                    show_pcs = 1;
                case 'cpoly_arg'
                    pcs_iarg = varargin{i+1};
                    i = i+1;
                case 'bicluster'
                    bicluster_memb = varargin{i+1};
                    i = i+1;
                case 'display_method'
                    disp_method = varargin{i+1};
                    i = i+1;
                case 'mrf'
                    source.dmesh = varargin{i+1};
                    i = i+1;
%                 case 'conn'
%                     conn = varargin{i+1};
%                     i = i+1;
                otherwise
                    error('Argument %s not known', varargin{i});
            end
            i = i+1;
        end
    end
    
    % detect MySQL or PgSQL and assign symbol query SQL code
    query_sym = [];
    if ~isempty(conn),       
        if isempty(strfind(conn.URL, 'mysql')),
            query_sym = pg_insitus('symbol');
        else
            query_sym = my_insitus('symbol');
        end
    end

    if min(size(clus_iassign)) == 1,
        clus_assign = clus_iassign;
    elseif clus_idx > 0,
        clus_assign = clus_iassign(:, clus_idx);
    else
        clus_assign = clus_iassign;
    end
    
    
    if show_clus > 0,
        
        if clus_idx <= 0,
            % if clus_idx is 0 or -1, display all the cluster centers

            x_pan = 10;
            y_pan = 8;
            panels = x_pan * y_pan;

            if min(size(clus_assign)) == 1,
                clus_memb = unique(clus_assign)';
                clus_assignc = clus_assign;
            else
                clus_memb = unique(clus_assign);
                % there may be 0s in there - get rid of them
                if ~isempty(find(clus_memb) == 0),
                    clus_memb = clus_memb(clus_memb ~= 0);
                end
            end

            res = zeros(size(emesh.stain, 1), clus_memb);
            
            j = 0;
            while j < length(clus_memb),

                figure;
                final_panel = j+panels;
                if final_panel > length(clus_memb),
                    final_panel = length(clus_memb);
                end
                for i=j+1:final_panel
                    
                    if min(size(clus_assign)) == 1,
                        clus_center = clus_memb(i);
                    else
                        clus_assignc = clus_assign(:, i);
                        clus_center = find(clus_assign(:, clus_memb(i)) == clus_memb(i));
                        clus_center = clus_center(1);
                    end
                    
                    % sym=pg_query(conn,0, pg_insitus('symbol'),emesh.files{source.idx(clus_center)}); 
                    if isempty(source.sym),
                        if isempty(conn)
                            sym{1} ='';
                        else
                            sym=pg_query(conn,0, query_sym,emesh.files{clus_center});
                        end
                    else
                        sym{1} = source.sym{clus_center};
                    end
                    fprintf('%d = %s\n', clus_memb(i), sym{1});
                    
                    displaymesh_adjust = 0;
                    
                    
                    if clus_idx < 0
                        curr_stain = [];
                        if ~isempty(source.dmesh),
                            if isempty(disp_method) || strcmp(disp_method, 'median')
                                curr_stain = median(source.dmesh(:, clus_assignc == clus_memb(i)), 2);
                                curr_stain = abs(255 - curr_stain .* 255);
                            elseif strcmp(disp_method, 'mesh_logo')
                                curr_stain = mesh_logo(emesh, source.dmesh(:, clus_assignc == clus_memb(i)), 0);
                                displaymesh_adjust = -1;
                            elseif strcmp(disp_method, 'svd')
                                [U, S, V] = svd(source.dmesh(:, clus_assignc == clus_memb(i)));
                                curr_stain = U(:,1);
                                displaymesh_adjust = -1;
                            end
                        else
                            if strcmp(disp_method, 'svd')
                                [U, S, V] = svd(source.stain(:, clus_assignc == clus_memb(i)));
                                curr_stain = U(:,1);
                                displaymesh_adjust = -1;
                            elseif strcmp(disp_method, 'median'),
                                curr_stain = median(source.stain(:, clus_assignc == clus_memb(i)), 2);
                            end
                        end
                        if isempty(curr_stain),
                            error('Method %s not known or need discretized mesh', disp_method);
                        end
                    else
                        curr_stain = source.stain(:,clus_center);
                    end

                    subplot(y_pan,x_pan,i-j); 
                    % displaymesh(source.stain(:,clus_memb(i)), emesh.p_scale, emesh.t);
                    displaymesh(curr_stain, emesh.p_scale, emesh.t, displaymesh_adjust);
                    if ~isempty(bicluster_memb)
                        overlay_biclusters(bicluster_memb(:, clus_memb(i)), emesh.p_scale, emesh.t);
                    end
                    title(sprintf('%s (%d)=%d', sym{1}, clus_memb(i), length(find(clus_assignc == clus_memb(i)))));
                    res(:,i) = curr_stain;
                end
                j = j+panels;

            end

        else
            % else display the contents of the cluster in question
            x_pan = 1;
            y_pan = 8;
            same_pan = 1;
            
            if show_img > 0,
                x_pan = x_pan + 1;
                same_pan = same_pan + 1;
            end
            if ~isempty(source.dmesh),
                x_pan = x_pan + 1;
                same_pan = same_pan + 1;
            end
            
            % if only meshes are displayed, display 5 in a row
            if x_pan == 1,
                x_pan = 5;
                panels = x_pan * y_pan;
            else
                panels = y_pan;
            end
            
            clus_memb = find(clus_assign == clus_idx);
            % clus_pan = length(clus_memb);
            
            res = zeros(size(emesh.stain,1), clus_memb);
            j = 0;
            while j < length(clus_memb),
                figure;
                
                final_panel = j+panels;
                if final_panel > length(clus_memb),
                    final_panel = length(clus_memb);
                end

                for i=j+1:final_panel
                    if isempty(source.sym),
                        if isempty(conn)
                            sym{1} ='';
                        else
                            sym=pg_query(conn,0, query_sym, emesh.files{source.idx(clus_memb(i))}); 
                        end
                    else
                        sym{1} = source.sym{source.idx(clus_memb(i))};
                    end
                    fprintf('%d = %s\n', clus_memb(i), sym{1});
                    
                    this_panel = i-j + ((i-j) - 1) * (same_pan - 1);
                    this_panel_inc = 1;
                    
                    subplot(y_pan,x_pan,this_panel); 
                    displaymesh(source.stain(:,clus_memb(i)), emesh.p_scale, emesh.t); 
                    title(sprintf('%s', sym{1}));
                    if ~isempty(bicluster_memb)
                        overlay_biclusters(bicluster_memb(:, clus_idx), emesh.p_scale, emesh.t);
                    end
                    res(:,i) = source.stain(:,clus_memb(i));

                    if ~isempty(source.dmesh),
                        subplot(y_pan,x_pan,this_panel+this_panel_inc);
                        displaymesh(abs(255 - source.dmesh(:,clus_memb(i)) .* 255), emesh.p_scale, emesh.t);
                        this_panel_inc = this_panel_inc + 1;
                    end
                    
                    if show_img > 0,
                        subplot(y_pan,x_pan,this_panel+this_panel_inc);
                        displayimage(conn, emesh.files{source.idx(clus_memb(i))});
                        title(sprintf('%s', emesh.files{source.idx(clus_memb(i))}));
                        this_panel_inc = this_panel_inc + 1;
                    end
                end
                j = j+panels;

            end
        end
    end
    
    if show_comp > 0,
        disp_mcomposites(emesh, source.idx(clus_idx), source.idx(clus_assign == clus_idx), ...
            2, conn, source.stain, source.idx)
    end
    
    if show_overl > 0 || show_pcs > 1,
        if isempty(source.dmesh),
            fprintf('No MRF meshes, running MRF discretization for selected cluster\n');
            dmesh = meshdiscretize_emrf(emesh, source.stain(:, clus_assign == clus_idx), 0);
            fprintf('Done.\n');
        else
            dmesh = source.dmesh(:, clus_assign == clus_idx);
        end
        res = dmesh;
    end
    
    if show_overl > 0,
        ol.idx = source.idx(clus_assign == clus_idx);
        ol.stain = source.stain(:, clus_assign == clus_idx);
        if show_res == 0,
            disp_moverlay(emesh, source.stain(:,clus_idx), ol, dmesh, conn);
        end
        res = median(dmesh, 2);
    end
    
    if show_svd > 0,
        if isempty(source.dmesh),
            [U, S, V] = svd(source.stain(:, clus_assign == clus_idx));
        else
            [U, S, V] = svd(source.dmesh(:, clus_assign == clus_idx));
        end
        SS = diag(S);
        figure; plot(SS, 'b-'); title('Scree plot of SVD');
        
        U_ct = length(find(SS > 0));
        
        x_pan = 10;
        y_pan = 8;
        panels = x_pan * y_pan;
        if panels > U_ct,
            x_pan = fix(U_ct / y_pan);
            if x_pan <= 1,
                x_pan = 1;
                y_pan = U_ct;
            end
            panels = x_pan * y_pan;
        elseif panels < U_ct,
            U_ct = panels;
        end
        
        figure;
        for i=1:panels,
            subplot(y_pan,x_pan, i);
            % adj_U = adj_range(U(:,i), 0, 255);
            displaymesh(U(:,i), emesh.p_scale, emesh.t, -1); 
            % displaymesh(U(i,:), emesh.p_scale, emesh.t, 1); 
        end
        
        res.U = U;
        res.S = S;
        res.V = V;
        
    end
    
    
    if show_pcs > 0,
        pcs_arg = {};
        if show_res > 0,
            pcs_arg = {'nogui'};
        end
        
        if show_pcs == 3,
            if isempty(pcs_iarg),
                pcs_arg = {pcs_arg{:}, 'metric', 'euclidean', 'max_slider', 3,  'dist_cut', 1.5};
            else
                pcs_arg = {pcs_arg{:}, pcs_iarg{:}};
            end
            res = polygon_clusters_slider2(emesh, dmesh, pcs_arg{:});
        elseif show_pcs == 2,
            if isempty(pcs_iarg),
                pcs_arg = {pcs_arg{:}, 'metric', 'euclidean', 'max_slider', 3,  'dist_cut', 1.5};
            else
                pcs_arg = {pcs_arg{:}, pcs_iarg{:}};
            end
            res = polygon_clusters_slider2(emesh, dmesh, pcs_arg{:});
        else
            mesh = source.stain(:, clus_assign == clus_idx);
            pcs_arg = {pcs_arg{:}, pcs_iarg{:}};
            res = polygon_clusters_slider2(emesh, mesh, pcs_arg{:});
            %    'metric', 'euclidean', 'max_slider', 3,  'dist_cut', 1.5);
        end
        
        pclus = res;
    end
    
    if show_res > 0,
        
        if isempty(pclus),
            error('Need to run pclus too!');
        end
        
        mag_sna_rest = find_restricted(emesh, allTF_meshu.idx, [], mag_sna_clus, 60, 'method', 'naivemrf', 'full');
    end
    
    
end


function overlay_biclusters(clusters, p_scale, t)

    hold on;
    for i = 1:size(t,1)
    % for i = 1:find(clusters == 1)',
        if clusters(i) ~= 0,
            roic = fix(p_scale(t(i,:),1)');
            roir = fix(p_scale(t(i,:),2)');
            roic = [roic, roic(1)];
            roir = [roir, roir(1)];
            plot(roic, roir, 'r-');
        end
    end;
    hold off;
end    

function displaybimesh(stain, p_scale, t, adjust, clusters)
% displaymesh(stain, p_scale, t, adjust)
% adjust is optional
% adjust:    0 -> no changes
%           +1 -> adjust scale of stain in positive correlation
%           -1 -> adjust scale of stain in negative correlation
%

    if exist('adjust', 'var'),
        if adjust > 0,
            stain = adj_mesh(stain, 0, 255);
        elseif adjust < 0,
            stain = adj_mesh(stain, 255, 0);
        end
    end

    colormap(gray); caxis([1 255])
    axis ij; axis off;
    hold on;
    for i = 1:size(t,1)
        co = stain(i) / 255;
        roic = fix(p_scale(t(i,:),1)');
        roir = fix(p_scale(t(i,:),2)');
        fill(roic,roir,stain(i),'EdgeColor', [co co co]);
        if clusters(i) ~= 0,
            plot(roic, roir, 'r-');
        end
    end;
    hold off;
end


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