# Algorithm Model Details ## GitHub Repository {% embed url="" %} {% hint style="info" %} **We developed Phases 1-4 using MATLAB and switched to Python for Phase 7 onward.** {% endhint %} | Table of Content | | ------------------------------------------------------------------------------------------------------------ | | **Programmed in MATLAB** | | [Phase 1: Static Model, Single Particle](#phase-1-static-model-single-particle) | | [Phase 2: Linear Velocity, Single Particle](#phase-2-linear-velocity-single-particle) | | [Phase 3: Constant Acceleration, Single Particle](#phase-3-constant-acceleration-single-particle) | | [Phase 4: Variable Acceleration, Single Particle](#phase-4-variable-acceleration-single-particles) | | [Phase 5: Experimental Data Testing for Single Particle](#phase-5-experimental-data-testing-single-particle) | | **Programmed in Python** | | [Phase 7: Variable Acceleration, Multiple Particles](#phase-7-variable-acceleration-multiple-particles) | ## Phase 1: Static Model, Single Particle The main test code, `Test_Cyrus` , performs a series of tests on various methods to compute the percent error in determining the position of a particle based on its projections. The methods are compared over a series of runs to determine their accuracy. Finally, the results are plotted for visual comparison. #### **Inputs** * `Cumulative_NOS`: Cumulative number of shots or projections taken. * `theta_degree`: Rotation angle between projections in degrees. * `noise_ratio`: Ratio of position value to the error for each projection. 1. **Original Method Test**: * This method tests the original approach using pairs of shots. * Results are stored in `Original_Result_PercentError`. 2. **Consecutive Number Average Method Test**: * This approach tests using consecutive shots. * Results are stored in `Consecutive_Result_PercentError`. 3. **Original Method with Bigger Matrix Test**: * This method uses the original approach but with a bigger matrix to account for the z-axis. * Results are stored in `OriginalBigger_Result_PercentError`. 4. **Big Matrix Method Test**: * This method tests using the `BigMatrix2` function. * Results are stored in `BigMatrix_Result_PercentError`. 5. **Consecutive Method with Kalman Filter Test**: * This method tests the consecutive approach but incorporates a Kalman filter. * Results are stored in `Kalman_Result_PercentError`. In each method file: **Inputs** * `true_position_3d`: A 1x3 vector containing the true x, y, z positions of the particle. * `noise`: The amount of random noise introduced in the projection. * `NOS`: Number of shots or projections taken. * `theta_degree`: Rotation angle between projections in degrees. **Output** * `percent_error`: The percent error between the predicted and true x-position of the particle. Testing with synthetic particle data shows that the Big Matrix method had the lowest percent error, so we moved forward with the linear least square approach. ## Phase 2: Linear Velocity, Single Particle `Phase2_pt_3d` aims to estimate the 3D position and velocities of a particle over time using a series of projections. The particle's position is rotated and translated over time. The function constructs a system of equations based on the particle's projected positions and solves them to estimate the particle's position and velocities. **Inputs** * `initial_position_3d`: Initial 3D position of the particle as a column vector `[x0; y0; z0]`. * `noise`: Deviation of the projection measurement from the perfect measurement (68% chance). * `delta_T`: Time interval between shots. * `NOS`: Number of shots or projections. * `theta_degree`: Rotation angle between projections in degrees. **Outputs** * `results`: Estimated 3D position and velocities. **Code Highlights** 1. **Initialization**: Set parameters, velocity function, and errors. ```matlab SRD = 1; % m, Source-Reference Distance RDD = 1; % m, Reference-Detector (screen) Distance method = 0; % 0 for least square, 1 for kalman v = @(t) [3;2;1]; % velocity function ``` 2. **Projection Generation**: Generate projections of the particle based on rotations and velocity function. ```matlab theta = theta_degree/180*pi; % % First shot: without rotation r0_0=initial_position_3d; real_Positions = initial_position_3d; M_p = (SRD+RDD)/(SRD+r0_0(2)); x_proj=M_p*r0_0(1)+randn*noise/2; z_proj=M_p*r0_0(3)+randn*noise/2; for k = 1:NOS-1 r0_k=r0_0+integral(v,0,delta_T*k,'ArrayValued', true); % true location in the original frame of reference real_Positions = [real_Positions; r0_0 + integral(v,0,delta_T*k,'ArrayValued', true)']; r_now=T(r0_k,theta*k); M_p = (SRD+RDD)/(SRD+r_now(2)); % magnification of particle x_proj=[x_proj;M_p*r_now(1)+randn*noise/2]; z_proj=[z_proj;M_p*r_now(3)+randn*noise/2]; end ``` 3. **Estimating True Position**: Estimate particle's position from the projections and derive velocity terms. ```matlab proj2r0(xz_proj,theta,SRD,RDD,delta_T) ``` 4. **Helper Functions**: * **Coordinate Rotation Function (T)**: Rotate a point in 3D space about the y-axis. 
function [r2]=T(r1,alpha) 
   %r1=[x1;y1;z1]and r2=[x2;y2;z2] describe the same physical pt but in 2 coordinate systems with
   % theta degress rotation CCW
   r2=[cos(-alpha) -sin(-alpha) 0; sin(-alpha) cos(-alpha) 0;0 0 1]*r1;
   % note that we are rotating the coordinate axes, not the location vector of the particle
