Predictive Tracking of Continuous Object Boundaries Using Sparse Local Estimates

Environmental hazards (wildfires, floods, oil spills) are often modeled as “continuous objects” which evolve in space and time taking irregular shapes. Tracking their boundaries and accurately predicting their spatiotemporal spreading patterns is of paramount importance to combat their often catastrophic consequences. Wireless Sensor Networks (WSN) have been very instrumental for this purpose. However, current WSN-based methods require a prohibitively large density of deployed sensors to achieve a reasonable boundary reconstruction accuracy because they are based on the explicit identification of nodes close to the boundary. We present a novel approach that can track and predict the global boundary using only a relatively small number of distributed local front estimates. Our approach first filters and fuses the available sparse set of local front estimates (e.g. vectors of local orientation, direction and speed) and then uses the resulting information to reconstruct a smooth curve prediction of the evolving object’s boundary at a future time. Moreover, since the uncertainty of the local front parameters is modeled, it can provide a heatmap representation of the evolving object, indicating the probability for each point in the area of interest to be reached at a future time by the spreading hazard. These predictive modeling capabilities when combined enable effective decision support for crisis management. We demonstrate that different types of continuous objects can be tracked with accuracy, even when only a relatively small number of noisy local front estimates is available. Our approach is practical since it can be applied in many situations where global boundary prediction updates are important to obtain by new sparsely distributed noisy local front estimates as soon as they reach a control center while the hazard is still progressing.


