Filling holes under non-linear constraints

In this paper we handle the problem of filling the hole in the graphic of a surface by means of a patch that joins the original surface with C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {C}^{1}$$\end{document}-smoothness and fulfills an additional non-linear geometrical constraint regarding its area or its mean curvature at some points. Furthermore, we develop a technique to estimate the optimum area that the filling patch is expected to have that will allow us to determine optimum filling patches by means of a system of linear and quadratic equations. We present several numerical and graphical examples showing the effectiveness of the proposed method.

CAGD, Earth Sciences or image reconstruction (Li and Li 2010;Ju 2009;Wang and Hung 2012), physics (Yang et al. 2017;Dong and Cao 2015), computer vision in robotics (Carr et al. 2001;Caselles et al. 2008), etc. A particular case of the above-referred research field is the one of filling the hole of a given surface, i.e., of determining a patch reconstructing a piece of a surface that for some reason is unknown or not properly defined. Several works regarding this specific topic define the filling patch to join the surface to be filled with some smoothness conditions and minimizing some quantity related to the geometrical shape or to a physical feature of the patch to, somehow, obtain a simple and fair filling [see e.g. (Caselles et al. 2008;Fortes et al. 2011;Zhong et al. 2021] and references therein). Nevertheless, these approaches based on the minimization of some quantity often lead to 'flat' patches, unable to retrieve complex shapes or even be faithful to the shape of the original surface to be filled. A simple example is the one of filling the hole on the top of a semisphere or a cone: depending on which the role of the filling patch is meant to be, we may want it to be somehow minimal, or we may want it to be faithful to the holed surface, or even more generally, to adapt to some specific prescribed shape or geometrical feature. Some papers in the literature develop techniques to overcome the problem of the flatness of the filling patches that some methods provide, like biharmonic optimization (Smurygin and Zhurbin 2015); transfinite interpolation (Dyken and Floater 2009); algorithms to handle weakly defined control points by means of B-spline surfaces (Weiss et al. 2002); or functionals involving some geometric features required for the filling patch, like Fortes et al. (2015) or Fortes and Medina (2022), where the patch is forced to somehow 'inherit' the shape of the surface to be filled, or Fortes et al. (2017), where a volume condition over the filling patch is imposed. Most of the additional geometrical constraints imposed with the aim that the filling patches be faithful to a prescribed shape or pattern are linear, in the sense that the associated problems lead to determinate linear systems. Nevertheless, filling patches fulfilling non-linear constraints are also useful in several researching fields like image processing, where curvature constraints are used to develop some models [see e.g. (Ambrosio and Masnou 2003;Goldluecke and Cremers 2011;Schoenemann et al. 2009) and references therein], or biological problems, where area constraints are considered in the frame of biological cell membranes [see e.g. (Colli and Laurençot 2012) and references therein].
In this paper, we consider the problem of defining fair filling patches under the two specific non-linear constraints pointed out above: having a prescribed area and having prescribed values of the mean curvature at some points. In both cases, we consider the filling patches to be quadratic Powell-Sabin splines, which allows us not only to join both the patch and the original holed surface with C 1 -smoothness, but also to simplify the expressions of the area and of the mean curvature of the filling patches insofar as both geometrical features can be expressed in terms of the coefficients of the expansion of the quadratic Powell-Sabin splines in the basis of the corresponding vectorial space. Another reason to consider quadratic splines is related to the fact that increasing the degree or the smoothness of the fitting splines implies handling much larger linear systems that are worse conditioned, adding further to the difficulty of solving them numerically. Besides, in the case of the area constraint, we consider Bézier techniques that will allow us to approximate the quadratic Powell-Sabin splines by means of triangular patches, leading to simpler expressions of the area. Moreover, we develop an ad-hoc algorithm to determine an optimum value of the area that the filling patch is expected to have when there is no prescribed value of the area. The joint implementation of both, Bézier approximation and determination of the optimum area, allows us to translate the filling problem with area constraint into a system of linear and quadratic equations, that have been already highly explored in the literature. Therefore, not only do we fill under an arbitrary  Powell-Sabin sub-triangulation (right) area constraint, but we also lean on results in quadratic equations to obtain area-constrained 'optimum' filling patches.
It is worth to mention that the systems of equations associated to the constraints considered in this paper are non-linear, infinite filling patches fulfilling such constraints can be found, but of course not all of them are fair enough. To get proper filling patches we will look for the ones minimizing a relative error balancing the bending energy of the filling patch -in order to get surfaces with no roughness or irregularities-and the non-linear constraint considered.
The outline of the paper is as follows: In Sect. 2 we give a brief description of all the preliminaries and basic concepts that will be used throughout the work and we describe the general frame of the problem to be considered. In Sect. 3 we introduce the main tools used to handle the non-linear systems arisen. In Sects. 4 and 5 we consider the non-linear area and mean curvature constraints, respectively. In both cases, we give numerical and graphical results and we include an analysis of the numerical methods leading to the best results. In Sect. 4 we also include the method to determine an optimum value of the area of the filling patch when there is no particular value to impose. Finally, we end by presenting some concluding remarks.

