Patent Publication Number: US-6906714-B2

Title: Accelerated visualization of surface light fields

Description:
TECHNICAL FIELD 
   This disclosure relates to computer graphics and, more particularly, to techniques for visualizing surface light fields representing an object from multiple view points under multiple lighting locations. 
   BACKGROUND 
   Image synthesis, which may be used, for example, when creating computer graphics for video games, web pages and the like, is typically carried out using techniques based either on analytical models or on acquired images of a physical object. Not only must the synthesized images accurately represent the shape and color of the physical object, the synthesized images must accurately represent the ways in which the physical object reflects light (e.g., the radiance of the object). As described below, techniques based on analytical models, as well as those based on acquired images, have shortcomings. 
   One technique that is based on analytical models is performed using global illumination techniques, which combine light transport simulation with analytic material reflectance models to achieve simple, compact and flexible representation of three-dimensional (3D) scenes or objects. Global illumination techniques are computationally efficient once the analytical models have been developed, and, therefore, are easily rendered to create real-time graphics. However, analytical material reflectance models and, therefore, global illumination imaging based on such reflectance models are costly to compute and difficult to develop. 
   Image-based techniques have recently emerged as an alternative approach to creating costly analytical models for realistic image synthesis. Image-based techniques synthesize images based on graphical information of an object. Such graphical information may include pixel information of the object from a digital camera or from a scanner. Image-based techniques represent radiance data directly in a sample-based format without using any analytical models and, therefore, have gained popularity because they promise simple acquisition of 3D models and an accurate portrayal of physical objects and the physical world. 
   Presently, image-based techniques process data from, for example, a number of digital camera acquisitions taken from different angles while a lighting source is located in a fixed position. For example, an object may be digitally photographed from thirty or more different angles or positions, while a light source illuminating the object is held in a fixed position. As will be appreciated, the information of an object developed by the digital camera may occupy several megabytes (MB) of storage space, making it difficult or impractical for a website or a computer graphics intensive game to render the object in real time. Additionally, the fixed lighting position leads to acquired images having fixed illumination thereon. Such images do not always realistically reproduce how the object would appear under diverse lighting conditions. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
       FIG. 1  illustrates an example object illuminated from a number of different positions and digitally captured from a number of different positions. 
       FIG. 2  is a diagram illustrating example data flow through the disclosed system. 
       FIG. 3  is a block diagram illustrating example hardware on which the disclosed system may be implemented. 
       FIG. 4  is a flow diagram illustrating an example image preparation routine that may be executed by the processing unit of FIG.  3 . 
       FIG. 5  is an illustration of an example object to be acquired and rendered. 
       FIG. 6  is an illustration of the object of  FIG. 5  having an example mesh topology applied thereto. 
       FIG. 7  is a flow diagram of an example partition the BSRFs routine. 
       FIGS. 8A-8D  are illustrations of an example partitioning of the object of FIG.  6 . 
       FIG. 9  is an illustration of an example reparameterization of one vertex of the object of FIG.  6 . 
       FIG. 10  is an illustration of an example approximate the partitioned BSRFs routine. 
       FIGS. 11 and 12  collectively form a flow diagram of an example image rendering routine. 
   

   DETAILED DESCRIPTION 
   Although the following discloses an example of one system including software executed on hardware, it should be noted that such a system is merely illustrative and should not be considered to be limiting. For example, it is contemplated that the disclosed functionality could be embodied exclusively in dedicated hardware, exclusively in software, exclusively in firmware or in some combination of hardware, firmware and/or software. Accordingly, while the following describes one example of a system, the example is not the only way to implement such a system. 
   Turning now to  FIG. 1 , an object  20  that is to be processed using image-based acquisition may be illuminated at any one time by any one of a number of light sources, three of which are shown at reference numerals  22 ,  24  and  26 . While the object  20  is illuminated by one of the light sources  22 - 26 , one or more digital cameras  30 ,  32  may be used to acquire images of the object. For example, the camera  30  may acquire one image while the light source  22  is powered and may likewise acquire single images while each one of the light sources  24  and  26  are individually powered. The camera  30  may then be moved to a different position, where the camera  30  will acquire three more images, one while each of the light sources  22 - 26  is powered. As will be readily appreciated, any suitable number of lighting positions and camera positions may be used. For example, thirty camera positions may be selected and thirty lighting positions may be selected, wherein a camera acquires one image under each light source for each of the thirty light source positions and each of the thirty camera positions. As will be further appreciated, such a process will result, for example, in the acquisition of 900 images, each of which may be one or more megabytes in size. 
   While images of the entire object  20  are acquired, for purposes of discussion particular attention will be paid to an example point  40  of the object  20 . Thus, the remaining disclosure herein pertains to the point  40  and the way with which image data for the point  40  is processed. For example, as shown in  FIG. 1 , the light sources  22 - 26  illuminate the point  40  with vectors l 1 -l 3 , respectively. Additionally, the point  40  is perceived by the cameras  30  and  32  with vectors v 1  and v 2 , respectively. 
   A bidirectional surface radiance function (BSRF) is used to represent the outgoing radiance of substantially every surface point s=(u, v) to be modeled on the object  20  in every light source direction l=(θ i , φ i ) and every viewing direction v=(θ o , φ o ) The BSRF function is a six-dimensional (6D) function that may be represented as a function of three two dimensional (2D) vectors, as shown in Equation 1.
 
