EFFICIENT MULTI-VIEW 3D TRACKING OF ARBITRARY ROCK FRAGMENTS UPON IMPACT

This paper presents a new methodology to accurately obtain 3D rotational velocities of blocks and fragments. Four high speed cameras are used to capture the scene. An additional two tilted mirrors are used to multiply the number of views. Hence, a total of six different viewing perspectives can be used to track translational and rotational velocities in 3D. The focus in the current work is on the rotational velocities, as tracking of the translation is generally straightforward. A common outline tracking algorithm based on the visual hull is adapted. The visual hull is further meshed using triangular elements to approximate the shape of the object. This 3D reconstruction is then used to track the 3D motion of the object. However, the accuracy of the results strongly depends on the accuracy of the 3D reconstruction which is mainly influenced by the number and position of the available views. In any case, the 3D reconstruction from the visual hull is only an approximation and significant errors can be introduced which influence the tracking accuracy. Hence, an inhouse post-processing algorithm based on the knowledge of the real geometry of the object, which can generally be accurately determined after a test, was developed. The improved performance of this new post-processing method is shown by controlled spinning tests. Finally, results of a real example of an impact fragmentation test are discussed.


INTRODUCTION
Rockfall is a natural hazard that requires a rigorous and comprehensive assessment to be successfully mitigated. Although rockfall trajectory modelling has advanced very rapidly in recent years, there are still many questions around how rocks behave upon impact. In particular, the fragmentation of blocks upon impact and the evolution of the resulting fragments is very poorly understood (Giacomini et al., 2009;Volkwein et al., 2011). It is necessary to investigate the fragmentation process in more detail to be able to accurately predict its likely outcomes such as number of fragments, volume of fragments, trajectories and residual energies. Images and videos recorded by high-speed cameras are commonly used for tracking rockfall trajectories during smallscale laboratory and large-scale field studies (e.g., Asteriou & Tsiambaos, 2018;Azzoni & De Freitas, 1995;Dorren et al., 2005;Giacomini et al., 2012;Giani et al., 2004;Glover, 2015;Spadari et al., 2012;Volkwein & Klette, 2014). Different algorithms can be used to track moving objects, such as feature tracking, circle tracking, centre of gravity tracking, corner contour tracking or outline tracking. Feature tracking or colour feature tracking is widely used because of its simplicity (Zhivkov & Kruse, 2004). Such tracking can be conducted from one or more viewpoints to obtain two-dimensional (2D) or threedimensional (3D) coordinates of the feature points. Translational velocities and accelerations of the colour feature or object can then be estimated. To obtain the rotational velocity, two or more features including the centre of gravity (CofG), have to be tracked (Shum & Komura, 2005). This allows estimation of the velocity around the axis of rotation. However, feature tracking is not always practical, especially for tracking 3D rotational velocities. An alternative is the more sophisticated and * Corresponding author computationally expensive outline tracking algorithm, which is based on the idea of reconstructing an approximated 3D shape of an object from silhouettes captured from different viewpoints. The 3D reconstructed geometry is referred to as a visual hull construction or shape-from-silhouette (Cheung et al., 2003). This algorithm has some advantages compared to feature tracking as it also provides an approximated 3D shape of the object in the form of a 3D mesh. The geometric properties of this mesh (e.g. volume, centre of mass, moment of inertia and principal axes) can easily be determined and used in the analysis. However, several well-positioned camera views (more than three) are generally required for a good approximation of the visual hull since the accuracy depends on the number of distinct silhouette images. This requirement may be a limiting factor in many situations because of financial (buying more cameras) or physical (system setup) reasons. Nevertheless, as pointed out by several researchers (Gluckman & Nayar, 2001;Murray, 1995;Teoh & Zhang, 1984;Fujigaki et al.;, it is possible to reconstruct a scene by imaging the scene reflection in a planar mirror. As such, the number of views can be increased without the need for additional cameras. Several algorithms were developed to track the 3D motion of an object from an array of cameras, among which is the visual hull approach where the shape of the object is approximated from silhouettes. Such approach is flexible but suffers from some limitations potentially leading to significant errors and inaccurate estimation of rotational velocities. Consequently, a postprocessing algorithm based on a more accurate fragment geometry, obtained from scanning post-testing, was developed and combined to the visual hull approach. This paper presents this new post-processing algorithm and its validation, which was achieved through controlled spinning tests and a drop test.

EXPERIMENTAL SETUP
A program of experimental testing in a controlled environment, with accurate tracking of 3D motion of high speed motion using high-speed cameras, is being undertaken to better understand the rockfall fragmentation phenomenon. A custom-made hexagonal fragmentation cell (Figure 1) was developed at the University of Newcastle (Australia) in which blocks can be dropped in a safe and controlled way onto a concrete slab (Guccione et al., 2019). The walls of the fragmentation cell alternate between plywood and transparent polycarbonate. This arrangement allows three cameras to be set up in front of the transparent polycarbonate panels with the opposite plywood panels acting as background. An additional camera is set up on an extended tripod looking down onto the slab where the block is impacting. Hence, a total of four high-speed cameras are installed. The cameras are all of the same type: Optronis CR600x2. The main characteristics and settings of the cameras are given in Table 1. A custom-built external synchronisation box is used to simultaneously trigger the four cameras. The camera on the extended tripod (Cam 4) uses a Nikkor AF 50mm f/1.8D lens whereas the three cameras behind the plywood (Cam 1, 2 and 3) use a Nikkor AF 35mm f/2D. A frame rate of 500 frames per second (fps) with a shutter speed of 1/3003 s was used. Adequate lighting is needed to reduce object blurring while retaining a reasonable depth of field and image quality. To achieve this, several LED strips and LED spotlights were installed on each side of the fragmentation cell. In addition, LED panels were mounted on each plywood panel opposite Cam 1, 2 and 3 allowing for a very good contrast between fragments and background. Initial tests by the authors suggested that four cameras are not enough to accurately track all major fragments upon impact as some fragments were not visible in all four views (Guccione et. al., 2019). Hence, it was decided to add mirrors to the setup to increase the number of views. Tilted mirrors were installed on the plywood panels opposite Cam 1 and Cam 2 ( Figure 1 and Figure 2), offering two additional views. A total of six views are thus available for the final image processing. The views (V1-V6) are named according to Figure 2b. It should be noted that a compromise had to be made when focusing Cam 1 and Cam 2 as the object distance in the mirrored views is much bigger than that in the real views ( Figure 2a). Also, in order to obtain V5 and V6 the images taken with Cam 1 and Cam 2 had to be flipped horizontally before processing. An example image and its flipped counterpart are shown in Figure 3a   The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition)

