[128] | 1 | function [kspace,data,msg_out]=epi_recon(kspace,Dat,procpar) |
---|
| 2 | % This is a custom VNMR k-space reconstruction code for EPI data used by |
---|
| 3 | % aedes_readvnmr. |
---|
| 4 | |
---|
| 5 | % If called without input arguments, return the sequence names that |
---|
| 6 | % this code reconstructs |
---|
| 7 | if nargin==0 |
---|
[170] | 8 | kspace = {'epi_se_rapid_sp3','epi','epi_fMRI','epip','epip_fMRI','epiT1rho'}; |
---|
[128] | 9 | return |
---|
| 10 | end |
---|
| 11 | |
---|
| 12 | data=[]; |
---|
| 13 | msg_out = ''; |
---|
| 14 | |
---|
[152] | 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 | % ===================================================== |
---|
[128] | 27 | if isfield(procpar,'readres') && isfield(procpar,'phaseres') |
---|
| 28 | |
---|
| 29 | % Number of slices |
---|
| 30 | tmp_ns=length(procpar.pss); |
---|
| 31 | |
---|
| 32 | if isfield(procpar,'navecho') && strcmpi(procpar.navecho{1},'y') |
---|
| 33 | tmp_nv = procpar.nv-procpar.nseg; |
---|
| 34 | else |
---|
| 35 | tmp_nv = procpar.nv; |
---|
| 36 | end |
---|
| 37 | kspace = reshape(kspace,[size(kspace,1) ... |
---|
| 38 | size(kspace,2)/tmp_nv tmp_nv]); |
---|
| 39 | kspace = permute(kspace,[1 3 2]); |
---|
| 40 | |
---|
| 41 | % Reshape to 4D matrix |
---|
| 42 | kspace = reshape(kspace,[size(kspace,1) size(kspace,2) ... |
---|
| 43 | tmp_ns size(kspace,3)/tmp_ns]); |
---|
[152] | 44 | |
---|
| 45 | |
---|
[153] | 46 | data = abs(fftshift(fftshift(fft(fft(kspace,[],1),[],2),1),2)); |
---|
| 47 | data = permute(data,[2 1 3 4]); |
---|
| 48 | data = flipdim(flipdim(data,1),2); |
---|
| 49 | data(:,:,:,1) = []; % Remove reference image |
---|
| 50 | |
---|
[152] | 51 | % =================================== |
---|
| 52 | % Reconstruct VNMRj and EPIP sequences |
---|
| 53 | % =================================== |
---|
[128] | 54 | else |
---|
[152] | 55 | |
---|
[169] | 56 | if isfield(procpar,'nphase') && ... |
---|
[152] | 57 | isfield(procpar,'nread') |
---|
| 58 | isEPIP = true; |
---|
| 59 | else |
---|
| 60 | isEPIP = false; |
---|
| 61 | end |
---|
| 62 | |
---|
| 63 | % Number of navigator echos |
---|
| 64 | if isfield(procpar,'navigator') && strcmpi(procpar.navigator,'y') |
---|
| 65 | nNav = length(procpar.nav_echo); |
---|
| 66 | else |
---|
| 67 | nNav = 0; |
---|
| 68 | end |
---|
| 69 | |
---|
[177] | 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') |
---|
[152] | 103 | nRef = 3; |
---|
[177] | 104 | elseif strcmpi(Dat.EPI_PC,'pointwise') |
---|
[152] | 105 | nRef = 1; |
---|
| 106 | end |
---|
| 107 | |
---|
[177] | 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 |
---|
| 115 | |
---|
[162] | 116 | % Number of segments |
---|
| 117 | if isfield(procpar,'nseg') && ~isempty(procpar.nseg) |
---|
| 118 | nSeg = procpar.nseg; |
---|
| 119 | else |
---|
| 120 | nSeg = 1; |
---|
| 121 | end |
---|
[177] | 122 | |
---|
| 123 | |
---|
[170] | 124 | %procpar.image = [1 1 0 -2 -1 1 1 1 1]; |
---|
| 125 | % Get index to first reference image |
---|
| 126 | ref_start = find(procpar.image == 0); |
---|
| 127 | ref_start = ref_start(1); |
---|
| 128 | nDummys = 0; |
---|
[177] | 129 | if isempty(ref_start) && doPhaseCorrection |
---|
[170] | 130 | kspace = []; |
---|
| 131 | data = []; |
---|
| 132 | msg_out = 'Reference image(s) not found.'; |
---|
| 133 | return |
---|
| 134 | elseif ref_start(1) > 1 |
---|
| 135 | nDummys = ref_start-1; |
---|
| 136 | warning('Rejecting %d volumes before first reference image as dummy scans...',nDummys); |
---|
| 137 | end |
---|
[162] | 138 | |
---|
[152] | 139 | % Number of volumes |
---|
[170] | 140 | nVols = length(procpar.image(nDummys+1:end))-nRef; |
---|
[152] | 141 | |
---|
| 142 | % Number of slices |
---|
| 143 | ns = procpar.ns; |
---|
| 144 | |
---|
| 145 | if isEPIP |
---|
[153] | 146 | nv = procpar.nphase; % Phase sampling |
---|
| 147 | np = procpar.nread/2; % Read sampling |
---|
[152] | 148 | else |
---|
| 149 | nv = procpar.nv; |
---|
| 150 | np = procpar.np/2; |
---|
| 151 | end |
---|
| 152 | |
---|
| 153 | % Get orientation |
---|
[153] | 154 | orient = procpar.orient{1}; |
---|
[152] | 155 | |
---|
[177] | 156 | % Detect partial fourier EPI |
---|
| 157 | if isfield(procpar,'fract_ky') && procpar.fract_ky~=nv/2 |
---|
| 158 | nv_full = nv; |
---|
| 159 | nv = nv/2+procpar.fract_ky; |
---|
[192] | 160 | elseif isEPIP && isfield(procpar,'kzero') && procpar.kzero~=0 |
---|
| 161 | nv_full = nv; |
---|
| 162 | nv = nv/2+procpar.kzero; |
---|
[177] | 163 | else |
---|
| 164 | nv_full = nv; |
---|
| 165 | end |
---|
| 166 | |
---|
[192] | 167 | if isEPIP && (nv_full==64 & np==64) |
---|
| 168 | % The information about navigator echoes in the procpar doesn't seem to |
---|
| 169 | % be consistent. This code tried to guess the number of navigators from |
---|
| 170 | % the k-space size |
---|
| 171 | nNav = size(kspace,1)/np-nv; |
---|
| 172 | end |
---|
| 173 | |
---|
[162] | 174 | % Calculate indexes for navigator and data echoes |
---|
| 175 | if nNav > 0 |
---|
| 176 | ind = reshape([1:(nv+nNav*nSeg)],[],nSeg); |
---|
| 177 | NavInd = ind(1:nNav,:); |
---|
| 178 | DataInd = ind(nNav+1:end,:); |
---|
| 179 | NavInd = NavInd(:).'; |
---|
| 180 | DataInd = DataInd(:).'; |
---|
[192] | 181 | |
---|
| 182 | % Warn that navigator echoes are not used... |
---|
| 183 | warning('Navigator correction has not been implemented. All navigator echoes are ignored.'); |
---|
[152] | 184 | else |
---|
[162] | 185 | DataInd = [1:nv]; |
---|
| 186 | NavInd = []; |
---|
[152] | 187 | end |
---|
| 188 | |
---|
| 189 | if Dat.EPIPhasedArrayData |
---|
[177] | 190 | data = zeros(nv_full,np,ns,nVols,nRcvrs,'single'); |
---|
[152] | 191 | else |
---|
[177] | 192 | data = zeros(nv_full,np,ns,nVols,'single'); |
---|
[152] | 193 | end |
---|
| 194 | |
---|
[162] | 195 | % Look for phase table |
---|
| 196 | pe_table = []; |
---|
| 197 | if isfield(Dat,'phasetable') && ~isempty(Dat.phasetable) |
---|
| 198 | pe_table = Dat.phasetable; |
---|
| 199 | end |
---|
| 200 | |
---|
[152] | 201 | % Reshape kspace and detach reference images |
---|
[170] | 202 | kspace = reshape(kspace,np,nv+nNav*nSeg,ns,nRcvrs,nVols+nRef+nDummys); |
---|
[152] | 203 | |
---|
[177] | 204 | |
---|
[170] | 205 | ref1_ind = find(procpar.image==0,1); |
---|
| 206 | ref2_ind = find(procpar.image==-2,1); |
---|
| 207 | im1_ind = find(procpar.image==-1,1); |
---|
| 208 | |
---|
[177] | 209 | if ( isempty(ref2_ind) | isempty(im1_ind) ) && nRef == 3 |
---|
| 210 | warning(['Triple reference phase correction was requested, but ',... |
---|
| 211 | ' < 3 reference images were found. Skipping phase correction.']) |
---|
| 212 | doPhaseCorrection = false; |
---|
[162] | 213 | end |
---|
[177] | 214 | |
---|
| 215 | if doPhaseCorrection |
---|
| 216 | |
---|
| 217 | ref_data = kspace(:,DataInd,:,:,[ref1_ind ref2_ind im1_ind]); |
---|
| 218 | if ~isempty(pe_table) |
---|
| 219 | ref_data(:,pe_table,:,:,:) = ref_data; |
---|
[152] | 220 | end |
---|
[177] | 221 | % Flip even kspace lines |
---|
| 222 | ref_data(:,2:2:end,:,:,:) = flipdim(ref_data(:,2:2:end,:,:,:),1); |
---|
| 223 | %kspace(:,:,:,:,1:nRef)=[]; |
---|
| 224 | |
---|
| 225 | % Calculate phase corrections for all receivers |
---|
| 226 | phase_e = ref_data(:,:,:,:,1); |
---|
| 227 | rev_phase = ref_data(:,:,:,:,1); |
---|
| 228 | for ii=1:nRcvrs |
---|
| 229 | if nRef==1 |
---|
| 230 | %ref_im = ref_data(:,:,:,ii,ref1_ind); |
---|
| 231 | ref_im = ref_data(:,:,:,ii,1); |
---|
| 232 | |
---|
| 233 | phase1 = exp(-sqrt(-1)*angle(fft(ref_im,[],1))); |
---|
| 234 | rev_phase = 1; |
---|
| 235 | phase_e(:,:,:,ii,:) = phase1; |
---|
| 236 | |
---|
| 237 | %tmp = fft(ref_im,[],1); |
---|
| 238 | %tmp_max = repmat(max(tmp),[size(ref_im,1),1,1]); |
---|
| 239 | %phase1 = conj(tmp./tmp_max); |
---|
| 240 | |
---|
| 241 | phase_e(:,:,:,ii,:) = phase1; |
---|
| 242 | elseif nRef==3 |
---|
| 243 | %ref1_ind = find(procpar.image==0,1); |
---|
| 244 | %ref1 = ref_data(:,:,:,ii,ref1_ind); |
---|
| 245 | ref1 = ref_data(:,:,:,ii,1); |
---|
| 246 | phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1))); |
---|
| 247 | |
---|
| 248 | %ref2_ind = find(procpar.image==-2,1); |
---|
| 249 | %ref2 = flipdim(ref_data(:,:,:,ii,ref2_ind),1); |
---|
| 250 | ref2 = flipdim(ref_data(:,:,:,ii,2),1); |
---|
| 251 | phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1))); |
---|
| 252 | |
---|
| 253 | %im1_ind = find(procpar.image==-1,1); |
---|
| 254 | %im1 = flipdim(ref_data(:,:,:,ii,im1_ind),1); |
---|
| 255 | im1 = flipdim(ref_data(:,:,:,ii,3),1); |
---|
| 256 | |
---|
| 257 | |
---|
| 258 | rev_phase(:,:,:,ii,:) = fft(im1,[],1).*phase2; |
---|
| 259 | phase_e(:,:,:,ii,:) = phase1; |
---|
| 260 | end |
---|
| 261 | end |
---|
[152] | 262 | end |
---|
| 263 | |
---|
| 264 | % Reconstruct in blocks |
---|
| 265 | kssz=size(kspace); |
---|
| 266 | blksz = Dat.EPIBlockSize; % Process EPI data in 100 volume blocks (default) |
---|
| 267 | nBlocks = ceil(nVols/blksz); |
---|
| 268 | lnum = length(num2str(nBlocks)); |
---|
| 269 | lnumstr = num2str(lnum); |
---|
| 270 | bsl = lnum*2+1; |
---|
| 271 | fprintf(1,'Processing data in blocks of %d volumes\n',blksz) |
---|
| 272 | fprintf(1,['Processing block...%0',lnumstr,'d/%0',lnumstr,'d'],1,nBlocks); |
---|
| 273 | |
---|
| 274 | for ii=1:nBlocks |
---|
| 275 | fprintf(1,repmat('\b',1,bsl)); |
---|
| 276 | fprintf(1,['%0',lnumstr,'d/%0',lnumstr,'d'],ii,nBlocks); |
---|
| 277 | tmp_data = []; |
---|
| 278 | if ii==nBlocks |
---|
[170] | 279 | vol_inds = (ii-1)*blksz+1+nRef+nDummys:nVols+nRef+nDummys; |
---|
[152] | 280 | else |
---|
[170] | 281 | vol_inds = (ii-1)*blksz+1+nRef+nDummys:ii*blksz+nRef+nDummys; |
---|
[152] | 282 | end |
---|
| 283 | |
---|
[162] | 284 | % At the moment navigator echoes are not used... |
---|
| 285 | %tmp_kspace = kspace(:,nNav+1:end,:,:,vol_inds); |
---|
| 286 | tmp_kspace = kspace(:,DataInd,:,:,vol_inds); |
---|
| 287 | if ~isempty(pe_table) |
---|
| 288 | tmp_kspace(:,pe_table,:,:,:) = tmp_kspace; |
---|
| 289 | end |
---|
| 290 | |
---|
[153] | 291 | % Flip even lines |
---|
| 292 | tmp_kspace(:,2:2:end,:,:,:) = flipdim(tmp_kspace(:,2:2:end,:,:,:),1); |
---|
[177] | 293 | |
---|
| 294 | if doPhaseCorrection |
---|
| 295 | for kk=1:nRcvrs |
---|
| 296 | if nRef==3 |
---|
| 297 | for tt=1:length(vol_inds) |
---|
| 298 | tmp_kspace(:,:,:,kk,tt) = ifft(rev_phase(:,:,:,kk,:)+fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:)); |
---|
| 299 | end |
---|
| 300 | else |
---|
| 301 | for tt=1:length(vol_inds) |
---|
| 302 | tmp_kspace(:,:,:,kk,tt) = ifft(fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:)); |
---|
| 303 | end |
---|
[152] | 304 | end |
---|
| 305 | end |
---|
| 306 | end |
---|
| 307 | |
---|
[177] | 308 | if nv~=nv_full |
---|
| 309 | tmp_kspace(:,nv_full,:,:,:) = single(0); |
---|
| 310 | end |
---|
| 311 | tmp_data = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2); |
---|
[152] | 312 | |
---|
| 313 | if Dat.EPIPhasedArrayData |
---|
| 314 | data_block = abs(tmp_data); |
---|
| 315 | % Permute to correct orientation |
---|
| 316 | if strcmpi(orient,'trans90') |
---|
[153] | 317 | data_block = permute(data_block,[2 1 3 5 4]); |
---|
| 318 | data_block = flipdim(data_block,2); |
---|
| 319 | else |
---|
| 320 | data_block = permute(data_block,[1 2 3 5 4]); |
---|
[152] | 321 | end |
---|
[158] | 322 | data(:,:,:,vol_inds-nRef,:) = data_block; |
---|
[152] | 323 | else |
---|
[162] | 324 | if nRcvrs == 1 |
---|
| 325 | data_block = abs(tmp_data); |
---|
| 326 | else |
---|
| 327 | data_block = sqrt(sum(tmp_data.*conj(tmp_data),4)); |
---|
| 328 | end |
---|
[152] | 329 | % Permute to correct orientation |
---|
| 330 | if strcmpi(orient,'trans90') |
---|
[153] | 331 | data_block = permute(data_block,[2 1 3 5 4]); |
---|
| 332 | data_block = flipdim(data_block,2); |
---|
[152] | 333 | end |
---|
[170] | 334 | data(:,:,:,vol_inds-nRef-nDummys) = squeeze(data_block); |
---|
[152] | 335 | end |
---|
| 336 | |
---|
| 337 | end |
---|
| 338 | fprintf(1,'\n') |
---|
| 339 | end |
---|
| 340 | |
---|
| 341 | |
---|
| 342 | |
---|
| 343 | |
---|
| 344 | |
---|
| 345 | |
---|
| 346 | |
---|
| 347 | |
---|
| 348 | |
---|
| 349 | |
---|
| 350 | |
---|
| 351 | |
---|