3D particle tracking velocimetry for the determination of temporally resolved particle trajectories within laser powder bed fusion of metals

Within this work, we present a system for the measurement of the three-dimensional (3D) trajectories of spatters and entrained particles during laser powder bed fusion (L-PBF) of metals. It is comprised of two ultrahigh-speed cameras and a reconstruction task specific processing reconstruction algorithm. The system enables an automated determination of 3D measures from the trajectories of a large number of tracked particles. Ambiguity evolving from an underdetermined geometrical situation induced by a two-camera setup is resolved within the tracking using a priori knowledge of L-PBF of metals. All processing steps were optimized to run on a graphics processing unit to allow the processing of large amounts of data within an appropriate time frame. The overall approach was validated by a comparison of the measurement results to synthetic images with a known 3D ground truth.


Introduction
Laser powder bed fusion (L-PBF) of metals distinguishes itself from conventional manufacturing processes by various aspects. Most importantly, it is characterized by a high degree of freedom of design and allows flexible manufacturing without cost-intensive tooling [1]. This makes the technology an important enabler for the production of highly optimized parts for applications pushing the envelope (e.g. harsh environments [2] and weight reduction [3]).
However, the production of defect-free and highly dense parts out of arbitrary metals today requires extensive parameter studies and process knowledge. Nevertheless, it has been shown that a more general process understanding can reduce the complexity of parameter determination [4].
However, within this approach, process dynamics are disregarded, even though they do have a large impact on the process result. These are mainly impacted by evaporation processes within the interaction zone (IZ) [5,6].
The evaporation effects within the IZ induce a shear flow in the surrounding atmosphere, applying a drag force to particles around the IZ. This leads to a local redistribution of the powder particles forming a denudation zone [5]. The evaporation direction and the speed of the vapor jet itself depend on the material, the process parameters, the inert gas, and the ambient pressure [7].
Particles entrained into the laser beam can be melted be molten, leading to the formation of spatters [6], or the inclusion in the melt pool, contributing to material buildup [5].
Both small and large spatters can cause process defects leading to a reduction in the quality of a part. However, large spatters cannot be fully melted by the laser beam in subsequent layers, leading to the formation of fusion defects within the processed layer [8]. These may result in voids which have a substantial negative impact on the mechanical properties of the part [9]. In contrast, depending on the material, small spatters can lead to defects caused by oxidation [10].
In addition, the stability of the spatter behavior correlates with the smoothness of the track morphology [11]. This is in accordance with the observation that the spatter characteristics highly coincide with the behavior of the plume and qualify, therefore, as an indicator of evaporation in the IZ [12].
The current process knowledge suggests a correlation between the process result and the spatter behavior. To investigate this effect, insight into and quantification of the spatter properties is crucial. However, velocity measurements presented in the published literature only investigate the velocity within the projection plane of the observing camera, resulting in a projection error [12,13]. In addition, the approaches are not applicable to large datasets.
These limitations can be overcome by the use of automated three-dimensional (3D) particle tracking velocimetry (3D PTV). 3D PTV is a Lagrangian measurement method that uses observations of single particles within a fluid to gain 3D velocity information with high spatial resolution [14]. It is based on the detection, matching, and tracking of individual particles within an imaged scene. Velocity information on the fluid flow field is made accessible by determining the 3D trajectories of tracer particles [15]. A 3D PTV system is comprised of four main parts: the illumination system, image recording devices, tracer particles, and image evaluation methods. The image evaluation consists of particle detection, correspondence analysis, 3D reconstruction, and a tracking algorithm. The tracking can be conducted in two dimensions or three [14]. 3D PTV is used in a variety of applications to non-intrusively identify flow fields of fluids.
The objective of this work is to introduce a method for gaining insight into the particle behavior within L-PBF based on 3D PTV. We present an experimental setup that allows spatial and temporal highly resolved process observation via a calibrated setup of two ultrahigh-speed cameras under realistic process conditions. We describe a specialized 3D PTV method, which makes 3D particle trajectories within L-PBF accessible. We show that ambiguity within the data of single time steps can be resolved by the incorporation of a priori process knowledge.
Our tracking method delivers the measured positions for each individual particle. This information makes characteristic 3D measures (e.g. the velocity, direction, and time of ejection) accessible. Furthermore, statistical measures over whole experiments can reveal trends within the data, which can provide a more general process understanding.

