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.

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; and the 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; and the 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.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

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

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

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

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

(8) 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;

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

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

(11) 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;

(12) 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;

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

(14) 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;

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

(16) 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;

(17) 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;

(18) 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;

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

(20) 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;

(21) 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;

(22) 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;

(23) 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;

(24) 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;

(25) 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;

(26) 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;

(27) 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;

(28) 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;

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

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

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

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

(33) 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;

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

(35) 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;

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

(37) 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;

(38) 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;

(39) 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;

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

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

(42) 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

(43) 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.

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

(45) 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 .fwdarw.. 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)}.sub.1 and {right arrow over (n)}.sub.2 in R.sup.3 is denoted by {right arrow over (n)}.sub.1{right arrow over (n)}.sub.2.

(46) General Overview

(47) 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.

(48) Problem Formulation

(49) 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.

(50) Triangular Mesh

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

(52) Assume that G=(V.sub.0,E.sub.0) is a given triangular mesh on R.sup.2 where V.sub.0 is the set of vertices and E.sub.0 is the set of (undirected) edges, i.e.,
V.sub.0:={p.sub.i=(p.sub.i1,p.sub.i2)R.sup.2|i{1,2, . . . ,n}},(1)
E.sub.0{p.sub.ip.sub.jp.sub.jp.sub.i|p.sub.i,p.sub.jV.sub.0}.(2)
From G, one can derive the set of triangles of the mesh as follows
T.sub.0:={=(p.sub.ip.sub.jp.sub.k)p|{p.sub.ip.sub.j,p.sub.jp.sub.k,p.sub.kp.sub.i}E.sub.0}.(3)

(53) One first aims to
find a vector z=(z.sub.1, . . . ,z.sub.n)X(4)
such that the points
{P.sub.i=(p.sub.i1,p.sub.i2,z.sub.i)}.sub.i{1, . . . ,n}satisfy a given set of constraints.(5)

(54) Clearly, the points (P.sub.i).sub.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.sup.2). FIG. 9 illustrates a simple triangular mesh with vertices
V.sub.0:={P.sub.1,P.sub.2,P.sub.3,P.sub.4},
edges
E.sub.0:={P.sub.1P.sub.4,P.sub.1P.sub.2,P.sub.2P.sub.4,P.sub.3P.sub.4,P.sub.2P.sub.4},
and the two triangles
T.sub.0:={(P.sub.1P.sub.2P.sub.4),(P.sub.2P.sub.3P.sub.4)}.
in accordance with one or more embodiments of the invention.

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

(56) Interval Constraint

(57) 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 iI, the entries z.sub.i must lie in a given interval of R.

(58) Edge Slope Constraint

(59) 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.sub.0, and for every edge p.sub.ip.sub.jE, the slope

(60) $\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.

(61) Edge Alignment Constraint

(62) 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.sub.ip.sub.j, p.sub.mp.sub.nE.sub.0, the edges must have the same slope s.sub.ij=s.sub.mn.

(63) Surface Alignment Constraint

(64) 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 .sub.i, .sub.jT.sub.0 the normal vectors {right arrow over (n)}.sub..sub.i and {right arrow over (n)}.sub..sub.j must be parallel.

(65) Surface Slope Constraint

(66) 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.sub.0 and for each triangle T, there is a given set of vectors Q.sub. R.sup.3 such that for each {right arrow over (q)}Q.sub., the angle between the normal vector {right arrow over (n)}.sub. and {right arrow over (q)} must lie in a given subset .sub.{right arrow over (q)} of [0, ], which is also referred to as a surface orientation constraint.

(67) Curvature Minimization

(68) 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 .sub.i, .sub.jT.sub.0 the normal vectors {right arrow over (n)}.sub..sub.i and {right arrow over (n)}.sub..sub.j we want to minimize the angle between the two normal vectors.

(69) Methods Overview

(70) 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.

(71) 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.

(72) 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.

(73) Projections onto Constraint Sets

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

