Toward a Real-Time Analysis of Column Height by Visible Cameras: An Example from Mt. Etna, in Italy

: Volcanic plume height is one the most important features of explosive activity; thus, it is a parameter of interest for volcanic monitoring that can be retrieved using different remote sensing techniques. Among them, calibrated visible cameras have demonstrated to be a promising alternative during daylight hours, mainly due to their low cost and low uncertainty in the results. However, currently these measurements are generally not fully automatic. In this paper, we present a new, interactive, open-source MATLAB tool, named ‘Plume Height Analyzer’ (PHA)


Introduction
Multiple volcanic hazards are associated with tephra dispersal [1][2][3], which encourages volcanological observatories to permanently improve their monitoring systems with the aim of tracking the main features of an explosive eruption [4][5][6][7]. Eruption column height is one of the most important source parameters for volcanic monitoring purposes [8,9]. In fact, this parameter is reported in the VONA (Volcano Observatory Notices for Aviation) messages issued in real-time by volcano observatories when an ash-producing event occurs and/or when there is a change in volcanic behavior [10]. The VONA messages are usually sent by fax or email by the observatory to the pertinent Area Control Centre, Meteorological Watch Office, and Volcanic Ash Advisory Centre [10]. Plume height estimation is also essential to evaluate the mass eruption rate of an explosive event [11][12][13][14], and represents a critical factor in forecasting volcanic ash dispersion [15][16][17] through numerical modeling [18][19][20][21][22]. Moreover, the level reached by the volcanic plume is essential in some gas plumes and aerosol retrievals by satellite [4]. Eruption column height can be obtained using different remote sensing techniques, including satellite [23][24][25], thermal or visible

Camera Network of INGV-OE and Dataset
The surveillance network of Etna, whose data are managed by the INGV-OE, includes seismic and infrasonic stations, tiltmeters, GPS (global positioning systems), strainmeters, UV scanners, and thermal and visible cameras [48], among other instruments installed for monitoring purposes. In particular, the visible camera network of INGV-OE includes two low-cost visible cameras: EBVH (Etna Bronte Visible High-definition camera) and ECV (Etna Catania Visible camera), which sit in the west flank of Etna near the town of Bronte and in the southern flank of Etna in Catania, respectively. These cameras (model VIVOTEK IP8172P) present a maximum resolution of 2560 × 1920, with a field-of-view of 33 • -93 • (vertical), 24 • -68 • (horizontal), and 40 • -119 • (diagonal). The cameras are calibrated in terms of location and orientation [49]. With this information, assuming that the plume has a negligible depth and is confined to a vertical plane that rotates according to the wind direction, it is possible to manually estimate the height associated with each pixel of the images and thus calculate the plume height from the record of visible cameras [4], with uncertainties of the order of 0.5 km. Additional details can be found in Scollo et al. [4,49]. Specifically, the dataset considered in this work consists of videos recorded by the ECV static camera for different explosive eruptions (Table 1). This dataset includes both eruptions with optimal visibility conditions and eruptions where column height is not detectable, even manually (Table 1). The images are particularly dark. Plume height is beyond the measurement limit during the whole video The video includes a period with plume height beyond the measurement limit. 1 Video employed to complement the set of calibration images used for the construction of an operative calibration file, but it was discarded for the analysis of plume height versus time because this parameter is not considered informative in this case. 2 Video employed to complement the set of calibration images used for the construction of an operative calibration file, but it was discarded for the analysis of plume height versus time due to visibility limitations.

PHA (Plume Height Analyzer)
is an open-source MATLAB tool that is able to analyze images and videos of volcanic eruptions (derived from visible cameras), with the objective of detecting the volcanic plume and calculation of the temporal evolution of plume height. The Image Processing Toolbox of MATLAB is required for launching the program PHA, whose graphical interface includes six sections (see Table 2). PHA is a self-customizing model, which means that the user must perform an iterative calibration procedure based on the analysis of images of previous eruptions before the model can be used to automatically detect plume height from new images. The underlying, final goal of this approach is to create a 'universal' functional algorithm that would automatically work for new eruptions (without further calibration), thus dealing with different eruption, atmospheric, and illumination conditions. Such a functionality would open the doors for us to implement real-time procedures to compute the column height of ongoing eruptions from the analysis of visible cameras.

Section Function Description
Fixed mask Load Load a previously created fixed mask.

Fixed mask Create
The program displays a graphical interface, where the user can draw a fixed mask on a reference image (selected by the user). This information is then saved in the folder MaskFiles.
Fixed mask Plot Plot the reference image and fixed mask.
Vent position Load Load a previously created vent position.

Vent position Create
The program displays a graphical interface, where the user can select the vent position on a reference image (selected by the user). This information is then saved in the folder VentFiles.

Vent position Plot
Plot the reference image, vent position and fixed mask.
Pixel to height Load Load a previously created pixel-height conversion matrix.

Pixel to height Create
Three modalities to create a pixel-height conversion matrix are available: -Constant, vertical gradient: the user is asked to indicate the height associated with the vent and with the top of the reference image.
-Bilinear interpolation: a graphical interface is displayed, and the user is asked to select a set of pixels of the image and indicate their heights. The resulting conversion matrix is computed by fitting this information as a function of pixel position, using a bilinear interpolation.
-Second-order interpolation: a graphical interface is displayed, and the user is asked to select a set of pixels of the image and indicate their heights. The resulting conversion matrix is computed by fitting this information as a function of pixel position, using a second-order interpolation. The resulting pixel-height conversion matrix is saved in the folder PixelHeightFiles.
Pixel to height Plot Plot the reference image and the isolines of height, derived from the pixel-height conversion matrix.
Calibration: Lab mask Load Load a previously created calibration function.