Experimental setup
For the experiments presented in this paper, 316L stainless steel (1.4404) was used as the powder and base material. The powder particle size was between 20 and 53 μm. It was spread with a scraper on a sandblasted specimen. In order to achieve reproducible layer heights, distance gauges were used. The layer height was set to 70 μm for all experiments and confirmed by laser scanning microscope measurements.
The specimen was placed inside a test chamber that was flooded with argon (Ar) inert gas. The chamber was attached to an axes system, which allowed a precise movement of the entire chamber through the focus of a stationary laser beam.
Moving the whole chamber relative to the laser beam and moving the beam relative to the work piece, as commonly done in L-PBF are equivalent and only differ in the reference coordinate system of the observer. A difference, which occurs regarding the process, is that the incidence angle is constant within our setup, whereas it changes depending on the position within the scan field for an industrial L-PBF machine. In the case of statistical measures, it is actually more desirable to investigate the influence in a set of experiments with a constant angle of incidence each, as is possible with the proposed setup.
A two-axis stage system was used. The coordinate system used within this paper was chosen relative to the camera setup, where the x-y plane describes the image plane and z is the depth along the observation axis. A fast stage can accelerate the whole setup up to usual scan speeds (Aerotech act115Dl-1000, = v 5 m s feed, max 1 ) along the x axis. A slower but more precise stage (PI M-605.2DD) was used for the realization of a defined hatch distance along the z axis.
The beam source used within the experiments was a Trumpf TruDisk 6001 Yb:YAG laser with a wavelength of 1030 nm and maximal output power of 6 kW. The end of the 100 μm fiber is projected by custom optics with an image ratio of 1:1 and a focal length of 200 mm, resulting in a 100 μm top hat spot on the specimen. A cross-section of the process chamber used is shown in figure 1.
The IZ was observed through two viewing windows. In order to adequately illuminate the particles above the IZ without illuminating resting particles within the powder bed, the process was illuminated using a flat angle. The optical axes of the two ultrahigh-speed cameras intersected within the process chamber, as shown in figure 2. Only trajectories of particles within the field of view (FOV) of both cameras can be reconstructed. Therefore, the measurement volume is formed by the intersection of the FOVs.
Two Phantom v1210 ultrahigh-speed cameras were used to observe the process at a frame rate of 60 kfps (@ 512× 256 px). The angle between the two camera axes was 30°, and both were inclined by 10°relative to the specimen plane.
The observed FOV is 20 mm×10 mm (compare figure 3), which results in a pixel equivalent of roughly 40 μm. Particles present within L-PBF can be smaller than this. Nevertheless, these particles are imaged on one pixel and are detectable under the premise that particle information within the image is higher than the image noise.
For illumination, a Cavilux HF illumination system was used. However, it was only used in every second frame causing an alternating illumination. This allowed the acquisition of both, illuminated and non-illuminated images from the same experiment with an effective frame rate of 30 kfps. This enables the classification of particles into hot and cold, since hot particles are also visible in non-illuminated images due to their emission of thermal radiation [16].
Ideally, the illumination is chosen to be in plane with the specimen, in order to reach particles above the specimen but not the powder layer itself. However, with the current setup, due to a high numerical aperture of our illumination optics, parts of the illumination light still reached the powder bed. Exemplary consecutive images with and without illumination from an experiment are shown in figure 3. One can see that some particles are visible in both images, identifying them as being hot particles, whereas a majority of the particles is only present in figure 3(b), classifying them as cold particles. Particles visible at the lower section of figure 3(b) are to avoid reflections from the powder bed.

