function artdetect4(varargin)
% FORMAT artdetect4
% 
%
% artdetect 4 is the latest version of a script intended to allow 
% visual inspection of MRI data and automated or by-hand removal of 
% 'outlier' scans, whose intensities are radically different from the
% mean of the timeseries.  "Radically different" is a user-set threshold
% of standard deviations from the mean which is set interactively in 
% the graphical viewing window.  The script offers several choices of
% removal methods - insertion of the global mean, interpolation from 
% surrounding scans or simple removal from the timeseries. 
% ----------------------------------------------------------------------
% Jeff Cooper, Nov. 2002
%
% enhanced from artdetect3.m, by Sue Whitfield + Max.


% -----------------------
% Initialize, begin loop
% -----------------------

pfig = [];
%intfig = [];
%do_more = 1;
%while do_more == 1;

% ------------------------
% Default values
% ------------------------

default_thresh = 2;
sd = default_thresh;


% ------------------------
% Collect files
% ------------------------

num_sess = spm_input('How many sessions?',1,'n',1,1);

M = [];
P = [];
for i = 1:num_sess
    P{i} = spm_get(Inf,'.img',['Select data images for session'  num2str(i) ':']);
    
    mvmt_file = spm_get(1,'.txt',['Select movement params file for session' num2str(i) ':']);
    load(mvmt_file);
    [mv_path,mv_name,mv_ext] = fileparts(mvmt_file);

    dot = findstr(mv_name, '.');
    if ~(isempty(dot))
       mv_name = mv_name(1:(dot - 1));
    end
    M{i} = eval(mv_name);
end

global_type_flag = spm_input('Which global mean to use?', 1, 'b', 'Regular | Every Vox | Masked', [1 2 3], 1);

if global_type_flag==3
    mask = spm_get(1, '.img', 'Select mask image in functional space');
    [maskY,maskXYZmm] = spm_read_vols(spm_vol(mask));
    maskXYZmm = maskXYZmm(:,find(maskY==max(max(max(maskY)))));
end

drop_flag = spm_input('Drop 1st scan of each session?', '+1', 'y/n', [1 0], 2);

if drop_flag
    for i = 1:num_sess
        currP = P{i};
        currP(1,:) = [];
        P{i} = currP;
        currM = M{i};
        currM(1,:) = [];
        M{i} = currM;
    end
end


P = char(P);
mv_data = [];
for i = 1:length(M)
    mv_data = vertcat(mv_data,M{i});
end



% -------------------------
% get file identifiers and Global values
% -------------------------

fprintf('%-4s: ','Mapping files...')                                  
VY     = spm_vol(P);
fprintf('%3s\n','...done')                                          

if any(any(diff(cat(1,VY.dim),1,1),1)&[1,1,1,0])
	error('images do not all have the same dimensions')
end

nscans = size(P,1);
%keyboard;
% ------------------------
% Compute Global variate
%--------------------------

%GM     = 100;
g      = zeros(nscans,1);

fprintf('%-4s: %3s','Calculating globals...',' ')
if global_type_flag==1  % regular mean
    for i  = 1:nscans  
	    g(i) = spm_global(VY(i));
    end
elseif global_type_flag==2  % every voxel
    for i = 1:nscans
        g(i) = mean(mean(mean(spm_read_vols(VY(i)))));
    end
else %masked mean
    [dummy, XYZmm] = spm_read_vols(VY(1));
    vinv = inv(VY(1).mat);
    [dummy, idx_to_mask] = intersect(XYZmm', maskXYZmm', 'rows');
    for i = 1:nscans
        Y = spm_read_vols(VY(i));
        Y(idx_to_mask) = [];
        g(i) = mean(Y);
    end
end
    
    
fprintf('%s%3s\n','...done\n')
if global_type_flag==3
    fprintf('\n%g voxels were in the mask.\n', length(idx_to_mask))
end


% ------------------------
% Compute default out indices
% ------------------------

out_idx = (find(abs(zscore(g)) > sd))';


% -----------------------
% Draw initial figure
% -----------------------


figure('Units', 'normalized', 'Position', [0.2 0.2 0.6 0.7]);
rng = range(g);
pfig = gcf;