Calibration: Lab mask Create
An interactive, iterative procedure is launched, which samples frames from a calibration dataset provided by the user (a single video, a folder containing videos, or a folder containing images) and shows the application of different threshold values in the Lab mask (see Section 2.2.2). The user is asked to indicate the best conservative threshold value. Once the iteration is finished by the user, the program computes the functions used to calculate the Lab mask threshold as a function of the image properties (see Section 2.2.2). The resulting function is saved in the folder CalibrationFiles/LabMask.

Calibration: Lab mask Improve
This routine reproduces the same process associated with the creation of a Lab mask calibration, but the information is added to an existent Lab mask calibration. Since the performance of this mask is strongly controlled by the amount of data that the calibration includes, this function allows to improve the program performance for a given static camera.
Calibration: Lab mask Merge This routine allows to merge the data contained in two or more existent Lab mask calibrations and creates a new, likely better calibration function.

Calibration: Lab mask Test
An iterative procedure is launched that shows the application of the Lab mask, whose threshold is computed with the loaded calibration function, on a set of frames selected by the user (a single video, a folder containing videos, or a folder containing images).

Section Function Description
Calibration: Lab mask Compare The user is asked to select two or more calibration functions (in order to compare them) and a set of frames (a single video or a folder of images). Two modalities are available: -Plot Threshold: the program calculates the Lab mask thresholds of the selected frames considering the different calibration functions. Then, this information is plotted.
-Show images: an iterative procedure is launched that shows the application of different Lab masks, whose thresholds are computed using the selected calibration functions, on the set of frames indicated by the user.

Calibration: Default Parameters Load
Load a previously created set of default parameters. Even though the results presented in this paper consider the same set of default parameters, this section allows to increase the applicability field of this code.

Calibration: Default Parameters Create
A window is displayed, where the user can change some of the constant parameters used in the code (e.g., maximum number of clusters in Lab mask function). Even though the results presented in this paper consider the same set of default parameters, this section allows to increase the applicability field of this code.

Analysis Single video
This function allows analyzing a single video. The user is asked to provide the frame step adopted to analyze the video and the time interval between two consecutives frames in the video. The output is a plot of the temporal evolution of plume height, and this information can be saved in the folder Results.

Analysis Folder with images
This function allows analyzing a folder containing images. The user is asked to provide the time interval between two consecutives images. The output is a plot of the temporal evolution of plume height, and this information can be saved in the folder Results.

Analysis
Analyze manually This function allows to select manually the pixel associated with the maximum height on a set of frames (a single video or a folder of images). The output is a plot of the temporal evolution of plume height, and this information can be saved in the folder Results.

Single plot
This function allows plotting the results of previously analyzed videos/folders with images. The input of this function is the output file saved by any of the three functions of the section Analysis.

Compare plots
This function allows comparing the results of previously analyzed videos/folders with images. The inputs of this function are the output files saved by any of the three functions of the section Analysis.
For each frame analyzed (e.g., Figure 1a), the general procedure to compute the plume height consists of a series of successive steps ( ing those in contact with the plume (Figure 1d-f). (d) A procedure to evaluate the internal variability of the non-masked zone of the images and eventually exclude the low-variability zones. (e) Finally, considering a pixel-to-height conversion matrix (Figure 2), the highest pixel belonging to the plume is identified.   Table 1) [4]. Winds blowing to the E (to the right in the images) translate into height isocurves dipping to the W, while when winds blow to the W (to the left in the images), the resulting height isocurves dip to the E. Data are presented in m a.s.l.
In this section, we present in detail each one of the described steps, together with all the interactive functions that PHA includes for their development (see Table 2).

Fixed Mask
Given a static camera, a fixed mask is introduced (Figure 1b) in order to discard all the pixels of the images that are associated with infrastructure and volcano topography. This is important to avoid processing pixels whose properties are not associated with the eruption and/or atmospheric characteristics. Since this process must be performed only once for each static camera, the best way to define this mask is manually. For this, PHA includes ad hoc functions that allow creating interactively, saving, loading, and plotting the fixed mask (Table 2).  Table 1) [4]. Winds blowing to the E (to the right in the images) translate into height isocurves dipping to the W, while when winds blow to the W (to the left in the images), the resulting height isocurves dip to the E. Data are presented in m a.s.l.

Lab Mask
In this section, we present in detail each one of the described steps, together with all the interactive functions that PHA includes for their development (see Table 2).

Fixed Mask
Given a static camera, a fixed mask is introduced (Figure 1b) in order to discard all the pixels of the images that are associated with infrastructure and volcano topography. This is important to avoid processing pixels whose properties are not associated with the eruption and/or atmospheric characteristics. Since this process must be performed only once for each static camera, the best way to define this mask is manually. For this, PHA includes ad hoc functions that allow creating interactively, saving, loading, and plotting the fixed mask ( Table 2).

