Continuous-Time Fast Motion of Explosion Fragments Estimated by Bundle Adjustment and Spline Representation Using HFR Cameras

: This paper introduces a continuous-time fast motion estimation framework using high frame-rate cameras. To recover the high-speed motions trajectory, we inherent the bundle adjustment using a different frame-rate strategy. Based on the optimized trajectory, a cubic B-spline representation was proposed to parameter the continuous-time position, velocity and acceleration during this fast motion. We designed a high-speed visual system consisting of the high frame-rate cameras and infrared cameras, which can capture the fast scattered motion of explosion fragments and evaluate our method. The experiments show that bundle adjustment can greatly improve the accuracy and stability of the trajectory estimation, and the B-spline representation of the high frame-rate can estimate the velocity, acceleration, momentum and force of each fragments at any given time during its motion. The related estimated result can achieve under 1% error. show that our adjusted bundle adjustment cooperated with B-spline representation, high accuracy and robust position, that the velocity and acceleration can be estimated at any time, and that the related error of the estimated motion can be controlled under 1% by using high frame-rate cameras.


Introduction
A continuous-time motion analysis has great engineering support for the damage assessment and cause analysis of an explosion [1]. Reconstructing a continuous-time trajectory of scattered fragments has great value. In this paper, we present a method of fast motion estimation of the explosion fragments. The difficulty of the fragment's motion estimation is randomly scattered. The fragments created by an explosion follow these four stages: 1. Radial expansion of the shell; 2. A crack appears somewhere in the shell; 3. Material leakage, resulting in pressure difference between the inside and outside of the shell; 4. The shell is broken randomly, forming irregular fragments which are scattered with a related high initial velocity.
To capture this fast, diffuse motion, the high frame-rate (HFR) cameras are deployed in this research. We used two HFR cameras to build a high-speed vision system capturing the high-speed motion.
The HFR camera has overcome the restrictions posed by conventional video signals. Many offline HFR cameras can be operated at over 1000 FPS (Frame Per Second) during a very short time [2][3][4]. HFR cameras can help us recognize the high-speed phenomena in the real world, which cannot be captured by conventional cameras (e.g., 30 FPS). Many

Bundle Adjustment Reconstructs Motion's Trajectory
The key idea of BA is that there exist bundles of light rays from the camera center and a projected position in images at each moment, and all parameters adjust simultaneously and iteratively together to obtain the optimal bundles of all times. Reprojection errors in bundle adjustment come from all deviations between the observed position in the projected images and the reprojected position at the same image at each time. The motion and the trajectory are estimated along with the image sequences whereas the cameras are static to the world coordinate. We define the reproject error as the following equation: where is the observed position, is the reprojected position, i is the identification of the recording camera, is the time stamp.
, , , , The reprojected function can be expressed in short as: where is the position of the world coordinate at the time stamp ; is the rotation vector presented in lie algebra so (3) from the world coordinate to camera coordinate; is the translation vector from the world to camera , presented in camera i coordinate; is the intrinsic parameter of camera = , , whereas is the focal length of the camera i; and are the distortion parameters of the camera. Equation (2) follows the pinhole camera model with the distortion expression. It can be detailed by the following equations: The optimized parameters can be expressed as, with m recording cameras and n captured positions, , , , , , , , All single reprojection errors form a × 1 reprojection vector 11 21 1 , ,..., ,... ,..., The optimization process for the bundle adjustment is finding the minimum of all reprojection errors. In this paper, we use the sum of the squares constructed by a large number of nonlinear functions as the cost function to form a nonlinear least-square problem.
To solve this unconstrained optimization Equation (4), basing on Taylor's theorem, a combination of linear search and trust region method Levenberg-Marquardt (LM) algorithm is deployed in this paper. For each iteration of the LM algorithm, the increment of optimized i can be calculated by the following equation: where is the Jacobian matrix, calculated analytically by the partial derivatives of the reprojection error; damping factor controls the LM algorithm search for the result following the linear search or trust region; D is the damping matrix. For the original LM algorithm, D is set as the identity matrix. In this paper, we set D as a diagonal matrix; each element at the diagonal is the eigenvalue of the matrix . The benefit of this design is that it can balance the magnitude gap between the position, orientation and intrinsic parameters. + λ is also denoted as the Hessian matrix. To solve the bundle adjustment, the key is the Jacobian matrix. From Equation (3), we give the analytically partial derivatives of the reprojection error to each component of by using rotation lie algebra so (3) properties and the chain rule of derivatives. For iw φ , the partial derivative can be expressed as: where the equation in that means the skew matrix of the related vector; the other two partial derivatives can be easily calculated by using Equation (3). is the left Jacobian for the BCH formula, with For translation vector , the partial derivate is expressed in Equation (7). Meanwhile, Equation (8) is for intrinsic , and Equation (9) is for position at time stamp j, . Equations (6)-(9) form the Jacobian matrix. We denote , , formed by Equations (6)-(8) as the camera's Jacobian for camera i at time stamp j, with formed by Equation (9) as the position's Jacobian observed by camera i at time stamp j. The Jacobian matrix will be detailed by the subsequent section which combines the HFR camera's properties and deploying strategy.

