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

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

Answer 1

[William Rose]

@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.

[William Rose]

@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]

@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]

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]

@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.

[William Rose]

@payam samadi,

The reproduction of the method of Poulsen et al 2009 is not easy! It seems far too difficult for a homework assignment. I am impressed that you are trying to do it. I recommend that you try to get it to work for one target.

[payam samadi]

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]

@payam samadi,

I inserted your corrected formulas for f and Xp into your script. Now the Xp;Yp values from your script match the values I get with projectionTest. The version of your script, as I have modified it so far, is untitled3b.m. It generates 1 projection and it has 6 targets. There is no provision for random motion of the targets.

The next task is to modifiy Step 1 of the script, to include random motion of the targets. Since you want to estimate 9 unknown parameters for each target, I suggest you have at least 45 observations: 5 observations per unknown parameter. Therefore use num_projections >= 45.

[William Rose]

@payam samadi,

Poulsen et al, 2009 recommend adjusting parameters

varLR, varCC, varAP, covLRCC, covLRAP, covCCAP

to minimize the negative log likelihood of the observed data (see eq.A.1 and text between eqs. A.2 and A.3).

I recommend fitting

sLR, sCC, sAP, rhoLRCC, rhoLRAP, rhoCCAP

(where s indicates standard deviation) instead. The reason is that if you fit the parameters recommended by Poulsen, the minimization routine may find a best-fit value for the covariances that is impossible for real or simulated data. For actual x-y data, cov(x,y)=rho(x,y)*s.d.(x)*s.d.(y). Therefore abs(cov(x,y)) must be <= sqrt(var(x)*var(y)). It may not be easy to enforce this constraint when doing the minimixzation, if we are adjusting cov(x,y), etc. But if we adjust rho(x,y) and rho(y,z) and rho(z,x), it will be easy to include the constraint -1<=rho(x,y)<=+1.

[William Rose]

@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]

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]

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

[payam samadi]

Dear William,

I made another attempt to implement the second part of Poulsen et al., 2009, but unfortunately, I encountered issues.

• In my first attempt, I used:

objective = @(L_Vec) -sum(log(P_xp_yp(L_Vec)));

as the objective function (Please see line 94 in computeError_a.m). However, this resulted in the following error:

Error using sqpInterface

Objective function is undefined at initial point. Fmincon cannot continue.

• In my second attempt, I modified the approach and defined:

objective = @(L_Vec) Objective(L_Vec, P_xp_yp);

(Please see line 95 in computeError_a.m). This time, I received the following error:

Objective evaluation: P_xp_yp = 0.000000

Error using Objective

P_xp_yp is invalid (zero, negative, or NaN). fmincon cannot continue.

I have attached my code (untitled4b.m, computeError_a.m, covarianceConstraint.m, and Objective.m). I would be truly grateful if you could take some time to review it and help me identify where I made a mistake.

Looking forward to your insights.

[Torsten]

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]

@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]

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]

@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]

@payam samadi,

In your new function computeError_c.m, I would move the code starting with

%% Add Equations from A.4 to A.11

and ending just before

%% Define the Function to Minimize

to Objective(). If you do not do that, then fmincon() cannot adjuist the value of FMLE as it tries different values for the unknown vector x. When you do this, you may need to adjust the code to correctly refer to the elements of the dummy variable x, or L_Vec, or whatever you call the dummy variable containing the unknown parameters, in Objective.m.

[payam samadi]

@Torsten and @William Rose

I am reconsidering whether starting with %% Add Equations from A.4 to A.11 is the best approach, as we need xpv for each marker in the subsequent steps. I am currently working on modifying the implementation to ensure it runs correctly.

In the meantime, if you have any suggestions on how to make it functional, I would greatly appreciate your insights. Please let me know if you come up with any ideas.

Thanks in advance for your time.

[William Rose]

@payam samadi,

The attached main program plus two functions runs without error. Foir this demonstration, I used 721 projections (every half degree), in hopes of attaining more accurate esitmates of the unknown parameters.

