Patent Publication Number: US-9846676-B2

Title: Method and apparatus for converting spherical harmonics representations of functions into multi-resolution representations

Description:
RELATED APPLICATIONS 
     This application is a continuation of and claims priority to U.S. patent application Ser. No. 12/354,583 filed Jan. 15, 2009 entitled “Method and Apparatus for Converting Spherical Harmonics Representations of Functions into Multi-Resolution Representations,” which claims priority to U.S. Provisional Patent Application Ser. No. 61/094,789 filed Sep. 5, 2008 entitled “Method and Apparatus for Importance Sampling Spherical Harmonics in Computer Graphics.” These disclosures are incorporated by reference herein in their entirety. 
    
    
     BACKGROUND 
     Illumination Functions 
     In 3-D computer graphics, techniques exist that attempt to capture the way in which light interacts in the real world. Such techniques may take into account, at locations on a surface in a scene, the light that comes directly from a light source (direct lighting, or direct illumination), and light rays from one or more light sources reflected by other surfaces in the scene (indirect illumination) to determine global illumination for the locations. A reflectance function may be applied to a surface to determine reflectance at locations of the surface; this reflectance may, for example, be used in determining indirect illumination at locations of the scene. Given a direction from a perspective point, and given incoming light from another direction or directions, a reflectance function may determine how much light is reflected from a surface towards the perspective point. An example of a reflectance function is the bidirectional reflectance distribution function (BRDF). 
     Environment Maps 
     In computer graphics, it is common to have functions that need to be represented or defined over a sphere, or a space. For example, a user might want to load an image and have it wrapped over a sphere, and assume that the image is lighting a scene within the sphere. A function may encode, at every pixel location in the scene, the intensity of light that is coming in from the environment represented by the image. The wrapped image may be referred to as an environment map. In this example, the environment map is used for lighting the scene within the sphere. However, environment maps may be used for other purposes than lighting in computer graphics, and other functions than a lighting function may be used with environment maps. 
     There may be various techniques to approximate these functions over a sphere. For example, the function may be given a direction, represented by the azimuth and polar angle, around the sphere. The function returns a value for the input direction. In the case of a lighting function, for a given azimuth and polar angle on the environment map, the function returns an intensity. In the case of a pixel representation with multiple channels, e.g. RGB pixels, the function may be applied separately to each pixel channel, e.g. to R, G and B separately for RGB pixels. 
     Work has been done on compressing, encoding, and representing these spherical signals or functions defined over a sphere. Two such representations are spherical harmonics and wavelets. 
     Spherical Harmonics 
     Spherical harmonics are a frequency-space basis for representing functions defined over a sphere. Spherical harmonics have been studied extensively, and have widespread applicability to many physical problems ranging from atomic electron configurations in physical chemistry, to the representation of gravitational and magnetic fields of planetary bodies in geodesy. Spherical harmonics also appear in quantum mechanics as the solutions of the Schrödinger equation in spherical coordinates. 
     In computer graphics, many quantities are naturally defined over the spherical or hemispherical domain, making spherical harmonics a natural basis for these computations. Furthermore, spherical harmonics have many convenient properties that lead to efficient implementations for convolutions, rotations, and double product integrals. Consequently, spherical harmonics have found recent popularity in computer graphics in the area of precomputed radiance transfer (PRT) techniques, and have also been used heavily in related domains such as radiative transfer. These properties and operations of spherical harmonics in the prior art are summarized in the table illustrated in  FIG. 1 , which summarizes and compares capabilities of spherical harmonics and Haar wavelets according to prior art. 
     Wavelets 
     Another class of basis functions that may be used in computer graphics is wavelets. Whereas spherical harmonics decompose functions into independent frequency bands, wavelets hierarchically decomposing functions into components that are localized in both frequency and space. Due to their simplicity, Haar wavelets have received much attention in computer graphics research. Haar wavelets, for example, have been successfully applied to PRT and radiosity computation. Furthermore, Haar wavelets support efficient sampling using a simple hierarchical warping technique. 
     Haar wavelets and spherical harmonics share many of the same strengths.  FIG. 1  shows a table that summarizes and compares capabilities of spherical harmonics and Haar wavelets according to the prior art. Haar wavelets can more efficiently approximate high-frequency functions; however, they do not support efficient rotations and may produce distracting aliasing artifacts. Spherical harmonics, on the other hand, support efficient rotation and are efficient for representing low-frequency functions. Spherical harmonics do not suffer from the same aliasing artifacts as Haar wavelets but instead exhibit error in the form of “ringing”. 
     Importance Sampling 
     Generating samples according to the intensity distribution of a function is a useful operation in many disciplines. In Monte Carlo techniques, this takes the form of importance sampling and can be used to significantly reduce variance in numerical integration. Importance sampling techniques also have widespread use outside of Monte Carlo rendering. For example, importance sampling techniques can be used in non-photorealistic rendering, artificial stippling, and procedural geometry placement. It may be beneficial to also enforce certain quality measures on the resulting sample distribution in order to reduce variance further or simply to obtain more visually pleasing distributions. 
     Mip Maps 
     In computer graphics, a mip map (also referred to as MIP map or mipmap) is a technique for storing multiple resolutions of a single image. Typically this includes the original high-resolution image, plus successive copies of the image reduced in resolution by a factor of 2 in each dimension. This collection of images forms an “image-pyramid” where the level in the pyramid determines the resolution and level-of-detail of the image. 
     SUMMARY 
     Embodiments of a method and apparatus for converting spherical harmonics representations of functions to multi-resolution representations, such as wavelet representations and mip map representations, are described. In computer graphics, reflectance functions on surfaces tend to be smooth and low-frequency; thus, these functions are effectively stored as spherical harmonic representations. On the other hand, wavelets are a good basis for environment maps because environment maps tend to be high-frequency. Embodiments may implement a method to efficiently convert a function that is stored in spherical harmonics into a function that is stored in a multi-resolution representation. Embodiments may implement a method to efficiently convert a function that is stored in spherical harmonics into a function that is stored in wavelets, for example Haar wavelets. Implementations of this method perform the conversion from a spherical harmonic representation to a wavelet representation efficiently and analytically. Implementations of this method may also perform conversion from a spherical harmonic representation to a mip map representation. 
     The analytic method for converting from a spherical harmonic representation to a wavelet or mip map representation enables embodiments to perform a variety of computer graphics tasks, and to perform computer graphics tasks more efficiently than could be done in conventional techniques. For example, one or the other of the environment map or the reflectance function may need to be rotated so that the two can be aligned. Since a spherical harmonic representation can be efficiently rotated, a function may be stored as a spherical harmonic representation and rotated in the spherical harmonic domain; the function can then be converted to a wavelet representation, for example a Haar wavelet representation if necessary. This may all be performed efficiently throughout the rendering pipeline. As another example, there is no conventional, efficient method for importance sampling spherical harmonics. The Haar wavelet domain does have efficient methods for importance sampling. Since there are methods to importance sample wavelets (e.g., Haar wavelets), embodiments may use an analytic conversion method to take the spherical harmonic function, convert the spherical harmonic function to Haar wavelets, and then apply an importance sampling technique for Haar wavelets to the Haar wavelet representation to generate a set of importance samples for the function. The end result is a set of sample rays in different directions over a sphere or hemisphere. 
     Thus, embodiments enable quick and efficient rotation since spherical harmonics representations can be efficiently rotated, provide a method to analytically convert from one representation to another on the fly, and enable importance sampling of spherical harmonics by efficiently converting from a spherical harmonic representation to a Haar wavelet representation and performing the importance sampling on the Haar wavelet representation. Generally speaking, embodiments may store different functions in different representations, and combine them together in different ways, with the analytic conversion method used to efficiently convert from one basis to the other. 
     In addition, the methods for conversion of function representations and for importance sampling may be used in embodiments to importance sample the product of two Haar wavelets. Some conventional techniques store both the lighting (environment map) and the reflectance in the Haar wavelet domain, and thus may importance sample the product of the two. However, these conventional techniques suffer from various limitations; for example, there is no efficient way to rotate a Haar wavelet representation. Embodiments may directly use the analytic conversion method in this case. One of the two functions, for example a reflectance function, may be represented in the spherical harmonic domain. The analytic conversion method may be applied to perform a fast conversion of the function into the Haar wavelet domain, and an importance sampling technique may then be applied which samples the product of the two functions in the Haar wavelet domain. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  shows a table that summarizes and compares capabilities of spherical harmonics and Haar wavelets according to prior art. 
         FIG. 2  shows a table that summarizes and compares capabilities of spherical harmonics and Haar wavelets according to embodiments. 
         FIGS. 3A through 3C  graphically illustrate aspects of importance sampling according to embodiments. 
         FIG. 3A  illustrates importance sampling a reflectance function according to one embodiment. 
         FIG. 3B  illustrates Haar wavelet importance sampling on an environment map according to embodiments. 
         FIG. 3C  illustrates importance sampling the product of a spherical harmonic function and a Haar wavelet function according to one embodiment. 
         FIG. 4  illustrates visualizations of the first nine spherical harmonic basis functions in Cartesian and spherical coordinates. 
         FIG. 5  illustrates a hierarchical sample warping method according to one embodiment. 
         FIG. 6  illustrates a transition between (b) and (c) of  FIG. 5  according to one embodiment. 
         FIG. 7  illustrates the separation of spherical harmonic basis functions. 
         FIGS. 8 and 9  illustrate examples of pseudo-code for recurrence procedures that may be used in embodiments. 
         FIGS. 10A and 10B  are flowcharts of methods in importance sampling spherical harmonics according to embodiments. 
         FIG. 10A  is a flowchart of a method of performing rotation and importance sampling according to one embodiment. 
         FIG. 10B  is a flowchart of a method of importance sampling the product of two functions according to one embodiment. 
         FIG. 11  illustrates a function representation conversion module in a rendering system according to one embodiment. 
         FIG. 12  illustrates an example of a computer system that may be used in embodiments. 
     
    
    
     While the invention is described herein by way of example for several embodiments and illustrative drawings, those skilled in the art will recognize that the invention is not limited to the embodiments or drawings described. It should be understood, that the drawings and detailed description thereto are not intended to limit the invention to the particular form disclosed, but on the contrary, the intention is to cover all modifications, equivalents and alternatives falling within the spirit and scope of the present invention. The headings used herein are for organizational purposes only and are not meant to be used to limit the scope of the description. As used throughout this application, the word “may” is used in a permissive sense (i.e., meaning having the potential to), rather than the mandatory sense (i.e., meaning must). Similarly, the words “include”, “including”, and “includes” mean including, but not limited to. 
     DETAILED DESCRIPTION OF EMBODIMENTS 
     In the following detailed description, numerous specific details are set forth to provide a thorough understanding of claimed subject matter. However, it will be understood by those skilled in the art that claimed subject matter may be practiced without these specific details. In other instances, methods, apparatuses or systems that would be known by one of ordinary skill have not been described in detail so as not to obscure claimed subject matter. 
     Some portions of the detailed description which follow are presented in terms of algorithms or symbolic representations of operations on binary digital signals stored within a memory of a specific apparatus or special purpose computing device or platform. In the context of this particular specification, the term specific apparatus or the like includes a general purpose computer once it is programmed to perform particular functions pursuant to instructions from program software. Algorithmic descriptions or symbolic representations are examples of techniques used by those of ordinary skill in the signal processing or related arts to convey the substance of their work to others skilled in the art. An algorithm is here, and is generally, considered to be a self-consistent sequence of operations or similar signal processing leading to a desired result. In this context, operations or processing involve physical manipulation of physical quantities. Typically, although not necessarily, such quantities may take the form of electrical or magnetic signals capable of being stored, transferred, combined, compared or otherwise manipulated. It has proven convenient at times, principally for reasons of common usage, to refer to such signals as bits, data, values, elements, symbols, characters, terms, numbers, numerals or the like. It should be understood, however, that all of these or similar terms are to be associated with appropriate physical quantities and are merely convenient labels. Unless specifically stated otherwise, as apparent from the following discussion, it is appreciated that throughout this specification discussions utilizing terms such as “processing,” “computing,” “calculating,” “determining” or the like refer to actions or processes of a specific apparatus, such as a special purpose computer or a similar special purpose electronic computing device. In the context of this specification, therefore, a special purpose computer or a similar special purpose electronic computing device is capable of manipulating or transforming signals, typically represented as physical electronic or magnetic quantities within memories, registers, or other information storage devices, transmission devices, or display devices of the special purpose computer or similar special purpose electronic computing device. 
     Various embodiments of a method and apparatus for converting spherical harmonics representations of functions to multi-resolution representations, such as wavelet representations and mip map representations, are described. Given a valid probability density function (PDF) represented as a vector of spherical harmonic coefficients, embodiments of the method may warp an input point set to match the target PDF using hierarchical sample warping. This approach is efficient and produces high quality sample distributions. As a by-product of the sampling procedure, embodiments may produce a multi-resolution representation of the spherical harmonic function as either a spherical mip map representation or wavelet representation (e.g., Haar wavelet representation). By exploiting this implicit conversion, the method may be extended in embodiments to distribute samples according to the product of a spherical harmonic function with a spherical mip map or a wavelet (e.g., a Haar wavelet). This generalization has immediate applicability in rendering, in, for example, importance sampling the product of a bidirectional reflectance distribution function (BRDF) and an environment map where the lighting (the environment map) is stored as a single high-resolution wavelet and the BRDF is represented in spherical harmonics. Since spherical harmonics can be efficiently rotated, this product can be computed on-the-fly even if the BRDF is stored in local-space. Thus, embodiments may seamlessly combine the strengths of spherical harmonics and wavelets, and, compared to conventional techniques, embodiments may significantly reduce precomputation time and storage requirements. 
     There are no conventional techniques for efficiently importance sampling spherical harmonics. Embodiments address this and other limitations of conventional techniques by implementing a number of techniques as described herein to bridge the gap between the capabilities of spherical harmonics and Haar wavelets. 
     Embodiments may provide a high-quality importance sampling method for spherical harmonics. As described herein, the importance sampling method for spherical harmonics can be used to sample the product of a spherical harmonic function and a Haar wavelet function, thus bridging the gap between spherical harmonics and Haar wavelets. This product can be evaluated efficiently on-the-fly with little or no precomputation, and allows for high-frequency content in one function, and efficient rotation of the other. Such a situation is common in compute graphics when, for example, convolving a smooth bidirectional reflectance distribution function (BRDF) with a high frequency environment map. Embodiments may provide one or more of, but not limited to, the following:
         An analytic conversion between spherical harmonics and multi-resolution representations (spherical mip maps/wavelets).   A method to rotate and align a function represented by spherical harmonics with, for example, a function represented by wavelets.   An importance sampling method for spherical harmonic functions.   A method for sampling the product of spherical harmonics and wavelets.       

     In computer graphics, the reflectance functions on surfaces tend to be smooth and low-frequency; thus, these functions are effectively stored as spherical harmonic representations. On the other hand, wavelets are a good basis for environment maps because environment maps tend to be high-frequency. 
     Embodiments may implement a method to efficiently convert a function that is stored in spherical harmonics into a function that is stored in wavelets, for example Haar wavelets. See, for example, equations 6 and 10 provided later in this document. Implementations of this method may perform the conversion from a spherical harmonic representation to a Haar wavelet representation efficiently and analytically. 
       FIG. 2  shows a table that summarizes and compares capabilities of spherical harmonics and Haar wavelets according to embodiments. Note that, in contrast to the table illustrated in  FIG. 1 ,  FIG. 2  shows importance sampling for spherical harmonics, a capability that is enabled by embodiments as described herein. The efficient analytic method for converting from a spherical harmonic representation to a Haar wavelet representation enables embodiments to perform a variety of computer graphics tasks, and to perform computer graphics tasks more efficiently than could be done in conventional techniques. For example, one or the other of the environment map or the reflectance function needs to be rotated so that the two can be aligned. Since a spherical harmonic representation can be efficiently rotated, a function may be stored as a spherical harmonic representation and rotated in the spherical harmonic domain; the function can then be converted to a Haar wavelet representation if necessary. This may all be performed efficiently throughout the rendering pipeline. As another example, there is no conventional, efficient method for importance sampling spherical harmonics. The Haar wavelet domain does have efficient methods for importance sampling. Since there are methods to importance sample Haar wavelets, embodiments may use the analytic conversion method to take the spherical harmonic function, convert the spherical harmonic function to Haar wavelets, and then apply an importance sampling technique for Haar wavelets to the Haar wavelet representation to generate a set of importance samples for the function. The end result is a set of sample rays in different directions over a sphere or hemisphere. 
     Thus, embodiments may enable quick and efficient rotation since spherical harmonics representations can be efficiently rotated, may provide a method to analytically convert from one representation to another on the fly, and may enable importance sampling of spherical harmonics by efficiently converting from a spherical harmonic representation to a Haar wavelet representation and performing the importance sampling on the Haar wavelet representation. Generally speaking, embodiments may store different functions in different representations, and combine them together in different ways, with the analytic conversion method used to efficiently convert from one basis to the other. 
     In addition, the methods conversion of function representations and for importance sampling may be used to importance sample the product of two Haar wavelets. Some conventional techniques store both the lighting (environment map) and the reflectance in the Haar wavelet domain, and thus may importance sample the product of the two. However, these conventional techniques suffer from various limitations; for example, there is no efficient way to rotate a Haar wavelet representation. Embodiments may directly use the analytic conversion method in this case. One of the two functions, for example a reflectance function, may be represented in the spherical harmonic domain. The analytic conversion method may be applied to perform a fast conversion of the function into the Haar wavelet domain, and an importance sampling technique may then be applied which samples the product of the two functions in the Haar wavelet domain. 
       FIGS. 3A through 3C  graphically illustrate aspects of importance sampling according to embodiments. Note that the white or gray dots on the images represent sample points. Embodiments provide a method for importance sampling spherical harmonic functions such as the bidirectional reflectance distribution function (BRDF) illustrated in  FIG. 3A .  FIG. 3B  illustrates Haar wavelet importance sampling on an environment map. Embodiments may also be combined with wavelets to importance sample the product of a spherical harmonic function and a Haar wavelet function, as illustrated in  FIG. 3C . 
     Conventional Sampling Techniques 
     The following briefly summarizes some conventional sampling techniques, focusing on importance sampling in the context of rendering, and contrasts the conventional techniques with embodiments as described herein. 
     Bidirectional Reflectance Distribution Function (BRDF) Sampling 
     A commonly used conventional technique for increasing the efficiency of ray tracing based algorithms is to importance sample the reflected rays based on the scattering function. Several BRDF models can be directly importance sampled, including the Phong model, the Ward model, the Lafortune model, and many others. Complex BRDF models such as the Torrance-Sparrow and Cook-Torrance models cannot be analytically inverted and require numerical approximations. However, parametric models which are directly invertible, such as the Ward, Lafortune, and half-way vector disk BRDFs have been successfully used to approximate complex measured BRDF data. Several researchers have addressed the problem of importance sampling measured BRDFs. Though spherical harmonics can be used to efficiently represent low frequency BRDFs and indirect lighting, no conventional BRDF sampling techniques have considered spherical harmonics since practical importance sampling techniques have not existed for this basis. 
     In radiative transfer within participating media, phase functions are generally smooth and commonly represented using spherical or zonal harmonics. Unfortunately, since spherical harmonics do not support importance sampling, most techniques rely on simpler, analytically invertible, models such as the Henyey-Greenstein and Schlick phase functions. Multi-lobe variants of these phase functions allow for very complex scattering distributions; however, though these phase functions can also be sampled, the required splitting process destroys stratification and the resulting distributions are of low quality. 
     Embodiments of the spherical harmonic importance sampling method may, for example, be used on arbitrarily complex distributions, and may preserve the variance reduction properties of stratification. Furthermore, embodiments are particularly well suited for phase functions, since the spherical harmonic importance sampling method gains further efficiency when using only the zonal subset of spherical harmonics. 
     Environment Map Sampling 
     Several conventional techniques have been developed for sampling complex distant lighting by placing pre-integrated direction lights at the brightest locations in environment maps. These methods may work well within diffuse scenes but become increasingly inefficient on glossy and specular (mirror-like) surfaces. This issue has been addressed by using spherical harmonics to pre-filter the environment map according to the BRDF. However, since efficient sampling techniques have not been available for spherical harmonics, these methods could not account for occlusions from other objects. 
     Embodiments of the spherical harmonic importance sampling method directly address this limitation of spherical harmonics. 
     Product Sampling 
     Conventional approaches have also been developed which attempt to use information about both the BRDF and the lighting in the sampling process. For example, multiple importance sampling can be used to optimally combine several sampling strategies. Approaches have been developed which try to distribute samples according to the product of the lighting and the BRDF. The simplest techniques first draw samples according to one of the terms and then adjust these samples to match the product distribution. This can be done using importance resampling, or by rejection sampling; however, these become ineffective if both functions have high peaks. Furthermore, none of these techniques can preserve stratification. 
     A wavelet importance sampling approach has been developed which can directly generate samples according to the product of two or more functions. This approach employs a sampling procedure that hierarchically warps a uniform input point set to match the wavelet product. Since warping is performed hierarchically and in small increments, stratification properties of the input set are preserved in the output distribution. Unfortunately, since the two functions need to lie in the same coordinate system and no efficient rotation exists for Haar wavelets, expensive pre-rotation is necessary. Due to memory constraints, this limits the method to relatively low-resolution lighting, and still may require hours of precomputation. Recent work has tried to address this limitation by either heuristically splitting the environment map based on peaks in the BRDF or by constructing a wavelet approximation of the BRDF on-the-fly. 
     Embodiments of the spherical harmonic importance sampling method rely on a sample warping procedure to efficiently produce high quality distributions for spherical harmonics. Embodiments may also be applied to product sampling, and address the limitations of wavelet importance sampling by relying on efficient rotations in the spherical harmonic domain. Thus, embodiments may seamlessly combine and exploit the complementary properties of spherical harmonics and Haar wavelets. 
     An algorithm for importance sampling the product of distant environment map lighting with an orientable clamped-cosine lobe has been developed. This algorithm creates a triangulated representation of the environment map and stores the illumination premultiplied by each of the first nine spherical harmonic basis functions at every vertex. This forms a steerable basis where the clamped-cosine can be efficiently rotated to any orientation. Sampling this triangulation is done using an approach much like a hierarchical sample warping technique. 
     Embodiments of the spherical harmonic importance sampling method have several advantages when compared to the above-described algorithm. For example, the algorithm for importance sampling the product of distant environment map lighting with an orientable clamped-cosine lobe only considers a circularly symmetric cosine lobe, which is approximated using order 3 spherical harmonics. In contrast, embodiments of the spherical harmonic importance sampling method work for arbitrary spherical harmonic functions of any order. Though generalizing the other approach to higher orders is possible, since a full 2D environment triangulation must be stored for each additional spherical harmonic basis, memory usage becomes an issue. In contrast, embodiments of the spherical harmonic importance sampling method require no additional storage if importance sampling a single spherical harmonic function, and only utilize a 1D mip map per spherical harmonic basis for product sampling. Furthermore, embodiments of the spherical harmonic importance sampling method can exploit the fast wavelet product to multiply two functions on-the-fly, while the conventional method requires a fixed, premultiplied environment map. 
     Basis Representations—Spherical Harmonics and Haar Wavelets 
     Embodiments of the spherical harmonic importance sampling method may rely on both spherical harmonics and hierarchical sample warping for Haar wavelets. Therefore, a brief introduction to these two basis representations is provided in this section. In subsequent sections, an embodiment of the spherical harmonic importance sampling method that leverages both spherical harmonics and hierarchical sample warping is described, and the description thereof builds upon and refers to portions of the discussion in this section. 
     Spherical Harmonics 
     Spherical harmonics are a frequency-space basis for representing functions defined over the sphere. If a direction vector {right arrow over (ω)} is represented using the standard spherical parameterization,
 
