THE BLAKE-ZISSERMAN MODEL FOR DIGITAL SURFACE MODELS SEGMENTATION

The Blake-Zisserman functional is a second-order variational model for data segmentation. The model is build up of several terms, the nature and the interaction of them allow to obtain a smooth approximation of the data that preserves the constant-gradient areas morphology, which are explicitly detected by partitioning the data with the graph of two special functions: the edge-detector function, which detects discontinuities of the datum, and the edge/crease-detector function, which also detects discontinuities of the gradient. First, the main features of the model are presented to justify the sense of the application of the model to DSMs. It is stressed the fact that the model can yield an almost piecewise-linear approximation of the data. This result is certainly of some interest for the specific application of the model to urban DSMs. Then, an example of its application is presented and the results are discussed to highlight how the features of the model affect the model outputs. The smooth approximation of the data produced by the model is thought to be a better candidate for further processing. In this sense, the application of the Blake-Zisserman model can be seen as a useful preprocessing step in the chain of DSMs processing. Eventually, some perspectives are presented to show some promising applications and developments of the Blake-Zisserman model.


INTRODUCTION
Variational models for data segmentation presented in (Mumford and Shah, 1989;Blake and Zisserman, 1987) allowed for a mathematical formulation of several significant problems such as signal and image segmentation (Vitti, 2012a,b), and material fractures analysis (Del Piero et al., 2007).We have to mention here that other classical variational models widely used to face the problem of segmentation are Total Variation Flow (e.g.Dogan et al. (2007)) and Anisotropic diffusion (e.g.Perona and Malik (1990)), which are essentially first order models.By segmentation we intend here the result of minimization of a specific energy, the choise of the energy components is properly made in order to obtain desired properties ot the minimizer (i.e. of the segmentation).In this work the Blake-Zisserman energy is used.An elliptic approximation of the energy gives rise to two auxiliary functions which represent, each one, an explicit partition of the image.By labelling those regions contoured by closed curves one can perform a segmentation in the classical way, this fact justify the terminology segmentation for a minimizer.In this work the results of the minimization procedure is presented and applied to a Digital Surface Model (DSM) of a mixed urban-agricultural area.By following the work by Zanetti (2013), the Blake-Zisserman variational model appears to be a proper model for the processing of an urban DSM.The features of the model we want to highlight are essentially two: first, the presence of both first and second order terms in the functional allow some kind of control in the smoothing of the data and in the edge/crease detection.Last but not least, an explicit detection of edges and creases of the data is possible due to the presence of the two special auxiliary functions.By describing the features of the functional and of its outputs, we try to show why the model can be considered a valuable tool for directly tackle the segmentation of urban DSMs but also as a preprocessing step useful to improve the results of further DSM processing.

