1 Introduction

The purpose of this paper is to describe a novel algorithm for computing a smooth surface from a given finite set of interlacing curves in R ³ satisfying a " geodecity condition” at the intersections. The surface computed is such that the curves are geodesics on it, and such that each cell in the net of geodesics is approximately a minimal Gaussian curvature (squared) surface [8]. It has been suggested and researched in the past that markings along geodesics on surfaces of objects are used by the human visual system to infer the type of object observed [11]. We also know that illustrators draw the projections of geodesic curves (which may intersect) in order to illustrate a specific object. The former study and the later observation have motivated the author to obtain the result in [8].

In [8] we deal with finding surfaces in R ³ , which are as close as possible to being flat and span a given contour such that the contour is a geodesic on the sought surface (formally on any surface containing it). We look for a surface which minimizes the total Gaussian curvature squared. We show that by a change of coordinates the curvature of the optimal surface is controlled by a PDE which can be reduced to the biharmonic equation with an easy-to-define Dirichlet boundary condition and Neumann boundary condition zero. We then state a system of PDEs for the function whose graph is the optimal surface. This result allows us to consider each cell in the net of geodesics separately. We regard the smallest net component corresponding to a " cell” as a contour. It is not smooth as we assumed in [8] , as it has corners, but we expect this gap not to have a considerable effect. Treating each cell separately allows a great simplification of the task of finding the complete smooth surface.

Works regarding the computation of surfaces with geodesic boundary curves appeared before in the context of computer-aided-design. The first work that we are aware of is [16], which had in mind a similar problem to the one that we are dealing with. The authors were looking to be able to process geodesic data acquired by a device that they are using. In this work they dealt with only two non-intersecting curves. This work was later followed by [4, 5, 6]. All of these works consider Coons patches and Hermite interpolation in order to interpolate an unknown surface bounded by geodesics. However, the outputted interpolated surfaces that they compute are not " consistent” surfaces in the sense that computing these surfaces for adjacent cells in a net of geodesic curves, will not necessarily result in a smooth surface, and this is where we find our approach significantly advantageous. The framework in which an energy is defined on the surface has the property that minimizing it on the possible surfaces constrained by the given geodesic curves, is equal to optimizing it over the possible patches in the net of the given curves. Intuitively, the energy used encodes the tension of the surface spanning the bounding contour, where minimizing it corresponds to the natural inchnation of many organic or synthetic materials which are inclined to stay as flat as possible, but allow stretching. We describe an algorithm to construct this approximately optimal smooth surface when given a net of curves satisfying the following condition at the intersections. The condition is that the principal normals (unit normals) of each two intersecting curves should agree modulo sign at the intersection, see [5, Section 4.3] for details. This is a necessary condition for the curves to be geodesics on a surface. For each cell, the algorithm finds an approximate solution for the PDE system presented in [8, Eq. 5]. This is done by considering a cell’s contour, heuristically fixing a plane such that the projection of the unknown minimal Gaussian curvature surface spanning the contour is with high probability one-to-one, solving a biharmonic equation for the conformal factor for the approximate optimal curvature by using a finite element method, and then solving the curvature Monge- Ampere equation for a function whose graph over the chosen affine plane has the prescribed curvature at each point and its trace equals the cell’s contour. Note that in [8] we work in isothermal coordinates (coordinates for which the surface’s metric can be written as ds ² = Φ· (dx ² + dy ² ), Φ a smooth real function on the parameterization domain) and we do not know how to compute the coordinate change maps from isothermal coordinates to " projection on a plane” coordinates. Inputting the approximate optimal curvature in Cartesian coordinates seems to work well, and most importantly, we are able to evaluate the error between the computed function and the solution of [8, Eq. 5]. This can be done by applying a discrete Laplace-Beltrami operator to evaluate the result. This will potentially allow refinement of the net of geodesics in regions where the Laplace-Beltrami operator applied to the coordinates of the projection on the plane function with respect to the computed surface, is significantly larger in absolute value from zero. In addition, we can check if near the given curves, at the seam of two cells, the error is low. We plan to to implement this in a subsequent work which will put the emphasis on the implementation side of the method.

The most important contribution of this work is its prospects in applications in visualization technology, architecture, sketch to 3d technology and computer vision. This was the initial motivation for this study. However, we are now exploring applications in physics and engineering.

2 The full PDE of a cell in the net

