Introduction

In the current Industry 4.0 framework, the use of image and video-image data in industrial process monitoring and quality modeling applications is becoming more and more common and widespread thanks to the quick development of discrete manufacturing process technologies and integrated sensing equipment (Barari et al., 2021; Colosimo, 2014, 2018a; Colosimo et al., 2018a, 2021; Megahed et al., 2011, 2012). In this framework, a major challenge regards how to make sense of big and complicated data to identify, model and analyze relevant information patterns. In many cases, the quantities of interest, also known as process signatures, exhibit spatial and/or spatio-temporal dependencies and variations whose proper characterization is the key to draw reliable conclusions about the process stability and the product quality.

In-line video-image analysis can be applied for different purposes. On the one hand, it can be used to make process optimization and process calibration more efficient. Indeed, in-line data may be used to reduce the experimental effort and the need for expensive and time-consuming post-process inspections, as response variables can be measured directly during the production process. On the other hand, in-line video-image data provide rich information for process monitoring, aiming at detecting process shifts and onsets of anomalies while the part is being produced.

In the recent years, many authors outlined the need for appropriate statistical modeling and monitoring methods for image data (Megahed et al., 2011, 2012; Qiu, 2005; Yan et al., 2017), but still limited attention has been devoted to video-image data, where the underlying evolution of the process in space and along time has to be properly captured and monitored. In this framework, we may identify two major streams of research. A first stream of methods consists of identifying regions of interest (ROIs) in each video frame and monitoring the evolution over time of their synthetic descriptors. This first route also represents the most common way to process video-image data in the industrial practice, especially when the engineering knowledge can guide the selection of the synthetic descriptors, i.e., the golden standard to be identified and counted in each frame. One main limitation of these approaches regards the information loss imposed by the selection of the quantities to be monitored. Another limitation is that, in some cases, the spatial spread of ROIs within the image can not be fully captured and modeled.

Another stream of research consists of modeling video-images as a sequence of frames, where each frame is represented as a matrix of pixel intensities. Following this route, Yan et al. (2015) compared different extensions of Principal Component Analysis (PCA) to detect shift of the spatial intensity pattern observed in a video sequence. A different method was proposed by Colosimo and Grasso (2018), where a so-called spatio-temporal PCA (ST-PCA) was used to enable a spatial localization of the occurred anomaly within the image, by combining a T-mode PCA implementation (Jolliffe, 2002) with a spatially weighted variance-covariance matrix definition. With the same aim of detecting an anomaly both in time and space, Yan et al. (2021) presented a penalized regression model, called spatio-temporal smooth sparse decomposition, combined with a control chart on the model parameters.

A limitation of this second family of methods is that they look at the frames of video-image data as matrices (or tensors) of pixel intensities, without frame pre-processing to look for specific elements or regions of the frame carrying detailed information about the process of interest.

This study focuses on modeling video-image data where the spatial signatures of interest can be represented as spatial point patterns (SPPs). An SPP data set is a nite collection of points specifying the locations of events within a given area of interest, which can be a 2-D, a 3-D space or even a manifold (Møller & Waagepetersen, 2017). Examples of SPPs can be found in various fields, from medical applications, e.g., spatial distributions of neurons (Jafari Mamaghani et al., 2010); to ecology, e.g., spatial distribution of trees for forest mapping (Velázquez et al., 2016); epidemiology (Diggle et al., 2005); astronomy (Stoica et al., 2007). In industrial applications, SPPs have been used to model the spatial distribution of fibers in composite or particles in nanocomposite manufacturing (Trias et al., 2007; Huang et al., 2015, 2017; Kam et al., 2013; Zhou et al., 2014; Dong et al., 2017a) or to study the porosity distribution in additively manufactured parts (Liu et al., 2019). Examples of industrial images exhibiting spatial point patterns are shown in Fig. 1. Despite a large literature on K-functions and SPP methods for image data, there is a lack of applications involving video-image data, although the dimensionality reduction enabled by K-functional representations makes such approach particularly attractive in video-image processing and analysis too.

Fig. 1
figure 1

Examples of industrial images where K-functions may be suitable to model spatial point patterns: SEM image of a metal matrix nanocomposite (left), TEM image of a polymer nanocomposite (centre) (Huang et al., 2015), and a spatial spread of pores in a metal specimen (right)

Full size image

In our study, the motivating application regards the statistical analysis of spattering in additive manufacturing, more specifically Laser Powder Bed Fusion (L-PBF) (Gibson et al., 2014).

L-PBF allows producing metal parts on a layer-by-layer basis by using a highly focused laser beam to locally melt a thin layer of metal powder. Different authors showed that spatters generated by the laser-material interaction enclose a rich information about the L-PBF process conditions and their stability over time (Zhang et al., 2019; Repossini et al., 2017; Andani et al., 2018; Bidare et al., 2018). Spatters are hot metal particles ejected from the melted region or the surrounding areas. Their spatial spread captured with a high-speed machine vision equipment can be regarded as an SPP. Indeed, at each given point in time, i.e., in each video-frame, the spatters can be represented as point particles whose spread in space is relevant to analyze the underlying process phenomena (more details are provided in the motivating example description in section “Motivating example”).

A technique suitable to synthesize the information enclosed in SPPs passing from 2-D spatial maps to 1-D curves is known as Ripley’s K-function (Ripley, 1977; Diggle et al., 2000). The K-function methodology allows one to describe the spatial spread of point processes and to distinguish between homogeneous, clustered and random spatial scatters of the analyzed quantities as they reflect into different shapes of the K-function curve. Ripley’s K-function has been applied to model nanoparticles (Dong et al., 2017a), The interpretability of resulting functional patterns represents a further benefit that makes Ripley’s K-function a suitable data modeling tool.