Lab Mask
We have observed that the third channel of an image in Lab scale (i.e., the color dimension b) generally discriminates well between clear sky pixels and other zones of the processed frames (e.g., volcanic plume and clouds; Figure 2c). The suitability of this channel to recognize sky pixels under optimal illumination conditions is similar to that observed by considering the subtraction of the red and the blue channels of the image in RGB scale, which is used in the software PlumeTraP [6], while it is slightly better during sunrise and sunset. For most frames of the studied videos, we have observed that a threshold value for the third channel of the image in Lab scale can be set in order to create a mask to split the image in two regions (hereafter, the Lab mask). However, it is not possible to set a single threshold value that works for a large set of images. Instead, it depends on specific image characteristics, which are in turn controlled by atmospheric factors as well as eruption and illumination conditions. Due to this, we developed an iterative procedure (hereafter, the Lab calibration) to calculate a function that, for each analyzed frame, provides the threshold (T) used to compute the Lab mask. This calibration can be performed using a reference video, a folder of videos, or a folder of images, and consists of an iterative procedure where different frames are sampled and analyzed, and the user is asked to indicate the best conservative mask (i.e., not masking the plume) from a choice of nine alternatives ( Figure 3). Let us consider a calibration procedure based on N frames. For each frame (i = 1, . . . , N), and considering the portion of the images not masked by the fixed mask, the algorithm saves the following information: on specific image characteristics, which are in turn controlled by atmospheric factors as well as eruption and illumination conditions. Due to this, we developed an iterative procedure (hereafter, the Lab calibration) to calculate a function that, for each analyzed frame, provides the threshold ( ) used to compute the Lab mask. This calibration can be performed using a reference video, a folder of videos, or a folder of images, and consists of an iterative procedure where different frames are sampled and analyzed, and the user is asked to indicate the best conservative mask (i.e., not masking the plume) from a choice of nine alternatives ( Figure 3). Let us consider a calibration procedure based on frames. For each frame ( = 1, … , ), and considering the portion of the images not masked by the fixed mask, the algorithm saves the following information:  Illustrative example of the iterative procedure used to define the threshold value of the Lab mask, necessary to recognize the pixels belonging to the volcanic plume as a function of the image properties. In these images, the lightened areas correspond to pixels potentially considered as part of the ash plume. In this case, the recommended choice is E or F.
The properties of the -th calibration image are thus given by = { , , , , , } and the information provided by the user, by selecting the best conservative mask, is given by . Using and distance-based criteria, calibration images are automatically clustered. The number of clusters ( ) clearly depends on : the larger the set of calibration . Illustrative example of the iterative procedure used to define the threshold value of the Lab mask, necessary to recognize the pixels belonging to the volcanic plume as a function of the image properties. In these images, the lightened areas correspond to pixels potentially considered as part of the ash plume. In this case, the recommended choice is E or F.
The properties of the i-th calibration image are thus given by and the information provided by the user, by selecting the best conservative mask, is given by T i . Using X i and distance-based criteria, calibration images are automatically clustered. The number of clusters (M) clearly depends on N: the larger the set of calibration images, the higher the number of clusters. The use of weighted, distance-based clustering techniques allows automatically generating classes of images characterized by common eruption, atmospheric, and illumination conditions. In this procedure, the distance between the different frames is computed by adopting the vector X i associated with each of them to define their positions. For each cluster, a specific function is derived to calculate the threshold (T) used in the Lab mask: where j is a subscript referring to the cluster (i.e., j = 1, . . . , M). F j is defined using a polynomial fit computed by adopting the subset of calibration images that defines the j-th cluster. The order of this polynomial fit depends on the amount of data available in this cluster (first order regression for clusters with 20 data or less, and second order regression for clusters with more than 20 data). We also assume that F j (X) is bounded by the minimum and maximum thresholds within the calibration images.
Thus, when a new image is analyzed, characterized by the vector X n , the code finds its cluster by computing the minimum distance between X n and X i (i = 1, . . . , N), and then it adopts the function F j associated with this specific cluster to calculate F j (X n ). On the other hand, a calibration function defined by the nearest calibration frame (in terms of the image characteristics, that is, in terms of the vector X i ) is also present in PHA. To present conservative results, in this work we use the more conservative choice between both the alternatives for each frame analyzed.
The number of iterations determines the spectrum of images that the code will be able to analyze correctly. Since this is a critical factor in the effectiveness of the presented procedure, PHA includes a set of functions that allows creating and interactively improving the calibration data, as well as for loading, comparing, and testing them ( Table 2). We remark that the only relevant correction needed by the Lab mask is associated with dark zones of the plume during sunrise and sunset, but this is automatically solved with no need of calibration. Additionally, note that this iterative procedure allows the user to indicate the images where the plume is not recognizable, which permits the code to automatically identify the images where the estimation of plume height is not possible using visible cameras (e.g., cloudy conditions, night images).

Cloud Identification
At this point, in general, the algorithm has masked all but the plume and clouds. Three successive procedures are then employed to discard the pixels associated with clouds: (a) We trace a large number of segments between border points of the image (above the vent position) including both horizontal and inclined segments. When a line intersecting a completely masked zone is identified (i.e., with no plume or clouds), the entire region above this line is masked, reducing the computation time and discarding pixels associated with clouds ( Figure 1d).
(b) Then, all the not masked pixels are clustered considering a distance-based criterion, and only the nearest cluster to the volcanic vent is conserved for the following steps ( Figure 1e). (c) Finally, the perimeter of the non-masked region is studied to identify lobe-like geometries. When the distance (calculated through a line) between two points in the non-masked region border is much lower than the distance calculated through the perimeter of this region, these points are assumed to define a lobe-like geometry, and this part of the non-masked zone is discarded. In this way, clouds superposed to the plume tend to be discarded (Figure 1f).
Since this procedure needs the volcanic vent position as an input parameter, PHA includes a set of functions that allow to interactively set the position of the vent (note that the definition of vent position must be performed only once for each static camera). The program allows plotting and loading this information as well ( Table 2).