The Powell-Sabin triangulation
Let D ⊂ R 2 be a polygonal domain (an open, non-empty, connected set) in such a way that D admits a 1 -type triangulation (see Fig. 1 left), defined as the ones induced by integer translates of x = 0, y = 0 and x + y = 0 [see e.g. (Davydov et al. 1998)]. Given a 1triangulation T of D, we will consider the associated Powell-Sabin triangulation T 6 of T [see e.g. (Laghchim-Lahlou and Sablonnière 1996)], which is obtained by joining an appropriate interior point T of each T ∈ T to the vertices of T and to the interior points T of the neighbouring triangles T ∈ T of T . When T has a side lying on the boundary of D, the point T is joined to the mid-point of this side, to the vertices of T and to the interior points T of the neighbouring triangles T ∈ T of T . Hence, the six micro-triangles inside any T ∈ T have the point T as a common vertex. There are several ways to consider appropriates points T (Powell and Sabin 1977), nevertheless, a good choice (Sablonnière 1987) is considering T to be the incenter of T , for all T ∈ T (see Fig. 1).
It is well known (Powell and Sabin 1977) that given the values of a function f (defined on D) and the ones of its first partial derivatives at all the knots of T , there exists a unique S in the spline space where P 2 stands for the space of bivariate polynomials of total degree at most two, such that the values of S and the ones of its first partial derivatives coincide with those of f at all the knots of T .