This paper presents a method which combines an image segmentation and processing step to transform video-image data of spatter ejections into Ripley’s K-function curves, and a curve fitting approach based on Functional Data Analysis (Ramsay, 2004), which enables the classification of different process behaviors via FANOVA. We show that by analyzing the statistical differences of K-functions associated to the spatial spread of spatters in L-PBF, it is possible to distinguish between process conditions that cause different types of volumetric defects (porosity) in the manufactured part. A comparison against benchmark methods that rely on computing synthetic spatter descriptors shows that the proposed approach is more effective in detecting actual departures from a stable and in-control process condition.

The paper is organized as follows: section “Motivating example” presents the motivating case study in L-PBF; section “Proposed methodology” describes the K-function methodology and our proposed FANOVA solution; section “Results” summarizes the main results in the spatter analysis application and section “Conclusions and future developments” concludes the paper.

Fig. 2
figure 2

Examples of cropped video-frames showing spatter ejections during the L-PBF process of maraging steel samples under lack-of-fusion (top panels) and over-melting (bottom panels) conditions

Full size image
Fig. 3
figure 3

Superimposition of centroid locations of all the spatters observed during the L-PBF production of one layer of maraging steel samples under lack-of-fusion (left panel) and over-melting (right panel) conditions

Full size image

Motivating example

The motivating example for the proposed approach regards the analysis of spatter ejections in L-PBF as a signature of the process behavior. Spatters, as many other relevant quantities, can be measured in-line and in-situ thanks to the layerwise paradigm that allows to observe (e.g., with high-speed cameras or thermal cameras) many process-related phenomena during the production of each layer (Colosimo, 2018a; Colosimo et al., 2018a; Grasso & Colosimo, 2017; Grasso et al., 2021; McCann et al., 2021). An increasing number of studies have been devoted in the last years to the possibility of using spatters and other L-PBF by-products (i.e., vapor gas and plasma emissions) for in-line and in-situ process monitoring (Repossini et al., 2017; Andani et al., 2018; Nassar et al., 2019; Zhang et al., 2018; Ly et al., 2017).

Spatter ejections are either caused by a vapor-driven entrapment of powder particles or by unstable solid-liquid transitions leading to molten material ejections (Liu et al., 2015; Khairallah et al., 2016; Ly et al., 2017). The interested reader is referred to (Grasso & Colosimo, 2017; Grasso et al., 2021; McCann et al., 2021; Kumar et al., 2022; Li et al., 2020; Tercan & Meisen, 2022; Liu et al., 2022; Zhang et al., 2022) for review studies devoted to these phenomena and to the development of process monitoring tools in L-PBF.

Various authors showed the effect of different process parameters on the spatter behavior (Repossini et al., 2017; Nassar et al., 2019; Zhang et al., 2018; Ly et al., 2017). They mainly used high-speed video imaging to measure spatter-related quantities and their evolution over time. The mainstream analytic approach consists of estimating, frame by frame, various synthetic descriptors, e.g., the number of ejected spatters, their distance from the melt pool, their size, velocity, etc., and relating them to the variation of process conditions. Some authors also proposed classification methods based on logistic regression (Repossini et al., 2017) or artificial neural networks (Zhang et al., 2018).

Figure 2 shows an example of high-speed (1000 fps) video frames acquired during the L-PBF of a maraging steel sample, showing the spatter ejected during the production of one single layer as bright (hot) particles on a dark background (more detailed information about the experimental settings are provided in Sect. 4). The laser beam is not visible because the laser wavelength is outside the camera’s sensitivity range. Top and bottom panels of Fig. 2 show two opposite out-of-control conditions, namely the spatters ejected when a too low energy input was provided to the material (top panel) and spatters ejected when a too high energy input was provided (bottom panel), respectively. The former condition is also known as lack-of-fusion, and the energy density is not sufficient to completely melt the metal powder, leading to large and irregular pores within the final part, and possible delaminations between adjacent layers. The latter condition is also known as over-melting, as the excessive energy density causes a local overheating and unstable solid-liquidus interfaces, with consequent formation of pores mainly characterized by a regular and spherical shape. The different effects of these two conditions on the final defectiveness of the part make their quick detection and identification extremely important to keep under control the process and aid the part qualification task. Figure 3 shows a superimposition of all the spatters detected in the two videos of Fig. 2 (lack-of-fusion condition on the left and over-melting condition on the right). Each spatter is identified by its centroid in Fig. 3. Figures 2 and 3 show that the spatial spread of spatters represents a suitable signature to identify different process states. Moreover, the evolution of spatters generated in L-PBF in space and along time can be regarded as a spatial point process, as the spatial mapping of the spatters within the image domain represents the feature of interest for process state classification. Therefore, SPP analysis can be particularly suitable to synthesize the relevant information content of the video-image stream, translating the image-based problem to a functional data analysis framework (Ramsay, 2004). It is worth noticing that this approach does not require any a-priori selection of synthetic descriptors, while focusing the attention on the actual spatial mapping of the spatter behavior.

Proposed methodology

Brief review of the Ripley’s K-function

The theoretical K-function can be defined as (Ripley, 1977; Diggle et al., 2000):

(1)

where \(\theta \) is the spatial density of points, i.e., the number of points per unit area. Therefore, \(\theta \, K(t)\) represents the expected number of points that are within a distance t around a randomly chosen point. This allows describing a point process at different distance scales.

Let U be an image of \(m \times p\) pixels, including n connected components that can be treated as points whose coordinates in the image domain are the coordinates of the connected component’s centroid. Then, an asymptotically unbiased estimator of the K-function K(t) is