Methodology-3D particle tracking velocimetry
In the following section, we present the algorithm which has been developed to calculate 3D spatter trajectories from stereoscopic images (e.g. figure 3) of the L-PBF. The flow scheme of our algorithm is displayed in figure 4. A similar approach was proposed in [15] for the tracking of particles in images generated by a three-or four-camera system. However, with a two-camera system and a high particle density, ambiguity occurs during the correspondence analysis, as two or more particle images can provide a potential match [17]. Several approaches are described in the literature to avoid or reduce ambiguity. One can use an additional camera in order to gain additional information for the correspondence analysis [17]. Alternatively, an increased frame rate leads to decreased particle movement within a time step, which reduces the complexity of the correspondence analysis within two consecutive measurements [14]. In addition, features of the particles can be used for their distinction [15].
Neither of the proposed approaches is applicable to our setup, due to the unavailability of another ultrahigh-speed camera, their limited frame rate, and the similarity of particles to be tracked. Therefore, in our case the algorithm is not only meant to reliably track the particles but also resolve ambiguity caused by the limited number of cameras.
The experimental setup used for the imaging was described in section 2. The pre-processing includes unpacking    of the packed data of the ultrahigh-speed cameras and noise reduction and low-pass filtering by means of a Gaussian filter, described in section 3.1. After this, the particle positions within the images are detected. By incorporating the information gained within the calibration of the camera system, a projection into 3D space of potential detected particles is possible. Within this step, ambiguity arises which is resolved in a following tracking step. After the tracking, the information gained about the particles is processed to acquire the desired measure (e.g. average particle velocity, number of particles, etc). The algorithm is described in detail in the following sections.

Pre-processing and detection of particles
In the first step, the particle positions need to be derived from the raw images (e.g. figure 3). The particle detector described in the following has been designed in accordance to the particle characteristics derived from the process observations distinguishing the particles from the background. Particles can be described as being spherical or close to spherical. In addition, particles can have a highly varying brightness, depending on the size, temperature, and reflectance, in the case of an illuminated scene. Furthermore, the effective particle size is, in all cases, in the range of several pixels. In the first step, the image noise and high frequency information within the raw image I raw is suppressed by the convolution of I raw with the 2D Gaussian-distribution G. The standard deviation σ can be estimated by the expected particle diameter in pixels. The filtered image, I, can then be written as shown in equation (3.1), with s and t being the coordinate relative to x and y: In the next step, the second derivative, in the form of the Hessian matrix, is calculated for each position (x, y) within the image. It describes the change of the slope in grayscale values in the pixel's neighborhood. ⎡ The local eigenvalues of the Hessian resemble the maximum and minimum of the second derivative at (x, y). Those can be calculated from H as follows: In the case of quasi-round particles, the second derivative needs to be uniform in all directions. In the present case of bright particles, the eigenvalues are negative [18].
The eigenvalues are used in an evaluation function delivering I Ev , whose value increases with the similarity of the image region to a round particle of size σ. By applying a threshold b thresh to the resulting image according to (3.6), each pixel is classified as being part of a particle or background (3.7).
each centroid of the formed regions is calculated. The result is a list of the detected particle coordinates from the left and right cameras. Figure 5 shows the result of the detector applied to the raw images shown in figure 3.
Within figure 5(b), a larger number of particles was detected compared to figure 5(a). A large number of detected particles seems to have rested within the powder bed. Nevertheless, the number of particles expected to be above the IZ is higher in figure 5(b), as expected.

Calibration
A photogrammetric approach is used for the calibration of the camera system. This means that the calibration is performed with a target of known geometry. We use a ceramic checkerboard pattern. The checkerboard pattern is moved by a selfbuilt calibration system and continuously monitored by the ultrahigh-speed camera. This results in very good coverage of the measurement volume with calibration points. The checkerboard with the detected pattern and the calibration system are shown in figure 6. The calibration was performed with the MATLAB ® 2018b built-in stereo calibration app, which is based on the methods described in [19,20]. With this method, an overall mean error of 0.15 pixels could be realized for our setup. The error describes the average deviation which occurs if a known checkerboard point is projected from one image into the other, using the calibration result.
In addition, the calibration provides the intrinsic parameters of each camera, which allows the compensation of image distortions caused by the lens or misalignment of the  image sensor to the optical axis. In addition to that, it delivers the extrinsic parameters, linking the cameras' coordinate systems to the world's coordinate system. The transformation of one coordinate system into the world coordinates can be fully described by a rotation matrix, R, a translation vector, t, and a scaling vector, s. By these, the two cameras' coordinate systems can be associated with each other, which is needed for the determination of 3D information, described in upcoming sections [19].

