1function [labels, aux] = autoLabel(seq, varargin)
2% Automatically find values of the entire sequence or its part.
3%
4% autoLabel(seqObj) automatically detects the label evolution from the
5% k-space analysis and applies it to the sequence. Return value of the
6% function is the structure 'labels' with fields named after the labels
7% used in the sequence containing the evolution of the labes for alls
8% ADCs. This structure is similar to that returned by evalLabels().
9%
10% autoLabel(...,'blockRange',[first last]) Evaluate label
11% values starting from the first specified block to the last
12% one.
13%
14% autoLabel(...,'reflect',[1,2]) Reflect k-space trajectories along
15% directions 1 and 2 (any number of directions between 0 and 3 can be
16% specified) prior to generating label ranges. Affects both Fourier
17% encoding dimensions (redout, phase and partition encoding) and the
18% slice ordering. If used in combination with 'reorder', it is performed
19% first.
20%
21% autoLabel(...,'reorder',[2,1]) Reorder axes of k-space trajectories
22% (in this example by swapping directions 1 and 2 prior to generating
23% label ranges (reordering vector can contain 2 or 3 entries). If used in
24% combination with 'reflect', it is performed last.
25%
26% autoLabel(...,'mirrorFourier',true) Mirrors all Fourier encoding
27% directions, e.g. if ifft() is used instead of fft() by the image
28% reconstruction algorithm. It essentially reverts all Fourier encoding
29% dimensions (readout, phase, partition encoding) but does not affect
30% slice ordering. Can be arbitrarily combined with 'reflect'.
31%
32% autoLabel(...,'useLabels',labels_struct) Skip label evolution
33% detection and only apply externally-calculated labels to the sequence.
34%
35% autoLabel(...,'useAux',aux_struct) use externally-calculated aux struct
36% to apply definitions the sequence. Une combination with 'useLabels'.
37%
38% autoLabel(...,'skipApply',true) Skip applying label evolution to the
39% sequence and onoly return the labaels struct.
40%
41% From the experience, Siemens scanners need 'mirrorFourier' if the "XYZ
42% in TRA" coordinate mapping (default) mode is applied, which transforms
43% the sequence by the following rotation matrix [0 -1 0; -1 0 0; 0 0 -1]
44% prior to passing it to IDEA; On Siemens 'sortSlices'='descending' is
45% optional, otherwise the interpreter will change the slice indexes.
46%
47% see evalLabels()
49validSliceSorting={'acquisition','ascending','descending'};
50persistent parser
51if isempty(parser)
52 parser = inputParser;
53 parser.FunctionName = 'autoLabel';
54 parser.addParamValue('blockRange',[1 inf],@(x)(isnumeric(x) && length(x)==2));
55 parser.addParamValue('useLabels',struct([]),@(x)(isempty(x) || isstruct(x)));
56 parser.addParamValue('useAux',struct([]),@(x)(isempty(x) || isstruct(x)));
57 parser.addParamValue('skipApply',false,@(x)(islogical(x) && isscalar(x)));
58 parser.addParamValue('mirrorFourier', false, @(x)(islogical(x) && isscalar(x)));
59 parser.addParamValue('reflect', [], @(x)(isnumeric(x) && numel(x)<4));
60 parser.addParamValue('reorder', [], @(x)(isnumeric(x) && numel(x)<4));
61 parser.addParamValue('sortSlices','acquisition',@(x) any(validatestring(x,validSliceSorting)));
62 parser.addParamValue('noPlots', false, @(x)(islogical(x) && isscalar(x)));
63end
64parse(parser,varargin{:});
65opt = parser.Results;
67if ~isempty(opt.useLabels) && (~isempty(opt.reflect) || ~isempty(opt.reorder))
68 error('Optional parameters ''reflect'' or ''reorder'' only effective for the detection part and cannot be used together with ''useLabels''');
69end
71if ~isempty(opt.reflect) && numel(opt.reflect)~=numel(unique(opt.reflect))
72 error('All indices in ''reflect'' must be unique');
73end
75if ~isempty(opt.reorder)
76 if numel(opt.reorder)~=numel(unique(opt.reorder))
77 error('All indices in ''reorder'' must be unique');
78 end
79 if numel(opt.reorder)~=3 && min(opt.reorder)==1 && max(opt.reorder)==2
80 error('If ''reorder'' contains two indices they must be [1 2] or [2 1]');
81 end
82end
84if ~isfinite(opt.blockRange(2))
85 opt.blockRange(2)=length(seq.blockEvents);
86end
88aux = struct();
90%% index ADCs, part 1
91blockStartTimes=cumsum([0 seq.blockDurations]);
92blockStartTimes=blockStartTimes(opt.blockRange(1):opt.blockRange(2));
93adcLengths=[];
94b_adc=[];
95for iB=opt.blockRange(1):opt.blockRange(2)
96 raw_block = seq.getRawBlockContentIDs(iB); % this is much faster than seq.getBlock(iB)
97 if ~isempty(raw_block.adc)
98 libData = seq.adcLibrary.data(raw_block.adc).array;
99 % from getBlock(): adc = cell2struct(num2cell(libData(1:end-1)), {'numSamples', 'dwell', 'delay', ...
100 adcLengths(end+1)=libData(1); % libData(1) is adc.numSamples
101 b_adc(end+1)=iB-opt.blockRange(1)+1;
102 end
103 %block = seq.getBlock(iB);
104 %if ~isempty(block.adc)
105 % adcLengths(end+1)=block.adc.numSamples;
106 % b_adc(end+1)=iB-opt.blockRange(1)+1;
107 %end
108end
109adcStartCounts=cumsum([1, adcLengths(1:end-1)]);
110% b_adc=zeros(size(t_adcStarts));
111% for i=1:numel(t_adcStarts)
112% b_adc(i)=find(blockStartTimes<t_adcStarts(i),1,'last');
113% end
115if isempty(opt.useLabels)
117 %% preparing calculations
118 [ktraj_adc, t_adc, ~, ~, t_excitation, ~, slicepos, t_slicepos, gw_pp] = seq.calculateKspacePP('blockRange',opt.blockRange);
120 t_adcStarts=t_adc(adcStartCounts);
121 firstNonNoiseAdc=find(t_adcStarts>t_excitation(1),1,'first');
122 firstNonNoiseAdcSampe=find(t_adc>t_excitation(1),1,'first');
123 nReadouts = numel(t_adcStarts)-firstNonNoiseAdc+1;
125 sliceGrads=zeros(size(slicepos));
126 for i=1:3
127 sliceGrads(i,:)=ppval(gw_pp{i},t_slicepos);
128 end
130 if opt.mirrorFourier
131 ktraj_adc=-ktraj_adc;
132 end
133 if ~isempty(opt.reflect)
134 ktraj_adc(opt.reflect,:)=-ktraj_adc(opt.reflect,:);
135 slicepos(opt.reflect,:)=-slicepos(opt.reflect,:);
136 sliceGrads(opt.reflect,:)=-sliceGrads(opt.reflect,:);
137 end
138 if ~isempty(opt.reorder)
139 ktraj_adc(1:numel(opt.reorder),:)=ktraj_adc(opt.reorder,:);
140 slicepos(1:numel(opt.reorder),:)=slicepos(opt.reorder,:);
141 sliceGrads(1:numel(opt.reorder),:)=sliceGrads(opt.reorder,:);
142 end
144 %% slice positions
145 [~,i]=max(abs(sliceGrads));
146 mainSliceGradSigns=sign(sliceGrads(sub2ind(size(sliceGrads),i,1:size(sliceGrads,2))));
147 sliceNormals=normalize(sliceGrads,'norm').*mainSliceGradSigns(ones(1,3),:);
148 sliceOffsets=dot(slicepos,sliceNormals);
149 sliceOffsets(~isfinite(sliceOffsets))=0;
150 [uniqueSlicePositions, I, sliceCountersAcquisitionOrder] = unique(sliceOffsets,'stable');
151 t_usp = t_slicepos(I);
152 b_usp = zeros(size(t_usp));
153 for i=1:numel(t_usp)
154 b_usp(i) = find(blockStartTimes<t_usp(i),1,'last');
155 end
156 [sortedSlicePositions, ~, sliceCountersSortedPositions] = unique(sliceOffsets);
158 % aux.SliceThickness
159 % aux.SliceGap
160 % todo: more checks, e.g. if all slice thickeness and gaps are the same...
161 GAs1=vecnorm(sliceGrads(:,1));
162 if GAs1>0
163 bSlice1=find(blockStartTimes<t_slicepos(1),1,'last');
164 b=seq.getBlock(bSlice1);
165 bw=mr.calcRfBandwidth(b.rf);
166 aux.SliceThickness=bw/GAs1;
167 if length(sortedSlicePositions)>1
168 aux.SliceGap=diff(sortedSlicePositions(1:2))-aux.SliceThickness;
169 end
170 end
171 if strcmp(opt.sortSlices,'descending')
172 sortedSlicePositions=sortedSlicePositions(end:-1:1);
173 sliceCountersSortedPositions = max(sliceCountersSortedPositions) + 1 - sliceCountersSortedPositions;
174 end
176 % useful results:
177 % uniqueSlicePositions : speaks for itself, sorted in the order of occurence
178 % sliceCountersAcquisitionOrder : slice cointers in the acquisition order
179 % b_usp: block indices containg RF pulses corresponding to the first
180 % occurence of each of the unique slice positions
181 % sortedSlicePositions : slice positions in the ascending or descending deorder
182 % sliceCountersSorted : slice counters corresponding to the sorted slice order
184 %% index ADCs, part 2
185 if ~strcmp(opt.sortSlices,'acquisition')
186 uniqueSlicePositions=sortedSlicePositions;
187 end
188 sliceCountersAdc=zeros(1,nReadouts);
189 for i=firstNonNoiseAdc:numel(t_adcStarts)
190 j=find(t_slicepos<t_adcStarts(i),1,'last');
191 sliceCountersAdc(i-firstNonNoiseAdc+1)=find(uniqueSlicePositions==sliceOffsets(j),1,'first');
192 end
193 % useful results:
194 % b_adc : vector containing block indices with ADC objects
195 % adcLengths : numbers of samples of each of the ADC objects
196 % t_adcStarts : time point of each ADCs first sample
197 % firstNonNoiseAdc : index of the first ADC after the noise scan (1 means no noise scan)
198 % firstNonNoiseAdcSample : index of the first ADC sample after the noise scan
199 % sliceCountersAdc : slice counters per ADC
201 %% readout analysis
202 cEchoPos = zeros(1,nReadouts);
203 tEcho = zeros(1,nReadouts);
204 kEcho = zeros(3,nReadouts);
205 gradReadout = zeros(3,nReadouts);
206 signReadout = zeros(1,nReadouts);
207 [~,kspaceCenterSample]=min(vecnorm(ktraj_adc(:,firstNonNoiseAdcSampe:end)));
208 kspaceCenterPoint=ktraj_adc(:,kspaceCenterSample);
209 isCartesianReadout=true;
210 for i=1:nReadouts
211 c1=adcStartCounts(i+firstNonNoiseAdc-1);
212 c2=c1+adcLengths(i+firstNonNoiseAdc-1)-1;
213 % find echo pos for each readout
214 [~,cEchoPos(i)]=min(vecnorm(ktraj_adc(:,c1:c2)-kspaceCenterPoint)); % for echos positioned between samples there can be some jitter... to reduce jitter we compare not to 0 but to the smallest absolute...
215 kEcho(:,i)=ktraj_adc(:,c1+cEchoPos(i)-1);
216 t_adcThisEcho=t_adc(c1+cEchoPos(i)-1);
217 if vecnorm(kspaceCenterPoint)>eps % this is more or less copied from seq.testRepoert()
218 % the actual echo might be between k-space samples, try to interpolate it
219 i2check=[];
220 % check if adc kspace trajectory has elements left and right to index_echo
221 if cEchoPos(i) > 1
222 i2check=c1+cEchoPos(i)-2;
223 end
224 if c1+cEchoPos(i)-1 < c2
225 i2check(end+1)=c1+cEchoPos(i);
226 end
227 for a=1:numel(i2check)
228 v_i_to_0=-kEcho(:,i);
229 v_i_to_t=ktraj_adc(:,i2check(a))-kEcho(:,i);
230 % project v_i_to_0 to v_o_to_t
231 p_vit=v_i_to_0'*v_i_to_t/(vecnorm(v_i_to_t)^2);
232 if p_vit>0
233 % we have forund a bracket for the echo and the proportionality coefficient is p_vit
234 t_adcThisEcho=t_adcThisEcho*(1-p_vit) + t_adc(i2check(a))*p_vit;
235 break;
236 end
237 end
238 end
240 i_sliseposThisEcho=find(t_slicepos<t_adcThisEcho,1,'last');
241 tEcho(i)=t_adcThisEcho-t_slicepos(i_sliseposThisEcho);
242 for j=1:3
243 gradReadout(j,i) = ppval(gw_pp{j},t_adcThisEcho); % this call is relatively expensive for long sequences, we should go away from it
244 end
245 if opt.mirrorFourier
246 gradReadout(:,i)=-gradReadout(:,i);
247 end
248 if ~isempty(opt.reflect)
249 gradReadout(opt.reflect,i)=-gradReadout(opt.reflect,i);
250 end
251 if ~isempty(opt.reorder)
252 gradReadout(1:numel(opt.reorder),i)=gradReadout(opt.reorder,i);
253 end
254 % store the most central readout trajectory as a reference
255 if c1<=kspaceCenterSample && c2>=kspaceCenterSample
256 cCentralReadout=i; % counter corresponding to the central readout
257 kCentralReadout=ktraj_adc(:,c1:c2);
258 [~,cCentralReadoutCenter]=min(vecnorm(kCentralReadout));
259 end
260 % detect readout alternation like for EPI
261 [~,j]=max(abs(gradReadout(:,i)));
262 signReadout(i)=sign(gradReadout(j,i));
263 % detect if the readout is sufficiently the same to see if the trajectory is Cartesian-like
264 if i>1 && vecnorm(gradReadout(:,i)*signReadout(i)-gradReadout(:,1)*signReadout(1))>1e-4 % why such threshold?
265 isCartesianReadout=false;
266 end
267 % % dictionary of readouts
268 % kReadoutKey = {ktraj_adc(:,c1:c2) - ktraj_adc(:,c1+cEchoPos(i)-1)};
269 % if i==1
270 % readoutDict=dictionary(kReadout,i);
271 % else
272 % if ~readoutDict.isKey(kReadout)
273 % readoutDict(kReadout)=i;
274 % end
275 % end
276 end
278 % more readout analysis
279 kCentralReadoutDirection=normalize(gradReadout(:,cCentralReadout),'norm');
280 kCentralReadoutProjection=dot(kCentralReadout,kCentralReadoutDirection(:,ones(1,size(kCentralReadout,2))));
281 dkR=median(diff(kCentralReadoutProjection));
283 % detect non-constant sampling
284 if any(abs(diff(kCentralReadoutProjection,2)/dkR*numel(kCentralReadoutProjection))>0.1)
285 % see if we can extraxt trapezoid resampling parameters
286 bCentralRO=find(blockStartTimes<t_adcStarts(firstNonNoiseAdcSampe+cCentralReadout-1),1,'last');
287 b=seq.getBlock(bCentralRO);
288 if ~isempty(opt.reorder) && opt.reorder(1)~=1
289 warning('EPI readout gridding parameter detection code is incompatible with the reorder option');
290 else
291 if isfield(b,'gx') && ~isempty(b.gx) && strcmp(b.gx.type,'trap') && isfield(b,'adc') && ~isempty(b.adc)
292 aux.TrapezoidGriddingParameters=[b.gx.riseTime, b.gx.flatTime, b.gx.fallTime, b.adc.delay-b.gx.delay b.adc.numSamples*b.adc.dwell];
293 % TODO: figure out a better way to find TargetGriddedSamples... In Siemens sequences it is often adc.numsamples
294 % we need kReadout to be able to do it, so we will patch this value below (in the trajectory analysis part)
295 aux.TargetGriddedSamples=b.adc.numSamples;
296 end
297 end
298 % TODO: figure out a better way to find TargetGriddedSamples... In Siemens sequences it is often adc.numsamples
299 end
301 % bipolar gradients: detect both and find the reflection point
302 if any(diff(signReadout)~=0)
303 % find an inverted sign readout close to k-space center
304 adcCenters=adcStartCounts(firstNonNoiseAdc:end)+cEchoPos-1;
305 isInverted=signReadout==-signReadout(cCentralReadout);
306 [~,cInvertedReadout]=min(vecnorm(ktraj_adc(:,adcCenters(isInverted))));
307 % here cInvertedReadout is the counter only amongst the inverted readouts, we need to restore its meaning for the total readouts vector
308 countInverted=cumsum(isInverted);
309 cInvertedReadout=find(cInvertedReadout==countInverted,1);
310 assert(isInverted(cInvertedReadout));
311 %
312 c1=adcStartCounts(firstNonNoiseAdc+cInvertedReadout-1);
313 c2=c1+adcLengths(firstNonNoiseAdc+firstNonNoiseAdc-1)-1;
314 kInvertedReadoutProjection=dot(ktraj_adc(:,c1:c2),kCentralReadoutDirection(:,ones(1,size(kCentralReadout,2))));
315 if signReadout(cCentralReadout)==1
316 aux.kReadout = [kCentralReadoutProjection; kInvertedReadoutProjection]; % 1: positive readout, 2: negative readout
317 else
318 aux.kReadout = [kInvertedReadoutProjection; kCentralReadoutProjection]; % 1: positive readout, 2: negative readout
319 end
320 if ~opt.noPlots
321 figure; plot(aux.kReadout(1,:),'.-'); hold on; plot(aux.kReadout(2,:),'.-');plot(aux.kReadout(2,end:-1:1),'x'); title('bipolar readout k-space alignment');
322 legend('forward','reverse','rev.refl'); % if crosses overlap the line then REV is symmetric
323 xlabel('sample number (acq. order)');
324 ylabel('k-space, 1/m');
325 end
326 else
327 aux.kReadout = kCentralReadoutProjection;
328 end
330 % useful results:
331 % isCartesianReadout : is this trajectory sufficiently Cartesian-like
332 % gradReadout : readout gradients as 3D vectors
333 % signReadout : sign of the dominant readout component
334 % kCentralReadout : the most central readout
335 % cCentralReadoutCenter : center of the above
336 % cEchoPos : position of the echo within each ADC vector for each readout
337 % tEcho : echo time for each readout
338 % kEcho : k-space location of each echo
340 %% detect navigators
341 isNavigator=false(1,nReadouts);
342 isNavigatorCandidate = vecnorm(kEcho-kspaceCenterPoint) < 1e-4; % why this threshold?
343 orderedReadoutIndicator=diff(kEcho-kspaceCenterPoint,1,2);
344 i=find(max(abs(orderedReadoutIndicator'))>1e-4);
345 orderedReadoutIndicator=orderedReadoutIndicator(i,:); % discard empty dimensions
346 %TODO: finish this detection some day... e.g. if all(isNavigatorCandidate(1:3)) && abs(orderedReadoutIndicator(1)-orderedReadoutIndicator(2))<1e-4 && max(abs(diff(orderedReadoutIndicator(4:x))))<1e-4
347 if nReadouts>=16 && all(isNavigatorCandidate(1:3)) && abs(orderedReadoutIndicator(1)-orderedReadoutIndicator(2))<1e-4 && ...
348 max(abs(diff(orderedReadoutIndicator(4:16))))<1e-4
349 aux.epiWithThreeEchoNavigator=true;
350 % this is a real hack to close those singular k=0 lines in the middle of readouts
351 isNavigator = (isNavigatorCandidate + circshift(isNavigatorCandidate,1) + circshift(isNavigatorCandidate,-1))>1.5;
352 % TODO: set SEG (010 10101...) and AVG (001 00000...)
353 end
355 %% Accelerate calculation for Cartesian-like sequences by only considering the echo position
356 if isCartesianReadout
357 ktraj_adc = ktraj_adc(:, adcStartCounts(firstNonNoiseAdc:end)+cEchoPos-1);
358 else
359 error('autoLabel() only supports sufficiently Cartesian sequences');
360 end
362 %% Analyze the trajectory data
364 dkR=median(diff(kCentralReadoutProjection));
366 k_extent=max(abs(ktraj_adc-kspaceCenterPoint),[],2);
367 k_scale=max(k_extent);
368 k_threshold=abs(dkR/50);
370 % calculate/update gridding samples
371 if any(any(abs(diff(aux.kReadout,2,2))>1e-6)) % threshold?
372 % calculate samples based on the largest k-space step, so that we never add samples in between that we didn't measure
373 aux.TargetGriddedSamples=floor((max(aux.kReadout(1,:))-min(aux.kReadout(1,:)))/max(abs(diff(aux.kReadout(1,:)))))+1;
374 end
376 % detect unused dimensions and delete them
377 if any(k_extent<k_threshold)
378 ktraj_adc(k_extent<k_threshold,:)=[]; % delete rows
379 kspaceCenterPoint(k_extent<k_threshold)=[];
380 k_extent(k_extent<k_threshold)=[];
381 end
383 % detect dK, k-space reordering and repetitions (or slices, etc)
384 kt_sorted=sort(ktraj_adc-kspaceCenterPoint,2);
385 dk_all=kt_sorted(:,2:end)-kt_sorted(:,1:(end-1));
386 dk_all(dk_all<k_threshold)=NaN;
387 dk_min=min(dk_all,[],2);
388 dk_max=max(dk_all,[],2);
389 dk_all(dk_all-dk_min(:,ones(1,size(dk_all,2)))>k_threshold)=NaN;
390 dk_all_cnt=sum(isfinite(dk_all),2);
391 dk_all(~isfinite(dk_all))=0;
392 dk=sum(dk_all,2)./dk_all_cnt;
393 dk(~isfinite(dk))=0; % dk vector ready
394 [~,k0_ind]=min(sum(ktraj_adc.^2,1)); % k-space center
395 kindex=round((ktraj_adc-ktraj_adc(:,k0_ind*ones(1,size(ktraj_adc,2))))./dk(:,ones(1,size(ktraj_adc,2))));
396 kindex(~isfinite(kindex))=0;
397 kindex_min=min(kindex,[],2);
398 kindex_mat=kindex-kindex_min(:,ones(1,size(ktraj_adc,2)))+1;
399 kindex_end=max(kindex_mat,[],2);
400 sampler=zeros(length(uniqueSlicePositions),prod(kindex_end));
401 repeat=zeros(1,size(ktraj_adc,2));
402 for i=1:size(kindex_mat,2)
403 if (~isNavigator(i))
404 % ignore navigators
405 switch size(kindex_mat,1)
406 case 3
407 ind=sub2ind(kindex_end,kindex_mat(1,i),kindex_mat(2,i),kindex_mat(3,i));
408 case 2
409 ind=sub2ind(kindex_end,kindex_mat(1,i),kindex_mat(2,i));
410 otherwise
411 ind=kindex_mat(1,i); %sub2ind(kindex_end,kindex_mat(1,i));
412 end
413 r=sampler(sliceCountersAdc(i),ind);
414 repeat(i)=r;
415 sampler(sliceCountersAdc(i),ind)=r+1;
416 end
417 end
418 % if (max(repeat(:))>0)
419 % kindex=[kindex;(repeat+1)];
420 % kindex_mat=[kindex_mat;(repeat+1)];
421 % kindex_end=max(kindex_mat,[],2);
422 % end
423 %figure; plot(kindex(1,:),kindex(2,:),'.-');
424 %figure; plot(kindex_mat(1,:),'.-');
426 % see if some repetitions are actully echoes/contrasts
427 nRep=max(repeat)+1;
428 skipECO=false;
429 if nRep>1
430 TE=zeros(1, nRep);
431 for i=1:nRep
432 try
433 if isvector(kindex)
434 TE(i)=tEcho(kindex==0 & sliceCountersAdc==1 & repeat==(i-1)); % we assume all slices have the same TEs and only check slice 1
435 else
436 TE(i)=tEcho(all(kindex==0) & sliceCountersAdc==1 & repeat==(i-1)); % we assume all slices have the same TEs and only check slice 1
437 end
438 catch
439 warning('unclear sequence structure, skipping TE & ECO counter detection');
440 skipECO=true;
441 break;
442 end
443 end
444 if ~skipECO
445 [TE_sorted, TE_order]=sort(TE);
446 TE_cluster=cumsum([1, diff(TE_sorted)>10e-6]);
447 unique_TE=zeros(1,max(TE_cluster));
448 for i=1:length(unique_TE)
449 unique_TE(i)=mean(TE_sorted(TE_cluster==i));
450 TE(TE_order(TE_cluster==i))=unique_TE(i);
451 end
452 aux.TE=unique_TE; % maybe we should fill this also for single-TE sequences?
453 echo=zeros(1,size(ktraj_adc,2));
454 echo_rep=zeros(1,nRep);
455 for i=1:nRep
456 cecho=find(abs(unique_TE-TE(i))<=1e-6);
457 echo(repeat==(i-1))=cecho;
458 repeat(repeat==(i-1))=sum(echo_rep==cecho);
459 echo_rep(i)=cecho;
460 end
461 end
462 end
464 %% create labels struct
465 labels=struct();
466 % NOISE
467 % SLC
468 % REV
469 % LIN,
470 % PAR
471 % ECO
472 % REP
474 % TODO: SEG,ECO,REF,IMA,AVG,SET
476 nADCs=nReadouts+firstNonNoiseAdc-1;
477 if firstNonNoiseAdc>1
478 labels.NOISE=zeros(1,nADCs);
479 labels.NOISE(1:(firstNonNoiseAdc-1))=1;
480 end
481 if any(sliceCountersAdc~=1)
482 if firstNonNoiseAdc>1
483 labels.SLC=zeros(1,nADCs);
484 labels.SLC(firstNonNoiseAdc:end)=sliceCountersAdc-1;
485 else
486 labels.SLC=sliceCountersAdc-1;
487 end
488 end
489 if any(signReadout<0)
490 if firstNonNoiseAdc>1
491 labels.REV=zeros(1,nADCs);
492 labels.REV(firstNonNoiseAdc:end)=(signReadout<0);
493 else
494 labels.REV=(signReadout<0);
495 end
496 end
497 if any(kindex_mat(1,:)~=1)
498 if firstNonNoiseAdc>1
499 labels.LIN=zeros(1,nADCs);
500 labels.LIN(firstNonNoiseAdc:end)=kindex_mat(1,:)-1;
501 else
502 labels.LIN=kindex_mat(1,:)-1;
503 end
504 end
505 if size(kindex_mat,1)>1 && any(kindex_mat(2,:)~=1)
506 if firstNonNoiseAdc>1
507 labels.PAR=zeros(1,nADCs);
508 labels.PAR(firstNonNoiseAdc:end)=kindex_mat(2,:)-1;
509 else
510 labels.PAR=kindex_mat(2,:)-1;
511 end
512 end
513 if exist('echo','var') && max(echo(:))>1
514 if firstNonNoiseAdc>1
515 labels.ECO=zeros(1,nADCs);
516 labels.ECO(firstNonNoiseAdc:end)=echo(:).'-1;
517 else
518 labels.ECO=echo(:).'-1;
519 end
520 end
521 if max(repeat(:))>0
522 if firstNonNoiseAdc>1
523 labels.REP=zeros(1,nADCs);
524 labels.REP(firstNonNoiseAdc:end)=repeat(:).';
525 else
526 labels.REP=repeat(:).';
527 end
528 end
529 if any(isNavigator)
530 labels.NAV=isNavigator;
531 end
533 if isfield(labels,'LIN'), aux.kSpaceCenterLine = labels.LIN(cCentralReadout); end
534 if isfield(labels,'PAR'), aux.kSpaceCenterPartition = labels.PAR(cCentralReadout); end
535 aux.kSpaceCenterSample = cCentralReadoutCenter-1;
536 if length(uniqueSlicePositions)>1
537 aux.SlicePositions = uniqueSlicePositions;
538 end
539 % aux.kSpacePhaseEncodingLines
540 % aux.PhaseResolution
541 % aux.ReadoutOversamplingFactor
542 % aux.SliceGap
543 % aux.SlicePositions
544 % aux.SliceThickness
545 % aux.TargetGriddedSamples
546 % aux.TrapezoidGriddingParameterss
547 % aux.AccelerationFactorPE
548 % aux.AccelerationFactor3D
549 % aux.FirstFourierLine
550 % aux.FirstRefLine
551 % aux.FirstFourier3D
552 % aux.FirstRef3D
554else
555 labels=opt.useLabels;
556end
558if ~isempty(opt.useAux)
559 aux=opt.useAux;
560end
562%% apply labels
564if ~opt.skipApply
565 lblNames=fieldnames(labels);
566 warnBcp = warning ('off','mr:fixmePreviousRotationExtension');
567 for i=1:numel(b_adc)
568 blkLabels={};
570 for j=1:length(lblNames)
571 if i==1 && labels.(lblNames{j})(1)~=0 || i>1 && labels.(lblNames{j})(i)~=labels.(lblNames{j})(i-1)
572 blkLabels{end+1}=mr.makeLabel('SET',lblNames{j},labels.(lblNames{j})(i));
573 end
574 end
576 if ~isempty(blkLabels)
577 iB=b_adc(i);
578 b=seq.getBlock(iB);
579 e=mr.block2events(b);
580 seq.setBlock(iB,e{:},blkLabels{:});
581 end
582 end
583 warning(warnBcp);
585 % convert usable fields of aux to sequence definitions
586 fieldsToExport={'kSpaceCenterLine','kSpaceCenterPartition','kSpaceCenterSample','kSpacePhaseEncodingLines',...
587 'PhaseResolution','ReadoutOversamplingFactor','SliceGap','SlicePositions','SliceThickness',...
588 'TargetGriddedSamples','TrapezoidGriddingParameters','AccelerationFactorPE','AccelerationFactor3D',...
589 'FirstFourierLine','FirstRefLine','FirstFourier3D','FirstRef3D'};
590 for i=1:length(fieldsToExport)
591 if isfield(aux,fieldsToExport{i})
592 prevDef=seq.getDefinition(fieldsToExport{i});
593 if ~isempty(prevDef)
594 warning('Overwriting existing sequence definition %s = %s', fieldsToExport{i}, num2str(prevDef));
595 end
596 seq.setDefinition(fieldsToExport{i},aux.(fieldsToExport{i}));
597 end
598 end
599end
601%% test/plot label settings
602if opt.noPlots
603 return;
604end
606if ~opt.skipApply
607 lbls=seq.evalLabels('evolution','adc');
608else
609 lbls=labels;
610end
611lbl_names=fieldnames(lbls);
612figure; hold on;
613for n=1:length(lbl_names)
614 plot(lbls.(lbl_names{n}));
615end
616legend(lbl_names(:));
617title('evolution of labels/counters/flags');
618xlabel('adc number');