$$\begin{aligned} K(t) = \frac{1}{n^2} \sum _{(x,y) \; \in \; U} w(x,y) \, {\mathbb {I}} \left( 0 < d(x,y) \le t \right) \end{aligned}$$
(2)

where \({\mathbb {I}} \left( \bullet \right) \) is the indicator function, d(x, y) is the Euclidean distance between points x and y, and w(x, y) is the edge correction factor.

The edge correction factor is applied to take into consideration the fact that circles of radius d(x, y) centered in x may be not fully included in U, depending on the closeness to the image border and the radius d(x, y). Ignoring edge effects may lead to a biased estimation of the K-function, especially for large values of t. The correction factor w(x, y) can be defined as \(w(x,y)= \frac{1}{P_{CIRC} (x, \, d(x,y) \, \mid \, U)}\), where \(P_{CIRC} (x, \, d(x,y) \, \mid \, U)\) is the proportion of the circumference of radius d(x, y) centered in x included into the image domain U. If the circle is completely inside the image domain, \(w(x,y)=1\). Different variants of the K(t) have been proposed, involving different edge corrections for boundary effects (Baddeley et al., 2006; Yamada & Rogerson, 2003). In this study, boundary effects do not represent a critical issue as the spatters are always located in a reduced area in the center of the image, making the edge correction not necessary. The interested reader may refer to Baddeley et al. (2006) and Yamada and Rogerson (2003) for a discussion on different applicable edge corrections.

Two other variants of the Ripley’s K-functions were proposed, involving a normalization operation that leads to a linear (L-function) or constant and equal to zero (H-function) expected value, respectively (Kiskowski et al., 2009; Baddeley et al., 2006). These two variants were mainly proposed to ease the qualitative analysis and interpretation of Ripley’s K-functions and to aid statistical hypothesis testing with respect to the null hypothesis of a homogeneous Poisson process.

K-functions are also connected to “pair correlation functions”, which may be obtained after differentiation and normalization of the K-function (Baddeley et al., 2006). However, the estimation of pair correlation functions involves the use of a kernel function, and the selection of the corresponding kernel bandwidth. As the K-functions avoid the need for this parameter selection, they are proposed in this study. However, the method can be easily extended to L-, H- and pair correlation functions.

K-functions have different properties that make them suitable to capture different patterns of spatial point processes. Indeed, the shape of K-function reflects not only the type of underlying distribution, but also the presence of clusters together with their within- and between-cluster distances at different scales. Figure 4, left panel, shows examples of SPPs obtained by generating 200 points whose coordinates were simulated as uncorrelated random numbers from a normal distribution with mean \(\upmu = \left( 0,0\right) ^{T}\) and \(\Sigma = \begin{bmatrix} 0.6 &{}\quad 0 \\ 0 &{}\quad 0.6 \end{bmatrix}\) or \(\Sigma = \begin{bmatrix} 1 &{}\quad 0 \\ 0 &{}\quad 1 \end{bmatrix}\), uncorrelated uniform numbers with parameters \((-1,1)\) and \((-2,2)\), random numbers from two overlapped clusters centered in (0, 0) but with different variances, and random numbers from two separate clusters. Figure 4, right panel, shows the corresponding estimated K-functions. The value of t at which the estimated K-function reaches the asymptote reflects the dispersion of the points, whereas the slope and changes of the first and second derivatives reflect whether the points are clustered together and their relative distances at different scales. In particular, the presence of clustered data leads to inflection points in the K-function. The K-functions also capture information about the number of points in a predefined area, since the asymptote is \(K(\infty ) = A - \frac{A}{n}\), being A the area of the image. Thus, being fixed A, the larger the number of points, the higher is the value of the asymptote. These features make the K-function suitable to capture salient patterns of SPPs including in one single curve information that could be captured via multiple simple statistics (e.g., number of points, average distance from the center, etc.), but also information that is difficult or impossible to fully describe via synthetic descriptors.

Fig. 4
figure 4

Left panel: examples of normal, uniform and clustered spatial distributions; right panel: corresponding estimated k-functions

Full size image

Methodology

A video can be represented as a temporal sequence of video frames \(U_1,U_2, \dots , U_j\) where the j-th video frame, \(U_j\), is an image of size \(A=m \times p\). The proposed approach involves an image pre-processing phase followed by the frame-by-frame estimation of the corresponding K-function. The pre-processing phase consists of (i) identifying, in each video frame, the \(n_j\) connected components representing the objects of interest and (ii) computing, for each of them, the centroid’s coordinates. We assume that the number of connected components, \(n_j\), may vary from one frame to another. Once the centroid’s coordinates of all connected components have been computed, the j-th K-function \(K_j(t)\) can be estimated using Eq. (2), for any pair of components (x, y).

Figure 5 shows an example of the sequence of steps to pass from one video frame where multiple connected components are present to the corresponding K-function estimation. The operation is then repeated for all video frames, as a parametric model of estimated functions can be fitted to enable statistical inference of K-functional patterns.

Fig. 5
figure 5

Flowchart of the proposed approach with an example of conversion from one original video frame to the K-function representation

Full size image

