From d8b1f0a91ce402df93b7f79b66a27906b24732eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jens=20Sj=C3=B6lund?= Date: Tue, 20 Oct 2020 09:32:23 +0200 Subject: [PATCH 1/2] Support for mixed linear- and nonlinear motion constraints --- NOW_GUI.m | 2 +- optimizationProblem.m | 18 ++++++++++++++++- private/nonlconAnalytic.m | 41 +++++++++++++++++++-------------------- private/optimize.m | 11 +++++------ scripted_NOW_Example.m | 11 ++++------- 5 files changed, 47 insertions(+), 36 deletions(-) diff --git a/NOW_GUI.m b/NOW_GUI.m index f77aded..9ada5ba 100644 --- a/NOW_GUI.m +++ b/NOW_GUI.m @@ -828,7 +828,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.linear = [true, true]; handles.problem = optimizationProblem(handles.problem); set(handles.motionCompensationOrder0, 'Value', false); set(handles.motionCompensationOrder1, 'Value', false); diff --git a/optimizationProblem.m b/optimizationProblem.m index 4194f0a..cd5876d 100644 --- a/optimizationProblem.m +++ b/optimizationProblem.m @@ -36,7 +36,7 @@ MaxwellIndex = 100; MaxFunEval = 1e5; MaxIter = 5e3; - motionCompensation = struct('order', [], 'maxMagnitude', [], 'linear', true) + motionCompensation = struct('order', [], 'maxMagnitude', [], 'linear', []) end properties (SetAccess = private) @@ -98,6 +98,22 @@ obj.signs = zeros(obj.N - 1,1)+eps; % setting to zero flips out due to sqrt(0)=complex (??) end + if ~isempty(obj.motionCompensation.order) + + if isempty(obj.motionCompensation.linear) + % Infer empty motionCompensation.linear from values of + % motionCompensation.maxMagnitude + obj.motionCompensation.linear = (obj.motionCompensation.maxMagnitude <= 0); + elseif length(obj.motionCompensation.linear) == 1 + % Convert scalar motionCompensation.linear to a vector + % with the same value. + obj.motionCompensation.linear = obj.motionCompensation.linear * ones(size(obj.motionCompensation.order)); + end + + if length(obj.motionCompensation.linear) ~= length(obj.motionCompensation.order) + error('motionCompensation.linear should be a vector of the same size as motionCompensation.order.') + end + end end end diff --git a/private/nonlconAnalytic.m b/private/nonlconAnalytic.m index a0ce2a0..1b65202 100644 --- a/private/nonlconAnalytic.m +++ b/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 = []; diff --git a/private/optimize.m b/private/optimize.m index 8cdb408..53c7348 100644 --- a/private/optimize.m +++ b/private/optimize.m @@ -163,7 +163,7 @@ [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; @@ -171,12 +171,11 @@ 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 diff --git a/scripted_NOW_Example.m b/scripted_NOW_Example.m index 0a250be..e5cd7b7 100644 --- a/scripted_NOW_Example.m +++ b/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. From a4cad27c1a36759999d5e4229c24f3cb83ce14cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jens=20Sj=C3=B6lund?= Date: Mon, 26 Oct 2020 10:55:48 +0100 Subject: [PATCH 2/2] Removed user option to manually specify linear or nonlinear motion-compensation. --- NOW_GUI.m | 5 +++-- optimizationProblem.m | 30 +++++++++++++----------------- 2 files changed, 16 insertions(+), 19 deletions(-) diff --git a/NOW_GUI.m b/NOW_GUI.m index 9ada5ba..20e5a07 100644 --- a/NOW_GUI.m +++ b/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, true]; +handles.problem.motionCompensation.maxMagnitude = zeros(size(order)); handles.problem = optimizationProblem(handles.problem); set(handles.motionCompensationOrder0, 'Value', false); set(handles.motionCompensationOrder1, 'Value', false); diff --git a/optimizationProblem.m b/optimizationProblem.m index cd5876d..f85085e 100644 --- a/optimizationProblem.m +++ b/optimizationProblem.m @@ -41,7 +41,7 @@ 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,24 +99,19 @@ obj.signs = zeros(obj.N - 1,1)+eps; % setting to zero flips out due to sqrt(0)=complex (??) end - if ~isempty(obj.motionCompensation.order) - - if isempty(obj.motionCompensation.linear) - % Infer empty motionCompensation.linear from values of - % motionCompensation.maxMagnitude - obj.motionCompensation.linear = (obj.motionCompensation.maxMagnitude <= 0); - elseif length(obj.motionCompensation.linear) == 1 - % Convert scalar motionCompensation.linear to a vector - % with the same value. - obj.motionCompensation.linear = obj.motionCompensation.linear * ones(size(obj.motionCompensation.order)); - end - - if length(obj.motionCompensation.linear) ~= length(obj.motionCompensation.order) - error('motionCompensation.linear should be a vector of the same size as motionCompensation.order.') - 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