An integer algorithm for rendering curved surfaces

An integer algorithm for rendering curved surfaces

Citation for published version (APA):

Overveld, van, C. W. A. M. (1987). An integer algorithm for rendering curved surfaces. (Computing science notes; Vol. 8718). Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/1987

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

An

integer algorithm for rendering

curved surfaces

by

C.W.A.M. van Overveld

87/18

11 -is is a series of notes of the Computing Science Section of the Department

of Mathematics and Computing Science of Eindhoven University of

Technol-ogy.

Since many of these notes are preliminary versions

or

may

be

published

else-where, they have a limited distribution only and

are

not for review.

Copies of

these

notes

are

available from the author or the editor.

Eindhoven University of Technology

Department of Mathematics and Computing Science

P.O. Box 513

5600 MB Eindhoven

The Netberlands

All rights reserved

·

.

An integer algorithm for rendering curved surfaces

C.W.A.M. van Overveld

Eindhoven University of Technology p.o. box 513

5600 MB Eindhoven The Nelherlands

ABSTRACT

In this I. ,:er we address an algorilhm for rendering free form surfaces, staning from boundary represenations. The melhod being used is a combination of bilinear blending and bicubic interpolation, combined in one recursion scheme in order to allow for on-the-fly texture mapping of colour texture, normal vector modulation and also some forms of geomeuic transformations. A correction against undersampling Ihe texture maps (anti-aliasing) takes place as a pre-process 10 avoid slowing down Ihe rendering process. Since Ihe bilinear part of Ihe algorilhm is based on the genemtion of piecewise cubic boundary curves, we derive, staning from some general syturnetry considemtions, a class of boundary curve algorithms which includes Ihe De Casteljau algorithm for four-point Bezier curves.

o.

Introduction

Wilhin Ihe context of computer aided geomeuical design of 3-D objects, surface modeling plays an important role. Notwilhstanding Ihe growing variety of basic 3-D elements !hat is becoming available in the realm of solid modeling [WU861, Ihe relatively simple data structures and rendering algorilhms, togelher with Ihe diversity of the achievable surfaces makes surface modeling prefemble in those cases where 'solid information' is not explicitly needed.

An often used approach is to describe a surface by first specifying one or more rectangular mauices of control points, which are located in 3-space, and next applying either an itemtive algorilhm (such as in [NEW 791, p. 328) or, preferably, a recursive algorithm (such as in [L1E87]) in order 10 render all points of Ihe associated Bezier- or B-spline surface.

However frequently used, Ihere seem 10 be some discrepancies between the representation of surfaces in terms of a matrix of control points and Ihe associated Bezier- or B-spline surface and Ihe representation of Ihe surface in the designer's mind. Indeed, when asking a designer for a hand dmft of the surface being designed, one most likely gets a dmwing of its boundaries instead of a set of contrOl points as an answer, sometimes wilh one or more additional curves 10 indicate the internal form of the surface. But in Ihe latter case Ihese additional curves usually are attached 10 the boundary, as a guide 10 the eye, so that in this case a boundary with one additional curve can be seen as a composite surface build from two surfaces specified by Iheir boundaries, provided a suitable continuity condition:

..

One can object that this is due to the limitations of a 2-dimensional projection, but there are indications that it is a more or less intrinsic property of human designers:

-firstly, the human perceptive system is organized such that boundaries of objects rather than their solid interior draw the first attention. It is probably for

this

reason that in the development of artificial vision systems much effort has been spent in boundary detection techniques;

-secondly, if a designer builds a physical model of a complex shaped surface he most likely also proceeds by shaping the (3-D) boundary curves first, say, of metal wire before filling with e.g. paper mache.

This relative importance of the boundary of a surface over the interior is by no means reflected in the matrix of control points: a nXln matrix contains (n-2)(m-2) internal points and only 2n+2m-4 boun-dary points, indicating that for large values of m and n much more effort is needed to specify the inte-rior region than the boundary.

Apart from the disproportion between boundary and interior in surface design using conventional Bezier· or B-spline techniques, there is what we will call the problem of interpolating vs.

approxi-mating the given control points. This problem rises from the fact that a substantial part of so-called free form designing is not as free as it seems. A surface being designed generally is part of a larger structure, and in order to have things fit nicely together it is necessary to define accurate interfaces between the several parts of this structure, these interfaces being again (among other things) the boun· daries of the designed surfaces. When using the control points technique for Bezier- or B-spline sur-faces, the consistency of this interface definition obliges the designer to use the same control points technique for all neighbour surfaces, and next for all their neighbour surfaces, and so on. In other words: using the technique of control points .and approximate surface definitions it is not well possible to combine different designing strategies for several components of one and the same object If interpo-lation would be used instead of approximation, or maybe even an analytic description of the boundary curves, one could much easier meet engineering requirements in the sense that certain parts of boun-daries should have exactly predefined forms such as straight line pieces and circular arcs, sine functions or otherwise specified forms.

The concept of linearly blended surfaces, as introduced by Coons [C0067] provides for an alternative methodology for defining surfaces. In his description,

\) only the boundary of the surface is being specified; 2) the surface passes exactly through this boundary;

Unfortunately, two adjacent Coons patches will be CO , even if their boundaries form a C1 network. In this paper we present a family of surfaces which can be seen as C1 generalisations of the bilinear Coons patch. We discuss a simple algorithm for rendering this family of surfaces which uses integer arithmetic and bitwise operations only and is therefore well-suited for machine code or hardware imple-mentation. It is based on a two-step approach:

3

-step I) Starting from rhe description of rhe boundaries of rhe surface, a set of 16-tuples of con-lrolpoinlS is generated using a recursive version of bilinear blending, each 16-tuple describing one bicubic Bezier patch.

step 2) This bicubic patch is rendered subsequently, using a recursive subdivision technique. The recursion schemes of rhe two steps are compatible in rhe sense rilat rbings can be arranged such rhat rhe algorithm of step 2) is been executed after rhe stopping criterion of rhe algorithm of step I) has been met. In orher words: step I) takes care of rhe topmost recursion levels and step 2) of rhe lower-most levels, rhus avoiding intermediate storage of all 16-tuples.

The subdivion scheme rhat is been used for step I) consists of repeated subdivision of rhe surface into four subsurfaces according to the following scheme:

For a reason to be explained in a minute, in step 2) we render a bicubic patch by. first replacing it by 4 equivalent biquadratic patches. These patches are equivalent in the sense

that

the same boundary condi-tions are being met. Next we render every biquadratic patch, but, for reasons of efficiency, we don't use subdivision into four subpatches, corresponding to simultanuous subdivision in !he directions of borh parameters, but into 2 subpatches, corresponding to subdivision in rhe direction of only one parameter. We continue rhis subdivision until the widLb of one subpatch becomes small enough to allow for a qua-dratic Bezier

curve

algoriLbm to do Lbe subdivisions in Lbe direction of Lbe remaining parameter.

The recursive nature of the total algorithm will have the following implications:

I) all boundary descriptions which are recursively computable could in principle be used; 2) in this way we avoid the integrity problem as addressed in [T AM85];

3) computing rhe surface normal vector in every point of the surface, necessary for implementing a Phong-like shading model is accomplished 'on rhe fly' without the loss of accuracy rhat one would expect on behalf of rhe diminishing size of the subsurfaces;

4) texture mapping, boLb color texture and normal vector modulation, and some additional geometric transformations can be supported.

