source: vnmr_recon/epi_recon.m @ 152

Last change on this file since 152 was 152, checked in by tjniskan, 8 years ago
  • aedes_readvnmr.m almost works with epi data...
  • The ROI averages of fMRI time series are now shown in one column instead of two in plugins/fmri_plugins/basic_fmri_analysis.m

M aedes_readvnmr.m
M plugins/fmri_plugins/basic_fmri_analysis.m
M aedes_revision.m
M vnmr_recon/epi_recon.m

File size: 5.1 KB
RevLine 
[128]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
[152]8  kspace = {'epi_se_rapid_sp3','epi','epi_fMRI','epip'};
[128]9  return
10end
11
12data=[];
13msg_out = '';
14
[152]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% =====================================================
[128]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]);
[152]44       
45       
46        % ===================================
47        % Reconstruct VNMRj and EPIP sequences
48        % ===================================
[128]49else
[152]50       
51        if strncmpi(procpar.seqfil,'epip',4) && isfield(procpar,'nphase') && ...
52                        isfield(procpar,'nread')
53                isEPIP = true;
54        else
55                isEPIP = false;
56        end
57       
58        % Number of navigator echos
59        if isfield(procpar,'navigator') && strcmpi(procpar.navigator,'y')
60                nNav = length(procpar.nav_echo);
61        else
62                nNav = 0;
63        end
64       
65        % Number of reference images
66        if isfield(procpar,'epiref_type') && ~isempty(procpar.epiref_type) && ...
67                        strcmpi(procpar.epiref_type,'triple')
68                nRef = 3;
69        else
70                nRef = 1;
71        end
72       
73        % Number of volumes
74        nVols = length(procpar.image)-nRef;
75       
76        % Number of slices
77        ns = procpar.ns;
78       
79        % Phase sampling
80        if isEPIP
81                nv = procpar.nphase;
82               
83                % Read sampling
84                np = procpar.nread/2;
85        else
86                nv = procpar.nv;
87                np = procpar.np/2;
88        end
89       
90        % Get orientation
91        orient = procpar.orient{1}
92       
93        % Remove navigators
94        if nNav~=0
95                nav_data = kspace(1:(np*nNav),:);
96                kspace(1:(np*nNav),:)=[];
97        else
98                nav_data = [];
99        end
100       
101        if Dat.EPIPhasedArrayData
102                data = zeros(nv,np,ns,nVols,nRcv,'single');
103        else
104                data = zeros(nv,np,ns,nVols,'single');
105        end
106       
107        % Reshape kspace and detach reference images
108        if isEPIP
109                kspace = reshape(kspace,np,nv,nRcvrs,ns,nVols+nRef);
110        else
111                kspace = reshape(kspace,np,nv,nRcvrs,ns,nVols+nRef);
112        end
113       
114        % Flip even kspace lines
115        if ~Dat.FastDataRead
116                for tt=2:2:size(kspace,2)
117                        kspace(:,tt,:,:,:) = flipdim(kspace(:,tt,:,:,:),1);
118                end
119        else
120                kspace(:,2:2:end,:,:,:) = flipdim(kspace(:,2:2:end,:,:,:),1);
121        end
122        ref_data = kspace(:,:,:,:,1:nRef);
123        kspace(:,:,:,:,1:nRef)=[];
124       
125        % Calculate phase corrections for all receivers
126        phase_e = ref_data(:,:,:,:,1);
127        rev_phase = ref_data(:,:,:,:,1);
128        for ii=1:nRcvrs
129                if nRef==1
130      ref_im = ref_data(:,:,ii,:,1);
131      phase1 = exp(-sqrt(-1)*angle(fft(ref_im,[],1)));
132                        rev_phase = 1;
133                        phase_e(:,:,ii,:,:) = phase1;
134                elseif nRef==3
135                        ref1_ind = find(procpar.image==0,1);
136                        ref1 = ref_data(:,:,ii,:,ref1_ind);
137                        phase1 = exp(-sqrt(-1)*angle(fft(ref1,[],1)));
138                       
139                        ref2_ind = find(procpar.image==-2,1);
140                        ref2 = flipdim(ref_data(:,:,ii,:,ref2_ind),1);
141                        phase2 = exp(-sqrt(-1)*angle(fft(ref2,[],1)));
142                       
143                        im1_ind = find(procpar.image==-1,1);
144                        im1 = flipdim(ref_data(:,:,ii,:,im1_ind),1);
145                       
146                        rev_phase(:,:,ii,:,:) = fft(im1,[],1).*phase2;
147                        phase_e(:,:,ii,:,:) = phase1;
148                end
149        end
150       
151        % Reconstruct in blocks
152        kssz=size(kspace);
153        blksz = Dat.EPIBlockSize; % Process EPI data in 100 volume blocks (default)
154        nBlocks = ceil(nVols/blksz);
155        lnum = length(num2str(nBlocks));
156        lnumstr = num2str(lnum);
157        bsl = lnum*2+1;
158        fprintf(1,'Processing data in blocks of %d volumes\n',blksz)
159        fprintf(1,['Processing block...%0',lnumstr,'d/%0',lnumstr,'d'],1,nBlocks);
160       
161        for ii=1:nBlocks
162                fprintf(1,repmat('\b',1,bsl));
163                fprintf(1,['%0',lnumstr,'d/%0',lnumstr,'d'],ii,nBlocks);
164                tmp_data = [];
165                if ii==nBlocks
166                        vol_inds = (ii-1)*blksz+1:nVols;
167                else
168                        vol_inds = (ii-1)*blksz+1:ii*blksz;
169                end
170                tmp_kspace = kspace(:,:,:,:,vol_inds);
171               
172                for kk=1:nRcvrs
173                        if nRef==3
174                                for tt=1:length(vol_inds)
175                                        tmp_kspace(:,:,kk,:,tt) = ifft(rev_phase(:,:,kk,:,:)+fft(kspace(:,:,kk,:,vol_inds(tt)),[],1).*phase_e(:,:,kk,:,:));
176                                end
177                        else
178                                for tt=1:length(vol_inds)
179                                        tmp_kspace(:,:,kk,:,tt) = ifft(fft(kspace(:,:,kk,:,vol_inds(tt)),[],1).*phase_e(:,:,kk,:,:));
180                                end
181                        end
182                        tmp_data = fftshift(fftshift(fft(fft(tmp_kspace,[],1),[],2),1),2);
183                end
184               
185               
186                if Dat.EPIPhasedArrayData
187                        data_block = abs(tmp_data);
188                        data_block = permute(data_block,[1 2 4 5 3]);
189                        % Permute to correct orientation
190                        if strcmpi(orient,'trans90')
191                                data_block = permute(data_block,[2 1 3 4 5]);
192                        end
193                        data(:,:,:,vol_inds,:) = data_block;
194                else
195                        data_block = sqrt(sum(tmp_data.*conj(tmp_data),3));
196                        % Permute to correct orientation
197                        if strcmpi(orient,'trans90')
198                                data_block = permute(data_block,[2 1 3 4 5]);
199                        end
200                        data(:,:,:,vol_inds) = squeeze(data_block);
201                end
202               
203        end
204        fprintf(1,'\n')
205       
206end
207       
208
209
210
211
212
213
214
215
216
217
218
Note: See TracBrowser for help on using the repository browser.

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