%program to isolate images which fall outside 2 std from the global mean.
%sue whitfield
%********Max modification***********
% modified by Max 10/1/02 to allow cutoff specification and interpolation.
% currently not recommended for multiple consecutive bad images. 
% Fixed plots to include useful information for interpolation decision making
% as well as the option to view a movie of a selected slice from images before repair.
% There is a second place, after selecting "fix it", to exit without interpolating
% any images by entering a 0 in the "edit cutoff vectors" box.
% suppressed unnecessary command line output.
infig = [];
pfig = [];
mfig = [];
do_more = 1;
while do_more == 1;
[cl, ypos] = spm_input_ui('Specify cutoff?(default = 2)','-1','y/n',[1,0],2);
infig = gcf;
      if cl == 1
         [sd, ypos] = spm_input_ui('Enter a value!(1-3)','1','r', ' ', 1, [1 3]);
      else
         sd = 2;
      end
%********Max modification 1 end************
P = spm_get(Inf,'.img','Select data images:');
% 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
%if any(any(any(diff(cat(3,VY.mat),1,3),3)))
%   error('images do not all have same orientation & voxel size'), end

nscan = size(P,1);
%-Compute Global variate
%-----------------------------------------------------------------------
GM     = 100;
q      = sum(nscan);
g      = zeros(q,1);
fprintf('%-4s: %3s','Calculating globals...',' ')                     %-#
for i  = 1:q
	%fprintf('%s%30s',sprintf('\b')*ones(1,30),sprintf('%4d/%-4d',i,q)) %-#
	g(i) = spm_global(VY(i));
end
fprintf('%s%3s\n','...done')               %-#

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

%plot motion params from realignment

sub = spm_get(1,'.txt',' Select movement params file:');
load(sub);
[fpath,fname,fext] = fileparts(sub);