I. INTRODUCTION
During the last decade, there is a growing interest in exploiting Wireless Sensor Networks (WSNs) and remote sensing technologies to track effectively diffusive hazardous phenomena (e.g. wildfires, oil spills, chemical leaks etc.) [1]- [17]. Predicting accurately the evolution of such ''continuous objects'' spreading over a large geographical area is of paramount importance, since it allows authorities to optimize their response (hazard suppression, evacuation plans etc.) The associate editor coordinating the review of this manuscript and approving it for publication was Gongbo Zhou .
to save lives and limit the scope of a hazard's catastrophic effects.
Several WSN-based methods have been proposed for tracking the boundaries of a continuous object [2]- [15]. Their key strategy is to identify over time the sensor nodes closest to the object's evolving front-line (called boundary nodes). More specifically, they either use techniques to reconstruct the continuous object's boundary (e.g. by determining the convex-hull polygon of the boundary nodes locations [15]) or rely on human operators to identify the boundary's shape using the boundary nodes locations. As a results, existing methods demand unrealistic sensor densities (i.e. thousands of sensors per km 2 ), a requirement that renders them VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ impractical for real-world environmental monitoring applications. In addition, they do not track the local evolution characteristics of different parts of the continuous object and thus cannot provide decision support using predictive modeling. Dynamic Data Driven Application Systems (DDDAS) [16]- [19] try to address some of these severe limitations. Their main objective is to continuously improve the quality of predictions produced by hazard-specific models (mathematical, semi-empirical, etc.) by recalibrating them periodically using available sensing data streams from the field. However, the development of such systems is a very difficult task since in most cases it is almost impossible to exploit raw sensor data to recalibrate the parameters of hazard-specific simulation models.
Computing the convex-hull, is a very well studied problem [20], [21] and it has been extensively used in many application domains (e.g. image analysis, epidemiology, etc.) to determine the boundary of a given set of points (e.g. pixels, patient locations, etc.). However, there are many cases where the convex-hull cannot accurately represent the shape of the area occupied by the set of points (e.g. when the points form a non-convex shape). To overcome this difficulty, many convex-hull generalization algorithms have been proposed [22]- [25]. Although, they can determine a ''tighter'' non-convex area occupied by the given set of points, their boundary reconstruction accuracy greatly depends on the selection of specific parameters which are very difficult or even impossible to estimate [25], [26]. Moreover, all these algorithms assume a dense set of points uniformly distributed along the object's boundary area, a fact that renders them impractical for scenarios where only a small number of irregularly distributed points is available. Finally, these algorithms cannot take into account boundary reconstruction constraints (like C1 and C2 as defined in Section IV-A) which are necessary to capture correctly the evolution dynamics of a continuous object.
The boundary (global front) of a continuous object can be approximated by a set of line segments (local fronts) and each segment's evolution can be represented by a model with a small number of parameters (e.g. orientation, direction, and speed), as shown in [27]- [29]. In our prior work we have introduced a probabilistic framework and associated distributed algorithms that can estimate the local evolution characteristics of a propagating hazard's front using WSNs of realistic density. As the hazard progresses, we have shown how the deployed sensor nodes can be dynamically re-organized into small-size ad-hoc clusters (node triplets) and compute, within their local area, estimates of the hazard's characteristics. Moreover, we have shown how the sensor nodes can update the local front model parameters using closed-form expressions (analytical solutions of a Bayesian parameter estimation problem) that are easy to implement using their commodity microprocessors [30]. Importantly, a unique feature of our WSN-based environmental monitoring methodology is that each sensor node can also estimate the uncertainty associated with the local front estimates it produces as the hazard progresses [29].
In this work (that complements [29]) we address the interesting problem of how a stream of sparse (in space and time) and inherently noisy local front estimates can be used (as they are becoming available to the decision center) to reconstruct accurately the global boundary of a continuous object. Our approach, which combines machine learning and computational geometry methods, offers the following novelties and unique advantages that are also very important from a practical standpoint for environmental monitoring applications: (i) It can reconstruct the boundary of an evolving continuous object using only a small number of sparse local front estimates; (ii) it is robust to local front parameter estimation errors; (iii) in contrast to other methods [2], [6], [8]- [14] it determines the global boundary explicitly and does not require a human in the loop to recognise it implicitly based on the locations of the identified boundary nodes; (iv) it supports a probabilistic notion of a continuous object that is now viewed as a field, indicating for each point of the considered area the probability to be reached by the hazard at any given time. The methods presented here are quite general as different sensor modalities (WSN's, smart cameras etc.) can be abstractly viewed as nodes of a heterogeneous monitoring system generating inherently noisy estimates of local scope that a data center has to reconcile and fuse continuously in order to predict with accuracy an evolving hazard's irregular global boundary, even in areas with a limited number of available sensors.
The rest of the paper is organized as follows: In Section II we introduce the notation and preliminaries needed to understand the paper and motivate the key elements of the proposed global boundary tracking and prediction approach. In Section III we explain how the available local front estimates can be filtered and fused in order to generate consistent information for feeding the boundary reconstruction algorithm, presented in Section IV. In Section V we present and discuss simulation results for different types of hazards. Finally, we summarize our findings and outline work in progress in Section VI. Illustrative video animations [31], [32] (provided also as supplementary material) are used to demonstrate the online reconstruction algorithm ''in action''. We urge the readers to view these videos as they study Section V.

A. LOCAL FRONT MODELS
The boundary of a continuous object can be approximated by a piecewise linear curve. We represent each segment of this curve (to be called a local front) by a simple model that captures locally (i.e. within a circular area C i of radius R i which depends on the sensor node type, see Fig. 1) its local evolution characteristics (e.g. direction, orientation and speed). So in our abstract model, a local front estimate m i is represented by the following parameters: • Location (P i = (x i , y i )): The physical coordinates of the line segment's middle point. Local front segments are assumed to be at close proximity to the surface of the earth (the z offset is ignored).
• Length (l i ): The length of the line segment that approximates locally (within a circle of radius R i = l i 2 ) the continuous object's boundary. The larger the radius R i the larger the scope of the local model. • Orientation parameters (φ i , σ i ): The angle i formed by the local segment and the horizontal axis (x-axis) is considered to be a random variable following a Normal distribution i ∼N (φ i , σ 2 i ). We will denote withφ i a realization (random sample) of i .
• Speed parameters (u i , s i ): The speed U i of the local front's line segment is considered to be a random variable that follows a Normal distribution U i ∼N (u i , s 2 i ). We will denote withû i a realization (random sample) of U i . • Local front model parameters estimation time (t i ), measured with respect to a global time reference. Using the center coordinates P i = (x i , y i ), the length l i , and the orientationφ i of a local front estimate m i , we can determine the coordinates of its end points P Ez We will use the notation m i (t j ) to denote the evolution of the local front model estimate m i generated at time t i to a future time instance t j (t j > t i ). In the text we may refer to m i (t j ) as the future ''projection'' of estimate m i (t i ). The midpoint, the end-points and the circular area of the space-time evolved local front m i (t j ) will be denoted as and C i (t j ) respectively (see Fig. 2). In Appendix B we provide the derivation of the closed form algebraic expressions used to calculate the location coordinates of a space-time evolved local front. Finally, we will use A i (t j ) to denote the ''area covered'' or ''swept'' by m i during its space-time evolution from time t i to time t j (see green shaded area in Fig. 2).

B. KEY ELEMENTS OF THE PROPOSED APPROACH
We assume that a distributed monitoring system (e.g. based on remote sensing, WSNs etc. [29]) provides to a decision center a stream of sparse in space and time local front estimates sampling the evolution behavior of a continuous object (hazard) at different locations and/or time instances, denoted by the black line segments in Fig. 3a. As soon as a sufficient number of local front estimates has become available, the proposed algorithm filters and fuses their information to determine the subset of them (see black segments in Fig. 3b) that can best (most consistently) represent the continuous object's evolution and will thus be used to reconstruct its boundary at a desirable time point (e.g. time instance t 12 in Fig. 3). Next, using the local fronts' evolution parameters, the algorithm predicts their projected locations at the desired reconstruction time (see Fig. 3c) and determines a polygon approximating the continuous object's boundary (see black dashed polygon in Fig. 3d). Subsequently, using uniform cubic B-splines the algorithm determines a ''smooth'' curve which is the reconstruction of the hazard's boundary curve (see red curve in Fig. 3d). Finally, based on the uncertainties associated with the local front estimates, the algorithm provides a probabilistic view of the continuous object, computing for each location and time instance the probability to be covered by the object (see Section V-D).

III. FILTERING AND FUSING LOCAL ESTIMATES
In this section we elaborate on the procedure used to determine the best subset of local front estimates providing a consistent view of the phenomenon and can be used to reconstruct a faithful representation of the continuous object's boundary.
Let's assume that we want to predict where the hazard's global boundary will be located at time t b . From the available local front estimates we select those with estimation times t i ≤ t b , and form a set M containing their parameters (see Section II-A). Our objective is to combine (i.e. filter and fuse) the spatio-temporal information of available local front estimates in M in order to determine the best subset of local front estimates that provide us with a consistent view of the hazard's front at the boundary reconstruction time t b . If M contains at least N (specified by the user) of local front estimates (i.e. (|M| ≥ N )), we initiate the following procedure. At this point we have to mention that in order to reconstruct a smooth boundary (see Section IV-B) we need at least 2 local front estimates (N ≥ 2). However, for higher values of N , we expect to have delayed (the algorithm waits until N local fronts become available) but more accurate reconstructions of the evolving boundary.
Step 1: We assume that two local front estimates m i and m j , generated at t i and t j respectively where t i ≤ t j , describe the same reality (same part of the evolving boundary), if the Euclidean distance of their mid-points (||P i , P j ||) is smaller or equal to the larger radius of their circular areas of scope (see Fig. 4).
For each pair of local front estimates in M (e.g. (m i , m j )) that satisfies condition (1), we keep only the estimate generated more recently (m j in Fig. 4, since we assume w.l.o.g that t i < t j ). This decision is justified based on the presumption that the more recent local front estimate captures better the current time varying evolution characteristics of the continuous object. The time complexity of Step 1 is bounded by O(|M| 2 ). At the end of this filtering step, we use the local fronts remaining in M to form a new set M f ⊆ M containing their parameters information. Next, we check if M f contains at least N local front estimates, and if so proceed to the next step.
Step 2: Reconstructing the boundary of a continuous object at time t b , based on the space-time evolved local front estimates in M f without considering possible intersections of their evolution paths usually leads to boundary shapes that deviate considerably from reality. Below we propose a method that first identifies and then handles ''events'' that may lead to such distortions.
Events Identification Procedure: For each pair of local front estimates in M f (e.g. (m i , m j ), where w.l.o.g. t i ≤ t j ), we check if any of the following mutually exclusive events may occur: Type-1 event -Evolved center P i approaches center P j : After projecting the local front estimate m i to the time t j of the local front estimate m j we check if {{P i (t j ) ∈ C j and P j ∈ A i (t j )} or {P j ∈ C i (t j ) and P j ∈ A i (t j )}}. If the above condition holds (see Fig. 5a), an event of type-1 is registered for time t j .
Type-2 event -Evolved m i area covers P j : After projecting the local front estimate m i to time t j we check if {P j ∈ A i (t j )}. If the above condition holds (see Fig. 5b), an event of type-2 is registered for time t j .
Type-3 event -Two evolving local front models block each other: If during the joint space-time evolution of two local fronts (m i , m j ) beyond time t j there exist a time instance t d (t i ≤ t j < t d ≤ t b ) such that the evolution path of one event blocks the evolution path of the other i.e.
, an event of type-3 is for time t d (see Fig. 5c).
Type-4 event -Two evolving local front models can be fused: If during the space-time evolution of the two local fronts (m i , m j ), there is a time instance t f such that t i ≤ t j < t f ≤ t b ) where the Euclidean distance of their mid-points becomes equal to the larger radius of their circular areas and none of them falls within the already covered area of the other, an event of type-4 is registered for time t f (see Fig. 5d). At this point we have to mention that all the registered event time instances (e.g. t j , t d and t f ) correspond at the starting time of the corresponding events. Moreover, the events identification procedure can be parallelized to increase the computational efficiency. At the end of the events identification procedure the algorithm finds the registered event with the earliest time and passes the involved pair of model estimates (e.g. {m i , m j }) to the events handler (see below) that implements the following procedure for each event type.
Events Handling Procedure: Handling type-1 events: In this case the two local front estimates m i (t j ) and m j describe at time t j the evolution behavior of the same part of the continuous object boundary (see Fig. 5a). Using the presumption that the more recently estimated local front (m j ) describes better the current evolution behavior of the continuous object, we keep in M f only the most recently updated local front (m j ).
Handling type-2 events: In this case the local estimate m j appears inside the area which has been already covered by m i (t j ) (see green shaded area in Fig. 5b). This implies that the evolved local front m i (t j ), failed to describe the evolution behavior of the continuous object boundary, since the more recent local front m j at time t j , should be located on the continuous object's boundary and not inside the area which has already been covered by the continuous object. Based on this reasoning we keep in M f only the local front estimate m j .
Handling type-3 events: In this case the evolution path of a local front estimate m j (t d ), ''blocks'' the evolution path of m i (t d ) or vice versa (see Fig. 5c). To handle this event type we keep in M f only the most recent local front estimate (m j (t d )).
Handling type-4 events: In this case two local front estimates (m i (t f ) and m j (t f ) in Fig. 5d) ''meet'' at some point during their space-time evolution without actually blocking each other. We consider this as an indication that two parts of the global boundary are ''merging'' at t f . To handle this event type properly we apply a technique that fuses the information of the two local front estimates {m i (t f ), m j (t f )} to produce a new local front estimate {m f (t f )} that better describes at time t f the local evolution behavior of the continuous object's boundary (see darkest local front m f in Fig. 5d). To estimate the parameters of the new local front model we apply the following equations: The equations in (2)  At the end of Step 2 we check if |M f | ≥ N : If the condition holds, the procedure finds the space-time evolutions of the local fronts in M f at the desired boundary reconstruction time t b , forms a new set M b that contains their information, and proceeds to the third and final step. The derivation of the equations used to calculate the location of projected local fronts is provided in Appendix B.
Step 3: This final correction step adjusts the length parameter of some projected local front estimates in M b in order to avoid boundary deformations during the final curve reconstruction (see Fig. 6a).
The procedure checks if there is (are) any pair(s) of local front estimates in M b for which their local circular areas overlap (e.g. see Fig. 6a). Next, it finds the local front with the earliest estimation time and reduces their lengths by an amount that is inversely proportional to the confidence ({w i , w j }) we have about their corresponding estimated evolution parameters, such that their circular areas do not overlap but osculate externally (see equation (3) and Fig. 6b). where {R i , R j } are the new radii of the local fronts' circular areas calculated using equations (4) and (5) below.
Equations (4) reduce the radii {R i , R j } by an amount that is inversely proportional to the corresponding local front weights {w i , w j }. Equation (5) is the solution for c of the system of equations (3) and (4).
Using the new radii R k where k = {i, j}, we calculate the end-point coordinates P

