- Email: [email protected]
Unified sun and sky illumination for shadows under trees
CVGIP:
GRAPHICAL
MODELS
AND
IMAGE
PROCESSING
Vol. 53, No. 3, May, pp. 223-230, 1991
Unified Sun and Sky Illumination NELSON Lawrence
Livermore
National
for Shadows under Trees’ MAX
Laboratory,
Livermore,
California
94550
Received July 27, 1987; accepted July 18, 1990
Assumethat the sun,sky, and cloudsare representedasa raster imageon a planeat infinity, and the shadowsof a tree canopy are representedas a raster transparency mask on a canopy plane. Then the illumination on the ground isthe correlationof thesetwo rasters,which can be rapidly computedby fast Fourier transform techniques. 0 1991 Academic Press, Inc.
INTRODUCTION
This paper presents an analytic method for computing sun and sky illumination under a canopy of trees. Recent illumination models have become increasingly accurate in accounting for light emanating from the whole environment, and have also become increasingly expensive to compute. Distributed ray tracing [l] requires enough rays at each pixel to adequately sample the illuminating environment. For a point light source, a single “shadow feeler” will suffice, but for area light sources, multiple rays must sample the source area, and for diffuse reflection of the environment, the rays must sample the whole hemisphere above the surface position [2]. Certain simplifications result if the environment is “at infinity,” so that it is independent of the viewer’s position [3-51. A clear or cloudy sky is, for practical purposes, at infinity. Other polygon-based algorithms generate the environment analytically, with a visible-surface computation from the point of view of each surface position. Finding
the visible surfaces of only the light sources is sufficient to render shadow penumbras from area light sources [8]. The “radiosity” method can handle diffuse [6] and specular [7] interreflection between surfaces. In order to be computationally practical, these analytic algorithms compute the environment at only a sampled collection of surface points, and render the image by interpolation. Since the fast Fourier transform techniques used here compute all the shadows simultaneously, interpolation is not required. Nishita and Nakamae [9] compute sky illumination using semicircular scan lines across a sky hemisphere. But it is impractical to make these scan lines fine enough to sample the sun, so a separate shadow computation is required to include direct sunlight. Takita et al. [IO] cast sun shadow penumbras from trees defined by several color and transparency rasters on intersecting planes. For every pixel, they intersect these rasters with a cone toward the sun. The method presented here treats the sky as a raster image, so it can include skylight, cloud glow, and sunlight in a uniform manner. The shadow casting objects are also represented in a raster image, allowing finely detailed leaf shadows. However, the method assumes all shadow casting objects lie near a canopy plane parallel to the ground, and only works for shadows cast onto the ground. Shadows cast onto other surfaces must be computed separately by some other method. SHADOWS AS CORRELATIONS
I This document was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor the University of California nor any of their employees make any warranty, expressed or implied, or assume any legal liability or responsibility for the accuracy, completeness, or usefulness or any information, apparatus, product, or process disclosed, or represent that its use would not infringe on privately owned rights. Reference herein to any specific commercial products, processes, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or the University of California. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or the University of California, and shall not be used for advertising or other product endorsement purposes.
Consider the shadows cast onto the ground plane G in Fig. 1, by a collection of polygonally defined leaves and branches lying near a canopy plane C parallel to the ground, from a point light source at infinity in the direction L. The shadows can be considered as a polygonally defined surface texture on the plane G. If this texture is translated in the direction L onto the plane C, it defines a collection Q of polygons in the plane C which would generate the same shadows as the actual leaf and branch polygons. We compute the sun, sky, and cloud illumination, as obscured by the polygons Q in C. Since rays from the sun are all nearly parallel to L, and the leaves all lie 223 1049-%52/91 $3.00 Copyright 0 1991by Academic Press, Inc. All rights of reproduction in any form reserved.
NELSON
MAX
and
Fig. 2. The unit hemisphere H of directions the plane S tangent to it at the north pole
above a surface point 0, N. An element do of solid
angle on H, with sides d0 and sin Od4, projects to an element of area dA on S with sides se&d8 and tan Od+, which in turn projects on the ground plane G to an element of area dB with sides cos OdO and sin Qd4. Fig. source
1. The canopy is in the direction
plane C parallel L. Polygons
and Sz on G, which translate in the shadow of P, from varies between
as the product L’ and L.
to the ground P,, and I’? near
plane G. The light C cast shadows S,
to Q, and Q? in C. The positioning error a direction L’, caused by using Q, instead,
of the
distance
d of P, from
C and
the
angle
e
IO = j,h( V, L)E(L)cos 8do
0
near C, the errors in penumbras cast by the leaves onto the ground will be of second order. (See the caption to figure 1.) The illumination from directions farther away from L is much less bright, so the errors in assuming all the leaves lie in C will not be very noticeable. The polygons in the set Q on the canopy plane C can be scan converted into an anti-ahased raster transparency mask, m(i, j), which is 1 for unobscured pixels, 0 for totally shadowed pixels, and some value between 0 and 1 for partially obscured pixels. Since it is this mask raster which will generate the shadows, polygons are not actually required. Any method which can generate a raster shadow texture on the ground plane (for example, a visible surface algorithm from the point of view of the sun) can be used to prepare the transparency mask. Let E(L) be the sky color in the direction L. Instead of projecting the sky environment onto a hemisphere at infinity, we project it onto a sky plane at infinity, parallel to the ground. The sky will appear the same from any point on the ground; its color will depend only on the direction of a ray, not the position. Thus we can think of the sky as an infinite plane S tangent to the unit hemisphere H of upward directions. As shown on Fig. 2, an element of solid angle dw = &(sin I%@) on the hemisphere H corresponds to an area dA = (sec2 O)(tan &f+) on the plane S, so dw = cos3fldA. Then, when there are no shadows, the unobscured intensity of the ground is (see 151)
= s &( V, L)E(L)cos4 BdA, i
(1)
where fr( V, L) is the bidirectional reflection function of the viewing direction V and the illumination direction L. In the images presented here, I used Lambert’s law, wheref,(V, L) is constant. But the method works whenever f,(V, L) is independent of V, so it is appropriate whenever the pyramid of view containing the possible viewing vectors V is small, and/or fr varies slowly. Because of the cos4 8 weighting factor, it is possible to get a good approximation to I0 by integrating only over a
Fig. 3. The coordinates on the canopy plane spond to those on the plane S of directions above ground G. The coordinates on the ground plane lated from those on C. When the hemisphere H tions are intersects
moved above C in the point
the pixel (k, I) on (i + k, j + I).
G, the
C are scaled to correthe origin (0, 0) on the G are vertically transand plane .S of direcray
in direction
(i, j)
225
UNIFIEDSUNANDSKYILLUMINATION
square region R of the infinite plane S, as long as R contains the sun. On R, we representf,(V, L)E(L)cos4 0 as an illumination raster e(i, j), defined for 0 I i, j 5 n - 1. Then Z0 = 2::; xyzd e(i, j), up to a normalizing factor which we neglect here for simplicity. Now we wish to include the shadows, and to compute an intensity raster g(i, j) at each point on the ground plane G. We organize the coordinate systems on the mask raster m(i, j), the illumination raster e(i, j), and the ground intensity raster g(i, j), as shown in Fig. 3. We assume the leaves are dark and opaque, and contribute no illumination of their own; they merely obscure light from the sky. A ray from (0, 0) on the ground, in the direction (i, 1) on the sky plane, will intersect the canopy plane at index (i,j) also, so g(0, 0) = ~~~~ xycd e(i,j) m(i j). Now extend m(i, j) and g(i, j) to be periodic functions in i andj, of period n. This is equivalent to assuming that the canopy of trees is periodic, which is reasonable for a section of a tightly packed forest. For an isolated clump of trees, we scale the rasters so there are gaps between the periodic clumps, and then artificially remove the extraneous shadows, by assuming points on the ground outside the first period of g(i,j) are completely unshadowed. We interpolate to I0 in a region near the edges of the first period to eliminate any shading discontinuity. Because of the cos4 8 factor in Eq. (1), the errors introduced by the periodicity rapidly become small if the periods are made large enough. As shown in Fig. 3, a ray from (k, 1) on the ground, in the direction (i,J, will intersect the canopy at (i + k, j + 0 so g(k, 1) = 2 c e(i,j)m(i i-0
+ k, j + 0.
one dimension while vectorizing over the other. In addition, since e and m are real, the FFT along the first index gives only n/2 + 1 independent complex numbers for each column of pixels, instead of n, resulting in a time saving factor of almost 2. The shadows of the leaves and branches on each other were computed by the method of Max [12] along radial scan lines. Each scan line is the intersection with the picture plane of a scan plane through the eye and light source. After all the shadow segments were accumulated in each scan plane, they were written out into a file. The scan planes also intersect the ground plane in a series of radial scan lines. Therefore, the shadows in the file can be scan converted into the mask raster m(i, j) shown in Fig. 4, using the anti-aliasing algorithm in Max 1131. The image of the cloudy sky shown in Fig. 5 was generated by the method of Gardner [14]. The sun was added, with the sun glow through the clouds calculated by the method of Blinn [15]. When multiplied by cos4 8, this resulted in the raster e(i,j) shown in Fig. 6. A raster g(i, j) for the intensity reflected from the environmental illumination by the ground, as shown in Fig. 7, was computed by the FFT correlation technique discussed above. Then the ground and sky images were resampled in perspective as a background raster in the picture plane, using standard area averaging techniques [16, 171. The visible segments of the trees were also scan converted as a foreground raster, shown in Fig. 8, using the methods of [13]. An (Ychannel, shown in Fig. 9, carrying
(2)
j=()
This means that g is the correlation of e and m (see [I 11). The summation in Eq. (2) requires n2 multiplications for reach pixel (k, Z), and there are n2 pixels, so the total computation is of order n 4. With fast Fourier transform (FFT) techniques, it is possible to compute the correlation at all pixels in time of order n* log n, as follows. Let G, E, and M be the Fourier transforms of g, e, and m, respectively. Then by the correlation theorem (see [ 1l]), G = Ea, where @is the complex conjugate of M. Therefore, to get g, we use the FFT algorithm to find E and M, multiply E by the complex conjugate of M at each of the n2 points in frequency space, and then use the FFT algorithm to find the inverse Fourier transform g of Ea. IMPLEMENTATION
AND
RESULTS
On the Cray-1 or the Cray XMP, I used an FFT algorithm for two-dimensional arrays, which transforms in
Fig. 4.
Transparency
mask
m(i. j) from
leaf shadows.
226
NELSON
Fig. 5.
Image
of cloudy
MAX
sky. Fig. 7.
Ground
intensity
raster
~(i, j).
pixel coverage information [IS], was produced at the same time. Then the background and foreground were combined into the final image, shown in Fig. 10: jinal(i, j) = forPground(i, j) + (1 - a(i,j))background(i,
j).
If atmospheric illumination is desired (see [ 12]), it can be scan converted as an extra channel, and added on top of the final image.
Fig. 6. Sky raster e(i, j). The small sun disk in the center one thousand times as bright as the brightest clouds.
is actually Fig. 9.
The (Y channel,
used for compositing.
UNIFIED
SUN AND SKY ILLUMINATION
227
Fig. 8. Trees, with shadows computed by the method of Max [121, as if the sun were a point source.
Fig. 10. The final image, made by resampling Figs. 5 and 7 in perspective, and compositing Fig. 8 on top.
Fig. 11. Similar to Fig. 10 but without sky illumination bras, using Fig. 4 instead of Fig. 7.
and penum-
Fig. 13. A composite image in a partial eclipse, using Fig. 12 as the ground raster. Due to the limitations in the color printing process, the subtle differences between Fig. 13 and Fig. 10 are not visible here.
Note the effects of the sky illumination in Figs. 7 and 10. The ground is darker in the center of the shadow of a large clump of trees than it is in shadows of branches or smaller clumps. Also, the ground is brighter far from the shadow of the trees than it is nearer to the shadow, but
still outside its penumbra from the sun. Compare this with Fig. 11, without these effects and without penmbras. Figures 12 and 13 show the effects of a partial eclipse on the penumbra of the sun. I have produced a computer animated film “Sun and
NELSON MAX
228
Fig. 14. The square R in the sky plane S projects radially to the spherical polygon, F = ABCD on the hemisphere H, which in turn projects vertically to the area K = A’B’C’D’.
Fig. 12. The ground raster for an eclipsed sun. Note crescent shaped spots of light.
Shade” [19], showing the illumination change in a grove of trees as the moon or clouds cover the sun. I used a 48 frame cycle of images of the trees, without ground or sky, showing the leaves fluttering in the breeze. I rendered the foreground images twice, once including point source illumination at the sun, and once with only diffuse illumination from the sky. I also computed the shadow mask raster m(i, j) for each frame of the cycle, and its Fourier transform. Then for each frame of the longer final sequence, the background image is found as above, using the precomputed Fourier transform M of the mask m. The fraction S of the direct sunlight which passes through the clouds is used to interpolate between the two foreground images, which are then combined with the background using the precomputed (Y channel. Finally, the atmospheric illumination is multiplied by S, and added to the image. The total cost on one processor of a Cray XMP was on the order of 1 min per 640 x 480 pixel frame, and was dominated by the anti-aliased resampling algorithm in [13]. This cost is amortized over the cyclic repeats, while the image interpolation and compositing costs are small. ERROR ANALYSIS
The integral of Eq. (1) is over an infinite plane S, while the sum of Eq. (2) approximates the integral only over a
finite square R in this plane. The size of this square is determined by the angle (Y between the vertical line ON, and the line OM from 0 to the center of a side of R, as shown in Fig. 14. To estimate the error introduced by neglecting the region of S outside R, we assume Lambet-t’s law, so thatf,(V, L) = 1, and that the sky illumination is also 1, independent of L. The integral (1) then takes the form Z,,(a) = j)os4 =
8dA
I F cos 8do
=I dB. K
Here F is the radial projection of R onto the unit hemisphere H, as shown in Fig. 14, and K is the vertical projection of F onto the ground plane G. The quadrilateral F is bounded by great circle arcs through the points A, B, C, and D, so K is bounded by ellipses through their projections A’, B’, C’, and D’. The last two integrals are equal because the area element dB on the ground plane is equal to cos 8dw. (See the caption to Fig. 2.) Now the last integral is simply the area of the region K. Figure 15 shows K divided into eight equal elliptical sectors like OPC’. The ellipses are tangent to the unit circle, with minor axis OP = sin (Y. Let U be the point on the unit circle above C’, and let 4 be the angle POW. The sector OPC’ is formed by vertically shrinking the circular sector OTU, of area 44, by a factor sin (Y, so this area is t+ sin (Y, and also QO = QC’ = QU sin (Y. Then, tan 4 = VU/V0 = QOIQU = QU sin alQU = sin cy, so $I = tan-‘(sin a). Thus the area of G is lo(a) = 8(&tan-‘(sin a)sin a) = 4sin (Y tan-‘(sin a). When (Y = 90”, Z,,(a) =
UNIFIED
229
SUN AND SKY ILLUMINATION
EXTENSIONS TO MULTIPLE SHADOW CASTING PLANES
Fig. 15. A top view of the area A’B’C’D’, which is bounded by four ellipses, and is divided by horizontal, vertical, and 45” diagonal lines into eight equal elliptical sectors like OPC’.
4tan-‘(I) = 7~, the area of the unit circle. Figure 16 shows a graph of the fraction Zo(a)17r of the unit circle area covered by K. A related curve appears as Fig. 3 in Reeker et al. [20]. In the examples in this paper, I used CY= 70”, and included 90% of a uniform sky illumination. For pixels below the center of R, this graph also measures the error from extending the mask m(i, j) periodically outside R, but for pixels farther out, the error from assuming periodicity is larger. When the sun is not obscured by clouds, and lies inside R, its energy outweighs that from the sky, and the errors discussed above become a smaller proportion of the total illumination on the ground.
Image correlation by fast Fourier transforms is an efficient technique for computing shadows on the ground from environmental illumination, provided the shadow casting objects lie in or near a single horizontal plane. Mask rasters ml and m2, in two planes at different heights, combine with the sky raster e in a more complex formula, for which there is no obvious FFT speedup. However, if we assume that the two masks are uncorrelated, we can approximate the composition of their two shadowing effects by computing the two ground intensity rasters g as above, multiplying them together, and dividing by the unobstructed intensity Z,. Takita et al. [IO] represent their trees by several horizontal raster images with color and transparency, and one such vertical image, which is always turned to face the eye. The cone toward the sun intersects the plane of one these images in an ellipse, and the fraction of pixels inside this ellipse which are transparent gives the penumbra illumination. If there are several planes, Takita et al. take the combined penumbra illumination to be the minimum of the results for the individual planes, while the uncorrelated assumption above would take the product. If the several rasters were instead sampled by a fixed or jittered grid of rays in the cone, the transparencies for each ray could be multiplied separately, giving a better estimate for the correlation. The cone computation can produce the variably sized penumbra along the trunk of the tree, as well as the correctly sized penumbras for each level of leaves. However, it requires extensive computations at each pixel. The algorithm here, applied to several leaf levels assumed to be uncorrelated, produces less accurate results much more rapidly. ACKNOWLEDGMENTS This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract Number W-7405ENG-48. The idea for this algorithm arose from discussions about shadows with Chuck Grant. I thank Brian Cabral, Becky Springmeyer, Chuck Grant, Michael Gwilliam, and Jeff Kallman for helpful suggestions on the exposition, Jules Bloomenthal for the tree trunk data and bark texture, Becky Springmeyer and Chuck Grant for the typing and layout, Ted Einwohner for Fig. 16, and Brian Cabral for providing the cloud texture in Fig. 5. The section on error analysis was recommended by an anonymous referee.
REFERENCES ,
20
40
60
80
Fig. 16. A graph of the function J,,((Y)/~, showing the fraction of a uniform sky’s energy included in a square subtending an angle of 2a.
I. R. L. Cook, T. Porter, and L. Carpenter, Distributed ray tracing, Compuf. Graphics 18, No. 3, 1984, 137-145. 2. J. T. Kajiya, The rendering equation, Comput Graphics 20, No. 4, 1986, 143-150.
230
NELSON
3. N. Greene, Applications of world projections, IEEE Cotnput. Graphics Appl. 6, No. 11, 1986, 21-29. 4. G. S. Miller and R. C. Hoffman, Illumination and Reflection Maps: Simulated Objects in Simulated and Real Environments, SIGGRAPH ‘84 Advanced Computer Graphics Animation Course Notes, July 1984. 5. B. Cabral, N. Max, and R. Springmeyer, Bidirectional reflection functions from surface bump maps, Comput. Graphics 24, No. 4, 1987, 273-281. 6. M. F. Cohen and D. P. Greenberg, A radiosity solution for complex environments, Comput. Graphics 19, No. 3, 1985, 31-40. 7. D. S. Immel, M. F. Cohen, and D. P. Greenberg, A radiosity method for non-diffuse environments, Comput. Graphics 20, No. 4, 1986, 133-142. 8. T. Nishita and E. Nakamae, Half tone representation of 3-D objects illuminated by area sources or polyhedron sources, IEEE Proc. COMSAC, 1983, 237-241. 9. T. Nishita and E. Nakamae, Continuous tone representation of three dimensional objects illuminated by sky light, Comput. Graphics 20, No. 4, 1986, 125-132. 10. S. Takita, K. Kaneda, T. Akinoba, H. Iriyama, E. Nakamae, and T. Nishita, A simple rendering for penumbra caused by sunshine, in Proceedings, CG International 90; New Advances in Computer Graphics, Springer-Verlag, New York, 1990.
MAX 11. R. C. Gonzalez and P. Wintz, Digital Image Processing, AddisonWesley, Reading, MA, 1977. 12. N. Max, Atmospheric illumination and shadows, Comput. Graphics 20, No. 4, 1986, 117-124. 13. N. Max, Anti-aliasing scan line data, IEEE Comput. Graphics Appl. 10, No. 1, 1990, 18-30. 14. G. Y. Gardner, Visual simulation of clouds, Comput. Graphics 19, No. 3, 1985, 297-303. 15. J. F. Blinn, Light reflection functions for simulation of clouds and dusty surfaces, Comput. Graphics 16, No. 3, 1982, 21-29. 16. E. Catmull and A. R. Smith, 3-D transformations in scanline order, Comput. Graphics, 14, No. 3, 1980, 279-285. 17. K. M. Fant, A nonaliasing, real time spatial transform technique, IEEE Computer Graphics Appl. 6, No. 1, 1986, 71-80 (see also Letters to the Editor in 6, Nos. 3 and 7). 18. T. Porter and T. Duff, Compositing digital images, Comput. Graphics 18, No. 3, 1984, 253-259. 19. N. Max, Sun and Shade, SIGGRAPH Video Review, Issue No. 36, Segment 8, 1988, ACM, available from First Priority, P.O. Box 576, Itasca, IL 60143. 20. R. J. Reeker, D. W. George, and D. P. Greenberg, Acceleration techniques for progressive refinement radiosity, Comput. Graphics 24, No. 2, 1990, 59-66.
GRAPHICAL
MODELS
AND
IMAGE
PROCESSING
Vol. 53, No. 3, May, pp. 223-230, 1991
Unified Sun and Sky Illumination NELSON Lawrence
Livermore
National
for Shadows under Trees’ MAX
Laboratory,
Livermore,
California
94550
Received July 27, 1987; accepted July 18, 1990
Assumethat the sun,sky, and cloudsare representedasa raster imageon a planeat infinity, and the shadowsof a tree canopy are representedas a raster transparency mask on a canopy plane. Then the illumination on the ground isthe correlationof thesetwo rasters,which can be rapidly computedby fast Fourier transform techniques. 0 1991 Academic Press, Inc.
INTRODUCTION
This paper presents an analytic method for computing sun and sky illumination under a canopy of trees. Recent illumination models have become increasingly accurate in accounting for light emanating from the whole environment, and have also become increasingly expensive to compute. Distributed ray tracing [l] requires enough rays at each pixel to adequately sample the illuminating environment. For a point light source, a single “shadow feeler” will suffice, but for area light sources, multiple rays must sample the source area, and for diffuse reflection of the environment, the rays must sample the whole hemisphere above the surface position [2]. Certain simplifications result if the environment is “at infinity,” so that it is independent of the viewer’s position [3-51. A clear or cloudy sky is, for practical purposes, at infinity. Other polygon-based algorithms generate the environment analytically, with a visible-surface computation from the point of view of each surface position. Finding
the visible surfaces of only the light sources is sufficient to render shadow penumbras from area light sources [8]. The “radiosity” method can handle diffuse [6] and specular [7] interreflection between surfaces. In order to be computationally practical, these analytic algorithms compute the environment at only a sampled collection of surface points, and render the image by interpolation. Since the fast Fourier transform techniques used here compute all the shadows simultaneously, interpolation is not required. Nishita and Nakamae [9] compute sky illumination using semicircular scan lines across a sky hemisphere. But it is impractical to make these scan lines fine enough to sample the sun, so a separate shadow computation is required to include direct sunlight. Takita et al. [IO] cast sun shadow penumbras from trees defined by several color and transparency rasters on intersecting planes. For every pixel, they intersect these rasters with a cone toward the sun. The method presented here treats the sky as a raster image, so it can include skylight, cloud glow, and sunlight in a uniform manner. The shadow casting objects are also represented in a raster image, allowing finely detailed leaf shadows. However, the method assumes all shadow casting objects lie near a canopy plane parallel to the ground, and only works for shadows cast onto the ground. Shadows cast onto other surfaces must be computed separately by some other method. SHADOWS AS CORRELATIONS
I This document was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor the University of California nor any of their employees make any warranty, expressed or implied, or assume any legal liability or responsibility for the accuracy, completeness, or usefulness or any information, apparatus, product, or process disclosed, or represent that its use would not infringe on privately owned rights. Reference herein to any specific commercial products, processes, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or the University of California. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or the University of California, and shall not be used for advertising or other product endorsement purposes.
Consider the shadows cast onto the ground plane G in Fig. 1, by a collection of polygonally defined leaves and branches lying near a canopy plane C parallel to the ground, from a point light source at infinity in the direction L. The shadows can be considered as a polygonally defined surface texture on the plane G. If this texture is translated in the direction L onto the plane C, it defines a collection Q of polygons in the plane C which would generate the same shadows as the actual leaf and branch polygons. We compute the sun, sky, and cloud illumination, as obscured by the polygons Q in C. Since rays from the sun are all nearly parallel to L, and the leaves all lie 223 1049-%52/91 $3.00 Copyright 0 1991by Academic Press, Inc. All rights of reproduction in any form reserved.
NELSON
MAX
and
Fig. 2. The unit hemisphere H of directions the plane S tangent to it at the north pole
above a surface point 0, N. An element do of solid
angle on H, with sides d0 and sin Od4, projects to an element of area dA on S with sides se&d8 and tan Od+, which in turn projects on the ground plane G to an element of area dB with sides cos OdO and sin Qd4. Fig. source
1. The canopy is in the direction
plane C parallel L. Polygons
and Sz on G, which translate in the shadow of P, from varies between
as the product L’ and L.
to the ground P,, and I’? near
plane G. The light C cast shadows S,
to Q, and Q? in C. The positioning error a direction L’, caused by using Q, instead,
of the
distance
d of P, from
C and
the
angle
e
IO = j,h( V, L)E(L)cos 8do
0
near C, the errors in penumbras cast by the leaves onto the ground will be of second order. (See the caption to figure 1.) The illumination from directions farther away from L is much less bright, so the errors in assuming all the leaves lie in C will not be very noticeable. The polygons in the set Q on the canopy plane C can be scan converted into an anti-ahased raster transparency mask, m(i, j), which is 1 for unobscured pixels, 0 for totally shadowed pixels, and some value between 0 and 1 for partially obscured pixels. Since it is this mask raster which will generate the shadows, polygons are not actually required. Any method which can generate a raster shadow texture on the ground plane (for example, a visible surface algorithm from the point of view of the sun) can be used to prepare the transparency mask. Let E(L) be the sky color in the direction L. Instead of projecting the sky environment onto a hemisphere at infinity, we project it onto a sky plane at infinity, parallel to the ground. The sky will appear the same from any point on the ground; its color will depend only on the direction of a ray, not the position. Thus we can think of the sky as an infinite plane S tangent to the unit hemisphere H of upward directions. As shown on Fig. 2, an element of solid angle dw = &(sin I%@) on the hemisphere H corresponds to an area dA = (sec2 O)(tan &f+) on the plane S, so dw = cos3fldA. Then, when there are no shadows, the unobscured intensity of the ground is (see 151)
= s &( V, L)E(L)cos4 BdA, i
(1)
where fr( V, L) is the bidirectional reflection function of the viewing direction V and the illumination direction L. In the images presented here, I used Lambert’s law, wheref,(V, L) is constant. But the method works whenever f,(V, L) is independent of V, so it is appropriate whenever the pyramid of view containing the possible viewing vectors V is small, and/or fr varies slowly. Because of the cos4 8 weighting factor, it is possible to get a good approximation to I0 by integrating only over a
Fig. 3. The coordinates on the canopy plane spond to those on the plane S of directions above ground G. The coordinates on the ground plane lated from those on C. When the hemisphere H tions are intersects
moved above C in the point
the pixel (k, I) on (i + k, j + I).
G, the
C are scaled to correthe origin (0, 0) on the G are vertically transand plane .S of direcray
in direction
(i, j)
225
UNIFIEDSUNANDSKYILLUMINATION
square region R of the infinite plane S, as long as R contains the sun. On R, we representf,(V, L)E(L)cos4 0 as an illumination raster e(i, j), defined for 0 I i, j 5 n - 1. Then Z0 = 2::; xyzd e(i, j), up to a normalizing factor which we neglect here for simplicity. Now we wish to include the shadows, and to compute an intensity raster g(i, j) at each point on the ground plane G. We organize the coordinate systems on the mask raster m(i, j), the illumination raster e(i, j), and the ground intensity raster g(i, j), as shown in Fig. 3. We assume the leaves are dark and opaque, and contribute no illumination of their own; they merely obscure light from the sky. A ray from (0, 0) on the ground, in the direction (i, 1) on the sky plane, will intersect the canopy plane at index (i,j) also, so g(0, 0) = ~~~~ xycd e(i,j) m(i j). Now extend m(i, j) and g(i, j) to be periodic functions in i andj, of period n. This is equivalent to assuming that the canopy of trees is periodic, which is reasonable for a section of a tightly packed forest. For an isolated clump of trees, we scale the rasters so there are gaps between the periodic clumps, and then artificially remove the extraneous shadows, by assuming points on the ground outside the first period of g(i,j) are completely unshadowed. We interpolate to I0 in a region near the edges of the first period to eliminate any shading discontinuity. Because of the cos4 8 factor in Eq. (1), the errors introduced by the periodicity rapidly become small if the periods are made large enough. As shown in Fig. 3, a ray from (k, 1) on the ground, in the direction (i,J, will intersect the canopy at (i + k, j + 0 so g(k, 1) = 2 c e(i,j)m(i i-0
+ k, j + 0.
one dimension while vectorizing over the other. In addition, since e and m are real, the FFT along the first index gives only n/2 + 1 independent complex numbers for each column of pixels, instead of n, resulting in a time saving factor of almost 2. The shadows of the leaves and branches on each other were computed by the method of Max [12] along radial scan lines. Each scan line is the intersection with the picture plane of a scan plane through the eye and light source. After all the shadow segments were accumulated in each scan plane, they were written out into a file. The scan planes also intersect the ground plane in a series of radial scan lines. Therefore, the shadows in the file can be scan converted into the mask raster m(i, j) shown in Fig. 4, using the anti-aliasing algorithm in Max 1131. The image of the cloudy sky shown in Fig. 5 was generated by the method of Gardner [14]. The sun was added, with the sun glow through the clouds calculated by the method of Blinn [15]. When multiplied by cos4 8, this resulted in the raster e(i,j) shown in Fig. 6. A raster g(i, j) for the intensity reflected from the environmental illumination by the ground, as shown in Fig. 7, was computed by the FFT correlation technique discussed above. Then the ground and sky images were resampled in perspective as a background raster in the picture plane, using standard area averaging techniques [16, 171. The visible segments of the trees were also scan converted as a foreground raster, shown in Fig. 8, using the methods of [13]. An (Ychannel, shown in Fig. 9, carrying
(2)
j=()
This means that g is the correlation of e and m (see [I 11). The summation in Eq. (2) requires n2 multiplications for reach pixel (k, Z), and there are n2 pixels, so the total computation is of order n 4. With fast Fourier transform (FFT) techniques, it is possible to compute the correlation at all pixels in time of order n* log n, as follows. Let G, E, and M be the Fourier transforms of g, e, and m, respectively. Then by the correlation theorem (see [ 1l]), G = Ea, where @is the complex conjugate of M. Therefore, to get g, we use the FFT algorithm to find E and M, multiply E by the complex conjugate of M at each of the n2 points in frequency space, and then use the FFT algorithm to find the inverse Fourier transform g of Ea. IMPLEMENTATION
AND
RESULTS
On the Cray-1 or the Cray XMP, I used an FFT algorithm for two-dimensional arrays, which transforms in
Fig. 4.
Transparency
mask
m(i. j) from
leaf shadows.
226
NELSON
Fig. 5.
Image
of cloudy
MAX
sky. Fig. 7.
Ground
intensity
raster
~(i, j).
pixel coverage information [IS], was produced at the same time. Then the background and foreground were combined into the final image, shown in Fig. 10: jinal(i, j) = forPground(i, j) + (1 - a(i,j))background(i,
j).
If atmospheric illumination is desired (see [ 12]), it can be scan converted as an extra channel, and added on top of the final image.
Fig. 6. Sky raster e(i, j). The small sun disk in the center one thousand times as bright as the brightest clouds.
is actually Fig. 9.
The (Y channel,
used for compositing.
UNIFIED
SUN AND SKY ILLUMINATION
227
Fig. 8. Trees, with shadows computed by the method of Max [121, as if the sun were a point source.
Fig. 10. The final image, made by resampling Figs. 5 and 7 in perspective, and compositing Fig. 8 on top.
Fig. 11. Similar to Fig. 10 but without sky illumination bras, using Fig. 4 instead of Fig. 7.
and penum-
Fig. 13. A composite image in a partial eclipse, using Fig. 12 as the ground raster. Due to the limitations in the color printing process, the subtle differences between Fig. 13 and Fig. 10 are not visible here.
Note the effects of the sky illumination in Figs. 7 and 10. The ground is darker in the center of the shadow of a large clump of trees than it is in shadows of branches or smaller clumps. Also, the ground is brighter far from the shadow of the trees than it is nearer to the shadow, but
still outside its penumbra from the sun. Compare this with Fig. 11, without these effects and without penmbras. Figures 12 and 13 show the effects of a partial eclipse on the penumbra of the sun. I have produced a computer animated film “Sun and
NELSON MAX
228
Fig. 14. The square R in the sky plane S projects radially to the spherical polygon, F = ABCD on the hemisphere H, which in turn projects vertically to the area K = A’B’C’D’.
Fig. 12. The ground raster for an eclipsed sun. Note crescent shaped spots of light.
Shade” [19], showing the illumination change in a grove of trees as the moon or clouds cover the sun. I used a 48 frame cycle of images of the trees, without ground or sky, showing the leaves fluttering in the breeze. I rendered the foreground images twice, once including point source illumination at the sun, and once with only diffuse illumination from the sky. I also computed the shadow mask raster m(i, j) for each frame of the cycle, and its Fourier transform. Then for each frame of the longer final sequence, the background image is found as above, using the precomputed Fourier transform M of the mask m. The fraction S of the direct sunlight which passes through the clouds is used to interpolate between the two foreground images, which are then combined with the background using the precomputed (Y channel. Finally, the atmospheric illumination is multiplied by S, and added to the image. The total cost on one processor of a Cray XMP was on the order of 1 min per 640 x 480 pixel frame, and was dominated by the anti-aliased resampling algorithm in [13]. This cost is amortized over the cyclic repeats, while the image interpolation and compositing costs are small. ERROR ANALYSIS
The integral of Eq. (1) is over an infinite plane S, while the sum of Eq. (2) approximates the integral only over a
finite square R in this plane. The size of this square is determined by the angle (Y between the vertical line ON, and the line OM from 0 to the center of a side of R, as shown in Fig. 14. To estimate the error introduced by neglecting the region of S outside R, we assume Lambet-t’s law, so thatf,(V, L) = 1, and that the sky illumination is also 1, independent of L. The integral (1) then takes the form Z,,(a) = j)os4 =
8dA
I F cos 8do
=I dB. K
Here F is the radial projection of R onto the unit hemisphere H, as shown in Fig. 14, and K is the vertical projection of F onto the ground plane G. The quadrilateral F is bounded by great circle arcs through the points A, B, C, and D, so K is bounded by ellipses through their projections A’, B’, C’, and D’. The last two integrals are equal because the area element dB on the ground plane is equal to cos 8dw. (See the caption to Fig. 2.) Now the last integral is simply the area of the region K. Figure 15 shows K divided into eight equal elliptical sectors like OPC’. The ellipses are tangent to the unit circle, with minor axis OP = sin (Y. Let U be the point on the unit circle above C’, and let 4 be the angle POW. The sector OPC’ is formed by vertically shrinking the circular sector OTU, of area 44, by a factor sin (Y, so this area is t+ sin (Y, and also QO = QC’ = QU sin (Y. Then, tan 4 = VU/V0 = QOIQU = QU sin alQU = sin cy, so $I = tan-‘(sin a). Thus the area of G is lo(a) = 8(&tan-‘(sin a)sin a) = 4sin (Y tan-‘(sin a). When (Y = 90”, Z,,(a) =
UNIFIED
229
SUN AND SKY ILLUMINATION
EXTENSIONS TO MULTIPLE SHADOW CASTING PLANES
Fig. 15. A top view of the area A’B’C’D’, which is bounded by four ellipses, and is divided by horizontal, vertical, and 45” diagonal lines into eight equal elliptical sectors like OPC’.
4tan-‘(I) = 7~, the area of the unit circle. Figure 16 shows a graph of the fraction Zo(a)17r of the unit circle area covered by K. A related curve appears as Fig. 3 in Reeker et al. [20]. In the examples in this paper, I used CY= 70”, and included 90% of a uniform sky illumination. For pixels below the center of R, this graph also measures the error from extending the mask m(i, j) periodically outside R, but for pixels farther out, the error from assuming periodicity is larger. When the sun is not obscured by clouds, and lies inside R, its energy outweighs that from the sky, and the errors discussed above become a smaller proportion of the total illumination on the ground.
Image correlation by fast Fourier transforms is an efficient technique for computing shadows on the ground from environmental illumination, provided the shadow casting objects lie in or near a single horizontal plane. Mask rasters ml and m2, in two planes at different heights, combine with the sky raster e in a more complex formula, for which there is no obvious FFT speedup. However, if we assume that the two masks are uncorrelated, we can approximate the composition of their two shadowing effects by computing the two ground intensity rasters g as above, multiplying them together, and dividing by the unobstructed intensity Z,. Takita et al. [IO] represent their trees by several horizontal raster images with color and transparency, and one such vertical image, which is always turned to face the eye. The cone toward the sun intersects the plane of one these images in an ellipse, and the fraction of pixels inside this ellipse which are transparent gives the penumbra illumination. If there are several planes, Takita et al. take the combined penumbra illumination to be the minimum of the results for the individual planes, while the uncorrelated assumption above would take the product. If the several rasters were instead sampled by a fixed or jittered grid of rays in the cone, the transparencies for each ray could be multiplied separately, giving a better estimate for the correlation. The cone computation can produce the variably sized penumbra along the trunk of the tree, as well as the correctly sized penumbras for each level of leaves. However, it requires extensive computations at each pixel. The algorithm here, applied to several leaf levels assumed to be uncorrelated, produces less accurate results much more rapidly. ACKNOWLEDGMENTS This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract Number W-7405ENG-48. The idea for this algorithm arose from discussions about shadows with Chuck Grant. I thank Brian Cabral, Becky Springmeyer, Chuck Grant, Michael Gwilliam, and Jeff Kallman for helpful suggestions on the exposition, Jules Bloomenthal for the tree trunk data and bark texture, Becky Springmeyer and Chuck Grant for the typing and layout, Ted Einwohner for Fig. 16, and Brian Cabral for providing the cloud texture in Fig. 5. The section on error analysis was recommended by an anonymous referee.
REFERENCES ,
20
40
60
80
Fig. 16. A graph of the function J,,((Y)/~, showing the fraction of a uniform sky’s energy included in a square subtending an angle of 2a.
I. R. L. Cook, T. Porter, and L. Carpenter, Distributed ray tracing, Compuf. Graphics 18, No. 3, 1984, 137-145. 2. J. T. Kajiya, The rendering equation, Comput Graphics 20, No. 4, 1986, 143-150.
230
NELSON
3. N. Greene, Applications of world projections, IEEE Cotnput. Graphics Appl. 6, No. 11, 1986, 21-29. 4. G. S. Miller and R. C. Hoffman, Illumination and Reflection Maps: Simulated Objects in Simulated and Real Environments, SIGGRAPH ‘84 Advanced Computer Graphics Animation Course Notes, July 1984. 5. B. Cabral, N. Max, and R. Springmeyer, Bidirectional reflection functions from surface bump maps, Comput. Graphics 24, No. 4, 1987, 273-281. 6. M. F. Cohen and D. P. Greenberg, A radiosity solution for complex environments, Comput. Graphics 19, No. 3, 1985, 31-40. 7. D. S. Immel, M. F. Cohen, and D. P. Greenberg, A radiosity method for non-diffuse environments, Comput. Graphics 20, No. 4, 1986, 133-142. 8. T. Nishita and E. Nakamae, Half tone representation of 3-D objects illuminated by area sources or polyhedron sources, IEEE Proc. COMSAC, 1983, 237-241. 9. T. Nishita and E. Nakamae, Continuous tone representation of three dimensional objects illuminated by sky light, Comput. Graphics 20, No. 4, 1986, 125-132. 10. S. Takita, K. Kaneda, T. Akinoba, H. Iriyama, E. Nakamae, and T. Nishita, A simple rendering for penumbra caused by sunshine, in Proceedings, CG International 90; New Advances in Computer Graphics, Springer-Verlag, New York, 1990.
MAX 11. R. C. Gonzalez and P. Wintz, Digital Image Processing, AddisonWesley, Reading, MA, 1977. 12. N. Max, Atmospheric illumination and shadows, Comput. Graphics 20, No. 4, 1986, 117-124. 13. N. Max, Anti-aliasing scan line data, IEEE Comput. Graphics Appl. 10, No. 1, 1990, 18-30. 14. G. Y. Gardner, Visual simulation of clouds, Comput. Graphics 19, No. 3, 1985, 297-303. 15. J. F. Blinn, Light reflection functions for simulation of clouds and dusty surfaces, Comput. Graphics 16, No. 3, 1982, 21-29. 16. E. Catmull and A. R. Smith, 3-D transformations in scanline order, Comput. Graphics, 14, No. 3, 1980, 279-285. 17. K. M. Fant, A nonaliasing, real time spatial transform technique, IEEE Computer Graphics Appl. 6, No. 1, 1986, 71-80 (see also Letters to the Editor in 6, Nos. 3 and 7). 18. T. Porter and T. Duff, Compositing digital images, Comput. Graphics 18, No. 3, 1984, 253-259. 19. N. Max, Sun and Shade, SIGGRAPH Video Review, Issue No. 36, Segment 8, 1988, ACM, available from First Priority, P.O. Box 576, Itasca, IL 60143. 20. R. J. Reeker, D. W. George, and D. P. Greenberg, Acceleration techniques for progressive refinement radiosity, Comput. Graphics 24, No. 2, 1990, 59-66.