Internal Variability of Non-Masked Zone
At this point, the algorithm has identified a set of pixels that represent a good candidate for the volcanic plume. However, sometimes the limits of the plume are diffuse, and thus we need an additional criterion. Since color variability tends to be higher in the plume, we adopted a threshold to consider the portion of the mask characterized by a large color variability. A fixed threshold value has works for all the videos studied here, and thus calibration is not needed. In any case, PHA includes proper functions to set different values of threshold when it is required (e.g., for other volcanoes or other static cameras).

Pixel to Height Conversion
In order to calculate plume height, we need a procedure that is able to relate pixels to plume height. This is represented by a conversion matrix with the same dimensions of the images, indicating the height associated with each pixel of the image. In this way, the model is able to evaluate the height of all the pixels belonging to the plume and determine the plume height (i.e., the maximum height associated with a plume pixel). While PHA offers different alternatives to create a conversion matrix (see Table 2), in this work we use specific conversion matrixes associated with the wind fields observed during some of the eruptions described here (see Section 2.1 and Figure 2). The frames studied present a maximum measurement height of the order of 10 km a.s.l. due to limitations in the field view, which is influenced by wind direction and intensity as well. Whereas winds blowing to the E (to the right in the images) translate into height isocurves dipping to the W, when winds blow to the W (to the left in the images), the height isocurves dip to the E. On the other hand, winds blowing to the south (i.e., directly towards the ECV camera) produce more spaced height isocurves (i.e., reduction in the measurement limit) and winds blowing to the north translate into less spaced height isocurves (i.e., increment of the measurement limit). Specific details about the geometric considerations needed to define these pixel-to-height conversion matrixes can be found in Scollo et al. [4,49].

Results
Once the information of each frame is computed, PHA plots the temporal evolution of height plume, excluding results considered outliers within the temporal series. Results are also saved in MATLAB files that can be imported, plotted, and compared later using the program PHA.

Internal Calibrations
As a first step, we show examples of internal calibrations, i.e., we use a few frames of single videos to create a calibration file and then we apply it to analyze the same videos. We focus on three short videos with different degrees of complexity (V18, V08 and V21; see Table 1) and a long-lasting video with significant changes in eruption and illumination conditions (V30; see Table 1). Note that all the calibration files and results described in this paper are included in the Git repository https://github.com/AlvaroAravena/PHA (accessed on 12 January 2023).
3.1.1. V18 (18 March 2012) This eruption is part of the 25 low-explosivity events observed at Etna between January 2011 and April 2012 [50]. V18 is characterized by extremely favorable illumination and atmospheric conditions. We performed three independent calibrations using only 5 frames out of 450 recorded over a 15 min period (calibration files V18_a-c).
In Figure 4, we show the evolution of the Lab threshold estimated using these three calibrations on the 450 frames of V18. We can note that the application of the different calibration files produces reasonably similar decreasing trends for the Lab threshold. Then, we used the same set of calibration files to analyze V18 with the pixel-to-height conversion matrix presented in Figure 2b (see Section 2.2.5). We highlight that, due to the optimal illumination and atmospheric conditions of V18, an even smaller number of calibration frames seems enough to produce reproducible numerical results and similar to manual estimates, as observed in Figure S1a  In Figure 4, we show the evolution of the Lab threshold estimated using these three calibrations on the 450 frames of V18. We can note that the application of the different calibration files produces reasonably similar decreasing trends for the Lab threshold. Then, we used the same set of calibration files to analyze V18 with the pixel-to-height conversion matrix presented in Figure 2b (see Section 2.2.5). We highlight that, due to the optimal illumination and atmospheric conditions of V18, an even smaller number of calibration frames seems enough to produce reproducible numerical results and similar to manual estimates, as observed in Figure S1a   . Threshold values for the Lab mask as a function of frame number of the video V18, computed using three different calibrations based on five frames extracted from V18 (see legends). In the different panels, we present the application of different fit strategies to define the threshold for the Lab mask. Results derived from the application of clustering and polynomial fit are presented in panel a, results in panel b are associated with a criterion of nearest value, and panel c presents, for each frame, the more conservative choice between panels a and b (i.e., the minimum value). Note that PHA (see Table 2) can generate this figure automatically.
The temporal evolutions of column height, calculated using the different calibration files, are strongly consistent between them (Figure 5a), showing percentage differences of the order of 0.2% during the video (maximum value of 1.0%). Note that similar average percentage differences are computed when numerical results are compared with manual estimates ( Figure S1a). The results indicate a continuous increase in column height from ~6700 m up to ~8200 during these 15 min, with oscillations with a period of about 5 min (see Supplementary Videos S1-S3 in the Supplementary Material). The regularity in the evolution of plume height highlights the precision of this program under favorable atmospheric and illumination conditions. We dispose of an additional set of 61 frames for a longer period (1 h) of the same eruption (V18b). Since a significant change in the eruption and illumination conditions occurs, the previous calibrations are not able to capture the characteristics of the complete video (note that it was recorded at 8 a.m.). In fact, these calibrations only work for the first ~35 frames of the video. Instead, we created three new, independent calibrations based on 5 frames out of 61 recorded over a period of 1 h (calibration files V18b_a-c). Results are consistent between them and show an increasing trend of plume height from ~6700 up to more than the measurement limit (Figure 5b and Supplementary Videos S4-S6). In this case, our calibrations produce average percentage differences of the order of 0.3% (0.5% before reaching the measurement limit), with a maximum value of 4.1%. . Threshold values for the Lab mask as a function of frame number of the video V18, computed using three different calibrations based on five frames extracted from V18 (see legends). In the different panels, we present the application of different fit strategies to define the threshold for the Lab mask. Results derived from the application of clustering and polynomial fit are presented in panel a, results in panel b are associated with a criterion of nearest value, and panel c presents, for each frame, the more conservative choice between panels a and b (i.e., the minimum value). Note that PHA (see Table 2) can generate this figure automatically.
The temporal evolutions of column height, calculated using the different calibration files, are strongly consistent between them (Figure 5a), showing percentage differences of the order of 0.2% during the video (maximum value of 1.0%). Note that similar average percentage differences are computed when numerical results are compared with manual estimates ( Figure S1a). The results indicate a continuous increase in column height from 6700 m up to~8200 during these 15 min, with oscillations with a period of about 5 min (see Supplementary Videos S1-S3 in the Supplementary Material). The regularity in the evolution of plume height highlights the precision of this program under favorable atmospheric and illumination conditions. We dispose of an additional set of 61 frames for a longer period (1 h) of the same eruption (V18b). Since a significant change in the eruption and illumination conditions occurs, the previous calibrations are not able to capture the characteristics of the complete video (note that it was recorded at 8 a.m.). In fact, these calibrations only work for the first 35 frames of the video. Instead, we created three new, independent calibrations based on 5 frames out of 61 recorded over a period of 1 h (calibration files V18b_a-c). Results are consistent between them and show an increasing trend of plume height from~6700 up to more than the measurement limit (Figure 5b and Supplementary Videos S4-S6). In this case, our calibrations produce average percentage differences of the order of 0.3% (0.5% before reaching the measurement limit), with a maximum value of 4.1%.
Remote Sens. 2023, 15, x FOR PEER REVIEW 12 of 22 Figure 5. Plume height as a function of time for some reference videos (see titles) of Mt. Etna (see Table 1) using different files of internal calibration (see legends and Section 3.1). The measurement limit in panel b is 9326 m a.s.l.

