/ concept-collection / seqlab
Sign in
concept-collection / seqlab
seqlab / src / engine / pulseq / +mr / traj2grad.m
78 lines · 2.6 KBBlameHistoryRaw
1function [g sr]=traj2grad(k,varargin)
2%traj2grad Convert a k-space trajectory to gradient waveform.
3% g=traj2grad(k) Convert k into gradient waveform g using finite
4% differences. The trajectory is in units of 1/m. The k-space points are
5% assumed to be sampled on the raster edges.
6% The size of k = [nChannel nTime].
7%
8% g=traj2grad(k,'RasterTime',T) Calculate gradient waveforms assuming
9% the given raster time.
11% See also Sequence.makeArbitraryGrad
13persistent parser
14if isempty(parser)
15 parser = inputParser;
16 parser.FunctionName = 'traj2grad';
17 parser.addRequired('k',@isnumeric);
18 parser.addParamValue('first',[],@isnumeric);
19 parser.addParamValue('firstGradStepHalfRaster',true,@islogical);
20 parser.addParamValue('conservativeSlewEstimate',false,@islogical);
21 parser.addParamValue('system',[],@isstruct);
22 parser.addParamValue('RasterTime',[],@isnumeric);
23end
24parse(parser,k,varargin{:});
25opt = parser.Results;
26if isempty(opt.system)
27 opt.system=mr.opts();
28end
29if isempty(opt.RasterTime)
30 opt.RasterTime=opt.system.gradRasterTime;
31end
32if isempty(opt.first)
33 opt.first=zeros(size(k,1),1); % QC: if the first gradient point is not given, set it to zero. 2025.01.03
34end
36% Compute finite difference for gradients in Hz/m
37%g=([k(:,2:end)-k(:,1:end-1) zeros(size(k,1),1)])/opt.RasterTime; % MZ: with zero-padding
38g=[(k(:,2:end)-k(:,1:end-1))/opt.RasterTime]; % MZ: no zero-padding!
40% Compute the slew rate (time derivative of the gradient)
41sr0=(g-[opt.first g(:,1:end-1)])/opt.RasterTime;
42if opt.firstGradStepHalfRaster
43 sr0(:,1)=sr0(:,1)*2; % account for the half-step in the beginning of the shape
44end
46% now we think how to post-process the results
47% gradient is now sampled between the k-points (on raster cell centers)
48% whilst the slew rate is between the gradient points, except of the first
49% point, which relies on the opt.first value (and may be a bit off anyway,
50% but this is the best estimate that we have)
51sr=zeros(size(sr0));
52sr(:,1)=sr0(:,1);
53if (opt.conservativeSlewEstimate)
54 if opt.firstGradStepHalfRaster
55 sr(:,2)=sr0(:,2);
56 sr(:,3:end)=max_abs(sr0(:,2:end-1),sr0(:,3:end));
57 else
58 sr(:,2:end)=max_abs(sr0(:,1:end-1),sr0(:,2:end));
59 end
60else
61 if opt.firstGradStepHalfRaster
62 sr(:,2)=sr0(:,2);
63 sr(:,3:end)=0.5*(sr0(:,2:end-1)+sr0(:,3:end));
64 else
65 sr(:,2:end)=0.5*(sr0(:,1:end-1)+sr0(:,2:end));
66 end
67end
69end
71function out=max_abs(in1, in2)
72 if size(in2)~=size(in2)
73 error('arrays of incompatible sizes');
74 end
75 abs1gtoe=abs(in1)>=abs(in2);
76 out=in1.*abs1gtoe+in2.*(~abs1gtoe);
77end
moveopenescclose