The hole-filling problem
To later present the techniques developed to fill holes by means of non-linear constraints, we first introduce the general notation that will be considered throughout the paper.
Let H (the hole) be a connected and nonempty subset of D (see Fig. 2, left) such that ∂ D ∩ ∂ H = ∅, where ∂ X stands for the boundary of the set X . If H was not connected, the techniques developed to fill one connected hole would be applied to each of the connected components of H . Let T be a 1 -type triangulation of D, with associated Powell-Sabin triangulation T 6 , and consider H * is a polygonal domain surrounding H (see Fig. 2, middle). The reason to extend the original hole H to the polygonal one H * is because the filling patches to be considered are splines defined over triangulations. Of course, H * tends to H as the triangulation T becomes finer, which, furthermore, must be fine enough to have ∂ D ∩∂ H * = ∅. Let us consider the 1 -type triangulation be the spaces of the Powell-Sabin splines of D− • H * and H * associated to the triangulations T D−H * and T H * , respectively.
The problem we are going to consider in this work is: We want to fill the hole in the graphic of f over H * by means of a C 1 -function in such a way that: i) f = s f be as close as possible to f over D − H * ; ii) f = σ s f fills the hole of f over H * with some desired geometric non-linear properties.
Function s f over D − H * is the only one [see (Barrera et al. 2008a), Proposition 1] minimizing the 'energy functional' defined on S 1 2 (D, T 6 ) by ; ρ is the evaluation operator ρ(v) = (v( p 1 ), . . . , v(p q )); P = {p 1 , . . . , p q } is a subset of points in D − H * ; τ 1 ≥ 0 and τ 2 > 0. The first term of J 1 measures how well v approximates the dataset { f ( p i )} q i=1 (in the least squares sense), while the second and the third ones represent the 'minimal energy condition'. In Barrera et al. (2008a) it is shown that s f can be expressed as On the other hand, the filling patch σ s f over H * will be defined to fulfill three conditions: First, it must join s f with class C 1 . To this end, σ s f must belong to the set S 1 2 (H * and {t 1 , . . . , t s } are the knots of T laying on the boundary of H * (see Fig. 2, right). Secondly, σ s f is required to minimize the functional J 2 : where |v| m,H * = |β|=m H * ∂ β u(x) 2 dx 1/2 and λ 0 ≥ 0, i.e., we want to control the bending energy of σ s f in such a way that its graphics does not have roughness or irregularities. In Theorem 1 of (Fortes et al. 2011) it is shown that there exists a unique σ s f ∈ S 1 2 (H * s f ) minimizing J 2 which has the expression where are the functions of the usual Hermite basis of S 1 2 (H * ) associated to the boundary knots {t i } s i=1 , i.e., the ones verifying {w 0 1 , . . . , w 0 3s , w 1 , . . . , w n } is the extension of B ∂ H * to the usual Hermite basis of S 1 2 (H * ) and the vector α = (α j ) n j=1 is the solution of the system of linear equations where for t = 1, . . . , n. For the sake of simplicity, we will also consider the matrix form of (5): where and T denotes the transposition operation. Finally, we will impose σ s f to fulfill N additional geometrical constraints, which will lead to the N additional non-linear equations To obtain the filling patch σ s f , we will handle simultaneously the equation systems (6)-(7). In Sects. 4 and 5 we consider particular cases for f n+i .

Solving non-linear equations systems with optimization techniques
The solution of a system of non-linear equations has close connections with non-linear optimization problems. In fact, it is well known that the solution of is the global minimum of g(x) , with g(x) = (g 1 (x), . . . , g m (x)).
With this motivation, to carry out the experiments reported in Sects. 4 and 5, we have considered different reformulations of the original non-linear system of equations as an optimization problem. In all cases, g i (α) = f i (α) − d i for all i = 1, . . . , n + N , and μ ∈ (0, 1).
• Reformulation O 1 : The most natural way of addressing the non-linear systems of equations is by solving it directly: Variants can consider the absolute value of components, (9) or penalize componentwise errors, The parameter μ allows to allocate different weights to the linear and to the nonlinear constraints. For solving these systems of non-linear equations, we considered the functions lsqnonlin and fsolve, from the Optimization toolbox of Matlab Matlab (2022), and knitro_nlneqs and knitro_nlnlsq, from Knitro Knitro (2022).
Inspired by the above-mentioned relation between the solution of systems of equations and optimization problems, the different equations can be scalarized into a single function that needs to be minimized: Penalizing deviations by squaring them also removes the non smoothness associated with the previous function, promoting the success of derivative-based solvers in addressing the problem: In this case, function knitro_nlp from Knitro Knitro (2022), GlobalSearch from the Global Optimization toolbox of Matlab Matlab (2022), and the solvers SID-PSM Custódio and Vicente (2007); Custódio et al. (2010) and GLODS Custódio and Madeira (2015) were considered. The last two are derivative-free optimization solvers, thus more likely to succeed in the formulation using the absolute value. Additionally, GLODS and GlobalSearch are global optimization solvers, allowing a more thorough exploration of the variable space.
The solution of the system of nonlinear equations can also be regarded as a multiobjective optimization problem: However, considering that n and N can be large, the linear and the nonlinear constraints were independently aggregated, thus generating the biobjective problem: A simple variant penalizes the errors, by squaring each one of the components: There are not too many solvers available to address biobjective optimization problems. We considered BoostDMS Brás and Custódio (2020), MultiGLODS Custódio and Madeira (2018), and paretosearch Matlab (2022), the last from the Global Optimization toolbox of Matlab. In general, rather than a single point, the solution of a multiobjective optimization problem is a set of points, which constitutes the Pareto front of the problem. If the solution of the non-linear system of equations is unique, then the Pareto front will be a singleton. Usually, that will not happen and a solution needs to be selected from the final approximation to the Pareto front. The parameter μ can be used for scalarization of the Pareto points, as in reformulation O 2 , allowing the selection of a single solution.
Finally, a Chebyshev formulation can be considered by minimizing the worst component of the system of non-linear equations: . As usual, deviations can be penalized: In this case, solvers SID-PSM Custódio and Vicente (2007); Custódio et al. (2010), GLODS Custódio and Madeira (2015), functions GlobalSearch and fminimax from Matlab Matlab (2022), and knitro_nlp from Knitro Knitro (2022) were attempted to solve the problem. Observe that, with exception of the multiobjective solvers (paretosearch, BoostDMS, and MultiGLODS), with the aim of solving the non-linear system (6)- (7), all solvers were provided with an initialization vectorα = (α 1 , . . . ,α n ). To this end, let {ξ i } n 3 i=1 be the set of knots of T laying in the interior of H * and let {w 3(i−1)+t } 3 t=1 be the usual Hermite basis functions of S 1 2 (H * ) associated to the knot ξ i , for i = 1, . . . , n 3 . Then, any On this basis, we have considered the initialization vectorα to be ⎧ ⎪ ⎨ ⎪ ⎩α whereĥ is the filling function of the test function h when considering the C 1 -basic filling method developed in Fortes et al. (2011), i.e., somehow we lean on the 'basic' j w j to obtain filling patches fulfilling the additional nonlinear constraints.
graphic of the filling σ s f defined in (1) to have a given area A over H * . The expression of the area Ar ea(σ ) of the graphic of a function σ ∈ S 1 2 (H * s f ) is not easy to handle as it implies to consider a great number of integrals of square roots of quadratic functions over triangles. Hence, we will consider the unique non-linear f n+1 in (7) to be a suitable approximation of Ar ea(σ ) and d n+1 to be the prescribed area (5) and d n+1 in (7) are strongly dependent insofar as high prescribed area values A will of course increase the bending energy of σ , moving the values . Anyway, fixed the prescribed area A, we are still interested in minimizing the bending energy of σ , i.e. in minimizing , in order to obtain a filling patch σ s f as smooth as possible, without roughness or irregularities. In this frame, the role of the parameter μ appearing in the objective functions O i introduced in Sect. 3, balancing the weight given to the linear part {g i (α)} n i=1 and to the non-linear one g n+1 (α) is quite trascendental as it should lead to a proper balanced pair bending energy-area.
For the sake of clarity, we will divide this section into several subsections: the first three handle the questions of estimating suitable values of μ, f n+1 and prescribed A, respectively. In the two final sections, we present the general setting under which the experiments have been carried out and the numerical results, respectively.

Estimation of the optimum parameter
To be able to estimate 'optimal' values of μ, we have considered 1600 where each α j has been randomly chosen in the interval [−2.5, 2.5]. For each one of these functions σ we have computed its bending energy J 2 (σ ) (3) and its area, and we have considered the point (Ar ea(σ ), J 2 (σ )). To shorten this procedure, the bending energy J 2 (σ ) has been computed by means of a numerical integration formula, exact for splines of order two, that evaluates in three points. In Fig. 3 we show the cloud of points obtained for the triangulation T and the polygonal H * shown in Fig. 2 right. These provide an insight of how the balance of the bending energy-area is. Next, given a prescribed value of the area A, we compute an estimation E(A) of which the bending energy associated with this area is expected to be by means of a quadratic regression function E (in red in Fig. 3). Then, the value of the parameter μ considered to be optimum is the one (5), regarding the bending energy, and the non-linear one g n+1 (α), regarding the area constraint, have the same weight in the objective functions O i considered, i.e.
In Tables 1, 2, 3, 4, 5 and 6 we give the values of the optimum parameters associated with the different experiments carried out. It can be observed that the higher the value of the prescribed area is, the less the optimum value of the parameter μ is. This is reasonable insofar as the bending energy of a function σ ∈ S 1 2 (H * s f ) grows faster than its area and, therefore, higher values of A require lower weights of the bending energy to achieve an equilibrium between them.

Approximation f n+1 of the area expression
To obtain a suitable approximation f n+1 of the area Ar ea(σ ) of a filling patch σ ∈ S 1 A thorough development of the control nets and of the properties that we will use in this work can be found in Farin (1986). For completeness, we briefly describe its construction. Given a triangle T , it is well-known that the Bernstein polynomials of degree m ≥ 1 over T , defined as where τ = (τ 1 , τ 2 , τ 3 ) are baricentric coordinates with respect to T (the multi-index notation λ = (λ 1 , λ 2 , λ 3 ) ∈ Z 3 is used, |λ| = λ 1 + λ 2 + λ 3 = m, and λ ≥ 0 indicates that λ i ≥ 0 for i = 1, 2, 3), form a partition of the unity over T and constitute a basis for P m (T ). For any p ∈ P m (T ), the unique representation is called the Bernstein-Bézier representation of p with respect toT . The coefficients b λ are called the Bézier ordinates of p.
, where ξ m λ are the Bézier triangle points of order m of T , defined as the ones having barycentric coordinates λ 1 m , λ 2 m , λ 3 m (see Fig. 4), are called the Bézier control points of p, and the linear interpolant B of the Bézier control points is called the Bézier control net of p (see Fig. 4 right). The values of the Bézier control net coincide with the ones of p at the vertices of T and, moreover, the control net is tangent to the graphic of p at these points (see Fig. 4 right).
If p ∈ P 2 , then p ∈ P 2+i for all i ≥ 1, and hence it can be written as The i th -elevation B (i) of the Bézier control net of p, defined as the linear interpolant of the control points {(ξ , lies in the convex hull of B (k) for all k = 0, . . . , i − 1 (we adopt the convention B (0) ≡ B), is tangent to p at the vertices of T , and verifies lim i→+∞ B (i) = p. Now, with the aim to obtain a filling patch over the polygonal hole H * whose graphic has a prescribed area A, we will consider the next process. Given a generic function σ = If for any elevation level k we consider a function σ k ∈ S 1 2 (H * s f ) verifying Ar ea(B σ k ) = A, then Ar ea(lim k→+∞ σ k ) = A. Therefore, for any elevation level k, we will consider i.e., σ k will be a function whose k th -elevation of its control net has area A. It is worth mentioning that the Bézier ordinates of a polynomial p ∈ P 2 (T ) can be obtained by means of the expressions This fact, together with (11), allows us to obtain the expression of the vertices of any of the triangular patches of B σ k in terms of the unknowns α. As a consequence, if we express the area of each triangular patch in terms of the coordinates of their vertices, we get that f n+1 where ≡ (k) is the number of triangular patches in B σ k .

Estimation of optimum area
All non-linear systems (6)-(7) associated to the experiments presented in this section have been carried out by means of the different formulations and solvers presented in Sect. 3. We have also considered different values of the prescribed area A to be set equal to d n+1 in (7). Nevertheless, we have found it appropriate to include a numerical procedure in order to estimate an optimum value of A when there is no particular value to impose. This numerical procedure consists of estimating the optimum value of the area that each of the triangles of the k-elevation of the control net of the filling patch σ H * for the triangulation T 5 . The procedure is as follows: (i) Let C be the centroid of the polygonal hole H * (in Fig. 5 right we show H * and its centroid for the hole defined in (17) and shown in Fig. 2). H * to C. In Fig. 6 we show, for k = 2 and for both test functions h 1 and h 2 , the cloud of points defined in ii) above (blue dots), and the points ( C T − C 2 , A opt T ) (red dots) for each triangle T ∈ N H * , where · 2 stands for the usual Euclidean norm in R 2 . The optimum value A opt In Tables 1 and 4 which can be equivalently written as where c ii if i = j. Problem (6)-(15) can now be restated as finding α that satifies This new explicit reformulation has theoretical and practical implications. On one hand, it allows for a better understanding of the underlying geometry since it describes the arrangement of intersecting lines and ellipses. More importantly, it sheds light on the nonconvexity of the problem due to the presence of these equality quadratic constraints. Also, given that the problem is inherently infeasible, solving (16) has to be thought in terms of minimizing an infeasibility measure of the constraints. Having a geometrical interpretation can help choosing an appropiate merit function that measures infeasibility and the best optimization solver to use. On the other hand, a nice property of explicit quadratic constrains is that solvers have access to first-and second-order derivatives practically at no computational cost, If q(x) = 1/2 x T Qx + c T x − d = 0 is a quadratic constraint, then ∇ x q(x) = Qx + c and ∇ 2 x x q(x) = Q are the first-and second-order derivatives respectively. Thus, dedicated quadratically-constrained solvers by design exploit the explicitly revealed structure in (16), as we will see in the experimental part.

Numerical and graphical examples settings
All examples of filling patches under area constraints presented in this section, except the last one, have been carried out over the domain D = (0, 1) × (0, 1); with the hole H defined implicitly by shown in Fig. 2 left, and with the 1 -type triangulations T m of D defined as the ones associated to uniform partitions of [0, 1] into m subintervals, for m = 5, 9 (T 5 is shown in Fig. 2). The smoothing parameters in the functional (2) have been chosen to be τ 1 = 10 −4 and τ 2 = 10 −5 . These values have been checked to lead to proper fitting functions s f over D − H * -a deep study on how to choose these parameters values was undertaken in Barrera et al. (2008b)-, while the smoothing parameter in functional (3) has been taken to be λ 0 = 10 in order to give the same weight to the semi norms | · | 1 and | · | 2 over D − H * and over H * in the functionals (2) and (3), respectively. We have considered two test functions, whose graphics are shown in Fig. 7: · Franke's function: · Semisphere function:

Numerical results
To check the feasibility of the proposed method, we will consider, for both test functions, three different values of the prescribed area: the optimal one A opt , obtained with the method previously developed (14), and another two values A + and A ++ , greater than A opt . All experiments have been carried out over the triangulations T 5 and T 9 and with the elevations k = 0, 4, 7 of the control net. We will denote by the best solution of (6)- (7) when handling the test function h i with prescribed area A and for the non-linear f n+1 as in (12), where the 'best solution'α = (α 1 , . . . ,α n ) has been considered to be the one minimizing the relative error where σ s f ∈ S 1 2 (H * s f ), defined in (4), is the unique spline in S 1 2 (H * s f ) minimizing the bending energy J 2 , i.e., σ s f is the C 1 -'minimum' filling patch of s f , in such a way that the first summand of (18) gives a measure of the relative bending energy of the spline 3s i=1 ϕ i (s f )w 0 i + n j=1 α j w j with respect to the minimum one, while the second summand is a measure of the relative error associated with the area constraint. Therefore, somehow we define the best solution to be the one exhibiting an equilibrium between a relative minimal bending energy and a relative minimal error with respect to the area A to be achieved. Of course, as expected, we obtain very bad results for values of A lower than a certain quantity (in fact, A cannot be less than the area of the polygonal H * , these experiments result in very high values of ( f n+1 (α)−A) 2 A 2 ). In Tables 1, 2 and 3, we show results for Franke's function h 1 , while in Tables 4, 5 and 6, we show results for semisphere's function h 2 . In all cases, we give the area A(σ Tables 1 and 4 we also include an estimation of the relative error where {a 1 , . . . , a 1500 } are random points in H * , when filling both test functions with the 'optimum' splines σ (k) h i ,A opt obtained with the optimal area constraint A opt , as well as the exact area of the test function A(h i ) over H * .
It is well known that the higher the elevation level of the control net of a bivariate quadratic polynomial is, the better its area approximates the area of the polynomial. All our experiments showcase this behaviour and can be observed that the relative error E decreases and the area constraint is better satisfied as the elevation level k increases.
We can also observe that the results obtained from triangulation T 9 are better than the ones associated to T 5 . Of course, this is also reasonable since disposing of a higher dimensional spline space where to look for filling patches must lead to better results. Although, of course, raising dimensions turns into higher computational costs.
In the case A = A opt , the best results, reported in Tables 1 and 4 and in Figs. 8 and 9 left, have been obtained when considering the non-linear system of separable equations (6)-(15).  Table 2 Numerical results for Franke's test function for A = A + = 0.8 h 1 ,0.8 ) = 0.792 h 1 ,0.8 ) = 0.798 Table 3 Numerical results for Franke's test function for A = A ++ = 1.8 h 1 ,1.8 ) = 1.781 h 1 ,1.8 ) = 1.787 Table 4 Numerical results for semisphere's test function for A = A opt k = 0 k = 4 k = 7 T 5 μ opt = 2.14 · 10 −2 A(h 2 ) = 0.3656 h 2 ,A opt ) = 0.2567 E = 9.87 · 10 −5 Table 5 Numerical results for semisphere's test function for A = A + = 0.7 h 2 ,0.7 ) = 0.679 h 2 ,0.7 ) = 0.689 Table 6 Numerical results for semisphere's test function for A = A ++ = 1.3 h 2 ,1.3 ) = 1.289 In practical terms, Problem (16) can be numerically solved using interior point (or barrier) methods (IPM) (Wächter and Biegler 2006;Byrd et al. 2000Byrd et al. , 1997Conn et al. 2000;Gould et al. 2001). IPM define a family of iterative methods that can solve convex or nonconvex nonlinear programming problems of the form where f (α) : R n → R is the objective function and c(α) : R n → R m describes a set of m ≥ 0 constraints. The algorithm consists of outer and inner loops. Given that (19) has no inequalities nor bound constraints, in the outer loop a sequence of (unconstrained) problems similar to are solved where λ is a Lagrange multiplier vector. These problems are generally solved using Newton's method and thus require to solve a linear system of equations involving the Hessian of the Lagrangian of (19) to find a suitable descent direction. For our case and since we are interested in finding a feasible point, there is no objective function and the constraints function c(α) = 0 is described in (16). Only linear and quadratic constraints are involved and the Hessian of the Lagrangian contains a block of constant elements. Quadratically constrained IPM methods exploit this fact and can keep a factorization of this block without the need to refactorize at each outer iteration. The inner loop consists of a line search using e.g. exact penalty functions or filter methods (Wächter and Biegler 2006). State-of-the-art IPM solvers such as NAG's e04st [40] or IpOpt Wächter and Biegler (2006) solve problem (16) by finding a feasible point. In this case, the inner iteration linesearch minimizes a merit function, f I (α), defined as a metric that measures infeasibility. Many metrics have been proposed, each with its own advantage. Sparsity-preserving 1 -norm, Euclidean 2 -norm or weighted variants of the previous two are all quite commonplace. We use an 1 -norm and the merit function chosen is f I (α) = i=1 |c i (α)|. This choice for f I (α) matches with the implemented one in the solvers e04st and IpOpt Wächter and Biegler (2006), with the former one used in the experimentation. It is worth mentioning that both implemententations achieve the best performance by making use of the first-and secondorder derivatives. Tables 1 and 4 report results for the numerical experimentation using the e04st solver.
On the other hand, regarding the case A = A opt , formulations O 1 , which directly attempt to solve the system of non-linear equations, are the most successful ones, when solved by fsolve, from Matlab, or knitro_nlneqs, from Knitro. It is worthy to mention that for high values of the prescribed area A, the minimal bending energy condition is much harder to fulfill, since the unique solution of the linear system (6) is naturally associated to a low value of A. From this point of view, we find it reasonable that the best results in the case , as it happens in the case A = A + . In fact, the results obtained for A = A + by means of the solution of the system = 0 are worse than the ones obtained with (8) and (9). Similarly, the results obtained for A = A ++ by solving the system are worse than the ones obtained with 10. These facts show that obtaining good results depends not only on the value of μ opt but also on the formulation considered.
In Fig. 8 above we show the graphics of the filling patches σ h 1 ,A of Franke's function obtained for triangulation T 9 and for the three different prescribed values of A considered in Tables 1, 2 and 3, respectively. In Fig. 8 below we show the graphics of the global reconstructed functions (1) all over D. Figure 9 follows the same pattern for the semisphere's function. In the case of the semisphere, due to the symmetry of its graphic around the top point (0.5, 0.5, h 2 (0.5, 0.5)), we have obtained, in all experiments, two different kinds of solutions: the ones leading to filling patches 'above' the graphic of the semisphere and the ones leading to filling patches 'below' the graphic. Since the components {α 3(i−1)+1 } n/3 i=1 of any solution α of (6)-(7) are exactly the values of the filling patch σ s f at the knots of the triangulation T lying in the interior of H * , it is easy to decide whether a solution is desirable or not. In the examples provided in this section we have chosen the solutions leading to filling patches 'above' the graphic of the semisphere to somehow keep its shape. Obviously, another solution leading to 'sunked' filling patches can be considered.
To end this section, we have considered an example over a hole H 0 with a more complex geometry. This hole, whose graphic is shown in Fig. 10 left, is inspired by the one considered in Chapter XI.4 of Arcangéli et al. (2004). This last example has also been carried out over the domain D = (0, 1) × (0, 1), for Franke's test function, and for the smoothing parameters τ 1 = 10 −4 and τ 2 = 10 −5 in the functional (2), and λ 0 = 10 in the functional (3). In order Fig. 9 Semisphere's filling patches with prescribed areas in T 9 for k = 7 h 1 ,A opt ) = 0.4334 E = 4.49 · 10 −4 to obtain a surrounding H * 0 close enough to H 0 , we have considered the triangulation T 25 of D. To perform this example, we have directly considered the non-linear system of separable equations (6)-(15), which in the previous ones led to the best filling patches. Numerical results of this experiment are reported in Table 7, while the filling patch is shown in Fig. 10 right. As in the previous examples, Table 7 shows that the relative error decreases and the area constraint are better satisfied as the elevation level k increases.

