Charged particle trajectory tracking and data acquisition method and system based on bidirectional vision

Abstract

The invention discloses a charged particle trajectory tracking and data acquisition method and system based on bidirectional vision, and relates to the technical field of machine vision, a charged particle two-dimensional image sequence is obtained through a plurality of image devices, a geometric transformation model containing lens distortion compensation and dynamic baseline correction is established for multi-view fusion, and a two-dimensional image sequence is obtained. Sub-pixel precision three-dimensional positioning is realized; time sequence correlation is carried out based on timestamps and spatial continuity to obtain trajectory data, trajectory quality is calculated according to motion parameters and field intensity information, and geometric parameters are reversely corrected to form a closed-loop optimization mechanism.

Description

Technical Field

[0001] This invention relates to machine vision technology, and more particularly to a method and system for tracking and acquiring data on the trajectory of charged particles based on bidirectional vision. Background Technology

[0002] Traditional charged particle trajectory tracking technologies mainly rely on single-view imaging and fixed-parameter image processing methods, which are insufficient to handle the tracking needs of high-speed particle motion in complex electric field environments. Existing technologies suffer from three main problems: First, single-view imaging lacks depth information, resulting in insufficient 3D spatial positioning accuracy, especially when charged particles move along the optical axis; second, the geometric transformation model does not consider changes in the relative position of the device caused by environmental disturbances, and geometric parameter drift during long-term operation leads to a significant decrease in positioning accuracy; third, there is a lack of organic integration between physical constraints and image processing, making it impossible to utilize the motion patterns of charged particles in electric fields for trajectory optimization and parameter correction. Summary of the Invention

[0003] This invention proposes a method and system for tracking and acquiring the trajectory of charged particles based on bidirectional vision, which can solve the problems in the prior art.

[0004] A first aspect of the present invention provides a method for tracking and acquiring the trajectory of charged particles based on bidirectional vision, comprising:

[0005] Two-dimensional image sequences of charged particles in an electric field environment are acquired synchronously by multiple image acquisition devices with spatial relationships.

[0006] Based on the spatial positional relationship of the image acquisition device, a geometric transformation model including lens distortion compensation and dynamic baseline correction is established to perform multi-view fusion transformation on the pixel coordinates of charged particles in the two-dimensional image sequence, thereby obtaining a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel accuracy;

[0007] Based on the timestamp information and spatial continuity constraints in the position coordinate sequence, the three-dimensional spatial position coordinates at different times are temporally correlated to obtain the three-dimensional trajectory data of the charged particle;

[0008] Based on the motion state parameters and electric field strength distribution information in the three-dimensional trajectory data, trajectory quality evaluation indicators are calculated through force balance relationships;

[0009] Based on the trajectory quality assessment index and the physical constraint relationship, the error of the geometric transformation parameters in the geometric transformation model is derived in reverse, and the geometric transformation parameters are corrected. The corrected geometric transformation parameters are then synchronously applied to the multi-view fusion transformation and the temporal correlation to update the position coordinate sequence, forming a closed-loop image processing optimization mechanism.

[0010] Optionally, establishing a geometric transformation model that includes lens distortion compensation and dynamic baseline correction includes:

[0011] By acquiring a sequence of calibration images containing known spatial location calibration points, the lens distortion parameters and relative position and attitude parameters of multiple image acquisition devices are simultaneously solved. The lens distortion parameters and the relative position and attitude parameters are then integrated into a unified parameter matrix to establish a geometric transformation model that couples lens distortion compensation and baseline correction.

[0012] The system monitors the pixel position changes of calibration points in multiple image acquisition devices in real time. When the detected pixel position change exceeds a preset change threshold, it determines that the relative positions of the multiple image acquisition devices have changed.

[0013] Based on the pixel position change, the reprojection error of the calibration point under different viewpoints is calculated. By minimizing the reprojection error, the lens distortion parameters and relative position and attitude parameters in the unified parameter matrix are jointly optimized and updated to obtain the updated geometric transformation model.

[0014] Optionally, performing multi-view fusion transformation on the pixel coordinates of charged particles in the two-dimensional image sequence to obtain a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel accuracy includes:

[0015] The geometric transformation model is used to perform distortion correction and viewpoint transformation on the pixel coordinates of charged particles in the two-dimensional image sequence, thereby obtaining the corrected pixel coordinates corresponding to each image acquisition device;

[0016] The electric field direction information of the electric field environment is obtained, the angle between the observation direction of each of the multiple image acquisition devices and the electric field direction is calculated, and a confidence weight is assigned to each image acquisition device according to the angle.

[0017] The corrected pixel coordinates are weighted and fused according to the confidence weight, and the candidate position coordinate set of charged particles in three-dimensional space is solved by multi-view geometric constraints;

[0018] Based on the physical constraints of the motion of charged particles under the influence of electric field forces in an electric field environment, the consistency metric between the direction of the electric field force and the direction of motion of the charged particles corresponding to each candidate position coordinate in the candidate position coordinate set is calculated. Candidate position coordinates with consistency metrics lower than the consistency threshold are eliminated, and candidate position coordinates that satisfy the force equilibrium condition are retained as the optimal position coordinates.

[0019] The optimal position coordinates are subjected to temporal interpolation and sub-pixel refinement to obtain a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel precision.

[0020] Optionally, the steps of acquiring the electric field direction information of the electric field environment, calculating the angle between the observation direction of each of the multiple image acquisition devices and the electric field direction, and assigning confidence weights to each image acquisition device according to the angles include:

[0021] The local field strength vector at the location of the charged particles in the electric field environment is obtained, and the direction of the local field strength vector is extracted as the field strength direction information; the angle between the optical axis direction of the multiple image acquisition devices and the field strength direction is calculated respectively.

[0022] Image quality parameters of the two-dimensional images acquired by each image acquisition device are obtained, including image sharpness and signal-to-noise ratio; an angle weighting coefficient is determined based on the numerical range of the included angle. When the included angle is greater than a preset vertical angle threshold, the angle weighting coefficient increases monotonically as the included angle increases; when the included angle is less than or equal to the preset vertical angle threshold, the angle weighting coefficient decreases monotonically as the included angle decreases.

[0023] The angle weight coefficient is combined with the image quality parameter to obtain the initial confidence weight; the initial confidence weight is then normalized to obtain the final confidence weight.

[0024] Optionally, the step of temporally associating the three-dimensional spatial position coordinates at different times based on the timestamp information and spatial continuity constraints in the position coordinate sequence to obtain the three-dimensional trajectory data of the charged particle includes:

[0025] Extract timestamp information from the location coordinate sequence and sort them according to the timestamp order to obtain a time-series arrangement of location coordinates;

[0026] Calculate the spatial distance and velocity change between adjacent position coordinates in the time-series position coordinate arrangement. When the spatial distance exceeds the position deviation threshold calculated based on the motion velocity or the velocity change exceeds the acceleration threshold calculated based on the direction of the electric field force, it is determined that there is a trajectory break.

[0027] For coordinates where there is a trajectory break, based on the velocity direction and electric force direction of the position coordinates before the break, the interpolated position coordinates between the break times are calculated, and the interpolated position coordinates are inserted into the time-series of position coordinates.

[0028] The similarity in spatial location distribution and motion direction between different trajectory segments in the repaired temporal arrangement position coordinate sequence is calculated. When the similarity is higher than the association threshold, the trajectory segments belonging to the same charged particle are spliced ​​together in the order of timestamps. The spliced ​​position coordinate sequence is then smoothed and filtered to obtain the three-dimensional trajectory data of the charged particle.

[0029] Optionally, the step of calculating the trajectory quality evaluation index based on the motion state parameters in the three-dimensional trajectory data and the field strength distribution information of the electric field environment through force balance relationships includes:

[0030] The position coordinates and timestamp information of charged particles at different times are extracted from the three-dimensional trajectory data. The instantaneous velocity of the charged particles is obtained by calculating the difference between the position coordinates of adjacent times, and the instantaneous acceleration is also calculated.

[0031] Obtain the field strength distribution information of the electric field environment, and query the local field strength vector corresponding to the position of the charged particle based on the position coordinate sequence;

[0032] Based on the force balance relationship of charged particles under the action of electric force, gravity and air resistance in an electric field environment, the theoretical acceleration of charged particles is calculated according to the local field strength vector.

[0033] The deviation between the instantaneous acceleration and the theoretical acceleration is calculated as a trajectory quality evaluation index.

[0034] Optionally, the steps of deriving the error of the geometric transformation parameters in the geometric transformation model based on the trajectory quality evaluation index and the physical constraint relationship, correcting the geometric transformation parameters, and synchronously applying the corrected geometric transformation parameters to the multi-view fusion transformation and the temporal correlation to update the position coordinate sequence, forming a closed-loop image processing optimization mechanism, include:

[0035] When the trajectory quality evaluation index exceeds the preset quality threshold, it is determined that there is an error in the geometric transformation parameters in the geometric transformation model;

[0036] A physical constraint equation is established based on the force equilibrium relationship of charged particles in an electric field environment. The physical constraint equation describes the relationship between the motion state parameters and position coordinates of charged particles.

[0037] Substitute the trajectory quality evaluation index into the physical constraint equation, calculate the error of the position coordinate sequence by reverse derivation, and derive the error of the geometric transformation parameters in the geometric transformation model based on the error of the position coordinate sequence.

[0038] The geometric transformation parameters are corrected based on the geometric transformation parameters to obtain the corrected geometric transformation parameters;

[0039] The corrected geometric transformation parameters are synchronously applied to the multi-view fusion transformation of the charged particle pixel coordinates in the two-dimensional image sequence, and the updated position coordinate sequence is recalculated. The updated position coordinate sequence is then applied to the temporal correlation to obtain the updated three-dimensional trajectory data.

[0040] Based on the updated 3D trajectory data, the trajectory quality evaluation index is recalculated to form a closed-loop image processing optimization mechanism.

[0041] Secondly, a bidirectional vision-based system for tracking and acquiring the trajectory of charged particles is provided, including:

[0042] The first unit is used to simultaneously acquire a two-dimensional image sequence of charged particles in an electric field environment through multiple image acquisition devices with spatial positional relationships;

[0043] The second unit is used to perform multi-view fusion transformation on the pixel coordinates of charged particles in the two-dimensional image sequence based on the spatial positional relationship of the image acquisition device, by establishing a geometric transformation model that includes lens distortion compensation and dynamic baseline correction, to obtain a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel accuracy;

[0044] The third unit is used to temporally correlate the three-dimensional spatial position coordinates at different times based on the timestamp information and spatial continuity constraints in the position coordinate sequence, to obtain the three-dimensional trajectory data of the charged particle;

[0045] The fourth unit is used to calculate the trajectory quality evaluation index based on the motion state parameters in the three-dimensional trajectory data and the field strength distribution information of the electric field environment through the force balance relationship;