{right arrow over (ω)}=( x,y,x )=(sin θ cos φ, sin θ sin φ, cos θ)
 
     then the real-valued spherical harmonic basis functions are defined as: 
                       y   l   m     ⁡     (     θ   ,   ϕ     )       =     {             K   l   m     ⁢       P   l        m          ⁡     (     cos   ⁢           ⁢   θ     )       ⁢     cos   ⁡     (          m        ⁢   ϕ     )                   if   ⁢           ⁢   m     ≥   0     ,                 K   l   m     ⁢       P   l        m          ⁡     (     cos   ⁢           ⁢   θ     )       ⁢     sin   ⁡     (          m        ⁢   ϕ     )                   if   ⁢           ⁢   m     &lt;   0     ,                     (   1   )               
where P l   m  are the associated Legendre polynomials and K l   m  are normalization constants. In one embodiment, the Condon-Shortley phase factor of (−1) m , which is sometimes included in the definition of P l   m  or y l   m , is omitted, since this simplifies the notation. Some sources include this factor, while others do not, so care should be taken in determining the convention used when referencing the literature.
 
     The basis functions are indexed according to two integer constants, the order or band, l, and the degree, m. These satisfy the constraints that l is a nonnegative integer, lεZ*, and −l≦m≦l. Therefore, 2l+1 basis functions exist at order l and n 2  total basis functions occupy the first n bands, 1=0 . . . (n−1). The basis functions can be visually arranged into a pyramid of increasing order, as shown in  FIG. 4 .  FIG. 4  illustrates visualizations of the first nine spherical harmonic basis functions in Cartesian and spherical coordinates, where the lighter shading indicates positive values and the darker shading indicates negative values. It may be convenient to reformulate the indexing scheme to use a single parameter:
 
 i=l ( l+ 1)+ m,  
 
which flattens the pyramid using a level-order traversal.
 
Orthogonality and Normalization.
 
     The spherical harmonics are orthonormal, satisfying the condition:
 
   y   i   y   i   =∫ Ω     4π     y   i ({right arrow over (ω)}) y   j ({right arrow over (ω)}) d{right arrow over (ω)}=δ   i,j   (2)
 
where δ ij  denotes the Kronecker delta function. The normalization constant that satisfies this condition is:
 
     
       
         
           
             
               
                 
                   
                     K 
                     l 
                     m 
                   
                   = 
                   
                     
                       
                         
                           
                             
                               ( 
                               
                                 2 
                                 - 
                                 
                                   δ 
                                   
                                     0 
                                     ⁢ 
                                     m 
                                   
                                 
                               
                               ) 
                             
                             ⁢ 
                             
                               ( 
                               
                                 
                                   2 
                                   ⁢ 
                                   l 
                                 
                                 + 
                                 1 
                               
                               ) 
                             
                           
                           
                             4 
                             ⁢ 
                             π 
                           
                         
                         ⁢ 
                         
                           
                             
                               ( 
                               
                                 l 
                                 - 
                                 
                                    
                                   m 
                                    
                                 
                               
                               ) 
                             
                             ! 
                           
                           
                             
                               ( 
                               
                                 l 
                                 + 
                                 
                                    
                                   m 
                                    
                                 
                               
                               ) 
                             
                             ! 
                           
                         
                       
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     Note that other normalizations are possible and may be used, but orthonormality is used because of its convenience in this context. 
     Projection and Reconstruction. 
     The spherical harmonic functions define a complete basis over the sphere. Thus, any real-valued spherical square-integrable function ƒ may be represented exactly as a linear combination of the infinite set of basis functions. By limiting the number of bands to a finite l=n−1, an nth order, band-limited approximation {tilde over (ƒ)} of the original function ƒ may be obtained as: 
                       f   ⁡     (     ω   _     )       ≈       f   ~     ⁢     (     ω   _     )         =       ∑     l   =   0       n   -   1       ⁢       ∑     m   =     -   l       l     ⁢         y   l   m     ⁡     (     ω   _     )       ⁢       f   l   m     .                   (   4   )               
where the coefficients ƒ l   m  are computed by projecting ƒ onto each basis function y l   m . Low-frequency functions can be well approximated using only a few bands, and as the number of coefficients increases, higher frequency signals can be approximated more accurately. Because the spherical harmonic basis is orthonormal, the projection coefficients that minimizes least-squares error can be computed by simply integrating the original function against each basis function:
 
ƒ l   m =∫ Ω     4π     y   l   m ({right arrow over (ω)})ƒ({right arrow over (ω)}) d{right arrow over (ω)}   (5)
 
Other Properties and Operations of Spherical Harmonics
 
     Spherical harmonics have other useful properties such as, for instance, rotation invariance, which guarantees alias-free projections. These properties lead to a number of efficient operations for computing rotations, convolutions, and double product integrals, all in the spherical harmonic domain. These properties of spherical harmonics are summarized in the table illustrated in  FIG. 2 . 
     Limitations of Spherical Harmonics 
     Spherical harmonics are a natural basis for many problems in computer graphics due to the prevalence of graphics problems posed in the spherical or hemispherical domain. However, conventionally, spherical harmonics do not support an effective importance sampling procedure. 
     Hierarchical Sample Warping for Haar wavelets 
     In 2D, both mip maps and Haar wavelets can be seen as a quad-tree hierarchy of coefficients. In a mip map, each tree node (past the root) points to four children and contains exactly four scaling coefficients, while in a Haar wavelet each node points to four children but contains at most three wavelet coefficients. If either type of tree is traversed in top-down order, at each node it is possible to extract the average value within the four child regions. This observation leads to an efficient hierarchical sample warping scheme that warps an initially uniform set of points by first computing the marginal probability of intensity in the top and bottom half of the node. The input point set is split according to this ratio and scaled to fit within the two halves. This process is repeated independently on each of the two rows, and the whole process repeats hierarchically on the four children. It is easy to generalize this procedure to mip maps or wavelets in any dimension. This process is illustrated in  FIGS. 5 and 6 . For simplicity, the mip map basis will be used in subsequent discussion, but a Haar basis could easily be used instead. 
       FIG. 5  illustrates a hierarchical sample warping method according to one embodiment. To perform hierarchical sample warping, the method starts with a uniform set of points (a) and integrates the spherical harmonic basis functions (b) to obtain the quadrant probabilities (c). The point set is first partitioned vertically according to the marginal row probabilities (d), and then each row is divided horizontally according to the probabilities of its cells (e). The whole process repeats on the four quadrants recursively. 
       FIG. 6  illustrates a transition between (b) and (c) of  FIG. 5  according to one embodiment. The spherical harmonic basis functions (b) are converted into mip map or Haar wavelet representation, as described elsewhere herein (see, e.g., equations 6 and 10) to obtain the quadrant probabilities (c). 
     Spherical Harmonic Importance Sampling Method 
     An insight behind the spherical harmonic importance sampling method as described herein is that hierarchical sample warping can be applied to spherical harmonics as long as there is a suitable mapping onto the sphere and as long as the average value of the spherical harmonic function can be computed at each level in the traversal. The second observation, which makes this idea practical, recognizes that since any spherical harmonic function {tilde over (ƒ)} is a weighted sum of basis functions y l   m , the integral of {tilde over (ƒ)} over any region Ω R  is simply the weighted sum of basis function integrals: 
     
       
         
           
             
               
                 
                   
                     
                       ∫ 
                       
                         Ω 
                         R 
                       
                     
                     ⁢ 
                     
                       
                         
                           f 
                           ~ 
                         
                         ⁡ 
                         
                           ( 
                           
                             ω 
                             → 
                           
                           ) 
                         
                       
                       ⁢ 
                       d 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         ω 
                         → 
                       
                     
                   
                   = 
                   
                     
                       
                         ∫ 
                         
                           Ω 
                           R 
                         
                       
                       ⁢ 
                       
                         
                           ( 
                           
                             
                               ∑ 
                               
                                 i 
                                 = 
                                 0 
                               
                               
                                 n 
                                 - 
                                 1 
                               
                             
                             ⁢ 
                             
                               
                                 ∑ 
                                 
                                   m 
                                   = 
                                   
                                     - 
                                     1 
                                   
                                 
                                 l 
                               
                               ⁢ 
                               
                                 
                                   f 
                                   l 
                                   m 
                                 
                                 ⁢ 
                                 
                                   
                                     y 
                                     l 
                                     m 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       ω 
                                       → 
                                     
                                     ) 
                                   
                                 
                               
                             
                           
                           ) 
                         
                         ⁢ 
                         d 
                         ⁢ 
                         
                           ω 
                           → 
                         
                       
                     
                     = 
                     
                       
                         ∑ 
                         
                           i 
                           = 
                           0 
                         
                         
                           n 
                           - 
                           1 
                         
                       
                       ⁢ 
                       
                         
                           ∑ 
                           
                             m 
                             = 
                             
                               - 
                               1 
                             
                           
                           l 
                         
                         ⁢ 
                         
                           
                             f 
                             l 
                             m 
                           
                           ⁢ 
                           
                             
                               
                                 ∫ 
                                 
                                   Ω 
                                   R 
                                 
                               
                               ⁢ 
                               
                                 
                                   
                                     y 
                                     l 
                                     m 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       ω 
                                       → 
                                     
                                     ) 
                                   
                                 
                                 ⁢ 
                                 d 
                                 ⁢ 
                                 
                                   ω 
                                   → 
                                 
                               
                             
                             
                               ︸ 
                               
                                 
                                   y 
                                   ^ 
                                 
                                 l 
                                 m 
                               
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   6 
                   ) 
                 
               
             
           
         
       
     
     Equation 6 says that if there is access to a collection of preintegrated basis functions ŷ l   m , then hierarchical sample warping may be applied to any arbitrary spherical harmonic function by performing a weighted sum with its projection coefficients. A naive application of this concept could be to create a 2D mip map or Haar wavelet representation of each basis function in a pre-process and then perform a weighted sum of these hierarchical representations on-the-fly to perform warping. This naive approach works, but embodiments may do significantly better by exploiting the structure of the basis functions. 
     Computing Basis Function Integrals 
     If a standard spherical mapping is used, the average value over a mip map node is proportional to the following integral in spherical coordinates:
 
 ŷ   l   m   =∫∫y   l   m (θ,φ)sin θ dθdφ,   (7)
 
     This expression can be simplified by observing that the spherical harmonic basis functions are separable. This can be expressed as:
 
 y   l   m (θ,φ)= K   l   m   P   l   |m| (cos θ)Φ m (φ),  (8)
 
