1function [ report ] = testReport( obj, varargin )
2%testReport Analyze the sequence and return a text report
3% Currently no parameters are required. In future versions it may be
4% possible to (de)select some tests. Optional parameter 'system' allows
5% to test the limits against the given MR system
6%
7% maxim.zaitsev@uniklinik-freiburg.de
9 function skey = makeSkey(iKey)
10 if iKey>=0
11 skey=['p' sprintf('%d',iKey)];
12 else
13 skey=['m' sprintf('%d',-iKey)];
14 end
15 end
17 function s = optionalOut(format, number)
18 if any(abs(number)>eps)
19 s=sprintf(format, number);
20 else
21 s='';
22 end
23 end
25persistent parser
26if isempty(parser)
27 parser = inputParser;
28 parser.FunctionName = 'testReport';
30 addParamValue(parser,'system',struct([]),@isstruct);
31end
32parse(parser,varargin{:});
33opt = parser.Results;
35% find the RF pulses and list flip angles
36flipAnglesDeg=[];
37for k=obj.rfLibrary.keys
38 libData=obj.rfLibrary.data(k).array;
39 if length(obj.rfLibrary.type)>=k
40 rf=obj.rfFromLibData(libData,obj.rfLibrary.type(k));
41 else
42 rf=obj.rfFromLibData(libData);
43 end
44 flipAnglesDeg=[flipAnglesDeg abs(sum(rf.signal(1:end-1).*(rf.t(2:end)-rf.t(1:end-1))))*360]; %we use rfex.t(1) in place of opt.system.rfRasterTime
45end
46flipAnglesDeg=unique(flipAnglesDeg);
48% calculate TE and TR
50[duration, numBlocks, eventCount]=obj.duration();
52[wnt.gw_data, wnt.tfp_excitation, wnt.tfp_refocusing, wnt.t_adc]=obj.waveforms_and_times();
53%[ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc] = obj.calculateKspace();
54%[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = obj.calculateKspacePP();
55[ktraj_adc, t_adc, ~, ~, t_excitation, ~] = obj.calculateKspacePP('externalWaveformsAndTimes',wnt);
57% remove all ADC events that come before the first RF event (noise scans or alike)
58if ~isempty(t_excitation)
59 ktraj_adc=ktraj_adc(:,t_adc > t_excitation(1));
60 t_adc = t_adc(t_adc > t_excitation(1));
61end
63% trajectory calculation will fail for spin-echoes if seq is loaded from a
64% file for the current file format revision (1.2.0) because we do not store
65% the use of the RF pulses. Read function has an option 'detectRFuse' which
66% may help...
68%
69kabs_adc=sum(ktraj_adc.^2,1).^0.5;
70[kabs_echo, index_echo]=min(kabs_adc);
71t_echo=t_adc(index_echo); % just a first estimate, see if we can improve it
72if kabs_echo>eps
73 i2check=[];
74 % check if adc kspace trajectory has elements left and right to index_echo
75 if index_echo > 1
76 i2check=[i2check (index_echo-1)];
77 end
78 if index_echo < length(kabs_adc)
79 i2check=[i2check (index_echo+1)];
80 end
81 for a=1:numel(i2check)
82 v_i_to_0=-ktraj_adc(:,index_echo);
83 v_i_to_t=ktraj_adc(:,i2check(a))-ktraj_adc(:,index_echo);
84 % project v_i_to_0 to v_o_to_t
85 p_vit=v_i_to_0'*v_i_to_t/(vecnorm(v_i_to_t)^2);
86 if p_vit>0
87 % we have forund a bracket for the echo and the proportionality
88 % coefficient is p_vit
89 t_echo=t_adc(index_echo)*(1-p_vit) + t_adc(i2check(a))*p_vit;
90 break;
91 end
92 end
93end
95if ~isempty(t_excitation)
96 t_ex_tmp=t_excitation(t_excitation<t_echo);
97 TE=t_echo-t_ex_tmp(end);
98 % TODO detect multiple TEs
99else
100 TE=NaN;
101end
103if (length(t_excitation)<2)
104 TR=duration; % best estimate for now
105else
106 t_ex_tmp1=t_excitation(t_excitation>t_echo);
107 if isempty(t_ex_tmp1)
108 TR=t_ex_tmp(end)-t_ex_tmp(end-1);
109 else
110 TR=t_ex_tmp1(1)-t_ex_tmp(end);
111 end
112 % TODO check frequency offset to detect multiple slices
113end
115% check sequence dimensionality and spatial resolution
116k_extent=max(abs(ktraj_adc),[],2);
117k_scale=max(k_extent);
118if (k_scale~=0)
119 k_bins=4e6; % this defines our ability to separate k-space samples.
120 % lower values give us imunity to rounding errors in k-space calculations
121 % current code below (2nd pass) however merges neighboring cells (+-1)
122 k_threshold=k_scale/k_bins;
124 % detect unused dimensions and delete them
125 if any(k_extent<k_threshold)
126 ktraj_adc(k_extent<k_threshold,:)=[]; % delete rows
127 k_extent(k_extent<k_threshold)=[];
128 end
130 % bin the k-space trajectory to detect repetitions / slices
131 k_len=size(ktraj_adc,2);
132 k_repeat=zeros(1,k_len);
133 k_storage=zeros(1,k_len);
134 k_storage_next=1;
135 % the fastest option should be dictionary(), if it is available
136 hasDict=false;
137 try
138 kmap=configureDictionary("string","double"); %is too new and only supported sice 2023; %containers.Map('KeyType','char','ValueType','double'); % works as well but is substantially slower
139 hasDict=true;
140 catch
141 end
142 if hasDict
143 for i=1:k_len
144 key = sprintf('a%d', int32(k_bins+round(ktraj_adc(:,i)/k_threshold)));
145 k_storage_ind=kmap.lookup(key,'FallbackValue',0);
146 if k_storage_ind==0
147 k_storage_ind=k_storage_next;
148 kmap(key)=k_storage_ind;
149 k_storage_next=k_storage_next+1;
150 end
151 k_storage(k_storage_ind)=k_storage(k_storage_ind)+1;
152 k_repeat(i) = k_storage(k_storage_ind);
153 end
154 else
155 % we use strings in combination with structs to go fast... (matlab is strange)
156 kmap = struct();
157 for i=1:k_len
158 key = sprintf('a%d', int32(k_bins+round(ktraj_adc(:,i)/k_threshold)));
159 assert(length(key)<=63);
160 if isfield(kmap,key)
161 k_storage_ind = kmap.(key);
162 else
163 k_storage_ind=k_storage_next;
164 kmap.(key)=k_storage_ind;
165 k_storage_next=k_storage_next+1;
166 end
167 k_storage(k_storage_ind)=k_storage(k_storage_ind)+1;
168 k_repeat(i) = k_storage(k_storage_ind);
169 end
170 end
171 % % containers.Map only supports string as a key... (in older matlabs)
172 % % in Octave we use the built-in, in Matlab the Java version
173 % if mr.aux.isOctave()
174 % kmap = containers.Map('KeyType', 'char', 'ValueType', 'int32');
175 % for i=1:k_len
176 % key_string = sprintf('%d ', int32(k_bins+round(ktraj_adc(:,i)/k_threshold)));
177 % % containers.Map does not have a proper find function so we use direct
178 % % access and catch the possible error
179 % try
180 % k_storage_ind = kmap(key_string);
181 % catch
182 % k_storage_ind=k_storage_next;
183 % kmap(key_string)=k_storage_ind;
184 % k_storage_next=k_storage_next+1;
185 % end
186 % k_storage(k_storage_ind)=k_storage(k_storage_ind)+1;
187 % k_repeat(i) = k_storage(k_storage_ind);
188 % end
189 % else
190 % kmap = java.util.HashMap;
191 % for i=1:k_len
192 % key_string = sprintf('%d ', int32(k_bins+round(ktraj_adc(:,i)/k_threshold)));
193 % k_storage_ind = kmap.get(key_string);
194 % if isempty(k_storage_ind)
195 % k_storage_ind=k_storage_next;
196 % kmap.put(key_string,k_storage_ind);
197 % k_storage_next=k_storage_next+1;
198 % end
199 % k_storage(k_storage_ind)=k_storage(k_storage_ind)+1;
200 % k_repeat(i) = k_storage(k_storage_ind);
201 % end
202 % end
203 % at this point k_storage(1:(k_storage_next-1)) is our visit frequency map
204 Repeats_max=max(k_storage(1:(k_storage_next-1)));
205 Repeats_min=min(k_storage(1:(k_storage_next-1)));
206 Repeats_median=median(k_storage(1:(k_storage_next-1)));
207 Repeats_unique=unique(k_storage(1:(k_storage_next-1)));
208 Counts_unique=zeros(size(Repeats_unique));
209 for i=1:numel(Repeats_unique)
210 Counts_unique(i)=sum(Repeats_unique(i)==k_storage(1:(k_storage_next-1)));
211 end
213 ktraj_rep1=ktraj_adc(:,k_repeat==1);
214 % TODO: think of something clever, e.g. detecting maximum delta-k
215 % if length(k_extent)==2
216 % dt = delaunayTriangulation(ktraj_rep1');
217 % k = convexHull(dt);
218 % figure; plot(dt.Points(:,1),dt.Points(:,2), '.', 'markersize',10); hold on;
219 % plot(dt.Points(k,1),dt.Points(k,2), 'r'); hold off;
220 % %[V,R] = voronoiDiagram(dt);
221 % figure; voronoi(dt);
222 % end
223 % try to detect k-space lines or columns.
224 k_counters=zeros(size(ktraj_rep1));
225 dims=size(ktraj_rep1,1);
226 %ordering=cell(1,dims);
227 if hasDict
228 for j=1:dims
229 kmap=configureDictionary("int32","int32");%dictionary(int32.empty(1,0),int32.empty(1,0));%configureDictionary("int32","int32"); % containers.Map('KeyType', 'int32', 'ValueType', 'int32'); % works as well but is substantially slower
230 k_storage=zeros(1,k_len);
231 k_storage_next=1;
232 for i=1:size(ktraj_rep1,2)
233 key=int32(round(ktraj_rep1(j,i)/k_threshold));
234 k_storage_ind = kmap.lookup(key,'FallbackValue',int32(0));
235 if k_storage_ind==0
236 % attempt to account for rounding errors
237 k_storage_ind = kmap.lookup(key+1,'FallbackValue',int32(0));
238 if k_storage_ind==0
239 k_storage_ind = kmap.lookup(key-1,'FallbackValue',int32(0));
240 % did not find anywhere...
241 if k_storage_ind==0
242 k_storage_ind=k_storage_next;
243 kmap(key)=k_storage_ind;
244 k_storage_next=k_storage_next+1;
245 k_storage(k_storage_ind)=ktraj_rep1(j,i);
246 %fprintf('%d:%d(%g) ',k_storage_ind,key,ktraj_rep1(j,i));
247 end
248 end
249 end
250 %assert(k_storage_ind==k_storage(k_storage_ind));
251 k_counters(j,i) = k_storage_ind;
252 end
253 %ordering{j}=cell2mat(kmap.values);
254 %fprintf('\n');
255 end
256 else
257 for j=1:dims
258 k_storage=zeros(1,k_len);
259 k_storage_next=1;
260 kmap = struct(); % use struct() as a replacement for dict()
261 for i=1:size(ktraj_rep1,2)
262 key=int32(round(ktraj_rep1(j,i)/k_threshold));
263 skey=makeSkey(key);
264 try
265 k_storage_ind = kmap.(skey);
266 catch
267 skey1=makeSkey(key+1); % attempt to account for rounding errors
268 try
269 k_storage_ind = kmap.(skey1);
270 catch
271 skey1=makeSkey(key-1); % attempt to account for rounding errors
272 try
273 k_storage_ind = kmap.(skey1);
274 catch
275 k_storage_ind=k_storage_next;
276 kmap.(skey)=k_storage_ind;
277 k_storage_next=k_storage_next+1;
278 k_storage(k_storage_ind)=ktraj_rep1(j,i);
279 %fprintf('%d:%d(%g) ',k_storage_ind,key,ktraj_rep1(j,i));
280 end
281 end
282 end
283 %assert(k_storage_ind==k_storage(k_storage_ind));
284 k_counters(j,i) = k_storage_ind;
285 end
286 % for i=1:size(ktraj_rep1,2)
287 % key=int32(round(ktraj_rep1(j,i)/k_threshold));
288 % skey=makeSkey(key);
289 % if isfield(kmap,skey)
290 % k_storage_ind = kmap.(skey);
291 % else
292 % skey1=makeSkey(key+1); % attempt to account for rounding errors
293 % if isfield(kmap,skey1)
294 % k_storage_ind = kmap.(skey1);
295 % else
296 % skey1=makeSkey(key-1); % attempt to account for rounding errors
297 % if isfield(kmap,skey1)
298 % k_storage_ind = kmap.(skey1);
299 % else
300 % k_storage_ind=k_storage_next;
301 % kmap.(skey)=k_storage_ind;
302 % k_storage_next=k_storage_next+1;
303 % k_storage(k_storage_ind)=ktraj_rep1(j,i);
304 % %fprintf('%d:%d(%g) ',k_storage_ind,key,ktraj_rep1(j,i));
305 % end
306 % end
307 % end
308 % %assert(k_storage_ind==k_storage(k_storage_ind));
309 % k_counters(j,i) = k_storage_ind;
310 % end
311 %ordering{j}=cell2mat(kmap.values);
312 %fprintf('\n');
313 end
314 end
315 unique_kpositions=max(k_counters,[],2);
316 isCartesian=(prod(unique_kpositions)==size(ktraj_rep1,2));
317else
318 unique_kpositions=1;
319end
320% check gradient amplitudes and slew rates
322% gradient waveform
323%gw_data=obj.waveforms_and_times(); % FIXME: avoid this second call for generating gradient shapes (1st one was inside of the k-space calculation routine)
324gws=cell(size(wnt.gw_data));
325ga=zeros(length(wnt.gw_data),1);
326gs=zeros(length(wnt.gw_data),1);
327% to calculate max absolute gradients and slew rates we have to play
328% tricks... we namely have to interpolate the data to the common time axis
329dim1ind = @(x, n) x(n,:);
330common_time=unique(dim1ind([wnt.gw_data{:}],1));
331gw_ct=zeros(length(wnt.gw_data),length(common_time));
332gs_ct=zeros(length(wnt.gw_data),length(common_time)-1);
333for gc=1:length(wnt.gw_data)
334 if size(wnt.gw_data{gc},2)>0
335 gws{gc}=(wnt.gw_data{gc}(2,2:end)-wnt.gw_data{gc}(2,1:end-1))./(wnt.gw_data{gc}(1,2:end)-wnt.gw_data{gc}(1,1:end-1)); % slew
336 % interpolate to the common time
337 gw_ct(gc,:)=interp1(wnt.gw_data{gc}(1,:),wnt.gw_data{gc}(2,:),common_time,'linear',0);
338 gs_ct(gc,:)=(gw_ct(gc,2:end)-gw_ct(gc,1:end-1))./(common_time(2:end)-common_time(1:end-1));
339 % max grad/slew per channel
340 ga(gc)=max(abs(wnt.gw_data{gc}(2,:)));
341 gs(gc)=max(abs(gws{gc}));
342 % TODO: calculate grad RMS values (this is an interesting task in the piece-wise-linear domain)
343 end
344end
346%figure; plot(common_time, gw_ct');
347%figure; plot(common_time, sum(gw_ct.^2,1).^0.5);
348%figure; plot(0.5*(common_time(1:end-1)+common_time(2:end)), sum(gs_ct.^2,1).^0.5);
350% max absolute value grad/slew -- check for a worst case upon rotation
351if ~isempty(gw_ct)
352 ga_abs=max(sum(gw_ct.^2,1).^0.5);
353else
354 ga_abs=0;
355end
356if ~isempty(gs_ct)
357 gs_abs=max(sum(gs_ct.^2,1).^0.5);
358else
359 gs_abs=0;
360end
362% check timing of blocks and delays (raster alignment)
363[timing_ok, timing_error_report] = obj.checkTiming();
365report = { sprintf('Number of blocks: %d\n',numBlocks),...
366 [ sprintf( 'Number of events:\n'),...
367 optionalOut(' RF: %6d\n',eventCount(2)),...
368 optionalOut(' Gx: %6d\n',eventCount(3)),...
369 optionalOut(' Gy: %6d\n',eventCount(4)),...
370 optionalOut(' Gz: %6d\n',eventCount(5)),...
371 optionalOut(' ADC: %6d\n',eventCount(6))],...
372 [ sprintf( 'Event library use:\n'),...
373 optionalOut(' RF: %6d\n',numel(obj.rfLibrary.keys)),...
374 optionalOut(' Grad: %6d\n',numel(obj.gradLibrary.keys)),...
375 optionalOut(' Shape: %6d\n',numel(obj.shapeLibrary.keys)),...
376 optionalOut(' ADC: %6d\n',numel(obj.adcLibrary.keys)),...
377 optionalOut(' Extn: %6d\n',numel(obj.extensionLibrary.keys)),...
378 optionalOut(' Trigg: %6d\n',numel(obj.trigLibrary.keys)),...
379 optionalOut(' Label: %6d\n',numel(obj.labelsetLibrary.keys)+numel(obj.labelincLibrary.keys)),...
380 optionalOut(' RfShm: %6d\n',numel(obj.rfShimLibrary.keys)),...
381 optionalOut(' Rot: %6d\n',numel(obj.rotationLibrary.keys)),...
382 optionalOut(' SoDel: %6d\n',numel(obj.softDelayLibrary.keys))],...
383 [ sprintf('Sequence duration: %.6fs\n',duration),...
384 sprintf('TE: %.6fs\n',TE),...
385 sprintf('TR: %.6fs\n',TR) ],...
386 sprintf('Flip angle: %.02f°\n', flipAnglesDeg),...
387 sprintf('Unique k-space positions (a.k.a. columns, rows, etc): %d\n', unique_kpositions)};
388if any(unique_kpositions>1)
389 report = { report{:},...
390 [ sprintf('Dimensions: %d\n', length(k_extent)),...
391 sprintf(' Spatial resolution: %.02f mm\n', 0.5./k_extent*1e3) ],...
392 sprintf('Repetitions/slices/contrasts: %.d range: [%.d %.d]\n', Repeats_median, Repeats_min, Repeats_max) };
393 report = { report{:},...
394 sprintf(' %d k-space position(s) repeated %d times\n', [Counts_unique;Repeats_unique])};
396 if isCartesian
397 report = { report{:}, sprintf('Grid-like/Cartesian encoding trajectory detected\n') };
398 else
399 report = { report{:}, sprintf('Non-Cartesian/irregular encoding trajectory detected (e.g. spiral, radial, some EPI, etc)\n') };
400 end
401end
402if (timing_ok)
403 report = { report{:}, sprintf('Block timing check passed successfully\n') };
404else
405 report = { report{:}, [ sprintf('Block timing check failed! Error listing follows:\n'),...
406 sprintf([timing_error_report{:}]) ] };
407end
408msg_ga='';
409if ~isempty(opt.system) && any(ga > opt.system.maxGrad)
410 msg_ga=' [some component EXCEEDED]';
411end
412msg_gs='';
413if ~isempty(opt.system) && any(gs > opt.system.maxSlew)
414 msg_gs=' [some component EXCEEDED]';
415end
416report = { report{:},...
417 sprintf(['Max. Gradient: %.0f Hz/m == %.02f mT/m' msg_ga '\n'], [ga mr.convert(ga,'Hz/m','mT/m')]'),...
418 sprintf(['Max. Slew Rate: %g Hz/m/s == %.02f T/m/s' msg_gs '\n'], [gs mr.convert(gs,'Hz/m/s','T/m/s')]') };
419msg_ga='';
420if ~isempty(opt.system) && ga_abs > opt.system.maxGrad
421 msg_ga=' [EXCEEDED]';
422end
423msg_gs='';
424if ~isempty(opt.system) && gs_abs > opt.system.maxSlew
425 msg_gs=' [EXCEEDED]';
426end
427report = { report{:},...
428 sprintf(['Max. Absolute Gradient: %.0f Hz/m == %.02f mT/m' msg_ga '\n'], [ga_abs mr.convert(ga_abs,'Hz/m','mT/m')]'),...
429 sprintf(['Max. Absolute Slew Rate: %g Hz/m/s == %.02f T/m/s' msg_gs '\n'], [gs_abs mr.convert(gs_abs,'Hz/m/s','T/m/s')]') };
431end