[0046] The fifth unit is used to deduce the error of the geometric transformation parameters in the geometric transformation model based on the trajectory quality evaluation index and the physical constraint relationship, and to correct the geometric transformation parameters. The corrected geometric transformation parameters are then applied synchronously to the multi-view fusion transformation and the temporal correlation to update the position coordinate sequence, forming a closed-loop image processing optimization mechanism.

[0047] Thirdly, a computer-readable storage medium is provided, having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0048] This invention achieves sub-pixel-level precision 3D spatial positioning through multi-view imaging and fusion transformation, solving the problem of insufficient depth information in traditional single-view methods. The geometric transformation model, including lens distortion compensation and dynamic baseline correction, effectively addresses device position changes caused by environmental disturbances, significantly improving system stability during long-term operation. It innovatively combines physical constraints with image processing, utilizing the force balance relationship of charged particles in an electric field environment to evaluate trajectory quality, which not only improves the accuracy of trajectory reconstruction but also provides a physical basis for parameter optimization; a complete closed-loop image processing optimization mechanism is established. Attached Figure Description

[0049] Figure 1 This is a flowchart illustrating the method for tracking and acquiring the trajectory of charged particles based on bidirectional vision, according to an embodiment of the present invention. Detailed Implementation

[0050] The invention will now be described in further detail with reference to the accompanying drawings. Figure 1 This is a flowchart illustrating the method for tracking and acquiring the trajectory of charged particles based on bidirectional vision according to the present invention. Figure 1 As shown, the method includes:

[0051] Two-dimensional image sequences of charged particles in an electric field environment are acquired synchronously by multiple image acquisition devices with spatial relationships.

[0052] Based on the spatial positional relationship of the image acquisition device, a geometric transformation model including lens distortion compensation and dynamic baseline correction is established to perform multi-view fusion transformation on the pixel coordinates of charged particles in the two-dimensional image sequence, thereby obtaining a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel accuracy;

[0053] Based on the timestamp information and spatial continuity constraints in the position coordinate sequence, the three-dimensional spatial position coordinates at different times are temporally correlated to obtain the three-dimensional trajectory data of the charged particle;

[0054] Based on the motion state parameters and electric field strength distribution information in the three-dimensional trajectory data, trajectory quality evaluation indicators are calculated through force balance relationships;

[0055] Based on the trajectory quality assessment index and the physical constraint relationship, the error of the geometric transformation parameters in the geometric transformation model is derived in reverse, and the geometric transformation parameters are corrected. The corrected geometric transformation parameters are then synchronously applied to the multi-view fusion transformation and the temporal correlation to update the position coordinate sequence, forming a closed-loop image processing optimization mechanism.

[0056] For example, multiple image acquisition devices are configured in the electric field environment. Each device uses an industrial camera with a fixed-focus lens, with a camera resolution of 1920×1080 pixels and a frame rate of 100 frames per second. The spatial relationship between the image acquisition devices is fixed by a rigid bracket made of a low-thermal-expansion-coefficient alloy (coefficient of thermal expansion less than 1×10⁻⁶ / ℃) to ensure stable spatial relationships during long-term operation. When three image acquisition devices are configured, their optical centers form an equilateral triangle with a side length of 1000 mm, and the optical axes of the three devices point towards the center of the electric field environment. A synchronization trigger signal is issued by the main controller. The trigger signal is a TTL level pulse; the rising edge of the pulse triggers synchronous exposure of all image acquisition devices. The trigger signal transmission delay is less than 1 μs, ensuring that the synchronization error is controlled within 1 μs. The image acquisition devices acquire two-dimensional images of charged particles in the electric field environment. The image format is grayscale, with pixel values ​​ranging from 0 to 255. The image storage format is lossless compression, and the file size of a single frame is approximately 2 MB. Two-dimensional image sequences are stored in order of timestamps with a timestamp precision of 1ms. The timestamp information is recorded in the metadata area of ​​the image file.

[0057] The geometric transformation model establishment process involves acquiring a sequence of calibration images, including those of the calibration board. The calibration board uses a checkerboard pattern with a corner spacing of 10mm. The corner position accuracy is calibrated using a laser rangefinder, achieving an accuracy of no less than 0.1mm. The image acquisition device simultaneously acquires images of the calibration board. Image processing algorithms extract the pixel coordinates of the checkerboard corners, and corner detection employs a sub-pixel-level positioning algorithm with an accuracy of 0.1 pixels. A projection relationship from 3D spatial coordinates to pixel coordinates is established, encompassing lens distortion and relative position / attitude. Lens distortion parameters include radial and tangential distortion coefficients. The radial distortion coefficient includes second-order and fourth-order terms, while the tangential distortion coefficient contains two independent components. Relative position / attitude parameters include a rotation matrix and a translation vector. The rotation matrix is ​​represented by rotation angles, and the translation vector represents the relative spatial position between the image acquisition devices. The optimization algorithm employs nonlinear least squares iterative solutions, aiming to minimize the reprojection error of all calibration points. The reprojection error is calculated as the Euclidean distance between the theoretical and actual pixel coordinates. The initial parameters of the iterative process are roughly estimated based on the physical installation location of the image acquisition device. The iteration terminates when the change in reprojection error is less than 0.001 pixels or the number of iterations exceeds 500. After optimization, the lens distortion parameters and relative position and attitude parameters of all image acquisition devices are obtained and integrated into a unified parameter matrix for storage. Dynamic baseline correction is achieved by real-time monitoring of the calibration point pixel position. The monitoring frequency is set to acquire a monitoring image containing the calibration board every 10 seconds. When a change in the calibration point pixel position is detected to exceed 0.5 pixels, the parameter update process is triggered to re-optimize the relative position and attitude parameters in the unified parameter matrix.

[0058] Multi-view fusion transformation processes the pixel coordinates of charged particles in a 2D image sequence. The pixel coordinates are obtained by detecting the center of the charged particle's light spot using image processing algorithms. Light spot detection employs threshold segmentation and centroid calculation, achieving a centroid calculation accuracy of 0.1 pixels. A geometric transformation model extracts corresponding lens distortion parameters based on the pixel coordinates and the image acquisition device number, performing distortion correction. Distortion correction converts the original pixel coordinates into normalized image coordinates. The normalization process uses inverse iterative solving, with the iteration terminating when the change in normalized coordinates is less than 0.0001 or the number of iterations exceeds 10. Normalized image coordinates are combined with rotation matrices and translation vectors to calculate 3D spatial rays. The ray origin is the optical center position of the image acquisition device, and the ray direction is the 3D direction vector obtained after the normalized image coordinates are transformed by the rotation matrix. Multi-view fusion calculates the nearest point for all ray pairs. The nearest point between two rays is obtained by solving for the parameter value that minimizes the distance between the two rays, and the set of nearest points constitutes candidate position coordinates. Based on the electric field direction information, a confidence weight is assigned to each image acquisition device. This confidence weight is calculated by combining the angle between the observation direction and the electric field direction, as well as image quality parameters; the weight is highest when the angle is close to 90°. Candidate position coordinates are weighted and fused according to their confidence weights, with the fusion weight being the product of the confidence weights of the two image acquisition devices. The fused 3D spatial positions are further filtered through physical constraints. These constraints calculate a consistency metric between the direction of charged particle motion and the direction of the electric field force; candidate position coordinates with a consistency metric below 0.5 are discarded. The retained candidate position coordinates are used as the optimal position coordinates. The 3D coordinate values ​​are then fine-tuned through iterative optimization, with the optimization objective being to minimize the reprojection error. The iteration terminates when the reprojection error is less than 0.05 pixels or the number of iterations exceeds 20. The optimal position coordinates at all times constitute a position coordinate sequence with sub-pixel accuracy.

[0059] The temporal correlation process extracts timestamp information from the position coordinate sequence and sorts the positions by timestamp values ​​in ascending order to obtain a temporally arranged position coordinate sequence. The spatial distance and velocity change between adjacent position coordinates are calculated. The spatial distance is the Euclidean distance between the two position coordinates, and the velocity change is the magnitude of the difference between adjacent instantaneous velocity vectors. The position deviation threshold is calculated based on the previous velocity magnitude × timestamp difference × deviation coefficient 1.5, and the acceleration threshold is calculated based on the theoretical acceleration derived from the electric field force × coefficient 2.0. When the spatial distance exceeds the position deviation threshold or the velocity change exceeds the acceleration threshold, a trajectory break is determined. The interpolation calculation for the trajectory break position uses a kinematic model, assuming the charged particle is subjected to a constant electric field force. The interpolated position coordinates are calculated as initial position + initial velocity × time interval + 0.5 × acceleration × time interval 2, applied to the X, Y, and Z directions respectively. The interpolated position coordinates are inserted into the position coordinate sequence in timestamp order to repair the trajectory break. The repaired position coordinate sequence contains multiple trajectory segments. The similarity between different trajectory segments is calculated, including spatial distribution similarity and motion direction similarity. When the overall similarity is higher than 0.7, the two trajectory segments are concatenated into a single trajectory. The concatenated position coordinate sequence is then smoothed and filtered using a moving average algorithm with a window length of 5 position coordinates. The filtered result is the 3D trajectory data.

[0060] The trajectory quality assessment index is calculated by extracting position coordinates and timestamps from the 3D trajectory data. Instantaneous velocity is calculated by the difference between adjacent position coordinates, and instantaneous acceleration is calculated by the difference between adjacent instantaneous velocities. The electric field strength distribution information is stored in the form of 3D spatial grid data with a grid node spacing of 10mm. The local field strength vector can be obtained by querying or interpolating based on the position coordinates. The theoretical acceleration is calculated based on the force balance relationship. The forces include electric force, gravity, and air resistance. The electric force is the field strength vector × charge 1.6 × 10⁻¹⁹ C, gravity is the mass 1 × 10⁻¹⁵ kg × gravitational acceleration 0.0098 mm / ms², and air resistance is the drag coefficient 3.4 × 10⁻⁷ N·s / m × velocity vector. The resultant force ÷ mass yields the theoretical acceleration. The deviation vector is instantaneous acceleration minus theoretical acceleration, and the root mean square value of the deviation vector magnitude is used as the trajectory quality assessment index.