ƒ( u,v,φ   i ,θ i ,φ o ,θ o )=ƒ( s,l,v )  Equation 1
 
   Complex global illumination effects, such as, for example, inter-reflections, translucency and self-shadowing, are inherent in the BSRF representation and further enhance the realism of 3D scenes rendered to a viewer. As described in detail below, rendering using BSRFs comprises modulating each sample of the BSRF function by the light source intensity E(l) as shown in Equation 2. 
               I   ⁡     (     s   ,   l   ,   v     )       =       ∫   v     ⁢         f   ⁡     (     s   ,   l   ,   v     )       ·     E   ⁡     (   l   )         ⁢     ⅆ   l                 Equation   ⁢           ⁢   2             
 
   Because under the BSRF formulation, each modeled surface point can have a unique reflectance property, the BSRF representation is ideally suited for modeling the reflectance properties of physical objects and scenes scanned through 3D photography. Scanned objects are often characterized by fine surface details, called micro-geometry, that are difficult or impossible to scan. Using BSRF, the lighting effects that the micro-geometry has on the object appearance can be represented. Accordingly, the BSRF has the effect of both simplifying the acquisition process and improving the visual fidelity of the scanned objects. 
   Because BSRFs are large and have high dimensionality, it is not possible to perform direct representation and manipulation of BSRFs. Accordingly, as described in detail below, in the illustrated system, BSRFs are separably factorized, or decomposed, and light field mapped. Techniques of separable factorization are described in detail in Kautz, J., “Interactive Rendering with Arbitrary BRDFs using Separable Approximations,” Eurographics Rendering Workshop 1999 (June 1999) and light field mapping techniques are addressed in Chen, “Light Field Mapping: Efficient Representation and Hardware Rendering of Surface Light Fields,” To Appear in the Proceedings of SIGGRAPH 2002 (July 2002). 
   More specifically, the disclosed system approximates the BSRF function ƒ as a sum of lower dimensional products, called maps. In particular, as shown in Equation 3, variables g, α, and β represent surface, light and view maps, respectively. 
               f   ⁡     (     s   ,   l   ,   v     )       ≅       ∑     i   =   1     k     ⁢           ⁢         g   i     ⁡     (   s   )       ⁢       α   i     ⁡     (     ω   1     )       ⁢       β   i     ⁡     (     ω   2     )                   Equation   ⁢           ⁢   3             
 