Outline tracking
The commercial software TEMA3D was used to process the images from the different views. For V5 and V6, the images of Cam 1 and Cam 2 were flipped horizontally and masked so they only show the relevant view (see e.g. Figure 3b). In each view, the object was tracked by a suitable outline tracker. An outline tracker finds the outline of the object in the image (Figure 4) by using a thresholding algorithm. In order for the thresholding algorithm to work properly, it is necessary to have a high contrast between the object and the background (achieved by the LED panels in the background, see Section 2). The outline is made up by four-connective chain code, which is a sequence of movements, up, down, left or right (Anliot, 2005). The chain code is complete when the sum of all movements ends up back at the starting pixel.

Visual hull
The outline of a view shows the silhouette of an object. The silhouette back-projected to the camera creates a 3D viewing cone within which the object resides. The intersection of all silhouette viewing cones, one from each camera view, gives the 3D visual hull of the object (Anliot, 2005), which is then approximated as a 3D triangular mesh. The physical attributes of the approximated 3D mesh (such as volume, principal axes, moment of inertia) can then be calculated. It is clear that the accuracy of the visual hull depends on the number of available views. Figures 5a and 5b, for example, show the concept of the visual hull of a fragment in 2D for a case with two and three cameras respectively. It can clearly be seen that the visual hull is only an approximation of the real shape and, the more views, the more accurate the visual hull. Another major limitation of the visual hull is that it cannot capture concave pits on an object. This means that the estimated volume is always overestimated.
The approximation of the shape using the visual hull algorithm is generally good enough to estimate the barycentre of the object, which is needed for estimating the translational velocities. However, estimation of the principal axes is generally not very accurate. This is a major issue, as the principal axes are needed to estimate the angular velocity. Another limitation of this approach resides in the difficulty to track rotation of symmetric objects, such as spheres or cubes. Although this is generally not a problem for fragments that tend to be irregular, it is problematic to track rotation of regular objects in free fall, prior to impact and fragmentation. In order to address some of the current limitations and improve accuracy of angular velocity tracking, a new postprocessing algorithm was developed. This is presented in Section 3.3.

