FIELD OF THE INVENTION

[0001]
The invention relates generally to the field of systems and computerimplemented methods for evaluating integrals, and more particularly to such systems and computerimplemented methods that evaluate integrals using Monte Carlo and quasiMonte Carlo methodologies. The invention particularly provides a new and improved system and computerimplemented method for evaluating integrals using a Monte Carlo methodology in which the integration domain is stratified using a stratification methodology in which the number of strata is independent of the number of dimensions in the integration domain. Specifically, the invention provides a methodology that makes use of stratification based on rank1 lattices. Systems and computerimplemented methods according to the invention find utility in a number of applications, including but not limited to computer graphics.
BACKGROUND OF THE INVENTION

[0002]
In computer graphics, a computer is used to generate digital data that represents the projection of surfaces of objects in, for example, a threedimensional scene, illuminated by one or more light sources, onto a twodimensional image plane, to simulate the recording of the image by, for example, a camera. The simulated camera may include a lens for projecting the image of the scene onto the image plane, or it may comprise a pinhole camera in which case no lens is used. The twodimensional image is in the form of an array of picture elements (which are variably referred to as “pixels” or “pels”), and the digital data generated for each pixel represents the color and luminance of the scene as projected onto the image plane at the point of the respective pixel in the image plane. The surfaces of the objects may have any of a number of characteristics, including shape, color, specularity, texture, and so forth, which are preferably rendered in the image as closely as possible, to provide a realisticlooking image.

