1function [B, m1, m2, m3] = calcMomentsBtensor(obj, varargin)
2%give input arguments of calcB,calcm1,calcm2,calcm3 with true or false and
3%Ndummy with 0,1,... as the Number of Dummy scans (scans that should be
4%skipped for the calculation)
6%m1/m1 are of the structure array(R,3), where R denotes the number of the
7%repetition readout (for single shot sequences) with the three orientations
8%x,y,z respectively
9%B is of a similar structure with array(R,3,3)
10%CONSIDER: This so far only works for one slice or rather one measurement PER TR
11defaultB = true;
12defaultM = false;
13defaultD = 0;
14%temporalTolerance=1e-9; % 1ns
16persistent p;
17if isempty(p)
18 p=mr.aux.InputParserCompat;
19 addRequired(p, 'obj');
21 % addParameter(p, 'calcB', defaultB);
22 % addParameter(p, 'calcm1', defaultM);
23 % addParameter(p, 'calcm2', defaultM);
24 % addParameter(p, 'calcm3', defaultM);
25 % addParameter(p, 'Ndummy', defaultD);
27 % Add optional arguments with default values
28 addOptional(p, 'calcB', defaultB);
29 addOptional(p, 'calcm1', defaultM);
30 addOptional(p, 'calcm2', defaultM);
31 addOptional(p, 'calcm3', defaultM);
32 addOptional(p, 'Ndummy', defaultD);
33end
35% Parse the input
36parse(p, obj, varargin{:});
38% Assign the results to variables
39calcB = p.Results.calcB;
40calcm1 = p.Results.calcm1;
41calcm2 = p.Results.calcm2;
42calcm3 = p.Results.calcm3;
43Ndummy = p.Results.Ndummy;
45%%
46% calcB = true;
47% calcm1 = true;
48% calcm2 = true;
49% calcm3 = true;
50% Ndummy = 0;
51%% first get gradient shape and timing properties out of the sequence
52warning('OFF', 'mr:restoreShape')
53[~, ~, ~, ~, t_excitation, t_refocusing, ~, ~, gw_pp] = obj.calculateKspacePP();
55R = length(t_excitation); %repetions defined over excitations
57gx_pp_l = gw_pp{1};
58gy_pp_l = gw_pp{2};
59gz_pp_l = gw_pp{3};
60%% splitting pp-functions into equal parts
61tSeq = [];
62t_echo = [];
63for i=1:R
64 t_echo(end+1)=(2*t_refocusing(i) - t_excitation(i)); % TODO: fixme for double-refocused sequences
65 tSeq = [tSeq [t_excitation(i), t_refocusing(i), t_echo(i)]];
66end
68t1 = gx_pp_l.breaks;
69t2 = gy_pp_l.breaks;
70t3 = gz_pp_l.breaks;
71tn = unique([t1,t2,t3,tSeq]);
73% % debuging / visualization
74% tnew = linspace(0,tn(end),10000);
75% figure; hold on;
76% plot(tnew,ppval(gx_pp_l,tnew));
77% plot(tnew,ppval(gy_pp_l,tnew));
78% plot(tnew,ppval(gz_pp_l,tnew));
80gx_pp_coefs=fillPpCoefs(gx_pp_l,tn);
81gx_pp_l = mkpp(tn,gx_pp_coefs);
82gy_pp_coefs=fillPpCoefs(gy_pp_l,tn);
83gy_pp_l = mkpp(tn,gy_pp_coefs);
84gz_pp_coefs=fillPpCoefs(gz_pp_l,tn);
85gz_pp_l = mkpp(tn,gz_pp_coefs);
87%% splitting pp at multiples of TR
88n = (1+Ndummy);
89gz_pp = cell(1,R);
90gy_pp = cell(1,R);
91gx_pp = cell(1,R);
92for i = 1:size(gz_pp_l.breaks,2)
93 if n == R
94 break
95 elseif abs(gz_pp_l.breaks(i) - t_excitation(n+1)) == 0 %skipping first excitation
96 gz_pp{n-Ndummy} = fnbrk(gz_pp_l,[t_excitation(n) t_excitation(n+1)]);
97 gy_pp{n-Ndummy} = fnbrk(gy_pp_l,[t_excitation(n) t_excitation(n+1)]);
98 gx_pp{n-Ndummy} = fnbrk(gx_pp_l,[t_excitation(n) t_excitation(n+1)]);
99 n = n +1;
100 end
101end
102%add the last block as well
103gx_pp{R}= fnbrk(gx_pp_l,[t_excitation(end) gx_pp_l.breaks(end)]);
104gy_pp{R}= fnbrk(gy_pp_l,[t_excitation(end) gy_pp_l.breaks(end)]);
105gz_pp{R}= fnbrk(gz_pp_l,[t_excitation(end) gz_pp_l.breaks(end)]);
108%% considering effects of the rf pulse -> effective gradients
109%i.e gradient becomes negative after
110%the refocusing pulse
111%as gradients are split at excitation, nothing necessary here
112for j=1:R
113 for i = 1:size(gx_pp{j}.coefs,1)
114 if gx_pp{j}.breaks(:,i) == t_refocusing(j+Ndummy)
115 gx_pp{j}.coefs(i:end,:) = -gx_pp{j}.coefs(i:end,:);
116 gy_pp{j}.coefs(i:end,:) = -gy_pp{j}.coefs(i:end,:);
117 gz_pp{j}.coefs(i:end,:) = -gz_pp{j}.coefs(i:end,:);
118 end
119 end
120end
121%% B-Tensor calculation
122B = zeros(R,3,3);
123if calcB
125 qz_pp = cell(1,R);
126 qy_pp = cell(1,R);
127 qx_pp = cell(1,R);
129 for i = 1:R
130 qx_pp{i}=compat_fnint(gx_pp{i});
131 qy_pp{i}=compat_fnint(gy_pp{i});
132 qz_pp{i}=compat_fnint(gz_pp{i});
134 qx_pp{i}.coefs = qx_pp{i}.coefs*2*pi;
135 qy_pp{i}.coefs = qy_pp{i}.coefs*2*pi;
136 qz_pp{i}.coefs = qz_pp{i}.coefs*2*pi;
137 end
139 % tnew = linspace(0,gx_pp_l.breaks(end),1000);
140 % figure;
141 % plot(tnew,ppval(qz_pp{2},tnew),'LineWidth',2,'Color','r');
142 % % hold on;
143 % % plot(tnew,ppval(qy_pp{1},tnew),'LineWidth',2,'Color','b');
144 % %plot(tnew,ppval(qz_pp{2},tnew),'LineWidth',2,'Color','y');
145 % hold off
146 % title('wave vector')
147 % legend('kx','ky','kz')
148 % xlabel('Time [ms]')
149 % ylabel('Wave Vector [1/m]')
150 % grid;
152 for m = 1:R
153 q_pp{1} = qx_pp{m};
154 q_pp{2} = qy_pp{m};
155 q_pp{3} = qz_pp{m};
156 for i = 1:3
157 for j = 1:3
158 coefs=zeros(size(qx_pp{m}.coefs,1),5);
159 for k=1:size(qx_pp{m}.coefs,1)
160 coefs(k,:)=conv(q_pp{i}.coefs(k,:),q_pp{j}.coefs(k,:));
161 end
162 Bpp_div = mkpp(qx_pp{m}.breaks,coefs);
163 Bpp = compat_fnint(Bpp_div);
164 %evaluate b-value at TE
165 B(m,i,j) = ppval(Bpp,tSeq(3 + 3*(m-1+Ndummy)));
166 end
167 end
168 end
169end
171%% calculate first moment m1
172m1 = zeros(R,3);
174if calcm1
176 m1z_pp = cell(1,R);
177 m1y_pp = cell(1,R);
178 m1x_pp = cell(1,R);
179 for i = 1:R
180 %defining a new piecewise polynomial aquivalent to g(t)=t
181 t_pp_coefs = zeros(length(gx_pp{i}.breaks)-1,2);
182 t_pp_coefs(:,1) = 1;
183 t_pp_coefs(:,2) = gx_pp{i}.breaks(1:(end-1))-gx_pp{i}.breaks(1);
185 new_cx = zeros(size(gx_pp{i}.coefs,1),size(gx_pp{i}.coefs,2)+1);
186 new_cy = zeros(size(gy_pp{i}.coefs,1),size(gy_pp{i}.coefs,2)+1);
187 new_cz = zeros(size(gz_pp{i}.coefs,1),size(gz_pp{i}.coefs,2)+1);
189 for k=1:size(t_pp_coefs,1)
190 new_cx(k,:)=conv(gx_pp{i}.coefs(k,:),t_pp_coefs(k,:));
191 new_cy(k,:)=conv(gy_pp{i}.coefs(k,:),t_pp_coefs(k,:));
192 new_cz(k,:)=conv(gz_pp{i}.coefs(k,:),t_pp_coefs(k,:));
193 end
195 tgx_pp = mkpp(gx_pp{i}.breaks, new_cx);
196 tgy_pp = mkpp(gy_pp{i}.breaks, new_cy);
197 tgz_pp = mkpp(gz_pp{i}.breaks, new_cz);
199 m1x_pp{i}=compat_fnint(tgx_pp);
200 m1y_pp{i}=compat_fnint(tgy_pp);
201 m1z_pp{i}=compat_fnint(tgz_pp);
202 m1x_pp{i}.coefs = m1x_pp{i}.coefs*2*pi;
203 m1y_pp{i}.coefs = m1y_pp{i}.coefs*2*pi;
204 m1z_pp{i}.coefs = m1z_pp{i}.coefs*2*pi;
205 end
207 %as an integration from 0 to TE, I will evaluate the function at TE to
208 %gain the first order vector
210 for m=1:R
211 m1(m,1) = ppval(m1x_pp{m},tSeq(3 + 3*(m-1+Ndummy)));%-ppval(m1x_pp{m},m1x_pp{m}.breaks(1));
212 m1(m,2) = ppval(m1y_pp{m},tSeq(3 + 3*(m-1+Ndummy)));%-ppval(m1y_pp{m},m1y_pp{m}.breaks(1));
213 m1(m,3) = ppval(m1z_pp{m},tSeq(3 + 3*(m-1+Ndummy)));%-ppval(m1z_pp{m},m1z_pp{m}.breaks(1));
214 end
216end
217%% calculate second moment m2
218m2 = zeros(R,3);
220if calcm2
222 m2z_pp = cell(1,R);
223 m2y_pp = cell(1,R);
224 m2x_pp = cell(1,R);
225 for i = 1:R
226 t2_pp_coefs = zeros(length(gx_pp{i}.breaks)-1,3);
227 t2_pp_coefs(:,1) = 1;
228 t2_pp_coefs(:,2) = 2*(gx_pp{i}.breaks(1:(end-1))-gx_pp{i}.breaks(1));
229 t2_pp_coefs(:,3) = (gx_pp{i}.breaks(1:(end-1))-gx_pp{i}.breaks(1)).^2;
231 new_cx = zeros(size(gx_pp{i}.coefs,1),size(gx_pp{i}.coefs,2)+2);
232 new_cy = zeros(size(gy_pp{i}.coefs,1),size(gy_pp{i}.coefs,2)+2);
233 new_cz = zeros(size(gz_pp{i}.coefs,1),size(gz_pp{i}.coefs,2)+2);
235 for k=1:size(t2_pp_coefs,1)
236 new_cx(k,:)=conv(gx_pp{i}.coefs(k,:),t2_pp_coefs(k,:));
237 new_cy(k,:)=conv(gy_pp{i}.coefs(k,:),t2_pp_coefs(k,:));
238 new_cz(k,:)=conv(gz_pp{i}.coefs(k,:),t2_pp_coefs(k,:));
239 end
241 tgx_pp = mkpp(gx_pp{i}.breaks, new_cx);
242 tgy_pp = mkpp(gy_pp{i}.breaks, new_cy);
243 tgz_pp = mkpp(gz_pp{i}.breaks, new_cz);
245 m2x_pp{i}=compat_fnint(tgx_pp);
246 m2y_pp{i}=compat_fnint(tgy_pp);
247 m2z_pp{i}=compat_fnint(tgz_pp);
248 m2x_pp{i}.coefs = m2x_pp{i}.coefs*2*pi;
249 m2y_pp{i}.coefs = m2y_pp{i}.coefs*2*pi;
250 m2z_pp{i}.coefs = m2z_pp{i}.coefs*2*pi;
251 end
253 %as an integration from 0 to TE, I will evaluate the function at TE to
254 %gain the first order vector
256 for m=1:R
257 m2(m,1) = ppval(m2x_pp{m},tSeq(3 + 3*(m-1+Ndummy)))-ppval(m2x_pp{m},m2x_pp{m}.breaks(1));
258 m2(m,2) = ppval(m2y_pp{m},tSeq(3 + 3*(m-1+Ndummy)))-ppval(m2y_pp{m},m2y_pp{m}.breaks(1));
259 m2(m,3) = ppval(m2z_pp{m},tSeq(3 + 3*(m-1+Ndummy)))-ppval(m2z_pp{m},m2z_pp{m}.breaks(1));
260 end
262end
265%% calc third moment m3:
266m3 = zeros(R,3);
268if calcm3
270 m3z_pp = cell(1,R);
271 m3y_pp = cell(1,R);
272 m3x_pp = cell(1,R);
273 for i = 1:R
274 t3_pp_coefs = zeros(length(gx_pp{i}.breaks)-1,4);
275 t3_pp_coefs(:,1) = 1;
276 t3_pp_coefs(:,2) = 3*(gx_pp{i}.breaks(1:(end-1))-gx_pp{i}.breaks(1));
277 t3_pp_coefs(:,3) = 3*(gx_pp{i}.breaks(1:(end-1))-gx_pp{i}.breaks(1)).^2;
278 t3_pp_coefs(:,4) = (gx_pp{i}.breaks(1:(end-1))-gx_pp{i}.breaks(1)).^3;
279 % t3_pp_breaks = gx_pp{i}.breaks;
280 % t3_pp = mkpp(t3_pp_breaks,t3_pp_coefs);
281 % tnew = linspace(0,gx_pp_l.breaks(end),1000);
282 % figure;
283 % plot(tnew,ppval(t3_pp,tnew),'LineWidth',2,'Color','r');
284 % hold on
285 % plot(tnew,(tnew-gx_pp{i}.breaks(1)).^3,'Color','b');
286 % hold off
288 new_cx = zeros(size(gx_pp{i}.coefs,1),size(gx_pp{i}.coefs,2)+3);
289 new_cy = zeros(size(gy_pp{i}.coefs,1),size(gy_pp{i}.coefs,2)+3);
290 new_cz = zeros(size(gz_pp{i}.coefs,1),size(gz_pp{i}.coefs,2)+3);
292 for k=1:size(t3_pp_coefs,1)
293 new_cx(k,:)=conv(gx_pp{i}.coefs(k,:),t3_pp_coefs(k,:));
294 new_cy(k,:)=conv(gy_pp{i}.coefs(k,:),t3_pp_coefs(k,:));
295 new_cz(k,:)=conv(gz_pp{i}.coefs(k,:),t3_pp_coefs(k,:));
296 end
298 tgx_pp = mkpp(gx_pp{i}.breaks, new_cx);
299 tgy_pp = mkpp(gy_pp{i}.breaks, new_cy);
300 tgz_pp = mkpp(gz_pp{i}.breaks, new_cz);
302 m3x_pp{i}=compat_fnint(tgx_pp);
303 m3y_pp{i}=compat_fnint(tgy_pp);
304 m3z_pp{i}=compat_fnint(tgz_pp);
305 m3x_pp{i}.coefs = m3x_pp{i}.coefs*2*pi;
306 m3y_pp{i}.coefs = m3y_pp{i}.coefs*2*pi;
307 m3z_pp{i}.coefs = m3z_pp{i}.coefs*2*pi;
308 end
310 %as an integration from 0 to TE, I will evaluate the function at TE to
311 %gain the first order vector
313 for m=1:R
314 m3(m,1) = ppval(m3x_pp{m},tSeq(3 + 3*(m-1+Ndummy)))-ppval(m3x_pp{m},m3x_pp{m}.breaks(1));
315 m3(m,2) = ppval(m3y_pp{m},tSeq(3 + 3*(m-1+Ndummy)))-ppval(m3y_pp{m},m3y_pp{m}.breaks(1));
316 m3(m,3) = ppval(m3z_pp{m},tSeq(3 + 3*(m-1+Ndummy)))-ppval(m3z_pp{m},m3z_pp{m}.breaks(1));
317 end
319end
320end
321%% functions
322function pp1_coefs=fillPpCoefs(pp1,xn)
323 idx1 = slookup(xn(1:end-1),pp1.breaks(1:end-1));
324 pp1_coefs = zeros(length(xn)-1,pp1.order);
325 for i=1:size(pp1_coefs,1)
326 if idx1(i)>0
327 % simple copy
328 pp1_coefs(i,:)=pp1.coefs(idx1(i),:);
329 elseif i > 1
330 % copy from the left with a reference shift (else leave at 0)
331 for k=0:(pp1.order-1)
332 for l=0:k
333 pp1_coefs(i,end-l) = pp1_coefs(i,end-l) + pp1_coefs(i-1,end-k)*nchoosek(k,l)*(xn(i)-xn(i-1))^(k-l);
334 end
335 end
336 end
337 end
338end
340function idx=slookup(what,where)
341% finds indices of values given by the sorted vector 'what' in the sorted vector 'where'
342% for failing searches indices of 0 are returned
343 idx=zeros(size(what));
344 wb=1; % where bound
345 for c=1:length(what)
346 i=find(what(c)==where(wb:end),1);
347 if isempty(i), continue; end
348 idx(c)=wb+i-1;
349 wb=wb+i;
350 end
351end
353function idx=sintlookup(what,where)
354% finds indices of intervals to which the sorted vector 'what' belongs in the sorted vector 'where'
355% for failing searches indices of 0 are returned
356 idx=zeros(size(what));
357 wb=1; % where bound
358 for c=1:length(what)
359 i=find(what(c)>=where(wb:end));
360 if isempty(i), continue; end
361 idx(c)=wb+i(end)-1;
362 wb=idx(c);
363 end
364end
366function pp_out = compat_fnint(pp_in)
367% Octave-compatible wrapper: use ppint (built-in) on Octave, fnint on MATLAB
368 if mr.aux.isOctave()
369 pp_out = ppint(pp_in);
370 else
371 pp_out = fnint(pp_in);
372 end
373end