Post-processing algorithm
The proposed post-processing algorithm relies on the real geometry of the object, which can generally be determined before and/or after the test by an accurate photogrammetric survey or a 3D scan. In this work, a structured light scanner (EinScan Pro 2X Plus) was used to reconstruct the real 3D geometry of the object. The objective is to align the mesh created from the scanned real geometry with the visual hull mesh exported from TEMA3D, at each time step in order to calculate the orientation of the principal axes ( Figure 6). The alignment process is performed using the iterative closest point (ICP) algorithm where points of the real geometry are automatically aligned to a sub-sampled point cloud of the visual hull. This step is executed within the open-source program CloudCompare whereas the calculation of the principal axes is performed with the open-source library trimesh (Dawson-Haggerty et al., 2019). This two-step process is repeated for each time step. Finally, the angular velocity is estimated based on the orientation of the new calculated principal axes. The procedure can be summarised as follows, where all steps, except step 1 and 2, have been implemented into a Python script: 1. Automatically export visual hull meshes for all time steps from TEMA3D (note that AutoHotkey was used to automatize the process of exporting the mesh, at each timestep). 2. Manually align the mesh of the real geometry with the first visual hull mesh ( = 1). This is necessary to ensure optimal convergence of the ICP. 3. Sub-sample all visual hull meshes to achieve a point density similar to that of the real geometry. 4. Automatically align the real geometry (manually aligned mesh in the first timestep =1 and mesh from the previous timestep − 1 thereafter) to the subsampled point clouds using the ICP.  The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) 5. Calculate the orientation of the principal axes using the aligned real geometry. This is done using the library trimesh. It should be noted that the calculated axes do not always represent a right-handed coordinate system. Hence, the third principal axis is always recalculated as cross product of the first two principal axes. 6. Check the orientation of the principal axes for timestep based on the orientation in the previous timestep − 1. This is necessary since the solution of the principal axes is not unique (four solutions are possible), i.e. the axes can be orientated in the positive or negative direction. Hence, changing the direction of the principal axes by 180 degrees is sometimes required. 7. Calculate the rotation between timesteps and − 1.
For each timestep, the rotational increment around each axis is calculated relative the axes' positions from the previous timestep. This increment is then accumulated with the increments from earlier times, to produce the total rotation around the principal axes since timestep = 1. 8. Fit linear trendline to accumulated increments to estimate the angular velocity around the principal axes.
The main steps of the post-processing algorithm are outlined in Figure 6. An example of a series of principal axes directly obtained from TEMA3D and after post-processing is shown in Figure 7a and 7b respectively.

EXPERIMENTAL PROCEDURE
A static calibration was performed before conducting the tests, by taking a picture of a calibration structure with 18 reference points ( Figure 8) simultaneously with all the cameras. The calibration was performed through TEMA3D by collimating at least 5 reference points. The object distance, the ground sampling distance (GSD) and the average root mean square error (RMSE) of all reference points after calibration are reported in Table 2. It can be seen that the RMSE of the reference points of the normal camera views (V1, V2, V3 and V4) is generally smaller than two times the GSD. For the mirrored views (V5 and V6), the RMSE is slightly bigger than two times the GSD. Two controlled spinning tests were conducted in order to evaluate the benefit of two additional mirrored views and to validate the proposed post-processing method. In particular, two different objects were rotated around an axis with a known rotational velocity. The first object used is a prismatic masonry brick (Figure 9a). Its dimensions are 7.5cm x 11cm x 22.9cm and its volume is equal to 1,872 cm 3 . The second object is a typical fragment resulting from a drop test of an artificial rock sphere with diameter of 100 mm (Figure 9b). Its shape can be approximated to a quarter of a sphere with two rough surfaces due to the crack propagation. Its volume is 131 cm 3 . Finally, a drop test with an artificial rock sphere was performed to assess the ability of the setup and the new post-processing algorithm to accurately capture the 3D motion of irregularly shaped fragments. The sphere (Figure 9c) has a diameter of 100 mm and its volume is 524 cm 3 . It was released from a drop height of about 3.05 m.

