function spike_repair
% FORMAT spike_repair 
% 
% SUMMARY
%   Reads a set of images and writes a new set of images 
%   after filtering the data for noise.  User is asked to choose the
%   repair method or methods. This program is best applied
%   to the raw images, so the cleaned output images can be fed into
%   slicetiming or realignment.
%
%   No original data is removed- the old images are all preserved.
%
% User is asked via GUI to supply:
%   - Set of images
%   - A repair method
% POSSIBLE REPAIRS
% 1. Filter all the images with a 3-point median filter in time. An 
%     excellent filter for TR=2 or less.  If TR > 2, this filter
%     may clobber activations in event-related experiments.
%     Adds a prefix "f" to the cleaned output images.
% 2. Detect and repair bad slices. Derives new values for the bad slice
%     using linear interpolation of the before and after volumes.
%     Bad slices are detected by the amount of data scattered outside
%     the head, relative to the usual amount for the slice. (The usual
%     amount differs for each slice, and is determined automatically.)
%     It removes outliers, but also reduces activations, so it's safest 
%     when fewer than 5% of the slices are cleaned up. 
%     Adds a prefix "g" to the cleaned output images.
% 3. Eliminate data outside the head. Generates a head mask automatically,
%     and writes it as ArtifactMask.img. Sets voxels outside the head 
%     to zero. This process may help realignment to be more accurate.
%     Adds a prefix "h" to the cleaned output images.
% 4. Combinations of methods 1 and 3, or 2 and 3.
%
% NOTES
%  1. Use biac_movie to compare the input and repaired images.
%  2. First and last volumes will not be filtered.
% Paul Mazaika - September 2004


% ----------------------
% Interrogate GUI
% ----------------------

% handles = guihandles;
% rng = str2num(get(handles.rangenum, 'String'));
% out_idx = round(str2num(get(handles.indexedit, 'String')));
% P = getappdata(gcbo, 'data');

% PARAMETERS
spm_defaults;
OUTSLICE =  5;  %  Threshold above sample means to filter slices
   %  15 is very visible on contrast image. 8 is slightly visible.
% THREEE CONTROL OPTIONS
%   Mask the images    - Set amask = 1.
%   Median filter all the images  - Set allfilt = 1.
%   Detect and repair bad slices  - Set slicefilt = 1.

Pimages  = spm_get(Inf,'.img','select images');
P = spm_vol(Pimages);

repair_flag = spm_input('Which repair methods to use?', 1, 'm', ...
    ' Repair Bad Slices ( best choice for TR > 2 sec ) | Median Filter All Data  ( best choice for TR = 2 sec or less ) | Eliminate data outside head | Repair Bad Slices and Eliminate data outside head | Median Filter All Data and Eliminate data outside head',...
              [1 2 3 4 5], 4);
if ( repair_flag == 3 ) amask = 1; allfilt = 0; slicefilt = 0; end
if ( repair_flag == 1 ) amask = 0; allfilt = 0; slicefilt = 1; end
if ( repair_flag == 2 ) amask = 0; allfilt = 1; slicefilt = 0; end
if ( repair_flag == 4 ) amask = 1; allfilt = 0; slicefilt = 1; end
if ( repair_flag == 5 ) amask = 1; allfilt = 1; slicefilt = 0; end

%  Set the prefix correctly on the written files
if  ( amask ==1 & allfilt == 0 & slicefilt == 0 ) prechar = 'h'; end
if  ( amask ==1 & allfilt == 1 & slicefilt == 0 ) prechar = 'fh'; end
if  ( amask ==1 & allfilt == 0 & slicefilt == 1 ) prechar = 'gh'; end
if  ( amask ==0 & allfilt == 1 & slicefilt == 0 ) prechar = 'f'; end
if  ( amask ==0 & allfilt == 0 & slicefilt == 1 ) prechar = 'g'; end

fprintf('\n NEW IMAGE FILES WILL BE CREATED');
fprintf('\n The filtered scan data will be saved in the same directory');
fprintf('\n with %s pre-pended to their filenames.\n',prechar);

if amask == 1 | slicefilt == 1 %  Automask options
    fprintf('\n Generated mask image is written to file ArtifactMask.img.\n');
    fprintf('\n');
    Pnames = P(1).fname;
    Automask = biac_automask(Pnames(1,:),-1,1);
    maskcount = sum(sum(sum(Automask)));  %  Number of voxels in mask.
    voxelcount = prod(size(Automask));    %  Number of voxels in 3D volume.
end

if amask == 1
    fprintf('\n ELIMINATE DATA OUTSIDE HEAD');
    fprintf('\n All scans will have voxels outside the head set to zero.\n');
end