Filling patches satisfying curvature constraints
In this section, we consider the problem of determining minimal energy-filling patches with prescribed values of the mean curvature at some interior points of H * . Let us recall that the mean curvature of a bivariate function at a point (x, y) is defined as Fig. 11 Franke's filling results with prescribed mean curvature at one point where ∂ x , ∂ y , ∂ x x , ∂ xy and ∂ yy represent the first and second orders partial derivatives operators. So, if {χ q } q=1 is a set of points in • H * in which we want the filling patch to have prescribed values K = {K q } q=1 of the mean curvature, we will have to handle the non-linear system (6)- (7) where, for each q = 1, . . . , , the function f n+q in (7) has an expression of the form where we adopt the convention α 0 = 1, and d n+q = K q . Contrary to what happens in the area constraint case, since mean curvature is a local geometrical concept, we cannot expect to obtain any regression function providing a kind of dependence between the values (5) and the prescribed {K q } q=1 . In fact, it is easy to define a function having zero mean curvature at one point with a very high bending energy. Bending energy is expected to take values much higher than the prescribed ones of the mean curvature and, therefore, the values of the parametes μ in the objective functions considered must necessarily be low to get a properly balanced bending energymean curvature.
We have carried out experiments with both test functions Franke (h 1 ) and semisphere (h 2 ) and, in all cases, the best solution σ h i ,K ∈ S 1 2 (H * s f ) associated to prescribed values of the mean curvature K has been considered to be the one minimizing the relative error where σ s f ∈ S 1 2 (H * s f ), defined in (4), is the unique spline in S 1 2 (H * s f ) minimizing the bending energy J 2 . All the experiments in this section have been performed with the parameters values μ = 10 −1 , 10 −4 , 10 −7 , being the best results in all cases the ones associated to μ = 10 −4 .
In Fig. 11 we present results for Franke's function h 1 over the triangulation T 5 when imposing mean curvature K 1 = 0 at the point χ 1 = (0.5, 0.5) -central point in the graphics-. We show the graphic of σ h 1 ,K 1 (left), the one of the global reconstructed function (right) and the values of the mean curvature at χ 1 of h 1 , of the C 1 -minimal filling of h 1 obtained just minimizing functional J 2 (4), and of the function σ h 1 ,K 1 .
In Fig. 11 we observe that, although the mean curvature of σ h 1 ,K 1 at (0.5, 0.5) is closer to the prescribed K 1 = 0 than the ones of Franke's function and of its C 1 -filling, it is still a little bit far from being 0. This is due to the fact that we have considered the coarse triangulation T 5 . To improve this result, in Fig. 12 we consider Franke's function over the finer triangulation  T 22 when imposing mean curvature K 1 = 0 and K 2 = 0 at the points χ 1 = (0.36, 0.5) and χ 2 = (0.64, 0.5).
Following the same pattern that for Franke's function, in Fig. 13 and in Table 8 we show the results for different prescribed values K 1 of the mean curvature for the semisphere at the point χ 1 = (0.5, 0.5), while in Fig. 14 and in Table 9 we show the results for different prescribed values K 1 and K 2 of the mean curvature for the semisphere at the points χ 1 = (0.36, 0.5) and χ 2 = (0.64, 0.5), respectively.
Once all experiments have been carried out, we have found that for prescribed values of the mean curvature closest to the exact ones of the functions to be filled (semisphere with prescribed K 1 = 0, K 1 = −5 and K 1 = K 2 = −5), the best results are obtained by means of the solver knitro-nlneqs, using reformulation O 1 , which treats directly the nonlinear system of equations, considering the absolute value of the different components (9). For prescribed values of the mean curvature farthest to the exact ones (the remaining experiments), the best results are obtained by means of GlobalSearch, when considering reformulation O 4 , again considering the absolute value of the different components. We point out that solver GlobalSearch aims at global optimization, allowing a better exploration of the variable space.

