Spatial constraint based triangular mesh operations in three dimensions

Abstract

A method and system provide the ability to design a (land) surface. A triangular surface mesh representative of an existing surface is obtained. The mesh includes triangles that are connected by vertices and edges. Design constraint sets are determined based design constraints. The design constraints include a maximum slope constraint for a first triangle of the two or more triangles in the triangular surface mesh. The maximum slope constraint is a maximum angle between a normal vector of the first triangle and a reference vector. Heights of the vertices of the first triangle are projected onto the design constraint sets such that the normal vector satisfies all of the design constraints. The projecting includes modifying the heights by a minimum Euclidian distance. A design of the surface represented by the triangular surface mesh is generated based on the projecting.

Description

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates generally to road design, and in particular, to a method, apparatus, system, and article of manufacture for designing a grading of a triangular mesh representing a physical site to accommodate proper drainage.

2. Description of the Related Art

(Note: This application references a number of different publications as indicated throughout the specification by reference numbers enclosed in brackets, e.g., [x]. A list of these different publications ordered according to these reference numbers can be found below in the section entitled “References.” Each of these publications is incorporated by reference herein.)

A common representation of three dimensional objects in computer applications such as graphics and design, are triangular meshes. In many instances, individual or groups of triangles need to satisfy spatial constraints that are imposed either by observation from the real world, or by concrete design specifications for the objects. As these problems tend to be of large scale, choosing a mathematical optimization approach can be particularly challenging. To better understand the problems of the prior art, an explanation of the design problems may be useful.

One issue of the prior art relates to how to solve design problems with concrete applications in the computational surface generation of triangular meshes. In many cases, these surfaces need to satisfy certain spatial constraints.

Digital Terrain Modeling in Computational Geometry

Digital terrain modeling may be done through a set of contours, or a collection of ordered elevation points at arbitrary datums in a landscape. These points can then be connected through a triangular mesh, which is usually called a Triangular Irregular Network (TIN) [22]. As pointed out in [5], traditional triangulation methods may result in terrain models that do not represent the true conditions in terms of drainage pools and channels, and may result in local minima and sharp slope changes that do not exist in reality. In [6], the authors apply heuristics to generate higher order Delaunay triangulations that attempt to mitigate these problems. A more deterministic approach would be to identify existing drainage channels and, within certain boundaries, adjust the triangle surfaces in order to smooth out grade change and create drainage channels.

Computer-Aided Design of Golf Course Features

Contour changes are an integral part of any golf course design. They may involve smaller terrain changes such as cut and fill of bunker walls for erosion control [13], or higher earthwork quantity features such as artificially embanked mounds, as well as differently sloped greens that require cut and fill in order to fit into the existing ground [13]. Tee flattening with constrained end slopes may require large fill quantities, or can be designed by balancing quantities through benched fills [15]. Hence, for each feature of a golf course, the architect needs to consider a variety of surface maximum and minimum slope constraints, surface alignment constraints, as well as other surface slope requirements that are not only part of the course difficulty, but also essential for drainage and environmental impact.

Mesh Conversion

A common problem in topology is the conversion of one surface representation into another. One such example is converting an existing triangular, irregular mesh into the closest regular rectangular mesh. In each context, the term closest could mean different things such as the vertical distance of the corner points of the rectangle to the triangular grid, or the minimum distance of the average of the triangle elevations to the average of the corner points, or the minimum of the volume between the triangles and the rectangular surface that encompasses them. Methods to convert two dimensional, triangular meshes into quadrilateral meshes while keeping some or all of the original vertices are shown in [17, 18, 24]. An extension of these methods to three dimensions would require that the triangles that are grouped into quadrilaterals to satisfy some surface alignment conditions. Furthermore, one would like to minimize the distance between the original and final elevations of the vertices, using either the 1-norm or the 2-norm. FIG. 1 shows an example of an existing triangular, irregular mesh at the bottom, and a quadrilateral grid on the top. The designer would like to convert the three dimensional triangular mesh 102 into a three dimensional quadrilateral mesh that aligns with the grid lines 104 shown on the top. FIG. 2 shows how the existing triangular mesh 102 is mapped in the quadrilateral grid 104. If a quadrilateral cuts an existing triangle, it may need to be retriangulated to form new triangles so that each quadrilateral is triangulated. The mentioned surface alignment conditions would then apply to triangles that are inside an individual quadrilateral. The design objective function is then to minimize the elevation difference between the triangle vertices in the original triangular mesh 102, and the triangle vertices in the new mesh that groups the triangles into quadrilaterals.

Computer-Aided Design for Architecture and Civil Engineering Structures

In Computer-Aided Design (CAD), triangular meshes are widely used in various engineering disciplines. Existing and finished ground surfaces in architecture and civil engineering are represented by triangulated irregular networks. Slopes are relevant in the context of drainage, in both, civil engineering (transportation structures), and architecture (roof designs), as well as in the context of surface alignments such as solar farms, embankment dams, and airport infrastructure layouts.

In architecture, the design of a complex roof requires slope constraints for drainage, but can also include orientation and maximum slope constraints. Consider the example in FIG. 3 of a roof seen from above. The roof polygons 302 with the light dashed lines are triangulated, in order to create a triangular mesh. The engineer may want to place multiple solar panels 304 on a roof. The solar panels 304 require a minimum inclination towards the position of the sun at noon, shown by the arrows 306. Also, the panel 304 should not exceed a certain angle, to preserve a maximum average efficiency. This requires a surface orientation constraint for some of the triangles.

A concrete problem that arises in civil engineering design is the grading of a parking lot. Within a given area, the engineer has to define grading slopes such that the parking lot fits within existing structures, the drainage requirements on the lot are met, and safety and comfort is taken into account. Besides these requirements, the engineer would like to change the existing surface as little as possible, in order to save on earthwork costs.

FIG. 4 depicts an example of a parking lot that has four horizontal parking rows. In FIG. 4, the designer wants water to flow in the direction of the arrows 402 towards and along the thick lines 404 into the inlets 406 that are represented as circles. Lines 408 are curbs that need to be aligned.

Drainage schemes can vary for each situation. The triangular mesh in FIG. 5 represents the existing ground of a planned parking lot. The engineer would like the water to drain away from the building, and along the drainage lines 502 into the four corners. Lines 504 indicate where the triangle edges need to be aligned. In FIG. 6, the engineer would like to study a different scheme, where the drainage 602 happens in parallel, on either side of the building. Lastly, FIG. 7 represents the mesh of an existing ground for a roundabout, where a minimum slope is required from the inner circle to the outer one, and water needs to drain from the road at the top along the outer circle to the roads on either side at the bottom (indicated by drainage lines 702).

SUMMARY OF THE INVENTION

Embodiments of the invention overcome the problems of the prior art based on geometric design constraints for triangles. In particular, a design constraint is defined that includes a maximum slope constraint for at least one of the triangles in a triangular surface mesh. The maximum slope constraint is a maximum angle between a normal vector of the triangle and a reference (e.g., Z) vector. Heights of the vertices of the triangle are projected onto design constraint sets such that the normal vector satisfies all of the design constraints (e.g., the maximum slope constraint). During the projection, the height(s) (e.g., the Z-coordinate [and not the x or y coordinate(s)]) of the triangle is moved by a minimum Euclidian distance. The projection is performed iteratively such that the triangular mesh converges to an optimal solution for a CAD design of the land surface. Accordingly, a design of the surface represented by the triangular surface mesh is generated based on the projecting.

BRIEF DESCRIPTION OF THE DRAWINGS

Referring now to the drawings in which like reference numbers represent corresponding parts throughout:

FIG. 1 shows an example of an existing triangular, irregular mesh at the bottom, and a quadrilateral grid on the top;

FIG. 2 shows how an existing triangular mesh is mapped in a quadrilateral grid;

FIG. 3 illustrates an exemplary roof structure view from above;

FIG. 4 depicts an example of a parking lot that has four horizontal parking rows;

FIG. 5 illustrates a triangle mesh representative of an existing ground of a planned parking lot;

FIG. 6 illustrates a triangular mesh for a planned parking lot where drainage happens in parallel, on either side of the building;

FIG. 7 represents the mesh of an existing ground for a roundabout, where water needs to drain from the road at the top along the outer circle to the roads on either side at the bottom;

FIG. 8 shows a triangular mesh in accordance with one or more embodiments of the invention;

FIG. 9 illustrates a simple triangular mesh with vertices, edges, and two triangles in accordance with one or more embodiments of the invention;

FIG. 10 illustrates two lines A and B that are shown as examples of two constraint sets in accordance with one or more embodiments of the invention;

FIG. 11 shows the iterates of the DR method on two constraint sets A and B in accordance with one or more embodiments of the invention;

FIG. 12 shows a triangular mesh from a birds-eye view in accordance with one or more embodiments of the invention;

FIG. 13 shows the birds-eye view of a feature line on a triangular mesh in accordance with one or more embodiments of the invention;

FIG. 14 illustrates a least squares problem in accordance with one or more embodiments of the invention;

FIG. 15 illustrates a maximum angle for a normal vector of a ground plane with respect to gravity in accordance with one or more embodiments of the invention;

FIG. 16 illustrates an optimal solution based on tangent points between a circle and level sets in accordance with one or more embodiments of the invention;

FIG. 17 illustrates a desired orientation ({right arrow over (q)}) of a normal vector ({right arrow over (n)}) for a triangle plane in accordance with one or more embodiments of the invention;

FIG. 18 illustrates a visualization of a surface minimum slope constraint in accordance with one or more embodiments of the invention;

FIG. 19 illustrates a feasible set region for a normal vector with respect to a Surface Oriented Minimum Slope Constraint in accordance with one or more embodiments of the invention;

FIG. 20 shows that there are exactly two tangent points (u, v) with positive Lagrange multipliers λ, one of which has u<0 and the other has u>0 in accordance with one or more embodiments of the invention;

FIG. 21 shows that there is exactly one tangent point (u, v) with a positive Lagrange multiplier λ and u>0, and at least one tangent point with a nonpositive Lagrange multiplier in accordance with one or more embodiments of the invention;

FIG. 22 illustrates a circular segment defining a feasible set for a maximum slope and an oriented minimum slope in accordance with one or more embodiments of the invention;

FIG. 23 illustrates a circle representing a feasible region for a maximum slope constraint and for two oriented minimum slope constraints in accordance with one or more embodiments of the invention;

FIG. 24 illustrates multiple surface oriented minimum slope constraints via the intersection of two corresponding circular segments in accordance with one or more embodiments of the invention;

FIG. 25 illustrates two adjacent triangles with surface normal vectors that are used to minimize the curvature in accordance with one or more embodiments of the invention;

FIG. 26 illustrates a hexagon that is used to replace a unit disk for a maximum slope constraint in accordance with one or more embodiments of the invention;

FIG. 27 shows how the original half cone is now replaced with the absolute value function in accordance with one or more embodiments of the invention;

FIG. 28 shows the starting situation for the problem shown in FIG. 7 in accordance with one or more embodiments of the invention;

FIG. 29 shows the mesh of FIG. 28 after 10 iterations in accordance with one or more embodiments of the invention;

FIG. 30 shows the mesh of FIG. 28 after 29 iterations in accordance with one or more embodiments of the invention;

FIG. 31 shows the final mesh of FIG. 28 after convergence in accordance with one or more embodiments of the invention;

FIG. 32 illustrates the mesh of FIG. 28 based on the problem shown in FIG. 8 after 15 iterations in accordance with one or more embodiments of the invention;

FIG. 33 illustrates the mesh of FIG. 28 based on the problem shown in FIG. 8 after 100 iterations;

FIG. 34 illustrates the final mesh of FIG. 28 based on the problem shown in FIG. 8 in accordance with one or more embodiments of the invention;

FIG. 35 is the starting situation mesh for the problem shown in FIG. 9 in accordance with one or more embodiments of the invention;

FIG. 36 illustrates the mesh for the problem shown in FIG. 9 after 10 iterations in accordance with one or more embodiments of the invention;

FIG. 37 illustrates the mesh for the problem shown in FIG. 9 after 20 iterations in accordance with one or more embodiments of the invention;

FIG. 38 is the final solution mesh for the problem shown in FIG. 9 after 20 iterations in accordance with one or more embodiments of the invention;

FIG. 39 illustrates the logical flow for designing a surface in accordance with one or more embodiments of the invention;

FIG. 40 is an exemplary hardware and software environment used to implement one or more embodiments of the invention; and

FIG. 41 schematically illustrates a typical distributed/cloud-based computer system using a network to connect client computers to server computers in accordance with one or more embodiments of the invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

In the following description, reference is made to the accompanying drawings which form a part hereof, and which is shown, by way of illustration, several embodiments of the present invention. It is understood that other embodiments may be utilized and structural changes may be made without departing from the scope of the present invention.

The figures herein are used to illustrate various geometric constraints on the triangles, design applications, and the way how iterative projections are performed.

Mathematical notation is utilized herein for a precise description of the geometric constraints and the corresponding projections. The mathematical notation is fairly standard and follows largely [2]. R denotes the set of real numbers. By “x:=y”, or equivalently “y=: x”, we mean that “x is defined by y”. The assignment operators are denoted by “δ” and “→”. The angle between two vectors {right arrow over (n)} and {right arrow over (q)} is denoted by ∠({right arrow over (n)},{right arrow over (q)}). The cross product of {right arrow over (n)}₁ and {right arrow over (n)}₂ in R³ is denoted by {right arrow over (n)}₁×{right arrow over (n)}₂.

General Overview

Embodiments of the invention provide a method and system for utilizing projections to various geometric design constraints for triangles in three dimensions, and proximity operators for design objective functions. The projections and proximity operators are executed on triangles that are part of an existing triangular mesh. Further, these projections and proximity operators are executed in an iterative way, such that the triangular mesh converges to an optimal solution for a computer-aided design in three dimensions.

Problem Formulation

A common representation of three dimensional objects in computer applications such as graphics and design, are triangular meshes. In many instances, individual or groups of triangles need to satisfy spatial constraints that are imposed either by observation from the real world, or by concrete design specifications for the objects. As these problems tend to be of large scale, choosing a mathematical optimization approach can be particularly challenging. For embodiments of the invention, various geometric constraints are defined as convex sets in Euclidean spaces, and the corresponding projections are found. Furthermore, proximity operators are introduced for certain objective functions. The projections and the proximity operators are used in an iterative method to find optimal design solutions.

Triangular Mesh

In computer-aided design, a triangular mesh is a surface that is represented by a set of connected triangles in three dimensional space R³. This form of mesh is also called a Triangular Irregular Network (TIN). FIG. 8 shows a triangular mesh in R³ with underlying contours on the R² plane in accordance with one or more embodiments of the invention. What follows is a precise mathematical definition of such a triangular mesh.

Assume that G=(V₀,E₀) is a given triangular mesh on R² where V₀ is the set of vertices and E₀ is the set of (undirected) edges, i.e.,V₀:={p_i=(p_i1,p_i2)∈R²|i∈{1,2, . . . ,n}},  (1)E₀⊂{p_ip_j≡p_jp_i|p_i,p_j∈V₀}.  (2)From G, one can derive the set of triangles of the mesh as followsT₀:={Δ=(p_ip_jp_k)p|{p_ip_j,p_jp_k,p_kp_i}⊂E₀}.  (3)

One first aims tofind a vector z=(z₁, . . . ,z_n)∈X  (4)such that the points{P_i=(p_i1,p_i2,z_i)}_{i∈{1, . . . ,n}}satisfy a given set of constraints.  (5)

Clearly, the points (P_i)_{i∈{1, . . . ,n}} also form a corresponding triangular mesh S in three dimensions, where the mesh does not contain triangles that have a vertical surface (since we first defined all the vertices as unique points in R²). FIG. 9 illustrates a simple triangular mesh with verticesV₀:={P₁,P₂,P₃,P₄},edgesE₀:={P₁P₄,P₁P₂,P₂P₄,P₃P₄,P₂P₄},and the two trianglesT₀:={(P₁P₂P₄),(P₂P₃P₄)}.in accordance with one or more embodiments of the invention.

Therefore, if there is no confusion, E₀, V₀, T₀ will be reused to denote the vertices, the edges, and the triangles of the three dimensional mesh.

Interval Constraint