(75) $\begin{array}{cc}\mathrm{find}a\mathrm{point}xC:=\underset{j\left\{1,2,\mathrm{.Math.},J\right\}}{.Math.}{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.

(76) 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

(77) $\begin{array}{cc}\min F\left(z\right)\mathrm{subject}\mathrm{to}zC:=\underset{j\left\{1,2,\mathrm{.Math.},J\right\}}{.Math.}{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.

(78) A constraint set CX is the collection of all feasible data points, i.e., points that satisfy some requirements. Suppose the given data point zX is not feasible (i.e., z.Math.C), one aims to modify z so that the newly obtained point x is feasible (i.e., xC ); 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.sub.Cz, is the solution of the optimization problem

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

(80) It is well known that if C is nonempty, closed, and convex.sup.1, then P.sub.Cz is singleton, see, e.g., [26, Theorem 2.10]. .sup.1C is convex if for all x, yC and [0,1], we have (1)x+yC

(81) 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.

(82) Iterative Projection Methods for Feasibility Problems

(83) Consider the feasibility problem

(84) $\begin{array}{cc}\mathrm{find}a\mathrm{point}xC:=\underset{j\left\{1,2,\mathrm{.Math.},J\right\}}{.Math.}{C}_j,& \left(10\right)\end{array}$

(85) 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.sub.C.sub.j onto each constraint, if the explicit formula is available. Then with an initial point, one can iteratively execute the projections P.sub.C.sub.j's to derive a solution for (10). Let z.sub.0 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).

(86) Cyclic Projections

(87) 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.sub.1 that is projected to the set A, resulting in the new point a.sub.0 that then gets projected to the set B. The obtained point b.sub.0 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.sub.k, one updates z.sub.k+1:=Tz.sub.k where
T:=P.sub.C.sub.JP.sub.C.sub.J-1 . . . P.sub.C.sub.2P.sub.C.sub.1.(11)

(88) Parallel Projections

(89) 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.sub.k, one updates z.sub.k+1:=Tz.sub.k where

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

(91) String-Average Projections

(92) The string-average projections is a combination of cyclic and parallel projections. One defines it as follows. Given z.sub.k, one updates z.sub.k+1:=Tz.sub.k where

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

(94) 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.

(95) Iterative Projection Methods for Optimization Problems

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

(97) $\begin{array}{cc}X.fwdarw.\left(-,+\right]\text{:}yf\left(y\right)+\frac{1}{2}{.Math..Math.x-y.Math..Math.}^2& \left(14\right)\end{array}$
has a unique minimizer, which will be denoted by P.sub. (x). The induced operator P.sub.: X.fwdarw.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.sub.C is defined by t.sub.C(x)=0 if xC and t.sub.C(x)=+ otherwise), then P.sub.=P.sub.C. Thus, proximity operators are generalizations of projections.

(98) The Douglas-Rachford Method

(99) 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.

(100) Using indicator functions, (8) is converted into
min F(x)+t.sub.C.sub.1+t.sub.C.sub.2+ . . . +t.sub.C.sub.J subject to xX.(15)

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

(102) $\begin{array}{cc}\min \underset{i=1}{\overset{m}{.Math.}}{f}_i\left(x\right)\mathrm{subject}\mathrm{to}xX,& \left(16\right)\end{array}$
where each .sub.i is a proper convex lower semicontinuous function on X. The DR operates in the product space X:=X.sup.m with inner product x,y:=.sub.i=1.sup.mx.sub.i, y.sub.i for x=(x.sub.i).sub.i=1.sup.m and y=(y.sub.i).sub.i=1.sup.m. Set the starting point x.sub.0=(z, . . . , z)X, where zX. Given x.sub.k=(x.sub.k,l, . . . , x.sub.k,m)X, one computes

(103) 0$\begin{array}{cc}{\stackrel{\_}{x}}_k:=\frac{1}{m}\underset{iK}{.Math.}{x}_{k,i},& \left(17\right)\end{array}$
i{1, . . . ,m}:x.sub.k+1,i:=x.sub.k,ix.sub.k+P.sub..sub.i(2x.sub.k,ix.sub.k,i),(18)
then update x.sub.k+1:=(x.sub.k+1,l, . . . ,x.sub.k+1,m).(19)
Then the monitored sequence (x.sub.k).sub.kN converges to a solution of (16).

(104) It is worth mentioning that when all .sub.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.

(105) Parallelization

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

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

(108) Suppose X=R.sup.n and the objective function : X.fwdarw.(, +] only depends on certain entries, i.e., (x)=((x.sub.i).sub.iI), where I{1, 2, . . . , n} is a given index set.

(109) Let z=(z.sub.1, . . . , z.sub.n)X and execute the proximity operator P.sub.Z. It is obvious that after the execution, the coordinates that possibly change are (z.sub.i).sub.iI; while all other coordinates (z.sub.i).sub.i.Math.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.sub.i).sub.iI in the subspace R.sup.|I| where |I| is the number of indices in I.

(110) Remark 2 (Parallelization)

(111) If two functions .sub.1, .sub.2: X.fwdarw.(,+] 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.

(112) Projections onto Linear Constraints

(113) 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.

(114) 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.

(115) Interval Constraint

(116) Similar to [3, Section 2.2], one may assume that I is a subset of {1, 2, . . . , n} and (l.sub.i).sub.iI, (u.sub.i).sub.iIR.sup.I are given. Set
Y:={z=(z.sub.1,z.sub.2, . . . ,z.sub.n)X|iI:l.sub.iz.sub.iu.sub.i}.(20)

(117) 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.

(118) Method Component 1 (Projection onto Interval Constraints)