The script untitled4c.m runs wqithout errors. It gives answers that are not crazy - but the answers are not correct. After running untitled4c, I get the console results shown below. The 4x9 array Final_Vec_opt contains the 9 estimated parameters for the each of the 4 markers. The mean values mx, my,mz (columns 1,2,3) for markers 2 and 4 don't look right (red, green boxes on screenshot). The standard deviations for all 4 markers should be about sdx=1, sdy=2, sdz=3 (blue box on screenshot). But the results are sdx~=sdy~=sdz=2 fror all four markers. And the correlations are also not as expected. But they are all in the range -1 to +1, as expected.

Even though the answers are not correct, the answers are also not crazy. I hope that Objective0.m will give you an idea of what you could do. Maybe you can find the cause of the errors.

[payam samadi]

Dear William,

Many thanks for your time and for providing such valuable code.

I have added the function Estimate_3D_coordinates_from_Final_Vec_opt.m to estimate the 3D tumor position from Final_Vec_opt and have updated the code to version untitled4d.m.

However, I have noticed that each time I run the code, I obtain a slightly different Final_Vec_opt. While the optimization converges, it does not do so exactly in the way we expect. I was wondering if you might have any insights into why this occurs.

In meantime, I will work on the Objective0.m to see what I can do.

[William Rose]

@payam samadi, I commented out the line rng(42) because I wanted to see if the results were due to an unlucky choice of the initial random seed. The results do not change a small amount from run to run, but not significantly.

[William Rose]

@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:

1. The x,y,z info for the simulated targets is not what we expect it to be, due to an error in untitled4c.m.
2. 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.
3. 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.
4. 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]

@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].)

More on Item 2 later.

[William Rose]

@payam samadi,

Earlier I developed my own equaiton for the coordinates of the projection of the target, expressed in 3D cooridnates and in image-plane coordinates. It is a different approach than the formula you use. In script projectionTest.m, I compute the projected coordinates Xp,Yp for six targets, for one angle. I do this calculatin using your formula and my formula. I also usi eq.A.4 to back-calculate the projected point coordinates in 3D from the 2D coordinates obtained with your formula. The distance between the 3D coords computed in twpo ways shoudl be the same, and the distance between the 2D coordinates (image plane coordinates) shoudl also be very small. The script computes and distays those distances. I tried angles from 0 to 359 degrees. In every case, the distance between the points is very small - on the order of 10^-14. This indicates that the projected points Xp,Yp are correct. Therefore I think Item 2 on the list above is not the cause of the unexpected results with untitled4c.m.

[payam samadi]

Dear William,

Many thanks for your time.

I agree that Items 1-2 are working correctly.

• If the target is small and close to the axis, then ||true_T()||<<1. In this case, it follows from the formula that:

f_theta ~= SAD/SID​

• For this part, please refer to Equation (6) in the attached article, Page 3, which states that "the coordinates of the target position X and Z are usually much smaller than the SAD."

Additionally, there is one important factor we forgot to include in our code:

• Since there are four fiducial markers inside the tumor (e.g., for a moving tumor), we should consider that each fiducial marker is moving over time in 3D space (x,y,z;t)(x, y, z; t)(x,y,z;t).
• The time parameter ttt should match the number of projections.
• Thus, in summary, we should have four targets, with their positions changing at each projection. However, the motion of these targets should not exceed 5 mm, meaning they should remain within 0-5 mm in three directions.

I believe that incorporating this consideration will improve our results. What do you think?

Looking forward to your thoughts.

[payam samadi]

Dear William,

I modified the FMLE = Objective 1(x, xp, yp, SAD, SAD, alpha) and uploaded a new version. In this version, I Included compute log-likelihood (Eq. A.12) into the for loop.

Also, do you think it's a good idea that we used B = pinv(A) instead of B = inv(A)?

Looking forward to your thoughts.

[William Rose]

@payam samadi, please send me a secure message so we can continue this dicussion offline. If you click on the "WR" circle next to my posts, you will see an envelope icon at top right of the pop-up window. Click it to send me a message. Thank you.

[payam samadi]

@William Rose

Dear William,

I sent you a message. If you got that, please let me know, any time.