A designer may require that some of the vertices in the triangular mesh cannot exceed a certain elevation, or cannot be lower than a certain elevation. Hence, the elevation of these vertices needs to lie within a given interval. This is referred to as the interval constraint and is defined as follows. For a given subset I of {1, 2, . . . , n}, for all i∈I, the entries z_i must lie in a given interval of R.

Edge Slope Constraint

A designer may require that the slope of a triangle edge cannot be steeper than a given value, or must be steeper than a given value. Hence, the edge slope needs to lie within a given interval. This is referred to as the edge slope constraint and is defined as follows. For a given subset E of the edges E₀, and for every edge p_ip_j∈E, the slope

$\begin{array}{cc}{s}_{\mathrm{ij}}:=\frac{z_i-z_j}{\sqrt{{\left(p_{i1}-p_{j1}\right)}^2+{\left(p_{i2}-p_{j2}\right)}^2}}& \left(6\right)\end{array}$ must lie in a given subset of R.

Edge Alignment Constraint

A designer may require that the edges of some of the triangles in the mesh align themselves. This is referred to as the edge alignment constraint and is defined as follows. For a given pair of edges p_ip_j, p_mp_n∈E₀, the edges must have the same slope s_ij=s_mn.

Surface Alignment Constraint

A designer may require that the surfaces of some of the triangles in the mesh align themselves. Aligning triangle surfaces can result in useful patches inside the mesh that may be non-triangular. This is referred to as the surface alignment constraint and is defined as follows. For a given pair of triangles Δ_i, Δ_j∈T₀ the normal vectors {right arrow over (n)}_Δi and {right arrow over (n)}_Δj must be parallel.

Surface Slope Constraint

A designer may require that one or multiple triangles cannot exceed a certain slope with respect to gravity. This is referred to as a maximum surface slope constraint. A designer may also require that some triangles need to have a certain minimum slope in the direction of a given vector. This is referred to as an oriented minimum surface slope constraint. Both surface slope constraints can be defined as a more general surface slope constraints that follow. For a given subset T of the triangles T₀ and for each triangle Δ∈T, there is a given set of vectors Q_Δ ⊂R³ such that for each {right arrow over (q)}∈Q_Δ, the angle between the normal vector {right arrow over (n)}_Δ and {right arrow over (q)} must lie in a given subset θ_{{right arrow over (q)}} of [0, π], which is also referred to as a surface orientation constraint.

Curvature Minimization

In some design problems, the designer may wish to construct a surface that is as “smooth” as possible. This problem is referred to as minimizing the curvature between adjacent triangles in the mesh. One way to minimize the curvature is to minimize the grade change between each pair of adjacent triangles in the mesh. One may define this as follows. For a given pair of triangles Δ_i, Δ_j∈T₀ the normal vectors {right arrow over (n)}_Δi and {right arrow over (n)}_Δj we want to minimize the angle between the two normal vectors.

Methods Overview

Embodiments of the invention are based on projections, and the use of projections in iterative ways to find solutions for feasibility or optimization problems. A projection of a point onto a convex set is the shortest distance of that point to the set. For example, if a set is a plane in three dimensions, then the projection of a point onto that set is its shadow on the plane. If the point to be projected is already on the plane, than the projection is the point itself.

Given multiple sets (e.g. multiple planes in three dimensions) and a starting point, iterative projection methods can be used to find a point in the intersection of all these sets that is closest from the starting point. These methods project iteratively onto the different sets. If the given sets are convex, then the iterative projections will converge to a solution.

In embodiments of the invention, the sets for each constraint are identified. Closed-form projections are then developed for each of the constraint sets. A proximity operator can be found for the curvature minimization. Iterative projection methods, such as the Douglas-Rachford splitting method, are then used to project iteratively to the constraint sets. This alters the elevations of the triangles in the triangular mesh until the mesh converges to a configuration that is desired by the engineer.

Projections onto Constraint Sets

Suppose there are J constraints imposed on a model. For j ∈{1, 2, . . . , J}, we denote by C_j the set of all points that satisfies the j-th constraint. Thus, (4) can be written in the mathematical form

$\begin{array}{cc}\mathrm{find}a\mathrm{point}x\in C:=\bigcap _{j\in \left\{1,2,\dots ,J\right\}}{C}_j,& \left(7\right)\end{array}$ which is referred to herein as a feasibility problem. Moreover, of the infinitude of possible solutions for (7), one may be particularly interested in those that are optimal in some sense. For instance, it could be desirable to find a solution that minimizes the slope change between adjacent triangles, a solution that minimizes the volume between the initial triangles and the triangles in the final solution, or variants and combinations thereof.

If more than one objective function is of interest, it is common to additively combine these functions, perhaps by scaling the functions based on their different levels of importance. In summary, one goal is to solve the problem

$\begin{array}{cc}\min F\left(z\right)\mathrm{subject}\mathrm{to}z\in C:=\bigcap _{j\in \left\{1,2,\dots ,J\right\}}{C}_j,& \left(8\right)\end{array}$ where F may be a sum of (scaled) objective functions. The function F is typically nonsmooth which prevents the use of standard optimization methods.

A constraint set C⊂X is the collection of all feasible data points, i.e., points that satisfy some requirements. Suppose the given data point z∈X is not feasible (i.e., z∉C), one aims to modify z so that the newly obtained point x is feasible (i.e., x∈C ); and one would like to do it with minimal adjustment on z. This task can be achieved by using the projection onto C. Recall that the projection of z onto C, denoted by P_Cz, is the solution of the optimization problem

$\begin{array}{cc}{P}_{C}z:=\underset{x\in C}{\arg \min }z-x=\{x\in C,x-z=\underset{y\in C}{\min }z-y\}.& \left(9\right)\end{array}$

It is well known that if C is nonempty, closed, and convex¹, then P_Cz is singleton, see, e.g., [26, Theorem 2.10]. ¹C is convex if for all x, y∈C and λ∈[0,1], we have (1−λ)x+λy∈C

It turns out that all spatial constraints encountered in various exemplary settings may be convex and closed. Hence, their projections always exist and are unique. Moreover, one can also successfully obtain explicit formulas for these projections.

Iterative Projection Methods for Feasibility Problems

Consider the feasibility problem

$\begin{array}{cc}\mathrm{find}a\mathrm{point}x\in C:=\bigcap _{j\in \left\{1,2,\dots ,J\right\}}{C}_j,& \left(10\right)\end{array}$

The intersection C is the feasible set containing the desired solutions. Clearly, it might be tempting to solve this problem by finding the projection onto C directly. However, this is often not possible due to the complexity of C. A workaround is to utilize the projection P_Cj onto each constraint, if the explicit formula is available. Then with an initial point, one can iteratively execute the projections P_Cj's to derive a solution for (10). Let z₀ ∈X be the initial data point. The following are just a few of many (iterative) projection methods (see [2, 10, 12] and the references therein).

Cyclic Projections

In this method, a point is projected in sequence, from one constraint set to the next. Once the projected point reaches the first constraint set again, it is referred to as one cycle or one iteration. FIG. 10 illustrates two lines A and B that are shown as examples of two constraint sets. The method starts with point b₋₁ that is projected to the set A, resulting in the new point a₀ that then gets projected to the set B. The obtained point b₀ completes one cycle. The process is repeated until the iterations converge to the intersection 1002 of the two sets A and B. The cyclic projections are defiend as follows. Given z_k, one updates z_k+1:=Tz_k whereT:=P_CJP_CJ-1 . . . P_C2P_C1.  (11)

Parallel Projections

In this method, a point is projected in parallel to all the constraint sets at once. The resulting points are then averaged to obtain a new point for the next iteration. One defines the parallel projections as follows. Given z_k, one updates z_k+1:=Tz_k where

$\begin{array}{cc}T:=\frac{1}{J}\left(P_{C_1}+P_{C_2}+\dots +P_{C_J}\right).& \left(12\right)\end{array}$

String-Average Projections

The string-average projections is a combination of cyclic and parallel projections. One defines it as follows. Given z_k, one updates z_k+1:=Tz_k where

$\begin{array}{cc}T:=\frac{1}{J}\left(P_{C_1}+P_{C_2}{P}_{C_1}+\dots +P_{C_J}{P}_{C_{J-1}}\dots P_{C_2}{P}_{C_1}\right).& \left(13\right)\end{array}$

In fact, these are probably the simplest methods for solving feasibility problems. In addition, when projection methods succeed (see [11] for some interesting examples), they have various attractive features: they are easy to understand, simple for implementation and maintenance, and sometime can be very fast. [1, 3, 10, 12] provide additional discussions on projection methods.

Iterative Projection Methods for Optimization Problems

Iterative optimization methods are often used for solving (8), which may require the computations of proximity operators for the functions involved. Suppose ƒ:X→(−∞,+∞] is a proper convex lower semicontinuous function (see [25] and [2] for a discussion of convex analysis) and fix x∈X. Then it is well known (see [2, Section 12.4]) that the function

$\begin{array}{cc}X\to \left(-\infty ,+\infty \right]\text{:}y\mapsto f\left(y\right)+\frac{1}{2}{x-y}^2& \left(14\right)\end{array}$ has a unique minimizer, which will be denoted by P_ƒ (x). The induced operator P_ƒ: X→X is called the proximal mapping or proximity operator of ƒ (see [23]). Note that if ƒ is the indicator function of a set C (the indicator function t_C is defined by t_C(x)=0 if x∈C and t_C(x)=+∞ otherwise), then P_ƒ=P_C. Thus, proximity operators are generalizations of projections.

The Douglas-Rachford Method

As proximity operators may be made available for several types of objective functions, any iterative optimization methods that utilize proximity operators (e.g., [7, 8, 21]) can be used to solve the corresponding problems. Let us describe one such method, the Douglas-Rachford (DR) splitting. DR emerged from the field of differential equations [14], and was later made widely applicable in other areas thanks to the seminal work [21]. FIG. 11 shows the iterates of the DR method on two constraint sets A and B in accordance with one or more embodiments of the invention.

Using indicator functions, (8) is converted intomin F(x)+t_C1+t_C2+ . . . +t_CJ subject to x∈X.  (15)

So, it suffices to present DR for the following general problem

$\begin{array}{cc}\min \sum _{i=1}^m{f}_i\left(x\right)\mathrm{subject}\mathrm{to}x\in X,& \left(16\right)\end{array}$ where each ƒ_i is a proper convex lower semicontinuous function on X. The DR operates in the product space X:=X^m with inner product x,y:=Σ_i=1^mx_i, y_i for x=(x_i)_i=1^m and y=(y_i)_i=1^m. Set the starting point x₀=(z, . . . , z)∈X, where z∈X. Given x_k=(x_k,l, . . . , x_k,m)∈X, one computes

$\begin{array}{cc}{\stackrel{\_}{x}}_k:=\frac{1}{m}\sum _{i\in K}{x}_{k,i},& \left(17\right)\end{array}$ ∀i∈{1, . . . ,m}:x_k+1,i:=x_k,i−x_k+P_γƒi(2x_k,i−x_k,i),  (18)then update x_k+1:=(x_k+1,l, . . . ,x_k+1,m).  (19)Then the monitored sequence (x_k)_k∈N converges to a solution of (16).

It is worth mentioning that when all ƒ_i's are indicator functions of the constraints (thus, all proximity operators become projections), then (16) becomes a pure feasibility problem. Therefore, all optimization methods for (16) can also be used for feasibility problems.

Parallelization

By using the indicator function, projections become a special case of proximity operators. Thus, the following remarks simply refer to proximity operators.

Remark 1 (Entries that Possibly Change after Executing a Proximity Operator)

Suppose X=Rⁿ and the objective function ƒ: X→(−∞, +∞] only depends on certain entries, i.e., ƒ(x)=ƒ((x_i)_i∈I), where I⊂{1, 2, . . . , n} is a given index set.

Let z=(z₁, . . . , z_n)∈X and execute the proximity operator P_ƒZ. It is obvious that after the execution, the coordinates that possibly change are (z_i)_i∈I; while all other coordinates (z_i)_i∉I are left unchanged, see, e.g., Method component 2 (projection onto an edge minimum slope constraint). Thus, for simplicity of the presentation, for each projection or proximity operator, one may often focus on the “possibly affected” coordinates (z_i)_i∈I in the subspace R^|I| where |I| is the number of indices in I.

Remark 2 (Parallelization)

If two functions ƒ₁, ƒ₂: X→(−∞,+∞] are independent in the sense that they depend on separate coordinates, then one can execute their respective proximity operators in parallel. Therefore, one can pack multiple independent proximity operators to accelerate the convergence of the method.

Projections onto Linear Constraints

The following describes the projections found for the linear constraints that were defined in the problem formulation above. These closed-form projections ensure that the method will converge to a solution in a reasonable amount of time, using an iterative projection method.

By linear constraints, one refers to any constraint on the triangular mesh that can be represented by a system of linear equalities and inequalities. Indeed, this class includes several important constraints in design problems. In the sections that follow, those constraints will be analyzed.

Interval Constraint

Similar to [3, Section 2.2], one may assume that I is a subset of {1, 2, . . . , n} and (l_i)_i∈I, (u_i)_i∈I∈R^I are given. SetY:={z=(z₁,z₂, . . . ,z_n)∈X|∀i∈I:l_i≤z_i≤u_i}.  (20)

The closed set Y is an affine subspace, thus, convex. The following explicit formula is for the projection onto Y, whose proof is straightforward, thus, omitted.

Method Component 1 (Projection onto Interval Constraints)