Conclusions
We have developed a method to reconstruct holes in a given surface or, more generally, in a given dataset. The filling patch joins the original surface with C 1 -smoothness and it is defined as a spline of bivariate polynomials of total degree at most two, leading to minimum possible computational costs insofar as no lower degree allows to obtain C 1 -splines. Moreover, the filling patch is required to be a proper approximation of a non-linear system composed of two types of equations: the first ones (linear) control the fairness of the filling patch, while the second ones control non-linear features, as the area or the mean curvature. To obtain the filling patch, some reformulations of the original non-linear system and some solvers to address biobjective optimization problems are considered. In the particular case of the area constraint, we develop a method to estimate optimum values of the weight to be given to each one of the geometric characteristic (fairness and non-linear constraint) to be handled, as well as a method to estimate the optimal area to be imposed. The numerical and graphical examples presented show the effectiveness of the numerical methods proposed. Regarding the area constraint, experiments carried out show that the filling patch associated with the optimum value of the area is better achieved when the original non-linear system is reformulated in terms of several disaggregated equations numerically solved by means of interior point methods. While, in the general case, reformulations directly attempting to solve the system of non-linear equations are the most successful ones, when solved by fsolve, from Matlab, or knitro_nlneqs, from Knitro. Regarding the mean curvature constraint we also found that solvers and reformulations leading to the best results depend on the value of the mean curvature imposed. More precisely, solver knitro-nlneqs and reformulation O 1 work better for prescribed values closer to the exact ones, while for farthest values the best results are obtained by means of GlobalSearch, when considering reformulation O 4 .