V21 (3 April 2013)
V21 is a video characterized by favorable illumination conditions and partially favorable atmospheric conditions due to the presence of several clouds that often cover part of the eruption plume. In this case as well, we performed three independent calibrations with only 5 frames out of 450 (calibration files V21_a-c), and then they were used to analyze the same video. The pixel-to-height conversion matrix is presented in Figure 2c (see Section 2.2.5). The evolution of plume height calculated using the different calibration files is similar (Figure 5c), showing oscillations between ~5800 and ~6800 m a.s.l. The only significant differences between numerical results are observed in a few, specific frames near the end of the video (differences of the order of 1000 m or less), when small clouds disturb the plume identification for two of the calibrations (see Supplementary Videos S7-S9). On average, the percentage differences of plume height associated with these calibrations are of the order of 1.3%. When compared with manual estimates, the average percentage differences of plume height are instead of the order of 3.0% ( Figure S1b). Interestingly, Figure  S1b shows a regular decrease in the average percentage difference with respect to manual estimates when the number of frames used for the Lab calibration increases, suggesting that a calibration based on >8 frames would produce a significantly better performance.
An additional video with 61 frames over a longer period of the same eruption was analyzed using three different calibrations based on only 5 frames (calibration files V21b_a-c). A general trend of plume height increase can be observed (from ~5300 to ~6900 m a.s.l.; Figure 5d). The main discrepancies between the curves are product of the interference of small clouds and the plume (instead, large clouds are typically recognized and discarded; see Supplementary Videos S10-S12). The percentage differences between the results obtained with these independent 5-frames-based calibrations are of the order of 2.8% (average value), with a maximum value of 11.5%. Interestingly, the calibrations constructed using a period of 15 min (i.e., calibration files V21_a-c) are able to capture well  Table 1) using different files of internal calibration (see legends and Section 3.1). The measurement limit in panel b is 9326 m a.s.l.

V21 (3 April 2013)
V21 is a video characterized by favorable illumination conditions and partially favorable atmospheric conditions due to the presence of several clouds that often cover part of the eruption plume. In this case as well, we performed three independent calibrations with only 5 frames out of 450 (calibration files V21_a-c), and then they were used to analyze the same video. The pixel-to-height conversion matrix is presented in Figure 2c (see Section 2.2.5). The evolution of plume height calculated using the different calibration files is similar (Figure 5c), showing oscillations between~5800 and~6800 m a.s.l. The only significant differences between numerical results are observed in a few, specific frames near the end of the video (differences of the order of 1000 m or less), when small clouds disturb the plume identification for two of the calibrations (see Supplementary Videos S7-S9). On average, the percentage differences of plume height associated with these calibrations are of the order of 1.3%. When compared with manual estimates, the average percentage differences of plume height are instead of the order of 3.0% ( Figure S1b). Interestingly, Figure  S1b shows a regular decrease in the average percentage difference with respect to manual estimates when the number of frames used for the Lab calibration increases, suggesting that a calibration based on >8 frames would produce a significantly better performance.
An additional video with 61 frames over a longer period of the same eruption was analyzed using three different calibrations based on only 5 frames (calibration files V21b_a-c). A general trend of plume height increase can be observed (from~5300 to~6900 m a.s.l.; Figure 5d). The main discrepancies between the curves are product of the interference of small clouds and the plume (instead, large clouds are typically recognized and discarded; see Supplementary Videos S10-S12). The percentage differences between the results obtained with these independent 5-frames-based calibrations are of the order of 2.8% (average value), with a maximum value of 11.5%. Interestingly, the calibrations constructed using a period of 15 min (i.e., calibration files V21_a-c) are able to capture well the characteristics of the complete video V21b (Supplementary Figure S2), reflecting that the illumination conditions are nearly constant during the 1 h video V21b (note that it was recorded at 1 p.m.).