Matching of particles
The particle detector delivers the points detected in the scene imaged by the two cameras. Particles within the measurement volume (see figure 2) are imaged on both cameras. The two projections of the individual particle within the camera images are matched.
The calibration results can be expressed in form of the fundamental matrix F . fund It allows to validate a matching set of projections, consisting of a point within the left (x, y) and the right image (x′, y′) in form of the condition [21]: The set of points within the image plane P' of the right camera matching to a point u in the right image plane P, forms the epipolar line l u [21] and vice versa.
The given formula, however, is only valid for an ideal calibration and measurement, without any error. By also including candidates within a specific distance to the epipolar line, the errors can be accounted for. The distance threshold should be chosen within the order of the mean calibration error. The distance d between the epipolar line l u and the point  u can be analytically calculated by: This approach reduces the search area from the whole image to a narrow line. If multiple particles are lying on the same epipolar line the matching is not unique but ambiguous, since there are multiple possible matches that fit the condition. In figure 7, the epipolar lines are drawn into the image shown in figure 5(a) and the corresponding image of the second camera.
Due to the uncertainty, the number of possible matches is higher than the actual number of particles, leading to the detection of particles that are not actually present. We refer to those detections as ghost particles. Within the non-illuminated observations (e.g. figure 7), the number of possible matches was 25% higher than the number of particles detected within the single images. For the illuminated scenes, roughly 90%-95% are false detections.
In accordance with the description used for the particle detection, the particles resemble each other; therefore, no features can be used to differentiate them, as stated in [21]. Due to this, it is not possible to identify false pairs at this step within the algorithm. The use of additional cameras would allow the identification of the crossover points of the epipolar lines, reducing the search along a line to a single point [17].
However, we show, in the following sections, that the ghost particles can be separated from the real particles by a tracking algorithm, which incorporates a priori knowledge of the real particle movement.

3D coordinate determination
From the matches, we derive 3D information by means of triangulation. For this purpose, the detected points of both cameras are first projected into one common plane, equivalent to the case where the optical axes of both cameras are parallel to each other.
This coordinate transform is expressed in the homography matrices H hom and ¢ H , hom which incorporate the calibration results. The projection of a detection, u, to the position, u * , in the common plane is expressed as: 3 . 1 3 hom * From the disparity, d, of the projected point coordinates of both projections, the z coordinate can be calculated by means of triangulation: The baseline, b, and focal length, f are available from the camera calibration.
This step delivers the 3D coordinates ( ) P x y z t , , , n k of each matched pair of particles at each given time step, t k . However, the earlier mentioned ghost particles are still present within the data. Figure 9 shows the projection of the particles shown in the images in figure 5 into 3D space. One can see that the number of projected particles within the illuminated case ( figure 8(b)) was much larger than in the non-illuminated case ( figure 8(a)). This is caused by the large number of ghost particles. In addition to that, the detections within the powder bed were also reconstructed.

Tracking of particles-a priori knowledge
In order to distinguish ghost particles from the real particles, a priori knowledge of the process and the particles is incorporated into the tracking routine. I. The expected particle speeds can be taken from the literature. Ly et al reported spatter speeds, measured at the process level, between 6 and 20 m s −1 [13]. For large ejecta, expected to cause defects within the build process, speeds of 1-3 m s −1 were reported [6,12].
II. Since all driving forces for the generation of spatters and particle movement originate from the IZ, the particle movement, from a macroscopic observation, may be expected to be directed away from the IZ. III. Individual particles have to be detected in multiple frames since there is no 'particle sink' which makes them disappear. IV. The particles follow a quasi-linear movement over several frames. V. It has to be expected that the particles cannot be detected in every frame.