%sub = input('which subject? ', 's');
%load (['c:\data\' sub '\realignment_params_aI_001.txt']);
data1 = eval(fname);

%*****************modification 2************
rng = range(g);
pfig = gcf;
subplot(4,1,1);
plot(g);
%plotliny (mean(g), 'k:');
plot(g);
xlabel(['artifact index list [' int2str(out_idx') ']'], 'FontSize', 8, 'Color','r');
ylabel(['Range = ' num2str(rng)], 'FontSize', 8);
title('Global Mean');

subplot(4,1,2);
plot(abs(zscore(g)));
%set(gca, 'XTick', [out_idx']);
%xlabel('Location of indices above threshold', 'FontSize', 8)
ylabel('std away from mean');
%title('Z-scored Deviation');
%plotliny (sd, 'k:');
%plot(abs(zscore(g)));

subplot(4,1,3);
%page = figure('NumberTitle', 'off', 'PaperOrientation', 'landscape', 'PaperPosition', [0 0 11 8.5], 'Units', 'inches', 'Position', [0 0 11 8.5]);
 plot(data1(:,1:3));

%plotline(14,'k:');
%plotline(320,'k:');
ylabel('movement in mm');
xlabel('scans');
h = gca;
set(h,'Ygrid','on');

subplot(4,1,4);
%page = figure('NumberTitle', 'off', 'PaperOrientation', 'landscape', 'PaperPosition', [0 0 11 8.5], 'Units', 'inches', 'Position', [0 0 11 8.5]);
 plot(data1(:,4:6));

%plotline(14,'k:');
%plotline(320,'k:');
ylabel('movement in rad');
xlabel('scans');
h = gca;
set(h,'Ygrid','on');

fix_it = spm_input_ui(['Range of data is ' num2str(rng) '. Fix it?'],'1','y/n',[1,0],2);
if fix_it == 1%chance to interpolate
    
    % this mod and the next - the movie - removed by Jeff 10/15/02 because, frankly,
    % they don't work.  At least without max_slice_means.
    % -----------------------------------------------------------------------------
    
    %bslice = [];
    %find_it = spm_input_ui(['Find the bad slices?'],'1','y/n',[1,0],2);
    %if find_it == 1%chance to find the bad slices within the offending imgs
    %    tot_Bs = [];
    %    sln = size(out_idx,1);
    %   %sample the offending images + a min subset of 10
    %       if sln < 5
    %           bad_idx = [out_idx; random('unid',254,10,1)];
    %       else
    %            bad_idx = [out_idx; random('unid',254,(size(out_idx,1)*2),1)]; 
    %       end
    %       sxln = size(bad_idx,1)
    %   for B_s = 1:sxln
    %       nV = spm_vol(P(bad_idx,:));
    %       Bs = max_slice_means(nV(B_s));
    %       tot_Bs = [tot_Bs; Bs];%img x slice matrix to pass to zscore
    %   end
    %    [bsr,bsc] = find(abs(zscore(tot_Bs)) > sd);
    %    
    %       bslice = mode(bsc);
       
       
    %end
    %mov = 1;
    %while mov == 1
    %      mov = spm_input_ui(['See a movie before editing vectors?'],'1','y/n',[1,0],2);
    %      if mov == 1
    %         m_type = spm_input_ui('Pick your poison!',1,'b','All Imgs|Pick a Set|The Bad',[1,2,3])
    %         if~isempty(bslice)
    %             slice = bslice;
    %         else
    %          slice = round(VY(1).dim(3)*.66);% calculates the axial slice at the top 3rd point
    %         end
    %     
    %          if m_type == 1
    %             mfig = spm_movie('Load',1,P,slice,0,1);% change the second arg for diff planes
    %                                                    % 1 = axial 2 = saggital 3 = coronal
    %          elseif m_type == 2
    %             P_2 = spm_get(Inf,'.img',[]);
    %             mfig = spm_movie('Load',1,P_2,slice,0,1);
    %          elseif m_type == 3
    %             P_3 =  P(out_idx,:);
    %             mfig = spm_movie('Load',1,P_3,slice,0,1);
    %          end
    %      end
    %end
    
    % ----------------------------------------------------------------------------------------------
    % Jeff's removal ends
    
    sd_mod = spm_input_ui(['Edit cutoff vectors.(0 for none)'],'1','e',int2str(out_idx'),inf);
    if sd_mod == 0
        do_more = spm_input_ui('No images to fix!',1,'b','Keep working|Quit Now',[1,2]);
        if do_more == 2
           close all
           break
        end
    else
        for j = 1:size(sd_mod, 1)
            if sd_mod(j,:) == 1 
                if (size(P,1) < 3)
                    error('Cannot interpolate first image: less than three scans.');
                end
                im_in = spm_vol([P(2,:);P(3,:)]);
                im_map = spm_vol(P(1,:));
                spm_imcalc(im_in,im_map,'(i1+i2)/2');
                continue;
            else 
                if sd_mod(j,:) == size(P,1)
                    im_in = spm_vol([P(sd_mod(j,:)-2,:);P(sd_mod(j,:)-1,:)]);
                    im_map = spm_vol(P(sd_mod(j,:),:));
                    spm_imcalc(im_in,im_map,'(i1+i2)/2');
                    continue;
                end
            end
            im_in = spm_vol([P(sd_mod(j,:)-1,:);P(sd_mod(j,:)+1,:)]);
            im_map = spm_vol(P(sd_mod(j,:),:));
            spm_imcalc(im_in, im_map, '(i1+i2)/2');
        end
    end
end
do_more = spm_input_ui('Session done!',1,'b','New Sess|Quit Now',[1,2]);
        if do_more == 1
            if ~isempty(pfig)
               close(pfig)
               pfig = [];
            end
            if ~isempty(mfig)
               close(mfig)
               mfig = [];
            end
            figure(infig)
        elseif do_more == 2
           if ~isempty(pfig)
               close(pfig)
           end
           if ~isempty(infig)
               close(infig)
           end
        end
end
%************Max modification 2 end*******************
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%PUT ME BACK IN~%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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'));