(119) $\begin{array}{cc}{P}_Y\text{:}X.fwdarw.X\text{:}\left(z_1,z_2,\mathrm{.Math.},z_n\right)\left(x_1,x_2,\mathrm{.Math.},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\}& iI,\\ z_i,& i\left\{1,2,\mathrm{.Math.},n\right\}\backslash I\end{array}.& \left(22\right)\end{array}$

(120) Edge Minimum Slope Constraint

(121) Let P.sub.i=(p.sub.i1, p.sub.i2, h.sub.i) and P.sub.j=(p.sub.j1, p.sub.j2, h.sub.j). The designer may require that the (directional) slope from P.sub.j to P.sub.i must be no less than a threshold level s.sub.ij. More specifically, since the value p.sub.i1, p.sub.i2, p.sub.j1, p.sub.j2 are fixed, one can write this constraint as
h.sub.ih.sub.j.sub.ij:=s.sub.ij{square root over ((p.sub.i1p.sub.j1).sup.2+(p.sub.i2q.sub.j2).sup.2)},(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.

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

(123) Let
E.sub.i,j={x=(x.sub.1,x.sub.2, . . . ,x.sub.n)R.sup.n|x.sub.ix.sub.j.sub.ij}(24)
be the minimum slope constraint on the edge P.sub.iP.sub.j. Let x:=(x.sub.1, . . . , x.sub.n):=P.sub.E.sub.i,jz be the projection of a point z=(z.sub.1, z.sub.2, . . . , z.sub.n) on to E.sub.i,j. Then

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

(125) Low Point Constraint

(126) A point P.sub.j=(p.sub.j1, p.sub.j2, z.sub.j) on the mesh is called a low point if each point P.sub.i=(p.sub.i1, p.sub.i2, z.sub.i) connected to P.sub.j satisfies the edge minimum slope constraint
z.sub.iz.sub.j.sub.i:=s.sub.i{square root over ((p.sub.i1p.sub.j1).sup.2+(p.sub.i2q.sub.j2).sup.2)},(28)
where all s.sub.iR are given.

(127) FIG. 12 shows a triangular mesh from a birds-eye view, where the low point P.sub.0 is connected to points P.sub.1, P.sub.2, P.sub.3, P.sub.4, P.sub.5, and P.sub.6. So in this picture, one requires that P.sub.0 is lower than all its connected points and that the slopes on each connecting edge i is at least s.sub.i.

(128) 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.

(129) Method Component 3 (Projection onto a Low Point Constraint)

(130) Let .sub.2, . . . , .sub.mR. Define the set
E:={(x.sub.1, . . . ,x.sub.m)R.sup.m|i{2, . . . ,m}:x.sub.ix.sub.1.sub.i}.(29)

(131) Let z=(z.sub.1, . . . , z.sub.m)R.sup.m. Let .sub.1=z.sub.1 and let .sub.2.sub.3 . . . .sub.m be the values {z.sub.i.sub.i}.sub.i{2, . . . , m} in nondecreasing order. Let k be the largest number in {1, . . . , m} such that .sub.k(.sub.1+ . . . +.sub.k)/k. Then the projection x:=(x.sub.1, x.sub.2, . . . , x.sub.m):=P.sub.EZ is given by
x.sub.1=(.sub.1+ . . . +.sub.k)/k and x.sub.i=max{x.sub.1,z.sub.i.sub.i}+.sub.i,i{2, . . . ,m}. (30)

(132) To show this formula, first, x is the solution of
min(x.sub.1z.sub.1).sup.2+ . . . +(x.sub.mz.sub.m).sup.2(31)
s.t x.sub.1x.sub.i+.sub.i0,i{2, . . . ,m}.(32)

(133) Let y:=(y.sub.1, . . . , y.sub.m) where y.sub.1=x.sub.1 and y.sub.i:=x.sub.i.sub.i for i{2, . . . , m}. Without loss of generality, one may assume .sub.i=z.sub.i.sub.i for i{2, . . . , m}. Then (4) becomes
min (y):=(y.sub.1.sub.1).sup.2+ . . . +(y.sub.m.sub.m).sup.2(33)
s.t g.sub.i(y):=y.sub.1y.sub.i0,i{2, . . . ,m}.(34)

(134) To find the (unique) solution, we use Lagrange multipliers for convex optimization: there exist .sub.i0, {2, . . . , m}, such that

(135) $\begin{array}{cc}\frac{1}{2}\left(y\right)+{}_2{g}_2\left(y\right)+\mathrm{.Math.}+{}_m{g}_m\left(y\right)=0& \left(35\right)\\ i\left\{2,\mathrm{.Math.},m\right\}\text{:}{}_i{g}_i\left(y\right)=0& \left(36\right)\\ i\left\{2,\mathrm{.Math.},m\right\}\text{:}{g}_i\left(y\right)0.& \left(37\right)\end{array}$

(136) It follows from (35) that
(y.sub.1.sub.1)+.sub.2+ . . . +.sub.m=0(38)
(y.sub.2.sub.2).sub.2=0(39)
. . .(40)
(y.sub.m.sub.m).sub.m=0.(41)

(137) Then y.sub.1.sub.1 because .sub.2, . . . , .sub.m0. By substitution, one gets
y.sub.1+y.sub.2+ . . . +y.sub.m=.sub.1+.sub.2+ . . . .sub.m.(42)

(138) Next, for each i{2, . . . , m}, one has y.sub.i.sub.i=.sub.i0, y.sub.iy.sub.10, and (36) implies,
0=.sub.ig.sub.i(y)=(y.sub.i.sub.i)(y.sub.1y.sub.i).(43)
This leads to either y.sub.i=.sub.iy.sub.1 or y.sub.1=y.sub.i.sub.i. That means
y.sub.i=max{y.sub.1,.sub.i} for all i{2, . . . ,m}.(44)

(139) So (42) becomes
.sub.1+.sub.2+ . . . +.sub.m=y.sub.1+max{y.sub.1,.sub.2}+ . . . +max{y.sub.1,.sub.m}.(45)
Next, suppose k is the largest number in {1, . . . , m} such that
.sub.k(.sub.1+ . . . +.sub.k)/k.(46)

(140) The choice of k is always possible since at least (46) is true for k=1. Now setting .sub.m+1=+, one claims that
(.sub.1+ . . . +.sub.k)/k<.sub.k+1.(47)

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

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

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

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

(145) Edge Alignment Constraint

(146) 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.sub.1, A.sub.2, . . . , A.sub.m on the triangular mesh where A.sub.1=(a.sub.i1, a.sub.i2, x.sub.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.sub.iA.sub.i+1, the length of its shadow on the xy-plane is
.sub.i:=(a.sub.i+1,1a.sub.i,1).sup.2+(a.sub.i+1,2a.sub.i,2).sup.2.(61)
Define also
t.sub.1:=0,t.sub.2:=.sub.1,t.sub.3:=.sub.1+.sub.2, . . . ,t.sub.m:=.sub.1+ . . . +.sub.m1.(62)
Then the alignment constraint is written as
i{1, . . . ,m2}:(x.sub.i+1x.sub.i)/(t.sub.i+1t.sub.i)=(x.sub.i+2x.sub.i+1)/(t.sub.i+2t.sub.i+1).(63)

(147) Method Component 4 (Projection onto an Edge Alignment Constraint)

(148) Suppose the points (a.sub.i1, a.sub.i2), i{1, . . . , m} forms a path in R.sup.2. Define (t.sub.1, . . . , t.sub.m) by (62) and define the convex set
F:={(x.sub.1, . . . ,x.sub.m)R.sup.m|((a.sub.i1,a.sub.i2,x.sub.i)).sub.i{i, . . . ,m}satisfies(63)}R.sup.m.(64)

(149) Let z=(z.sub.1, . . . , z.sub.m)R.sup.m. Then the projection P.sub.Fz is given by
i {1, . . . ,m}:(P.sub.FZ).sub.i=(t.sub.i).(65)
where (t):=+t is the least squares line for the points (t.sub.i, z.sub.i).sub.i{1, . . . m}. More specifically, (, ) is obtained from the linear system

(150) $\begin{array}{cc}\left[\begin{array}{cc}m& \underset{i=1}{\overset{m}{.Math.}}{t}_i\\ \underset{i=1}{\overset{m}{.Math.}}{t}_i& \underset{i=1}{\overset{m}{.Math.}}{t}_i^2\end{array}\right]\left[\begin{array}{c}\\ \end{array}\right]=\left[\begin{array}{c}\underset{i=1}{\overset{m}{.Math.}}{z}_i\\ \underset{i=1}{\overset{m}{.Math.}}{t}_i{z}_i\end{array}\right]& \left(66\right)\end{array}$

(151) To show this, first, F is convex since all constraints in (63) are linear. Next, one considers the points (t.sub.i, z.sub.i) in R.sup.2. The problem is to find (x.sub.1, . . . , x.sub.m) such that the points (t.sub.i, x.sub.i) are aligned and
xz.sup.2=E.sub.i=1.sup.m(x.sub.iz.sub.i).sup.2 is minimized.(67)

(152) This is the least squares problem for the points (t.sub.i, z.sub.i), shown in FIG. 14. So the least squares line is (t)=+t where (, ) satisfies (66). Thus, the projection coordinates are (P.sub.F (z)).sub.i=+t.sub.i for i{1, . . . , m}.

(153) Surface Alignment Constraint

(154) 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.

(155) Method Component 5 (Projection onto a Surface Alignment Constraint)

(156) Let (a.sub.i1, a.sub.i2), i{1, . . . , m} be a collection of points in R.sup.2 that are not on the same line. Define
F:={(x.sub.1, . . . ,x.sub.m)R.sup.m|{(a.sub.i1,a.sub.i2,x.sub.i)}.sub.i{1, . . . ,m}lie on the same plane}. (69)
Let z=(z.sub.1, . . . , z.sub.m)R.sup.m. Then the projection P.sub.Fz is given by
i{1, . . . ,m}:(P.sub.Fz).sub.i=(a.sub.i1,a.sub.i2).(70)
where (t.sub.1, t.sub.2):=+t.sub.1+t.sub.2 is the least squares plane for the points (a.sub.i1, a.sub.i2, z.sub.i), i{1, . . . , m}. More specifically, (, , ) is the solution of the system

(157) $\begin{array}{cc}\left[\begin{array}{ccc}m& .Math.e,a_1.Math.& .Math.e,a_1.Math.\\ .Math.e,a_1.Math.& {.Math.a_1.Math.}^2& .Math.a_1,a_2.Math.\\ .Math.e,a_2.Math.& .Math.a_1,a_2.Math.& {.Math.a_2.Math.}^2\end{array}\right]\left[\begin{array}{c}\\ \\ \end{array}\right]=\left[\begin{array}{c}.Math.e,z.Math.\\ .Math.a_1,z.Math.\\ .Math.a_2,z.Math.\end{array}\right].& \left(71\right)\end{array}$
where e=(1, 1, . . . , 1), a.sub.1=(a.sub.11, a.sub.21, . . . , a.sub.m1), and a.sub.2=(a.sub.12 a.sub.22, . . . , a.sub.m2).

(158) To show this, first, F is a convex set. Let x=(x.sub.1, . . . , x.sub.m)=P.sub.FZ, then x minimizes
xz.sup.2=.sub.i=1.sup.m(x.sub.iz.sub.i).sup.2(72)
subject to the constraint that {(a.sub.i1, a.sub.i2, x.sub.i)}.sub.i{1, . . . , m} lie on the same plane in R.sup.3. Thus, it is equivalent to find the least squares plane (t.sub.1, t.sub.2)=+t.sub.1+t.sub.2 that fits {(a.sub.i1, a.sub.i2, z.sub.i)}.sub.i{1, . . . , m}. One defines e, a.sub.1, a.sub.2 as above, then (, , ) is the solution of (71). Since {(a.sub.i1, a.sub.i2)}.sub.i{1, . . . , m} are not on the same line, (71) has a unique solution (, , ). Finally, one computes x.sub.i=+a.sub.i1+a.sub.i2 for i{1, . . . , m}.
Projections onto General Surface Slope Constraints

(159) 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.

(160) Normal Vector and Surface Slope Constraint

(161) Let (e.sub.1, e.sub.2, e.sub.3) be the standard basis of R.sup.3. Given three points P.sub.1=(P.sub.11, P.sub.12, h.sub.1), P.sub.2=(P.sub.21, p.sub.22, h.sub.2), and P.sub.3=(p.sub.31, p.sub.32, h.sub.3) in R.sup.3, the normal vector to the plane P.sub.1P.sub.2P.sub.3 is the cross product

(162) $\begin{array}{cc}\stackrel{.fwdarw.}{n}=\left[\begin{array}{c}{p}_{11}-p_{31}\\ p_{12}-p_{32}\\ h_1-h_3\end{array}\right]\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]\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.sub.1, P.sub.2, P.sub.3 on xy-plane are (p.sub.11, p.sub.12), (p.sub.21, p.sub.22), (p.sub.31, p.sub.32), which are not on the same line. This implies a.sub.1b.sub.2a.sub.2b.sub.10. So one can rescale and use

(163) $\begin{array}{cc}\stackrel{.fwdarw.}{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

(164) 0$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},$
then
t.sub.1=a.sub.1u+b.sub.1v,t.sub.2=a.sub.2u+b.sub.2v, 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 form
g(u,v)0.(76)
An important case of such constraints is the surface orientation constraint below.

(165) Definition 1 (surface orientation constraint) Let be a triangle with normal vector {right arrow over (n)}=(u,v,1)R.sup.3 as above. Let {right arrow over (q)}=(q.sub.1, q.sub.2, q.sub.3)R.sup.3 \{0} be a unit vector and be an angle in [0, ], the constraint

(166) $\begin{array}{cc}\left(\stackrel{.fwdarw.}{n},\stackrel{.fwdarw.}{q}\right),\mathrm{orequivalently},\cos \left(\stackrel{.fwdarw.}{n},\stackrel{.fwdarw.}{q}\right)=\frac{-{\mathrm{uq}}_1-{\mathrm{vq}}_2+q_3}{\sqrt{u^2+v^2+1}}\cos ,& \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.

(167) Projection onto a General Surface Slope Constraint

(168) Consider a single triangle determined by three points Q.sub.1=(p.sub.11, p.sub.12, w.sub.1), Q.sub.2=(p.sub.21, p.sub.22, w.sub.2), and Q.sub.3=(p.sub.31, p.sub.32, w.sub.3) in R.sup.3. Without loss of generality, one may assume
w.sub.1+w.sub.2+w.sub.3=0.(78)
Projecting onto an orientation constraint is to find the new heights h.sub.1, h.sub.2, h.sub.3 that minimizes
:=(h.sub.1w.sub.1).sup.2+(h.sub.2w.sub.2).sup.2+(h.sub.3w.sub.3).sup.2(79)
such that the new points P.sub.1=(p.sub.11, p.sub.12, h.sub.1), P.sub.2=(p.sub.11, p.sub.12, h.sub.2), and P.sub.3=(p.sub.11, p.sub.12, h.sub.3) satisfy the given orientation constraint. Define t.sub.1:=h.sub.1h.sub.3 and t.sub.2:=h.sub.2h.sub.3, then
=(t.sub.1+h.sub.3w.sub.1).sup.2+(t.sub.2+h.sub.3w.sub.2).sup.2+(h.sub.3w.sub.3).sup.2.(80)
Using the change of variables (75), one has:
=(a.sub.1u+b.sub.1v+h.sub.3w.sub.1).sup.2+(a.sub.2u+b.sub.2v+h.sub.3w.sub.2).sup.2+(h.sub.3w.sub.3).sup.2.(81)
The new triangle P.sub.1P.sub.2P.sub.3 needs to satisfy some given orientation constraint g(u, v)0. So the projection reduces to

(169) $\begin{array}{cc}\underset{\left(u,v,h_3\right)R^3}{\min }={\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)0.& \left(83\right)\end{array}$

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

(171) $\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}$