As far as the spatter analysis in L-PBF is concerned, each video frame consists of a dark background and different bright areas in foreground. The foreground areas correspond either to spatters ejected by the laser-material interaction or to the laser heated zone. The laser heated zone is usually the largest bright area in the image, due to the high radiation emission from the melt pool and the surrounding heated material. Each spatter can be described as a distinct connected component by converting the original gray-scale video frames into a binary format. To this aim, the approach proposed by Repossini et al. (2017), which foresees two steps, was applied. First, the Otsu’s thresholding algorithm (Chaki et al., 2014) was applied to binarize the images, such that the foreground pixels in each connected component have intensity equal to 1 and background pixels have intensity equal to 0. The optimal threshold value was identified by minimizing the intra-class intensity variance within the class of dark (background) pixels and the one of bright (foreground) pixels. In order to keep the same threshold value for all the frames, and to reduce the in-line computational effort, the optimal Otsu’s threshold was identified during a calibration phase, choosing the sample mean of threshold values identified in calibration frames. Such calibration was carried out by using a video-image data gathered in one of initial layers, before starting the data acquisition for the actual test phase. In principle, historical data from previous builds can be used to calibrate the binarization settings as well. Once binary frames have been obtained, the largest connected component in correspondence of the melted region was labeled as laser heated zone and not considered in next steps, to focus the analysis only on the spatters and their spatial spread. The spatter’s centroid was computed as the center of mass of the connected component. In this study a connected component labeling was applied with connectivity 8, i.e., such that each connected components includes pixels that are connected if their edges or their corners touch. All pre-processing operations were carried out using the Image Processing Toolbox of Matlab.

The segmented centers of mass were then used to estimate the Ripley’s K-function corresponding to a specific frame. It is also worth noticing that the K-function estimation and its use in our proposed approach do not rely on distributional assumptions. One of the properties of a K-function is that it has to be non decreasing, so a fitting method able to preserve this peculiarity must be used. As suggested by Ramsay (2004), a non decreasing function can be estimated using the equation:

$$\begin{aligned} K(t) = \beta _0 + \beta _1 \exp \left\{ \int \limits _{t_0}^{t} W(u) \, du \right\} \end{aligned}$$
(3)

where \( W(t) = \varvec{f}^{T}(t) \, \varvec{\alpha } \) is the usual regression term, whereas \(\beta _0\) and \(\beta _1\) are the constant and the tangent term of the model. Equation (3) describes a general set of curves with the constraint that they have to be non decreasing, i.e. the derivative approaches 0 if \(t \rightarrow \infty \). It should be noted that the terms \(\beta _0\) and \(\beta _1\) have the same meaning that they have in simple linear regression: \(\beta _0\) is the value that the function assumes when t is equal to 0 and \(\beta _1\) is the tangent at the same value. The model has to describe a generic function, so the basis functions of \(\varvec{f}(t)\) must be suitable to achieve a good curve reconstruction. To this aim, we advocate the use of a third degree B-splines basis functions. Once the basis in Eq. (3) has been selected, the parameters’ vectors \(\varvec{\beta }\) and \(\varvec{\alpha }\) are estimated by minimizing the sum of squares:

$$\begin{aligned} \min _{\varvec{\beta }, \varvec{\alpha }} \sum _{i=1}^{n} \left[ k_i - K(t_i) \right] ^2 + \lambda \int \limits _{T} \frac{d^2 K(t)}{dt^2} \, dt. \end{aligned}$$
(4)

The parameter \(\lambda \) controls the so-called smoothing penalty. If \(\lambda \) is equal to 0, the function in Eq. (4) is the classic sum of squares error function. In B-splines regression, a large number of basis function are commonly needed to estimate the underling curve, leading to possible undesired effects like undulation of the fitted curve or an ill-posed system of equations. To avoid these effects a small value of \(\lambda > 0\) should be set. It is also worth noticing that a curve of degree three or higher must be used to compute the smoothing penalty in Eq. (4). Figure 6 shows an example of estimated K-functions corresponding to each video frame during the production of one layer of one specimen in our real case study, and the corresponding fitted average curve \({\widehat{K}}(t)\). A four B-splines basis function with five equally distant internal knots and \(\lambda \) equal to \(10^{-5}\) was used.

Fig. 6
figure 6

Example of measured K-functions during all frames of a video and fitted mean function

Full size image

Once all the functions are fitted, they can be used to identify different conditions of the underlying process. Let \(\upmu _1(t), \upmu _2(t), \dots , \upmu _l(t)\) be the mean K-functions fitted starting from video-image data gathered in l different video sequences. The sample average curve associated to the i-th video, with \(i=1, \dots , l\) can be estimated as follows:

$$\begin{aligned} {\widehat{\upmu }}_i(t) = \frac{1}{n_i} \sum _{j=1}^{n_i} {\widehat{K}}_{ij}(t) \end{aligned}$$
(5)

where \(n_i\) is the number of frames in the i-th video sequence. In addition to the average curve, the corresponding \(95\%\) confidence band can be estimated as:

$$\begin{aligned} \begin{aligned} {L}_i(t)&= {\widehat{\upmu }}_i(t) - z_{0.975} \, {\widehat{\sigma }}_i (t) \\ {U}_i(t)&= {\widehat{\upmu }}_i(t) + z_{0.975} \, {\widehat{\sigma }}_i (t) \end{aligned} \end{aligned}$$
(6)

where \(z_{0.975} \approx 1.96\) is the 97.5% quantile of the Gaussian distribution and \({\widehat{\sigma }}_i (t) = \sqrt{{\widehat{\sigma }}_i^2 (t)}\) is the estimation of the standard deviation of the i-th function, with:

$$\begin{aligned} {\widehat{\sigma }}_i^2 (t) = \frac{1}{n_i - 1} \sum _{j = 1}^{n_i} \left( K_{ij}(t) - {\widehat{\upmu }}_i \right) ^2. \end{aligned}$$
(7)
Table 1 L-PBF process parameters
Full size table