In Equation 3, ω 1  and ω 2  are the reparameterized vectors l and v, respectively, which as described below, are reparameterized based on the portioning used for breaking down the BSRF. In Equation 3, the variable k represents a small integer that is the number of sums, or approximation terms, used to represent the BSRF. In general, BSRF approximations are compact, accurate approximations that are formed by partitioning the BSRF data across small surface primitives and building the approximations for each primitive independently. The BSRF partitioning is carried out in such a way as to ensure substantially continuous approximations across neighboring surface elements. As described below, a reparameterization scheme may be used to improve further the separability of BSRF data.
 
   In general, as shown in  FIG. 2 , raw BSRF data is provided to a decomposer  50 , which, as described in detail below, decomposes the BSRF data into a sum of products and passes the result of the decomposition to a quantizer  52 . The decomposition of the BSRF data may alone lead to a 100:1 data compression. The vector quantizer  52 , which may be any standard run-time memory reduction technique, may further compress the data by a factor of 10:1. Further, hardware texture compression techniques (block  54 ) may also be used to provide additional compression on the order of 6:1. As will be appreciated, significant compression, which may be on the order of 6000:1 or higher, may be achieved. This compression is significant given the large amount of data required to represent a complete BSRF without decomposition. 
   After hardware texture compression is performed at block  54 , the results thereof may be passed to online texture memory  56  for storage or may be passed to further disk storage reduction techniques  58  before being written to disk storage  60 . When the time is appropriate, a renderor  62  accesses the disk storage  60  to recall the decomposed and compressed graphics information stored therein generates graphics for display therefrom. 
   Although  FIG. 2  shows a number of structures linked together in one example manner, it will be appreciated by persons of ordinary skill in the art that some of the structures may reside apart from the remaining structures of FIG.  2 . For example, the renderor  62  may reside separately from the decomposer  50  and vector quantizer  52 . 
   As shown in  FIG. 3 , an example hardware system  70  on which the system of  FIG. 2  may be implemented includes a processing unit  72 , a memory  74  and an input/output circuit  76  that may be coupled via a first communication bus  78 . The input/output circuit  76  is further coupled, via a second communication bus  80 , to a number of other components that may include, for example, a network  82 , an external memory  84 , an image acquisition device  86 , graphics hardware  87 , a display device  88  and a user interface device  90 . 
   In practice, the components  72 - 76  may be implemented by a conventional personal computer, a workstation or in any other suitable computing device. The network  82  may be implemented by an intranet, the Internet or in any other network that allows the exchange of data across distances. The external memory  84  may be implemented by any suitable magnetic, electrical or optical media-based device. For example, the external memory  84  may be implemented by a floppy disk, a hard disk, a read/write compact disk (CD), a zip disk or any other suitable device or media. The image acquisition device  86  may be implemented by a flatbed or any other type of scanner and/or any digital camera. The graphics hardware  87  may be implemented by most consumer level add-on AGP graphics cards produced today. The display device  88  may be implemented by a conventional cathode ray tube (CRT) display, a liquid crystal display (LCD) or in any other suitable display device that is adapted to present graphics to a user of the system  70 . The user interface device  90  may be implemented by a keyboard, a touchscreen, a mouse or any other suitable device that allows a user to input information into the processing unit  72  via the input/output circuit  76 . 
   Routines (described in detail below) that may be stored in memory  74  or external memory  84  and executed by the processing unit  72 , are used to compress and store graphics information for later use. Additionally, the graphics hardware  87  of the system  70  is used to decompress and render graphics information for display to a user. These two general modes are referred to herein as image preparation and image rendering, respectively. 
   Generally, during image preparation, graphics information is provided to the processing unit  72  via the image acquisition device  86  and the input/output circuit  76 . Alternatively, the graphics information may be previously stored in either the memory  74  or the external memory  84 . Either way, the processing unit  72  accesses the graphics information, applies a mesh topology to the graphics information and calculates a BSRF for each vertex of the mesh topology. The processing unit  72  then partitions and approximates each vertex BSRF, before storing the partitioned and approximated BSRFs. As described below, partitioning may be carried out using incident/reflected parameterization or Gram-Schmidt parameterization. 
   As described in detail below, when graphics information is to be rendered, the graphics hardware  87  cooperates with the processing unit  72  to retrieve the partitioned and approximated graphics information and to reassemble such information into graphics for display on the display device  88 . In such an arrangement, the graphics may be used for a computer game resident on the system  70  or for any other application in which graphics are used. Alternatively, the system  70  may not have the compressed and approximated graphics information stored thereon. In such cases, the compressed and approximated graphics information may be stored on a remote system (not shown) to which the system  70  may have access via the network  82 . For example, the processing unit  72  may access compressed and approximated graphics information from another system over an intranet, the Intranet or over any other suitable network. Such an arrangement would be desirable in, for example, applications in which Internet web pages use graphics that, without being compressed and approximated, would require an enormous amount of bandwidth and time to display to a client. 
   Turning now to  FIG. 4 , an example image preparation routine  100 , which may be stored in either the memory  74  or the external memory  84  for execution by the processing unit  72  is shown in representative flow diagram format. The routine  100  begins execution by acquiring images for each desired viewing angle and each desired lighting position (block  102 ). For example, referring back to  FIG. 1 , images of the object  20  may be acquired from each camera position  30 ,  32  for each lighting position  22 - 26 . In practice, the images of the object  20  may be acquired from, for example, thirty positions under, for example, thirty different lightings, for a total of 900 images. Each of the 900 images, which may occupy several megabytes or more of storage space, may be stored in either the memory  74  or the external memory  84 . One example view of an acquired object  104  is shown in FIG.  5 . For clarity and illustrative purposes, the remainder of the discussion will focus on the single view of the object  104  shown in FIG.  5 . 
   Returning to  FIG. 4 , after the images are acquired (block  102 ), a mesh topology is input for the object  104  (block  106 ). The mesh topology may be constructed before the images are acquired at block  102 . Alternatively, one or more of the images acquired at block  102  may serve as a starting point for constructing the mesh topology. As a result, each of the views of the object  104  has a mesh topology  108  overlaid thereon, as shown in FIG.  6 . The mesh topology may be, for example, a number of points, three of which are shown in  FIG. 6  at reference numerals  110 A,  110 B and  110 C, that are connected with vectors or lines, three of which are shown at reference numerals  112 A,  112 B and  112 C, to form triangles, one of which is shown at reference numeral  114 A (see FIG.  6 ). Of course, as will be readily appreciated, any suitable polygon other than a triangle may be used to form a mesh topology and the use of triangularly shaped polygons is merely an example of one such polygon topology. 
   After the mesh topology is applied to the acquired images, the processor  72  calculates a full BSRF for each sample (block  120 ). The process by which BSRFs are calculated. is known, but the full BSRF representation entails a vast quantity of data. For example, a resampled vertex BSRF, which is described below, may require on the order of several hundred megabytes of storage space. Accordingly, direct representation and manipulation of BSRFs are not possible because of the 6D nature of BSRFs. As described in detail below, separable factorization of BSRFs and light field mapping are used to compress the information contained in the BSRFs. 
   Because not all triangles are visible from all directions, after the mesh topology is applied (block  120 ), the visibility of each triangle (e.g., the triangle  114 A, see  FIG. 6 ) in the mesh topology  108  is determined from each viewpoint by comparing the images acquired in block  102  to the camera locations used to acquire the images (block  121 ). For example, the triangle  114 A may not be visible in images acquired from a camera angle on the opposite side of the object. 
   The largest reprojection size of each triangle, among all acquired images of that triangle (e.g., the largest size that the triangle appears in any of the acquired images), is used to sample the triangle from all visible acquired images (block  122 ). This process results in, for each triangle, a list of 2D images, each of which has an identical size and shape and is associated with a light source and camera position. The BSRF having the sampled normalized triangles is referred to as a renormalized triangle BSRF, which is saved to the memory  74  or the external memory  84  before resampling is carried out. 
   Resampling of the renormalized triangle BSRFs is then carried out in the lighting and viewing directions, referred to with vectors l and v in  FIG. 1  (block  123 ). For each particular vertex, all of the renormalized triangle BSRF surrounding that vertex are loaded into memory  74 . Triangle views that are not commonly visible from a particular view for the ring of triangles around the vertex are then discarded, which ensures that visual continuity is preserved about each vertex  110  of the mesh topology  108 . For certain places on the object  104  where the curvature around a vertex  110  is high, the constraint of discarding triangles that are not commonly visible may be relieved. 
   If the BSRF function is to be manipulated in incident/reflected parameterization, then regularly resampling on the lighting and viewing directions through linear interpolation of the three closest ring triangle view is carried out. This dictates that the BSRF function is evaluated at the point (s,v,l), where variables (v,l) take values from the set V n×n ×L n×n , which is shown in Equation 4 below.
 
 V   n×n ={( x,y )ε G|x   2   +y   2 ≦1},
 
 L   n×n ={( x,y )ε G|x   2   +y   2 ≦1}
 
