How to extract 3D position from the 2D projection using (an optimzation method + 3D Gaussian PDF)

Question

payam samadi 2025-3-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2174843-how-to-extract-3d-position-from-the-2d-projection-using-an-optimzation-method-3d-gaussian-pdf

编辑： payam samadi 2025-3-13

Hi there,

I am currently working on a Matlab code for my homework (the attached PDF file). As a beginner with Matlab, I initially get help from ChatGPT, but unfortunately, it hasn't been very effective for me.

Here are some quick guidelines to clarify my approach:

In Step 1 of the code, I define a few system parameters—nothing more. I have assumed that I already have 3D datasets (true_T) and have subsequently used the code to compute 2D projections (lines 20-23) to estimate the projections ([xp, yp]).

In Step 2, I attempted to use the lsqnonlin method, along with the detailed information provided in the PDF file, to obtain the final covariance matrix (A). However, this part is not functioning as expected.

Any assistance or alternative approaches to resolve this issue would be greatly appreciated.

Thank in advance!

clc; 
clear; 
close all;
%% **Step 1: Load 2D Projection Data**
% Imaging system parameters
SAD = 100;  
SID = 150;  
num_projections = 15;  
% Generate projection angles (e.g., uniformly spaced)
degrees = 30;                                  
alpha = linspace(0, degrees, num_projections)';  
theta_rad = deg2rad(alpha);  
% Generate synthetic ground truth 3D position
rng(42); % Fix random seed for reproducibility
true_T = randn(3, num_projections);  
% Projection Model: Compute 2D projections
f_theta = (SAD - (true_T(1, :) .* cos(theta_rad) + true_T(3, :) .* sin(theta_rad))) ./ SID;
Xp = (true_T(1, :) .* sin(theta_rad) - true_T(3, :) .* cos(theta_rad)) ./ f_theta;                    
Yp = true_T(2, :) ./ f_theta;                    
% Initialize 3D position estimates
f_theta = SAD ./ SID;  
y_est = f_theta .* Yp;  
% Initial guess for estimated target positions
x_est = zeros(1, num_projections);
z_est = zeros(1, num_projections);
%% **Step 2: Compute Probability for Each Projection & Optimize B**
% Initial Cholesky decomposition for B
initial_L = ones(3);  
initial_params = ones(3); 
% Define optimization settings for lsqnonlin
opts = optimoptions('lsqnonlin', 'Display', 'iter', 'TolFun', 1e-6);
lb = [];  % No lower bounds
ub = [];  % No upper bounds
% Optimize B using lsqnonlin
fMLE = @(L_vec) computeError(L_vec, Xp, Yp, alpha, SAD, SID);
opt_L_vec = lsqnonlin(fMLE, initial_params, lb, ub, opts);
Unable to perform assignment because the left and right sides have a different number of elements.

Error in solution>computeError (line 122)
    L(tril(true)) = L_vec;

Error in solution>@(L_vec)computeError(L_vec,Xp,Yp,alpha,SAD,SID) (line 44)
fMLE = @(L_vec) computeError(L_vec, Xp, Yp, alpha, SAD, SID);

Error in lsqnonlin (line 221)
            initVals.F = feval(funfcn{3},xCurrent,varargin{:});

Caused by:
    Failure in initial objective function evaluation. LSQNONLIN cannot continue.
