function [xXa,Sessa,Ka,Pa,nscana,rowa] = spm_bch_tsampl3(xX,Sess,K,P,nscan,row)
% SPM batch system: bch mode (only) non regular temporal sampling
% FORMAT [xXa,Sessa,Ka,Pa,nscana,rowa] = spm_bch_tsampl3(xX,Sess,K,P,nscan,row)
%
% xX     - 
% Sess   - 
% K      - 
% P      - 
% nscan  - 
% row    - 
% xXa    - 
% Sessa  - 
% Ka     - 
% nscana - 
% rowa   - 
%_______________________________________________________________________
% @(#)spm_bch_tsampl.m	2.8 Stephanie Rouquette 00/01/25
%       spm_bch_tsampl3 modified by Jeff Cooper from spm_bch_tsampl2  10/4/02
% 
% -----------------------------------------------------------------
% -----------------------------------------------------------------
% Notes on operation
% -----------------------------------------------------------------
% -----------------------------------------------------------------
% For most purposes, spm_bch_tsampl3 operates just like its predecessors,
% spm_bch_tsampl and spm_bch_tsampl2.  It's called from spm_fmri_spm_ui.m,
% and it functions to reshape the design matrix according to some 
% specifications.
%
% If it's called in non-batch mode, there is no difference from
% spm_bch_tsampl2 - the code requests a vector of scans to omit for each
% session, by hand.  If called in batch mode, the code examines the 
% 'model' batch structure (where things like nscans, nsessions, filters, etc.
% are specified) for a field called 'omitted'; each entry of omitted (each omitted(i)) is
% a cell containing a vector of scans to omit for that entry's session (the ith session).
% In other words, if your experiment had five sessions, and you 
% wanted to throw out the 1st and 5th in each session, you would edit
% your batch file containing your model structure so it looked like this:
%
% model(1) = struct( ...
% 'types',          4, ...  		% 1-'spec'; 2-'review'; 3-'est'; 4-'spec & est'
% 'global_effects', {'None'}, ...   	% 'Scaling|None'
%  .
%  .
%  etc.
%  .
%  .
% 'omitted',        {{[1 5] [1 5] [1 5] [1 5] [1 5]}}, ...
%  .
%  .
%  etc.
% )
%
% Note the double curly brace surrounding the list of vectors; this 
% distributes the vectors each to their own cell.
%
% If your model doesn't contain this field, the default is that
% no scans will be thrown out.

%=======================================================================
% Programmers Guide
%=======================================================================
% Batch system implemented on this routine. See spm_bch.man
% If inputs are modified in this routine, try to modify spm_bch.man
% and spm_bch_bchmat (if necessary) accordingly. 
% Calls to spm_input in this routine use the BCH gobal variable.  
%    BCH.bch_mat 
%    BCH.index0  = {'normalize',index_of_Analysis};
% or 
%    BCH.index0  = {'normalisation',index_of_Analysis}; (when
%                   spm_spn3d is launched for edit_defaults 
%modified 9/26/02 -sue

global BCH;

[xXa,Sessa,Ka,Pa,nscana,rowa] = deal(xX,Sess,K,P,nscan,row);

if ~isempty(BCH) | spm_input('Specify omitted scans?','+1','y/n')=='y'

	xXa.X = [];
	Pa = '';
	compt1 = 0;
	compt2 = 0;
    omitted = cell(1,length(Sess));
	for s = 1:length(Sess)
		if ~isempty(BCH)
                    omitted{s} = spm_input('batch',{},'omitted',s);
                    if isempty(omitted{s})
                        sample = spm_input('batch',{},'remain',s);
                    else
                        sample = 1:length(Sess{s}.row);
                        sample(omitted{s}) = [];
                    end
                else
                   omitted{s} = spm_input('Enter omitted scans','+1','n');
                   sample      = 1:length(Sess{s}.row);
                   sample(omitted{s}) = [];
                end
		nscana(s) = length(sample);
		xXa.X = [xXa.X;xX.X(compt1+sample,:)];
		Pa = [Pa;P(compt1+sample,:)];
		rowa{s} = (compt2+(1:nscana(s)))';
		Sessa{s}.row = rowa{s};
        % next, repair pst - very naive way, but this is a
        % crude tool to start with.  We don't want to make
        % it seem as though these events didn't happen -
        % they've just had chunks taken out of them.
        for i = 1:length(Sessa{s}.pst)
            Sessa{s}.pst{i}(omitted{s}) = [];
        end
        
        % repair sf - delete missing bins
        nbins = length(Sess{s}.sf{1}) / nscan(s);  % number of bins per scan
        omitted{s} = sort(omitted{s});
        for i = 1:length(Sessa{s}.sf)
            omitted_to_toss = omitted{s};
            for j = 1:length(omitted_to_toss)
                miss_bin_start = ((omitted_to_toss(j)-1)*nbins) + 1;
                miss_bin_end = (miss_bin_start + nbins) - 1;
                Sessa{s}.sf{i}(miss_bin_start:miss_bin_end) = [];
                % this next line is to re-index omitted -
                % now that we've taken out nbins from sf, we
                % need to adjust omitted to reflect the fact
                % that we have one less scan.
                omitted_to_toss = omitted_to_toss - 1; % can trash omitted because this is the last thing to fix
            end
        end
		Ka{s}.row = rowa{s};	
                KL = K{s}.KL(sample,sample);

                % Renormalize
                k = length(sample);
                KL = spdiags(1./sum(KL')',0,k,k)*KL;

                Ka{s}.KL = KL;
                Ka{s}.KH = Ka{s}.KH(sample,:);

		compt1 = compt1+nscan(s);
		compt2 = compt2+nscana(s);
	end % for s = 1:length(Sess)
end