IV. GLOBAL BOUNDARY RECONSTRUCTION
We present now the algorithm which reconstructs the continuous object's boundary using the information of the selected local front estimates in M b . The algorithm has to two main phases: Phase 1: It forms a polygonal approximation of the boundary.
Phase 2: It constructs a curve that smoothly approximates the polygonal boundary.
To better explain the two phases and without loss of generality, we will use a running example: Let's assume that M b contains after filtering and fusing 14 local front estimates (see Fig. 7). In Fig. 7 the black arrows indicate the evolution directions of the local front segments and the black dots their mid-and end-points.

A. PHASE 1: POLYGONAL APPROXIMATION
This phase has two Steps: Step 1 determines the order in which the local front segments in M b have to be ''stitched'' in order to form the polygonal approximation of the continuous object's boundary. Using the local fronts' connection order determined in Step 1 and their evolution direction parameters, Step 2 constructs the polygonal approximation.
For the approximation to be realistic, the formed polygon should satisfy the following two conditions: • C1: Be simple (no intersections between its edges). • C2: Local front evolution direction vectors should point towards outside the polygon's area. Before initiating Step 1, we form the following sets: C: Contains the mid-points of all local front estimates in M b .
V: Contains the vertices of the convex hull of points in set C (in our example, V = {P 3 , P 4 , P 7 , P 8 , P 11 , P 13 , P 14 }, see Fig. 7a) in a connection order that forms the convex hull polygon. To find the convex hull polygon we use the Graham scan algorithm. The time complexity of this algorithm is O(|C| · log(|C|)). It holds that V ⊆ C. Q = C ∩ V c : Contains the points in C that do not belong to the convex hull polygon. In Fig. 7a, Q = {P 1 , P 2 , P 5 , P 6 , P 9 , P 10 , P 12 }. Note that V and Q are partition subsets of C.
The time complexity for the formation of the sets C and Q is bounded by O(|M b |).

STEP 1
At this step the algorithm appropriately places the points of Q in V. At the end of this step, the ordering of points in V, indicates the order in which the local fronts in M b have to be connected to form the polygon that approximates the continuous object's boundary.
We check the set Q : • If {Q = ∅} : We continue to Step 2, since all local front mid-points are already included in V.
• If {Q = ∅} (the general case): We repeat the following procedure until Q becomes empty: Using the polygon determined by V (e.g. see convex polygon in Fig. 7a) we partition Q into the following subsets: Q on : Contains the points of Q located on the polygon's edges.
Q out : Contains the points of Q located outside the polygon's area.
Q in : Contains the points of Q located inside the polygon's area.
Then we check: If {Q on = ∅}: We select a point from Q on (e.g. P 2 in Fig. 7c), determine the polygon's edge on which it lies (e.g. P 1 P 3 in Fig. 7c) and apply an ''edge split'' operation (i.e. the point becomes a vertex of the polygon). An edge split occurs when a point from Q (P 2 in Fig. 7c) becomes a polygon vertex (and is inserted into V, between points P 1 , P 3 which define the edge on which it lies). After an edge split we remove the ''new'' polygon vertex (P 2 ) from Q. The aforementioned procedure is repeated until Q on becomes empty.
If {Q on = ∅ and Q out = ∅}: We use the points in Q out and form the following sets: Q I out : Contains the points of Q out that are centers of local front segments intersecting at  least one edge of the polygon (e.g. P 9 in Fig. 7d). Q T out : Contains the points of Q out that are centers of local front segments located totally outside (do not intersect) the polygon (e.g. P 5 in Fig. 7f). Q I out and Q T out are partition subsets of Q out .
After forming these subsets we check: If {Q I out = ∅}: We find the point in Q I out with the smallest distance to the polygon's edge that it intersects (e.g. P 9 in Fig. 7d). Using this point, we split the intersected edge (edge P 8 P 10 is split to P 8 P 9 and P 9 P 10 , see Fig. 7e) and update the sets V and Q. After an edge split we repeat Step 1 from the beginning.
If {Q T out = ∅}: We find the point in Q T out (e.g. P 5 in Fig. 7f) with the smallest distance to a polygon's edge and apply an edge split (e.g. edge P 4 P 6 is split into P 4 P 5 and P 5 P 6 , in Fig. 7g). After an edge split, we update the sets V and Q and repeat Step 1 from the beginning.
If {Q on = ∅ and Q out = ∅ and Q in = ∅}: We partition Q in into two subsets: Q I in : Contains the points of Q in that are centers of local front segments intersecting at least one edge of the polygon (e.g. {P 1 , P 6 , P 10 } in Fig. 7a). Q T in : Contains the points of Q in that are centers to local front segments located totally inside the polygon (e.g. {P 2 , P 5 , P 9 , P 12 } in Fig. 7a). Q I in and Q T in are partition subsets of Q in . VOLUME 8, 2020 Then we check: We find the point in Q I in (e.g. P 1 , in Fig. 7a) with the smallest distance from the polygon's edge that it intersects. Using this point we split the intersected edge (P 14 , P 3 , in Fig. 7a), and update the sets V and Q. After an edge split we repeat Step 1 from the beginning.
If {Q T in = ∅}: We derive for each point in Q T in the equation of the line that emanates from it and is perpendicular to the local front segment (e.g. line ε in Fig. 7h for point P 12 ). For each such line, we determine its intersection points with the polygon. We select the intersection point that is consistent with the local front's evolution direction (point A in Fig. 7h). Next, we find the point in Q T in with the smallest distance from the corresponding intersection point and use it to split the corresponding polygon's edge (e.g. P 11 P 13 is split into P 11 P 12 and P 12 P 13 , in Fig. 7h). After the edge split we update the sets V, and Q and repeat Step 1 from the beginning. The time complexity for the formation of the sets Q on , Q out , Q in , where |E| is the number of the polygon's edges.
Step 1 is repeated |Q| times. Thus the complexity of Step 1 is bounded by O(|Q| 2 · |E|). At the end of Step 1, the ordered set V uniquely determines a polygon (see Fig. 7i) which indicates the connection order of the local fronts.