V08 (20 August 2011)
V08 is also included in the sequence of 25 low-explosivity events observed at Etna between January 2011 and April 2012 [50]. This video presents partially unfavorable illumination and atmospheric conditions, with permanent presence of a strata of lowaltitude clouds. We constructed three independent calibrations with 5 frames each (out of 900 frames), which were used to analyze the same video (calibration files V08_a-c). On the other hand, the matrix of pixel-to-height conversion of this video is presented in Figure 2a (see Section 2.2.5). The results produced by PHA, highly consistent between them, are characterized by a continuous increase in plume height up to reach the image top and thus exceed the measurement limit (Figure 6a). The increase in plume height occurs rapidly, at a rate of~600 m/min, from~4200 to~9800 m a.s.l. Some of the calibrations differ in specific frames during the waxing phase due to the diffuse limits of the plume in this period (see Supplementary Videos S13-S15), while they capture well the general tendency of plume height. In this case, on average, the percentage differences of plume height are of the order of 0.3% (1.1% before reaching the measurement limit), with a maximum value of 11.5%. With respect to manual estimates, the average percentage differences are of the order of 3.0%, and we also observe a regular decrease in the average percentage difference with respect to manual estimates when the number of frames used for the Lab calibration increases ( Figure S1c). Results suggest that, for this video, an internal calibration based on 5 frames is enough to produce reliable data of column height. V08 is also included in the sequence of 25 low-explosivity events observed at Etna between January 2011 and April 2012 [50]. This video presents partially unfavorable illumination and atmospheric conditions, with permanent presence of a strata of low-altitude clouds. We constructed three independent calibrations with 5 frames each (out of 900 frames), which were used to analyze the same video (calibration files V08_a-c). On the other hand, the matrix of pixel-to-height conversion of this video is presented in Figure  2a (see Section 2.2.5). The results produced by PHA, highly consistent between them, are characterized by a continuous increase in plume height up to reach the image top and thus exceed the measurement limit (Figure 6a). The increase in plume height occurs rapidly, at a rate of ~600 m/min, from ~4200 to ~9800 m a.s.l. Some of the calibrations differ in specific frames during the waxing phase due to the diffuse limits of the plume in this period (see Supplementary Videos S13-S15), while they capture well the general tendency of plume height. In this case, on average, the percentage differences of plume height are of the order of 0.3% (1.1% before reaching the measurement limit), with a maximum value of 11.5%. With respect to manual estimates, the average percentage differences are of the order of 3.0%, and we also observe a regular decrease in the average percentage difference with respect to manual estimates when the number of frames used for the Lab calibration increases ( Figure S1c). Results suggest that, for this video, an internal calibration based on 5 frames is enough to produce reliable data of column height.  Table 1) using different files of internal calibration (see legends and Section 3.1). The measurement limit in panels a and b is 9774 m a.s.l. and that of panel c is 9517 m a.s.l.
The same set of calibrations was able to analyze the frames of the video V08b (Table  1), which comprises a longer period of the same eruption (Figure 6b). Results associated Figure 6. Plume height as a function of time for some reference videos (see titles) of Mt. Etna (see Table 1) using different files of internal calibration (see legends and Section 3.1). The measurement limit in panels a and b is 9774 m a.s.l. and that of panel c is 9517 m a.s.l.
The same set of calibrations was able to analyze the frames of the video V08b (Table 1), which comprises a longer period of the same eruption (Figure 6b). Results associated with the different calibration files (see Supplementary Videos S16-S18) are strongly consistent and capture well the characteristics of the waxing phase described in the previous paragraph.

V30 (12 March 2021)
V30 is a~6-h-lasting video with strong changes in the illumination and eruption conditions, comprising well-defined waxing and waning phases. Consequently, a larger number of calibration frames must be considered in order to create a functional calibration file. We constructed three independent 10-frames-based calibrations (out of 361 frames), which were used to analyze the same video (Figure 6c and Supplementary Videos S19-S21). Our results, obtained by adopting the pixel-to-height conversion matrix presented in Figure 2d, show a waxing phase from a plume height of~5200 m a.s.l. to a value beyond the measurement limit, with an average rate of plume height increase in the order of 40 m/min, much smaller than that observed in V08. Note, however, that the waxing phase can be divided in two steps with different variation rates of plume height (Figure 6c). The waning phase, from beyond the measurement limit up to a plume height of~5500 m a.s.l., occurred at a rate of~120 m/min and, after that, plume height decreases slowly up to~4700 m a.s.l. The results obtained using the different calibrations are strongly similar, with average differences of the order of 2.0%, while average percentage differences of the order of 7.0% are obtained when numerical results are compared with manual estimates ( Figure S1d). Note that at least 10 calibration frames are needed to obtain average percentage differences (with respect to manual measurements) below 10% ( Figure S1d).