Thanks to the computation of the confidence interval for each sample mean of fitted K-functions, it is possible to determine whether process conditions captured in different portions of the same process, and/or during the production of different parts, produced a statically different mean pattern of spatter ejections. It is also possible to compare the K-function representation of spatter patterns in any new layer or any new build against a reference pattern, representative of a stable and in-control process state. Once a statistically significant difference from the reference confidence band is observed, a departure from desired process conditions can be detected. Moreover, the shape of the sample mean of fitted K-functions in the current process state provides a rich and interpretable information about the type of occurred deviation. Because of this, the proposed approach can be used for in-line and in-situ monitoring of the L-PBF process, taking advantage of a reference pattern estimate carried out during a training phase. It can be also used to support process characterization and optimization relying on in-line measurement of spatter ejections as a signature of the process.

Results

Experimental settings

The experimental case study involves the production of simple specimens via L-PBF of 18Ni (300) maraging steel. By changing the process parameters, different melting conditions were obtained, corresponding to different severities and types of volumetric defects in the final parts. A high-speed camera was used to monitor the spatter ejections during the production of the specimens in different layers. At the end of the process, two non-destructive methods were used to estimate (i) the density of the samples (Archimede’s measurement method) and (ii) the nature of pores, especially their shape and location within the sample (X-ray computed tomography, CT). The experimentation was carried out by using an industrial L-PBF machine, namely a Renishaw AM250, equipped with a single laser scanning system that enables a point-wise scanning method. Along each scanned track, the laser moves from one point to another spaced a distance \(d_p\) apart with an exposure time \(\tau \) on each point. The purpose of the study was to study the capability of the proposed approach to properly identify statistically significant differences in spatter patterns induced by different energy inputs provided to the material. Six treatments, corresponding to six different energy density levels, were obtained by varying two controllable parameters, namely \(\tau \) and \(d_p\). The rationale behind the variation of both \(\tau \) and \(d_p\) is that acting on these two parameters allowed exploring a wider range of energy densities than controlling one single parameters, while preserving safe operating modes. Indeed, due to the pulsed melting mode implemented in the Renishaw AM250, varying both the exposure time and the distance between exposed points allowed inducing treatments ranging from severe under-melting (also known as lack of fusion) to over-melting condition. Each treatment was replicated three times by printing three specimens for each energy density level, leading to a total number of 18 specimens. All other controllable process parameters were kept fixed, apart from the scan direction that was varied layerwise according to the standard industrial practice. A stripe scanning strategy was applied, with a 67\(^{\circ }\) scan direction rotation every layer. The response variable was the spatial spread of the spatters captured by the high speed camera, which was then translated into the 1D K-function form. Table 1 shows the process parameters corresponding to each energy density level and Fig. 7 shows the specimen allocation within the build area. All the specimens had equal geometry (i.e., parallelepipeds of size 5 \(\times \) 5 \(\times \) 12 mm). The material used was a gas atomized powder with average particle size of 25–35 \(\upmu \)m.

Fig. 7
figure 7

Scheme of the specimen allocation within the build area; numbers indicate the energy density level for each specimen

Full size image

The specimens were placed in the middle of the build platform (Fig. 7), at a sufficient distance from each other to avoid thermal interference. The intermediate energy density level (namely, level 3) corresponds to the default set of process parameters suggested by the L-PBF system developer. Figure 8 shows the in-situ sensing setup used in this study. It consists of a high-speed camera (Olympus I-speed 3 with CMOS sensor) in the visible range (about 400 to 700 nm) placed outside the protective window of the build chamber. Video-images were acquired at 1000 frames per second (fps), in order to capture the spattering behaviour with a sufficient temporal resolution. A chessboard camera calibration was performed to estimate the spatial resolution in the powder bed plane, which was about 250 \(\upmu \)m/pixel. A total of 240 layers were printed, but in-situ videos were acquired only during six non-consecutive layers spanning the entire duration of the process. In each monitored layer, a video of the L-PBF of all 18 specimens was acquired.

Fig. 8
figure 8

High-speed video imaging setup

Full size image

As-built quality characterization

The effects of different process parameters on the volumetric defects of the specimens were investigated by measuring the overall density of each specimen and its porosity structure. The final part density was evaluated by means of the Archimede’s test and the results are shown in Fig. 9.

Fig. 9
figure 9

Density values measured by using the Archimede’s test

Full size image
Fig. 10
figure 10

Top view of the segmented pores using the CT results

Full size image

Figure 9 shows that energy density level 1 (30 kJ/cm) caused a severe lack-of-fusion condition, with a low density in the order of \(90\%\). Energy density levels 2, 3 and 4 (50, 80, 100 kJ/cm) produced the highest density. Finally, levels 5 and 6 (115, 130 kJ/cm) caused a higher part-to-part variability and an average density slightly lower than 98\(\%\). To gather additional information about the structure of the internal porosity, an X-ray CT scan of the specimens was performed, with a voxel size of 10\(\upmu \)m. Top views of the segmented pores are shown in Fig. 10.

To characterize the different porosity structures, a further descriptor was computed for each pore, namely the sphericity index. It was computed as:

$$\begin{aligned} \phi = \frac{\root 3 \of {36 \, \pi \, V^2}}{S} \end{aligned}$$
(8)

where V and S are the volume and the sphericity of the pore respectively. This index ranges from 0 to 1, where 1 indicates a perfect sphere. Figure 11 shows the \(95\%\) confidence intervals of the sphericity index for all the energy density levels. Specimens printed with energy density level 1 (30 kJ/cm) were characterized by large pores that were connected to each other. Due to the impossibility of separating connected pores in a reliable way, these pores were removed from the estimate of confidence intervals shown in Fig. 11. Moreover, one of the three samples printed with energy density equal to 30 kJ/cm detached from the baseplate during the process because of a severe delamination. Thus, the non-destructive characterization was performed only for two specimens at energy density level 1.

Fig. 11
figure 11

\(95\%\) confidence intervals for the sphericity index of each pore

Full size image

