Changeset 152 for vnmr_recon/epi_recon.m


Ignore:
Timestamp:
Jan 18, 2011, 1:04:28 PM (8 years ago)
Author:
tjniskan
Message:
  • aedes_readvnmr.m almost works with epi data...
  • The ROI averages of fMRI time series are now shown in one column instead of two in plugins/fmri_plugins/basic_fmri_analysis.m

M aedes_readvnmr.m
M plugins/fmri_plugins/basic_fmri_analysis.m
M aedes_revision.m
M vnmr_recon/epi_recon.m

File:
1 edited

Legend:

Unmodified
Added
Removed
  • vnmr_recon/epi_recon.m

    r128 r152  
    66% this code reconstructs
    77if nargin==0
    8   kspace = {'epi_se_rapid_sp3','epi','epi_fMRI'};
     8  kspace = {'epi_se_rapid_sp3','epi','epi_fMRI','epip'};
    99  return
    1010end
     
    1313msg_out = '';
    1414
    15 % EPI data is measured with INOVA system and sorted using CORREPI
     15% Number of receivers
     16if isfield(procpar,'rcvrs') && ~isempty(procpar.rcvrs)
     17        nRcvrs = length(procpar.rcvrs{1}=='y');
     18else
     19        nRcvrs = 1;
     20end
     21
     22
     23
     24% =====================================================
     25% Reconstruct EPI sequences meaured with the old INOVA
     26% =====================================================
    1627if isfield(procpar,'readres') && isfield(procpar,'phaseres')
    1728 
     
    3142  kspace = reshape(kspace,[size(kspace,1) size(kspace,2) ...
    3243    tmp_ns size(kspace,3)/tmp_ns]);
     44       
     45       
     46        % ===================================
     47        % Reconstruct VNMRj and EPIP sequences
     48        % ===================================
    3349else
    34  
     50       
     51        if strncmpi(procpar.seqfil,'epip',4) && isfield(procpar,'nphase') && ...
     52                        isfield(procpar,'nread')
     53                isEPIP = true;
     54        else
     55                isEPIP = false;
     56        end
     57       
     58        % Number of navigator echos
     59        if isfield(procpar,'navigator') && strcmpi(procpar.navigator,'y')
     60                nNav = length(procpar.nav_echo);
     61        else
     62                nNav = 0;
     63        end
     64       
     65        % Number of reference images
     66        if isfield(procpar,'epiref_type') && ~isempty(procpar.epiref_type) && ...
     67                        strcmpi(procpar.epiref_type,'triple')
     68                nRef = 3;
     69        else
     70                nRef = 1;
     71        end
     72       
     73        % Number of volumes
     74        nVols = length(procpar.image)-nRef;
     75       
     76        % Number of slices
     77        ns = procpar.ns;
     78       
     79        % Phase sampling
     80        if isEPIP
     81                nv = procpar.nphase;
     82               
     83                % Read sampling
     84                np = procpar.nread/2;
     85        else
     86                nv = procpar.nv;
     87                np = procpar.np/2;
     88        end
     89       
     90        % Get orientation
     91        orient = procpar.orient{1}
     92       
     93        % Remove navigators
     94        if nNav~=0
     95                nav_data = kspace(1:(np*nNav),:);
     96                kspace(1:(np*nNav),:)=[];
     97        else
     98                nav_data = [];
     99        end
     100       
     101        if Dat.EPIPhasedArrayData
     102                data = zeros(nv,np,ns,nVols,nRcv,'single');
     103        else
     104                data = zeros(nv,np,ns,nVols,'single');
     105        end
     106       
     107        % Reshape kspace and detach reference images
     108        if isEPIP
     109                kspace = reshape(kspace,np,nv,nRcvrs,ns,nVols+nRef);
     110        else
     111                kspace = reshape(kspace,np,nv,nRcvrs,ns,nVols+nRef);
     112        end
     113       
     114        % Flip even kspace lines
     115        if ~Dat.FastDataRead
     116                for tt=2:2:size(kspace,2)
     117                        kspace(:,tt,:,:,:) = flipdim(kspace(:,tt,:,:,:),1);
     118                end
     119        else
     120                kspace(:,2:2:end,:,:,:) = flipdim(kspace(:,2:2:end,:,:,:),1);
     121        end
     122        ref_data = kspace(:,:,:,:,1:nRef);
     123        kspace(:,:,:,:,1:nRef)=[];
     124       
     125        % Calculate phase corrections for all receivers
     126        phase_e = ref_data(:,:,:,:,1);
     127        rev_phase = ref_data(:,:,:,:,1);
     128        for ii=1:nRcvrs
     129                if nRef==1
     130      ref_im = ref_data(:,:,ii,:,1);
     131      phase1 = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
     132                        rev_phase = 1;
     133                        phase_e(:,:,ii,:,:) = phase1;
     134                elseif nRef==3
     135                        ref1_ind = find(procpar.image==0,1);
     136                        ref1 = ref_data(:,:,ii,:,ref1_ind);
     137                        phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
     138                       
     139                        ref2_ind = find(procpar.image==-2,1);
     140                        ref2 = flipdim(ref_data(:,:,ii,:,ref2_ind),1);
     141                        phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
     142                       
     143                        im1_ind = find(procpar.image==-1,1);
     144                        im1 = flipdim(ref_data(:,:,ii,:,im1_ind),1);
     145                       
     146                        rev_phase(:,:,ii,:,:) = fft(im1,[],1).*phase2;
     147                        phase_e(:,:,ii,:,:) = phase1;
     148                end
     149        end
     150       
     151        % Reconstruct in blocks
     152        kssz=size(kspace);
     153        blksz = Dat.EPIBlockSize; % Process EPI data in 100 volume blocks (default)
     154        nBlocks = ceil(nVols/blksz);
     155        lnum = length(num2str(nBlocks));
     156        lnumstr = num2str(lnum);
     157        bsl = lnum*2+1;
     158        fprintf(1,'Processing data in blocks of %d volumes\n',blksz)
     159        fprintf(1,['Processing block...%0',lnumstr,'d/%0',lnumstr,'d'],1,nBlocks);
     160       
     161        for ii=1:nBlocks
     162                fprintf(1,repmat('\b',1,bsl));
     163                fprintf(1,['%0',lnumstr,'d/%0',lnumstr,'d'],ii,nBlocks);
     164                tmp_data = [];
     165                if ii==nBlocks
     166                        vol_inds = (ii-1)*blksz+1:nVols;
     167                else
     168                        vol_inds = (ii-1)*blksz+1:ii*blksz;
     169                end
     170                tmp_kspace = kspace(:,:,:,:,vol_inds);
     171               
     172                for kk=1:nRcvrs
     173                        if nRef==3
     174                                for tt=1:length(vol_inds)
     175                                        tmp_kspace(:,:,kk,:,tt) = ifft(rev_phase(:,:,kk,:,:)+fft(kspace(:,:,kk,:,vol_inds(tt)),[],1).*phase_e(:,:,kk,:,:));
     176                                end
     177                        else
     178                                for tt=1:length(vol_inds)
     179                                        tmp_kspace(:,:,kk,:,tt) = ifft(fft(kspace(:,:,kk,:,vol_inds(tt)),[],1).*phase_e(:,:,kk,:,:));
     180                                end
     181                        end
     182                        tmp_data = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
     183                end
     184               
     185               
     186                if Dat.EPIPhasedArrayData
     187                        data_block = abs(tmp_data);
     188                        data_block = permute(data_block,[1 2 4 5 3]);
     189                        % Permute to correct orientation
     190                        if strcmpi(orient,'trans90')
     191                                data_block = permute(data_block,[2 1 3 4 5]);
     192                        end
     193                        data(:,:,:,vol_inds,:) = data_block;
     194                else
     195                        data_block = sqrt(sum(tmp_data.*conj(tmp_data),3));
     196                        % Permute to correct orientation
     197                        if strcmpi(orient,'trans90')
     198                                data_block = permute(data_block,[2 1 3 4 5]);
     199                        end
     200                        data(:,:,:,vol_inds) = squeeze(data_block);
     201                end
     202               
     203        end
     204        fprintf(1,'\n')
     205       
    35206end
     207       
     208
     209
     210
     211
     212
     213
     214
     215
     216
     217
     218
Note: See TracChangeset for help on using the changeset viewer.

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