(172) One may also define

(173) $\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

(174) $\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

(175) $\begin{array}{cc}\underset{\left(\hat{u},\hat{v},h_3\right)R^3}{\min }={\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)0,& \left(88\right)\end{array}$
Therefore, one can treat (87)-(88) as (82)-(83) with new respective coefficients .sub.i, {circumflex over (b)}.sub.i

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

(177) $\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

(178) $\begin{array}{cc}9={\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 defining
A=2(a.sub.1.sup.2+a.sub.2.sup.2a.sub.1a.sub.2)>0,(97)
B=2(b.sub.1.sup.2+b.sub.2.sup.2b.sub.1b.sub.2)>0,(98)
C=2a.sub.1b.sub.1+2a.sub.2b.sub.2a.sub.1b.sub.2a.sub.2b.sub.1,(99)
w.sub.a=3(a.sub.1w.sub.1+a.sub.2w.sub.2),(100)
w.sub.b=3(b.sub.1w.sub.1+b.sub.2w.sub.2),(101)
one has 9=3Au.sup.2+3Bv.sup.2+6Cuv6w.sub.au6w.sub.bv+Z. Thus, (82)-(83) reduces to

(179) $\begin{array}{cc}\underset{\left(u,v\right)R^2}{\min }\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)0.& \left(103\right)\end{array}$
As long as one can find the solution (u, v), one can obtain the projection (h.sub.1, h.sub.2, h.sub.3) from (90), (75), and (73). More specifically,