We consider a single cell in the net of geodesics and its contou Γr in R ³ . We need to heuristically fix an affine plane for which the projection of the minimal Gaussian curvature surface spanning Γ would be one-to-one. We then have a PDE system for the functions h, x, y, real functions defined on the region on that plane that is bounded by the projection of F. is the approximate optimal curvature that we compute after choosing an affine plane, its domain is the same as of h, x, y. We are specifically interested in the function h, whose graph over its domain is the minimal Gaussian curvature surface spanning F. In particular, its trace is F. The PDE system is stated in [8] and is the following:

Note that E, F, G are the elements of the first fundamental form of the surface we are looking for and therefore depend on h. Ω is the domain of h, x, y and as mentioned, is a subset of some chosen affine plane. Note that the projection of Γ will be different for different choices of planes, and so will be the boundary conditions in the next section, however the approximate computed optimal curvature will be very similar due to the theory in [8]. is the approximate optimal curvature in the case of this paper. It is a solution for the linear part of the Euler-Lagrange equation for the integral of the square of the Gauss curvature on the surface. It is constructed by solving a biharmonic equation Δ ² ƒ = 0 with appropriate boundary conditions. We refer to ƒ as the " conformal factor” of the corresponding Riemaniann metric, e ^2ƒ (da: ² + dy ² ). We are prescribing the Gauss curvature given by

The analysis which yields the above derivation relies heavily on the usage of isothermal coordinates on surfaces. For related theory refer to [15] and [17] .

3 Solving the biharmonic equation

We consider a single cell in the net of curves. One needs to guess an affine plane for which the minimal Gaussian curvature surface spanning the cell’s contour, can be projected one-to-one on. Heuristically, one can define the affine plane by any choice of three corner points of the cell. Let be the projection of the cell’s contour on the chosen affine plane, and let Γ be the cell’s contour. Let Γ be a parameterization of the cell’s contour which is the inverse of the projection on the affine plane. As shown in [8], the Dirichlet boundary condition of the biharmonic equation for the optimal curvature is given by at point for a parameterization (arc-length or not) of

For experimenting, We considered contours constructed by modulating a periodic function, such a sine or a cosine, on the unit circle in the x-y plane. For these type of curves, one can compute the derivative and therefore have the boundary condition for the biharmonic equation for the conformal factor. (In Figure 1 in Drawings), under each sub-figure we state the function modulated on the circle, which was used to form the contour.) We therefore had a well-defined biharmonic problem in the plane for each curve, that we had to solve. For each periodic function z(t), modulated on the unit circle, the Dirichlet boundary condition is given by | g ^ The Neumann boundary condition is zero, corresponding to the demand that the contour will be a geodesic.

In order to solve such an equation we used FreeFEM++ [10]. One needs to put the PDE (biharmonic, in our case) in a variational formulation. We looked for a formluation which would also yield the Laplacian for the solution, as both the solution and its Laplacian are needed in order to compute the curvature. Such a formulation was given in the Ciarlet-Raviart approach to biharmonic problems [3], which is also due to Bertrand Mercier [12].

They both dealt with homogeneous boundary conditions, but by applying an elementary modification to their approach, one can consider the two unknowns u and ω (that will come out be equal to — Δu). We are looking for and for such that:

Generallƒy can be a non-zero function when looking for Δ ² u = ƒ· In our case, ƒ = 0. Note that now the particular solution (say, u _g ) needed to construct is easy to generate (just assigning the value of g at the boundary nodes and zero at internal nodes). 4 The curvature Monge-Ampere equation

In the previous section we showed how to obtain the approximate optimal curvature of the minimal Gauss curvature surface spanning a geodesic contour. We use the term " approximate” because the biharmonic equation we solved is the linear part of the Euler-Lagrange equation for the considered energy (the square of the Gauss curvature) . The boundary conditions provided for the solver are exact, for the case that the contour is a geodesic on an ambient surface containing the surface that we are searching. We also remark that we input the curvature as if it were given in Cartesian coordinates, which may also lead to inaccuracies, however we can evaluate the result. We elaborate on this in the end of this section and in the next section.

We input the curvature and the contour to a finite-difference solver, which solves for a function which satisfies the Monge-Ampere curvature equation, and its trace is the inputted contour. The solver we used was generously made available to us by Brittany Froese Hemfeldt. The method used by the solver is described in [7] in general context and also utilized in Gaussian curvature context in [9]. Qualitatively, it provided very good results, in addition to being meshfree - there are no structural constraints on the data points. A survey of methods for solving numerically Monge-Ampere equations can be found in the introduction section of [13].