if slicefilt == 1  % Prepare some thresholds for slice testing.
    fprintf('\n INTERPOLATE THROUGH BAD SLICES\n');
    % Find the slice orientation used to collect the data
    [ vx, vy, vz ] = size(Automask);
    orient = 0;
    if ( vx < vy & vx < vz ) orient = 1; disp(' Remove bad Sagittal slices'); end
    if ( vy < vx & vy < vz ) orient = 2; disp(' Remove bad Coronal slices'); end
    if ( vz < vx & vz < vy ) orient = 3; disp(' Remove bad Axial slices'); end
    if ( orient == 0 ) 
        disp('Cannot determine slice orientation for bad slice filtering.')
        return; 
    end
    % Find 3 samples of slice baseline activity outside the head.
    p = zeros(3,min([vx vy vz]));
    for i = 1:3
        Y1 = spm_read_vols(P(i));
        Y1 = ( 1 - Automask ).*Y1;
        % Get the plane sums perpendicular to orient direction
        if ( orient == 1 ) p(i,:) = mean(mean(Y1,3),2); end 
        if ( orient == 2 ) p(i,:) = mean(mean(Y1,1),3); end 
        if ( orient == 3 ) p(i,:) = mean(mean(Y1,1),2); end 
    end
    % Select a low value for each plane, and set threshold a bit above it.
    pq = 0.5*( min(p) + median(p,1));
    pq = pq + OUTSLICE;  % pq array is the threshold test for bad slices.
    fprintf('\n Interpolating new values for bad slices when ');
    fprintf('\n average value outside head is %3d counts over baseline.\n',OUTSLICE);
    fprintf('\n');
end
     
if allfilt == 1
    fprintf('\n MEDIAN FILTER ALL THE DATA');
	fprintf('\n All scans are being filtered by median 3-point time filter.');
	fprintf('\n This is safe to do when TR = 2. Otherwise, use the');
	fprintf('\n    Remove Bad Slices option instead.\n');
    fprintf('\n');
end



% Main Loop - Filter everything but the first and last volumes.
nscans = size(P,1);
Y4(1,:,:,:) = spm_read_vols(P(1));  % rows vary fastest
Y4(2,:,:,:) = spm_read_vols(P(2));
for i = 3:nscans
        Y4(3,:,:,:) = spm_read_vols(P(i));
        if allfilt == 1  % Filter all data by median.
            Yn = median(Y4,1);
            Yn2 = squeeze(Yn);
        end
        if slicefilt == 1  % Repair bad slices by linear interpolation.
            % Check if outside head value is over pq. If so, filter.
            Yn2 = squeeze(Y4(2,:,:,:));
            Y = ( 1 - Automask).*Yn2;
            if ( orient == 1 ) 
                py = mean(mean(Y,3),2); 
                for j = 1:vx
                    if py(j) > pq(j)  % slice needs to be filtered.
                        Yn2(j,:,:) = squeeze((Y4(1,j,:,:) + Y4(3,j,:,:))/2.0);
                        fprintf('   Fixed sagittal. Vol. %d, Slice %d.\n',i-1,j);
                    end  
                end 
            end
            if ( orient == 2 )
                py = mean(mean(Y,1),3); 
                for j = 1:vy
                    if py(j) > pq(j)  % slice needs to be filtered.
                        Yn2(:,j,:) = squeeze((Y4(1,:,j,:) + Y4(3,:,j,:))/2.0);
                        fprintf('   Fixed coronal. Vol %d, Slice %d.\n',i-1,j);
                    end  
                end 
            end 
            if ( orient == 3 )
                py = mean(mean(Y1,1),2);
                for j = 1:vz
                    if py(j) > pq(j)  % slice needs to be filtered.
                        Yn2(:,:,j) = squeeze((Y4(1,:,:,j) + Y4(3,:,:,j))/2.0);
                        fprintf('   Fixed axial. Vol %d, Slice %d.\n',i-1,j);
                    end  
                end 
            end    
        end
        if  ( allfilt == 0 & slicefilt == 0 )  % Head mask only.
            Yn2 = squeeze(Y4(2,:,:,:));
        end
        % Yn is a 4D file including the smoothed Y2 now.
        if ( amask ) Yn2 = Yn2.*Automask; end
        % Prepare the header for the filtered volume.
           V = spm_vol(P(i-1).fname);
           v = V;
           [dirname, sname, sext ] = fileparts(V.fname);
           sfname = [ prechar, sname ];
           filtname = fullfile(dirname,[sfname sext]);
           v.fname = filtname;
        spm_write_vol(v,Yn2);  
        showprog = [' Corrected volume   ', sname, sext ];
        disp(showprog); 
        if i == 3   % Write unfiltered first scan
            Yn1 = squeeze(Y4(1,:,:,:)); 
            if ( amask ) Yn1 = Yn1.*Automask; end
            [dirname, sname, sext ] = fileparts(P(1).fname);
            sfname = [ prechar, sname ];
            filtname = fullfile(dirname,[sfname sext]);
            v.fname = filtname;
            spm_write_vol(v,Yn1); 
        end
        if i == nscans  % Write unfiltered last scan.
            Yn1 = squeeze(Y4(3,:,:,:));
            if ( amask ) Yn1 = Yn1.*Automask; end
            Vlast = spm_vol(P(nscans).fname);
            [dirname, sname, sext ] = fileparts(P(nscans).fname);
            sfname = [ prechar, sname ];
            filtname = fullfile(dirname,[sfname sext]);
            v.fname = filtname;
            spm_write_vol(v,Yn1); 
        end
        % Slide the read volumes window up.
        Y4(1,:,:,:) = Y4(2,:,:,:);
        Y4(2,:,:,:) = Y4(3,:,:,:); 
end
disp('Done!')