(180) 0$\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.sub.i, b.sub.i, u, v) is replaced by (.sub.i, {circumflex over (b)}.sub.i, , {circumflex over (v)}).
Projection onto Surface Maximum Slope Constraint

(181) Let P.sub.1P.sub.2P.sub.3 be a triangle in R.sup.3 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.sub.1P.sub.2P.sub.3 represents the desired ground, that cannot be too steep with respect to gravity, i.e., the slope of the plane P.sub.1P.sub.2P.sub.3 must not exceed a threshold s:=s.sub.maxR.sub.+. Then the angle between {circumflex over (n)} and e.sub.3=(0,0,1) must satisfy ({circumflex over (n)},e.sub.3)tan.sup.1 (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.sub.3)tan.sup.1 (s) is equivalent to
cos ({right arrow over (n)},e.sub.3)={right arrow over (n)},e.sub.3/{right arrow over (n)}cos(tan.sup.1(s))=1/{square root over (1+s.sup.2)}.(107)

(182) 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.sub.3, tan.sup.1 (s)). From (107), one may conclude

(183) $\begin{array}{cc}\frac{1}{\sqrt{u^2+v^2+1}}\frac{1}{\sqrt{1+s^2}}{u}^2+v^2-s^{2}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

(184) $\begin{array}{cc}\underset{\left(u,v\right)R^2}{\min }=\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}0,& \left(110\right)\end{array}$
Rescaling by

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

