Algorithm Model Details

GitHub Repository

GitHub - yuxinliuhubert/Fluid-tracking: fluid tracking projectGitHub

We developed Phases 1-4 using MATLAB and switched to Python for Phase 7 onward.

Table of Content
Programmed in MATLAB
Phase 1: Static Model, Single Particle
Phase 2: Linear Velocity, Single Particle
Phase 3: Constant Acceleration, Single Particle
Phase 4: Variable Acceleration, Single Particle
Phase 5: Experimental Data Testing for Single Particle
Programmed in Python
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.

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

1. Projection Generation: Generate projections of the particle based on rotations and velocity function.

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

1. Estimating True Position: Estimate particle's position from the projections and derive velocity terms.

proj2r0(xz_proj,theta,SRD,RDD,delta_T)

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

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.

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.

plot3(real_Positions(:, 1), real_Positions(:, 2), real_Positions(:, 3), 'r', 'LineWidth', 2);
plot3(estimatedPositions(:, 1), estimatedPositions(:, 2), estimatedPositions(:

proj2r0 Upgraded to Include Acceleration

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

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