Merge pull request #39 from filip-szczepankiewicz/fsz_dev_crosscomp

New feature: cross-term compensation

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,:);

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
+

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;
+

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

optimizationProblem.m

@@ -37,6 +37,8 @@
         MaxFunEval = 1e5;
         MaxIter    = 5e3;
         motionCompensation = struct('order', [], 'maxMagnitude', [], 'linear', [])
+        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
 
     properties (SetAccess = private)
@@ -77,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
+            %% 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)
@@ -111,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

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;
@@ -85,23 +85,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
@@ -110,6 +115,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);
@@ -120,23 +126,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];

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,22 +165,36 @@
 
 [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,:);
 
 % Motion compensation
 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 > 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)

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);