Merge pull request #35 from jsjol/mixed-motion-comp

Support for mixed linear- and nonlinear motion constraints

NOW_GUI.m

@@ -794,6 +794,7 @@ function motionCompensationOrder0_Callback(hObject, eventdata, handles)
 % Hint: get(hObject,'Value') returns toggle state of motionCompensationOrder0
 order = []; % order zero is always ensured by the spin echo condition
 handles.problem.motionCompensation.order = order;
+handles.problem.motionCompensation.maxMagnitude = zeros(size(order));
 handles.problem = optimizationProblem(handles.problem);
 set(handles.motionCompensationOrder0, 'Value', true);
 set(handles.motionCompensationOrder1, 'Value', false);
@@ -811,7 +812,7 @@ function motionCompensationOrder1_Callback(hObject, eventdata, handles)
 % Hint: get(hObject,'Value') returns toggle state of motionCompensationOrder1
 order = 1; % order zero is always ensured by the spin echo condition
 handles.problem.motionCompensation.order = order;
-handles.problem.motionCompensation.linear = true;
+handles.problem.motionCompensation.maxMagnitude = zeros(size(order));
 handles.problem = optimizationProblem(handles.problem);
 set(handles.motionCompensationOrder0, 'Value', false);
 set(handles.motionCompensationOrder1, 'Value', true);
@@ -828,7 +829,7 @@ function motionCompensationOrder2_Callback(hObject, eventdata, handles)
 % Hint: get(hObject,'Value') returns toggle state of motionCompensationOrder2
 order = [1, 2]; % order zero is always ensured by the spin echo condition
 handles.problem.motionCompensation.order = order;
-handles.problem.motionCompensation.linear = true;
+handles.problem.motionCompensation.maxMagnitude = zeros(size(order));
 handles.problem = optimizationProblem(handles.problem);
 set(handles.motionCompensationOrder0, 'Value', false);
 set(handles.motionCompensationOrder1, 'Value', false);

optimizationProblem.m

@@ -36,12 +36,12 @@
         MaxwellIndex = 100;
         MaxFunEval = 1e5;
         MaxIter    = 5e3;
-        motionCompensation = struct('order', [], 'maxMagnitude', [], 'linear', true)
+        motionCompensation = struct('order', [], 'maxMagnitude', [], 'linear', [])
     end
 
     properties (SetAccess = private)
         zeroGradientAtIndex = [];
-        tolIsotropy = .5e-2; %before 1e-4
+        tolIsotropy = .5e-2;
         tolMaxwell
         signs
         tolSlew
@@ -78,6 +78,7 @@
             obj.sMaxConstraint = obj.sMax*obj.dt^2;
             obj.integralConstraint = obj.eta*obj.gMaxConstraint^2*obj.totalTimeActual/obj.dt;
 
+            %% Maxwell compensation
             obj.tolMaxwell = obj.MaxwellIndex/obj.dt; %
 
             if ~isempty(obj.zeroGradientAtIndex) && obj.doMaxwellComp
@@ -98,8 +99,19 @@
                 obj.signs = zeros(obj.N - 1,1)+eps; % setting to zero flips out due to sqrt(0)=complex (??)
             end
 
+            %% Motion compensation
+            if length(obj.motionCompensation.maxMagnitude) ~= length(obj.motionCompensation.order)
+                error('motionCompensation.maxMagnitude must have the same size as motionCompensation.order.')                
+            end
+
+            if isempty(obj.motionCompensation.maxMagnitude)
+                obj.motionCompensation.linear = [];
+            else
+                % Infer empty motionCompensation.linear from values of
+                % motionCompensation.maxMagnitude
+                obj.motionCompensation.linear = (obj.motionCompensation.maxMagnitude <= 0);
+            end
         end
     end
-
 end
 

private/nonlconAnalytic.m

@@ -106,31 +106,30 @@
 % Magn Reson Med. 2020
 
 % Non-linear constraints for motion compensation
-if ~motionCompensation.linear && ~isempty(motionCompensation.order)
-    c7 = zeros(1, length(motionCompensation.order));
-    dc7_dx = zeros(length(motionCompensation.order), 3*N + 1);
-
-    % 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.
-    for i = 1:length(motionCompensation.order)
-        order = motionCompensation.order(i);
-        moment_weighting = - order * dt * t.^(order-1);
-        moment_vector = moment_weighting * Q;
-        c7(i) = sum(moment_vector.^2) - (motionCompensation.maxMagnitude(i) * 1000^order / (gamma * 1e-6))^2;
-        dc7_dx(i, 1:(3*N)) = 2 * kron(moment_vector, moment_weighting);
-    end
-
-else
+nonlinear_ind = find(~motionCompensation.linear);
+if isempty(nonlinear_ind)
     c7 = [];
     dc7_dx = [];
-end
+else
+    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 
+        % 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
 
 ceq     = [];
 gradceq = [];

private/optimize.m

@@ -163,20 +163,19 @@
 [firstDerivativeMatrix, ~] = getDerivativeMatrices(problem);
 
 Aeq = zeros(2 + length(problem.zeroGradientAtIndex) + ...
-    problem.motionCompensation.linear * length(problem.motionCompensation.order),problem.N);
+    nnz(problem.motionCompensation.linear), 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
-if problem.motionCompensation.linear
-    t = ((1:problem.N)-1/2) * problem.dt;
-    for i = 1:length(problem.motionCompensation.order)
-        order = problem.motionCompensation.order(i);
+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);
-    end
 end
 
 % Enforce symmetry about zero gradient interval

scripted_NOW_Example.m

@@ -51,14 +51,11 @@
 % to set this parameter.
 problem.MaxwellIndex = 100; %In units of (mT/m)^2 ms
 
-% Set the desired orders of motion compensation, and whether to use linear
-% constraints (requiring moments of the desired orders to be zero) or 
-% nonlinear constraints with corresponding thresholds for allowed
-% deviations. Please see Szczepankiewicz et al., MRM, 2020 for further
-% details. maxMagnitude in units s^order / m.
+% Set the desired orders of motion compensation and corresponding
+% thresholds for allowed deviations. See Szczepankiewicz et al., MRM, 2020
+% for details. maxMagnitude in units s^order / m.
 problem.motionCompensation.order = [1, 2];
-problem.motionCompensation.linear = false;
-problem.motionCompensation.maxMagnitude = [1e-4, 1e-4];
+problem.motionCompensation.maxMagnitude = [0, 1e-4];
 
 % Make a new optimizationProblem object using the updated specifications.
 % This explicit call is necessary to update all private variables.