Explicitly, the solver solves numerically the foil wing equation:

The solver’s input is the curvature on different nodes inside a domain, and the trace of h, e.g. its values on points on the boundary of the domain. Error estimates and empirical tests for the solver can be found in [7].

Note that in running the solver we inputted the curvature function found in Section 3 with no coordinate transformation. Thus we regarded the curvature as given in Cartesian coordinates. This is not compatible with the theory in [8], which is also mentioned in Section 2, but is used to compute an approximation of the target surface.

5 Applying the Laplace-Beltrami operator

We regard the approximate optimal curvature found in Section 3 to be in Cartesian coordinates and to correspond to the curvature of the graph of the height function at each point on the plane. Following that, we use the Laplace-Beltrami operator as described in [1] to obtain a value at each vertex of the mesh (a triangulation of the graph). The method described in [1] is suitable for closed discrete surfaces or for surfaces with boundary, whose boundary is a geodesic. It is therefore suitable in our case. The recent work by Burman et al. [2], which is a finite element approximation of the Laplace-Beltrami operator on a surface with boundary, can be used alternatively.

The definition of the discrete Laplace-Beltrami that we used, is Definition 16 in [1]. (The weights defined are the same as the original weights defined in the classic work by Pinkall and Polthier [14], but assumes a preprocessing step of an intrinsic Delaunay triangulation of the mesh). The Laplace-Beltrami operator is given by (the set of vertices of the mesh): where E _D is the edge set of a Delaunay triangulation of S and the weights are given by

The angles α _ij and α _ji are the α- angles of an internal edge. Refer to Figure 1 in [1] . The metric g is implicit, and only used for notation purposes.

The wishful result would be that when applying the discrete Laplace-Beltrami operator with respect to the constructed surface on coordinate projection functions π _x , π _y , we would observe uniformly close to zero values on the vertices consisting of the discrete surface. If the surface constructed is an approximate minimal Gauss curvature surface, this would be consistent with the defining PDE system, Equations 1. This is unrealistic to expect in all cases, as that would imply that a coordinates chart given by a projection map of a surface on the plane, is always an isothermal coordinates chart (see [17, 15] for a reference regarding the theory of isothermal coordinates).

Recall that we regard each cell separately. We use the values of a Laplace-Beltrami operator as an indication of how close our constructed surface is to the minimal Gauss curvature surface that spans the contour (with the boundary conditions being that the contour is a geodesic). In the case that there is a large deviation from zero in the values of the Laplace-Beltrami operator on the coordinate functions, we do the following: We compute a geodesic on the constructed surface emanating from two points of maximal distance on the cell’s contour. This geodesic defines two new cells’ contours. We can then construct surfaces for the two new cells, and discard the original constructed surface. The assumption is that the new geodesic that we find, is close to a geodesic on the true minimal Gauss curvature surface spanning the original contour. See Figure 2 in Drawings for computations of minimal Gaussian surfaces for two nets of curves, each consisting of two cells. Note that the " geodicity condition” does not hold, but it is a nice qualitative demonstration of the visually smoothness of the generated surface, even when the condition does not hold.

6 Method Summary

The intention of this section is to describe a complete potential method to create smooth surfaces from nets of geodesic curves. The theory and basic implementations of small scale examples presented in this paper, imply that such a system is feasible.

Given a net of curves in (smooth or discretized), we treat each cell in the net independently. We perform the following process for each cell in the net. Looking at a specific cell, the method finds an affine plane such that the projection of the cell’s contour is one-to-one on that affine plane. Let γ(t) be an arc- length parameterization of the projection of the cell’s contour. Let z(·) be the height function on the image of γ, such that h(t) := z(γ(·)) is the contour. One computes numerically and uses it to define the Dirchlet boundary conditions for the curvature’s biharmonic equation [8]. Explicitly the Dirichlet boundary condition is at γ(t). These values are computed numerically on nodes along the boundary γ(·). Alternatively, if one has an analytic expression for the velocity of the cell’s contour segments, this can be used instead of the numeric computation. Since the cell’s contour is not entirely smooth but has corners, a value corresponding to either one of the meeting segments should be assigned to the corner, or one can take the average of the two values.

The biharmonic equation described in [8] is solved by using the Ciarlet-Raviart/Mercier approach with a finite elements method. The Dirichlet boundary conditions are as described in the above paragraph. The Neumann boundary condition is zero everywhere. The solution to the biharmonic equation is the " conformal factor” aƒnd the Ciar let-Raviart/Mercier approach yields the Laplacian of ƒ as well. We then have the approximate optimal curvature of the minimal Gauss curvature surface spanning the cell’s contour - it is given by K(x, y) = —e ^-2ƒ(x,y) Δ(x, y).