[0061] The closed-loop image processing optimization mechanism compares the trajectory quality assessment index with a preset quality threshold of 0.5 mm / ms². If the index exceeds this threshold, errors in the geometric transformation parameters are identified. The physical constraint equations are established based on force equilibrium relationships, and the residuals of these equations represent the trajectory quality assessment index. Backward derivation calculates the impact of position coordinate errors on the residuals through sensitivity analysis. The sensitivity matrix is ​​calculated using numerical differentiation with a differentiation step size of 0.01 mm. The position coordinate error is obtained by solving the linear equations between the sensitivity matrix and the residual vector using the least squares method. The geometric transformation parameter error is derived from the position coordinate error. The derivation process calculates the partial derivatives of the position coordinates with respect to the geometric transformation parameters, constructs a linear relationship matrix, and solves for the parameter error using a pseudo-inverse method. Geometric transformation parameter correction subtracts the error from the current parameter values. The correction magnitude is constrained: the radial distortion coefficient change cannot exceed 10% of the current value, the rotation angle change cannot exceed 0.05 rad, and the translation vector change cannot exceed 5 mm. The corrected geometric transformation parameters replace the original parameters in the unified parameter matrix and are simultaneously applied to multi-view fusion transformation and temporal correlation. Updated position coordinate sequences and 3D trajectory data are recalculated, and trajectory quality evaluation metrics are recalculated. When the metrics fall below a threshold, optimization is considered converged, forming a closed-loop optimization. Closed-loop optimization is periodically triggered, once every 1000 frames processed, with a maximum iteration limit of 10 times.

[0062] Optionally, establishing a geometric transformation model that includes lens distortion compensation and dynamic baseline correction includes:

[0063] By acquiring a sequence of calibration images containing known spatial location calibration points, the lens distortion parameters and relative position and attitude parameters of multiple image acquisition devices are simultaneously solved. The lens distortion parameters and the relative position and attitude parameters are then integrated into a unified parameter matrix to establish a geometric transformation model that couples lens distortion compensation and baseline correction.

[0064] The system monitors the pixel position changes of calibration points in multiple image acquisition devices in real time. When the detected pixel position change exceeds a preset change threshold, it determines that the relative positions of the multiple image acquisition devices have changed.

[0065] Based on the pixel position change, the reprojection error of the calibration point under different viewpoints is calculated. By minimizing the reprojection error, the lens distortion parameters and relative position and attitude parameters in the unified parameter matrix are jointly optimized and updated to obtain the updated geometric transformation model.

[0066] For example, when acquiring images of the calibration board, the calibration board uses a checkerboard pattern, with the corner points of the checkerboard serving as spatial positioning calibration points. The spacing between the corner points is controlled within 0.01 mm. The calibration board is fixed at a known three-dimensional coordinate position in the electric field environment. The three-dimensional coordinates are obtained by measuring with a laser rangefinder, with a measurement accuracy of no less than 0.1 mm. Multiple image acquisition devices are simultaneously triggered for exposure, with the synchronization error controlled within 1 microsecond. The resolution of the acquired calibration image sequence is set to 1920×1080 pixels, and the frame rate is set to 100 frames per second. The pixel coordinates of the checkerboard corner points in the calibration image sequence are extracted. The corner point detection uses a sub-pixel level corner point positioning algorithm, achieving a positioning accuracy of 0.1 pixels.

[0067] The 3D spatial coordinates of the calibration points are known, and the corresponding pixel coordinates are obtained through corner detection. A projection calculation path from 3D space to the pixel plane is established. The projection process consists of two stages. In the first stage, the 3D spatial coordinates are transformed to the camera coordinate system of the image acquisition device through rotation and translation transformations. The rotation transformation is described by three independent rotation angles, and the translation transformation is described by a 3D translation vector. In the second stage, the 3D coordinates in the camera coordinate system are projected onto the 2D pixel plane. The projection calculation includes intrinsic parameters such as focal length, principal point position, and pixel scaling ratio, and a distortion correction stage is introduced. Radial distortion correction uses second-order and fourth-order distortion coefficients. The correction calculation nonlinearly adjusts the ideal pixel coordinates according to the radial distance from the principal point; the larger the radial distance, the more significant the distortion correction. Tangential distortion correction uses two independent coefficients. The correction calculation superimposes a tangential offset component on top of the radial correction.

[0068] Substituting the known 3D coordinates of all calibration points into the projection calculation path yields the theoretical pixel coordinates for each calibration point. The theoretical pixel coordinates differ from the actual pixel coordinates obtained from corner detection; this difference is the reprojection error. The optimization objective is to adjust parameters such as rotation angle, translation vector, focal length, principal point position, and distortion coefficient to minimize the sum of squared reprojection errors of all calibration points across all image acquisition devices. The optimization algorithm employs a nonlinear least squares iterative solution. Initially, the rotation angle is set to zero degrees, the translation vector is roughly estimated based on the physical installation position of the image acquisition device, the focal length is set according to the lens's nominal focal length, the principal point position is set to the image center, and the initial distortion coefficient is set to zero. During iteration, the total reprojection error under the current parameter combination is calculated, and the partial derivatives of the error with respect to each parameter are calculated. The parameter values ​​are adjusted based on the direction of the partial derivatives. An adaptive strategy is used for parameter adjustment step size. The initial step size is set to one-hundredth of the parameter's magnitude. When the error decreases, the step size increases by a factor of 1.2; when the error increases, the step size decreases by a factor of 0.5, and the parameters are reverted to the previous iteration value. The iteration termination condition is set to the total reprojection error change being less than 0.001 pixels or the number of iterations exceeding 500. After the iteration is completed, the rotation angle and translation vector obtained are the relative position and attitude parameters, and the obtained distortion coefficients are the lens distortion parameters.

[0069] A unified parameter matrix is ​​constructed to store the lens distortion parameters and relative position and pose parameters of all image acquisition devices. Assuming three image acquisition devices are configured, each device corresponds to a set of lens distortion parameters, including two radial distortion coefficients and two tangential distortion coefficients, for a total of four parameters. The relative position and pose parameters are based on the first image acquisition device as the reference frame. The rotation angle and translation vector of the first image acquisition device are set to zero. The second and third image acquisition devices each have three rotation angles and three translation vector components, resulting in six relative pose parameters for each device. There are a total of 12 lens distortion parameters and 12 relative pose parameters across the three image acquisition devices. Adding the intrinsic parameters such as focal length and principal point coordinates for each device, a total of 21 parameters are generated, resulting in a unified parameter matrix with a total dimension of 45. The parameter matrix is ​​stored in row-major order, and the parameter precision is represented as double-precision floating-point numbers, with each parameter occupying 8 bytes of storage space.

[0070] The geometric transformation model performs coordinate transformation based on a unified parameter matrix. The input is any two-dimensional pixel coordinates and the corresponding image acquisition device number; the output is three-dimensional spatial coordinates. The transformation process first extracts the corresponding distortion coefficients based on the image acquisition device number, and then applies the inverse operation of distortion correction to convert the pixel coordinates into distortion-free normalized image coordinates. The normalized image coordinates, combined with focal length and principal point position parameters, are converted into direction vectors in the camera coordinate system. These direction vectors are then combined with the rotation angle and translation vector parameters of the image acquisition device, and transformed into a unified world coordinate system through inverse rotation and inverse translation transformations to obtain the three-dimensional spatial ray direction. Multiple image acquisition devices calculate their own three-dimensional spatial rays for the same charged particle; the intersection of these rays represents the three-dimensional spatial coordinates of the charged particle.

[0071] During actual operation in the electric field environment, monitoring images, including those of the calibration board, are continuously acquired at a frequency of once every 10 seconds to avoid frequent calculations that consume computational resources. The current pixel position of the calibration point in the monitoring image is extracted and compared with the reference pixel position recorded during initial calibration. The change in pixel position is calculated as the square root of the sum of the squares of the differences between the current pixel position and the reference pixel position in the horizontal and vertical directions, with the change measured in pixels. A preset change threshold is set to 0.5 pixels. When the change in pixel position of any calibration point exceeds 0.5 pixels, it is determined that the relative positions between the multiple image acquisition devices have changed, and the geometric transformation parameter update process is initiated.

[0072] The parameter update process uses the current unified parameter matrix to project the known 3D spatial coordinates of the calibration points onto the image planes of each image acquisition device to obtain the theoretical pixel positions. The difference between the theoretical pixel positions and the actual pixel positions detected in the monitored images constitutes the reprojection error vector, which contains two components: horizontal and vertical. The reprojection error vectors of all calibration points and all image acquisition devices are statistically analyzed, and the total reprojection error is calculated as the sum of the squares of the horizontal and vertical components of all error vectors. The optimization algorithm adjusts the lens distortion coefficients and relative pose parameters in the unified parameter matrix, while keeping the 3D spatial coordinates of the calibration points unchanged. The focal length and principal point position in the intrinsic parameters are usually kept unchanged to improve optimization stability. The optimization iteration process calculates the total reprojection error under the current parameters, calculates the gradient of the error with respect to the lens distortion coefficients and relative pose parameters, and updates the parameter values ​​according to the gradient direction. The parameter update step size adopts an adaptive adjustment strategy, with an initial step size set to 0.01. When the total reprojection error decreases, the step size increases by a factor of 1.2; when the total reprojection error increases, the step size decreases by a factor of 0.5, and the parameters are rolled back to the previous iteration value. The iteration termination condition is set to a total reprojection error of less than 0.1 pixels or an iteration count exceeding 300. After optimization, an updated unified parameter matrix is ​​obtained. The updated matrix replaces the original matrix, and subsequent calculations of the 3D coordinates of charged particles are performed using the updated parameters and geometric transformations.

[0073] This method triggers joint optimization and update of geometric transformation parameters by real-time monitoring of changes in the position of calibration point pixels, adapts to changes in the relative position of the image acquisition device, and ensures the long-term accuracy and stability of multi-view fusion transformation.

[0074] Optionally, performing multi-view fusion transformation on the pixel coordinates of charged particles in the two-dimensional image sequence to obtain a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel accuracy includes:

[0075] The geometric transformation model is used to perform distortion correction and viewpoint transformation on the pixel coordinates of charged particles in the two-dimensional image sequence, thereby obtaining the corrected pixel coordinates corresponding to each image acquisition device;

[0076] The electric field direction information of the electric field environment is obtained, the angle between the observation direction of each of the multiple image acquisition devices and the electric field direction is calculated, and a confidence weight is assigned to each image acquisition device according to the angle.

[0077] The corrected pixel coordinates are weighted and fused according to the confidence weight, and the candidate position coordinate set of charged particles in three-dimensional space is solved by multi-view geometric constraints;

[0078] Based on the physical constraints of the motion of charged particles under the influence of electric field forces in an electric field environment, the consistency metric between the direction of the electric field force and the direction of motion of the charged particles corresponding to each candidate position coordinate in the candidate position coordinate set is calculated. Candidate position coordinates with consistency metrics lower than the consistency threshold are eliminated, and candidate position coordinates that satisfy the force equilibrium condition are retained as the optimal position coordinates.

[0079] The optimal position coordinates are subjected to temporal interpolation and sub-pixel refinement to obtain a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel precision.