STEP 2
This step uses the connection order and evolution direction parameters of the local fronts and first produces a refined polygonal approximation of the continuous object's boundary.
First, we sort the information in M b according to the ordering of the mid-points of local fronts in V. To sort the information in M b we used the Heapsort algorithm that has time complexity O(|M b | · log(|M b |)). Next, using the coordinates of the end-points of the first two local front segments in M b (e.g. Fig. 8a) we find the minimum distance pair of end-points (excluding pairs that belong to the same segment), and connect them (see red dashed edge P E 1 1 P E 1 2 in Fig. 8a). Then, using the non-connected end-point of the second local front segment (P E 2 2 ) and the end-points of the third local front segment (P E 1 3 , P E 2 3 ), we find the minimum distance pair and connect these points (see edge P E 2 2 P E 1 3 in Fig. 8a). Applying the aforementioned sequential procedure successively to all local front segments in M b , results to a refined polygon which has as vertices the midand the end-points of the local front segments (see Fig. 8b).
However, the aforementioned construction procedure does not guarantee that the refined polygon will satisfy conditions C1 and C2 discussed in Phase I for polygonal approximations of continuous object boundaries. To ensure that these conditions hold we do the following: For each local front segment we check: a) if the edges that connect it with its adjacent local front segments intersect with other edge(s) (non-simple polygon), and b) if its evolution direction vector points towards inside the polygon's area. If at least one of these conditions hold (e.g. in Fig. 8b the direction vectors of m 1 and m 10 point towards inside the polygon), we interchange the corresponding local front's end-point connections to its neighbors (see Fig. 8c). If after a connections interchange, either one of the aforementioned conditions still holds true, we drop the ''problematic'' local front estimate from M b and form a new polygon by connecting the end-points of its adjacent remaining local front segments. The time complexity for the formation of the refined polygon is bounded by O(|M b | 2 ).

B. PHASE 2: SMOOTH BOUNDARY RECONSTRUCTION
In this phase we use the vertices of the refined polygon in conjunction with uniform cubic B-splines to generate a curve that smoothly approximates the continuous object's boundary. In nature, the boundaries of phenomena that can be modeled as continuous objects (e.g. wildfires, oil spills, etc.) tend to form smooth shapes and not polygonal. Inspired by this observation we used the popular cubic B-splines [33]- [36] to produce a smooth approximation of the continuous object polygonal boundary. Cubic B-splines use 3rd degree polynomials to represent each curve segment, and constraint the points that join the curve segments to meet C 2 continuity (no polarity changes in slope at the transition between segment i and segment i + 1) which guaranties a smooth curve boundary representation [37].
In our case the cubic spline curve uses as first control point the end-point of a local front segment (P E 1 1 , see Fig. 8a). This selection guaranties that the formed cubic spline will be tangent to mid-points of the local front segments. The red curve in Fig. 8d corresponds to the formed B-spline curve approximating the refined polygon approximation of the continuous object. For the smooth boundary reconstruction we applied the Carl R. de Boor's algorithm which requires O(ξ 2 ) + O(ξ ) operations where ξ is the degree of the polynomial functions (in our case ξ = 3). The time complexity for reconstructing a smooth boundary using cubic B-splines is bounded by O(|M b | · ξ 2 ) [37].

V. EVALUATION RESULTS
We present next simulation results demonstrating the ability of the proposed algorithm to track and predict with accuracy the boundary of continuous objects under different conditions.

A. THE SIMULATION ENVIRONMENT
For the evaluation of the boundary tracking and prediction algorithm we developed a flexible simulator in MATLAB which allows us to generate scenarios with different: a) propagating continuous object properties (shape, speed, acceleration), b) placement strategies (number and locations) of distributed sensors in the continuous object's propagation area, and c) induced estimation errors of local fronts' evolution parameters.
Before initiating a simulation, a MATLAB procedure takes as input the continuous object's propagation properties and the location coordinates of the assumed available sensors that when reached generate the local front estimates. Using this information, the MATLAB procedure computes the local front evolution characteristics (speed, orientation and direction) when the hazard reached a sensor, thus simulating the local sensing process (see Appendix D for details), and store them in a structure containing for each local front the values of its parameters, as defined in Section II. The generated structure is then passed as input, along with the continuous object's propagation properties, to another MATLAB simulator which implements the proposed boundary reconstruction algorithm and evaluates its accuracy.

B. EVALUATION METRICS
To evaluate the accuracy of the proposed algorithm, we consider the continuous object's area as a grid of square cells, and based on the real and reconstructed boundaries we classify the cells as True Positive, False Positives or False Negatives. To assess the boundaries similarity of the true and the reconstructed continuous object we use the F1-Score which is the harmonic mean of precision and recall (see Appendix E for details).

