function xjview(varargin)
% xjview, version 4
%
% usage 1: xjview (no argument)
%           for displaying a result img file, or multiple image files,
%           (which will be loaded later) and changing p-value or t/f-value
% usage 2: xjview(imagefilename)
%           for displaying the result img file and changing p-value or
%           t-value
%           Example: xjView spmT_0002.img
%                    xjView('spmT_0002.img')
%                    xjView mymask.img
% usage 3: xjview(imagefilename1, imagefilename2, ...)
%           for displaying the result img files and changing p-value or
%           t/f-value
%           Example: xjView spmT_0002.img spmT_0005.img spmT_0007.img
%                    xjView('spmT_0002.img', 'spmT_0003.img', 'spmT_0006.img')
%                    xjView myMask1.img myMask2.img myMask3.img
% usage 4: xjview(mnicoord, intensity)
%           for displaying where are the mni coordinates
%           mnicoord: Nx3 matrix of mni coordinates
%           intensity: (optional) Nx1 matrix, usually t values of the
%           corresponding voxels
%           Example: xjView([20 10 1; -10 2 5],[1;2])
%                    xjView([20 10 1; -10 2 5])
%           Note: to use xjview this way, you may need to modify the value
%           of M and DIM in the begining of xjview.m
%
% http://people.hnl.bcm.tmc.edu/cuixu/xjView
%
% by Xu Cui and Jian Li 2/21/2005
% last modified: 11/20/2007 (new database, all database from wfu_pickatlas, including aal
% last modified: 06/01/2007 (correct FDR corrected p-value list, change intensity to handles.intensity{1})
% last modified: 02/18/2007 (add colorbar max control)
% last modified: 11/16/2006 (keyboard shortcut for open image and open roi file)
% last modified: 06/16/2006 (spm5 compatible)
% last modified: 05/30/2006 (left/right flip, path of mask.img and templateFile.img)
% last modified: 05/08/2006 (debug CallBack_volumePush function, change handles.intensity{1} to intensity)
% last modified: 04/03/2006 (modify tr)
% last modified: 12/28/2005 (modify SPM process)
%
% Thank Sergey Pakhomov for sharing his database (MNI Space Utility).
% Thank Yuval Cohen for the maximize figure function (maximize.m)
%

warnstate = warning;
warning off;

% pre-set values
% important! you need compare the display of xjview and spm. If you find
% xjview flipped the left/right, you need to set leftrightflip = 1;
% otherwise leave it to 0.
leftrightflip = 0; 

% You only need to change M and DIM when you want to use xjview under
% 'usage 4'.
M = ...
    [-4     0     0    84;...
     0     4     0  -116;...
     0     0     4   -56;...
     0     0     0     1];
DIM = [41 48 35]';
TR = 2; % you don't need to set TR is you only use the viewing part of xjview

% system settings
try
    spmdir = spm('dir');
    %spm('defaults', 'fmri');
catch
    disp('Please add spm path.');
    warning(warnstate(1).state);
    return
end

try
    spm('defaults', 'fmri');
catch
	[];
end


if ispc
    os = 'windows';
elseif isunix
    os = 'linux';
else
    warndlg('I don''t know what kind of computer you are using. I assumed it is unix.', 'What computer are you using?');
    os = 'linux';
end
screenResolution = get(0,'ScreenSize');

xjviewpath = fileparts(which('xjview'));

% pre-set values
pValue = 0.001;
intensityThreshold = 0;
clusterSizeThreshold = 5;


% Appearance Settings
figurePosition =                    [0.100,   0.050,    0.660,    0.880];
sectionViewPosition =               [0.5,0.61,0.45,0.45];
glassViewAxesPosition =             [0.000,   0.600,    0.464,    0.400];
if screenResolution(3) <= 1024
    figurePosition =                    [0.100,   0.050,    0.700,    0.900];
    sectionViewPosition =               [0.5,       0.61,   0.45,   0.46];
    glassViewAxesPosition =             [0.000,   0.600,    0.464,    0.400];    
end

left =                              0.01;
editBoxHeight =                     0.05;
editBoxWidth =                      0.200;
editBoxLeft =                       0.100;

controlPanelPosition =              [left,     0.080,  0.500,      0.500];
stretchMatrix =                     diag([controlPanelPosition(3),controlPanelPosition(4),controlPanelPosition(3),controlPanelPosition(4)]);
controlPanelOffset =                controlPanelPosition' .* [1,1,0,0]';
heightUnit  =                       0.055;

sliderPosition =                    stretchMatrix*[0.000,   0*heightUnit,    1.000,           editBoxHeight]' + controlPanelOffset;
pValueTextPosition =                stretchMatrix*[0.000,   1*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
pValueEditPosition =                stretchMatrix*[0.110,   1*heightUnit,    editBoxWidth*4/3,    editBoxHeight]' + controlPanelOffset;
intensityThresholdTextPosition =    stretchMatrix*[0.400,   1*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
intensityThresholdEditPosition =    stretchMatrix*[0.520,   1*heightUnit,    editBoxWidth*3/3,    editBoxHeight]' + controlPanelOffset;
dfTextPosition =                    stretchMatrix*[0.740,   1*heightUnit,    0.8-0.74,    editBoxHeight]' + controlPanelOffset;
dfEditPosition =                    stretchMatrix*[0.800,   1*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
clusterSizeThresholdTextPosition =  stretchMatrix*[0.000,   2*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
clusterSizeThresholdEditPosition =  stretchMatrix*[0.150,   2*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
pickThisClusterPushPosition =       stretchMatrix*[0.400,   2*heightUnit,    editBoxWidth*1,    editBoxHeight]' + controlPanelOffset;
selectThisClusterPushPosition =     stretchMatrix*[0.600,   2*heightUnit,    editBoxWidth*1,    editBoxHeight]' + controlPanelOffset;
clearSelectedClusterPushPosition =   stretchMatrix*[0.800,   2*heightUnit,    editBoxWidth*1,    editBoxHeight]' + controlPanelOffset;
thisClusterSizeTextPosition =       stretchMatrix*[0.600,   2*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
thisClusterSizeEditPosition =       stretchMatrix*[0.800,   2*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
loadImagePushPosition =             stretchMatrix*[0.000,   3*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
imageFileEditPosition =             stretchMatrix*[0.200,   3*heightUnit,    1-editBoxWidth,    editBoxHeight]' + controlPanelOffset;
saveImagePushPosition =             stretchMatrix*[0.000,   4*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
saveImageFileEditPosition =         stretchMatrix*[0.200,   4*heightUnit,    1-editBoxWidth,    editBoxHeight]' + controlPanelOffset;
saveResultPSPushPosition =          stretchMatrix*[0.000,   5*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
resultPSFileEditPosition =          stretchMatrix*[0.200,   5*heightUnit,    1-editBoxWidth,    editBoxHeight]' + controlPanelOffset;
displayIntensityTextPosition =      stretchMatrix*[0.000,   3*heightUnit,    editBoxWidth+0.1,    editBoxHeight]' + controlPanelOffset;
allIntensityRadioPosition =         stretchMatrix*[0.250,   3*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
positiveIntensityRadioPosition =    stretchMatrix*[0.400,   3*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
negativeIntensityRadioPosition =    stretchMatrix*[0.550,   3*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
renderViewCheckPosition =           stretchMatrix*[0.780,   3*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
hideControlPushPosition =           stretchMatrix*[0.200,   4*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
volumePushPosition =                stretchMatrix*[0.000,   4*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
commonRegionPushPosition =          stretchMatrix*[0.200,   4*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
displayPushPosition =               stretchMatrix*[0.200,   4*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
allinonePushPosition =              stretchMatrix*[0.400,   4*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
searchPushPosition =                stretchMatrix*[0.000,   6*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
searchContentEditPosition =         stretchMatrix*[0.200,   6*heightUnit,    editBoxWidth*2,    editBoxHeight]' + controlPanelOffset;
searchTextPosition =                stretchMatrix*[0.600,   6*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
searchEnginePopPosition =           stretchMatrix*[0.600,   6*heightUnit,    1-0.6,    editBoxHeight]' + controlPanelOffset;
overlayPushPosition =               stretchMatrix*[0.000,   5*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
overlayEditPosition =               stretchMatrix*[0.200,   5*heightUnit,    0.6-editBoxWidth,    editBoxHeight]' + controlPanelOffset;
overlayPopPosition =                stretchMatrix*[0.600,   5*heightUnit,    1-0.6,    editBoxHeight]' + controlPanelOffset;
helpPosition =                      stretchMatrix*[0.800,   5*heightUnit,    editBoxWidth,    editBoxHeight]' + controlPanelOffset;
infoTextBoxPosition =               stretchMatrix*[0.000,   8*heightUnit,    1,           editBoxHeight*9]' + controlPanelOffset;
%xjViewPosition =                    stretchMatrix*[0.400,   13*heightUnit,    editBoxWidth*2.5,    editBoxHeight*3]' + controlPanelOffset;

sectionViewListboxPosition =        [sectionViewPosition(1)+0.4, sectionViewPosition(2)+0.02, 0.1, 0.14];
sectionViewMoreTargetPushPosition = [sectionViewListboxPosition(1),sectionViewListboxPosition(2)-0.02,0.10,0.02];
xHairCheckPosition =                [sectionViewListboxPosition(1),sectionViewListboxPosition(2)+0.14,0.15,0.02];
setTRangeEditPosition =             [sectionViewListboxPosition(1),sectionViewListboxPosition(2)-0.06,0.10,0.02];
setTRangeTextPosition =             [sectionViewListboxPosition(1),sectionViewListboxPosition(2)-0.04,0.10,0.02];

getStructurePushPosition =          [glassViewAxesPosition(1), glassViewAxesPosition(2)-0.06, editBoxWidth/2, editBoxHeight/2];
structureEditPosition =             [getStructurePushPosition(1), getStructurePushPosition(2), 1, getStructurePushPosition(4)];
framePosition =                     (controlPanelPosition - controlPanelOffset')*1.05 + 0.95*controlPanelOffset';

% draw figure and control
figureBKcolor=[176/255 252/255 188/255];
figureBKcolor=get(0,'Defaultuicontrolbackgroundcolor');
f = figure('unit','normalized','position',figurePosition,'Color',figureBKcolor,'defaultuicontrolBackgroundColor', figureBKcolor,...
        'Name','xjView', 'NumberTitle','off','resize','off','CloseRequestFcn', {@CallBack_quit, warnstate(1).state}, 'visible','off');
handles = guihandles(f);
 
% databases
try
    X = load('TDdatabase');
    handles.DB = X.DB;
    handles.wholeMaskMNIAll = X.wholeMaskMNIAll;
catch
    errordlg('I can''t find TDdatabase.mat','TDdatabase not found');
end

handles.figure = f;
handles.frame = uicontrol(handles.figure,'style','frame',... 
        'unit','normalized',...
        'position',framePosition,...
        'Visible','off');  
handles.slider = uicontrol(handles.figure,'style','slider',... 
        'unit','normalized',...
        'position',sliderPosition,...
        'max',1,'min',0,...
        'sliderstep',[0.01,0.10],...
        'callback',@CallBack_slider,...
        'value',0,'Visible','on');
handles.pValueTextPosition = uicontrol(handles.figure,'style','text',...
        'unit','normalized','position',pValueTextPosition,...
        'string','pValue=','horizontal','left');    
handles.pValueEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',pValueEditPosition,...
        'horizontal','left',...
        'String', num2str(pValue),...
        'BackgroundColor', 'w',...
        'callback',@CallBack_pValueEdit);
handles.intensityThresholdText = uicontrol(handles.figure,'style','text',...
        'unit','normalized','position',intensityThresholdTextPosition,...
        'string',' intensity=','horizontal','left');        
handles.intensityThresholdEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',intensityThresholdEditPosition,...
        'horizontal','left',...
        'BackgroundColor', 'w',...
        'String', num2str(intensityThreshold),...
        'callback',@CallBack_intensityThresholdEdit);
handles.dfText = uicontrol(handles.figure,'style','text',...
        'unit','normalized','position',dfTextPosition,...
        'string','df= ','horizontal','right');        
handles.dfEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',dfEditPosition,...
        'horizontal','left',...
        'BackgroundColor', 'w',...
        'String', '',...
        'callback',@CallBack_dfEdit);    
handles.clusterSizeThresholdText = uicontrol(handles.figure,'style','text',...
        'unit','normalized','position',clusterSizeThresholdTextPosition,...
        'string','cluster size >=','horizontal','left');
handles.clusterSizeThresholdEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',clusterSizeThresholdEditPosition,...
        'horizontal','left',...
        'BackgroundColor', 'w',...
        'String', num2str(clusterSizeThreshold),...
        'callback',@CallBack_clusterSizeThresholdEdit);
handles.thisClusterSizeText = uicontrol(handles.figure,'style','text',...
        'unit','normalized','position',thisClusterSizeTextPosition,...
        'string','size= ','horizontal','right', 'visible', 'off');
handles.thisClusterSizeEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',thisClusterSizeEditPosition,...
        'horizontal','left',...
        'Enable', 'inactive',...        
        'String', '','visible','off');   

handles.imageFileEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',imageFileEditPosition,...
        'horizontal','left',...
        'String', '',...
        'BackgroundColor', 'w',...
        'callback',@CallBack_imageFileEdit,...
        'visible','off');
handles.saveImageFileEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',saveImageFileEditPosition,...
        'horizontal','left',...
        'BackgroundColor', 'w',...
        'String', 'myMask.img',...
        'callback',@CallBack_saveImageFileEdit,...
        'visible','off');    
handles.saveResultPSEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',resultPSFileEditPosition,...
        'horizontal','left',...
        'BackgroundColor', 'w',...
        'String', 'myResult.ps',...
        'callback',@CallBack_saveResultPSEdit,...
        'visible','off');    
handles.loadImagePush = uicontrol(handles.figure,'style','push',...
        'unit','normalized','position',loadImagePushPosition,...
        'string','Load Image','callback',@CallBack_loadImagePush,...
        'visible','off');    
handles.saveImagePush = uicontrol(handles.figure,'style','push',...
        'unit','normalized','position',saveImagePushPosition,...
        'string','Save Image','callback',@CallBack_saveImagePush,...
        'visible','off');  
handles.saveResultPSPush = uicontrol(handles.figure,'style','push',...
        'unit','normalized','position',saveResultPSPushPosition,...
        'string','Save Result','callback',@CallBack_saveResultPSPush,...
        'visible','off');      
handles.getStructurePush = uicontrol(handles.figure,'style','push',...
        'unit','normalized','position',getStructurePushPosition,...
        'string','Get Structure','callback',@CallBack_getStructurePush,'visible','off');      
handles.structureEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',structureEditPosition,...
        'horizontal','center',...
        'enable', 'on',...
        'UserData',struct(...
		'hReg',	[],...
		'M',	M,...
		'D',	DIM,...
		'xyz',	[0 0 0]	));
handles.pickThisClusterPush = uicontrol(handles.figure,'style','push',...
        'unit','normalized','position',pickThisClusterPushPosition,...
        'string','Pick Cluster/Info','callback',@CallBack_pickThisClusterPush);          
handles.selectThisClusterPush = uicontrol(handles.figure,'style','push',...
        'unit','normalized','position',selectThisClusterPushPosition,...
        'string','Select Cluster','callback',@CallBack_selectThisClusterPush);          
handles.clearSelectedClusterPush = uicontrol(handles.figure,'style','push',...
        'unit','normalized','position',clearSelectedClusterPushPosition,...
        'string','Clear Selection','callback',@CallBack_clearSelectedClusterPush);          
    
handles.displayIntensityText = uicontrol(handles.figure,'style','text',...
        'unit','normalized','position',displayIntensityTextPosition,...
        'string','display intensity','horizontal','left'); 
handles.allIntensityRadio = uicontrol(handles.figure,'style','radio',...
        'unit','normalized',...
        'string','All',...
        'position',allIntensityRadioPosition,...
        'value', 1,...
        'callback',@CallBack_allIntensityRadio);
handles.positiveIntensityRadio = uicontrol(handles.figure,'style','radio',...
        'unit','normalized',...
        'string','Only +',...
        'position',positiveIntensityRadioPosition,...
        'callback',{@CallBack_allIntensityRadio,'+'});
handles.negativeIntensityRadio = uicontrol(handles.figure,'style','radio',...
        'unit','normalized',...
        'string','Only -',...
        'position',negativeIntensityRadioPosition,...
        'callback',{@CallBack_allIntensityRadio,'-'});
handles.renderViewCheck = uicontrol(handles.figure,'style','checkbox',... 
        'unit','normalized',...
        'string','Render View' ,...
        'horizontal', 'right',...
        'position',renderViewCheckPosition,...
        'callback', @CallBack_renderViewCheck);   
    
handles.sectionViewListbox = uicontrol(handles.figure,'style','listbox',...
        'unit','normalized',...
        'String', {'single T1','avg152PD','avg152T1','avg152T2','avg305T1','ch2','ch2bet','aal','brodmann'}, ...
        'value',3,...
        'position',sectionViewListboxPosition,...
        'callback',@CallBack_sectionViewListbox);
    
handles.xHairCheck = uicontrol(handles.figure,'style','checkbox',... 
        'unit','normalized',...
        'string','XHairs Off' ,...
        'horizontal','left',...
        'position',xHairCheckPosition,...
        'callback',@CallBack_xHairCheck);  
handles.sectionViewMoreTargetPush = uicontrol(handles.figure,'style','push',...
        'unit','normalized','position',sectionViewMoreTargetPushPosition,...
        'string','other ...','callback',@CallBack_sectionViewMoreTargetPush);   
handles.setTRangeEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',setTRangeEditPosition,'BackgroundColor', 'w',...
        'string','auto','callback',@CallBack_setTRangeEdit);
handles.setTRangeText = uicontrol(handles.figure,'style','text',...
        'unit','normalized','position',setTRangeTextPosition,...
        'string','colorbar max');    
handles.hideControlPush = uicontrol(handles.figure, 'style', 'push',...
        'unit','normalized',...
        'String', '<', 'position', hideControlPushPosition,...
        'visible','off');
handles.volumePush = uicontrol(handles.figure, 'style', 'push',...
        'unit','normalized',...
        'String', 'volume', ...
        'position', volumePushPosition,...
        'callback', @CallBack_volumePush);    
handles.commonRegionPush = uicontrol(handles.figure, 'style', 'push',...
        'unit','normalized',...
        'String', 'common region', ...
        'position', commonRegionPushPosition,...
        'callback', @CallBack_commonRegionPush);       
handles.displayPush = uicontrol(handles.figure, 'style', 'push',...
        'unit','normalized',...
        'String', 'display', ...
        'position', displayPushPosition,...
        'callback', @CallBack_displayPush,...
        'visible','off');    
handles.allinonePush = uicontrol(handles.figure, 'style', 'push',...
        'unit','normalized',...
        'String', 'all in one', ...
        'position', allinonePushPosition,...
        'callback', @CallBack_allinonePush,...
        'visible','off');        
handles.searchPush = uicontrol(handles.figure, 'style','push',...
        'unit','normalized','position',searchPushPosition,...
        'String', 'search','callback',@CallBack_searchPush, 'ForeGroundColor',[0 0 1]);
handles.searchContentEdit = uicontrol(handles.figure, 'style','edit',...
        'unit','normalized','position',searchContentEditPosition,...
        'ForeGroundColor',[0 0 1],...
        'BackgroundColor', 'w',...
        'horizontal','left');
handles.searchText = uicontrol(handles.figure, 'style','text',...
        'unit','normalized','position',searchTextPosition,...
        'string', '  in',...
        'horizontal','left',...
        'visible','off');
handles.searchEnginePop = uicontrol(...
		'Units','normalized', ...
		'ListboxTop',0, ...
		'Position',searchEnginePopPosition, ...
		'String',{'in xBrain.org';'in google scholar';'in pubmed';'in wiki'}, ...
		'Style','popupmenu', ...
		'value',1);    
handles.overlayPush = uicontrol(handles.figure, 'style','push',...
        'unit','normalized','position',overlayPushPosition,...
        'String', 'overlay','callback',@CallBack_overlayPush);
handles.overlayEdit = uicontrol(handles.figure, 'style','edit',...
        'unit','normalized','position',overlayEditPosition,...
        'BackgroundColor', 'w',...
        'horizontal','left',...
        'callback', @CallBack_overlayEdit);
handles.overlayPop = uicontrol(handles.figure, 'style','popupmenu',...
        'unit','normalized','position',overlayPopPosition,...
        'string', sort(fieldnames(handles.wholeMaskMNIAll)),...
        'horizontal','left',...
        'callback', @CallBack_overlayPop);
    
handles.helpPush = uicontrol(handles.figure, 'style','push',...
        'unit','normalized','position',helpPosition,...
        'String', 'help','callback','web http://people.hnl.bcm.tmc.edu/cuixu/xjView','ForeGroundColor',[0 0 1],...
        'horizontal','left', ...
        'visible','off');    
handles.infoTextBox = uicontrol(handles.figure, 'style','edit',...
        'unit','normalized','position',infoTextBoxPosition,...
        'String', 'Welcome to xjView 4','ForeGroundColor','k', 'BackgroundColor', 'w',...
        'horizontal','left', ...'fontname','times',...
        'max',2, 'min',0); 
try
    urlread('http://people.hnl.bcm.edu/cuixu/xjView/guestbook/stat.php');
	s = urlread('http://people.hnl.bcm.tmc.edu/cuixu/xjView/toUser.txt');
    set(handles.infoTextBox, 'String', s);
end
 
handles.glassViewAxes = axes('unit','normalized','position',glassViewAxesPosition,'XTick',[],'YTick',[],'visible','off');


handles.testEdit = uicontrol(handles.figure,'style','edit',...
        'unit','normalized','position',[0.1 0.4 0.2 0.05],...
        'horizontal','left',...
        'callback',@test,...
        'visible','off');

% menu
cSHH = get(0,'ShowHiddenHandles');
set(0,'ShowHiddenHandles','on')
hMenuFile = findobj(get(handles.figure,'Children'),'flat','Label','&File');
if ~isempty(hMenuFile)
	hMenuFileOpen = findobj(get(handles.figure,'Children'),'Label','&Open...');
	set(hMenuFileOpen, 'label', 'Open Figure...');
	hMenuFileSave = findobj(get(handles.figure,'Children'),'Label','&Save');
	set(hMenuFileSave, 'label', 'Save Figure ...');
	hMenuFileSaveAs = findobj(get(handles.figure,'Children'),'Label','Save &As...');
	set(hMenuFileSaveAs, 'label', 'Save Figure As ...');
else
    hMenuFile = uimenu(handles.figure, 'label', '&File');
end

set(hMenuFile,'ForegroundColor',[0 0 1]);
set(findobj(hMenuFile,'Position',1),'Separator','on');
uimenu('Parent',hMenuFile,'Position',1,'ForegroundColor',[0 0 1],...
	'Label','Open Images (*.img) ...',...
	'CallBack',@CallBack_loadImagePush, 'Accelerator', 'o');
uimenu('Parent',hMenuFile,'Position',2,'ForegroundColor',[0 0 1],...
	'Label','Save Current Image (*.img) ...',...
	'CallBack',{@CallBack_saveImagePush, '', 0});
uimenu('Parent',hMenuFile,'Position',3,'ForegroundColor',[0 0 1],...
	'Label','Save Current Image as Mask (*.img) ...',...
	'CallBack',{@CallBack_saveImagePush, '', 1});
uimenu('Parent',hMenuFile,'Position',4,'ForegroundColor',[0 0 1],...
	'Label','Save Result (*.ps/pdf) ...',...
	'CallBack',@CallBack_saveResultPSPush);
    
hMenuHelp = findobj(get(handles.figure,'Children'),'flat','Label','&Help');
if isempty(hMenuHelp)
    hMenuHelp = uimenu(handles.figure, 'label', 'xjView &Help');
end
set(hMenuHelp,'ForegroundColor',[0 0 1]);
uimenu('Parent',hMenuHelp,'Position',1,'ForegroundColor',[0 0 1],...
	'Label','xBrain.org: brain mapping database',...
	'CallBack','web http://www.xbrain.org -browser');
uimenu('Parent',hMenuHelp,'Position',2,...
	'Label','xjview help','ForegroundColor',[0 0 1],...
	'CallBack','web http://people.hnl.bcm.tmc.edu/cuixu/xjView');
set(findobj(hMenuHelp,'Position',3),'Separator','on');
set(0,'ShowHiddenHandles',cSHH)

if exist('cuixuBOLDretrieve')
    hMenuAnalyze = uimenu('label','&Analyze','ForegroundColor',[0 0 1],'visible','on');
    hMenuPreprocess = uimenu(hMenuAnalyze,'label','Preprocess','ForegroundColor',[0 0 1],'callback', @CallBack_preprocess);
    hMenuProcess = uimenu(hMenuAnalyze,'label','Process (GLM estimation)','ForegroundColor',[0 0 1],'callback',@CallBack_process);
    hMenuSPMProcess = uimenu(hMenuAnalyze,'label','SPMProcess (GLM using SPM)','ForegroundColor',[0 0 1],'callback',@CallBack_SPMProcess);
    hMenuGLMPeak = uimenu(hMenuAnalyze,'label','GLM on peak BOLD','ForegroundColor',[0 0 1],'callback',@CallBack_GLMPeak);
    hMenuContrast = uimenu(hMenuAnalyze,'label','Contrast','ForegroundColor',[0 0 1],'callback',@CallBack_contrast);
    hMenuFDR = uimenu(hMenuAnalyze,'label','FDR','ForegroundColor',[0 0 1],'callback',@CallBack_fdr);
    hMenuROI = uimenu(hMenuAnalyze,'label','ROI: retrieve signal','ForegroundColor',[0 0 1],'callback',@CallBack_timeSeries);
    hMenuROIPlot = uimenu(hMenuAnalyze,'label','ROI: plot','ForegroundColor',[0 0 1],'callback',@CallBack_plotROI, 'Accelerator', 'm');
    hMenuROIIndividualPlot = uimenu(hMenuAnalyze,'label','ROI: individual plot','ForegroundColor',[0 0 1],'callback',@CallBack_plotIndividualROI);   
    hMenuROIIndividualPlotWithBehavior = uimenu(hMenuAnalyze,'label','ROI: individual plot with behavior','ForegroundColor',[0 0 1],'callback',@CallBack_plotIndividualROIWithBehavior);
    hMenuROICorrelationPlot = uimenu(hMenuAnalyze,'label','ROI: correlation plot','ForegroundColor',[0 0 1],'callback',@CallBack_plotCorrelationROI);
    %hMenuWholeBrainCorrelation = uimenu(hMenuAnalyze,'label','Whole brain correlation','ForegroundColor',[0 0 1],'callback',@CallBack_wholeBrainCorrelation);
    hMenuBehaviorAnalysis = uimenu(hMenuAnalyze,'label','Behavior analsyis','ForegroundColor',[0 0 1],'callback',@CallBack_behaviorAnalysis);
	hMenuHeadMovementAnalysis = uimenu(hMenuAnalyze,'label','Head movement analsyis','ForegroundColor',[0 0 1],'callback',@CallBack_headMovementAnalysis);
    hMenuModelComparison = uimenu(hMenuAnalyze,'label','Linear model comparison','ForegroundColor',[0 0 1],'callback',@CallBack_modelComparison);
    
    hMenuHNLOnly = uimenu('label','For H&NL Only','ForegroundColor',[0 0 1],'visible','on');
    hMenuNew2Old = uimenu(hMenuHNLOnly,'label','Format convert','ForegroundColor',[0 0 1],'callback',@CallBack_new2old);
    hMenuPreprocessCluster = uimenu(hMenuHNLOnly,'label','Preprocess (on cluster)','ForegroundColor',[0 0 1],'callback',{@CallBack_preprocess, 'cluster'});
    hMenuProcessCluster = uimenu(hMenuHNLOnly,'label','Process (GLM, on cluster)','ForegroundColor',[0 0 1],'callback',{@CallBack_process, 'cluster'});
    hMenuSPMProcessCluster = uimenu(hMenuHNLOnly,'label','SPM Process (GLM using SPM, on cluster)','ForegroundColor',[0 0 1],'callback',{@CallBack_SPMProcess, 'cluster'});
    hMenuGLMPeakCluster = uimenu(hMenuHNLOnly,'label','GLM on peak BOLD (on cluster)','ForegroundColor',[0 0 1],'callback',{@CallBack_GLMPeak, 'cluster'});
end

figurecm = uicontextmenu;
uimenu(figurecm,'label','Red','callback','set(gcf,''color'',''r'')');
uimenu(figurecm,'label','White','callback','set(gcf,''color'',''w'')');
uimenu(figurecm,'label','Gray','callback','set(gcf,''color'',[0.925, 0.914,  0.847])');
%set(handles.figure,'uicontextmenu',figurecm);
set(handles.figure,'WindowButtonDownFcn',@figureMouseUpFcn);

set(handles.figure,'visible','on');

% save pre-set values
handles.system = os;
handles.spmdir = spmdir;
handles.screenResolution = screenResolution;
handles.pValue = pValue;
handles.intensityThreshold = intensityThreshold;
handles.clusterSizeThreshold = clusterSizeThreshold;
handles.sectionViewPosition = sectionViewPosition;
handles.sectionViewTargetFile = getSectionViewTargetFile(spmdir, 'avg152T1');

guidata(f, handles);

% global variables for rotation matrix M and dimension
global M_;
global DIM_;
global TR_;
global LEFTRIGHTFLIP_;
global TMAX_; % colorbar max to display in section view
M_ = M;
DIM_ = DIM;
TR_ = TR;
LEFTRIGHTFLIP_ = leftrightflip;
TMAX_ = 'auto';

% check input arguments
if length(varargin) == 0
    [];
elseif isstr(varargin{1})
    CallBack_loadImagePush(handles.loadImagePush, [], varargin);
else
    mniCoord = varargin{1};
    if length(varargin) < 2
        intensity = ones(size(mniCoord,1),1);
    else
        intensity = varargin{2};
    end
    thisStruct.mni = mniCoord;
    thisStruct.intensity = intensity;
    thisStruct.M = M;
    thisStruct.DIM = DIM';    
    CallBack_loadImagePush(handles.loadImagePush, [], thisStruct);
end



function test(hObject, eventdata)
handles = guidata(gcbo);
set(hObject, 'String', num2str(handles.pValue));
vars = evalin('base','who');
vars
x = evalin('base',vars{1});
x
handles.DIM


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% FDR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_fdr(hObject, eventdata)
handles = guidata(hObject);
set(handles.infoTextBox, 'string', {'Right now FDR only works for T-test image.'});
if length(handles.imageFileName) > 1
    set(handles.infoTextBox, 'string', {'I can only work for a single image file. You opened multiple files.'});
    return
end

q = get(handles.pValueEdit,'String');
q = str2num(q);
if get(handles.allIntensityRadio, 'Value')
    positive = 1;
elseif get(handles.positiveIntensityRadio, 'Value')
    positive = 1;
elseif get(handles.negativeIntensityRadio, 'Value')
    positive = -1;
end
xjviewpath = fileparts(which('xjview'));
maskImageFile = fullfile(xjviewpath, 'mask.img');
[threshold, pvalue] = fdr(handles.imageFileName{1}, q, positive, maskImageFile);
set(handles.pValueEdit,'string',pvalue);
CallBack_pValueEdit(handles.pValueEdit, eventdata);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% model comparison
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_modelComparison(hObject, eventdata)
answer = getVariable({'y','x (full model)', 'x (reduced model)'});

if isempty(answer)
    return;
end

if ~isempty(answer{1})
    y = evalin('base',answer{1});
else
    return;
end
if ~isempty(answer{2})
    xf = evalin('base',answer{2}); % full model
else
    return;
end
if ~isempty(answer{3})
    xr = evalin('base',answer{3}); % reduced model
else
    xr = [];
end

% make vector column vector
[r, c] = size(y);
if r==1; y = y'; end
[r, c] = size(xf);
if r<c; xf = xf'; end
[r, c] = size(xr);
if r<c; xr = xr'; end

format long;
disp('------------------------------------------')
TotalVariance = var(y)
n = size(y, 1);
dfTotal = n - 1

disp('------------------------------------------')
disp('Full Model:');
beta = linearregression(y,xf)
predictedy = [xf ones(size(xf, 1),1)] * beta;
residual = y - predictedy;
r2 = 1 - (var(residual) / var(y))
ResidualVarianceFullModel =  var(residual)
ModelVarianceFullModel =  var(predictedy);
dfResidualFull = n - size(xf, 2) - 1

if isempty(xr)
    return;
end

disp('------------------------------------------')
disp('Reduced Model:');
beta = linearregression(y,xr)
predictedy = [xr ones(size(xr, 1),1)] * beta;
residual = y - predictedy;
r2 = 1 - (var(residual) / var(y))
ResidualVarianceReducedModel =  var(residual)
ModelVarianceReducedModel =  var(predictedy);
dfResidualReduced = n - size(xr, 2) - 1

disp('------------------------------------------')
disp('Model Comparison (is the full model significantly better than the reduced model?):');
dfDenominator = dfResidualFull;
dfNumerator = dfResidualReduced - dfResidualFull;
F = (ResidualVarianceReducedModel - ResidualVarianceFullModel)/dfNumerator / (ResidualVarianceFullModel / dfDenominator);
disp(['f(' num2str(dfNumerator) ',' num2str(dfDenominator) ')=' num2str(F)])
pvalue = 1-spm_Fcdf(F, [dfNumerator dfDenominator])
% 
% F = ModelVarianceFullModel/size(xf,2) / (ModelVarianceReducedModel/size(xr,2));
% disp(['f(' num2str(size(xf,2)) ',' num2str(size(xr,2)) ')=' num2str(F)])
% pvalue = 1-spm_Fcdf(F, [size(xf,2) size(xr,2)])

format short;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% behavior analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_behaviorAnalysis(hObject, eventdata)
global TR_;
answer = getVariable({'event data file','correlator', 'head movement files'});

if isempty(answer)
    return;
end

if ~isempty(answer{1})
    Pe = evalin('base',answer{1}); % event files
    numsubj = size(Pe, 1);
else
    Pe = [];
end
if ~isempty(answer{2})
    Pc = evalin('base',answer{2}); % correlator files
    numsubj = size(Pc, 1);
else
    Pc = [];
end
if ~isempty(answer{3})
    P = evalin('base',answer{3}); % headmovment files
    numsubj = size(P, 1);
else
    P = [];
end

if isempty(Pe) && isempty(Pc) && ~isempty(P)
    colors = 'rgbcmk';
    figure;
    Maximize(gcf);
    for ii=1:size(P,1)

        kk = mod(ii,8);
        if kk == 0
            kk = 8;
        end
        subplot(4,2,kk)

        tmp = load(deblank(P(ii,:)));
        tmp2 = fieldnames(tmp);
        hm = eval(['tmp.' tmp2{1}]);
        hmcell = struct2cell(hm);
        hmname = fieldnames(hm);
        for jj=1:length(hmname)
            plot(hmcell{jj}, colors(jj));
            hold on
        end
        legend(hmname)
        title(['subject ' num2str(ii)])

        if mod(ii,8) == 0
            figure;
            Maximize(gcf);
        end    

    end
    return
end

colors=[1 0 0;
    0 0 1;
    0 0 0;
    1 0 1;
    0 1 1;
    0 1 0;
    1 1 0;
    1 0.5 0.5;
    .5 .5 1;
    .5 .5 .5;
    1 .5 0;
    1 0 .5;
    0 1 .5;
    .5 1 0;
    0 .5 1;
    .5 0 1];
colors = repmat(colors, 10, 1);

% allow user to select which subject to plot
prompt=['I find you have ' num2str(numsubj) ' subjects.  Please let me know which subject(s) you want me to use to generate the plots.'];
name='Which subject(s) to plot?';
numlines=1;
defaultanswer={['[1:'  num2str(numsubj) ']']};
 
whichtoplot=inputdlg(prompt,name,numlines,defaultanswer);
if isempty(whichtoplot)
    return;
end
whichtoplot = eval(whichtoplot{1});


%% event
if ~isempty(Pe)
    figure;
    Maximize(gcf);
    for ii=whichtoplot
        kk = mod(ii,4);
        if kk == 0
            kk = 4;
        end
        subplot(4,2,2*kk-1)



        tmp = load(deblank(Pe(ii,:)));
        tmp2 = fieldnames(tmp);
        hm = eval(['tmp.' tmp2{1}]);
        hmcell = struct2cell(hm);
        hmname = fieldnames(hm);
        for jj=1:length(hmname)
            toplot1 = [];
            toplot2 = [];
            for mm = 1:length(hmcell{jj})
                toplot1 = [toplot1 hmcell{jj}(mm) hmcell{jj}(mm) nan];
                toplot2 = [toplot2 0 1 nan];
            end
            plot(toplot1, toplot2, 'color', colors(jj,:));
            ylim([0 4])
            hold on
        end
        legend(hmname)
        title(['subject ' num2str(ii)])
        xlabel('time in s');

        %% prepare interval data for later plot
        alltimes{ii} = [];
        for jj=1:length(hmname)
            alltimes{ii} = union(alltimes{ii}, hmcell{jj});
        end

        subplot(4,2,2*kk)
        for jj=1:length(hmname)
            roundhmcell{jj} = round(hmcell{jj}*1);
            signal = zeros(1, max(roundhmcell{jj})+1);
            signal(roundhmcell{jj}+1) = 1;
            signal = signal - mean(signal);
            y = fft(signal);
            plotlength = round(length(y))/1;
            plot([1:plotlength]/plotlength, abs(y(1:plotlength)).^2, 'color', colors(jj,:));
            hold on
        end
        legend(hmname)
        xlabel('frequency');
        ylabel('fft power');
        title(['subject ' num2str(ii)])

        subplot(4,2,2*kk-1)
        if ~isempty(P) % if there is head movement data
            tmp = load(deblank(P(ii,:)));
            tmp2 = fieldnames(tmp);
            hm = eval(['tmp.' tmp2{1}]);
            hmcell = struct2cell(hm);
            hmname = fieldnames(hm);
            for jj=1:length(hmname)
                plot(TR_*[0:length(hmcell{jj})-1],hmcell{jj},'color', [.5 .5 .5]);
                hold on
            end        
        end

        if mod(ii,4) == 0 && ii~=size(Pe, 1)
            figure;
            Maximize(gcf);
        end     
    end

    %% now I plot correlogram of intervals
    figure;
    kk = 1;
    for ii=whichtoplot
        alltimes{ii} = sort(alltimes{ii});
        intervals = diff(alltimes{ii});
        corrlag(intervals,intervals,[max(-10, -length(intervals)+2) : min(10, length(intervals)-2)], colors(ii,:));
        hold on;
        legendname{kk} = ['sub ' num2str(ii)];
        kk = kk+1;
    end
    legend(legendname)
    xlabel('lag');
    ylabel('correlation coefficient');
    title('correlogram of intervals')
end

%% correlator

if ~isempty(Pc)
    
	tmp = load(deblank(Pc(1,:)));
	tmp2 = fieldnames(tmp);
	hm = eval(['tmp.' tmp2{1}]);

    C = hm;
    cenames = fields(C);
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            eval(['C.' cenames{jj} '.' correlatornames{jj}{kk} '= [];']);
        end
    end
    
    %list all plot names in command window for selection
    promt{1} = 'Please select which two correlators do you want to plot. The first one will be x and the second will be y.  The two have to be of the same length.';
    promt{2} = 'index   event    correlator';
    count = 1;
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            promt{count+2} = sprintf('%d         %s         %s', count,  cenames{jj}, correlatornames{jj}{kk});
            count = count + 1;
        end
    end
    promt{count+2} = '';

    if count == 1
        return
    end
    promt = char(promt);
    name='Please select which two correlator do you want to plot. The first one will be x and the second will be y.';
    numlines=1;
    if count == 1
        defaultanswer={['You don''t have more than one correlators. Click Cancel']};
    else
        defaultanswer={['[1 '  num2str(count-1) ']']};
    end

    whichplot=inputdlg(promt,name,numlines,defaultanswer);
    if isempty(whichplot)
        return;
    end    
    whichplot = eval(whichplot{1});
    
    figure;
    Maximize(gcf);
    for ii=whichtoplot
        tmp = load(deblank(Pc(ii,:)));
    	tmp2 = fieldnames(tmp);
        hm = eval(['tmp.' tmp2{1}]);
        count = 1;
        for jj=1:length(cenames)
            eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
            for kk=1:length(correlatornames{jj})
                tmp2 = eval(['hm.' cenames{jj} '.' correlatornames{jj}{kk} ';']);
                if size(tmp2,1) == 1
                    tmp2 = tmp2';
                end
                eval(['C.' cenames{jj} '.' correlatornames{jj}{kk} '= [C.' cenames{jj} '.' correlatornames{jj}{kk} '; tmp2];']);
                
                if whichplot(1) == count
                    x = tmp2;
                    xlabeltext = [correlatornames{jj}{kk} ' at ' cenames{jj}];
                end
                if whichplot(2) == count
                    y = tmp2;
                    ylabeltext = [correlatornames{jj}{kk} ' at ' cenames{jj}];
                end
                count = count + 1;
            end
        end
        
        kk = mod(ii,16);
        if kk == 0
            kk = 16;
        end
        subplot(4,4,kk)
        
        plot2(x,y, 'b', 1);
        xlabel(xlabeltext)
        ylabel(ylabeltext)
        title(['subject ' num2str(ii)])    
        try
            % convert to colum vector
            if size(x,1)==1
                x = x';
            end
            if size(y,1)==1
                y = y';
            end

            % remove NaN
            tmpposx = find(isnan(x) | isinf(x));
            tmpposy = find(isnan(y) | isinf(y));
            x([tmpposx; tmpposy]) = [];
            y([tmpposx; tmpposy]) = [];

            b = linearregression(y,x);

            totalvar = var(y);
            residvar = var(y - x*b(1) - b(2));
            modelvar = totalvar - residvar;
            dfm = 1;
            dft = length(x) - 1;
            dfr = dft - dfm;
            mmodelvar = modelvar/dfm;
            mresidvar = residvar/dfr;
            F = mmodelvar/mresidvar;
            pvalue = 1-spm_Fcdf(F, [dfm dfr]);

            title(sprintf('sub %d, pValue=%g', ii, pvalue))
            xx=[min(x):(max(x)-min(x))/10:max(x)];
            yy=b(1) * xx + b(2);
            hold on;
            plot(xx,yy,'g');
        end
        
        if mod(ii,16) == 0 && ii~=size(Pc, 1)
            figure;
            Maximize(gcf);
        end    
    end

    % plot the average
    count = 1;
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            if whichplot(1) == count
                eval(['x = C.' cenames{jj} '.' correlatornames{jj}{kk} ';']);
                xlabeltext = [correlatornames{jj}{kk} ' at ' cenames{jj}];
            end
            if whichplot(2) == count
                eval(['y = C.' cenames{jj} '.' correlatornames{jj}{kk} ';']);
                ylabeltext = [correlatornames{jj}{kk} ' at ' cenames{jj}];
            end
            count = count + 1;
        end
    end
    
    figure
    plot2(x,y, 'b', 1);
    xlabel(xlabeltext)
    ylabel(ylabeltext)
    try
        % convert to colum vector
        if size(x,1)==1
            x = x';
        end
        if size(y,1)==1
            y = y';
        end

        % remove NaN
        tmpposx = find(isnan(x) | isinf(x));
        tmpposy = find(isnan(y) | isinf(y));
        x([tmpposx; tmpposy]) = [];
        y([tmpposx; tmpposy]) = [];

        b = linearregression(y,x);

        totalvar = var(y);
        residvar = var(y - x*b(1) - b(2));
        modelvar = totalvar - residvar;
        dfm = 1;
        dft = length(x) - 1;
        dfr = dft - dfm;
        mmodelvar = modelvar/dfm;
        mresidvar = residvar/dfr;
        F = mmodelvar/mresidvar;
        pvalue = 1-spm_Fcdf(F, [dfm dfr]);

        title(sprintf('pValue=%g; slope=%g; constant=%g', pvalue, b(1), b(2)))
        xx=[min(x):(max(x)-min(x))/10:max(x)];
        yy=b(1) * xx + b(2);
        hold on;
        plot(xx,yy,'g');
    end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% head movement analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_headMovementAnalysis(hObject, eventdata)
answer = getVariable({'head movement files'});

if isempty(answer)
    return;
end
if isempty(answer{1})
    return;
end

P = evalin('base',answer{1}); % headmovment files

colors = 'rgbcmk';
figure;
Maximize(gcf);
for ii=1:size(P,1)
    
    kk = mod(ii,8);
    if kk == 0
        kk = 8;
    end
    subplot(4,2,kk)
    
    tmp = load(deblank(P(ii,:)));
    tmp2 = fieldnames(tmp);
    hm = eval(['tmp.' tmp2{1}]);
    hmcell = struct2cell(hm);
    hmname = fieldnames(hm);
    for jj=1:length(hmname)
        plot(hmcell{jj}, colors(jj));
        hold on
    end
    legend(hmname)
    title(['subject ' num2str(ii)])

    if mod(ii,8) == 0
        figure;
        Maximize(gcf);
    end    
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% contrast
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_contrast(hObject, eventdata)

answer = getVariable({'beta files', 'contrast vector'});

if isempty(answer)
    return;
end
if isempty(answer{1}) | isempty(answer{2}) 
    return;
end

[filename, pathname] = uiputfile('*.img', 'Save t-test image as', '');
if isequal(filename,0) | isequal(pathname,0)
   return
else
   thisfilename = fullfile(pathname, filename);
end


P1 = evalin('base',answer{1}); % beta image files
c = evalin('base',answer{2}); % contrast vector

if iscell(P1)
    [];
else
    disp('format wrong');
    return;
end


h = waitbar(0,'Please wait...');
V=spm_vol(deblank(P1{1}(1,:)));
M1={spm_read_vols(V)};
contr = zeros(size(M1{1},1), size(M1{1},2), size(M1{1},3), size(P1{1},1));

for ii=1:size(P1{1},1)
    for jj=1:length(P1)
        if c(jj)~=0
            V=spm_vol(deblank(P1{jj}(ii,:)));
            M1{jj}=spm_read_vols(V);
            contr(:,:,:,ii) = contr(:,:,:,ii)+c(jj)*double(M1{jj});
        end
    end

    waitbar(ii/size(P1{1},1),h,['finished ' num2str(ii) ' of ' num2str(size(P1{1},1))]);
end
close(h);


df = size(P1{1},1)-1;

T = squeeze(t(permute(contr,[4 1 2 3])));   % return a 41*48*35 matrix
xjviewpath = fileparts(which('xjview'));
V=spm_vol(fullfile(xjviewpath, 'mask.img'));
m=spm_read_vols(V);
T = T.*m;

targetfilename = thisfilename;
global M_;
global DIM_;
V.mat = M_;
V.dim = [DIM_(1) DIM_(2) DIM_(3) 16];
V.fname = targetfilename;
V.descrip = ['SPM{T_[' num2str(df) '.0]} ' num2str(c)];
V = spm_write_vol(V,T);

handles = guidata(hObject);
[tmp,tmp,ext] = fileparts(thisfilename);
if isempty(ext)
    ext = '.img';
elseif isequal(lower(ext), '.img')
    ext = '';
end
CallBack_loadImagePush(handles.loadImagePush, [], {[thisfilename ext]});

return;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% SPM process (use SPM processing)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_SPMProcess(hObject, eventdata, singleOrCluster)

answer = getVariable({'Image files', 'Event time data', 'Modulator','Other regressor (e.g. headmovement)'});

if isempty(answer)
    return;
end
if isempty(answer{1}) | isempty(answer{2})
    return;
end

P = evalin('base',answer{1}); % subject directories
E = evalin('base',answer{2}); % event time
try
    Modulator = evalin('base',answer{3}); % modulator
catch
    Modulator = [];
end
try
    headmovement = evalin('base',answer{4}); % other regressor such as headmovement
catch
    headmovement = [];
end

if nargin < 3
    singleOrCluster = 'single';
end


prompt='I will create a folder under each subject''s directory and put the SPM output files (including beta images, SPM.mat etc) for that subject into that folder. Please give a name to the folder (example: FixedEffect). No space or strange characters allowed.';
name='Give a name to SPM output folder';
numlines=1;
defaultanswer={'FixedEffect'};
 
wheretosave=inputdlg(prompt,name,numlines,defaultanswer);

if isempty(wheretosave)
    return
else
    wheretosave = wheretosave{1};
end

currentdir = pwd;

if strcmp(singleOrCluster, 'cluster')    
    handles = guidata(hObject);
    set(handles.infoTextBox, 'string', {'This will do SPM GLM estimation using cluster.', 'Check out xjviewtmp directory'});
    
    mkdir('xjviewtmp');
    cd('xjviewtmp');
    system('rm *');

    fid = fopen('processALL.pbs','wt');
    fprintf(fid,'#!/bin/bash\n\n');
    fprintf(fid,'for ((s=1;s<=%d;s++))\n', size(P,1));
    fprintf(fid,'do\n\tsleep 3.5\n\tqsub %s/process.pbs -v "s=$s" -N SPM_$s\ndone\n', pwd);
    fclose(fid);

    save P P;
    save E E;
    save Modulator Modulator;
    save headmovement headmovement;
    fid = fopen('process.pbs','wt');
    fprintf(fid,'#!/bin/bash\n\n');
    fprintf(fid,'#PBS -l nodes=1:ppn=1,walltime=04:00:00\n');
    fprintf(fid,'#PBS -N process\n');
    fprintf(fid,'#PBS -e %s\n', pwd);
    fprintf(fid,'#PBS -o %s\n\n', pwd);
    fprintf(fid,'sleep 3.5\n\n');
    fprintf(fid,'PATH=/usr/local/bin:$PATH\n\n');
    fprintf(fid,'matlab -nojvm -nodisplay \\\n');
    fprintf(fid,'-logfile "%s/process.$s.$PBS_JOBID.log" \\\n', pwd);


    fprintf(fid,['-r "addpath ~/gang/xjview;'...
        'wheretosave = ''%s'';'...
        'cd %s; '...
        'load P; '...
        'load E; '...
        'load Modulator;'... 
        'load headmovement;'...
        'tmp=load(deblank(P($s,:)));'...
        'fd=fields(tmp);'...
        'eval([''Pi=tmp.'' fd{1} '';'']);'...
        'N = size(Pi,1);'...
        'tmp=load(deblank(E($s,:))); '...
        'fd=fields(tmp); '...
        'eval([''EE=tmp.'' fd{1} '';'']); '...
        'thisModulator = [];'...
        'if ~isempty(Modulator);'...
            'tmp = load(deblank(Modulator($s,:)));  '...
            'fd=fields(tmp); '...
            'eval([''thisModulator=tmp.'' fd{1} '';'']);'...
        'end;'...
        '[otherregressor{1}.C, otherregressor{1}.name]=struct2matrix(N, EE, thisModulator);'...
        'if ~isempty(headmovement);'...
            'tmp = load(deblank(headmovement($s,:)));  '...
            'fd=fields(tmp); '...
            'eval([''thisheadmovement=tmp.'' fd{1} '';'']);'...
            'tmphd = struct2cell(thisheadmovement);'...
            'tmpname = fieldnames(thisheadmovement);'...
            'for kk=1:length(tmphd);'...
                '[r,c] = size(tmphd{kk});'...
                'if r==1; tmphd{kk} = tmphd{kk}''; end;'...
                'otherregressor{1}.C = [otherregressor{1}.C tmphd{kk}];'...
                'otherregressor{1}.name = [otherregressor{1}.name {tmpname{kk}}];'...
            'end;'...
        'end;'...
        'subDir = fileparts(deblank(Pi(1,:))); '...
        'for kk=1:length(subDir); '...
            'if subDir(end-kk+1)==filesep;'...
                'pos=kk;'...
                'break;'...
            'end;'...
        'end; '...
        'subDir(end-pos+1:end)=[]; '...
        'cd(subDir); '...
        'try;' ...
            'rmdir(wheretosave, ''s'');'...
        'end;'...
        'mkdir(wheretosave);'...
        'cd(wheretosave);'...
        'spm_defaults;'...
        'disp(subDir);'...
        'cuixuprocess(N, [0], {}, {}, {}, [], otherregressor, Pi, 16);'...
        'disp(''done!'');'...
        'exit;"'], wheretosave, pwd);
    
    fclose(fid);

    system(['ssh cluster.hnl.bcm.tmc.edu  PATH=$PATH:/usr/local/pbs/i686/bin bash ' pwd '/processALL.pbs']);
    cd ..
elseif  strcmp(singleOrCluster, 'single')
    handles = guidata(hObject);
    set(handles.infoTextBox, 'string', {'This will do SPM GLM estimation.'});
    for ii=1:size(P,1)
        tmp=load(deblank(P(ii,:)));
        fd=fields(tmp);
        eval(['Pi=tmp.' fd{1} ';']);
        N = size(Pi,1);
        tmp=load(deblank(E(ii,:))); 
        fd=fields(tmp); 
        eval(['EE=tmp.' fd{1} ';']); 
        thisModulator = [];
        if ~isempty(Modulator)
            tmp = load(deblank(Modulator(ii,:)));  
            fd=fields(tmp); 
            eval(['thisModulator=tmp.' fd{1} ';']);
        end
        [otherregressor{1}.C, otherregressor{1}.name]=struct2matrix(N, EE, thisModulator);
        if ~isempty(headmovement)
            tmp = load(deblank(headmovement(ii,:)));  
            fd=fields(tmp); 
            eval(['thisheadmovement=tmp.' fd{1} ';']);
            
            tmphd = struct2cell(thisheadmovement);
            tmpname = fieldnames(thisheadmovement);
            for kk=1:length(tmphd)
                [r,c] = size(tmphd{kk});
                if r==1; tmphd{kk} = tmphd{kk}'; end
                otherregressor{1}.C = [otherregressor{1}.C tmphd{kk}];
                otherregressor{1}.name = [otherregressor{1}.name {tmpname{kk}}];
            end
        end
        subDir = fileparts(deblank(Pi(1,:))); 
        for kk=1:length(subDir); 
            if subDir(end-kk+1)==filesep;
                pos=kk;
                break;
            end;
        end; 
        subDir(end-pos+1:end)=[]; 
        cd(subDir); 
        try
            rmdir(wheretosave, 's');
        end
        mkdir(wheretosave);
        cd(wheretosave);

        spm_defaults;
        disp(subDir);
        cuixuprocess(N, [0], {}, {}, {}, [], otherregressor, Pi, 16);
        
        tmp = get(handles.infoTextBox, 'string');
        set(handles.infoTextBox, 'string', [tmp; {['SPM Processing subject ' num2str(ii) ' / ' num2str(size(P,1)) ' ...']}]); 
        disp(['----finished ' num2str(ii) ' out of ' num2str(size(P,1)) ' subjects----']);
    end
	tmp = get(handles.infoTextBox, 'string');
	set(handles.infoTextBox, 'string', [tmp; {'Done!'}]); 
end

cd(currentdir)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% GLMPeak (only do regression on peak BOLD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_GLMPeak(hObject, eventdata, singleOrCluster)

answer = getVariable({'Image data', 'Event time data', 'Modulator'});

if isempty(answer)
    return;
end
if isempty(answer{1}) | isempty(answer{2}) | isempty(answer{3})
    return;
end

directoryname = uigetdir(pwd, 'Pick a Directory to Save Your Beta Files...');
if isequal(directoryname, 0)
    return;
end

peakpoint = inputdlg({'What points in time do you considered as peak BOLD?'}, 'Peak point',1,{'[4 6]'});
if isempty(peakpoint); return; end
peakpointstr = peakpoint{1};
peakpoint = str2num(peakpoint{1});
    
    
P = evalin('base',answer{1}); % subject directories
E = evalin('base',answer{2}); % event time
Modulator = evalin('base',answer{3}); % modulator


if nargin < 3
    singleOrCluster = 'single';
end

if strcmp(singleOrCluster, 'cluster')
    mkdir('xjviewtmp');
    cd('xjviewtmp');
    system('rm *');

    fid = fopen('processALL.pbs','wt');
    fprintf(fid,'#!/bin/bash\n\n');
    fprintf(fid,'for ((s=1;s<=%d;s++))\n', size(P,1));
    fprintf(fid,'do\n\tsleep 3.5\n\tqsub %s/process.pbs -v "s=$s" -N pro_$s\ndone\n', pwd);
    fclose(fid);

    save P P;
    save E E;
    save Modulator Modulator;
    fid = fopen('process.pbs','wt');
    fprintf(fid,'#!/bin/bash\n\n');
    fprintf(fid,'#PBS -l nodes=1:ppn=1,walltime=08:00:00\n');
    fprintf(fid,'#PBS -N process\n');
    fprintf(fid,'#PBS -e %s\n', pwd);
    fprintf(fid,'#PBS -o %s\n\n', pwd);
    fprintf(fid,'sleep 3.5\n\n');
    fprintf(fid,'PATH=/usr/local/bin:$PATH\n\n');
    fprintf(fid,'matlab -nojvm -nodisplay \\\n');
    fprintf(fid,'-logfile "%s/process.$s.$PBS_JOBID.log" \\\n', pwd);

    fprintf(fid,['-r "addpath ~/gang/xjview; '...
            'cd %s; '...
            'load P; '...
            'load E; '...
            'load Modulator; '...
            'tmp=load(deblank(P($s,:))); '...
            'fd=fields(tmp); '...
            'eval([''SS=tmp.'' fd{1} '';'']); '...
            'tmp=load(deblank(E($s,:))); '...
            'fd=fields(tmp); '...
            'eval([''EE=tmp.'' fd{1} '';'']); '...
            'tmp = load(deblank(Modulator($s,:)));  fd=fields(tmp);  eval([''MM=tmp.'' fd{1} '';'']);'...
                    'names = fieldnames(EE);'...
            'for jj=1:length(names);'...
                'eval([''MM.'' names{jj} ''.time = EE.'' names{jj} '';'']);'...
            'end;'...
            'xjviewpath = fileparts(which(''xjview''));'...
            'cuixuGLMpeak(MM, SS, ''%s'', [''subject'' num2str($s)], fullfile(xjviewpath, ''mask.img''), %s, 100, 2, fullfile(xjviewpath, ''templateFile.img''));'...
            'exit;"'], pwd, directoryname, peakpointstr);

    fclose(fid);

    system(['ssh cluster.hnl.bcm.tmc.edu  PATH=$PATH:/usr/local/pbs/i686/bin bash ' pwd '/processALL.pbs']);
    cd ..
elseif strcmp(singleOrCluster, 'single')
    xjviewpath = fileparts(which('xjview'));
    handles = guidata(hObject);
    set(handles.infoTextBox, 'string', {'This will do simple GLM estimation on peak BOLD.'});
    for ii = 1:size(P,1)
        tmp = get(handles.infoTextBox, 'string');
        set(handles.infoTextBox, 'string', [tmp; {['Processing subject ' num2str(ii) ' / ' num2str(size(P,1)) ' ...']}]); 
        tmp=load(deblank(P(ii,:)));
        fd=fields(tmp);
        eval(['SS=tmp.' fd{1} ';']);
        tmp=load(deblank(E(ii,:)));
        fd=fields(tmp);
        eval(['EE=tmp.' fd{1} ';']);
        tmp=load(deblank(Modulator(ii,:)));
        fd=fields(tmp);
        eval(['MM=tmp.' fd{1} ';']);
        % convert Modulator and E to be recognized by cuixuGLMpeak
        names = fieldnames(EE);
        for jj=1:length(names)
            eval(['MM.' names{jj} '.time = EE.' names{jj} ';']);
        end
        cuixuGLMpeak(MM, SS, directoryname, ['subject' num2str(ii)], fullfile(xjviewpath, 'mask.img'), peakpoint, 100, 2, fullfile(xjviewpath, 'templateFile.img'));
        disp(['----finished ' num2str(ii) ' out of ' num2str(size(P,1)) ' subjects----']);
    end
    tmp = get(handles.infoTextBox, 'string');
    set(handles.infoTextBox, 'string', [tmp; {'Done!'}]); 
end
  
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% process (use my own process code, cuixuAGLM, estimation on detrended relative BOLD)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_process(hObject, eventdata, singleOrCluster)

answer = getVariable({'Image data', 'Event time data', 'Modulator'});

if isempty(answer)
    return;
end
if isempty(answer{1}) | isempty(answer{2})
    return;
end

directoryname = uigetdir(pwd, 'Pick a Directory to Save Your Beta Files...');
if isequal(directoryname, 0)
    return;
end

P = evalin('base',answer{1}); % subject directories
E = evalin('base',answer{2}); % event time
try
    Modulator = evalin('base',answer{3}); % modulator
catch
    Modulator = [];
end


if nargin < 3
    singleOrCluster = 'single';
end

if strcmp(singleOrCluster, 'cluster')
    mkdir('xjviewtmp');
    cd('xjviewtmp');
    system('rm *');

    fid = fopen('processALL.pbs','wt');
    fprintf(fid,'#!/bin/bash\n\n');
    fprintf(fid,'for ((s=1;s<=%d;s++))\n', size(P,1));
    fprintf(fid,'do\n\tsleep 3.5\n\tqsub %s/process.pbs -v "s=$s" -N pro_$s\ndone\n', pwd);
    fclose(fid);

    save P P;
    save E E;
    save Modulator Modulator;
    fid = fopen('process.pbs','wt');
    fprintf(fid,'#!/bin/bash\n\n');
    fprintf(fid,'#PBS -l nodes=1:ppn=1,walltime=08:00:00\n');
    fprintf(fid,'#PBS -N process\n');
    fprintf(fid,'#PBS -e %s\n', pwd);
    fprintf(fid,'#PBS -o %s\n\n', pwd);
    fprintf(fid,'sleep 3.5\n\n');
    fprintf(fid,'PATH=/usr/local/bin:$PATH\n\n');
    fprintf(fid,'matlab -nojvm -nodisplay \\\n');
    fprintf(fid,'-logfile "%s/process.$s.$PBS_JOBID.log" \\\n', pwd);

    fprintf(fid,['-r "addpath ~/gang/xjview; '...
            'cd %s; '...
            'load P; '...
            'load E; '...
            'load Modulator; '...
            'tmp=load(deblank(P($s,:))); '...
            'fd=fields(tmp); '...
            'eval([''M=tmp.'' fd{1} '';'']); '...
            'if isnumeric(M); M=double(M); elseif isstr(M); V=spm_vol(M); M =spm_read_vols(V); else; error(''I don not understand the fileformat''); end;' ...
            'tmp=load(deblank(E($s,:))); '...
            'fd=fields(tmp); '...
            'eval([''EE=tmp.'' fd{1} '';'']); '...
            'if ~isempty(Modulator); tmp = load(deblank(Modulator($s,:)));  fd=fields(tmp);  eval([''thisModulator=tmp.'' fd{1} '';'']); else; thisModulator=[]; end;'...
            'M = permute(M,[4 1 2 3]); '...
            '[tmp,tmp,M]=mAveDetrend(M,100); '...
            'M = permute(M,[2 3 4 1]); '...
            '[a,b,c]=fileparts(deblank(P($s,:)));'...
            'cuixuAGLM(M,EE,{''%s'', b}, thisModulator); '...
            'exit;"'], pwd, directoryname);

    fclose(fid);

    system(['ssh cluster.hnl.bcm.tmc.edu  PATH=$PATH:/usr/local/pbs/i686/bin bash ' pwd '/processALL.pbs']);
    cd ..
elseif strcmp(singleOrCluster, 'single')
    handles = guidata(hObject);
    set(handles.infoTextBox, 'string', {'This will do simple GLM estimation.'});
    for ii=1:size(P,1)
        tmp = get(handles.infoTextBox, 'string');
        set(handles.infoTextBox, 'string', [tmp; {['Processing subject ' num2str(ii) ' / ' num2str(size(P,1)) ' ...']}]); 
        tmp=load(deblank(P(ii,:)));
        fd=fields(tmp);
        eval(['M=tmp.' fd{1} ';']);
        if isnumeric(M); M=double(M); elseif isstr(M); V=spm_vol(M); M =spm_read_vols(V); else; error('I don''t understand the fileformat'); end;
        tmp=load(deblank(E(ii,:)));
        fd=fields(tmp);
        eval(['EE=tmp.' fd{1} ';']);
        if ~isempty(Modulator); tmp = load(deblank(Modulator(ii,:)));  fd=fields(tmp);  eval(['thisModulator=tmp.' fd{1} ';']); else; thisModulator=[]; end;
        M = permute(M,[4 1 2 3]);
        [tmp,tmp,M]=mAveDetrend(M,100);
        M = permute(M,[2 3 4 1]);
        [a,b,c]=fileparts(deblank(P(ii,:)));
        cuixuAGLM(M,EE,{directoryname, b}, thisModulator);
        disp(['----finished ' num2str(ii) ' out of ' num2str(size(P,1)) ' subjects----']);
    end
	tmp = get(handles.infoTextBox, 'string');
	set(handles.infoTextBox, 'string', [tmp; {'Done!'}]); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% convert new data format to old
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_new2old(hObject, eventdata)
answer = getVariable({'subjects directories'});
if isempty(answer)
    return;
end
if isempty(answer{1})
    return;
end
P = evalin('base',answer{1}); % subject directories
new2old(P);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% pre-process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_preprocess(hObject, eventdata, singleOrCluster)
answer = getVariable({'subjects directories'});
if isempty(answer)
    return;
end
if isempty(answer{1})
    return;
end
P = evalin('base',answer{1}); % subject directories

if nargin < 3
    singleOrCluster = 'single';
end

if strcmp(singleOrCluster, 'cluster')
    mkdir('xjviewtmp');
    cd('xjviewtmp');
    system('rm *');

    fid = fopen('preprocessALL.pbs','wt');
    fprintf(fid,'#!/bin/bash\n\n');
    fprintf(fid,'for ((s=1;s<=%d;s++))\n', size(P,1));
    fprintf(fid,'do\n\tsleep 3.5\n\tqsub %s/preprocess.pbs -v "s=$s" -N pre_$s\ndone\n', pwd);
    fclose(fid);

    save P P;
    fid = fopen('preprocess.pbs','wt');
    fprintf(fid,'#!/bin/bash\n\n');
    fprintf(fid,'#PBS -l nodes=1:ppn=1,walltime=08:00:00\n');
    fprintf(fid,'#PBS -N preprocess\n');
    fprintf(fid,'#PBS -e %s\n', pwd);
    fprintf(fid,'#PBS -o %s\n\n', pwd);
    fprintf(fid,'sleep 3.5\n\n');
    fprintf(fid,'PATH=/usr/local/bin:$PATH\n\n');
    fprintf(fid,'matlab -nojvm -nodisplay \\\n');
    fprintf(fid,'-logfile "%s/preprocess.$s.$PBS_JOBID.log" \\\n', pwd);
    fprintf(fid,['-r "addpath ~/gang/xjview; cd ' pwd '; load P; cuixuSmartPreprocess(P, $s);exit;"']);
    fclose(fid);

    system(['ssh cluster.hnl.bcm.tmc.edu  PATH=$PATH:/usr/local/pbs/i686/bin bash ' pwd '/preprocessALL.pbs']);
    cd ..
elseif strcmp(singleOrCluster, 'single')
    handles = guidata(hObject);
    set(handles.infoTextBox, 'string', {'What will be done in preprocessing is:',...
        '1. Converting DICOM files',...
        '2. Re-align',...
        '3. Coregister',...
        '4. Slice timing',...
        '5. Normalize',...
        '6. Smooth','',...
        'This function creates a directory ''preprocessing'' under your homeDir and put all preprocessing files there. It also saves the head movement matrix in ''headmovement.mat''. The intermediate files are removed',''});
    for ii=1:size(P,1)
        tmp = get(handles.infoTextBox, 'string');
        set(handles.infoTextBox, 'string', [tmp; {['Preprocessing subject ' num2str(ii) ' / ' num2str(size(P,1)) ' ...']}]); 
        cuixuSmartPreprocess(P, ii);
    end
	tmp = get(handles.infoTextBox, 'string');
	set(handles.infoTextBox, 'string', [tmp; {'Done!'}]); 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ROI: TimeSeries
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_timeSeries(hObject, eventdata)

global TR_;
handles = guidata(hObject);
if ~isfield(handles,'currentDisplayMNI') | isempty(cell2mat(handles.currentDisplayMNI'))
    %errordlg('No voxels found.','error');
	set(handles.infoTextBox, 'string', 'No voxels found.'); 
    beep
    return;
end

answer = getVariable({'Image data', 'Event time data', 'Correlator', 'Effects to be removed (e.g. headmovement)'});
%answer=inputdlg('Input the name of the variable (in the base workspace) which contains the image file name list, or the 4-D data matrix. If you wish to select the image files directly, click cancel.', 'Variable Name?', 1, {'P'});
if isempty(answer)
    return;
end
if isempty(answer{1})
    return;
end

searchContent = get(handles.searchContentEdit,'string');
if isempty(searchContent)
    searchContent = num2str(coord(1,3));
end
searchContent = ['ROI_' searchContent];
[filename, pathname] = uiputfile('*.mat', 'Save result to', searchContent);
if isequal(filename,0) | isequal(pathname,0)
   return
else
   thisfilename = fullfile(pathname, filename);
end


P = evalin('base',answer{1}); % data
try
    E = evalin('base',answer{2}); % event
catch
    E = [];
end
try
    C = evalin('base',answer{3}); % correlator
catch
    C = [];
end
try
    headmovement = evalin('base',answer{4}); % headmovement
catch
    headmovement = [];
end

mni = cell2mat(handles.currentDisplayMNI');
coord = mni2cor(mni, handles.M{1});


if ischar(P)
    h = waitbar(0,'Please wait...');
    for ii=1:size(P,1)
        tmp = load(deblank(P(ii,:)));
        tmp1 = fields(tmp);
        eval(['M = tmp.' tmp1{1} ';']);
        if isnumeric(M)
            M = double(M);
        end
        if ~isempty(E)
            tmp = load(deblank(E(ii,:)));
            tmp1 = fields(tmp);
            eval(['event{ii} = tmp.' tmp1{1} ';']);        
            if ~isstruct(event{ii})
                error('I need you specify the event in a structure, not a cell array!');
                return
            end
        else
            event = {};
        end
     
        if ~isempty(headmovement) % if headmovement is present (or any other effects to be removed), then get rid of it first
            H = [];
            tmp = load(deblank(headmovement(ii,:)));  
            fd=fields(tmp); 
            eval(['thisheadmovement=tmp.' fd{1} ';']);
            
            tmphd = struct2cell(thisheadmovement);
            for kk=1:length(tmphd)
                [r,c] = size(tmphd{kk});
                if r==1; tmphd{kk} = tmphd{kk}'; end
                H = [H  tmphd{kk}];
            end    
            for kk=1:size(coord,1)
                b = linearregression(squeeze(M(coord(kk,1), coord(kk,2), coord(kk,3), :)),H);
                M(coord(kk,1), coord(kk,2), coord(kk,3), :) = squeeze(M(coord(kk,1), coord(kk,2), coord(kk,3), :)) - H*b(1:end-1);
            end
        end
        
        [eventResponse{ii}, time{ii}, wholeOriginal{ii}, wholeAbsolute{ii}, wholeBaseline{ii}, removedEvent{ii}, remainedEvent{ii}] = cuixuBOLDretrieve(M, coord, event{ii}, 50, -20, 100, 0, TR_);

        if ~isempty(C)
            tmp = load(deblank(C(ii,:)));
            tmp1 = fields(tmp);
            eval(['correlator{ii} = tmp.' tmp1{1} ';']);        
            if ~isstruct(correlator{ii})
                error('I need you specify the correlator in a structure, not a cell array!');
                return
            end

            names = fields(correlator{ii});
            for kk=1:length(names)
                eval(['names2 = fields(correlator{ii}.' names{kk} ');']);
                for mm=1:length(names2)
                    eval(['correlator{ii}.' names{kk} '.' names2{mm} ' = correlator{ii}.' names{kk} '.' names2{mm} '(remainedEvent{ii}.' names{kk} ');']);;
                end
            end
        end
        waitbar(ii/size(P,1),h,['finished ' num2str(ii) ' of ' num2str(size(P,1))]);
        disp(['finished ' num2str(ii) ' of ' num2str(size(P,1))]);
    end
    close(h)    
elseif isnumeric(P) % never used!
    P = double(P);
	if ~isempty(E)
        tmp = load(deblank(E(ii,:)));
        tmp1 = fields(tmp);
        eval(['event{ii} = tmp.' tmp1{1} ';']);        
	else
        event = {};
	end
    [eventResponse{ii}, time{ii}, wholeOriginal{ii}, wholeAbsolute{ii}, wholeBaseline{ii}, removedEvent{ii}, remainedEvent{ii}] = cuixuBOLDretrieve(P, coord, event{ii}, 50, -20, 100, 0, TR_);
end

if ~exist('correlator')
    correlator = {};
end

mat = handles.M{1};
save(thisfilename, 'eventResponse', 'time', 'wholeOriginal', 'wholeAbsolute', 'wholeBaseline', 'correlator', 'coord', 'mni', 'removedEvent','remainedEvent', 'event', 'mat', 'P', 'headmovement');
CallBack_plotROI(hObject, eventdata, thisfilename);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% plot ROI
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_plotROI(hObject, eventdata, filename)
handles = guidata(hObject);

if ~exist('filename')
    if findstr('SPM2',spm('ver'))
        P = spm_get([0:1],'*.mat','Select ROI result file');
    elseif findstr('SPM5',spm('ver'))
        P = spm_select([0:1],'mat','Select ROI result file');
    end
else
    P = filename;
end

if isempty(P)
    return;
end
load(deblank(P));

if ~exist('mni')
    set(handles.infoTextBox, 'string', {'Oops! I donnt find the variable mni. Please open a proper mat file.'}); 
    return
end

delete(gcf)
xjview(mni);

% allow user to select which subject to plot
prompt=['I find you have ' num2str(length(eventResponse)) ' subjects.  Please let me know which subjects you want me to use to generate the plots.'];
name='Which subject(s) to plot?';
numlines=1;
defaultanswer={['[1:'  num2str(length(eventResponse)) ']']};
 
whichtoplot=inputdlg(prompt,name,numlines,defaultanswer);
if isempty(whichtoplot)
    return;
end
whichtoplot = eval(whichtoplot{1});


response = eventResponse{1};

eventisstruct = 0;
if isstruct(response)
    eventisstruct = 1;
	f = fieldnames(response);
    response = struct2cell(response); 
else
    f = num2cell([1:length(response)]);
    for ii=1:length(response)
        f{ii} = num2str(f{ii});
    end
end
for jj=1:length(response)
    response{jj} = [];
end

if exist('correlator') & ~isempty(correlator) & eventisstruct==1
    C = correlator{1};
    cenames = fields(C);
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            eval(['C.' cenames{jj} '.' correlatornames{jj}{kk} '= [];']);
        end
    end
    for ii=whichtoplot%1:length(correlator)
        for jj=1:length(cenames)
            eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
            for kk=1:length(correlatornames{jj})
                tmp2 = eval(['correlator{ii}.' cenames{jj} '.' correlatornames{jj}{kk} ';']);
                if size(tmp2,1) == 1
                    tmp2 = tmp2';
                end
                eval(['C.' cenames{jj} '.' correlatornames{jj}{kk} '= [C.' cenames{jj} '.' correlatornames{jj}{kk} '; tmp2];']);
            end
        end
    end
end

if eventisstruct == 1
    for ii=1:length(eventResponse)
        for kk=1:length(f)
            eval(['tmptmp{kk} = eventResponse{ii}.' f{kk} ';']);
        end
        eventResponse{ii} = tmptmp;
    end
end

for ii=whichtoplot%1:length(eventResponse)
    for jj=1:length(eventResponse{1})
        response{jj} = [response{jj}; eventResponse{ii}{jj}];
    end
end

if isempty(response)
    msgbox('Nothing to plot');
    return
end




v = version;
figure;
%colors='rbkmcgyrbkmcgyrbkmcgyrbkmcgyrbkmcgyrbkmcgyrbkmcgy';
colors=[1 0 0;
    0 0 1;
    0 0 0;
    1 0 1;
    0 1 1;
    0 1 0;
    1 1 0;
    1 0.5 0.5;
    .5 .5 1;
    .5 .5 .5;
    1 .5 0;
    1 0 .5;
    0 1 .5;
    .5 1 0;
    0 .5 1;
    .5 0 1];
colors = repmat(colors, 10, 1);
thislabel = {};
for ii=1:length(response)
    if v(1)=='6'
        ff = errorbar(time{1}, meannan(response{ii}), stdnan(response{ii})/sqrt(size(response{ii},1)));
        set(ff,'color',colors(ii, :));
        thislabel = [thislabel {''} f(ii)];
    elseif v(1)=='7'
        errorbar(time{1}, meannan(response{ii}), stdnan(response{ii})/sqrt(size(response{ii},1)), 'color', colors(ii, :));
        thislabel = [thislabel f(ii)];
    end
    hold on;
end
xlabel('time (s)');
ylabel('relative BOLD');
legend(thislabel)
hold off;

[row,col]  = size(response);
if row == 1; response = response'; end;
ResponseForPlot = cell2struct(response, f, 1);

if exist('correlator') & ~isempty(correlator)
    %list all plot names in command window for selection
    promt{1} = 'I will plot BOLD amplitude vs correlator. Here are all the correlators:';
    promt{2} = 'index   event    correlator';
    count = 1;
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            promt{count+2} = sprintf('%d         %s         %s', count,  cenames{jj}, correlatornames{jj}{kk});
            count = count + 1;
        end
    end
    promt{count+2} = 'Which plots do you want me to show?';

    promt = char(promt);
    name='Please select which plots to show';
    numlines=1;
    defaultanswer={['[1:'  num2str(count-1) ']']};

    whichplot=inputdlg(promt,name,numlines,defaultanswer);
    if isempty(whichplot)
        return;
    end    
    whichplot = eval(whichplot{1});
    
    count = 1;
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            if ~ismember(count, whichplot)
                count = count + 1;
                continue;
            end
            
            eval(['x = C.' cenames{jj} '.' correlatornames{jj}{kk} ';']);
            
            uniquex = unique(x);
            pos = find(isnan(uniquex));
            uniquex(pos) = [];
            
            if(length(uniquex)<10)  % only if the correlator is discrete, then plot BOLD vs time for each correlator
                figure;
                thislabel = {};
                for ii=1:length(uniquex)
                    pos = find(x==uniquex(ii));
                    eval(['z = meannan(ResponseForPlot.' cenames{jj} '(pos,:));']);
                    eval(['stdz = stdnan(ResponseForPlot.' cenames{jj} '(pos,:))/sqrt(length(pos));']);

                    if v(1)=='6'
                        ff = errorbar(time{1}, z, stdz);
                        set(ff,'color',colors(ii, :));
                        thislabel = [thislabel {''} {(num2str(uniquex(ii)))}];
                    elseif v(1)=='7'
                        errorbar(time{1}, z, stdz, 'color', colors(ii, :));
                        thislabel = [thislabel {(num2str(uniquex(ii)))}];
                    end
                    hold on;
                end
                legend(thislabel);
                xlabel('time in second');
                ylabel('BOLD');
                title([cenames{jj} ' ' correlatornames{jj}{kk}])
            end
            
            count = count + 1;
        end
    end

    peakpoint = inputdlg({'I will plot peak BOLD. What points in time do you considered as peak BOLD?'}, 'Peak point',1,{'4 6'});
    if isempty(peakpoint); return; end
    peakpoint = str2num(peakpoint{1});
    
    
	%list all plot names in command window for selection, effect(s) to be
	%removed
    clear promt;
    promt{1} = 'Which effect(s) do you want me to remove first? (If none to remove, leave it blank)';
    promt{2} = 'index   event    correlator';
    count = 1;
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            promt{count+2} = sprintf('%d         %s         %s', count,  cenames{jj}, correlatornames{jj}{kk});
            count = count + 1;
        end
    end
    %promt{count+2} = 'Which effect(s) do you want me to remove first?';

    promt = char(promt);
    name='Please select which effect to remove';
    numlines=1;
    defaultanswer={''};

	whichremovefirst=inputdlg(promt,name,numlines,defaultanswer);
    if isempty(whichremovefirst)
        whichremovefirst = [];
    elseif isempty(deblank(whichremovefirst{1}))
        whichremovefirst = [];
    else
        whichremovefirst = eval(whichremovefirst{1});
    end
    
   
%% remove effect first
    tmpResponseForPlot = ResponseForPlot;
    for jj=1:length(cenames)
        eval(['tmpResponseForPlot.' cenames{jj} ' = mean(tmpResponseForPlot.' cenames{jj} '(:,findind(time{1},peakpoint)), 2);']);
    end
    
    if ~isempty(whichremovefirst)
        count = 1;
        for jj=1:length(cenames)
            eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
            eval(['y = tmpResponseForPlot.' cenames{jj} ';']);
            x = [];
            if size(y,1)==1
                y = y';
            end
            for kk=1:length(correlatornames{jj})
                if ~ismember(count, whichremovefirst)
                    count = count + 1;
                    continue;
                end

                eval(['tmpx = C.' cenames{jj} '.' correlatornames{jj}{kk} ';']);

                % convert to colum vector
                if size(tmpx,1)==1
                    tmpx = tmpx';
                end
                x = [x tmpx];



                count = count + 1;
            end

            % remove NaN
            tmpposx = find(isnan(mean(x,2)) | isinf(mean(x,2)));
            tmpposy = find(isnan(mean(y,2)) | isinf(mean(y,2)));
            x([tmpposx; tmpposy],:) = [];
            y([tmpposx; tmpposy],:) = [];

            b = linearregression(y,x);
            y = y - x*b(1:end-1) - b(end);
            eval(['tmpResponseForPlot.' cenames{jj} ' = y;']);
        end
    end

%% now plot the residual amplitude against correlator
    count = 1;
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            if ~ismember(count, whichplot)
                count = count + 1;
                continue;
            end

            eval(['y = tmpResponseForPlot.' cenames{jj} ';']);
            eval(['x = C.' cenames{jj} '.' correlatornames{jj}{kk} ';']);
            
           
            figure
            plot(x,y, 'sb');
            hold on;            
            plot2(x,y,'r',1);
            
            xlabel(correlatornames{jj}{kk})
            ylabel('peak relative BOLD');
            legend(cenames{jj})
            try
                % convert to colum vector
                if size(x,1)==1
                    x = x';
                end
                if size(y,1)==1
                    y = y';
                end
                
                % remove NaN
                tmpposx = find(isnan(x) | isinf(x));
                tmpposy = find(isnan(y) | isinf(y));
                x([tmpposx; tmpposy]) = [];
                y([tmpposx; tmpposy]) = [];

                b = linearregression(y,x);
                
                totalvar = var(y);
                residvar = var(y - x*b(1) - b(2));
                modelvar = totalvar - residvar;
                dfm = 1;
                dft = length(x) - 1;
                dfr = dft - dfm;
                mmodelvar = modelvar/dfm;
                mresidvar = residvar/dfr;
                F = mmodelvar/mresidvar;
                pvalue = 1-spm_Fcdf(F, [dfm dfr]);
                
                title(sprintf('pValue=%s; slope=%g; constant=%g', pvalue, b(1), b(2)))
                xx=[min(x):(max(x)-min(x))/10:max(x)];
                yy=b(1) * xx + b(2);
                hold on;
                plot(xx,yy,'g');
            end
            count = count + 1;
        end
    end
end

function index = findind(vector1, vector2)
kk = 1;
for ii=1:length(vector2)
    tmp = find(vector1 == vector2(ii));
    if ~isempty(tmp)
        index(kk) = tmp;
        kk = kk+1;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% plot ROI, individual by individual, with behavior timing on it
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_plotIndividualROIWithBehavior(hObject, eventdata)
global TR_;
handles = guidata(hObject);

if findstr('SPM2',spm('ver'))
    P = spm_get([0:1],'*.mat','Select ROI result file');
elseif findstr('SPM5',spm('ver'))
    P = spm_select([0:1],'mat','Select ROI result file');
end

if isempty(P)
    return;
end
load(P);

answer = getVariable({'headmovement data (optional)'});


try
    Ph = evalin('base',answer{1}); % headmovement
catch
    Ph = [];
end

colors=[1 0 0;
    0 0 1;
    0 0 0;
    1 0 1;
    0 1 1;
    0 1 0;
    1 1 0;
    1 0.5 0.5;
    .5 .5 1;
    .5 .5 .5;
    1 .5 0;
    1 0 .5;
    0 1 .5;
    .5 1 0;
    0 .5 1;
    .5 0 1];
colors = repmat(colors, 10, 1);


figure;
Maximize(gcf);
for ii=1:length(event)
    kk = mod(ii,8);
    if kk == 0
        kk = 8;
    end
    subplot(4,2,kk)

    m = mean(wholeOriginal{ii});
    s = std(wholeOriginal{ii});
    mx = max(wholeOriginal{ii});
    mn = min(wholeOriginal{ii});
    ylowerlimit = mn;
    
    hm = event{ii};
    hmcell = struct2cell(hm);
    hmname = fieldnames(hm);
    for jj=1:length(hmname)
        toplot1 = [];
        toplot2 = [];
        for mm = 1:length(hmcell{jj})
            toplot1 = [toplot1 hmcell{jj}(mm) hmcell{jj}(mm) nan];
            toplot2 = [toplot2 m-s m+s nan];
        end
        plot(toplot1, toplot2, 'color', colors(jj,:));
        hold on
    end
    legend(hmname)
    title(['subject ' num2str(ii)])
    xlabel('time in s');

    if ~isempty(Ph) % if there is head movement data
        tmp = load(deblank(Ph(ii,:)));
        tmp2 = fieldnames(tmp);
        hm = eval(['tmp.' tmp2{1}]);
        hmcell = struct2cell(hm);
        hmname = fieldnames(hm);
        for jj=1:length(hmname)
            plot(TR_*[0:length(hmcell{jj})-1],mn-s+s*hmcell{jj}/max(abs(hmcell{jj})),'color', [.5 .5 .5]);
            ylowerlimit = min(ylowerlimit, min(mn-s+s*hmcell{jj}/max(abs(hmcell{jj}))));
            hold on
        end        
    end

    % plot imaging data
    plot(TR_*[0:length(wholeOriginal{ii})-1], wholeOriginal{ii}, 'k');
    hold on;
    ylim([ylowerlimit mx])

    if mod(ii,8) == 0 && ii~=length(event)
        figure;
        Maximize(gcf);
    end     
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% plot ROI, individual by individual
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_plotIndividualROI(hObject, eventdata, filename)
handles = guidata(hObject);

if ~exist('filename')
    if findstr('SPM2',spm('ver'))
        P = spm_get([0:1],'*.mat','Select ROI result file');
    elseif findstr('SPM5',spm('ver'))
        P = spm_select([0:1],'mat','Select ROI result file');
    end
else
    P = filename;
end

if isempty(P)
    return;
end
load(deblank(P));

if ~exist('mni')
    set(handles.infoTextBox, 'string', {'Oops! I donnt find the variable mni. Please open a proper mat file.'}); 
    return
end

delete(gcf)
xjview(mni);

response = eventResponse{1};

eventisstruct = 0;
if isstruct(response)
    eventisstruct = 1;
	f = fieldnames(response);
    response = struct2cell(response); 
else
    f = num2cell([1:length(response)]);
    for ii=1:length(response)
        f{ii} = num2str(f{ii});
    end
end
for jj=1:length(response)
    response{jj} = [];
end

if exist('correlator') & ~isempty(correlator) & eventisstruct==1
    C = correlator{1};
    cenames = fields(C);
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            eval(['C.' cenames{jj} '.' correlatornames{jj}{kk} '= [];']);
        end
    end
    for ii=1:length(correlator)
        for jj=1:length(cenames)
            eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
            for kk=1:length(correlatornames{jj})
                tmp2 = eval(['correlator{ii}.' cenames{jj} '.' correlatornames{jj}{kk} ';']);
                if size(tmp2,1) == 1
                    tmp2 = tmp2';
                end
                eval(['C.' cenames{jj} '.' correlatornames{jj}{kk} '= [C.' cenames{jj} '.' correlatornames{jj}{kk} '; tmp2];']);
            end
        end
    end
end

if eventisstruct == 1
    for ii=1:length(eventResponse)
        for kk=1:length(f)
            eval(['tmptmp{kk} = eventResponse{ii}.' f{kk} ';']);
        end
        eventResponse{ii} = tmptmp;
    end
end

for ii=1:length(eventResponse)
    for jj=1:length(eventResponse{1})
        response{jj} = [response{jj}; eventResponse{ii}{jj}];
    end
end

if isempty(response)
    msgbox('Nothing to plot');
    return
end




v = version;
figure;
%colors='rbkmcgyrbkmcgyrbkmcgyrbkmcgyrbkmcgyrbkmcgyrbkmcgy';
colors=[1 0 0;
    0 0 1;
    0 0 0;
    1 0 1;
    0 1 1;
    0 1 0;
    1 1 0;
    1 0.5 0.5;
    .5 .5 1;
    .5 .5 .5;
    1 .5 0;
    1 0 .5;
    0 1 .5;
    .5 1 0;
    0 .5 1;
    .5 0 1];
colors = repmat(colors, 10, 1);
thislabel = {};
for ii=1:length(response)
    if v(1)=='6'
        ff = errorbar(time{1}, meannan(response{ii}), stdnan(response{ii})/sqrt(size(response{ii},1)));
        set(ff,'color',colors(ii, :));
        thislabel = [thislabel {''} f(ii)];
    elseif v(1)=='7'
        errorbar(time{1}, meannan(response{ii}), stdnan(response{ii})/sqrt(size(response{ii},1)), 'color', colors(ii, :));
        thislabel = [thislabel f(ii)];
    end
    hold on;
end
xlabel('time (s)');
ylabel('relative BOLD');
legend(thislabel)
hold off;

figure;
Maximize(gcf);
for jj=1:length(eventResponse)
    wheretoplot = mod(jj, 16);
    if wheretoplot == 0; 
        wheretoplot = 16; 
        subplot(4,4,wheretoplot)
    else
        subplot(4,4,wheretoplot)
    end
    for ii=1:length(response)
        if v(1)=='6'
            ff = errorbar(time{1}, meannan(eventResponse{jj}{ii}), stdnan(eventResponse{jj}{ii})/sqrt(size(eventResponse{jj}{ii},1)));
            set(ff,'color',colors(ii, :));
            %thislabel = [thislabel {''} f(ii)];
        elseif v(1)=='7'
            errorbar(time{1}, meannan(eventResponse{jj}{ii}), stdnan(eventResponse{jj}{ii})/sqrt(size(eventResponse{jj}{ii},1)), 'color', colors(ii, :));
            %errorbar(time{1}, meannan(response{ii}), stdnan(response{ii})/sqrt(size(response{ii},1)), 'color', colors(ii, :));
            %thislabel = [thislabel f(ii)];
        end
        hold on;
    end
    title(['subject ' num2str(jj)])
    
    if mod(jj, 16)==0 && jj~=length(eventResponse)
        figure;
        Maximize(gcf);
    end
end



[row,col]  = size(response);
if row == 1; response = response'; end;
for jj=1:length(eventResponse)
    ResponseForPlot{jj} = cell2struct(eventResponse{jj}, f, 2);
end

if exist('correlator') & ~isempty(correlator)
    %list all plot names in command window for selection
    promt{1} = 'I will plot BOLD amplitude vs correlator. Here are all the correlators:';
    promt{2} = 'index   event    correlator';
    count = 1;
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            promt{count+2} = sprintf('%d         %s         %s', count,  cenames{jj}, correlatornames{jj}{kk});
            count = count + 1;
        end
    end
    promt{count+2} = 'Which plots do you want me to show?';

    promt = char(promt);
    name='Please select which plots to show';
    numlines=1;
    defaultanswer={['[1:'  num2str(count-1) ']']};

    whichplot=inputdlg(promt,name,numlines,defaultanswer);
    if isempty(whichplot)
        return;
    end    
    whichplot = eval(whichplot{1});

    peakpoint = inputdlg({'What points in time do you considered as peak BOLD?'}, 'Peak point',1,{'4 6'});
    if isempty(peakpoint); return; end
    peakpoint = str2num(peakpoint{1});   
    
    count = 1;
    for jj=1:length(cenames)
        eval(['correlatornames{jj} = fields(C.' cenames{jj} ');']);
        for kk=1:length(correlatornames{jj})
            if ~ismember(count, whichplot)
                count = count + 1;
                continue;
            end
            figure;
            Maximize(gcf);
            for ii=1:length(eventResponse)
                wheretoplot = mod(ii, 16);
                if wheretoplot == 0;
                    wheretoplot = 16;
                    subplot(4,4,wheretoplot)
                else
                    subplot(4,4,wheretoplot)
                end


                eval(['y = mean(ResponseForPlot{ii}.' cenames{jj} '(:,findind(time{1},peakpoint)), 2);']);
                eval(['x = correlator{ii}.' cenames{jj} '.' correlatornames{jj}{kk} ';']);

                %plot(x,y, 'sb');
                %hold on;            
                pos = find(isnan(x));
                x(pos) = [];
                y(pos) = [];
                pos = find(isnan(y));
                x(pos) = [];
                y(pos) = [];                
                if length(x) < 1
                    continue;
                end
                plot2(x,y,'r',1);

                xlabel(correlatornames{jj}{kk})
                ylabel(['peak BOLD at ' cenames{jj}]);
                %legend(cenames{jj})
                try
                    % convert to colum vector
                    if size(x,1)==1
                        x = x';
                    end
                    if size(y,1)==1
                        y = y';
                    end

                    % remove NaN
                    tmpposx = find(isnan(x) | isinf(x));
                    tmpposy = find(isnan(y) | isinf(y));
                    x([tmpposx; tmpposy]) = [];
                    y([tmpposx; tmpposy]) = [];

                    b = linearregression(y,x);

                    totalvar = var(y);
                    residvar = var(y - x*b(1) - b(2));
                    modelvar = totalvar - residvar;
                    dfm = 1;
                    dft = length(x) - 1;
                    dfr = dft - dfm;
                    mmodelvar = modelvar/dfm;
                    mresidvar = residvar/dfr;
                    F = mmodelvar/mresidvar;
                    pvalue = 1-spm_Fcdf(F, [dfm dfr]);

                    title(sprintf('subject %d, pValue=%s', ii, pvalue))
                    xx=[min(x):(max(x)-min(x))/10:max(x)];
                    yy=b(1) * xx + b(2);
                    hold on;
                    plot(xx,yy,'g');
                end
                if mod(ii, 16)==0 && ii~=length(eventResponse)
                    figure;
                    Maximize(gcf);
                end
            end % end of loop subject
            count = count + 1;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% plot ROI: correlation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_plotCorrelationROI(hObject, eventdata)

if findstr('SPM2',spm('ver'))
    P = spm_get([0:100],'*.mat','Select ROI result file(s)');
elseif findstr('SPM5',spm('ver'))
    P = spm_select([0:100],'mat','Select ROI result file(s)');
end

if isempty(P)
    return;
end

N = size(P,1);
t = [-20:20];

v = version;

if N == 1
    load(deblank(P));
    tmp = [];
    for jj=1:length(wholeAbsolute)
        tmp = [tmp; corrlag(wholeAbsolute{jj},wholeAbsolute{jj},t)];
    end
    figure;
    errorbar(t, mean(tmp,1), std(tmp)/sqrt(length(wholeAbsolute)));
    xlabel('lag in scan number');
    ylabel('auto correlation');
    title(deblank(P))
    return
else
    figure;
    colors='rbkmcgyrbkmcgyrbkmcgyrbkmcgyrbkmcgyrbkmcgyrbkmcgy';
    legendlabel = {};
    totalPermute = nchoosek([1:N],2);
    for ii=1:size(totalPermute,1)
        tmp = [];
        s1 = load(deblank(P(totalPermute(ii,1),:)));
        s2 = load(deblank(P(totalPermute(ii,2),:)));        
        for jj=1:length(s1.wholeAbsolute)
            tmp = [tmp; corrlag(s1.wholeAbsolute{jj},s2.wholeAbsolute{jj},t)];
        end
        errorbar(t, mean(tmp,1), std(tmp)/sqrt(length(s1.wholeAbsolute)), colors(ii));
        xlabel('lag in scan number');
        ylabel('cross correlation');
        hold on;
        if v(1) == '6'
            legendlabel = [legendlabel, {'', P(totalPermute(ii,:),:)}];
        elseif v(1) == '7'
            legendlabel = [legendlabel, {P(totalPermute(ii,:),:)}];
        end
    end
    legend(legendlabel);
    hold off;
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% whole brain correlation analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_wholeBrainCorrelation(hObject, eventdata)
% answer = getVariable({'Which data file?'});
% if isempty(answer)
%     return;
% end
% if isempty(answer{1})
%     return;
% end
% P = evalin('base',answer{1}); % subject directories
if findstr('SPM2',spm('ver'))
    P = spm_get([0:1],'*.mat','Whole Brain Signal file');
elseif findstr('SPM5',spm('ver'))
    P = spm_select([0:1],'mat','Whole Brain Signal file');
end

if isempty(P)
    return
end

load(deblank(P));
M = double(M);

N=zeros(size(M,4), size(M,1) * size(M,2) * size(M,3));

mm = 1;
for ii=1:2:size(M,1)
    for jj=1:3:size(M,2)
        for kk=1:2:size(M,3)
            N(:, mm) = squeeze(M(ii,jj,kk,:));
            mm = mm + 1;
            cor(mm,:) = [ii jj kk];
        end
    end
end

N(:, mm:end) = [];

C=corrcoef(N);
C(find(isnan(C))) = 0;

for ii=1:4
    CC{ii} = abs(C)>(0.5+ii/10);
end

K = {zeros(1, size(C,1)),zeros(1, size(C,1)),zeros(1, size(C,1)),zeros(1, size(C,1))};
for ii=1:size(C,1)
    for jj=1:4
        K{jj}(ii) = sum(CC{jj}(ii,:)) - 1;
    end
end

x=[0:size(C,1)];
for jj=1:4
    y{jj}=zeros(size(x));
end

for ii=1:length(x)
    for jj=1:4
        y{jj}(ii) = sum(K{jj} == x(ii));
    end
end

figure;
loglog(x,y{1}, x, y{2}, x, y{3}, x, y{4});
xlabel('degree');
ylabel('counts');
legend('0.6', '0.7', '0.8', '0.9');
axis equal


z = sum(abs(C))-1;
zx = 0:max(z);
for ii=1:length(zx)-1
    pos = find(z>zx(ii) & z<=zx(ii+1));
    zy(ii) = length(pos);
end
figure;
loglog(zx(1:length(zy)),zy);
xlabel('weighted degree');
ylabel('counts');
axis equal

xjview(cor2mni(coord, handles.M{1}), z);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% get/select variable names in base workspace
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function vars = getVariable(titles)
vars = evalin('base','who');
height = 0.07;
f = dialog('unit','normalized', 'menubar','none', 'position', [0.3 0.2 0.4 0.4], 'name', 'pick variables', 'NumberTitle','off');
uicontrol(f,'style','text',...
        'unit','normalized',...
        'String', 'Available variables', ...
        'position',[0 0.9 0.5 0.1/2]);    
variableListbox = uicontrol(f,'style','listbox','tag','variableListbox',...
        'unit','normalized',...
        'String', vars, ...
        'value',1,...
        'position',[0 0 0.5 0.9]);
okPush = uicontrol(f,'style','push',...
        'unit','normalized',...
        'String', 'OK', ...
        'position',[0.55 0.05 0.2 height],...
        'callback','uiresume');
cancelPush = uicontrol(f,'style','push',...
        'unit','normalized',...
        'String', 'Cancel', ...
        'position',[0.75 0.05 0.2 height],...
        'callback','delete(gcf)');  

uicontrol(f,'style','text',...
        'unit','normalized',...
        'String', titles{1}, ...
        'position',[0.5 0.9 0.5 0.1/2]);   
imageDataPush = uicontrol(f,'style','push',...
        'unit','normalized',...
        'String', '->', ...
        'position',[0.5 0.8 0.1 height],...
        'callback', [...
            'xyouneverfindme=findobj(get(gcf,''Children''),''flat'',''tag'',''imageDataEdit'');'...
            'yyouneverfindme=findobj(get(gcf,''Children''),''flat'',''tag'',''variableListbox'');'...
            'allyouneverfindme=get(yyouneverfindme,''string'');'...
            'selectedyouneverfindme=allyouneverfindme{get(yyouneverfindme,''value'')};'...
            'set(xyouneverfindme,''string'',selectedyouneverfindme);'...
            'clear allyouneverfindme selectedyouneverfindme xyouneverfindme yyouneverfindme;']);      
imageDataEdit = uicontrol(f,'style','edit','tag','imageDataEdit',...
        'unit','normalized',...
        'String', '', ...
        'BackgroundColor', 'w',...
        'position',[0.6 0.8 0.3 height]);

if length(titles)>=2
uicontrol(f,'style','text',...
        'unit','normalized',...
        'String', titles{2}, ...
        'position',[0.5 0.7 0.5 0.1/2]);   
eventDataPush = uicontrol(f,'style','push',...
        'unit','normalized',...
        'String', '->', ...
        'position',[0.5 0.6 0.1 height],...
        'callback', [...
            'xyouneverfindme=findobj(get(gcf,''Children''),''flat'',''tag'',''eventDataEdit'');'...
            'yyouneverfindme=findobj(get(gcf,''Children''),''flat'',''tag'',''variableListbox'');'...
            'allyouneverfindme=get(yyouneverfindme,''string'');'...
            'selectedyouneverfindme=allyouneverfindme{get(yyouneverfindme,''value'')};'...
            'set(xyouneverfindme,''string'',selectedyouneverfindme);'...
            'clear allyouneverfindme selectedyouneverfindme xyouneverfindme yyouneverfindme;']);           

eventDataEdit = uicontrol(f,'style','edit','tag','eventDataEdit',...
        'unit','normalized',...
        'String', '', ...
        'BackgroundColor', 'w',...
        'position',[0.6 0.6 0.3 height]);    
end
if length(titles)>=3
uicontrol(f,'style','text',...
        'unit','normalized',...
        'String', titles{3}, ...
        'position',[0.5 0.5 0.5 0.1/2]);   
correlatorDataPush = uicontrol(f,'style','push',...
        'unit','normalized',...
        'String', '->', ...
        'position',[0.5 0.4 0.1 height],...
        'callback', [...
            'xyouneverfindme=findobj(get(gcf,''Children''),''flat'',''tag'',''correlatorDataEdit'');'...
            'yyouneverfindme=findobj(get(gcf,''Children''),''flat'',''tag'',''variableListbox'');'...
            'allyouneverfindme=get(yyouneverfindme,''string'');'...
            'selectedyouneverfindme=allyouneverfindme{get(yyouneverfindme,''value'')};'...
            'set(xyouneverfindme,''string'',selectedyouneverfindme);'...
            'clear allyouneverfindme selectedyouneverfindme xyouneverfindme yyouneverfindme;']);           
    
correlatorDataEdit = uicontrol(f,'style','edit','tag','correlatorDataEdit',...
        'unit','normalized',...
        'String', '', ...
        'BackgroundColor', 'w',...
        'position',[0.6 0.4 0.3 height]); 
end
if length(titles)>=4
uicontrol(f,'style','text',...
        'unit','normalized',...
        'String', titles{4}, ...
        'position',[0.5 0.3 0.5 0.1/2]);   
otherDataPush = uicontrol(f,'style','push',...
        'unit','normalized',...
        'String', '->', ...
        'position',[0.5 0.2 0.1 height],...
        'callback', [...
            'xyouneverfindme=findobj(get(gcf,''Children''),''flat'',''tag'',''otherDataEdit'');'...
            'yyouneverfindme=findobj(get(gcf,''Children''),''flat'',''tag'',''variableListbox'');'...
            'allyouneverfindme=get(yyouneverfindme,''string'');'...
            'selectedyouneverfindme=allyouneverfindme{get(yyouneverfindme,''value'')};'...
            'set(xyouneverfindme,''string'',selectedyouneverfindme);'...
            'clear allyouneverfindme selectedyouneverfindme xyouneverfindme yyouneverfindme;']);           
   
otherDataEdit = uicontrol(f,'style','edit','tag','otherDataEdit',...
        'unit','normalized',...
        'String', '', ...
        'BackgroundColor', 'w',...
        'position',[0.6 0.2 0.3 height]); 
end

uiwait(f);
try
    var{1} = get(imageDataEdit,'string');
    if length(titles)>=2
        var{2} = get(eventDataEdit,'string');
    end
    if length(titles)>=3
        var{3} = get(correlatorDataEdit,'string');    
    end
    if length(titles)>=4
        var{4} = get(otherDataEdit,'string');    
    end
    vars = var;
    delete(f);
catch
    vars = {};
    try    
        delete(f);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% quit xjview
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_quit(hObject, eventdata, warnstate)
warning(warnstate)
% try
%     rmdir('xjviewtmp');
% end
delete(gcf);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% mouse double click
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function figureMouseUpFcn(hObject, eventdata)
status = get(hObject, 'SelectionType');

% double click
if strcmp(status, 'open')
    handles = guidata(hObject);
    CallBack_loadImagePush(handles.loadImagePush, eventdata);
else
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% change edit image file name
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_imageFileEdit(hObject, eventdata)
handles = guidata(hObject);
filename = get(hObject, 'String');
filename = str2cell(filename);
CallBack_loadImagePush(handles.loadImagePush, eventdata, filename);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% click load image file button
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_loadImagePush(hObject, eventdata, thisfilename)
handles = guidata(hObject);

handles.imageFileName=[]; handles.M=[]; handles.DIM=[]; handles.TF=[]; handles.df=[]; 
handles.mni=[]; handles.intensity=[]; handles.currentmni=[]; handles.currentintensity=[]; handles.currentDisplayMNI=[]; handles.currentDisplayIntensity=[];

if ~exist('thisfilename')
    thisfilename = ''; 
end

if isstruct(thisfilename)
    handles.imageFileName = {''};
    handles.mni = {thisfilename.mni};
	handles.intensity = {thisfilename.intensity};
    handles.M = {thisfilename.M};
    handles.DIM = {thisfilename.DIM};        
    handles.TF = {'S'};
    handles.df = {1e6};
    handles.pValue = 1;
    set(handles.pValueEdit, 'string', '1');
    handles.clusterSizeThreshold = 0;
    set(handles.clusterSizeThresholdEdit, 'String', '0');
else
    [handles.imageFileName, handles.M, handles.DIM, handles.TF, handles.df, handles.mni, handles.intensity] = getImageFile(thisfilename);
end

if isempty(handles.imageFileName)
    return
end

different = 0; % image files are same or different?
if length(handles.TF)>1
    for ii=1:length(handles.TF)
        if strcmp(handles.TF{ii},handles.TF{1}) & isequal(handles.df{ii},handles.df{1}) & isequal(handles.M{ii},handles.M{1}) & isequal(handles.DIM{ii},handles.DIM{1})
            continue;
        else
            warndlg('Images are from different statistics or sources.', 'Warning');
            %set(handles.infoTextBox, 'string', 'Images are from different statistics or sources.'); 
            beep;
            different = 1;
            break;
        end
    end
end

% reset files with empty df/TF to df=1 and TF='S'. 'S'=='T' but has a tag
% meaning it is changed.
resetTF = 0;
for ii=1:length(handles.TF)
    if isempty(handles.TF{ii}) 
        handles.TF{ii} = 'S';
        resetTF = 1;
    end
    if isempty(handles.df{ii}) | isequal(handles.df{ii},0)
        handles.df{ii} = 1e6;
        resetTF = 1;
    end
end
    
handles.currentmni = handles.mni;
handles.currentintensity = handles.intensity;

set(handles.dfEdit, 'String', cell2str(handles.df));
set(handles.imageFileEdit, 'String', cell2str(handles.imageFileName)); % s=-log10(p)
maxs = maxcell(t2s(cellmax(handles.intensity,'abs'),handles.df, handles.TF),'abs'); 
if isinf(maxs); maxs = 20; end
set(handles.slider, 'Max', maxs, 'Min', 0, 'sliderstep',[min(maxs/100,0.05),min(maxs/100,0.05)]);
if handles.TF{1}=='T' & different == 0
    str = [blanks(length(' intensit')) 'T=']; 
elseif handles.TF{1}=='F' & different == 0
    str = [blanks(length(' intensit')) 'F=']; 
else
    str=' intensity='; 
end
set(handles.intensityThresholdText, 'String', str);
set(handles.figure,'Name',['xjView: ' cell2str(handles.imageFileName)]);
try
    for ii=1:length(handles.hLegend)
        delete(handles.hLegend{ii});
    end
end
if length(handles.TF)>1
    colours = repmat([1 0 0;1 1 0;0 1 0;0 1 1;0 0 1;1 0 1],100,1);
    for ii=1:min(length(handles.TF),10)
        [tmp,filename] = fileparts(handles.imageFileName{ii});
        if ii == 10; filename = '......'; end
        pos0 = handles.sectionViewPosition;
        pos(1) = pos0(1)+pos0(3)/2+0.03;
        pos(2) = pos0(2)+ii/50-0.03;
        pos(3) = 0.12;
        pos(4) = 0.02;
        handles.hLegend{ii}=uicontrol(handles.figure, 'style','text',...
        'unit','normalized','position',pos,...
        'string', filename,...
        'horizontal','left',...
        'fontweight','bold',...
        'ForeGroundColor',colours(ii,:));
    end
	c_names = {'red';'yellow';'green';'cyan';'blue';'magenta'};
end
guidata(hObject, handles);

global M_;
global DIM_;
M_ = handles.M{1};
DIM_ = handles.DIM{1};

if resetTF==1
    set(handles.pValueEdit,'string',1);
end
CallBack_pValueEdit(handles.pValueEdit, eventdata);

% display info
try
    urlread('http://people.hnl.bcm.edu/cuixu/xjView/guestbook/stat.php');
	s = urlread('http://people.hnl.bcm.tmc.edu/cuixu/xjView/toUser.txt');
    report{1} = s;
catch
    report{1} = 'Welcome to xjView 4';
end
for jj=1:length(handles.imageFileName)
    report{2+6*(jj-1)} = cell2str(handles.imageFileName(jj));
    if handles.TF{jj} == 'T' | handles.TF{jj} == 'F'
        report{3+6*(jj-1)} = ['This is a ' handles.TF{jj} ' test image.'];
    else
        report{3+6*(jj-1)} = '';%['I don''t know what test this image came from.'];
    end
    report{4+6*(jj-1)} = 'mat = ';
    report{5+6*(jj-1)} = num2str(handles.M{jj});
    report{6+6*(jj-1)} = 'dimension = ';
    report{7+6*(jj-1)} = num2str(handles.DIM{jj});    
end
set(handles.infoTextBox, 'string', report);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% pValueEdit
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_pValueEdit(hObject, eventdata)
handles = guidata(hObject);
tmp = str2double(get(hObject, 'String'));
tmps = -log10(tmp);
if isnan(tmp) | tmp < 0 | tmp > 1
    errordlg('I don''t understand the input.','error');
    set(hObject, 'String', handles.pValue);
    return
end

if tmps>get(handles.slider,'max') | tmps<get(handles.slider,'min')
    errordlg('pValue is too small. No suprathreshold voxels.','error');
    set(hObject, 'String', handles.pValue);
    return
end

if isempty(handles.df) | handles.df{1}==0
    set(handles.pValueEdit,'String', 'NaN');
end

handles.pValue = tmp;
handles.intensityThreshold = p2t(num2cell(handles.pValue*ones(1,length(handles.TF))), handles.df, handles.TF);

set(handles.pValueEdit,'string', num2str(handles.pValue));
guidata(hObject, handles);
CallBack_slider(hObject, eventdata, -log10(handles.pValue));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% intensity threshold edit
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_intensityThresholdEdit(hObject, eventdata)
handles = guidata(hObject);
tmp = str2double(get(hObject, 'String'));
if isnan(tmp) | tmp<0
    errordlg('Please input a single valid number bigger than 0.','error');
    return
end

if tmp > maxcell(cellmax(handles.intensity,'abs'))
    tmp = maxcell(cellmax(handles.intensity,'abs'));
end

handles.pValue = t2p(tmp, handles.df{1}, handles.TF{1});
handles.intensityThreshold = p2t(num2cell(handles.pValue*ones(1,length(handles.TF))), handles.df, handles.TF);
set(handles.slider,'Value', -log10(handles.pValue));
set(handles.intensityThresholdEdit, 'String', cell2str(handles.intensityThreshold));
guidata(hObject, handles);
CallBack_slider(hObject, eventdata, -log10(handles.pValue));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% slider bar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_slider(hObject, eventdata, value)

handles = guidata(hObject);

if exist('value')
    s = value;
    if s > get(handles.slider,'max')
        s = get(handles.slider,'max')*0.99;
    end
    if s < get(handles.slider,'min')
        s = get(handles.slider,'min');
    end
else
    s = get(hObject,'Value');
end

set(handles.slider, 'value', s);
pvalue = 10^(-s);
set(handles.pValueEdit,'String',num2str(pvalue));
t = p2t(num2cell(pvalue*ones(1,length(handles.TF))), handles.df, handles.TF);
handles.intensityThreshold = t;
set(handles.intensityThresholdEdit, 'String', cell2str(handles.intensityThreshold));
handles.pValue = pvalue;
for ii=1:length(handles.TF)
    pos{ii} = find(abs(handles.intensity{ii})>=t{ii});
    handles.currentintensity{ii} = handles.intensity{ii}(pos{ii});
    handles.currentmni{ii} = handles.mni{ii}(pos{ii},:);
end

guidata(hObject,handles);

if get(handles.allIntensityRadio, 'Value')
    CallBack_allIntensityRadio(handles.allIntensityRadio, eventdata);
elseif get(handles.positiveIntensityRadio, 'Value')
    CallBack_allIntensityRadio(handles.positiveIntensityRadio, eventdata, '+');
elseif get(handles.negativeIntensityRadio, 'Value')
    CallBack_allIntensityRadio(handles.negativeIntensityRadio, eventdata, '-');
end

set(handles.infoTextBox, 'string', {'Don''t drag the slider bar too fast. Release your mouse button at least 1 second later.', ...
    'This sounds stupid. But there is a bug (probably MatLab bug) which I can'' fix right now.', ...
    'I suggest you confirm the correctness of the current display by press Enter in the pValue edit box.'}); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% display intensity all+- radios
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_allIntensityRadio(hObject, eventdata, pnall)
% pnall = '+', '-', or 'c'. c means current (simply update drawing)
%
handles = guidata(hObject);

currentselect = [];
if get(handles.allIntensityRadio, 'Value'); currentselect = handles.allIntensityRadio; thispnall = 'a'; end
if get(handles.positiveIntensityRadio, 'Value'); currentselect = handles.positiveIntensityRadio; thispnall = '+';  end
if get(handles.negativeIntensityRadio, 'Value'); currentselect = handles.negativeIntensityRadio; thispnall = '-';  end

set(handles.allIntensityRadio, 'Value', 0);
set(handles.positiveIntensityRadio, 'Value', 0);
set(handles.negativeIntensityRadio, 'Value', 0);

if exist('pnall')
    if pnall=='c'
        hObject = currentselect;
        pnall = thispnall;
    end
end

set(hObject, 'Value', 1);

if ~isfield(handles,'currentintensity')
    return
end
for ii=1:length(handles.TF)
	if exist('pnall')
        if pnall == '-'
            pos{ii} = find(handles.currentintensity{ii} < 0);
        elseif pnall == '+'
            pos{ii} = find(handles.currentintensity{ii} > 0);
        elseif pnall == 'a'
            pos{ii} = 1:length(handles.currentintensity{ii});
        end
	else
        pos{ii} = 1:length(handles.currentintensity{ii});
	end
    intensity{ii} = handles.currentintensity{ii}(pos{ii});
    mni{ii} = handles.currentmni{ii}(pos{ii},:);
    cor{ii} = mni2cor(mni{ii}, handles.M{ii});

	if ~isempty(cor{ii})
		A = spm_clusters(cor{ii}');
		pos0 = [];
		for kk = 1:max(A)
            jj = find(A == kk);
            if length(jj) >= handles.clusterSizeThreshold; pos0 = [pos0 jj]; end
		end
		handles.currentDisplayMNI{ii} = mni{ii}(pos0,:);
		handles.currentDisplayIntensity{ii} = intensity{ii}(pos0);
	else
		handles.currentDisplayMNI{ii} = mni{ii}([],:);
		handles.currentDisplayIntensity{ii} = intensity{ii}([]);    
	end
end
[handles.hReg, handles.hSection, handles.hcolorbar] = Draw(handles.currentDisplayMNI, handles.currentDisplayIntensity, hObject, handles);
if get(handles.allIntensityRadio, 'Value') & max(handles.currentDisplayIntensity{1}) < 0
    warndlg('No supra-threshold positive intensity. Only negative intensity is displayed.');
end

try
	set(handles.figure,'currentaxes', handles.glassViewAxes);
	xrange = xlim;
	yrange = ylim;
	try
        delete(handles.hGlassText)
	end
    if ~isempty(handles.selectedCluster,1)
    	%handles.hGlassText = text(xrange(1)+diff(xrange)*0.6, yrange(1)+diff(yrange)*0.9, [num2str(size(handles.selectedCluster,1)) ' clusters selected']);
        set(handles.infoTextBox, 'string', [num2str(size(handles.selectedCluster,1)) ' clusters selected']);
    end
end
guidata(hObject, handles);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% render view check?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_renderViewCheck(hObject, eventdata)
check = get(hObject, 'Value');
if check
    CallBack_allIntensityRadio(hObject, eventdata, 'c');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% which section view target file? list box
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_sectionViewListbox(hObject, eventdata)
handles = guidata(hObject);
contents = get(handles.sectionViewListbox,'String');
currentsel = contents{get(handles.sectionViewListbox,'Value')};
handles.sectionViewTargetFile = getSectionViewTargetFile(handles.spmdir, currentsel);
guidata(hObject, handles);
CallBack_allIntensityRadio(hObject, eventdata, 'c');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% get sectionviewtargetfile
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function sectionViewTargetFile = getSectionViewTargetFile(spmdir, selectedcontent)
if findstr('SPM2',spm('ver'))
    fileext = 'mnc';
elseif findstr('SPM5',spm('ver'))
    fileext = 'nii';
end
currentsel = selectedcontent;
if ~isempty(strfind(currentsel, 'single'))
    sectionViewTargetFile = fullfile(spmdir, 'canonical', ['single_subj_T1.' fileext]);
elseif ~isempty(strfind(currentsel, '152PD'))
    sectionViewTargetFile = fullfile(spmdir, 'canonical', ['avg152PD.' fileext]);
elseif ~isempty(strfind(currentsel, '152T1'))
    sectionViewTargetFile = fullfile(spmdir, 'canonical', ['avg152T1.' fileext]);    
elseif ~isempty(strfind(currentsel, '152T2'))
    sectionViewTargetFile = fullfile(spmdir, 'canonical', ['avg152T2.' fileext]);
elseif ~isempty(strfind(currentsel, '305T1'))
    sectionViewTargetFile = fullfile(spmdir, 'canonical', ['avg305T1.' fileext]);
elseif strcmp(currentsel, 'ch2')
    sectionViewTargetFile = fullfile(spmdir, 'canonical', 'ch2.img');    
elseif strcmp(currentsel, 'ch2bet')
    sectionViewTargetFile = fullfile(spmdir, 'canonical', 'ch2bet.img');      
elseif strcmp(currentsel, 'aal')
    sectionViewTargetFile = fullfile(spmdir, 'canonical', 'aal.img');      
elseif strcmp(currentsel, 'brodmann')
    sectionViewTargetFile = fullfile(spmdir, 'canonical', 'brodmann.img');      
end
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% xhairs in section view?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_xHairCheck(hObject, eventdata)
check = get(hObject, 'Value');
if check
    spm_orthviews('Xhairs','off');
else
    spm_orthviews('Xhairs','on');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% other target file for section view, push button
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_sectionViewMoreTargetPush(hObject, eventdata)    
handles = guidata(hObject);
[filename, pathname, filterindex] = uigetfile('*', 'Pick an target file');
if isequal(filename,0) | isequal(pathname,0)
   return;
end
handles.sectionViewTargetFile = fullfile(pathname, filename);
guidata(hObject, handles);
CallBack_allIntensityRadio(hObject, eventdata, 'c');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% set colorbar range
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_setTRangeEdit(hObject, eventdata)    
global TMAX_;
handles = guidata(hObject);
TMAX_ = get(hObject, 'String');
if isempty(str2num(TMAX_)) && ~strcmp(TMAX_, 'auto')
    return;
end
guidata(hObject, handles);
CallBack_allIntensityRadio(hObject, eventdata, 'c');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% change degree of freedome 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_dfEdit(hObject, eventdata) 
handles = guidata(hObject);
tmp = str2cell(get(hObject, 'String'));
if iscellstr(tmp)
    errordlg('Please input a valid number.','error');
    try
        set(hObject, 'String', handles.df);
    catch
        set(hObject, 'String', '');
    end
    return
end
handles.df=tmp;


if isfield(handles,'TF')
    t = p2t(mat2cell(handles.pValue*ones(1,length(handles.TF))), handles.df, handles.TF); 
    handles.intensityThreshold = t;
    set(handles.intensityThresholdEdit, 'String', cell2str(t));    
end

guidata(hObject, handles);
CallBack_pValueEdit(handles.pValueEdit, eventdata);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% get structure push
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_getStructurePush(hObject, eventdata) 
handles = guidata(hObject);
xyz = spm_XYZreg('GetCoords',handles.hReg);
tmp_coor = cuixuFindStructure(xyz', handles.DB);
set(handles.structureEdit,'String', tmp_coor{1});
handles.currentxyz = xyz';
guidata(hObject, handles);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% structure edit
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_structureEdit(action, xyz, thisfun, hReg)

try
    handles=guidata(hReg);
catch
    return;
end

handles.currentxyz = xyz';
guidata(hReg, handles);

try
    [tmp_coor, cellstructure] = cuixuFindStructure(xyz', handles.DB);
catch
    return;
end

set(handles.structureEdit,'String', tmp_coor{1});

for ii=[5 3 2 6 1 4]
    if strfind('undefined', cellstructure{ii})
        continue;
    else
        set(handles.searchContentEdit,'string', trimStructureStr(cellstructure{ii}));
        return;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% trim str
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = trimStructureStr(str)
pos = findstr('(', str);
if ~isempty(pos)
    str(pos-1:end)=[];
end
out = str;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% cluster size threshold edit
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_clusterSizeThresholdEdit(hObject, eventdata) 
handles = guidata(hObject);
tmp = str2double(get(hObject, 'String'));
if isnan(tmp) | tmp<0
    errordlg('Please input a valid number.','error');
    try
        set(hObject, 'String', handles.clusterSizeThreshold);
    catch
        set(hObject, 'String', '5');
    end
    return
end
handles.clusterSizeThreshold=tmp;

guidata(hObject, handles);
CallBack_allIntensityRadio(hObject, eventdata, 'c');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% save image push button
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_saveImagePush(hObject, eventdata, thisfilename, isMask) 
handles = guidata(hObject);

if exist('thisfilename')
    if ~strcmp(deblank(thisfilename), '')
        if isfield(handles,'imageFileName')
            if ~isempty(handles.imageFileName)
                mni2mask(cell2mat(handles.currentDisplayMNI'), thisfilename, cell2mat(handles.currentDisplayIntensity), handles.M{1}, handles.DIM{1}, handles.imageFileName{1});
                return;
            end
        end
        mni2mask(cell2mat(handles.currentDisplayMNI'), thisfilename, cell2mat(handles.currentDisplayIntensity), handles.M{1}, handles.DIM{1});
        return;
    end
end

[filename, pathname] = uiputfile('*.img', 'Save image file as', get(handles.saveImageFileEdit, 'string'));
if isequal(filename,0) | isequal(pathname,0)
   return
else
   thisfilename = fullfile(pathname, filename);
end

if isfield(handles,'imageFileName')
    if ~isempty(handles.imageFileName{1})
        mni2mask(cell2mat(handles.currentDisplayMNI'), thisfilename, cell2mat(handles.currentDisplayIntensity'), handles.M{1}, handles.DIM{1}, handles.imageFileName{1}, isMask);
        return;
    end
end

mni2mask(cell2mat(handles.currentDisplayMNI'), thisfilename, cell2mat(handles.currentDisplayIntensity), handles.M{1}, handles.DIM{1}, '', isMask);

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% save image edit
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_saveImageFileEdit(hObject, eventdata) 
handles = guidata(hObject);
CallBack_saveImagePush(handles.saveImagePush, eventdata, get(hObject,'string'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% save result push
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_saveResultPSPush(hObject, eventdata, thisfilename) 

handles = guidata(hObject);

if exist('thisfilename')
    if ~strcmp(deblank(thisfilename), '')
        spm_print(handles.figure);
        if strcmp(handles.system, 'linux')
            system(['ps2pdf spm2.ps']);
            system(['mv spm2.ps ' thisfilename]);
            [p,f,ext]=fileparts(thisfilename);
            system(['mv spm2.pdf ' fullfile(p,f) '.pdf']);            
        elseif strcmp(handles.system, 'windows')
            system(['move spm2.ps ' '"' thisfilename '"']);
		end        
        return;
    end
end

[filename, pathname] = uiputfile('*.ps', 'Save result as', get(handles.saveResultPSEdit, 'string'));
if isequal(filename,0) | isequal(pathname,0)
   return
else
   thisfilename = fullfile(pathname, filename);
end


% print
H  = findobj(get(handles.figure,'Children'),'flat','Type','axes');
un = cellstr(get(H,'Units'));
pos = get(H,'position');
index = [];

for ii=1:length(H)
    if findstr('pixels', un{ii})
        continue;
    end
    if pos{ii}(1)>0.4 & pos{ii}(2) > 0.5
        set(H(ii),'position',[pos{ii}(1), pos{ii}(2), pos{ii}(3)*3/4, pos{ii}(4)]);
        index = [index ii];
    end
end


spm_print(handles.figure);
if strcmp(handles.system, 'linux')
    system(['ps2pdf spm2.ps']);
    system(['mv spm2.ps ' thisfilename]);
    [p,f,ext]=fileparts(thisfilename);
    system(['mv spm2.pdf ' fullfile(p,f) '.pdf']);            
elseif strcmp(handles.system, 'windows')
    system(['move spm2.ps ' '"' thisfilename '"']);
end 

% set the position back
set(H(index), {'position'}, pos(index)); 


% printstr = ['print -dpsc2 -painters -noui ' '''' thisfilename ''''];
% try
%     orient portrait
%     eval(printstr);
%     printsuccess = 1;
% catch
%     errordlg('Print to ps file failed', 'print error');
%     printsuccess = 0;
% end
%set(H,{'Units'},un);
% if strcmp(handles.system, 'linux') & printsuccess == 1
%     system(['ps2pdf ' '''' thisfilename '''']);
% end

return;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% save result edit
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_saveResultPSEdit(hObject, eventdata) 
handles = guidata(hObject);
CallBack_saveResultPSPush(hObject, eventdata, get(handles.saveResultPSEdit,'string'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% select a cluster
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_selectThisClusterPush(hObject, eventdata) 
handles = guidata(hObject);

try
    xyz = handles.currentxyz';
catch
    xyz = spm_XYZreg('GetCoords',handles.hReg);
end
try
    handles.selectedCluster = [handles.selectedCluster; xyz'];
catch
    handles.selectedCluster = xyz';
end

set(handles.figure,'currentaxes', handles.glassViewAxes);
xrange = xlim;
yrange = ylim;
try
    delete(handles.hGlassText)
end
%handles.hGlassText = text(xrange(1)+diff(xrange)*0.6, yrange(1)+diff(yrange)*0.9, [num2str(size(handles.selectedCluster,1)) ' clusters selected']);
set(handles.infoTextBox, 'string', [num2str(size(handles.selectedCluster,1)) ' clusters selected']);
guidata(hObject, handles);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% unselect a cluster
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CallBack_clearSelectedClusterPush(hObject, eventdata) 
handles = guidata(hObject);
try
    handles = rmfield(handles,'selectedCluster');
end
set(handles.figure,'currentaxes', handles.glassViewAxes);
xrange = xlim;
yrange = ylim;
try
    %delete(handles.hGlassText)
    set(handles.infoTextBox, 'string', ['No clusters selected']);
end

guidata(hObject, handles);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% pick a cluster
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%