[0080] For example, pixel coordinates of charged particles are extracted from a two-dimensional image sequence. Pixel coordinates include two components: horizontal and vertical coordinates. The origin of the coordinates is located at the upper left corner of the image, with the horizontal direction to the right being positive and the vertical direction downwards being positive. The geometric transformation model receives the pixel coordinates and the corresponding image acquisition device number as input, and extracts the corresponding lens distortion parameters from the unified parameter matrix based on the image acquisition device number. Distortion correction calculation converts the original pixel coordinates into normalized image coordinates. The normalized image coordinates have the principal point position as the origin, and the coordinate values ​​are divided by the focal length to obtain dimensionless coordinates. Distortion correction uses an inverse iterative solution. The initial normalized coordinates are set as the original pixel coordinates minus the principal point position divided by the focal length. The radial distance corresponding to this normalized coordinate is calculated as the square root of the sum of the squares of the horizontal and vertical coordinates. The radial distortion correction is calculated as radial distance 2 × second-order radial distortion coefficient + radial distance 4 × fourth-order radial distortion coefficient. The corrected radial coordinates are the original radial coordinates × (distortion correction + 1). The tangential distortion correction is calculated separately in the horizontal and vertical directions. The horizontal correction is calculated as: first tangential distortion coefficient × radial distance 2 × 2 + second tangential distortion coefficient × radial distance 2 × horizontal coordinate. The vertical correction is calculated as: second tangential distortion coefficient × radial distance 2 × 2 + first tangential distortion coefficient × radial distance 2 × vertical coordinate. During the iteration process, the corrected normalized coordinates are substituted into the distortion calculation to derive the original normalized coordinates in reverse. The iteration terminates when the difference between the normalized coordinates of two consecutive iterations is less than 0.0001 or the number of iterations exceeds 10. After iteration, the distortion-free normalized image coordinates are obtained. The corrected pixel coordinates are obtained by multiplying the normalized image coordinates by the focal length and adding the principal point position.

[0081] The electric field strength distribution information of the electric field environment is obtained in advance through electric field simulation calculations or actual measurements and stored as three-dimensional spatial grid data with a grid resolution of 10 mm. Each grid node records the field strength vector at that location, and the field strength vector contains three directional components. The field strength vector of the nearest grid node is queried based on the coarse position estimate of the charged particles; the field strength direction is the normalized direction vector of the field strength vector. The observation direction of the image acquisition device is the optical axis direction of the device. The optical axis direction is stored in the relative position and attitude parameters of the unified parameter matrix as a rotation matrix, and the third column vector of the rotation matrix is ​​the optical axis direction vector. The included angle is calculated as the dot product of the observation direction vector and the field strength direction vector; the inverse cosine of the dot product is the included angle, with the unit being radians and a value ranging from 0 to π.

[0082] The confidence weight allocation determines the angle weight coefficient based on the included angle value, with a preset vertical angle threshold of π / 2, or 90 degrees. When the included angle is greater than π / 2, the angle between the observation direction and the field strength direction is relatively large, and the angle weight coefficient is calculated as (including angle - π / 2) ÷ (π / 2), with a value range of 0 to 1; the larger the included angle, the greater the weight coefficient. When the included angle is less than or equal to π / 2, the angle between the observation direction and the field strength direction is relatively small, and the angle weight coefficient is calculated as (π / 2 - included angle) ÷ (π / 2), with a value range of 0 to 1; the smaller the included angle, the smaller the weight coefficient. This weight allocation strategy ensures that the highest weight is obtained when the observation direction is nearly perpendicular to the field strength direction, because at this time, the projected displacement of the charged particle moving along the field strength direction on the image plane is the largest, resulting in the highest positioning accuracy. Simultaneously, image quality parameters of the two-dimensional images acquired by each image acquisition device are obtained. Image sharpness is evaluated by calculating the statistical measure of the image gradient magnitude. The gradient magnitude statistic uses the peak position of the gradient magnitude histogram; the higher the peak position, the sharper the image. Signal-to-noise ratio (SNR) is evaluated by calculating the ratio of image signal power to noise power. Signal power is the variance of image pixel grayscale values, and noise power is the variance of pixel grayscale values ​​in dark areas of the image. An initial confidence weight is obtained by combining the angle weight coefficient with image quality parameters. The combination operation is: angle weight coefficient × normalized image sharpness value × normalized SNR value. Normalization maps each parameter value to the range of 0 to 1. The initial confidence weight is summed over all image acquisition devices and then normalized. The normalization process divides the initial confidence weight of each image acquisition device by the sum of the initial confidence weights of all image acquisition devices, resulting in a normalized sum of confidence weights of 1.

[0083] The weighted fusion process uses the corrected pixel coordinates to calculate the 3D spatial rays corresponding to each image acquisition device. The ray origin is the optical center position of the image acquisition device, and the ray direction is the 3D direction vector obtained by transforming the normalized image coordinates corresponding to the corrected pixel coordinates using a rotation matrix. When solving for the 3D spatial position of charged particles using multi-view geometric constraints, the closest point between all ray pairs is calculated. The closest point between two rays is obtained by solving for the parameter value that minimizes the distance between the two rays. A ray is represented as origin + parameter × direction vector. The square of the distance between two rays is the difference vector between the two origins + the square of the magnitude of the difference vector between the two parameters × the corresponding direction vectors. The partial derivatives of the square of the distance with respect to the two parameters are calculated, and the partial derivatives are set to 0 to obtain the two parameter values. Substituting these parameter values ​​into the ray expression yields the closest point between the two rays. The candidate position coordinate set contains the closest points of all ray pairs. If three image acquisition devices are configured, the number of ray pairs is 3, and the number of candidate position coordinates is 3. The weighted fusion calculation is to calculate the weighted average of all candidate position coordinates according to the confidence weights of the two image acquisition devices involved in the corresponding ray pair. The weight of the weighted average is the product of the confidence weights of the two image acquisition devices. The weighted sum of all candidate position coordinates divided by the total weights gives the fused three-dimensional spatial position.

[0084] The direction of motion of charged particles is calculated by the difference between their position coordinates at adjacent time points. If the current time is frame t and the previous time was frame (t-1), the motion direction vector is the position coordinate of frame t minus the position coordinate of frame (t-1), and the motion direction is normalized to a unit vector. The direction of the electric force is determined based on the local field strength vector at the location of the charged particle. The direction of the electric force for a positive charge is the same as the field strength direction, while the direction of the electric force for a negative charge is opposite to the field strength direction. The consistency metric is calculated as the dot product of the electric force direction vector and the motion direction vector. The dot product ranges from -1 to +1; the closer the dot product is to +1, the more consistent the electric force direction and the motion direction are. The consistency threshold is set to 0.5. When the consistency metric corresponding to a candidate position coordinate is lower than 0.5, the candidate position coordinate is determined not to meet the force equilibrium condition and is discarded. The retained candidate position coordinates are used as the optimal position coordinates. If all candidate position coordinates are discarded, the consistency threshold is lowered to 0.3 for re-selection. If no coordinates are still retained, the fused 3D spatial position is used as the optimal position coordinate.

[0085] Temporal interpolation addresses uneven inter-frame time intervals or partial frame loss by detecting the timestamp difference between adjacent optimal position coordinates. A lost frame is identified when the timestamp difference exceeds twice the standard frame interval. The interpolation calculation employs a linear interpolation strategy: the interpolated position coordinate is the previous optimal position coordinate plus a time ratio multiplied by the position coordinate difference. The time ratio is calculated as (interpolation timestamp - previous timestamp) ÷ (difference between the two optimal position coordinate timestamps). Sub-pixel-level refinement fine-tunes the optimal position coordinates through iterative optimization. The optimization objective is to minimize the reprojection error between the optimal position coordinates projected onto the image plane of each image acquisition device and the detected charged particle pixel coordinates. The reprojection error is calculated as the Euclidean distance between the theoretical and actual pixel coordinates. The theoretical pixel coordinates are obtained by projecting the optimal position coordinates onto the image plane using a geometric transformation model. The optimization algorithm adjusts the three spatial components of the optimal position coordinates. The initial iteration step size is set to 0.01 mm, increasing by 1.1 when the reprojection error decreases and by 0.5 when the reprojection error increases. The iteration terminates when the reprojection error is less than 0.05 pixels or the number of iterations exceeds 20. After thinning, position coordinates with sub-pixel accuracy are obtained, and the position coordinates at all times are arranged in timestamp order to form a position coordinate sequence.

[0086] This method eliminates the nonlinear effects of the lens through distortion correction, optimizes the multi-view fusion accuracy based on the confidence weight of field strength direction and image quality, eliminates unreasonable candidate positions by combining motion physics constraints, and improves the spatiotemporal continuity and accuracy of position coordinate sequence through temporal interpolation and sub-pixel level refinement.

[0087] Optionally, the steps of acquiring the electric field direction information of the electric field environment, calculating the angle between the observation direction of each of the multiple image acquisition devices and the electric field direction, and assigning confidence weights to each image acquisition device according to the angles include:

[0088] The local field strength vector at the location of the charged particles in the electric field environment is obtained, and the direction of the local field strength vector is extracted as the field strength direction information; the angle between the optical axis direction of the multiple image acquisition devices and the field strength direction is calculated respectively, and the optical axis direction is used as the observation direction of the corresponding image acquisition device;

[0089] Image quality parameters of the two-dimensional images acquired by each image acquisition device are obtained, including image sharpness and signal-to-noise ratio; an angle weighting coefficient is determined based on the numerical range of the included angle. When the included angle is greater than a preset vertical angle threshold, the angle weighting coefficient increases monotonically as the included angle increases; when the included angle is less than or equal to the preset vertical angle threshold, the angle weighting coefficient decreases monotonically as the included angle decreases.

[0090] The angle weight coefficient is combined with the image quality parameter to obtain the initial confidence weight; the initial confidence weight is then normalized to obtain the final confidence weight.

[0091] For example, the electric field strength distribution information of the electric field environment is pre-calculated using electric field simulation software or obtained through actual measurement and stored as a three-dimensional spatial grid data structure. The grid data structure divides the electric field environment into regular cubic grids with a grid resolution of 10 mm × 10 mm × 10 mm. Each grid node stores the field strength vector at that location. The field strength vector includes X-direction, Y-direction, and Z-direction components, with component values ​​in volts per meter and numerical precision expressed as single-precision floating-point numbers. The field strength distribution data is queried based on the three-dimensional coordinates of the charged particle's location. The query process uses a nearest-neighbor interpolation strategy, calculating the Euclidean distance between the charged particle's location coordinates and the coordinates of each grid node, and selecting the grid node with the smallest distance as the query result. The queried field strength vector is the local field strength vector at the location of the charged particle, and the three components of the local field strength vector are denoted as FX, FY, and FZ.

[0092] Field strength direction information is extracted by normalizing the local field strength vector. The normalization process calculates the magnitude of the field strength vector as the square root of (FX² + FY² + FZ²). The three components of the field strength direction are FX ÷ magnitude, FY ÷ magnitude, and FZ ÷ magnitude, respectively. The normalized field strength direction vector has a magnitude of 1, and its direction is consistent with the original field strength vector. When the magnitude of the field strength vector is less than 0.1 volts per meter, the field strength at that location is considered close to zero, and the field strength direction information is set to the default value (0, 0, -1), with the default direction pointing vertically downwards.