end
* **proj2r0 Function**: Estimate the original position and velocities of the particle using projections. ```matlab function [r0]= proj2r0(proj,theta,SRD,RDD,delta_T) % This function estimates the original 3D position of a particle % and its velocities using a set of projections. % Initialization NOS=height(proj); SDD=(SRD+RDD); row_number_A = round(2*NOS + 2*(NOS-1), 0); col_number_A = round(1 + 2*NOS, 0); A = zeros(row_number_A, col_number_A); b = zeros(row_number_A, 1); % Constructing Equations from Magnification for j = 1:NOS ... end % Transformation Equations for k = 2:NOS ... end % Expanding Matrix to Account for Velocity A = [A, zeros(height(A), 3)]; new_col_num = size(A, 2); % Incorporating Velocity Coefficients for j = 1:NOS A(2*j-1, new_col_num) = delta_T*(j-1); end % Adding Equations for Velocities for k = 2:NOS ... end % Solving the System of Equations x = (A\b); r0 = [x(2), x(3), x(1), x(new_col_num-2), x(new_col_num-1), x(new_col_num)]; end ``` ## Phase 3: Constant Acceleration, Single Particle In this phase, we sought to deal with situations involving a constant acceleration. The `Phase3_pt_3d_success` function builds upon Phase 2 by incorporating an additional movement model that considers constant acceleration. It then visualizes the real and estimated positions in 3D space, offering a tool for assessing the accuracy of the estimation process. ```matlab function results = Phase3_pt_3d_success(initial_position_3d, noise, delta_T, NOS, theta_degree, ``` **Inputs** * `initial_position_3d`: Initial position of the particle in 3D space. * `noise`: The deviation of the projection measurement from a perfect measurement. * `delta_T`: Time interval between successive camera shots. * `NOS`: Number of shots or projections taken. * `theta_degree`: Clockwise rotation angle of the camera before each shot. * `v`: Velocity function representing the particle's movement in 3D space. **Visualization** This segment plots the real positions (red) and estimated positions (blue) of the particle in a 3D space. ```matlab plot3(real_Positions(:, 1), real_Positions(:, 2), real_Positions(:, 3), 'r', 'LineWidth', 2); plot3(estimatedPositions(:, 1), estimatedPositions(:, 2), estimatedPositions(: ``` **`proj2r0` Upgraded to Include Acceleration** ```matlab function [r0]= proj2r0(proj,theta,SRD,RDD,delta_T) % This function aims to predict the original 3D coordinates of a particle % when the camera angle is theta=0 (i.e., before the camera starts rotating). % Initialization: NOS=height(proj); % Number of shots SDD=(SRD+RDD); % Sum of the Source-Reference and Reference-Detector Distances % Pre-allocation for efficiency: ... % Matrix Construction for Magnification: for j = 1:(NOS) ... % This loop constructs equations arising solely from magnification. % For every shot, two new variables are introduced leading to two new equations. end % Matrix Construction for Transformation: for k = 2:(NOS) ... % This loop constructs transformation equations. % These equations account for the rotation of the camera between shots. end % Incorporating New Variables: A=[A,zeros(height(A),6)]; % The matrix A is expanded to incorporate new variables % representing velocity and acceleration components. % Coefficients Related to Velocity and Acceleration: for j = 1:(NOS) ... % Updating the matrix A with coefficients related to velocity and acceleration. end for k=2:NOS ... % Adding the effects of velocity and acceleration to the transformation equations. end % Solving the Equations: x=(A\b); % Solving the system of equations. % Output: r0=[...]; % Contains the estimated 3D coordinates, velocity, and acceleration of the particle. end ``` **Compared to Phase 2**: * **Similarities**: * The core structure of the `proj2r0` function remains the same, with equations being formulated based on magnification and transformation. * The function still aims to predict the 3D coordinates of the particle based on its projections. * **Differences**: * In this phase (Phase 3), the function has been expanded to consider not just initial position and velocity but also constant acceleration. This is reflected in the introduction of acceleration components $$ax​, ay​, az​$$ in the matrix `A`. * The matrix `A` in Phase 3 is expanded further to accommodate six new variables (3 for velocity and 3 for acceleration) compared to Phase 2. * The equations in this phase account for the effect of constant acceleration on the particle's position over time. ## Phase 4: Variable Acceleration, Single Particles The critical aspect of developing this phase is to generalize the movement of a particle, which could be stationary, moving with a linear velocity, with a constant acceleration, or changing accelerations. To address that, Phase 3's code is used over a section of the particle's movement, to estimate its path over a smaller time range, with the assumption that in a smaller time interval, the movement can be characterized by a simple acceleration model. This assumption is valid because the experimental data mostly focuses on particles moving in a slow flow field. Before introducing the strategies, we introduce a variable named `NOS_per_section`. If there exists n NOS (number of snapshots), NOS\_per\_section dictates the number of snapshots that will be fed into the **`proj2r0`** algorithm for an estimated path. In other words, it represents the size of a batch of data that will be analyzed. Two strategies are impleme[^1]nted to divide the data, discussed as below. #### Strategy 1: Serial Given `NOS`, and `NOS_per_section` = M, this strategy breaks the entire snapshot data set into M chunks, each containing a unique set of snapshots. Each chunk is then analyzed through the linear big matrix to get the estimated 3D positions of each time step, termed "subpath." After all subpaths are generated, they are connected in a head-to-tail fashion to form the final estimated 3D path for the particle. The main advantage of the strategy is speed, as every snapshot is only evaluated once. However, it does not work well when the particle is relatively stationary, since fewer iterations on the same time step do not have the resilience against noise. **Strategy 2: Overlap** Following the methodology of Strategy 1, Strategy 2 adopts an overlapping approach to analyze the snapshot data set. In this strategy, the data is divided into overlapping chunks, where each chunk contains a sequence of snapshots from the 1st to the k th, followed by the 2nd to the (k + 1) th, and so on. This creates a series of overlapping sections, each containing a distinct but partially shared set of snapshots. The linear big matrix is then applied to each chunk to estimate the 3D positions for each time step, termed a "subpath." Unlike Strategy 1, where subpaths are connected directly, Strategy 2 involves averaging the 3D position estimates of points that are evaluated multiple times due to the overlap. This results in a more refined and potentially accurate final estimated 3D path for the particle, as it incorporates repeated assessments of the same points, reducing uncertainty and improving precision in the trajectory estimation. It also requires less NOS to produce the same number of batches. However, this strategy does not work well when the angle between each shot is large. *** Here is a diagram for visualizing and comparing the strategies:
## Phase 5: Experimental Data Testing (Single Particle) ## Phase 7: Variable Acceleration, Multiple Particles In this phase, there are two main goals: 1. Convert all code to Python 2. Develop the sorting algorithm that can take a 2D snapshot (image) with the 2D coordinates of the particles, and match each position with a particle. The idea is to classify the positions in each image into shadows of different particles. Then, when all particles receive a time-serial chain of 2D positions, this chain will be fed into the 2D to 3D big matrix to generate the 3D paths. We first converted all code to Python. Take **`projr20`**, one of the fundamental blocks of the algorithm, as an example. Its Python form is as follows: ```python def proj2r0_acc(xz_proj, theta, SOD, ODD, dt): NOS = xz_proj.shape[0] SDD = SOD + ODD row_number_A = int(2 * NOS + 2 * (np.math.factorial(NOS) / (np.math.factorial(NOS - 2) * 2))) col_number_A = int(1 + 2 * NOS) A = np.zeros((row_number_A, col_number_A)) b = np.zeros((row_number_A, 1)) for j in range(NOS): xi_j, zi_j = xz_proj[j] A[2 * j, 0] = 1 A[2 * j, 2 * j +2] = -zi_j / SDD b[2 * j] = zi_j * SOD / SDD A[2 * j + 1, 2 * j+1] = -1 A[2 * j + 1, 2 * j + 2] = xi_j / SDD b[2 * j + 1] = -xi_j * SOD / SDD IoR = 2 * NOS for k in range(1, NOS): trans_count = k A[IoR:IoR + 2 * trans_count, 2 * k+1:2 * k + 3] = np.tile(np.array([[-1, 0], [0, -1]]), (trans_count, 1)) for i in range(trans_count): delta_theta = theta * (i + 1) A[IoR + 2 * i:IoR + 2 * i + 2, 2 * (k - 1) - 2 * i+1:2 * (k - 1) - 2 * i + 2+1] = np.array([[np.cos(delta_theta), np.sin(delta_theta)], [-np.sin(delta_theta), np.cos(delta_theta)]]) IoR += 2 * trans_count A = np.hstack([A, np.zeros((A.shape[0], 6))]) new_col_num = A.shape[1] u_ind, v_ind, w_ind, ax_ind, ay_ind, az_ind = new_col_num - 6, new_col_num - 5, new_col_num - 4, new_col_num - 3, new_col_num - 2, new_col_num - 1 for j in range(NOS): A[2 * j, w_ind] = dt * j A[2 * j, az_ind] = 0.5 * (dt * j)**2 IoR = 2 * NOS for k in range(1, NOS): trans_count = k T2 = dt * k for i in range(trans_count): delta_theta = theta * (i + 1) T1 = T2 - (i+1) * dt theta_prime = theta * (k - i - 1) A[IoR, u_ind] = (np.cos(theta_prime) * np.cos(delta_theta) - np.sin(theta_prime) * np.sin(delta_theta)) * (T2 - T1) A[IoR + 1, u_ind] = (-np.cos(theta_prime) * np.sin(delta_theta) - np.sin(theta_prime) * np.cos(delta_theta)) * (T2 - T1) A[IoR, v_ind] = (np.sin(theta_prime) * np.cos(delta_theta) + np.cos(theta_prime) * np.sin(delta_theta)) * (T2 - T1) A[IoR + 1, v_ind] = (-np.sin(theta_prime) * np.sin(delta_theta) + np.cos(theta_prime) * np.cos(delta_theta)) * (T2 - T1) A[IoR, ax_ind] = 0.5 * (np.cos(theta_prime) * np.cos(delta_theta) - np.sin(theta_prime) * np.sin(delta_theta)) * (T2**2 - T1**2) A[IoR + 1, ax_ind] = 0.5 * (-np.cos(theta_prime) * np.sin(delta_theta) - np.sin(theta_prime) * np.cos(delta_theta)) * (T2**2 - T1**2) A[IoR, ay_ind] = 0.5 * (np.sin(theta_prime) * np.cos(delta_theta) + np.cos(theta_prime) * np.sin(delta_theta)) * (T2**2 - T1**2) A[IoR + 1, ay_ind] = 0.5 * (-np.sin(theta_prime) * np.sin(delta_theta) + np.cos(theta_prime) * np.cos(delta_theta)) * (T2**2 - T1**2) IoR += 2 x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None) r0 = [x[1], x[2], x[0], x[u_ind], x[v_ind], x[w_ind], x[ax_ind], x[ay_ind], x[az_ind]] return r0 ``` Then, I tackled the 2D sorting section, which turned out to be the most difficult part of the algorithm development. The main challenge is that [^1]: --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://www.hubertyl.com/hubert-lius-porfolio/engineering-related/temporal-super-resolution-particle-and-feature-tracking/algorithm-model-details.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.