C. EXPERIMENTAL SETUP
A notable advantage of our proposed boundary reconstruction algorithm is that it can track with accuracy continuous objects using a small number of local front estimates (sparsity). To demonstrate this novel feature we have used local front estimate densities that are considered low for environmental monitoring applications, specifically {10 −5 , 2 × 10 −5 , 3 × 10 −5 , 4 × 10 −5 } estimates m 2 , corresponding to 10, 20, 30 and 40 sensors available within a 1km 2 square area respectively. For each density value we used a large number of randomly drawn sensor placements and demonstrate how the algorithm performs under different densities as well as errors of local front evolution parameters (speed and direction). For all experiments the local front segments length was set to l = 100m (i.e. assuming sensors of sensing area with radius R = 50m). The evolution parameters (speed, orientation angle and direction) of the local fronts were set equal to the true values that the spreading hazard had when it reached the corresponding sensor that provided the measurement (perfect estimation, see assignment method details in Appendix D). However, to evaluate how the algorithm performs under local front orientation and speed estimation errors, we have repeated 6 times each experiment and tested the accuracy of the boundary construction. Specifically the data values used to model the orientation and percent speed errors were {10 o , 20 o and 30 o } degrees and {10%, 20% and 30%} deviation from the ground truth boundary's speed respectively. For all experiments the minimum number of local front estimates required for the initiation of the boundary reconstruction algorithm was set to N = 3.

D. DIFFUSIVE HAZARD WITH REGULAR SHAPE
In Experiment 1 the diffusive hazard is modeled as a circular continuous object of fixed center but increasing radius, with the center located at the middle point of a 2km × 2km square area. Considering the area's bottom left corner as being at the origin, the circle is centered at point (1000m, 1000m) and has initial radius equal to 1m. The speed at which the radius is increasing is described by a triangular function (see Fig. 9) with initial value 2.5m/min. The speed increases with constant rate (0.0265m/min 2 ) until it reaches its maximum value (5m/min). Beyond that time point the speed starts decreasing at the same rate until it returns back to its initial value. At the end of the simulation the circle encloses the 1km 2 local fronts'  placement square area centered at point (1000m, 1000m) (see red circle in Fig. 10f). Modeling propagating hazards with circular shape is justified as an approximation because Fick's second law indicates that the diffusion of a substance emanating from a single point source covers a circular area whose size is increasing at a rate indicated by the diffusion coefficient [39]. FIGURE 10. Video frame snapshots: Panels (a) through (f) show 6 snapshots of the simulation of Experiment 1. Line segments represent the available sparse local front estimates. With black color we mark the local front estimates selected to participate to the boundary reconstruction algorithm at each t b . The blue curve represents the reconstructed boundary. The full video is available at [31].
We have created a MATLAB animation [31] (provided online and as supplementary material) that demonstrates the ability of the proposed algorithm to track the boundary of the evolving continuous object. Fig. 10a shows the circular front (red circle), its center (red dot), the 20 randomly placed local front estimates (magenta line segments) and their placement area (dashed square area of 1km 2 ). We assume that a local front estimate starts participating to the boundary reconstruction algorithm when the front reaches the mid-point of the corresponding line segment (assumed measurement time at the sensor located at the mid-point). With black color, we indicate the subset of the local front segments (see Fig. 10b -10f) that are selected to contribute to the approximation of the boundary curve, based on the selection method described in Section III. Using the characteristics of the corresponding local fronts, we determine their space-time evolutions (see Appendix B) and based on their ''new'' (predicted) locations (black segments on blue curve) we derive the ''smooth'' B-spline curve (blue curve) that approximates the circular continuous object, based on the algorithm described in Section IV.
To evaluate the boundary tracking accuracy of the proposed algorithm we attempt to reconstruct the boundary of the diffusive phenomenon 20 times during its evolution at equally spaced (every 188/20 min) time points. We have to note that at each boundary reconstruction time point (t b ), we use only the local front estimates at locations reached by the diffusive phenomenon up to this time. The boxplots in Fig. 11a summarize the distribution of the boundary tracking accuracy considering 1000 random estimate placements (1000 runs) per local front estimates' density value. For the generation of each boxplot we used as sample points the F1-scores computed at all the time points a boundary reconstruction was attempted. As we observe from the boxplots the boundary reconstruction accuracy increases with the local front estimates density. Moreover the variability of the boundary reconstruction accuracy decreases as the density of local estimates increases. This implies that as the number of local front estimates increases we become more certain that the proposed algorithm will reconstruct with accuracy the diffusive object's boundary.   11b shows the average boundary reconstruction accuracy (F1-score) as a function of time for different local front estimate densities. We observe that the accuracy of the boundary tracking drops slightly for a certain time interval and then recovers again. The time extent and experienced accuracy dip increases as the density of available local estimates decreases. To explain this behavior we have to consider the following: At the first time steps the size of the continuous object is small and therefore a small number of local estimates suffices to produce an accurate representation of its boundary. As the size of the diffusive phenomenon gets bigger with time the accuracy of the boundary reconstruction drops until a sufficient number of local estimates becomes available (see also Fig. 11c) to participate in the boundary reconstruction algorithm.
In Fig. 11c we see that as the local front estimates' density increases, the mean number of the available local estimates that participate to the boundary reconstruction algorithm at each time step of the diffusive phenomenon evolution also increases. This larger number of available local estimates explains why the time interval where the accuracy drops is smaller at larger densities (see Fig. 11b). The slightly decreasing trend of the accuracy observed beyond 150 min is explained if we consider that up to that time point most local estimates have become available (see Fig. 11c) and therefore the reconstruction of the boundary beyond that point is based mostly on ''outdated'' local front estimate projections.
Another interesting observation is that beyond t = 150 min the boundary reconstruction accuracy reaches almost the same high level for all considered local front estimate densities (see Fig. 11b). This behavior can be explained if we consider that the proposed algorithm selects only a small number (the most suitable subset) of the available local front estimates (see Section III) to reconstruct the object's boundary, which was about the same for all the local front densities. This shows that the prediction ability of our algorithm is almost insensitive to the local front estimates' density, a fact that further supports our claim that it can be used to accurately track and predict the evolving boundary of a continuous object using sparse local front estimates. Fig. 11d shows the probability to have a successful boundary reconstruction event at each one of the considered time instances. We observe that as the local estimates density increases, the probability for a successful boundary reconstruction event also increases for all reconstruction time instances (dots on the curves of the Fig. 11d). This behavior is justified if we consider that the higher local estimates density implies more available local estimates at each reconstruction time instance which in turn increases the probability to have a sufficient number of local front estimates in M b (|M b | ≥ N see Section III) which is necessary condition for a boundary reconstruction to be attempted.