[0093] The optical axis direction of the image acquisition device is extracted from the unified parameter matrix of the geometric transformation model. This unified parameter matrix stores the relative position and attitude parameters of each image acquisition device. The relative position and attitude parameters include a rotation matrix and a translation vector. The rotation matrix is ​​a 3×3 matrix, and its elements describe the rotational transformation relationship between the image acquisition device's coordinate system and the world coordinate system. The third column vector of the rotation matrix represents the direction of the Z-axis of the image acquisition device's coordinate system in the world coordinate system; this direction is the optical axis direction. The optical axis direction vector contains three components, denoted as OX, OY, and OZ, with a magnitude of 1. By iterating through all image acquisition devices, the optical axis direction vector corresponding to each device is extracted as the observation direction for that device.

[0094] The angle is calculated using a vector dot product operation. The dot product of the optical axis direction vector and the field intensity direction vector is equal to OX×FX÷magnitude + OY×FY÷magnitude + OZ×FZ÷magnitude. The dot product value ranges from -1 to +1, and the inverse cosine function value of the dot product is the angle, with the unit being radians and a value range from 0 to π. The inverse cosine function is approximated using a lookup table or Taylor series, with the calculation accuracy controlled within 0.001 radians. The angle between the optical axis direction and the field intensity direction of each image acquisition device is calculated, and the angle value is stored as a floating-point number.

[0095] Image quality parameters are extracted from two-dimensional images acquired by various image acquisition devices. These images are in grayscale format, with pixel values ​​ranging from 0 to 255. Image sharpness is evaluated by calculating the image gradient magnitude using the Sobel operator, which performs convolution operations on the image horizontally and vertically. The horizontal gradient convolution kernel is a 3×3 matrix, with the center column being the transpose of (+2, 0, -2), the left column the transpose of (+1, 0, -1), and the right column the transpose of (+1, 0, -1). The vertical gradient convolution kernel is the transpose of the horizontal kernel. The gradient magnitude is calculated as the square root of the sum of the horizontal and vertical gradient values, and the image sharpness is the average of the gradient magnitudes of all pixels. Signal-to-noise ratio (SNR) calculation requires separating image signal power from noise power. Signal power is represented by the variance of the image pixel grayscale values, calculated as the sum of the squares of the differences between all pixel grayscale values ​​and the average grayscale value, divided by the total number of pixels. Noise power is represented by the variance of pixel gray values ​​in areas where the gray value is 20% lower than the average gray value in the image. This area is typically a dark background region, and gray value fluctuations are mainly caused by noise. Signal-to-noise ratio (SNR) is calculated as signal power ÷ noise power; a higher SNR indicates better image quality. Image sharpness and SNR are both normalized, mapping the values ​​to the range of 0 to 1 using the formula (current value - minimum value) ÷ (maximum value - minimum value). The minimum and maximum values ​​are determined by statistically analyzing historical image data. The minimum image sharpness is set to 0, and the maximum to 50; the minimum SNR is set to 1, and the maximum to 20. When the actual value exceeds the set range, the normalization result is truncated to 0 or 1.

[0096] The angle weighting coefficient is calculated segmentally based on the range of included angle values, with a preset vertical angle threshold of π / 2 radians (90 degrees). The included angle value is compared to the vertical angle threshold. When the included angle is greater than π / 2, the angle between the observation direction and the field strength direction is too large, and the angle weighting coefficient is calculated using a linearly increasing function: (including angle - π / 2) ÷ (π / 2). The angle weighting coefficient is 0 when the included angle is π / 2 and 1 when the included angle is π, with the coefficient monotonically increasing within the range of π / 2 to π. When the included angle is less than or equal to π / 2, the angle between the observation direction and the field strength direction is too small, and the angle weighting coefficient is calculated using a linearly decreasing function: (π / 2 - included angle) ÷ (π / 2). The angle weighting coefficient is 1 when the included angle is 0 and 0 when the included angle is π / 2, with the coefficient monotonically decreasing within the range of 0 to π / 2. This weighting strategy ensures that the angle weighting coefficient is 0 when the observation direction is perpendicular to the field strength direction, and 1 when the observation direction is parallel or antiparallel to the field strength direction. The physical basis for this is that when a charged particle moves along the field strength direction, an image acquisition device with its observation direction perpendicular to the field strength direction can completely observe the displacement, resulting in the highest positioning accuracy. Conversely, an image acquisition device with its observation direction parallel to the field strength direction only observes the projected component of the motion, resulting in lower positioning accuracy.

[0097] The combination operation multiplies the angle weight coefficient by the image quality parameter to obtain the initial confidence weight. The calculation formula is: angle weight coefficient × normalized image sharpness value × normalized signal-to-noise ratio value. The initial confidence weight ranges from 0 to 1; a larger value indicates higher reliability of the image acquisition device in calculating the 3D position of charged particles. The initial confidence weights of all image acquisition devices are calculated, and all initial confidence weights are summed to obtain the total weight. Normalization is performed by dividing the initial confidence weight of each image acquisition device by the total weight; the sum of the normalized confidence weights is 1. When the initial confidence weight of an image acquisition device is 0, the normalized confidence weight is also 0, and that image acquisition device does not participate in subsequent multi-view fusion calculations. When the initial confidence weights of all image acquisition devices are 0, the confidence weights of all image acquisition devices are set to an equal value of 1 ÷ the number of image acquisition devices to ensure that multi-view fusion calculations can continue.

[0098] This method assigns confidence weights to each image acquisition device based on the angle between the observation direction and the field strength direction, as well as image quality parameters. This ensures that image acquisition devices with observation directions that are close to or opposite to the field strength direction and with high image quality receive higher weights, thereby improving the positioning accuracy and robustness of multi-view fusion.

[0099] Optionally, the step of temporally associating the three-dimensional spatial position coordinates at different times based on the timestamp information and spatial continuity constraints in the position coordinate sequence to obtain the three-dimensional trajectory data of the charged particle includes:

[0100] Extract timestamp information from the location coordinate sequence and sort them according to the timestamp order to obtain a time-series arrangement of location coordinates;

[0101] Calculate the spatial distance and velocity change between adjacent position coordinates in the time-series position coordinate arrangement. When the spatial distance exceeds the position deviation threshold calculated based on the motion velocity or the velocity change exceeds the acceleration threshold calculated based on the direction of the electric field force, it is determined that there is a trajectory break.

[0102] For coordinates where there is a trajectory break, based on the velocity direction and electric force direction of the position coordinates before the break, the interpolated position coordinates between the break times are calculated, and the interpolated position coordinates are inserted into the time-series of position coordinates.

[0103] The similarity in spatial location distribution and motion direction between different trajectory segments in the repaired temporal arrangement position coordinate sequence is calculated. When the similarity is higher than the association threshold, the trajectory segments belonging to the same charged particle are spliced ​​together in the order of timestamps. The spliced ​​position coordinate sequence is then smoothed and filtered to obtain the three-dimensional trajectory data of the charged particle.

[0104] For example, each position coordinate in the position coordinate sequence contains a three-dimensional spatial coordinate and a corresponding timestamp. The three-dimensional spatial coordinate consists of X, Y, and Z coordinates, with units of millimeters. The timestamp records the moment the position coordinate was acquired, with units of milliseconds and a precision of 1 millisecond. The position coordinate sequence is traversed to extract all timestamp information, which is stored as a 64-bit integer to avoid loss of floating-point precision. A quicksort algorithm is used, with the comparison function arranging the coordinates according to their timestamp values ​​from smallest to largest. Position coordinates with the same timestamp maintain their original order. After sorting, a temporally ordered position coordinate sequence is obtained, where the position coordinates are arranged sequentially according to time, and the sequence index corresponds to a monotonically increasing timestamp.

[0105] Spatial distance calculation is performed on two adjacent position coordinates in a temporally arranged sequence of position coordinates. Adjacent position coordinates refer to two position coordinates whose sequence indices differ by 1. The spatial distance is calculated as the Euclidean distance between the three-dimensional coordinate differences of the two position coordinates, and the formula is: √((X2-X1)²+(Y2-Y1)²+(Z2-Z1)²), where X1, Y1, and Z1 are the three-dimensional components of the previous position coordinate, and X2, Y2, and Z2 are the three-dimensional components of the next position coordinate. Velocity change calculation requires first calculating the instantaneous velocity corresponding to the adjacent position coordinates. The instantaneous velocity is the position coordinate difference divided by the timestamp difference. The velocity vector includes velocity components in the X, Y, and Z directions. The velocity component calculation formula is (X2-X1)÷(T2-T1), where T1 and T2 are the corresponding timestamps in milliseconds, and the velocity unit is millimeters per millisecond. The change in velocity is the magnitude of the difference between adjacent instantaneous velocity vectors. The formula is: take the square root of ((VX2-VX1)2+(VY2-VY1)2+(VZ2-VZ1)2), where VX1, VY1, and VZ1 are the velocity components of the previous segment, and VX2, VY2, and VZ2 are the velocity components of the next segment.

[0106] The position deviation threshold is calculated based on the velocity. The calculation method is to multiply the velocity magnitude of the previous segment by the timestamp difference, and then multiply by the deviation coefficient. The velocity magnitude is calculated as the square root of (VX12 + VY12 + VZ12). The deviation coefficient is set to 1.5, which physically means that the allowable position deviation is 1.5 times the theoretical displacement. This coefficient takes into account the accelerated motion of charged particles under the influence of the electric field force. When the actual spatial distance exceeds the position deviation threshold, a trajectory break is determined. The acceleration threshold is calculated based on the direction of the electric field force. The calculation method is to multiply the field strength magnitude in the direction of the electric field force by the unit charge, divide by the unit mass, and then multiply by the acceleration coefficient. The field strength magnitude is obtained from the electric field distribution information of the electric field environment. The unit charge is set to 1.6 × 10⁻¹⁹ coulombs, the unit mass is set to 1 × 10⁻¹⁵ kg, and the acceleration coefficient is set to 2.0, which physically means that the allowable acceleration is twice the theoretical electric field acceleration. When the actual velocity change divided by the timestamp difference exceeds the acceleration threshold, a trajectory break is determined. The system iterates through the temporally ordered sequence of position coordinates, calculating the spatial distance and velocity change for each pair of adjacent coordinates. If any breakage condition is met, the system marks that a trajectory break exists between the coordinate pairs.

