from the modified spotlight to all other elements. The delta
`emission is determined according to the change in the direc-
`tionality parameter, in this case, E;.Q'—"2i12(cos""(¢v,‘+Av)-—
`cos""(¢.,k)) where duh is the angle between the original di-
`rection vector of the light and the direction of the element
`and ¢vk+A.,
`is the angle between the new spotlight direc-
`tion vector and the direction of the element. Subsequent
`shooting steps proceed in the normal fashion in order to
`handle multiple bounce effects. The same technique can be
`used when the distribution pattern parameter m. is changed
`as illustrated in figure 3.
`In this case the radiosity cast is
`Ek( n -l-fin-l-1 c0sn,‘(¢vk)__ n2-l-1 c0sn,,(¢vk))_
`The cost functions that measure patterns of light or sub-
`jective impressions are defined in terms of perception. The
`partial derivatives of the functions examining lighting pat-
`terns with respect to light emission Eh are:
`
`afbrightnesa = _Zi
`3E1:
`2.. A.‘
`
`a_fnon—unI'form
`
`3E1;
`
`Z: (P‘“’9r' _
`
`2:‘. A.‘
`(BISQE:
`Z.‘ A‘
`
`' —-
`
`afperipheral
`BE).
`
`__ E: 38119:,
`2:‘./1.‘
`
`' _
`
`A]
`X2111,‘
`
`The partials of the subjective impressions are just a lin-
`ear combination of the partial derivatives of f1,r.'gm,.”,,
`_fnon—un:'form, and _fpert'pheral-
`The partials 8P,-/3E1; are derived by differentiating equa-
`tion 2 giving,
`
`BP,‘ __
`
`35. ‘ 1°
`
`.,
`
`aa 83,’
`
`[3, an. + an
`
`Ba
`
`(5)
`
`the adaption level which can be approx-
`is
`where or
`Iog1o(B.'/10,000)A.')/ ::‘.A.-,
`7 is aa a:
`imated by
`Iog1o(B.'/10,000) + bb,
`and c is 0.4l0g1o(B.'/10, 000) —
`In(10)(0.8a + 2.6).
`If the the adaptation level is assumed constant with re-
`spect to a change in emission Eh, 80//6E). = 0, otherwise
`
`A,‘
`__
`30/
`6E1. " B,zn(1o)§j,(A.-) .95..
`
`3.4 Optimization
`
`The optimization process uses the BFGS algorithm, which
`evaluates the objective function and gradient at a current
`step in the design space in order to compute a search di-
`rection. Once a search direction is derived, a line search
`is performed in this direction. Each step in the line search
`involves a reevaluation of the objective function, hence a
`reevaluation of the element radiosities which are displayed,
`allowing the user to watch the progress of the optimization.
`This process is repeated until the system has converged to
`a minimum.
`
`BUNGIE - EXHIBIT 1018 - PART 2 OF 4
`
`SIGGRAPH 93, Anaheim, California, 1-6 August_1993
`
`
`
`“#24 1~:..<oos'=« ¢v....- °°S"“¢v.)
`
`\
`
`\
`
`'
`‘
`
`'
`
`L
`
`v,_ v.+Av\
`
`A3;
`
`IIIIIIIIIIIIIIII4
`.l
`
`Figure 2: Estimation of 831/ 8Vk by shooting a delta emis-
`sion from source k.
`
`1;
`
`n +1
`"!2_"'
`1
`A
`(2x:2_n+_
`
`Dk¢vk
`
`,...,+A..,, _g.3~_1c°S,,M
`
`\
`
`AB,-
`
`IIIIIIJIIIIIII
`
`1'
`
`
`vk‘
`'
`fig. wsn.+An¢vx
`
`
`
`I l
`
`Figure 3: Estimation of BB;/Bnk by shooting a delta emis-
`sion from source k.
`
`by finite differences. A small “delta” emission, AE, is shot
`from the variable emission light source as indicated in fig-
`ure 1 and allowed to interreflect. The iterative shooting
`operations are very rapid since the links representing the
`form factors are precomputed during the baseline rendering.
`The result of shooting a small amount of energy through
`the network of links results in an effect on each element
`radiosity, 133;, thus providing all the derivative estimates
`AB,/AE. If the only free variables in the optimization are
`light emissions,
`these influence factors need only be eval-
`uated once, due to linearity. On the other hand,
`if any
`spotlight directionality or element reflectance is allowed to
`be variable, light emission influence factors must be updated
`each iteration.
`
`The partial derivative of the objective with respect to a
`variable element reflectivity is handled in a similar fashion.
`The element reflectivity p;. is adjusted by a small delta Ap.
`The effect on all other elements can be evaluated by “shoot-
`ing” the unshot radiosity due to the change in reflectivity:
`B;,Ap. As with light sources, several shooting iterations may
`be necessary to account for multiple bounce effects. Once
`convergence has been achieved, the effect of Ap on element
`radiosity AB,‘ is available and the influence factor estimate
`AB,‘/Ap can be recorded.
`Influence factors for spotlight directionality variables, V;,
`and n;,, are also approximated through finite differences.
`For example, a small change, AV, can be made to the di-
`rection vector V1. and the effect on each element radiosity
`can be determined by a series of shooting steps. The first
`shooting step, illustrated in figure 2, shoots a delta emission
`
`152
`
`BUNGIE - EXHIBIT 1018 - PART 2 OF 4
`
`

`
`
`
`COMPUTER GRAPHICS Proceedings, Annual Conference Series, 1993
`
`
`
`4 Experiences and Results
`
`The first implementation of the Radioptimization system al-
`lowed an objective function based only on photometric mea-
`sures and did not take into account the psychophysical prop-
`erties of lighting. The system could successfully optimize
`fighting but required quite a bit of unintuitive “tweaking”
`of the objective function weights in order to achieve light-
`ing that had the right subjective appearance. These early
`experiences led to the investigation of the psychophysical
`objective functions.
`Figure 6 shows the effects that the subjective impressions
`have on an optimization. The top image constrains the table
`to have a small amount of illumination while conserving en-
`ergy and creating an overall impression of visual clarity. To
`improve efficiency the optimization was run at a low resolu-
`tion on a simplified model, without the chairs and television
`set. The optimization process took 1 minute and 21 seconds
`on an IBM Model 550 RISC System 6000. The bottom im-
`age has the same design goals as the top image except that
`it tries to elicit an impression of privateness. This optimiza-
`tion took 2 minutes and 11 seconds.
`It took two or three hours of performing design iterations
`before developing an intuitive “feel” for the optimization
`process and the effects of the weights on the objective func-
`tion. One of the problems with the design cycle is that there
`may be local minima of the specified objective that are vi-
`sually unattractive. For example, in addition to the design
`goals mentioned above for figure 6, we needed to add an
`additional constraint limiting the illumination of the ceil-
`ing because pointing the lights directly at the ceiling was an
`optimal way of increasing the overall brightness of the room.
`One drawback of the system at this point is that it is not
`fast enough to allow a highly interactive feedback cycle for
`complex models. However since the system allows a designer
`to think in terms of their own design goals, it requires fewer
`design iterations to achieve the desired result.
`
`5 Conclusions
`
`This paper has presented a new method of designing illumi-
`nation in a computer simulated environment, based on goal
`directed modeling. A library of functions were developed
`that approximate a room’s success in meeting certain light-
`in g design goals such as minimizing energy or evoking an im-
`pression of privacy. The objective functions were developed
`through an experiment in which subjects ordered a set of im-
`ages according to a particular impression. Processing this
`data with INDSCAL, showed a correlation between quanti-
`tative lighting patterns and subjective measures of visually
`clarity, pleasantness, and privacy. Once the lighting design
`goals have been set, the software system searches the space
`of lighting configurations for the illumination pattern that
`“best” meets the design specifications. The system absorbs
`much of the burden for searching the design space allow-
`ing the user to focus on the goals of the illumination design
`rather than the intricate details of a complete illumination
`specification.
`.
`The radioptimization system explores only one possi-
`ble path in the application of optimization techniques to
`image synthesis design problems. Constrained optimiza-
`tion techniques may be more suitable than the uncon-
`strained penalty method technique used here when the de-
`sign goals must be satisfied precisely. Discrete optimization
`methods may be appropriate in some instances, for exam-
`ple when emissivities are constrained to a finite set, e.g.
`
`Geometric properties of the
`{60 Watts ,10O Watts
`model, such as the position of the lights or the size and po-
`sition of the windows, could be allowed as free variables.
`More general image synthesis methods could be applied to ,
`account for non-diffuse effects such as glare.
`
`Acknowledgments
`Pat Hanrahan, Larry Aupperle and David Salzman pro-
`vided the radiosity software used as a basis for this work.
`Greg Ward offered suggestions on useful objective functions
`for lighting design. Shinichi Kasahara participated in dis-
`cussions about the work as it developed.
`The first and second author’s work was supported by NFS
`(CCR-9210587). All opinions, findings, conclusions, or rec-
`ommendations expressed in this document are those of the
`authors and do not necessarily reflect the views of the spon-
`soring agencies.
`
`...,,,....._..m:.
`
`References
`
`[1] J. J. Chang and J. D. Carroll. How to use INDSCAL:
`a computer program for canonical decomposition of N-
`way tables and individual differences in multidimen-
`sional scaling. Technical report, Bell Telephone Labo-
`ratories, 1972.
`
`[2] M. F. Cohen, S. E. Chen, J. R. Wallace, and D. P.
`Greenberg. A progressive refinement approach to fast
`radiosity image generation. Computer Graphics {SIG-
`GRAPH ’88 Proceedings}, 22(4):75—82, July 1988.
`
`[3] M. F. Cohen and D. P. Greenberg. The hemi-cube: A
`radiosity for complex environments. Computer Graph-
`ics (SIGGRAPH ’85 Proceedings), 19(3):31—40, July
`1985.
`
`[4] J. E. Flynn. A study of subjective responses to low en-
`ergy and nonuniform lighting systems. Lighting Design
`and Application, Feb. 1977.
`
`[5] J. E. Flynn, C. Hendrick, T. J. Spencer, and 0. Mar-
`tyniuk. A guide to methodology procedures for mea-
`suring subjective impressions in lighting.
`Journal of
`the IE5, Jan. 1979.
`
`[6] J. E. Flynn, T. J. Spencer, 0. Martyniuk, and C. Hen-
`drick. Interim study of procedures for investigating the
`effect of light on impression and behavior. Journal of
`the IE5, Oct. 1973.
`
`[7] P. E. Green, F. J. Carmone, Jr., and S. M. Smith.
`Multidimensional Scaling Concepts and Applications.
`Smith, Allyn, and Bacon, 1989.
`
`[8] P. Hanrahan, D. Salzman, and L. Aupperle. A Rapid
`Hierarchical Radiosity Algorithm. Computer Graph-
`ics (SIGGRAPH '91 Proceedings}, 25(4):197—206, July
`1991.
`C
`
`[9] J. T. Kajiya. The rendering equation. Computer Graph-
`ics (SIGGRAPH ’86 Proceedings}, 20(4):143—150, Aug.
`1986.
`
`[10] H. N. McKay. Energy optimization and quality lighting
`design. Lighting Design and Application, Mar. 1986.
`
`[11] P. Y. Papalambros and D. J. Wilde. Principles of Op-
`timal Design. Cambridge University Press, Cambridge,
`England, 1988.
`
`
`
`
`
`153
`
`
`
`

`
`
`
`8_|§-‘:(_3_FiAPH 93, Anaheim, _Ca|if_ornia, 1-6 August 1993
`
`[19]
`
`[13]
`
`[14]
`
`[15]
`
`P. Poulin and A. Fournier. Lights from highlights and
`shadows. 1992 Symposium on Interactive 3D Graphics,
`pages 31-38, Mar. 1992.
`
`W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T.
`Vetterling. Numerical Recipes. Cambridge University
`Press, New York, 1986.
`
`J. B. Rosen. The gradient projection method for non-
`linear programming, part i: Linear constraints. SIAM,
`s:131—217, 1930.
`
`J. B. Rosen. The gradient projection method for non-
`linear programming, part ii: Non-linear constraints.
`SIAM, 9514-532, 1961.
`
`llfil
`
`P. C. Sorcar. Architectural Lighting for Commercial
`Interiors. John Wiley and Sons Inc., 1987.
`
`[171
`
`[13]
`
`S. S. Stevens and J. C. Stevens. Brightness function:
`Effects of adaptation. Journal of the Optical Society of
`America, 53(3), Mar. 1963.
`
`J. Turnblin and H. Rushmeier. Tone reproductions for
`realistic computer generated images. Technical Report
`GIT—GVU-91-13, Graphics, Visualization, and Usabil-
`ity Center, Georgia Institute of Technology, July 1991.
`
`.I—_. .-jazuunnuuu
`
`nu-u-u rlxunuluual
`an-u .pqnunuuIunI.|Iu
`an-5 .-.n-nuns uunuln
`
`Figure 5: Sample interface which allows the user to set the
`weights of the objective and [or specify constraints.
`
`Figure 4: Computer generated rooms used to test subjects
`on which illumination patterns illicit particular subjective
`impressions.
`
`Figure 6: The top image constrains the table to have a small
`amount of illumination while preserving energy and creat-
`ing an overall impression of visual clarity. The bottom image
`also constrains the table to have a small amount of illumi-
`nation while preserving energy. In addition, it trys to create
`a feeling of privacy.
`
`154
`
`

`
`
`
`COMPUTER GRAPHICS Proceedings, Annual C’pnference_$eries, 1993
`
`A Hierarchical Illumination Algorithm for Surfaces
`with Glossy Reflection
`
`Larry Aupperle
`
`Pat Hanrahan
`
`Department of Computer Science
`Princeton University
`
`
`
`Abstract
`We develop a radiance formulation for discrete three point transport,
`and a new measure and description of reflectance: area reflectance.
`This formulation and associated reflectance allow an estimate of er-
`ror in the computation of radiance across triples of surface elements,
`and lead directly to a hierarchical refinement algorithm for global
`illumination.
`
`We have implemented and analyzed this algorithm over surfaces
`exhibiting glossy specular and diffuse reflection. Theoretical growth
`in light transport computation is shown to be O(n+k3 ) for sufficient
`refinement, where n is the number of elements at the finest level
`of subdivision over an environment consisting of k input polygonal
`patches —- this growth is exhibited in experimental trials. Naive
`application of three point transport would require computation over
`O(n3) element-triple interactions.
`CR Categories and Subject Descriptors: 1.3.7 [Computer Graph-
`ics]: Three—Dimensional Graphics and Realism.
`Key Words: adaptive meshing, global illumination, radiosity, ray
`tracing.
`
`1
`
`Introduction
`
`A major open problem in image synthesis is the efficient solution of
`the rendering equation. Radiosity methods have been quite success-
`ful over environments containing surfaces that exhibit only diffuse
`reflection. Unfortunately, very few materials are purely Lamber-
`tian reflectors, and efficient solution techniques have not yet been
`developed for more general specular or glossy reflection functions.
`The rendering equation is an integral equation, and the solutions
`to complicated integral equations are generally obtained using either
`Monte Carlo or finite element techniques. Monte Carlo algorithms
`sometimes go under the name of distributed or stochastic ray tracing
`and are the most commonly employed in computer graphics (e.g.
`see [4, 5, 9, 12, 16]). Monte Carlo techniques have the advantage
`that they are easy to implement and can be used for complicated
`geometries and reflection functions. Unfortunately, their disadvan-
`tage is that they are notoriously inefficient. The second approach,
`the finite element method, has been very successfully applied to
`the rendering equation under the radiosity assumption, but has only
`begun to be employed in the general case, and with limited success.
`For example, lmmel et al. [8] discretized radiance into a lattice of
`cubical environment maps, and solved the resulting system. More
`recently, Sillion et al. [13] used a mesh of spherical harmonic func-
`tions to represent radiance, and solved the resulting system using a
`shooting algorithm.
`
`Permission to copy without fee all or part of this material is granted
`provided that the copies are not made or distributed for direct
`commercial advantage, the ACM copyright notice and the title of the
`publication and its date appear, and notice is given that copying is by
`permission of the Association for Computing Machinery. To copy
`otherwise, or to republish, requires a fee and/or specific permission.
`
`© 1993
`
`ACM-0-89791-601-8/93/008/0155
`
`$0150
`
`
`
`
`
`
`
`There are many ways to parameterize the rendering equation, and
`each leads to a different choice of basis functions. In the transport
`theory community two techniques are common: directional sub-
`division (the method of discrete ordinates or SN), and sphericad
`harmonics (PN). These two techniques roughly correspond to the
`methods of lmmel et al. and Sillion et al., although many interest-
`ing variations are possible. Our approach is somewhat different, and
`based on Kajiya’s original formulation of the rendering equation [9].
`Under this formulation, the rendering equation is expressed in terms
`of three point transport. That is, the kernel of the integral expresses
`the transport of light from a point on the source to a point on the
`receiver, via a point on a reflector. Given this formulation, the three
`point rendering equation can be discretized over pairs of elements
`to form a linear system of equations. Solving this system yields the
`radiance transported between elements. Note that this approach is
`very similar to the radiosity formulation.
`The problem with finite element methods is that the matrix of
`interactions is very large for interesting environments. For a given
`environment of k input polygonal patches containing 11. elements
`at the finest level of refinement, the three point discretization that
`we are proposing generates an 71.3 matrix of interactions. However,
`in this paper we show that we can accurately approximate the n3
`reflectance matrix with O(n + k3) blocks, in a way very similar
`to our recent hierarchical radiosity algorithm [7]. In that paper we
`showed how the n2 form factor matrix could be approximated with
`O(n + k2) blocks, resulting in a very efficient algorithm in both
`space and time. Although the results presented in this paper are
`preliminary, we believe a hierarchical finite element approach along
`these lines will ultimately lead to a fast, efficient algorithm.
`In the following section we describe our application of the fi-
`nite element element method to the three point rendering equation,
`yielding a radiance formulation for discrete transport. In Section 3
`we present a simple adaptive refinement algorithm for computation
`over this formulation, and the iterative solution technique employed
`for the actual calculation of transport. In Section 4 we discuss our
`implementation of the algorithm over glossy reflection, and in Sec-
`tion 5 we present some experiments and results. An appendix to
`this paper contains details of our error analysis for discrete transport
`under the glossy model.
`
`2 Discrete Three Point Transport
`
`The algorithm presented in this paper operates through two func-
`tions: refinement of the environment to form a hierarchy of discrete
`interactions, patches and elements, and the actual computation of
`illumination over this hierarchy.
`In this section we develop the basis for both discretization and
`transport. We derive a radiance formulation for three point transport,
`and a new measure and description of reflectance, area reflectance.
`This radiance formulation and associated reflectance provide a natu-
`ral criterion for discretization under illumination and reflection, and
`allow both the computation of radiance across triples of individual
`surface elements, and the expression and computation of all light
`transport over all surfaces.
`
`155
`
`

`
`
`
`
`
`kxjxi
`
`ixjxl
`
`kxjxl
`
`SIGGRAPH 93, Anaheim, California, 1-6 August 1993
`
`
`
`Figure 1: Geometry of Reflection
`
`Figure 2: A Reflection Product
`
`2.1 A Radiance Formulation for Three Point
`
`‘Transport
`When computing and imaging illumination within an environment,
`we are interested in the transport of light from surface to surface —
`it is this interaction of surfaces that characterizes illumination, in the
`absence of participatory media. Reflection within an environment
`may thus be naturally expressed over triples of surfaces. Consider
`surfaces A, A’ , and A" (Figure 1) — we will examine the transport
`of light incident at A’ originating at A and reflected toward A".
`Let Lu: and to; be the solid angles subtended at point :5’ by A and
`A”, respectively. Consider differential solid angles at 63$ and 655. —
`by definition of the bidirectional reflectancerdistribution function
`(BRDF), f, [11], the radiance L(cI3:.) along 631. due to illumination
`through solid angle wfi is:
`
`L(w:)= / rr<a2,w:)L(w:>cos0£dw£
`
`cu’.
`
`Integrating this expression over wfi , and introducing cos 0; , we have:
`
`[ L(tZ2'f.) cos0',dw', =
`Ir
`I I f,(632,<23’,)L(d3£) cos 9: cos 9',dw£dw',
`
`by definition of the diffuse form factor, F,-;..
`We may similarly discretize A as A.-, and rewrite the right side
`of equation (1) as:
`
`L,~J-‘f I
`
`A.’ Aj Ab
`
`f,(a:,:5',a:")G(a:,:5')G(:5',:5")da:"d:1;'d:5 =
`
`Z7TLijAiFijRijk
`
`where R;,1. is defined such that
`
`7|'AiF:'jR«z'jk =
`
`‘/l I
`
`A;
`
`A1‘ A],
`
`fr(:5, :5’, :5")G(:5,z')G(:5', a:")d:5"d:5'd:5
`
`Note that, by the symmetry of fr and G:
`
`AiFijRijk = AkFlcjRkji
`
`We thus have:
`
`7rL,~,.A,-F,-,. = 7r EL.-,~A,.F,.,-R,.,«.~
`
`7' Z Lo‘Ai1“3'I=Rkji
`
`We may then reparameterize over A and A" to yield:
`
`by the reciprocity of form factors, and thus:
`
`Ljk = ZLajRkji
`
`The three dimensional character of Rkji over indices 1‘, j, It leads
`naturally to a three dimensional matrix formulation for the above
`system. Consider a product over an n x n x n R;,,~,' “matrix” and
`an n x n X 1 L.-,- matrix producing an n x n x 1 matrix of reflected
`radiances, as shown in Figure 2. Note that the R1.,-,- matrix is of size
`O(n3) — the hierarchical method discussed in subsequent sections
`of this paper addresses more tractable representation of this matrix.
`Taking into account emission, we have derived a radiance for-
`mulation for three point transport:
`
`Ljk = Ejk + Z LijRkj£
`
`(2)
`
`This formulation states that:
`
`The radiance atArea j in the direction ofA rea k is equal
`to the radiance emitted by j in the direction of is, plus,
`for every Area i, the radiance at i in the direction of j
`multiplied by the area refiectance Rt,-.‘.
`
`Note that equation (2) is very similar to the radiosity formulation:
`
`Bi = E: +/7:‘ Z325‘-'
`
`All
`
`[ L(xI, 1/_//)G(zI, z//)d$I/ =
`[
`f,(:5, :5’, :5")L(a:, :5')G(:5, :5')G(:5', a:")d:5"d:5
`
`A A”
`
`where
`
`I
`’ _
`cos 0, cos 9,-
`G(z7$ ) — '1: __$,'2 v(z,z')
`where v(:5,:5') is 1 if points :5,
`:5’, are mutually visible, and 0
`otherwise. Note that G is very similar to a differential form factor.
`We integrate over A’, thus introducing all three areas into the
`formulation:
`
`A] All
`
`(1)
`[ [ L(:5',a:’')G(:5',:5'’)d:5"d:5'=
`ff
`f,(:5, :5’, :5")L(:1:,:I;')G(:5,a:')G(:5', z")d:5"da:'da:
`
`A AI All
`
`We may now rewrite the equation in discrete form. Let A,- and At
`be subareas of A’ and A" such that L(:5’ , :5”) is nearly constant over
`their surfaces. The left side of equation (1) may then be rewritten,
`bringing radiance out of the integral as L,-;,:
`
`
`
`Ljk/ [ G'(:l,‘/,:I:")d:I."’d:I:, = 7rLJ'kA,'FJ~;¢
`
`A,‘
`
`AL.
`
`156
`
`
`
`

`
`
`
`COMPUTER GRAPHICS Proceedings, Annual Conference Series, 1993
`
`_-scam.»
`
`2.2 Area Reflectance
`
`The quantity Rt,-, has a natural and satisfying physical significance
`_— it is an expression of reflectance over areas A.-, A,~, and At.
`Consider the fraction of the radiant flux transported from A,
`incident to A,- that is reflected in the direction of area At:
`
`fm fAJ_ flak f,.(:1:,2:’,:1:")L(:1:,:1:')G(:I:,:z:')G(:1:',:1:")d:1:"d:I:'d:I:
`I
`.7
`IA. fA_ L(:1:,:1:’)G(1:,:1:’)d2:’d:1:
`
`In general, and as is shown for glossy reflection in Section 4,
`the accuracy of estimators for F,-.~ and Sr,-; is dependent on the
`size of the patches over which reflectance is computed, relative
`to their distance apart. As relative size decreases, so does error in
`computation, leading directly to the adaptive refinement strategy for
`illumination presented in Section 3 below.
`
`3 Algorithms for Three Point Transport
`
`If we assume that incident radiance is uniform and isotropic over
`both wfi (as induced by Ar) and A,-, we may divide through by
`L(a:, as’ ), yielding:
`
`3.1
`
`Introduction
`
`Recall equation (2):
`
`p(A,‘, Aj,Ak) E
`fit; fAJ_ flak f,(:c,:v’,:t:")G(:c,:z:')G(:c’,a:")d:c"d:c’d:z:
`fA‘_ ff“ G(:c,:c’)da:’d:z:
`
`We define p(A.-,A,-,A;,) to be area reflectance. Note that area
`reflectance is similar to biconical reflectance [11], save that it is also
`integrated over the reflecting surface.
`By definition of Rijki
`
`R-we = Pi/L, AjaAk)
`
`Conservation of energy over reflection, and the reciprocity rela-
`tion derived for R,-,-;, above, constitute fundamental properties of
`area reflectance:
`
`1. 2:12,,-,, g 1,for fixed i, j.
`k
`
`2. A.'F.'jR.'jk = AkFkjRkji-
`
`where equality is achieved in property 1 over complete enclosures
`and perfect reflectivity.
`
`2.3 Evaluation of R“,
`
`In this section we examine the evaluation of Rkji over given patches
`A,, A,- , At.
`Recall:
`
`Rkji =
`
`J
`I
`flak fit’ fA_ f,(:z:",:z:’,:c)G(:z:",:1:')G(a:',:z:)d:cd:z:'d:z:"
`J
`IA,‘ IA‘ G(:1:”,:z:’)d:1:’d:z:”
`
`We assume that discrete areas A,-, Aj, A;. are of small enough
`scale that f,. and G are relatively constant over their surfaces. Then:
`
`Rkji
`
`Skj-.'GkjGjiAlcAjAi
`G)¢_,'AkA_,'
`= SkjiGjiAi
`
`where S is the discretized value of f,, S,-,-;. = S;,,-.- = S,,c$J.,_..
`is
`Note that the average value of G(:c',a:) over A, and A,-
`7rF,~,-/A,- — we thus estimate G,-;A,- by 1rF,-,-, and compute R;.,-.-
`as:
`
`Rkji = 7TFj.'Skjt
`
`In practice, it will not be possible to compute the exact values
`of Fj; and S;,,-,- over A,-, A,-, A;.. We assume that we are able to
`estimate these values, along with error bounds for each estimation.
`Let AF}; and AS;,,-.- be error estimates for computed F‘,-,~ and SW,
`pespectively. We then have an estimate for area reflectance in the
`orm:
`
`Rkji
`
`1r(F,-.- + AF}'.')(Skj.' + A-Skji)
`7r(Fj;Skj; + AFj,‘Sk_,'-,‘ + ASkjiFji + AF_,‘,;ASk_,',~)
`7F(Fj{Skji + AFj.'Skj.' + ASkjiFji)
`
`=2:
`
`Assuming AFj.' < Fji. AS;,,-.- < Skji, we have neglected the last
`term and estimate the error in Rug as 7r(AF};S;,,-,- + AS;,J~,-Fj,-).
`
`.
`
`.,
`
`,,,,.,_-L:-_y-.—»w...~.~.n-s~a.:)..:r,.....».w..«.m...,
`
`Li), = E,-1; + E L.'_,'Rk_,-,~
`
`This equation suggests both a solution strategy for radiance under
`three point transport, and a natural representation for illumination
`within the solution system.
`We may interpret equation (2) as a gathering iteration similar to
`that employed for radiosity under diffuse reflection: the radiance
`L,-k at patch A,- in the direction of patch A;, is found by gathering
`radiances L,-,- in the direction of A,» at patches A,-. We may solve
`for transport by gathering radiance for each L,-;., and successively
`iterating to capture all significant re-reflection.
`We are left with the question of what structure we are gathering
`over and iterating upon. Note that all illumination is expressed
`as the radiance at a given patch in the direction of another — it
`is these patch-patch interactions that form the primary structure
`within the solution system. All operation is over interactions: both
`the representation and transport of radiance, and the iteration and
`solution for illumination.
`Consider the following structure:
`
`typedef struct _interaction (
`Patch * from.-
`Patch ‘to;
`Color
`L;
`Color
`Lg,-
`List:
`*gat:her;
`struct _irrt:eraction *nw,
`) Interaction;
`
`‘sw, *se,
`
`‘me;
`
`A given interaction ij is defined by two patches i j —>from and
`ij —>to, and represents the radiance at from in the direction of
`to. This radiance is stored within the interaction as attribute L.
`Lg is radiance gathered during the current solution iteration from
`interactions contained in the list gather. Subinteractions nw, sw,
`se, ne are the children of z’j , induced by subdivision over either
`from or to. The structure assumes quadtree refinement, leaving
`northwest, southwest, southeast, and northeast descendants.
`In the following sections we will present an algorithm for the
`refinement and computation of illumination over a hierarchy of
`interactions. The algorithm will operate by refining pairs of inter-
`actions ij, jk (such that ij —>to == jk—>from), to ensure that
`computed reflectance across the interaction pairs, and associated
`patch triples, satisfies user specified error bounds. If a given in-
`teraction pair ij, jk is satisfactory, the interactions are linked to
`record that radiance may be gathered from ij to jk, otherwise one
`or both interactions are subdivided and refinement applied to their
`descendants.
`
`After refinement, a gathering iteration may be carried out, each
`interaction gathering radiance from interactions to which it has been
`linked. The gathered radiances are then distributed within each re-
`ceiving interaction hierarchy, and subsequent iterations computed
`until satisfactory convergence has been achieved.
`Note that, within this system, the eye may be regarded as simply
`another object with which patches may interact. The radiance along
`interactions to the eye provides the resulting view.
`
`157
`
`
`
`

`
`
`
`
`
`Figure 3: Refinement and Subdivision
`
`3.3 Gathering Radiance
`
`Gathering radiance over interactions may be written as a simple
`procedure:
`
`%ather(Interaction *jk)
`Interaction *ij;
`if (jk)
`(
`jk—>Lg = O;
`jk—>gather)
`ForAllElements(ij,
`jk—>Lg += ij—>L * Reflectance(ij,
`Gather(jk—>nw);
`Gather(jk—>sw);
`Gather(jk—>se);
`
`jk);
`
`) Gather(jk—>ne);
`
`l W
`
`e gather radiance into jk—->Lg rather than directly into jk—>L
`to avoid the necessity of a push/pull with every invocation of the
`procedure (see Section 3.4). The solution method is thus simple
`Jacobi iteration, as opposed to Gauss-Seidel, as the hierarchical
`structure imposes simultaneous rather than successive displacement.
`
`3.4 Radiance within a Hierarchy
`A gathering iteration results in received radiance scattered through-
`out each interaction hierarchy. This gathered radiance must be dis-
`tributed and accounted for over all ancestors and descendants of
`each receiving interaction, in order to maintain the consistency and
`correctness of the hierarchical representation of radiance between
`patches.
`We employ a distribution algorithm similar to that presented in
`[7] for radiosity over patch/element hierarchies: gathered radiance
`is “pushed” to the leaf interactions within each hierarchy to ensure
`propagation to all descendants, and then “pulled” and distributed
`back up from the leaves through all higher level interactions to
`their common ancestor at the root. As is shown in [2], radiance
`may be pushed unchanged within the interaction hierarchy, and area
`averaged as it is pulled from child to parent.
`
`4 Application over Glossy Reflection
`
`In this section we discuss our implementation of the above algo-
`rithms over glossy reflection.
`
`4.1 The Reflection Function
`
`We employ a highly simplified Torrance-Sparrow [15] model for
`our glossy reflection function:
`
`f(d)._d),)_n+2 cos"l9,,,
`9
`" T _ 87r cos0.-cos0,
`This function incorporates the facet distribution function cos‘ Om
`developed by Blinn [3], normalized for projected facet area under
`
`sh(9;,0r)
`
`
`
`SIGGRAPH 93, Anaheim, California, 1-6 August 1993
`
`3.2 Adaptive Refinement
`
`Consider the following procedure:
`ReEine(Interaction *ij, Interaction *jk,
`(
`float Feps, float Seps, float Aeps)
`float feps, seps;
`feps
`GeometryErrorEstimate(ij);
`seps
`ReflectionErrorEstimate(ij,
`if (feps < Feps && seps < Seps)
`Link(ij:
`jk);
`else if (seps >= seps)
`switch(SubdivS(ij,
`case PATCH_I:
`
`(
`jk, Aeps))
`
`C
`
`jk);
`
`1 1 E i i l
`
`jk, seps, Feps, Aeps);
`jk, seps, Feps, Aeps);
`jk, Seps, Feps, Aeps);
`jk, seps, Feps, Aeps);
`
`Refine(ij—>nw,
`Refine(ij—>sw,
`Refine(ij—>se,
`Refine(ij—>ne,
`break;
`case PATCH_J:
`/* refine over children of ij and jk */
`case PATCH_K:
`/* refine over children of jk
`case NONE:
`
`*/
`
`)
`
`Link(ij,
`
`jk);
`
`lse (
`
`/* feps >= Feps */
`(
`switch(SubdivG(ij, jk, Aeps))
`/* refine over children, or link, as
`/* directed by PATCH_I, J, K, or NONE.
`l
`
`*/
`*/
`
`l e
`
`)
`
`} T
`
`his procedure computes over pairs of interactions, and associated
`patch triples, subdividing and recursively refining if estimated error
`exceeds user specified bounds, linking the interactions for gathering
`if the bounds are satisfied, or if no further subdivision is possible.
`Feps and seps are the bounds for geometric and reflection error,
`respectively; Aeps specifies the minimum area a patch may possess
`and still be subdivided. GeometryErrorEstimate and Re-
`f lec tionErrorEs timate provide estimations for 7rAF,-,~.S'k_,-.-
`and‘nZSSkj;P}a
`Subdivs and SubdivG control refinement for reflection and
`geometry error, respectively. Both routines select a patch for refine-
`ment, subdividing the patch and associated interaction(s) if required.
`An identifier for the selected patch is returned — if no patch may
`be subdivided, then NONE is passed back. Note that a given in-
`teraction/patch may be refined against many different interactions
`within the system, and thus may have already been subdivided when
`selected by a Subdiv routine — in this case, the routine simply
`returns the proper identifier.
`The Subdiv routines should select for refinement patches that
`are of large size relative to their distance from their partner(s) in the
`transport triple. Form fact