Where G is a grid of points defined by 
             G   =     {       (       k     n   -   1       ,     l     n   -   1         )     ,     k   =       0   ⁢           ⁢   …   ⁢           ⁢   n     -   1       ,     l   =       0   ⁢           ⁢   …   ⁢           ⁢   n     -   1         }             Equation   ⁢           ⁢   4             
 
   In the general case, the resampling is a four dimensional (4D) interpolation problem. However, because of the data acquisition process, this problem can be reduced to a series of 2D interpolation problems. If the lighting direction is denoted as l i  and for each lighting direction l i , a set of images are taken at viewing directions v i   j , all vectors l i  may be collected and projected into 2D using an XY-map, and a 2D Delaunay triangulation may be constructed. For a fixed i, all vectors v i   j  are collected and a 2D Delaunay triangulation for each of the groups is generated. These triangulations establish interpolation functions that allow the closest lighting and viewing directions to be queried together with the interpolation weights shown in Equation 5 below.
 
 q   l   :l →{( w   k     1     ,l   k     1   ),( w   k     2     ,l   k     2   ),( w   k     3     ,l   k     3   )}
 
 q   v   i   :v →{( w   i   m     1     ,v   i   m     1   ),( w   i   m     2     ,v   i   m     2   ),( w   i   m     3     ,v   i   m     3   )}  Equation 5
 
   For incident/reflected parameterization, which is discussed below, the BSRF function may be resampled as shown in Equation 6. 
                     f   ⁡     (     s   ,   v   ,   l     )       =       ∑     i   =   1     3     ⁢           ⁢       w     k   i       ⁢     f   ⁡     (     s   ,   v   ,     l     k   i         )                       =       ∑     i   =   1     3     ⁢           ⁢       w     k   i       ⁢       ∑     j   =   1     3     ⁢           ⁢       ω     m   j       k   i       ⁢     f   ⁡     (     s   ,     v     m   j       k   i       ,     l     k   i         )                             Equation   ⁢           ⁢   6             
 
   In the case of half-vector parameterization, which is also discussed below, the resampling process is similar. First grids for variables h and d are formed as shown in Equation 7.
 
 H   n×n ={( x,y )ε G|x   2   +y   2 ≦1},
 
 D   n×n ={( x,y )ε G|x   2   +y   2 ≦1}  Equation 7
 