Theoretical Background
An introduction to variational methods and a 1st order model.In the first order variational formulation for data segmentation proposed by Mumford and Shah, a function g : Ω → [0, 1], with Ω ⊂ R 2 , is the representation of the data, and the results of the segmentation process are a function u and a set K ⊂ Ω.The former is a regularized version of g and the latter represents the edges of the distinguishable objects in g.The basic idea is to identify the smooth approximation u and the edges of the objects of g by minimizing the functional: among any closed set K ⊂ Ω and any function u ∈ C 1 (Ω \ K).Parameters α, µ > 0 are introduced in order to properly adjust the characteristics of the minimizers, i.e., u, K. H 1 is the 1dimensional Hausdorff measure.
In the minimization process, the term Ω |u − g| 2 forces u to be close to g.On the other hand, the term containing |∇u| 2 induces the smoothness of the solution u on set Ω \ K, by neglecting small discontinuities (considered as noise) of g and allowing discontinuities only on K, which represents the edges of objects in g.The size of the set K is controlled by the term H 1 (K ∩ Ω) and the parameter α.
In optimization problems of this kind, both bulk and surface energies are minimized simultaneously and their supports are also undefined, being unknowns of the problem themselves.Consequently, the discontinuity set K is not known a priori, thus leading to the expression free discontinuity problems (De Giorgi, 1991;Ambrosio et al., 2000).
In practice, classical results of Calculus of Variations cannot be applied in this framework, since the formulation (1) turns out to be too strong.In order to prove the existence of minima a relaxation of the functional is necessary, so the problem is moved where now u ∈ SBV (Ω) and Su is the discontinuity set of u.
Because of the term H 1 (Su ∩ Ω), the functional is not differentiable, hence classical gradient descent methods cannot be directly applied for computing the minima.
In order to numerically compute a minimum, a Γ-convergence approximation (Braides, 2002b) of the functional with elliptical functionals defined on proper Sobolev spaces (Ambrosio and Tortorelli, 1992) is given.The properties of the Γ-convergence ensure that the obtained minima converge to a minimum of the functional F itself.
A 2nd order model: the Blake-Zisserman functional.
The Mumford-Shah functional, being a first-order model, presents some drawbacks (Blake and Zisserman, 1987;Vitti, 2012b).Firstly, some features of the data, such as creases (gradient discontinuities) are not sensed.Secondly, the gradient term leads to the so-called steep gradient over-segmentation, that is regions with very steep gradient are approximated by step functions, thus decreasing the definition of the solution u.In addition, the model leads to discontinuity sets composed of unions of C 1 arcs with at most 3-points junctions (in this case arcs are 2/3π wide).
In order to overcome these limitations, Blake and Zisserman (1987) proposed a second-order variational model: where As for the Mumford-Shah model, a weak formulation of ( 3) is necessary to prove the minima existence.The weak formulation of the Blake-Zisserman functional, given in the space of the generalized special functions of bounded variation, GSBV (Ω), is: A correspondence with minima for the strong formulation has been proved by operating a suitable identification (Carriero et al., 1997).Again, in order to compute the minima, a Γ-convergence approximation is necessary.By properly adapting the techniques developed for the Mumford-Shah case, Ambrosio et al. (2001) proposed the following functionals: defined on proper Sobolev spaces.The quantities α, β, µ are positive parameters, ξ , o are infinitesimals faster than , the Γconvergence parameter, and s, z : Ω → [0, 1].We see that in the formulation (5) every term of the functionals is an integral over the domain Ω, and that the lengths of the sets Su and S∇u \ Su are now approximated by the integrals involving two new auxiliary functions: s and z.In practice, the function u results to be a smooth approximation of the data, with the effect of the smoothing concentrated only on homogeneous portions of the data.On the other hand, the graph of each auxiliary function s and z forms a partition of the data.Let (u , s , z ) be a minimizer of (5).By letting → 0, Γ -convergence properties ensure u to be sufficiently close to a minimizer of (4).Moreover, thanks to the minimization (Ambrosio and Tortorelli, 1992)  Let us take a qualitative look to the choice of the functional parameters.Increasing values of µ will correspond to different minimizers where the closeness to the original datum g is forced, indeed high values of µ penalize the term squared distance from g.This is not a desiderable effect in noise-removal applications for which low values of that parameter should be used.On the other hand, in processing synthetic or low-noised images, high values of µ allow to an high-fidelity description of data features.
An high penalization, with parameter δ, of the second order term -the integral of the squared norm of the Hessian -locally force z to detect an information even when that norm is small compared to the other terms in the functional.As a consequence we have an high detection of second order informations of the image, such as creases and shadows.A similar approach on the parameters α and β will produce a control on the length of the countours of the constant gradient homogeneous areas (parameter α − β) and of the second order objects detected (parameter β).We can see that the segmentation features are strictly related to the choise of these parameters, which have to be chosen keeping into account that any parameter can't be stressed so much, in order to have a right and predictable behaviour of the result.

Proposed Algorithm
By the properties of Γ-convergence and by the regularity results described in the previous section, we know that, for small enough, a minimizer of the functional ( 5) is sufficiently close to a minimizer of F .The elliptic functional F is differentiable, furthermore we easily observe that it is quadratic with respect to the single variables u, s, z: So we expect, after a discretization of the energy, to obtain three symmetric linear systems associated to the partial gradients in the variables u, s, z.
Discretization of the model.
In order to minimize the functional we define the discrete representation of the rectangle Ω ⊂ R 2 with a lattice of coordinates.
The differential operators of the first and second orders in dimension two can be approximated using classic difference-schemes discretization.
Because of the cross terms z 2 |∇ 2 u| 2 and (s 2 + o )|∇u| 2 a classical difference-schemes method cannot be used.Anyway, the functional F is strictly convex in the directions (•, s, z), (u, •, z) and (u, s, •), so we propose the following subsequent minimization procedure: given an initial triplet (u 0 , s 0 , z 0 ) we compute subsequently, until small variations of the functional, according to a given threshold, are achieved.
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey The gradient of the discrete approximation of ( 5) can be splitted in three components and written in a matrix form: The square matrices Au, As, Az are symmetric and positive-definite, hence, stationary points of the three functionals given above can be computed by solving the following linear systems with a preconditioned conjugate gradient method.
The first initial values for the minimization procedure are set to Investigations that motivate this procedure can be found in a forthcoming work by Zanetti and Ruggiero.
We just mention here that other strategies can been adopted to face the variational approximation of the so called free discontinuity problems (see Braides, 1998) and the numerical implementation of the Mumford-Shah and Blake-Zisserman functionals (Carriero et al., 2002;March and Dozio, 1997;Chambolle, 1999;Bellettini and Coscia, 1994;Ambrosio et al., 2001;Braides, 2002a).