$\begin{array}{cc}{P}_Y\text{:}X\to X\text{:}\left(z_1,z_2,\dots ,z_n\right)\mapsto \left(x_1,x_2,\dots ,x_n\right),& \left(21\right)\\ \mathrm{where}{x}_i=\{\begin{array}{cc}\max \left\{l_i,\min \left\{u_i,z_i\right\}\right\}& i\in I,\\ z_i,& i\in \left\{1,2,\dots ,n\right\}\backslash I\end{array}.& \left(22\right)\end{array}$

Edge Minimum Slope Constraint

Let P_i=(p_i1, p_i2, h_i) and P_j=(p_j1, p_j2, h_j). The designer may require that the (directional) slope from P_j to P_i must be no less than a threshold level s_ij. More specifically, since the value p_i1, p_i2, p_j1, p_j2 are fixed, one can write this constraint ash_i−h_j≥α_ij:=s_ij√{square root over ((p_i1−p_j1)²+(p_i2−q_j2)²)},  (23)which is called the edge minimum slope constraint. This constraint is a linear inequality, thus, convex. Projection on to this type of constraint is derived analogously to [3, Section 2.3], thus, its proof is also omitted.

Method Component 2 (Projection onto an Edge Minimum Slope Constraint)

LetE_i,j={x=(x₁,x₂, . . . ,x_n)∈Rⁿ|x_i−x_j≥α_ij}  (24)be the minimum slope constraint on the edge P_iP_j. Let x:=(x₁, . . . , x_n):=P_Ei,jz be the projection of a point z=(z₁, z₂, . . . , z_n) on to E_i,j. Then

$\begin{array}{cc}{x}_i=\frac{1}{2}\left(z_i+z_j+\max \left\{z_i-z_j,{\alpha }_{\mathrm{ij}}\right\}\right),& \left(25\right)\\ x_j=\frac{1}{2}\left(z_i+z_j-\max \left\{z_i-z_j,{\alpha }_{\mathrm{ij}}\right\}\right),\mathrm{and}& \left(26\right)\\ x_k=z_k\mathrm{for}\mathrm{all}k\notin \left\{i,j\right\}.& \left(27\right)\end{array}$

Low Point Constraint

A point P_j=(p_j1, p_j2, z_j) on the mesh is called a low point if each point P_i=(p_i1, p_i2, z_i) connected to P_j satisfies the edge minimum slope constraintz_i−z_j≥α_i:=s_i√{square root over ((p_i1−p_j1)²+(p_i2−q_j2)²)},  (28)where all s_i∈R are given.

FIG. 12 shows a triangular mesh from a birds-eye view, where the low point P₀ is connected to points P₁, P₂, P₃, P₄, P₅, and P₆. So in this picture, one requires that P₀ is lower than all its connected points and that the slopes on each connecting edge i is at least s_i.

One can treat the low point constraint as a single constraint even though it is the intersection of finitely many edge slope constraints. The following result shows how to project onto this constraint.

Method Component 3 (Projection onto a Low Point Constraint)

Let α₂, . . . , α_m∈R. Define the setE:={(x₁, . . . ,x_m)∈R^m|∀i∈{2, . . . ,m}:x_i−x₁≥α_i}.  (29)

Let z=(z₁, . . . , z_m)∈R^m. Let δ₁=z₁ and let δ₂≤δ₃≤ . . . ≤δ_m be the values {z_i−α_i}_{i∈{2, . . . , m}} in nondecreasing order. Let k be the largest number in {1, . . . , m} such that δ_k≤(δ₁+ . . . +δ_k)/k. Then the projection x:=(x₁, x₂, . . . , x_m):=P_EZ is given byx₁=(δ₁+ . . . +δ_k)/k and x_i=max{x₁,z_i−α_i}+α_i,∀i∈{2, . . . ,m}.   (30)

To show this formula, first, x is the solution ofmin(x₁−z₁)²+ . . . +(x_m−z_m)²  (31)s.t x₁−x_i+α_i≤0,i∈{2, . . . ,m}.  (32)

Let y:=(y₁, . . . , y_m) where y₁=x₁ and y_i:=x_i−α_i for i∈{2, . . . , m}. Without loss of generality, one may assume δ_i=z_i−α_i for i∈{2, . . . , m}. Then (4) becomesmin ϕ(y):=(y₁−δ₁)²+ . . . +(y_m−δ_m)²  (33)s.t g_i(y):=y₁−y_i≤0,i∈{2, . . . ,m}.  (34)

To find the (unique) solution, we use Lagrange multipliers for convex optimization: there exist λ_i≥0, λ∈{2, . . . , m}, such that

$\begin{array}{cc}\frac{1}{2}\nabla \varphi \left(y\right)+{\lambda }_2\nabla g_2\left(y\right)+\dots +{\lambda }_m\nabla g_m\left(y\right)=0& \left(35\right)\\ \forall i\in \left\{2,\dots ,m\right\}\text{:}{\lambda }_i{g}_i\left(y\right)=0& \left(36\right)\\ \forall i\in \left\{2,\dots ,m\right\}\text{:}{g}_i\left(y\right)\le 0.& \left(37\right)\end{array}$

It follows from (35) that(y₁−δ₁)+λ₂+ . . . +λ_m=0  (38)(y₂−δ₂)−λ₂=0  (39). . .  (40)(y_m−δ_m)−λ_m=0.  (41)

Then y₁≤δ₁ because λ₂, . . . , λ_m≥0. By substitution, one getsy₁+y₂+ . . . +y_m=δ₁+δ₂+ . . . δ_m.  (42)

Next, for each i∈{2, . . . , m}, one has y_i−δ_i=λ_i≥0, y_i−y₁≥0, and (36) implies,0=λ_ig_i(y)=(y_i−δ_i)(y₁−y_i).  (43)This leads to either y_i=δ_i≥y₁ or y₁=y_i≥δ_i. That meansy_i=max{y₁,δ_i} for all i∈{2, . . . ,m}.  (44)

So (42) becomesδ₁+δ₂+ . . . +δ_m=y₁+max{y₁,δ₂}+ . . . +max{y₁,δ_m}.  (45)Next, suppose k is the largest number in {1, . . . , m} such thatδ_k≤(δ₁+ . . . +δ_k)/k.  (46)

The choice of k is always possible since at least (46) is true for k=1. Now setting δ_m+1=+∞, one claims that(δ₁+ . . . +δ_k)/k<δ_k+1.  (47)

In fact, if k=m, then (47) clearly holds. If k<m, then by the choice of k, one has δ₁+ . . . +δ_k+δ_k+1<(k+1)δ_k+1, which implies δ₁+ . . . +δ_k<kδ_k+1. So (47) holds. Next, one can show thatδ_k≤y₁≤δ_k+1.  (48)

Indeed, on the one hand, (46) and (45) imply

$\begin{array}{cc}k{\delta }_k+{\delta }_{k+1}+\dots +{\delta }_m\le \left({\delta }_1+\dots +{\delta }_k\right)+{\delta }_{k+1}+\dots +{\delta }_m& \left(49\right)\\ =y_1+\max \left\{y_1,{\delta }_2\right\}+\dots +\max \left\{y_1,{\delta }_m\right\}& \left(50\right)\\ \le y_1+\underset{\underset{k-1\mathrm{terms}}{︸}}{\max \left\{y_1,{\delta }_k\right\}+\dots +\max \left\{y_1,{\delta }_k\right\}}& \left(51\right)\\ +\max \left\{y_1,{\delta }_{k+1}\right\}+\dots +\max \left\{y_1,{\delta }_m\right\}& \left(52\right)\\ =y_1+\left(k-1\right)\max \left\{y_1,{\delta }_k\right)+\max \left\{y_1,{\delta }_{k+1}\right\}+\dots +\max \left\{y_1,{\delta }_m\right\}.& \left(53\right)\end{array}$ If y₁<δ_k, then (4) implieskδ_k=y₁+(k−1)max{y₁,δ_k}<kδ_k,  (54)which is a contradiction. Thus, y₁≥δ_k. On the other hand, (47) and (45) implies

$\begin{array}{cc}k{\delta }_{k+1}+{\delta }_{k+1}+\dots +{\delta }_m\ge \left({\delta }_1+\dots +{\delta }_k\right)+{\delta }_{k+1}+\dots +{\delta }_m& \left(55\right)\\ =y_1+\max \left\{y_1,{\delta }_2\right\}+\dots +\max \left\{y_1,{\delta }_m\right\}& \left(56\right)\\ =y_1+\underset{\underset{k-1\mathrm{terms}}{︸}}{y_1+\dots +y_1}+\max \left\{y_1,{\delta }_{k+1}\right\}+\dots +\max \left\{y_1,{\delta }_m\right\}& \left(57\right)\\ \ge {\mathrm{ky}}_1+\max \left\{y_1,{\delta }_{k+1}\right\}+{\delta }_{k+2}+\dots +{\delta }_m.& \left(58\right)\end{array}$ If y₁>δ_k+1, then (55)-(58) implies(k+1)δ_k+1≥ky₁+max{y₁,δ_k+1}>(k+1)δ_k+1,  (59)which is a contradiction. Thus, y₁≤δ_k+1. Therefore, the claim (48) is true. Then (45) impliesδ₁+δ₂+ . . . +δ_m=y₁+max{y₁,δ₂}+ . . . +max{y₁,δ_m}=ky₁+δ_k+1+ . . . +δ_m,   (60)which leads to y₁=(δ₁+ . . . +δ_k)/k. Combining with (44) completes the proof.

Edge Alignment Constraint

On the triangular mesh, the designer may want a constant slope on a particular path, in which case one can say the path is “aligned”. Such a path is sometimes called a feature line. To formulate this constraint, suppose the feature line is given by adjacent points A₁, A₂, . . . , A_m on the triangular mesh where A₁=(a_i1, a_i2, x_i). FIG. 13 shows the birds-eye view of a feature line on a triangular mesh in accordance with one or more embodiments of the invention. For A_iA_i+1, the length of its “shadow” on the xy-plane isδ_i:=(a_i+1,1−a_i,1)²+(a_i+1,2−a_i,2)².  (61)Define alsot₁:=0,t₂:=δ₁,t₃:=δ₁+δ₂, . . . ,t_m:=δ₁+ . . . +δ_m−1.  (62)Then the alignment constraint is written as∀i∈{1, . . . ,m−2}:(x_i+1−x_i)/(t_i+1−t_i)=(x_i+2−x_i+1)/(t_i+2−t_i+1).  (63)

Method Component 4 (Projection onto an Edge Alignment Constraint)

Suppose the points (a_i1, a_i2), i∈{1, . . . , m} forms a path in R². Define (t₁, . . . , t_m) by (62) and define the convex setF:={(x₁, . . . ,x_m)∈R^m|((a_i1,a_i2,x_i))_{i∈{i, . . . ,m}}satisfies(63)}⊂R^m.  (64)

Let z=(z₁, . . . , z_m)∈R^m. Then the projection P_Fz is given by∀i ∈{1, . . . ,m}:(P_FZ)_i=ƒ(t_i).  (65)where ƒ(t):=α+βt is the least squares line for the points (t_i, z_i)_{i∈{1, . . . m}}. More specifically, (α, β) is obtained from the linear system

$\begin{array}{cc}\left[\begin{array}{cc}m& \sum _{i=1}^m{t}_i\\ \sum _{i=1}^m{t}_i& \sum _{i=1}^m{t}_i^2\end{array}\right]\left[\begin{array}{c}\alpha \\ \beta \end{array}\right]=\left[\begin{array}{c}\sum _{i=1}^m{z}_i\\ \sum _{i=1}^m{t}_i{z}_i\end{array}\right]& \left(66\right)\end{array}$

To show this, first, F is convex since all constraints in (63) are linear. Next, one considers the points (t_i, z_i) in R². The problem is to find (x₁, . . . , x_m) such that the points (t_i, x_i) are aligned and∥x−z∥²=E_i=1^m(x_i−z_i)² is minimized.  (67)

This is the least squares problem for the points (t_i, z_i), shown in FIG. 14. So the least squares line is ƒ(t)=α+βt where (α, β) satisfies (66). Thus, the projection coordinates are (P_F (z))_i=α+βt_i for i∈{1, . . . , m}.

Surface Alignment Constraint

The designer may want to patch several adjacent triangles on the mesh into a single polygon, in which case one may say that these triangles are “aligned”. This is equivalent to require all vertices of the triangles to lie on the same plane. So the following may result.

Method Component 5 (Projection onto a Surface Alignment Constraint)

Let (a_i1, a_i2), i∈{1, . . . , m} be a collection of points in R² that are not on the same line. DefineF:={(x₁, . . . ,x_m)∈R^m|{(a_i1,a_i2,x_i)}_{i∈{1, . . . ,m}}lie on the same plane}.   (69)Let z=(z₁, . . . , z_m)∈R^m. Then the projection P_Fz is given by∀i∈{1, . . . ,m}:(P_Fz)_i=ƒ(a_i1,a_i2).  (70)where ƒ(t₁, t₂):=α+∈t₁+γt₂ is the least squares plane for the points (a_i1, a_i2, z_i), i∈{1, . . . , m}. More specifically, (α, β, γ) is the solution of the system

$\begin{array}{cc}\left[\begin{array}{ccc}m& \langle e,a_1\rangle & \langle e,a_1\rangle \\ \langle e,a_1\rangle & {a_1}^2& \langle a_1,a_2\rangle \\ \langle e,a_2\rangle & \langle a_1,a_2\rangle & {a_2}^2\end{array}\right]\left[\begin{array}{c}\alpha \\ \beta \\ \gamma \end{array}\right]=\left[\begin{array}{c}\langle e,z\rangle \\ \langle a_1,z\rangle \\ \langle a_2,z\rangle \end{array}\right].& \left(71\right)\end{array}$ where e=(1, 1, . . . , 1), a₁=(a₁₁, a₂₁, . . . , a_m1), and a₂=(a₁₂ a₂₂, . . . , a_m2).

To show this, first, F is a convex set. Let x=(x₁, . . . , x_m)=P_FZ, then x minimizes∥x−z∥²=Σ_i=1^m(x_i−z_i)²  (72)subject to the constraint that {(a_i1, a_i2, x_i)}_{i∈{1, . . . , m}} lie on the same plane in R³. Thus, it is equivalent to find the least squares plane ƒ(t₁, t₂)=α+βt₁+γt₂ that fits {(a_i1, a_i2, z_i)}_{i∈{1, . . . , m}}. One defines e, a₁, a₂ as above, then (α, β, γ) is the solution of (71). Since {(a_i1, a_i2)}_{i∈{1, . . . , m}} are not on the same line, (71) has a unique solution (α, β, γ). Finally, one computes x_i=α+βa_i1+γa_i2 for i∈{1, . . . , m}.Projections onto General Surface Slope Constraints

Surface slope constraints are any requirement imposed on the slope of a triangle. In this section, a general setup is provided for projections onto such constraints.

Normal Vector and Surface Slope Constraint

Let (e₁, e₂, e₃) be the standard basis of R³. Given three points P₁=(P₁₁, P₁₂, h₁), P₂=(P₂₁, p₂₂, h₂), and P₃=(p₃₁, p₃₂, h₃) in R³, the normal vector to the plane P₁P₂P₃ is the cross product

$\begin{array}{cc}\overrightarrow{n}=\left[\begin{array}{c}{p}_{11}-p_{31}\\ p_{12}-p_{32}\\ h_1-h_3\end{array}\right]\times \left[\begin{array}{c}{p}_{21}-p_{31}\\ p_{22}-p_{32}\\ h_2-h_3\end{array}\right]=\text{:}\left[\begin{array}{c}{a}_1\\ b_1\\ t_1\end{array}\right]\times \left[\begin{array}{c}{a}_2\\ b_2\\ t_2\end{array}\right]=\left[\begin{array}{c}{b}_1{t}_2-b_2{t}_1\\ a_2{t}_1-a_1{t}_2\\ a_1{b}_2-a_2{b}_1\end{array}\right]& \left(73\right)\end{array}$ The “shadows” of P₁, P₂, P₃ on xy-plane are (p₁₁, p₁₂), (p₂₁, p₂₂), (p₃₁, p₃₂), which are not on the same line. This implies a₁b₂−a₂b₁≠0. So one can rescale and use

$\begin{array}{cc}\overrightarrow{n}=\left(\frac{b_1{t}_2-b_2{t}_1}{a_1{b}_2-a_2{b}_1},-\frac{a_1{t}_2-a_2{t}_1}{a_1{b}_2-a_2{b}_1},1\right).& \left(74\right)\end{array}$ If one defines

$u=\frac{b_2{t}_1-b_1{t}_2}{a_1{b}_2-a_2{b}_1}\mathrm{and}v=\frac{a_1{t}_2-a_2{t}_1}{a_1{b}_2-a_2{b}_1},$ thent₁=a₁u+b₁v,t₂=a₂u+b₂v, and {right arrow over (n)}=(−u,−v,1).  (75)Obviously, surface slope constraints depend solely on the normal vector {right arrow over (n)}. Thus, one can represent a general surface slope constraint in the formg(u,v)≤0.  (76)An important case of such constraints is the surface orientation constraint below.

Definition 1 (surface orientation constraint) Let Δ be a triangle with normal vector {right arrow over (n)}=(−u,−v,1)∈R³ as above. Let {right arrow over (q)}=(q₁, q₂, q₃)∈R³ \{0} be a unit vector and θ be an angle in [0, π], the constraint

$\begin{array}{cc}\angle \left(\overrightarrow{n},\overrightarrow{q}\right)\le \theta ,\mathrm{orequivalently},\cos \angle \left(\overrightarrow{n},\overrightarrow{q}\right)=\frac{-{\mathrm{uq}}_1-{\mathrm{vq}}_2+q_3}{\sqrt{u^2+v^2+1}}\ge \cos \theta ,& \left(77\right)\end{array}$ is called the surface orientation constraint specified by the pair ({right arrow over (q)}, θ). Below is a description of the general framework for projection onto general surface slope constraints, followed by considerations of two special surface orientation constraints: the surface maximum slope constraint and the surface oriented minimum slope constraint.

Projection onto a General Surface Slope Constraint

Consider a single triangle determined by three points Q₁=(p₁₁, p₁₂, w₁), Q₂=(p₂₁, p₂₂, w₂), and Q₃=(p₃₁, p₃₂, w₃) in R³. Without loss of generality, one may assumew₁+w₂+w₃=0.  (78)Projecting onto an orientation constraint is to find the new heights h₁, h₂, h₃ that minimizesϕ:=(h₁−w₁)²+(h₂−w₂)²+(h₃−w₃)²  (79)such that the new points P₁=(p₁₁, p₁₂, h₁), P₂=(p₁₁, p₁₂, h₂), and P₃=(p₁₁, p₁₂, h₃) satisfy the given orientation constraint. Define t₁:=h₁−h₃ and t₂:=h₂−h₃, thenϕ=(t₁+h₃−w₁)²+(t₂+h₃−w₂)²+(h₃−w₃)².  (80)Using the change of variables (75), one has:ϕ=(a₁u+b₁v+h₃−w₁)²+(a₂u+b₂v+h₃−w₂)²+(h₃−w₃)².  (81)The new triangle P₁P₂P₃ needs to satisfy some given orientation constraint g(u, v)≤0. So the projection reduces to

$\begin{array}{cc}\underset{\left(u,v,h_3\right)\in R^3}{\min }\varphi ={\left(a_{1}u+b_{1}v+h_3-w_1\right)}^2+{\left(a_{2}u+b_{2}v+h_3-w_2\right)}^2+{\left(h_3-w_3\right)}^2& \left(82\right)\\ \mathrm{subject}\mathrm{to}g\left(u,v\right)\le 0.& \left(83\right)\end{array}$

Next, suppose that changing variables is necessary, for instance, (u, v) is replaced by the new variables (û, {circumflex over (v)}) under the linear transformation

$\begin{array}{cc}\left[\begin{array}{c}\hat{u}\\ \hat{v}\end{array}\right]:=Q\left[\begin{array}{c}u\\ v\end{array}\right],\mathrm{where}Q\mathrm{is}\mathrm{an}\mathrm{invertible}\mathrm{matrix}.& \left(84\right)\end{array}$

One may also define

$\left[\begin{array}{cc}{\hat{a}}_1& {\hat{b}}_1\\ {\hat{a}}_2& {\hat{b}}_2\end{array}\right]:=\left[\begin{array}{cc}{a}_1& b_1\\ a_2& b_2\end{array}\right]Q^{-1}.$ Then

$\begin{array}{cc}{\hat{a}}_1\hat{u}+{\hat{b}}_1\hat{v}=\left[\begin{array}{cc}{\hat{a}}_1& {\hat{b}}_1\end{array}\right]\left[\begin{array}{c}\hat{u}\\ \hat{v}\end{array}\right]=\left[\begin{array}{cc}{a}_1& b_1\end{array}\right]Q^{-1}Q\left[\begin{array}{c}u\\ v\end{array}\right]=a_{1}u+b_{1}v,\mathrm{and}& \left(85\right)\\ {\hat{a}}_2\hat{u}+{\hat{b}}_2\hat{v}=a_{2}u+b_{2}v.& \left(86\right)\end{array}$ So (82)-(83) becomes

$\begin{array}{cc}\underset{\left(\hat{u},\hat{v},h_3\right)\in R^3}{\min }\varphi ={\left({\hat{a}}_1\hat{u}+{\hat{b}}_1\hat{v}+h_3-w_1\right)}^2+{\left({\hat{a}}_2\hat{u}+{\hat{b}}_2\hat{v}+h_3-w_2\right)}^2+{\left(h_3-w_3\right)}^2& \left(87\right)\\ \mathrm{subject}\mathrm{to}\hat{g}\left(\hat{u},\hat{v}\right):=g\left(u,v\right)\le 0,& \left(88\right)\end{array}$ Therefore, one can treat (87)-(88) as (82)-(83) with new respective coefficients â_i, {circumflex over (b)}_i

Next, one can simplify the model (82)-(83). Since (83) does not involve h₃, one can convert problem (82)-(83) into two variables (u, v) as follows: first, set the derivative of ϕ with respect to h₃ to zero∇_h3ϕ(u,v,h₃)=2(a₁u+b₁v+h₃−w₁)+2(a₂u+b₂v+h₃−w₂)+2(h₃−w₃)=0.(89)Then one obtains

$\begin{array}{cc}{h}_3=-\frac{1}{3}\left[\left(a_1+a_2\right)u+\left(b_1+b_2\right)v\right].& \left(90\right)\end{array}$ Substituting into (82) and multiplying by 9, one obtains

$\begin{array}{cc}9\varphi ={\left[\left(2a_1-a_2\right)u+\left(2b_1-b_2\right)v-3w_1\right]}^2& \left(91\right)\\ +{\left[\left(2a_2-a_1\right)u+\left(2b_2-b_1\right)v-3w_2\right]}^2& \left(92\right)\\ +{\left[\left(a_1+a_2\right)u+\left(b_1+b_2\right)v-3\left(w_1+w_2\right)\right]}^2& \left(93\right)\\ =\text{:}\left(6a_1^2+6a_2^2-6a_1{a}_2\right)u^2+\left(6b_1^2+6b_2^2-6b_1{b}_2\right)v^2& \left(94\right)\\ +6\left(2a_1{b}_1+2a_2{b}_2-a_1{b}_2-a_2{b}_1\right)\mathrm{uv}& \left(95\right)\\ -18\left(a_1{w}_1+a_2{w}_2\right)u-18\left(b_1{w}_1+b_2{w}_2\right)v+Z.& \left(96\right)\end{array}$ where Z is some constant. So, by definingA=2(a₁²+a₂²−a₁a₂)>0,  (97)B=2(b₁²+b₂²−b₁b₂)>0,  (98)C=2a₁b₁+2a₂b₂−a₁b₂−a₂b₁,  (99)w_a=3(a₁w₁+a₂w₂),  (100)w_b=3(b₁w₁+b₂w₂),  (101)one has 9ϕ=3Au²+3Bv²+6Cuv−6w_au−6w_bv+Z. Thus, (82)-(83) reduces to

$\begin{array}{cc}\underset{\left(u,v\right)\in R^2}{\min }\varphi \left(u,v\right)=\frac{A}{2}{u}^2+\mathrm{Cuv}+\frac{B}{2}{v}^2-w_{a}u-w_{b}v& \left(102\right)\\ \mathrm{subjectto}g\left(u,v\right)\le 0.& \left(103\right)\end{array}$ As long as one can find the solution (u, v), one can obtain the projection (h₁, h₂, h₃) from (90), (75), and (73). More specifically,

$\begin{array}{cc}{h}_3=-\frac{1}{3}\left[\left(a_1+a_2\right)u+\left(b_1+b_2\right)v\right],& \left(104\right)\\ h_1=a_{1}u+b_{1}v+h_3,& \left(105\right)\\ h_2=a_{2}u+b_{2}v+h_3.& \left(106\right)\end{array}$ Due to (85)-(86), if new variables (û, {circumflex over (v)}) are used, then one also has (104)-(106) in which (a_i, b_i, u, v) is replaced by (â_i, {circumflex over (b)}_i, û, {circumflex over (v)}).Projection onto Surface Maximum Slope Constraint

Let P₁P₂P₃ be a triangle in R³ with normal vector {right arrow over (n)}=(−u,−v,1) obtained as in the normal vector and surface slope constraint section above. In certain cases, it is required that the angle between {right arrow over (n)} and a given direction {right arrow over (q)} must not exceed a given threshold. For example, suppose P₁P₂P₃ represents the desired ground, that cannot be too steep with respect to gravity, i.e., the slope of the plane P₁P₂P₃ must not exceed a threshold s:=s_max∈R₊. Then the angle between {circumflex over (n)} and e₃=(0,0,1) must satisfy ∠({circumflex over (n)},e₃)≤tan⁻¹ (s), which is shown as cone 1502 in FIG. 15. Accordingly, FIG. 15 illustrates a maximum angle for a normal vector of a ground plane with respect to gravity in accordance with one or more embodiments of the invention. Thus, whenever the normal vector {right arrow over (n)} is outside the cone 1502, the constraint is violated. The inequality ∠({right arrow over (n)}, e₃)≤tan⁻¹ (s) is equivalent tocos ∠({right arrow over (n)},e₃)={right arrow over (n)},e₃/∥{right arrow over (n)}∥≥cos(tan⁻¹(s))=1/√{square root over (1+s²)}.  (107)

Definition 2 (surface maximum slope constraint) One calls (107) the surface maximum slope constraint with maximum slope s. This is a special case of surface orientation constraint where ({right arrow over (q)}, θ)=(e₃, tan⁻¹ (s)). From (107), one may conclude

$\begin{array}{cc}\frac{1}{\sqrt{u^2+v^2+1}}\ge \frac{1}{\sqrt{1+s^2}}\iff u^2+v^2-s^2\le 0.& \left(108\right)\end{array}$ Using the general model (102)-(103) with the constants defined by (97)-(101), projection onto the surface maximum slope constraint reduces to

$\begin{array}{cc}\underset{\left(u,v\right)\in R^2}{\min }\varphi =\frac{A}{2}{u}^2+\mathrm{Cuv}+\frac{B}{2}{v}^2-w_{a}u-w_{b}v& \left(109\right)\\ \mathrm{subject}\mathrm{to}{u}^2+v^2-s^2\le 0,& \left(110\right)\end{array}$ Rescaling by

$\left(u,v,w_a,w_b\right)\leftarrow \left(\frac{u}{s},\frac{v}{s},\frac{w_a}{s},\frac{w_b}{s}\right),$ one can obtain an equivalent problem

$\begin{array}{cc}\underset{\left(u,v\right)\in R^2}{\min }\varphi \left(u,v\right)=\frac{A}{2}{u}^2+\mathrm{Cuv}+\frac{B}{2}{v}^2-w_{a}u-w_{b}u& \left(111\right)\\ \mathrm{subject}\mathrm{to}g\left(u,v\right):=u^2+v^2-1\le 0.& \left(112\right)\end{array}$

Hence, the projection to the surface maximum slope constraint will find new elevations for the triangle vertices, such that the angle between the z-axis and the normal vector of the triangle surface does not exceed tan⁻¹ (s) (i.e. the normal vector lays in the cone 1502 depicted in FIG. 15), and that equation (79) is minimized.

This problem may be solved directly by several ways including numerical methods. One exemplary method is presented herein, in which the well-known Lagrange multipliers, also known as Karush-Kuhn-Tucker (KKT) conditions [19, 20], and Ferrari's method for quartic equations [16] are used. First, one observes the following:

1. ϕ(u, v) is a strictly convex quadratic surface.

2. For each value η≥0, the level set ϕ(u, v)=η is an ellipse in R² whose center is the vertex (u₀, v₀) of ϕ(u, v), which satisfiesAu₀+Cv₀=w_a  (113)Cu₀+Bv₀=w_b.  (114)

3. This is minimizing a strictly convex function over a closed, bounded, convex set. Thus, there exists a unique solution.

One may consider two cases:

Case 1: g(u₀, v₀)≤0. Then (u₀, v₀) is the solution of (111)-(112).

Case 2: g(u₀, v₀)>0. Then the solution of (111)-(112) is the tangent point of the circle g(u, v)=0 and the smallest elliptical level set ϕ(u, v)=r that intersects the circle (see FIG. 16). Accordingly, FIG. 16 illustrates an optimal solution based on tangent points between a circle and level sets. To find the maximizers/minimizers, one considers the following Lagrange multipliers system

$\begin{array}{cc}{\nabla }_u\varphi \left(u,v\right)+\frac{\lambda }{2}{\nabla }_{u}g\left(u,v\right)=0,& \left(115\right)\\ {\nabla }_u\varphi \left(u,v\right)+\frac{\lambda }{2}{\nabla }_{v}g\left(u,v\right)=0,& \left(116\right)\\ u^2+v^2-1=0.& \left(117\right)\end{array}$

Using the geometry of R² (see FIG. 16) and the theory of Lagrange multipliers, one can obtain the below properties.

Consider problem (111)-(112) and its Lagrange multipliers system (115)-(117). Suppose the solution (u₀, v₀) of (113)-(114) satisfies g(u₀, v₀)>0. Then

1. (115)-(117) yields a unique solution (u₁, v₁, λ₁) such that λ₁>0; and in this case, (u₁, v₁) is the unique solution of (111)-(112).

2. (115)-(117) yields at least one solution (u₂, v₂, λ₂) such that λ₂<0.

Therefore, one seeks a solution (u, v, λ) of (115)-(117) with λ>0. First, (115) and (116) yield(A+λ)u+Cv=w_a  (118)Cu+(B+λ)v=w_b.  (119)

Using the Cauchy-Schwarz inequality, one hasAB=[a₁²+a₂²+(a₁−a₂)²][b₁²+b₂²+(b₁−b₂)²](120)≥[a₁b₁+a₂b₂+(a₁−a₂)(b₁−b₂)]²=C².  (121)

Since A, B, λ>0, the determinant of (118)-(119) is(A+λ)(B+λ)−C²>AB−C²>0.  (122)This implies the system (118)-(119) has a unique solution

$\begin{array}{cc}u=\frac{w_a\left(B+\lambda \right)-w_{b}C}{\left(A+\lambda \right)\left(B+\lambda \right)-C^2}=\frac{\lambda w_a+w_{a}B-w_{b}C}{{\lambda }^2+\left(A+B\right)\lambda +\mathrm{AB}-C^2},& \left(123\right)\\ v=\frac{w_b\left(A+\lambda \right)-w_{a}C}{\left(A+\lambda \right)\left(B+\lambda \right)-C^2}=\frac{\lambda w_b+w_{b}A-w_{a}C}{{\lambda }^2+\left(A+B\right)\lambda +\mathrm{AB}-C^2},& \left(124\right)\end{array}$

Next, substituting u and v into (117) yields

$\begin{array}{cc}{\left[{\lambda }^2+\left(A+B\right)\lambda +\mathrm{AB}-C^2\right]}^2& \left(125\right)\\ ={\left[w_a\lambda +\left(w_{a}B-w_{b}C\right)\right]}^2+{\left[w_b\lambda +\left(w_{b}A-w_{a}C\right)\right]}^2& \left(126\right)\\ =\left(w_a^2+w_b^2\right){\lambda }^2+2\left(w_a^{2}B+w_b^{2}A-2w_a{w}_{b}C\right)\lambda & \left(127\right)\\ +{\left(w_{a}B-w_{b}C\right)}^2+{\left(w_{b}A-w_{a}C\right)}^2.& \left(128\right)\end{array}$

Defining the constants accordingly, one may rewrite as(λ²+R₁λ+R₂)²=R₃λ²+2R₄λ+R₅.  (129)

One can simply solve this equation by the classic Ferrari's method for quartic equations [9]. However, since there is a unique positive λ, one will find its explicit formula following Ferrari's technique. A real variable y is introduced(λ²+R₁λ+R₂+y)²=R₃λ²+2R₄λ+R₅+2(λ²+R₁λ+R₂)y+y²  (130)=(R₃+2y)λ²+2(R₄+R₁y)+R₅+2R₂y+y².  (131)

One can choose y such that the right hand side is a perfect square in λ. Thus, the right hand side must have zero discriminant, i.e.,[R₄+R₁y]²−(R₃+2y)(R₅+2R₂y+y²)=0  (132)−2y³+[R₁²−R₃−4R₂]y²+[2R₁R₄−2R₂R₃−2R₅]y+R₄²−R₃R₅=0.   (133)

This is a cubic equation in y, so Cardano's method [9] can be used to find one real solution y₀. Then (130)-(131) becomes

$\begin{array}{cc}{\left({\lambda }^2+R_1\lambda +R_2+y_0\right)}^2=\left(R_3+2y_0\right){\left(\lambda +\frac{R_4+R_1{y}_0}{R_3+2y_0}\right)}^2.& \left(134\right)\end{array}$ From the discussion above, this equation must have at least one positive and one negative solutions. Next, one solves for the (unique) positive Δ:

Case 2a: R₃+2y₀<0. Then

$\lambda =-\frac{R_4+R_1{y}_0}{R_3+2y_0}$ is the unique value of λ, which is a contradiction. Thus, this case cannot happen.

Case 2b. R₃+²y₀=0. Then

$y_0=-\frac{R_3}{2}$ andλ²+R₁λ+R₂+y₀=0.  (135)So one has

$\begin{array}{cc}\lambda =\frac{1}{2}\left(-R_1\pm \sqrt{R_1^2-4\left(R_2-\frac{R_3}{2}\right)}\right)=\frac{1}{2}\left(-R_1\pm \sqrt{R_1^2-4R_2+2R_3}\right).& \left(136\right)\end{array}$ Since there is only one positive λ, one takes

$\begin{array}{cc}\lambda =\frac{1}{2}\left(-R_1+\sqrt{R_1^2-4R_2+2R_3}\right).& \left(137\right)\end{array}$

Case 2c: R₃+2y₀>0. Define r:=√{square root over (R₃+2y₀)}. From (134),

$\begin{array}{cc}{\left({\lambda }^2+R_1\lambda +R_2+y_0\right)}^2={\left(\lambda \sqrt{R_3+2y_0}+\frac{R_4+R_1{y}_0}{\sqrt{R_3+2y_0}}\right)}^2& \left(138\right)\\ ={\left(\lambda r+\frac{R_4+R_1{y}_0}{r}\right)}^2.& \left(139\right)\end{array}$ This leads to two quadratic equations in λ

$\begin{array}{cc}{\lambda }^2+\left(R_1\pm r\right)\lambda +\left(R_2+y_0\pm \frac{R_4+R_1{y}_0}{r}\right)=0.& \left(140\right)\end{array}$

Since there is only one positive λ, one must have R₄+R₁y₀≠0. Otherwise, (140) consists of two quadratic equations that have constant term R₂+y₀; thus, will yield an even number (possibly none) of positive solutions, which is a contradiction. Moreover, one must take the equation with negative constant term. So one sets r←r·sgn(R₄+R₁y₀) and takes only the equation

$\begin{array}{cc}{\lambda }^2+\left(R_1-r\right)\lambda +\left(R_2+y_0-\frac{R_4+R_1{y}_0}{r}\right)=0.& \left(141\right))\end{array}$ It follows that the positive Δ is

