From 6da6171ca32d180ebe9ed22c9e0395fa036ef9a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jens=20Sj=C3=B6lund?= Date: Wed, 11 Nov 2020 09:22:25 +0100 Subject: [PATCH 1/9] Initial implementation of background compensation as a linear equality constraint. --- optimizationProblem.m | 3 ++- private/optimize.m | 15 ++++++++++++++- scripted_NOW_Example.m | 4 ++++ 3 files changed, 20 insertions(+), 2 deletions(-) diff --git a/optimizationProblem.m b/optimizationProblem.m index f85085e..cae83f5 100644 --- a/optimizationProblem.m +++ b/optimizationProblem.m @@ -37,6 +37,7 @@ MaxFunEval = 1e5; MaxIter = 5e3; motionCompensation = struct('order', [], 'maxMagnitude', [], 'linear', []) + doBackgroundCompensation = false end properties (SetAccess = private) @@ -81,7 +82,7 @@ %% Maxwell compensation obj.tolMaxwell = obj.MaxwellIndex/obj.dt; % - if ~isempty(obj.zeroGradientAtIndex) && obj.doMaxwellComp + if ~isempty(obj.zeroGradientAtIndex) && (obj.doMaxwellComp || obj.doBackgroundCompensation) signs = ones(obj.N - 1,1); % Ghost points excluded during opt signs(obj.zeroGradientAtIndex(end) + 1:end) = -1; obj.signs = signs; diff --git a/private/optimize.m b/private/optimize.m index 53c7348..7e7d329 100644 --- a/private/optimize.m +++ b/private/optimize.m @@ -162,11 +162,15 @@ [firstDerivativeMatrix, ~] = getDerivativeMatrices(problem); +% Allocate matrix Aeq that will act on a single component Aeq = zeros(2 + length(problem.zeroGradientAtIndex) + ... - nnz(problem.motionCompensation.linear), problem.N); + nnz(problem.motionCompensation.linear) + ... + 1 * problem.doBackgroundCompensation, problem.N); + % Require start and end in q-space origin (echo condition) Aeq(1,1)=1; Aeq(2,problem.N)=1; + % Require zero gradient at the specified indices Aeq(2+(1:length(problem.zeroGradientAtIndex)),:) = firstDerivativeMatrix(problem.zeroGradientAtIndex,:); @@ -177,6 +181,15 @@ order = problem.motionCompensation.order(linear_ind(i)); Aeq(2+length(problem.zeroGradientAtIndex)+i,:) = - order * problem.dt * t.^(order-1); end + +% Background compensation +if problem.doBackgroundCompensation == true + H = cumtrapz([1; problem.signs])'; % Safe to ignore dt because we'll equate to zero. + % Compute the integral of q*H using the trapezoidal rule (ignoring dt again): + Aeq(2 + length(problem.zeroGradientAtIndex) + length(linear_ind) + 1, 1) = H(1)/2; + Aeq(2 + length(problem.zeroGradientAtIndex) + length(linear_ind) + 1, 2:(end-1)) = H(2:(end-1)); + Aeq(2 + length(problem.zeroGradientAtIndex) + length(linear_ind) + 1, end) = H(end)/2; +end % Enforce symmetry about zero gradient interval if problem.enforceSymmetry == true diff --git a/scripted_NOW_Example.m b/scripted_NOW_Example.m index e5cd7b7..ed8bc50 100644 --- a/scripted_NOW_Example.m +++ b/scripted_NOW_Example.m @@ -16,6 +16,7 @@ % doMaxwellComp = true; % MaxwellIndex = 100; % Motion compensation: disabled +% Background compensation: disabled % % Written by Jens Sjölund and Filip Szczepankiewicz @@ -57,6 +58,9 @@ problem.motionCompensation.order = [1, 2]; problem.motionCompensation.maxMagnitude = [0, 1e-4]; +% Toggle compensation for background gradients +problem.doBackgroundCompensation = true; + % Make a new optimizationProblem object using the updated specifications. % This explicit call is necessary to update all private variables. problem = optimizationProblem(problem); From a9f39ca7851b2120d63f064caf2e2aa161dd5d98 Mon Sep 17 00:00:00 2001 From: Filip Szczepankiewicz Date: Thu, 4 Mar 2021 13:49:41 +0100 Subject: [PATCH 2/9] Create now_cross_term_sensitivity.m --- now_cross_term_sensitivity.m | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 now_cross_term_sensitivity.m diff --git a/now_cross_term_sensitivity.m b/now_cross_term_sensitivity.m new file mode 100644 index 0000000..673343a --- /dev/null +++ b/now_cross_term_sensitivity.m @@ -0,0 +1,12 @@ +function [n, nt, H] = now_cross_term_sensitivity(result) +% function [n, nt, H] = now_cross_term_sensitivity(result) + +gwf = result.g; +rf = result.rf; +dt = result.dt; + +qt = now_gwf_to_q(gwf, dt); + +H = now_gamma * cumsum(rf)*dt; +nt = cumsum(qt.*H)*dt; +n = nt(end,:); From 87f97b87cd58bfb6812322bb24c1352c78719ac2 Mon Sep 17 00:00:00 2001 From: Filip Szczepankiewicz Date: Thu, 4 Mar 2021 13:50:01 +0100 Subject: [PATCH 3/9] Create now_gwf_to_q.m --- now_gwf_to_q.m | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 now_gwf_to_q.m diff --git a/now_gwf_to_q.m b/now_gwf_to_q.m new file mode 100644 index 0000000..b412edd --- /dev/null +++ b/now_gwf_to_q.m @@ -0,0 +1,7 @@ +function q = now_gwf_to_q(gwf, dt) +% function q = now_gwf_to_q(gwf, dt) +% +% Assume gwf is given as the effective waveform. + +q = 2 * pi * now_gamma * cumsum(gwf, 1) * dt; + From d7e7610e90cd6affdbdb02608d3685016ec148f6 Mon Sep 17 00:00:00 2001 From: Filip Szczepankiewicz Date: Thu, 4 Mar 2021 13:50:16 +0100 Subject: [PATCH 4/9] Create now_gamma.m --- now_gamma.m | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 now_gamma.m diff --git a/now_gamma.m b/now_gamma.m new file mode 100644 index 0000000..223d8a5 --- /dev/null +++ b/now_gamma.m @@ -0,0 +1,18 @@ +function gamma = now_gamma(nuc) +% function gamma = now_gamma(nuc) +% +% Returns the gyromagnetic constant for given nucleus in units of 1/s/T. + +if nargin < 1 + nuc = 'H'; +end + + +switch nuc + case 'H' + gamma = 42.6e6; + + otherwise + error('Nucleus not recognized') +end + From 3e43411536c67df8a6ca7bb9b1e85c4940c80421 Mon Sep 17 00:00:00 2001 From: Filip Szczepankiewicz Date: Thu, 4 Mar 2021 13:55:58 +0100 Subject: [PATCH 5/9] Update optimizationProblem.m Added specific behavior for specific and general optimization Updated calculation of sign vector to have reversal in the middle of pause time. --- optimizationProblem.m | 91 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 73 insertions(+), 18 deletions(-) diff --git a/optimizationProblem.m b/optimizationProblem.m index cae83f5..997d3f2 100644 --- a/optimizationProblem.m +++ b/optimizationProblem.m @@ -37,7 +37,8 @@ MaxFunEval = 1e5; MaxIter = 5e3; motionCompensation = struct('order', [], 'maxMagnitude', [], 'linear', []) - doBackgroundCompensation = false + doBackgroundCompensation = 2; % 0 = off; 1 = general timing cond.; 2 = specific timing cond. + startTime = 0; % Time from the excitataion (t=0) to the first gradietn waveform sample in ms. end properties (SetAccess = private) @@ -78,31 +79,39 @@ obj.gMaxConstraint = obj.gMax*obj.dt; obj.sMaxConstraint = obj.sMax*obj.dt^2; obj.integralConstraint = obj.eta*obj.gMaxConstraint^2*obj.totalTimeActual/obj.dt; + obj.tolMaxwell = obj.MaxwellIndex/obj.dt; + + % Turn Maxwell compensation off if requested + % Can also be written as a more confusing and compact form: + % obj.tolMaxwell = obj.MaxwellIndex/obj.dt / (obj.doMaxwell>0); + % But doMaxwell can also be entirely replaced by + % obj.MaxwellIndex where < inf indicates that its ON but this will + % complicate the GUI. + if ~obj.doMaxwellComp + obj.tolMaxwell = inf; + end - %% Maxwell compensation - obj.tolMaxwell = obj.MaxwellIndex/obj.dt; % - if ~isempty(obj.zeroGradientAtIndex) && (obj.doMaxwellComp || obj.doBackgroundCompensation) + %% Create spin dephasing direction vector + if ~isempty(obj.zeroGradientAtIndex) signs = ones(obj.N - 1,1); % Ghost points excluded during opt - signs(obj.zeroGradientAtIndex(end) + 1:end) = -1; - obj.signs = signs; - else - % Maxwell terms cannot be compensated if no 180 pulses are - % used. Normally this kind of optimization is intended for - % a repetition of two identical self-balanced waveforms, - % but it may be necessary to warn users that single-sided - % experiments will always incurr some error due to - % concomitant fields. In practice the "weight" of the - % optimization with respect to Maxwell terms is removed by - % setting the sign vector to all zeros. - obj.doMaxwellComp = false; - obj.signs = zeros(obj.N - 1,1)+eps; % setting to zero flips out due to sqrt(0)=complex (??) + % Assume that sign change happens in the middle of the pause + mi = median(obj.zeroGradientAtIndex); + + signs(round(mi):end) = -1; + + if (mi==round(mi)) + signs(round(mi)) = 0; + end + + obj.signs = signs; end + %% Motion compensation if length(obj.motionCompensation.maxMagnitude) ~= length(obj.motionCompensation.order) - error('motionCompensation.maxMagnitude must have the same size as motionCompensation.order.') + error('motionCompensation.maxMagnitude must have the same size as motionCompensation.order.') end if isempty(obj.motionCompensation.maxMagnitude) @@ -112,6 +121,52 @@ % motionCompensation.maxMagnitude obj.motionCompensation.linear = (obj.motionCompensation.maxMagnitude <= 0); end + + %% Cross-term-compensation + switch obj.doBackgroundCompensation + case 0 % No compensation + % Do nothing + + case 1 % General timing condition + % This case requires velocity compensation. + + ind_velo = find(obj.motionCompensation.order==1, 1); + + % Check if velocity compensation is requested at all. + % If not, force linear constraint. + % Else, check that requested velo compensation is + % linear. If not, throw error since instructions cannot + % be fulfilled. + + if isempty(ind_velo) + ind = numel(obj.motionCompensation.order)+1; + obj.motionCompensation.order(ind) = 1; + obj.motionCompensation.linear(ind) = 1; + obj.motionCompensation.maxMagnitude(ind) = 0; + else + if ~obj.motionCompensation.linear(ind_velo) + error('Cross-term-compensation for a general timing requires velocity compensation!') + end + end + + case 2 % Specific timing condition + % Specific timing condition requires that the start + % time is set. The start time can be zero, but it is + % unlikely that this is an accurate setup, so we throw + % a warning. + + if obj.startTime == 0 + warning('Start time for waveform is t = 0, which is an unlikely setting. Please check!') + end + + if obj.startTime < 0 + error('Start time cannot be smaller than zero!') + end + + otherwise + error('Selection for Cross-term-compensation not recognized! Use value 0, 1 or 2.') + end + end end end From 31dd940a8f0e8c7a208683750453e9371ce79b00 Mon Sep 17 00:00:00 2001 From: Filip Szczepankiewicz Date: Thu, 4 Mar 2021 14:01:30 +0100 Subject: [PATCH 6/9] Update optimize.m Updated export of sign function (rf) to the mddmri format. Added support for arbitrary starting time. --- private/optimize.m | 40 ++++++++++++++++++++++------------------ 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/private/optimize.m b/private/optimize.m index 7e7d329..4d4efa9 100644 --- a/private/optimize.m +++ b/private/optimize.m @@ -81,7 +81,7 @@ %% Evaluate and store results [firstDerivativeMatrix, secondDerivativeMatrix] = getDerivativeMatrices(problem); -gamma = 42.6e6*2*pi; +gamma = now_gamma*2*pi; q = reshape(x(1:3*problem.N),[problem.N,3]); % remove noise from output - disabled awaiting improved handling of @@ -112,17 +112,20 @@ result.optimizationTime = optimizationTime; result.iter = iter-1; result.rawq = x(1:(end-1)); +result.zind = problem.zeroGradientAtIndex; -if problem.doMaxwellComp - rf = problem.signs; -else - rf = ones(size(problem.signs)); + +% Create output that is compatible with mddMRI framework +rf = ones(size(g,1), 1); + +if ~isempty(problem.zeroGradientAtIndex) + rf = [problem.signs(1); problem.signs; problem.signs(end)]; + rf = rf/max(abs(rf)); % This is dirty, but signs is ill used when mxwl and cross comp is off. end -result.zind = problem.zeroGradientAtIndex; % Keep this info for save function -result.rf = [rf(1); rf; rf(end)]; % Spin direction -result.gwf = bsxfun(@times, result.rf, g/1000); % T/m -result.dt = problem.dt/1000; % s +result.rf = rf; % Spin dephasing direction +result.gwf = (g/1000) .* repmat(rf, 1, 3); % T/m +result.dt = problem.dt/1000; % s end @@ -162,10 +165,10 @@ [firstDerivativeMatrix, ~] = getDerivativeMatrices(problem); -% Allocate matrix Aeq that will act on a single component +% Allocate matrix Aeq that will act on a single component Aeq = zeros(2 + length(problem.zeroGradientAtIndex) + ... - nnz(problem.motionCompensation.linear) + ... - 1 * problem.doBackgroundCompensation, problem.N); + nnz(problem.motionCompensation.linear) + ... + 1 * problem.doBackgroundCompensation, problem.N); % Require start and end in q-space origin (echo condition) Aeq(1,1)=1; @@ -178,19 +181,20 @@ t = ((1:problem.N)-1/2) * problem.dt; linear_ind = find(problem.motionCompensation.linear); for i = 1:length(linear_ind) - order = problem.motionCompensation.order(linear_ind(i)); - Aeq(2+length(problem.zeroGradientAtIndex)+i,:) = - order * problem.dt * t.^(order-1); + order = problem.motionCompensation.order(linear_ind(i)); + Aeq(2+length(problem.zeroGradientAtIndex)+i,:) = - order * problem.dt * t.^(order-1); end % Background compensation -if problem.doBackgroundCompensation == true - H = cumtrapz([1; problem.signs])'; % Safe to ignore dt because we'll equate to zero. - % Compute the integral of q*H using the trapezoidal rule (ignoring dt again): +if problem.doBackgroundCompensation > 0 + s = problem.startTime; % ms + H = cumsum([1; problem.signs])' * problem.dt + s; % Safe to ignore dt because we'll equate to zero. + % Compute the integral of q*H using the trapezoidal rule (ignoring dt again): Aeq(2 + length(problem.zeroGradientAtIndex) + length(linear_ind) + 1, 1) = H(1)/2; Aeq(2 + length(problem.zeroGradientAtIndex) + length(linear_ind) + 1, 2:(end-1)) = H(2:(end-1)); Aeq(2 + length(problem.zeroGradientAtIndex) + length(linear_ind) + 1, end) = H(end)/2; end - + % Enforce symmetry about zero gradient interval if problem.enforceSymmetry == true if isempty(problem.zeroGradientAtIndex) From 7e752a0ea960bc2351bcc5bbd4aa5d1acaf06c9d Mon Sep 17 00:00:00 2001 From: Filip Szczepankiewicz Date: Thu, 4 Mar 2021 14:02:01 +0100 Subject: [PATCH 7/9] Update now_maxwell_coeff.m --- now_maxwell_coeff.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/now_maxwell_coeff.m b/now_maxwell_coeff.m index 59cc073..85d6ed8 100644 --- a/now_maxwell_coeff.m +++ b/now_maxwell_coeff.m @@ -42,7 +42,7 @@ end function k = get_k_matrix(g, rf, B0, dt) -gammaOver2Pi = 42.6e6; % 1 / (T s) +gammaOver2Pi = now_gamma; % 1 / (T s) % Waveform moment vector scale matrix % k-vector along direction [x y z] for one part of the waveform (pre or From db82da0426451ce2401681cc59e0a31e0a51834e Mon Sep 17 00:00:00 2001 From: Filip Szczepankiewicz Date: Thu, 4 Mar 2021 14:03:18 +0100 Subject: [PATCH 8/9] Update nonlconAnalytic.m Updated how the maxwell condition is fed to optimizer Minor update to reference gamma factor needs to be updated throughout --- private/nonlconAnalytic.m | 57 ++++++++++++++++++++++----------------- 1 file changed, 32 insertions(+), 25 deletions(-) diff --git a/private/nonlconAnalytic.m b/private/nonlconAnalytic.m index 1b65202..e6af640 100644 --- a/private/nonlconAnalytic.m +++ b/private/nonlconAnalytic.m @@ -8,7 +8,7 @@ % Remove action on axes that are turned off - disabled awaiting further % tests or improved handling of disabled axes -% Q = Q .* ((diag(targetTensor))>0)'; +Q = Q .* ((diag(targetTensor))>0)'; integrationWeights = ones(N,1); integrationWeights(1) = 0.5; @@ -79,23 +79,28 @@ % Magn Reson Med. 2019, Vol. 82, Issue 4, p.1424-1437. % https://doi.org/10.1002/mrm.27828 -signedg = bsxfun(@times, g, signs); -M = g'*signedg; -m = sqrt(trace(M'*M)); % 'Maxwell index' -c6 = m - tolMaxwell; % Check whether to square or not (as in ref) - -dc6_dM = 1/m * M; -% M is a "matrix quadratic form", so we use the same procedure as above, -% but for performance reasons we keep it inline instead of -% defining a function -weightedQ = firstDerivativeMatrix'*signedg; -firstTerm = kron(eye(3), weightedQ'); -secondTerm = reshape(firstTerm, [3,3, 3*N]); -secondTerm = permute(secondTerm, [2, 1, 3]); -secondTerm = reshape(secondTerm, [9, 3*N]); -dM_dq = firstTerm + secondTerm; -dc6_dq = reshape(dc6_dM, [1, 9]) * dM_dq; -dc6_dx = [dc6_dq, 0]; +if isinf(tolMaxwell) + c6 = []; + dc6_dx = []; +else + signedg = bsxfun(@times, g, signs); + M = g'*signedg; + m = sqrt(trace(M'*M)); % 'Maxwell index' + c6 = m - tolMaxwell; % Check whether to square or not (as in ref) + + dc6_dM = 1/m * M; + % M is a "matrix quadratic form", so we use the same procedure as above, + % but for performance reasons we keep it inline instead of + % defining a function + weightedQ = firstDerivativeMatrix'*signedg; + firstTerm = kron(eye(3), weightedQ'); + secondTerm = reshape(firstTerm, [3,3, 3*N]); + secondTerm = permute(secondTerm, [2, 1, 3]); + secondTerm = reshape(secondTerm, [9, 3*N]); + dM_dq = firstTerm + secondTerm; + dc6_dq = reshape(dc6_dM, [1, 9]) * dM_dq; + dc6_dx = [dc6_dq, 0]; +end % Constraint for ballistic motion encoding @@ -104,6 +109,7 @@ % Motion-compensated gradient waveforms for tensor-valued % diffusion encoding by constrained numerical optimization % Magn Reson Med. 2020 +% https://doi.org/10.1002/mrm.28551 % Non-linear constraints for motion compensation nonlinear_ind = find(~motionCompensation.linear); @@ -114,23 +120,24 @@ for i = 1:length(nonlinear_ind) c7 = zeros(1, length(nonlinear_ind)); dc7_dx = zeros(length(nonlinear_ind), 3*N + 1); - - % dt is in ms; % Behaves better if calculation uses ms but - % requested tolerance has some strange units. Fixed this by + + % dt is in ms; % Behaves better if calculation uses ms but + % requested tolerance has some strange units. Fixed this by % rescaling the maxMagnitude by 1000^order. - + t = ((1:N)-1/2) * dt; - + gamma = 2.6751e+08; % radians / T / s for hydrogen. - + order = motionCompensation.order(nonlinear_ind(i)); moment_weighting = - order * dt * t.^(order-1); moment_vector = moment_weighting * Q; c7(i) = sum(moment_vector.^2) - (motionCompensation.maxMagnitude(nonlinear_ind(i)) * 1000^order / (gamma * 1e-6))^2; dc7_dx(i, 1:(3*N)) = 2 * kron(moment_vector, moment_weighting); - end + end end + ceq = []; gradceq = []; c = [c1 c2 c3 c4 c5 c6 c7]; From 9049835fec8bba6c06ae5d66a077aff6d5008892 Mon Sep 17 00:00:00 2001 From: Filip Szczepankiewicz Date: Sat, 6 Mar 2021 15:54:38 +0100 Subject: [PATCH 9/9] Reset default to non-bck-comp --- optimizationProblem.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/optimizationProblem.m b/optimizationProblem.m index 997d3f2..1327970 100644 --- a/optimizationProblem.m +++ b/optimizationProblem.m @@ -37,7 +37,7 @@ MaxFunEval = 1e5; MaxIter = 5e3; motionCompensation = struct('order', [], 'maxMagnitude', [], 'linear', []) - doBackgroundCompensation = 2; % 0 = off; 1 = general timing cond.; 2 = specific timing cond. + doBackgroundCompensation = 0; % 0 = off; 1 = general timing cond.; 2 = specific timing cond. startTime = 0; % Time from the excitataion (t=0) to the first gradietn waveform sample in ms. end