Tracking of particles-Kalman filtering
The a priori knowledge of the particle movement is incorporated in the form of a Kalman filter, first described in [22]. It allows the tracking of particles by means of a motion model (see IV). The Kalman filter is used to estimate the current state, x k (e.g. velocity, position) of an object based on previous measurements of the state x k−1 . The transition is expected to be described by the state transition matrix, F, and the process noise, v k . Consequently, the filter can be applied if the state model follows the equation: The available measurement can be related to the current state with the measurement matrix, H.
Since the process noise cannot be determined, the particle state,x k is estimated at the current point in time, t k , with consideration of the previous measurement and inaccuracies, all of which causing deviation to the real state, x k . With the previous estimated state at time, t k−1 , the prediction of the particle state,ˆ| -x , k k 1 within the current time step, t k is possible considering the state-transition matrix, F. The transition matrix F is specified by the applied motion model, which is, in our case, a movement with constant velocity. It is important to note that the motion model is only an approximation of the real motion. This means that particles not exactly following the motion model can also be tracked and described by:ˆ( The covariance, P k , of the prediction with respect to the expected measurement noise (measurement noise covariance R) and the expected accuracy of the motion model (process noise covariance Q) can be calculated as follows: The prediction delivers a state estimate, including the position of the particle at the time step, t k , and a measure for the expected accuracy of the estimation. If a measured particle can be found close to the estimated position at t k , the particle stateˆ| x k k can be updated with respect to the information gained by the measured values of the particle Z k at t k . The Kalman gain K k can be interpreted as the sensitivity of the Figure 8. Projection of an exemplary particle from 3D (purple) to the detection position (green) on the image plane (P P, P P¢) and the common plane (P P P P ¢ , * * ). update to the new measurement.
With this recursive approach, the precision of the particle state increases over time; and the predictions of the filter become more accurate.

Tracking of particles
From II, we can derive that the density of detected particles decreases with distance to the IZ. As discussed in section 3.3, it is obvious that the number of particles close to the epipolar line rises and more ghost particles appear with increased particle density. According to III, the routine needs to detect the particles until they leave the observation volume.
We consider the scene to play backwards, meaning that particles are flying towards the IZ. This means that their first occurrence is not within the IZ but close to the boundaries of the observation volume. There, the particle density is relatively low, which is the result of II. With decreasing distance to the IZ, the particle density rises. However, due to simultaneously increasing the number of included measurements over multiple time steps, the predictions of the Kalman filter become more accurate and it is less prone to misdetections.
The flow chart of the resulting algorithm is shown in figure 10. It is applied to the 3D particle detections (purple) from section 3.4.
First, the tracks are initialized. In our case, all combinations of particle detections over three consecutive time steps are identified as fulfilling I, II, and IV, forming multiple triplets of points. For this case only, particles with a certain distance to the IZ are considered, as their first occurrence is close to the boundaries of the observation volume. The data of each triplet is used to initialize a Kalman filter with state x 1 . An initial estimation for the covariance, P 1 , is done according to the user-driven parameterization of the algorithm.
Since we track backwards, the position in the previous frame is predicted for every track. Matching particles are appended, and the tracks are updated.
The seeding of new tracks is done for every time step. However, only particles which were not found to be part of previously initialized tracks are used for the seeding of new tracks.
The tracking is continued until no more valid particle positions can be measured for a given track over multiple time steps. Only tracks of a certain length are considered to be real tracks due to III, other tracks are disregarded as misdetections.
Particles which were not appended to a track are considered to be ghost particles.
At the end of the tracking procedure, we distinguish ghost particles from real particles due to their appearance within a valid track. Figure 11 shows an exemplary tracking result in which the IZ represents the origin of the reference coordinate system.

Speeding up processing steps by means of GPU processing
In order to analyze trends within the process behavior over a large parameter set, it is crucial to be able to process the raw data within an acceptable amount of time.
Therefore, all processing steps were optimized to run on the graphics processing unit (GPU) or in parallel tasks on the CPU. The system we used is equipped with 128 GB of RAM, two Nvidia Titan Xp GPUs, and a 16 Core AMD Ryzen Threadripper 2950X CPU.
We were able to speed up the particle detection within raw images from 20 to 500 fps by means of parallelization, with common MATLAB code and up to 6000 fps with an optimized GPU code. The 3D reconstruction achieved up to 2500 fps by means of parallelization in MATLAB. The tracking algorithm is not yet fully optimized. The main issue is the cubical rise of the number of seeding triplets with the number of particle detections, which makes the planning of memory allocations harder. However, within most experiments we can achieve 250 fps.

Validation
In order to validate our approach, we created a set of synthetic datasets with a known ground truth of particle trajectories. A specific number of individual particle tracks was placed within a 3D volume. The trajectories followed all criteria given in section 3.5 (I-V).
For each set of synthetic images, a mean velocity and a standard deviation along x, y, and z was set. According to the resulting normal distribution, the velocity components of each particle were randomly assigned. In addition to the abovementioned criteria, only particle velocities with a positive y component were considered as valid.
Based on the previous sections, we can conclude that the number of ghost particles and, therefore, potential misdetections increases with the number of particles within the image I of a single time step. Two main causes can be identified for this. First, a low average speed of the particles leads to a longer time of the particles being within the observation volume. Second, the number of particles ejected per time step is increased.
Even though this approach is not directly based on real experiment data but synthetic images, the particle properties are within the range of particles measured manually within real experiments and in accordance with the literature (I). This allows us to test the performance of our algorithm on a data set, which covers a large variety of particle properties. This cannot be achieved with manually labeled data. The overall number of synthetic images used for our validation was 1500 000 images, which is equivalent to 400 GB of experimental data.
Both cases were considered within the validation. First we increased the number of individual particles within the single experiment, keeping the particle velocities and standard deviations constant. Within the second test, we continuously increased the velocity. Figure 12 shows the validation results for particle numbers, as seen within our initial experiments. For the purpose of comparability to the results from the previous sections, the average number of particles per frame is also plotted (blue). The number of particles within one frame is an equivalent to the number of detections shown in figure 5.
In figure 12, four measures of the validation result are presented. The number of correctly detected particles (green) can be read as a measure for the information actually measured. The number is set relative to the overall number of real particles within the known ground truth. The number of ghost particles is a measure for the measurement noise.
The average ratio between ghost particles introducing noise to the measurement and real particles contributing to the measurement is 5%. However, looking at the deviation of the average speed taken from the measurement relative to the average speed within the ground truth (brown), the ghost particles do not have a high impact (below 1.5%).
However, by increasing the number of particles further, it is visible that the measurement error increased with the number of particles within the image, as expected. The results for a larger number of particles are shown in figure 13.
From figure 13, it follows that the accuracy will drop with an increasing number of particles. For the illuminated scene presented in figure 5(b), the to-be-expected ratio between ghost particles and real particles rose to 19%. This underlines the necessity of an illumination setup optimized towards the illumination of to-be-detected particles only.
The validation results for an increase in velocity are shown in figure 14. For low velocities, the number of visible particles rose cubically, as expected. A large number of particles also led to an increased probability of particles traveling with the same velocity. This led to the detection of ghost particles because the particles have similar epipolar lines over multiple time steps. However, the number of correctly detected particles rose for expected average particle velocities of 3-5 m s −1 .
Increasing the velocity further, the number of correct tracks decreased. However, ghost particles were not present. Both can be explained by the particles leaving the measurement volume in shorter time and, consequently, within less   frames. As a result, fewer particles per frame were visible, which decreased ambiguity. Since the particles were only visible within fewer frames, the algorithm was more prone to missing detections, leading to non-tracked particles. This can be counteracted by increasing the frame rate or the use of both illuminated and non-illuminated images for the detection of particles with a velocity larger than 10 m s −1 .
However, in both validation cases, the comparison of the detected average velocity with the actual average velocity shows good accuracy. In all presented validation experiments applicable to our measurement tasks, the error was below 5%.

Conclusion
Within our work. we have shown and validated a new approach to measuring 3D trajectories of hot and cold particles and their properties in the context of L-PBF by employing a 3D PTV approach. To the authors' knowledge, this is a methodology which has not yet been introduced in the field of laser materials processing.
Our validation has shown that the system delivers good results within the specified measurement task. However, the validation also revealed that the system is prone to measurement errors for velocities and particle numbers, exceeding those present in L-PBF. Nevertheless, the former can be avoided by an increased number of measurement points; and the latter can be resolved by an improved illumination setup, if needed. All in all, it can be concluded that the system is applicable to measuring particle properties within L-PBF.
With the availability of this new measurement technique, future work will focus on the particle and, especially, the spatter behavior within different process situations and materials. The results may then be used to quantify the evaporation behavior and its influence on the process. Furthermore, we plan to investigate the correlation between process regime and spatter characteristics in order to assess the applicability of the approach for quality control within industrial machines.
The system allows us to get a detailed view into the process behavior with a statistically sufficient sample size, literally accessing an additional dimension within the process observation of L-PBF. Increased process knowledge is needed to realize a faster process, higher quality, more resilient parts-pushing the envelope further.