As shown in Fig. 11, all specimens printed with energy density levels 1 and 2 (i.e., energy density lower than 80 kJ/cm) exhibited a porosity structure characterized by irregularly shaped pores. This structure is typical of lack-of-fusion conditions, where insufficient energy input is provided to the material. Although the number of pores at level 1 was much higher than in level 2 (as shown also by the Archimede’s test), specimens produced with these two energy density levels shared the same internal porosity structure. On the other hand, all specimens produced with energy density levels 3, 4, 5 and 6 (i.e., energy density \(\le \) 80 kJ/cm) exhibited near-spherical pores mainly placed in the corners and along the borders of the specimen. A higher concentration of nearly spherical pores is typical of local over-melting conditions. Pores are mainly located along the contours of the part because of the inversion of laser motion at the end of each scan track, which causes locally slower scan speed and higher heat accumulation. This becomes even more severe in the corners where the length of adjacent scan tracks quickly drops leading to a shorter time between two laser scans in the same area and consequent over-heating effects.

Figure 11 clearly shows that the different energy density levels tested in this study produced two different kinds of porosity structures. In particular, there is a clear transition between the lack-of-fusion condition and a “plateau” characterized by a high density and a high concentration of near-spherical pores along the borders of the specimen. Optimal process parameters—in terms of volumetric defects—can be found at the beginning of such plateau, i.e. in correspondence of the energy density level 3. The energy density level 3 produced not only the lowest porosity in the manufactured specimens but also pore that were more spherical than all other tested conditions. On the one hand, the capability to determine the transition between the lack-of-fusion condition and the following plateau represents a key issue to predict the final quality of the part relying on in-line and in-situ data. On the other hand, the capability to detect statistically significant differences from a reference (stable and in-control) process behavior is necessary to design in-line and in-situ process monitoring tools. The analysis discussed in the following demonstrates the suitability of the proposed approach to this aim, where the spatter behavior observed with a volumetric energy density of 80 kJ/cm was used as a signature of reference process performance.

Spatter analysis via K-function modeling

In the present case study, the dataset consists of one video per monitored layer. Each video was split into six sub-portions, each of them capturing the spatial spread of spatters during the production of specimens with a given energy density. An example of the spatial spread of spatters in one layer for different energy densities is shown in Fig. 12.

Fig. 12
figure 12

Example of spatial spread of spatters in one layer for different volumetric energy density levels

Full size image

K-functions estimated in different layers produced with a given energy density can be considered as replicates of K-functions generated under the same process conditions. Figure 13 shows the fitted average K-functions for each layer and each energy density level (left panel) together with the \(95\%\) confidence band of the fitted grand mean K-function for each energy density level (right panel).

Figure 13 shows that energy density level 1 (30 kJ/cm) was characterized by a K-function pattern considerably different from the one observed at other energy density levels. At the lowest energy density, the grand mean of the K-functions was much faster than all other function in reaching the asymptote, which means that the spatial spread of spatters was smaller than in other cases. Moreover, the asymptote was significantly lower than in other K-functions, because of a much lower number of spatters. The grand mean of the K-function for energy density level 2 (50 kJ/cm) exhibited an intermediate pattern between the one at energy density level 1 and the one of all other energy density levels. The grand means of K-functions at energy density equal to or larger than 50 kJ/cm were characterized by overlapped confidence bands, which is consistent with the plateau condition discussed in Sect. 4.2.

Fig. 13
figure 13

Sample mean K-functions fitted in different layers and for different energy densities (a) and the corresponding grand mean K-functions with 95% confidence intervals (b)

Full size image

To better highlight and interpret the different patterns of the analyzed K-functions, the derivatives of the fitted grand mean K-functions were computed as:

$$\begin{aligned} K^{\prime }(t) = \frac{dK(t)}{dt} = \beta _1 \, W(t) \, \exp \left\{ \int _{t_0}^{t} W(u) \, du \right\} \end{aligned}$$
(9)

The curve in Eq. (9) is the derivative of the K-function estimated using Eq. (3). The additional analysis of the K-function derivative allows one to better highlight some local differences in the K-function patterns. Figure 14 shows the derivative of the fitted average K-functions for each layer and each energy density level (left panel) together with the corresponding \(95\%\) confidence band (right panel). Figure 14 confirms the existence of three distinct patterns, two corresponding to energy density levels 1 and 2, and another one consisting of a partial overlap of confidence bands for all other energy densities. Figure 14 also shows an inflection point of the derivatives of the grand means of K-functions. At energy density levels 1 and 2, such inflection point is more evident and it occurs at values t in the order \(t=25\text {-}30\) pixels. At higher energy density levels the inflection point is present but less evident, and if occurs at larger values of t, i.e., \(t=75\text {-}80\) pixels. This is caused by the fact that the spatter spread typically consists of most spatters close together in the neighborhood of the melt pool, and a smaller number of spatters with a kinetic energy sufficient to move away at higher distances from the melt pool. The value t at which the derivatives of the K-function exhibit an inflection point is representative of the scale at which the transition occurs between a majority of slower spatters and a minority of faster ones able to travel a longer distance.

Fig. 14
figure 14

Derivative of sample mean K-functions fitted in different layers and for different energy densities (a) and the corresponding grand mean derivative K-functions with 95% confidence intervals (b)

Full size image
Fig. 15
figure 15

Contrast plots with \(95\%\) confidence interval using as a reference the K-function pattern at energy density level 3, i.e., 80 kJ/cm

Full size image

In order to have a quantitative evaluation of the statistical differences among the K-functions associated to different energy densities, a functional variant of the analysis of variance can be applied, i.e., a hypothesis testing where the response variable is a function. It consists of testing the following hypothesis:

