- Email: [email protected]
Visibility of surfaces via differential geometry
353
Computer Aided Geometric Design 7 (1990) 353-373 North-Holland
Visibility of surfaces via differential geometry Rolf WALTER Institut
ftir Mathematik,
Universitiit Dortmund,
D-4600 Dortmund
50, FRG
Presented at Oberwolfach October 1988 and April 1989 Received June 1989; revised November 1989 Abstract. Using the means of differential geometry, we develop a theory for rendering surfaces with complete visibility claryfying. The results are based on the new concept of ‘sight indices’ which control the change of coverings dynamically along arbitrary given curves. The method is entirely curve oriented, so suitable for most graphic outputs, in particular for plotters. Important attributes of surfaces, like silhouette, selfintersection, and boundary curves will be generated and, in fact, be viewed as constitutive for the whole process. Keywords. Surface rendering, visibility, hidden lines, silhouettes, surface intersections, boundary curves
1. Introduction The handling of surfaces including their attributes like boundaries, silhouettes and intersections is a central subject of CAGD. In this paper, we shall attack these questions using the means of differential geometry. For the present, we will concentrate on surface rendering and visibility clarifying with special attention paid to these attributes. Geometric surface design has developed in various branches, especially in the very successful field of approximation theory. We refer to the excellent survey articles of [Sutherland, Sproull & Schumaker ‘741 and [Boehm, Farin & Kahmann ‘841, the book of [Newman & Sproull ‘841, and the vast literature given there. However, the visibility problem is still a challenge. The traditional procedures start, in one way or the other, with discretizations of the surface at an early stage. The present method may be considered as an alternative and supplement to this, in particular when the surfaces come in complicated analytic representations. Its main feature is that it does not use any patches or pixels at all: it is completely curve oriented. The important additives of a surface like boundaries, silhouettes and selfintersections will not only be generated but considered as consstitutive for the whole process. Additional curves, like parameter lines, niveau lines, geodesics, curvature lines, etc. can be added at any time. The visible (or covered) parts of curves will be drawn in one stroke which is especially important for controlling plotters. The growth of data with respect to refinements is roughly linear compared to the quadratic ascent of patches or pixels. No sorting of big fields is needed. On the other hand, due to the complexity of differentialgeometric surfaces, one encounters situations which are not present in finite straight models. So, for higher singularities, further investigations may be necessary and seem possible along similar lines. 2. The main ideas Visibility is a global problem, because decisions have to be made for points far apart on the surface. It is in the spirit of differential geometry to relate global questions to suitable local or pointwise invariants. 0167-8396/90/$03.50
0 1990 - Elsevier Science Publishers B.V. (North-Holland)
R. Walter / Visibility of surfaces
354
Figs. l-3.
Rendering procedure for a bounded piece of a rotational torus. In Fig. 2, the sight indices L$ and D4 refer to the silhouette, while the sight index D5 refers to the boundary as the current curve.
As to the projection scheme, the method here applies to all types as well (parallel or central). In principle, the observer rays have to be cut first with the object in 3-space and then with the image plane. For simplicity of description, we assume that the object is a single surface and consider only curves on this surface. But these restrictions are not really necessary. The observer rays shall be called rays, for short. The first idea is to avoid the calculation of the intersecting points of the rays with the surface but rather to keep a balance of the number of these points along a current curve under investigation. The number of points on a ray, common with the surface and before the current point towards the observer, will be called the sight number. So a sight number of zero means that the current point on the surface is visible. Given a certain projection scheme, the sight number is a globally defined function on the surface. It is constant during transversal intersections and away from boundaries and selfintersections. Changes of the sight number along a current curve happen only occasionally, but these jumps are decisive for our visibility method. We call each change of this type a sight index. The salient point is that these sight indices can be computed from local differential-geometric quantities like the second fundamental form of the surface. The derivation of such index formulas for the most important events is the main theoretical concern here. The actual procedure in the computer is then as follows: (A) Generate the distinguished curves, namely the silhouettes, the selfintersections and eventually the boundaries as point sequences. The mathematical principles for this part are essentially known. We will describe a specific arrangement which has proven very effective. The distinguished curves are generated in full, no sight indices enter this part (see Fig. 1). (B) Clarify the visibility of any surface curve according to the following rules: (B.l) If the sight number is known for one point of a given curve, then it is known for each point on the curve by calculation of the sight indices and the places where these occur. So the curves on the surface do the job of ‘transporting’ the sight number from one surface point to another. (B.2) Starting with one visible point on the surface (or with one point whose sight number is known), the visibility of each surface curve can be clarified individually, together with its plotting. This is in particular possible for each of the distinguished curves like the silhouettes, the selfintersections and the boundaries (see Fig. 2, where also some sight indices are marked, and Fig. 3).
3. Differential-geometric The most general mapping
foundations
surface here will not only be viewed
as a point
set in R3 but rather
x: M+R3, where
M is a two-dimensional,
as a
(3.1) compact,
oriented,
and connected
Hausdorff
manifold
with
R. Walter / Visibility of surfaces
355
boundary aM, and x is immersive, i.e. of rank 2. For simplicity, the differentiability is assumed to be C”, but this can be relaxed considerably, mostly to C* or C3, depending on the problem. For many questions, M may be viewed as a subset of Iw2 and x as a parametrized surface. However, the concept of M as a manifold does not stick to a specific parametrization and also comprehends more complicated topological types. In any case, the set M will be addressed to as the domain of the surface. In the following, square brackets [ ] denote the determinant, angle brackets ( ) denote the scalar product, and the lying cross X denotes the vector product in Iw3. The field of unit normal vectors corresponding to the given orientation is viewed as a mapping: N: M+R3.
(3.2)
The first and second fundamental I := (dx, II:= (dx,
form of the surface
are:
dx),
(3.3)
dN).
(3.4)
The evaluation of such forms, e.g. of ZZ, in p E M for two tangent vectors p, v of M in p will be denoted in full by II( p; p, v) := ((dx),(p), (dN),(v)) or, for short, by JZ(p, v). By the theorema egregium, the expression for the Gaussian curvature is: (3.5) where CL,v are two linearly independent vectors at the domain point p considered. References for further details on differential geometry are e.g. [Hicks ‘651 and [Walter ‘891. Among the properties of the projection scheme, only the directions of the observer rays along the surface play a role in our index formulas. These direction vectors may be described for parallel or central projection in a common way, as follows. In the first case, let a be the constant direction vector of the rays towards the observer, and in the second case, let a be the constant position vector of the center. The rays through the surface points x then have the direction vector L=a-8.x,
(3.6)
where 6 = 0 in case of parallel projection and S = 1 in case of central projection. The observer ray through the surface point x(p) is thus given by the parameter tion t++x(p)+t*L(p),
tER.
representa-
(3.7)
One of its points, say x(p) + t1 . L(p), lies before another, say x(p) + t, . L(p), if t, > t,. As the general symbol for the projection mapping of Iw3 (or lR3\ { a }, in case of central projection) to the image plane B, we use the letter P. The explicit formulas for P are well known but we don’t have to use them in our theory. However, in case of central projection, we assume for simplicity that no ray through a surface point is parallel to the image plane, and that the center never lies on the surface nor between two surface points on the same ray. For any domain point p E M, x(p) is the corresponding object point in Iw3, and P(x( p)) is the corresponding image point in B. (As to the computer handling, there are three corresponding choices for the storage of points which should be used flexibly according to practical requirements.) We will utilize the sight function u defined by the scalar product: a:=(L,
N):
M-l!%.
(3.8)
356
R. Walter / Visibility of surfaces
The zeros of the sight function u yield those points on the surface whose rays are tangential to the surface there. These points are the so called silhouette points. The silhouette C is the following closed subset of M: c:=
{pEM~a(p)=O}.
We will decompose L = (dx)(
the direction
(3.9) vector L of rays in a tangent
and normal
part, according
X) + u. N.
to:
(3.10)
So X is the vector field on M, corresponding to the tangential part of L along the surface in space. (3.10) is evident, since L - u *N is orthogonal to N and thus is of the form (dx)( X). By forming the differential of (3.6) and (3.8) and insertion of (3.10), we obtain do=
(dL,
thus da = ((dx)(X), (de)(p)
N)+ dN)
(L,
dN)
=(
-6.dx,
or, for any p tangent
N) + ((dx)(X)
+u.N,
dN),
(3.11)
to M:
= LL(X, CL),
(3.12)
expressing the differential of the sight function by the second fundamental form. If the mapping x is not injective, then the surface has selfintersections. We define selfintersection (set) S by: S:={p~MIthereexistsq~Mwithpfqandx(p)=x(q)}.
the
(3.13)
It is a consequence of the local injectivity of the map x that S is a closed subset of M. A point of S will be called simple if x(p) has exactly two pre-images under x. Also, the boundary aM is known to form a closed part of M. To distinguish from the three subsets C, S, and EM4 of M are their images in R3 under x, which we shall supply with the additive ‘in the object space’, and the images of those in B under x, which we shall supply with the additive ‘in the image plane ‘. Thus, e.g., C is the silhouette, x(C) is the silhouette in the object space, and P(x(C)) is the silhouette in the image plane. All nine sets C, S, aM, x(C), x(S), x(aM), P(x(C)), P(x(S)), P(x(aM)) are compact. Remark. In classical descriptive geometry, the adjectives instead. So, e.g., x(C) is then called the true silhouette
‘true’, resp. ‘apparent ’ are often used and P(x(C)) the apparent silhouette.
Theorem 3.1. Assume, for a fixed p0 E M, that the ray through x( pO) does not hit the silhouette, the selfintersection, and the boundary in the object space before and including x( p,,); i.e. wheneuer, for some t, 2 0 and somep, E M, we haue x( pO) + t,L( p,,) = x( pl), then p1 6E C U S u 3M. Then the following holds: (i) There exist only finitely many intersection points of this ray with the surface x(M) before and including x( p,,), say n points (ii) There is a neighborhood U of p0 in M such that, for any p E U, the ray through x(p) has exactly n intersection points with the surface x(M) before and including x(p), and again does not hit the silhouette, the selfintersection and the boundary in the object space before and including x(P)* Proof. First, an intersecting point p1 of the given fixed ray with x(M) is isolated on this ray, because otherwise an easy limit argument would show that L( pI) is tangent to the surface at x( pl). Since the intersection of the fixed ray with x(M) is compact, there are only finitely many such pl. Now observe that, generally, two surface points x(q,) and x(q,) are situated on the same ray if and only if q1 and q2 have the same image in the image plane B under the composition P 0 x. And: the points where this map P 0 x: M + B has a vanishing Jacobian are exactly those of the silhouette C.
R. Walter / Visibility
of surfaces
357
This reduces the question, how many points q of M have images x(q) on a ray, to the problem of the cardinality of pre-images of a given fixed point in B under P 0x.There are well known arguments showing that, if the Jacobian does not vanish, this number is constant in a neighborhood of a point where it is finite (cf. [Walter ‘681, e.g.). A little care has to be taken to stay outside C, S, and EM, which is possible by the compactness of these sets. Also, one has to focus on points before x(p), but this does not disturb the local constant argument. So, for p in this neighborhood, any intersection point of the ray through x(p) with x(M) before and including x(p) has only one pre-image in M under x. 0 The general definition of the sight number will be put in the following way: For any p E M, it is the cardinality of the intersection of the halfray t - x(p) + t L(p), t > 0, with the point set x(M) in R3. Remark. Alternatively, one could include the number of pre-images of any such intersection point under x in case that the intersection point is on the selfintersection set x(S) in the object space. But we will not do so. So, the meaning of Theorem 3.1 is that the sight number is finite and constant around any p. of M such that the ray through x( po) does not hit x(C) U x(S) U x( CIM) before and including x( PO). As usual, for local calculations, representation of the form
the mapping x may be replaced by a Gaussian parameter
X1(% u> x:A*R3,
(u, U)_X(U,
u)=
x2@,
u)
(3.14)
,
i X3(% u) 1 where the domain A is now contained in lR*. Any domain point p is then represented by its Gaussian coordinates u, u [within the validity of (3.14)]. We presuppose a little familiarity in passing over from (3.1) to (3.14). [The formal scheme is that the x in (3.14) is in fact the composition of the inverse of a chart of M with the x in (3.1).] The field of unit normal vectors is then correspondingly represented by:
N: A-+R3,
(u, u) *N(u,
u) =
N’(u,
u)
NZ(u,
u)
1 N3(U,
.
(3.15)
u> 1
The orientation of the Gaussian coordinates is always such that [XU, X”,
Nl ’ 0.
(3.16)
Here and in the sequel, indices like U, u denote partial derivatives. In later parts, the visibility shall always be clarified ‘dynamically’ for points running along a given surface curve under consideration, the so called current curue. Any surface curve can locally be described by a parameter representation of the form: w c, x(r+),
U(W))>
(3.17.a)
where the curve parameter w runs through a real interval. The point c(w) := (u(w),
U(W))
(3.17.b)
describes the corresponding curve in the parameter domain A. Its component functions U(W) and u(w) are assumed to be of class C” (with similar possibilities for relaxing as for x). The
358
R. Walter /
restriction of any surface quantity, defined (3.17) will be denoted by a bar. Thus, e.g. x(w)
:=x(c(w)),
Visibility
ofsurfaces
on the domain
L(w):=L(c(w)),
If (3.17) describes the current curve then, sometimes, and the ray through X(w) as the current ray.
A, to the points
of a surface
etc. X(w) is addressed
curve (3.18)
to as the current point,
4. The sight indices The sight number, which is globally defined on M or on x(M), will now be investigated along a given surface curve (3.17) the so called current curve. We call a surface curve distinguished if all its points c(w) lie on one of the sets: silhouette, selfintersection or boundary, and we call it general is c(w) hits those sets only at isolated points of the parameter interval. So, we will have to consider the cases where the current curve is a general one or a distinguished one (‘mixed’ cases are excluded here). We assume the parameter interval of the current curve to be open. Of course, there is no problem at a fixed parameter point w0 of the current curve when the sight number is constant around w,, i.e. in an open interval containing w,. In the contrary case, the sight index at w,, of the current curve exists, by definition, if the sight number on the surface points X(w) has a left and right (proper) limit N,, N_ at w,, and its value is then the difference N+- N_ of these limits. Observe that this implies that the sight number is constant near left and near right of w,,. If the sight number is a piecewise constant function along the current curve, then the sight numbers in the constancy intervals can be calculated from the sight indices at the jump points and a single value of the sight number in one of the constancy intervals. (The sight numbers at the jump points are practically of no interest.) By Theorem 3.1, the sight number changes at most at those points w, of the current curve where the ray through x( wO) hits one of the three sets x(C), x(S), x(aM) before and including x(wa). Any circumstance which causes a change of the sight number along the current curve will be called an euent, for short. For the most important events, we shall calculate the sight index under certain genericity assumptions. General assumption. For simplicity of formulation, we describe each event individually, assuming that no other events take place at the point considered on the ray. In practice, if more than one event is happening at the same parameter value w, but at most one at each point on the corresponding ray and each of them obeys the genericity assumptions individually, then the sight index at w, is gained by addition of the individual ones at w,. There are two different types of events possible because either x( wa) itself may meet the sets x(C), x(S), x( aM> or else a point on the ray through x( w,,) before x( wO) may meet those sets. We say that the event is of point type, resp. of ray type. In the arrangement of results, we follow the degree of complexity of the surface; so we start with the theorems which would suffice in case S=fl, aM=,0. The normal case Theorem 4.1. For w = w,, let the current curve w ++X(w) = x(c( w)) meet the silhouette object space, i.e. a(c(w,,)) = 0. Then the sight index for this event is given by:
~q+o>; 4wo), mwm
D1= -SignIl(c(wo); X(c(w,)),
X(c(wo))) ’
in the
(4.1)
359
R. Walter / Visibility of surfaces
where one has to assume, nonzero.
in addition,
that the numerator
and the denominator
in (4.1)
are
The general assumption implies here that X(wO) is not contained in x( S U aikf ), and similarly in the following theorems. By (3.12), the last assumptions mean intuitively that neither the current curve not the sight vector ‘touch’ the silhouette in the object space at the point considered. Theorem 4.2. For w = w,, let the ray through the point X(wO) of the current curve w c* X(w) = x( c( w)) hit the silhouette in the object space before X( wO), say at x( pl) with p, E C. Then the sight index of this event is given by
(~‘(WIA N D2= -2siWIl(p,; where one has to assume, nonzero.
h(p,), in addition,
PI)> h(p,))’ that the numerator
(4.2) and the denominator
in (4.2)
are
The last assumptions mean intuitively that neither the ray through X(wO) is tangent to the silhouette in the object space at x( p,) nor the ruled surface formed by the rays through the current points X(w) touches the surface x(M) at x( p,). The silhouette
itself
Now, consider the case that our current curve is fully contained in the silhouette. This implies that the numerator in (4.1) is identically 0. Genericly, the denominator in (4.1) will only vanish at isolated parameter values (our assumptions shall guarantee this and, moreover, that the current curve locally exhausts C). If this denominator vanishes at w = w,, then the ray through X( wO) touches the silhouette in the object space at X(wO). This is an interesting phenomenon which is responsible for the seemingly spontaneous stopping of the silhouette in the image plane, known e.g. from the torus (see Fig. 2 at D3 and also the torus picture in Fig. 4).
In the following theorem, there enters the covariant derivative of the second fundamental form which is a Codazzi tensor and therefore will be abbreviated by Cod := 011. For differentialgeometric
(4.3) details, see e.g. [Walter ‘851.
Theorem 4.3. Let the current curve w * c(w) u( c( w)) = 0,
be contained
in the silhouette,
identically.
Fig. 4. Rotational
i.e (4.4)
torus with a torus knot as current
curve.
R. Walter / Visibility of surfaces
360
If, for w = w,, II( c( wO); A( c( wO)), A( c( wO))) # 0, then no point type sight index occurs at w,. If II( c( wO); A( c( w,,)), A( c( wO))) = 0, then this causes a point type sight index of value D3 = -sign(X’(w,,),z(w,)), where, in addition,
(4.5)
it is assumed that
Rk(w0))
(4.6)
f 0,
Cod(c(w,);
h(c(wo)),
+(wo)),
+(w,)))
~0.
(4.7)
Under these assumptions, the current curve is locally open as a map w - c(w) into C, and it can be reparumetrized near w, in such a way that X’(wO) = p - z(w,,) for some p # 0. The boundary case For the discussion of the boundary behaviour, one should note first that the immersion x : M --$ R 3 can be continued locally over any boundary point as an immersion. Let R, resp. E, denote the field of unit tangent, resp. normal vectors of the boundary along the restriction x 1i3M, viewed as mappings of aA4 into R 3. Both fields are tangent to x(M), i.e. (R, N) = (E, N) = 0 on M. The orientation of R is chosen canonically (‘interior to the left’), and that of E according to [E, R, N] > 0. Of course, there is no sight index at w, of point type if c( wO)E aM\ S. Theorem 4.4. For w = w,,, let the ray through the point X( wO) of the current curue w - X(w) = x(c(w)) hit the boundary in the object space before X(wO), say at x( p,) with p1 E 3M. Then the sight index of this event is given by
D = sign[+%)~ w%)~ NPdl 4 m%),
where one has to assume, nonzero.
in addition,
(4.8)
’
NP,))
that the numerator
and the denominator
in (4.8)
are
Now, consider the case that the current curve is a piece of the boundary. If the ray through T(w,,) hits the sets x(C), x(S), and x(&M) before X( wO), then there is no difference to the same question for a general surface curve. Theorem 4.5. Let the current curve w ++c(w) c(w)
E aM,
D5= -sipIl(c(wo); h(c(w,)),
The selfintersection
i.e
identically.
If, for w = w,, this curve meets the silhouette, is given by:
where one has to assume, are nonzero and (E(w,),
describe a piece of the boundary,
(4.9) i.e. a( c( wO)) = 0, then the sight index for this event
1 --skn(E(2w,), z(w,)) X(c(w,)))
in addition, that the numerator z(w,)) z 0.
and the denominator
(4.10)
of the first quotient
case
Here the sight indices will be considered which are connected with a transversal selfintersection whose points are simple, i.e. have exactly two pre-images. So, for x( pl) = x( p2), p1 # p2,
R. Walter / Visibilityof surfaces
361
we presuppose that N( pl), N( p2) are linearly independent. Then the selfintersection consists of curves. If the ray through the point X( wO)of a general current curve meets the selfintersection set in the object space before X(wO), then the corresponding sight index is zero. (Observe that our general assumption excludes the case that this ray point is also on the silhouette or on the boundary in the object space.) Since the cases where the ray meets the silhouette or the boundary in the object space before X(wO) are already covered by the foregoing indices, there remains one further important situation: Theorem 4.6. For w = w,,, let the current curve w ++X(w) = x(c(w)) meet a simple transversal selfintersection in the object space and let p1 be the pre-image of X( wO) under x, different from c(wO). Then the sight index for this event is given by:
06=
. -sign
e+i3)~
N(P& N(p,))
(L(w,),
where one has to assume, nonzero.
in addition,
(4.11)
’ that the numerator
and the denominator
in (4.11)
are
Along a current curve which is a transversal selfintersection, a change of the sight number happens at most when the ray through the current point X( wO) meets the set x(C U S U EM) before X( wO)or at X( w,,). Then the sight index has to be calculated as for a general curve where analogous assumptions have to be made. Exceptional
cases
There are exceptional cases which are not covered by the theorems above. Such exceptional situations may occur because some of the assumptions are violated or by a coincidence of more than one of the ‘pure’ events considered above. Normally, these exceptions are nongeneric, i.e. they vanish if the situation is disturbed slightly. For example, one drastic exception is that the silhouette contains inner points, i.e. whole two-dimensional parts of the surface have tangent sight vector. Nevertheless, this can be removed by a slight change of the projection parameters. In other, less drastic situations, there exists the possibility to derive substitutive index formulas. As a model for such a remedy, we give the following theorem which concerns the case of Theorem 4.1 when the denominator there vanishes: Theorem 4.7. For w = w,, let the current curve w ++X(w) = x( c(w)) meet the silhouette in the object space at a point where the sight vector is tangent to the silhouette, i.e. a(c(w,,)) = 0 and ZZ( c( w,,); h( c( w,,)), A( c( wO)) = 0. Then, instead of (4.1), the sight index for this event is given by:
Z~Wo>; +d~ WWJ)
D’= -SignCod(c(wO); where one has to assume, nonzero.
X(c(w,J,
in addition,
Another exceptional situation curve in R3. For this, an index instead of D,. The procedure for traced out roughly by the proofs complicated then).
A(c(w,,),
that the numerator
X(c(w,)))
’
and the denominator
(4.12) in (4.12)
are
would be, that the boundary itself has selfintersections as a formula could be derived similar to (4.10) where D6 enters establishing index formulas for such exceptional situations are to be given in Section 7 (though the details may become more
362
R. Walter/ Visibilityof surfaces
5. Generation of the distinguished curves The actual computation of the sight indices requires the intersection of the current ray with the silhouette, selfintersection and boundary in the object space. So the data of these three sets must be produced and kept available in the computer. For the present, assume that each of these sets, as viewed in the domain, consists of finitely many curves (components) which are either closed or run from boundary to boundary. Each such component of a distinguished curve can be represented in the computer in form of a finite sequence of points, a polygon. These sequences have to be fine and accurate enough, because they enter subsequent calculations. The boundary components play a somewhat separate role, because normally they are given explicitly by Gaussian parameter representations. Then, of course, a storage as point sequences would be excessive. If, however, the boundary is defined implicitly, e.g. by cutting out ‘windows’ of the surface, then the generation of such boundary components is subject to the same process as for the silhouette and the selfintersection components. Remark. In principle, ‘parametric’ manner manner in form of distinguished curves, as least temporarily,
there are always two possibilities for storing surface points, namely the in form of domain points p (as Gaussian coordinates) and the ‘spatial the points x(p) in object space (as suitable coordinates there). For the points should be stored parametrically, but it may be helpful to use them, in spatial form.
The numerical means for this part are essentially nent is obtained by solving a system of (nonlinear) degree 1, i.e. of the form
known because any distinguished compoequations, which is underdetermined with
fi(UI,..., 24,)=o )...) fn-l(ul )..., UJ = 0.
(5.1)
For the silhouette, we have n = 2, and (u,, u2) = (u, u) are the Gaussian surface coordinates, and fi is the sight function cr. For the intersection of two surfaces, we have n = 4, and (u,, z.+) = (u, u) are the Gaussian coordinates of the first surface, and (u,, uq) = (ii, 6) those of the second, while the f, are the components of the difference of the position vector functions of the two surfaces. As to the numerical solution of the system (5.1) we prefer a generalized Newton’s method which applies as well to all types of distinguished curves and has the advantage that it can also be used for initializing the solutions, i.e. for the global solution part. The local background will only be outlined briefly (without discussion of converges): We start with an arbitrary point a,) and try to find a solution point which is as near to a as possible, say a + h, u=(u,,..., h:=(h,,...,h,). The requirement fi(a+h)=O for iE{l,...,n-1) is, in first order, equivalent to: f,(a)
+ t
f,&)h,
= 07
iE {l,...,
(5 4
n - l},
j=l
where the commas in the index indicate partial derivatives. Among system (5.2) we are searching one for which the squared norm g:=$
i:
(Q2
the solutions
h of the linear
(5.3)
j=l
becomes
minimal,
that
is we like to come near a solution with a small change of a. The of g under the side condition (5.2) leads to the following
Lugrungiun method for the extremum
R Walter / Visibility ofsurfaces
363
equations for the Lagrangian multipliers X,: (54
and for the associated extreme point: n-1 h,
=
c
fk,jWL.
(5.5)
k=l
If the system (5.1) has maximal rank n - 1 at a, then the coefficient matrix of (5.4) is of the Gram type and regular, and the unique inversion of (5.4) leads with (5.5) to an exact solution a# :=a+h
(5.6)
of (5.2) and thus to another approximative solution of (5.1) which usually is better than the original one. The iteration a, a#, a##, . . . yields, under favourable circumstances, an exact solution a * of (5.1) in form of a limit. In practice, it may be useful to replace the partial derivatives occuring in (5.2) by suitable difference quotients. Though this procedure does not exactly generate the solution point nearest to the arbitrary starting point, it is helpful to almost think in this terms. So, the assignment a e a * will be called a quasi-footpoint-projection, or QF-projection, for short. Of course, it may happen that the shows irregularities, nonconvergence or exit of the domain. This is sequence a, a#, a#*,... not as dramatic as it looks at first sight, because the starting point a is arbitrary and may be varied until the convergence becomes satisfactory. Anyway, a first solution point of (5.1) may thus be found. Then one can generate a second one, for example by QF-projection of a suitable neighboring point. Once there are two neighbouring solution points, then the chord of them can be prolonged over the newer point in a small increment, and QF-projection again produces from this an improved third solution. This is of iterative character, always producing from the last two points the next one. Thus, a whole solution component may be generated. This generation process of a solution component must be stopped at the boundary or if the component closes up. For the recognition of this second possibility, it proved useful to introduce a tube around an initial piece of the component and to test the entrance in that tube tangentially, i.e. using two or three successive points instead of one. (Tubes, or parallel sets, are familiar and very useful tools in geometry; see e.g. [Walter ‘811.) In order to generate all solution components, the above process “arbitrary starting point * solution component” should be applied on a sufficiently large set of starting points, which can be obtained by subdividing the domain after some suitable rule (e.g. by halfing). We would like to emphasize that ‘starting point’ here does not mean a first point of solution, but an arbitrary point on which the QF-projection may be applied to produce a first solution. This specific arrangement contributes to the efficiency of this part. After running over to the next starting point, one has to decide whether the QF-projection ‘lands’ on a solution component already generated. This again can be done by a tube construction.
6. Program structure The actual programming is not the purpose of the present paper. So we will indicate only the rough structure in order to show that the ideas, developed here, fit together to a reasonable algorithm. The whole method is composed of two main acts (A) and (B) which were described in Section 2. Act (A) is essentially numerical (see Section 5).
R. Walter / Visibility of surfaces
364
In the first part (B.l) of act (B) one moves along a given current curve in suitably fine steps, intersects the current ray in each step with the distinguished curves and at each event calculates the corresponding sight index, according to the theorems of Section 4. It is important that outside the current point, this intersection has actually to be done very rarely, because a ray whose trace point in the image plane does not meet the distinguished curves there, does not hit them in object space. Now the second part (B.2) of act (B) is performed as follows: With (B.l) at hand, and knowing the sight number of one specific point poo on the surface, the sight number along any current curve can be calculated. In fact, for one specific initial point on the given curve, one can ‘fetch’ the sight number by connecting it through any curve with poo, and then, parallel to the generation of the current curve points, one can update the sight indices by the procedure mentioned before. It is useful to do this for one specific pair of curves, which intersect in poo and to store the visibility data (sight indices and occuring places) for these two ‘base curves’. The initialization of the sight number for any other current curve may then be done by fetching it from the nearest point on the base curves. The role of the current curve can be taken by the distinguished ones or by any general curve, and there is no order, what has to be done first. Only the index formulas to be used depend on the type of the current curve. Since the sight number is a globally defined function on the surface, there are many possibilities to test the sight number of a point calculated by sight indices along curves: Any curve reaching the point considered must deliver the same sight number there! If the sight number for one specific point should not be known a priori, one may find a visible point by a (rough) preliminary search or, more sophisticated, by assigning a fictive sight number to some arbitrary point and renormalizing the lowest sight number found on the silhouette and boundary to zero. Surface
input
The standard way of input of the given surface is by Gaussian parameter representations and their derivatives up to the order 2. (As to the Codazzi tensor vZZ, only sign decisions at the zeros of the quadratic form of ZZ are involved which are numerically possible using ZZ alone.) In principle, any input equivalent to this or capable of producing this, is also possible. A conceivable goal could be to input just the coefficient functions of the derivation equations of Gauss and Weingarten, and to include their integration in some way, because the necessary derivatives then come out as by-products. Runtime acceleration
We mention some means to reduce the runtime: The intersection of current rays with the distinguished components occurs very rarely, as we saw above. A further reduction is possible by first doing the intersection for the projections in the image plane and by a sort of metric preliminary test: if a point has larger distance from the vertices of a polygon as is given by the maximal distance of two consecutive vertices, then the point is not on the polygon (and similar for the intersection of a straight interval with a polygon). Furthermore, the image plane may be divided in a small number of rectangles. This induces a corresponding division of the polygons which represent the projections of the distinguished components. Search through the whole polygon may then be replaced by search through some few of these subpolygons, because the others are too far away. As far as operations in the image plane are necessary, it is useful to replace (at least temporarily) Cartesian space coordinates by triplets whose first pair are coordinates in the image plane and whose third component is the distance from the projection
365
R. Walter / Visibility of surfaces
plane measured on the ray, the height. Of course, there are easy and quick transformations between these two types of space coordinates.
7. Proof of the index formulas In the whole proof we will use, in the neighborhood surface representation of the form F(z)
of a point zi of x(M),
= 0.
an implicit (7.1)
By this, the intersection problem of the rays with the surface may be viewed under the bifurcation aspect of a single scalar equation, depending on one additional parameter. Finally, the result can be re-interpreted for the Gaussian parameter representation of the surface. In (7.1), F is a real P-function of rank 1, defined on a spatial neighborhood U of zi such that x(V) n U is the zero set of F (for some open parameter domain I’). The insertion of the points of the current rays along a curve X(w) = x(c(w)) under consideration yields a real function of two variables: G(t,
w) := F@(w)
(7.2)
+ t.z(w)).
If the ray point X(w) + t. z(w) is contained in U, then the condition G(t, w) = 0 expresses that this ray point is on the surface part x(V). So, we shall have to study the zero behaviour of the real function t - G( t, w) in dependence on the parameter w. For this, the following two lemmas are basic. These lemmas refer, quite generally, to a real P-function (t, w) - G(t, w), defined in a neighborhood Q of a fixed point (t*, w*) of R*. Let the corresponding zero set be denoted by: w)~QandG(t,
Z(w):={t~Rl(t,
w)=O}.
(7.3)
Of course, the differentiability can again be relaxed considerably, as the proofs will show. Without loss of generality, neighborhoods in R shall be intervals, and for any neighborhood W of w,, let W+, resp. W- denote the left, resp. right part of W (including w, itself). Lemma 7.1. Assume:
G(t,,
w*> =0,
al := G,(t,,
we) # 0,
b:= G,(t,,
w*) # 0.
Then there exist neighborhoods T* oft, and W, of w, such that: (i) For any w E W,, Z(w) n T* consists of exactly one number t = f(w). (ii) f(w*) = t,. (iii) The function f is in W, of class C”, and strictly increasing if a,/b decreasing if a/b > 0.
(7.4) (7.5) (7.6)
< 0, and strictly
Lemma 7.2. Assume:
G(t,,
w*> =0,
(7.7)
G,(t,,
w*) = 0,
(7.8)
a2 := G,,(t*, b := G,,,(t,,
we) + 0, w*) # 0.
(7.9) (7.10)
366
R. Walter / Visibility of surfaces
Then there exist neighborhoodrr T* of t, and W, of w, such that, in case a/b < 0, we have: (i) For any w E Wi \ { we}, the set Z(w) n T* is empty. (ii) Z(w*) n T* = {t*}. (iii) For any w E Wz \ { w*}, the set Z(w) n T* consists of exactly two points on different sides oft*, and these points converge, for w + w,, (in W,‘), to t,. In case a,/b > 0 the same is true with interchanged roles of Wi and Wz. For the proof, it is very useful, to utilize generalized Taylor formulas, which expand along increasing powers of t - t,, but have main coefficients depending on the parameter w; cf. [Whitney ‘431. By this, we have for fixed p > 0:
5 $Gq(w)(t-t.)‘+ (p:l) ~p+l(t>w)(t - t*)p+l,
G(t, w)=
!
(7.11)
q=o
where VG
Gq(W)=F(‘*,
w)
(7.12)
forO
but where also the GP+i are of class C”. Proofs. Lemma 7.1 is a direct consequence of the standard theorem on implicit functions. We displayed it here because, in connection with certain factor decompositions, also higher cases can be dealt with it. Lemma 7.2 will be treated, assuming a2 > 0,
(7.13)
b
the other cases run analogously. First of all, the function t - G(t, w) is near (t*, w*) strictly convex: By a2 > 0, there are neighborhoods T* of t* and W, of w, such that G,,(t, w) > 0 in T* X W,. By (7.7) and (7.8), this implies the convexity. In the sequel, we shall have to reduce the neighborhoods T* and W, again. For 7.2 (i), (ii): In the expansion (7.11) for p = 2, we have with positive constants ci, c2, c3: Go(w)=,Iw-w,l IG,(w)
IG&,
I
~c,Iw--w*l
w) I dc3
in W;,
(7.14)
in W,,
(7.15)
in W, x T*.
(7.16)
Here, (7.14) is based on (7.10) and (7.13). With this, (7.11) yields in WC X T* the estimate: G(t,
w)a
Iw-w&c,--+%I)+
lt--t,12(~G~(w)--31f--t*l).
(7.17)
By c1 > 0 and G,(w) > 0, this implies G(t, w) > 0 in (W; \ { w,}) X T*. For 7.2(iii): By (7.7), (7.10) and (7.13), G(t,, w) < 0 in Wz \ { w*}. Furthermore, for the two end points t, and t, of T* and all w E Wz, we may achieve: G(t,, w) > 0 and G(t,, w) > 0: For w = w,, this is true by the convexity mentioned before. The convexity then further implies that the function t - G( t, w) has, for any w E Wz \ { we}, exactly one zero in T; \ { t*} and Tz \ { t*}, resp. These two zeros converge, for w + w, (in W,’ \ { we}), to t*, since otherwise there would exist a further zero of t c) G,(t, w*) in T,, besides t,. 0 The index computations to follow are based on Lemmas 7.1 and 7.2. In every case, we then have to do the re-interpretation from the implicit to the parameter representation. For this, we compile here the most important exchange formulas for further use:
ofsurfaces
R. Walter / Visibility
367
Lemma 7.3. If, as in (7.1): F(z) = 0, is an implicit local representation of our surface which has, on the other hand, the parameter representation x, as in (3.141, then the unit normal vector, the second fundamental form and its covariant derivative can be expressed by means of F in the following way: N=s-‘*(grad II(q
j3) =s-‘*((Hess
(vII)(cx,
Here,
(7.18)
F)ox, F)ox)(A,
(7.19)
B),
p, y) =s-‘.((E;,~,~ox)A’~‘C~-_I~(~,
P)((Hess
-II(P,
y)((Hess
F)ox)(N,
A)
-II(y,
a)((Hess
F)ox)(N,
B)).
F)ox)(N,
C)
(7.20)
the abbreviation s:=c](grad
(7.21)
F)oxj
is used, where c is a fixed sign, defined by the equation (7.18). Moreover, a, p, y denote tangent vectors in the domain and A, B, C their images in R3 under dx. Remark. In convention.
differential-geometric
formulas
like
(7.20)
Proof. The formulas (7.18) and (7.19) are to be found obtained by another differentiation of (7.19). q In the following
proofs
of our theorems,
Proof of Theorem 4.1. The Taylor form: G(t,
w) = G,(w)
Here, the zero term is not present G,(w,)
(7.11) takes,
t* + &(t,
use
in [Dombrowski
G always is the function
formula
t + fG2(w)
we sometimes
Einstein’s
sum
‘681. Eq. (7.20)
is
(7.2).
for p = 2 and
t, := 0, w, = w,, the
w) t3.
(7.22)
by F(0, w) = 0. Moreover: (7.23)
= 0;
for we have: G,(w,,)
= G,(O, wO) = ((grad
We make the additional
F)(Z(w,)),
z(wO))
= 5(w0)
= 0.
(7.24)
assumptions:
G&O)
+ 0,
(7.25)
G;(w,)
+ 0,
(7.26)
whose meaning G(t,
will be clarified w) = t.H(t,
with the complementary H(t, Besides the H. Now we function t one, if al/b
in a moment.
By (7.22), we have the factor
w)
decomposition: (7.27a)
factor:
w) := G,(w)
+ $G,(w)
t + $e3(t,
w) t*.
(7.27b)
trivial solution t = 0, which exists for any w, the function G has the same zeros as apply Lemma 7.1 to H. The assumptions made there are fulfilled. Consequently the H(t, w) games a positive zero when w passes through w, if a,/b < 0, and loses > 0. This implies:
D, = - sign( al/b).
(7.28)
368
R. Walter / Visibility of surfaces
It remains to calculate the quantities a, and b, which are associated to H, but may be finally expressed by (7.27), (7.22), (7.2) and Lemma 7.3 by quantities referring to x. In addition, one has to observe that the vectors C(w,), h(c(w,)) of the domain are mapped under dx to the vectors X’(wO), z(w,,) of R3 [cf. (3.10) and (7.24)]: a, =H,(O,
wa) = :G,(w,)
= &I(grad b = &(O,
= :G,,(O,
F)(x(wa))
Y,) = G;(w,)
=&,,(X(wa))
X”(Wo)
w,,) = :&(X(wO))~(wO)
I I~(h(cbd
Ii(,) (7.29)
X(~bd))~
= G,,(O, ~0) z’(wo) +F@(we))
P(w,).
(7.30)
Here, the second term on the RHS vanishes, because it is proportional to (F( wa), z’( w,,)), and because z’ is by (3.6) a multiple of X’. Thus: b = c I(grad F)(x(wa))
I I~(~(wo),
(7.31)
X(C(WO)).
The assumptions (7.25), (7.26) are thus explained with surface data, and the insertion of (7.29), (7.31) into (7.28) yields the assertion (4.1). 0 Proof of Theorem 4.2. Here, the function G has to be studied near the place (to, wa) with t, > 0. The corresponding Taylor formula (7.11) reads for p = 2: G(t,
w) =G,,(w)+G,(w)(t-t,)
+ fG2(~)(t-r,,)2+$G3(?,
w)(~-z~)~. (7.32)
We have: G&a)
= G,(w,)
(7.33)
= 0,
because G(t, wa) has as function of t a double zero at I,, as is seen exactly as for (7.24). We make here the additional assumptions: G2bd
(7.34)
#03
(7.35)
G;( wa) Z 0. This time, Lemma 7.2 can be applied directly to G. This yields D2 = 2, if a,/b D2 = -2, if a,/b > 0, hence: D, = -2
< 0, and (7.36)
sign(a,/b).
Now follows the calculation of a2 und b, analogous to (7.29)-(7.31): u2=G,,(t,,,
~,)=I;;,~(~(w,,)+t~~(w,,))
Here, z(w,) = c3 L(u,, there on the centre (a,,
zj(w,).
(7.37)
ui) with a positive c3, as follows from (3.6) and the assumption made ur are the Gaussian coordinates of pl). Thus
a2 = c3e Ibad b = G,(r,,
L’(w,)
F)(xh,
d)
I I~(%,
wa) = ((grad F)(x(u,,
ud, %,
ul)),
ud),
X’(wa) + raE’(wa)).
(7.38) (7.39)
By (3.6), z’( w,,) = -8 . X’( w,,) with 6 = 0 resp. = 1 in case of parallel resp. central projection. Since in both cases 1 - 6. t, > 0, we get with a positive cq: b = c4c
Ibad
f’)(x(uI,
ul)) I(N(~I,
(7.40)
UI), X’(e)).
Insertion of (7.38), (7.40) into (7.36) gives the assertion (4.2).
0
369
R. Walter / Visibility of surfaces
Proof of Theorem further: G(t,
4.7. Instead
w) = G,(w)
of (7.22),
t + ;G,(w)
the Taylor
t2 + ;G3(w)
expansion
has to be pushed
t3 + &&&,
one term (7.41)
w) t4.
Still we have: G,(w,)
= 0,
(7.42)
G2(w,)
=O,
(7.43)
but also
which springs
from (7.29). We assume
G3(wo)
G;(w,)
+O,
f 0.
(7.44)
Again, we can split th factor t in (7.41) and apply Lemma This gives, with the quantities a,, b, referring to H:
7.2 to the complementary
factor H.
= - sign( a,/b),
0:
(7.45)
however without the factor 2, since the two zeros, which are referred 7.2, have different signs. Calculation of a2: a2=fG3(wJ
= fGt,,(O, wO>
=&I(grad
=f4,j,k(x(wCl))
F)(x(w,))I(vZZ)(A(c(w,)),
L'(wO)
to in part (iii) of Lemma
L'(wCJ)
~(c(wo>)~
G(t,
w) = +G2(w)
expansion
t2 + ;G3(w)
X(c(w)))
= 0,
ZZ(A(c(w,)),
X(c(w,)))
used here is: t3+ &G4(t,
(7.47)
w) t4.
The terms of order zero and one vanish identically due to the assumption this vanishing follows by another derivation for w, using (7.19): ZZ(e(w),
(7.46)
~(cbd)).
At the end, we must apply (7.20), observing that by (7.43) und (7.29): = 0. The calculation of b runs exactly as in (7.29) with the same result. 0 Thus, with (7.45) and (7.46), the assertion (4.12) is proved. Proof of Theorem 4.3. The Taylor
Ek(W13)
(4.4); cf. (7.24). From (7.48)
identically.
If, first, G2(w0) # 0 then the map t - G(t, w) has be (7.47) the exact double zero 0 for any w in a neighborhood of w,. Thus, along X(w), there is not change of the sight number. If however G2(w0) = 0, i.e. by (7.29): ZZ(X(c(w,,)), X(c(w,))) = 0, then we assume G&O)
(7.49)
G;( wa) # 0.
+ 0,
Again, the factor t2 splits in (7.47) and, on the complementary applied. This yields with the quantities a,, b, referring to H:
factor
H, Lemma
7.1 can be
D3 = - sign( q/b)
(7.50)
and: a, =
$3(w)
b = ;G;(w,) = ;~,j,,(x(w,))
=&,(O,
wO)=i4,j,k(x(wlJ))
= :G,,,(O,
z'(wO)
Lj(wO)
zk(WCJ),
(7.51)
w,,)
P(Wa)
E’(w,)
P(W,)
+ 2 q,j(x(wo))
L’(w,)
Z”(,). (7.52)
370
R. Walter / Visibility of surfaces
Since, by (3.6) z’ observing (7.48):
is proportional
b = fQ&(w,))
If the vectors c(w,,) would be:
E(w,),
to X’, the second
Zi’(wO)
zj(,)
term
on the RHS
of (7.52) vanishes,
Zk(wO).
h(c( wO)) were linearly
(7.53)
independent,
then
the Gaussian
curvature
in
(7.54) contradicting
the assumption
(4.6). Thus for a p E Iw :
x’(w()) = p Z(,)
(7.55)
and by (7.50), (7.51), (7.53): Ds=sign(--p)=sign(-(?(w,), as asserted in (4.5). Under the additional because the differential u*
z(,))),
(7.56)
assumptions, the silhouette is a regular curve in a neighborhood of the sight function at w, is by (3.12):
(dc)(c(w,);u)
=Il(X(c(wo),
of w,, (7.57)
u),
and this is not the zero form since A( c( wa)) + 0 and K( c( w,,)) # 0. Hence, such a way that X’(wO) # 0. 0
X can be chosen
in
Proof of Theorem 4.4. By the additional assumptions, the ray through X( wa) cuts the surface at x( q, ut) transversally. Thus, in a neighborhood of this point, there arises a regular intersection curve of the rays through X(w) with the surface, which allows a parameter representation of the form y(w) For its tangent
:=x(w)
+ f(w)
L(w).
(7.58)
vector, we have:
y’ = x’ + t’E 4 tZ’ = (1 - qx
+ t’Z,
where (3.6) has been applied. Since y’(w,,) ( y ‘( wO), N( ut, q)) = 0, which implies:
(7.59)
is tangent
to the surface
at x(u,,
ut), we have
(7.60) The ray through X(w) gains resp. loses at w, an intersection to the sign of (y’(wO), E( ui, u,)): D4 = -sign(y’(we), The conversion
E(u,,
of the scalar product
point
with the surface,
according (7.61)
ui)). in (7.61) using (7.59), (7.60) yields:
(7.62) One has to use that R = N X E along the boundary. (7.62) is positive, we get the assertion from (7.61). Proof of Theorem 4.5. In a neighborhood F(z)
= 0,
B(z)
< 0.
Since the factor in front of the fraction 0
of X(wo), let the bounded
in
surface be given by (7.63)
371
R. Walter/ Visibilityof surfaces Table 1 section with cont. surf
section with bound. surf
sightnumb.
w.r.t.
index w.r.t. bound. surf.
bound. surf. fd ’ 0
c; > 0 c, c 0
fd (0
c; > 0 cg < 0
p: a: p: a:
1 1 1 1
behind before behind before
p: 1 before
remains dropped dropped remains
0 0 0 1
a:
1 behind
remains dropped
1 0
p: a:
1 before 1 behind
dropped remains
0 0
0 1
-1 0
Here, B is a suitable P-function on the neighborhood such that the gradients of F and B are linearly independent there. If we think the surface locally continued over the boundary, then the intersection behaviour with the rays through X(w) near w0 would be as in Theorem 4.1. However, certain of the intersection points calculated in this way will fulfill B > 0 and thus do not belong to the surface with boundary. This changes the counting of the intersection points. The additional intersection point X(w) +f(w)z(w) with the continued surface, calculated by (7.27), lies on the bounded part if C(w) := B(x(w) +f(w)Z(w)) d 0, and otherwise outside. Since the function C has a zero at wO,its sign behaviour is dictated by the derivative: C’(%>
= ((grad B)@(%)),
f’(w&(wO))
= c&%+&
%~J)f’(w,).
(7.64)
Here, c5 is a positive constant, because E und grad B are positively proportional along the boundary. According to the sign of f ‘( wo)=:f, und C’( wO)=: Cd, one has the sectional behaviour given in Table 1. In this table, ‘p’ means previously (w -Z wO), ‘a’ means afterwards (w > wO), and ‘remains’ and ‘dropped’ refer to the passage from the continued surface to the bounded one. (The smallest sight number is normed to 0.) From this table, one reads off the assertion (4.10). n Proof of Theorem 4.6. By the transversality of all curves and surfaces in X(wa), the section problem may here completely be linearized. We only have to determine how the rays along the tangent to X( wO)cut the tangent plane of the ‘other’ surface part. Without restriction, let X( w,,) be the origin. The ray through the point w * X’(wo) of the tangent has the parameter representation t c) z = we X’(wO) + t z(w,). The tangent plane of the other surface part in X( w,,) has the equation (z, N( ui, ui)) = 0. Insertion gives: (7.65) When w runs through w,, the ray through X(w) loses an intersection point before X(w), if t/w < 0, and it gains one, if t/w > 0. This together with (7.65) proves the assertion (4.11). EI
8. First experience and concluding remarks The pictures (Figs. 1-8) were produced by the present methods. Figs. l-3 follow the program structure in the case of a bounded piece of a rotational torus. Fig. 4 shows an example of a general surface curve, a torus knot. The last four pictures are surfaces with selfintersec-
R. Walter / Visibility of surfaces
372
Figs. 5-6. Canal surface with selfintersection.
Figs. 7-8. Constant
mean curvature
torus. Selfintersection
(left), exterior
view (right).
tions, namely a canal surface (Figs. 5, 6) and a torus of constant mean curvature H = 1 (Figs. 7, 8). This H-torus is a real analytic immersion whose Gaussian parameter representation is rather complicated, including several integrals; see [Walter ‘871. Numerically, it must be viewed as almost singular, because its Gaussian curvature varies between (roughly) plus 1 and minus 11 hundred! Nevertheless, its silhouette curves and selfintersection components came out very good. (Some parameter lines showed numerical instabilities which had to be corrected half-automatically. So the complete automatization for such highly complex surfaces should be developed further.) Running times (except for Figs. 7, 8) are about 5 to 15 minutes after surface input and adaptions of some accuracy bounds (on a PC). For simplicity of description, we have made several assumptions which are not really necessary. For example, the same methods can be applied to more than one surface, and for curves outside of them. Also boundaries with comers may be dealt with, etc. We intend to come back to these and other extensions on a subsequent opportunity.
Note added in proof (January 15, 1990)
In a late referee statement, the author was which is different from the present procedure, author), with eventual contact to the material curves may require additional measures in considered here, have something to do with author’s book [Walter ‘891, new edition.
informed on the screen-oriented scan line algorithm of [Lane et al. ‘801, and on a talk of W. Thurston (1985, no title, no place; unknown to the here. The author agrees with this referee that the generation of solution practice. He disagrees with him in that discretizations of curves, as in the discretizations of the surface. - Figs. 6-8 were also published
R. Walter / Visibility of surfaces
313
References Boehm, W., Farin, G. and Kahmann, J. (1984), A survey of curve and surface methods in CAGD, Computer Aided Geometric Design 1, l-60. Dombrowski, P. (1968) Krtimmungsgriissen gleichungsdefinierter Untermanmgfabigkeiten Riemannscher Mannigfaltigkeiten. Math. Nachr. 38, 133-180. Hicks, N.J. (1965) Notes on Differential Geometry. Van Nostrand Reinhold, New York. Lane, J.M., Carpenter, L.C., Whitted, T. and Blinn, J.F. (1980), Scan line methods for displaying parametrically defined surfaces, Comm. ACM 23, 23-34. Newman, W.M. and Sproull, R.F. (1984) Principles of Znteracfioe Computer Graphics. McGraw-Hill, New York, 4 editions 1973 to 1984. Sutherland, I.E., SprouII, R.F. and Schumaker, R.A. (1974) A characterization of ten hidden-surface algorithms, Comput. Surv. 6, l-55. Walter, R. (1968), Uber Grenziibergange in der Differentialgeometrie. Math. 2. 104, 337-348. Walter, R. (1981), Konvexitlt in riemannschen Mamrigfahigkeiten, J. Ber. d. Dt. Math.-Verein. 83, l-31. Walter, R. (1985), Compact hypersurfaces with a constant higher mean curvature function. Math. Ann. 270, 125-145. Walter, R. (1987), Explicit examples to the H-problem of Heinz Hopf, Geom. Ded. 23, 187-213. Walter, R. (1989) Differentialgeometrie, Bibl. Inst. Mannheim/Wien/Ztirich, B.I.-Wissenschaftsverlag, 1. Aufl. 1978, 2. Aufl. 1989. Whitney, H. (1943) Differentiability of the remainder term in Taylor’s formula. Duke Math. J., 10, 153-158.
Computer Aided Geometric Design 7 (1990) 353-373 North-Holland
Visibility of surfaces via differential geometry Rolf WALTER Institut
ftir Mathematik,
Universitiit Dortmund,
D-4600 Dortmund
50, FRG
Presented at Oberwolfach October 1988 and April 1989 Received June 1989; revised November 1989 Abstract. Using the means of differential geometry, we develop a theory for rendering surfaces with complete visibility claryfying. The results are based on the new concept of ‘sight indices’ which control the change of coverings dynamically along arbitrary given curves. The method is entirely curve oriented, so suitable for most graphic outputs, in particular for plotters. Important attributes of surfaces, like silhouette, selfintersection, and boundary curves will be generated and, in fact, be viewed as constitutive for the whole process. Keywords. Surface rendering, visibility, hidden lines, silhouettes, surface intersections, boundary curves
1. Introduction The handling of surfaces including their attributes like boundaries, silhouettes and intersections is a central subject of CAGD. In this paper, we shall attack these questions using the means of differential geometry. For the present, we will concentrate on surface rendering and visibility clarifying with special attention paid to these attributes. Geometric surface design has developed in various branches, especially in the very successful field of approximation theory. We refer to the excellent survey articles of [Sutherland, Sproull & Schumaker ‘741 and [Boehm, Farin & Kahmann ‘841, the book of [Newman & Sproull ‘841, and the vast literature given there. However, the visibility problem is still a challenge. The traditional procedures start, in one way or the other, with discretizations of the surface at an early stage. The present method may be considered as an alternative and supplement to this, in particular when the surfaces come in complicated analytic representations. Its main feature is that it does not use any patches or pixels at all: it is completely curve oriented. The important additives of a surface like boundaries, silhouettes and selfintersections will not only be generated but considered as consstitutive for the whole process. Additional curves, like parameter lines, niveau lines, geodesics, curvature lines, etc. can be added at any time. The visible (or covered) parts of curves will be drawn in one stroke which is especially important for controlling plotters. The growth of data with respect to refinements is roughly linear compared to the quadratic ascent of patches or pixels. No sorting of big fields is needed. On the other hand, due to the complexity of differentialgeometric surfaces, one encounters situations which are not present in finite straight models. So, for higher singularities, further investigations may be necessary and seem possible along similar lines. 2. The main ideas Visibility is a global problem, because decisions have to be made for points far apart on the surface. It is in the spirit of differential geometry to relate global questions to suitable local or pointwise invariants. 0167-8396/90/$03.50
0 1990 - Elsevier Science Publishers B.V. (North-Holland)
R. Walter / Visibility of surfaces
354
Figs. l-3.
Rendering procedure for a bounded piece of a rotational torus. In Fig. 2, the sight indices L$ and D4 refer to the silhouette, while the sight index D5 refers to the boundary as the current curve.
As to the projection scheme, the method here applies to all types as well (parallel or central). In principle, the observer rays have to be cut first with the object in 3-space and then with the image plane. For simplicity of description, we assume that the object is a single surface and consider only curves on this surface. But these restrictions are not really necessary. The observer rays shall be called rays, for short. The first idea is to avoid the calculation of the intersecting points of the rays with the surface but rather to keep a balance of the number of these points along a current curve under investigation. The number of points on a ray, common with the surface and before the current point towards the observer, will be called the sight number. So a sight number of zero means that the current point on the surface is visible. Given a certain projection scheme, the sight number is a globally defined function on the surface. It is constant during transversal intersections and away from boundaries and selfintersections. Changes of the sight number along a current curve happen only occasionally, but these jumps are decisive for our visibility method. We call each change of this type a sight index. The salient point is that these sight indices can be computed from local differential-geometric quantities like the second fundamental form of the surface. The derivation of such index formulas for the most important events is the main theoretical concern here. The actual procedure in the computer is then as follows: (A) Generate the distinguished curves, namely the silhouettes, the selfintersections and eventually the boundaries as point sequences. The mathematical principles for this part are essentially known. We will describe a specific arrangement which has proven very effective. The distinguished curves are generated in full, no sight indices enter this part (see Fig. 1). (B) Clarify the visibility of any surface curve according to the following rules: (B.l) If the sight number is known for one point of a given curve, then it is known for each point on the curve by calculation of the sight indices and the places where these occur. So the curves on the surface do the job of ‘transporting’ the sight number from one surface point to another. (B.2) Starting with one visible point on the surface (or with one point whose sight number is known), the visibility of each surface curve can be clarified individually, together with its plotting. This is in particular possible for each of the distinguished curves like the silhouettes, the selfintersections and the boundaries (see Fig. 2, where also some sight indices are marked, and Fig. 3).
3. Differential-geometric The most general mapping
foundations
surface here will not only be viewed
as a point
set in R3 but rather
x: M+R3, where
M is a two-dimensional,
as a
(3.1) compact,
oriented,
and connected
Hausdorff
manifold
with
R. Walter / Visibility of surfaces
355
boundary aM, and x is immersive, i.e. of rank 2. For simplicity, the differentiability is assumed to be C”, but this can be relaxed considerably, mostly to C* or C3, depending on the problem. For many questions, M may be viewed as a subset of Iw2 and x as a parametrized surface. However, the concept of M as a manifold does not stick to a specific parametrization and also comprehends more complicated topological types. In any case, the set M will be addressed to as the domain of the surface. In the following, square brackets [ ] denote the determinant, angle brackets ( ) denote the scalar product, and the lying cross X denotes the vector product in Iw3. The field of unit normal vectors corresponding to the given orientation is viewed as a mapping: N: M+R3.
(3.2)
The first and second fundamental I := (dx, II:= (dx,
form of the surface
are:
dx),
(3.3)
dN).
(3.4)
The evaluation of such forms, e.g. of ZZ, in p E M for two tangent vectors p, v of M in p will be denoted in full by II( p; p, v) := ((dx),(p), (dN),(v)) or, for short, by JZ(p, v). By the theorema egregium, the expression for the Gaussian curvature is: (3.5) where CL,v are two linearly independent vectors at the domain point p considered. References for further details on differential geometry are e.g. [Hicks ‘651 and [Walter ‘891. Among the properties of the projection scheme, only the directions of the observer rays along the surface play a role in our index formulas. These direction vectors may be described for parallel or central projection in a common way, as follows. In the first case, let a be the constant direction vector of the rays towards the observer, and in the second case, let a be the constant position vector of the center. The rays through the surface points x then have the direction vector L=a-8.x,
(3.6)
where 6 = 0 in case of parallel projection and S = 1 in case of central projection. The observer ray through the surface point x(p) is thus given by the parameter tion t++x(p)+t*L(p),
tER.
representa-
(3.7)
One of its points, say x(p) + t1 . L(p), lies before another, say x(p) + t, . L(p), if t, > t,. As the general symbol for the projection mapping of Iw3 (or lR3\ { a }, in case of central projection) to the image plane B, we use the letter P. The explicit formulas for P are well known but we don’t have to use them in our theory. However, in case of central projection, we assume for simplicity that no ray through a surface point is parallel to the image plane, and that the center never lies on the surface nor between two surface points on the same ray. For any domain point p E M, x(p) is the corresponding object point in Iw3, and P(x( p)) is the corresponding image point in B. (As to the computer handling, there are three corresponding choices for the storage of points which should be used flexibly according to practical requirements.) We will utilize the sight function u defined by the scalar product: a:=(L,
N):
M-l!%.
(3.8)
356
R. Walter / Visibility of surfaces
The zeros of the sight function u yield those points on the surface whose rays are tangential to the surface there. These points are the so called silhouette points. The silhouette C is the following closed subset of M: c:=
{pEM~a(p)=O}.
We will decompose L = (dx)(
the direction
(3.9) vector L of rays in a tangent
and normal
part, according
X) + u. N.
to:
(3.10)
So X is the vector field on M, corresponding to the tangential part of L along the surface in space. (3.10) is evident, since L - u *N is orthogonal to N and thus is of the form (dx)( X). By forming the differential of (3.6) and (3.8) and insertion of (3.10), we obtain do=
(dL,
thus da = ((dx)(X), (de)(p)
N)+ dN)
(L,
dN)
=(
-6.dx,
or, for any p tangent
N) + ((dx)(X)
+u.N,
dN),
(3.11)
to M:
= LL(X, CL),
(3.12)
expressing the differential of the sight function by the second fundamental form. If the mapping x is not injective, then the surface has selfintersections. We define selfintersection (set) S by: S:={p~MIthereexistsq~Mwithpfqandx(p)=x(q)}.
the
(3.13)
It is a consequence of the local injectivity of the map x that S is a closed subset of M. A point of S will be called simple if x(p) has exactly two pre-images under x. Also, the boundary aM is known to form a closed part of M. To distinguish from the three subsets C, S, and EM4 of M are their images in R3 under x, which we shall supply with the additive ‘in the object space’, and the images of those in B under x, which we shall supply with the additive ‘in the image plane ‘. Thus, e.g., C is the silhouette, x(C) is the silhouette in the object space, and P(x(C)) is the silhouette in the image plane. All nine sets C, S, aM, x(C), x(S), x(aM), P(x(C)), P(x(S)), P(x(aM)) are compact. Remark. In classical descriptive geometry, the adjectives instead. So, e.g., x(C) is then called the true silhouette
‘true’, resp. ‘apparent ’ are often used and P(x(C)) the apparent silhouette.
Theorem 3.1. Assume, for a fixed p0 E M, that the ray through x( pO) does not hit the silhouette, the selfintersection, and the boundary in the object space before and including x( p,,); i.e. wheneuer, for some t, 2 0 and somep, E M, we haue x( pO) + t,L( p,,) = x( pl), then p1 6E C U S u 3M. Then the following holds: (i) There exist only finitely many intersection points of this ray with the surface x(M) before and including x( p,,), say n points (ii) There is a neighborhood U of p0 in M such that, for any p E U, the ray through x(p) has exactly n intersection points with the surface x(M) before and including x(p), and again does not hit the silhouette, the selfintersection and the boundary in the object space before and including x(P)* Proof. First, an intersecting point p1 of the given fixed ray with x(M) is isolated on this ray, because otherwise an easy limit argument would show that L( pI) is tangent to the surface at x( pl). Since the intersection of the fixed ray with x(M) is compact, there are only finitely many such pl. Now observe that, generally, two surface points x(q,) and x(q,) are situated on the same ray if and only if q1 and q2 have the same image in the image plane B under the composition P 0 x. And: the points where this map P 0 x: M + B has a vanishing Jacobian are exactly those of the silhouette C.
R. Walter / Visibility
of surfaces
357
This reduces the question, how many points q of M have images x(q) on a ray, to the problem of the cardinality of pre-images of a given fixed point in B under P 0x.There are well known arguments showing that, if the Jacobian does not vanish, this number is constant in a neighborhood of a point where it is finite (cf. [Walter ‘681, e.g.). A little care has to be taken to stay outside C, S, and EM, which is possible by the compactness of these sets. Also, one has to focus on points before x(p), but this does not disturb the local constant argument. So, for p in this neighborhood, any intersection point of the ray through x(p) with x(M) before and including x(p) has only one pre-image in M under x. 0 The general definition of the sight number will be put in the following way: For any p E M, it is the cardinality of the intersection of the halfray t - x(p) + t L(p), t > 0, with the point set x(M) in R3. Remark. Alternatively, one could include the number of pre-images of any such intersection point under x in case that the intersection point is on the selfintersection set x(S) in the object space. But we will not do so. So, the meaning of Theorem 3.1 is that the sight number is finite and constant around any p. of M such that the ray through x( po) does not hit x(C) U x(S) U x( CIM) before and including x( PO). As usual, for local calculations, representation of the form
the mapping x may be replaced by a Gaussian parameter
X1(% u> x:A*R3,
(u, U)_X(U,
u)=
x2@,
u)
(3.14)
,
i X3(% u) 1 where the domain A is now contained in lR*. Any domain point p is then represented by its Gaussian coordinates u, u [within the validity of (3.14)]. We presuppose a little familiarity in passing over from (3.1) to (3.14). [The formal scheme is that the x in (3.14) is in fact the composition of the inverse of a chart of M with the x in (3.1).] The field of unit normal vectors is then correspondingly represented by:
N: A-+R3,
(u, u) *N(u,
u) =
N’(u,
u)
NZ(u,
u)
1 N3(U,
.
(3.15)
u> 1
The orientation of the Gaussian coordinates is always such that [XU, X”,
Nl ’ 0.
(3.16)
Here and in the sequel, indices like U, u denote partial derivatives. In later parts, the visibility shall always be clarified ‘dynamically’ for points running along a given surface curve under consideration, the so called current curue. Any surface curve can locally be described by a parameter representation of the form: w c, x(r+),
U(W))>
(3.17.a)
where the curve parameter w runs through a real interval. The point c(w) := (u(w),
U(W))
(3.17.b)
describes the corresponding curve in the parameter domain A. Its component functions U(W) and u(w) are assumed to be of class C” (with similar possibilities for relaxing as for x). The
358
R. Walter /
restriction of any surface quantity, defined (3.17) will be denoted by a bar. Thus, e.g. x(w)
:=x(c(w)),
Visibility
ofsurfaces
on the domain
L(w):=L(c(w)),
If (3.17) describes the current curve then, sometimes, and the ray through X(w) as the current ray.
A, to the points
of a surface
etc. X(w) is addressed
curve (3.18)
to as the current point,
4. The sight indices The sight number, which is globally defined on M or on x(M), will now be investigated along a given surface curve (3.17) the so called current curve. We call a surface curve distinguished if all its points c(w) lie on one of the sets: silhouette, selfintersection or boundary, and we call it general is c(w) hits those sets only at isolated points of the parameter interval. So, we will have to consider the cases where the current curve is a general one or a distinguished one (‘mixed’ cases are excluded here). We assume the parameter interval of the current curve to be open. Of course, there is no problem at a fixed parameter point w0 of the current curve when the sight number is constant around w,, i.e. in an open interval containing w,. In the contrary case, the sight index at w,, of the current curve exists, by definition, if the sight number on the surface points X(w) has a left and right (proper) limit N,, N_ at w,, and its value is then the difference N+- N_ of these limits. Observe that this implies that the sight number is constant near left and near right of w,,. If the sight number is a piecewise constant function along the current curve, then the sight numbers in the constancy intervals can be calculated from the sight indices at the jump points and a single value of the sight number in one of the constancy intervals. (The sight numbers at the jump points are practically of no interest.) By Theorem 3.1, the sight number changes at most at those points w, of the current curve where the ray through x( wO) hits one of the three sets x(C), x(S), x(aM) before and including x(wa). Any circumstance which causes a change of the sight number along the current curve will be called an euent, for short. For the most important events, we shall calculate the sight index under certain genericity assumptions. General assumption. For simplicity of formulation, we describe each event individually, assuming that no other events take place at the point considered on the ray. In practice, if more than one event is happening at the same parameter value w, but at most one at each point on the corresponding ray and each of them obeys the genericity assumptions individually, then the sight index at w, is gained by addition of the individual ones at w,. There are two different types of events possible because either x( wa) itself may meet the sets x(C), x(S), x( aM> or else a point on the ray through x( w,,) before x( wO) may meet those sets. We say that the event is of point type, resp. of ray type. In the arrangement of results, we follow the degree of complexity of the surface; so we start with the theorems which would suffice in case S=fl, aM=,0. The normal case Theorem 4.1. For w = w,, let the current curve w ++X(w) = x(c( w)) meet the silhouette object space, i.e. a(c(w,,)) = 0. Then the sight index for this event is given by:
~q+o>; 4wo), mwm
D1= -SignIl(c(wo); X(c(w,)),
X(c(wo))) ’
in the
(4.1)
359
R. Walter / Visibility of surfaces
where one has to assume, nonzero.
in addition,
that the numerator
and the denominator
in (4.1)
are
The general assumption implies here that X(wO) is not contained in x( S U aikf ), and similarly in the following theorems. By (3.12), the last assumptions mean intuitively that neither the current curve not the sight vector ‘touch’ the silhouette in the object space at the point considered. Theorem 4.2. For w = w,, let the ray through the point X(wO) of the current curve w c* X(w) = x( c( w)) hit the silhouette in the object space before X( wO), say at x( pl) with p, E C. Then the sight index of this event is given by
(~‘(WIA N D2= -2siWIl(p,; where one has to assume, nonzero.
h(p,), in addition,
PI)> h(p,))’ that the numerator
(4.2) and the denominator
in (4.2)
are
The last assumptions mean intuitively that neither the ray through X(wO) is tangent to the silhouette in the object space at x( p,) nor the ruled surface formed by the rays through the current points X(w) touches the surface x(M) at x( p,). The silhouette
itself
Now, consider the case that our current curve is fully contained in the silhouette. This implies that the numerator in (4.1) is identically 0. Genericly, the denominator in (4.1) will only vanish at isolated parameter values (our assumptions shall guarantee this and, moreover, that the current curve locally exhausts C). If this denominator vanishes at w = w,, then the ray through X( wO) touches the silhouette in the object space at X(wO). This is an interesting phenomenon which is responsible for the seemingly spontaneous stopping of the silhouette in the image plane, known e.g. from the torus (see Fig. 2 at D3 and also the torus picture in Fig. 4).
In the following theorem, there enters the covariant derivative of the second fundamental form which is a Codazzi tensor and therefore will be abbreviated by Cod := 011. For differentialgeometric
(4.3) details, see e.g. [Walter ‘851.
Theorem 4.3. Let the current curve w * c(w) u( c( w)) = 0,
be contained
in the silhouette,
identically.
Fig. 4. Rotational
i.e (4.4)
torus with a torus knot as current
curve.
R. Walter / Visibility of surfaces
360
If, for w = w,, II( c( wO); A( c( wO)), A( c( wO))) # 0, then no point type sight index occurs at w,. If II( c( wO); A( c( w,,)), A( c( wO))) = 0, then this causes a point type sight index of value D3 = -sign(X’(w,,),z(w,)), where, in addition,
(4.5)
it is assumed that
Rk(w0))
(4.6)
f 0,
Cod(c(w,);
h(c(wo)),
+(wo)),
+(w,)))
~0.
(4.7)
Under these assumptions, the current curve is locally open as a map w - c(w) into C, and it can be reparumetrized near w, in such a way that X’(wO) = p - z(w,,) for some p # 0. The boundary case For the discussion of the boundary behaviour, one should note first that the immersion x : M --$ R 3 can be continued locally over any boundary point as an immersion. Let R, resp. E, denote the field of unit tangent, resp. normal vectors of the boundary along the restriction x 1i3M, viewed as mappings of aA4 into R 3. Both fields are tangent to x(M), i.e. (R, N) = (E, N) = 0 on M. The orientation of R is chosen canonically (‘interior to the left’), and that of E according to [E, R, N] > 0. Of course, there is no sight index at w, of point type if c( wO)E aM\ S. Theorem 4.4. For w = w,,, let the ray through the point X( wO) of the current curue w - X(w) = x(c(w)) hit the boundary in the object space before X(wO), say at x( p,) with p1 E 3M. Then the sight index of this event is given by
D = sign[+%)~ w%)~ NPdl 4 m%),
where one has to assume, nonzero.
in addition,
(4.8)
’
NP,))
that the numerator
and the denominator
in (4.8)
are
Now, consider the case that the current curve is a piece of the boundary. If the ray through T(w,,) hits the sets x(C), x(S), and x(&M) before X( wO), then there is no difference to the same question for a general surface curve. Theorem 4.5. Let the current curve w ++c(w) c(w)
E aM,
D5= -sipIl(c(wo); h(c(w,)),
The selfintersection
i.e
identically.
If, for w = w,, this curve meets the silhouette, is given by:
where one has to assume, are nonzero and (E(w,),
describe a piece of the boundary,
(4.9) i.e. a( c( wO)) = 0, then the sight index for this event
1 --skn(E(2w,), z(w,)) X(c(w,)))
in addition, that the numerator z(w,)) z 0.
and the denominator
(4.10)
of the first quotient
case
Here the sight indices will be considered which are connected with a transversal selfintersection whose points are simple, i.e. have exactly two pre-images. So, for x( pl) = x( p2), p1 # p2,
R. Walter / Visibilityof surfaces
361
we presuppose that N( pl), N( p2) are linearly independent. Then the selfintersection consists of curves. If the ray through the point X( wO)of a general current curve meets the selfintersection set in the object space before X(wO), then the corresponding sight index is zero. (Observe that our general assumption excludes the case that this ray point is also on the silhouette or on the boundary in the object space.) Since the cases where the ray meets the silhouette or the boundary in the object space before X(wO) are already covered by the foregoing indices, there remains one further important situation: Theorem 4.6. For w = w,,, let the current curve w ++X(w) = x(c(w)) meet a simple transversal selfintersection in the object space and let p1 be the pre-image of X( wO) under x, different from c(wO). Then the sight index for this event is given by:
06=
. -sign
e+i3)~
N(P& N(p,))
(L(w,),
where one has to assume, nonzero.
in addition,
(4.11)
’ that the numerator
and the denominator
in (4.11)
are
Along a current curve which is a transversal selfintersection, a change of the sight number happens at most when the ray through the current point X( wO) meets the set x(C U S U EM) before X( wO)or at X( w,,). Then the sight index has to be calculated as for a general curve where analogous assumptions have to be made. Exceptional
cases
There are exceptional cases which are not covered by the theorems above. Such exceptional situations may occur because some of the assumptions are violated or by a coincidence of more than one of the ‘pure’ events considered above. Normally, these exceptions are nongeneric, i.e. they vanish if the situation is disturbed slightly. For example, one drastic exception is that the silhouette contains inner points, i.e. whole two-dimensional parts of the surface have tangent sight vector. Nevertheless, this can be removed by a slight change of the projection parameters. In other, less drastic situations, there exists the possibility to derive substitutive index formulas. As a model for such a remedy, we give the following theorem which concerns the case of Theorem 4.1 when the denominator there vanishes: Theorem 4.7. For w = w,, let the current curve w ++X(w) = x( c(w)) meet the silhouette in the object space at a point where the sight vector is tangent to the silhouette, i.e. a(c(w,,)) = 0 and ZZ( c( w,,); h( c( w,,)), A( c( wO)) = 0. Then, instead of (4.1), the sight index for this event is given by:
Z~Wo>; +d~ WWJ)
D’= -SignCod(c(wO); where one has to assume, nonzero.
X(c(w,J,
in addition,
Another exceptional situation curve in R3. For this, an index instead of D,. The procedure for traced out roughly by the proofs complicated then).
A(c(w,,),
that the numerator
X(c(w,)))
’
and the denominator
(4.12) in (4.12)
are
would be, that the boundary itself has selfintersections as a formula could be derived similar to (4.10) where D6 enters establishing index formulas for such exceptional situations are to be given in Section 7 (though the details may become more
362
R. Walter/ Visibilityof surfaces
5. Generation of the distinguished curves The actual computation of the sight indices requires the intersection of the current ray with the silhouette, selfintersection and boundary in the object space. So the data of these three sets must be produced and kept available in the computer. For the present, assume that each of these sets, as viewed in the domain, consists of finitely many curves (components) which are either closed or run from boundary to boundary. Each such component of a distinguished curve can be represented in the computer in form of a finite sequence of points, a polygon. These sequences have to be fine and accurate enough, because they enter subsequent calculations. The boundary components play a somewhat separate role, because normally they are given explicitly by Gaussian parameter representations. Then, of course, a storage as point sequences would be excessive. If, however, the boundary is defined implicitly, e.g. by cutting out ‘windows’ of the surface, then the generation of such boundary components is subject to the same process as for the silhouette and the selfintersection components. Remark. In principle, ‘parametric’ manner manner in form of distinguished curves, as least temporarily,
there are always two possibilities for storing surface points, namely the in form of domain points p (as Gaussian coordinates) and the ‘spatial the points x(p) in object space (as suitable coordinates there). For the points should be stored parametrically, but it may be helpful to use them, in spatial form.
The numerical means for this part are essentially nent is obtained by solving a system of (nonlinear) degree 1, i.e. of the form
known because any distinguished compoequations, which is underdetermined with
fi(UI,..., 24,)=o )...) fn-l(ul )..., UJ = 0.
(5.1)
For the silhouette, we have n = 2, and (u,, u2) = (u, u) are the Gaussian surface coordinates, and fi is the sight function cr. For the intersection of two surfaces, we have n = 4, and (u,, z.+) = (u, u) are the Gaussian coordinates of the first surface, and (u,, uq) = (ii, 6) those of the second, while the f, are the components of the difference of the position vector functions of the two surfaces. As to the numerical solution of the system (5.1) we prefer a generalized Newton’s method which applies as well to all types of distinguished curves and has the advantage that it can also be used for initializing the solutions, i.e. for the global solution part. The local background will only be outlined briefly (without discussion of converges): We start with an arbitrary point a,) and try to find a solution point which is as near to a as possible, say a + h, u=(u,,..., h:=(h,,...,h,). The requirement fi(a+h)=O for iE{l,...,n-1) is, in first order, equivalent to: f,(a)
+ t
f,&)h,
= 07
iE {l,...,
(5 4
n - l},
j=l
where the commas in the index indicate partial derivatives. Among system (5.2) we are searching one for which the squared norm g:=$
i:
(Q2
the solutions
h of the linear
(5.3)
j=l
becomes
minimal,
that
is we like to come near a solution with a small change of a. The of g under the side condition (5.2) leads to the following
Lugrungiun method for the extremum
R Walter / Visibility ofsurfaces
363
equations for the Lagrangian multipliers X,: (54
and for the associated extreme point: n-1 h,
=
c
fk,jWL.
(5.5)
k=l
If the system (5.1) has maximal rank n - 1 at a, then the coefficient matrix of (5.4) is of the Gram type and regular, and the unique inversion of (5.4) leads with (5.5) to an exact solution a# :=a+h
(5.6)
of (5.2) and thus to another approximative solution of (5.1) which usually is better than the original one. The iteration a, a#, a##, . . . yields, under favourable circumstances, an exact solution a * of (5.1) in form of a limit. In practice, it may be useful to replace the partial derivatives occuring in (5.2) by suitable difference quotients. Though this procedure does not exactly generate the solution point nearest to the arbitrary starting point, it is helpful to almost think in this terms. So, the assignment a e a * will be called a quasi-footpoint-projection, or QF-projection, for short. Of course, it may happen that the shows irregularities, nonconvergence or exit of the domain. This is sequence a, a#, a#*,... not as dramatic as it looks at first sight, because the starting point a is arbitrary and may be varied until the convergence becomes satisfactory. Anyway, a first solution point of (5.1) may thus be found. Then one can generate a second one, for example by QF-projection of a suitable neighboring point. Once there are two neighbouring solution points, then the chord of them can be prolonged over the newer point in a small increment, and QF-projection again produces from this an improved third solution. This is of iterative character, always producing from the last two points the next one. Thus, a whole solution component may be generated. This generation process of a solution component must be stopped at the boundary or if the component closes up. For the recognition of this second possibility, it proved useful to introduce a tube around an initial piece of the component and to test the entrance in that tube tangentially, i.e. using two or three successive points instead of one. (Tubes, or parallel sets, are familiar and very useful tools in geometry; see e.g. [Walter ‘811.) In order to generate all solution components, the above process “arbitrary starting point * solution component” should be applied on a sufficiently large set of starting points, which can be obtained by subdividing the domain after some suitable rule (e.g. by halfing). We would like to emphasize that ‘starting point’ here does not mean a first point of solution, but an arbitrary point on which the QF-projection may be applied to produce a first solution. This specific arrangement contributes to the efficiency of this part. After running over to the next starting point, one has to decide whether the QF-projection ‘lands’ on a solution component already generated. This again can be done by a tube construction.
6. Program structure The actual programming is not the purpose of the present paper. So we will indicate only the rough structure in order to show that the ideas, developed here, fit together to a reasonable algorithm. The whole method is composed of two main acts (A) and (B) which were described in Section 2. Act (A) is essentially numerical (see Section 5).
R. Walter / Visibility of surfaces
364
In the first part (B.l) of act (B) one moves along a given current curve in suitably fine steps, intersects the current ray in each step with the distinguished curves and at each event calculates the corresponding sight index, according to the theorems of Section 4. It is important that outside the current point, this intersection has actually to be done very rarely, because a ray whose trace point in the image plane does not meet the distinguished curves there, does not hit them in object space. Now the second part (B.2) of act (B) is performed as follows: With (B.l) at hand, and knowing the sight number of one specific point poo on the surface, the sight number along any current curve can be calculated. In fact, for one specific initial point on the given curve, one can ‘fetch’ the sight number by connecting it through any curve with poo, and then, parallel to the generation of the current curve points, one can update the sight indices by the procedure mentioned before. It is useful to do this for one specific pair of curves, which intersect in poo and to store the visibility data (sight indices and occuring places) for these two ‘base curves’. The initialization of the sight number for any other current curve may then be done by fetching it from the nearest point on the base curves. The role of the current curve can be taken by the distinguished ones or by any general curve, and there is no order, what has to be done first. Only the index formulas to be used depend on the type of the current curve. Since the sight number is a globally defined function on the surface, there are many possibilities to test the sight number of a point calculated by sight indices along curves: Any curve reaching the point considered must deliver the same sight number there! If the sight number for one specific point should not be known a priori, one may find a visible point by a (rough) preliminary search or, more sophisticated, by assigning a fictive sight number to some arbitrary point and renormalizing the lowest sight number found on the silhouette and boundary to zero. Surface
input
The standard way of input of the given surface is by Gaussian parameter representations and their derivatives up to the order 2. (As to the Codazzi tensor vZZ, only sign decisions at the zeros of the quadratic form of ZZ are involved which are numerically possible using ZZ alone.) In principle, any input equivalent to this or capable of producing this, is also possible. A conceivable goal could be to input just the coefficient functions of the derivation equations of Gauss and Weingarten, and to include their integration in some way, because the necessary derivatives then come out as by-products. Runtime acceleration
We mention some means to reduce the runtime: The intersection of current rays with the distinguished components occurs very rarely, as we saw above. A further reduction is possible by first doing the intersection for the projections in the image plane and by a sort of metric preliminary test: if a point has larger distance from the vertices of a polygon as is given by the maximal distance of two consecutive vertices, then the point is not on the polygon (and similar for the intersection of a straight interval with a polygon). Furthermore, the image plane may be divided in a small number of rectangles. This induces a corresponding division of the polygons which represent the projections of the distinguished components. Search through the whole polygon may then be replaced by search through some few of these subpolygons, because the others are too far away. As far as operations in the image plane are necessary, it is useful to replace (at least temporarily) Cartesian space coordinates by triplets whose first pair are coordinates in the image plane and whose third component is the distance from the projection
365
R. Walter / Visibility of surfaces
plane measured on the ray, the height. Of course, there are easy and quick transformations between these two types of space coordinates.
7. Proof of the index formulas In the whole proof we will use, in the neighborhood surface representation of the form F(z)
of a point zi of x(M),
= 0.
an implicit (7.1)
By this, the intersection problem of the rays with the surface may be viewed under the bifurcation aspect of a single scalar equation, depending on one additional parameter. Finally, the result can be re-interpreted for the Gaussian parameter representation of the surface. In (7.1), F is a real P-function of rank 1, defined on a spatial neighborhood U of zi such that x(V) n U is the zero set of F (for some open parameter domain I’). The insertion of the points of the current rays along a curve X(w) = x(c(w)) under consideration yields a real function of two variables: G(t,
w) := F@(w)
(7.2)
+ t.z(w)).
If the ray point X(w) + t. z(w) is contained in U, then the condition G(t, w) = 0 expresses that this ray point is on the surface part x(V). So, we shall have to study the zero behaviour of the real function t - G( t, w) in dependence on the parameter w. For this, the following two lemmas are basic. These lemmas refer, quite generally, to a real P-function (t, w) - G(t, w), defined in a neighborhood Q of a fixed point (t*, w*) of R*. Let the corresponding zero set be denoted by: w)~QandG(t,
Z(w):={t~Rl(t,
w)=O}.
(7.3)
Of course, the differentiability can again be relaxed considerably, as the proofs will show. Without loss of generality, neighborhoods in R shall be intervals, and for any neighborhood W of w,, let W+, resp. W- denote the left, resp. right part of W (including w, itself). Lemma 7.1. Assume:
G(t,,
w*> =0,
al := G,(t,,
we) # 0,
b:= G,(t,,
w*) # 0.
Then there exist neighborhoods T* oft, and W, of w, such that: (i) For any w E W,, Z(w) n T* consists of exactly one number t = f(w). (ii) f(w*) = t,. (iii) The function f is in W, of class C”, and strictly increasing if a,/b decreasing if a/b > 0.
(7.4) (7.5) (7.6)
< 0, and strictly
Lemma 7.2. Assume:
G(t,,
w*> =0,
(7.7)
G,(t,,
w*) = 0,
(7.8)
a2 := G,,(t*, b := G,,,(t,,
we) + 0, w*) # 0.
(7.9) (7.10)
366
R. Walter / Visibility of surfaces
Then there exist neighborhoodrr T* of t, and W, of w, such that, in case a/b < 0, we have: (i) For any w E Wi \ { we}, the set Z(w) n T* is empty. (ii) Z(w*) n T* = {t*}. (iii) For any w E Wz \ { w*}, the set Z(w) n T* consists of exactly two points on different sides oft*, and these points converge, for w + w,, (in W,‘), to t,. In case a,/b > 0 the same is true with interchanged roles of Wi and Wz. For the proof, it is very useful, to utilize generalized Taylor formulas, which expand along increasing powers of t - t,, but have main coefficients depending on the parameter w; cf. [Whitney ‘431. By this, we have for fixed p > 0:
5 $Gq(w)(t-t.)‘+ (p:l) ~p+l(t>w)(t - t*)p+l,
G(t, w)=
!
(7.11)
q=o
where VG
Gq(W)=F(‘*,
w)
(7.12)
forO
but where also the GP+i are of class C”. Proofs. Lemma 7.1 is a direct consequence of the standard theorem on implicit functions. We displayed it here because, in connection with certain factor decompositions, also higher cases can be dealt with it. Lemma 7.2 will be treated, assuming a2 > 0,
(7.13)
b
the other cases run analogously. First of all, the function t - G(t, w) is near (t*, w*) strictly convex: By a2 > 0, there are neighborhoods T* of t* and W, of w, such that G,,(t, w) > 0 in T* X W,. By (7.7) and (7.8), this implies the convexity. In the sequel, we shall have to reduce the neighborhoods T* and W, again. For 7.2 (i), (ii): In the expansion (7.11) for p = 2, we have with positive constants ci, c2, c3: Go(w)=,Iw-w,l IG,(w)
IG&,
I
~c,Iw--w*l
w) I dc3
in W;,
(7.14)
in W,,
(7.15)
in W, x T*.
(7.16)
Here, (7.14) is based on (7.10) and (7.13). With this, (7.11) yields in WC X T* the estimate: G(t,
w)a
Iw-w&c,--+%I)+
lt--t,12(~G~(w)--31f--t*l).
(7.17)
By c1 > 0 and G,(w) > 0, this implies G(t, w) > 0 in (W; \ { w,}) X T*. For 7.2(iii): By (7.7), (7.10) and (7.13), G(t,, w) < 0 in Wz \ { w*}. Furthermore, for the two end points t, and t, of T* and all w E Wz, we may achieve: G(t,, w) > 0 and G(t,, w) > 0: For w = w,, this is true by the convexity mentioned before. The convexity then further implies that the function t - G( t, w) has, for any w E Wz \ { we}, exactly one zero in T; \ { t*} and Tz \ { t*}, resp. These two zeros converge, for w + w, (in W,’ \ { we}), to t*, since otherwise there would exist a further zero of t c) G,(t, w*) in T,, besides t,. 0 The index computations to follow are based on Lemmas 7.1 and 7.2. In every case, we then have to do the re-interpretation from the implicit to the parameter representation. For this, we compile here the most important exchange formulas for further use:
ofsurfaces
R. Walter / Visibility
367
Lemma 7.3. If, as in (7.1): F(z) = 0, is an implicit local representation of our surface which has, on the other hand, the parameter representation x, as in (3.141, then the unit normal vector, the second fundamental form and its covariant derivative can be expressed by means of F in the following way: N=s-‘*(grad II(q
j3) =s-‘*((Hess
(vII)(cx,
Here,
(7.18)
F)ox, F)ox)(A,
(7.19)
B),
p, y) =s-‘.((E;,~,~ox)A’~‘C~-_I~(~,
P)((Hess
-II(P,
y)((Hess
F)ox)(N,
A)
-II(y,
a)((Hess
F)ox)(N,
B)).
F)ox)(N,
C)
(7.20)
the abbreviation s:=c](grad
(7.21)
F)oxj
is used, where c is a fixed sign, defined by the equation (7.18). Moreover, a, p, y denote tangent vectors in the domain and A, B, C their images in R3 under dx. Remark. In convention.
differential-geometric
formulas
like
(7.20)
Proof. The formulas (7.18) and (7.19) are to be found obtained by another differentiation of (7.19). q In the following
proofs
of our theorems,
Proof of Theorem 4.1. The Taylor form: G(t,
w) = G,(w)
Here, the zero term is not present G,(w,)
(7.11) takes,
t* + &(t,
use
in [Dombrowski
G always is the function
formula
t + fG2(w)
we sometimes
Einstein’s
sum
‘681. Eq. (7.20)
is
(7.2).
for p = 2 and
t, := 0, w, = w,, the
w) t3.
(7.22)
by F(0, w) = 0. Moreover: (7.23)
= 0;
for we have: G,(w,,)
= G,(O, wO) = ((grad
We make the additional
F)(Z(w,)),
z(wO))
= 5(w0)
= 0.
(7.24)
assumptions:
G&O)
+ 0,
(7.25)
G;(w,)
+ 0,
(7.26)
whose meaning G(t,
will be clarified w) = t.H(t,
with the complementary H(t, Besides the H. Now we function t one, if al/b
in a moment.
By (7.22), we have the factor
w)
decomposition: (7.27a)
factor:
w) := G,(w)
+ $G,(w)
t + $e3(t,
w) t*.
(7.27b)
trivial solution t = 0, which exists for any w, the function G has the same zeros as apply Lemma 7.1 to H. The assumptions made there are fulfilled. Consequently the H(t, w) games a positive zero when w passes through w, if a,/b < 0, and loses > 0. This implies:
D, = - sign( al/b).
(7.28)
368
R. Walter / Visibility of surfaces
It remains to calculate the quantities a, and b, which are associated to H, but may be finally expressed by (7.27), (7.22), (7.2) and Lemma 7.3 by quantities referring to x. In addition, one has to observe that the vectors C(w,), h(c(w,)) of the domain are mapped under dx to the vectors X’(wO), z(w,,) of R3 [cf. (3.10) and (7.24)]: a, =H,(O,
wa) = :G,(w,)
= &I(grad b = &(O,
= :G,,(O,
F)(x(wa))
Y,) = G;(w,)
=&,,(X(wa))
X”(Wo)
w,,) = :&(X(wO))~(wO)
I I~(h(cbd
Ii(,) (7.29)
X(~bd))~
= G,,(O, ~0) z’(wo) +F@(we))
P(w,).
(7.30)
Here, the second term on the RHS vanishes, because it is proportional to (F( wa), z’( w,,)), and because z’ is by (3.6) a multiple of X’. Thus: b = c I(grad F)(x(wa))
I I~(~(wo),
(7.31)
X(C(WO)).
The assumptions (7.25), (7.26) are thus explained with surface data, and the insertion of (7.29), (7.31) into (7.28) yields the assertion (4.1). 0 Proof of Theorem 4.2. Here, the function G has to be studied near the place (to, wa) with t, > 0. The corresponding Taylor formula (7.11) reads for p = 2: G(t,
w) =G,,(w)+G,(w)(t-t,)
+ fG2(~)(t-r,,)2+$G3(?,
w)(~-z~)~. (7.32)
We have: G&a)
= G,(w,)
(7.33)
= 0,
because G(t, wa) has as function of t a double zero at I,, as is seen exactly as for (7.24). We make here the additional assumptions: G2bd
(7.34)
#03
(7.35)
G;( wa) Z 0. This time, Lemma 7.2 can be applied directly to G. This yields D2 = 2, if a,/b D2 = -2, if a,/b > 0, hence: D, = -2
< 0, and (7.36)
sign(a,/b).
Now follows the calculation of a2 und b, analogous to (7.29)-(7.31): u2=G,,(t,,,
~,)=I;;,~(~(w,,)+t~~(w,,))
Here, z(w,) = c3 L(u,, there on the centre (a,,
zj(w,).
(7.37)
ui) with a positive c3, as follows from (3.6) and the assumption made ur are the Gaussian coordinates of pl). Thus
a2 = c3e Ibad b = G,(r,,
L’(w,)
F)(xh,
d)
I I~(%,
wa) = ((grad F)(x(u,,
ud, %,
ul)),
ud),
X’(wa) + raE’(wa)).
(7.38) (7.39)
By (3.6), z’( w,,) = -8 . X’( w,,) with 6 = 0 resp. = 1 in case of parallel resp. central projection. Since in both cases 1 - 6. t, > 0, we get with a positive cq: b = c4c
Ibad
f’)(x(uI,
ul)) I(N(~I,
(7.40)
UI), X’(e)).
Insertion of (7.38), (7.40) into (7.36) gives the assertion (4.2).
0
369
R. Walter / Visibility of surfaces
Proof of Theorem further: G(t,
4.7. Instead
w) = G,(w)
of (7.22),
t + ;G,(w)
the Taylor
t2 + ;G3(w)
expansion
has to be pushed
t3 + &&&,
one term (7.41)
w) t4.
Still we have: G,(w,)
= 0,
(7.42)
G2(w,)
=O,
(7.43)
but also
which springs
from (7.29). We assume
G3(wo)
G;(w,)
+O,
f 0.
(7.44)
Again, we can split th factor t in (7.41) and apply Lemma This gives, with the quantities a,, b, referring to H:
7.2 to the complementary
factor H.
= - sign( a,/b),
0:
(7.45)
however without the factor 2, since the two zeros, which are referred 7.2, have different signs. Calculation of a2: a2=fG3(wJ
= fGt,,(O, wO>
=&I(grad
=f4,j,k(x(wCl))
F)(x(w,))I(vZZ)(A(c(w,)),
L'(wO)
to in part (iii) of Lemma
L'(wCJ)
~(c(wo>)~
G(t,
w) = +G2(w)
expansion
t2 + ;G3(w)
X(c(w)))
= 0,
ZZ(A(c(w,)),
X(c(w,)))
used here is: t3+ &G4(t,
(7.47)
w) t4.
The terms of order zero and one vanish identically due to the assumption this vanishing follows by another derivation for w, using (7.19): ZZ(e(w),
(7.46)
~(cbd)).
At the end, we must apply (7.20), observing that by (7.43) und (7.29): = 0. The calculation of b runs exactly as in (7.29) with the same result. 0 Thus, with (7.45) and (7.46), the assertion (4.12) is proved. Proof of Theorem 4.3. The Taylor
Ek(W13)
(4.4); cf. (7.24). From (7.48)
identically.
If, first, G2(w0) # 0 then the map t - G(t, w) has be (7.47) the exact double zero 0 for any w in a neighborhood of w,. Thus, along X(w), there is not change of the sight number. If however G2(w0) = 0, i.e. by (7.29): ZZ(X(c(w,,)), X(c(w,))) = 0, then we assume G&O)
(7.49)
G;( wa) # 0.
+ 0,
Again, the factor t2 splits in (7.47) and, on the complementary applied. This yields with the quantities a,, b, referring to H:
factor
H, Lemma
7.1 can be
D3 = - sign( q/b)
(7.50)
and: a, =
$3(w)
b = ;G;(w,) = ;~,j,,(x(w,))
=&,(O,
wO)=i4,j,k(x(wlJ))
= :G,,,(O,
z'(wO)
Lj(wO)
zk(WCJ),
(7.51)
w,,)
P(Wa)
E’(w,)
P(W,)
+ 2 q,j(x(wo))
L’(w,)
Z”(,). (7.52)
370
R. Walter / Visibility of surfaces
Since, by (3.6) z’ observing (7.48):
is proportional
b = fQ&(w,))
If the vectors c(w,,) would be:
E(w,),
to X’, the second
Zi’(wO)
zj(,)
term
on the RHS
of (7.52) vanishes,
Zk(wO).
h(c( wO)) were linearly
(7.53)
independent,
then
the Gaussian
curvature
in
(7.54) contradicting
the assumption
(4.6). Thus for a p E Iw :
x’(w()) = p Z(,)
(7.55)
and by (7.50), (7.51), (7.53): Ds=sign(--p)=sign(-(?(w,), as asserted in (4.5). Under the additional because the differential u*
z(,))),
(7.56)
assumptions, the silhouette is a regular curve in a neighborhood of the sight function at w, is by (3.12):
(dc)(c(w,);u)
=Il(X(c(wo),
of w,, (7.57)
u),
and this is not the zero form since A( c( wa)) + 0 and K( c( w,,)) # 0. Hence, such a way that X’(wO) # 0. 0
X can be chosen
in
Proof of Theorem 4.4. By the additional assumptions, the ray through X( wa) cuts the surface at x( q, ut) transversally. Thus, in a neighborhood of this point, there arises a regular intersection curve of the rays through X(w) with the surface, which allows a parameter representation of the form y(w) For its tangent
:=x(w)
+ f(w)
L(w).
(7.58)
vector, we have:
y’ = x’ + t’E 4 tZ’ = (1 - qx
+ t’Z,
where (3.6) has been applied. Since y’(w,,) ( y ‘( wO), N( ut, q)) = 0, which implies:
(7.59)
is tangent
to the surface
at x(u,,
ut), we have
(7.60) The ray through X(w) gains resp. loses at w, an intersection to the sign of (y’(wO), E( ui, u,)): D4 = -sign(y’(we), The conversion
E(u,,
of the scalar product
point
with the surface,
according (7.61)
ui)). in (7.61) using (7.59), (7.60) yields:
(7.62) One has to use that R = N X E along the boundary. (7.62) is positive, we get the assertion from (7.61). Proof of Theorem 4.5. In a neighborhood F(z)
= 0,
B(z)
< 0.
Since the factor in front of the fraction 0
of X(wo), let the bounded
in
surface be given by (7.63)
371
R. Walter/ Visibilityof surfaces Table 1 section with cont. surf
section with bound. surf
sightnumb.
w.r.t.
index w.r.t. bound. surf.
bound. surf. fd ’ 0
c; > 0 c, c 0
fd (0
c; > 0 cg < 0
p: a: p: a:
1 1 1 1
behind before behind before
p: 1 before
remains dropped dropped remains
0 0 0 1
a:
1 behind
remains dropped
1 0
p: a:
1 before 1 behind
dropped remains
0 0
0 1
-1 0
Here, B is a suitable P-function on the neighborhood such that the gradients of F and B are linearly independent there. If we think the surface locally continued over the boundary, then the intersection behaviour with the rays through X(w) near w0 would be as in Theorem 4.1. However, certain of the intersection points calculated in this way will fulfill B > 0 and thus do not belong to the surface with boundary. This changes the counting of the intersection points. The additional intersection point X(w) +f(w)z(w) with the continued surface, calculated by (7.27), lies on the bounded part if C(w) := B(x(w) +f(w)Z(w)) d 0, and otherwise outside. Since the function C has a zero at wO,its sign behaviour is dictated by the derivative: C’(%>
= ((grad B)@(%)),
f’(w&(wO))
= c&%+&
%~J)f’(w,).
(7.64)
Here, c5 is a positive constant, because E und grad B are positively proportional along the boundary. According to the sign of f ‘( wo)=:f, und C’( wO)=: Cd, one has the sectional behaviour given in Table 1. In this table, ‘p’ means previously (w -Z wO), ‘a’ means afterwards (w > wO), and ‘remains’ and ‘dropped’ refer to the passage from the continued surface to the bounded one. (The smallest sight number is normed to 0.) From this table, one reads off the assertion (4.10). n Proof of Theorem 4.6. By the transversality of all curves and surfaces in X(wa), the section problem may here completely be linearized. We only have to determine how the rays along the tangent to X( wO)cut the tangent plane of the ‘other’ surface part. Without restriction, let X( w,,) be the origin. The ray through the point w * X’(wo) of the tangent has the parameter representation t c) z = we X’(wO) + t z(w,). The tangent plane of the other surface part in X( w,,) has the equation (z, N( ui, ui)) = 0. Insertion gives: (7.65) When w runs through w,, the ray through X(w) loses an intersection point before X(w), if t/w < 0, and it gains one, if t/w > 0. This together with (7.65) proves the assertion (4.11). EI
8. First experience and concluding remarks The pictures (Figs. 1-8) were produced by the present methods. Figs. l-3 follow the program structure in the case of a bounded piece of a rotational torus. Fig. 4 shows an example of a general surface curve, a torus knot. The last four pictures are surfaces with selfintersec-
R. Walter / Visibility of surfaces
372
Figs. 5-6. Canal surface with selfintersection.
Figs. 7-8. Constant
mean curvature
torus. Selfintersection
(left), exterior
view (right).
tions, namely a canal surface (Figs. 5, 6) and a torus of constant mean curvature H = 1 (Figs. 7, 8). This H-torus is a real analytic immersion whose Gaussian parameter representation is rather complicated, including several integrals; see [Walter ‘871. Numerically, it must be viewed as almost singular, because its Gaussian curvature varies between (roughly) plus 1 and minus 11 hundred! Nevertheless, its silhouette curves and selfintersection components came out very good. (Some parameter lines showed numerical instabilities which had to be corrected half-automatically. So the complete automatization for such highly complex surfaces should be developed further.) Running times (except for Figs. 7, 8) are about 5 to 15 minutes after surface input and adaptions of some accuracy bounds (on a PC). For simplicity of description, we have made several assumptions which are not really necessary. For example, the same methods can be applied to more than one surface, and for curves outside of them. Also boundaries with comers may be dealt with, etc. We intend to come back to these and other extensions on a subsequent opportunity.
Note added in proof (January 15, 1990)
In a late referee statement, the author was which is different from the present procedure, author), with eventual contact to the material curves may require additional measures in considered here, have something to do with author’s book [Walter ‘891, new edition.
informed on the screen-oriented scan line algorithm of [Lane et al. ‘801, and on a talk of W. Thurston (1985, no title, no place; unknown to the here. The author agrees with this referee that the generation of solution practice. He disagrees with him in that discretizations of curves, as in the discretizations of the surface. - Figs. 6-8 were also published
R. Walter / Visibility of surfaces
313
References Boehm, W., Farin, G. and Kahmann, J. (1984), A survey of curve and surface methods in CAGD, Computer Aided Geometric Design 1, l-60. Dombrowski, P. (1968) Krtimmungsgriissen gleichungsdefinierter Untermanmgfabigkeiten Riemannscher Mannigfaltigkeiten. Math. Nachr. 38, 133-180. Hicks, N.J. (1965) Notes on Differential Geometry. Van Nostrand Reinhold, New York. Lane, J.M., Carpenter, L.C., Whitted, T. and Blinn, J.F. (1980), Scan line methods for displaying parametrically defined surfaces, Comm. ACM 23, 23-34. Newman, W.M. and Sproull, R.F. (1984) Principles of Znteracfioe Computer Graphics. McGraw-Hill, New York, 4 editions 1973 to 1984. Sutherland, I.E., SprouII, R.F. and Schumaker, R.A. (1974) A characterization of ten hidden-surface algorithms, Comput. Surv. 6, l-55. Walter, R. (1968), Uber Grenziibergange in der Differentialgeometrie. Math. 2. 104, 337-348. Walter, R. (1981), Konvexitlt in riemannschen Mamrigfahigkeiten, J. Ber. d. Dt. Math.-Verein. 83, l-31. Walter, R. (1985), Compact hypersurfaces with a constant higher mean curvature function. Math. Ann. 270, 125-145. Walter, R. (1987), Explicit examples to the H-problem of Heinz Hopf, Geom. Ded. 23, 187-213. Walter, R. (1989) Differentialgeometrie, Bibl. Inst. Mannheim/Wien/Ztirich, B.I.-Wissenschaftsverlag, 1. Aufl. 1978, 2. Aufl. 1989. Whitney, H. (1943) Differentiability of the remainder term in Taylor’s formula. Duke Math. J., 10, 153-158.