We then input the curvature we found and the cell’s contour to a solver that solves for a function that its graph is a surface with the specified curvature and specified trace. We regard the curvature as given in Cartesian coordinates.

We then apply a discrete Laplace-Beltrami operator with regards to the computed cell’s surface on the coordinates projection functions π _x and π _y . If the application of the operator on the two functions, yields uniformly close to zero values, then we are done. If this is not the case, we split the cell into two cells. This is done by finding a geodesic on the computed surface (running from two maximal distance points on the original cell’s contour, for example), which together with the original contour defines two new contours. We discard the original surface (or just the " bad half’ of it) and compute new approximately optimal surfaces for the new contours.

7 Acknowledgments

I wish to thank Prof. Michel Bercovier from the Hebrew University for introducing me to FreeFEM++, and discussing different related topics with me.

I wish to thank Prof. Franco Brezzi for suggeting the Raviart-Ciarlet/Mercier formulation for the biharmonic equation, and its adaptation to non-homogeneous Dirichlet boundary condition.

I am deeply grateful and indebted to Prof. Brittany Froese Hemfeldt for her interest in the project, and most importantly her willingness to share her Matlab code with me.

References

[1] A. I. Bobenko and B.A. Springborn. A Discrete Laplace-Beltrami Operator for Simplicial Surfaces. Discrete Comput. Geom., 38:740-756, 2007.

[2] E. Burman, P. Hansbo, M.G. Larson, K. Larsson, and A. Massing. Finite element ap- proximation of the Laplace-Beltrami operator on a surface with boundary. Numer. Math., 141:141-172, 2019.

[3] P.G. Ciarlet and RA. Raviart. A Mixed Finite Element Method for the Biharmonic Equa- tion. In Carl de Boor, editor, Mathematical Aspects of Finite Elements in Partial Differential Equations, pages 125-145. Academic Press, 1974.

[4] R.T. Farouki, N. Szafran, and L. Biard. Construction of bezier surface patches with bezier curves as geodesic boundaries. Computer-Aided, Design, 41(11):772-781, 2009.

[5] R.T. Farouki, N. Szafran, and L. Biard. Existence conditions for coons patches interpolating geodesic boundary curves. Computer Aided Geometric Design, 26(5):599-614, 2009.

[6] R.T. Farouki, N. Szafran, and L. Biard. Construction and smoothing of triangular coons patches with geodesic boundary curves. Computer Aided Geometric Design, 27(4):301-312, 2010.

[7] B.D. Froese. Meshfree finite difference approximations for functions of the eigenvalues of the hessian. Num er. Math., 138:75-99, 2018.

[8] T. Gilat. Minimal Gaussian Curvature Surface. arXiv preprint: https://arxiv. org/ abs/2101.06673, 2021. [9] B. Froese Hamfeldt. Convergent approximation of non-continuous surfaces of prescribed gaussian curvature. Communications on Pure and Applied Analysis, 17(2):671-707, 2018.

[10] F. Hecht. New development in freefem++. J. Numer. Math., 20(3-4) :251-265, 2012.

[11] D.C. Knill. Perception of surface contours and surface shape: from computation to psy- chophysics. J. Opt. Soc. Am. A, 9(9):1449-1464, Sep 1992.

[12] B. Mercier. Numerical solution of the biharmonic problem by mixed finite elements of class C ⁰ . Boll. Un. Mat. Ital. (4), 10:133-149, 1974.

[13] M. Neilan. Finite element methods for fully nonlinear second order pdes based on a discrete hessian with applications to the monge-ampere equation. J. Comput. Appt. Math., 263:351- 369, 2014.

[14] U. Pinkali and K. Polthier. Computing discrete minimal surfaces and their conjugates. Experiment. Math., 2(1):15— 36, 1993.

[15] W. Schlag. A Course in Complex Analysis and Riemann Surfaces. Graduate Studies in Mathematics 154. American Mathematical Society, 2014.

[16] N. Sprynski, N. Szafran, B. Lacolle, and L. Biard. Surface reconstruction via geodesic interpolation. Comput. Aided Des., 40(4):480-492, April 2008.

[17] M.E. Taylor. Partial Differential Equations I: Basic Theory. Applied Mathematical Sciences 115. Springer- Verlag New York, 2nd edition, 2011.