source: vnmr_recon/epi_recon.m @ 157

Last change on this file since 157 was 157, checked in by tjniskan, 8 years ago
  • Tried to fix some navigator echo issues in reconstructing EPIP data

M aedes_revision.m
M vnmr_recon/epi_recon.m

File size: 5.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'};
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 strncmpi(procpar.seqfil,'epip',4) && 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 volumes
79        nVols = length(procpar.image)-nRef;
80       
81        % Number of slices
82        ns = procpar.ns;
83       
84       
85        if isEPIP
86                nv = procpar.nphase; % Phase sampling
87                np = procpar.nread/2; % Read sampling
88        else
89                nv = procpar.nv;
90                np = procpar.np/2;
91        end
92       
93        % Get orientation
94        orient = procpar.orient{1};
95       
96        % Remove navigators
97        if isEPIP && (nv==64 & np==64)
98                % The information about navigator echoes in the procpar doesn't seem to
99                % be consistent. This code tried to guess the number of navigators from
100                % the k-space size
101                nNav = size(kspace,1)/np-nv;
102                if nNav~=0
103                        nav_data = kspace(1:(np*nNav),:);
104                else
105                        nav_data = [];
106                end
107        else
108                if nNav~=0
109                        nav_data = kspace(1:(np*nNav),:);
110                        %kspace(1:(np*nNav),:)=[];
111                else
112                        nav_data = [];
113                end
114        end
115       
116        if Dat.EPIPhasedArrayData
117                data = zeros(nv,np,ns,nVols,nRcvrs,'single');
118        else
119                data = zeros(nv,np,ns,nVols,'single');
120        end
121       
122        % Reshape kspace and detach reference images
123        kspace = reshape(kspace,np,nv+nNav,ns,nRcvrs,nVols+nRef);
124       
125        % Flip even kspace lines
126%       if ~Dat.FastDataRead
127%               for tt=2+nNav:2:size(kspace,2)
128%                       kspace(:,tt,:,:,:) = flipdim(kspace(:,tt,:,:,:),1);
129%               end
130%       else
131%               for tt=2+nNav:2:size(kspace,2)
132%                       kspace(:,tt,:,:,:) = flipdim(kspace(:,tt,:,:,:),1);
133%               end
134%               %kspace(:,2+nNav:2:end,:,:,:) = flipdim(kspace(:,2+nNav:2:end,:,:,:),1);
135%       end
136        ref_data = kspace(:,nNav+1:end,:,:,1:nRef);
137        ref_data(:,2:2:end,:,:,:) = flipdim(ref_data(:,2:2:end,:,:,:),1);
138        %kspace(:,:,:,:,1:nRef)=[];
139       
140        % Calculate phase corrections for all receivers
141        phase_e = ref_data(:,:,:,:,1);
142        rev_phase = ref_data(:,:,:,:,1);
143        for ii=1:nRcvrs
144                if nRef==1
145      ref_im = ref_data(:,:,:,ii,1);
146      phase1 = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
147                        rev_phase = 1;
148                        phase_e(:,:,:,ii,:) = phase1;
149                elseif nRef==3
150                        ref1_ind = find(procpar.image==0,1);
151                        ref1 = ref_data(:,:,:,ii,ref1_ind);
152                        phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
153                       
154                        ref2_ind = find(procpar.image==-2,1);
155                        ref2 = flipdim(ref_data(:,:,:,ii,ref2_ind),1);
156                        phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
157                       
158                        im1_ind = find(procpar.image==-1,1);
159                        im1 = flipdim(ref_data(:,:,:,ii,im1_ind),1);
160                       
161                        rev_phase(:,:,:,ii,:) = fft(im1,[],1).*phase2;
162                        phase_e(:,:,:,ii,:) = phase1;
163                end
164        end
165       
166        % Reconstruct in blocks
167        kssz=size(kspace);
168        blksz = Dat.EPIBlockSize; % Process EPI data in 100 volume blocks (default)
169        nBlocks = ceil(nVols/blksz);
170        lnum = length(num2str(nBlocks));
171        lnumstr = num2str(lnum);
172        bsl = lnum*2+1;
173        fprintf(1,'Processing data in blocks of %d volumes\n',blksz)
174        fprintf(1,['Processing block...%0',lnumstr,'d/%0',lnumstr,'d'],1,nBlocks);
175       
176        for ii=1:nBlocks
177                fprintf(1,repmat('\b',1,bsl));
178                fprintf(1,['%0',lnumstr,'d/%0',lnumstr,'d'],ii,nBlocks);
179                tmp_data = [];
180                if ii==nBlocks
181                        vol_inds = (ii-1)*blksz+1+nRef:nVols+nRef;
182                else
183                        vol_inds = (ii-1)*blksz+1+nRef:ii*blksz+nRef;
184                end
185                tmp_kspace = kspace(:,nNav+1:end,:,:,vol_inds);
186               
187                % Flip even lines
188                tmp_kspace(:,2:2:end,:,:,:) = flipdim(tmp_kspace(:,2:2:end,:,:,:),1);
189
190                for kk=1:nRcvrs
191                        if nRef==3
192                                for tt=1:length(vol_inds)
193                                        tmp_kspace(:,:,:,kk,tt) = ifft(rev_phase(:,:,:,kk,:)+fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:));
194                                end
195                        else
196                                for tt=1:length(vol_inds)
197                                        tmp_kspace(:,:,:,kk,tt) = ifft(fft(tmp_kspace(:,:,:,kk,tt),[],1).*phase_e(:,:,:,kk,:));
198                                end
199                        end
200                        tmp_data = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
201                end
202               
203               
204                if Dat.EPIPhasedArrayData
205                        data_block = abs(tmp_data);
206                        % Permute to correct orientation
207                        if strcmpi(orient,'trans90')
208                                data_block = permute(data_block,[2 1 3 5 4]);
209                                data_block = flipdim(data_block,2);
210                        else
211                                data_block = permute(data_block,[1 2 3 5 4]);
212                        end
213                        data(:,:,:,:,vol_inds-nRef) = data_block;
214                else
215                        data_block = sqrt(sum(tmp_data.*conj(tmp_data),4));
216                        % Permute to correct orientation
217                        if strcmpi(orient,'trans90')
218                                data_block = permute(data_block,[2 1 3 5 4]);
219                                data_block = flipdim(data_block,2);
220                        end
221                        data(:,:,:,vol_inds-nRef) = squeeze(data_block);
222                end
223               
224        end
225        fprintf(1,'\n')
226end
227       
228
229
230
231
232
233
234
235
236
237
238
Note: See TracBrowser for help on using the repository browser.

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