The resampled function ƒ(s, h, d) is then computed by first transforming the vector (h, d) into the incident/reflected parameterization (l, v) and applying Equation 6 for resampling.
 
   In BSRF representation, the visibility of the light sources with respect to a surface do not need to be calculated, which implies that an image of a surface back-facing a light source should be incorporated in the BSRF. However, because the radiances of these images are typically quite small, the disclosed system parameterizes the light source direction on an upper hemisphere with respect to the surface, and discards the data from backward-facing light sources. 
   After resampling is carried out (block  123 ), a routine  124  is executed by the processor  72  to partitions the BSRFs. The routine  124  is shown in further detail in FIG.  7 . In general, to retain scalability of the system, the BSRFs are partitioned into local functions around vertices  110  of the mesh topology  108  (see  FIG. 6 ) (block  126 ). The partitioning used is similar to that disclosed in Chen, “Light Field Mapping: Efficient Representation and Hardware Rendering of Surface Light Fields,” to Appear in the Proceedings of SIGGRAPH 2002 (July 2002). The partitioning process divides the large BSRF function ƒ into multiple and separable vertex BSRFs ƒ v     j   , where ƒ v     j    is defined in Equation 8. 
    ƒ v     j   ( s,l,v )=Λ v     j   ( s )ƒ( s,l,v )  Equation 8 
   In Equation 8, Λ v     j    is referred to as a “hat function” that is defined as the barycentric weight around vertex v j , as illustrated in  FIGS. 8A-8D . Hat functions are calculated for each vertex (block  128 ). Generally, as shown in  FIGS. 8A-8D , a hat function is calculated for each vertex (e.g.,  110 A- 100 C) defining the corners of a triangle of interest (e.g.,  114 A). The vertex hat functions represent the radiance contributions of mesh components surrounding each vertex of the triangle of interest. Because the hat functions are barycentric, the sum of the vertex hat function&#39;s surrounding a particular triangle is unitary. 
   In particular, referring to  FIG. 8A , when determining the radiance contribution of vertex  110 A to the total radiance of the triangle  114 A, triangles  114 A- 114 F surrounding vertex  110 A are considered and the hat function Λ v     1    is calculated based thereon. Similarly, as shown in  FIG. 8B , the contribution by vertex  110 C is determined by considering contributions from triangles  114 A,  114 B and  114 G- 114 I and calculating Λ v     2    therefrom. The contribution due to the vertex  110 B, as shown in  FIG. 8C , is determined by considering the contributions from triangles  114 A,  114 F,  114 G,  114 J and  114 K and calculating Λ v3 . As shown in  FIG. 8D , the radiance of triangle  114 A is approximated as a sum of the hat functions of the surrounding vertices  110 A- 110 C. Within each triangle Δ i , the sum of all three vertex BSRFs thereof is equal to the original function ƒ. Because the hat functions are C 0  continuous, the set {f vj }}, j=1,2, . . . N defines a C 0  partitioning of the target BSRF over vertices V j , j=1,2, . . . N. 
   The hat functions calculated at block  128  are used to approximate each vertex BSRF independently without producing visible discontinuity artifacts (block  130 ). If a(i 1 , i 2 , i 3 ) is denoted as vector a expressed in the coordinate system spanned by (i 1 , i 2 , i 3 ), a partitioned vertex BSRF is expressed as ƒ v     j   (s, l(x, y, z), v(x, y, z)), where (x, y, z) is the global coordinate system. The vertex BSRFs are determined by multiplying the full BSRF by the hat functions of each vertex. The result is a number of vertex BSRFs that correspond to each vertex  110  in the mesh topology  108 . 
   After individual vertex BSRFs ƒ v     j   , are determined at block  130 , each vertex BSRF is reparameterized to a local coordinate system (block  132 ). Reparameterization may be carried out using incident/reflected parameterization or by using the Gram-Schmidt Half-vector, each of which is discussed below. 
   Incident/reflected parameterization entails reparameterizing functions ƒ v     j    on a local coordinate system for each vertex v j  to obtain the classic incident/reflected parameterization with respect to the newly assigned local coordinate system. Each vertex coordinate system is defined by the orthonormal vectors (n, s, t), where n is a vector that is normal to the surfaces of the triangles surrounding the vertex under consideration and s and t are the remaining orthonormal components. The surface parameters s are parameterized using the barycentric coordinate around vertex v j . Accordingly, the reparameterized functions ƒ v     j    may be defined ƒ v     j   (s, l(n, s, t), v(n, s, t)). The vector mathematics associated with incident/reflected parameterization will be readily apparent to those having ordinary skill in the art. 
   Alternatively, parameterization using the Gram-Schmidt half vector may be carried out to further improve the separability of the vertex BSRF function. Such a parameterization is described in Kautz, J., “Interactive Rendering with Arbitrary BRDFs using Separable Approximations,” Eurographics Rendering Workshop 1999 (June 1999). A graphical representation of the process of Gram-Schmidt half vector reparameterization is shown in  FIG. 9  with respect to vertex  110 C. First, a half vector h 1 , which has a direction equivalent to l 1 +v 1 , where l 1  is the lighting vector and v 1 , is the view vector, is calculated. The vector h 1  is referred to as the half vector because, as shown in  FIG. 9 , it is half way between the l 1  and v 1  vectors. As will be readily appreciated, each vertex of the mesh topology will have a different half vector because the vectors between a given light source and a given viewing location change as the location of the vertex being considered changes. 
   After the vector h is determined for a vertex, a new coordinate system is defined as (h, s′, t′) by the Gram-Schmidt orthonormalization of vectors (h, s, t). The reparameterized vertex BSRF is then expressed as ƒ v     j   (s, h(n, s, t), d(h, s′, t′)). After each vertex BSRF has been reparameterized, control returns to routine  140  of  FIG. 4 , which approximates the partitioned BSRFs. 
   As shown in  FIG. 10 , at block  142  the decomposer  50  decomposes the partitioned 6D BSRFs into 2D and 4D representations. As will be readily appreciated by those having ordinary skill in the art and as described below, a 6D BSRF may be represented as a 3D matrix of 2D functions. After reparameterization, the vertex BSRFs are available as a discrete function through resampling. These functions are represented as ƒ v     j   (s, ω 1 , ω 2 ), where (ω 1 , ω 2 )≡(l(n, s, t), v(n, s, t)) for the classic incident/reflected parameterization, and (ω 1 , ω 2 )≡(h(n, s, t), d(h, s′, t′)) for the Gram-Schmidt half vector parameterization. As a first step, the discrete function is transformed into a 3D matrix F={F l,m,n }, where l, m, n are combined indices of (u, v), h and d, respectively. The components of the 3D matrix F are approximated as shown in Equation 9. 
               F     l   ,   m   ,   n       ≅       ∑     i   =   0       k   -   1       ⁢           ⁢       u   l   i     ⁢     v   m   i     ⁢     w   n   i                 Equation   ⁢           ⁢   9             
 