Boundary tracking evaluation under orientation and speed estimation errors
In this experiment we evaluate the boundary tracking accuracy of the proposed algorithm under different local front orientation and percent speed error values. For each scenario (specific density and error values) the measured evolution parameters (orientation or speed) of each local front estimate are distorted by applying the selected error values (discussed in Section V-C). The distortion is implemented by adding to, or subtracting from, each estimate (randomly selected action) the chosen error value for the corresponding parameter. Fig. 12a and 12b show for each angle and percent speed error value used, the average boundary tracking accuracy for all the local estimates density scenarios considered. For the 20 o angle error and the 20% speed error we also provide for each local estimates density value the corresponding standard deviations (error bar). The standard deviations for the other angle (10 o , 30 o ) and percent speed errors (10%, 30%), are similar to these shown for the 20 o angle and 20% speed error respectively. For each case (e.g. 20 local estimates, 10 degrees angle error) the results are generated using as sample points the F1-score values estimated from all the boundary reconstructions that occurred using 1000 different (randomly selected) local estimates placements. The line plots in Fig. 12a and 12b indicate that as the angle and speed error increase, the boundary reconstruction accuracy decreases. Another observation is that the difference, between the boundary reconstruction accuracies achieved for each angle or percent speed error value, gets smaller as the number of local estimates increases. This behavior is justified if we consider that for higher densities, the mean distance that the local front segments (which participate in the boundary construction) have to ''travel'' as we project them to the future in order to predict the diffusive objects boundary is smaller (they continuously get replaced by more recent estimates). Traveling smaller distances with erroneous evolution parameters implies smaller boundary location errors and therefore higher boundary reconstruction accuracy. This also explains why the standard deviation drops as the local fronts' density increases for all the considered angle and speed error cases.

Probabilistic boundary tracking of a continuous object
As discussed in Section II, the evolution characteristics (orientation and speed) of the local fronts {m i } that participate in the boundary's reconstruction algorithm are described by scalar values {φ i ,û i } randomly sampled from the corresponding Normal distributions (for i and U i ). This random sampling implies that if we run the boundary reconstruction algorithm several times, we will get similar but not identical boundary shapes. Based on this observation we developed a technique which allows us to consider the area occupied by a continuous object as a field of probabilities where every point has a probability to be affected by the continuous object. To calculate the probability field we do the following: We consider the continuous object's evolution area as a grid of square cells. A cell is assumed to be affected by the evolving continuous object when its center is located inside the object's area. To calculate the probability for a cell C to belong to the continuous object at a specific time instance t b , we run the boundary reconstruction algorithm n times and count the number n C of simulation scenarios that cell C is located inside the object's area. The probability for the continuous object to arrive at cell C can then be estimated as:   13 shows for experiment 1 the probability field at t b = 188min that produced after running 1000 times the boundary reconstruction algorithm of a specific local fronts' placement scenario with 20 local fronts per km 2 . As observed, the proposed boundary reconstruction algorithm exploits the uncertainties of evolution parameters of local fronts to characterize the uncertainty about the continuous object's boundary location. In this example the used angle and speed uncertainties of the local front segments in i and U i were set to σ i = 10 o degrees and s i = 10% respectively.

