source: vnmr_recon/epi_recon.m @ 170

Last change on this file since 170 was 170, checked in by tjniskan, 8 years ago
  • Epi reconstruction code now consideres all images in the data prior to image=0 as dummy scans

M aedes_revision.m
M vnmr_recon/epi_recon.m

File size: 6.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        % Number of reference images
71        if isfield(procpar,'epiref_type') && ~isempty(procpar.epiref_type) && ...
72                        strcmpi(procpar.epiref_type,'triple')
73                nRef = 3;
74        else
75                nRef = 1;
76        end
77       
78        % Number of segments
79        if isfield(procpar,'nseg') && ~isempty(procpar.nseg)
80                nSeg = procpar.nseg;
81        else
82                nSeg = 1;
83        end
84        %procpar.image = [1 1 0 -2 -1 1 1 1 1];
85        % Get index to first reference image
86        ref_start = find(procpar.image == 0);
87        ref_start = ref_start(1);
88        nDummys = 0;
89        if isempty(ref_start)
90                kspace = [];
91                data = [];
92                msg_out = 'Reference image(s) not found.';
93                return
94        elseif ref_start(1) > 1
95                nDummys = ref_start-1;
96                warning('Rejecting %d volumes before first reference image as dummy scans...',nDummys);
97        end
98       
99        % Number of volumes
100        nVols = length(procpar.image(nDummys+1:end))-nRef;
101       
102        % Number of slices
103        ns = procpar.ns;
104       
105        if isEPIP
106                nv = procpar.nphase; % Phase sampling
107                np = procpar.nread/2; % Read sampling
108        else
109                nv = procpar.nv;
110                np = procpar.np/2;
111        end
112       
113        % Get orientation
114        orient = procpar.orient{1};
115       
116        if isEPIP && (nv==64 & np==64)
117                % The information about navigator echoes in the procpar doesn't seem to
118                % be consistent. This code tried to guess the number of navigators from
119                % the k-space size
120                nNav = size(kspace,1)/np-nv;
121        end
122       
123        % Calculate indexes for navigator and data echoes
124        if nNav > 0
125                ind = reshape([1:(nv+nNav*nSeg)],[],nSeg);
126                NavInd = ind(1:nNav,:);
127                DataInd = ind(nNav+1:end,:);
128                NavInd = NavInd(:).';
129                DataInd = DataInd(:).';
130        else
131                DataInd = [1:nv];
132                NavInd = [];
133        end
134       
135        if Dat.EPIPhasedArrayData
136                data = zeros(nv,np,ns,nVols,nRcvrs,'single');
137        else
138                data = zeros(nv,np,ns,nVols,'single');
139        end
140       
141        % Look for phase table
142        pe_table = [];
143        if isfield(Dat,'phasetable') && ~isempty(Dat.phasetable)
144                pe_table = Dat.phasetable;
145        end
146       
147        % Reshape kspace and detach reference images
148        kspace = reshape(kspace,np,nv+nNav*nSeg,ns,nRcvrs,nVols+nRef+nDummys);
149       
150        ref1_ind = find(procpar.image==0,1);
151        ref2_ind = find(procpar.image==-2,1);
152        im1_ind = find(procpar.image==-1,1);
153       
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;
190                end
191        end
192       
193        % Reconstruct in blocks
194        kssz=size(kspace);
195        blksz = Dat.EPIBlockSize; % Process EPI data in 100 volume blocks (default)
196        nBlocks = ceil(nVols/blksz);
197        lnum = length(num2str(nBlocks));
198        lnumstr = num2str(lnum);
199        bsl = lnum*2+1;
200        fprintf(1,'Processing data in blocks of %d volumes\n',blksz)
201        fprintf(1,['Processing block...%0',lnumstr,'d/%0',lnumstr,'d'],1,nBlocks);
202       
203        for ii=1:nBlocks
204                fprintf(1,repmat('\b',1,bsl));
205                fprintf(1,['%0',lnumstr,'d/%0',lnumstr,'d'],ii,nBlocks);
206                tmp_data = [];
207                if ii==nBlocks
208                        vol_inds = (ii-1)*blksz+1+nRef+nDummys:nVols+nRef+nDummys;
209                else
210                        vol_inds = (ii-1)*blksz+1+nRef+nDummys:ii*blksz+nRef+nDummys;
211                end
212               
213                % At the moment navigator echoes are not used...
214                %tmp_kspace = kspace(:,nNav+1:end,:,:,vol_inds);
215                tmp_kspace = kspace(:,DataInd,:,:,vol_inds);
216                if ~isempty(pe_table)
217                        tmp_kspace(:,pe_table,:,:,:) = tmp_kspace;
218                end
219               
220                % Flip even lines
221                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,:));
227                                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               
236               
237                if Dat.EPIPhasedArrayData
238                        data_block = abs(tmp_data);
239                        % Permute to correct orientation
240                        if strcmpi(orient,'trans90')
241                                data_block = permute(data_block,[2 1 3 5 4]);
242                                data_block = flipdim(data_block,2);
243                        else
244                                data_block = permute(data_block,[1 2 3 5 4]);
245                        end
246                        data(:,:,:,vol_inds-nRef,:) = data_block;
247                else
248                        if nRcvrs == 1
249                                data_block = abs(tmp_data);
250                        else
251                                data_block = sqrt(sum(tmp_data.*conj(tmp_data),4));
252                        end
253                        % Permute to correct orientation
254                        if strcmpi(orient,'trans90')
255                                data_block = permute(data_block,[2 1 3 5 4]);
256                                data_block = flipdim(data_block,2);
257                        end
258                        data(:,:,:,vol_inds-nRef-nDummys) = squeeze(data_block);
259                end
260               
261        end
262        fprintf(1,'\n')
263end
264       
265
266
267
268
269
270
271
272
273
274
275
Note: See TracBrowser for help on using the repository browser.

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