(186) $\begin{array}{cc}\underset{\left(u,v\right)R^2}{\min }\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-10.& \left(112\right)\end{array}$

(187) 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.sup.1 (s) (i.e. the normal vector lays in the cone 1502 depicted in FIG. 15), and that equation (79) is minimized.

(188) 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:

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

(190) 2. For each value 0, the level set (u, v)= is an ellipse in R.sup.2 whose center is the vertex (u.sub.0, v.sub.0) of (u, v), which satisfies
Au.sub.0+Cv.sub.0=w.sub.a(113)
Cu.sub.0+Bv.sub.0=w.sub.b.(114)

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

(192) One may consider two cases:

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

(194) Case 2: g(u.sub.0, v.sub.0)>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

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

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

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

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

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

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

(201) Using the Cauchy-Schwarz inequality, one has
AB=[a.sub.1.sup.2+a.sub.2.sup.2+(a.sub.1a.sub.2).sup.2][b.sub.1.sup.2+b.sub.2.sup.2+(b.sub.1b.sub.2).sup.2](120)
[a.sub.1b.sub.1+a.sub.2b.sub.2+(a.sub.1a.sub.2)(b.sub.1b.sub.2)].sup.2=C.sup.2.(121)

(202) Since A, B, >0, the determinant of (118)-(119) is
(A+)(B+)C.sup.2>ABC.sup.2>0.(122)
This implies the system (118)-(119) has a unique solution

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