$$\begin{aligned} {\left\{ \begin{array}{ll} H_0: \upmu _1(t) = \upmu _2(t) = \dots = \upmu _6(t), \; t \in T \\ H_1: \exists \; (i,j) : \upmu _i(t) \ne \upmu _j(t), \; t \in T \end{array}\right. } \end{aligned}$$
(10)

where \({\widehat{\upmu }}_1(t),\ldots ,{\widehat{\upmu }}_6(t)\) are the grand mean K-functions associated to the six different energy density levels tested in this study. The reader is referred to Ramsden et al. (2007) for an overview of the theoretical background of statistical tests in functional data analysis. The resulting p-value of the hypothesis test in Eq. (10) was close to zero (\(p<0.001\)), confirming the statistical difference among the K-functions associated to different energy density levels. A further analysis consists of determining the differences between each energy density level and a reference condition, i.e., the condition that produced the lowest porosity in the manufactured specimens. In this study, the reference conditions consists of level 3, i.e., 80 kJ/cm. Such differences (also known as contrasts) can be estimated as

$$\begin{aligned} c_{ik} = \upmu _i (t) - \upmu _k (t) \end{aligned}$$
(11)

The contrast \(c_{ik}\) represents the difference between the grand mean of K-functions corresponding to the i-th energy density level and k-th level used as a reference, where \(k=3\). Figure 15 shows the contrasts and their corresponding \(95\%\) confidence bands. If the confidence band of a contrast in Fig. 15 includes the 0 for all values of t, it means that there is no statistical difference at the \(95\%\) confidence between the K-functions associated to the i-th and k-th energy density levels within the whole domain. On the contrary, if the confidence band is entirely on the positive or the negative side of the y-axis, it means that the K-functions associated to the i-th and k-th energy density levels are statistical different at the \(95\%\) confidence within the whole domain.

Figure 15a and b show that for energy densities levels 1 and 2, the contrast with respect to the reference level 3 is well above 0 at lowest distances, t, and below 0 at highest distances. This clearly reflects the type of departure from the reference spatter behavior at 80 kJ/cm under lack-of-fusion conditions. Compared to the the reference spatter behavior, most spatters produced with too low energy density concentrated at smaller distances above the heat affected zone (positive contrast at low t), and a only a lower number of spatters reached higher distances from the heat affected zone. Figure 15 also shows the different spatter behavior in the over-melting conditions (energy density levels from 4 to 6). In this case, the confidence bands of contrasts with respect to the reference level are below 0, and the higher is the energy density, the more evident is the deviation of the corresponding contrast towards negative values. At energy density level 4 (Fig. 15c), the confidence band of the contrast includes 0 in most of the domain, which indicates the lack of statistical difference with respect to the reference condition. At energy density level 5 (Fig. 15d), the upper bound of the contrast’s confidence band is below 0 for most values of t although being very close to 0. At energy density level 6 (Fig. 15e), eventually, the confidence band of the contrast is entirely below 0 and it shows a significant departure from the reference condition. This result highlights a gradual and significant departure of the spatter behavior from the reference pattern as the energy density increases. Such increasing deviations in the over-melting state correspond to a higher porosity of the manufactured parts. Therefore, the K-function patterns capture a clear transition from lack-of-fusion conditions and a plateau condition starting from 80 kJ/cm. In addition, they also capture a gradual and significant departure from the reference level 80 kJ/cm within the plateau as higher energy densities are provided to the material, especially at the highest energy density level, 130 kJ/cm. Therefore, the proposed method allows synthesizing the spatter behavior captured in video-image data is a way that is suitable to capture the major differences between different process states and their consequent effects on the final part porosity. Because of this, the proposed approach is applicable to detect departures from a reference (stable and in-control) process state in the framework of in-line and in-situ process monitoring via video-image data. The interpretability of the fitted K-function patterns represents another advantage of the proposed approach, as it allows the characterization of the salient differences among spatial spreads of spatters in every monitored condition, enclosing in one single curve multiple properties of the analyzed phenomenon.

Sensitivity analysis

One parameter that may affect the performance of the K-function representation of spatter patterns is the Otsu’s threshold used to binarize the original video frames. To determine the robustness of the proposed approach to this parameter, a sensitivity analysis was performed. All video frames acquired during the production of one whole layer were used to test different thresholds, ranging from 0 to 1. Figure 16 shows \(95\%\) confidence intervals of the mean number of spatters identified as connected regions in video frames binarized using different thresholds. Figure 16 shows that too small values of the threshold introduce an anomalous peak in the number of detected components caused by the fact that several dark pixels of the background regions are wrongly assigned to the foreground class. For threshold values larger than 0.2, an increase of the threshold has a minor effect on the number of identified spatters, and the number drops only for very high thresholds, namely above 0.98. By testing our proposed approach with different threshold values, we observed that no significant change in the ANOVA results occurred within a wide range of values, roughly between 0.23 and 0.73. The optimal threshold used in this study, as a result of the calibration phase, was 0.42, which is well within this range. The amplitude of the range makes the proposed approach robust to uncertainty in video-image data pre-processing settings. Generally speaking, we advocate a calibration phase as the one described above to identify image binarization settings that are suitable to generate robust and consistent results.

Fig. 16
figure 16

\(95\%\) confidence intervals of the mean number of spatters identified as connected regions in video frames binarized using different thresholds from 0 to 1

Full size image
Fig. 17
figure 17

95% confidence interval of contrasts with respect to energy density level 3: number of spatter (left panel) and area of the convex hull including all spatters (right panel)

Full size image

Comparison analysis

A comparison analysis was carried out to investigate the benefits of the proposed approach for spatter analysis against benchmark methods mainly used in the literature. As discussed in Sect. 2, the most common approach adopted in the literature for the analysis of spatter ejections in L-PBF and for process state classification using spatter-related information involved the computation of synthetic descriptors, i.e., one or multiple scalar values associated to each spatter. It is also worth noticing that currently there is no industrial solution suitable to analyze spatter patterns. Indeed, most industrial L-PBF systems are equipped with sensors (including machine vision devices) for in-situ monitoring, but the focus is on the powder bed homogeneity and on the stability of melt pool emissions. Despite of the different target of the analysis, the best industrial practice for in-situ image data analysis in L-PBF consists of computing either pixel-wise variation metrics or synthetic descriptors associated to regions of interest. Because of this, the competitor approach considered in this study is not only representative of the mainstream literature on spatter analysis and monitoring in AM, but also of the best industrial practice applied to off-the-shelf monitoring toolkits. Two synthetic indexes used by different authors are the number of spatters and the area of the convex hull including all the spatters in one frame. The area of the convex hull is a synthetic measurement of the spatial spread of the spatters. More details about their computation can be found in Repossini et al. (2017). The two indexes were computed in every frame and their average values were then used as an alternative to the computation of sample mean K-functions. Figure 17 shows the \(95\%\) confidence intervals for the contrasts between the number of spatters and the convex hull area computed under different energy density levels and the reference level 3. Thus, Fig. 17 represents the analogous of Fig. 15. Figure 17 shows that both the synthetic indexes allow detecting a statistically significant difference between the two lack-of-fusion conditions (energy density level 1 and 2) and the reference level. However, they fail in detecting a statistical difference in the spatter behavior between overmelting conditions (energy density levels 4 to 6) and the reference level 3. The adoption of the K-function to characterize the spatial spread of spatters, instead, highlighted a statistically significant difference between the most severe over-melting state and the reference, which was also confirmed by a higher porosity in the part.

Figure 17 shows that synthetic indexes may be effective in detecting most severe departures from a reference condition, like the two lack-of-fusion states observed in this study, but they are not enough sensitive to less severe deviations. Moreover, they capture a poorer information content compared to K-functions, which results into a potentially lower interpretation capability about the monitored phenomena. More specifically, the K-function methodology allows capturing variations of the underlying pattern at different scales, i.e., at different distances t within the image, in a way that is different or even impossible to represent by means of discrete synthetic indexes. These results highlight the added potential enabled by the proposed method for in-line and in-situ analysis and monitoring of spatial point patterns in video-image data, and its possible application to L-PBF.

Conclusions and future developments

The increasing availability of image and video-image data for statistical quality modeling and monitoring applications is making more and more urgent the need for novel methods suitable to deal with real industrial problems. In this study, we presented a methodology that allows summarizing the spatial spread of multiple regions of interest in video-image data by passing from a 2-D image to a 1-D profile data format relying on the Ripley’s K-function approach. More specifically, this study addressed the family of problems where the time-varying spatial signatures of interest can be represented as SPPs, where each region of interest can be treated as a point object. The application of the proposed approach in L-PBF showed that the data synthesis based on the K-functions representing the spatial spread of spatters allowed distinguishing between different process states that led to volumetric defects in the final part. In-situ video-image data analysis via the Ripley’s K-function can be used to detect process shifts associated to a change in the spatter ejection behavior, which can be a possible driver of defects in the part. This approach also represents a viable way to reduce the dimensionality of big data gathered during the L-PBF process, with consequent benefits in terms of data storage and management requirements. Indeed, it allows passing from terabytes per build, needed to store high-speed videos gathered in every layer, to few megabytes needed to store K-function data. We showed that it is more effective than competitor methods commonly used in the literature, which consist of analyzing pre-defined synthetic indexes. The method is also suited to be used in a process optimization framework, where in-line spatter measurements may support the identification of optimal processing windows, with the further advantage of reducing the experimental effort and the costs associated to post-process inspections. In summary, the major strength points of the proposed methodology can be synthesized as follows: (i) enhanced in-situ detection of unstable process conditions: spatters are known to be a relevant signature for in-situ monitoring of the L-PBF process, and the proposed approach enables a more effective detection of departures from a stable process state than competing methods adopted in the mainstream literature; (ii) enhanced in-line data modeling for efficient process optimization, by replacing post-process inspection with in-line measurements, (iii) reduction of data handling and storage needs thanks to the massive data reduction enabled by the K-function modeling, moving from high-speed/high-resolution video-image streams to 1D functional representation of salient information. Two possible extensions of the proposed approach can be envisaged. One regards the extension from modeling spatial patterns enclosed in the data stream to modeling the spatio-temporal patterns. The use of K-functions to this aim was discussed by some authors (Hohl et al., 2017). It includes a temporal scale in addition to the spatial scale in the definition and estimation of the K-function. Such extension may be suitable to model time-varying spatial patterns.

Another extension regards the possibility to include additional descriptors of interest associated to each particle, e.g., area, shape, etc., into the SPP analysis. Variants of the K-function and point process methodologies have been proposed to this aim (Comas et al., 2011). They can be considered as a valuable extension when the description of the phenomenon requires other features of interest associated to each spatial location.

It is worth noticing that although spatters move in a three-dimensional space, a monocular vision like the one used in this study allows detecting only their apparent location in the 2-D image plane. Some recent papers (Eschner et al., 2019; Barrett et al., 2018) demonstrated the feasibility of a 3-D spatial spatter localization in L-PBF, although it requires more complex machine vision equipment and measurement settings. Albeit the 2-D spatial mapping of spatters, as discussed in this paper, has been shown to be sufficient for the classification of different energy density conditions, the proposed methodology can be further extended to 3-D spatial maps possibly gathered through high-speed stereo vision systems for a deeper analysis of process stability.