function plot_Z_aal(p)
% This function makes slice overlays of the groups found by IRM/BCD/IHW/IDM
% and prints to files for easy access. It also calls the make_html to make
% some simple html files, which can be used to inspect the results.
%
% Input: p.samplefile, file written by inference part
%        p.A, adjacency matrices
%        p.roifiles, filenames for the ROIs used
%        p.outputdir, output directory
%--------------------------------------------------------------------------
% Written by Kasper Winther Andersen, 2013
%

% use export_fig to print figures
addpath export_fig\

% file with sample.MAP structure
load(p.samplefile)

roifiles=p.roifiles;
n_roi=size(roifiles,1);
A=p.A;
%--------------------------------------------------------------------------

% structural scan to underlay
tmp=which('spm');
[spmpath]=fileparts(tmp);
underlay_file=[spmpath filesep 'canonical' filesep 'single_subj_T1.nii'];

% directories for html and plots
htmldir = [p.outputdir filesep 'html'];
plotdir = [p.outputdir filesep 'plot'];
niidir  = [p.outputdir filesep 'nii'];

if ~exist(htmldir)
    system(['mkdir ' htmldir]);
end
if ~exist(plotdir)
    system(['mkdir ' plotdir]);
end
if ~exist(niidir)
    system(['mkdir ' niidir]);
end

% get MAP solution
Z=sample.MAP.Z;
eta=sample.MAP.eta;
noc=size(Z,1);

for i=1:size(roifiles,1)
    [dummy tmp_name]=fileparts(char(roifiles(i,:)));
    roitmp{i}=tmp_name;
end

med_eta=median(eta,3);

%--------------------------------------------------------------------------
h=figure;
[A_s,Z_s,eta_s,perm_s]=sortGraphUnipartite(A{1},Z,eta);
mySpyPlot(A_s,2,Z_s,Z_s,eta_s)
print(h,'-djpeg100', [plotdir filesep 'sorted_Z_A_MAP.jpg']);
%--------------------------------------------------------------------------

indx=[];
roiindx=[];
for i=1:n_roi
    v=spm_vol(roifiles(i,:));
    [tmp xyz]=spm_read_vols(v);
    indx = [indx; (find(tmp>0.8))];
    roiindx = [roiindx; ones(length(find(tmp>0.8)),1)*i];
end

%
load blobs_so
SO.figure=spm_figure;
set(SO.figure,'color',[0 0 0]); % set black background
SO.img(2).hold   = 0;
SO.img(2).range  = [0 64];
SO.img(2).cmap   = lines;
SO.slices        = 76:-4:-40;
SO.cbar          = [];
%SO.area.position = [0 0 1 1];
SO.slicedef      = [-80 1 80; -110 1 75];
SO.img(3)=[];

for comp=1:noc
    y2=zeros(size(tmp));
    colcount=0.5;
    
    % generate .nii file for each cluster
    if true % set this to false if you already have generated the .nii files
        for i=1:size(Z,2)
            if Z(comp,i)
                y2(indx(find(roiindx==i)))=colcount;
                colcount=colcount+1;
            end
        end
        
        y2(find((y2==0)))=NaN;
        y_tmp=y2;
        y_tmp(find((y_tmp==0)))=NaN;
        
        v_out       = v;
        v_out.dt    = [16 0];
        v_out.fname = [niidir filesep 'cluster_' myint2str(comp) '.nii'];
        spm_write_vol(v_out, y_tmp);
    end
    %----------------------------------------------------------------------
    
    % set the files
    pth = pwd;
    SO.img(1).vol = spm_vol(underlay_file);
    SO.img(2).vol = spm_vol([niidir filesep 'cluster_' myint2str(comp) '.nii']);
    
    % make the plot
    paint(SO);
    
    % print component number
    ann=annotation(SO.figure,'textbox',[0.04 0.8 0.2 0.15],...
        'String',['COMP ' myint2str(comp)],...
        'color','w',...
        'LineStyle','none',...
        'FontWeight','bold',...
        'FitBoxToText','off');
    
    % save figure
    filename=[plotdir filesep 'cluster_' myint2str(comp) '.jpg'];
    %print(SO.figure, '-djpeg100', filename);
    export_fig(filename,'-r200')
    
    % write roi names
    if true
        s2=spm_figure;
        col=lines(100);
        indx2=find(Z(comp,:));
        for ele=1:length(indx2)
            if 0.015+0.025*(ele-1)<1
                %ann=annotation(SO.figure,'textbox',[0.6 0.015+0.025*(ele-1) 0.09692 0.03],...
                ann=annotation(s2,'textbox',[0.1 0.015+0.025*(ele-1) 0.4 0.03],...
                    'String',strrep(roitmp(indx2(ele)),'_',' '),...
                    'color',col(ele,:),...
                    'LineStyle','none',...
                    'fontsize',12,...
                    'FitBoxToText','off');
            end
        end
        
        % save figure
        filename=[plotdir filesep 'string_cluster_' myint2str(comp) '.jpg'];
        %print(SO.figure, '-djpeg100', filename);
        export_fig(filename)
        close(s2);
    end    
end

% make html files 
if true
    p.img_width = 600;
    p.sep = '_';
    
    p.plotdir = plotdir;
    p.htmldir = htmldir;
    p.Z       = Z;
    p.eta     = med_eta;
    
    make_html(p)
    disp(['done! open ' htmldir filesep 'cluster01.html to browse results'])
end