(204) Next, substituting u and v into (117) yields

(205) $\begin{array}{cc}{\left[{}^2+\left(A+B\right)+\mathrm{AB}-C^2\right]}^2& \left(125\right)\\ ={\left[w_a+\left(w_{a}B-w_{b}C\right)\right]}^2+{\left[w_b+\left(w_{b}A-w_{a}C\right)\right]}^2& \left(126\right)\\ =\left(w_a^2+w_b^2\right){}^2+2\left(w_a^{2}B+w_b^{2}A-2w_a{w}_{b}C\right)& \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}$

(206) Defining the constants accordingly, one may rewrite as
(.sup.2+R.sub.1+R.sub.2).sup.2=R.sub.3.sup.2+2R.sub.4+R.sub.5.(129)

(207) 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
(.sup.2+R.sub.1+R.sub.2+y).sup.2=R.sub.3.sup.2+2R.sub.4+R.sub.5+2(.sup.2+R.sub.1+R.sub.2)y+y.sup.2(130)
=(R.sub.3+2y).sup.2+2(R.sub.4+R.sub.1y)+R.sub.5+2R.sub.2y+y.sup.2.(131)

(208) 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.sub.4+R.sub.1y].sup.2(R.sub.3+2y)(R.sub.5+2R.sub.2y+y.sup.2)=0(132)
2y.sup.3+[R.sub.1.sup.2R.sub.34R.sub.2]y.sup.2+[2R.sub.1R.sub.42R.sub.2R.sub.32R.sub.5]y+R.sub.4.sup.2R.sub.3R.sub.5=0. (133)

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

(210) $\begin{array}{cc}{\left({}^2+R_1+R_2+y_0\right)}^2=\left(R_3+2y_0\right){\left(+\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 :

(211) Case 2a: R.sub.3+2y.sub.0<0. Then

(212) $=-\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.

(213) Case 2b. R.sub.3+.sup.2y.sub.0=0. Then

(214) 0$y_0=-\frac{R_3}{2}$
and
.sup.2+R.sub.1+R.sub.2+y.sub.0=0.(135)
So one has

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

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

(217) Case 2c: R.sub.3+2y.sub.0>0. Define r:={square root over (R.sub.3+2y.sub.0)}. From (134),

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

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

(220) Since there is only one positive , one must have R.sub.4+R.sub.1y.sub.00. Otherwise, (140) consists of two quadratic equations that have constant term R.sub.2+y.sub.0; 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 rr.Math.sgn(R.sub.4+R.sub.1y.sub.0) and takes only the equation

(221) $\begin{array}{cc}{}^2+\left(R_1-r\right)+\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

(222) $\begin{array}{cc}=\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}$

(223) 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.sub.1, h.sub.2, h.sub.3).

(224) Projection onto Surface Oriented Minimum Slope Constraint

(225) Nonconvex Surface Minimum Constraint

(226) Let P.sub.1P.sub.2P.sub.3 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 .

(227) This happens, for example, if P.sub.1P.sub.2P.sub.3 must have a slope at least s:=s.sub.min R.sub.+. In this case, the angle between {right arrow over (n)} and e.sub.3=(0,0,1) satisfies
({right arrow over (n)},e.sub.3)tan.sup.1(s).(143)

(228) 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.sub.3)tan.sup.1 (s) is equivalent to

(229) $\begin{array}{cc}\cos \left(\stackrel{->}{n},e_3\right)=\frac{.Math.\stackrel{->}{n},e_3.Math.}{.Math.\stackrel{->}{n}.Math.}\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

(230) $\begin{array}{cc}\frac{1}{\sqrt{u^2+v^2+1}}\frac{1}{\sqrt{1+s^2}}{u}^2+v^2-s^{2}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).

(231) The Surface Oriented Minimum Slope Constraint