Jacobian Matrix for Stereo HFR Cameras
Stereo vision is widely used for non-contact measurement, which is suitable for the physical evaluation of an explosion. In this section, we use the synchronized stereo HFR camera's captured motion. For each fragment captured by these stereo HFR cameras, during its movement, the formation of the Jacobian matrix is shown in Figure 2. Figure 2 presents the HFR cameras with 540 FPS of the first 20 samples (e.g., for the first 37 ms). In Figure 2, each 4 × 18 matrix block of yellow and purple presents the Jacobian matrix of both cameras at time stamp j, and the matrix block in yellow are the . The matrix block in light blue is at time stamp j. The matrix block in purple from Figure 2 are zeros. Both the Jacobian matrix and Hessian Matrix are sparse matrices. Therefore, we use the Shur Complement to solve Equation (5). Equation (5) can be specified as: Figure 2. (a) Jacobian matrix for stereo high frame-rate (HFR) cameras; (b) related Hessian Matrix of these two HFR cameras.
Since the dimension of the cameras is related less to the position, is the diagonal block matrix which makes the + λ full rank. We use the Shur complement in the following: Equation (11) efficiently reduces the complexity of the calculation. From Figure 2b we can see that since matrix P is a block diagonal matrix, its inversion is relatively simple. The original Equation (5) can be easily solved by Equations (11) and (12), especially since the number of captured images is large.

Jacobian Matrix for HFR Cameras with Frame-Rate Multiplication
For HFR cameras, the captured frame-rate is highly reliant on the resolution of the frame, the exposure time, the bandwidth of the hardware, etc. In general, the resolution affects the frame-rate the most. the resolution influences the monitoring field of view (FOV). For fragments produced by the explosion, we want to use high frame-rate (over 500 FPS) to capture the moment of the explosion, and we also want to record the scattered fragments as long as possible.
To solve this problem, we use two HFR cameras with frame-rate multiplication, e.g., one HFR camera under 540 FPS, and another HFR camera under 1080 FPS. The frame-rate of the second HFR camera is twice that of the first camera. This is the same as Figure 2 for the first 37 ms, i.e., 20 frames for the first camera and 40 frames for the second camera. The Jacobian matrix of these two HFR cameras are presented in Figure 3. Each 6 × 18 matrix block of yellow and purple presents the Jacobian matrix of both cameras at time stamp j and j + 1, with 4 × 18 for both cameras at time j, 2 × 18 for camera two at time j+1. The matrix block in yellow is the . The matrix block in light blue is at time stamp j and for time j+1. The Jacobian matrix and Hessian matrix are also the sparse matrices, and with Equations (11) and (12), the iteration step can be solved efficiently. The interval between captured frames and the synchronization of our HFR cameras is strictly controlled by hardware instruments. This gives us the benefits that: (1) the Jacobian matrix is formatted strictly along the time axis with the same interval; (2) the Jacobian matrix presented in this paper is blocked diagonally, which reduces the complexity of the iteration. These are clearly demonstrated in the examples of Figures 2 and 3.
To improve the accuracy and reduce the iteration step, the initialization of i is important. For synchronization stereo HFR cameras, we can initial i by a stereo calculation at each time stamp. In this section, for a certain time stamp, there is only one camera observing the fragment. To solve this problem, we use pose averaging to initialize the position which is observed by only one camera, i.e., ( ) where k is the time interval to the last time stamp where at least two cameras observed, and ' k is the time interval to the follow up time stamp where two cameras captured the fragment. Based on our initialization strategy, the iteration step can be reduced to less than 300 steps for both frame synchronization stereo cameras and the frame-rate multiplication where we set , and the norm of all reprojection errors is less than 10 pixels.