% Convert optimized L vector back to matrix form
L_opt = zeros(3,3);
L_opt(tril(true)) = opt_L_vec;
B_opt = L_opt * L_opt';  
disp('Optimized Covariance Matrix B:');
disp(B_opt);
% Compute probability and estimate 3D positions
detB = det(B_opt);
for i = 1:num_projections
    opt_L_vec = lsqnonlin(fMLE, initial_params, lb, ub, opts);
    % Extract new optimzation values
    var_LR  = opt_L_vec(1);
    var_CC  = opt_L_vec(5);
    var_AP  = opt_L_vec(9);
    cov_LR_CC  = opt_L_vec(2);
    cov_LR_AP  = opt_L_vec(3);
    cov_CC_AP  = opt_L_vec(6);
    A_opt = [ var_LR,  cov_LR_CC, cov_LR_AP;
              cov_LR_CC, var_CC, cov_CC_AP;
              cov_LR_AP, cov_CC_AP, var_AP];
    B_opt = inv(A_opt);  % Inverse covariance matrix
    detB_opt = det(B_opt);
    % Projection Point on Imager (Eq. A.4)
    p = [Xp(i) * cos(alpha(i)) - (SID - SAD) * sin(alpha(i));
         Yp(i);
        -Xp(i) * sin(alpha(i)) - (SID - SAD) * cos(alpha(i))];
    % Focus Point (Eq. A.5)
    f = [SAD * sin(alpha(i)); 0; SAD * cos(alpha(i))];
    % Compute Unit Vector Along the Ray (Eq. A.6)
    e = (f - p) / norm(f - p);
    % Compute Standard Deviation Along the Ray (Eq. A.8)
    sigma_ray = 1 / sqrt(e' * B_opt * e);
    % Compute Mean Position Along the Ray (Eq. A.9)
    mu_ray = ((p' * B_opt * e) / (e' * B_opt * e));
    % Compute Scaling Factor (Eq. A.10)
    K = sqrt(detB_opt / (2 * pi)^3) * exp(-0.5 * (p + e * mu_ray)' * B_opt * (p + e * mu_ray));
    % Compute Probability (Eq. A.11)
    P_xp_yp = K * sqrt(2 * pi) * sigma_ray;
    % Estimate 3D coordinates from projection
    X_est(i) = p(1) + e(1) * mu_ray;
    Y_est(i) = p(2) + e(2) * mu_ray;
    Z_est(i) = p(3) + e(3) * mu_ray;
end
%% **Step 3: Plot the 3D Reconstruction**
figure;
plot3(true_T(1,:), true_T(2,:), true_T(3,:), 'ro-', 'LineWidth', 2); hold on;
plot3(X_est, Y_est, Z_est, 'bo-', 'LineWidth', 2);
grid on;
xlabel('Left-Right (LR) Direction (cm)');
ylabel('Superior-Inferior (SI) Direction (cm)');
zlabel('Anterior-Posterior (AP) Direction (cm)');
title('True vs. Estimated 3D Tumor Position');
legend('True 3D Position', 'Estimated 3D Position');
view(3);
%% **computeError Function (lsqnonlin)**
function residuals = computeError(L_vec, Xp, Yp, alpha, SAD, SID)
    % Convert L vector to lower triangular matrix
    L = zeros(3,3);
    L(tril(true)) = L_vec;
    B = L * L';  % Ensure positive definite B
    detB = det(B);
    num_projections = length(alpha);
    residuals = zeros(num_projections, 1);  % Return vector for lsqnonlin
    for i = 1:num_projections
        % Projection Point on Imager
        p = [Xp(i) * cos(alpha(i)) - (SID - SAD) * sin(alpha(i));
             Yp(i);
            -Xp(i) * sin(alpha(i)) - (SID - SAD) * cos(alpha(i))];
        % Focus Point
        f = [SAD * sin(alpha(i)); 0; SAD * cos(alpha(i))];
        % Compute Unit Vector Along the Ray
        e = (f - p) / norm(f - p);
        % Compute Standard Deviation Along the Ray
        sigma_ray = 1 / sqrt(e' * B * e);
        if sigma_ray <= 0
            residuals(i) = 1e6;  % Large penalty
            continue;
        end
        % Compute Mean Position Along the Ray
        mu_ray = ((p' * B * e) / (e' * B * e));
        % Compute Scaling Factor
        K = sqrt(detB / (2 * pi)^3) * exp(-0.5 * (p + e * mu_ray)' * B * (p + e * mu_ray));
        % Compute Probability
        P_xp_yp = K * sqrt(2 * pi) * sigma_ray;
        if P_xp_yp <= 0
            residuals(i) = 1e6;  % Large penalty for invalid values
            continue;
        end
        % Compute residual as the difference from expected probability
        residuals(i) = log(P_xp_yp);
    end
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

William Rose 2025-3-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2174843-how-to-extract-3d-position-from-the-2d-projection-using-an-optimzation-method-3d-gaussian-pdf#answer_1561284

在 MATLAB Online 中打开

@payam samadi,

Good job by you to make a solid effort before asking for help on your homework. You are a good example.

The portion of the Appendix which you attached does not define AP, CC, LR. I assume those stand for anterior-posterior, caudal-cranial, left-right views. Is each one a vector or array or scalar?

Are the lines below doing what you want? My guess is you want views from angles 0, 30, 60,,...,330 degrees. If so, the lines below needs fixing. It's homework, so I'll let you figure it out.

num_projections = 15; 
degrees = 30;                                  
alpha = linspace(0, degrees, num_projections)';  

Instead of using random data, I suggest you use exact data (like points on the corners of a cube) to test your algorithm. Make sure the projections of the cube points are giving the results you expect.

More later.

28 个评论
显示 26更早的评论隐藏 26更早的评论

William Rose 2025-3-6

@payam samadi,

The Appendix says the patient coordinate system is x,y,z, where +x=left, +y=cranial, +z=anterior. I condisider this to be the fixed coordinate system. The source (of x-rays) is at f. The source and the image plane are on opposite sides of the y-axis. The source and the image plane rotate about the y (cranial-caudal) axis. The y-coordinate of f is always 0. SAD and SDD are distances from f to y-axis and from f to image plane. Source-imager rotation angle alpha is measured in the x-z plane, CW from the +z axis, in Figure 9.

The image plane coordinate system is (xp,yp), where +yp is into the page, in Fig.9.

I mis-interpreted equation A.4 when I first read it. I assume you did not make the mistake I made. Eq.A.4 gives the cooridinates of the image point, in the patient coordinate system, in terms of (xp,yp). Eq. A.4 does not tell us what (xp,yp) are, for a given target point.

I recommend you introduce a new constant for the number of target points to be identified. It appears in your script untitled3.m that the number of target points exactly equations the number of projections (num_projections). That seems like an odd choice to me.

William Rose 2025-3-6

编辑：William Rose 2025-3-6

@payam samadi,

The first few lines of your script untitled3.m compute the coordinates Xpand Yp, which are the image plane coordinates for the projections of the target points. I worote a script to do the same. The script also plots the source, the targt points, the image plane, and the projections of the target points opnto the image plane. My script gets different answer for Xp and Yp than your script. I recommend you compare the your code and mine to figure out what acocunts for the difference. If your script has a mistake in the calculation of Xp and Yp, then the later code, to reconstruct the 3-D positions, is not likely to succeed.

I am attaching a modified version of your script: untitled3a.m. It uses 6 target points, and just one angle, 30 degrees. It displays Xp and Yp on the console. I am also attaching my script. It uses the same target points and the same angle and the same SAD and (SID or SDD)*. My script also makes a plot of the source, the targets, the image plane, and the projected points. You can see how your script and my script give different results for [Xp; Yp].

*SDD in the Appendix you provided and in my script; SID in your code.

payam samadi 2025-3-7

编辑：payam samadi 2025-3-7

Dear William,

Many thanks for your time and response.

The portion of the Appendix which you attached does not define AP, CC, LR. I assume those stand for anterior-posterior, caudal-cranial, left-right views. Is each one a vector, or array, or scalar?

Ans: Yes, Indeed. The AP, CC, and LR are anterior-posterior, caudal-cranial, and left-right directions.

Instead of using random data, I suggest you use exact data (like points on the corners of a cube) to test your algorithm.

Ans: Unfortunately, in this stage, I don't have real datasets. So, I decided to create some datasets.

The Appendix says the patient coordinate system is x,y,z, where +x=left, +y=cranial, +z=anterior. I consider this to be the fixed coordinate system. The source (of x-rays) is at f. The source and the image plane are on opposite sides of the y-axis. The source and the image plane rotate about the y (cranial-caudal) axis. The y-coordinate of f is always 0. SAD and SDD are distances from f to the y-axis and from f to the image plane. Source-imager rotation angle alpha is measured in the x-z plane, CW from the +z axis, in Figure 9.

Ans: Correct. Also, for more clarity, I upload the original File of the article.

I checked my code, and I found that I made a mistake in my f_that and Xp definition. Here is the correct version.

f_theta = (SAD - (true_T(1, :) .* sin(theta_rad) + true_T(3, :) .* cos(theta_rad))) ./ SID;

Xp = (true_T(1, :) .* cos(theta_rad) - true_T(3, :) .* sin(theta_rad)) ./ f_theta;

Yp = true_T(2, :) ./ f_theta;

William Rose 2025-3-7

prostatePositionEstimation2009.pdf

@payam samadi,

I downloaded the full paper (P R Poulsen et al., 2009) rather than just the Appendix.

My interpretation of Poulsen 2009 follows:

A. A point-like target in 3D space has random Gaussian motion about its mean position. We obtain 2D projection views of the target, each view from a different angle. Each view provides a 2-vector [Xp;Yp] which are the coordinates of the target projection onto the image plane. The objective is to estimate nine parameters for the target: the location of the mean in 3D space (3 parameters), the variances of its movement distribution along the principal directions (3 parameters), and the covariances between the principal directions (3 parameters). Poulsen provides a series of fomulas in the Appendix of the 2009 paper, to estimate the parameters by minimizing the negative log likelihood.

Your script untitled3.m generates "true" locations, in 3D, for num_projections points, then computes Xp,Yp coordinates for each point, and does it for num_projections different views. You do not add any random motion to the 3D position at each frame. (There are no "true values" for the variances and covariances in your code.) Therefore your Xp,Yp values do not include the effects of random movement about the mean. Therefore you will not be able to estimate the variances and covariances of the target motion, with the simulated data generated by script untilted3.m.

If you want to implement Poulsen's method, then you need to add random mvement to the target. The random mobevement should follow a 3D Gaussian distribution with nine specified true values for the mean, variance, and covariance.

If you want to estimate the target location, and you do not care about variances or covariances of its motion, then you can use singular value decomposition. See https://www.cs.cmu.edu/~16385/s17/Slides/11.4_Triangulation.pdf. Or google for other explanations.

Poulsen et al. 2009 does not provide a formula for [Xp;Yp], given the 3D location [x;y;z] and the projection angle α. Your script untitled3.m includes a formula for this purpose. I derived my own formula, which gives different results from yours. I discussed the difference, and provided code to illustrate the difference, in an earlier comment.

payam samadi 2025-3-7

Dear William,

Many thanks for your quick response.

If you want to implement Poulsen's method, then you need to add random movement to the target. The random movement should follow a 3D Gaussian distribution with nine specified true values for the mean, variance, and covariance.

Ans: Yes, Indeed. My task is to implement Poulsen's method. I can use this method in Matlab (awgn(x,10,'measured') to add white Gaussian noise. I wonder if can you help me in the implementation of Step 2 of this code ( untitled3.m)? The implementation of the formula in the Appendix?

If you want to estimate the target location, and you do not care about variances or covariances of its motion, then you can use singular value decomposition.

Ans: For this purpose, I also reached the following article [1] which used the triangulation method to estimate 3D tumor location [That will be my second task :) ]. Also, I found there might be some errors in their article regarding equations 1 and 2. I sent them an email to ask for more details and waited for their response.

Poulsen et al. 2009 do not provide a formula for [Xp; Yp], given the 3D location [x;y;z] and the projection angle α. Your script untitled3.m includes a formula for this purpose.

Ans: Correct. The Poulsen et al. 2009 study does not provide a formula for [Xp, Yp] calculation. I get them from this article [1, 2].

1. Subsequent kilovoltage imaging-based fiducial triangulation for intra-fractional prostate motion monitoring and correction

2. Reconstruction of implanted marker trajectories from cone-beam CT projection images using interdimensional correlation modeling

William Rose 2025-3-8

untitled4a.m

@payam samadi,

I modified seciton 1 of your scrit. It now adds random motion to each target on each projection. It includes 4 targets. For each target, 9 parameters are specified: the mean posiion (3 parameters), the standard deivaitons along x,y,z directions, and the correlations between variables. The variances and covariances in Poulsen's eq. A.1. can be computed directly and easily from the standard deviations and the correlations. The 3-D positions are saved in array true_T, which has dimensions 3 x num_projections x N. The 2-D projection coordinates are saved in 2D arrays Xp and Yp, which have dimensions num_projections x N. The script displays the theoretical covariance matrix (eq. A.1) for each of the four targets. See comments in the code for more information.

The screenshot below shows part of the console output when I ran the script, with num_projections=10^4. The screenshot also shows how I computed and displayed the sample mean position and the sample covariance matrix for target 4, after the script finished. The sample mean is very close to the theoretical mean for target 4 (0,0,1). The sample covariance is also close to the theoretical covariance (botha re shown in the screenshot). This gives me confidence that Step 1 of the script works as desired.

You now have a script that simulates target movement with the desired statistical properties. The script generates Xp and Yp data for the markers, for each different angle of view. (You will want to change num_projecitons from 1e4 to a value such as 46 or 91.) I recommend that you now make another attempt to implement the minimization algorithm of Poulsen et al., 2009 (step 2 of your script). Specific recommendations for step 2:

Modify function compute_error() so that it returns

You will need to use most of the equations in the Appendix to compute F_MLE. The arguments to function compute_error() will be the arguments you currently show:

computeError(L_vec, Xp, Yp, alpha, SAD, SID)

but L_vec will be a nine-element vector containing

.

Make an outer loop, for i=1:N (where N=number of targets). On each loop pass, do an optimization to find the values of L_vec that minimize F_MLE for target i. Use fmincon() instead of lsqnonlin(), because we are not minimizing a sum of squares. When doing the minimization, use appropriate bounds for the mean values and (positive) standard deviations and correlations (rho values). The correlation bounds should be +-1. A reasonable initial guess for L_vec on each pass is Lvecinit=[0,0,0,1,1,1,0,0,0];.

Good luck.

payam samadi 2025-3-8

Dear William,

Many thanks for your time and effort—I deeply appreciate it.

I reviewed the code (untitled4a.m) and agree with you that both the sample mean and sample covariance are related to the theoretical covariance. Additionally, I came across another study where they used only one target by mathematically converting the positions of four fiducial markers inside the tumor into a single target. However, our code (untitled4a.m) allows flexibility in changing the number of targets. If you agree, I would prefer to proceed with the implementation using four targets for step 2.

I will do my best to finalize this step and will update you on the results ASAP.

Regarding the number of projections, I found that the maximum could be 360 in 1 min (corresponding to a full CBCT rotation). However, in some cases, a lower number of projections is used. While more projections improve optimization, we should also consider the ALARA (As Low As Reasonably Achievable) principle.

Additionally, I found that the inverse covariance target matrix must be positive definite during MLP optimization to ensure its invertibility. By the way, I will do my best to implement the minimization algorithm of Poulsen et al., 2009 (step 2).

William Rose 2025-3-8

在 MATLAB Online 中打开

@payam samadi,

You are welcome.

You are right to ALARA principle. Where I work, all experimental studies must be reviewed and approved beforehand. The review committees at my university and hospital would almost certainly not approve a study that used any more ionizing radiation than what the radiologists require. I predict that even with simulated data, we will find rather large error bars on our estimates of the variances and covariances, when we use 360 projections or fewer.

If you fit the standard deviations and the correlations, and you impose the constraint that -1<rho<+1, then you will be guaranteed to have covariance matrices that are positive definite. The requirement for invertibility is another reason why I think it is a good idea to fit the standard deviations and correlations rather than variances and covariances.

I agree that you should fit each of the N targets independently. That is why I recommended an outer for loop:

for i=1:N
    % Estimate L_vec(i,:) by finding the values for L_vec(i,:) that
    % minimize F_MLE for target i.
end

In my earlier post, I wrote

"The arguments to function compute_error() will be the arguments you currently show:

computeError(L_vec, Xp, Yp, alpha, SAD, SID)

but L_vec will be a nine-element vector containing

."

That was not exactly correct. Matrices Xp and Yp contain projected positions for all N targets. Since we are fitting one target at a time, we should only pass the values of Xp and Yp for the current target. Also, since we have N targets, L_vec should be a N-by-9 array, not a vector. Therefore we should pass variables as follows:

computeError(L_vec(i,:), Xp(:,i), Yp(:,i), alpha, SAD, SID)

Torsten 2025-3-9

编辑：Torsten 2025-3-9

在 MATLAB Online 中打开

Substantial changes had to be made to your objective function (see below).

Further, your initial guess vector for L_vec has length 12, but your vectors lb and ub have length 6. Why ?

Further, your objective function returns NaN because P_xp_yp equals 0.

Further, you overwrite p here for each value of j and you leave the loop for p as set for j = num_projections. You must decide which j to take:

    % Projection Point (Eq. A.4)
    for j = 1:num_projections
        p = [XP(j) * cos(alpha(j)) - (SID - SAD) * sin(alpha(j));
             YP(j);
            -XP(j) * sin(alpha(j)) - (SID - SAD) * cos(alpha(j))];
    end

% untitled4a.m
% Posted by @payam samadi, 2025-03-06. Modified by W. Rose.
% v.3  Original post by Payam Samadi.
% v.3a Adjust SAD, SID, num_projections, true_T to match projectionTest.m.
%      With these changes we see that Xp,Yp from this code do not match 
%      xp,yp from projectionTest.  They should match.
% v.3b Use new formulas for f and Xp, posted by Payam 20250307. Now Xp;Yp
%      from this script matches xp;yp from projectionTest.m.  This version
%      has 1 projection and 6 targets.  This version (like v3 and v3a)
%      does not include random motion of the targets.
% v.4a Add random motion of each target. true_T = 3D array of target 
%      position data, with dimensions 3 x num_projections x N,
%      where N=number of targets.
%      Xp, Yp are 2D arrays with dimensions num_projections x N.
%      Script returns before attempting optimization.
%      theta_rad is now a row vector.
% clc; 
clear; 
% close all;
%% **Step 1: Load 2D Projection Data**
% Imaging system parameters
SAD = 100;              % source-axis distance (cm) 
SID = 150;              % source-image plane distance (cm)
num_projections = 5;   % number of different views
% Generate projection angles
alpha = linspace(0, 360, num_projections);
% alpha=30;
theta_rad = deg2rad(alpha);  % view angles (radians)
% Simulate 3D positions, including random target motion
rng(42); % Fix random seed for reproducibility
% Specify nine parameters for each target  
mx=[0 1 0 0];           % mean x for targets
my=[0 0 1 0];           % mean y for targets
mz=[0 0 0 1];           % mean z for targets
sx=[1.0 1.0 1.0 1.0];   % s.d.(x)=sigma(LR) for each target
sy=[2.0 2.0 2.0 2.0];   % s.d.(y)=sigma(CC) for each target
sz=[3.0 3.0 3.0 3.0];   % s.d.(z)=sigma(AP) for each target
rxy=[0.0 0.0 0.0 0.7];  % corr(x,y) for each target
ryz=[0.0 0.0 0.7 0.7];  % corr(y,z) for each target
rzx=[0.0 0.7 0.7 0.7];  % corr(z,x) for each target
N=length(mx);           % number of targets
mu=[mx;my;mz];          % mean target positions
% Use sigma and rho values to create covariance matrix for each target
sigma=ones(3,3,N);      % 3D array of covariance matrices
for i=1:N
    sigma(:,:,i)=[sx(i)^2, rxy(i)*sx(i)*sy(i), rzx(i)*sz(i)*sx(i);
                  rxy(i)*sx(i)*sy(i), sy(i)^2, ryz(i)*sy(i)*sz(i);
                  rzx(i)*sz(i)*sx(i), ryz(i)*sy(i)*sz(i), sz(i)^2];
end
for i=1:N
    fprintf('sigma(%d)=theoretical covariance matrix for target %d=\n',i,i); 
    disp(sigma(:,:,i))
end
sigma(1)=theoretical covariance matrix for target 1=
     1     0     0
     0     4     0
     0     0     9
sigma(2)=theoretical covariance matrix for target 2=
    1.0000         0    2.1000
         0    4.0000         0
    2.1000         0    9.0000
sigma(3)=theoretical covariance matrix for target 3=
    1.0000         0    2.1000
         0    4.0000    4.2000
    2.1000    4.2000    9.0000
sigma(4)=theoretical covariance matrix for target 4=
    1.0000    1.4000    2.1000
    1.4000    4.0000    4.2000
    2.1000    4.2000    9.0000
% Compute 3D target locations, including random motion 
true_T=zeros(3,num_projections,N);
for i=1:N
    % true_T(:,:,i)=mu(:,i);      % no random motion
    true_T(:,:,i)=mu(:,i)+chol(sigma(:,:,i))'*randn(3,num_projections);
end
% Projection Model: Compute 2D projections
Xp=zeros(num_projections,N);    % allocate array
Yp=zeros(num_projections,N); 
for i=1:N
   f_theta = (SAD - (true_T(1,:,i) .* sin(theta_rad) + true_T(3,:,i) .* cos(theta_rad))) ./ SID;
   Xp(:,i) = (true_T(1,:,i) .* cos(theta_rad) - true_T(3,:,i) .* sin(theta_rad)) ./ f_theta;
   Yp(:,i) = true_T(2,:,i) ./ f_theta;                    
end
fprintf('Xp;Yp for targets, 1st projection:\n')
Xp;Yp for targets, 1st projection:
disp([Xp,Yp]);
   -0.8317    0.4713   -1.4145   -0.7970    2.6802    0.0419    0.5126   -1.0280
    2.3624    3.4083   -7.4665   -5.6539   -2.9984    0.9210    4.0433    3.9104
    1.9359    0.4351    1.3004   -1.6909    2.7377   -0.8796    2.7239    6.6020
    8.2375    3.3786   -4.8335    0.1689   -1.7612   -2.3568    1.4537   -0.7362
    2.6375    1.9386    1.6906   -3.7699   -0.3595   -6.4939    0.1094   -6.8973
%% **Step 2: Compute Probability for Each Projection & Optimize B**
% Initial guess for A Matrix
L_vec = rand(size(sigma), 'like', sigma);
Final_Vec_opt = computeError_a(L_vec, Xp, Yp, sigma, SAD, SID, alpha, num_projections);
Initial_guess = 1×12
    0.2713    0.2809    0.8022    0.7722    0.8155    0.7713    0.1159    0.3309    0.3252    0.8872    0.7132    0.7710
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
Warning: Length of lower bounds is < length(x); filling in missing lower bounds with -Inf.
Warning: Length of upper bounds is < length(x); filling in missing upper bounds with +Inf.
Objective evaluation: P_xp_yp = 0.000000
Error using solution>Objective (line 201)
P_xp_yp is invalid (zero, negative, or NaN). fmincon cannot continue.

Error in solution>@(L_Vec)Objective(L_Vec,i,XP,YP,SAD,SID,alpha,num_projections) (line 115)
    L_Vec_opt = fmincon(@(L_Vec)Objective(L_Vec,i,XP,YP,SAD,SID,alpha,num_projections), Initial_guess, [], [], [], [], lb, ub, @covarianceConstraint, options);        

Error in objfunEvaluator (line 5)
fval = feval(Objfun, x, self.FunArgs.AdditionalParameters{:});

Error in OptimFunctions/objectiveFirstEval (line 645)
                [fval, grad, hess] = self.ObjectiveFunAndGrad(self,self.FunFcn{3},...

Error in fmincon (line 501)
[initVals.f,initVals.g,HESSIAN,funObj] = funObj.objectiveFirstEval(X);

Error in solution>computeError_a (line 115)
    L_Vec_opt = fmincon(@(L_Vec)Objective(L_Vec,i,XP,YP,SAD,SID,alpha,num_projections), Initial_guess, [], [], [], [], lb, ub, @covarianceConstraint, options);        

Caused by:
    Failure in initial objective function evaluation. FMINCON cannot continue.
%% Step 3: Estimate 3D coordinates from projection
% X_est(i) = p(1) + e(1) * mu_ray;
% Y_est(i) = p(2) + e(2) * mu_ray;
% Z_est(i) = p(3) + e(3) * mu_ray;
%% **Step 3: Plot the 3D Reconstruction**
% figure;
% plot3(true_T(1,:), true_T(2,:), true_T(3,:), 'ro-', 'LineWidth', 2); hold on;
% plot3(X_est, Y_est, Z_est, 'bo-', 'LineWidth', 2);
% grid on;
% xlabel('Left-Right (LR) Direction (cm)');
% ylabel('Superior-Inferior (SI) Direction (cm)');
% zlabel('Anterior-Posterior (AP) Direction (cm)');
% title('True vs. Estimated 3D Tumor Position');
% legend('True 3D Position', 'Estimated 3D Position');
% view(3);
% Maximum Likelihood Estimation (MLE) approach
function Final_Vec_opt = computeError_a(L_vec_array, Xp_array, Yp_array, sigma, SAD, SID, alpha, num_projections)
Final_Vec_opt = zeros(size(sigma),'like',sigma);
N = size(L_vec_array, 1);  % Number of targets
for i = 1:N
    L_Vec = L_vec_array(i, :);
    XP = Xp_array(:, i);
    YP = Yp_array(:, i);
    
    %% Optimization with fmincon
    Initial_guess = L_Vec
    lb = [0, 0, 0, -1, -1, -1];     % Lower bounds
    ub = [Inf, Inf, Inf, 1, 1, 1];  % Upper bounds
    options = optimoptions('fmincon', 'Algorithm', 'sqp', 'Display', 'iter');
    L_Vec_opt = fmincon(@(L_Vec)Objective(L_Vec,i,XP,YP,SAD,SID,alpha,num_projections), Initial_guess, [], [], [], [], lb, ub, @covarianceConstraint, options);        
   
    % Final Optimized Vector
    Final_Vec_opt(i,:) = L_Vec_opt;
end
end
function val = Objective(L_Vec,i,XP,YP,SAD,SID,alpha,num_projections)
    sLR  = L_Vec(1);
    sCC  = L_Vec(2);
    sAP  = L_Vec(3);
    rhoLRCC  = L_Vec(4);
    rhoLRAP  = L_Vec(5);
    rhoCCAP  = L_Vec(6);
    % Compute covariance values
    cov_LR_CC  = rhoLRCC * sLR * sCC;
    cov_LR_AP  = rhoLRAP * sLR * sAP;
    cov_CC_AP  = rhoCCAP * sCC * sAP;
    % Covariance matrix A
    A_opt = [sLR^2,  cov_LR_CC, cov_LR_AP;
             cov_LR_CC, sCC^2, cov_CC_AP;
             cov_LR_AP, cov_CC_AP, sAP^2];
    % Inverse covariance matrix
    B_opt = inv(A_opt);  
    detB_opt = det(B_opt);
    
    % Projection Point (Eq. A.4)
    for j = 1:num_projections
        p = [XP(j) * cos(alpha(j)) - (SID - SAD) * sin(alpha(j));
             YP(j);
            -XP(j) * sin(alpha(j)) - (SID - SAD) * cos(alpha(j))];
    end
    
    % Focus Point (Eq. A.5)
    f = [SAD * sin(alpha(i)); 0; SAD * cos(alpha(i))];
    % Compute Unit Vector (Eq. A.6)
    e = (f - p) / norm(f - p);
    
    % Compute σ, μ, K (Eqs. A.8-A.10)
    sigma_ray = 1 / sqrt(e' * B_opt * e);
    mu_ray = ((p' * B_opt * e) / (e' * B_opt * e));
    K = sqrt(detB_opt / (2 * pi)^3) * exp(-0.5 * (p + e * mu_ray)' * B_opt * (p + e * mu_ray));
    % Compute Probability (Eq. A.11)
    % Function to minimize
    P_xp_yp =  K * sqrt(2 * pi) * sigma_ray;  
    %% Debugging Each Component of P_xp_yp
    % % Compute Standard Deviation Along the Ray (Eq. A.8)
    % denominator = e' * B_opt * e;
    % 
    % if denominator <= 0 || isnan(denominator)
    %     error('Invalid denominator in sigma_ray computation: e'' * B_opt * e = %f', denominator);
    % end
    % 
    % sigma_ray = 1 / sqrt(denominator);
    % fprintf('sigma_ray = %f\n', sigma_ray);
    % 
    % % Compute Scaling Factor (Eq. A.10)
    % exponent = -0.5 * (p + e * mu_ray)' * B_opt * (p + e * mu_ray);
    % 
    % % Limit exponent range to prevent NaN/Inf
    % exponent = max(min(exponent, 1), -1);
    % K = sqrt(detB_opt / (2 * pi)^3) * exp(exponent);
    % 
    % fprintf('K = %f\n', K);
    % 
    % % Compute Probability (Eq. A.11)
    % P_xp_yp = K * sqrt(2 * pi) * sigma_ray;
    % 
    % fprintf('P_xp_yp = %f\n', P_xp_yp);
    % 
    % if isnan(P_xp_yp) || P_xp_yp <= 0
    %     error('P_xp_yp is NaN or invalid! Check determinant or B_opt inversion.');
    % end
    % Debugging printouts
    prob = P_xp_yp;
    fprintf('Objective evaluation: P_xp_yp = %f\n', prob);
    
    % Ensure P_xp_yp > 0
    if prob <= 0 || isnan(prob) || isinf(prob)
        error('P_xp_yp is invalid (zero, negative, or NaN). fmincon cannot continue.');
    end
    % Objective function (negative log probability for MLE)
    val = -sum(log(P_xp_yp));
end
function [c, ceq] = covarianceConstraint(L_vec)
    sLR = L_Vec(1);
    sCC = L_Vec(2);
    sAP = L_Vec(3);
    rhoLRCC = L_Vec(4);
    rhoLRAP = L_Vec(5);
    rhoCCAP = L_Vec(6);
    % Compute covariance values
    cov_LR_CC  = rhoLRCC * sLR * sCC;
    cov_LR_AP  = rhoLRAP * sLR * sAP;
    cov_CC_AP  = rhoCCAP * sCC * sAP;
    % Define covariance matrix
    A = [sLR^2,  cov_LR_CC, cov_LR_AP;
         cov_LR_CC, sCC^2, cov_CC_AP;
         cov_LR_AP, cov_CC_AP, sAP^2];
    % Constraint: det(A) must be positive
    detA = det(A);
    
    fprintf('Det(A) = %f\n', detA); % Debugging
    
    if detA <= 0 || isnan(detA)
        warning('det(A) is non-positive or NaN. Adjusting constraints.');
    end
    c = -detA; % Ensure det(A) > 0
    ceq = [];
end

William Rose 2025-3-9

编辑：William Rose 2025-3-9

在 MATLAB Online 中打开

@payam samadi,

[Edit: fix spelling errors. No changes to code.]

[Edit 2: remove line fun=@(x)... from outside the for loop, since it is superceded by a similar line inside the loop.]

Your new script, untitled4b.m, begins with

% untitled4a.m
% Posted by @payam samadi, 2025-03-06. Modified by W. Rose.

It does not include comments explaining how 4b differs from 4a. When you make a new version, change the name in the comment on the first line. Include detailed comments as you modifiy your code. It will help you you remember what you were thinking, if you look at your code in a month or two. And it will help others help you.

In the main program, v.4b, Step 2 has only two lines of code, whereas in v.4a and v.3, Step 2 was many lines. This design change is fine. With this design change, you do not need to initialize L_vec in the main program, and you don't need to pass it to computeError_a(). Then Step 2 will be a single line in th4 main program: a call to computeError_a().

File computeError_a.m begins as follows:

% Maximum Likelihood Estimation (MLE) approach
function Final_Vec_opt = computeError_a(L_vec, Xp, Yp, sigma, SAD, SID, alpha, num_projections)
Final_Vec_opt = zeros(size(sigma),'like',sigma);
N = size(L_vec, 1);  % Number of targets
for i = 1:N
    L_Vec = L_vec(i, :);
    XP = Xp(:, i);
    YP = Yp(:, i);

I encourage students to follow the Matlab default format for a script. This format appears when you click on New > Function:

function [outputArg1,outputArg2] = untitled(inputArg1,inputArg2)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
outputArg1 = inputArg1;
outputArg2 = inputArg2;
end

Include a description of each input and output variable, with units of measure and dimesions, if applicable. In your case, I would not include num_projections as an input, since you can determine it from Xp. With your design change, you do not need or want L_vec as an input to computeError(). Therefore your function computeError_b() could look like this. Note the detailed comments at the start.

function LvecArray = computeError_b(Xp, Yp, SAD, SID, alpha)
%COMPUTEERROR_B Estimate unknown parameters for markers, from projection data.
% For each 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.
% There are N markers 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, N x 9
%   LvecArray(i,:)=[mux,muy,muz,sdx,sdy,sdz,pxy,ryz,rzx] for marker i
%   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.
% Payam Samadi and William Rose, 2025-03-09.
[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=[-1e3 -1e3 -1e3 0 0 0 -1 -1 -1];     % lower bounds 
ub=[1e3 1e3 1e3 1e3 1e3 1e3 1 1 1];     % upper bounds
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
    Lvec=fmincon(fun,LvecInit,[],[],[],[],lb,ub);
    LvecArray(i,:)=Lvec;                % save Lvec in row i of LvecArray
end
end

The function above calls Objective(x, xpv, ypv, SAD, SID, alpha), where xpv and ypv and alpha are vectors with length Np. x is the vector of unknown parameters, length=9. Objective() should return F_MLE=negative log likelihood, which is a scalar.

The function computeError_b(), above, may contain errors. I have not tested it. You must make your own version, or test the version above carefully.

When writing complicated functions, I find it useful to create scripts to test the functions.

payam samadi 2025-3-10

computeError_c.m

Dear William,

Many thanks for your response and valuable tips. I will ensure that in my future codes, I include all necessary comments for clarity.

I have a question regarding the computeError_b.m function that I hope you can help clarify:

Matrix Dimension Issue:

In line 34, the code initializes LvecArray as a (N × 9) matrix: LvecArray = zeros(N, 9); However, shouldn’t it be a (3 × 3 × N) matrix instead? While both structures contain 36 values in total, wouldn’t the current shape cause issues when computing LvecArray in subsequent calculations?

I have incorporated Equations (A.4 - A.11) into the function and made some modifications:

Added optimization options:

options = optimoptions('fmincon', 'Algorithm', 'sqp', 'Display', 'iter');

Integrated Equations (A.5 - A.11) inside fun = @(x) Objective(x, xpv, ypv, SAD, SID, alpha);

Modified the initial guess for Lvec to a random positive number to avoid Inf values in detB_opt.

Despite these changes, the function still encounters an issue where P_xp_yp = Inf. Upon further inspection, I found that:

sigma_ray is producing complex values (x + jy), which propagates into K, resulting in an invalid P_xp_yp.

At this point, I am a bit confused about where the mistake might be. I would greatly appreciate your guidance on this issue.

William Rose 2025-3-10

编辑：William Rose 2025-3-10

@payam samadi,

[Edit: fix spelling errors.]

My matrix LvecArray is N by 9 because I regard my nine unknowns as a vector. My unknowns are 3 means, 3 standard deviations, and 3 correlations: [mx, my, mz, sx, sy, sz, rxy, ryz, rzx]. You can compute the covariance matrix from the standard deviations and the correlations. If your Lvec is 3x3, then I assume it represents the covariance matrix. If so, then where are the unknown mean values?

"Modified the initial guess for Lvec to a random positive number to avoid Inf values in detB_opt."

If you define Lvec as I do, then the initial value I propose should work fine. If you initialize the covariance array to a set of random positive numbers, then it will probably not be a valid covariance matrix. Even if you make it symmetric (for example by adding it to its transpose), it may not be a valid covariance matrix. I have not had time to look at your Objective() funciton. I do not know what sigma_ray is, or why it is complex when you do not expect that.

@Torsten gave a vey thoughtful analysis of your objective funciton above. He is better at math and Matlab than me, so I always value his suggestions.

William Rose 2025-3-11

在 MATLAB Online 中打开

@payam samadi,

As I mentioned before, untitled4c.m, calling Objective0(), gives unexpected results. Specifically, the estimated mean positions (mx,my,mz) for each target are not correct for targets 2 and 4. The estimated variances sx and sz are not correct for any of the targets. The estimated correlations rxy, ryz, rzx are not correct for targets 2, 3, 4.

Possible causes of the unexpected results include:

The x,y,z info for the simulated targets is not what we expect it to be, due to an error in untitled4c.m.
The x,y,z info for the targets is correct, but the Xp,Yp data (projection data) is not what it should be, due to an error in untitled4c.m.
The x,y,z data is OK, and the Xp,Yp data and the alpha data is correct, but there is a mistake in computeError_b.m, which causes incorrect data to be passed to Objective0(), or computeError_b.m calls fmincon() in an incorrect way.
The x,y,z data is OK and Xp and Yp are OK, and computeError_b is OK, but Objective0 contains 1 or more errors.

Item 4 above is probably the most interesting option, but we should try to rule out possibilities 1, 2, 3 before we work on 4. Unless you see something obviouslty wrong with Objective0.m.

To evaluate 1. above, we should compare the statistics for the x,y,z data in true_T() to the expected values. The expected values are specified in untitled4c.m (each column is a different target; there are 4 targets):

% Specify nine parameters for each target  
mx=[0 1 0 0];           % mean x for targets
my=[0 0 1 0];           % mean y for targets
mz=[0 0 0 1];           % mean z for targets
sx=[1.0 1.0 1.0 1.0];   % s.d.(x)=sigma(LR) for each target
sy=[2.0 2.0 2.0 2.0];   % s.d.(y)=sigma(CC) for each target
sz=[3.0 3.0 3.0 3.0];   % s.d.(z)=sigma(AP) for each target
rxy=[0.0 0.0 0.0 0.7];  % corr(x,y) for each target
ryz=[0.0 0.0 0.7 0.7];  % corr(y,z) for each target
rzx=[0.0 0.7 0.7 0.7];  % corr(z,x) for each target

After running untitled4c.m (with 721 projections), we compute and display the sample mean values for the targets, and compare to mx, my, mz in the code above:

>> disp(squeeze(mean(true_T,2)))

-0.0436 1.0871 -0.0073 0.0151

0.0633 0.0248 1.0595 -0.0811

0.1383 0.0842 0.0273 0.9178

The sample means above have the expected values: approximately [0,1,0,0; 0,0,1,0; 0,0,0,1]. I conclude that untitled4c.m creates x,y,z data with the desired mean values.

Compute and display the sample standard deviations for the targets, and compare to sx, sy, sz in the code above:

>> disp(squeeze(std(true_T,0,2)))

1.0007 1.0112 0.9905 1.0047

2.0109 1.9885 1.9583 2.0111

2.9800 2.9357 3.0572 3.0048

The sample st.devs above have the expected values: approximately [1,1,1,1; 2,2,2,2; 3,3,3,3]. I conclude that untitled4c.m creates x,y,z data with the desired standard deviations.

Compute and display correlations between x,y,z for targets 1-4.

>> sampleCorr=zeros(3,4);

>> for i=1:4,CorrEst=corr(true_T(:,:,i)'); sampleCorr(1,i)=CorrEst(1,2); sampleCorr(2,i)=CorrEst(2,3); sampleCorr(3,i)=CorrEst(1,3); end

>> disp(sampleCorr)

0.0995 -0.0751 0.0650 0.7078

0.0669 -0.0175 0.7203 0.7072

-0.0184 0.6987 0.7248 0.6918

Expected values (see code above) is [0,0,0,0.7; 0,0,0.7,0.7; 0,0.7,0.7,0.7]. I conclude that untitled4c.m creates x,y,z data with the desired correlations.

Since the sample means, standard deviations, and correlations are consistent with the expected values, I rule out item 1 above as a likely cause of the unepected results.

I still have items 2,3,4 to check.

William Rose 2025-3-11

在 MATLAB Online 中打开

@payam samadi,

Item 2 in the list in my previous post is to check for erroprs in Xp,Yp. Poulsen et al., 2009, do NOT give a formula for Xp,Yp in terms of x,y,z. First, we inspect the code in untitled4c.m, to see if it agrees with Poulsen et al., 2009.

for i=1:N
    f_theta = (SAD - (true_T(1,:,i) .* sin(theta_rad) + true_T(3,:,i) .* cos(theta_rad))) ./ SID;
    Xp(:,i) = (true_T(1,:,i) .* cos(theta_rad) - true_T(3,:,i) .* sin(theta_rad)) ./ f_theta;
    Yp(:,i) = true_T(2,:,i) ./ f_theta;                    
end

In the code above, the indices k,j,i of true_T(k,j,i) correspond to:

k=1,2,3 <=> Tx,Ty,Tz. j=projection number. i=target number.

f_theta in the loop above does NOT correspond to f in eq.A.5 of Poulsen. f_theta is a vector of length Nproj, while f is a vector if length 3. f_theta is an inverse magnification factor in the x-z plane. If the target is small and close to the axis, then ||true_T()||<<1. In this case, it follows from the formula above that

f_theta ~= SAD/SID.

If SID=150 and SAD=100, then the image plane is 1.5 times farther from the source than the target, so the magnification is 1.5, for a very small target volume.

When the target volume is not small, the magnification is slightly different. The formula for f_theta, above, computes (inverse) magnification in such cases. The inverse magnification magnification looks reasonable, for some special cases. (Such as true_T(:,1,1)=[1,0,1] and alpha=[0, pi/4, pi/2].)

How to extract 3D position from the 2D projection using (an optimzation method + 3D Gaussian PDF)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

28 个评论
显示 26更早的评论隐藏 26更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

How to extract 3D position from the 2D projection using (an optimzation method + 3D Gaussian PDF)

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

28 个评论 显示 26更早的评论隐藏 26更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

28 个评论
显示 26更早的评论隐藏 26更早的评论