Validation
The brick (Figure 9a) was rotated around an axis of rotation very close to the third principal axis at a velocity of 11.99 rad/s. Results showing the influence of the number of views and the comparison between the original method implemented in TEMA3D and the new post-processing (PP) method are plotted in Figure 10. The figure shows the cumulative rotation angles around the three principal axes calculated with both methods over time. Fitting a linear trend through the data points provides an estimate of the angular velocity around the principal axes. In the case of a regular and simple shape, such as the brick, there is not much difference between four and six views. Indeed, Guccione at al. (2019) showed that in such cases even two views, if well positioned, are enough to capture the main component of the rotational velocity. When comparing the new post-processing method with the original results from TEMA3D, it can be seen that TEMA3D shows fluctuations of the angle of the principal axes were the angular velocity is low. The proposed postprocessing method seems to provide smoother and more accurate results. Table 3 summarises the computed rotational velocities around the main axes ( 1 , 2 , 3 ), the absolute rotational velocity ( ) and the relative error with the reference value (1 − / ). The new method shows results closer to the reference for both, four views and six views. The fragment (Figure 9b) was also rotated around an arbitrary axis of rotation close to the third axis at a velocity of 37.46 rad/s. The fragment has an arbitrary shape and it can be seen that in such a case, the more views, the better the result (Figure 11). The benefit of using the new post-processing algorithm is clearly visible when using an irregular shaped fragment: the error drops from around 30% to around 8% when using 4 views; and from 5% to 1% when using 6 views (    The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition)

Application
Figures 12a and 12b show images of the free-falling sphere before impact and after impact, respectively, taken from view V3. It can be seen that upon impact the sphere broke into three major fragments. A few very small fragments are also seen but fragments with mass smaller than 5 g were not tracked. Figures  12c and 12d show the corresponding tracking results from TEMA3D. Six views were used to track the sphere before the impact and the visual hull appears to be very close to a sphere. The volume of the visual hull is 539 cm 3 whereas the correct value is 524 cm 3 , which is very close. The fragments after impact are also relatively well approximated by the visual hull.   Table 5 is relative to the instant that tracking of each fragment started. This is generally when the outline of the fragment to be tracked is well defined and no overlaps with other fragments are seen from at least four views. Views where fragment overlap exists were discarded. This means that the number of views used in the image processing varies (for F1 and F2 between 4 and 5, for F3 between 5 and 6). Figures 12c and 12d show that the trajectory of the sphere in free fall follows a straight line while the trajectory of fragments follows a parabola. These two observations are consistent with laws of physics. The measured accumulated angle around the principal axes and the estimated rotational velocities of the fragments before and after post-processing are illustrated in Figure 13. The vertical dotted lines indicate the timestep where the number of views is changing according to Table 5. It can clearly be seen how the number of views influences the accuracy of the results. The original results of TEMA3D show quite a large scatter ( Figure  13b) and sometimes even a change in gradient (Figure 13a). From physics, and ignoring air resistance, we know that the rotational velocity should be constant. A change in gradient (i.e. a kink in the rotational velocity) is hence physically not possible unless the change comes from external factors such as the collision with another fragment. The fragments did not interact after separation and, hence, there is no reason for the rotational velocity not to be constant. It can be seen that the new post-processing algorithm resolves these issues and noticeably improves the estimation of the rotational velocity. Only for fragment F3 (Figure 13c), the change from six to five views introduces a very small jump and a slight change in the gradient for the second axis.  Table 5).

CONCLUSIONS
In this paper, a new post-processing algorithm for the accurate estimation of 3D rotational velocities is presented. Four synchronised high-speed cameras plus two additional views, obtained by installing tilted mirrors, were used to record drop tests and controlled spinning tests. Tracking of fragments was performed using outline tracking visual hull algorithms in TEMA3D. However, the initial results were found to be inaccurate for generally shaped objects and a novel postprocessing algorithm was developed to improve accuracy and solve issues associated to the identification of the rotational velocity. This new algorithm relies on a more accurate representation of the objects, obtained by scanning before or after the test. The proposed post-processing algorithm was implemented into a Python script and used to re-analyse the data obtained with TEMA3D. Validation tests and drop tests showed that the proposed algorithm significantly improves the estimation of rotational velocities. It was also shown that, with this new algorithm, less views can be used to accurately predict rotational velocities.