E. DIFFUSIVE HAZARDS WITH IRREGULAR SHAPE
In Experiment 2 we evaluate the proposed boundary reconstruction algorithm using diffusive hazards with irregular evolution patterns (e.g. non-geometric front shapes, large propagation speed variations). To simulate such hazards with realistic characteristics we exploited FLogA [38], a webbased interactive tool which allows to draw a forest area anywhere in Europe over Google Earth [40] insert fire ignition points (''hotspots''), define wind direction and speed scenarios, and then simulate and geo-animate the evolving wildfire.
Using FLogA [38] we defined a square forest area of 3km × 3km at Hymettus mountain in Attica Greece and simulated a wildfire scenario. The fire was initiated from a single ignition point which was placed at the middle of the considered 9km 2 forest area. The wind orientation and speed parameters were fixed within the forest area and their values were set to 0 o (with respect to the x-axis) and 2m/s respectively. Fig. 14 show four snapshots that help us visualize on Google Earth the ignition point as well as the corresponding part of the forest area that has been affected by the fire (continuous object is shown with red color) 50, 100 and 150 minutes after the ignition respectively. In Fig. 15 and 16 we present a summary of the results for experiment 2 following the same format as for Fig. 11 and 12 respectively.
Similarly to Experiment 1, we evaluated the proposed boundary reconstruction algorithm under different local front: a) densities ({10, 20, 30, 40} estimates per 1km 2 ), b) placements (1000 random placements), and c) orientation and percent speed error scenarios ({10 o , 20 o and 30 o } and {10%, 20% and 30%} respectively. In summary, in this experiment the mean boundary reconstruction accuracy (F1-scores) was on average smaller by 3.91% (standard deviation 1.15%) relatively to Experiment 1 considering the results from all simulation scenarios. This accuracy drop is justified considering the complex evolution behavior of the wildfire's front line (irregular accelerations/decelerations) which makes the boundary tracking a lot more difficult.
Similarly to Experiment 1 we have created a MATLAB animation [32] that is provided online and as supplementary material that shows how the proposed  algorithm can track accurately evolving continuous objects with complex boundary shapes and evolution behavior (Fig. 17). In Fig. 17a we can see the area affected by an evolving wildfire (pink area), the wildfire's ignition point (red dot), the 20 randomly placed local front estimates (magenta line segments) and their placement area (dashed square area of 1km 2 ). The simulation terminates when the wildfire covers the local fronts' placement area (152min after the ignition).
A local front model estimate gets activated and participates in the boundary reconstruction when its mid-point is reached by the wildfire's front line. With black color (see Fig. 17b -17f) we depict the set of the local front segments that participate (after the selection/fusion procedure, see Section III) to the boundary reconstruction algorithm. Based on local fronts evolution characteristics, we determine their space-time propagation (black segments on the blue curve, see also Appendix B) and using the proposed boundary reconstruction algorithm (see Section IV) we reconstruct the ''smooth'' B-spline curve (blue curve) that approximates the boundary of the wildfire. Fig. 17b -17f, demonstrate the ability of the proposed algorithm to track continuously with accuracy the irregular shapes of the boundaries which result from wildfire's complex evolution behavior, despite the very small number of the local fronts estimates used (only 20 in 9km 2 ). Notice that the boundary tracking algorithm underestimates the wildfire's boundary only in regions where no local front estimates are available (see Fig. 17b -17e). Furthermore, in Fig. 17f the underestimation of the wildfire's boundary is expected if we observe that the wildfire's boundary is located outside the sensors placement area and therefore its tracking is based mostly on ''outdated'' and far away local front estimate projections.
In real hazard tracking applications the topography at the hazard's site can skew the placement of the local front estimates in the field. To evaluate the accuracy of our boundary reconstruction algorithm under this assumption, we repeated Experiment 2 using highly skewed placements of the local front estimates in the 1 km 2 square area (dashed square in Fig. 17) generated as follows: We partition the 1 km 2 square deployment area into 4 equal square areas (0.25 km 2 each) by connecting the middle points of its opposite sides. Next we randomly select 2 of the 4 formed square areas and we randomly deploy there 80% of the local front estimates used in the specific local front estimate density scenario. The rest 20% of the local front estimates are randomly deployed in the other 2 non-selected squares. This procedure allows us to generate highly skewed placements of the local front estimates in the original deployment area (1 km 2 square area).
Using the placements we evaluated the proposed boundary reconstruction algorithm under different local front densities ({10, 20, 30, 40} estimates per 1km 2 ). As expected, the boundary reconstruction accuracy was higher in the square areas where the 80% of the local fronts were placed. However, an interesting finding is that the average boundary reconstruction accuracies (F1-scores) did not significantly drop compared to the nominal Experiment 2. More specifically, for each local front estimate densities scenario (i.e. {10, 20, 30, 40}) the percentage reduction of the average boundary reconstruction accuracy was {5.98%, 4.66%, 2.8%, 2.62%} respectively. The results clearly show that our boundary reconstruction algorithm is almost insensitive to the skewed placements of the local front estimates, and its sensitivity is further reduced as the density of the local front estimates increases.

F. COMPARISON TO CODA
As mentioned in the introduction Section, most of the reported WSN-based schemes use the locations of the boundary sensor nodes to implicitly determine the boundary of a continuous object. Therefore these schemes rely on a human expert in order to reconstruct the boundary curve of a continuous object. One notable exception is the work in [15] where the authors propose CODA -a Continuous Object Detection and Tracking Algorithm which approximates a continuous object's boundary using the convex hull polygon of the boundary node locations. For completeness we compared CODA to our method assuming the local front estimates density scenarios (10, 20, 30 and 40 sensor nodes deployed within 1km 2 and 9km 2 square areas for Experiment 1 and Experiment 2 respectively). For all density scenarios, CODA failed to track with accuracy the continuous object's boundary. The accuracy (F1-score) of CODA for Experiment 1 and Experiment 2 were on average smaller by 50.71% and 59.36% respectively, relatively to the proposed boundary reconstruction scheme. We remark that these results were extracted using for each density a large number (1000) of randomly drawn sensor node deployments.
We also repeated the aforementioned experiments with very high density sensor networks with 5 × 10 −4 , 10 −3 , 2 × 10 −3 and 3 × 10 −3 sensors km 2 corresponding to 500, 1000, 2000 and 3000 sensor nodes respectively deployed within a square area of 1km 2 . Even for these very high densities, the accuracy (F1-score) achieved by CODA for Experiment 1 and Experiment 2 were on average smaller by 1.55% and 7.79% percent, relatively to the proposed boundary reconstruction scheme. The smaller boundary reconstruction accuracy of CODA in this case can be explained if we consider its inability to predict the location of the boundary when the object moves beyond the sensor nodes deployment area (see animations [31], [32] and Fig. 10f, 17f). Moreover, the larger decrease of the boundary reconstruction accuracy for Experiment 2, can be justified if we consider that CODA can form only convex boundary shapes. Therefore, the more realistic non-convex shapes formed for example as a wildfire evolves (see Fig. 17f and animation [32]) cannot be tracked accurately using CODA.

VI. CONCLUSION
We have presented a novel approach which allows us to track accurately the boundary of an evolving continuous object using only sparse in space and time estimates of its local front characteristics provided by distributed sensors. Our continuous object tracking scheme has the following unique characteristics relatively to prior schemes: a) It first filters and fuses the information of the available local front estimates, and b) then uses the surviving subset of estimates to reconstruct a smooth curve approximating the evolving object's boundary; c) it requires a much smaller number of sensors to provide local front estimates than other methods; d) it supports predictive modeling of the boundary even in areas with limited number of sensors; e) it provides higher accuracy and remains robust to local front estimate errors. Moreover, by exploiting the inherent uncertainties associated with the local front estimates, it can provide a field representation of the continuous object, indicating for every location the probability to be reached by the continuous object at any given time. Realistic simulations demonstrate that the proposed method can track accurately the evolving boundaries of different types of continuous objects (e.g. with time-varying evolution characteristics and/or irregular shapes), using a realistic number of possibly noisy local estimates. Moreover, global boundary reconstruction remains robust to local front estimation errors.
We are currently investigating a method that will allow us to determine the number of continuous objects that may evolve simultaneously in an area of interest. Furthermore, we are developing methods to probabilistically associate local front estimates to their originating continuous object. Solving these challenging problems can help us extend the current framework allowing us to track the boundaries of multiple co-evolving hazards. Finally, we are working on the development of a novel method able to estimate for every boundary point (using the evolution characteristics of the sparse local front estimates) the following parameters i.e. latitude and longitude coordinates, speed, direction and acceleration. This information can be directly exploited by hazard-specific predictive models in order to re-calibrate their parameters and to improve their prediction accuracy.

APPENDIX A COORDINATES OF LOCAL FRONT END POINTS
Using the mid-point coordinates P i = (x i , y i ), the length l i = 2R i , and the orientationφ i of a local front estimate m i , we can calculate the coordinates of its end-points, , by solving the following system of equations: The first equation in (7) defines a line (ε i in Fig. 18) on which the local front's line segment lies. The second equation defines a circle (C i ) of radius R i = l i 2 , centered at the local front's segment mid-point P i the assumed sensor location. The solution of (7) provides the following closed form algebraic expressions which are used for the computation of the local front's end-point coordinates (

APPENDIX B SPACE-TIME EVOLUTION OF A LOCAL FRONT ESTIMATE
Let's assume that we want to determine the location of a local front m i (t i ) (with parameters estimated at time t i ) at a future time instance t j , where t j > t i . Using the speed realizationû i , we calculate the distanced ij ''traveled'' by m i if we assume that it ''moves'' for time interval equal to t ij = t j − t i .
Next, by solving the system of equations in (10) below, we can determine the mid-point's coordinates P i (t j ) = (x i (t j ), y i (t j )) of the space time projected local front estimate m i (t j ).
The first equation in (10) defines a line (ε i in Fig. 19) which is perpendicular to the local front and emanates form its middle point P i = (x i , y i ). The second equation defines a circle of radiusd ij centered at the local front's segment middle point P i (see Fig. 19).

FIGURE 19.
Space-time evolution of the local front m i (t i ) from time t i of its generation to some future time t j (t j > t i ).
Solving this system of equations provides two closed form expressions for computing the intersection points of line ε i and circle C i (see Fig. 19). From these solutions we accept only the point that lies in the half plane indicated by the direction parameter δ i (green point in Fig. 19) , y i (t j )), is the new location of the mid-point of m i (t j ). For the calculation of m i (t j ) end-point coordinates, , 2}), we apply the equations (8) of Appendix A.