Operational Calibration
We have shown that PHA permits to analyze videos of volcanic eruptions by means of a fast calibration process. Although this may be useful for research and operative purposes by itself, an autonomous operative use of PHA requires that the program can recognize the characteristics of new images and estimate proper calibration parameters. In other words, in order to advance toward the implementation of a real-time column height estimation procedure using visible cameras on erupting volcanoes, we need a calibration file that considers a large set of eruption, illumination, and atmospheric conditions. To do this, PHA includes a set of functions that allows to create, merge, and improve calibrations by adding new data, permitting to refine the performance of the program as new data becomes available for the calibration as well (see Table 2).
To show the capability of the program to deal with a large calibration file and to recognize different conditions simultaneously, we merged a set of calibration files constructed using different videos (Table 1). Our dataset includes: (a) 5 frames for short-lasting videos with a capture period of 1 min (V08b, V18b, and V21b), (b) 10 frames for long-lasting videos with a capture period of 1 min (V30), (c) 40 frames for long-lasting videos with a capture period of 2 s (V26 and V29), and (d) 10 frames for short-lasting videos with a capture period of 2 s (all the other videos), totaling 375 frames. The application of a common calibration file for a large set of videos has shown that: (a) For eruptions with favorable atmospheric conditions and when the outline of the plume is well defined (V02, V10, V12, V15, V18, V18b, V23, V28, V29, and V30), PHA is able to trace accurately plume height and in some cases small-scale oscillations of this parameter can be identified as well (Figure 7). Comparisons with manual estimates of plume height are presented in Figure 7, where we can observe remarkably consistent trends. (b) For eruptions with unfavorable atmospheric conditions (e.g., small clouds interfering the visual field; V06, V08, V08b, V14, V21, and V21b), the program is able to recognize well the range of values of plume height and general tendencies, but small-scale oscillations are not traced and occasional mistakes in punctual frames are observed. However, we stress that they are typically below the intrinsic uncertainty of plume height estimations based on visible cameras [4], as observed in Figure 8, where we present comparisons with manual estimates. (c) Finally, for eruptions with plumes characterized by diffuse outlines (V01, V03, V04, V05, V09, V11 and V26), PHA is able to trace well the range of values of plume height and recognizes large-scale tendencies. However, the results present a typically oscillating behavior around the manual estimates ( Figure 9) and occasional mistakes in specific frames are observed as well. height estimations based on visible cameras [4], as observed in Figure 8, where we present comparisons with manual estimates. (c) Finally, for eruptions with plumes characterized by diffuse outlines (V01, V03, V04, V05, V09, V11 and V26), PHA is able to trace well the range of values of plume height and recognizes large-scale tendencies. However, the results present a typically oscillating behavior around the manual estimates ( Figure 9) and occasional mistakes in specific frames are observed as well.  Table 1), using a common calibration file (see Section 3.2). These videos are characterized by favorable atmospheric and illumination conditions, and the plumes present a well-defined outline.  Table 1), using a common calibration file (see Section 3.2). These videos are characterized by favorable atmospheric and illumination conditions, and the plumes present a well-defined outline.   Table 1), using a common calibration file (see Section 3.2). These videos are characterized by unfavorable atmospheric conditions (e.g., presence of clouds interfering with the visual field), and the plumes present a well-defined outline.
In general terms, results are consistent with manual measurements, showing that PHA, when an enough large number of frames are present in the calibration file, is able to deal with different eruption, atmospheric, and illumination conditions. The average percentage difference between manual and automatic measurements of column height is 2.70%, with a median of 0.59% and 90th and 95th percentiles of 8.55% and 13.73%, respectively. Regarding the absolute difference between the two estimates of column height, the average value is 166.8 m, with a median of 34.5 m and with 90th and 95th percentiles of 548.0 m and 925.1 m, respectively.  Table 1), using a common calibration file (see Section 3.2). These videos are characterized by unfavorable atmospheric conditions (e.g., presence of clouds interfering with the visual field), and the plumes present a well-defined outline.
In general terms, results are consistent with manual measurements, showing that PHA, when an enough large number of frames are present in the calibration file, is able to deal with different eruption, atmospheric, and illumination conditions. The average percentage difference between manual and automatic measurements of column height is 2.70%, with a median of 0.59% and 90th and 95th percentiles of 8.55% and 13.73%, respectively. Regarding the absolute difference between the two estimates of column height, the average value is 166.  Table 1), using a common calibration file (see Section 3.2). These videos are characterized by plumes with diffuse outlines.

The Program PHA
In this paper, we have presented the program PHA, which is a novel, open-source MATLAB tool designed to analyze images from visible cameras of volcanic plumes, with the aim of automatically recognizing the plume and estimate its maximum height as a function of time. Due to the intrinsic variability of the images captured during volcanic eruptions, which is a consequence of changes in eruption, illumination, and atmospheric  Table 1), using a common calibration file (see Section 3.2). These videos are characterized by plumes with diffuse outlines.