$\begin{array}{cc}\lambda =\frac{1}{2}\left[r-R_1+\sqrt{{\left(r-R_1\right)}^2-4\left(R_2+y_0-\frac{R_4+R_1{y}_0}{r}\right)}\right].& \left(142\right)\end{array}$

Next, one obtains u and v from (123)-(124) and then rescales variables (u, v)←(su, sv). Finally, one can use (104)-(106) to obtain (h₁, h₂, h₃).

Projection onto Surface Oriented Minimum Slope Constraint

Nonconvex Surface Minimum Constraint

Let P₁P₂P₃ be a triangle with the normal vector {right arrow over (n)}=(−u, −v, 1) as described above. In some cases, it is required that the angle ∠({right arrow over (n)}, {right arrow over (q)})≥α for some given vector {right arrow over (q)} and number α.

This happens, for example, if P₁P₂P₃ must have a slope at least s:=s_min ∈R₊. In this case, the angle between {right arrow over (n)} and e₃=(0,0,1) satisfies∠({right arrow over (n)},e₃)≥tan⁻¹(s).  (143)

If one considers FIG. 15 again, in this case, the normal vector needs to be outside of the cone 1502. The inequality ∠({right arrow over (n)}, e₃)≥tan⁻¹ (s) is equivalent to

$\begin{array}{cc}\cos \angle \left(\stackrel{->}{n},e_3\right)=\frac{\langle \stackrel{->}{n},e_3\rangle }{\stackrel{->}{n}}\le \cos \left({\tan }^{-1}\left(s\right)\right)=\frac{1}{\sqrt{1+s^2}}.& \left(144\right)\end{array}$ One can refer to (143) or (144) as the surface minimum slope constraint with minimum slope s. From (144), one has