APPENDIX C LOCAL FRONTS INFORMATION FUSION
We present the proposed information fusion technique applied to handle events of type 4 (see Section III, Step 2). This technique uses the information of two close by local front estimates m i (t f ) and m j (t f ) in Fig. 20 that satisfy the condition of event 4 to calculate the parameters of a new local front m f capturing locally the evolution characteristics of the continuous object's boundary at time t f .
The information fusion technique calculates for each local front estimate m k (t f ) (where k ∈ {i, j}), the time differences t k = t f − t k . Next, using the parameters of the speed distributions U k , it determines for each local front estimate m k the distance D k that it will cover if it moves for time interval equal to t k . Since the speed U k is described by a Normal distribution N (u k , s 2 k ), the distance D k will also be described by Normal distribution N (d k , ω 2 k ) with parameters d k = u k t k ω k = s k t k .
It is worth mentioning that as the time difference t k increases, the uncertainty (ω k ) about the local front's (m k ) VOLUME 8, 2020 traveled distance will also increase. Next, we calculate for each local front m k the distanced k , that it will cover if it moves with speedû k for time interval t k d k =û k t k .
Using the distribution for the orientation k and distance D k , and the realization valuesφ k andd k , we calculate the weights (see equation (14)) that will be used for the calculation of the fused model m f parameters .
For the calculation of the weights we consider that the local front estimate with the larger likelihood values ( k (φ k ), D k (d k )) should be trusted more.
Using the weights w k k ∈ {i, j}, we calculate the parameters of the new local front estimate m f as follows To estimate the evolution direction parameter δ f we use the local fronts' orientations (φ k ), evolution directions (δ k ) and the weights (w k ) and determine two vectors which: a) are perpendicular to the corresponding local front segments, b) point towards the corresponding local fronts' evolution directions and c) the ratio of their lengths is equal to the corresponding ratio of their weights w i w j . Next, using these vectors we calculate their resultant which determines the evolution direction of the new local front estimate m f . Using the equation of the line where the new local front estimate m f lies, we determine in which half plane (positive or negative) the vector of the resultant points to and assign the corresponding value (+1 or −1) to the direction parameter δ f .
To calculate the speed (orientation) distribution of m f we apply the following procedure: Using the weights w k k ∈ {i, j} and the corresponding speed (orientation) distributions U k , ( k ) of the local front estimates, we find their Gaussian mixture (see Fig. 21): Next, by applying variational calculus we approximate the Gaussian mixture in (16) by the Normal distribution which minimizes the Kullback-Leibler (KL) divergence (maximizes the similarity) from the Gaussian mixture (see [41], [42]). The equations which can be used to compute the parameters of these Normal distributions are: Speed parameters Orientation parameters The parameters calculated from the above equations govern m f 's speed (U f ∼N (u f , s 2 f )) and orientation angle ( f ∼N (φ f , σ 2 f )) normal distributions after the fusion.

APPENDIX D INITIALIZATION OF LOCAL FRONT PARAMETERS
We explain here how we initialize the local front evolution parameters for the experiments presented in Section V-C, to simulate the measurements process of distributed sensors producing the local front estimates in the absence of errors (perfect sensing).

A. EXPERIMENT 1: CIRCULAR FRONT
Let's assume that the circular boundary of a continuous object reaches the mid-point (sensor location) of a local front segment (m i ) at time instance t R (see Fig. 22a). We set the speed (u i ) parameter of m i to the speed of the circle's radius at time t R , and its orientation (φ i ) parameter to the angle formed between the x-axis and the circle's tangent line (ε i in Fig. 22a) that osculates from the middle point of m i .

B. EXPERIMENT 2: IRREGULAR FRONT
FLogA (Fire Logic Animator) is a web-based software tool which allows the user to define (draw) a forest area on Google Earth anywhere in Europe, insert interactively fire ignition points, simulate and geo-animate the behavior of the evolving  fire line under different weather conditions [38]. FLogA considers the forest area as a grid of square cells, with every cell experiencing different topography and weather conditions. The dimensions of the cells are defined by the user. To simulate the evolution behavior of a wildfire, FLogA accepts as input a set of raster ASCII files that contain information about the forest's topographic layers (slope, aspect, fuel model, fuel moisture), the prevailing weather conditions (wind speed and wind direction) in the forest's area, as well as the number and locations of the fire ignition points (''hotspots''). FLogA uses cellular automata like algorithms to predict for each cell of the grid information such as the expected time of fire's arrival, the fire line's speed, the fire line's evolution direction etc. (see Fig. 22b). We set the speed (u i ) and orientation (φ i ) parameters of a local front estimate m i at a sensor location P i , to be equal to the speed and evolution direction values of the hazard for the forest cell containing the mid-point of the line segment of m i (sensor location).

APPENDIX E ASSESSING SIMILARITY OF CONTINUOUS OBJECTS
To evaluate the accuracy of the proposed boundary reconstruction algorithm, we compare the similarity of the areas occupied by the real and the reconstructed continuous object as described below: We consider the continuous object's evolution area as a grid of square cells (see Fig. 23). For all conducted experiments we set the size of the square cells to be 20m × 20m. Based on the real and reconstructed boundary of the continuous object, we classify the square cells to the following categories: • True Positives (TP): Cells located inside the real and the reconstructed object.
• False Positives (FP): Cells located outside the real but inside the reconstructed object.
• False Negatives (FN): Cells located inside the real but outside the reconstructed object. To estimate the boundary reconstruction accuracy we used the F1-Score which is the harmonic mean of precision and recall F 1 = 2 · precision · recall precision + recall where precision = TP TP + FP , recall = TP TP + FN . The i-th local front model.

APPENDIX F TABLE OF SYMBOLS
The physical coordinates of the i-th line segment's mid-point. The circular local area of model m i at time t j .

R i
The radius of the circular area C i .

R i
The adjusted radius of m i circular area.
The area covered by m i during its space time evolution from time t i to time t j . t i The time instance of m i parameters estimation by a sensor located at P i . The time instance of a boundary reconstruction.
The model for distance that m k (t f ) will cover if it moves for time interval equal to t k . d k The distance that will be covered by m k if it moves with the sample speedû k for time interval t k . t k The time difference between t f and t k . w k The weight used for the information fusion of the local front models m i and m j .
The set of m f parameters which result after fusing the information of m i and m j . M A set that contains the parameters of all the local front models.

M b
A set that contains the parameters of all the local front models selected to be used to reconstruct the boundary at time t b . C A set that contains the mid-points of all local front estimates in M b . V An ordered set that contains the vertices of the convex hull polygon of C. VOLUME 8, 2020 Q A set that contains the mid-points in C that do not belong to the convex hull polygon. Q on A set that contains the points of Q located on the polygon's edges. Q out A set that contains the points of Q located outside the polygon's area. Q I out A set that contains the points of Q out that are centers of local front segments intersecting at least one edge of the polygon. Q T out A set that contains the points of Q out that are centers of local front segments located totally outside (do not intersect) the polygon. Q in A set that contains the points of Q located inside the polygon's area. Q I in A set that contains the points of Q in that are centers of local front segments intersecting at least one edge of the polygon. Q T in A set that contains the points of Q in that are centers of local front segments located totally inside (do not intersect) the polygon.

P(C)
The probability for the continuous object to arrive at cell C. n C The number of simulation scenarios that cell C is located inside the object's area. n Total number of simulation scenarios. N The minimum number of local front estimates required for launching a boundary reconstruction.