[0107] The interpolation position coordinate calculation is performed on adjacent position coordinate pairs where a trajectory break exists. The position coordinates before the break refer to the earlier timestamp, and the position coordinates after the break refer to the later timestamp. The velocity direction is extracted by finding the velocity vector between the pre-break position coordinate and its preceding position coordinate, and normalized to a unit vector. The direction of the electric force is obtained by querying the local field strength vector at the pre-break position coordinate from the electric field strength distribution information of the electric field environment, and normalized to a unit vector. The timestamps of the interpolated position coordinates are evenly distributed between the pre-break and post-break timestamps, with the interval set as the standard frame interval. The standard frame interval is obtained by statistically analyzing the median difference between all adjacent timestamps in the time-series of position coordinates. The three-dimensional coordinate calculation of the interpolated position coordinates uses a kinematic model, assuming that the charged particle undergoes uniformly accelerated motion under the influence of a constant electric field force between the break and the time of breakage. The initial position is the pre-break position coordinate, the initial velocity is the pre-break velocity vector, and the acceleration is the electric field force direction multiplied by the acceleration amplitude. The acceleration amplitude is calculated as field strength magnitude × unit charge ÷ unit mass. The time interval between the interpolation time t and the time before the fracture is denoted as Δt. The formula for calculating the interpolation position coordinates is: initial position + initial velocity × Δt + 0.5 × acceleration × Δt². This formula is applied to the X, Y, and Z components respectively. The interpolation position coordinates corresponding to all interpolation times are calculated, and the interpolation position coordinates are inserted into the temporally arranged position coordinate sequence according to the timestamp order. The insertion operation maintains the sequence in a monotonically increasing order according to the timestamp.

[0108] The repaired temporal sequence of position coordinates contains multiple trajectory segments, separated by trajectory break points. The start and end coordinates of all trajectory segments are identified. The start position of a trajectory segment is either the beginning of the sequence or the coordinate after a break, while the end position is either the coordinate before a break or the end of the sequence. Similarity calculations are performed between different trajectory segments, comprising two components: spatial distribution similarity and motion direction similarity. Spatial distribution similarity is calculated as the reciprocal of the spatial distance between the end and start coordinates of two trajectory segments; the smaller the distance, the higher the similarity. Motion direction similarity is calculated as the dot product of the velocity directions at the end and start of two trajectory segments, with a dot product value ranging from -1 to +1; the closer the dot product value is to +1, the more consistent the motion directions. The overall similarity is calculated as spatial distribution similarity × weighting coefficient + motion direction similarity × weighting coefficient, with the two weighting coefficients set to 0.4 and 0.6 respectively. The sum of similarities is normalized to the range of 0 to 1. The association threshold is set to 0.7; when the overall similarity is higher than 0.7, the two trajectory segments are considered to belong to the same charged particle. Traverse all trajectory segment pairs, calculate the comprehensive similarity and compare it with the association threshold. Trajectory segments that meet the condition are marked as associated trajectory segments. The splicing operation merges the associated trajectory segments into a single trajectory segment according to the timestamp order. During splicing, all position coordinates of the two trajectory segments are preserved, and the merged trajectory segment is rearranged according to the timestamp order.

[0109] The smoothing filtering process employs a moving average filtering algorithm, with a filtering window length set to 5 position coordinates. The window slides across the concatenated position coordinate sequence, with a sliding step of 1 position coordinate. The filtering calculation averages the three-dimensional coordinates of the 5 position coordinates within the window. The average X-coordinate is calculated as the sum of the 5 X-coordinates within the window divided by 5; the average Y-coordinate and Z-coordinate are calculated using the same method. The filtered position coordinates replace the window center coordinates, with the window center being the 3rd position coordinate within the window. Since the start and end positions of the sequence cannot form a complete window, the processing strategy is to keep the original values ​​of the first two and last two position coordinates without filtering. After filtering, a smoothed position coordinate sequence is obtained, which represents the three-dimensional trajectory data of the charged particle. This three-dimensional trajectory data contains position coordinates and timestamp information at all times, maintaining continuity in both time and space.

[0110] This method establishes a temporal correlation basis by sorting by timestamps, determines trajectory breaks based on spatial continuity and physical constraints of motion and performs interpolation repair, associates different trajectory segments through similarity calculation, and obtains spatiotemporally continuous three-dimensional trajectory data by combining smoothing filtering.

[0111] Optionally, the step of calculating the trajectory quality evaluation index based on the motion state parameters in the three-dimensional trajectory data and the field strength distribution information of the electric field environment through force balance relationships includes:

[0112] The position coordinates and timestamp information of charged particles at different times are extracted from the three-dimensional trajectory data. The instantaneous velocity of the charged particles is obtained by calculating the difference between the position coordinates of adjacent times, and the instantaneous acceleration is also calculated.

[0113] Obtain the field strength distribution information of the electric field environment, and query the local field strength vector corresponding to the position of the charged particle based on the position coordinate sequence;

[0114] Based on the force balance relationship of charged particles under the action of electric force, gravity and air resistance in an electric field environment, the theoretical acceleration of charged particles is calculated according to the local field strength vector.

[0115] The deviation between the instantaneous acceleration and the theoretical acceleration is calculated as a trajectory quality evaluation index.

[0116] For example, the 3D trajectory data is stored as a time-series data structure. Each data point contains position coordinates and a timestamp. The position coordinates consist of X, Y, and Z coordinates, with units of millimeters, and the timestamps are in milliseconds. The 3D trajectory data is traversed to extract all position coordinates and their corresponding timestamps. The position coordinates are stored as floating-point arrays, and the timestamps are stored as integer arrays. The array indices are consistent with the time order of the trajectory data points. Instantaneous velocity calculation involves a difference operation on the position coordinates at adjacent moments, where adjacent moments refer to two data points whose timestamp indices differ by 1. The X-direction component of the velocity vector is calculated as (X2-X1)÷(T2-T1), where X1 and X2 are the X-components of adjacent position coordinates, and T1 and T2 are the corresponding timestamps. The velocity unit is millimeters per millisecond. The Y-direction and Z-direction components are calculated using the same method as the X-direction component. The velocity vector contains three directional components stored as floating-point arrays. The timestamp of the instantaneous velocity data point is set as the average of the timestamps of two adjacent position coordinates, calculated as (T1+T2)÷2. This timestamp indicates that the moment corresponding to the instantaneous velocity is located between the two position coordinates.

[0117] Instantaneous acceleration is calculated by performing a difference operation on the instantaneous velocities at adjacent moments. Adjacent instantaneous velocities refer to two velocity vectors whose velocity data point indices differ by 1. The X-direction component of the acceleration vector is calculated as (VX2 - VX1) ÷ (TV2 - TV1), where VX1 and VX2 are the X-direction components of adjacent instantaneous velocities, and TV1 and TV2 are the timestamps of the corresponding velocity data points. The acceleration unit is millimeters per millisecond squared. The Y-direction and Z-direction components are calculated in the same way as the X-direction component. The acceleration vector contains three directional components and is stored as a floating-point array. The timestamp of the instantaneous acceleration data point is set as the average of the timestamps of two adjacent velocity data points, calculated as (TV1 + TV2) ÷ 2. Since the velocity data point timestamp is already the average of the position coordinate timestamps, the acceleration data point timestamp corresponds to the vicinity of the middle moment of the original position coordinate sequence. Instantaneous velocities are calculated for all adjacent position coordinates in the 3D trajectory data, and instantaneous accelerations are calculated for all adjacent instantaneous velocities, resulting in a complete instantaneous acceleration sequence.

[0118] The electric field distribution information of the electric field environment is stored in the form of three-dimensional spatial grid data. The grid node spacing is set to 10 mm. Each grid node records the field strength vector at that location. The field strength vector includes X-direction, Y-direction, and Z-direction components, and the unit is volts per meter. When querying the local field strength vector based on the location coordinate sequence, the Euclidean distance between the location coordinate and each grid node is calculated, and the grid node with the smallest distance is selected as the query result. The distance is calculated as the square root of ((X-XN)² + (Y-YN)² + (Z-ZN)²), where X, Y, and Z are the three components of the location coordinate, and XN, YN, and ZN are the three components of the grid node coordinate. The field strength vector obtained by querying is the local field strength vector corresponding to the location of the charged particle. When the location coordinate is between two grid nodes, a trilinear interpolation strategy is used to calculate the interpolated field strength vector. The interpolation process performs linear interpolation in the X, Y, and Z directions respectively, and the interpolation weight is determined according to the ratio of the distance from the location coordinate to the adjacent grid node. Trilinear interpolation requires querying the field intensity vectors of the eight grid nodes surrounding the location coordinates. The accuracy of the interpolated field intensity vector is higher than that obtained by the nearest neighbor query method. For each location coordinate in the location coordinate sequence, the corresponding local field intensity vector is obtained by querying or interpolating. The local field intensity vector sequence has the same index as the location coordinate sequence.

[0119] The theoretical acceleration is calculated based on the force equilibrium relationship of charged particles in an electric field environment. The forces acting on the charged particles include three components: electric force, gravity, and air resistance. The electric force is calculated as the local field strength vector multiplied by the charge of the charged particle, which is set to 1.6 × 10⁻¹⁹ coulombs, corresponding to a single elementary charge. The X-direction component of the electric force vector is the electric field strength vector X component multiplied by the charge; the Y-direction components are calculated using the same method as the Z-direction components. Gravity has only a component in the Z-direction. The magnitude of gravity is the charged particle's mass multiplied by gravitational acceleration, where the charged particle's mass is set to 1 × 10⁻¹⁵ kg, and the gravitational acceleration is taken as 9.8 m / s² (0.0098 mm / ms²). Gravity is downward, and the Z-direction component is negative. The gravity vector is (0, 0, -mass × 0.0098). Air resistance is calculated using the Stokes drag model, where the magnitude of the drag is proportional to the velocity, with a proportionality constant of 6 × π × aerodynamic viscosity × particle radius. The aerodynamic viscosity is set to 1.8 × 10⁻⁵ Pascals per second, and the radius of charged particles is set to 1 micrometer (0.001 millimeters). The drag coefficient is calculated as 6 × 3.14159 × 1.8 × 10⁻⁵ × 0.001 = 3.4 × 10⁻⁷ Newtons per second per meter. The drag vector is in the opposite direction to the velocity vector. The X-component of the drag vector is -drag coefficient × instantaneous velocity X-component, while the Y-component is calculated in the same way as the Z-component.

[0120] The resultant force vector is the vector sum of the electric force vector, the gravity vector, and the drag vector. The X-direction component of the resultant force is the sum of the X-components of the electric force, gravity, and drag. The Y-direction components are calculated using the same method as the Z-direction components. The theoretical acceleration is calculated according to Newton's second law. The theoretical acceleration vector is the resultant force vector divided by the mass of the charged particle. The X-direction component of the theoretical acceleration is the sum of the X-components of the resultant force divided by the mass. The Y-direction components are calculated using the same method as the Z-direction components. The unit of theoretical acceleration is millimeters per millisecond squared. For each position coordinate, the local field strength vector and instantaneous velocity are used to calculate the electric force, gravity, and air drag. The sum of these is divided by the mass to obtain the theoretical acceleration. The theoretical acceleration sequence has the same index as the position coordinate sequence.