(232) (143) implies the angle between {right arrow over (n)} and the xy-plane satisfies
({right arrow over (n)},xyplane)(/2)tan.sup.1(s).(146)
In some cases, it is not unreasonable to align the plane P.sub.1P.sub.2P.sub.3 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.sub.1P.sub.2P.sub.3 to incline toward a unit direction {right arrow over (q)}=(q.sub.1, q.sub.2, 0)R.sup.3 (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

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

(234) 0$\begin{array}{cc}\cos \left(\stackrel{->}{n},\stackrel{->}{q}\right)=\frac{.Math.\stackrel{->}{n},\stackrel{->}{q}.Math.}{.Math.\stackrel{->}{n}.Math.}\cos \left(\frac{}{2}-{\tan }^{-1}\left(s\right)\right)=\frac{s}{\sqrt{1+s^2}}& \left(147\right)\end{array}$
and arrive at the following definition.

(235) 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

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

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

(238) $\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}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.

(239) Projection onto Surface Oriented Minimum Slope Constraint

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

(241) $\begin{array}{cc}\underset{\left(u,v\right)R^2}{\min }\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}$

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

(243) $Q:=\left[\begin{array}{cc}{q}_1& q_2\\ q_2& -q_1\end{array}\right],$
then Q.sup.2=Id since{right arrow over (q)}=q.sub.1.sup.2+q.sub.2.sup.2=1. So Q=Q.sup.T=Q.sup.1. Next, the variables are changed

(244) $\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]\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

(245) $\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+10.& \left(154\right)\end{array}$
Thus, the constraint (149)-(150) becomes

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

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

(248) $\begin{array}{cc}\left[\begin{array}{cc}{a}_1& a_2\\ b_1& b_2\end{array}\right]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

(249) 0$\begin{array}{cc}\underset{\left(\hat{u},\hat{v}\right)R^2}{\min }\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}+10\end{array}.& \left(159\right)\end{array}$
where A, B, C, w.sub.a, w.sub.b are defined below using the new a.sub.1, a.sub.2, b.sub.1, b.sub.2 from (157)
A=2(a.sub.1.sup.2a.sub.2.sup.2a.sub.1a.sub.2)>0,(160)
B=2(b.sub.1.sup.2+b.sub.2.sup.2b.sub.1b.sub.2)>0,(161)
C=2a.sub.1b.sub.1+2a.sub.2b.sub.2a.sub.1b.sub.2a.sub.2b.sub.1,(162)
w.sub.a=3(a.sub.1w.sub.1+a.sub.2w.sub.2),(163)
w.sub.b=3(b.sub.1w.sub.1+b.sub.2w.sub.2).(164)
Changing variables by

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

(251) $\begin{array}{cc}\underset{\left(\hat{u},\hat{v}\right)R^2}{\min }\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+10\end{array}.& \left(166\right)\end{array}$

(252) 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.sup.1(s) and the angle between the normal vector and the directional vector q is at most

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

(254) Before solving the above problem, one can observe that

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

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

(257) 3. For each value 0, the level set (u, v)= is an ellipse in R.sup.2 whose center is the vertex (u.sub.0, v.sub.0) of (u, v), which satisfies
Au.sub.0+Cv.sub.0=w.sub.a(167)
Cu.sub.0+Bv.sub.0=w.sub.b.(168)

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

(259) Case 1: (u.sub.0, v.sub.0) is feasible, i.e., it lies inside or on the left branch {(u, v)R.sup.2|g.sub.2(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.

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

(261) $\begin{array}{cc}{}_u\left(u,v\right)+\frac{}{2}{}_u{g}_2\left(u,v\right)=0,& \left(169\right)\\ {}_v\left(u,v\right)+\frac{}{2}{}_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}$

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

(263) 1. If g.sub.2(u.sub.0, v.sub.0)>0, then (169)-(171) has a unique solution (u.sub.1, v.sub.1, .sub.1) with .sub.1>0 and u.sub.1<0, and a unique solution (u.sub.2, v.sub.2, .sub.2) with .sub.2>0 and u.sub.2>0.

(264) 2. If g.sub.2(u.sub.0, v.sub.0)0 and g.sub.1(u.sub.0, v.sub.0)=u>0, then (169)-(171) has a unique solution (u.sub.1, v.sub.1, .sub.1) with .sub.1>0. Thus, u.sub.1<0. In addition, (169)-(171) also has at least one solution (u.sub.2, v.sub.2, .sub.2) with .sub.2<0.

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

(266) 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.sup.2u.sup.2+1=0 and an elliptical level set of (u, v).

(267) 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.

(268) 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 .

(269) 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.

(270) Next, one can show how to find the solution in Case 2. First, write (169)-(170) as
(A)u+Cv=w.sub.a,(173)
Cu+(B+)v=w.sub.b.(174)
Suppose for now that the determinant (A)(B+)C.sup.20. Then

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

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