Equation 9 denotes the k-term approximation of matrix F, in which the u, v and w components may be referred to as surface, view and light maps, respectively. Although Principal Component Analysis (PCA) algorithms can be used to calculate 2D matrix factorization, these algorithms need to be modified to accommodate for 3D matrix factorization. As described below, a two-step process is used that first initializes the principal vectors and then refines them to obtain better approximation.
 
   To carry out the initialization process, a modified and efficient variant of principal component analysis (PCA), called EMPCA is used. EMPCA was proposed by Roweis in “EM Algorithms for PCA and SPCA,” in the book entitled “Advances in Neural Information Processing Systems,” vol. 10, which was published by The MIT Press in 1998. EMPCA is efficient at decomposing large matrices when only a few eigenvectors corresponding to the largest eigenvalues are needed. The primary advantage of EMPCA is its online mode, which allows calculation of PCA without simultaneously loading the entire matrix into memory  74 . Briefly, EMPCA uses interpolated data generated in the resampling stage as online training vectors for the decomposition algorithm. The decomposition process may be implemented using the Math Kernal Library from Intel Corporation, which contains an optimized BLAS library for the PC. 
   The initialization process operates as described below. First, the 3D matrix F is rearranged to from dimensions of L×M×N into a 2D matrix Γ={Γ l,q }. In this process, the angular parameter dimensions are combined into one index q, and the matrix Γ becomes a L×MN matrix. Next a one-term PCA of Γ, denoted as u 0 (γ 0 ) T , where u 0  is a left eigenvector of Γ, which is referred to as the surface map and where γ 0  is a combination of the view and light maps. At this point, the decomposer  50  stores the surface map u 0  in either the memory  74  or the external memory  84  (block  144 ). 
   The decomposer  50  then, decomposes the remaining 4D representation into two, 2D representations by rearranging γ 0  into a 2D matrix and taking a one-term PCA thereof (block  146 ). The one-term PCA is denoted as v 0 (w 0 ) and this estimation of γ 0  is denoted as γ 0′ . The terms v 0  and w 0  represent the view and light maps, respectively. After developing the 2D representations, (block  146 ), the view and light maps (v 0  and w 0 , respectively) are stored (block  148 ). 
   After the view and light maps have been stored (block  148 ), the decomposer  50  determines if more approximations are desired. Typically, for PCA, three approximations are sufficient, but, as will be readily appreciated, any other number of suitable approximations may be used. If more approximations are desired, control passes from block  150  back to block  142 . 
   In one implementation, XY-maps may be used to generate texture maps for α i  and β i  in the decomposed approximations. An XY-map texture map reprojects a 3D unit vector onto a 2D surface by discarding the Z component of the 3D vector. Because the Z components of both the incident/reflected parameterization and the Gram-Schmidt parameterization are positive, XY-mapping does not present a problem. 
   The full parameterization on lighting directions will require one additional XY-map for either of the vectors l and h. Such an implementation requires the rendering process to be slightly modified to bind the correct texture map before rendering. 
   To improve rendering efficiency, the texture maps α i  and β i  are collected into larger texture maps. The collection process is performed to improve the access coherency of texture maps. The collection process is carried out by dividing the model into server connected partitions. The texture maps of g i , α i  and β i  corresponding to the same partition and same approximation term are then collected to form larger texture maps. The ultimate size of the collected texture maps is a trade-off between graphics memory bus utilization efficiency and overhead associated with texture switching. 
   If more approximations are desired, the error approximation shown in Equation 10 is used as the input matrix for finding the next approximation terms, which are shown in Equation 11. 
                     F     l   ,   m   ,   n     1     =       F     l   ,   m   ,   n       -       F   ~       l   ,   m   ,   n     0                   =       F     l   ,   m   ,   n       -       u   l   0     ⁢     v   m   0     ⁢     w   n   0                       Equation   ⁢           ⁢   10             
  {tilde over (F)}   i   l,m,n,   =u   l   i   v   m   i   w   n   i ,
 
 F   l,m,n   i+1   =F   l,m,n,   i   −{tilde over (F)}   i   l,m,n   Equation 11
 
The process proceeds until a desirable error threshold or number of approximation terms is reached.
 
   To prove the convergence of the above process, Equation 11 can be rewritten in 2D form as Equation 12, shown below.
 