[0121] The deviation is calculated as the difference between the instantaneous acceleration and the theoretical acceleration. The X-direction component of the deviation vector is the difference between the instantaneous acceleration X-component and the theoretical acceleration X-component. The Y-direction components are calculated using the same method as the Z-direction components. The magnitude of the deviation vector is calculated as the square root of (2x² of deviation X-component + 2x² of deviation Y-component + 2x² of deviation Z-component). The magnitude value characterizes the overall deviation between the instantaneous acceleration and the theoretical acceleration. The trajectory quality assessment index uses a statistical measure of the deviation vector magnitude, which includes both root mean square (RMS) and maximum values. The RMS value is calculated by taking the square root of the sum of the squares of all deviation vector magnitudes divided by the number of data points. The formula is the square root of (Σ(2x² of deviation magnitude) ÷ N), where N is the number of data points, and Σ represents the summation over all data points. The maximum value is the largest value among all deviation vector magnitudes. The RMS value reflects the overall deviation level of the trajectory, while the maximum value reflects the largest local deviation within the trajectory. Calculate the deviation vector sequence between the instantaneous acceleration sequence and the theoretical acceleration sequence, and statistically analyze the root mean square value and maximum value of the deviation vector magnitude. Use the root mean square value as the main trajectory quality evaluation index and the maximum value as the auxiliary evaluation index.

[0122] This method extracts motion state parameters through differential calculation, calculates theoretical acceleration by combining the force balance relationship between electric field force, gravity and air resistance, and evaluates trajectory quality by the statistical deviation between instantaneous acceleration and theoretical acceleration, quantitatively reflecting the physical consistency and measurement accuracy of three-dimensional trajectory data.

[0123] Optionally, the steps of deriving the error of the geometric transformation parameters in the geometric transformation model based on the trajectory quality evaluation index and the physical constraint relationship, correcting the geometric transformation parameters, and synchronously applying the corrected geometric transformation parameters to the multi-view fusion transformation and the temporal correlation to update the position coordinate sequence, forming a closed-loop image processing optimization mechanism, include:

[0124] When the trajectory quality evaluation index exceeds the preset quality threshold, it is determined that there is an error in the geometric transformation parameters in the geometric transformation model;

[0125] A physical constraint equation is established based on the force equilibrium relationship of charged particles in an electric field environment. The physical constraint equation describes the relationship between the motion state parameters and position coordinates of charged particles.

[0126] Substitute the trajectory quality evaluation index into the physical constraint equation, calculate the error of the position coordinate sequence by reverse derivation, and derive the error of the geometric transformation parameters in the geometric transformation model based on the error of the position coordinate sequence.

[0127] The geometric transformation parameters are corrected based on the geometric transformation parameters to obtain the corrected geometric transformation parameters;

[0128] The corrected geometric transformation parameters are synchronously applied to the multi-view fusion transformation of the charged particle pixel coordinates in the two-dimensional image sequence, and the updated position coordinate sequence is recalculated. The updated position coordinate sequence is then applied to the temporal correlation to obtain the updated three-dimensional trajectory data.

[0129] Based on the updated 3D trajectory data, the trajectory quality evaluation index is recalculated to form a closed-loop image processing optimization mechanism.

[0130] For example, the trajectory quality assessment index is numerically compared with a preset quality threshold, which is set to 0.5 mm / s². This threshold is determined based on the typical motion characteristics of charged particles in an electric field environment and the measurement accuracy of image acquisition. The trajectory quality assessment index is represented by the root mean square value of the deviation between instantaneous acceleration and theoretical acceleration. When the root mean square value exceeds 0.5 mm / s², it is determined that there is an error in the geometric transformation parameters in the geometric transformation model, triggering the parameter correction process. The geometric transformation parameters include lens distortion parameters and relative position and attitude parameters. Lens distortion parameters include radial distortion coefficients and tangential distortion coefficients, and relative position and attitude parameters include rotation matrices and translation vectors. Sources of parameter errors include mechanical vibration of the image acquisition device, lens deformation caused by temperature changes, and installation position drift after long-term operation.

[0131] The physical constraint equations are based on the force equilibrium relationship of a charged particle in an electric field environment. This force equilibrium relationship describes the vector sum of the electric force, gravity, and air resistance acting on the charged particle as equal to mass multiplied by acceleration. The physical constraint equations relate instantaneous acceleration, instantaneous velocity, position coordinates, and the local field strength vector. The equation form is: Mass × Instantaneous Acceleration = Charge × Local Field Strength Vector + Mass × Gravitational Acceleration Vector - Drag Coefficient × Instantaneous Velocity. This equation establishes independent constraint relationships in the X, Y, and Z directions. The constraint equation in the X direction is: Mass × Instantaneous Acceleration X Component = Charge × Field Strength X Component - Drag Coefficient × Velocity X Component. The constraint equations in the Y and Z directions have the same form, with the Z direction additionally including a gravity term. The instantaneous acceleration and instantaneous velocity in the physical constraint equations are calculated using the first and second time derivatives of the position coordinates. Therefore, errors in the position coordinates will propagate to the instantaneous acceleration and instantaneous velocity, thus affecting the balance of the physical constraint equations.

[0132] The process of substituting the trajectory quality assessment index into the physical constraint equation is achieved through residual analysis. The residual is defined as the difference between the left and right sides of the physical constraint equation, and the root mean square value of the residual is the trajectory quality assessment index. The reverse derivation process assumes that the physical constraint equation itself holds, and that the residual is entirely caused by position coordinate errors. Position coordinate errors are calculated through sensitivity analysis of the physical constraint equation, which measures the impact of small changes in position coordinates on the residual. The impact of the error in the X-component of the position coordinate on the residual is obtained by calculating the partial derivative of the residual with respect to the X-component. The partial derivative is calculated using a numerical differentiation method with a step size of 0.01 mm. The numerical differentiation calculation involves increasing the position coordinate X-component by 0.01 mm and recalculating the residual; the change in residual is divided by 0.01 mm to obtain the partial derivative value. The calculation method for the partial derivatives of the Y-component and Z-component is the same. The sensitivity matrix is ​​a matrix composed of the partial derivatives of all position coordinates with respect to all residuals, and its dimension is three times the number of residuals multiplied by the number of position coordinates. The error of the position coordinate sequence is obtained by solving a system of linear equations, which is: sensitivity matrix × position coordinate error vector = residual vector. The algorithm used to solve the equations is the least squares method, which calculates the position coordinate error vector that minimizes the sum of squared residuals.

[0133] The error quantities of geometric transformation parameters are derived based on the error quantities of the position coordinate sequence, utilizing the backpropagation relationship of the geometric transformation model. The geometric transformation model maps pixel coordinates to three-dimensional spatial position coordinates, and backpropagation calculates the dependence of position coordinate errors on geometric transformation parameters. The sensitivity of position coordinate errors to rotation matrix parameters is obtained by calculating the partial derivatives of position coordinates with respect to rotation angle. For small changes in rotation angle, the change in position coordinates is approximately equal to the cross product of the position coordinate vector and the rotation axis direction vector, multiplied by the change in rotation angle. The sensitivity of position coordinate errors to translation vectors is the identity matrix; the error in translation vectors is directly equal to the error in position coordinates. The sensitivity of position coordinate errors to lens distortion parameters is calculated using the partial derivatives of the distortion correction formula. Changes in the radial distortion coefficient lead to changes in the radial direction of the position coordinates, with the change proportional to the square or fourth power of the radial distance. A linear relationship matrix between geometric transformation parameter errors and position coordinate errors is constructed. The matrix is ​​solved using a pseudo-inverse method, calculated as matrix transpose × matrix inverse × matrix transpose, yielding estimates of the geometric transformation parameter errors.

[0134] The geometric transformation parameter correction subtracts the estimated error from the current parameter values. The corrected radial distortion coefficient is the current radial distortion coefficient minus the radial distortion coefficient error; the corrected rotation angle is the current rotation angle minus the rotation angle error; and the corrected translation vector is the current translation vector minus the translation vector error. The correction process is performed synchronously on all geometric transformation parameters to ensure consistent coupling between them. The corrected geometric transformation parameters are stored in a unified parameter matrix, replacing the original parameter values. The magnitude of parameter correction is constrained: the change in radial distortion coefficient in a single correction cannot exceed 10% of the current value, the change in rotation angle cannot exceed 0.05 radians, and the change in translation vector cannot exceed 5 millimeters. These constraints prevent over-correction that could lead to instability.

[0135] When the corrected geometric transformation parameters are synchronously applied to multi-view fusion transformation, the pixel coordinates of charged particles in the original 2D image sequence are reread, and distortion correction and view transformation are performed using the corrected geometric transformation parameters. The distortion correction process recalculates the corrected pixel coordinates based on the corrected radial and tangential distortion coefficients, while the view transformation process recalculates the 3D spatial rays based on the corrected rotation matrix and translation vector. The multi-view fusion calculation re-solves for the closest points of ray pairs, and weighted fusion yields the updated position coordinate sequence. Compared to the original position coordinate sequence, the position values ​​in the updated position coordinate sequence change, and the amount of change is related to the correction amount of the geometric transformation parameters. When the updated position coordinate sequence is applied to temporal correlation, timestamp sorting, trajectory breakage detection, interpolation repair, trajectory segment stitching, and smoothing filtering are re-executed to obtain updated 3D trajectory data. The position deviation threshold and acceleration threshold remain unchanged during the temporal correlation process to ensure the consistency of the temporal correlation logic.

[0136] The updated 3D trajectory data is used to recalculate the trajectory quality assessment index. The calculation process extracts the instantaneous velocity and acceleration from the updated position coordinates, queries the corresponding local field strength vector, calculates the theoretical acceleration, and calculates the root mean square value of the deviation between the instantaneous acceleration and the theoretical acceleration. The recalculated trajectory quality assessment index is compared with a preset quality threshold. If it still exceeds the threshold, the parameter correction process is repeated; if it is below the threshold, the optimization is considered converged, and the closed-loop image processing optimization mechanism completes one iteration. The trajectory quality assessment index value for each iteration is recorded. When the index change is less than 0.01 mm / millisecond squared for three consecutive iterations, convergence is determined, and the iteration process terminates. The maximum number of iterations is limited to 10 to avoid infinite loops. The closed-loop optimization mechanism is periodically triggered during operation, with a trigger cycle set to perform optimization once every 1000 frames of images processed, ensuring that the geometric transformation parameters remain accurate over time.

[0137] This method triggers parameter correction through trajectory quality assessment indicators, derives geometric transformation parameter errors in reverse based on physical constraint equations, and applies the corrected parameters synchronously to multi-view fusion and temporal correlation to form a closed-loop optimization mechanism, thereby continuously improving the accuracy and long-term stability of 3D trajectory data.

