Changeset 177 for vnmr_recon/epi_recon.m


Ignore:
Timestamp:
Sep 15, 2011, 10:13:39 AM (8 years ago)
Author:
tjniskan
Message:
  • Added support for reading complex data from bruker 2dseq files.
  • Fixed scaling of 3D data from 2dseq files.
  • Aedes now prompts for handling of complex data (absolute, real

or imag parts) if only complex data is read from a file.

  • Added an option for EPI phase correction in aedes_readvnmr and epi_recon

M aedes_readbruker.m
M aedes_readvnmr.m
M aedes.m
M aedes_revision.m
M vnmr_recon/epi_recon.m

File:
1 edited

Legend:

Unmodified
Added
Removed
  • vnmr_recon/epi_recon.m

    r170 r177  
    6868        end
    6969       
    70         % Number of reference images
    71         if isfield(procpar,'epiref_type') && ~isempty(procpar.epiref_type) && ...
    72                         strcmpi(procpar.epiref_type,'triple')
     70        % EPI phase correction
     71        doPhaseCorrection = true;
     72        if strcmpi(Dat.EPI_PC,'auto')
     73                % Select phase correction using procpar
     74                if isfield(procpar,'epi_pc')
     75                        switch lower(procpar.epi_pc{1})
     76                                case 'triple_ref'
     77                                        nRef = 3;
     78                                case 'scaled_triple'
     79                                        % Issue a warning because unscaled triple ref is used anyway...
     80                                        warning(['Using unscaled TRIPLE REFERENCE phase correction although "%s" ',...
     81                                                'is indicated in PROCPAR.EPI_PC.'],procpar.epi_pc{1});
     82                                        nRef = 3;
     83                                case {'linear','quadratic',...
     84                                                'center_pair','pairwise','first_pair'}
     85                                        % Issue a warning because pointise correction is used anyway...
     86                                        warning(['Using POINTWISE phase correction although "%s" ',...
     87                                                'is indicated in PROCPAR.EPI_PC.'],procpar.epi_pc{1});
     88                                        nRef = 1;
     89                                case 'pointwise'
     90                                        nRef = 1;
     91                                otherwise
     92                                        warning('Unknown EPI_PC "%s". Skipping phase correction,',procpar.epi_pc{1})
     93                        end
     94                else
     95                        % Don't use phase correction
     96                        warning('Cannot find EPI_PC parameter in PROCPAR. Skipping phase correction...');
     97                        doPhaseCorrection = false;
     98                end
     99        elseif strcmpi(Dat.EPI_PC,'off')
     100                nRef = 0;
     101                doPhaseCorrection = false;
     102        elseif strcmpi(Dat.EPI_PC,'triple')
    73103                nRef = 3;
    74         else
     104        elseif strcmpi(Dat.EPI_PC,'pointwise')
    75105                nRef = 1;
    76106        end
     107       
     108%       % Number of reference images
     109%       if isfield(procpar,'epiref_type') && ~isempty(procpar.epiref_type) && ...
     110%                       strcmpi(procpar.epiref_type,'triple')
     111%               nRef = 3;
     112%       else
     113%               nRef = 1;
     114%       end
    77115       
    78116        % Number of segments
     
    82120                nSeg = 1;
    83121        end
     122               
     123       
    84124        %procpar.image = [1 1 0 -2 -1 1 1 1 1];
    85125        % Get index to first reference image
     
    87127        ref_start = ref_start(1);
    88128        nDummys = 0;
    89         if isempty(ref_start)
     129        if isempty(ref_start) && doPhaseCorrection
    90130                kspace = [];
    91131                data = [];
     
    121161        end
    122162       
     163        % Detect partial fourier EPI
     164        if isfield(procpar,'fract_ky') && procpar.fract_ky~=nv/2
     165                nv_full = nv;
     166                nv = nv/2+procpar.fract_ky;
     167        else
     168                nv_full = nv;
     169        end
     170       
    123171        % Calculate indexes for navigator and data echoes
    124172        if nNav > 0
     
    134182       
    135183        if Dat.EPIPhasedArrayData
    136                 data = zeros(nv,np,ns,nVols,nRcvrs,'single');
    137         else
    138                 data = zeros(nv,np,ns,nVols,'single');
     184                data = zeros(nv_full,np,ns,nVols,nRcvrs,'single');
     185        else
     186                data = zeros(nv_full,np,ns,nVols,'single');
    139187        end
    140188       
     
    148196        kspace = reshape(kspace,np,nv+nNav*nSeg,ns,nRcvrs,nVols+nRef+nDummys);
    149197       
     198       
    150199        ref1_ind = find(procpar.image==0,1);
    151200        ref2_ind = find(procpar.image==-2,1);
    152201        im1_ind = find(procpar.image==-1,1);
    153202       
    154         ref_data = kspace(:,DataInd,:,:,[ref1_ind ref2_ind im1_ind]);
    155         if ~isempty(pe_table)
    156                 ref_data(:,pe_table,:,:,:) = ref_data;
    157         end
    158         % Flip even kspace lines
    159         ref_data(:,2:2:end,:,:,:) = flipdim(ref_data(:,2:2:end,:,:,:),1);
    160         %kspace(:,:,:,:,1:nRef)=[];
    161        
    162         % Calculate phase corrections for all receivers
    163         phase_e = ref_data(:,:,:,:,1);
    164         rev_phase = ref_data(:,:,:,:,1);
    165         for ii=1:nRcvrs
    166                 if nRef==1
    167       %ref_im = ref_data(:,:,:,ii,ref1_ind);
    168                         ref_im = ref_data(:,:,:,ii,1);
    169       phase1 = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
    170                         rev_phase = 1;
    171                         phase_e(:,:,:,ii,:) = phase1;
    172                 elseif nRef==3
    173                         %ref1_ind = find(procpar.image==0,1);
    174                         %ref1 = ref_data(:,:,:,ii,ref1_ind);
    175                         ref1 = ref_data(:,:,:,ii,1);
    176                         phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
    177                        
    178                         %ref2_ind = find(procpar.image==-2,1);
    179                         %ref2 = flipdim(ref_data(:,:,:,ii,ref2_ind),1);
    180                         ref2 = flipdim(ref_data(:,:,:,ii,2),1);
    181                         phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
    182                        
    183                         %im1_ind = find(procpar.image==-1,1);
    184                         %im1 = flipdim(ref_data(:,:,:,ii,im1_ind),1);
    185                         im1 = flipdim(ref_data(:,:,:,ii,3),1);
    186                        
    187                        
    188                         rev_phase(:,:,:,ii,:) = fft(im1,[],1).*phase2;
    189                         phase_e(:,:,:,ii,:) = phase1;
     203        if ( isempty(ref2_ind) | isempty(im1_ind) ) && nRef == 3
     204                warning(['Triple reference phase correction was requested, but ',...
     205                        ' < 3 reference images were found. Skipping phase correction.'])
     206                doPhaseCorrection = false;
     207        end
     208               
     209        if doPhaseCorrection
     210
     211                ref_data = kspace(:,DataInd,:,:,[ref1_ind ref2_ind im1_ind]);
     212                if ~isempty(pe_table)
     213                        ref_data(:,pe_table,:,:,:) = ref_data;
     214                end
     215                % Flip even kspace lines
     216                ref_data(:,2:2:end,:,:,:) = flipdim(ref_data(:,2:2:end,:,:,:),1);
     217                %kspace(:,:,:,:,1:nRef)=[];
     218               
     219                % Calculate phase corrections for all receivers
     220                phase_e = ref_data(:,:,:,:,1);
     221                rev_phase = ref_data(:,:,:,:,1);
     222                for ii=1:nRcvrs
     223                        if nRef==1
     224                                %ref_im = ref_data(:,:,:,ii,ref1_ind);
     225                                ref_im = ref_data(:,:,:,ii,1);
     226                               
     227                                phase1 = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
     228                                rev_phase = 1;
     229                                phase_e(:,:,:,ii,:) = phase1;
     230                               
     231                                %tmp = fft(ref_im,[],1);
     232                                %tmp_max = repmat(max(tmp),[size(ref_im,1),1,1]);
     233                                %phase1 = conj(tmp./tmp_max);
     234                               
     235                                phase_e(:,:,:,ii,:) = phase1;
     236                        elseif nRef==3
     237                                %ref1_ind = find(procpar.image==0,1);
     238                                %ref1 = ref_data(:,:,:,ii,ref1_ind);
     239                                ref1 = ref_data(:,:,:,ii,1);
     240                                phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
     241                               
     242                                %ref2_ind = find(procpar.image==-2,1);
     243                                %ref2 = flipdim(ref_data(:,:,:,ii,ref2_ind),1);
     244                                ref2 = flipdim(ref_data(:,:,:,ii,2),1);
     245                                phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
     246                               
     247                                %im1_ind = find(procpar.image==-1,1);
     248                                %im1 = flipdim(ref_data(:,:,:,ii,im1_ind),1);
     249                                im1 = flipdim(ref_data(:,:,:,ii,3),1);
     250                               
     251                               
     252                                rev_phase(:,:,:,ii,:) = fft(im1,[],1).*phase2;
     253                                phase_e(:,:,:,ii,:) = phase1;
     254                        end
    190255                end
    191256        end
     
    220285                % Flip even lines
    221286                tmp_kspace(:,2:2:end,:,:,:) = flipdim(tmp_kspace(:,2:2:end,:,:,:),1);
    222 
    223                 for kk=1:nRcvrs
    224                         if nRef==3
    225                                 for tt=1:length(vol_inds)
    226                                         tmp_kspace(:,:,:,kk,tt) = ifft(rev_phase(:,:,:,kk,:)+fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:));
     287               
     288                if doPhaseCorrection
     289                        for kk=1:nRcvrs
     290                                if nRef==3
     291                                        for tt=1:length(vol_inds)
     292                                                tmp_kspace(:,:,:,kk,tt) = ifft(rev_phase(:,:,:,kk,:)+fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:));
     293                                        end
     294                                else
     295                                        for tt=1:length(vol_inds)
     296                                                tmp_kspace(:,:,:,kk,tt) = ifft(fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:));
     297                                        end
    227298                                end
    228                         else
    229                                 for tt=1:length(vol_inds)
    230                                         tmp_kspace(:,:,:,kk,tt) = ifft(fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:));
    231                                 end
    232                         end
    233                         tmp_data = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
    234                 end
    235                
     299                        end
     300                end
     301               
     302                if nv~=nv_full
     303                        tmp_kspace(:,nv_full,:,:,:) = single(0);
     304                end
     305                tmp_data = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
    236306               
    237307                if Dat.EPIPhasedArrayData
Note: See TracChangeset for help on using the changeset viewer.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008