Optimization method based on the Nonlinear least squares as well as fmincon for the defined objective function (Non-linear function).
Show older comments
Hi there,
I am currently working on a Matlab code for my task (optimization task), and I tried two different methods;
Method A: the fmincon was used to minimize the objective function in which different algorithms, such as interior-point, sqp, sqp-legacy, interior-point, and active-set, were applied.
Method B: Nonlinear least squares optimization was combined with three algorithms, levenberg-marquardt, and trust-region-reflective.
The best results I got now are (from the Nonlinear least squares optimization + 'interior-point'):
RMSE for X, Y, Z Direction(mm) : [1.18, 0.16, 1.33] mm; however, I want to reduce it as far as it is possible.
@ A quick guide for the code, and how it works.
1. I do have 3D data sets, and I used some mathematical methods to convert these 3D data into 2D data (Step 3).
2. The goal is to use some inverse method to get the same 3D datasets; therefore, after defining the objective function, I applied different optimization methods to address this task.
Any assistance or alternative approaches would be greatly appreciated.
Thank you in advance!
% @ Payam Samadi (2025.3.15) NLSP (Nonlinear least squares optimization)
% In step 3, the 2D projection data was extracted from the 3D data at
% different projections.
% In step4, the objective function is defined.
% In step 5, the Nonlinear least squares optimization is used to minimize
% the objective function.
% Algorithm:
% 1. 'levenberg-marquardt' : [1.21, 0.17, 1.39] mm,
% 2. 'trust-region-reflective' : [1.21, 0.17, 1.39] mm,
% 3. 'interior-point' : [1.18, 0.16, 1.33] mm
% Note: It is important to carefully consider the initial guesses for
% target locations, as well as the lower bounds (lb) and upper bounds (ub).
tic;
clc;
clear;
close all;
%% Step 1: Load xls file (3D Tumor position)
Load_Data = xlsread('Sample 1.xlsx'); % Data include 12 (sec)
% Load_Data = xlsread('Sample 2.xlsx'); % Data include 60 (sec)
Time = Load_Data(:,1); % Time (s)
xt=Load_Data(:,2); % xt for targets
yt=Load_Data(:,3); % yt for targets
zt=Load_Data(:,4); % zt for targets
true_T=[xt,yt,zt]; % target locations
[~,numT]=size(true_T); % number of targets
%% Step 2: Define imaging system parameters
SAD = 100; % source-axis distance (cm)
SID = 150; % source-image plane distance (cm)
num_projections = length(Time); % number of different views
% Generate projection angles
alpha = linspace(0, 360, num_projections);
% alpha=30;
theta_rad = deg2rad(alpha); % view angles (radians)
%% Step 3: Projection Model: Compute 2D projections
Xp=zeros(1,num_projections); % allocate array
Yp=zeros(1,num_projections);
for i=1:num_projections
f_theta = (SAD - (true_T(i,1) .* sin(theta_rad(i)) + true_T(i,3) .* cos(theta_rad(i)))) ./ SID;
Xp(1, i) = (true_T(i,1) .* cos(theta_rad(i)) - true_T(i,3) .* sin(theta_rad(i))) ./ f_theta;
Yp(1, i) = true_T(i,2) ./ f_theta;
end
xp = Xp';
yp = Yp';
%% Step 4: Define the objective function
% Define the objective function
objFun = @(T) computeProjectionError(T, xp, yp, theta_rad, SAD, SID, num_projections);
%% Step 5: Minimize the objective function using lsqnonlin
% Algorithm:
% 1. 'levenberg-marquardt' : [1.21, 0.17, 1.39] mm,
% 2. 'trust-region-reflective' : [1.21, 0.17, 1.39] mm,
% 3. 'interior-point' : [1.18, 0.16, 1.33] mm
% Define bounds
lb = []; % Lower bounds
ub = []; % Upper bounds
% Define optimization options
options = optimoptions('lsqnonlin', ...
'Algorithm', 'interior-point', ...
'Display', 'iter', ...
'MaxIterations', 300, ...
'MaxFunctionEvaluations', 20000, ...
'TolFun', 1e-10, ...
'TolX', 1e-10, ...
'StepTolerance', 1e-8, ...
'FiniteDifferenceStepSize', 1e-3, ...
'UseParallel', true);
% Generate a better initial guess
T0 = ones(3, num_projections);
% Run optimization for each projection
estimated_T = zeros(3, num_projections);
for i = 1:num_projections
objFun_i = @(T) computeProjectionError(T, xp(i), yp(i), theta_rad(i), SAD, SID, 1);
estimated_T(:, i) = lsqnonlin(objFun_i, T0(:, i), lb, ub, options);
end
%% Step 6: Calculate RMSE value
% Estimated Tumor Position in X, Y, and Z Direction
Estimate_3D_X = estimated_T(1,:);
Estimate_3D_Y = estimated_T(2,:);
Estimate_3D_Z = estimated_T(3,:);
RMSE_x = Calculate_RMSR (xt,Estimate_3D_X);
RMSE_y = Calculate_RMSR (yt,Estimate_3D_Y);
RMSE_z = Calculate_RMSR (zt,Estimate_3D_Z);
fprintf(['RMSE for X, Y, Z Direction(mm): ' ...
'[%.2f, %.2f, %.2f] mm\n'], RMSE_x, RMSE_y, RMSE_z);
%% Step 7: Plot Estimated vs Real Data
Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z);
toc;
%% Objective function: Computes error between observed and estimated projections
function error = computeProjectionError(T, xp, yp, theta_rad, SAD, SID, num_projections)
for i=1:num_projections
x_est = T(1,:);
y_est = T(2,:);
z_est = T(3,:);
% Compute estimated projection using the perspective transformation
f_theta = (SAD - (x_est(i) .* sin(theta_rad(i)) + z_est(i) .* cos(theta_rad(i))))./ SID;
xp_est(i) = (x_est(i) .* cos(theta_rad(i)) - z_est(i) .* sin(theta_rad(i))) ./ f_theta;
yp_est(i) = y_est(i) ./ f_theta;
end
%% Compute residual error (Method 1)
% Compute the residual error
error_x = xp - xp_est;
error_y = yp - yp_est;
% Return residuals for least squares minimization
error = [error_x; error_y];
end
function RMSE = Calculate_RMSR (True_value,Estimated_value)
% Calculate RMSE
e1 = (True_value(:)-Estimated_value(:)).^2;
RMSE = sqrt(mean(e1,'omitnan'));
end
function Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z)
figure('WindowState', 'maximized')
subplot(3,1,1);
plot(Time, Estimate_3D_X,'r');
hold on;
plot(Time, xt,'b');
legend ('Estimate', 'Real');
title ('X Dierction')
xlabel ('Time (s)');
ylabel ('Motion (mm)')
hold off;
subplot(3,1,2);
plot(Time, Estimate_3D_Y,'r');
hold on;
plot(Time, yt,'b');
legend ('Estimate', 'Real');
title ('Y Dierction')
xlabel ('Time (s)');
ylabel ('Motion (mm)')
hold off;
subplot(3,1,3);
plot(Time, Estimate_3D_Z,'r');
hold on;
plot(Time, zt,'b');
legend ('Estimate', 'Real');
title ('Z Dierction')
xlabel ('Time (s)');
ylabel ('Motion (mm)')
hold off;
end
3 Comments
As far as I understand, in the loop
% Run optimization for each projection
estimated_T = zeros(3, num_projections);
for i = 1:num_projections
objFun_i = @(T) computeProjectionError(T, xp(i), yp(i), theta_rad(i), SAD, SID, 1);
estimated_T(:, i) = lsqnonlin(objFun_i, T0(:, i), lb, ub, options);
end
you solve two (after rewriting: linear) equations in three unknowns for each value of i.
Is that really what you want ? Usually, there will be infinitely many solutions for each value of i.
Matt J
on 17 Mar 2025
Since you have no constraints, you may as well use fminunc instead of fmincon, or fsolve instead of lsqnonlin.
payam samadi
on 18 Mar 2025
Accepted Answer
The goal is to use some inverse method to get the same 3D datasets; therefore, after defining the objective function, I applied different optimization methods to address this task.
But for every time Time(i), you appear to have only one projection view of the tumor position. This is not enough, I'm afraid. As noted by Torsten as well, it will be an underdetermined problem. However, since this appears to be a radiotherapy scenario, and since radiotherapy treatment systems typically have two x-ray imagers, you should be able to get two projections per tumor position. If we adapt your code to this scenario (see below), we can see the inverse problem is easily solved:
%% Step 1: Load xls file (3D Tumor position)
Load_Data = xlsread('Sample 1.xlsx'); % Data include 12 (sec)
% Load_Data = xlsread('Sample 2.xlsx'); % Data include 60 (sec)
Time = Load_Data(:,1); % Time (s)
xt=Load_Data(:,2); % xt for targets
yt=Load_Data(:,3); % yt for targets
zt=Load_Data(:,4); % zt for targets
true_T=[xt,yt,zt]; % target locations
numT=height(true_T); % number of targets
%% Step 2: Define imaging system parameters
SAD = 100; % source-axis distance (cm)
SID = 150; % source-image plane distance (cm)
num_projections = length(Time); % number of different views
% Generate projection angles
alpha = linspace(0, 360, num_projections)';
theta_rad = deg2rad(alpha)+[0,pi/2];
num_detectors=width(theta_rad);
%% Step 3: Projection Model: Compute 2D projections
[Xp,Yp]=deal(nan(num_projections,num_detectors)); % allocate arrays
for i=1:numT
for j=1:num_detectors
f_theta = (SAD - (true_T(i,1) .* sin(theta_rad(i,j)) + true_T(i,3) .* cos(theta_rad(i,j)))) ./ SID;
Xp(i,j) = (true_T(i,1) .* cos(theta_rad(i,j)) - true_T(i,3) .* sin(theta_rad(i,j))) ./ f_theta;
Yp(i,j) = true_T(i,2) ./ f_theta;
end
end
xp=Xp; yp=Yp;
%% Step 4: Define the objective function
tic
T = optimvar('T',[numT,3]);
[eqn1,eqn2]=deal(cell(num_projections,1));
for i=1:numT
for j=1:num_detectors
f_theta = (SAD - (T(i,1) .* sin(theta_rad(i,j)) + T(i,3) .* cos(theta_rad(i,j)))) ./ SID;
eqn1{i,j}=xp(i,j).*f_theta == (T(i,1) .* cos(theta_rad(i,j)) - T(i,3) .* sin(theta_rad(i,j))) ;
eqn2{i,j}=yp(i,j).*f_theta == T(i,2) ;
end
end
prob=eqnproblem();
prob.Equations.eqn1=vertcat(eqn1{:});
prob.Equations.eqn2=vertcat(eqn2{:});
s=prob2struct(prob);
estimated_T=reshape(s.C\s.d,[],3);
toc
%% Step 6: Calculate RMSE value
% Estimated Tumor Position in X, Y, and Z Direction
Estimate_3D_X = estimated_T(:,1);
Estimate_3D_Y = estimated_T(:,2);
Estimate_3D_Z = estimated_T(:,3);
RMSE_x = Calculate_RMSR (xt,Estimate_3D_X);
RMSE_y = Calculate_RMSR (yt,Estimate_3D_Y);
RMSE_z = Calculate_RMSR (zt,Estimate_3D_Z);
fprintf(['RMSE for X, Y, Z Direction(mm): ' ...
'[%.2f, %.2f, %.2f] mm\n'], RMSE_x, RMSE_y, RMSE_z);
%% Step 7: Plot Estimated vs Real Data
Plot_results_Matt(Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z);
34 Comments
Matt J
on 18 Mar 2025
Note, the approach above reformulates the problem as a set of linear equations, as we have been discussing, and estimates T non-iteratively. This may not give you the most optimal triangulation, as the stochastic errors in xp,yp will be weighted in a weird way. However, it will give you a better initial guess for an iterative optimization, if you wish to go back and refine it.
payam samadi
on 18 Mar 2025
Your link doesn't lead anywhere, so I don't know what Poulsen did. However, it has already been pointed out to you that it is absolutely not possible with only 1 projection and no other constraints. The problem is underdetermined in that case.
Perhaps Poulsen et al. apply special constraints or structure on the problem, but you have not applied any special tricks in your code. You are using a very basic unconstrained model. So you won't be able to get much from it.
payam samadi
on 18 Mar 2025
payam samadi
on 18 Mar 2025
Torsten
on 18 Mar 2025
I think the people with a tumor wouldn't be amused to hear that you localized it with fake projections.
payam samadi
on 18 Mar 2025
If we generate a fake projection from the original one (which I have already done), could this serve as a potential solution for our task? I have attached the results for your review.
No, if the synthetic ("fake") projection is derived from the original projection, it cannot contain any new information. I get very poor results when I run your Method_1.m. Your Results.docx attachment contains no results, only a description of your idea. So, I'm not sure what we were supposed to see.
how can we use the current code (the optimized one with the initial datasets; ex 12-second data)
If you mean, how do you use the original, nonlinear equation model, you can adapt your main loops as below,
T0=____; % Nx3 matrix of initial guesses. Get this from the linear method
estimated_T=nan(numT,3);
options = optimoptions('fsolve', 'Algorithm', 'levenberg-marquardt', 'Display', 'off');
for i=1:numT
i,
T = optimvar('T',[1,3]);
[eqn1,eqn2]=deal(cell(2,1));
for j=1:num_detectors
f_theta = (SAD - (T(1) .* sin(theta_rad(i,j)) + T(3) .* cos(theta_rad(i,j)))) ./ SID;
eqn1{j}=xp(i,j) == (T(1) .* cos(theta_rad(i,j)) - T(3) .* sin(theta_rad(i,j)))./f_theta ;
eqn2{j}=yp(i,j) == T(2)./f_theta ;
end
prob=eqnproblem();
prob.Equations.eqn1=vertcat(eqn1{:});
prob.Equations.eqn2=vertcat(eqn2{:});
sol0.T=T0(i,:);
estimated_T(i,:) = solve(prob, sol0, 'Options', options).T;
end
@Torsten If there is no measurement noise in xp,yp then either modeling approach will give the same result. However, since there is generally additve stochastic noise in the measurements, cross-multiplying by f_theta creates a non-additive noise-model, and that is usually not ideal.
Although, as I've said, the linear solution is still useful as an initializer for the nonlinear solution.
payam samadi
on 19 Mar 2025
Edited: payam samadi
on 21 Mar 2025
payam samadi
on 19 Mar 2025
Matt J
on 19 Mar 2025
Your Xp, Yp aren't dervied from real measurements. You are still simulating them from a set of known T values with no errors. If you add errors to your simulation, like below, you will see that you don't get a perfect re-estimate of T.
sig=0.3; %simulated error standard deviation
for i=1:numT
for j=1:num_detectors
f_theta = (SAD - (true_T(i,1) .* sin(theta_rad(i,j)) + true_T(i,3) .* cos(theta_rad(i,j)))) ./ SID;
Xp(i,j) = (true_T(i,1) .* cos(theta_rad(i,j)) - true_T(i,3) .* sin(theta_rad(i,j))) ./ f_theta;
Yp(i,j) = true_T(i,2) ./ f_theta;
Xp(i,j)=Xp(i,j) + sig*randn;
Yp(i,j)=Yp(i,j) + sig*randn;
end
end
Torsten
on 19 Mar 2025
The equations you define and solve for T just recover exactly the true_T values you used in the loop before. This means that you get RMSE values of 0.
payam samadi
on 19 Mar 2025
Even if so, they are working with a different model for the tumor position, one that may have fewer degrees of freedom than the one in your code. In your code, you give the tumor position a full 3 degrees of freedom (x,y,z) per time point. It is easy to see that this is under-determined when you have only 1 projection.
payam samadi
on 19 Mar 2025
Moved: Matt J
on 19 Mar 2025
payam samadi
on 19 Mar 2025
Edited: payam samadi
on 19 Mar 2025
William and I are both working on it together.
If you give the link to that thread, Torsten and/or I and/or others might join it to see what we can contribute. However, as for the fmincon error message, it looks pretty simple. It is telling you that evaluating your objective function at x0 is returning nan or inf. So, you need to finish debugging your objective function before trying to minimize it. Usually you would investigate the error by calling the objective function on x0 separately with the debugger active and see what is going on.
payam samadi
on 20 Mar 2025
Dear Matt,
Thanks for your clue. I finished debugging the code. It gave me some initial results, but it was too far from my expectations (Please see the difference between the covariance matrix from the Real data and Estimation).
I attached the code here and I think there might be two possibilities for it.
1. The code may be getting trapped in local minima (Please take a moment to review the notes at the beginning of the code). This could be due to an incorrect initial guess for LvecInit in the computeError function, improper selection of the lower (LB) and upper (UB) bounds, or the use of an unsuitable optimization method for this task. Also, three methods to optimize the objective function, including fmincon (results in no convergence since the matrix B (inverse covariance matrix) is not positive definite), elder-Mead (simple) as well as Nelder-Mead (with constraint) [both run and get some results; but it so far from my expectations].
2. The code functions correctly, but step 3 (Projection Model: Compute 2D projections) employs an incorrect method for projecting the 3D data into a 2D projection.
I also attached my notes from the article that are related to this work (I used this information during implementation), as I thought they would be helpful.
Please let me know if you come up with any ideas.
Thanks in advance for your time.
% Payam Samadi and William Rose, 2025-03-19.
% After loading the data and defining the system parameters in
% steps 1 and 2, the 2D projections are estimated from the original datasets.
% A simulated error standard deviation (sigma) is then added to
% the 2D projection datasets (Step 3).
% Step 4: Estimate distribution parameters (Appendix of Poulsen et al.,
% 2009).
% The code is not working; the problem comes from FMLE (Nan or INF).
%%
% New Update (2025.3.21): The lb and ub parameters have been replaced
% with LvecInit, which serves as the initial guess for Lvec.
% The code has been executed accordingly.
% Step 5 has also been included to evaluate the results
% from the estimated Final_Vec_opt against the real data.
% Three methods were utilized to optimize the objective function (Objective.m):
% 1. Method 1: fmincon
% 2. Method 2: Simple Nelder-Mead
% 3. Method 3: Nelder-Mead with constraints
% Two possibilities have been identified for this part:
% 1. The code may be getting trapped in local minima [The LvecInit
% (initial guess for Lvec) in computer error is wrong, or
% the lower bounds (LB), as well as upper bounds (UP), should be modified].
% 2. The code functions correctly, but step 2 employs
% an incorrect method for projecting the 3D data into a 2D projection.
%%
tic
clc
clear
close all
%% Step 1: Load xls file (3D Tumor position)
Load_Data = xlsread('Sample 1.xlsx'); % Data include 12 (sec)
Time = Load_Data(:,1); % Time (s)
xt=Load_Data(:,2); % xt for targets
yt=Load_Data(:,3); % yt for targets
zt=Load_Data(:,4); % zt for targets
true_T=[xt,yt,zt]; % target locations
numT=height(true_T); % number of targets
%% Step 2: Define imaging system parameters
SAD = 100; % source-axis distance (cm)
SID = 150; % source-image plane distance (cm)
num_projections = length(Time); % number of different views
% Generate projection angles
alpha = linspace(0, 360, num_projections)';
theta_rad = deg2rad(alpha);
num_detectors= 1;
%% Step 3: Projection Model: Compute 2D projections
% Add simulated error standard deviation (sig) to the projection
sig=0.05; %simulated error standard deviation
[Xp,Yp]=deal(nan(num_projections,num_detectors)); % allocate arrays
for i=1:numT
for j=1:num_detectors
f_theta = (SAD - (true_T(i,1) .* sin(theta_rad(i,j)) + true_T(i,3) .* cos(theta_rad(i,j)))) ./ SID;
Xp(i,j) = (true_T(i,1) .* cos(theta_rad(i,j)) - true_T(i,3) .* sin(theta_rad(i,j))) ./ f_theta;
Yp(i,j) = true_T(i,2) ./ f_theta;
Xp(i,j)=Xp(i,j) + sig*randn;
Yp(i,j)=Yp(i,j) + sig*randn;
end
end
%% Step 4: Estimate distribution parameters (Appendix of Poulsen et al., 2009).
% Matrix LvecArray is nine unknowns vector.
% My unknowns are 3 means, 3 standard deviations, and 3 correlations:
% [mx, my, mz, sx, sy, sz, rxy, ryz, rzx].
Final_Vec_opt = computeError(Xp, Yp, SAD, SID, alpha);
%% Step 5: Evalute the Results from the Estimated Final_Vec_opt and Real Data
A_matrix = covariance_matrix(xt,yt,zt);
Estimated_A_matrix = A_matrix_From_Estimation (Final_Vec_opt);
%% Step 6: Estimate 3D coordinates from projection
[Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z] = Estimate_3D_coordinates_from_Final_Vec_opt (Xp, Yp, ...
SID, SAD, alpha, Final_Vec_opt, num_projections);
%% Step 7: Calculate RMSE value in each direction
RMSE_x = Calculate_RMSR (xt,Estimate_3D_X);
RMSE_y = Calculate_RMSR (yt,Estimate_3D_Y);
RMSE_z = Calculate_RMSR (zt,Estimate_3D_Z);
fprintf(['RMSE for X, Y, Z Direction(mm): ' ...
'[%.2f, %.2f, %.2f] mm\n'], RMSE_x, RMSE_y, RMSE_z);
%% Step 7: Plot Estimated vs Real Data
Plot_results(Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z);
toc
function Estimated_A_matrix = A_matrix_From_Estimation (Final_Vec_opt)
% This function calculate the A_matrix (covariance matrix) from the
% estimation value of the Final_Vec_opt.
% Given Final_Vec_opt = [mx, my, mz, sx, sy, sz, rxy, ryz, rzx]
mx = Final_Vec_opt(1); % Mean value of the X direction
my = Final_Vec_opt(2); % Mean value of the Y direction
mz = Final_Vec_opt(3); % Mean value of the Z direction
sx = Final_Vec_opt(4); % Standard deviation of X
sy = Final_Vec_opt(5); % Standard deviation of Y
sz = Final_Vec_opt(6); % Standard deviation of Z
rxy = Final_Vec_opt(7); % Correlation between X and Y
ryz = Final_Vec_opt(8); % Correlation between Y and Z
rzx = Final_Vec_opt(9); % Correlation between Z and X
% Compute covariance values
Cov_x_y = rxy * sx * sy; % Covariance between X and Y Direction
Cov_y_z = ryz * sy * sz; % Covariance between Y and Z Direction
Cov_z_x = rzx * sz * sx; % Covariance between X and Z Direction
% Construct covariance matrix A
Estimated_A_matrix = [sx^2, Cov_x_y, Cov_z_x;
Cov_x_y, sy^2, Cov_y_z;
Cov_z_x, Cov_y_z, sz^2];
% Display matrix
disp('Covariance Matrix A_matrix From Estimation:');
disp(Estimated_A_matrix);
end
function RMSE = Calculate_RMSR (True_value,Estimated_value)
% Calculate RMSE
e1 = (True_value(:)-Estimated_value(:)).^2;
RMSE = sqrt(mean(e1,'omitnan'));
end
function LvecArray = computeError(Xp, Yp, SAD, SID, alpha)
% "computeError" Estimate unknown parameters for markers, from projection data.
% For one marker, estimate mean position and standard deviations and
% correlations of movement, using method of Poulsen et al., 2009.
% Input data is projections of 3D marker positions onto 2D image planes.
% First scenario: There is one marker (N = 1) and Npos positions of the image plane.
% The image plane rotates about the y-axis. See Fig.9 in Poulsen et al., 2009.
% Inputs
% Xp=marker X-coordinates on image plane, Npos x N (cm)
% Yp=marker Y-coordinates on image plane, Npos x N (cm)
% SAD=source-axis distance (cm)
% SID=source-image plane distance (cm)
% alpha=projection angles, 1 x Npos (radians)
% Output
% LvecArray=estimated marker parameters.
% LvecArray=[mux,muy,muz,sdx,sdy,sdz,pxy,ryz,rzx]
% where mux,y,z=estim.mean position (cm), sdx,y,z=estim. st.dev. (cm), and
% rxy,etc=estim.correlation between instantaneous x,y position, etc.
% Function(s) called
% Objective(x, Xp, Yp, SAD, SID, alpha) returns F_MLE.
[Npos,N]=size(Xp); % determine number of targets and projections from Xp
% error checking
if size(Xp)~=size(Yp)
error('Error: size(Yp) must = size(Xp).');
end
if length(alpha)~=Npos
error('Error: length(alpha) must = number of rows in Xp.')
end
LvecArray=zeros(N,9); % allocate array for results
% LvecINit and lb and ub are the same on every loop pass.
% lb, ub are bounds for [mux,muy,muz,sdx,sdy,sdz,pxy,ryz,rzx].
LvecInit=[0,0,0,1,1,1,0,0,0]; % initial guess for Lvec
lb=[-10 -10 -10 -10 -10 -10 -10 -10 -10]; % lower bounds
ub=[10 10 10 10 10 10 10 10 10]; % upper bounds
% For Methods 2 and 3 (Nelder-Mead optimization [simple and with constraint])
% Set optimization options to improve convergence
options = optimset('fminsearch');
options.MaxIter = 2000; % Increase max iterations
options.MaxFunEvals = 4000; % Increase function evaluations
options.TolFun = 1e-10; % Reduce function tolerance
options.TolX = 1e-10; % Reduce solution tolerance
options.Display = 'iter'; % Show iterations
for i=1:N
% Find 1x9 vector Lvec that minimizes F_MLE for marker i.
xpv = Xp(:,i); % column vector, length Np
ypv = Yp(:,i); % column vector, length Np
% fun=@(x) Objective(x, xpv, ypv, SAD, SID, alpha); % function to
% minimize (Not included matrix B is positive definite).
fun=@(x) Objective_1(x, xpv, ypv, SAD, SID, alpha); % function to
% minimize (Included matrix B is positive definite).
%% Methods to Optimize Objective Function
% Lvec = fmincon(fun,LvecInit,[],[],[],[],lb,ub); % Method 1: Use the fmincon to calculate the Lvec
% Lvec = fminsearch(fun, LvecInit, options); % Method 2: Perform Nelder-Mead (simple) optimization using fminsearch
% Method 3: Perform Nelder-Mead (with constraint) optimization using
% fminsearch + add constraint to store best solution found.
% Multi-start strategy: Run multiple times with different initial guesses
num_restarts = 5; % Number of restarts
best_Lvec = [];
best_fval = Inf;
for j = 1:num_restarts
% Run Nelder-Mead optimization
[Lvec, fval] = fminsearch(fun, LvecInit, options);
% Store best solution found
if fval < best_fval
best_fval = fval;
best_Lvec = Lvec;
end
end
LvecArray(i,:) = Lvec; % save Lvec in row i of LvecArray
end
end
function A_matrix = covariance_matrix(xt,yt,zt)
% This function calculate the A_matrix (covariance matrix) from the real
% motion data.
% Calculate the Variance
Var_x = var(xt); % X Direction
Var_y = var(yt); % Y Direction
Var_z = var(zt); % Z Direction
% Calculate the Covariance
Cov_x_y = cov(xt,yt);
Cov_x_y = Cov_x_y(1,2); % Covariance between X and Y Direction
Cov_x_z = cov(xt,zt);
Cov_x_z = Cov_x_z(1,2); % Covariance between X and Z Direction
Cov_y_z = cov(yt,zt);
Cov_y_z = Cov_y_z(1,2); % Covariance between Y and Z Direction
% Construct covariance matrix (A_matrix)
A_matrix = [Var_x, Cov_x_y, Cov_x_z;
Cov_x_y, Var_y, Cov_y_z;
Cov_x_z, Cov_y_z, Var_z];
% Display matrix
disp('Covariance Matrix A_matrix From Real Data:');
disp(A_matrix);
end
function [Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z] = Estimate_3D_coordinates_from_Final_Vec_opt (xp, yp, ...
SDD, SAD, alpha, Final_Vec_opt, num_projections)
% @ Payam Samadi (2025.3.21). Estimate_3D_coordinates_from_Final_Vec_opt is
% used to estimate 3D coordinates based on the Poulsen et al., 2009 (See page 4035).
% The Final_Vec_opt is [mux,muy,muz] for eahc marker (i).
% Np refers to the number of projections.
Number_of_target=length(xp); % number of targets
Estimate_3D_X = zeros(num_projections,1);
Estimate_3D_Y = zeros(num_projections,1);
Estimate_3D_Z = zeros(num_projections,1);
for i=1:num_projections
p=[xp(i)'.*cos(alpha(i))-(SDD-SAD)*sin(alpha(i));...
yp(i)';...
-xp(i)'.*sin(alpha(i))-(SDD-SAD)*cos(alpha(i))]; % Projection Point (Eq. A.4)
f=SAD.*[sin(alpha(i)); 0; cos(alpha(i))]; % Focal point (Eq. A.5)
e_hat=(f-p)/norm(f-p); % Unit vector "e_hat" (Eq. A.6)
% Estimate 3D Position (Eq. A-7)
Estimate_3D_X (i) = p(1) + e_hat(1) * Final_Vec_opt (1);
Estimate_3D_Y (i) = p(2) + e_hat(2) * Final_Vec_opt (2);
Estimate_3D_Z (i) = p(3) + e_hat(3) * Final_Vec_opt (3);
end
end
function FMLE = Objective(x, xp, yp, SAD, SID, alpha)
% "Objective" Compute FMLE for observations of one marker.
% Compute FMLE=neg.log likelihood of the observed marker positions, for
% a marker whose observed positions xp,yp are random samples from the
% multivariate normal distribution specified by the values in x.
% The marker is in 3D space. It is observed Np times, each time
% projected onto a 2D image plane. The x-ray source and image plane
% rotate about the y axis (Poulsen Fig.9).
% See Poulsen et al., 2009, especially Appendix.
% Inputs
% x=[mx,my,mz,sdx,sdy,sdz,rxy,ryz,rzx]=1x9 vector with parameters of
% the marker's statistical distribution. The means and SDs are in
% cm; the correlations rxy,... are dimensionless.
% xp=observed marker x-values (Np x 1) (cm)
% yp=observed marker y-values (Np x 1) (cm)
% SAD=source-axis distance (cm)
% SID=source-image plane distance (cm)
% alpha=orientation of source and image plane (1 x Np) (radian)
% Output
% FMLE=sum(fMLE)
mx=x(1); my=x(2); mz=x(3); % extract values from x
sdx=x(4); sdy=x(5); sdz=x(6);
rxy=x(7); ryz=x(8); rzx=x(9);
%% Calculate FMLE
r0=[mx;my;mz]; % theoretical mean position
Np=length(xp); % number of projections
A=[sdx^2, rxy*sdx*sdy, rzx*sdz*sdx;... % theo. covariance matrix
rxy*sdx*sdy, sdy^2, ryz*sdy*sdz;...
rzx*sdz*sdx, ryz*sdy*sdz, sdz^2];
B=inv(A); % inverse covariance matrix
K=zeros(Np,1); % allocate K
sigma=zeros(Np,1); % allocate sigma
for i=1:Np
p=[xp(i)*cos(alpha(i))-(SID-SAD)*sin(alpha(i));... % p=3D coords of(Xp,Yp)
yp(i);...
-xp(i)*sin(alpha(i))-(SID-SAD)*cos(alpha(i))]; % A.4
f=SAD*[sin(alpha(i)); 0; cos(alpha(i))]; % A.5 focal point
ehat=(f-p)/norm(f-p); % A.6 unit vector
sigma(i)=1/sqrt(ehat'*B*ehat); % A.8
mu=((r0-p)'*B*ehat)/(ehat'*B*ehat); % A.9
K(i)=((sqrt(det(B))/((2*pi)^(3/2))))*...
exp(((-1/2)*(p+ehat*mu-r0)'*B*(p+ehat*mu-r0))); % A.10
end
fMLE=-log(sqrt(2*pi)*K.*sigma); % A.12
FMLE=sum(fMLE); % A.12
end
function FMLE = Objective_1(x, xp, yp, SAD, SID, alpha)
% "Objective_1" Compute FMLE for observations of one marker.
% Compute FMLE=neg.log likelihood of the observed marker positions, for
% a marker whose observed positions xp,yp are random samples from the
% multivariate normal distribution specified by the values in x.
% The marker is in 3D space. It is observed Np times, each time
% projected onto a 2D image plane. The x-ray source and image plane
% rotate about the y axis (Poulsen Fig.9).
% See Poulsen et al., 2009, especially Appendix.
% Inputs
% x=[mx,my,mz,sdx,sdy,sdz,rxy,ryz,rzx]=1x9 vector with parameters of
% the marker's statistical distribution. The means and SDs are in
% cm; the correlations rxy,... are dimensionless.
% xp=observed marker x-values (Np x 1) (cm)
% yp=observed marker y-values (Np x 1) (cm)
% SAD=source-axis distance (cm)
% SID=source-image plane distance (cm)
% alpha=orientation of source and image plane (1 x Np) (radian)
% Output
% FMLE=sum(fMLE)
% Ensures that matrix B (inverse covariance matrix) is positive definite.
mx=x(1); my=x(2); mz=x(3); % extract values from x
sdx=x(4); sdy=x(5); sdz=x(6);
rxy=x(7); ryz=x(8); rzx=x(9);
%% Calculate FMLE
r0=[mx;my;mz]; % theoretical mean position
Np=length(xp); % number of projections
% Construct the covariance matrix A
A=[sdx^2, rxy*sdx*sdy, rzx*sdz*sdx;...
rxy*sdx*sdy, sdy^2, ryz*sdy*sdz;...
rzx*sdz*sdx, ryz*sdy*sdz, sdz^2];
% Ensure positive definiteness using regularization
epsilon = 1e-10; % Small positive value
A = A + epsilon * eye(3); % Adds a small value to the diagonal to ensure positive definiteness
% Compute inverse covariance matrix B
B=inv(A);
% Verify B is positive definite
if any(eig(B) <= 0)
error('Matrix B is not positive definite. Adjusting A may be needed.');
end
K=zeros(Np,1); % allocate K
sigma=zeros(Np,1); % allocate sigma
for i=1:Np
p=[xp(i)*cos(alpha(i))-(SID-SAD)*sin(alpha(i));... % p=3D coords of (Xp,Yp)
yp(i);...
-xp(i)*sin(alpha(i))-(SID-SAD)*cos(alpha(i))]; % A.4
f=SAD*[sin(alpha(i)); 0; cos(alpha(i))]; % A.5 focal point
ehat=(f-p)/norm(f-p); % A.6 unit vector
sigma(i)=1/sqrt(ehat'*B*ehat); % A.8
mu=((r0-p)'*B*ehat)/(ehat'*B*ehat); % A.9
K(i)=((sqrt(det(B))/((2*pi)^(3/2))))*...
exp(((-1/2)*(p+ehat*mu-r0)'*B*(p+ehat*mu-r0))); % A.10
end
fMLE=-log(sqrt(2*pi)*K.*sigma); % A.12
FMLE=sum(fMLE); % A.12
end
function Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z)
figure('WindowState', 'maximized')
subplot(3,1,1);
plot(Time, Estimate_3D_X,'r-','LineWidth',2);
hold on;
plot(Time, xt,'bo');
legend ('Estimate', 'Real');
title ('X Direction')
xlabel ('Time (s)');
ylabel ('Motion (mm)')
hold off;
subplot(3,1,2);
plot(Time, Estimate_3D_Y,'r-','LineWidth',2);
hold on;
plot(Time, yt,'bo');
legend ('Estimate', 'Real');
title ('Y Direction')
xlabel ('Time (s)');
ylabel ('Motion (mm)')
hold off;
subplot(3,1,3);
plot(Time, Estimate_3D_Z,'r-','LineWidth',2);
hold on;
plot(Time, zt,'bo');
legend ('Estimate', 'Real');
title ('Z Direction')
xlabel ('Time (s)');
ylabel ('Motion (mm)')
hold off;
end
One thing to test for in regads to local minimum is to run the optimization with a known global solution, with the intial x0 chosen at the known solution and again also with x0 near the known solution. When initialized at the known solution, it should stop right away. When initialized near it, it should quickly move to the solution, without getting stuck.
payam samadi
on 23 Mar 2025
Please let me know if you come up with any ideas.
It is not clear to me whether the advice from my last comment was followed. Do any of these choices of Lvecinit correspond to the known ground truth solution of simulated data, as I advised?
A few random observations though:
(1) fminsearch will usually perform poorly for a problem with 9 unknown variables. There is no theoretical gaurantee that it will converge for more than 1 variable. Empirically, it does okay for about 6, but 9 is pushing it.
(2) Your expression for the covariance matrix requires rxy,ryz,rzx to all satisfy abs( r )<=1, if it is supposed to give you a positive semidefinite matrix. It doesn't appear that you have chosen lb,ub to ensure this. That is in addition to the fact that you are using fminsearch, which does not enforce bounds at all.
payam samadi
on 23 Mar 2025
The condition is that I used very close LvecInit, and as you predicted, it quickly ran and gave me results close to it.
My advice, though, was to choose Lvecinit to be the known global minimum. In other words, I am asking you to choose an LvecTrue and use it to generate hypothetical projection data xp,yp. When you run the optimizer with Lvecinit=LvecTrue, the optimizer should be able to immediately recognize LvecInit as a solution and stop.
The values of LvecInit that you have selected in your first test do not seem to be the global minimum because they produce inf values. You suppressed the infs with sum(fMLE(~isinf(fMLE))), but why would there be any Infs to suppress if LvecInit was the global minimum?
Could you please provide more specific details that I can try?
I have provided specific details. I specifically told you not to use fminsearch. Your code is already set up to use fmincon instead, and that is what you should do. Also, I told you that you need to set your lb, ub to bound to ensure that -1<=[rxy,ryz,rzx]<=1.
Also, remove the line,
fMLE=sum(fMLE(~isinf(fMLE)))
It makes your objective discontinuous, which violates the assumptions of fmincon.
Another possible way to deal with the covariance being non-positive definite is to reparametrize as below. With this method, you do not need to put bounds on x(4:9) or even to construct A at all.
mx=x(1); my=x(2); mz=x(3); %Extract mean values
L=zeros(3);
L([1,2,3,5,6,9]0=x(4:9); %Form a lower triangular matrix
B=L*L'; % Construct inverse covariance matrix directly
I would also recommend some changes to the main loop that computes the FMLE. In particular, I think you should compute negative loglikelihood terms directly,
%% Calculate FMLE
r0=[mx;my;mz]; % theoretical mean position
Np=length(xp); % number of projections
mahalanobisDist=zeros(Np,1); % allocate K
sigma=zeros(Np,1); % allocate sigma
for i=1:Np
p=[xp(i)*cos(alpha(i))-(SID-SAD)*sin(alpha(i));... % p=3D coords of (Xp,Yp)
yp(i);...
-xp(i)*sin(alpha(i))-(SID-SAD)*cos(alpha(i))]; % A.4
f=SAD*[sin(alpha(i)); 0; cos(alpha(i))]; % A.5 focal point
ehat=(f-p)/norm(f-p); % A.6 unit vector
sigma(i)=1/sqrt(ehat'*B*ehat); % A.8
mu=((r0-p)'*B*ehat)/(ehat'*B*ehat); % A.9
mahalanobisDist(i) = (p+ehat*mu-r0)'*B*(p+ehat*mu-r0) ; % A.10
end
FMLE = -Np*log(det(B)) + sum(mahalanobisDist); % A.12
This avoids potential numerical problems with evaluating log(exp(___)).
payam samadi
on 24 Mar 2025
Do you think other power optimization methods such as GA and gravity search could handle and optimize these nine parameters?
It might help with the local minima, but I feel like this problem should be easy to come up with a good Lvecinit. I don't know how you are deriving your Lvecinit currently.
Still, the results are so far from our expectations.
In Objective_2 you have imposed bounds on x(4:9) which I said you should not have. You can impose bounds on x(1:3), but I don't know where you are getting them. In your test data, it looks like the target point drifts more than a centimeter from isocenter, so +/-10 is the least I would choose.
Aside from all that, you still don't appear to have run the experiment initializing at ground truth. I will take a break from this thread until you do.
payam samadi
on 25 Mar 2025
Edited: Walter Roberson
on 29 Mar 2025
More Answers (2)
You can add a multiple of the vector you obtain by the below command "double(null(A))" to "estimated_T", and you will still get a solution to your linear system of two equations. The system is underdetermined and there exist infinitly many solutions for it.
% @ Payam Samadi (2025.3.15) NLSP (Nonlinear least squares optimization)
% In step 3, the 2D projection data was extracted from the 3D data at
% different projections.
% In step4, the objective function is defined.
% In step 5, the Nonlinear least squares optimization is used to minimize
% the objective function.
% Algorithm:
% 1. 'levenberg-marquardt' : [1.21, 0.17, 1.39] mm,
% 2. 'trust-region-reflective' : [1.21, 0.17, 1.39] mm,
% 3. 'interior-point' : [1.18, 0.16, 1.33] mm
% Note: It is important to carefully consider the initial guesses for
% target locations, as well as the lower bounds (lb) and upper bounds (ub).
tic;
clc;
clear;
close all;
%% Step 1: Load xls file (3D Tumor position)
Load_Data = xlsread('Sample 1.xlsx'); % Data include 12 (sec)
% Load_Data = xlsread('Sample 2.xlsx'); % Data include 60 (sec)
Time = Load_Data(:,1); % Time (s)
xt=Load_Data(:,2); % xt for targets
yt=Load_Data(:,3); % yt for targets
zt=Load_Data(:,4); % zt for targets
true_T=[xt,yt,zt]; % target locations
[~,numT]=size(true_T); % number of targets
%% Step 2: Define imaging system parameters
SAD = 100; % source-axis distance (cm)
SID = 150; % source-image plane distance (cm)
num_projections = length(Time); % number of different views
% Generate projection angles
alpha = linspace(0, 360, num_projections);
% alpha=30;
theta_rad = deg2rad(alpha); % view angles (radians)
%% Step 3: Projection Model: Compute 2D projections
Xp=zeros(1,num_projections); % allocate array
Yp=zeros(1,num_projections);
for i=1:num_projections
f_theta = (SAD - (true_T(i,1) .* sin(theta_rad(i)) + true_T(i,3) .* cos(theta_rad(i)))) ./ SID;
Xp(1, i) = (true_T(i,1) .* cos(theta_rad(i)) - true_T(i,3) .* sin(theta_rad(i))) ./ f_theta;
Yp(1, i) = true_T(i,2) ./ f_theta;
end
xp = Xp';
yp = Yp';
% Run optimization for each projection
estimated_T = zeros(3, num_projections);
T = sym('T',[3 1]);
for i = 1:num_projections
f_theta = (SAD - (T(1) .* sin(theta_rad(i)) + T(3) .* cos(theta_rad(i))))./ SID;
eqn1 = xp(i)*f_theta - (T(1) .* cos(theta_rad(i)) - T(3) .* sin(theta_rad(i)))==0;
eqn2 = yp(i)*f_theta - T(2)==0;
[A,b] = equationsToMatrix([eqn1,eqn2]);
double(null(A))
estimated_T(:, i) = double(A)\double(b);
end
%% Step 6: Calculate RMSE value
% Estimated Tumor Position in X, Y, and Z Direction
Estimate_3D_X = estimated_T(1,:);
Estimate_3D_Y = estimated_T(2,:);
Estimate_3D_Z = estimated_T(3,:);
RMSE_x = Calculate_RMSR (xt,Estimate_3D_X);
RMSE_y = Calculate_RMSR (yt,Estimate_3D_Y);
RMSE_z = Calculate_RMSR (zt,Estimate_3D_Z);
fprintf(['RMSE for X, Y, Z Direction(mm): ' ...
'[%.2f, %.2f, %.2f] mm\n'], RMSE_x, RMSE_y, RMSE_z);
%% Step 7: Plot Estimated vs Real Data
Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z);
toc;
function RMSE = Calculate_RMSR (True_value,Estimated_value)
% Calculate RMSE
e1 = (True_value(:)-Estimated_value(:)).^2;
RMSE = sqrt(mean(e1,'omitnan'));
end
function Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z)
figure('WindowState', 'maximized')
subplot(3,1,1);
plot(Time, Estimate_3D_X,'r');
hold on;
plot(Time, xt,'b');
legend ('Estimate', 'Real');
title ('X Dierction')
xlabel ('Time (s)');
ylabel ('Motion (mm)')
hold off;
subplot(3,1,2);
plot(Time, Estimate_3D_Y,'r');
hold on;
plot(Time, yt,'b');
legend ('Estimate', 'Real');
title ('Y Dierction')
xlabel ('Time (s)');
ylabel ('Motion (mm)')
hold off;
subplot(3,1,3);
plot(Time, Estimate_3D_Z,'r');
hold on;
plot(Time, zt,'b');
legend ('Estimate', 'Real');
title ('Z Dierction')
xlabel ('Time (s)');
ylabel ('Motion (mm)')
hold off;
end
6 Comments
Matt J
on 17 Mar 2025
The system is underdetermined and there exist infinitly many solutions for it.
Only because equationsToMatrix is being used here on each projection individually. When you combine the equations from multiple, widely separated projection angles, you get a well-determined system.
When you combine the equations from multiple, widely separated projection angles, you get a well-determined system.
Yes, but the original code also solves for each projection individually. The variable "num_projections" in function "computeProjectionError" is set to 1.
Matt J
on 17 Mar 2025
Yes, perhaps. I'm just saying, I think the OP took a wrong turn with that.
payam samadi
on 18 Mar 2025
I noticed that I made a mistake and the number "1" should be "num_projections".
After replacing 1 by num_projections, your code throws an error. In function "computeProjectionError", you try to access theta_rad(i) for 1 <= i <= num_projections in your loop, but you only pass theta_rad(i) to the function which is a single scalar value.
payam samadi
on 18 Mar 2025
payam samadi
on 28 Mar 2025
Edited: payam samadi
on 28 Mar 2025
0 votes
Dear Matt,
I have carefully reviewed the code and made several refinements. While I observed noticeable improvements compared to the previous version, the results are still not fully meeting my expectations. Below is a summary of my findings and ongoing challenges:
1. Ensuring Matrix B is Positive Definite:
- I tested 24 different methods to enforce positive definiteness on Matrix B. However, only methods 9 and 10 produced meaningful results (see Objective_3.m, lines 57-213).
- To validate this, I had a colleague test the code using a Gravitational Search Algorithm (GSA) for optimization. Interestingly, modifying the method for ensuring B's positive definiteness led to wildly different results (Test 6-GSA.zip file).
- Although multiple methods exist (to make sure Matrix B is positive definite), methods 9 and 10 performed best (though not optimally). However, since they rely on L=zeros(3) or L = tril(randn(3)), enforcing proper bounds on correlation coefficients (ρxy, ρyz, ρzx) becomes ineffective.
2. Improved Mahalanobis-Based Approach (Method 2):
- To address numerical instability due to high condition numbers (nearly singular B), I implemented an improved Mahalanobis-based approach (see Objective_3.m, lines 236-239).
- This modification aims to stabilize computations and enhance reliability.
- Initial Guess for LvecInit:
- Currently, I am using [0,0,0,1,1,1,-0.9,-0.9,0.9] as the initial guess across different samples [1-5].
- While this yields reasonable results, further refinement may improve convergence and accuracy.
3. Bounding Mean Tumor Motion (mx, my, mz):
- I incorporated constraints into Objective_3.m to ensure that mean tumor motion (mx, my, mz) remains within the range -10 ≤ m ≤ 10.
- This helps maintain realistic motion estimations while improving numerical stability.
4. To avoid getting into the local minima:
- Two functions computeError and computeError_1 are proposed to avoid getting into the local minima. The computeError_1 is better, and included two different optimzation; strating with fminsearch and then for the next iteration fmincon is used.
Overall, these updates have resulted in improved stability, but there is still room for this code. Let me know your thoughts on the next steps or if you have any suggestions.
Here is results of the code
Categories
Find more on Optimization Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