subplot(5,1,1);
plot(g);
%xlabel(['artifact index list [' int2str(out_idx') ']'], 'FontSize', 8, 'Color','r');
ylabel(['Range = ' num2str(rng)], 'FontSize', 8);
title('Global Mean');

subplot(5,1,2);
%thresh_axes = gca;
%set(gca, 'Tag', 'threshaxes');
plot(abs(zscore(g)));
ylabel('std away from mean');

thresh_x = 1:nscans;
thresh_y = sd*ones(1,nscans);
line(thresh_x, thresh_y, 'Color', 'r');

axes_lim = get(gca, 'YLim');
axes_height = axes_lim(1):1:axes_lim(2);
for i = 1:length(out_idx)
    line((out_idx(i)*ones(1, length(axes_height))), axes_height, 'Color', 'r');
end

subplot(5,1,3);
plot(mv_data(:,1:3));
ylabel('movement in mm');
xlabel('scans');
legend('x mvmt', 'y mvmt', 'z mvmt',0);
h = gca;
set(h,'Ygrid','on');

subplot(5,1,4);
plot(mv_data(:,4:6));
ylabel('movement in rad');
xlabel('scans');
y_lim = get(gca, 'YLim');
legend('pitch', 'roll', 'yaw', ['max:' num2str(y_lim(2))], ['min:' num2str(y_lim(1))],0);
h = gca;
set(h,'Ygrid','on');

%keyboard;
h_rangetext = uicontrol(gcf, 'Units', 'characters', 'Position', [10 10 18 2],...
        'String', 'Range of data is: ', 'Style', 'text', ...
        'HorizontalAlignment', 'left', ...
        'BackgroundColor', [0.8 0.8 0.8]);
h_rangenum = uicontrol(gcf, 'Units', 'characters', 'Position', [29 10 10 2], ...
        'String', num2str(rng), 'Style', 'text', ...
        'HorizontalAlignment', 'left',...
        'Tag', 'rangenum',...
        'BackgroundColor', [0.8 0.8 0.8]);
h_threshtext = uicontrol(gcf, 'Units', 'characters', 'Position', [25 8 18 2],...
        'String', 'Current threshold is: ', 'Style', 'text', ...
        'HorizontalAlignment', 'left', ...
        'BackgroundColor', [0.8 0.8 0.8]);
h_threshnum = uicontrol(gcf, 'Units', 'characters', 'Position', [44 8 10 2],...
        'String', num2str(sd), 'Style', 'text', ...
        'HorizontalAlignment', 'left', ...
        'BackgroundColor', [0.8 0.8 0.8],...
        'Tag', 'threshnum');
h_indextext = uicontrol(gcf, 'Units', 'characters', 'Position', [10 3 15 2],...
        'String', 'Outlier indices: ', 'Style', 'text', ...
        'HorizontalAlignment', 'left', ...
        'BackgroundColor', [0.8 0.8 0.8], ...
        'ForegroundColor', 'r');
h_indexedit = uicontrol(gcf, 'Units', 'characters', 'Position', [25 3.25 40 2],...
        'String', int2str(out_idx), 'Style', 'edit', ...
        'HorizontalAlignment', 'left', ...
        'Callback', 'artdetect_edit',...
        'BackgroundColor', [0.8 0.8 0.8],...
        'Tag', 'indexedit');
h_indexinst = uicontrol(gcf, 'Units', 'characters', 'Position', [66 3 40 2],...
        'String', '[Hit return to update after editing]', 'Style', 'text',...
        'HorizontalAlignment', 'left', ...
        'BackgroundColor', [0.8 0.8 0.8]);
h_repair = uicontrol(gcf, 'Units', 'characters', 'Position', [10 1 10 2],...
        'String', 'Repair', 'Style', 'pushbutton', ...
        'Tooltipstring', 'Repair data by eliminating outliers', ...
        'Callback', 'artdetect_repair');
h_up = uicontrol(gcf, 'Units', 'characters', 'Position', [10 8 10 2],...
        'String', 'Up', 'Style', 'pushbutton', ...
        'TooltipString', 'Raise threshold for outliers', ...
        'Callback', 'thresh_up');
h_down = uicontrol(gcf, 'Units', 'characters', 'Position', [10 6 10 2],...
        'String', 'Down', 'Style', 'pushbutton', ...
        'TooltipString', 'Lower threshold for outliers', ...
        'Callback', 'thresh_down');

guidata(gcf, g);
setappdata(h_repair,'data',P);




% --------------------------
% Finish loop, close up figures
% --------------------------

%keyboard;
%do_more = spm_input('Session done!',1,'b','New Sess|Quit Now',[1,2]);

%if do_more==2
%    if ~isempty(intfig)
%        close(intfig);
%    end
%end

%if ~isempty(pfig)
%    close(pfig);
%    pfig = [];
%end

%end






% Things to add
% - scaling (see code below)?
% - masking (see code below)?


% scale if specified (otherwise session specific grand mean scaling)
%-----------------------------------------------------------------------
%gSF     = GM./g;
%if strcmp(Global,'None')
	%for i = 1:nsess
	%	j      = row{i};
	%	gSF(j) = GM./mean(g(j));
	%end
%end

%-Apply gSF to memory-mapped scalefactors to implement scaling
%-----------------------------------------------------------------------
%for  i = 1:q, VY(i).pinfo(1:2,:) = VY(i).pinfo(1:2,:)*gSF(i); end


%-Masking structure
%-----------------------------------------------------------------------
%xM     = struct('T',	ones(q,1),...
		%'TH',	g.*gSF,...
	%	'I',	0,...
	%	'VM',	{[]},...
	%	'xs',	struct('Masking','analysis threshold'));