Continuous-Time Spline Representation for Fast Motion Estimation
With the bundle adjustment under HFR cameras, each fragment's movement can be optimized to a trajectory sampled by the high frame-rate. For continuous-time motion estimation, the velocity and the acceleration cannot be calculated preciously by the difference in the bundle adjustment's results directly. To solve this problem, we use B-spline for continuous-time motion estimation, i.e., position, velocity and acceleration estimation based on bundle adjustment optimized results.

B-Spline Interpolation
Splines are an excellent choice for representing a continuous trajectory because their derivatives can be computed analytically. B-splines are chosen in this paper as they are a well-known representation for trajectories in R 3 [22][23][24]. B-splines use a series of control points (also denoted as knots) to define a linear combination of continuous function which is constrained by the knots. The interpolation between a series of 3D positions (knots) presented in the B-spline is as follows: where , ( ) is the De Boor-Cox recursive formula, the fractional coefficient for each given control point . The formula of , ( ) is:

Cubic B-Spline Representation for Continuous-Time Motion Estimation
Estimated velocity and acceleration need a C 2 continuity, a good approximation of minimal torque trajectories and a parameterization of rigid-body motion devoid of singularities [25,26]. A cubic B-spline representation can satisfy these conditions. We update the 3D position to SE(3), i.e., where is the rotation between the world coordinate and the fragment coordinate at time stamp j. When we treat the fragment as a single point, We define time fraction as: . By using lie algebra of SE(3) and the related propriety, we give the analytical result for motion estimation of the position (Equation (18)), the velocity (Equation (19)), and the acceleration (Equation (20)) at time u(t) during time and .

u t T A A A A A A A A A A A A A A A A A A
We define , the parameters are under these formulas:

Cubic B-Spline Representation for Continuous-Time Motion Estimation
Using Equations (18)- (20), the continuous-time's position, velocity and acceleration of a moving object can be estimated. The sampling frequency of knots (i.e., the frame-rate of cameras) influences the accuracy of the estimated result. We use HFR cameras with a different frame-rate to capture the known motion of a moving object. Using cubic B-spline motion estimation, the distribution of the position error, velocity error and acceleration error in three axes of the world coordinate are shown in Figure 4.
With the initial velocity at around 30 m/s, with a captured frame-rate of over 500 FPS, the velocity error of the B-splines estimation can be controlled under 16 cm/s. The percentage error of the position, velocity and acceleration are under 1%.

Experiments Setup of Explosion Simulation
We used two balloons with coffee beans installed between them to simulate the explosion fragments. The reason why we installed coffee beans between two balloons is because the broken shell's initial velocity is slow. With around 0.1 g of coffee beans, a high speed scattered motion can be created. Similar to the regular explosion presented in the introduction, our explosion followed these four similar stages: 1. The balloons are inflated in the inner shell and the inflation behavior is mainly radial; 2. The elastic limit causes breakage of the balloons; 3. Because of the elastic contraction of the balloon, the pressure difference between the inside and outside is constructed at the same time; 4. The shells of the balloons are randomly broken, and the coffee beans are distributed randomly during the inflation. After the balloons are broken, the coffee beans between these two balloons are scattered with a high initial velocity due to the elasticity and pressure difference. We used two Emergent Vision Technology (EVT) HFR cameras (10 GigE, HT-2000). When using these two cameras for experiments, the frame-rate can be set from 1080 FPS to 91 FPS, where the resolution of these cameras are from 512 × 320 to 2048 × 1088, respectively. HFR cameras captured the scattered motion of the explosion. We also deployed a monitoring module composed of two long-wavelength infrared (LWIR) cameras and one large FOV visible camera under the regular frame-rate (60 FPS). The covering FOV of this module is 120 degrees in both the visible band and LWIR band. The monitoring module helped us (1) monitorthe explosion's environment for safety concerns, (2) control the stop time of the HFR cameras to save storage space, and (3) capture the possible landing position of each fragment for further research. Our instruments and constructed balloons are shown in Figure 5.