In section 1 we address some mathematical properties of bilinear blended surfaces. The algorithm for their recursive generation (the algorithm for step 1) ) is derived in section 2, making use of the tech-nique of invariants [OVE86]. The algorithm for step 2) is dealt with in section 3. In section 4 we dis-cuss the computation of normal vectors and texture mapping. Some

aspects

of the computation of the boundary curves

are

dealt with in section 5, and we conclude with a discussion, some examples and a summary of our conclusions in section 6.

1. Recursively defined surfaces

We will describe the surface using a function

1

1:(u,v)IO;;u';;I" O;;v';;I)-+{(x,y,z)lx,y,zeR) subject to the boundary conditions:

1 (O,v);go(v) where 1 (I,v);g,(v) 1 (u ,O);ho(u) l(u,l);h,(u) It 0(0); g 0(0) ho(l);g ,(0) h,(O);go(I) h ,(I);g ,(I) (0) (Ia) (Ib) (Ic) (Id) (2a) (2b) (2c) (2d) In order to relate the values of

1

in the interior of a region to the boundary of

this

region, we demand that for any u, uo, u,: OSu~u';;u,SI and v, vo, v,: OSv~v';;vpl:

vo+v, , 1 (u

'-2->=}

(f (u ,vol+/(u ,v,) (3a) (u,-u) vo+v, + ( ) (2/ (uo.-2-)--1 (u 0 •• 0)--/ (uo.v,» uI-U O (u-uo) vo+v, I ( ) (21 (u "-2-}-1 (u "vo}-I (u"v,))) UI-UO UO+Ul 1 1 (-2-'V>=} (f (uo.v)+/(u"v) (3b) (v ,-v) uo+u,

+ (

)

(21 (-2-,vo}-1 (uo.vo)--I (u"vo»

VI-Va

(v-vo) uo+u,

+ (

)

(21 (-2-,v,}-1 (uo,v,}-I (u"v

,»)

VI-VO

5

-uo+u, vo+v,

These definitions have been chosen such that 1(-2-'-2-) is unique:

UO+UI vO+,vl

1(-2-'-2-); (4)

,uo+u, uo+u, vo+v, vo+v,

Zif(-2-,v o)+/(-2-,v,)+/(u o'-2-)+/(u"-2-»

-t

if (uo.vo)+1 (uo,v,)+1 (u

"v

0>+1 (u"v,»

In the special case that I (uo,v) happens to be linear in v for v<l>v5v" we note

that

(3a) reduces to:

vo+v,

I (u ' - 2 -

_Ft

if (u ,vo)+/(u,v,) (3a')

(u-uo) vo+v,

I (21 (u "-2-}-1 (u"vo}-I (u"v

,»).

(u,-uo)

Similar simplifications can be made if

I

is linear for other parameter values.

Let us examine the function

I

a bit closer. Firstly, we note that

I

is a solution of the partial differential equation:

a

2

a

2 _{•• ••}

(au'

+ av2)f (u,v);(I-u)go+ug,+(I-v)ho+vh, (5)

subject to boundary conditions (I). Indeed:

a

22 I (u ,v);lim

S-V

(u +2h ,v)+1 (u ,v}-21 (u +h ,v»

au

&-.0

=Iim S-V(u+2h ,v)+1 (u ,v}-21 (u+h ,v»

&....0

=lim

S-V

(u+2h,v)+1 (u,v}-21 (U+2:+u ,v»

&....0

( using (3b) and setting uo=u; u,=u+2h; vo=O; v,=1 )

=lim S-2«I-v _{)(2ho(u+h }-ho(u }-h o(u+2h»+} &-.0

v (2h ,(u+h}-h ,(u}-h ,(u+2h»)

.. ..

=(I-v )ho(u )+vh,(u),

and similar for

a

2

2/ (U,V). Note, by the way, that the factors (I-v) and v are irrelevant lor the above

av

proof; they might be replaced by any arbitrary functions of v. The solution of the partial differential equation is, in matrix notation:

f

(u ,v)= -

[I~U]T rgO~V) ~~i~? ~',i~?]

[i=!v]

u

l~,(V)

ho(l) h,(I) v

(6)

which is exactly the definition of the bilinear blended surface according to Coons as described in

[BOE84].

Now consider the special case that ho(u) and h,(u) have the same shape and that also go(v) and g,(v)

have the same shape, that is:

(Au :05u5 I :h,(u );ho(u)+h ,(O}-ho(O» "

(Av :05v5 I:g, (u );go(u )+h o(l}-h 0(0» Then

f

(u ,v) takes the simple form:

(8) In othe.r words: the generated surface then equals the Minkowski sum of the two boundary curves ho(u)

and go(v)-ho(O):

surJace=(ho(u)+go(v)-ho(O)IQ,;uSI" Q,;vSI} (9)

which is a well suited description for a large class of sweep objects:

2. The bilinear blending algorithm

Now we arrive at the actual algorithm. Its aim is to generate a set of bicubic Bezier patches, together forming the surface

if

(u ,v) I Q,;u';; 1 " Q,;v';; 1], starting from the four boundary functions. The lauer can be thought of as being represented by a series of mB control points (Pi ,Q;), (),;;j <mB of cubic Bezier splines, specifying a piecewise cubic approximation of the boundary curve, together with mB

additional control points (Ri ) that specify the cross-boundary derivatives as depicted in the following

figure: . / , / , / \ , /

\

\

We assume that Ai :Q,;i <mB-l : _Q'i=2Pi+₁-Qi+" thus guaranteeing G 1 continuity for the boundaries, and we allow for different values of

m.

for every boundary curve.

We present several versions of the a]gorithm; every next version can be seen as a transformation of the

previous one. In the algorithms we use the indices II, 1m, Ir, ml, mm, mr, ul, urn, ur for lower left,

lower middle, lower right etcetera; the index

m

stands for 'middle' in Vm and

um •

These notions refer

7

-ul um ur

ml mm mr

11 1m lr

Furthennore we will use the addition and (scalar) multiplication in the algorithm both for numbers and for vectors; there will be no difference in notation for the two. Vectors in 3-space are denoted by the type pointJd , and tuples of 4 points together fonning one comer of a bicubic Bezier patch are denoted by the type point3d4 In the example of the above figure, the points Po, Qo, Ro and Qa+Ro-Po together fonn one pointJd4 , as do PI> QI> R I> R ,+Q,-P, and P z, Qz, R z, Rz+Qz-P z. Geometrical operations that take place on all 4 points of one pointJd4 are written as single expressions, thus indicating that the bilinear blending is applied on all control points simultanuously. For this reason the function

f

(u ,v) should be seen from now on as a function that returns for a given parameter pair (u ,v) nOl only a 3-dimensional point but also 3 additional points specifying a tangent plane to the mentioned surface, interpreted as the four control points that detennine one of the four comers of a bicubic Bezier patch.

In the first few versions of the algorithm all points consist of a record of three reals; we will mention explicitly when they should be considered as a record of three integers.

The booleans low JIOP, UP_SlOP, etcetera indicate whether the computation of the corresponding boun-daries has exhausted all input infonnation (that is the complete set (Pi, Qi' R;) for

0Si

<WIB for that boundary) and therefore whether any additional points of that boundary, needed in the computation of the blending. can be obtained by applying the De Casteljau algorithm On one of the pieces cubic Bezier spline given by one set (Pi, Qi' Pi+1> Q'i)

procedure squarefilC 0(f,,1 ITI MIl ... : point3d4; uo.u"v",v,:real; low_stop ,up_stop .Iefestop .rightjtop: boolean;); (f,,;f (uo.vo) 1\ f ,,;f (u,.vo) 1\

f MI;f (uo.v,) 1\ f ... ;f (u"v,)

Note thaI the above expressions equate four 4-tuples of 3-dimensionaJ points 10 four parameters of type pointJd4 )

var f _ IIm,f ml,f _,f _: point3d4; _um.vm: real; hor _stop .vertjtop: boolean;

begin

ir low_stop and up_stop and leftjtop and right_stop tben render _bicuhicJXJtch (f" I",f

MI! ... )

else begin

end;

u .. :;(uo+u ,)12; Vm :;(v o+v ,)12;

f .. :;(I-vm)o _ho(uo}+vm

*

_{h ,(uo}+(I-uo)o g o(vm}+U(jl< g, (vm)}

-(I-uo)o (I-v .. )o ho(O}-(I-uo)o v ..

*

h, (O}-U(jl< (I-v .. )o ho(1}-u(jI< V ..

*

h ,(I); Ie f estop :;stop (g 0);

(f .. ;f (uo.v .. ))

( ... likewise for the points with indices mr. 1m and urn ...• )

f _:;(f Im+f _+f .. +f _)12-(f MI+f .... fu.f ,,)14; (f .... _{;f (um.vm))}

vert_stop :;righestop 1\ lefestop; horjtop:;low_stop "up_SlOP;

squarefilCO(f",fImI .. 1 _,U"'Um,vO,vm' lowjtop ,hor jtop .Ieft_stop ,verestop); squarefilCO(f 1m ,fl,,f _,f _.um,u"v",vm'

low jlOp ,hor jtop ,vert_stop ,right_stop ); squaref ilC O(f .. ,f .... ,f MI,f _ ,"o,U .. ,v ... v"

hor _stop ,up_stop ,1efestop .verestop); squarefilCO(f .... ,f _,f _ I .... U .. ,U "v .. ,v"

hor _stop ,up_stop .vert_stop ,righL_stop);

end; {squarefilC O}

begin

end.

... {the conlrolpoints for the functions

f

o ....

g, gel defined} ...

9

-The function render _bicubicJXllch renders a bicubic Bezier patch given its four 4-ruples of control points. As we will see in section 3, this procedure also operates according 10 the successive subdivisions scheme.

The function stop rerums a boolean indicating whelber there are yet unused inpUlpoints on the argu-ment boundary.

The above version of the algorithm takes 576 floating point multiplications for one recursion step. This can be reduced 10 192 by observing that Ibe terms in the second lines of the expressions for

f

mJ ,

f

~

,

[1m, [_

are bilinear in u and v. This means that these terms can be computed by successive subdivisions. Let gil. g". gul. g ... represent Ibe values of the respective terms in

the four

comers of Ibe

u -v -

rectangle that corresponds wilb Ibe surface being drawn. Then we have for the values in Ibe mid-points of Ibe edges and Ibe midpoint of the rectangle:

g"", :=(g ... +gul )/2; gmJ :=(gul +gll )/2; glm :=(gIl+g,,)/2; g~ :=(g"+g,,, )/2;

g_ :=(gIm+g"",+gmJ+g~ )/4;

Now if we set [.=k.-g. for a=ul.um.ur,ml,mm.mrJI,lmJr, then Ibe next version of the algorithm (square[ilC I) is equivalent with Ibe previous version:

procedure square[ilC I(kll.gu);",g .. ,k.,.,,gul.k.u,g,,,: point3d4; uo.U"VO,Vl: real; lOW_SlOP .up_SIOP ,le[cstop ,righcslop: boolean;);

(ku-g ll=[ (uo.v0) "

gll=(l-uo)(1-v olh o(O)+(1-u o)voh I(O)+(U o(l-v olho(l)+u ovoh 1(1) etcetera} var k"",.g ... );Im,glm);mJ,gmJ);"",g"");"",,,g_: point3d4; u..,v .. : real;

begin

ir

up_stop and low_stop and Ie/estop and righestop tben

render _bicubic...palch (ku-gll );,,-KIT );"rg .. );.,-g",) else

begin

end;

.... :=(u o+u ,)/2;

v .. :=(vo+v,)/2;

k..J

:=(1-v ..

>0

hcf.u,,)+v ..

*

h, (uo)+(1-u

0>0

g o(v .. )+u"" g ,(v .. ); Ie/estop :=stop (g 0);

gmI:=(g .. +gll)/2; {kmI-gmI=/ (uo,v .. ) /\

gmI =(I-u o)(1-v .. )h o(O)+(I-uo)v .. h, (O)+(I-v .. )uoho( I)+u oV .. h, (I)) ( ... .Iikewise for the points with indices

mr. 1m

and

urn ... )

k....

:=(k .. +k .... +k"" +k_ )12-(/C .. +k", +ku+k" )/4; g""":=(g .. +g .... +gmI+g_)/4;

{k"""-g",,,,=/ (u .. ,v .. ) /\

g""" =(I-u .. )(l-v .. )ho(O)+(I-u .. )v .. h, (O)+(I-v .. )u .. hcf.l)+u .. v .. h ,(I)) verestop :=le/estop /\ righestop;

hor_stop:=low_stop /\ up_stop;

square/ilC I(kll,gll);" .g");mI,g"");"",, .g""" .uo.u ... vo.v .. , low_stop ,hor _stop .Ie/estop .verestop );

square/ilC I(k ... ,g ... );" 08");....08 .... );_08_"' ... u"v .. v .. , low_stop ,hor _stop .vert_stop .righestop);

square/ilC l(kmI ,gmI);""" 08 .... ); .. ,g .. ); ... g ... "'o.U ... v .. ,v" hor _stop ,up_stop Je/t_stop .verestop);

square/ilC I(k .... 08"",,);_ ,g_); .... ,g .... );., ,g", "' .. ,u,.v ... v" hor _stop ,up_stop ,verestop ,righljtop);

end; {square/ilC I}

begin

~nd.

... {the controlpoints for the functions h ... g, get defined} ....

square/ilC I (2h o(O),h o(O),2h o(l).hcf.I),2h ,(O),h, (O),2h, (I),h ,(1),0,1,0,1 / alse J alse J alse J alse );

11

--The multiplications with (I-uo), (I-u,), Uo etcetera to compute the k 's;

-The evaluation of the boundary functions ho, hI, go and g, with non-integer arguments.

Moreover, it occurs that all boundary functions have to be evaluated for the same arguments many times. We now propose a data structure both for the

u's

and

v's

and for the h 's and

g's

which solves all three problems:

First, replace the u 's and

v's

by records, consisting of a bitstring and two integers:

type param = record

B: string

of

bit; I ,k: integer; end; (param

I

,'ar VO,Ut.VO.V1.U""V". : param;

We introduce a function BVAL as follows: BV AL :param-4R

lEN(U.B)

BVAL(V)=

1:

V.B [i] .2''; ;=1

where LEN (A) is the length of the bitstring A. We will maintain as an invariant P:

P: UtFBVAL (Vo) " .. " v .. =BVAL(V .. ) " LEN(Vo.B)= .. =LEN(V .. .B)=n+I"

BV AL (V ,)-BV AL (V 0)=2-' " BV AL (V ,)-BV AL (V 0)=2-- " Ai:V • .i<i:V • .B[i]=Oforae(O,I,ml "

Ai:V •. I<i:V • .B [i]=O for ae (O,I,m

I

Due to the latter two conjuncts we can also write:

U.l

BVAL (V)=1:V.B [i].2''; ;=1

(10)

(11)

Since in squarefill_O and in squarefill_1 we had initially

uo=O ,

u,=I,

vo=O,

v,=1 we have to initialize: V oB :=O;V 0.1 :=0; V ,.B :=I;V ,.1:=1; V o.B :=0; V

0.1

:=0; V ,.B :=I;V ,.1:=1; n:=O; (PI

Within the procedure body we get (

e

stands for the cOIlcatenation of bitsttings ): {PJ u~ :=(uo+u,)I2; V~.B :=Vo.B

e

1; V o.B :=V 0.B eO; V ,.B :=V ,.B eO; V~.I:=n+2;

(uo=BVAL(V o) " u,=BVAL(V,) " ~=BVAL(V~)"

LEN (V 0.B )=LEN (V ,.B )=LEN(V~.B)J

vm:=(vo+v,)/2; Vm.B:=Vo.Bel; Vo·B:=Vo.BeO; V,.B:=v,.BeO; Vm·I:=n+2

(vo=BVAL(Vo) " v ,=BVAL (V,) "vm=BVAL(V .. )" LEN (V 0.B )=LEN (V

,.fJ

)=LEN (V m.B) J

n:=n+1;

[BVAL(V,rBVAL(Vo)=2'-o "LEN(V .. .B)=LEN(VoB)=LEN(V,.B)=n+1" BVAL(V,rBVAL (V 0)=2'-' " LEN (V .. .B)=LEN(V oB)=LEN(V ,.B )=n+1 " BVAL(V,rBVAL(Vm)=2-o "BVAL(V .. rBVAL(Vo)=2-o "

BVAL(V,rBVAL(V .. )=2-' "BVAL(VmrBVAL (V o)=2-o " Ai:V •. ki:V • .B[i]=OforaE{O.l.mJ "

Ai :V • . ki :V • .B [i]=O for aE {0.1.m J J

After the invocation of square fill with VI. V .. • respectively V .. • Vo as arguments instead of V" Vo. we have restored P again. Moreover. we note that n equals the recursion level; also note that we have not yet specified the meaning of the record fields V.k . In our next version of the algorithm. we will not use "0 .. v, as program variables any more. With these changes. we have (leaving out the

assertions of squaref

ilC

0 and squarefilC I ):

procedure squarefilC2(kl/.gl/.k".g".k.J .gw},k,...g.,.: point3d4; Vo.V,,vo.V,:param;

n: integer;

low_stop .uP_stop ,lefcstop .right_slop: boolean;); {P ,.., BVAL (U ,rBVAL (V o)=2->O " BVAL (V,rBVAL(Vo)=r' " LEN (V 0.B )= .. =LEN (V,.B )=n+I)

var k~.g_.klm .g'm .kml.gml.k~ .g~.kmm .gmm: point3d4; U"..Vm:param;

begin

13

-if up_stop and low_stop and left_stop and rightJtop

tben

render _bicubiC...JXlIch (ku -gu .I<"-g,, .k",-goJ .k.-g",)

else begin

end;

UM.B:=Uo-B GlI; VM.B:=VO.B GlI; U 0.B :=U 0.B Gl 0; Vo-B :=Vo-B Gl 0; UI.B:=UI.BGlO; VI.B:=VI.BGlO; VM.I:=n+2; UM.I:=n+2;

.,:=n+l; _{{BVAL(UM)-BVAL(Uo)=2-" ... "BVAL(VI)-BVAL(VM)=2- "}

Ai:U • .I<i:U • .B{i]=Oforae{O,I,m} "

Ai:V •. l<i:V • .B[i]=Oforae{O,l,m]]

""" :=k30mpute (n ,U o,Vo.); left_stop :=stop (g 0); g"" :=(goJ +gll )12;

( ... .likewise for the points with indices mr, 1m and urn .... )

k_ :=(k ... +k_ +k"" +k"" )/2-(koJ +k,.,. +kll +k" )/4; g_

:=<8 ...

+g_+g",,+g..,.)/4;

verestop :=left_stop " righestop; hor_stop:=low_stop " up_stop;

squarefilC2(ku ,gll.k ... ,g ... .k"",g"" .I<_,g_,U o,Uo. ,V .. Vo. ,n, low_stop ,hor _stop Jefestop ,vert_stop);

squarefilC2(k ... ,g ... .k" ,g".k- ,g_.k..,.,g..,.,U .. ,U "V .. V .. ,n, low_stop ,hor _stop ,verestop ,right_stop);

squarefilC 2(k"",g"".k_,g_.koJ,goJ.k_ ,g_ ,U o,Uo. ,V .. ,V lon, hor _stop ,up_stop Jefestop ,verestop);

squarefilC 2(k_ ,g_.k"" ,g..,..k-,g ...

.k.

,g", ,Uo. ,U I ,V .. , V l.n , hor _stop ,up_stop ,vert_stop .rightJtop);

begin

end.

... (the functions ho .. g, get defined) ...

U 0.B :=O;Uo.I:=O;

U ,.B :=I;U

,.I

:=1;

VoJJ :=0; V

0.1

:=0;

V ,.B :=I;V ,.1:=1;

squareJilC 2(2h o(0).h o(0).2ho(l).ho(l).2h I (O).h, (0).2h, (1).h ,(l).U o.U 1. V .. V 1.0. J alse

J

alse

J

alse

J

alse );

It is straightforward to see that the function kJompute should be: function kJompUie (n :integer; U.V :param ;):point3d4;

(n+!=LEN(U.B)=LEN(V.B) " AJ:V.I<i:V.B [i]=O" Ai:u'/<i:U.B [i]=O) var K _.hd.gd: point3d4;

i :integer; begin

K :=ho(BVAL (U))+go(BVAL (V)); hd :=h ,(BVAL (U)}-ho(BVAL (U)); gd:=g, (BVAL (V)}-go(BVAL (V)); i :=1;

wbile i!>V.land hdoO do begin

end;

if V.B [i ]=1 tben K :=K +hd;

h_d:=h_dI2;

i:=i+1;

(K =(l-BV AL (V))ho(BVAL (U))+BVAL (V)h,(BVAL (U))+go(BVAL (V))))

i:=l;

while i!>U.land gdoO do begin

end;

if U.B [i]=! then K:=k+gd;

gd:=gd12;

;:=i+l;

(K =(l-BVAL (V))ho(BVAL (U))+BVAL (V)h,(BVAL (U)) +(l-BVAL (U))go(BVAL (V))+BVAL (U)g ,(BVAL (V)))

k_compute :=K;

end; (kJompute )

Next we are in need for an implementation of the function calls like ho(BVAL (U)). We do so by intro-ducing data structures for the boundary functions h.,..g, which are well suited to supPOrt the recursive

15

-nature of our algorithm and which can contain a sufficient amount of values of the boundary functions in order to avoid a too large loss of efficiency due to their repeated computing. Such a data structure will be:

where

type boundary=array[O .. 2" lor point3d4;(p'" I]

var Ho.H"Go,G,:boundory;

H o[O]=ho(O); Ho[l]=ho(l);

Aq:15qSp:(Ai:15iS2Q

- ' : H

o

[2Q-'+i]=h

o(

2~~1»

and similar for H, .. G,. Suppose that initially rnB points of the boundary function ho have been given. These must occupy the first rnB elements of the array H 0, and therefore they get assigned the parameter values u=O, I, 0.5, 0.25, 0.75, 0.125, 0.375, etcetera. The points must be ordered accordingly. All residual elements of Ho get some non-occurring value, say XXX. For clarity, we depict the situation as follows:

In this representation, it is clear that all elements apan from the elements 0 and I have two upper neighbors. If we think of these elements somewhere along the boundary curve ho, than every element is located between its two upper neighbors. But this indicates in what way we have to proceed to find points of the boundary,

that

is 10 implement the function calls hr/..BVAL (U» etcetera, thereby either locating the proper elements in the array H 0 or assigning values 10 yet unused elements:

-Firs~ we have to translate the argument value of the functions ho(BVALO), h,(BVALO), g.,(BVALO)

and g, (BV AL 0), i.e. the param -s U and V , into the corresponding array indices of the arrays H o .. G ,. We will do so by making use of the yet unused record fields U.k and V.k. We state an additional

invariant:

Q. (U):H o[U.k]=ho(BVAL (U»" H ,[U.k]=h,(BVAL (U» Q,(V):Go[V.k]=go(BVAL(V» " G,[V.k]=g,(BVAL(V»

(12a)

Q

is initially fulfilled by setting

Vo-k:=O; V,.k:=I; Vo.k:=O; V,.k:=I;

(Q.(V

o) "

Q.(V,) " Q,(Vo) " Q,(V,)}

As caD be verified in the above scheme. we have to set in the procedure body: (Q.(Vo)" Q.(V,) "

Q,(V

_{o) "}

Q,(V,)}

VM .8 :=Vo.8

e

I; VM .8 :=Vo.B Ell I; V 0.8 :=V 0.8 eO; V 0.8 :=V 0.8 eO;

V,.8:=v,.8eo; V,.8:=v,.8eo; Vm.l:=n+2; Vm.l:=n+2

n:==n+l;

if _{V o.8=O} and V,.8=! then _Vm.k:=2 else

begin

end;

if Vo.8=O and V,.8=! then VM.k:=2

else begin

end;

(Q.(Vm) " Q,(Vm)}

-Next. test if the element H o[k

1

has a proper value. that is. if it corresponds to a real point of the boun-dary curve. If so. we are ready to return this poinL Otherwise. we have to identify the two upper neigh-bors of Ho[kl. say Ho[kt! and Ho[k,.l (due 10 the order of the recursion we are sure that they do have proper values) and then compute Ho[kl starting from these. In section 4 we will address some methods for recursive computation of points of the bour,dary curves; for now we can

think

of just applying the De Castelau algorithm.

17

-In order 10 compute k, and k" given k, we use the following procedure:

procedure upper_neighbor (k :integer;var k,,k,. :integer); begin if k mod 2=0 then begin end else begin end k,:=k div 2; k,:=k,;

while k, mod 2=0 do k, :=k, div 2;

k,:=(k,+I) div2

k,:=(k+l) div 2;

k,:=k,;

while k, mod 2=1 do k,:=(k,+I) div 2;

k,:=k, div2

if k,=1 then k,:~

end; (upper_neighbor )

The correctness of this procedure can be seen as follows: map the interval 0 .. 1 onto a circle, thereby identifying the points 0 and 1. Every next 'layer' in the array H 0 as depicted above then corresponds to what we will call an arc refinement:

0=1 0=1 0=1 0=1

3

3 4

2 2 2

If we have taken r refinement steps, then we have introduced 2' circular arcs of length

±

21t . Now for

2'

2,0'<kS:2'. we are asked 10 find the clockwise (left) and counterclockwise (right) neighbors of k. First we note that if

k

is even then the left neighbor is

k12;

if

k

is odd then the right neighbor is

(k+I)12

(This is caused by the construction of the arc refinements). Concerning the right neighbor in the case of even k, (and likewise the left neighbor in the case of odd k), we assume that the while- loop has taken

s steps. This means that the IOtai arc length that has been traversed equals

I 1 I

-21t(-+--+ ... + - - )

2' 2,-1 2'-.1'+1

1 1

=-21t( 2'0'

-2<)

k,:=(k,+I) div2

which is equivalent to an additional rotation over an arc length of

+2~.

The total arc length thus

T""' equals _27«_' ___ 1 )+2~ 2"- 2"

2"-=2.2!..

2' .

which is exactly the distance between two neighbor points. The same reasoning bolds for the case of odd k. Finally we should remember that the problem was not defined on a circle but on

a

straight line piece. which causes the point I to be in fact two coinciding points (I and 0) indicating the necessity of the final if-statement.

Resuming. we can write down the function definition HoBVALO. which yields the value of ho(BVALO):

function H oBV AL (U :param ;):point3d4;

var k, ,k, :integer;

begin

if H o[U.k ]oXXX then H oBV AL :=H o[U.k]

else begin end; upper_neighbor (U.k,k,,k,); H o[U.k ]:=hO-curveJompule (k, ,k,); H oBV AL :=H o[U.k]; end; [HoBVAL ]

Possible implementations of the function hO-curveJompule (k,,k,). which computes a point of the

curve in between the points with indices k, and k,. will be discussed in section 4. In the same way we can define functions H,BVALO. GoBVALO. G,BVALO.

Apart from the functions render _bicubic..palch

O.

SlOp

0

and the functions hO-curve_compule ...

g I_curve_compute. we now have given a complete derivation of the algorithm. It only contains addi~

tions. subtractions and divisions by 2 and 4. and when being used as a rendering algorithm for raster devices it is therefore easily transformed into an integer algorithm. The only change is to redefine the type point3d4 as a record consisting of four times three integers and to replace the

/2

and /4 operations by shifts. Note that at every recursion step we loose about one bit of accuracy. Therefore it is recom-mended to usc a longer word length than that corresponds to the screen resolution and reserve the right-most p bits as a fractional part where p equals the expected maximal recursion depth.

19

-3. The rendering of the bicuhic patches.

The problem of rendering a bicubic patch has been addressed by many authors. In [LIE87] some often used methods have been reponed. The method we use is a recursive subdivision method; it operates as follows:

a) replace the bicubic patch to be rendered by two square-cubic patches; b) replace every square-cubic patch by two bisquare patches;

c) divide every bisquare patch recursively into twO bisquare patches until the resulting bisquare patches are narrow enough to be treated as quadratic Bezier splines (parabolas);

d) divide every quadratic Bezier spline into two quadratic splines until the resulting quadratic splines have sub-pixel dimensions;

e) the resulting pixel receives a colour using a texture table lookup and an intensity using an intensity table lookup. The latter table contains the precomputed values of Phong shading intensities as a func-tion of the x and y components of the surface normal vector; these components therefore serve as table

indices.

We will now allucidate the above steps a) to e).

Ad a) and b).

The reason why we do the recursive subdivisions for biquadratic patches in stead of bicubic patches will become clear in the next section where we deal with the normal vectors; for now it suffices to state that the expressions become somewhat simpler. We have the obligation to show, however, that the sul>-divisions from steps a) and b) do not destroy the C' continuity. As follows:

In the left

pan

of the above figure we show the bicubic subdivision according to the De Casteljau algo-rithm on four controlpoints _{Po ... P 3; since t(Po+2P,+pv.}

~(PO+3P,+3Pzi'P3)

and t(P,+2P zi'P'} are colinear this subdivision preserves

C'

continuity.

In the right

pan

we show our reduction from cubic to square as it takes place twice (once in each parameter direction. in steps a) and b}. respectively}: the input again is the four-tuple Po ... P 3 and the

1 1

outputs are the two three-tuples Po.

t

(Po+3P,) • g(Po+3P,+3Pzi'P3}. and g(Po+3P,+3P2+P3}. t(P,+3P:0 • P 3. Since t(Po+3P,} •

~(PO+3P,+3Pzi'P,)

and t(P,+3P:0 are colinear. too. the C'

continuity is preserved. Ad c} and d}.

In recursive subdivision schemes, an important issue is the stopping criterion. Stopping too early with the recursion causes gaps (pixels that receive no shading). Stupping too late waists processortime. In [LIE87]. the stopping problem is dealt with by maintaining the maximal orthogonal distance between two adjacent Bezier splines between 0.5 and 1.0 times the interpixel distance. In order to compute the spacing ou in between the current curve

f

(u =Ui .v) and the next curve. Lien et al. run a series of test-ing curves in the orthogonal direction (i.e.

u

direction) at

v=O.O.

0.125 ... 1.0 and examine the step size used by these curves at the position U =Ui •

In order to reduce the amount of calculational effort. involved in this approach. we use

a

somewhat different method. We start by observing that all points from a Beziercurve are completely determined by its controlpoints. From this and from the way in which Bezier curves are being constructed by the De Casteljau algorithm. it follows that two sets of controlpoints. poop .... ,P ~ and Q 0. Q ... Q~ with the property that

(Ai: ()o;i <m :d ~(Pi .Qi )<1) where for

a

and b

e

R3

,

d~(a,b)=lim

(

L

(Iai-bi IL})L

L_oo i=x,y,z

determine two Bezier curves B p and BQ with

(Ax:x eBp:{fu>:y eBQ :d~(x;y )<1»

and, conversely.

(Ax:x eBQ:(Ey :yeBp :d~(x,y }<I»

This means that we only have to check whether the control points are close enough in order to test if the resulting Beziers as a whole are close enough. The above conclusions can be extended to hold for the case we are working in fixed point representation too. by taking enough "spare" bits to avoid a loss of significant bits at the end of the recursion. Since the maximal recursion depth is the logarithm of the screen resolution. we know in advance the appropriate number of additional bits.

Using the fixed point representation. the stopping criterion for step d) is formulated analogously. Ad e).

In order to render one pixel. we have to compute its colour and we have to examine its visibility. The latter is most easily been done using the Z-buffer algorithm. The former. since it is really in the heart of the entire rendering algorithm should be as efficient as possible. This means that we should avoid computations as much as possible. We succeeded in implementing Phong shading. texture mapping, sur-face normal vector modulation. and optional trimming and dispacement under texture-control (the latter could be useful in computer animation as a special effect). with only a few 2-dimensional table lookups (to be implemented as double-indirect addressing. using a one- dimensional auxi1iary table containing the row addresses of the main 2-dimensional table) and some additions. To do so. we arranged the colour lookup table in such a way that the 5 least significant bits represent the intensity and the (in our case 3) most significant bits represent the hue and saturation. A result from this arrangment is. that the value from the Phong-table. laying between 0 and 31. simply can be added to the value from the texture

21

-map. laying between 0 and 224 in steps of 32. The following piece of code controls rendering one pixel:

const MaxX = 1023; MaxY = 1023; SpareBits = 10; MaxMaplndex = 127; MaxIntensitylndex = 99;

type punt = record or x

.y.': integer; end;

{Punt)

type map: array[O.MaxTablelndex )[O.MaxTablelnde.t) or integer; type inttab: array[O.MaxIntensitylndex)[O.MaxIntensitylndex) or integer; extern n 1': punt;

[n.x ,fI.y =x.y-component of nonnalvector. scaled and shifted to within the range 5 •..• 94;

p.x I'.y 1'.' =x.y,z-components of pixel and associated z-vaIue. scaled by a factor 1024;

)

extern TextureMap. NormvecModulationMapX • NormvecModulationMapY.

TrimMap. DisplacementMapX. DisplacementMapY. DisplacementMapZ: map; (TextureMap contains values 0,32.64.96 •... ,224;

NormvecModulationMapX contains values -5.-4 •.. 5; NormvecModulationMapY contains values -5.-4 ... 5; TrimMap contains values 0 or 1;

DisplacementMapX .Y;Z contain arbitrary integers; )

extern I ntTab: inttab;

(IntTab contains values 0 .. 31;

extern Zbuf: array [O.MaxX )[O.MaxY) or integer; (Zbuf is the Z-buffer;

)

extern I 1,12: integer;

(11 and t2 are within the range 0 .. MaxMaplndex; they serve as indices within the several maps; in section 4 we discuss their computation.

if (TrimMap [11][12]=1) then begin

p.x :=p.x arithmeticshift]ight SpareBils; p.y:=p.y arithmetic_shifCright SpareBils; p.z :=p.z arithmeticshifCright SpareBils; p.x :=p.x+DisplacememMapX [11][12]; p.y :=p.y+DisplacemneIMapY [11][12]; p.z :=p.z +DisplacementMapZ [11][12];

if(P.x~ and p.xgfax){ and p.y~ and p.ygfaxY) then

if(p.z~uf fp.x]fp.y)) then

hegin

draw-pixel(p.x /l.y ;I"exlureMap [11][12]+IntTab [n.x+NormvecModulalionMapX [11][12]]

[n.y+NorrnvecModulationMapY [11][12]

end end

);

Zbuf fp.x]fp.y]:=p.z;

4. The computation of the normal vectors and texure mapping. The computation of the normal vectors.

We will denote the normal vector in the point f(u,v) by n(u,v). Since normal vectors sbould have

normalized lengths, it is unpractical to represent them by coordinates, since in general no three coordi-nates n"

n,., "',

represented with a finite machine precision will fulfill n.2+n,.2+n.'=a2 for some fixed a.

Instead we will represent normal vectors by two angles

a.

and a, with

n,

(Xl; =arclan- ;

n,

a, =arclan.!2..

n,

(13) (14)

Now we will discuss a method for approximating the

a.

and a, on the surface as defined above. It

consists of simply averaging the

a,

and a, at each subdivision of the surface, assuming that the

a,

and

a, in the given points of the boundaries are known:

. "'~~~ . n,.~~~

(0.;11 approxunates arclan ( ) ;a,;11 approxunates arclan (

"'u~o "'u~~

similar for (Xl:;ul .CI.y;JJ ,az;lr.o, ;lr .Ox;IU'.a,;1O" ) a,;"":=(a,,wl+a,,11 )/2; a,;"" :=(a,;u1+a,;ll )/2; a,;~ :=(a,;~ +a.;" )/2; a,;mr :=(a,,,,,+a,;/, )/2;

a,;~ :=(a,;",+a"w )/2; a,;um :=(a,;u1+a"",)/2; a,;1m :=(a,;11+0.;/,)l2; a,,1m :=(a,;11+a,;/,)/2;

23

-. nx(uo,v .. ) . n,(uo,V .. )

(0,;"" approxlInates arctan ( ) ;Ily"", approxunates arctan ( )

n. uo.v".

n,

Uo,V".

similar for ax;"" ,Ily;_,cx,_ ,Ily_ ,0,;", ,Ily;", ,ax"",,, ,Ily,....

J

This program fragment is used both in the bilinear blending pan of the algorithm and in the bicubic patch rendering parI.

It is not necessary in general to represent the

a

's in degrees or radians; since in rendering they only serve as indices into a precomputed 2-dimensional array of Phong-shading values, any proportional scale will do.

This method of course is just an approximation: it would only be exact if the equi-

u

and equi- v curves would be large circles on a spherical surface or straight lines on a fiat surface. Nevertheless, due to its simplicity and the possibility to be transformed into an integer algorithm, it can be a very useful approach if only shading information is to be computed. Indeed, human observers are generally unable to distinguish "exact" shading models, based on the exact normalvectorfield, and approximations thereof, as long as the most striking qualitative features, such as smoothness and the presence of highligths, are present. A nice example of this fact is illustrated in the wod< of Shantz et al. , where figures 3 to 6 from [SHA87] all are very convincing images of one and the same object, be it that the underlaying computational models differ dramatically.

There is, however, one aspect that should be taken good care of in using approximations for the normal vector field, as we do. The normal vector, namely, not only serves as the argument for the Phong shad-ing function, but it also distshad-inguishes between the back and the front of a surface. In situations like the figure below, we have at the viSibility contour a change of the sign of the z-component of the normal vector, and, accordingly, a discontinuous change in the observed intensity.

,

_,

,

,

,

•

Suppose that we apply the method as described above to compute the normal vector components for the

patch in the above figure, which is assumed not to be spherical. Then the change of the sign of the

z-component most likely will occur on the wrong place, thus resulting in annoying errors in the picture. The way in which we solve this problem is the following. We start by observing that within one bicubic patch sign changes of the z-component of the normal vector may occur zero times, once or twice along one equi-parameter curve. This means that inspecting the normal vectors in the corner points does not reveal the structure of the normal vector field with respect to its z-component within the patch. More-over, since the components of the normal vectors are quintic functions of the parameters, it is very hard to compute the location of the viSibility contour in advance. Things are much more optimistic, however, in the case of biquadratic patches. First of all, the computation of the visibility contour could now in principle be done. But there is another circumstance which allows for easy detection of sign changes within a quadratic patch. Consider the situation below:

P,

Here we have the 9 control points PI' .. P 9 that detennine a biquadratic patch. This patch is to be ili'/ided into two patches by the introduction of the points PI 3, p. 6 and P 7 9' From some elementary algebra it follows that the nonnal vectors in the laUer three pOints can be expressed in tenns of the nor-mal vectors in the points PI' P 3, P 7 and P 9, say n" n3, n7 and n9, and some additional cross products.

Let

Then:

n4=(P_{2 -} P 4)X(P_{3 -} ~ .),

n6=(P_{I -} P₆)X(P₂ -P 6) n,=(P, - P,)x(P 3 - P,).

nl_391roportional to(n4+n<+nl+n3+2n,). Similar expressions hold for n4_6 and n7_9'

This means that if all vectors in the right hand sides of these expressions have the same sign for their z<omponents, then so wiIl, according to the recursive subdivision scheme, all the nonnal vectors within the biquadratic patch. In practice, we only examine the vectors nl, n3, n7 and n9 because in most patches being used in 3-dimensional modelling the equality of the signs of their z-components implies the equality of the signs of the other vectors involved as well, In case not all the examined components have the same sign, we divide the patch into 4 subpatches and repeat the test, This means that we not only use the recursive subdivision scheme to render the biquadratic patches, but also to find, on the fly, the visibility contour. In practice this method works very elegant, very fast and very precise, Note that norma!isation of the normal vectors never occurs; moreover the (integer) multiplications in the cross products are the only multiplications in the entire algorithm and they only take place in

the

direct neighbourhood of the visibility contours. Typically, the number of multiplications in the entire algo-rithm is little over 4 times the square root of the number of rendered biquadratic patches.

Texture mapping,

We assume that the texture, or the nonnalvector modulation, or the trimfunction, or whatever we want to be mapped onto the surface, can be written as a function T:

T: {(p ,q)E R2 p:S: I " Q;; q:S: I J --+something (15)

where T (p ,q) is the value of some function we want to define in the point (p ,q). We define a series of finite resolutiDn representations of T, Tm,m" as follows:

25

-PI; ql

where A is an operator which returns a function value that corresponds widl the "average" value

P=Po;q=qo

on die domain {(P.q)lprfopS,P1" qrfoqS,qd. The T",l"" should be precomputed and stored for mlo m2E {O.I •.. ,m""',}. Note dlat the amount of storage required is only a factor 4 more

than

die storage of 1"'_M_ only.

A rnce way to represent . all th e ... ₁

,M,.

I S '

[O .. 2M_+1_IJ[O .. 2M~ +1__{1] as follows:}

to arrange them

in

one 2-dimensional array

2"'IIIalI_I

TM_M_(p .Q)

o

-o

With every subpatch to be processed. a pair of values of the pair of indices I I and 12 exists that point to the T -value that should be used when rendering the lower left corner of this subpatch. For the algo-rithm, we proceed as follows. Initially. 11 and 12 both equal 2",_+1__2. In case the entire original patch should be of subpixelsize. we render the entire patch (=one pixel) using one T value. namely Too(O.O). This is the upper righunost entry in the table. and it is exactly the value where the pointers 11 and 12 are pointing to. In case the patch should be split into two halves in the direction of the first parameter. the left and right subpatch both receive a new value for 11. say lefl (11) and right (r 1). The

value of r2 reamains unchanged during the spit Conversely. if the patch is split in the direction of the second parameter, then I 1 remains die same for the two subpatches and 12 is replaced by two new values for the upper and lower subpatch. say up(r 2) and low(12). Symmetry considerations lead to the

conclusion that the function Ie/tO is identical to the function upO and rightO is identical to lowO.

From inspection of the above figure we learn that

m +1 m +1

lefr (2 - -2)=2 - -4

and

m +1 Maax+1

righr(2 - -2)=2 -3. From inspection it is also clear that

Closer inspection reveals that 1ft '" +1 (A I :2 ~:s; I < 2 - -2:lefl(I)=M(I}+2(I-K(I») where (A I :2m-:s; I

<:

2m-+'-2:right(I)=lefl(I}+I) M(I)=C!nax.j: kj+l:S; I: kj); K(I)=(max

j:

kj:S; I: _{kj );} kj=(:E i: O:s; i

<j:

2m_ - \

As an example we give the functions leflO and righlO for mmax=3: lef I (I) righl (I)

o

0

0

I I I 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 0 I 9 2 3 10 4 5 11 6 7 12

8

9 13 10 11 14 12 13 15

Note that with increasing values of the recursion level, the values of I I and 12 decrease until they end up inside the lower leftmost subsquare of the table. The simplest way to implement this way of texture mapping is to precompute the functions leflO and righlO and store the values in two arrays; in recur-sive subdivion calls the indices I I and 12 can be passed by either replacing I I by lefl (I I) and right (I I) for the two subpatches or replacing 12 by lefl(12) and right(12), depending on whether the subdivision takes place into one parameter or the other parameter. The code fragment of the previous section assumed the values of I 1 and 12 to be computed as described above; note that it is in fact hardly "computing": it is merely a repeated function application, coded as repeated table lookups.

5. Recursively defined boundary functions

In their survey paper [BOE84], Boehm, Farin and Kahmann present an overview of the most frequently used curve description methods in CAGD. Since most of these are recursively computable, it would suffice in principle to refer to this paper for an idea of the implementation of

the functions hO-curve_compule .... g,_curve_compuJe. We prefer, however, to give our own derivations for the following reasons:

-consistency of representation and notation;

-a somewhat different set of specifications to start with; -the need for integer algorithms.

Our approach therefore appears to be a bit upside-down: we start by giving a form of an algorithm we think the curve generator should have in order to be easily integer-implementable and next bring in the desired (symmetry) properties that the curves should posses in order to set constraints to the algorithm;

27

-it is interesting to note that eventually we end up w-ith a family of algor-ithms for smooth interpolation which includes the well-known De Casteljau algorithms [BOE84] for four-point and three-point Bezier cwves but which have one degree of freedom more than the De Casteljau algorithms.

Due to the structure of the surface algorithm of the previous section, we are in need for functions ho....curvejompute and so on , in this section to be referred to as curvejompUle, that generate a point of a cwve gi ven its two endpoints. Moreover, we allow for some additional information, but then again this should be provided as attributes with the two endpoints.

An

obvious choice

seems

to be to give the two endpoints together with two additional points that indicate the tangents of the curve in the two end-points, as in the Bezier definition:

Let uS denote the two endpoints by

x

0 and

x,

and the two vectors, pointing from the endpoints to the additional points by

ao

and

a,.

Note that in our convention both

ao

and

al

are parallel to the tangents, whereas in the Bezier-configuration one is parallel and one is anti-parallel. The cwve segment will be given by

c(xO,x,.aO.al):

c :R'xR3XR3XR3~p

(R3), (18)

where

PO

denotes the powerset.

We demand some obvious propenies:

XO,x,E C(;t:O'>:haO,a,) (c should contain its endpoints)

c (xo,x,.aO,a,Fc (x ,,xo,-a ,,-a

0)

(c should posses 'exchange symmetry'; note the signs!)

c (Mx oNx ,Na oNa ,FMc (XO,xl.aO.a ,)

(c should be invariant under a linear transformation M)

c

(xo+v,>: l+v .ao.a lFC (XO,xl.aO,al)+v

(c should be invariant under a translation v)

c (xo+(xo,e)v ,x ,+(x"e

)v

.ao+(ao,e)v.a ,+(a he)V F

c (xo,x "a o,a

l)+(C

(x

0'>:

,.a o,a 1),e)v

(c should be invariant under parallel projection with vectors a and v)

(19)

(20)

(21)

(22)

(23)

We want an algorithm that produces an element of

c(xo,x"ao.a,),

say

x .. ,

and its associate vector

a ..

c

(xo,X,.a."a,)=c (xo,x~ .a."a~) u c(x~,x,.a~ .a')

(Note the advantage of our sign definition for

ao

and

a,

above the Bezier convention.) A simple class of algorithms that fulfill both (21) and (23) is based on:

~Hm

(24)

(25a)

In order to fulfill (22), we have to demand that P oo+P 0,=1 and P llr~P 11=0; note that no similar con-straints exist for the other coefficients for reasons of spatial symmetry, since there are no simple spatial transformations that influence the tangent vectors and leave the points themselves W1Changed, in a way as the translation transformation influences the points whiie leaving the tangent vectors unchanged. In order to guarantee a correct behavior of the sign of a~ under exchange as expressed in (20) we have to demand:

p~;~l=rl o)pl;~l

(26)

ao

Lo

-I

-a,

ad

-ao

We find from (26) as additional constraints: P crr+P03=O, P 01=

~

and P ,z=P 13 , thus our matrix P

looks:

P _

r

112 112 -PooP 03] (27)

-l-P

Il

P

Il P'2 P'2

In order to study the nature of the curve that is generated with lItese matrices, we assume it to be parametrized as follows:

X:{IERIQ5I:S;ll~R3 (28)

a

(I )=4>(I)i (t); (j>(I) is a scalar function x (O)=xo;

x(l)=x,;

From all possible pararnetrisations, we will only consider those that exhibit the following property: let,

at a given recursion level p,

XIP]O=X(lo) 1\ xIP],=x(t,) =>

xIP]~=x(I~I').

This means that we assume that the pararnetrisation is such that in every recursion step the curve

[x (r) I 10:S; I :s; I,) is split into two halves:

10+1, 10+1,

{X(I)llo:S; I:S; Itl = [x(I)llo:S; I:S; - 2 - l u

[X(I)I-2-:S; I:S;

ttl

where the two halves thus correspond with equal intervals in parameter space. In order to do a follow-ing recursion step we have to replace I by

t

(1-10> and by

t

(I

1~ll),

respectively, irrespective of the values of 10 and I" in order to set up for the next two halves. A consequence of this is that we have to scale the tangent vectors (that are proportional to the derivatives of lite curve in lite endpoints of lite considered piece) accordingly.

29

-The total computation of Xm and aM' reckoning with subsequent recursive calls, therefore should read:

[:+m

09)

ao:=o,j2;

a

1:=0/2;

am

:=a",.I2;

From this chosen parametrisation we learn that for every recursion level p there exists at least one 1 such that:

x fp 1m =x (I); xfplo=x (I -li); X fp II=x (1+8); (30)

a fp1m =2&p(1 )x (I); a fplo=2&P(1 -li).i (1-<'5); afpl,=2&p(I+8)X(I+O),

where 5=TP, and in the first step we have p=1 and therefore 1=

~

as only value for I; in the second step we have p =2 and possible values for 1 are

!,

~

and so on, The factor

28

in the expressions for a accounts for the above mentioned scaling; note that

20=1

when p=L

We can perform Taylor expansion round

8=0

for (30) and next substitute into (29). Equating all terms with equal powers of

Ii

yields the following results:

( ;1 )" x (1)=0 for n

>3,

thus the resulting curve is a cubic polynomial; cjl does not depend on 1 ;

Pll= cjl(I-2Pd;

I

Pm= 8cjl; P ll=

-64>P

12;

Combining constraints (3Ic) and (3Ie) we find:

I P

lr

-'4;

3<1> Pll

=2';

(3Ia) (3Ib) (3Ic) (3Id) (3Ie) (313) (32b) This means that we still have one degree of freedom left: we can choose an arbitrary value for cjl , thus defining an entire class of algorithms. The four-point De Casteljau algorithm, with an opposite sign for a I for computing Bezier curves [BOE84] is a member of this class: it has cjl=

~

and thus

3

Pm=g; I Pll

='2'

(33)

as is easily verifIed from the geometrical constructions which form the base for the De Casteljau algo-rithms.

(34) which displays nicely the degree of freedom. We observe that in the special case that

(o,+oo)=2$(x,-xo) the curves reduce to parabolas. In this case we can drop consb'aint (3Ie). since this constraint originated from an expression involving the third derivative of X(I) only. The well-known three-point De Casteljau algorithm is an example of the latter case: it has

I

~;

~

I P o

,;"2;

P,:!",O; I P11

;"2;

I Poo

;'4;

but again it is only one member of a class of algorithms, distinguished by different values of

$.

6 discussion and conclusiol)s

By dealing with a swface algorithm and a boundary algorithm in one unified approach we have increased the applicability and versatility of bilinear blended swfaces. Because of the fact that the boundary curves can be arbitrarily shaped, we have the impression that relatively complex objects can be described using relatively few swfaces. This means that the problem of the continuous transition between one swface and its neighbors looses a bit of its importance, because there are simply less of these transitions. Still, by using a linear blending scheme also for the cross-boundary contrOl points and using a bicubic patch rendering algorithm for the lowermost recursion levels, we succeeded in deriving an integer algorithm for a true

C'

scheme.

Some other topics for future research include:

-Development of a more general theory concerning recursively defined swfaces and curves; -An error-analysis of the integer versions of the algorithms.

Acknowledgement

The author wishes to thank Mr. Anton Smeets of the department of Mathematics and Computing Sci-ence of the EUT for the proof of the correctness of the algorithm upper_neighbor.