Γ i+1 =Γ i   −u   i (γ i′ ) T   Equation 12
 
   The singular-value decomposition of Γ is denoted as UV T  where the singular values are absorbed into V for simplicity. The vector u i  is the singular vector corresponding to the largest singular value of Γ and is thus equal to the first column of U. The variable v j  is denoted as the j th  column of V. Equation 12 above now takes the form of Equation 13. 
               Γ     i   +   1       =       U   ⁢           ⁢     V   T       -         u   i     ⁡     (     γ     i   ′       )       T                     =     U   ⁡     (     V   ′     )         ,             
 
where
 
 V ′=[( v   0   −g   i ),  v   1   , . . . , v   MN ]  Equation 13
 
The vector (v 0 −g i )/|v 0 −g i | is a right eigenvector of Γ i+1  corresponding to the left eigenvector u i . Because g i  is a linearized one-term approximation of v 0 , |v 0 −g i |&lt;|v 0 |. Accordingly, |v 0 | and |v 0 −g i | are singular values of Γ i  and Γ 1+1 , respectively, and all singular values of Γ i+1  are otherwise equivalent to the ones in Γ i . Because the squared sum of elements in a matrix is equivalent to the squared sum of all its singular values, the squared sum of Γ i+1  is thus lower than Γ i , and, therefore, Equation 11 converges. Alternatively, if no more approximations are desired, the routine  140  returns control to the routine  100  of FIG.  4 .
 
   As an alternative to the foregoing described process, to further enhance the results of decomposition, a pseudo-single value decomposition (SVD) technique has been proposed independently by Carroll and Chang in “Analysis of Individual Differences in Multidimensional Scaling Via an N-Way Generalization of Eckart-Young Decomposition,” published in Psychometrica 35, 3 (1970), pp 283-319. A similar approach was also suggested by Harshman in “Foundations of Parafac Procedure: Models and Conditions for an Explanatory Multi-Modal Factor Analysis,” published in the UCLA Working Papers in Phonetics 16 (1970), pp. 1-84. For a budget of k-term approximation, these processes find a local minimum by fixing two out sets of vectors and linearly solving for the third set. The process proceeds as follows. First, an error matrix, as shown in Equation 14 is defined. 
               E     l   ,   m   ,   n     k     =            F     l   ,   m   ,   n       -       ∑     i   =   0       k   -   1       ⁢           ⁢       u   l   i     ⁢     v   m   i     ⁢     w   n   i                        Equation   ⁢           ⁢   14             
 
   Next, a 3D matrix is defined as A l,m,i =u l   i v m   i . The matrix F is then reorganized into a 2D matrix F′ by combining indices (l, m) into one index. The matrix A is also reorganized into a 2D matrix A′ in a similar fashion. The above error equation 14 can be written in matrix notation as shown below in Equation 15.
 
 E   k   =∥F′   T   −A′w∥   Equation 15
 