The Program PHA
In this paper, we have presented the program PHA, which is a novel, open-source MATLAB tool designed to analyze images from visible cameras of volcanic plumes, with the aim of automatically recognizing the plume and estimate its maximum height as a function of time. Due to the intrinsic variability of the images captured during volcanic eruptions, which is a consequence of changes in eruption, illumination, and atmospheric conditions, PHA was conceived as a self-customizing tool. This means that, before operational use, an iterative calibration procedure must be performed, based on the analysis of previous eruptions of the volcanic system of interest, possibly occurred under different environmental and volcanic conditions. The algorithm created to identify the volcanic plume largely relies on the analysis of the third channel of the images in Lab scale, which allows recognizing and discarding pixels associated with the sky, and the application of a series of procedures aimed at discarding clouds and determining high-color gradient zones in the plume.
By means of a series of illustrative applications of PHA on some events at Mt. Etna, Italy, we showed that a bounded number of frames can be used to calibrate the model and create a functional tool that is able to process data from past eruptions and trace the plume height automatically, with reproducible and accurate results when compared with manual measurements (see Section 3.1). However, to create a program that is able to recognize in real-time the characteristics of new images and estimate proper calibration parameters, a calibration with a large set of images with different characteristics in terms of illumination, atmospheric conditions, and eruption parameters, is needed. PHA includes a large set of interactive functionalities to facilitate the construction of a truly operative tool in the context of volcano monitoring, with functions to create, improve, and merge the data of different calibration files. In this paper, these functions were used to create a calibration file based on 375 images captured by the ECV visible camera of Etna. This calibration file has been shown to be useful to analyze videos of 23 events of Mt. Etna. Results are remarkably consistent with manual estimations when illumination and atmospheric conditions are favorable, while some occasional mistakes are still present when illumination and atmospheric conditions are not optimal. However, even in this case, the program permits to recognize large-scale, time-dependent tendencies of the eruptions, with differences between PHA and manual estimates typically below the intrinsic uncertainty of these measurements [4]. This first attempt shows that PHA is potentially useful to construct a tool that is able to analyze automatically visible camera images of Mt. Etna, with important applications for monitoring purposes, and may represent a significant approximation toward a real-time analysis of column height using visible cameras on erupting volcanoes.

Limitations, Strengths, and Future Advances of PHA
This tool presents the intrinsic limitations of images of visible cameras: reliable results can be obtained exclusively during daylight hours and with a small to moderate presence of clouds. They also present a bounded visual field, which translated into the presence of measurement limits as those observed in the ECV visible camera of Mt. Etna (see . On the other hand, due to the frequent presence of clouds in the summit of stratovolcanoes, this tool cannot be directly used to recognize the onset of an explosive eruption. In fact, the current version of the code does not include an automatic procedure to detect the absence of an ash plume and, in such a case, it will only indicate that the plume height, if present, is less than the minimum measurement limit.
Additionally, a calibration step is strictly needed before its use and, even if this tool can be set to analyze a posteriori the record of visible cameras from any explosive volcanic eruption, the availability of large datasets of past eruptions is necessary to construct an operative tool in real-time for monitoring purposes. Thus, this type of application is only feasible in volcanoes with frequent and well-monitored volcanic activity such as Etna, where the need of having a rapid, reliable tool to detect and measure volcanic plumes represents an impelling necessity.
Even though PHA is still not an operative instrument at INGV-OE, the results presented in this paper are encouraging in terms of the applicability of a customizable tool to estimate plume height as a function of time for monitoring purposes. To further improve the performance of PHA, the inclusion of additional videos of past eruptions is needed, as well as more code reliability testing and analysis of videos from other volcanoes with enough eruptions recorded. Additionally, the structure of the code allows for the refine-ment of the model by including new variables that are able to characterize the calibration frames (i.e., increasing the dimensions of the vector X i ; see Section 2.2.2); thus improving the characteristics of the fit used to calculate the calibration-derived inputs. We emphasize that this code can also be adapted to analyze images from thermal cameras, which uncovers additional development opportunities for this code in the future.

Conclusive Remarks
A new open-source MATLAB tool named Plume Height Analyzer' (PHA) able to analyze images from visible cameras of volcanic plumes and automatically estimate its maximum height as a function of time is presented in this paper. Although the tool uses an iterative calibration procedure based on the analysis of previous eruptions of a given volcano and should be tested under different environmental and volcanic conditions, it may represent a first approximation toward a real-time analysis of column height using visible cameras on erupting volcanoes.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs15102595/s1. Figure S1: Average percentage difference with respect to manual measurements of column height for a set of results obtained with the program PHA, as a function of the number of frames used for the construction of the calibration files adopted (internal calibration). For each video (see titles), we created fifteen independent calibration files with different numbers of calibration frames (from 2 to 10 in panels a-c, from 4 to 20 in panel d). In each panel, we present the mean value (filled circles) and standard deviation (bars) associated with the application of three independent calibration files, created using different numbers of calibration frames (see x-axis). Figure S2: Plume height as a function of time for video V21b of Etna (see Table 1) using different files of internal calibration (see legends and Section 3.1). Video S1: Video V18 analyzed using the calibration file V18_a. Video S2: Video V18 analyzed using the calibration file V18_b. Video S3: Video V18 analyzed using the calibration file V18_c. Video S4: Video V18b analyzed using the calibration file V18b_a. Video S5: Video V18b analyzed using the calibration file V18b_b. Video S6: Video V18b analyzed using the calibration file V18b_c. Video S7: Video V21 analyzed using the calibration file V21_a. Video S8: Video V21 analyzed using the calibration file V21_b. Video S9: Video V21 analyzed using the calibration file V21_c. Video S10: Video V21b analyzed using the calibration file V21b_a. Video S11: Video V21b analyzed using the calibration file V21b_b. Video S12: Video V21b analyzed using the calibration file V21b_c. Video S13: Video V08 analyzed using the calibration file V08_a. Video S14: Video V08 analyzed using the calibration file V08_b. Video S15: Video V08 analyzed using the calibration file V08_c. Video S16: Video V08b analyzed using the calibration file V08_a. Video S17: Video V08b analyzed using the calibration file V08_b. Video S18: Video V08b analyzed using the calibration file V08_c. Video S19: Video V30 analyzed using the calibration file V30_a. Video S20: Video V30 analyzed using the calibration file V30_b. Video S21: Video V30 analyzed using the calibration file V30_c.

Data Availability Statement:
The data used in this paper belong to INGV and can be made available upon request to the authors.