[0003]
Generally, the contributions of the light reflected from the various points in the scene to the pixel value representing the color and intensity of a particular pixel are expressed in the form of the one or more integrals of relatively complicated functions. Since the integrals used in computer graphics generally do not have a closedform solution, numerical methods must be used to evaluate them to generate the pixel value. Typically, a methodology referred to as the “Monte Carlo” method has been used in computer graphics to numerically evaluate the integrals. Generally, in the Monte Carlo method, to evaluate an integral
<ƒ>=∫_{[0,1[} _{ s }ƒ(x)dx (1)
where f(x) is a real function defined on an integration domain comprising the “s”dimensional unit cube [0,1)^{s }(that is, an sdimensional cube each dimension of which includes “zero,” and excludes “one”),a random point set P_{r}=(ξ_{1}. . . ξ_{r})⊂[0,1)^{s }of “r” independently realized points ξ_{i}, i=1, . . . , r, is generated over the integration domain. The random points ξ_{i }are used as sample points for which sample values f(ξ_{i}) are generated for the function f(x). An estimate {overscore (ƒ)} for the value of the integral (equation (1)) is generated as
$\begin{array}{cc}\langle f\rangle \approx \stackrel{\_}{f}=\frac{1}{r}\sum _{i=1}^{r}f\left({\xi}_{i}\right)& \left(2\right)\end{array}$
As the number “r” of sample points that are used in generating the sample values f(ξ_{i}) increases, the value of the estimate {overscore (ƒ)} will converge toward the actual value of the integral <ƒ>. Generally, the distribution of values of the estimate {overscore (ƒ)} that will be generated for various values of “r,” that is, for various numbers of sample points ξ_{i }are distributed around the actual value <ƒ> with a standard deviation σ which can be estimated by
$\begin{array}{cc}\sigma =\sqrt{\frac{1}{N1}\left(\stackrel{\_}{{f}^{2}}{\stackrel{\_}{f}}^{2}\right)}& \left(3\right)\end{array}$
if the sample points ξ_{i }used to generate the sample values f(ξ_{i}) are statistically independent, that is, if the points ξ_{i }are truly positioned at random in the integration domain. This yields the probabilistic error bound
$\begin{array}{cc}\mathrm{Prob}\left(\uf603{\int}_{{\left[0,1\right)}^{s}}f\left(x\right)dx\frac{1}{r}\sum _{k=1}^{r}f\left({\xi}_{k}\right)\uf604\le \frac{3\sigma \left(f\right)}{\sqrt{r}}\right)\approx 0.997.& \left(4\right)\end{array}$
The rate of convergence of
$O\left(\frac{1}{\sqrt{r}}\right)$
is independent of the smoothness of the integrand f(x) and the number “s” of dimensions comprising the integration domain.

[0004]
Several strategies can be used to improve convergence of the estimate {overscore (ƒ)} to the actual value <ƒ> of the integral. One strategy, which is described in the aforementioned Grabenstein patent application, is to use a lowdiscrepancy, strictly deterministic point methodology to generate the set of sample points that will be used in generating the sample values f(ξ_{i}) for the estimate. A lowdiscrepancy, strictly deterministic methodology will ensure that the sample points are distributed throughout the integration domain [0,)^{s }without clumping in particular regions, which can occur if the sample points are randomly distributed.

[0005]
Another strategy to improve the convergence is to stratify the integration domain [0,1)^{s }into a number “K” of disjoint strata A_{j}, such that [0,1)^{2}=∪_{j+1} ^{K}A_{j}, and separately evaluate the resulting integrals
$\begin{array}{cc}{\int}_{{\left[0,1\right)}^{s}}f\left(x\right)dx=\sum _{j=1}^{K}{\int}_{{A}_{j}}f\left(x\right)dx& \left(5\right)\end{array}$
using a sample point ξ_{i }obtained from each stratum A_{j}. Stratifying the integration domain [0,1)^{s }can also be helpful in ensuring that sample points are reasonably well distributed throughout the integration domain [0,1)^{s}, and that they do not clump in the integration domain, which can occur if the locations of the sample points ξ_{i }are randomly located in the integration domain.

[0006]
Several stratification strategies have been proposed. Typically, previouslyproposed stratification strategies, particularly stratification strategies that are based on axisaligned grids, suffer exponential growth in the number of strata with increasing numbers of dimensions comprising the integration domain. Since integrands ƒ that are used in fields such as, for example, computer graphics, often have very highdimensional integration domains, stratifying domains using such stratification strategies can be computationally quite intensive.
SUMMARY OF THE INVENTION

[0007]
The invention provides a new and improved system and computerimplemented method for evaluating integrals using a Monte Carlo methodology in which the integration domain is stratified using a stratification methodology in which the number of strata is independent of the number of dimensions in the integration domain. Specifically, the invention provides a methodology that makes use of stratification based on rank1 lattices.

[0008]
In brief summary, in one aspect, the invention provides a system for numerically evaluating an integral of a function over an sdimensional integration domain, the system comprising a sample point generator, a function value generator and an integral value estimate generator. The sample point generator is configured to generate a selected number of sample points over the integration domain, the sample points being generated such that there is at least one sample point in each of a plurality of strata distributed over the integration domain, the strata being defined by a rank1 lattice. The function value generator is configured to, for respective ones of the sample points, generate a value for the function at the respective sample point. The integral value estimate generator is configured to use the function values generated by the function value generator at the respective sample points in generating an estimate for the value of the integral.

[0009]
In another aspect, the invention provides a computer implemented method of numerically evaluating an integral of a function over an sdimensional integration domain, the method comprising a sample point generation step, a function value generation step and an integral value estimate generation step. During the sample point generation step, a selected number of sample points are generated over the integration domain, the sample points being generated such that there is at least one sample point in each of a plurality of strata distributed over the integration domain, the strata being defined by a rank1 lattice. During the function value generation step, for respective ones of the sample points, a value for the function at the respective sample point is generated. During the integral value estimate generation step, the function values generated during the function value generation step at the respective sample points are used in generating an estimate for the value of the integral.

[0010]
In another aspect, the invention provides computer program product for use with a computer to provide a system for numerically evaluating an integral of a function over an sdimensional integration domain. The computer program product comprises a computer readable medium having encoded thereon a sample point generator module, a function value generator module and an integral value estimate generator module. The sample point generator module is configured to enable the computer to generate a selected number of sample points over the integration domain, the sample points being generated such that there is at least one sample point in each of a plurality of strata distributed over the integration domain, the strata being defined by a rank1 lattice. The function value generator module is configured to enable the computer to, for respective ones of the sample points, generate a value for the function at the respective sample point. The integral value estimate generator module is configured to enable the computer to use the function values generated by the function value generator at the respective sample points in generating an estimate for the value of the integral.
BRIEF DESCRIPTION OF THE DRAWINGS

[0011]
This invention is pointed out with particularity in the appended claims. The above and further advantages of this invention may be better understood by referring to the following description taken in conjunction with the accompanying drawings, in which:

[0012]
FIG. 1 depicts an illustrative computer graphics system that evaluates integrals using a Monte Carlo methodology in which the integration domain is stratified using a stratification methodology in which the number of strata is independent of the number of dimensions in the integration domain;

[0013]
FIG. 2 depicts a table of information that is useful in understanding the operation of the computer graphics system depicted in FIG. 1.
DETAILED DESCRIPTION OF AN ILLUSTRATIVE EMBODIMENT

[0014]
The invention provides a computer graphics system and method for generating pixel values for pixels in an image of a simulated scene that makes use of a Monte Carlo methodology in which an integral or integrals is evaluated of function(s) by using sample points determined by stratifying the integration domain using rank1 lattices. The function(s) represent the contributions of the light that is emitted from simulated light sources and reflected from the various simulated surfaces in the scene directed towards a simulated camera, and the integral(s) provide the pixel values for the respective pixels in the image. Stratifying the integration domain in this manner provides that the number of strata will not grow exponentially with the number of dimensions in the integration domain, which can occur with, for example, axisaligned strata.

[0015]
FIG. 1 attached hereto depicts an illustrative computer system 10 that makes use of such a stratification methodology. With reference to FIG. 1, the computer system 10 in one embodiment includes a processor module 11 and operator interface elements comprising operator input components such as a keyboard 12A and/or a mouse 12B (generally identified as operator input element(s) 12) and an operator output element such as a video display device 13. The illustrative computer system 10 is of the conventional storedprogram computer architecture. The processor module 11 includes, for example, one or more processor, memory and mass storage devices, such as disk and/or tape storage elements (not separately shown), which perform processing and storage operations in connection with digital data provided thereto. The operator input element(s) 12 are provided to permit an operator to input information for processing. The video display device 13 is provided to display output information generated by the processor module 11 on a screen 14 to the operator, including data that the operator may input for processing, information that the operator may input to control processing, as well as information generated during processing. The processor module 11 generates information for display by the video display device 13 using a socalled “graphical user interface” (“GUI”), in which information for various applications programs is displayed using various “windows.” Although the computer system 10 is shown as comprising particular components, such as the keyboard 12A and mouse 12B for receiving input information from an operator, and a video display device 13 for displaying output information to the operator, it will be appreciated that the computer system 10 may include a variety of components in addition to or instead of those depicted in FIG. 1.

[0016]
In addition, the processor module 11 includes one or more network ports, generally identified by reference numeral 14, which are connected to communication links which connect the computer system 10 in a computer network. The network ports enable the computer system 10 to transmit information to, and receive information from, other computer systems and other devices in the network. In a typical network organized according to, for example, the clientserver paradigm, certain computer systems in the network are designated as servers, which store data and programs (generally, “information”) for processing by the other, client computer systems, thereby to enable the client computer systems to conveniently share the information. A client computer system which needs access to information maintained by a particular server will enable the server to download the information to it over the network. After processing the data, the client computer system may also return the processed data to the server for storage. In addition to computer systems (including the abovedescribed servers and clients), a network may also include, for example, printers and facsimile devices, digital audio or video storage and distribution devices, and the like, which may be shared among the various computer systems connected in the network. The communication links interconnecting the computer systems in the network may, as is conventional, comprise any convenient informationcarrying medium, including wires, optical fibers or other media for carrying signals among the computer systems. Computer systems transfer information over the network by means of messages transferred over the communication links, with each message including information and an identifier identifying the device to receive the message.

[0017]
It will be helpful to initially provide some background on operations performed by the computer graphics system 10 in generating an image. Generally, the computer graphics system 10 generates an image that attempts to simulate an image of a threedimensional scene as would be recorded by a camera. In that operation, for each pixel of the camera's image plane, the computer graphics system 10 numerically evaluates one or more integrals of functions that represent light that is directed toward the respective pixel from the scene. An illustrative function will be described below. The functions typically contain a number of terms, including one or more terms that represent light that is directly incident on the respective pixel from one or more light source sources, as well as one or more terms that represent light that, after having been emitted by a light source, is reflected by one or more surfaces of objects in the scene and ultimately directed toward the pixel.

[0018]
To numerically evaluate the integral, the computer graphics system 10 shoots, or traces, simulated rays from a plurality of points in the pixel into the scene. At points in the scene at which the shot rays intersect surfaces of objects in the scene, the computer graphics system determines the contributions of light that are directed from those intersection points in the scene toward the respective points in the pixel from which the rays were shot. The computer graphics system can determine the contributions that are directed to the respective points in the pixel in a number of ways. In a path tracing methodology, which is described in more detail below, the computer graphics system 10 determines the illumination at the point in the scene at which the ray shot from the pixel intersects an object in the scene by shooting one or more rays from the intersection point in a plurality of directions, and determining the extent to which that point is illuminated directly by light sources, and indirectly by light that is reflected off of other surfaces of objects in the scene. After the computer graphics system has determined the extent to which the intersection point is illuminated, the direction or directions from which the point is illuminated, and other information that will be known by those skilled in the art, it generates a value that represents the light reflected toward the point in the pixel from which the ray was shot. The computer graphics system will repeat these operations for the various points in the pixel from which rays are shot, and will determine a pixel value for the entire pixel from the various values that are generated therefor. The computer graphics system also repeats the operations for all of the pixels to generate the final image. In a bidirectional path tracing methodology, the computer graphics system also shoots rays from light source or sources into the scenes and determines the pixel values in relation to points in at which those rays intersect points on surfaces of objects in the scene.

[0019]
In a photon map methodology, prior to shooting rays from the pixels in the image plane into the scene, the computer graphics system 10 first generates a photon map representing photons that are incident on the surfaces of objects in the scene directly from the light sources and indirectly by photons that are reflected from various surfaces in the scene. After the computer graphics system has generated the photon map, it will shoot rays from the points on the pixels into the scene, and generate the pixel value in relation to the photons, if any, that are proximate the point of intersection of those rays with points on surfaces of the objects in the scene.

[0020]
In both methodologies, as well as other methodologies as will be appreciated by those skilled in the art, the computer graphics system 10 uses the coordinates of the sample points that are obtained from the integration domain of the integral to control various aspects regarding the geometry while generating the image. These aspects may include, for example, determining the locations of the points on the various pixels from which rays are shot, determining the locations on light sources from which rays representing photons are shot, determining directions of reflection, determining whether to terminate a ray at a surface or to allow the ray to be reflected, and other aspects as will be aware to those skilled in the art.

[0021]
As noted above, formally, the operations performed by the computer graphics system in generating pixel values for an image comprise the numerical evaluation of an integral over an integration domain that is defined over an sdimensional unit cube [0,1)^{s}. In numerically evaluating the integral, the computer graphics system partitions the integration domain [0,1)^{s }into a number “K” of disjoint strata A_{j}, such that [0,1)^{s}∪=_{j=1} ^{K}A_{j}, and separately evaluates the resulting integrals
$\begin{array}{cc}{\int}_{{\left[0,1\right)}^{s}}f\left(x\right)dx=\sum _{j=1}^{K}{\int}_{{A}_{j}}f\left(x\right)dx& \left(6\right)\end{array}$
using a sample point obtained from each stratum A_{j}. The strata A_{4 }are defined by an “s” dimensional lattice. The lattice divides the integration domain into a plurality of strata A_{j }of equal size, and the computer graphics system 10 obtains a sample point from each stratum A_{j }for use in generating the approximation for the integral.

[0022]
A lattice is a discrete set of points L:=P_{n}+Z^{s}⊂R^{s }that is closed under addition and subtraction, where P_{n}:={z_{0}, . . . , z_{n1}l }, “Z” is an integer, “R” is a real number, and superscript “s” corresponds to the dimension of the lattice. In the case of a lattice defined by an axisaligned regular grid, the dimension “s” is the number of basis vectors that, in linear combinations, yield the lattice points. According to another methodology, and in accordance with the invention, the lattice is defined using a single generator vector g ε N^{s}, where “N” is a member of the natural numbers, with the coordinates of the lattice points being defined by
$\begin{array}{cc}{z}_{j}:=\frac{j}{n}g,& \left(7\right)\end{array}$
where j=0, . . . ,n1. The lattice is referred to as a “rank1” lattice, because it is generated using a single generator vector g.

[0023]
An example of a rank1 lattice that will be used herein is a lattice that is generated based on the Fibonacci sequence. The Fibonacci sequence {F_{k}} is defined as F_{k}:=F_{k−1}+F_{k−2}, with F_{1}=F_{2}=1 and index “k” an integer. The twodimensional Fibonacci lattice, which is based on the Fibonacci sequence, for n=F_{k}, k≧3, lattice vertices has the generator vector
g=(1, F _{k−1}) (8)
It will be apparent from equations (7) and (8) that the coordinates (α,β) of the vertices z_{j }comprising a twodimensional Fibonacci lattice on the unit square [0,1)^{2 }can be readily generated in a “for”loop, in which the first component “α” is
$\u201cj\frac{1}{{F}_{k}}\u201d$
and the second component “β” is
$\u201cj\frac{{F}_{k1}}{{F}_{k}}\mathrm{mod}\text{\hspace{1em}}1\u201d,$
since, as noted above, n=F_{k}, the “kth” element of the Fibonacci sequence.
Fibonacci lattices of dimension “s” higher than “two” can be defined in a similar manner using an “s”dimensional generator vector. For example, the generator vector
$g=\left(1,{F}_{k1},{F}_{k1}^{2},\dots \text{\hspace{1em}},{F}_{k1}^{s1}\right)$
can be used, which would provide a lattice in Korobov form.

[0024]
As noted above, the computer graphics system 10 obtains sample points from strata A_{j }that are defined by the lattice. The strata relate to unit cells that are associated with the respective vertices z_{j}. The unit cell of a rank1 lattice is not uniquely specified, but a suitable unit cell, which is referred to as a Delaunay unit cell, can be specified such that the volume of an sdimensional sphere that is inscribed in the unit cell is maximized. The unit cells that result on the rank1 lattice are all of similar shape and can be mapped from one to another by means of rigid body transformations that basically perform a translation and/or rotation to move them from vertex to vertex.

[0025]
Each unit cell is identified by the Voronoi diagram of the rank1 lattice, whose nerve is constructed out of the “s” shortest vectors v_{1}, . . . , v_{s }from the origin of the coordinate system defining the lattice, to points of the lattice, where the value of the first component of each vector v_{j} ^{(1) }is greater than zero. This implies that, for all vectors, v_{j}, v_{j}≠0, and further that the value of “n,” the number of vertices of the lattice, is greater than “one.” The vectors span an sdimensional parallelopiped
B:=(v_{1 }. . . v_{s}) (9)
whose volume
$\uf603B\uf604:=\uf603\mathrm{det}\text{\hspace{1em}}B\uf604=\frac{1}{n}.$
Each of the “n” sdimensional parallelopipeds in a rank1 lattice restricted to the unit cube [0,1)^{s }is a stratum A_{j }that is induced by the bijection
$\begin{array}{cc}{R}_{j}:{\left[0,1\right)}^{s}>{A}_{j}\text{}x\mapsto \left(\frac{j}{n}\xb7g+B\text{\hspace{1em}}x\right)\text{\hspace{1em}}{\mathrm{mod}\text{\hspace{1em}}\left[0,1\right)}^{s}& \left(10\right)\end{array}$
where the jth lattice vertex
${z}_{j}=\frac{j}{n}\xb7g$
is used as an offset into the “s”dimensional unit cube [0,1)^{s }and “B,” which is used as a matrix with the vectors v_{j }comprising the columns of the matrix, defines a linear transformation that represents a rotation and/or shear. Equation (10) defines what can be referred to as a “reduced CranleyPatterson rotation,” which determines the location of the “jth” sample point at which the function ƒ will be evaluated for the contribution to the estimate {overscore (ƒ)} of the value of the integral. Generally, the first term in the second line of equation (10), namely,
$\u201c\frac{j}{n}\xb7g\u201d,$
j=0, . . . ,n1, identifies the locations of the vertices for the basic “unrotated” rank1 lattice. The second term in the second line of equation (9), namely, B x, essentially provides a displacement from locations of the respective unrotated vertices to the particular sample point that will be used in evaluating the function. As will be described below, the computer graphics system 10 may use the same value for “x” for some or all of the vertices (reference the discussion below regarding correlated sampling), or it may use different values (reference the discussion below regarding jittered sampling). In equation (10), if different values are to be used for “x,” “x” may be jth element of an sdimensional point set P_{r}=(ξ_{1}, . . . ξ_{r})⊂[0,1)^{s}, which may be generated using a random methodology, a strictlydeterministic lowdiscrepancy methodology, or other methodologies that will be apparent to those skilled in the art.

[0026]
The particular parallelopiped that is defined by B will depend on the particular form of the generator vector g that is used to define the lattice, the number “n” of vertices that are to comprise the lattice, the number of dimensions, and other criteria as will be appreciated by those skilled in the art. For the twodimensional Fibonacci lattice with n=F
_{k}>1 vertices, the matrix B is given by
$\begin{array}{cc}B=\frac{1}{{F}_{k}}\left(\begin{array}{cc}{j}_{1}& {j}_{2}\\ {j}_{1}{F}_{k1}& {j}_{2}{F}_{k1}\end{array}\right)& \left(11\right)\end{array}$
where
${j}_{1}={F}_{2\xb7\lfloor \frac{k1}{4}\rfloor +1}\text{\hspace{1em}}\mathrm{and}\text{\hspace{1em}}{j}_{2}={F}_{2\xb7\lfloor \frac{k+1}{4}\rfloor}.$
FIG. 2 depicts a table that, for twodimensional Fibonacci lattices of various numbers “n” of vertices, gives the generator vector g and values of j
_{1 }and j
_{2 }for use in connection with equation (11). As noted above, in the Fibonacci sequence, F
_{1}=F
_{2}=1. Accordingly, the third element of the sequence, F
_{3}, is the sum of the two preceding elements, and so F
_{3}=F
_{2}+F
_{1}=2, which, in turn, corresponds to the number “n” of vertices in the lattice, as shown in the first line of the table depicted in
FIG. 2. For the twodimensional lattice associated with the F
_{3 }element of the Fibonacci sequence

 (i) the generator vector g=(1,F_{k−1})=(1,F_{2})=(1,1);
$\left(\mathrm{ii}\right)\text{\hspace{1em}}{j}_{1}={F}_{2\xb7\lfloor \frac{31}{4}\rfloor +1}={F}_{2\xb70+1}={F}_{1}=1;\text{\hspace{1em}}$
$\left(\mathrm{iii}\right)\text{\hspace{1em}}{j}_{2}={F}_{2\xb7\lfloor \frac{3+1}{4}\rfloor}={F}_{2\xb71}={F}_{2}=1;$
$\mathrm{and}\text{\hspace{1em}}\mathrm{so}\text{}\left(\mathrm{iv}\right)\text{\hspace{1em}}B=\frac{1}{2}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right).\text{\hspace{1em}}$
The values for items (i) through (iii) are depicted in the second, third and fourth columns, respectively, of the table depicted in FIG. 2. It will be apparent that the term
$\frac{j}{n}\xb7g$
in the second line of equation (9), which defines the unrotated lattice, corresponds to vertices having coordinates (0,0) for j=0 and (1/2,1/2) for j=1. The matrix B, in turn, defines strata that are on the order of 1/2by1/2 in size. That is,

[0028]
(a) the stratum that is associated with the vertex having coordinates (0,0) comprises the quadrilateral bounded by line segments between points (0,0), (1/2,0), (1/2,1/2) and (0,1/2), with the quadrilateral including the line segments between points (0,0)(1/2,0) and (0,0)(0,1/2) (except for points (1/2,0) and (0,1/2)), and not including the line segments between points (1/2,0)(1/2,1/2) and (0,1/2)(1/2,1/2), and

[0029]
(b) the stratum that is associated with the vertex having coordinates (1/2,1/2) comprises the quadrilateral bounded by line segments between points (1/2,1/2), (1,1/2), (1,1) and (1/2,1), with the quadrilateral including the line segments between points (1/2,1/2)(1,1/2) and (1/2,1/2)(1/2,1) (except for points (1,1/2) and (1/2,1)), and not including the line segments through points (1,0)(1,1) and (0,1)(1,1).

[0030]
As noted above, the strata associated with the respective vertices essentially represent a rigid translation, and it is apparent that the stratum associated with the vertex having coordinates (1/2,1/2) has the same size and shape as the stratum associated with the vertex having coordinates (0,0). The particular location within each stratum at which the vertex of the rotated lattice will be determined by the value of “x” (reference the second line of equation (10)).

[0031]
As another illustration of the generation of values for g, j
_{1}, j
_{1}, j
_{2 }and B, for the twodimensional 19 lattice associated with the ninth element of the Fibonacci sequence, F
_{9}, there are F
_{9}=F
_{8}+F
_{7}=21+13=34 vertices and, for the lattice that is associated with th

 (i) the vector generator g=(1,F_{k−1})=(1,F_{8})=(1,21);
$\left(\mathrm{ii}\right)\text{\hspace{1em}}{j}_{1}={F}_{2\xb7\lfloor \frac{91}{4}\rfloor +1}={F}_{2\xb7\lfloor \frac{8}{4}\rfloor +1}={F}_{2\xb72+1}={F}_{5}=5;\text{\hspace{1em}}$
$\left(\mathrm{iii}\right)\text{\hspace{1em}}{j}_{2}={F}_{2\xb7\lfloor \frac{9+1}{4}\rfloor}={F}_{2\xb7\lfloor \frac{10}{4}\rfloor}={F}_{2\xb72}={F}_{4}=3;\mathrm{and}\text{\hspace{1em}}\mathrm{so}\text{}\left(\mathrm{iv}\right)\text{\hspace{1em}}B=\frac{1}{34}\left(\begin{array}{cc}5& 3\\ 5\xb721& 3\xb721\end{array}\right)=\frac{1}{34}\left(\begin{array}{cc}5& 3\\ 105& 63\end{array}\right).\text{\hspace{1em}}$
It will be apparent that the term
$\frac{j}{n}\xb7g$
in the second line of equation (10), which defines the unrotated lattice, corresponds to vertices having coordinates (0,0) for j=0, (1/34,21/34) for j=1,(2/34,8/34) for j=2, (3/34,29/34) for j=3, (4/34,16/34) for j=4, (5/34,3/34) for j=5, and so forth. The matrix B, in turn, defines for the vertex associated with j=0, namely, the vertex having coordinates (0,0), the stratum comprising the quadrilateral bounded by line segments between points (0,0), (5/34,3/34), (8/34,8/34) and (3/34,5/34), with the quadrilateral including the line segments between points (0,0)(5/34,3/34) and (0,0)(3/34,5/34) (except for points (5/34,3/34) and (3/34,5/34), and not including the line segments between points (5/34,3/34)(8/34,8/34) and (3/34,5/34)(8134,8/34). The strata that are associated with the other vertices can be determined in a similar manner. As with the strata associated with the twovertex lattice (F_{3}) described above, in the case of the thirtyfour vertex (F_{9}) lattice, the strata associated with all of the vertices have the same sizes and shapes

[0033]
Other lines of the table depicted in FIG. 2 depict similar information for twodimensional Fibonacci lattices having different numbers of vertices.

[0034]
Benefits of stratified sampling using rankone lattices will be illustrated using several examples. Generally, in contrast to stratification using gridaligned lattices, stratifying an “s”dimensional domain using rank1 lattices provides that the number of strata will not grow exponentially with the number of dimensions in the integration domain. This will be illustrated in connection with jittered sampling, correlating sampling and correlated trajectory sampling methodologies.

[0035]
Jittered Sampling by Rank1 Lattices

[0036]
Generally, jittered sampling is achieved by placing the samples points used for evaluating the integrand ƒ in random positions within respective strata. Using the stratification that is induced by rank1 lattices (equation (10)), the value of an integral can be estimated as follows:
$\begin{array}{cc}\begin{array}{c}{\int}_{{\left[0,1\right)}^{s}}f\left(x\right)dx=\sum _{j=0}^{n1}{\int}_{{A}_{j}}f\left(x\right)dx\\ =\frac{1}{n}\sum _{j=0}^{n1}{\int}_{{\left[0,1\right)}^{s}}f\left({R}_{j}\left(x\right)\right)dx\\ \approx \frac{1}{n}\sum _{j=0}^{n1}f\left({R}_{j}\left({\xi}_{j}\right)\right)\end{array}& \left(12\right)\end{array}$
In equation (12), the first line reflects the fact that the integral over the entire integration domain, represented by the left hand side, corresponds to the sum of the integrals over the various strata A_{j}, j=0, . . . ,n1. The second line of equation (12) follows from the fact that, since, if “x” is a continuous variable spanning the integration domain, R_{j}(x) will also be continuous and span the same entire integration domain. In that case, if the integral is summed “n” times, the value of the integral can be obtained by dividing the sum by “n.” Finally, the third line reflects that the value of the “jth” integral in the sum in the second line can be estimated by evaluating the function ƒ using the “jth” sample point, which corresponds to R_{j}(ξ_{j}).

[0037]
In determining the value for the estimate {overscore (ƒ)}, which corresponds to the last line of equation (12), the computer graphics system 10 will, for each value of index “j,” make use of the generator vector g to generate the
$\frac{j}{n}\xb7g$
term of R_{j}(ξ_{j}) (reference equation (10)), and the “jth” element ξ_{j }of an “s”dimensional random sequence
${\left\{{\xi}_{j}\right\}}_{j=0}^{n1}$
to generate the Bx term, where the “jth” element ξ_{j }is used as vector “x” in generating R_{j}(ξ_{j}) (equation (10)). The sum Rj(ξ_{j}) is used as the sample point in evaluating the function ƒ for the respective “jth” sample point. Since different elements ξ_{j }of the random sequence
${\left\{{\xi}_{j}\right\}}_{j=0}^{n1}$
are used in generating R_{j}(ξ_{j}) for different values of index “j,” the sample points ξ_{j }are likely to have different displacements within the various strata A_{j}. However since there is at least one sample point R_{j}(ξ_{j}) in each stratum, and since the strata are uniformly distributed over the integration domain, the sample points will be reasonably well distributed around the integration domain. It will be appreciated that several elements ξ_{j }of the random sequence
${\left\{{\xi}_{j}\right\}}_{j=0}^{n1}$
may clump, and may even have the same value, but the stratification induced by the
$\frac{j}{n}\xb7g$
term of equation (10) will ensure that the sample points R_{j}(ξ_{j}) that are used in evaluating the function ƒ are well distributed over the integration domain. A more uniform distribution may be achieved if the elements ξ_{j }of the random sequence
${\left\{{\xi}_{j}\right\}}_{j=0}^{n1}$
are generated using, for example, a strictlydeterministic, lowdiscrepancy methodology.

[0038]
Using rank1 lattices, the computer graphics system 10 can readily partition the sdimensional unit cube [0,1)^{s }into a number “n” of strata, with the number “n” being independent of the number “s” of dimensions. The estimator for the value of the integral, which is described in equation (12), and particularly in the last line of that equation, is unbiased and can also be used in connection with generating estimates of the values of integrals of functions related to fields other than computer graphics.

[0039]
Correlated Sampling by Rank1 Lattices

[0040]
As with jittered sampling, correlated sampling also involves evaluating a sum of integrals of a function ƒ over the sdimensional unit cube. However, unlike jittered sampling, in which elements of a random “s”dimensional sequence
${\left\{{\xi}_{j}\right\}}_{j=0}^{n1}$
is used in generating the sample points at which the function ƒ is evaluated (equation (12)), in correlated sampling one randomly generated element ξ is used in generating all of the sample points (equation (10)) that are used in evaluating the function ƒ for use in the sum. This is derived from equation (12) by, with reference to the second line of equation (12), interchanging the finite sum and the integral. That is,
$\begin{array}{cc}\begin{array}{c}\frac{1}{n}\sum _{j=0}^{n1}{\int}_{{\left[0,1\right)}^{s}}f\left({R}_{j}\left(x\right)\right)dx=\frac{1}{n}{\int}_{{\left[0,1\right)}^{s}}\sum _{j=0}^{n1}f\left({R}_{j}\left(x\right)\right)dx\\ \approx \frac{1}{n}\sum _{j=0}^{n1}f\left({R}_{j}\left(\xi \right)\right)\end{array}& \left(13\right)\end{array}$
It will be appreciated that, in equation (13), the left hand side of the first line of the equation corresponds to the right hand side of the second line of equation (12). In generating the estimate {overscore (ƒ)}, which corresponds to the last line of equation (13), the computer graphics system 10 will, for each value of “j,” make use of the generator vector g in generating the
$\frac{j}{n}\xb7g$
term of R_{j}(ξ), and the same randomlygenerated vector “ξ” as “x” in the Bx term of the sample point R_{j}(ξ). As noted above, in generating the image, the computer graphics system 10 will separately evaluate the integral for each of the pixels in the image; and in that case the computer graphics system 10 may use the same element ξ for all of the integrals associated with the various pixels, or, alternatively, the computer graphics system 10 may use different random elements . . . ξ_{i−1}, ξ_{i}, ξ_{1+1}, . . for the integral associated with the . . . “i1st,” “ith”, “i+1st” pixel, respectively. One advantage of using different elements for the various integrals is that that may result in the generation of noise in the image, which can obscure aliasing artifacts that might otherwise be present in an image.

[0041]
Since stratified samples are independent, there is no correlation between them; in particular, there is no minimum distance property. Therefore, no variance reduction will be guaranteed. Minimum distance constraints can be imposed by, for example, applying Lloyd's relaxation algorithm in a conventional manner. Combined with modified sample weights by poststratification, the estimate {overscore (ƒ)}converges to the actual value <ƒ> of the integral. For Riemannintegrable functions the unbiased estimator (reference the second line of equation (13)) has a reduced variance of order
$O\left(\frac{1}{{n}^{2}}{\mathrm{ln}}^{2s}n\right),$
if the set of sample points
${\left({R}_{j}\left(\xi \right)\right)}_{j=0}^{n1}$
is of sufficiently low discrepancy.

[0042]
The estimator (reference the second line of equation (13)) that is used for an integral can be realized by generating one random sdimensional vector ξ ε [0,1)^{s }and replicating it using equation (10), that is, using the product of B and the random element ξ as an offset for each of the unrotated lattice vertices, which, as noted above, are defined
$\frac{j}{n}\xb7g.$
The resulting correlated samples have a guaranteed minimum distance and good uniformity properties, if the generator vector g is suitably chosen. If values for “x” are uniformly distributed the transformation B can be omitted, in which case the bijection (reference equation (10)) reduces to
$\begin{array}{cc}{R}_{j}^{\mathrm{CP}}:{\left[0,1\right)}^{s}>{\left[0,1\right)}^{s}\text{}x\mapsto \left(\frac{j}{n}\xb7g+x\right)\text{\hspace{1em}}{\mathrm{mod}\text{\hspace{1em}}\left[0,1\right)}^{s}& \left(14\right)\end{array}$
which corresponds to a CranleyPatterson rotation on a rank1 lattice.

[0043]
Correlated Trajectory Splitting by Rank1 Lattices

[0044]
Trajectory splitting can improve efficiency depending on the correlation of the integrand with respect to its lower dimensional projections. Trajectory splitting is used in a number of types of simulations, including but not limited to computer graphics. Trajectory splitting can be used in several ways in computer graphics, one of which will be described here and another below. One use of trajectory splitting in connection with computer graphics is tracing multiple shadow rays to an light source that has a nonzero area. As noted above, in ray tracing, path tracing and similar methodologies, to determine a pixel value for a pixel in an image, the computer graphics system 10 traces a plurality of rays from respective points in the pixel into a simulated scene. For each such ray, if the ray intersects a surface of an object in the scene, the computer graphics system 10 can trace one or more rays, which are referred to as shadow rays, from the point of intersection of the shot ray to various points on the area of the light source. If a shadow ray intersects another object in the scene, the point of intersection of the shot ray is in the other object's shadow for that shadow ray.

[0045]
On the other hand, if the shadow ray does not intersect another object in the scene, the light source is deemed to contribute to the direct illumination of the point in the scene for the pixel from which the ray was shot. The shadow rays that are traced from the point of intersection to the light source are distributed over the surface of the light source, and so some of the shadow rays may intersect other objects in the scene, whereas others of the shadow rays may not. The extent to which the point of intersection of the shot ray is illuminated by the light source is a function of the number of shadow rays from that point on the object that do not intersect other objects. Generally, if a scene contains multiple light sources, similar operations will be performed for each light source to determine the extent to which the point at which the ray that was shot from the pixel intersects the object, is illuminated by the various light sources.

[0046]
As noted above, a pixel also has an area, and the computer graphics system will typically shoot simulated rays into the scene from several points within the area of the pixel, and trace shadow rays from points at which they intersect the various objects. Generally, the variance of the integrand with respect to the area of the pixel is less than the variance with respect to the area of the light source.

[0047]
In trajectory splitting, the samples (illustratively, the shadow rays) that are split off of one instance (illustratively, the ray that is traced from the pixel into the scene) do not need to be independent. Continuing with the computer graphics shadow ray example, using a vector “x” to refer to a location in the pixel, and a vector “y” to refer to a location on the light source, the function ƒ in the integral is a function of both vector “x” and vector “y,” that is f(x,y), with vector “x” being associated with s_{1 }dimensions and vector “y” being associated with s_{2 }dimensions. Accordingly, the portion of the integration domain that is associated with vector “x” is the s_{1}dimensional unit cube [0,1)^{s} ^{ 1 }and the portion of the integration domain that is associated with the vector “y” is the s_{2 }dimensional unit cube [0,1)_{2} ^{ 1 }. In that case, if s_{1}+s_{2}=s, which will be the case if the vectors “x” and “y” are associated with all of the dimensions comprising the integration domain, the integral of ƒ is evaluated as follows
$\begin{array}{cc}\begin{array}{c}{\int}_{{\left[0,1\right)}^{{s}_{1}}}{\int}_{{\left[0,1\right)}^{{s}_{2}}}f\left(x,y\right)dydx={\int}_{{\left[0,1\right)}^{{s}_{1}}}\sum _{j=0}^{n1}{\int}_{{A}_{j}}f\left(x,y\right)dy\text{\hspace{1em}}dx\\ ={\int}_{{\left[0,1\right)}^{{s}_{1}}}\frac{1}{n}{\int}_{{\left[0,1\right)}^{{s}_{2}}}\sum _{j=0}^{n1}f\left(x,{R}_{j}\left(y\right)\right)dy\text{\hspace{1em}}dx\\ \approx \frac{1}{r}\sum _{i=0}^{r1}\frac{1}{n}\sum _{j=0}^{n1}f\left({\xi}_{i},{R}_{j}\left({\zeta}_{i}\right)\right)\end{array}& \left(15\right)\end{array}$
where ξ_{i }and ζ_{i }are the “ith” elements of two random sequences, namely, an s_{1}dimensional random sequence
${\left\{{\xi}_{i}\right\}}_{i=0}^{r1}$
and an s_{2}dimensional random sequence
${\left\{{\zeta}_{i}\right\}}_{i=0}^{r1}.$
Equation (15) essentially reflects the fact that the rank1 lattice, instead of being used to stratify the entire “s”dimensional integration domain, can be used to stratify a subset s_{2}<s of the dimensions comprising the integration domain. Continuing with the shadowray example, the element ξ_{i }is used to determine the location of a point within the pixel from which the simulated ray is to be shot into the scene. On the other hand, R_{j}(ζ_{i}) (reference equation (10)) is used to determine the location on the area light source toward which a shadow ray will be traced from the point at which the ray shot from the pixel intersects an object in the scene.

[0048]
It will be apparent that, in evaluating the integral according to the last line of equation (15), the sampling methodology essentially corresponds to subdividing the area light source into a number j=0, . . . ,n1 of light sources that are induced by the rank1 lattice stratification (reference equation (10)). For each of the “r” rays that are shot into the scene from the point in the pixel that is associated with the “ith” element ξ_{i }of the random sequence
${\left\{{\xi}_{i}\right\}}_{i=0}^{r1},$
if the shot ray intersects an object in the scene, one shadow ray is traced from the point at which the shot ray intersects an object in the scene, to one point within each of the j=1, . . . ,n strata on the light source, the point in the stratum that is determined by the “ith” element of the random sequence
${\left\{{\zeta}_{i}\right\}}_{i=0}^{r1}.$
Since for each element ξ_{i }of the random sequence
${\left\{{\xi}_{i}\right\}}_{i=0}^{r1}$
(which, as noted above, is used in determining the location within the pixel from which the ray shot will be shot into the scene), the same element ζ_{i }of the random sequence
${\left\{{\zeta}_{i}\right\}}_{i=0}^{r1}$
will be used in generating values for R_{j}(ζ_{i}) (which is used in determining the locations in the various strata on the light source towards which the “j” respective shadow rays will be traced) for various values of “j,” the shadow rays that are traced from the point of incidence of the shot ray will have the same displacements in the various strata that are induced on the light source.

[0049]
On the other hand, for the “r” different points in the pixel from which rays are shot, locations 10 of which are determined using the different elements . . . ,ξ_{i−1}, ξ_{i}, ξ_{i+1 }. . . of the random sequence
${\left\{{\xi}_{i}\right\}}_{i=0}^{r1},$
a different element . . . ,ζ_{i−1}, ζ_{i}, ζ_{i+1 }. . . of the random sequence
${\left\{{\zeta}_{i}\right\}}_{i=0}^{r1}$
will be used in determining the displacements within the various strata on the light source toward which the shadow rays will be traced. Accordingly, the displacements in the respective strata on the light source towards which shadow rays will be shot will usually differ for the different rays that are shot in the scene from the respective pixel. It will be appreciated, however, that, since the elements ζ_{i }are from a random “s_{2}”dimensional sequence, it is possible that two more elements of the sequence
${\left\{{\zeta}_{i}\right\}}_{i=0}^{r1}$
may have the same or close to the same value, in which case, shadow rays that are traced from the intersection points in the scene will be traced to the same or almost the same points on the light source.

[0050]
The jittered sampling, correlated sampling and trajectory splitting methodologies as described above will be further described in connection with three particular applications that are used in connection with generating images in computer graphics, namely, generating photon maps, path tracing and distribution ray tracing.

[0051]
Photon Map Generation

[0052]
As noted above, the photon map methodology generally includes two phases, namely, a first “photon map generation” phase followed by a “gather” phase. In the photon map generation phase, the computer graphics system 10 shoots rays that represent paths taken by simulated photons emitted by one or more simulated light sources into a simulated scene. Each rays is in the form of a random walk that is traced through at least a portion of the simulated scene. The computer graphics system determines whether a ray intersects a surface of an object in the scene, and, if so, stores information relating to the photon associated with the ray in a database. The information which may include, for example, the photon's color, its direction of incidence, the location of the intersection point, and other information that will be apparent to those skilled in the art, in a database. Depending on a number of criteria, including characteristics of the surface of the simulated object at the point at which the photon intersects the object, the number of times (if any) the photon has been previously reflected, and other criteria that will be apparent to those skilled in the art, the computer graphics system may terminate the photon at the point at which it intersects the object, or it may allow the photon to be reflected along another ray. If the computer graphics system determines that the photon is to be reflected, the computer graphics system determines the direction of reflection using a number of criteria, including, among other things, the angle of incidence and characteristics of the simulated object's surface, such as whether the surface is specular, glossy or diffuse. The computer graphics system can make use of random factors in several ways, including, for example, determining whether to terminate a ray at a surface that the ray intersects, or to allow the ray it to be reflected. The computer graphics system may also use random factors to jitter the angle with which a ray is to reflected, and as well as in other ways as will be apparent to those skilled in the art.

[0053]
After the computer graphics system has shot a number of rays into the scene, it initiates the second “gather” phase of the photon map methodology. In that phase, as described above, the computer graphics system 10 shoots simulated rays from a simulated image plane of a simulated camera into the scene. The rays that are shot from the image plane into the scene during the gather phase essentially correspond to the directions into the scene that are visible from the respective points on the pixel from which they are shot. If a gather ray intersects the surface of an object in the scene, the computer graphics system makes use of the information regarding the photons, if any, that intersected the surface proximate the point at which the ray shot from the camera intersected the surface, including the directions of incidence of the respective photons, the direction of incidence of the gather ray, and the various characteristics of the surface on which the gather ray is incident. The computer graphics system uses that information to determine the brightness and color of the scene at that point in the image.

[0054]
As noted above, for rays that are shot into the simulated scene, random factors can be used in several ways, including determining whether to terminate a ray at a surface that the ray intersects or whether to allow the ray to be reflected, determining the direction of reflection of the photons that are shot into the scene during the photon map generation phase, and other ways that will be appreciated by those skilled in the art. The random factors that are used to jitter each ray are collected into an “s”dimensional vector ξ ε[0,1)^{s}, where “s” is related to the maximum number of reflections or refractions that a photon may have before it is terminated. Problems have arisen in the past, since stratification of highdimensional integration domains has not been available. However, stratification using rank1 lattices provides a solution to this problem. In particular, the computer graphics system 10 makes use of a sequence of random sdimensional vectors
${\left\{{\xi}_{j}\right\}}_{j=0}^{r1},$
ξ_{j }ε[0,1)^{s}, and uses the “jth” element ξ_{j }of the sequence in generating R_{j}(ξ_{j}) (reference equation (10)) for the “jth” random walk. By jittering the random walk using a random number set, the computer graphics system 10 can reduce the likelihood that correlation patterns will develop in the traces that the computer graphics system shoots into the scene.

[0055]
Path Tracing

[0056]
In path tracing, the computer graphics system 10 evaluates a rendering equation that is of the form
L({right arrow over (x)},{right arrow over (ω)})+∫_{A} L({right arrow over (x)}′,{right arrow over (ω)}({right arrow over (x)}′,{right arrow over (x)})ƒ_{r}({right arrow over (x)},{right arrow over (ω)}){right arrow over (x)}′,{right arrow over (x)}))G({right arrow over (x)},{right arrow over (x)}ω′)V({right arrow over (x)},{right arrow over (x)}′)dA′ (16)
where

[0057]
(i) L({right arrow over (x)},{right arrow over (ω)}) on the lefthand side of equation (16) represents the luminance at point {right arrow over (x)} in direction {right arrow over (ω)};

[0058]
(ii) L_{e}({right arrow over (x)}, {right arrow over (ω)}) on the righthand side is an emission term representing the illumination emitted from point {right arrow over (x)} in the direction {right arrow over (ω)};

[0059]
(iii) L({right arrow over (x)}′, ({right arrow over (x)}′,{right arrow over (x)})) on the righthand side represents the luminance at another point {right arrow over (x)}′ in the scene that is reflected in the direction ω({right arrow over (x)}′,{right arrow over (x)}) toward point {right arrow over (x)};

[0060]
(iv) ƒ_{r}({right arrow over (x)},ω({right arrow over (x)}′,{right arrow over (x)})) on the righthand side is a bidirectional scattering distribution function that describes how much of the light incident on from the point {right arrow over (x)}′ is reflected, refracted or otherwise scattered in the direction ω({right arrow over (x)}′,{right arrow over (x)}) toward the point {right arrow over (x)}, and is generally the sum of a diffuse component, a glossy component and a specular component;

[0061]
(v) G({right arrow over (x)},{right arrow over (x)}′) on the righthand side is a geometric term that is a function of the orientation of the surfaces that contain points {right arrow over (x)}′ and {right arrow over (x)} relative to one another, specifically,
$\begin{array}{cc}G\left(\stackrel{>}{x},{\stackrel{>}{x}}^{\prime}\right)=\frac{\mathrm{cos}\text{\hspace{1em}}\theta \text{\hspace{1em}}\mathrm{cos}\text{\hspace{1em}}{\theta}^{\prime}}{{\uf603\stackrel{>}{x}{\stackrel{>}{x}}^{\prime}\uf604}^{2}}& \left(17\right)\end{array}$
where θ and θ′ are angles of the line between points {right arrow over (x)} and {right arrow over (x)}′, respectively relative to the normals of the surfaces that contain points {right arrow over (x)} and {right arrow over (x)}′, and (vi) V({right arrow over (x)},{right arrow over (x)}′) on the righthand side of equation (16) is a visibility function that equals the value one if the point {right arrow over (x)}′ is visible from the point {right arrow over (x)} and zero if the point {right arrow over (x)}′ is not visible from the point {right arrow over (x)}.
It will be appreciated that point {right arrow over (x)} represents the point at which the luminance is to be determined and point {right arrow over (x)}′ represents another point in the scene from which light may be emitted (in the case of a light source) or reflected toward the point {right arrow over (x)}. Essentially equation (16) reflects the fact that the luminance at point {right arrow over (x)} in the direction {right arrow over (ω)} is generally the sum of two components. One component, which is represented by item (ii) above (L_{e}({right arrow over (x)},{right arrow over (ω)})), represents the light (if any) that is emitted from the point {right arrow over (x)} in the direction {right arrow over (ω)}. This component will be nonzero if the point {right arrow over (x)} is on a light source.

[0062]
The second component, which is represented by the integral on the right hand side of equation (16), represents the light (if any) that is directed toward point {right arrow over (x)} from elsewhere in the scene, including, for example, light sources and other surfaces, and that is directed to or reflected or otherwise scattered towards the point {right arrow over (x)} from other points {right arrow over (x)}′ in the scene that are visible from point {right arrow over (x)}. The integral may be nonzero for points on light sources as well as for points on other objects in the scene. The integral in equation (16) is evaluated over the area A′ of the surface of a hemisphere that is formed from a sphere that is centered on the point {right arrow over (x)} and generally bisected by a plane that is tangent to the surface at point {right arrow over (x)}.

[0063]
The computer graphics system 10 numerically evaluates the integral in equation (16) by shooting a plurality of rays from respective points on the pixel into the scene. For each ray shot from the pixel, if the path of the ray intersects the surface of an simulated object in the scene, in addition to shooting shadow rays to light sources as described above to determine the extent to the point on the object is illuminated by a light source, the computer graphics system 10 can also determine the extent to which the point on the object is illuminated by light reflected from elsewhere in the scene. In the latter operation, the computer graphics system 10 can extend the ray by tracing a path from the point on the object that the path of the ray as shot from the pixel intersected, in a direction that is generally a reflection of the path of incidence. The computer graphics system can repeat these operations through a plurality of iterations, or until the path of the ray does not intersect an object in the scene. In addition, the computer graphics system 10 will perform similar operations for all of the rays that are shot from the pixel, and the combine the results for the respective rays to determine the pixel value for the pixel.

[0064]
It will be appreciated that the computer graphics system 10 makes use of a large number of paths in determining the pixel value, it can use correlated sampling to determine the locations of the points of the pixel. Since in correlated sampling the same value ξ is used for all of the sample points for the respective pixel, the computer graphics system 10 can generate the sample points R_{j}(ξ) for the pixel by use of a rigid transformation, which is more efficient than, for example, jittered sampling.

[0065]
Distribution Ray Tracing

[0066]
The distribution ray tracing methodology extends the path tracing methodology as described above, by adding trajectory splitting. In the distribution ray tracing methodology, if a ray shot from the simulated image plane into the scene intersects the simulated surface of a simulated object in a scene, instead of using a single ray reflected from the surface, the computer graphics system reflects a plurality of rays along a like plurality of paths from the intersection point on the surface. The computer graphics system divides the hemisphere into a plurality of strata, one for each path that is to be reflected from the surface. The computer graphics system directs a ray through each stratum, and can make use of trajectory splitting as described above to determine the specific direction through the respective strata. Accordingly, the computer graphics system will trace a plurality of paths from the point at which the ray's path intersects the point on the object. The computer graphics system can perform trajectory splitting at any point at which a ray's path intersects a point on an object.

[0067]
The invention provides a number of advantages. In particular, the invention provides a new and improved system and computerimplemented method for evaluating integrals using a Monte Carlo methodology in which the integration domain is partitioned into a plurality of strata using a stratification methodology in which the number of strata is independent of the number of dimensions comprising the integration domain. Specifically, the invention provides a methodology that makes use of stratification based on rank1 lattices. Systems and computerimplemented methods according to the invention find utility in a number of applications, including but not limited to computer graphics.

[0068]
The invention is also beneficial in connection with one problem that can arise in computer graphics, namely, aliasing, which gives rise to artifacts in the form of, for example, jagged edges. A computer graphics system 10 in accordance with the invention does not eliminate aliasing artifacts, but it does, however, mask them by noise that can arise from randomization. To accomplish this, the sample points are preferably random as between different pixels; however, the samples can be correlated within a pixel. The estimator described above in connection with the last line of equation (13) matches these requirements. Within a pixel, the samples are a shifted lattice, with the amount of shift being determined by the value of the random element ξ that is used for the pixel, providing good uniformity of distribution throughout the area of the pixel, which, in turn, provides for fast convergence. For different pixels, the various values of elements . . . ξ_{i−1},ξ_{i},ξ_{i+1 }. . . that are used provides noise that can mask the aliasing artifacts in the image that is generated.

[0069]
It will be appreciated that a number of changes and modifications may be made to the invention as described herein. For example, although the invention has been described in the context of a computer graphics system, it will be appreciated that the invention will find utility in other applications in which integrals are numerically evaluated using the Monte Carlo methodology.

[0070]
In addition, although the invention has been described in connection with a rank1 lattice in the form of a lattice based on the Fibonacci sequence, it will be appreciated that other forms of rank1 lattices may be used.

[0071]
It will be appreciated that a system in accordance with the invention can be constructed in whole or in part from special purpose hardware or a general purpose computer system, or any combination thereof, any portion of which may be controlled by a suitable program. Any program may in whole or in part comprise part of or be stored on the system in a conventional manner, or it may in whole or in part be provided in to the system over a network or other mechanism for transferring information in a conventional manner. In addition, it will be appreciated that the system may be operated and/or otherwise controlled by means of information provided by an operator using operator input elements (not shown) which may be connected directly to the system or which may transfer the information to the system over a network or other mechanism for transferring information in a conventional manner.

[0072]
The foregoing description has been limited to a specific embodiment of this invention. It will be apparent, however, that various variations and modifications may be made to the invention, with the attainment of some or all of the advantages of the invention. It is the object of the appended claims to cover these and such other variations and modifications as come within the true spirit and scope of the invention.