Wherein, w is an N×k matrix consisting of vectors w i . Given fixed F′ and A′, error E k  is minimized by setting w=pinv(A′)F′ T , where pinv(·) denotes the pseudo-inverse of a matrix. The process then performs the same steps for u, v and w repeatedly until there is convergence.
 
   The refinement stage can be performed relatively efficiently. Because matrix A′ is full-rank, the pseudo-inverse of pinv(A′) can be calculated without calculating SVD. The bottleneck of multiplying pinv(A′) and F′ T  can be carried out in online mode without bringing the entire matrix F′ into the processing unit  72 . 
   Rendering using the factorized BSRF representation developed above includes manipulation of three group functions (g i , α i  and β i ) per approximation term i for each vertex of each triangle in the mesh topology  108 . The BSRF function is reconstructed by evaluating each group of functions using Equation 3 for each surface element, adding the results together, drawing the triangles with correct visibility and displaying the graphics to the user on the display device  88 . 
   Returning to  FIG. 4 , after the routine  140  has completed execution, block  142  may further compress the light, view and surface maps. This can be done using standard image compression techniques that allow for fast decompression in graphics hardware. One example of such techniques is vector quantization (VQ), which groups light, view and surface maps based on their appearance similarity and then represents each group by the most representative individual of the group. Another example is the use of hardware compression techniques that are commonly implemented on graphics hardware that is available today, such as, for example, S3TC. 
     FIG. 11  illustrates an example image rendering routine  160 , which may be carried out by the graphics hardware  87 , the processing unit  72  or any combination thereof. The routine  160  begins operation at block  162 , where a frame buffer and depth buffer, each of which may be in the graphics hardware  87 , are cleared. Control then passes to block  164 , which selects the next mesh component (e.g., the next mesh polygon, which in the foregoing exemplary description, is a triangle), before passing control to block  166 , which selects the next approximation of the BSRF that is to be rendered. The graphics hardware  87  loads the light, surface and view maps for the next vertex of the selected approximation and the selected mesh component (block  168 ). 
   At block  170 , the light, surface and view maps are decoded for the selected mesh component and the selected approximation and the results of the decoding are added to the frame buffer with the correct visibility. Vertices are colored using the intended light source irradiance. The decoded BSRF is obtained by modulating the vertex color with the product of texture maps g, l and h. 
   Block  172  determines if there are more vertices in the mesh component and if there are, passes control back to block  168 . If, however, there are no more vertices to be considered in the selected mesh component, control passes from block  172  to block  174 , which determines if more approximations are desired. If more approximations are desired, control passes back to block  166 . Conversely, if no more approximations are desired, control passes from block  174  to block  176 , which determines if there are more mesh components to consider. If there are more mesh components to consider, control passes back to block  164 . Alternatively, if all vertices of all mesh components have been considered for the appropriate number of approximations, the routine  160  ends and the frame buffer is displayed to the user. 
   The rendering process described above may be implemented in efficient BSRF renderers that exploit commodity graphics hardware features. First of all, decomposed BSRF functions are 2D functions and can be stored in the graphics hardware  87  as textures. The decoding of BSRFs (performed at block  170 ) requires pixel-by-pixel multiplication, which is supported in the graphics hardware  87  via multi-texturing operations. Multi-term approximations are rendered using multiple passes of the routine  160  and on each pass information is written into the frame buffer. 
   While the foregoing description of  FIG. 11  addresses rendering in a general manner, the following description provides further detail on two cases of hardware-accelerated rendering, namely point light source and environment mapping. Because each rendering pass for each triangle proceeds identically, for simplicity, the following focuses on rendering one term approximation for one vertex BSRF within a triangle. The corresponding decomposed BSRF maps are denoted as simply as g, l, and h. 
   When rendering a decomposed BSRF scene with one single point light source, three texture units are used in the hardware, one for each of the texture maps g(s), α(ω 1 ) and β(ω 2 ). For all three vertices within the triangle, we calculate the texture coordinates for l i  and h i  based on the light source position and viewing direction. The function g i  is parameterized on surfaces and thus does not need to be calculated on a per-frame basis. The texture coordinate within the triangle are linearly interpolated, which works well for both incident/reflected and Gram-Schmidt Half-vector parameterizations. In cases involving multiple light sources, the rendering procedure may simply render multiple passes using the process described in  FIGS. 11 and 12 . Accordingly, the rendering cost grows linearly with the number of light sources. 
   When the scene include only directional light sources, and when the incident/reflected parameterization is used, the number of texture units can be reduced by one, because the texture coordinate for α(ω 1 ) is identical everywhere within the triangle. Two aspects of the rendering process may be modified as follows. First, a light source direction ω 0   1  is calculated with respect to the local coordinate on the vertex, and α 0 =α(ω 0   1 ) is evaluated. The term α 0  is them multiplied by light source intensity I and the vertex color is set as I·l 0 . Then the vertex color is modulated with the product of the texture map g and β to obtain the decoded BSRF. Second, because a distant viewer will have almost a parallel projection, the texture coordinates of β may be identical for all vertices and, therefore, the number of required texture units may be reduced to one. 
   To render environments, the rendering equation for the BSRF representation in Equation 2 can be rewritten in discrete, decomposed format, as shown in Equation 16. 
                     I   ⁡     (     s   ,   v     )       =       ∑   l             ⁢           ⁢       f   ⁡     (     s   ,   l   ,   v     )       ⁢     E   ⁡     [   l   ]                       =       ∑   l             ⁢       ∑     i   =   0       k   -   1       ⁢           ⁢         g   i     ⁡     (   s   )       ⁢       α   i     ⁡     (     ω   1     )       ⁢       β   i     ⁡     (     ω   2     )       ⁢     E   ⁡     [   l   ]                           Equation   ⁢           ⁢   16             
 
Wherein E[l] is a subsampled environment map.
 
   In the case of incident/reflected parameterization, the above equation can be written as shown in Equation 17. 
               I   ⁡     (     s   ,   v     )       =       ∑     i   =   0       k   -   1       ⁢           ⁢         g   i     ⁡     (   s   )       ⁢       β   i     ⁡     (   v   )       ⁢     (       ∑   l             ⁢         α   i     ⁡     (   v   )       ⁢     E   ⁡     [   l   ]           )                 Equation   ⁢           ⁢   17             
 
In Equation 17, the inner summation is independent of the output parameters (s, v), and once this summation is evaluated, it may be assigned in the vertex color and two texture units may be used for each of g and β to calculate the final rendering.
 
   The inner summation of Equation 17 represents a vector-vector inner product. If functions α i  are parameterized for all vertices in the global coordinate system, the inner product may be calculated efficiently for all vertices using a matrix-vector multiplication operation. For a 3-term decomposition model with 3000 vertices and a 16×16 size of α i , this multiplication requires 2.3 million MPEG audio decoder (MAD) operations. 
   Although certain apparatus constructed in accordance with the teachings of the invention have been described herein, the scope of coverage of this patent is not limited thereto. On the contrary, this patent covers all embodiments of the teachings of the invention fairly falling within the scope of the appended claims either literally or under the doctrine of equivalents.