par;
defaults=spm_get_defaults;
%spm_defaults;
flags = defaults.coreg;
resFlags = struct(...
        'interp', 1,...                       % trilinear interpolation
        'wrap', flags.write.wrap,...           % wrapping info (ignore...)
        'mask', flags.write.mask,...           % masking (see spm_reslice)
        'which',1,...                         % write reslice time series for later use, don't write the first 1
        'mean',0);                            % do write mean image
par;
for sb =1:length(PAR.subjects) % for each subject 
    clear P;
    clear SPM;
    sprintf('Registering T1 and segmented WM CSF to BOLD data space for subject #%g''s data...%g subjects left...',sb,length(PAR.subjects)-sb)
    
    
        t1img=spm_select('FPList', PAR.structdir{sb}, ['^' PAR.subjects{sb} '.*_t1_struct\.img$']);
        if isempty(t1img)
            t1img=spm_select('FPList', PAR.structdir{sb}, ['^' PAR.subjects{sb} '.*_t1_struct\.nii$']);
        end
        v1=spm_vol(t1img);
        dat=spm_read_vols(v1);
        segt1=spm_select('FPList', PAR.structdir{sb}, ['^c\w*' PAR.subjects{sb} '.*_t1_struct\.nii$']);
        fname=fullfile(spm_str_manip(t1img, 'h'), 't1_noskull.nii');
        if ~exist(fname, 'file')
            vs=spm_vol(segt1(1:3,:));
            msk=spm_read_vols(vs);
            msk=squeeze(sum(msk,4));
            msk=msk>0.3;
            dat=dat.*msk;
            vo=v1;
            vo.fname=fname;
            vo=spm_write_vol(vo,dat);
        end
        PF=fullfile(spm_str_manip(t1img, 'h'), 't1_noskull.nii');
        P=[];
        for c=1:PAR.ncond
            if isempty(PAR.condirs{sb,c}), continue; end
            PO=spm_select('FPList', PAR.condirs{sb,c} , ['^rbk_c.*$']);
            if ~isempty(PO), 
                fprintf('rbk_c images have been found in %s\n', PAR.condirs{sb,c});
                continue; 
            end;
            PF1=fullfile(PAR.condirs{sb,c}, ['bk_' spm_str_manip(PF(1,:),'t')]);
            cpstr=['!cp '  PF ' ' PF1  ];
            eval(cpstr);    % make a copy of the T1
            PF=PF1;
            if size(PF,1)>1, PG=PF(1,:);end
            VF = spm_vol(PF);
            PO=segt1;
            % make a copy of the segmented maps since we will do registration of T1 to mean ASL image for each of the 2 PASL and each of the pCASL sequence separately.
            if size(PO,1)<3, 
               fprintf('Only %d segmented maps were found for %s ', size(PO,1), PAR.structdir{sb,nses});
               continue;
            end
            for i=1:3
                 cpstr=['!cp ' PO(i,:) ' ' fullfile(PAR.condirs{sb,c} , ['bk_' spm_str_manip(PO(i,:),'t')]) ];
                 eval(cpstr);
            end
            PO=spm_select('FPList', PAR.condirs{sb,c} , ['^bk_c.*\.nii$']);
            PG=spm_select('FPList', PAR.condirs{sb,c} , ['^meana' PAR.subjects{sb} '.*img$']);
            if isempty(PG)
                PG=spm_select('FPList', PAR.condirs{sb,c} , ['^meana' PAR.subjects{sb} '.*nii$']);
            end
%             PG =spm_select('FPList', PAR.condirs{sb,s,c}, ['^EPI_noskull\.nii$']);
            
            if size(PG,1)>2,PG=PG(size(PG,1),:);end
       
      
            if isempty(PG), 
                fprintf('No mean images found for %s\n', PAR.condirs{sb,c});
                continue; 
            end
            PG=PG(1,:);
            VG = spm_vol(PG(1,:));
            if size(VG,1)>1, VG=VG(2);end;
            x  = spm_coreg(VG, VF,flags.estimate);

            %get the transformation to be applied with parameters 'x'
            M  = inv(spm_matrix(x));
            %transform the mean image
            %spm_get_space(deblank(PG),M);
            %in MM we put the transformations for the 'other' images
            MM = zeros(4,4,size(PO,1));
            for j=1:size(PO,1),
                %get the transformation matrix for every image
                MM(:,:,j) = spm_get_space(deblank(PO(j,:)));
            end
            for j=1:size(PO,1),
                %write the transformation applied to every image
                %filename: deblank(PO(j,:))
                spm_get_space(deblank(PO(j,:)), M*MM(:,:,j));
            end
            P=strvcat(PG,PO);
            spm_reslice_flt(P,resFlags);
            eval(['!rm -f '    '/bk_*.nii']);
    
        end
   
end
