1classdef TransformFOV < handle
3 properties (Access = public)
4 rotation=[]; % 3x3 rotation matrix or empty matrix if none
5 translation=[]; % 1x3 vector or empty matrix if none
6 scale=[]; % 1x3 vector or empty if none
7 end
9 properties (Access = private)
10 use_rotation_extension=false; % whether rotations should be explicitly applied to the gradients or implicitly by means of RotationExtension
11 prior_phase_cycle=0;
12 rotation_quaternion=[]; % rotation quaternion version of the provided rotation matrix or empty matrix if none
13 system;
14 % high_accuracy;
15 labels=struct('NOPOS',0,'NOROT',0,'NOSCL',0);
16 end
18 methods
20 function obj = TransformFOV(varargin)
21 % Creates an instance of the TransformFOV object, which can
22 % then be applied to shift, rotate or scale the imaging volume
23 % either of a single block, range of blocks or the entire
24 % sequence. The optional input parameters are:
25 % rotation: 3x3 rotation matrix
26 % translation: 1x3 translation vector in the original
27 % (non-rotated) logical Pulseq coordinates
28 % scale: 1x3 scailing vector in the original
29 % (non-rotated) logical Pulseq coordinates. The
30 % provided vector is the *gradient* scailing
31 % vector, so the FOV size is *divided* by the
32 % provided values. Scaling of 0 on one or
33 % several axes sets gradients on these axes to 0
34 % transform: 4x4 homogeneous transform matrix, used as an
35 % alternative to a combination of 'rotation' and
36 % 'translation'. If 'tansform' is specified
37 % neither 'translation' nor 'rotation' can be
38 % used. However, 'scale' parameter can be used
39 % in combination with 'transform'. In contrast
40 % to 'translation', the translation part of the
41 % 'transform' matrix is in the world/lab (e.g.
42 % new, rotated) coordinates in accordance with
43 % the standard homogeneous transform definition
44 % use_rotation_extension: defines whether the gradient events
45 % are rotated immediately by the function, or
46 % whether the rotation extension should be used
47 % and the actual rotation will then be applied
48 % by the interpreter. If the input sequence
49 % readily uses rotation extension and
50 % 'use_rotation_extension' is False, the
51 % gradient events will be transformed by the
52 % stored rotation extension information prior to
53 % further calculations
54 % prior_phase_cycle: allows to specify previously accumulated
55 % phase, e.g. for combining different instances
56 % of the FOV position classes
57 %
58 % The order of transformations is as follows:
59 % scale, translation, rotation. Translation and rotation can be
60 % mixed in into the homogeneous transform matrix, but that
61 % must exclude scaling, which then may still be provided
62 % separately
63 %
65 persistent parser
66 if isempty(parser)
67 parser = inputParser;
68 parser.FunctionName = 'TransformFOV_f';
70 addParameter(parser, 'rotation', [], @(m) (isnumeric(m) && all(size(m)==[3 3])));
71 addParameter(parser, 'translation', [], @(m) (isnumeric(m) && size(m,2)==3 && length(m)==3));
72 addParameter(parser, 'scale', [], @(m) (isnumeric(m) && size(m,2)==3 && length(m)==3));
73 addParameter(parser, 'transform', [], @(m) (isnumeric(m) && all(size(m)==[4 4])));
74 addParameter(parser, 'use_rotation_extension', false, @(m) (islogical(m) && numel(m)==1));
75 addParameter(parser, 'prior_phase_cycle', 0, @(m) (isnumeric(m) && numel(m)==1));
76 addParameter(parser, 'system', [], @isstruct);
77 end
79 parse(parser, varargin{:});
80 opt = parser.Results;
82 if ~isempty(opt.transform)
83 if ~isempty(opt.rotation) || ~isempty(opt.translation)
84 error('Neither ''translation'' nor ''rotation'' can be provided in combination with the ''transfrom'' option');
85 end
86 opt.rotation=opt.transform(1:3,1:3);
87 off=opt.transform(4,1:3); % TODO: check whether this is indeed the column
88 opt.translation=opt.rotation*off'; % TODO: check the direction of the rotation (or inverse)
89 %off=opt.transform(1:3,4); % alternative - translation column of the standard 4x4 homogeneous matrix
90 %opt.translation=(opt.rotation*off)'; % TODO: check the direction of the rotation (or inverse)
91 elseif isempty(opt.rotation) && isempty(opt.translation) && isempty(opt.scale)
92 error('At least one transforming parameter needs to be provided');
93 end
95 obj.rotation=opt.rotation;
96 obj.translation=opt.translation;
97 obj.scale=opt.scale;
98 obj.prior_phase_cycle=opt.prior_phase_cycle;
99 obj.use_rotation_extension=opt.use_rotation_extension;
101 if obj.use_rotation_extension && ~isempty(obj.rotation)
102 obj.rotation_quaternion=mr.aux.quat.fromRotMat(obj.rotation);
103 end
105 if isempty(opt.system)
106 obj.system=mr.opts();
107 else
108 obj.system=opt.system;
109 end
110 end
112 function out=applyToBlock(obj, varargin)
114 % convert the input into a plain cell array of events
115 if ~any(iscell(varargin))
116 block_events=varargin;
117 else
118 block_events={};
119 for i=varargin
120 if iscell(i)
121 block_events=[block_events i];
122 else
123 block_events{end+1}=i;
124 end
125 end
126 end
128 % see if we've got a block as a single struct, then convert it to a cell array
129 if (isstruct(block_events) && isfield(block_events, 'blockDuration')) || ...
130 (iscell(block_events) && ~isempty(block_events) && isstruct(block_events{1}) && isfield(block_events{1}, 'blockDuration'))
131 block_events=mr.block2events(block_events);
132 end
134 % extract various mr events including {rf,adc,gx,gy,gz} from "block_events" input
135 rf = [];
136 adc = [];
137 grads = cell(1,3);
138 other = {};
139 rotExtQuaternion = [];
140 for i = 1:length(block_events)
141 e=block_events{i};
142 if length(e)==1 && isstruct(e) && isfield(e, 'type')
143 switch e.type
144 case 'rf'
145 rf = e; % save rf event
146 case 'adc'
147 adc = e; % save adc event
148 case {'trap', 'grad'} % if gradient event, check 'channel'
149 if ~isfield(e, 'channel')
150 error('unspecified gradient channel for the gradient object');
151 end
152 switch e.channel
153 case 'x'
154 grads{1} = e;
155 case 'y'
156 grads{2} = e;
157 case 'z'
158 grads{3} = e;
159 otherwise
160 error('unsupported gradient channel %s for the gradient object', e.channel);
161 end
162 case 'labelset' %{'labelset', 'labelinc'} % we dont really need 'labelinc', as all labels that are important for us are flags and have no 'inc'
163 for j=1:length(e)
164 switch e(j).label % this switch has only one case on purpose, it is just a lazy way of checking that we deal with a relevant label setting
165 case {'NOPOS','NOROT','NOSCL'}
166 obj.labels(e(j).label)=e(j).value;
167 end
168 end
169 case 'rot3D'
170 rotExtQuaternion=e.rotQuaternion;
171 otherwise
172 other{end+1}=e;
173 end
174 else
175 other{end+1}=e;
176 %other=[other num2cell(e)]; % e can be an array of structs, and funny enough, num2cell can convert it to a cell array...
177 end
178 end
180 gradRasterTime = obj.system.gradRasterTime;
182 %% scale (apply it first prior to any other transformation)
183 if ~isempty(obj.scale)
184 %channel2index=struct('x',1,'y',2,'z',3);
185 for i=1:length(grads)
186 %grads{i}=mr.scaleGrad(grads{i},obj.scale(channel2index.(grads{i}.channel)),obj.system);
187 if ~isempty(grads{i})
188 if (isfield(grads{i},'id'))
189 grads{i}=mr.scaleGrad(rmfield(grads{i},'id'),obj.scale(i),obj.system);
190 else
191 grads{i}=mr.scaleGrad(grads{i},obj.scale(i),obj.system);
192 end
193 end
194 end
195 end
197 %% translation
198 if ~isempty(obj.translation)
199 % big picture of the algorithm
201 % if ~isempty(obj.translation)
202 % phase_cycle_this_block = 0
203 % if NOPOS==0
204 % apply prior_phase_cycle to rf-adc
205 % end
206 % for i=1:3
207 % g = grad{i};
208 % if ~isempty(g)
209 % if translation(i)~=0
210 % generate piecewise polynomial of the current gradient
211 % if NOPOS==0
212 % apply phase and freq offsets to rf-adc by the current gradient
213 % end
214 % update phase_cycle_this_block by the current gradient
215 % end
216 % end
217 % end
218 % update prior_phase_cycle =+ phase_cycle_this_block
219 % end
221 % extract the first and the last time points of possible rf or adc in the block
222 if isempty(rf)
223 if isempty(adc)
224 % both ADC and RF are not defined
225 t_end = [];
226 t_start = [];
227 else
228 % only ADC is defined
229 [t_start,t_end] = extract_time(adc);
230 end
231 else % RF is defined
232 if isempty(adc)
233 % only RF is defined
234 [t_start,t_end] = extract_time(rf);
235 else
236 % both ADC and RF are defined
237 [t_start_adc,t_end_adc] = extract_time(adc);
238 [t_start_rf,t_end_rf] = extract_time(rf);
239 t_start = min(t_start_adc,t_start_rf);
240 t_end = max(t_end_adc,t_end_rf);
241 end
242 end
244 % remove IDs because we will change the objects below
245 if ~isempty(rf) && isfield(rf,'id')
246 rf=rmfield(rf,'id');
247 end
248 if ~isempty(adc) && isfield(adc,'id')
249 adc=rmfield(adc,'id');
250 end
252 % if the current block uses rotation extension and the
253 % TransformFOV object is configured not to use the rotation
254 % extension then we apply the rotation to the gradients now
255 if ~obj.use_rotation_extension && ~isempty(rotExtQuaternion)
256 grads=mr.rotate3D(rotExtQuaternion,grads,'system',obj.sys);
257 rotExtQuaternion=[]; % now that we have applied the current rotation, we can discard it
258 end
260 % MZ: I think we have to rotate the gradient "backwards" if
261 % this block has 'NOROT'. We restore the grads object below
262 % WARNING: I don't think this is compatible with obj.use_rotation_extension
263 if obj.labels.NOROT
264 grads_backup=grads;
265 % MZ: HA! we could rotate obj.translation (or it's copy) in the opposite direction instead
266 % MZ: and, we could use the same mechanism to handle the rotation extention
267 grads=mr.rotate3D(obj.rotation',grads,'system',obj.sys); % MZ: I guess we have to rotate the gradients "back" because we are normally in local logical coordinates, which would be "rotated" if there were NOROT flag
268 % MZ: please check if the above point is correct
269 end
271 % define a temporary parameter for the phase cycle of the current block
272 phase_cycle_this_block=0;
274 % apply prior_phase_cycle to rf-adc only if NOPOS==0
275 if ~obj.labels.NOPOS % do fov positioning for the current block
276 % if there is rf in the current block, adjust its
277 % phase-offset
278 if ~isempty(rf)
279 rf.phaseOffset = rf.phaseOffset + 2*pi * obj.prior_phase_cycle;
280 end
281 % if there is adc in the current block, adjust its
282 % phase-offset
283 if ~isempty(adc)
284 adc.phaseOffset = adc.phaseOffset + 2*pi * obj.prior_phase_cycle;
285 end
286 end
288 % 1- take all gradintes in the current block and make
289 % piecewise polynimials for them
290 % 2- only if NOPOS==0 , calculate rf and adc phase and frequency offsets based on
291 % the single gradient
292 % 3- update the value of 'phase_cycle_this_block'
293 % 4- update the value of 'prior_phase_cycle'
294 for i=1:3
295 if abs(obj.translation(i))>eps % check whether there is shift in the specific direction (x or y or z)
296 g = grads{i};
297 if ~isempty(g) % check whether there is corresponding gradient in the direction of the shift (gx or gy or gz)
298 % 1- make pp
299 if strcmp(g.type , 'trap') % if g is a simple trapezoid
300 if (abs(g.flatTime)>eps) % interp1 gets confused by triangular gradients (repeating sample)
301 tt = g.delay+cumsum([0 g.riseTime g.flatTime g.fallTime]);
302 waveform = g.amplitude*[0 1 1 0];
303 else
304 if (abs(g.riseTime)>eps && abs(g.fallTime)>eps) % we skip 'empty' gradients
305 tt = g.delay+cumsum([0 g.riseTime g.fallTime]);
306 waveform = g.amplitude*[0 1 0];
307 else
308 if abs(g.amplitude)>eps
309 warning('''empty'' gradient with non-zero magnitude detected');
310 end
311 end
312 end
313 else % if g is a extended trapezoid or arbitrary gradient
314 tt = g.delay + g.tt;
315 waveform = g.waveform;
316 end
318 % generate breaks and coefs of the gradient required for
319 % making piecewise polynomial
320 [breaks, coefs, tt_extended, waveform_extended] = generate_breaks_coefs(g, tt, waveform, gradRasterTime, t_start, t_end);
322 % make piecewise polynomial for the gradient
323 f_pp = mkpp(breaks, coefs);
324 % integrate it analytically
325 if mr.aux.isOctave()
326 fi_pp = ppint(f_pp);
327 else
328 fi_pp = fnint(f_pp); % MZ: TODO: check whether we need more accurate functions here
329 end
332 if ~obj.labels.NOPOS
333 % apply adc or rf phase and freq offset contribution by only one (current) gradient
334 % (we separate all gradients in the block and apply effect of each one to the offsets)
335 event = {rf,adc};
336 for j=1:length(event)
337 e = event{j};
338 if ~isempty(e)
339 [t_s, t_e] = extract_time(e); % find the first and last time point of the rf or adc event
340 is_const = is_grad_const(tt_extended, waveform_extended, t_s, t_e); % check whether the gradient is constant during the rf or adc event: 1 means constant
341 if is_const % in case of constant gradient, we can easily adjust the frequency and phase offset of the rf or ADC event
342 freq = obj.translation(i) * ppval(f_pp, t_s);
343 if isfield(e,'t') %e=="rf"
344 rf.freqOffset = rf.freqOffset + freq;
345 phase_cycle = local_frac( accurate_mod_pp(fi_pp.breaks, fi_pp.coefs, t_s, obj.translation(i)) - freq * (t_s-rf.delay) );
346 rf.phaseOffset = rf.phaseOffset + 2*pi*phase_cycle;
347 else %if e=="adc"
348 adc.freqOffset = adc.freqOffset + freq;
349 phase_cycle = local_frac( accurate_mod_pp(fi_pp.breaks, fi_pp.coefs, t_s, obj.translation(i)) - freq * (t_s-adc.delay) );
350 adc.phaseOffset = adc.phaseOffset + 2*pi* phase_cycle;
351 end
352 else
353 % in case of non-constant gradient, we should calculate a phase vector for rf or ADC. for rf event we can easily add
354 % the phase vector to rf.signal and for ADC event we store the phase vector and use it in image reconstruction
355 if isfield(e,'t') %e=="rf"
356 % calculate the frequency at the center
357 ppval_f_center=ppval(f_pp, rf.delay+rf.center);
358 freq = obj.translation(i) * ppval_f_center;
359 rf.freqOffset = rf.freqOffset + freq;
360 % ppval_fi_center = ppval(fi_pp, rf.delay+rf.center);
361 ppval_fi_center_shift = accurate_mod_pp(fi_pp.breaks, fi_pp.coefs, rf.delay+rf.center, obj.translation(i));
362 phase_cycle = local_frac(ppval_fi_center_shift - freq*rf.center);
363 rf.phaseOffset = rf.phaseOffset + 2*pi*phase_cycle;
364 phase_cycle_vector_tmp = accurate_mod_pp(fi_pp.breaks, fi_pp.coefs, rf.t+rf.delay, obj.translation(i));
365 phase_cycle_vector = local_frac( phase_cycle_vector_tmp -ppval_f_center*(rf.t-rf.center) * obj.translation(i) - ppval_fi_center_shift );
366 rf.signal = rf.signal .* exp(1i*2*pi*phase_cycle_vector);
367 else % if e=="adc"
368 % calculate the frequency at the center
369 adc_center=0.5*adc.dwell*adc.numSamples;
370 ppval_f_center=ppval(f_pp, adc.delay+adc_center);
371 freq = obj.translation(i) * ppval_f_center;
372 adc.freqOffset = adc.freqOffset + freq;
373 % ppval_fi_center=ppval(fi_pp, adc.delay+adc_center);
374 ppval_fi_center_shift = accurate_mod_pp(fi_pp.breaks, fi_pp.coefs, adc.delay+adc_center, obj.translation(i));
375 % these -0.5 and +0.5 are needed to avoid unnecessay jumps for values that are very close to 0s
376 phase_cycle = local_frac( -0.5 + local_frac(0.5 + ppval_fi_center_shift - freq*adc_center) );
377 adc.phaseOffset = adc.phaseOffset + 2*pi*phase_cycle;
378 adc_t=adc.dwell*(0.5:adc.numSamples-0.5);
379 % these -0.5 and +0.5 are needed to avoid unnecessay jumps for values that are very close to 0s
380 phase_cycle_vector_tmp = accurate_mod_pp(fi_pp.breaks, fi_pp.coefs, adc_t+adc.delay, obj.translation(i));
381 phase_cycle_vector = local_frac( ( -0.5 + local_frac(0.5 - ppval_f_center*(adc_t-adc_center) ) ) * obj.translation(i) + ppval_fi_center_shift + phase_cycle_vector_tmp );
382 % phase_cycle_vector = accurate_mod_pp(fi_pp.breaks, fi_pp.coefs, adc.dwell*(0:adc.numSamples-1)+adc.delay, obj.translation(i));
383 if isempty(adc.phaseModulation)
384 adc.phaseModulation=2*pi* phase_cycle_vector(:); % store residual adc phase for image reconstruction
385 else
386 adc.phaseModulation=adc.phaseModulation+2*pi* phase_cycle_vector(:); % store residual adc phase for image reconstruction
387 end
388 end
389 end
390 end
391 end
392 end
394 % update the phase_cycle_this_block by only the current gradient
395 % phase_cycle_this_block = local_frac( phase_cycle_this_block + local_frac( g.area * obj.translation(i)));
396 phase_cycle_this_block = local_frac( phase_cycle_this_block + accurate_mod_pp(fi_pp.breaks, fi_pp.coefs, tt_extended(end), obj.translation(i)) );
397 end
398 end
399 end
401 % MZ: now restore the grads object
402 if obj.labels.NOROT
403 grads=grads_backup;
404 end
406 % now update the phase stored for the next block
407 obj.prior_phase_cycle = local_frac( obj.prior_phase_cycle + phase_cycle_this_block );
408 end
410 %% rotation
411 if ~isempty(obj.rotation) && ~obj.labels.NOROT
412 if obj.use_rotation_extension
413 if isempty(rotExtQuaternion)
414 rotExtQuaternion=obj.rotation_quaternion;
415 else
416 rotExtQuaternion=mr.aux.quat.multiply(rotExtQuaternion, obj.rotation_quaternion); % TODO: check left or right rotation
417 end
418 else
419 grads = mr.rotate3D(obj.rotation,grads,'system',obj.system);
420 end
421 end
423 %% rotation extension support
424 if ~isempty(rotExtQuaternion)
425 other{end+1}=mr.makeRotation(rotExtQuaternion);
426 end
428 %% output
429 out=[{rf} {adc} grads(:)' other(:)'];
430 out=out(~cellfun(@isempty,out)); % clean up empty
431 end
433 function seq2 = applyToSeq(obj, seq, varargin)
435 persistent parser
436 if isempty(parser)
437 parser = inputParser;
438 parser.FunctionName = 'applyToSeq';
440 %parser.addRequired('seq');
441 parser.addParamValue('sameSeq', false, @islogical); % TODO: add another option for an in-place transform
442 parser.addParamValue('blockRange',[1 inf],@(x)(isnumeric(x) && length(x)==2));
443 end
445 parse(parser, varargin{:});
446 opt = parser.Results;
448 if ~isfinite(opt.blockRange(2))
449 opt.blockRange(2)=length(seq.blockDurations);
450 end
452 if opt.sameSeq
453 seq2=seq;
454 else
455 seq2 = mr.Sequence(seq.sys);
456 seq2.copyDefinitions(seq); % copy definitions from the source sequence
457 end
459 obj.labels=struct('NOPOS',0,'NOROT',0,'NOSCL',0);
461 for iB=opt.blockRange(1):opt.blockRange(2)
462 B=seq.getBlock(iB,opt.sameSeq); % second parameter means 'addIDs'
463 B2=obj.applyToBlock(B);
464 seq2.addBlock(B2);
465 end
466 end
467 end
468end
470%% Local Functions
472% find the first and last time point of the rf or adc input in the block
473function [t_s, t_e] = extract_time(event)
474 if strcmp(event.type,'adc')
475 t_s = event.delay + event.dwell * 0.5;
476 t_e = event.delay + event.dwell * (event.numSamples-0.5);
477 elseif strcmp(event.type,'rf')
478 t_s = event.delay + event.t(1);
479 t_e = event.delay + event.t(end);
480 end
481end
483% check whether gradient is constant during rf or adc time points
485% inputs are: t=tt and amp=waveform of the gradient
486% and t_start and t_end are the first and last time point of the rf or adc event
487% flag=1 means the gradient is constant
488function is_const = is_grad_const(t, amp, t_start, t_end)
489 index_s = find(t <= t_start, 1, 'last');
490 index_e = find(t >= t_end, 1, 'first');
491 if isempty(index_s)
492 index_s=1;
493 end
494 if isempty(index_e)
495 index_e=length(t);
496 end
497 is_const = all(abs(amp(index_s:index_e) - amp(index_s)) <= 1e-10); % MZ: why this threshold? CHECKME
498end
501% generate breaks and coefficients for a gradient that is required for
502% making piecwise polynomial
503function [b, c, tt_extended, waveform_extended] = generate_breaks_coefs(g, tt, waveform, gradRasterTime, t_start, t_end)
504if g.type=='grad'
505 if abs(g.tt(1)-gradRasterTime/2)<eps % check whether gradient is arbitrary gradient or extended trapezoid
506 tt = [tt(1)-gradRasterTime/2; tt(:); tt(end)+gradRasterTime/2]; % if it's arbitrary gradient, we add first and last gradient amplitudes to the waveform
507 waveform = [g.first; waveform(:); g.last];
508 end
509end
510tt_extended = tt(:)';
511waveform_extended = waveform(:)';
513% generate initial breaks and coefs
514b = tt_extended;
515c = zeros(length(tt)-1,2);
516c(:,1) = diff(waveform) ./ diff(tt);
517c(find(isnan(c))) = 0;
518c(find(isinf(c))) = 0;
519c(:,2) = (waveform(1:end-1))';
521% MZ: question: why not adding (t_start,0) or (t_end,0) to tt/waveform? %
522% MS: if we add 0 at the begining/end of the waveform and then calculate the slope of the first/last time slot, we will get a non-zero slope (if g.first or g.last are not zero).
523% This is because the 0 value is connected to non-zero g.first or g.last and it means there is real gradient between them but should not exist.
524% so it's better to add [0,0] at the begining or end of the coef matrix later on.
525% modify breaks and coefs if the rf or adc event(s) is/are started (or ended)
526% before (or after) the first (or last) time point of the gradient
527if t_start<tt(1)
528 b = [t_start, b];
529 c = [0,0;c];
530 tt_extended = [t_start, tt_extended];
531 waveform_extended = [0, waveform_extended];
532end
533if t_end>tt(end)
534 b = [b, t_end];
535 c = [c;0,0];
536 tt_extended = [tt_extended, t_end];
537 waveform_extended = [waveform_extended, 0];
538end
539end
542% % make piecewise polynomial based on time points and amplitudes of gradient event
543% function f_pp = mkpp_linear(t,f)
544% breaks = t;
545% coefs = zeros(length(t)-1,2);
546% coefs(:,1) = diff(f) ./ diff(t);
547% coefs(find(isnan(coefs))) = 0;
548% coefs(find(isinf(coefs))) = 0;
549% coefs(:,2) = (f(1:end-1))';
550% f_pp = mkpp(breaks,coefs);
551% end
554% local function returning the fractional part
555function out=local_frac(in)
556 out = in-floor(in);
557end
559% % local mod function
560% function out=local_mod(in,m)
561% out = zeros(size(in));
562% for i=1:length(in)
563% if in(i)>0.5*m
564% z=floor(in(i)/m+0.5);
565% out(i)=in(i)-z*m;
566% elseif in(i)<=-0.5*m
567% z=floor(-in(i)/m+0.5);
568% out(i)=in(i)+z*m;
569% else
570% out(i)=in(i);
571% end
572% end
573% end
576%% Accurate mathematical operations using high resolution times (sum , multiply , mod)
577% in this function, we take the time breaks and coefficients of the
578% piecewise polynomial of the gradient and desired time point of the rf or ADC event (t) and
579% amount of the fov shift. and accurately calculate the area under each time
580% slot of the gradient and multiply it by shift and finally calculate the
581% mod of that by 1. (phase cycle)
582function Mod = accurate_mod_pp(breaks, coefs, t, shift)
583Mod = zeros(size(t));
584A = [];
585i_breaks = 0;
586for c=1:length(t)
587 % c
588 index = find(breaks <= t(c), 1, 'last');
589 t0 = t(c);
590 area = 0;
592 while index > (i_breaks + 1)
593 i_breaks = i_breaks +1;
594 AB = coefs(i_breaks,1) * shift * ( (breaks(i_breaks+1)-breaks(i_breaks))^2-(breaks(i_breaks)-breaks(i_breaks))^2 );
595 CD = coefs(i_breaks,2) * shift * ( (breaks(i_breaks+1)-breaks(i_breaks))-(breaks(i_breaks)-breaks(i_breaks)) );
596 A = [A,local_frac( local_frac(AB) + local_frac(CD) )];
597 end
599 if t0==breaks(i_breaks+1)
600 area = local_frac(sum(A));
601 else
602 AB_n = coefs(i_breaks+1,1) * shift * ( (t0-breaks(i_breaks+1))^2-(breaks(i_breaks+1)-breaks(i_breaks+1))^2 );
603 CD_n = coefs(i_breaks+1,2) * shift * ( (t0-breaks(i_breaks+1))-(breaks(i_breaks+1)-breaks(i_breaks+1)) );
604 area = local_frac( local_frac(sum(A)) + local_frac(AB_n) + local_frac(CD_n) );
605 end
606 Mod(c) = area;
607end
608end