source: vnmr_recon/epi_recon.m @ 177

Last change on this file since 177 was 177, checked in by tjniskan, 8 years ago
  • 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 size: 8.9 KB
Line 
1function [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
7if nargin==0
8  kspace = {'epi_se_rapid_sp3','epi','epi_fMRI','epip','epip_fMRI','epiT1rho'};
9  return
10end
11
12data=[];
13msg_out = '';
14
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% =====================================================
27if 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]);
44       
45       
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       
51        % ===================================
52        % Reconstruct VNMRj and EPIP sequences
53        % ===================================
54else
55       
56        if isfield(procpar,'nphase') && ...
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       
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')
103                nRef = 3;
104        elseif strcmpi(Dat.EPI_PC,'pointwise')
105                nRef = 1;
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
115       
116        % Number of segments
117        if isfield(procpar,'nseg') && ~isempty(procpar.nseg)
118                nSeg = procpar.nseg;
119        else
120                nSeg = 1;
121        end
122               
123       
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;
129        if isempty(ref_start) && doPhaseCorrection
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
138       
139        % Number of volumes
140        nVols = length(procpar.image(nDummys+1:end))-nRef;
141       
142        % Number of slices
143        ns = procpar.ns;
144       
145        if isEPIP
146                nv = procpar.nphase; % Phase sampling
147                np = procpar.nread/2; % Read sampling
148        else
149                nv = procpar.nv;
150                np = procpar.np/2;
151        end
152       
153        % Get orientation
154        orient = procpar.orient{1};
155       
156        if isEPIP && (nv==64 & np==64)
157                % The information about navigator echoes in the procpar doesn't seem to
158                % be consistent. This code tried to guess the number of navigators from
159                % the k-space size
160                nNav = size(kspace,1)/np-nv;
161        end
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       
171        % Calculate indexes for navigator and data echoes
172        if nNav > 0
173                ind = reshape([1:(nv+nNav*nSeg)],[],nSeg);
174                NavInd = ind(1:nNav,:);
175                DataInd = ind(nNav+1:end,:);
176                NavInd = NavInd(:).';
177                DataInd = DataInd(:).';
178        else
179                DataInd = [1:nv];
180                NavInd = [];
181        end
182       
183        if Dat.EPIPhasedArrayData
184                data = zeros(nv_full,np,ns,nVols,nRcvrs,'single');
185        else
186                data = zeros(nv_full,np,ns,nVols,'single');
187        end
188       
189        % Look for phase table
190        pe_table = [];
191        if isfield(Dat,'phasetable') && ~isempty(Dat.phasetable)
192                pe_table = Dat.phasetable;
193        end
194       
195        % Reshape kspace and detach reference images
196        kspace = reshape(kspace,np,nv+nNav*nSeg,ns,nRcvrs,nVols+nRef+nDummys);
197       
198       
199        ref1_ind = find(procpar.image==0,1);
200        ref2_ind = find(procpar.image==-2,1);
201        im1_ind = find(procpar.image==-1,1);
202       
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
255                end
256        end
257       
258        % Reconstruct in blocks
259        kssz=size(kspace);
260        blksz = Dat.EPIBlockSize; % Process EPI data in 100 volume blocks (default)
261        nBlocks = ceil(nVols/blksz);
262        lnum = length(num2str(nBlocks));
263        lnumstr = num2str(lnum);
264        bsl = lnum*2+1;
265        fprintf(1,'Processing data in blocks of %d volumes\n',blksz)
266        fprintf(1,['Processing block...%0',lnumstr,'d/%0',lnumstr,'d'],1,nBlocks);
267       
268        for ii=1:nBlocks
269                fprintf(1,repmat('\b',1,bsl));
270                fprintf(1,['%0',lnumstr,'d/%0',lnumstr,'d'],ii,nBlocks);
271                tmp_data = [];
272                if ii==nBlocks
273                        vol_inds = (ii-1)*blksz+1+nRef+nDummys:nVols+nRef+nDummys;
274                else
275                        vol_inds = (ii-1)*blksz+1+nRef+nDummys:ii*blksz+nRef+nDummys;
276                end
277               
278                % At the moment navigator echoes are not used...
279                %tmp_kspace = kspace(:,nNav+1:end,:,:,vol_inds);
280                tmp_kspace = kspace(:,DataInd,:,:,vol_inds);
281                if ~isempty(pe_table)
282                        tmp_kspace(:,pe_table,:,:,:) = tmp_kspace;
283                end
284               
285                % Flip even lines
286                tmp_kspace(:,2:2:end,:,:,:) = flipdim(tmp_kspace(:,2:2:end,:,:,:),1);
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
298                                end
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);
306               
307                if Dat.EPIPhasedArrayData
308                        data_block = abs(tmp_data);
309                        % Permute to correct orientation
310                        if strcmpi(orient,'trans90')
311                                data_block = permute(data_block,[2 1 3 5 4]);
312                                data_block = flipdim(data_block,2);
313                        else
314                                data_block = permute(data_block,[1 2 3 5 4]);
315                        end
316                        data(:,:,:,vol_inds-nRef,:) = data_block;
317                else
318                        if nRcvrs == 1
319                                data_block = abs(tmp_data);
320                        else
321                                data_block = sqrt(sum(tmp_data.*conj(tmp_data),4));
322                        end
323                        % Permute to correct orientation
324                        if strcmpi(orient,'trans90')
325                                data_block = permute(data_block,[2 1 3 5 4]);
326                                data_block = flipdim(data_block,2);
327                        end
328                        data(:,:,:,vol_inds-nRef-nDummys) = squeeze(data_block);
329                end
330               
331        end
332        fprintf(1,'\n')
333end
334       
335
336
337
338
339
340
341
342
343
344
345
Note: See TracBrowser for help on using the repository browser.

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