[0138] Secondly, a bidirectional vision-based system for tracking and acquiring the trajectory of charged particles is provided, including:

[0139] The first unit is used to simultaneously acquire a two-dimensional image sequence of charged particles in an electric field environment through multiple image acquisition devices with spatial positional relationships;

[0140] The second unit is used to perform multi-view fusion transformation on the pixel coordinates of charged particles in the two-dimensional image sequence based on the spatial positional relationship of the image acquisition device, by establishing a geometric transformation model that includes lens distortion compensation and dynamic baseline correction, to obtain a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel accuracy;

[0141] The third unit is used to temporally correlate the three-dimensional spatial position coordinates at different times based on the timestamp information and spatial continuity constraints in the position coordinate sequence, to obtain the three-dimensional trajectory data of the charged particle;

[0142] The fourth unit is used to calculate the trajectory quality evaluation index based on the motion state parameters in the three-dimensional trajectory data and the field strength distribution information of the electric field environment through the force balance relationship;

[0143] The fifth unit is used to deduce the error of the geometric transformation parameters in the geometric transformation model based on the trajectory quality evaluation index and the physical constraint relationship, and to correct the geometric transformation parameters. The corrected geometric transformation parameters are then applied synchronously to the multi-view fusion transformation and the temporal correlation to update the position coordinate sequence, forming a closed-loop image processing optimization mechanism.

[0144] A third aspect of the present invention provides an electronic device, comprising:

[0145] processor;

[0146] Memory used to store processor-executable instructions;

[0147] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0148] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0149] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0150] Thirdly, a computer-readable storage medium is provided, having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

Claims

1. A method for tracking and acquiring the trajectory of charged particles based on bidirectional vision, characterized in that, include: Two-dimensional image sequences of charged particles in an electric field environment are acquired synchronously by multiple image acquisition devices with spatial relationships. Based on the spatial positional relationship of the image acquisition device, a geometric transformation model including lens distortion compensation and dynamic baseline correction is established to perform multi-view fusion transformation on the pixel coordinates of charged particles in the two-dimensional image sequence, thereby obtaining a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel accuracy; Based on the timestamp information and spatial continuity constraints in the position coordinate sequence, the three-dimensional spatial position coordinates at different times are temporally correlated to obtain the three-dimensional trajectory data of the charged particle; Based on the motion state parameters and electric field strength distribution information in the three-dimensional trajectory data, trajectory quality evaluation indicators are calculated through force balance relationships; Based on the trajectory quality assessment index and the physical constraint relationship, the error of the geometric transformation parameters in the geometric transformation model is derived in reverse, and the geometric transformation parameters are corrected. The corrected geometric transformation parameters are then synchronously applied to the multi-view fusion transformation and the temporal correlation to update the position coordinate sequence, forming a closed-loop image processing optimization mechanism.

2. The method according to claim 1, characterized in that, The geometric transformation model that includes lens distortion compensation and dynamic baseline correction includes: By acquiring a sequence of calibration images containing known spatial location calibration points, the lens distortion parameters and relative position and attitude parameters of multiple image acquisition devices are simultaneously solved. The lens distortion parameters and the relative position and attitude parameters are then integrated into a unified parameter matrix to establish a geometric transformation model that couples lens distortion compensation and baseline correction. The system monitors the pixel position changes of calibration points in multiple image acquisition devices in real time. When the detected pixel position change exceeds a preset change threshold, it determines that the relative positions of the multiple image acquisition devices have changed. Based on the pixel position change, the reprojection error of the calibration point under different viewpoints is calculated. By minimizing the reprojection error, the lens distortion parameters and relative position and attitude parameters in the unified parameter matrix are jointly optimized and updated to obtain the updated geometric transformation model.

3. The method according to claim 1, characterized in that, Performing a multi-view fusion transformation on the pixel coordinates of charged particles in the two-dimensional image sequence to obtain a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel accuracy includes: The geometric transformation model is used to perform distortion correction and viewpoint transformation on the pixel coordinates of charged particles in the two-dimensional image sequence, thereby obtaining the corrected pixel coordinates corresponding to each image acquisition device; The electric field direction information of the electric field environment is obtained, the angle between the observation direction of each of the multiple image acquisition devices and the electric field direction is calculated, and a confidence weight is assigned to each image acquisition device according to the angle. The corrected pixel coordinates are weighted and fused according to the confidence weight, and the candidate position coordinate set of charged particles in three-dimensional space is solved by multi-view geometric constraints; Based on the physical constraints of the motion of charged particles under the influence of electric field forces in an electric field environment, the consistency metric between the direction of the electric field force and the direction of motion of the charged particles corresponding to each candidate position coordinate in the candidate position coordinate set is calculated. Candidate position coordinates with consistency metrics lower than the consistency threshold are eliminated, and candidate position coordinates that satisfy the force equilibrium condition are retained as the optimal position coordinates. The optimal position coordinates are subjected to temporal interpolation and sub-pixel refinement to obtain a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel precision.

4. The method according to claim 3, characterized in that, The steps of acquiring electric field strength direction information of the electric field environment, calculating the angle between the observation direction of each of the multiple image acquisition devices and the electric field strength direction, and assigning confidence weights to each image acquisition device according to the angles include: The local field strength vector at the location of the charged particles in the electric field environment is obtained, and the direction of the local field strength vector is extracted as the field strength direction information; the angle between the optical axis direction of the multiple image acquisition devices and the field strength direction is calculated respectively. Image quality parameters of the two-dimensional images acquired by each image acquisition device are obtained, including image sharpness and signal-to-noise ratio; an angle weighting coefficient is determined based on the numerical range of the included angle. When the included angle is greater than a preset vertical angle threshold, the angle weighting coefficient increases monotonically as the included angle increases; when the included angle is less than or equal to the preset vertical angle threshold, the angle weighting coefficient decreases monotonically as the included angle decreases. The angle weight coefficient and the image quality parameter are combined to obtain the initial confidence weight; the initial confidence weight is then normalized to obtain the final confidence weight.

5. The method according to claim 1, characterized in that, The steps for obtaining the three-dimensional trajectory data of the charged particle by temporally associating the three-dimensional spatial position coordinates at different times based on the timestamp information and spatial continuity constraints in the position coordinate sequence include: Extract timestamp information from the location coordinate sequence, sort them according to the timestamp order, and obtain a time-ordered location coordinate sequence; Calculate the spatial distance and velocity change between adjacent position coordinates in the time-series position coordinate arrangement. When the spatial distance exceeds the position deviation threshold calculated based on the motion velocity or the velocity change exceeds the acceleration threshold calculated based on the direction of the electric field force, it is determined that there is a trajectory break. For coordinates where there is a trajectory break, based on the velocity direction and electric force direction of the position coordinates before the break, the interpolated position coordinates between the break times are calculated, and the interpolated position coordinates are inserted into the time-series of position coordinates. The similarity in spatial location distribution and motion direction between different trajectory segments in the repaired temporal arrangement position coordinate sequence is calculated. When the similarity is higher than the association threshold, the trajectory segments belonging to the same charged particle are spliced ​​together in the order of timestamps. The spliced ​​position coordinate sequence is then smoothed and filtered to obtain the three-dimensional trajectory data of the charged particle.

6. The method according to claim 1, characterized in that, The steps for calculating trajectory quality evaluation indicators based on the motion state parameters and electric field strength distribution information in the three-dimensional trajectory data through force balance relationships include: The position coordinates and timestamp information of charged particles at different times are extracted from the three-dimensional trajectory data. The instantaneous velocity of the charged particles is obtained by calculating the difference between the position coordinates of adjacent times, and the instantaneous acceleration is also calculated. Obtain the field strength distribution information of the electric field environment, and query the local field strength vector corresponding to the position of the charged particle based on the position coordinate sequence; Based on the force balance relationship of charged particles under the action of electric force, gravity and air resistance in an electric field environment, the theoretical acceleration of charged particles is calculated according to the local field strength vector. The deviation between the instantaneous acceleration and the theoretical acceleration is calculated as a trajectory quality evaluation index.

7. The method according to claim 1, characterized in that, The steps of deriving the error of the geometric transformation parameters in the geometric transformation model based on the trajectory quality assessment index and the physical constraint relationship, correcting the geometric transformation parameters, and synchronously applying the corrected geometric transformation parameters to the multi-view fusion transformation and the temporal correlation to update the position coordinate sequence, forming a closed-loop image processing optimization mechanism, include: When the trajectory quality evaluation index exceeds the preset quality threshold, it is determined that there is an error in the geometric transformation parameters in the geometric transformation model; A physical constraint equation is established based on the force equilibrium relationship of charged particles in an electric field environment. The physical constraint equation describes the relationship between the motion state parameters and position coordinates of charged particles. Substitute the trajectory quality evaluation index into the physical constraint equation, calculate the error of the position coordinate sequence by reverse derivation, and derive the error of the geometric transformation parameters in the geometric transformation model based on the error of the position coordinate sequence. The geometric transformation parameters are corrected based on the geometric transformation parameters to obtain the corrected geometric transformation parameters; The corrected geometric transformation parameters are synchronously applied to the multi-view fusion transformation of the charged particle pixel coordinates in the two-dimensional image sequence, and the updated position coordinate sequence is recalculated. The updated position coordinate sequence is then applied to the temporal correlation to obtain the updated three-dimensional trajectory data. Based on the updated 3D trajectory data, the trajectory quality evaluation index is recalculated to form a closed-loop image processing optimization mechanism.

8. A charged particle trajectory tracking and data acquisition system based on bidirectional vision, used to implement the method of any one of claims 1-7, characterized in that, include: The first unit is used to simultaneously acquire a two-dimensional image sequence of charged particles in an electric field environment through multiple image acquisition devices with spatial positional relationships; The second unit is used to perform multi-view fusion transformation on the pixel coordinates of charged particles in the two-dimensional image sequence based on the spatial positional relationship of the image acquisition device, by establishing a geometric transformation model that includes lens distortion compensation and dynamic baseline correction, to obtain a sequence of position coordinates of charged particles in three-dimensional space with sub-pixel accuracy; The third unit is used to temporally correlate the three-dimensional spatial position coordinates at different times based on the timestamp information and spatial continuity constraints in the position coordinate sequence, to obtain the three-dimensional trajectory data of the charged particle; The fourth unit is used to calculate the trajectory quality evaluation index based on the motion state parameters in the three-dimensional trajectory data and the field strength distribution information of the electric field environment through the force balance relationship; The fifth unit is used to deduce the error of the geometric transformation parameters in the geometric transformation model based on the trajectory quality evaluation index and the physical constraint relationship, and to correct the geometric transformation parameters. The corrected geometric transformation parameters are then applied synchronously to the multi-view fusion transformation and the temporal correlation to update the position coordinate sequence, forming a closed-loop image processing optimization mechanism.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images