Changeset 152 for vnmr_recon/epi_recon.m
 Timestamp:
 Jan 18, 2011, 1:04:28 PM (8 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

vnmr_recon/epi_recon.m
r128 r152 6 6 % this code reconstructs 7 7 if nargin==0 8 kspace = {'epi_se_rapid_sp3','epi','epi_fMRI' };8 kspace = {'epi_se_rapid_sp3','epi','epi_fMRI','epip'}; 9 9 return 10 10 end … … 13 13 msg_out = ''; 14 14 15 % EPI data is measured with INOVA system and sorted using CORREPI 15 % Number of receivers 16 if isfield(procpar,'rcvrs') && ~isempty(procpar.rcvrs) 17 nRcvrs = length(procpar.rcvrs{1}=='y'); 18 else 19 nRcvrs = 1; 20 end 21 22 23 24 % ===================================================== 25 % Reconstruct EPI sequences meaured with the old INOVA 26 % ===================================================== 16 27 if isfield(procpar,'readres') && isfield(procpar,'phaseres') 17 28 … … 31 42 kspace = reshape(kspace,[size(kspace,1) size(kspace,2) ... 32 43 tmp_ns size(kspace,3)/tmp_ns]); 44 45 46 % =================================== 47 % Reconstruct VNMRj and EPIP sequences 48 % =================================== 33 49 else 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 = (ii1)*blksz+1:nVols; 167 else 168 vol_inds = (ii1)*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 35 206 end 207 208 209 210 211 212 213 214 215 216 217 218
Note: See TracChangeset
for help on using the changeset viewer.