Changeset 177 for vnmr_recon/epi_recon.m
 Timestamp:
 Sep 15, 2011, 10:13:39 AM (8 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

vnmr_recon/epi_recon.m
r170 r177 68 68 end 69 69 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') 73 103 nRef = 3; 74 else 104 elseif strcmpi(Dat.EPI_PC,'pointwise') 75 105 nRef = 1; 76 106 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 77 115 78 116 % Number of segments … … 82 120 nSeg = 1; 83 121 end 122 123 84 124 %procpar.image = [1 1 0 2 1 1 1 1 1]; 85 125 % Get index to first reference image … … 87 127 ref_start = ref_start(1); 88 128 nDummys = 0; 89 if isempty(ref_start) 129 if isempty(ref_start) && doPhaseCorrection 90 130 kspace = []; 91 131 data = []; … … 121 161 end 122 162 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 123 171 % Calculate indexes for navigator and data echoes 124 172 if nNav > 0 … … 134 182 135 183 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'); 139 187 end 140 188 … … 148 196 kspace = reshape(kspace,np,nv+nNav*nSeg,ns,nRcvrs,nVols+nRef+nDummys); 149 197 198 150 199 ref1_ind = find(procpar.image==0,1); 151 200 ref2_ind = find(procpar.image==2,1); 152 201 im1_ind = find(procpar.image==1,1); 153 202 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 190 255 end 191 256 end … … 220 285 % Flip even lines 221 286 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 227 298 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); 236 306 237 307 if Dat.EPIPhasedArrayData
Note: See TracChangeset
for help on using the changeset viewer.