$\begin{array}{cc}\frac{1}{\sqrt{u^2+v^2+1}}\le \frac{1}{\sqrt{1+s^2}}\iff u^2+v^2-s^2\ge 0.& \left(145\right)\end{array}$ This is clearly a nonconvex constraint in (u, v). Despite the projection onto this constraint that is still available, nonconvexity may prevent iterative methods from convergence. Therefore, one will find a convex replacement for this constraint which can still guarantee (143)-(144).

The Surface Oriented Minimum Slope Constraint

(143) implies the angle between {right arrow over (n)} and the xy-plane satisfies∠({right arrow over (n)},xy−plane)≤(π/2)−tan⁻¹(s).  (146)In some cases, it is not unreasonable to align the plane P₁P₂P₃ toward a given location/direction. For instance, in civil engineering drainage, the designer may want to direct the water to certain drains or inlets. Suppose one wants P₁P₂P₃ to incline toward a unit direction {right arrow over (q)}=(q₁, q₂, 0)∈R³ (see FIG. 17). In this regard, FIG. 17 illustrates a desired orientation ({right arrow over (q)}) of a normal vector (n) for a triangle plane. Then one can fulfill (146) by requiring

$\angle \left(\stackrel{->}{n},\stackrel{->}{q}\right)\le \frac{\pi }{2}-{\tan }^{-1}\left(s\right),$ i.e.,

$\begin{array}{cc}\cos \angle \left(\stackrel{->}{n},\stackrel{->}{q}\right)=\frac{\langle \stackrel{->}{n},\stackrel{->}{q}\rangle }{\stackrel{->}{n}}\ge \cos \left(\frac{\pi }{2}-{\tan }^{-1}\left(s\right)\right)=\frac{s}{\sqrt{1+s^2}}& \left(147\right)\end{array}$ and arrive at the following definition.

Definition 3 [surface oriented minimum slope constraint]. One calls (147) the surface oriented minimum slope constraint specified by ({right arrow over (q)}, s), the unit vector {right arrow over (q)} and the minimum slope s. This is a special case of the surface orientation constraint where

$\begin{array}{cc}\left(\stackrel{->}{q},\theta \right)=\left(\stackrel{->}{q},\frac{\pi }{2}-{\tan }^{-1}\left(s\right)\right).& \end{array}$ Substituting {right arrow over (n)}=(−u, −v, 1) and {right arrow over (q)}=(q₁, q₂, 0) into (147), one has

$\begin{array}{cc}\frac{-q_{1}u-q_{2}v}{\sqrt{u^2+v^2+1}}\ge \frac{s}{\sqrt{1+s^2}},& \left(148\right)\end{array}$ which is equivalent to

$\begin{array}{cc}{q}_{1}u+q_{2}v<0,& \left(149\right)\\ u^2+v^2+1-\frac{1+s^2}{s^2}{\left(q_{1}u+q_{2}v\right)}^2\le 0.& \left(150\right)\end{array}$ This constraint not only guarantees surface minimum slope, but also generates a convex set in (u, v). A visualization of this finding is given in FIG. 18. In other words, FIG. 18 illustrates a visualization of a surface minimum slope constraint in accordance with one or more embodiments of the invention. The constraint requires that the normal vector lays on or inside the half cone 1802, whereas the direction vector q is the axis of the full cone.

Projection onto Surface Oriented Minimum Slope Constraint

By the model (82)-(83), the projection onto the surface oriented minimum slope constraint is

$\begin{array}{cc}\underset{\left(u,v\right)\in R^2}{\min }\varphi \left(u,v,h_3\right)={\left(a_{1}u+b_{1}v+h_3-w_1\right)}^2+{\left(a_{2}u+b_{2}v+h_3-w_2\right)}^2+{\left(h_3-w_3\right)}^2& \left(151\right)\\ \mathrm{subject}\mathrm{to}\left(149\right)-\left(150\right).& \left(152\right)\end{array}$

Again, one solves (151)-(152) directly using Lagrange multipliers and Ferrari's method for quartic equations. First, define

$Q:=\left[\begin{array}{cc}{q}_1& q_2\\ q_2& -q_1\end{array}\right],$ then Q²=Id since∥{right arrow over (q)}∥=q₁²+q₂²=1. So Q=Q^T=Q⁻¹. Next, the variables are changed

$\begin{array}{cc}\left[\begin{array}{c}\hat{u}\\ \hat{v}\end{array}\right]:=Q\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{cc}{q}_{1}u& q_{2}v\\ q_{2}u& -q_{1}v\end{array}\right]\iff \left[\begin{array}{c}u\\ v\end{array}\right]=Q\left[\begin{array}{c}\hat{u}\\ \hat{v}\end{array}\right]=\left[\begin{array}{cc}{q}_1\hat{u}& q_2\hat{v}\\ q_2\hat{u}& -q_1\hat{v}\end{array}\right].& \left(153\right)\end{array}$ Then (150) becomes

$\begin{array}{cc}{\left(q_1\hat{u}+q_2\hat{v}\right)}^2+{\left(q_2\hat{u}-q_1\hat{v}\right)}^2+1-\left(1+\frac{1}{s^2}\right){\hat{u}}^2=\left(-\frac{1}{s^2}\right){\hat{u}}^2+{\hat{v}}^2+1\le 0.& \left(154\right)\end{array}$ Thus, the constraint (149)-(150) becomes

$\begin{array}{cc}\hat{u}<0,& \left(155\right)\\ {\hat{v}}^2-\frac{{\hat{u}}^2}{s^2}+1\le 0,& \left(156\right)\end{array}$

Following the projections onto a general surface slope constraint described above, the coefficients are changed (without relabeling)

$\begin{array}{cc}\left[\begin{array}{cc}{a}_1& a_2\\ b_1& b_2\end{array}\right]\leftarrow Q\left[\begin{array}{cc}{a}_1& a_2\\ b_1& b_2\end{array}\right]& \left(157\right)\end{array}$ and the following problem is obtained