RESULTS AND DISCUSSION
In the following, we present an exploratory application of the Blake-Zisserman model to a DSM (spatial resolution of 1m), of a portion (650m x 1250m) of a small town and of its surrounding landscape nearby the city of Trento, Italy.
An aerial orthophoto of the study area is given in Figure 1.In Figure 2 the differences between the data g and its smooth approximation u produced by the Blake-Zisserman model are plotted.Statistics of the differences: n. of pixels: 812500; minimum: -0.734; maximum: 0.963; range: 1.697; mean: 0.00610718; standard deviation: 0.0831212; mean of absolute values: 0.0569721.
In Figure 3 and 4 the graphs of the auxiliary functions s and z are given.We recall here that the function s maps the discontinuity of u and the function z map the discontinuities and the discontinuities of the gradient of u, the smooth approximation of the given data.The values of the functions s and z are 1 (white) on homogeneous areas where the function u is smoothing the data, the values are 0 (black) in correspondence of the edges and creases of the objects in the scene.Grey areas on the graph of the function s correspond to those areas where the smoothing effect is less, the amplitude of the un-smoothed variations depends on the ratio between the values of the functional parameters.
In Figures 6, 7, 8, and 9 the same information (i.e.aerial orthophoto, map of the differences between g and u, and graphs of the auxiliary functions s and z) refers to a small portion of the study area.Observing Figures 2 and 7, it is possible to see how the data has been smoothed in correspondence of the raws of the trees cultivated on the surrounding of the small town.Somehow, a similar result can be obtained by applying the Mumford-Shah model directly to the orthophoto, (see Vitti, 2012b).
In Figure 5 a cross section of the data g and of its smooth approximation u is plotted along with the auxiliary functions s and z.The white line in the upper left corner of Figure 6 represents the trace of the cross section.It is possible to observe how the method smooths the data while respecting the discontinuities of    We remark that the output of the Blake-Zissermamn model is somehow twofold.On one hand, the model produces an approximation u of the data g that is smooth just within regions presenting a certain level of homogeneity, given by the nature and interaction of the terms in the functional.The discontinuities, and hence the region boundaries, are therefore preserved from being smoothed out.On the other hand the method, by means of the two auxiliary functions s and z, directly detects discontinuities such as region boundaries, edges, and creases.From this point of view, the Blake-Zisserman model produces outputs of the same type of those given in general by either region or edge based methods.

CONCLUSIONS AND PERSPECTIVES
The multiple outputs of the Blake-Zisserman relaxed model, namely the functions u, s, z, represent a particular segmentation of a given DSM.In u the noise is filtered by preserving constant-gradient areas morphology, while a map of the edges and creases in the data are explicitly detected by the functions s, z.The second order variational model by Blake-Zisserman overcomes the Mumford-Shah model in many ways that are of particular relevance when such functionals are applied to DSM rather than to images.In this context it must be noted that while the Mumford-Shah functional, being a first-order model, can lead to a piece-wise constant approximation of the data, wereas the Blake-Zisserman functional can lead to a piece-wise linear approximation of the data.Moreover, the Blake-Zisserman model, by means of the function z, explicitly detects the discontinuities of the gradient of the data.These two facts motivate the authors to develop a currently ongoing work on the use of the Blake-Zisserman model to face the problem of roof-detection by exploiting not only the features of the solution u but also those of the function s and z, to detect and distinguish each single roof pitch.Certainly, a systematic assessment of the performance of the Blake-Zisseman model needs to be carried out.The advantages related to the second order nature of the functional could be also evaluated in a quantitative way by comparing the results that can be obtained by further processing the outputs of the Blake-Zisserman model with higher level methods such as those used in building reconstruction models or in DTM generation from DSM models.Following a work on the implementation of the Mumford-Shah model by using a finite-element approximation over an unstructured mesh (Bourdin and Chambolle, 2000), it could be interesting to implement the Blake-Zisserman model in a similar way for the direct processing of ALS data.Eventually, we mention that recently (Carriero et al., 2012) have presented a second-order variational model in the framework of image inpainting for the recovery of locally damaged images.This application suggest to investigate the feasibility of the application of such a model for the processing of LiDAR data, for example for the recovery of surfaces degraded by object-removal operations (e.g., ground surface and trees removal).

Figure 1 :
Figure 1: Aerial orthophoto of the study area (spatial resolution of 0.5m)

Figure 3 :Figure 4 :
Figure 3: The graph of the auxiliary function s mapping the discontinuities of the smooth approximation u of the data [white=1, black=0]

Figure 5 :
Figure5: Cross section graphs from the particular in the Figure6the smooth approximation u and of its gradient.The loci of such discontinuities are clearly detected by the functions s and z.In Figures10 and 11a 3D view of a portion of the data g and of its smooth approximation u is given.It is possible to see how the noise has been filtered out and how the action of the smoothing has not affected the edges of the objects corresponding to the buildings facades and even the creases of the data corresponding to the roof ridges.

Figure 6 :
Figure 6: Aerial orthophoto of a portion of the the study area, the white curve is the trace of the cross section the term 1/4 (s − 1) 2 induces s to be 1 almost everywhere over Ω, except where big values of |∇u | 2 are achieved (i.e., in presence of a discontinuity), in this case s 2 (hence s ) is forced to be 0. The term |∇s | avoids big oscillations of the function s .A similar argument stands for z , which goes to 0 in correspondence of big values of |∇ 2 u | 2 (i.e., in presence of a discontinuity of the gradient).