Evaluation
The procedure of the continuous-time fast motion estimation is detailed as follows, shown in Figure 6, where: (a) Explosion fragments' motion are captured by the HFR cameras with a related framerate setting strategy (synchronization or multiplication); the large FOV LWIR cameras monitor the explosion (shown in Figure 6a) (b) With recording images from both HFR cameras, we use the dense optical flow: the Farneback method [27] is used to extract the motion area of the images at each time stamp j. The pixel motion of different fragments are analyzed (shown in Figure 6b). (c) The bundle adjustment is optimized based on the HFR cameras' working strategy, which has been detailed in the previous section. Figure   The accuracy of our BA optimization with different frame-rate strategy is presented in this section. We compare synchronization frame-rate bundle adjustment (denoted as SBA) and bundle adjustment with frame-rate magnitude (denoted as BAM) to stereo calculation directly. Euclidean distance is deployed to evaluate the position error, while the norm of the reprojection error from each camera evaluates the accuracy of the proposed methods. We use boxplot to demonstrate the distribution of the position error, shown in Figure 7a, and the reprojection error as shown in Figure 7b (Where -Ci means reprojection in camera i).  Table 1.  Figure 7 and Table 1 demonstrate that by using bundle adjustment optimization, the accuracy and stability of the trajectory estimation is highly improved, especially using the frame synchronization strategy (SBA). BAM can efficiently increase the motion's recording time. Comparing to SBA, the lack of observed motion from camera one at certain time stamps for BAM causes roughly 5 mm of accuracy degradation of the position, 0.001 pixel in MESRE. This is why the outlier of BAM-C2 is significantly high. When captured motion is around the camera's center area (around 45 pixels), the distortion affection is small, and the position calculated by stereo directly can achieve related high accuracy. Bundle adjustment optimization updates the intrinsic parameters instead of the position around these areas. This is why the minimum position error of stereo calculation directly is lower than BA's optimization. This is the only drawback of trajectory estimation using bundle adjustment.
A representative of one fragment's motion estimated by SBA and BAM with uncertainty bounds (3σ) is shown in Figure 8. We plot the whole 419 ms from the explosion to being scattered out of the cameras' FOV. The sample frame-rate for SBA of both HFR cameras is 540 fps, with 1024 × 480 resolution. The frame-rate for BAM is: one camera under 540 fps, with 1024 × 480 resolution, another under 1080 fps, with 512 × 320 resolution. The sampling data calculated by B-spline representation is six times the higher framerate, e.g., 6480 FPS, which exceeds the limits of our HFR cameras. The uncertainty data of five fragments scattered in this explosion are shown in Tables 2 and 3. All these supplementary data are calculated by B-spline representation. Based on estimated velocity and acceleration, the momentum and force of each fragment can also be provided. These physical parameters are very important for damage assessment, etc.

Discussion and Conclusions
In this paper, a continuous-time explosion fragments motion estimation method is proposed. We balanced the FOV and the frame-rate of the HFR camera, designed a highspeed vision system, and captured the moment of the explosion and the motion of the high-speed fragments. Basing on our frame-rate strategy, an adjusted bundle adjustment, i.e., SBA and BMA are proposed for trajectory estimation of high speed. With the optimized trajectory, a cubic B-spline representation of the continuous-time position, velocity and acceleration is proposed for this fast motion. Furthermore, the momentum and force can also be estimated at any given time. With these quantization parameters, the damage assessment of fragments can be evaluated.
The experiments show that our adjusted bundle adjustment cooperated with B-spline representation, high accuracy and robust position, that the velocity and acceleration can be estimated at any given time, and that the related error of the estimated motion can be controlled under 1% by using high frame-rate cameras.

Conflicts of Interest:
The authors declare no conflict of interest.