$\begin{array}{cc}\underset{\left(\hat{u},\hat{v}\right)\in R^2}{\min }\varphi \left(\hat{u},\hat{v}\right)=\frac{A}{2}{\hat{u}}^2+C\hat{u}\hat{v}+\frac{B}{2}{\hat{v}}^2-w_a\hat{u}-w_b\hat{v}& \left(158\right)\\ \mathrm{subject}\mathrm{to}\{\begin{array}{c}\hat{u}<0,\\ {\hat{v}}^2-\frac{{\hat{u}}^2}{s^2}+1\le 0\end{array}.& \left(159\right)\end{array}$ where A, B, C, w_a, w_b are defined below using the new a₁, a₂, b₁, b₂ from (157)A=2(a₁²a₂²−a₁a₂)>0,  (160)B=2(b₁²+b₂²−b₁b₂)>0,  (161)C=2a₁b₁+2a₂b₂−a₁b₂−a₂b₁,  (162)w_a=3(a₁w₁+a₂w₂),  (163)w_b=3(b₁w₁+b₂w₂).  (164)Changing variables by

$\left(u,v,A,B,w_b\right)\leftarrow \left(\frac{\hat{u}}{s},\hat{v},\mathrm{sA},\frac{B}{s},\frac{w_b}{s}\right),$ then (158)-(159) becomes

$\begin{array}{cc}\underset{\left(\hat{u},\hat{v}\right)\in R^2}{\min }\varphi \left(u,v\right)=\frac{A}{2}{u}^2+\mathrm{Cuv}+\frac{B}{2}{v}^2-w_{a}u-w_{b}v& \left(165\right)\\ \mathrm{subjectto}\{\begin{array}{c}{g}_1\left(u,v\right):=u<0,\\ g_2\left(u,v\right):=v^2-u^2+1\le 0\end{array}.& \left(166\right)\end{array}$

Hence, the projection to the surface oriented minimum slope constraint will find new elevations for the triangle vertices, such that the angle between the normal vector of the triangle and the z-axis is at least tan⁻¹(s) and the angle between the normal vector and the directional vector q is at most

$\frac{\pi }{2}-{\tan }^{-1}\left(s\right),$ and that equation (151) is minimized.

Before solving the above problem, one can observe that

1. The constraint (166) is one side of a hyperbola, thus, convex.

2. The objective ϕ(u, v) in (165) is a nonnegative strictly convex quadratic surface.

3. For each value η≥0, the level set ϕ(u, v)=η is an ellipse in R² whose center is the vertex (u₀, v₀) of ϕ(u, v), which satisfiesAu₀+Cv₀=w_a  (167)Cu₀+Bv₀=w_b.  (168)

Therefore, (165)-(166) has a unique solution. To solve it, two cases are considered:

Case 1: (u₀, v₀) is feasible, i.e., it lies inside or on the left branch {(u, v)∈R²|g₂(u, v)=0, u<}. Then it is the unique solution of (165)-(166) because the optimal objective value is zero. FIG. 19 illustrates a feasible set region for a normal vector with respect to a Surface Oriented Minimum Slope Constraint in accordance with one or more embodiments of the invention.

Case 2: (u₀, v₀) is not feasible. Let (u, v) be a solution to (165)-(166). By Lagrange multipliers, there exists λ>0 such that

$\begin{array}{cc}{\nabla }_u\varphi \left(u,v\right)+\frac{\lambda }{2}{\nabla }_u{g}_2\left(u,v\right)=0,& \left(169\right)\\ {\nabla }_v\varphi \left(u,v\right)+\frac{\lambda }{2}{\nabla }_v{g}_2\left(u,v\right)=0,& \left(170\right)\\ g_2\left(u,v\right)=v^2-u^2+1=0,& \left(171\right)\\ g_1\left(u,v\right)=u<0.& \left(172\right)\end{array}$

Consider the system of three equations (169)-(171). Using the geometry of R² (see FIGS. 20 and 21), one may conclude that it has at least two solutions. To describe its solutions, let (u₀, v₀) be defined by (167)-(168), then

1. If g₂(u₀, v₀)>0, then (169)-(171) has a unique solution (u₁, v₁, λ₁) with λ₁>0 and u₁<0, and a unique solution (u₂, v₂, λ₂) with λ₂>0 and u₂>0.

2. If g₂(u₀, v₀)≤0 and g₁(u₀, v₀)=u>0, then (169)-(171) has a unique solution (u₁, v₁, λ₁) with λ₁>0. Thus, u₁<0. In addition, (169)-(171) also has at least one solution (u₂, v₂, λ₂) with λ₂<0.

3. In both cases 1 and 2, (u₁, v₁, λ₁) is the unique solution of the whole Lagrange system (169)-(172), thus, (u₁, v₁) is the unique solution of (165)-(166).

To prove the above conclusion, suppose (u, v, λ) is a solution of (169)-(171), then (u, v) is the tangent point of the hyperbola v²−u²+1=0 and an elliptical level set of ϕ(u, v).

FIG. 20 shows that there are exactly two tangent points (u, v) with positive Lagrange multipliers λ, one of which has u<0 and the other has u>0 in accordance with one or more embodiments of the invention. This justifies conclusion 1.

FIG. 21 shows that there is exactly one tangent point (u, v) with a positive Lagrange multiplier λ and u>0, and at least one tangent point with a nonpositive Lagrange multiplier in accordance with one or more embodiments of the invention. This justifies conclusion λ.

Now if (u, v) is a solution to (165)-(166), then there exists λ>0 such that (u, v, λ) is a solution to (169)-(171). So, conclusion 3 follows from conclusions 1 and 2.

Next, one can show how to find the solution in Case 2. First, write (169)-(170) as(A−λ)u+Cv=w_a,  (173)Cu+(B+λ)v=w_b.  (174)Suppose for now that the determinant (A−λ)(B+λ)−C²≠0. Then

$\begin{array}{cc}u=\frac{w_a\left(B+\lambda \right)-w_{b}C}{\left(A-\lambda \right)\left(B+\lambda \right)-C^2}\hspace{1em}=\frac{\left(\lambda w_a+w_{a}B-w_{b}C\right)}{-{\lambda }^2+\left(A-B\right)\lambda +\left(\mathrm{AB}-C^2\right)},& \left(175\right)\\ v=\frac{w_b\left(A-\lambda \right)-w_{a}C}{\left(A-\lambda \right)\left(B+\lambda \right)-C^2}\hspace{1em}=\frac{-\lambda w_b+\left(w_{b}A-w_{a}C\right)}{-{\lambda }^2+\left(A-B\right)\lambda +\left(\mathrm{AB}-C^2\right)},& \left(176\right)\end{array}$ Substituting into (171), one can derive

$\begin{array}{cc}1={\left(\frac{\lambda w_a+w_{a}B-w_{b}C}{{\lambda }^2+\left(B-A\right)\lambda +\left(C^2-\mathrm{AB}\right)}\right)}^2-{\left(\frac{\lambda w_b+\left(w_{a}C-w_{b}A\right)}{{\lambda }^2+\left(B-A\right)\lambda +\left(C^2-\mathrm{AB}\right)}\right)}^2.& \left(177\right)\end{array}$ This implies

$\begin{array}{cc}\begin{array}{c}{\left[{\lambda }^2+\left(B-A\right)\lambda +\left(C^2-\mathrm{AB}\right)\right]}^2\\ ={\left(\lambda w_a+w_{a}B-w_{b}C\right)}^2-{\left[\lambda w_b+\left(w_{a}C-w_{b}A\right)\right]}^2\\ =\left(w_a^2-w_b^2\right){\lambda }^2+2\left(w_a^{2}B+w_b^{2}A-2w_a{w}_{b}C\right)\lambda \\ +{\left(w_{a}B-w_{b}C\right)}^2-{\left(w_{a}C-w_{b}A\right)}^2.\end{array}& \begin{array}{c}\left(178\right)\\ \left(179\right)\\ \left(180\right)\\ \begin{array}{c}\\ \left(181\right)\end{array}\end{array}\end{array}$

By defining the constants accordingly, one rewrites as(λ²+R₁λ+R₂)²=R₃λ²+2R₄λ+R₅.  (182)

Again, one can analyze this equation analogously to the projection onto the surface maximum slope constraint described above. However, complications arise since there are possibly more than one positive λ's. Instead, one can expand (182)λ⁴+2R₁λ³+(R₁²+2R₂−R₃)λ²+2(R₁R₂−R₄)Δ+R₂²−R₅=0,  (183)and use the classical Ferrari's method for quartic equation. Then for each λ>0, defineD=(A−λ)(B+λ)−C²,  (184)D_u=w_a(B+λ)−w_bC,  (185)D_v=w_b(A−λ)−w_aC.  (186)

Case 2a: D=D_u=D_v=0. Then one has the systemCu+(B+λ)v=w_b,  (187)v²−u²+1=0.  (188)

Since B>0 and λ>0, the first equation implies v=(w_b−Cu)/(B+λ). Then the second one becomes:(w_b−Cu)²−(B+λ)²u²+(B+λ)²=0,i.e.,  (189)[C²−(B+λ)²]u²−2w_bCu+[w_b²+(B+λ)²]≥0.  (190)If the determinant(w_bC)²−(C²−(B+λ)²)(w_b²+(B+λ)²)≥0,  (191)then two solutions u are obtained. Since there is a unique pair (u, v) with u<0, the quadratic equation (190) must yield one positive and one negative solutions

$\begin{array}{cc}u=\frac{w_{b}C\pm \sqrt{{\left(w_{b}C\right)}^2-\left(C^2-{\left(B+\lambda \right)}^2\right)\left(w_b^2+{\left(B+\lambda \right)}^2\right)}}{C^2-{\left(B+\lambda \right)}^2},& \left(192\right)\end{array}$ and one must also have (C²−(B+λ)²)(w_b²+(B+λ)²)<0. So, the negative solution u is

$\begin{array}{cc}u=\frac{w_{b}C+\sqrt{{\left(w_{b}C\right)}^2-\left(C^2-{\left(B+\lambda \right)}^2\right)\left(w_b^2+{\left(B+\lambda \right)}^2\right)}}{C^2-{\left(B+\lambda \right)}^2}.& \left(193\right)\end{array}$

Case 2b: D≠0. Then by (175)-(176), u=D_u/D and v=D_v/D.

Next, among all (u, v)'s, choose the pair with u<0. Then rescale variables (u, v)←(su, v). Finally, one can use (104)-(106) to obtain (h₁, h₂, h₃).

Multiple Surface Orientation Constraints

It is possible to impose multiple surface orientation constraints on one triangle. However, it might cause infeasibility. For instance, it is easy to see that one cannot “tilt” the normal vector toward opposite directions. Therefore, care must be taken when enforcing multiple orientation constraints. Below is a description of the feasibility in such cases.

Feasibility Criterion for One Maximum and One Oriented Minimum Slope Constraints

Given a triangle with unit normal vector {right arrow over (n)}. Suppose constraints are enforced on {right arrow over (n)} as follows:cos ∠({right arrow over (n)},e₃)≥1/√{square root over (1+t²)}(maximum slope)  (194)cos ∠({right arrow over (n)},{right arrow over (q)})≥s/√{square root over (1+s²)}(oriented minimum slope)  (195)where e₃=(0,0,1); {right arrow over (q)}=(q₁, q₂, 0) ≡(q₁, q₂) is a given unit vector; and t, s∈R₊ are also given. Then the problem is feasible if and only if t≥s.

To verify this, assume {right arrow over (n)}=(u, v, w) is a unit vector (w>0), and its shadow on xy-plane is {right arrow over (n)}₀=(u, v). Then (194) implies

$\begin{array}{cc}w\ge \frac{1}{\sqrt{1+t^2}}\iff u^2+v^2=1-w^2\le 1-\frac{1}{1+t^2}=\frac{t^2}{1+t^2},& \left(196\right)\end{array}$ which is a disk on R² centered at the origin with radius t/√{square root over (1+t²)}, and (195) implies

$\begin{array}{cc}{\mathrm{uq}}_1+{\mathrm{vq}}_2\ge \frac{s}{\sqrt{1+s^2}},& \left(197\right)\end{array}$ which is a circular segment 2202 with height

$\frac{t}{\sqrt{1+t^2}}-\frac{s}{\sqrt{1+s^2}}$ (see FIG. 22). Thus, FIG. 22 illustrates a circular segment defining a feasible set for a maximum slope and an oriented minimum slope in accordance with one or more embodiments of the invention.

So {right arrow over (n)}₀=(u, v) must lie in the area 2202, which is nonempty if and only if

$\frac{t}{\sqrt{1+t^2}}\ge \frac{s}{\sqrt{1+s^2}}\iff \frac{t^2}{1+t^2}\ge \frac{s^2}{1+s^2}\iff t\ge s.$

Feasibility Criterion for One Maximum and Two Oriented Minimum Slope Constraints

Given a triangle with unit normal vector {right arrow over (n)}. Suppose {right arrow over (n)} satisfies:

$\begin{array}{cc}\cos \angle \left(\overrightarrow{n},e_3\right)\ge \frac{1}{\sqrt{1+t^2}},\cos \angle \left(\overrightarrow{n},{\overrightarrow{q}}_1\right)\ge \frac{s_1}{\sqrt{1+s_1^2}},\mathrm{and}\cos \angle \left(\overrightarrow{n},{\overrightarrow{q}}_2\right)\ge \frac{s_2}{\sqrt{1+s_2^2}}.& \left(198\right)\end{array}$ where e₃=(0,0,1), {right arrow over (q)}_i=(q_i1, q_i2, 0)≡(q_i1 q_i2) are unit vectors, and t, s₁, s₂ ∈R₊. Then the problem is feasible if and only if

$\hspace{1em}\begin{array}{cc}\left({\overrightarrow{q}}_1,{\overrightarrow{q}}_2\right)\le {\cos }^{-1}\left(\frac{s_1\sqrt{1+t^2}}{t\sqrt{1+s_1^2}}\right)+{\cos }^{-1}\left(\frac{s_2\sqrt{1+t^2}}{t\sqrt{1+s_2^2}}\right).& \left(199\right)\end{array}$

To see this, the (feasible) disk and circular segments (i.e., regions 2302 and 2304) are illustrated in FIG. 23, where

$r:=\frac{t}{\sqrt{1+t^2}},r_1:=\frac{s_1}{\sqrt{1+s_1^2}},\mathrm{and}{r}_2:=\frac{s_2}{\sqrt{1+s_2^2}}.$ Thus, the intersection is nonempty if and only if

$\angle \left({\stackrel{->}{q}}_1,{\stackrel{->}{q}}_2\right)\le {\cos }^{-1}\left(\frac{r_1}{r}\right)+{\cos }^{-1}\left(\frac{r_2}{r}\right).$ Accordingly, FIG. 23 illustrates a circle representing a feasible region for a maximum slope constraint and regions 2302 and 2304 show the feasible region for two oriented minimum slope constraints.

Feasibility Criterion for Multiple Surface Oriented Minimum Slope Constraints:

Given a triangle with unit normal vector {right arrow over (n)}, s, t∈ R, and a collection of unit vectors Q:={{right arrow over (q)}_i}_i∈1 ⊂R². Let t≥s>0. Suppose the constraints are enforced on n as follows:cos ∠({right arrow over (n)},e₃)≥1/√{square root over (1+t²)},  (200)∀i∈I:cos ∠({right arrow over (n)},{right arrow over (q)}_i)≥s/√{square root over (1+s²)}.  (201)Then it suffices to enforce 2 constraints in (201) since the rest will be automatically fulfilled.

To verify this, one can assume {{right arrow over (q)}₁, {right arrow over (q)}₂}⊂Q is the pair that make the largest angle (see FIG. 24) and the area S is the intersection of two corresponding circular segments. Accordingly, FIG. 24 illustrates multiple surface oriented minimum slope constraints via the intersection of two corresponding circular segments in accordance with one or more embodiments of the invention.

Similar to previous cases, one can illustrate the disk of radius

$\frac{t}{\sqrt{1+t^2}}$ and several circular segments of the same height in FIG. 24

$\left(\mathrm{where}r=\frac{s}{\sqrt{1+s^2}}\right).$ One can argue that all other {right arrow over (q)}_i's lie in between {right arrow over (q)}₁ and {right arrow over (q)}₂ (otherwise, the intersection with S is empty). Notice that the distances from the origin to all circular segments are equal. Thus, one can deduce that S satisfies all other oriented minimum slope constraint {right arrow over (q)}_i. Thus, it suffices to enforce at most two oriented minimum slope constraints.

Further, one can associate each {right arrow over (q)}_i with different minimum slope s_i and then obtain analogous feasibility criterion.

Curvature Minimization

In some design problems, the designer may wish to construct a surface that is as “smooth” as possible. This problem is referred to as minimizing the curvature between adjacent triangles in the mesh.

Given the points P_i=(p_i1, p_i2, h_i), i=1, 2, 3, 4, so that they form two adjacent triangles Δ₁=P₁P₂P₄ and Δ₂=P₂P₃P₄ (see FIG. 25). In this regard, FIG. 25 illustrates two adjacent triangles with surface normal vectors that are used to minimize the curvature in accordance with one or more embodiments of the invention.

Define a_i:=p_i1−p₄₁ and b_i:=p_i2−p₄₂, for i=1, 2, 3. Then the respective normal vectors are computed by the cross products

$\begin{array}{cc}{\overrightarrow{n}}_1=\left[\begin{array}{c}{a}_1\\ b_1\\ h_1-h_4\end{array}\right]\times \left[\begin{array}{c}{a}_2\\ b_2\\ h_2-h_4\end{array}\right]=\left[\begin{array}{c}-b_2{h}_1+b_1{h}_2+\left(b_2-b_1\right)h_4\\ a_2{h}_1-a_1{h}_2+\left(a_1-a_1\right)h_4\\ a_1{b}_2-a_2{b}_1\end{array}\right]& \left(202\right)\\ {\overrightarrow{n}}_2=\left[\begin{array}{c}{a}_2\\ b_2\\ h_2-h_4\end{array}\right]\times \left[\begin{array}{c}{a}_3\\ b_3\\ h_3-h_4\end{array}\right]=\left[\begin{array}{c}-b_3{h}_2+b_2{h}_3+\left(b_3-b_2\right)h_4\\ a_3{h}_2-a_2{h}_3+\left(a_2-a_3\right)h_4\\ a_2{b}_3-a_3{b}_2\end{array}\right]& \left(203\right)\end{array}$ One can rescale and obtain

$\begin{array}{cc}{\overrightarrow{n}}_1=\left(\frac{-b_2{h}_1+b_1{h}_2+\left(b_2-b_1\right)h_4}{a_1{b}_2-a_2{b}_1},\frac{a_2{h}_1-a_1{h}_2+\left(a_1-a_2\right)h_4}{a_1{b}_2-a_2{b}_1},1\right),& \left(204\right)\\ {\overrightarrow{n}}_2=\left(\frac{-b_3{h}_2+b_1{h}_3+\left(b_3-b_2\right)h_4}{a_2{b}_3-a_3{b}_2},\frac{a_3{h}_2-a_2{h}_3+\left(a_2-a_3\right)h_4}{a_2{b}_3-a_3{b}_2},1\right).& \left(205\right)\end{array}$ The curvature can be represented by the difference between {right arrow over (n)}₁ and {right arrow over (n)}₂, i.e.,

$\begin{array}{cc}{\overrightarrow{\delta }}_{12}:={\overrightarrow{\delta }}_{{\Delta }_1{\Delta }_2}={\overrightarrow{n}}_1-{\overrightarrow{n}}_2=\text{:}\left[\begin{array}{c}\langle u,h\rangle \\ \langle v,h\rangle \end{array}\right]& \left(206\right)\end{array}$ where

$\begin{array}{cc}u:=\left(u_1,u_2,u_3,u_4\right),v:=\left(v_1,v_2,v_3,v_4\right),h:=\left(h_1,h_2,h_3,h_4\right),& \left(207\right)\\ u_1=\frac{-b_2}{a_1{b}_2-a_2{b}_1},u_2=\frac{b_1}{a_1{b}_2-a_2{b}_1}+\frac{b_3}{a_2{b}_3-a_3{b}_2},& \left(208\right)\\ u_3=\frac{-b_2}{a_2{b}_3-a_3{b}_2},u_4=-\left(u_1+u_2+u_3\right),& \left(209\right)\\ v_1=\frac{a_2}{a_1{b}_2-a_2{b}_1},v_2=\frac{-a_1}{a_1{b}_2-a_2{b}_1}+\frac{-a_3}{a_2{b}_3-a_3{b}_2},& \left(210\right)\\ v_3=\frac{a_2}{a_2{b}_3-a_3{b}_2},v_4=-\left(v_1+v_2+v_3\right).& \left(211\right)\end{array}$

Now let be the set of all pairs of adjacent triangles (Δ_i, Δ_j) for which it is desirable to minimize the curvature. Then for each pair (Δ_i, Δ_j)∈, one finds the corresponding vectors h_ij=(h₁^ij, h₂^ij, h₃^ij, h₄^ij)∈R⁴, u_ij, v_ij∈R⁴, and computes the corresponding curvature δ_ij=(u_ij, h_ij, v_ij, h_ij). One can then aim to minimize all the computed curvatures. Thus, one can arrive at the objective

$\begin{array}{cc}{G}_{1,*}\left(x\right):=\sum _{\left({\Delta }_i,{\Delta }_j\right)\in 𝒵}{{\overrightarrow{\delta }}_{\mathrm{ij}}}_{*}& \left(212\right)\end{array}$ where ∥⋅∥_* can be either 1-norm or max-norm in R². For 1-norm, the objective is

$\begin{array}{cc}{G}_{1,1}\left(x\right):=\sum _{\left({\Delta }_i,{\Delta }_j\right)\in 𝒵}\langle u_{\mathrm{ij}},h_{\mathrm{ij}}\rangle +\langle v_{\mathrm{ij}},h_{\mathrm{ij}}\rangle ,& \left(213\right)\end{array}$ and for max-norm, the objective is

$\begin{array}{cc}{G}_{1,\infty }\left(x\right):=\sum _{\left({\Delta }_i,{\Delta }_j\right)\in 𝒵}\max \left\{\langle u_{\mathrm{ij}},h_{\mathrm{ij}}\rangle +\langle v_{\mathrm{ij}},h_{\mathrm{ij}}\rangle \right\}.& \left(214\right)\end{array}$

Simplified Computations for Symmetric Cases.

Suppose the two dimensional mesh satisfies the following symmetry: for every pair of adjacent triangles (P₁P₂P₄, P₂P₃P₄)∈, there exists t∈R such that{right arrow over (P₄P₁)}+{right arrow over (P₄P₃)}=t{right arrow over (P₄P₂)}.  (215)Then it follows that (a₁, b₁₎+(a₃, b₃)=t(a₂, b₂). So one can deduce a₁b₂−a₂b₁=a₂b₃−a₃b₂. Using (207)-(211), one obtains

$\begin{array}{cc}u=\frac{-b_2}{a_1{b}_2-a_2{b}_1}\left(1,-t,1,t-2\right)\mathrm{and}v=\frac{a_2}{a_1{b}_2-a_2{b}_1}\left(1,-t,1,t-2\right).& \left(216\right)\end{array}$ That means u and v are parallel. This simplifies the computations for (213) and (214).

Proximity Operators

Since optimization methods can require proximity operators, one will derive the necessary formulas. It is sometimes convenient to compute the proximity operator P_ƒ via the proximity operator of its Fenchel conjugate ƒ*, which is defined by ƒ*:X→R:x*sup_x∈X(x*, x−ƒ(x)). Indeed, if γ>0, then (see, e.g., [2, Theorem 14.3(ii)])∀x∈X:x=P_γƒ(x)+γP_γ−1ƒ*(γ⁻¹x).  (217)One may also recall a useful formula from [2, Lemma 2.3]: if ƒ:X→R is convex and positively homogeneous, α>0, w∈X, and h:X→R:xαƒ(x−w), then

$\begin{array}{cc}{P}_h\left(x\right)=w+\alpha P_f\left(\frac{x-w}{\alpha }\right)=x-\alpha P_{f^{*}}\left(\frac{x-w}{\alpha }\right).& \left(218\right)\end{array}$

Formula for Proximity Operators:

Let {u_i}_i∈I be a system of finitely many vectors in Rⁿ, and

$\begin{array}{cc}f\text{:}{R}^{″}\to R\text{:}x\to \underset{i\in I}{\max }\left\{\langle u_i,x\rangle \right\}.& \left(219\right)\end{array}$

Then ƒ*=t_D where D:=convU_i∈I{u_i,−u_i}. Consequently, the proximity operator of ƒ can be computed by P_ƒ=Id−P_D, where P_D is the projection onto the set D.

To see this, first, suppose x*∈D, then one can express

$\begin{array}{cc}{x}^{*}=\sum _{i\in I}{\lambda }_i{u}_i\mathrm{where}\sum _{i\in I}{\lambda }_i\le 1.& \left(220\right)\end{array}$ It follows that for all x∈X,

$\begin{array}{cc}\langle x^{*},x\rangle -f\left(x\right)=\sum _{i\in I}{\lambda }_i\langle u_i,x\rangle -\underset{i\in I}{\max }\langle u_i,x\rangle \le \left(\sum _{i\in I}{\lambda }_i-1\right)\underset{i\in I}{\max }\langle u_i,x\rangle \le 0.& \left(221\right)\end{array}$ So, ƒ* (x*)=sup_x∈X[x*, x−ƒ(x)]≤0. Notice that equality happens if one sets x=0. Thus, ƒ*(x*)=0.

Now suppose x*∉D. Since D is nonempty, closed, and convex, the classic separation theorem implies that there exists x∈X such thatx*,x>u,x for all u∈D  (222)This leads to x*, x−ƒ(x)>0. Since ƒ is homogeneous, one can haveƒ*(x*)≥x*,λx−ƒ(λx)→+∞ as λ→+∞.  (223)So ƒ*(x*)=+∞. Therefore, one concludes that ƒ*=t_D. Now employing formula (218) and noting that P_ƒ*=P_D, the projection onto D (see the beginning of the iterative projection methods optimization problems description above), one obtains the proximity operator for ƒ.

Next, one presents two special cases of the above proximity formula, which will be used in curvature minimization.

Example 1

Given α>0, a vector u ∈Rⁿ \{0} and the functionƒ:Rⁿ→R: x α|u,x|.  (224)Then the proximity operator of ƒ isP_ƒ=Id−P_D where D:=[−αu,αu].  (225)The explicit projection onto D is given below (see also, [4, Theorem 2.7]).

$\begin{array}{cc}{P}_{D}x=\min \left\{1,\max \left\{-1,\frac{\langle \alpha u,x\rangle }{{\alpha u}^2}\right\}\right\}\alpha u=\min \left\{\alpha ,\max \left\{-\alpha ,\frac{\langle u,x\rangle }{{u}^2}\right\}\right\}u.& \left(226\right)\end{array}$

Example 2

Given two vectors u, v∈Rⁿ andƒ:Rⁿ→R:xmax{|u,x|,|v,x|}.  (227)Then the proximity operator of ƒ isP_ƒ=Id−P_D where D:=conv{u,v,−u,−v}.  (228)In general, P_D is the projection onto a parallelogram in Rⁿ.Approximation of Surface Slope Constraint Projections

In this section, the problem of surface slope constraint projections are considered. Such projections may be written in the form

$\begin{array}{cc}\underset{u\in R^2}{\min }\phi \left(u\right):=\frac{A}{2}{u}_1^2+{\mathrm{Cu}}_1{u}_2+\frac{B}{2}{u}_2^2-w_a{u}_1-w_b{u}_2& \left(229\right)\end{array}$ subject to u∈Ω:={u∈R²|g(u)≤0}.  (230)

If one defines

$M:=\left[\begin{array}{cc}A& C\\ C& B\end{array}\right]\mathrm{and}w:=\left[\begin{array}{c}{W}_a\\ W_b\end{array}\right],$ then

$\begin{array}{cc}\phi \left(u\right)=\frac{1}{2}{u}^T\mathrm{Mu}-w^{T}u.& \left(231\right)\end{array}$ In practice, solving (229)-(230) directly as described above (e.g., as described above in the “Projection onto the surface maximum slope constraint” section and the “Projection onto Surface Oriented Minimum Slope Constraint” section) might not be straightforward as one may encounter numerical instability. Therefore, one can use an alternative approach where the constraint Ω is approximated by a regular polygon.

Approximation for Maximum Slope Constraint

Consider the projection onto maximum slope constraint written in the form

$\begin{array}{cc}\underset{u\in R^2}{\min }\phi \left(u\right):=\frac{1}{2}{u}^T\mathrm{Mu}-w^{T}u& \left(232\right)\\ \mathrm{subject}\mathrm{to}u\in \Omega :=\{\left(u_1,u_2\right)\in R^2u_1^2+u_2^2\le 1\}.& \left(233\right)\end{array}$ The unit disk Ω may be replaced by a convex (or simple) regular polygon that is contained in Ω. For example, the unit disk may be replaced by a hexagon (see FIG. 26). In this regard, FIG. 26 illustrates a hexagon that is used to replace a unit disk for a maximum slope constraint in accordance with one or more embodiments of the invention. More specifically, let m≥3 be an integer and 0≤θ₁≤θ₂≤ . . . ≤θ_m<2π. Define the points on the unit circle p_k:=(cos θ_k, sin θ_k) for k=1, 2, . . . , m, and set also p₀:=p_m. Now define the convex m-polygon U:=conv{p_k|k=0, 1, . . . , m}. One then arrives at the approximate problem

$\begin{array}{cc}\underset{u\in R^2}{\min }\phi \left(u\right):=\frac{1}{2}{u}^T\mathrm{Mu}-w^{T}u\mathrm{subject}\mathrm{to}u\in U.& \left(234\right)\end{array}$ This problem may be solved as follows. Let u₀:=M⁻¹w be the vertex of the paraboloid φ(u).

Case 1: u₀ ∈Ω. Then u₀ is even the solution of (232)-(233), thus, the problem is solved.

Case 2: u₀ ∉C. So one also has u₀∉U. The solution of (234) may be found by the following steps.

Step 2a: Find the index setK:={k|0 and u₀ are on opposite sides of the line p_kp_k+1}.  (235)

Step 2b: Minimize φ(u) on each edge p_kp_k+1 for k∈K. Consider φ(u) on an edge p_kp_k+1. Define Δp_k=p_k+1−p_k. Take a point u∈p_kp_k+1, thenu=p_k+tΔp_k for some t∈[0,1].  (236)So

$\begin{array}{cc}\psi \left(t\right):=\phi \left(u\right)=\frac{1}{2}{\left(p_k+t\Delta p_k\right)}^{T}M\left(p_k+t\Delta p_k\right)-w^T\left(p_k+t\Delta p_k\right)& \mathrm{ }\left(237\right)\\ =\frac{t^2}{2}\Delta p_k^{T}M\Delta p_k-t\left(w^T\Delta p_k-p_k^{T}M\Delta p_k\right)+\frac{1}{2}{p}_k^T{\mathrm{Mp}}_k-& \left(238\right)\\ w^T{p}_k.& \end{array}$ The vertex of ψ(t) is

$\stackrel{\_}{t}:=\frac{w^T\Delta p_k-p_k^{T}M\Delta p_k}{\Delta p_k^{T}M\Delta p_k},$ so the minimizer of ψ(t) on [0,1] ist_k:=min{1,max{0,t}}.  (239)

Step 2c: If there exists k∈K such that 0<t_k<1, then u=p_k+t_kΔp_k is the solution of (234). Otherwise, the solution of (234) is

$\begin{array}{cc}u=p_{k_0}+t_{k_0}\Delta p_{k_0}\mathrm{where}{k}_0:=\underset{k\in K}{\mathrm{argmin}}\phi \left(p_k+t_k\Delta p_k\right).& \left(240\right)\end{array}$

Hence, the projection to the surface maximum slope constraint is approximated by using a regular polygon instead of a unit disk.

Approximation for Oriented Minimum Slope Constraint

Consider the projection onto oriented minimum slope constraint written in the form

$\begin{array}{cc}\underset{u\in R^2}{\min }\phi \left(u\right):=\frac{1}{2}{u}^T\mathrm{Mu}-w^{T}u& \left(241\right)\\ \mathrm{subject}\mathrm{to}u\in \Omega :=\left\{\left(u_1,u_2\right)\in R^2|u_1<0,u_2^2-\frac{u_1^2}{s^2}+1\le 0\right\}.& \left(242\right)\end{array}$ The constraint Ω may be replaced byU:={(u₁,u₂)∈R²|u₁≤−s(|u₂|+1)}.  (243)FIG. 27 shows how the original half cone is now replaced with the absolute value function in accordance with one or more embodiments of the invention. One can easily check that U⊂Ω. So one arrives at an approximate problem

$\begin{array}{cc}\underset{u\in R^2}{\min }\phi \left(u\right):=\frac{1}{2}{u}^T\mathrm{Mu}-w^{T}u\mathrm{subject}\mathrm{to}u\in U.& \left(244\right)\end{array}$ Let u₀:=M⁻¹w be the vertex of the paraboloid φ(u).

Case 1: u₀ ∈C. Then u₀ is even the solution of (241)-(242), thus, the problem is solved.

Case 2: u₀ ∈C. So u₀ ∈U. One solves (244) by first finding the minimizer of φ(u) on the boundary of U, which consists of two half-lines{p+tΔp₁|t∈R₊} and {p+tΔp₂|tΣR₊},  (245)where p:=(−s, 0), Δp₁:=(−s, 1), and Δp₂:=(−s, −1).Let u:=p+tΔp for t∈[0, +∞). Similar to the computation above, one has

$\begin{array}{cc}\psi \left(t\right):=\phi \left(u\right)=\frac{t^2}{2}\Delta p^{T}M\Delta p-t\left(w^T\Delta p-p^{T}M\Delta p\right)+\frac{1}{2}{p}^T\mathrm{Mp}-w^{T}p.& \left(246\right)\end{array}$ So the minimizer of φ(u) on [0, +∞) is

$\begin{array}{cc}t:=\max \left\{0,\frac{w^T\Delta p-p^{T}M\Delta p}{\Delta p^{T}M\Delta p}\right\}.& \left(247\right)\end{array}$ Second, one sets Δp=Δp₁ and Δp=Δp₂ in the above formula to obtain t₁ and t₂, respectively. Finally, comparing the values φ(p+t₁Δp₁) and φ(p+t₂Δp₂), one obtains the minimizer of (244).

Visualization Examples

The closed-form proximity and projection operators can be used in the DR method to solve various design problems described above.

FIG. 28 shows the starting situation for the problem shown in FIG. 7 in accordance with one or more embodiments of the invention. Triangles 2802 represent triangles that violate the design constraints, and triangles 2902 (from FIG. 29) represent triangles that satisfy the design constraints. FIG. 29 shows the mesh after 10 iterations in accordance with one or more embodiments of the invention. FIG. 30 shows the mesh after 29 iterations, and FIG. 31 shows the final mesh after convergence in accordance with one or more embodiments of the invention.

Consider FIG. 28 again, this time as a starting situation for the problem shown in FIG. 8. FIG. 32 is after 15 iterations, FIG. 33 is after 100 iterations, and FIG. 34 is the final solution.

FIG. 35 is the starting situation for the problem shown in FIG. 9. FIG. 36 is after 10 iterations, FIG. 37 is after 20 iterations, and FIG. 38 is the final solution.

Logical Flow

FIG. 39 illustrates the logical flow for designing a surface in accordance with one or more embodiments of the invention.

At step 3902, a triangular surface mesh representative of an existing surface is obtained. The triangular surface mesh includes two or more triangles that are connected by vertices and edges. In embodiments of the invention, the surface may be a land surface (e.g., that is actually constructed based on the design [e.g., a computer-aided design—CAD]).

At step 3904, one or more design constraint sets are determined based on one or more design constraints. The design constraints can include various constraints including a maximum slope constraint for a first triangle of the two or more triangles in the triangular surface mesh. The maximum slope constraint is a maximum angle between a normal vector of the first triangle and a reference vector (e.g., a Z-vector representative of gravity).

An additional design constraint may be a surface oriented minimum slope constraint for the first triangle. The surface oriented minimum slope constraint is a minimum angle between the normal vector of the first triangle and a directional vector.

At step 3906, one or more heights of the vertices of the first triangle are projected onto the one or more design constraint sets such that the normal vector satisfies all of the one or more design constraints. The projecting includes modifying the one or more heights by a minimum Euclidian distance. Further, the projecting may be performed iteratively until all of the one or more design constraints are satisfied. In addition, the projecting onto the design constraint set for the maximum slope constraint may include moving the normal vector onto a cone defined by an origin of the reference vector on the triangular surface mesh and a unit circle around the reference vector that is defined by the maximum angle. Alternatively, or in addition, the projecting onto to the design constraint set for the maximum slope constraint may include moving the normal vector onto a polygonal shaped pyramid defined by the reference vector on the triangular surface mesh and a regular polygon around the reference vector, wherein edges of the polygonal shaped pyramid are at the maximum angle with the reference vector.

At step 3908, a design of the surface represented by the triangular surface mesh is generated (based on the projecting).

Hardware Environment

FIG. 40 is an exemplary hardware and software environment 4000 used to implement one or more embodiments of the invention. The hardware and software environment includes a computer 4002 and may include peripherals. Computer 4002 may be a user/client computer, server computer, or may be a database computer. The computer 4002 comprises a general purpose hardware processor 4004A and/or a special purpose hardware processor 4004B (hereinafter alternatively collectively referred to as processor 4004) and a memory 4006, such as random access memory (RAM). The computer 4002 may be coupled to, and/or integrated with, other devices, including input/output (I/O) devices such as a keyboard 4014, a cursor control device 4016 (e.g., a mouse, a pointing device, pen and tablet, touch screen, multi-touch device, etc.) and a printer 4028. In one or more embodiments, computer 4002 may be coupled to, or may comprise, a portable or media viewing/listening device 4032 (e.g., an MP3 player, IPOD, NOOK, portable digital video player, cellular device, personal digital assistant, etc.). In yet another embodiment, the computer 4002 may comprise a multi-touch device, mobile phone, gaming system, internet enabled television, television set top box, or other internet enabled device executing on various platforms and operating systems.

In one embodiment, the computer 4002 operates by the general purpose processor 4004A performing instructions defined by the computer program 4010 under control of an operating system 4008. The computer program 4010 and/or the operating system 4008 may be stored in the memory 4006 and may interface with the user and/or other devices to accept input and commands and, based on such input and commands and the instructions defined by the computer program 4010 and operating system 4008, to provide output and results.

Output/results may be presented on the display 4022 or provided to another device for presentation or further processing or action. In one embodiment, the display 4022 comprises a liquid crystal display (LCD) having a plurality of separately addressable liquid crystals. Alternatively, the display 4022 may comprise a light emitting diode (LED) display having clusters of red, green and blue diodes driven together to form full-color pixels. Each liquid crystal or pixel of the display 4022 changes to an opaque or translucent state to form a part of the image on the display in response to the data or information generated by the processor 4004 from the application of the instructions of the computer program 4010 and/or operating system 4008 to the input and commands. The image may be provided through a graphical user interface (GUI) module 4018. Although the GUI module 4018 is depicted as a separate module, the instructions performing the GUI functions can be resident or distributed in the operating system 4008, the computer program 4010, or implemented with special purpose memory and processors.

In one or more embodiments, the display 4022 is integrated with/into the computer 4002 and comprises a multi-touch device having a touch sensing surface (e.g., track pod or touch screen) with the ability to recognize the presence of two or more points of contact with the surface. Examples of multi-touch devices include mobile devices (e.g., IPHONE, NEXUS S, DROID devices, etc.), tablet computers (e.g., IPAD, HP TOUCHPAD), portable/handheld game/music/video player/console devices (e.g., IPOD TOUCH, MP3 players, NINTENDO 3DS, PLAYSTATION PORTABLE, etc.), touch tables, and walls (e.g., where an image is projected through acrylic and/or glass, and the image is then backlit with LEDs).

Some or all of the operations performed by the computer 4002 according to the computer program 4010 instructions may be implemented in a special purpose processor 4004B. In this embodiment, some or all of the computer program 4010 instructions may be implemented via firmware instructions stored in a read only memory (ROM), a programmable read only memory (PROM) or flash memory within the special purpose processor 4004B or in memory 4006. The special purpose processor 4004B may also be hardwired through circuit design to perform some or all of the operations to implement the present invention. Further, the special purpose processor 4004B may be a hybrid processor, which includes dedicated circuitry for performing a subset of functions, and other circuits for performing more general functions such as responding to computer program 4010 instructions. In one embodiment, the special purpose processor 4004B is an application specific integrated circuit (ASIC).

The computer 4002 may also implement a compiler 4012 that allows an application or computer program 4010 written in a programming language such as C, C++, Assembly, SQL, PYTHON, PROLOG, MATLAB, RUBY, RAILS, HASKELL, or other language to be translated into processor 4004 readable code. Alternatively, the compiler 4012 may be an interpreter that executes instructions/source code directly, translates source code into an intermediate representation that is executed, or that executes stored precompiled code. Such source code may be written in a variety of programming languages such as JAVA, JAVASCRIPT, PERL, BASIC, etc. After completion, the application or computer program 4010 accesses and manipulates data accepted from I/O devices and stored in the memory 4006 of the computer 4002 using the relationships and logic that were generated using the compiler 4012.

The computer 4002 also optionally comprises an external communication device such as a modem, satellite link, Ethernet card, or other device for accepting input from, and providing output to, other computers 4002.

In one embodiment, instructions implementing the operating system 4008, the computer program 4010, and the compiler 4012 are tangibly embodied in a non-transitory computer-readable medium, e.g., data storage device 4020, which could include one or more fixed or removable data storage devices, such as a zip drive, floppy disc drive 4024, hard drive, CD-ROM drive, tape drive, etc. Further, the operating system 4008 and the computer program 4010 are comprised of computer program 4010 instructions which, when accessed, read and executed by the computer 4002, cause the computer 4002 to perform the steps necessary to implement and/or use the present invention or to load the program of instructions into a memory 4006, thus creating a special purpose data structure causing the computer 4002 to operate as a specially programmed computer executing the method steps described herein. Computer program 4010 and/or operating instructions may also be tangibly embodied in memory 4006 and/or data communications devices 4030, thereby making a computer program product or article of manufacture according to the invention. As such, the terms “article of manufacture,” “program storage device,” and “computer program product,” as used herein, are intended to encompass a computer program accessible from any computer readable device or media.

Of course, those skilled in the art will recognize that any combination of the above components, or any number of different components, peripherals, and other devices, may be used with the computer 4002.

FIG. 41 schematically illustrates a typical distributed/cloud-based computer system 4100 using a network 4104 to connect client computers 4102 to server computers 4106. A typical combination of resources may include a network 4104 comprising the Internet, LANs (local area networks), WANs (wide area networks), SNA (systems network architecture) networks, or the like, clients 4102 that are personal computers or workstations (as set forth in FIG. 40), and servers 4106 that are personal computers, workstations, minicomputers, or mainframes (as set forth in FIG. 40). However, it may be noted that different networks such as a cellular network (e.g., GSM [global system for mobile communications] or otherwise), a satellite based network, or any other type of network may be used to connect clients 4102 and servers 4106 in accordance with embodiments of the invention.

A network 4104 such as the Internet connects clients 4102 to server computers 4106. Network 4104 may utilize ethernet, coaxial cable, wireless communications, radio frequency (RF), etc. to connect and provide the communication between clients 4102 and servers 4106. Further, in a cloud-based computing system, resources (e.g., storage, processors, applications, memory, infrastructure, etc.) in clients 4102 and server computers 4106 may be shared by clients 4102, server computers 4106, and users across one or more networks. Resources may be shared by multiple users and can be dynamically reallocated per demand. In this regard, cloud computing may be referred to as a model for enabling access to a shared pool of configurable computing resources.

Clients 4102 may execute a client application or web browser and communicate with server computers 4106 executing web servers 4110. Such a web browser is typically a program such as MICROSOFT INTERNET EXPLORER, MOZILLA FIREFOX, OPERA, APPLE SAFARI, GOOGLE CHROME, etc. Further, the software executing on clients 4102 may be downloaded from server computer 4106 to client computers 4102 and installed as a plug-in or ACTIVEX control of a web browser. Accordingly, clients 4102 may utilize ACTIVEX components/component object model (COM) or distributed COM (DCOM) components to provide a user interface on a display of client 4102. The web server 4110 is typically a program such as MICROSOFT'S INTERNET INFORMATION SERVER.

Web server 4110 may host an Active Server Page (ASP) or Internet Server Application Programming Interface (ISAPI) application 4112, which may be executing scripts. The scripts invoke objects that execute business logic (referred to as business objects). The business objects then manipulate data in database 4116 through a database management system (DBMS) 4114. Alternatively, database 4116 may be part of, or connected directly to, client 4102 instead of communicating/obtaining the information from database 4116 across network 4104. When a developer encapsulates the business functionality into objects, the system may be referred to as a component object model (COM) system. Accordingly, the scripts executing on web server 4110 (and/or application 4112) invoke COM objects that implement the business logic. Further, server 4106 may utilize MICROSOFT'S TRANSACTION SERVER (MTS) to access required data stored in database 4116 via an interface such as ADO (Active Data Objects), OLE DB (Object Linking and Embedding DataBase), or ODBC (Open DataBase Connectivity).

Generally, these components 4100-4116 all comprise logic and/or data that is embodied in/or retrievable from device, medium, signal, or carrier, e.g., a data storage device, a data communications device, a remote computer or device coupled to the computer via a network or via another data communications device, etc. Moreover, this logic and/or data, when read, executed, and/or interpreted, results in the steps necessary to implement and/or use the present invention being performed.

Although the terms “user computer”, “client computer”, and/or “server computer” are referred to herein, it is understood that such computers 4102 and 4106 may be interchangeable and may further include thin client devices with limited or full processing capabilities, portable devices such as cell phones, notebook computers, pocket computers, multi-touch devices, and/or any other devices with suitable processing, communication, and input/output capability.

Of course, those skilled in the art will recognize that any combination of the above components, or any number of different components, peripherals, and other devices, may be used with computers 4102 and 4106. Further, embodiments of the invention are implemented as a software application on a client 4102 or server computer 4106. In addition, as described above, the client 4102 or server computer 4106 may comprise a thin client device or a portable device that has a multi-touch-based display.

CONCLUSION

This concludes the description of the preferred embodiment of the invention. The following describes some alternative embodiments for accomplishing the present invention. For example, any type of computer, such as a mainframe, minicomputer, or personal computer, or computer configuration, such as a timesharing mainframe, local area network, or standalone personal computer, could be used with the present invention.

The foregoing description of the preferred embodiment of the invention has been presented for the purposes of illustration and description. It is not intended to be exhaustive or to limit the invention to the precise form disclosed. Many modifications and variations are possible in light of the above teaching. It is intended that the scope of the invention be limited not by this detailed description, but rather by the claims appended hereto.

REFERENCES

• [1] BAUSCHKE, H. H., AND BORWEIN, J. M. On projection algorithms for solving convex feasibility problems. SIAM Rev. 38, 3 (1996), 367-426.
• [2] BAUSCHKE, H. H., AND COMBETTES, P. L. Convex analysis and monotone operator theory in Hilbert spaces. CMS Books in Mathematics/Ouvrages de Mathématiques de la SMC. Springer, New York, 2011. With a foreword by Hédy Attouch.
• [3] BAUSCHKE, H. H., AND KOCH, V. R. Projection methods: Swiss army knives for solving feasibility and best approximation problems with halfspaces. In Infinite products of operators and their applications, vol. 636 of Contemp. Math. Amer. Math. Soc., Providence, R.I., 2015, pp. 1-40.
• [4] BAUSCHKE, H. H., KOCH, V. R., AND PHAN, H. M. Stadium norm and Douglas-Rachford splitting: a new approach to road design optimization. Oper. Res. 64, 1 (2016), 201-218.
• [5] BINIAZ, A., AND DASTGHAIBYFARD, G. Drainage reality in terrains with higher-order Delaunay triangulations. Springer Berlin Heidelberg, Berlin, Heidelberg, 2008, pp. 199-211.
• [6] BINIAZ, A., AND DASTGHAIBYFARD, G. Slope fidelity in terrains with higher-order delaunay triangulations. In 16th International Conference in Central Europe on Computer Graphics, Visualization, and Computer Vision (Plzen, Czech Republic, 2008), Václav Skala-UNION Agency, pp. 17-23.
• [7] BOT, R. I., CSETNEK, E. R., AND HEINRICH, A. A primal-dual splitting algorithm for finding zeros of sums of maximal monotone operators. SIAM J. Optim. 23, 4 (2013), 2011-2036.
• [8] BRICEÑO ARIAS, L. M., AND COMBETTES, P. L. A monotone+skew splitting model for composite monotone inclusions in duality. SIAM J. Optim. 21, 4 (2011), 1230-1250.
• [9] CARDANO, G. The great art or The rules of algebra. The M.I.T. Press, Cambridge, Mass.-London, 1968. Translated from the Latin and edited by T. Richard Witmer, With a foreword by Oystein Ore.
• [10] CEGIELSKI, A. Iterative methods for fixed point problems in Hilbert spaces, vol. 2057 of Lecture Notes in Mathematics. Springer, Heidelberg, 2012.
• [11] CENSOR, Y., CHEN, W., COMBETTES, P. L., DAVIDI, R., AND HERMAN, G. T. On the effectiveness of projection methods for convex feasibility problems with linear inequality constraints. Comput. Optim. Appl. 51, 3 (2012), 1065-1088.
• [12] CENSOR, Y., AND ZENIOS, S. A. Parallel optimization. Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 1997. Theory, algorithms, and applications, With a foreword by George B. Dantzig.
• [13] CORNISH, G., AND MUIR GRAVES, R. Classic Golf Hole Design. John Wiley & Sons, Inc., Hoboken DC, 2002.
• [14] DOUGLAS, JR., J., AND RACHFORD, JR., H. H. On the numerical solution of heat conduction problems in two and three space variables. Trans. Amer. Math. Soc. 82 (1956), 421-439.
• [15] HAWTREE, F. W. The Golf Course, Planning design, construction and maintenance. E & FN Spon, and imprint of Chapman & Hall, London, U K, 2005.
• [16] IRVING, R. Beyond the Quadratic Formula. The Mathematical Association of America, Washington D.C., 2010.
• [17] ITOH, T., AND SHIMADA, K. Automatic conversion of triangular meshes into quadrilateral meshes with directionality. International Journal of CADCAM 1, 1 (2002), 11-21.
• [18] JOHNSTON, B. P., SULLIVAN, J. M., AND KWASNIK, A. Automatic conversion of triangular finite element meshes to quadrilateral elements. International Journal for Numerical Methods in Engineering 31, 1 (1991), 67-84.
• [19] KARUSH, W. MINIMA OF FUNCTIONS OF SEVERAL VARIABLES WITH INEQUALITIES AS SIDE CONDITIONS. ProQuest LLC, Ann Arbor, Mich., 1939. Thesis (SM)—The University of Chicago.
• [20] KUHN, H. W., AND TUCKER, A. W. Nonlinear programming. In Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, 1950 (1951), University of California Press, Berkeley and Los Angeles, pp. 481-492.
• [21] LIONS, P.-L., AND MERCIER, B. Splitting algorithms for the sum of two nonlinear operators. SIAM J. Numer. Anal. 16, 6 (1979), 964-979.
• [22] MOORE, I.D., GRAYSON, R., AND LADSON, A. Digital terrain modelling: a review of hydrological, geomorphological, and biological applications. Hydrological processes 5, 1 (1991), 3-30.
• [23] MOREAU, J.-J. Proximité et dualité dans un espace hilbertien. Bull. Soc. Math. France 93 (1965), 273-299.
• [24] RANK, E., SCHWEINGRUBER, M., AND SOMMER, M. Adaptive mesh generation and transformation of triangular to quadrilateral meshes. Communications in Numerical Methods in Engineering 9, 2 (1993), 121-129.
• [25] ROCKAFELLAR, R. T. Convex analysis. Princeton Mathematical Series, No. 28. Princeton University Press, Princeton, N.J., 1970.
• [26] RUSZCZYŃSKI, A. Nonlinear optimization. Princeton University Press, Princeton, N.J., 2006.

Claims

1. A computer-implemented method for designing a surface, comprising: (a) obtaining, in a computer, a triangular surface mesh representative of an existing physical surface, wherein the triangular surface mesh comprises two or more triangles that are connected by vertices and edges;(b) determining, in the computer, one or more design constraint sets based on one or more design constraints, wherein at least one of the one or more design constraints comprise: a maximum slope constraint for a first triangle of the two or more triangles in the triangular surface mesh, wherein the maximum slope constraint comprises a maximum angle between a normal vector of the first triangle and a reference vector;(c) projecting, in the computer, one or more heights of the vertices of the first triangle onto the one or more design constraint sets such that the normal vector satisfies all of the one or more design constraints, wherein the projecting comprises modifying the one or more heights by a minimum Euclidian distance, and wherein the modifying changes the triangular surface mesh into a converged triangular surface mesh;(d) generating, in the computer, a design of the surface represented by the converged triangular surface mesh, wherein the generating is based on the projecting; and(e) rendering the design of the surface in a computer-aided design (CAD) application, wherein the surface is constructed based on the design.

2. The computer-implemented method of claim 1, wherein the one or more design constraints further comprise: a surface oriented minimum slope constraint for the first triangle, wherein the surface oriented minimum slope constraint comprises a minimum angle between the normal vector of the first triangle and a directional vector.

3. The computer-implemented method of claim 1, wherein: the surface comprises a land surface; andthe method further comprises constructing the land surface based on the design.

4. The computer-implemented method of claim 3, wherein: the design comprises a computer-aided design (CAD) for the land surface.

5. The computer-implemented method of claim 1, wherein the reference vector comprises a Z-vector representative of gravity.

6. The computer-implemented method of claim 1, wherein: the projecting is performed iteratively until all of the one or more design constraints are satisfied.

7. The computer-implemented method of claim 1, wherein: the projecting onto the design constraint set for the maximum slope constraint comprises moving the normal vector onto a cone defined by an origin of the reference vector on the triangular surface mesh and a unit circle around the reference vector that is defined by the maximum angle.

8. The computer-implemented method of claim 1, wherein: the projecting onto to the design constraint set for the maximum slope constraint comprises moving the normal vector onto a polygonal shaped pyramid defined by the reference vector on the triangular surface mesh and a regular polygon around the reference vector, wherein edges of the polygonal shaped pyramid are at the maximum angle with the reference vector.

9. A system for designing a surface, comprising: (a) a computer having a memory and a processor;(b) a computer aided design (CAD) application executed by the processor;(c) a triangular surface mesh, obtained in the CAD application, wherein the triangular surface mesh is representative of an existing physical surface and comprises two or more triangles that are connected by vertices and edges;(d) one or more design constraint sets, determined by the CAD application, based on one or more design constraints, wherein at least one of the one or more design constraints comprise: a maximum slope constraint for a first triangle of the two or more triangles in the triangular surface mesh, wherein the maximum slope constraint comprises a maximum angle between a normal vector of the first triangle and a reference vector;(e) the CAD application projecting, in the computer, one or more heights of the vertices of the first triangle onto the one or more design constraint sets such that the normal vector satisfies all of the one or more design constraints, wherein the projecting comprises modifying the one or more heights by a minimum Euclidian distance, and wherein the modifying changes the triangular surface mesh into a converged triangular surface mesh; and(d) a design of the surface represented by the converged triangular surface mesh, wherein the design is generated and rendered in the CAD application based on the projecting, wherein the surface is constructed based on the design.

10. The system of claim 9, wherein the one or more design constraints further comprise: a surface oriented minimum slope constraint for the first triangle, wherein the surface oriented minimum slope constraint comprises a minimum angle between the normal vector of the first triangle and a directional vector.

11. The system of claim 9, wherein: the surface comprises a land surface; andthe land surface is constructed based on the design.

12. The system of claim 11, wherein: the design comprises a computer-aided design (CAD) for the land surface.

13. The system of claim 9, wherein the reference vector comprises a Z-vector representative of gravity.

14. The system of claim 9, wherein: the application performs the projecting iteratively until all of the one or more design constraints are satisfied.

15. The system of claim 9, wherein: the application projects onto the design constraint set for the maximum slope constraint by moving the normal vector onto a cone defined by an origin of the reference vector on the triangular surface mesh and a unit circle around the reference vector that is defined by the maximum angle.

16. The system of claim 9, wherein: the application projects onto to the design constraint set for the maximum slope constraint by moving the normal vector onto a polygonal shaped pyramid defined by the reference vector on the triangular surface mesh and a regular polygon around the reference vector, wherein edges of the polygonal shaped pyramid are at the maximum angle with the reference vector.