Features
FeaturesHere's a quick way to display a textured surface that nevertheless holds its own against more time-consuming techniques.
Introduction
Lambert Shading [1] is a technique traditionally applied to regular solids (cylinders, spheres, tetrahedrons, etc.) or surfaces composed of collections of regular solids (e.g. triangular meshed landscapes) to quickly generate images of these objects. Lambert Shading (LS) is among the quickest and simplest methods used for rendering simple 3-D objects into images. A Lambertian surface has the property that the intensity at any point depends only on the direction of the illumination and the surface slope. The intensity is independent of the observer's location. More sophisticated methods such as ray tracing [1] are usually employed for rendering complex objects or scenes containing elements such as reflections and shadows. The disadvantage of ray-tracing and similar techniques is the time required to render a scene. For complicated images, this time can be on the order of hours or days, even when the fastest machines are employed.
Lambert Shading can provide a quick and very satisfactory alternative to the more sophisticated rendering techniques in the case of large, complex, irregular surfaces. The time required to render such a surface is on the order of seconds or minutes. To utilize LS the surface should ideally be large and digitized. A surface can be digitized if it is adequately represented by sampling and storing its elevations at some (usually) regular interval. This surface is considered large if its sampled dimensions are on the order of the number of pixels used to display its rendered image. For such surfaces, the results of LS rendering can be quite dramatic, yielding near-picture quality images.
A good example of surfaces that are irregular, large, and digitized is Digital Terrain Elevation Data (DTED). Digital terrain data is available for free from the United States Geological Survey (USGS) via their ftp site: edcftp.cr.usgs.gov. USGS has digitized and stored almost all the United States; each DTED image represents a one-square degree area of the earth's surface sampled at 3 arc-seconds (about 100 meters). Each file is compressed in gzip format. When uncompressed these ASCII formatted files expand to about 10MB. Each database contains 1201 X 1201 terrain elevations measured in meters relative to mean sea level.
Figure 1 shows an LS rendered image of the DTED whose southwest corner is located at 34º N Latitude, 120º W Longitude. This area contains Camarillo California, the Oxnard Plain (lower right), and part of the Channel Islands (lower left). The Pacific Ocean is colored black to provide sharper definition of the California coast line.
Lambert Shading Fundamentals
Before presenting the coded algorithm, We'll review the fundamentals of LS.
Lambert Shading uses surface gradient to determine the pixel shading at each point. For digital surfaces, each elevation, z, is mapped to a pixel. In the DTED case, the resulting image is therefore 1,201 X 1,201 pixels. The pixel shading or intensity, I, at point (x,y) is calculated as the dot product of two vectors: i, the illumination vector and n(x,y), the local surface normal (gradient):
I(x,y) = n(x,y) * i (Equation 1)The source of illumination is assumed to be at infinity, so the direction to the source is the same anywhere on the surface. This assumption allows the source vector to be treated as constant in the calculation, significantly increasing the rendering speed. This illumination vector is expressed in terms of an elevation angle, l, and azimuth, u, which represent the direction of the vector relative to the x-y plane. It may be helpful to think of the light source as positioned on the surface of an infinitely large sphere, at elevation l and azimuth u. The light source points toward the center of the sphere, where the x-y plane resides. A normal through the origin of the x-y plane passes through (0, 0).
Calculating the Dot Product
Calculating the dot product in Equation 1 requires that i and n(x,y) be expressed in the same terms. The following equations express the illumination and surface-normal vectors in terms of the Cartesian frame of elevations, where e1, e2, and e3 represent unit vectors in the x, y, and z directions respectively.
i = e1os(l)cos(u) + e2cos(l)sin(u) + e3sin(l) (Equation 2) n(x,y) = e1GRADxz(x,y) + e2GRADyz(x,y) - e3 (Equation 3)The symbol GRAD is the gradient operator with the x and y subscripts indicating direction.
Applying Equation 1 is now straightforward. The angles (l, u) are often chosen as 45º, representing illumination from the northeast (upper right). The coordinates (x,y) are discrete and mapped directly to pixels. In the DTED example, 0 < x,y < 1200. The gradients GRADxz(x,y) and GRADyz(x,y) are approximated as the difference between adjacent elevations. At this point, we introduce two scaling factors, Io and Smax, to take maximal advantage of the range of display intensities (usually 0 - 255). The following equation is equation 1 restated in terms of equations 2 and 3. It shows how to compute the intensity at each pixel:
I(x,y) = Ix[z(x+1,y) - z(x-1,y)] + Iy[z(x,y+1) - z(x, y-1)] + Io (Equation 4)where
Ix = (Imax - Io)Smax20.5cos(l)cos(u), Iy = (Imax-Io)Smax20.5cos(l)sin(u), Io = Imax20.5sin(l)/2,Imaxis the maximum intensity, and the factor 20.5 insures Io = Imax/2 when l = 45° (see below).
Equation 4 can be directly implemented. All that remains is to determine Smax and Io. Since, on average, the slope of the surface is expected to be neutral, we select Io as the mid-range intensity (usually 127). Computing Smax is more involved, since it depends on the magnitude of the slopes. For a smooth surface, the slope magnitude range is small and intensities cluster around a midpoint — the lowest and highest intensities remain unused. For rough surfaces the opposite is true; slope magnitude variation is large and pixel intensities tend to the extremes of the range. Therefore, Smax should be chosen to evenly spread pixel intensity through the available range of display intensities.
Smax will yield an even distribution if it is taken as the slope of the line that best fits the normalized cumulative histogram of surface slopes. (In this case, "best fits" means as computed by least squares, and the surface slopes are basically elevation differences.) The normalized cumulative histogram is represented by the symbol Hz.
Surface slopes are rarely greater than 45°; this makes it simple to choose the histogram bin size and maximum number of bins. We've found that it usually works well to choose a bin size (D) equal to elevation resolution (1 meter), and a maximum number of bins (Nz) equal to twice the range of the elevation spacing (100 meters). For different surfaces you may need to adjust these values to better represent the data.
Figure 2 illustrates how to calculate the best-fit line to the normalized cumulative histogram; we can then use the slope of this line to determine the best value for Smax. The histogram in Figure 2 was computed from the elevations in Figure 1. In calculating this best-fit line, we apply the constraint that it must pass through the center of the histogram (bin index 0) at a normalized elevation of 0.5. The equation for this line then becomes:
y(x) = (1 + sx)/2 (Equation 5)where s/2 is the line's slope and 0 < x < Nz. The line that best fits the histogram has a slope such that
S ((x) - Hz[x])2is a minimum.
The minimum of this expression is easily found by taking its derivative with respect to s, setting the derivative equal to zero, and solving for s/2. This is the desired value of Smax. The details of computing Smax appear in Listing 2, in the function surfaceFindSmax.
Implementation
Listing 1 contains the header for the CSurface base class used to render an irregular surface using Lambert Shading. The class contains data members describing the dimensions of the surface, pointers to the surface data, and the bitmap image created by the LS process. CSurface also contains the cumulative elevation difference histogram and Smax as data members. We use typdefs to define the surface elevation and bitmap data types; changing these typedefs allows support for a variety of data types. For the DTED image, the surface elevations are stored as short integers. Bitmaps are almost always stored as characters or unsigned characters.
The public member functions include a constructor for initializing data members to zero, a destructor for freeing allocated memory, and virtual read and draw functions. The read function takes as a single argument the name of the file containing the surface data. The base-class implementation of this function does nothing; you must override it in a derived class to read in the elevations of the desired surface. At a minimum, this function should allocate memory for the elevations and read the data into this storage. In addition it should initialize the data dimensions and spacing members, allocate memory for the bitmap, and initialize the bitmap buffer size. We've also made the draw member function virtual, since it is obviously implementation dependent. We've typedef'd the arguments to this function to support MS-Windows, X-Windows, and Macintosh development. For MS-Windows development redefine _HANDLE as CDC*. For X-Windows _HANDLE becomes Widget and _POINTER becomes XtPointer.
In 16-bit Windows environments you should #define _HUGE as the keyword huge to support large surfaces. In dealing with DTEDs this is a must since the elevation and bitmap pointers refer to memory blocks well in excess of the 65K far pointer limit.
The private member functions perform the Lambert Shading of the surface elevations. Function surfaceLSRender calls surfaceFindSmax to compute Smax and render the surface. Upon completion of surfaceLSRender the image bitmap is filled and ready for display. The member function surfaceFindSmax calls surfaceHistogram to compute the cumulative elevation difference histogram required to compute Smax. A more efficient approach would be to compute this histogram while reading in the elevation data.
After computing Smax, the intensity components Ix, Iy, and Io are computed. These constants are scaled to allow integer arithmetic in the nested for loops used to calculate the pixel intensities of the image. This scaling operation greatly increases performance when the elevations are integers. The while loop in the following code fragment determines the integer scaling factor, FAC:
z = ( x < y ? x : y); FAC = 1; while( z < 1.0) { z = 10.0*z; FAC = 10*FAC; }In the above snippet the local variable z is set to the minimum of Ix and Iy (called x and y at this point) and the scale factor FAC is set to one. The while loop increases both z and FAC by a factor of ten until z is greater than one.
The above procedure insures that precision is maintained by the scale factor. The intensity components are then scaled by FAC. In calculating pixel intensities, the algorithm effectively replaces the otherwise floating-point multiplication of Ix, Iy, and Io by integer multiplication and an integer divide, which can be significantly faster.
Upon finishing the calculation of the intensity components, the algorithm sets pointers to the first and third rows and columns of the elevations and enters nested for loops to compute the pixel intensities of surface elevations. Intensities are computed for all elevations except those on the edge; these pixels are filled in later by copying the values of adjacent interior pixels. Edge points are excluded since they lack surrounding elevation on either side to compute a gradient. The gradient calculation could be modified and intensities computed for the edge, but the result is not worth the effort since copying the adjacent pixels is unnoticeable in the final image.
How to Apply Lambert Shading
In implementing and applying the LS algorithm, you have a wide range of choices in terms of what you select for parameter values, number of bits, etc. It's worthwhile going over some of the choices and their implications.
First, a few words on the general characteristics of LS images. Flat regions in the image end up with intensity Io; very smooth surfaces will tend to have the majority of their intensities near this value also. The contrast in the resulting image will be poor. Rough surfaces tend to provide crisper and sharper images, since the resulting intensities are naturally better distributed. Low contrast effects are not always undesirable however. Since Lambert Shading is essentially a derivative process and derivatives are very sensitive to noise, "defects" or anomalies in the elevations become exaggerated by the scale factor Smax and are readily apparent in the final image. For instance, USGS originally created DTEDs by digitizing paper contour maps. Much of this data has subsequently been replaced with more accurate satellite measurements. However, in some of these older paper map-derived databases the contours are visible in the rendered image. Step-like areas show up especially in smooth terrain. In some cases, quality control was poor and a column of elevation data was missing (set to zero). This missing column is easily visible in the LS image.
Manipulating the illumination angles (l,u) creates both interesting and subtle effects on the final image. The elevation angle, l, has little effect on the image except at very low (0°) or very high (90°) incidences. At 90°, for example, Ix and Iy become zero and the image appears as a smooth flat surface. For this reason, l = 45° is almost always the best choice; both the cos(l) and sin(l) terms in Ix, Iy, and Io become 2-0.5.who apply LS usually drop the cos(l) and sin(l) terms from Ix, Iy, and Io.
The longitude angle, u, creates more subtle effects in the image. If two images are generated by illuminating the surface from opposite directions (us 180° apart) the results appear almost identical. However, small concavities in the surface appear as bumps when the illumination is from the upper right (u = 45°) and as depressions when the illumination is from the lower left (u = 225°). This effect is an optical illusion; it reflects the manner in which the mind has been conditioned to view reality and is not an inherent flaw in the LS process.
Although the code presented here utilizes the full range of intensity available, you don't need to exploit the full intensity to produce good imagery. The 256 or 8-bit intensity range is more than adequate. In fact, you can generate good images with as little as 5 bits. Figure 1 is an example of 5-bit imagery. The remaining intensities can be utilized to color the image. For example, in the case of terrain, you can color the image by elevation, thereby overlaying contours on the surface in a natural, familiar, and visual pleasing way. o
Reference
[1] C. Pokorny. Computer Graphics, an Object-Oriented Approach to the Art and Science (Franklin, Beedle & Associates Inc., Wilsonville, Oregon, 1994).
Ronald E. Huss is a consultant and former employee of Hughes Aircraft Company in El Seguendo, CA. He has a Ph.D. in bio-medical engineering. Ron is the resident thaumaturge at Hughes. His research interests include solving complex optimization problems involving the exploitation of digital terrain data.
Mark A. Pumar is a Senior Engineer/Programmer at COMPTEK Federal System in Camarillo, CA. He has MS degrees in Physics and Electrical Engineering from UCLA. Mark is a programmer by circumstance instead of choice, but has managed through luck and good fortune to bluff his way, thus far, into remaining employed as such. With the initials MAP he should, and is, interested in digital mapping. Mark can be reached via e-mail at: MAPar@aol.com or MPUMAR@COMPTEK.com or by phone at (805) 987-9739.