where:
 
                       Φ   m     ⁡     (   ϕ   )       =     {           cos   ⁡     (     m   ⁢           ⁢   ϕ     )               if   ⁢           ⁢   m     ≥   0               sin   ⁡     (          m        ⁢   ϕ     )               if   ⁢           ⁢   m     &lt;   0                     (   9   )               
and is illustrated in  FIG. 7 .  FIG. 7  illustrates the separation of spherical harmonic basis functions. In  FIG. 7 , the spherical harmonics (top) are separable into a θ-dependent component (middle) and a φ-dependent component (bottom). Therefore, it is possible to transform the two-dimensional integral in equation (7) into two one-dimensional integrals:
 
 ŷ   l   m   =∫∫K   l   m   P   l   |m| (cos θ)Φ m (φ)sin θ dθdφ=K   l   m (∫ P   l   |m| (cos θ)sin θ d θ)(∫Φ m (φ) d φ)  (10)
 
     The first term is simply the normalization constant defined in equation 3. The following subsections describe methods for computing the two integrals according to some embodiments. 
     Integrals of Associated Legendre Functions 
     The middle term in equation 10 involves computing the integral of an expression containing two trigonometric functions and an associated Legendre function. Solving this integral directly is difficult. However, the problem can be simplified by performing a change of variables to cylindrical coordinates and observing that,
 
∫ θ     −     θ     +     P   l   |m| (cos θ)sin θ dθ=∫   z     −     z     +     P   l   |m| ( z ) dz={circumflex over (P)}   l   |m| ,  (11)
 
where:
 
     z=cos θ 
     z − =cos(θ − ), 
     z + =cos(θ + ). 
     Therefore, a method is needed to compute the integral of the associated Legendre functions, {circumflex over (P)} l   |m| . 
     Since associated Legendre functions are just polynomials in z, the functions could be symbolically integrated using a computer algebra system, such as the Maple™ computer algebra system or the MATLAB® computer algebra system, and the functions could thus be hard-coded up to a sufficiently large order. Maple™ is a trademark of Waterloo Maple, Inc. MATLAB® is a trademark of The MathWorks, Inc. A similar approach has been taken for integrating the spherical harmonics over a spherical triangle. 
     Though this approach may be exact, and may be used in some embodiments, it has drawbacks. Firstly, a set of hard-coded functions requires a predetermined maximum number of bands in the spherical harmonic approximation. This may be sufficient for some applications, but reduces the generality of the approach. More importantly, in embodiments, when evaluating the integral of a spherical harmonic function using equation 6, it may be necessary to compute the integrals of all basis functions up to some band n, which would be very inefficient with this approach. 
     A more efficient approach that may be used in some embodiments is to derive recurrence relations that compute the integral {circumflex over (P)} l   |m|  directly from the integrals at lower values of l and m. Two such methods have been independently developed in the geodesy field. One of these approaches is more efficient and compact, and thus some embodiments may use this approach, although the other approach may also be used. The approach used involves two recurrence relations. The first is a recurrence on l: 
                         P   ^     l   m     =               (     l   -   2     )     ⁢     (     l   -   1   +   m     )           (     l   +   1     )     ⁢     (     l   -   m     )         ⁢       P   ^       l   -   2     m       -           2   ⁢   l     -   1         (     l   +   1     )     ⁢     (     l   -   m     )         ⁢     (     1   -     z   2       )     ⁢       P     l   -   1     m     ⁡     (   z   )           ⁢     |     z   -       z   +           ,           (   12   )               
which is valid for 0≦m&lt;l. When m=l−1, the first term above drops out since {circumflex over (P)} l-2   l-1  is implicitly zero. This relation has a singularity when m=1, so a second recurrence may be used for this case:
 
     
       
         
           
             
               
                 
                   
                     
                       P 
                       ^ 
                     
                     l 
                     l 
                   
                   = 
                   
                     
                       
                         1 
                         
                           l 
                           + 
                           1 
                         
                       
                       ⁡ 
                       
                         [ 
                         
                           
                             
                               
                                 l 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     
                                       2 
                                       ⁢ 
                                       l 
                                     
                                     - 
                                     3 
                                   
                                   ) 
                                 
                               
                               ⁢ 
                               
                                 ( 
                                 
                                   
                                     2 
                                     ⁢ 
                                     l 
                                   
                                   - 
                                   1 
                                 
                                 ) 
                               
                               ⁢ 
                               
                                 
                                   P 
                                   ^ 
                                 
                                 
                                   l 
                                   - 
                                   2 
                                 
                                 
                                   l 
                                   - 
                                   2 
                                 
                               
                             
                             + 
                             
                               z 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 
                                   P 
                                   l 
                                   l 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   z 
                                   ) 
                                 
                               
                             
                           
                           ⁢ 
                           
                             | 
                             
                               z 
                               - 
                             
                             
                               z 
                               + 
                             
                           
                         
                         ] 
                       
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
     To start the recurrence, the following may be used in some embodiments: 
     
       
         
           
             
               
                 
                   
                     
                       
                         P 
                         ^ 
                       
                       0 
                       0 
                     
                     = 
                     
                       
                         z 
                         + 
                       
                       - 
                       
                         z 
                         - 
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                         P 
                         ^ 
                       
                       1 
                       0 
                     
                     = 
                     
                       
                         
                           z 
                           + 
                           2 
                         
                         - 
                         
                           z 
                           - 
                           2 
                         
                       
                       2 
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       P 
                       ^ 
                     
                     1 
                     1 
                   
                   = 
                   
                     
                       
                         1 
                         2 
                       
                       ⁡ 
                       
                         [ 
                         
                           
                             z 
                             ⁢ 
                             
                               
                                 1 
                                 - 
                                 
                                   z 
                                   2 
                                 
                               
                             
                           
                           + 
                           
                             
                               sin 
                               
                                 - 
                                 1 
                               
                             
                             ⁡ 
                             
                               ( 
                               z 
                               ) 
                             
                           
                         
                         ] 
                       
                     
                     ⁢ 
                     
                       | 
                       
                         z 
                         - 
                       
                       
                         z 
                         + 
                       
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   16 
                   ) 
                 
               
             
           
         
       
     
     These two recurrence relations are numerically robust and allow embodiments to compute all the integrals up to order n in O(n 2 ) time (note that this is actually linear time since there are n 2  basis functions up to order n). In contrast, symbolic evaluation would require O(n 4 ) time and can suffer from catastrophic floating-point cancellation when implemented naïvely. An example of pseudo-code for this recurrence procedure is provided in Algorithm 1, shown in  FIG. 8 . 
     Integral Along Polar Direction 
     The last term in equation 10 is trivial to compute since it is the integral of a trigonometric function: 
     
       
         
           
             
               
                 
                   
                     
                       
                         Φ 
                         ^ 
                       
                       m 
                     
                     ⁡ 
                     
                       ( 
                       ϕ 
                       ) 
                     
                   
                   = 
                   
                     
                       ∫ 
                       
                         
                           
                             Φ 
                             m 
                           
                           ⁡ 
                           
                             ( 
                             ϕ 
                             ) 
                           
                         
                         ⁢ 
                         d 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         ϕ 
                       
                     
                     = 
                     
                       
                         1 
                         m 
                       
                       ⁢ 
                       
                         { 
                         
                           
                             
                               
                                 sin 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     m 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     ϕ 
                                   
                                   ) 
                                 
                               
                             
                             
                               
                                 
                                   
                                     if 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     m 
                                   
                                   ≥ 
                                   0 
                                 
                                 , 
                               
                             
                           
                           
                             
                               
                                 cos 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     m 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     ϕ 
                                   
                                   ) 
                                 
                               
                             
                             
                               
                                 
                                   if 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   m 
                                 
                                 &lt; 
                                 0. 
                               
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   17 
                   ) 
                 
               
             
           
         
       
     
     This can be evaluated directly; however, to compute integrals for all basis functions up to some band in, this requires O(n) calls to expensive transcendental functions. This may be made more efficient by using the following multiple-angle recurrence relations:
 
sin( m φ)=2 sin(( m− 1)φ)cos(φ)−sin(( m− 2)φ)  (18)
 
cos( m φ)=2 cos(( m− 1)φ)cos(φ)−cos(( m− 2)φ).  (19)
 
     This recurrence can be initialized with: 
     m=0: cos(0)=1, sin(0)=0, 
     m=1: cos(φ)=cos(φ), sin(φ)=sin(φ), 
     which requires only a single evaluation of the two trigonometric functions for any value of n. An example of pseudo-code for this recurrence method is provided in Algorithm 2, shown in  FIG. 9 . 
     Spherical Harmonic Importance Sampling 
     In embodiments, the concepts and results from the previous section may be applied into a practical approach for importance sampling spherical harmonics. In embodiments, standard cylindrical coordinates may be used to map the 2D mip map or Haar wavelets onto the spherical domain. Spherical coordinates may also be used in embodiments, but cylindrical coordinates are area-preserving, which simplifies the procedure. Given this mapping, hierarchical warping may need to evaluate equation 6 over the area of each node in the traversal. An embodiment may therefore need to compute the integral, ŷ l   m , of the spherical harmonic basis functions, which, in cylindrical coordinates, becomes:
 
 ŷ   l   m   =K   l   m (∫ z−   z+   P   l   |m| ( z ) dz )(∫ φ−   φ+ Φ m (φ) d φ).  (20)
 
     These two integrals can be efficiently evaluated for all values of l and m, for example using example algorithm 1 and example algorithm 2 illustrated in  FIG. 8  and  FIG. 9 , respectively. 
     To perform warping, an embodiment may start at the root and evaluate the definite integrals of all basis functions, ŷ l   m , over the four quadrants of the sphere. The weighted sum of these integrals according to equation 6 gives the exact integral of the spherical harmonic function for the four quadrants. An embodiment may warp the points according to this 2×2 importance function, and repeat hierarchically for each of the four child nodes. At each stage of warping, an embodiment may evaluate equation 20 with integration bounds defined by the extents of the four child nodes. At the end of the procedure the points may be distributed with a probability density function (PDF) proportional to the spherical harmonic function. In one embodiment, the PDF of each sample is simply the integral of the containing node divided by the integral of the whole function. 
     Probability Density Function (PDF) 
     In embodiments, a probability density function (PDF) may be used for sampling the integrand. Rather than sampling the integrand uniformly in an unbiased manner, embodiments may purposely sample the integrand more in regions that will have more impact on the result of the integral. To remove bias, embodiments may re-weight the samples that have been taken, since more samples are placed in some areas than in others. 
     Suppose the following 1D integral is to be sampled: 
     
       
         
           
             F 
             = 
             
               
                 ∫ 
                 a 
                 b 
               
               ⁢ 
               
                 
                   f 
                   ⁡ 
                   
                     ( 
                     x 
                     ) 
                   
                 
                 ⁢ 
                 dx 
               
             
           
         
       
     
     n samples can be uniformly generated to get an approximate value: 
     
       
         
           
             
               F 
               approx 
             
             = 
             
               
                 
                   b 
                   - 
                   c 
                 
                 n 
               
               * 
               
                 
                   ∑ 
                   1 
                   n 
                 
                 ⁢ 
                 
                   f 
                   ⁡ 
                   
                     ( 
                     x 
                     ) 
                   
                 
               
             
           
         
       
     
     If there is a probability density function (PDF) and samples are drawn non-uniformly from some distribution, then: 
     
       
         
           
             
               F 
               approx 
             
             = 
             
               
                 1 
                 n 
               
               * 
               
                 
                   ∑ 
                   1 
                   n 
                 
                 ⁢ 
                 
                   ( 
                   
                     
                       f 
                       ⁡ 
                       
                         ( 
                         x 
                         ) 
                       
                     
                     
                       p 
                       ⁡ 
                       
                         ( 
                         x 
                         ) 
                       
                     
                   
                   ) 
                 
               
             
           
         
       
     
     The PDF here is given by p(x). Non-uniform sampling may yield faster convergence when a good PDF is found. 
     Using 1D Mip Maps 
     Embodiments of the warping technique described above require no precomputation and may have infinite precision, since an embodiment can continue warping to an arbitrary depth. However, infinite precision may not be necessary, especially for low-order spherical harmonic functions that vary smoothly. To improve performance, a cylindrical mip map representation of all the needed basis functions may be computed, and then warping may be performed on the weighted sum of these mip maps. However, equations 10 and 20 show that, due to separability, the 2D integral of the spherical harmonic bases may be computed by only performing two 1D integrals. Thus, embodiments may implicitly construct each of these 2D mip maps by storing only two 1D mip maps. This significantly reduces memory requirements and also improves performance due to better cache coherence. This leads to embodiments that implement an efficient importance sampling scheme for spherical harmonics which only requires 2n 2  1D mip maps for order n spherical harmonic functions. Note that these mip maps need to be created only once, since they are not function-dependent and can be used for any spherical harmonic function by weighting using a different set of coefficients ƒ l   m . 
     Product Sampling 
     There are techniques to directly importance sample the product of two functions stored as Haar wavelets by exploiting a fast wavelet product. Embodiments of the spherical harmonic sampling method described herein may implicitly create a mip map or Haar wavelet representation of the function. Embodiments can therefore be used to exploit the fast wavelet product to importance sample the product of a spherical harmonic function and a mip map or Haar wavelet on-the-fly. Thus, embodiments may bridge the gap between the spherical harmonic and Haar wavelet bases, allowing embodiments to use spherical harmonics for smooth functions such as BRDFs and Haar wavelets for high-frequency environment maps. Furthermore, the spherical harmonic function can be efficiently rotated in embodiments, thus allowing embodiments to use a local coordinate system and a single, high-resolution environment map. 
     Flowcharts 
       FIGS. 10A and 10B  are flowcharts of methods in importance sampling spherical harmonics according to embodiments.  FIG. 10A  is a flowchart of a method of performing rotation and importance sampling according to one embodiment. As indicated at  800 , a reflectance function may be stored as a spherical harmonic representation, for example as a BRDF. As indicated at  802 , the spherical harmonic representation of the reflectance function may be rotated to align with an environment map. The environment map may be stored as a Haar wavelet representation. As indicated at  804 , the spherical harmonic representation of the reflectance function may be converted to a Haar wavelet representation using an efficient algorithmic conversion method as described elsewhere herein. As indicated at  806 , the Haar wavelet representation of the reflectance function can then be importance sampled. 
       FIG. 10B  is a flowchart of a method of importance sampling the product of two functions according to one embodiment. As indicated at  820 , an environment map may be stored as a Haar wavelet representation. As indicated at  822 , a corresponding reflectance function may be stored as a spherical harmonic representation, for example as a BRDF. As indicated at  824 , the spherical harmonic representation of the reflectance function may be rotated to align with an environment map. As indicated at  826 , the spherical harmonic representation of the reflectance function may be converted to a Haar wavelet representation using an efficient algorithmic conversion method as described elsewhere herein. As indicated at  828 , the product of the Haar wavelet representation of the reflectance function and the Haar wavelet representation of the environment map may be generated, and the product may be importance sampled. 
     Example Applications and Implementations of Embodiments 
     An example of an application of embodiments is in computer graphics, for example in 3-D rendering of synthetic scenes, or photorealistic rendering. In this application, there is generally an environment map, and a set of one or more reflectance functions defined on objects in a scene. A 3-D rendering process is used to compute a high-quality, illuminated result or results from the environment map and reflectance function(s). Embodiments may enable a 3-D rendering process to more quickly and efficiently generate quality results, i.e. output images, than could be done using conventional techniques. Embodiments may be applied to other areas or problems in computer graphics, such as artificial stippling, and procedural geometry placement, and may also be extended into areas outside of computer graphics. 
     Embodiments may be implemented as a stand-alone program, a module or plug-in, a library function, or generally as any sort of software entity. For simplicity, an implementation of an embodiment of the method described herein may be referred to as a function representation conversion module. 
       FIG. 11  illustrates a function representation conversion module performing importance sampling in a rendering system according to one embodiment. The function representation conversion module  910  may be implemented in, called by, or plugged into a rendering system  900  that traces rays via a ray tracing function and performs other rendering tasks to render an output synthetic image or images from environment map(s) and reflectance function(s) of surfaces.  FIG. 12  illustrates an example of a system in which a rendering system  900  may be implemented, e.g. as or in graphic application  1220 . Returning again to  FIG. 11 , in an application for 3-D rendering of synthetic scenes, function representation conversion module  910  may receive, as input, lighting (e.g., environment map  902 ) and reflectance (e.g., reflectance function  904 , for example a BRDF). Environment map  902  may, for example, be stored as a Haar wavelet representation, and reflectance  904  may be stored as a spherical harmonic representation. Function representation conversion module  910  may generate samples  906  according to the methods described herein. Samples  906 , in this example, represent directions over a sphere, hemisphere, or other shape that indicate where rays are to be traced. The output from the function representation conversion module  910  may be provided to a rendering module  920  that performs the ray tracing, and so on, for use in rendering an output image or images  930 . Output image(s)  930  may be stored to a storage device, transmitted over a network, printed, or otherwise used in a computer graphics application. The rendering system  900  may use the function representation conversion module  910  to perform rotations, to perform importance sampling of a reflectance function stored as a spherical harmonic representation, and to perform importance sampling of the product of the reflectance function and the environment map. In so doing, function representation conversion module  910  may convert a function from one representation (spherical harmonics) into another representation (e.g., Haar wavelets or mip map). Rotations may be performed using a spherical harmonic representation, and importance sampling may be performed using a Haar wavelet representation. 
     Example System 
     The methods described herein including those for converting spherical harmonics representations of functions to multi-resolution representations may be implemented, for example as a function representation conversion module  910  as illustrated in  FIG. 11 , by a computer system configured to provide the functionality described.  FIG. 12  is a block diagram illustrating one embodiment of a computer system  1200  suitable for implementing embodiments of a function representation conversion module  910 . A graphics application such as graphics application  1220  may be configured to perform various image processing functions and to render new images accordingly. In some embodiments, a user may invoke operations to perform various functions as described for embodiments via a user interface of graphics application  1220 . Graphics application  1220  may be configured to perform these operations using the methods described herein for importance sampling spherical harmonics. Graphics application  1220  may, for example, be configured to render an output synthetic image using embodiments of a function representation conversion module  910  to generate samples  906 . 
     Graphics application  1220  may represent various types of graphics applications, such as painting, publishing, photography, games, animation, and other applications. Additionally, graphics application  1220  may utilize a graphics processor  1240  when rendering or displaying images according to various embodiments. A graphics processing unit or GPU may be considered a dedicated graphics-rendering device for a personal computer, workstation, game console or other computer system. Modern GPUs may be very efficient at manipulating and displaying computer graphics and their highly parallel structure may make them more effective than typical CPUs for a range of complex graphical algorithms. For example, graphics processor  1240  may implement a number of graphics primitive operations in a way that makes executing them much faster than drawing directly to the screen with a host central processing unit (CPU), such as CPU  1230 . In various embodiments, the methods disclosed herein for may be implemented by program instructions configured for parallel execution on two or more such GPUs. The GPU  1200  may implement one or more application programmer interfaces (APIs) that permit programmers to invoke the functionality of the GPU. Suitable GPUs may be commercially available from vendors such as NVIDIA Corporation, ATI Technologies, and others. 
     Please note that functionality and/or features described herein as being part of, or performed by, graphics application  1220  may, in some embodiments, be part of, or performed by, one or more graphics processors, such as graphics processor  1240 . As described above, in some embodiments graphics application  1220  may be configured to render modified and/or reconstructed images into a different window than input images. 
     The various functions described herein for embodiments of a function representation conversion module may be implemented on various types of computer systems. Referring again to  FIG. 12 , computer system  1200  may be any of various types of devices, including, but not limited to, a personal computer system, desktop computer, laptop or notebook computer, mainframe computer system, handheld computer, workstation, network computer, a consumer device, video game console, handheld video game device, application server, storage device, a peripheral device such as a switch, modem, router, or in general any type of computing device. 
     Graphics application  1220 , which may be configured to implement a function representation conversion module  910  as described herein, may be provided as a computer program product, or software, that may include a computer-readable storage medium having stored thereon instructions, which may be used to program a computer system (or other electronic devices) to implement computer graphics functions such as rendering synthetic images using the techniques described herein. A computer-readable storage medium may include any mechanism for storing information in a form (e.g., software, processing application) readable by a machine (e.g., a computer). The machine-readable storage medium may include, but is not limited to, magnetic storage medium (e.g., floppy diskette); optical storage medium (e.g., CD-ROM); magneto optical storage medium; read only memory (ROM); random access memory (RAM); erasable programmable memory (e.g., EPROM and EEPROM); flash memory; electrical, or other types of medium suitable for storing program instructions. In addition, program instructions may be communicated using optical, acoustical or other form of propagated signal (e.g., carrier waves, infrared signals, digital signals, or other types of signals or mediums.). 
     As illustrated in  FIG. 12 , computer system  1200  may include one or more processor units (CPUs)  1230 . Processors  1230  may be implemented using any desired architecture or chip set, such as the SPARC™ architecture, an x86-compatible architecture from Intel Corporation or Advanced Micro Devices, or another architecture or chipset capable of processing data, and may in various embodiments include multiple processors, a single threaded processor, a multi-threaded processor, a multi-core processor, or any other type of general-purpose or special-purpose processor. Any desired operating system(s) may be run on computer system  1200 , such as various versions of Unix, Linux, Windows™ from Microsoft Corporation, MacOS™ from Apple Corporation, or any other operating system that enables the operation of software on a hardware platform. 
     The computer system  1200  may also include one or more system memories  1210  (e.g., one or more of cache, SRAM, DRAM, RDRAM, EDO RAM, DDR RAM, SDRAM, Rambus RAM, EEPROM, or other memory type), or other types of RAM or ROM) coupled to other components of computer system  1200  via interconnect  1260 . Memory  1210  may include other types of memory as well, or combinations thereof. One or more of memories  1210  may include program instructions  1215  executable by one or more of processors  1230  to implement a function representation conversion module  910  as described herein. Program instructions  1215 , which may include program instructions configured to implement graphics application  1220 , may be partly or fully resident within the memory  1210  of computer system  1200  at any point in time. Alternatively, program instructions  1215  may be provided to GPU  1240  for performing various computer graphics operations (or portions thereof) on GPU  1240  using one or more of the techniques described herein. In some embodiments, the techniques described herein may be implemented by a combination of program instructions  1215  executed on one or more processors  1230  and one or more GPUs  1240 , respectively. Program instructions  1215  may also be stored on an external storage device (not shown) accessible by the processor(s)  1230  and/or GPU  1240 , in some embodiments. Any of a variety of such storage devices may be used to store the program instructions  1215  in different embodiments, including any desired type of persistent and/or volatile storage devices, such as individual disks, disk arrays, optical devices (e.g., CD-ROMs, CD-RW drives, DVD-ROMs, DVD-RW drives), flash memory devices, various types of RAM, holographic storage, etc. The storage devices may be coupled to the processor(s)  1230  and/or GPU  1240  through one or more storage or I/O interfaces including, but not limited to, interconnect  1260  or network interface  1250 , as described herein. In some embodiments, the program instructions  1215  may be provided to the computer system  1200  via any suitable computer-readable storage medium including memory  1210  and/or external storage devices described above. Memory  1210  may also be configured to implement one or more data structures  1225 , such as one or more data structures for rendering synthetic images using a function representation conversion module  910  to generate samples  906 . Data structures  1225  may be accessible by processor(s)  1230  and/or GPU  1240  when executing graphics application  1220  or other program instructions  1215 . 
     As shown in  FIG. 12 , processor(s)  1230  may be coupled to one or more of the other illustrated components by at least one communications bus, such as interconnect  1260  (e.g., a system bus, LDT, PCI, ISA, or other communication bus type), and a network interface  1250  (e.g., an ATM interface, an Ethernet interface, a Frame Relay interface, or other interface). The CPU  1230 , the network interface  1250 , and the memory  1210  may be coupled to the interconnect  1260 . It should also be noted that one or more components of system  1200  might be located remotely and accessed via a network. 
     As noted above, in some embodiments, memory  1210  may include program instructions  1215 , comprising program instructions configured to implement graphics application  1220 , as described herein. Graphics application  1220  may be implemented in various embodiments using any desired programming language, scripting language, or combination of programming languages and/or scripting languages, e.g., C, C++, C#, Java™, Perl, etc. For example, in one embodiment, graphics application  1220  may be JAVA based, while in another embodiments, it may be implemented using the C or C++ programming languages. In other embodiments, graphics application  1220  may be implemented using specific graphic languages specifically for developing programs executed by specialize graphics hardware, such as GPU  1240 . In addition, graphics application  1220  may be embodied on memory specifically allocated for use by graphics processor(s)  1240 , such as memory on a graphics board including graphics processor(s)  1240 . Thus, memory  1210  may represent dedicated graphics memory as well as general-purpose system RAM, in various embodiments. Memory  1210  may in some embodiments also include a data store configured to store image data for one or more input images or other graphics information and/or output images, or other information for use in computer graphics such as environment maps and reflectance functions, in various embodiments. Other information not described herein may be included in memory  1210  and may be used to implement the methods described herein and/or other functionality of computer system  1200 . 
     Network interface  1250  may be configured to enable computer system  1200  to communicate with other computers, systems or machines, such as across a network. Network interface  1250  may use standard communications technologies and/or protocols, and may utilize links using technologies such as Ethernet, 802.11, integrated services digital network (ISDN), digital subscriber line (DSL), and asynchronous transfer mode (ATM) as well as other communications technologies. Similarly, the networking protocols used on a network to which computer system  1200  is interconnected may include multiprotocol label switching (MPLS), the transmission control protocol/Internet protocol (TCP/IP), the User Datagram Protocol (UDP), the hypertext transport protocol (HTTP), the simple mail transfer protocol (SMTP), and the file transfer protocol (FTP), among other network protocols. The data exchanged over such a network by network interface  1250  may be represented using technologies, languages, and/or formats, such as the hypertext markup language (HTML), the extensible markup language (XML), and the simple object access protocol (SOAP) among other data representation technologies. Additionally, all or some of the links or data may be encrypted using any suitable encryption technologies, such as the secure sockets layer (SSL), Secure HTTP and/or virtual private networks (VPNs), the international data encryption standard (DES or IDEA), triple DES, Blowfish, RC2, RC4, RC5, RC6, as well as other data encryption standards and protocols. In other embodiments, custom and/or dedicated data communications, representation, and encryption technologies and/or protocols may be used instead of, or in addition to, the particular ones described above. 
     GPUs, such as GPU  1240  may be implemented in a number of different physical forms. For example, GPU  1240  may take the form of a dedicated graphics card, an integrated graphics solution and/or a hybrid solution. GPU  1240  may interface with the motherboard by means of an expansion slot such as PCI Express Graphics or Accelerated Graphics Port (AGP) and thus may be replaced or upgraded with relative ease, assuming the motherboard is capable of supporting the upgrade. However, a dedicated GPU is not necessarily removable, nor does it necessarily interface the motherboard in a standard fashion. The term “dedicated” refers to the fact that hardware graphics solution may have RAM that is dedicated for graphics use, not to whether the graphics solution is removable or replaceable. Dedicated GPUs for portable computers may be interfaced through a non-standard and often proprietary slot due to size and weight constraints. Such ports may still be considered AGP or PCI express, even if they are not physically interchangeable with their counterparts. As illustrated in  FIG. 12 , memory  1210  may represent any of various types and arrangements of memory, including general-purpose system RAM and/or dedication graphics or video memory. 
     Integrated graphics solutions, or shared graphics solutions are graphics processors that utilize a portion of a computer&#39;s system RAM rather than dedicated graphics memory. For instance, modern desktop motherboards normally include an integrated graphics solution and have expansion slots available to add a dedicated graphics card later. As a GPU may be extremely memory intensive, an integrated solution finds itself competing for the already slow system RAM with the CPU as the integrated solution has no dedicated video memory. For instance, system RAM may experience a bandwidth between 2 GB/s and 8 GB/s, while most dedicated GPUs enjoy from 15 GB/s to 30 GB/s of bandwidth. Hybrid solutions may also share memory with the system memory, but may have a smaller amount of memory on-board than discrete or dedicated graphics cards to make up for the high latency of system RAM. Data communicated between the graphics processing unit  1240  and the rest of the computer system  1200  may travel through a graphics card slot or other interface, such as interconnect  1260  of  FIG. 12 . 
     Computer system  1200  may also include one or more additional I/O interfaces, such as interfaces for one or more user input devices  1270 , or such devices may be coupled to computer system  1200  via network interface  1250 . For example, computer system  1200  may include interfaces to a keyboard, a mouse or other cursor control device, a joystick, or other user input devices  1270 , in various embodiments. Additionally, the computer system  1200  may include one or more displays (not shown), coupled to processors  1230  and/or other components via interconnect  1260  or network interface  1250 . Such input/output devices may be configured to allow a user to interact with graphics application  1220  to request or invoke various computer graphics operations and/or to specify various parameters, thresholds, and/or other configurable options available to the user when performing computer graphics tasks while executing graphic application  1220 . It will be apparent to those having ordinary skill in the art that computer system  1200  may also include numerous other elements not shown in  FIG. 12 . 
     Note that program instructions  1215  may be configured to implement a graphic application  1220  as a stand-alone application, or as a module of another graphics application or graphics library, in various embodiments. For example, in one embodiment program instructions  1215  may be configured to implement graphics applications such as painting, publishing, photography, games, animation, and/or other applications, and may be configured to render and/or edit synthetic images as part of one or more of these graphics applications. In another embodiment, program instructions  1215  may be configured to implement the techniques of various embodiments as described herein in one or more functions called by another graphics application executed on GPU  1240  and/or processor(s)  1230 . Program instructions  1215  may also be configured to render images and present them on one or more displays as the output of an computer graphics operation and/or to store image data for synthesized, modified and/or reconstructed images in memory  1210  and/or an external storage device(s), in various embodiments. For example, a graphics application  1220  included in program instructions  1215  may utilize GPU  1240  when modifying, rendering, or displaying images in some embodiments. 
     Those skilled in the art will appreciate that computer system  1200  is merely illustrative and is not intended to limit the scope of a method and apparatus for importance sampling spherical harmonics as described herein. In particular, the computer system and devices may include any combination of hardware or software that can perform the indicated functions, including computers, network devices, internet appliances, PDAs, wireless phones, pagers, etc. Computer system  1200  may also be connected to other devices that are not illustrated, or instead may operate as a stand-alone system. In addition, the functionality provided by the illustrated components may in some embodiments be combined in fewer components or distributed in additional components. Similarly, in some embodiments, the functionality of some of the illustrated components may not be provided and/or other additional functionality may be available. 
     Those skilled in the art will also appreciate that, while various items are illustrated as being stored in memory or on storage while being used, these items or portions of them may be transferred between memory and other storage devices for purposes of memory management and data integrity. Alternatively, in other embodiments some or all of the software components may execute in memory on another device and communicate with the illustrated computer system via inter-computer communication. Some or all of the system components or data structures may also be stored (e.g., as instructions or structured data) on a computer-accessible medium or a portable article to be read by an appropriate drive, various examples of which are described above. In some embodiments, instructions stored on a computer-accessible medium separate from computer system  1200  may be transmitted to computer system  1200  via transmission media or signals such as electrical, electromagnetic, or digital signals, conveyed via a communication medium such as a network and/or a wireless link. Various embodiments may further include receiving, sending or storing instructions and/or data implemented in accordance with the foregoing description upon a computer-accessible medium. Accordingly, the present invention may be practiced with other computer system configurations. 
     Various embodiments may further include receiving, sending or storing instructions and/or data implemented in accordance with the foregoing description upon a computer-accessible medium. Generally speaking, a computer-accessible medium may include storage media or memory media such as magnetic or optical media, e.g., disk or DVD/CD-ROM, volatile or non-volatile media such as RAM (e.g., SDRAM, DDR, RDRAM, SRAM, etc.), ROM, etc., as well as transmission media or signals such as electrical, electromagnetic, or digital signals, conveyed via a communication medium such as network and/or a wireless link. 
     The various methods as illustrated in the Figures and described herein represent examples of embodiments of methods. The methods may be implemented in software, hardware, or a combination thereof. The order of method may be changed, and various elements may be added, reordered, combined, omitted, modified, etc. 
     Various modifications and changes may be made as would be obvious to a person skilled in the art having the benefit of this disclosure. It is intended that the invention embrace all such modifications and changes and, accordingly, the above description to be regarded in an illustrative rather than a restrictive sense.