JP4947394B2 - Simultaneous simulation of Markov chains using the quasi-Monte Carlo method - Google Patents

Description

Cross-reference of related applications

  This patent application (Docket MENT-101-B) claims the benefit of the priority of US Provisional Patent Application No. 60 / 822,417, filed August 15, 2006. This application is also a continuation-in-part of US Patent Application No. 11 / 474,517 (reference number MENT-101-US) filed on June 23, 2006. Claims the benefit of US Provisional Patent Application No. 60 / 693,231 (Docket No. MENT-101-PR) filed on Nov. 19, 2002, and also filed on Nov. 19, 2002. It is also a partial continuation application of No. 574 (reference number MENT-075).

  US patent application Ser. No. 10 / 299,574 (Docket No. MENT-075) is a continuation of US patent application No. 09 / 884,861 (Docket No. MENT-061) filed Jun. 19, 2001. This application is a U.S. Provisional Patent Application No. 60 / 265,934 filed on Feb. 1, 2001, and U.S. Provisional Patent Application No. 60/212, filed Jun. 19, 2000. Insist on the benefits of 286 priority.

  Each of the patent applications listed above, including, but not limited to, docket numbers MENT-061, MENT-075, MENT-101-PR, and MENT-101-US The entirety of which is hereby incorporated by reference as if set forth in full herein.

Browse computer program list

  A source code list is attached to the present application as follows. The source code list is referred to herein as a “computer program list” and is organized into sections identified by a three-digit reference number of the form “1.1.1”.

Field of Invention

  The present invention relates generally to a method and system for image rendering in and by a digital computing system, such as for moving pictures and other applications, and in particular, for real-time high precision ray tracing. The present invention relates to a method, a system, a device, and computer software.

Background of the Invention

  The term “ray tracing” refers to a technique for synthesizing photorealistic images by identifying all the light paths that link a light source with a camera and summing these contributions. The simulation traces rays along the line of sight to determine visibility and traces rays from the light source to determine illuminance.

  Ray tracing has become the mainstream in video and other applications. However, current ray tracing techniques involve many known issues, including numerical problems, limitations in the ability to process dynamic scenes, slow acceleration data structure settings, and large memory footprints. It is affected by the limitations and weaknesses. Therefore, current ray tracing techniques lack the ability to efficiently process scenes with overall movement, such as forests that flutter in the wind and human hair. Overcoming the limitations of current raytracing systems also allows for the rendering of higher quality motion blur, for example in movie production.

  Current attempts to improve the performance of ray tracing systems are inadequate for a number of reasons. For example, current real-time ray tracing systems typically use a 3D tree based on an axis-aligned binary space partition as its acceleration structure. Since the main emphasis of these systems is on rendering static scenes, these systems are generally quite long, needed to build the necessary data structures for an overall moving scene The set time cannot be dealt with. Along this line, certain manufacturers have improved real-time ray tracing by developing techniques that can build efficient 3D trees and reduce the time required to traverse the tree. However, it can be shown that as the number of objects to be ray-traced increases, the expected memory requirement of the system increases secondarily.

  Another manufacturer has designed a ray tracing integrated circuit that uses a bounding volume hierarchy to improve system performance. However, the performance of the architecture has been found to be impaired if too many incoherent secondary rays are traced.

  In addition, attempts have been made to improve system performance by implementing 3D tree traversal techniques using field programmable gate arrays (FPGAs). The major increase in processing speed of these systems is obtained by taking advantage of the FPGA's ability to trace coherent beam bundles and perform fast hardwired computations.

  The construction of the acceleration structure has not yet been implemented in hardware. FPGA implementations generally use a floating point technique with reduced precision.

  A number of different techniques can be used to simulate photon trajectories involving Markov chains. The use of the quasi-Monte Carlo method in simulating a Markov chain is described, for example, by L'Ecuyer et al., Randomized Quasi-Monte Carlo Simulation of Markov Chain Sce. Ecuyer "), which is incorporated herein by reference as if set forth in its entirety.

  L'Ecuyer describes the use of a randomized quasi-Monte-Carlo method to estimate the state distribution at each step of a Markov chain with a total order (discrete or continuous) state space. . The number of steps in the chain can be random and unlimited. This method should be used to obtain a low-variance unbalanced estimator with an expected total cost up to some random downtime, especially when state-dependent costs are paid at each step Can do. The paper provides a numerical explanation, which has a large variance reduction compared to standard Monte Carlo.

  L'Ecuyer is a deterministic quasi-Monte Carlo method for estimating a transient measure over a fixed number of steps for discrete-time and discrete-state Markov chains with full ordered state space. method) is proposed and studied in the following paper, which is also incorporated herein by reference in its entirety:

  Lecot et al., “Quasi-Random Walk Methods”, Monte Carlo and Quasi-Monte Carlo Methods 2000, pp. 63-85, Springer-Verlag (Berlin, 2002)

  Lecot et al., “Quasi-Monte Carlo Methods for Estimating Transient Measures of Discrete Time Markov Chains”, Monte Carlo and Qualiz-Monte Carlo Me. 329-43, Springer-Verlag (Berlin, 2004) ("Lecot 2004")

L'Ecuyer is radix 2 (0,2) - using the column, which describes techniques to simulate the n = ^{2 k} pieces copies of the chain in parallel over the same number of steps. In step j of the chain, the technique reorders the n copies according to their state, and replaces the elements from nj to nj + n-1 in the (0,2) -sequence instead of a uniform random number for simulation. By using it, a transition (ie, the next state) is simulated for n copies. The technique assumes that the simulation of each transition in the chain requires a single uniform random variable. In Lecot 2004, convergence to the correct value was proved under certain conditions on the structure of the Markov chain transition matrix.

  L'Ecuyer has attempted to generalize the above technique to a Markov chain with a random and unlimited number of steps and having a continuous state space (such a technique is especially intended to handle regeneration simulations). In contrast, the number of uniform random variables d required to generate the next state in one step of a Markov chain can be greater than one.

  However, in order to simulate photon orbits, various researchers assume the use of Markov chains, but conventional techniques that use quasi-Monte Carlo methods to simulate Markov chains, for example, photon orbits. Provided only a slight improvement over random sampling. However, the QMC method according to the present invention described below can take advantage of the inherent low-dimensional structure of the solution operator, resulting in a much more efficient technique.

  Thus, prior to the present invention, there was still a need for more efficient methods, systems, apparatus, devices, and software program / code products for simultaneously simulating Markov chains.

  Accordingly, it is an object to provide an improved and more efficient method, system, apparatus, device, and software program / code product for simulating a Markov chain using quasi-Monte Carlo methods and sorting strategies. desirable.

  The present invention provides a method, system, apparatus, and computer software / computer code product suitable for use in a computer graphics system for generating pixel values of pixels in an image, wherein the pixel values are in a scene. The point is represented as recorded on the image plane of the simulated camera, and the computer graphics system simulates at least one ray emitted from the pixel toward the scene along the selected direction. The method is configured to generate a pixel value of the image using a selected ray-tracing method, the method comprising: calculating the intersection of the ray with the surface of an object in the scene; and Simulating the trajectory of rays that illuminate an object inside Look, step to simulate the track, using the Markov chain.

  In one aspect of the invention, the method, system, apparatus, and computer software / computer code product according to the invention simulate (and / or means for simulating) a Markov chain using quasi-Monte Carlo methods. And simulating the Markov chain includes sorting states, and the sorting step includes proximity sorting.

  Further aspects of the invention include simulating a plurality of Markov chains simultaneously, the step of simultaneously simulating a plurality of Markov chains comprising a Halton sequence, a Halton sequence modification, a lattice sequence, or (t, s) −. Utilizing a quasi-Monte Carlo point that can be selected using any of the columns. The step of simultaneously simulating multiple Markov chains can utilize any number of chains.

  In a further aspect of the invention, simulating multiple Markov chains simultaneously includes adding additional trajectories on the fly, such as by trajectory splitting. The simulating step can also include utilizing a technique for simulating particle absorption, and can include using a Russian Roulette technique.

In one aspect of the invention, simultaneously simulating multiple Markov chains for an s-dimensional problem includes the following technique.
n: = 0
X _{0, l} is initialized using the quasi-Monte Carlo point _xl .
loop:
-Sort the state vector using the appropriate order.
-N: = n + 1
-Continue the chain by calculating _{Xn, l} using the subsequent quasi-Monte Carlo points _xl .

  The present invention can also include ordering states by proximity in state space and randomization techniques. Any selected point set or point sequence with any selected randomization can be used in conjunction with the method of the present invention.

In another aspect, the present invention further includes the following techniques.
n: = 0, point sequence _xl is selected.
Instantiate randomization.
rand (x _l ) is used to initialize X _{0, l} .
Loop • Sort the state vector using some order.
・ N: = n + 1
• Instantiate randomization.
Continue the chain by calculating X _{n, l} using rand (x _l ).

  The randomization technique can include randomizing a quasi-Monte Carlo point at each time step n, and more specifically, the randomizing step is more specifically a radix b = 2 (from which subsequent samples are taken. selecting a t, s) -column and subjecting the sample to an XOR operation with an s-dimensional random vector, wherein a random vector is generated after each transition step.

  Another aspect of the invention uses a bucket enumerated by sorting by brightness or state norm, sorting using a spatial hierarchy, and at least one space filling curve. And sorting. Sorting using a spatial hierarchy can include utilizing a binary hierarchy, where the spatial hierarchy is either a BSP tree, a kD tree, or any other axis-parallel subset of a BSP tree structure, or a boundary volume. It can contain hierarchies or regular voxels. The binary hierarchy is recursively subdivided using the plane selected by the chosen heuristics, and then traversing the hierarchy in order to enumerate the leaves of the hierarchy in order of proximity Can be built. More efficient traversal can be provided by left-balancing the tree and storing the tree in an array data structure. Sorting using the spatial hierarchy may also include utilizing bucket sorting and a selected space filling curve and, in one aspect of the invention, using regular voxels.

In another aspect of the invention, the simulation and / or sorting process comprises
Delimiting an object to be rendered by an axis-aligned bounding box;
By applying a split plane that passes through a selected point of the object, the bounding box is split recursively or iteratively into successive left and right branches, each ending with a leaf, thereby Generating a leaf of
Ordering the leaves hierarchically according to their brightness;
Performing a bucket sort on a matrix of a selected size divided into a selected number of buckets;
Using a selected space-filling curve to traverse the matrix, wherein within each bucket a smaller space-filling curve is used to traverse within an individual cell within the matrix; and including.

  The simulation process of the present invention can be adapted for construction to solve integral equations. The simulation is operable to allow a light transport simulation to synthesize a real-looking image. According to this aspect of the invention, the integral equation underlying the light transport simulation can be reformulated as a path integral so that sampling the path space corresponds to simulating a Markov chain, In doing so, the path is established by ray tracing and scattering events, the initial distribution is determined by the emission characteristics of the light source or sensor, and the transition probability is a bi-directional reflection distribution function for the surface in the image. function) and bi-directional subsurface scattering distribution function. The transmission effect can be handled by utilizing path integrals, initial distributions, and transition probabilities. In a related aspect of the invention, the path integral can be solved by utilizing high-dimensional quasi-Monte Carlo points or by padding low-dimensional quasi-Monte Carlo points, and the integral equation underlying the light transport simulation is: It is regarded as a Fredholm integral equation or a Volterra integral equation.

  Additional details, embodiments, practices, and examples are set forth in the following detailed description, read in conjunction with the accompanying drawings, or will be readily apparent to those skilled in the art from reading the detailed description in conjunction with the drawings. Will.

1 is a schematic diagram of a conventional digital processing system in which the present invention may be deployed. 1 is a schematic diagram of a conventional personal computer or similar computing device in which the present invention may be deployed. FIG. 2 shows the overall method according to the first aspect of the present invention. FIG. 6 is a ray tracing procedure illustrating the self-intersection problem. FIG. 6 is an elevational view of a split axis parallel bounding box used as an acceleration data structure according to a further aspect of the invention. FIG. 6 is an isometric view of the axis parallel bounding box shown in FIG. 5 showing the partitioning of the bounding box using the L and R planes. FIG. 6 is an isometric view of the axis parallel bounding box shown in FIG. 5 showing the partitioning of the bounding box using the L and R planes. FIG. 6 is an isometric view of the axis parallel bounding box shown in FIG. 5 showing the partitioning of the bounding box using the L and R planes. 4 is a flow chart of a ray tracing method according to a further aspect of the present invention. 4 is a flow chart of a ray tracing method according to a further aspect of the present invention. It is a figure illustrating shift Φ _S (n). FIG. 5 illustrates an exemplary Halton column with its hierarchical characteristics. It is a figure which shows sampling transport route space by bidirectional route tracing. FIG. 6 is a rendered image showing the difference between a simulation using the Monte Carlo method and a simulation using a sequence randomized quasi-Monte Carlo method. FIG. 6 is a rendered image showing the difference between a simulation using the Monte Carlo method and a simulation using a sequence randomized quasi-Monte Carlo method. FIG. 6 is a rendered image showing the difference between a simulation using the Monte Carlo method and a simulation using a sequence randomized quasi-Monte Carlo method. FIG. 6 is a rendered image showing the difference between a simulation using the Monte Carlo method and a simulation using a sequence randomized quasi-Monte Carlo method. It is a figure which shows the photon orbit from the starting point on a light source. It is a figure which shows the photon orbit from the starting point on a light source. It is a figure which shows the photon orbit from the starting point on a light source. It is a figure which shows the photon orbit from the starting point on a light source. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. FIG. 4 shows the steps of defining the order of proximity by classifying states into leaves of a spatial hierarchy, thereby traversing the hierarchy in order. It is a figure of the test scene "Shirley 6". It is a graph which shows the average RMS error with respect to a master solution. It is a graph which shows global dispersion | distribution. 6 is a graph showing pixel-based dispersion. It is a figure of a test scene "Invisible Date". It is a graph which shows the average RMS error with respect to a master solution. It is a graph which shows global dispersion | distribution. 6 is a graph showing pixel-based dispersion. FIG. 6 is a visual comparison of a test scene “Invisible Date” rendered using Monte Carlo. FIG. 6 is a visual comparison of a test scene “Invisible Date” rendered using randomized quasi-Monte Carlo. FIG. 7 is a visual comparison of a test scene “Invisible Date” rendered using randomized quasi-Monte Carlo with sorting. It is the schematic of a labyrinth test scene. Figure 6 is a graph showing the average amount of total radiation received by the camera when using each technique. It is a figure of the enlarged part of a graph. It is a table | surface which displays the deviation from a master solution. 4 is a flowchart of techniques and sub-techniques according to the described aspects of the invention. 4 is a flowchart of techniques and sub-techniques according to the described aspects of the invention. 4 is a flowchart of techniques and sub-techniques according to the described aspects of the invention. 4 is a flowchart of techniques and sub-techniques according to the described aspects of the invention. 4 is a flowchart of techniques and sub-techniques according to the described aspects of the invention. 4 is a flowchart of techniques and sub-techniques according to the described aspects of the invention. 4 is a flowchart of techniques and sub-techniques according to the described aspects of the invention.

  The present invention provides an improved technique for the efficient construction of acceleration data structures required for ray tracing and for fast ray tracing. The following description describes methods, structures, and systems according to these techniques.

  The described methods and systems operate in accordance with a conventional operating system such as Microsoft Windows, Linux, or Unix, or emulate a conventional operating system, in a stand-alone configuration or distributed over a network (PC) Those skilled in the art will appreciate that conventional computer equipment, such as) or equivalent devices, may be implemented in software, hardware, or a combination of software and hardware. The various processing and computing means described below and recited in the claims are thus implemented in software and / or hardware elements of a suitably configured digital processing device or network of devices. Good. Processing may be performed sequentially or in parallel, and may be performed using dedicated or reconfigurable hardware.

  A method, device, or software product according to the present invention is illustrated by way of example in FIG. 1 (eg, network system 100), whether it is stand-alone, networked, portable, or fixed. Can operate on any of a wide range of conventional computing devices and systems, such as a conventional PC 102, laptop 104, handheld or mobile computer 106, or Connected via the Internet or other network 108, the Internet or other network 108 may similarly include a server 110 and storage 112.

  In line with conventional computer software and hardware implementations, a software application configured in accordance with the present invention can run in, for example, a PC 102 as shown in FIG. 2, where program instructions are stored in a CD-ROM 116. Can be read from a magnetic disk or other storage 120 and loaded into the RAM 114 for execution by the CPU 118. Data can be entered into the system via any known device or means, including a conventional keyboard, scanner, mouse, or other element 103.

This description of the invention is divided into two large sections.
I. Real-time high-accuracy ray tracing II. A sorting strategy for rays simulated using Markov chains and quasi-Monte Carlo methods

I. Real-Time Precision Ray Tracing FIG. 3 illustrates an overall method 200 according to one aspect of the present invention. The method is implemented in connection with a computer graphics system, where a pixel value is generated for each pixel in the image. Each generated pixel value represents a point in the scene as recorded on the image plane of the simulated camera. The computer graphics system is configured to generate pixel values for the image using the selected ray tracing method. The selected ray tracing method includes the use of a ray tree that includes at least one ray emitted from the pixel toward the scene along the selected direction, and includes the ray and the object in the scene (and / or Or the calculation of the intersection with the surface of the object).

  In the method 200 of FIG. 3, a boundary volume hierarchy is used to calculate the intersection of the ray with the surface in the scene. In step 201, the bounding box of the scene is calculated. In step 202, it is determined whether predetermined termination criteria are met. If not, in step 203 the axis parallel bounding box is refined. The process continues recursively until termination criteria are met. According to one aspect of the invention, the termination criterion is defined as the condition that the bounding box coordinates differ from the floating point representation of the ray / surface intersection by one unit of resolution. However, the scope of the present invention extends to other termination criteria.

  The use of a border volume hierarchy as an acceleration structure is advantageous for a number of reasons. Memory requirements for the boundary volume hierarchy can be linearly limited to the number of objects that are ray-traced. Also, as will be explained below, the boundary volume hierarchy can be built much more efficiently than a 3D tree, which is an amortized analysis required for an overall moving scene. analysis) makes the boundary volume hierarchy very suitable.

  The following description describes in detail certain problems in ray tracing methods and specific aspects of the present invention that address those problems. FIG. 4 is a diagram illustrating the “self-intersection” problem. FIG. 4 shows a ray tracing procedure 300 that includes an image surface 302, an observation point 304, and a light source 306. In order to synthesize the surface image, a series of calculations are performed to find rays extending between the observation point 304 and the surface 302. FIG. 4 shows one such ray 308. Ideally, in that case, the exact intersection 310 of the ray 308 and the surface 302 is calculated.

  However, due to floating point arithmetic calculations in the computer, the calculated ray / surface intersection 312 may sometimes be different from the actual intersection 310. Further, as illustrated in FIG. 4, the calculated point 312 can be found on the “wrong” side of the surface 302. In that case, when calculations are performed to find a second ray 314 that extends from the calculated ray / surface intersection 312 to the light source 306, these calculations result in the second ray 314 reaching the light source 306. Instead of extending directly, it shows hitting the surface 302 at the second intersection point 316, thus leading to an image error.

  One known solution to the self-intersection problem is to start each ray 308 from a safe distance from the surface 302. This safe distance is generally expressed as a global floating point ε. However, the determination of the global floating point ε is highly dependent on the scene and the specific location within the scene itself where the image is synthesized.

  One aspect of the present invention provides a more accurate alternative. After reaching the calculated ray / surface intersection 312, the direction of the calculated point 312 and ray 308 is used to recalculate the intersection closer to the actual intersection 310. This recalculation of intersection points is incorporated into the ray tracing method as an iteration that increases accuracy. If an iteratively calculated intersection is found to be on the “wrong” side of the surface 302, the intersection is moved to the “correct” side of the surface 302. The iteratively calculated intersection can move along the surface normal or along an axis determined by the longest component of the normal. Instead of using the global floating point ε, the decimal point is moved by the integer ε toward the last bit of the floating point mantissa.

  The described procedure has the advantage of avoiding double precision calculations and implicitly adapting to the floating point scale determined by the exponent of the floating point number. Thus, in this implementation, all secondary rays start directly from these modified points that eliminate the ε offset. Thus, during the intersection calculation, it can be assumed that the effective ray section starts at 0 instead of some offset.

  The modification of the integer representation of the mantissa also avoids numerical problems when crossing a triangle and a plane to determine which point is on which side.

  Utilizing the convex hull nature of convex coupling, the intersection of a ray and a free-form surface can be found by refining the axis-parallel bounding box containing the intersection closest to the source of the ray. This refinement can continue until a floating point resolution is reached, i.e., the bounding box coordinates differ from the floating point representation by one unit of resolution. In that case, the self-intersection problem is avoided by selecting the corner of the bounding box that is closest to the surface normal at the center of the bounding box. This corner point is then used to start the second ray. This “ray object intersection test” is very efficient and benefits from avoiding the self-intersection problem.

  After building the acceleration data structure, the triangle is deformed in-plane. The new representation encodes the degenerate triangle so that the intersection test can handle the degenerate triangle without any special effort. The new representation can of course simply prevent the degenerate triangle from entering the graphics pipeline. Sections 2.2.1 and 2.2.2 of the computer program listing show a list of source code according to this aspect of the invention.

  The test first determines the intersection of the ray with the triangular plane, then the effective interval on the ray] 0, result. tfar] is excluded. This is achieved by a single integer test. Note that +0 is excluded from the valid interval. This is important when denormalized floating point numbers are not supported. If this initial decision is successful, the test proceeds by calculating the barycentric coordinates of the intersection. Note that in order to perform a full inclusion test, again only an integer test is needed, i.e. more specifically only a test of two bits. Therefore, the number of branches is minimal. In order to allow this efficient testing, triangle edges and normals are appropriately scaled in the deformation step.

  The accuracy of the test is sufficient to avoid false or lost ray intersections. However, during traversal situations may arise where it is appropriate to enlarge the triangles for robust intersection testing. This can be done before deforming the triangle. Since the triangle is projected along the axis identified by the longest component of its normal, this projection case must be preserved. This is accomplished by a counter at the leaf node of the acceleration data structure. Triangle references are sorted by projection case, and the leaf contains bytes for the number of triangles in each class.

  A further aspect of the invention provides an improved approach for building an acceleration data structure for ray tracing. Compared to conventional software implementations following many different optimizations, the approach described herein yields a significantly flat tree that provides excellent ray tracing performance.

  The split plane candidates are given by the coordinates of the triangle vertices inside the axis-parallel bounding box to be split. Note that this includes vertices that are actually outside the bounding box, but have at least one coordinate contained in one of the three intervals defined by the bounding box. From these candidates, the plane closest to the center of the longest side of the current axis-parallel bounding box is selected. Further optimization selects only those triangle candidates whose longest surface normal component matches the normal of the potential splitting plane. This procedure yields a much flatter tree because placing the splitting plane through the triangle vertices implicitly reduces the number of triangles split by the splitting plane. In addition, the surface is closely approximated and the empty space is maximized. If the number of triangles exceeds the specified threshold and there are no more split plane candidates, the box is split in the middle along its longest side. This avoids the inefficiencies of other approaches, including, for example, the use of long diagonal objects.

  The recursive procedure for determining which triangles belong to the left and right children of a node in the hierarchy generally required extensive bookkeeping and memory allocation. There is a much simpler approach that does not work only for exception cases. Only two arrays of references to the ray-tracing object are assigned. The first array is initialized with an object reference. During the recursive spatial partitioning, the left element stack grows from the beginning of the array, while the elements classified on the right are maintained on the stack growing from the end of the array toward the center. The second array maintains their stack so that the elements that intersect the split plane, i.e., the elements that are both left and right, can be quickly recovered. Thus, backtracking is efficient and simple.

  Instead of pruning the branches of the tree by using a surface area heuristic, the depth of the tree is pruned by approximately left equilibrating binary spatial partitioning starting from a fixed depth. As observed by exhaustive experiments, global fixed depth parameters can be specified across a wide variety of scenes. This can be understood by observing that a relatively flat connected element usually remains in space after a certain number of binary space divisions. Section 2.3.1 of the computer program listing shows a list of source code according to this aspect of the invention.

  Using the boundary volume hierarchy, each ray-tracing object is referenced exactly once. As a result, in contrast to 3D trees, a mailbox mechanism is not required to prevent multiple intersections of objects and rays during traversal of the hierarchy. This is a great advantage from a system performance perspective and makes it much easier to implement on a shared memory system. The second important result is that there can be no more internal nodes in the boundary volume hierarchy tree than the total number of ray-traced objects. Thus, the memory footprint of the acceleration data structure can be linearly limited to the number of objects prior to construction. Such pre-limits are not available in the construction of 3D trees, where it is expected that memory complexity will increase quadratically with the number of objects that are ray-traced.

  Thus, there is a new concept of boundary volume hierarchies that is significantly faster than current 3D tree ray tracing methods, and that memory requirements increase linearly rather than secondarily expected with the number of objects to be ray-traced. Will be explained from now on. The core concept that allows the boundary volume hierarchy to outperform the 3D tree is to focus on how the space can be divided instead of focusing on the boundary volume itself.

  In a 3D tree, the bounding box is divided by a single plane. In accordance with this aspect of the invention, two parallel planes are used to define two axis parallel bounding boxes. FIG. 5 is a diagram showing a main data structure 400.

  FIG. 5 shows the axis parallel bounding box 402 in elevation. The L plane 402 and the R plane 404 that are axially parallel and parallel to each other are used to divide the bounding box 402 into a left-axis parallel bounding box and a right-axis parallel bounding box. The left bounding box extends from the left wall 406 of the original bounding box 402 to the L plane 402. The right bounding box extends from the R plane 404 to the right wall 408 of the original bounding box 402. Therefore, the left bounding box and the right bounding box may overlap each other. The traversal of ray 412 is determined by the position of intersection of L plane 402 and R plane 404 with respect to the valid section [N, F] 412 of ray 410.

  In the data structure 400 of FIG. 5, the L plane 402 and the R plane 404 are positioned relative to each other so as to divide the set of objects contained within the original bounding box 402, rather than the space contained within the bounding box 402. It is done. In contrast to 3D tree partitioning, having two planes offers the possibility of maximizing the empty space between the two planes. As a result, the scene boundaries can be approximated much faster.

  FIGS. 6-8 are a series of three-dimensional diagrams that further illustrate the data structure 400. FIG. 6 shows a diagram of the bounding box 402. For illustrative purposes, a virtual object is shown in the bounding box 402 as an abstract circle 416. As shown in FIGS. 7 and 8, the L plane 404 and the R plane 406 are then used to divide the bounding box 402 into a left bounding box 402a and a right bounding box 402b. The L and R planes are chosen so that the empty space between them is maximized. Each virtual object 416 eventually enters the left bounding box 402a or the right bounding box 402b. As shown in the lower part of FIG. 8, the virtual object 416 is divided into a “left” object 416a and a “right” object 416b. Each of the resulting bounding boxes 402a, 402b is itself split and so on until the exit criteria are met.

  FIG. 9 is a flowchart of the method 500 described. In step 501, the bounding box of the scene is calculated. In step 502, the L and R planes are used to divide the axis parallel bounding box into left and right axis parallel bounding boxes that may overlap. In step 503, the left and right bounding boxes are used to divide the set of virtual objects contained in the original axis parallel bounding box into a set of left objects and a set of right objects. In step 504, the left and right objects are processed recursively until the exit criteria are met.

  Instead of the one splitting parameter used in the previous implementation, two splitting parameters are stored in the node. Since the number of nodes is linearly limited by the number of ray-tracing objects, the array of all nodes can be assigned only once. This eliminates the need for expensive memory management of the 3D tree under construction.

  The construction technique is much simpler than the analog for 3D tree construction and is easily implemented in a recursive manner or by using iterative versions and stacks. Given a list of objects and an axis-parallel bounding box, the L and R planes are determined and the set of objects is determined accordingly. The left and right objects are then processed recursively until some termination criterion is met. Since the number of internal nodes is limited, it is safe to rely on termination when only one object remains.

  Note that segmentation relies solely on classifying objects along planes perpendicular to the x, y, and z axes, which is very efficient and numerically absolutely stable. I want. In contrast to the 3D tree, the exact intersection of the object and the splitting plane does not need to be calculated, which is more expensive and difficult to achieve in a numerically robust manner. Numerical problems with 3D trees, such as the loss of triangles along vertices and edges, can be avoided by expanding the triangles before building the boundary volume hierarchy. Also, in a 3D tree, overlapping objects must be classified into both a left-axis parallel bounding box and a right-axis parallel bounding box, thereby causing the expected secondary growth of the tree.

  Various techniques may be used to determine the L and R planes, and thus the actual tree layout. Returning to FIGS. 6-8, one technique determines the plane M418 using the 3D tree construction technique described above, and the resulting overlap of the L and R planes of the new axis parallel bounding box has been proposed. The object is divided so as to overlap with the dividing plane M418 the smallest. The resulting tree is very similar to the corresponding 3D tree, but the resulting tree is much flatter because the object set is split instead of space. Another approach is to select the R and L planes in such a way that child box overlap is minimized and the empty space is maximized as much as possible. Note that for some objects, the axis parallel bounding box is inefficient. One example of such a situation is a long cylinder with a small radius on the diagonal of an axis-parallel bounding box.

  FIG. 10 is a flowchart of a method 600 according to this aspect of the invention. In step 601, the bounding box of the scene is calculated. In step 602, a 3D tree construction is performed to determine the split plane M. In step 603, the parallel L and R planes are used to divide the axis parallel bounding box into a left axis parallel bounding box and a right axis parallel bounding box that overlap the partition plane M the smallest. In step 604, the left and right bounding boxes are used to divide the set of virtual objects contained in the original axis parallel bounding box into a set of left objects and a set of right objects. In step 605, the left and right objects are processed recursively until the exit criteria are met. The method 600 described in FIG. 10 is similar to the method 200 described in FIG. 3 with other techniques described herein, including techniques related to 3D tree construction, real-time processing, bucket sorting, self-intersection, and the like. Note that may be combined with:

  In the case of a 3D tree, spatial subdivision is continued to clip away the empty portion of the space around the object. For the boundary volume hierarchy described, splitting such objects into smaller objects yields similar behavior. In order to maintain predictability of memory requirements, a maximum bounding box size is defined. All objects with a size that exceeds the maximum bounding box size are divided into smaller parts to meet the requirement. The maximum allowable size can be found by scanning the data set for the smallest size of all objects.

  The data structures described herein allow fast 3D tree traversal principles to be ported to the boundary volume hierarchy. (1) Left child only, (2) Right child only, (3) After left child, right child, (4) After right child, left child, or (5) Rays between split planes (ie empty The case of crossing that is in the space is the same. Since one node in the described technique is divided by two parallel planes, the order in which the box is traversed is determined by the direction of the rays. FIG. 6 illustrates a source code listing that includes the techniques described above.

  Previous boundary volume hierarchy techniques have not been able to efficiently determine in what order to traverse additional work required such as updating child nodes or heap data structures. In addition, although the entire bounding volume had to be loaded and tested for rays, this approach only requires two plane distances. However, checking the rays against two planes in software seems more expensive. Traversal is a bottleneck in the 3D tree, and doing a little more computation here better hides memory access latency. In addition, the boundary volume hierarchy tree tends to be much smaller than the corresponding 3D tree with the same performance.

  Although a new bounding volume hierarchy is described herein, there is a strong association with 3D tree traversal. Setting L = R acquires traditional binary space partitioning, and the traversal technique collapses into a 3D tree traversal technique.

  The described boundary volume hierarchy can also be applied to efficiently find ray free surface intersections by subdividing the free surface. Doing so allows for the intersection of a free-form surface with convex hullability, and a subdivision technique that is efficiently computed to floating-point precision depending on the actual floating-point computation. The subdivision step is performed on, for example, a polynomial curved surface, a rational curved surface, and an approximate subdivision surface. For each axis in space, bounding boxes that may overlap are determined as described above. For binary subdivisions, L box intersection and RL box intersection for a new bounding box of the new mesh. Now that the spatial order of the boxes is known, the above traversal can be performed efficiently. Instead of precomputing the boundary volume hierarchy, it can be computed on the fly. This procedure is efficient for freeform surfaces and allows saving memory for acceleration data structures, which is replaced by a small stack of boundary volumes that must be traversed by backtracking. The subdivision continues until the ray surface intersection enters a boundary volume that collapses into a point of floating point precision or a small size interval. Section 2.1.1 of the computer program listing shows a code listing according to this aspect of the invention.

Although the use of regular grids as an acceleration data structure in ray tracing is straightforward, efficiency is affected by a lack of spatial adaptability and subsequent traversal of many empty grid cells. Hierarchical regular grids can improve the situation, but are still inferior compared to boundary volume hierarchies and 3D trees. However, regular grids can be used to improve the speed of building acceleration data structures. The technique for building the acceleration data structure is similar to quick sorting and is expected to work with ◯ (nlogn). Improvements can be obtained by applying bucket sorting that operates in linear time. Accordingly, the object axis-aligned bounding box is divided into n _x × n _y × n _z number of axis-aligned box. Each object is then accurately classified into one of these boxes by one selected point, for example, the centroid or the first vertex of each triangle can be used. Thereafter, the actual axial parallel bounding box of the object in each grid cell is determined. These axis-parallel bounding boxes are used in place of the objects they contain unless the box intersects one of the split planes. If they intersect, the box is opened and instead the object in the box is used directly. This procedure can also be applied recursively, eliminating many comparisons and memory accesses, significantly improving the constants of the order of construction techniques. The above technique is particularly attractive for hardware implementations because it can be realized by processing a stream of objects.

  An acceleration data structure can be constructed on demand, that is, when a ray traverses a particular axis-parallel bounding box with an object. In that case, on the other hand, the acceleration data structure will not be elaborated in areas of light that do not shine, and the cache will not be polluted by data that will never be touched. On the other hand, after elaboration, objects that can intersect the light beam already exist in the cache.

  From the above description, the present invention addresses a long known problem in ray tracing and improves the accuracy, overall speed, and memory footprint of acceleration data structures, and techniques for ray tracing. You can see that The improvement in numerical accuracy is not only a conversion to other numerical systems, but also a conversion to a logarithmic system used in, for example, the hardware of an ART ray tracing chip. Note that the specific implementation of the IEEE floating point standard on a processor or dedicated hardware can have a severe impact on performance. For example, on a Pentium 4 chip, a denormalized number can degrade performance to 1/100 or less. As explained above, implementation of the present invention avoids these exceptions. The boundary volume hierarchy concept described herein makes the boundary volume hierarchy suitable for real-time ray tracing. In amortization analysis, the described technique outperforms previous situations in the art, and thus, for example, in production settings, for example, motion blur in motion scenes. Allows the use of more accurate techniques, such as techniques for computing. It is clear from the above description that the described boundary volume hierarchy has significant advantages when compared to 3D trees and other techniques, especially in hardware implementations and for large scenes. In the amortization analysis, the described boundary volume hierarchy is at least twice as good as the current 3D tree. In addition, the memory footprint can be predetermined and is linear with the number of objects.

II. Sorting Strategy for Rays Simulated Using Markov Chains and Quasi-Monte Carlo Methods Aspects of the invention described below use pixels in an image that use ray-tracing methods to calculate ray surface intersections. A method, technique, and method for efficient photorealistic simulation of light transport, suitable for use in computer graphics and other applications where a computer or other digital processor is used to calculate a pixel value of Provide a system.

  In particular, the present invention provides a method for making ray surface intersection testing more efficient by applying a sorting strategy that allows for fast processing of objects in a scene illuminated by a packet or unit of light such as photons. Provides improvements over conventional ray-tracing methods, including: optical packet or unit trajectories are simulated using Markov chains, which are simulated using appropriate quasi-Monte Carlo (QMC) methods Is done. The sorting strategy includes proximity sorting, and the present invention advantageously exploits a low-dimensional QMC method that uses such proximity sorting.

  The described techniques and systems can be used to solve integral equations in near real-time, a task that presents great difficulty with conventional systems.

  Previously, the use of quasi-Monte Carlo methods to simulate Markov chains, such as photon trajectories, provided only a slight improvement over random sampling. However, the QMC method according to the present invention can take advantage of the intrinsic low-dimensional structure of the solution operator, resulting in a much more efficient technique.

The description of this aspect of the invention is organized as follows.
1. 1. Markov chain and quasi-Monte Carlo method Simultaneous simulation of Markov chains 2.1 Analysis of deterministic techniques in one dimension 2.1.1 Simplification of techniques 2.1.2 When and why it works 2.2 Simplification in s dimensions Technique 2.3 Randomization 2.4 Sequence Randomization Quasi-Monte Carlo 3. Application to light transport simulation 3.1 Fredholm integration and Volterra integration Sorting Techniques and Systems 4.1 State Norm 4.2 Spatial Hierarchy 4.3 Bucket Sorting and Space Filling Curves 5. Result 6. Conclusion 7. General technique

1. Markov Chain and Quasi-Monte Carlo Method Mathematically, a Markov chain is a discrete-time probability sequence of system states. A Markov chain is a “memoryless” process in which, given the knowledge of the current state, the previous state is irrelevant in predicting the probability of the subsequent state. Markov chains can be used to model and analyze many different phenomena. For example, English was modeled using a Markov chain based on a calculation of the relative frequency of one word following another. Using these transition probabilities, English sentences can be generated. Many other evolutionary processes can also be described using Markov chains, including, for example, queuing, Brownian motion, particle transport analysis, and light transport simulations that are particularly relevant to the present description.

  More specifically, Markov chains can be used to simulate the trajectory of photons that illuminate objects in the scene of the simulated image. Light transport simulation is characterized by complex regions and discontinuities and reduction to higher dimensional integrals. In addition, light transport simulation is characterized by a low dimensional structure by repeated application of the same operator.

  The Markov chain itself can also be simulated using the quasi-Monte Carlo (QMC) method. The general idea behind the QMC method is that the known decision to provide a very uniform coverage of the selected interval instead of trying to adjust the number of randomly generated in some way. It starts with a theoretical line. In randomized QMC, the number has a random appearance, but the sequence is randomized in such a way that the desired property of covering the selected interval in a regular manner is retained as a group.

  There are many approaches that can be used to simulate Markov chains using QMC. One approach is to apply a high-dimensional, low-discrepancy point set. However, this approach offers significant advantages only for the middle dimension. Another approach is padding of a low-dimensional randomized point set. This second approach is easier to implement and achieves almost good results. The third technique combines state sorting and padding. This third approach results in an increase in performance where QMC capabilities are visible and can be used in conjunction with both deterministic and randomization techniques. Implementation of this technique is relatively easy.

  The characteristics of the process can be estimated by averaging the contributions of multiple Markov chains. Instead of simulating each trajectory individually, it was found that convergence can be significantly improved by simultaneously simulating Markov chains using correlation samples and sorting a group of states after each transition step. It was.

  Given the order on the state space, these approaches can be improved by sorting a group of Markov chains. In accordance with aspects of the present invention described below, the use of deterministic techniques provides simplification in the described techniques, and beneficial indications as to when and why the described sorting results in increased performance. I will provide a. As described further below, some or all of the described sorting strategies can be advantageously used to increase the efficiency of light transport simulations.

2. Simultaneous simulation of Markov chains To increase efficiency, Markov chains can be simulated simultaneously. The idea of improving the simultaneous simulation of Markov chains using the quasi-Monte Carlo method by an intermediate sorting step was first introduced by Lecot in a series of papers dealing with the Boltzmann equation and later dealing with the thermal equation. This idea was then used and refined by Morokoff and Caflisch to solve the heat equation on the grid, and recently L 'to incorporate a randomized version of the technique and partitioning for rare event simulation. Extended by Ecuyer, Tuffin, Demers et al. Independent research has been conducted in the field of computer graphics in a manner related to the above approach.

  In the following description, techniques for implementing improvements over the above scheme are described. The background of these improved techniques is described in H.W. Edited by Niederreiter, Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing 2002, pp. 329-44, Springer Verlag (2004), described in Lecot et al., “Quasi-Monte Carlo Methods for Estimating Transient Measurements of Discrete Time Markov Chains, partly by the literature,” It is incorporated herein by reference as if set forth in its entirety.

  The following description provides a conceptual framework for when and why the described techniques will be superior to conventional approaches. An advantageous implementation of the described technique is also described.

及び遷移行列

を有する離散状態空間Ｅ上で、マルコフ連鎖を同時にシミュレートする。基数ｂの（ｔ，２）−列

を使用して、Ｎ＝ｂ^ｍ個の連鎖が、並列にシミュレートされ、ここで、Ｘ_ｎ，ｌは、時間ステップｎにおける連鎖ｌの状態であり、

は、ゼロを含む自然数の集合である。さらに、当該技法は、ｍ＞ｔについて、Ｎ＝ｂ^ｍ個の後続要素ｘ_i，１が、基数ｂの（０，ｍ，１）−ネットを形成することを要求する。Ｌｅｃｏｔに示されるように、以下の基準が満たされる場合、当該技法は収束する。

2.1 Analysis of deterministic techniques in one dimension Deterministic techniques for one-dimensional problems will now be explained. The technique presented in Lecot uses an initial distribution

And transition matrix

A Markov chain is simultaneously simulated on a discrete state space E having (T, 2) -sequence of radix b

N = b ^m chains are simulated in parallel, where X _{n, l} is the state of chain l at time step n,

Is a set of natural numbers including zero. Further, the technique requires that N = b ^m subsequent elements x _{i, 1} form a (0, m, 1) -net of radix b for m> t. As shown in Lecot, the technique converges if the following criteria are met:

を初期化する。
ループ
何らかの順序を使用して、状態ベクトルをソートする。
ｎ：＝ｎ＋１
０≦ｌ＜Ｎについて、

を使用して、

を計算することによって、連鎖を継続させる。 X _{n, l} is the state of chain 1 at time step n. Therefore, the simulation can be performed as follows.
n: = 0
For 0 ≦ l <N, use x _{l, 2}

Is initialized.
Loop Sort the state vector using some order.
n: = n + 1
For 0 ≦ l <N,

using,

The chain is continued by computing

For N = ^{b m} pieces of chain, for m> t, ^{b m-number} of subsequent elements _{x i, 1} is, (0, m, 1) - so as to form a net, radix b (t, 2 )-Column x ⁱ is selected.

について、

となる。 It can be seen that the Sobol ′ (0,2) -column (φ ₂ (l), φ _S (l)) of the radix b = 2 satisfies the following requirements.
For m = 3> t = 0, N = 2 ³ = 8,

about,

It becomes.

  All substitutions are identical.

2.1.1 Simplification of technique It can be seen that stable sorting does not change the order of the states, assuming a uniform distribution. In fact, since the convergence condition remains unchanged, all substitutions can be omitted. Thus, the technique described above can be simplified.

について検討する。元のＳｏｂｏｌ’列は、Ｉ．Ｓｏｂｏｌ’、「Ｏｎ ｔｈｅ Ｄｉｓｔｒｉｂｕｔｉｏｎ ｏｆ Ｐｏｉｎｔｓ ｉｎ ａ Ｃｕｂｅ ａｎｄ ｔｈｅ Ａｐｐｒｏｘｉｍａｔｅ Ｅｖａｌｕａｔｉｏｎ ｏｆ Ｉｎｔｅｇｒａｌｓ」、Ｚｈ．Ｖｙｃｈｉｓｌ．Ｍａｔ．ｉＭａｔ．Ｆｉｚ．、７（４）：７８４−８０２（１９６７）において定義されており、同文献は、その全体があたかも説明されているかのように、参照により本明細書に組み込まれる。

についての簡単なコード例は、Ｔ．Ｋｏｌｌｉｇ及びＡ．Ｋｅｌｌｅｒ、「Ｅｆｆｉｃｉｅｎｔ Ｍｕｌｔｉｄｉｍｅｎｓｉｏｎａｌ Ｓａｍｐｌｉｎｇ」、Ｃｏｍｐｕｔｅｒ Ｇｒａｐｈｉｃｓ Ｆｏｒｕｍ、２１（３）：５５７〜５６３（２００２年９月）において提供されており、同文献は、その全体があたかも説明されているかのように、参照により本明細書に組み込まれる。 A Sobol 'sequence that is a (0,2) -sequence of radix b = 2 and satisfies the assumptions needed to maintain convergence conditions

To consider. The original Sobol 'column contains I.D. Sobol ', "On the Distribution of Points in a Cube and the Apploximation Evaluation of Integrals", Zh. Vychisl. Mat. iMat. Fiz. 7 (4): 784-802 (1967), which is incorporated herein by reference as if set forth in its entirety.

A simple code example for T. Kollig and A.K. Keller, “Efficient Multidimensional Sampling”, Computer Graphics Forum, 21 (3): 557-563 (September 2002), which is incorporated by reference as if set forth in its entirety. Incorporated herein.

を初期化する。その後、当該技法は、状態をソートすること（以下のセクション３において詳細に説明される）によって続行され、Ｐに従った遷移を実現するため、０≦ｌ≦Ｎについて、

を使用して、

を計算することによって、連鎖を継続させる。遷移の次の状態を選択するためのインデックス

は、実際には、基数２のｖａｎ ｄｅｒ Ｃｏｒｐｕｔ列Φ_２を使用し、これは、（０，１）−列であり、したがって、（０，ｍ，１）−ネットの系列である。例えば、ｍ＝３＞ｔ＝０を選択すると、Ｎ＝２^３＝８となり、

について、

となる。したがって、異なる時間ステップｎの間に使用されるすべてのインデックスは、実際には、すべてのｍについて同一である。 The simulation itself starts at time step n = 0 and uses x _{l, 2} for the realization of μ and states for 0 ≦ l ≦ N

Is initialized. The technique then continues by sorting the states (described in detail in Section 3 below), and for achieving a transition according to P, for 0 ≦ l ≦ N,

using,

The chain is continued by computing An index to select the next state of the transition

Actually uses a radix-2 van der Corp sequence [Phi] ₂ , which is a (0,1) -sequence and therefore a (0, m, 1) -net sequence. For example, if m = 3> t = 0, N = 2 ³ = 8,

about,

It becomes. Thus, all indexes used during different time steps n are actually the same for all m.

を仮定すると、収束定理（ｃｏｎｖｅｒｇｅｎｃｅ ｔｈｅｏｒｅｍ）が依然として適用されるが、より重要な安定なソーティングは、状態の順序を変化させない。したがって、実際には、インデックス置換は、収束条件に触れることなく、恒等置換（ｉｄｅｎｔｉｔｙ）として選択できることになる。同じことが、初期状態Ｘ_０，ｌを選択するために適用され、それが、（０，２）−列 ｘ_ｌの第２の要素だけを使用する、以下の簡略化された、しかし等価の技法をもたらす。
・ ｎ：＝０
・ ０≦ｌ＜２^ｍについて、ｘ_ｌ，２を使用して、Ｘ_０，ｌを初期化する。
・ ループ
− 適切な順序を使用して、状態ベクトルをソートする。
− ｎ：＝ｎ＋１
− ０≦ｌ＜２^ｍについて、Φ_Ｓ（ｌ＋ｎ・２^ｍ）を使用して、Ｘ_ｎ，ｌを計算することによって、連鎖を継続させる。 Uniform probability

, The convergence theorem still applies, but the more important stable sorting does not change the order of the states. Therefore, in practice, index replacement can be selected as identity replacement without touching the convergence condition. The same applies to select the initial state X _{0, l} , which uses only the second element of (0,2) -column x _l , the following simplified but equivalent Bring technique.
・ N: = 0
Initialize X _0,1 using x _{l, 2} for 0 ≦ l <2 ^m .
• Loop-Sort the state vector using the appropriate order.
-N: = n + 1
-For 0 ≦ l <2 ^m , continue the chain by calculating X _{n, l} using Φ _S (l + n · 2 ^m ).

によって理解することができ、Φ_Ｓ（ｌ）は、状態数ｌに依存する間隔１／（２^ｍ）を有するオフセットであり、一方、シフトΦ_Ｓ（ｎ）は、時間ステップｎにおけるすべての区間について同一である。図１１は、シフトΦ_Ｓ（ｎ）を図説する２次元配列７００を示している。 2.1.2 When and why it works The improved convergence of the scheme, observed in many applications, is used to achieve a transition of X _{n, l} according to P Sample Φ _S (l + n · 2 ^m )
It can be seen that it is caused by the structure of This is a decomposition of the radical inverse that exposes an implicit hierarchical structure.

Φ _S (l) is an offset with an interval 1 / (2 ^m ) depending on the number of states l, while the shift Φ _S (n) is the entire interval at time step n Are the same. FIG. 11 shows a two-dimensional array 700 illustrating the shift Φ _S (n).

Here, l = 0,. . . , 2 ^m −1, the whole Φ _S (l) must actually be a (0, m, 1) -net, and therefore must be an equidistant set of samples, so the low dimension The setting may lead to a misleading interpretation that the samples are shifted grids or layered samples. However, good performance is that of (t, m ′, s) -net for any m ′ where Φ _S (l) is in the (t, s) -sequence and therefore t ≦ m ′ ≦ m. Arises from the property of being a series. This is the same in the state space, and therefore means that b ^{m ′} states that are linked by post-sorting order sample transitions by the (t, m ′, s) -net, which is good Guarantees a discrete density approximation. Maximum improvement is obtained when all 2 ^m chains are in the same state. The better the chain states are in the state space, the smaller the performance improvement.

2.2 Simplification technique in s dimension Using the (t, s) -sequence of base b, which is a sequence of (t, m, s) -net, the scheme also works in s dimension. A Markov chain whose states are similar after sorting is guaranteed to sample transition probabilities with low discrepancy samples. The simplification technique in the s dimension now has the following aspect.
・ N: = 0
Initialize X _{0, l} using the quasi-Monte Carlo point _xl .
• Loop-Sort the state vector using the appropriate order.
-N: = n + 1
- (t, s) - using subsequent samples _{x l} from _{column, X n,} by calculating the _l, to continue the chain.

  Some simulations require trajectory splitting to obtain some local subtle effects. This has already been addressed using a more complex approach, but in practice it is simpler for the state to be split by simply taking more samples from the (t, s) -column. Can be achieved in a way.

This is a result of the fact that it is not necessary to simulate exactly b ^m chains simultaneously. To enable maximum performance gain, take subsequent samples from the (t, s) ^-sequence and minimize the number of points b ^m in the subsequent (t, m, s) -net Only is important. However, the selection of the (t, s) -column is limited by the condition that the (0, s) -column exists only if b ≧ s and m> t.

Note that other radical inversion-based point sets, such as the Halton sequence or its scrambled variants, satisfy similar characteristics as the (t, s) -sequence and yield similar performance gains. . FIG. 12 shows a graph 720 illustrating the Halton row, the point x _l = (φ ₂ (l), φ ₃ (l)). In the Halton column shown, the root inverse transform implies a hierarchical structure with subsequent points falling into the largest “hole”. This property is suitable for any starting point, and the (t, s) -column has the same property.

2.3 Randomization Although there is currently no general proof about the convergence of deterministic techniques in higher dimensions, the technique is unbiased by newly randomizing quasi-Monte Carlo points at each time step n. Become. In practice, this is the case for padded replica sampling, so the discussion of unbiasedness is simpler. However, randomization deserves special attention.

  The most efficient implementation consists of selecting the (t, s) -column with basis b = 2, then taking the subsequent sample, and subjecting that sample to an XOR operation with an s-dimensional random vector. This random vector is newly extracted after each transition step. However, random scrambling changes the order in which the points are listed, so the local characteristics of the (t, m, s) -net sequence are also changed.

  This result can be used as an explanation for some of the effects seen using previous approaches. The Sobol and Korovov points used in array (R) QMC simulations are worse than their transformed equivalents (Gray code for Sobol, Baker transform for Korovov) up to 10 times the variance reduction. An explanation for this is found in the dot structure. The (t, m, s) -net sequence extracted from the Sobol sequence is locally worse than its Gray code deformation. The same is true for the Korovov lattice and its transformed version.

2.4 Array Randomization Quasi-Monte Carlo The following is a randomization technique for the s-dimensional problem.
n: = 0, point sequence _xl is selected.
Instantiate randomization.
Initialize X _{0, l} using rand (x _l ) for 0 ≦ l <N.
Loop • Sort the state vector using the appropriate order.
・ N: = n + 1
• Instantiate randomization.
Continue the chain by computing X _{n, l} using rand (x _l ) for 0 ≦ l <N.

The technique has many interesting properties. First, the technique, one set was summarized in the table _{x l} for (one tabulated set), and a unbiased for continuous stream of _{x l.} Second, splitting and absorption are essential when using streams. This is true in Russian roulette system, it is also true in the case of dividing by taking more consecutive samples from x _l.

3. Light Transport Simulation For numerical evidence, the technique described above is applied to light transport simulation to synthesize realistic images. The underlying integral equation can be reformulated as a path integral.

  FIG. 13 is a diagram showing a sampling transport route space 740 by bidirectional route tracing. The trajectory from eye 742 to scene object 746 and the trajectory from light source 744 to scene object 746 are generated by a Markov chain and connected to determine the amount of light transported. The sampling path space shown in FIG. 13 corresponds to simulating a Markov chain, where the path is established by ray tracing and scattering events. The initial distribution is determined by the emission characteristics of the light source 744, and the transition probability is given by the bidirectional reflection distribution function on the impact surface.

  To solve the path integral, there are at least two basic strategies that use high-dimensional low discrepancy points or pad low-dimensional low discrepancy points. The latter approach is compatible with the above discovery, where high-dimensional events are configured as subsequent transitions of Markov chains. When measured, the use of high-dimensional sequences for Markov chain simulations is small enough that there is no difference, but the use of padded sequences is computationally faster and requires less implementation. Only effort is required. It is even easier for practitioners in the rendering industry.

  In addition, since the hierarchical properties of the (t, s) -sequence or Halton sequence and its scrambled variants are much better when the dimension is small, as explained above, the low-dimensional approach is Allows much better results.

  14A-B and FIGS. 15A-B show rendered images 760, 770, 780 showing the difference between the simulation using the Monte Carlo method and the simulation using the array randomized quasi-Monte Carlo (A-RQMC) method. 790. Simulations using the Monte Carlo method are shown in images 760 and 780 in FIGS. 14A and 15A. Simulations using the array randomized quasi-Monte Carlo method are shown in images 770 and 790 in FIGS. 14B and 15B. Arrows 762, 772, 782, 792 at the top of each rendered image indicate the emission position and direction. 14A-B and 15A-B, in these examples, the optical path simulated using the array randomized QMC method is much more uniform than the optical path simulated using the Monte Carlo method. It can be seen that they are distributed.

3.1 Fredholm integral and Volterra integral The integral equation underlying the light transport can be regarded as a Fredholm integral equation or a Volterra integral equation, which is important for implementation.

は、点ｘにおいて方向ωの方に表面から反射された放射であり、積分演算子の定義域

は、ｘにおける法線ｎ（ｘ）と軸平行な半球である。

Is the radiation reflected from the surface in the direction ω at the point x, and the domain of the integral operator

Is a hemisphere that is axially parallel to the normal n (x) at x.

  16A-D are a series of diagrams 800, 810, 820, 830 showing photon trajectories 802, 812, 822, 832 starting from a plurality of starting points 804, 814, 824, 834 on the light source. After tracing the ray, when it strikes the surface 806, 816, 826, 836, the bi-directional reflection distribution function is sampled to determine the direction of scattering and continue the path.

をもたらし、他方で、ｘとは独立に単位球Ｓ^２の全方向にわたって積分が実行されるような、第２種のＦｒｅｄｈｏｌｍ積分方程式をもたらす。 _fr is a bi-directional reflection distribution function describing the optical surface property and L is the incident radiation. The use of this integral operator is the second kind of Volterra integral equation where the integral domain depends on x.

On the other hand, a Fredholm integral equation of the second kind is obtained such that the integration is performed in all directions of the unit sphere S ² independently of x.

への写像は、当技術分野において知られている。 Although the use of the latter approach of generating all directions and eliminating directions with negative scalar products relative to the surface normal n (x) is computationally attractive, the first approach is locally Requires the generation of directions in the hemisphere, which must be converted to a smooth surface frame, which is more costly. ^{S 2} or from the unit sphere

The mapping to is known in the art.

A much more important argument for generating omnidirectional is related to the approach to be described. By sorting, it can happen that two nearby surface points with different surface normals, for example at the corners, use subsequent samples in the (t, s) -sequence to generate the scattering direction. Although omnidirectional generation now works well, directional generation in two different local frames using subsequent samples breaks down the low discrepancy. These discontinuities become clearly visible. However, the use of all directions, the sampling of importance by _{f r} (importance), or does not allow the cosine term, it is often inconvenient, deserves further study.

4). Sorting Techniques and Systems Sorting can be used to increase the efficiency of the techniques described above. As described herein, a (t, s) -sequence is a sequence of (t, m, s) -nets of base b and is similar to the Halton sequence and its variants. A similar state X _{n, l} uses successive points. Thus, the states should be ordered by proximity in the state space. The closer the state is in the state space, the more uniformly the state space is searched.

  To make the states as close as possible to each other, the states must be listed in an order that minimizes the sum of the distances of neighboring states. This is actually related to the traveling salesman problem, where each city is exactly 1 for a given set of cities and the cost of traveling from one city to another. After the visit, the cheapest patrol trip is required to return to the departure city. This problem has already been studied by Euler in the early 18th century and unfortunately can be shown to be NP-hard.

Our problem is very similar except that it is not necessary to return from the last state of the route to the first state. Several techniques are already available to efficiently compute a solution to the traveling salesman problem that is an approximation but close to optimal. However, since the execution time of these techniques shows the amount of calculation of O (N ² ) only by calculating the distance matrix, it is not acceptable in the simulation, but the amount of calculation of O (N) in the classical Monte Carlo method It is desirable to maintain the technique as close as possible.

  In the following, a number of sequences according to aspects of the present invention to achieve fast sorting for a high dimensional state space will be described.

が、状態空間上で順序を定義するために使用される。しかし、類似のノルムが、必ずしも状態空間における近接性を示すわけではない。これについての簡単な例は、輸送シミュレーションにおける、空間内で遠く離れて配置された粒子の類似のエネルギーである。 4.1 Norm of states The average complexity of quicksort is O (NlogN), but in some scenarios there are even O (N) techniques, such as radix sort, but it adds Requires temporary memory. In order to use these one-dimensional sorting techniques, the multidimensional state must be reduced to one dimension. Of many choices, often some norm

Are used to define the order on the state space. However, similar norms do not necessarily indicate proximity in state space. A simple example of this is the similar energy of particles located far away in space in a transport simulation.

4.2 Spatial Hierarchy A second possibility to enable multidimensional sorting is the use of a spatial hierarchy to define order on the state. An efficient data structure for this purpose is a BSP tree, a special axis-parallel subset thereof, a kD tree, or a boundary volume hierarchy. Building such a binary hierarchy is simple. The space is recursively subdivided using a plane selected by some heuristic method. The construction is performed with O (NlogN) on average. Hierarchical traversal enumerates leaves in order of proximity. This traversal is trivial if the tree is left-balanced so that it can be stored in an array.

  If the spatial hierarchy has to be used anyway, for example to accelerate ray tracing, there is no additional construction time for the hierarchy. In that case, the particles are stored as a linked list associated with the leaves of the hierarchy.

  FIGS. 17A-G, 18A-B, and 19A-C show that the classification of states into leaves of a spatial hierarchy defines the order of proximity by traversing the hierarchy in order. FIG.

  FIG. 17A shows an exemplary scene 840 that includes a light source 850 and an object 860 to be rendered. The object 860 has a star shape with five vertices. In FIG. 17B, the object 860 is bounded by an axis parallel bounding box 870. 17C-G, the bounding box 870 is recursively subdivided into branches. Each branch ends with a leaf.

  In FIG. 17C, the bounding box 870 is divided into a first generation left branch L and a first generation right branch R by a dividing plane 871 passing through the first vertex 861 of the object 860.

  In FIG. 17D, the left branch L is split by a second split plane 872 that passes through the second vertex 862 of the object 860, resulting in a second generation left branch LL and a second generation right branch LR.

  In FIG. 17E, the third dividing plane 873 passes through the third vertex 863 of the object 860 and divides the branch LR into the third generation left branch LRL and right branch LRR.

  In FIG. 17F, the fourth and fifth split planes 874 and 875 pass through the fourth and fifth vertices 864 and 865 of the object 860. As shown in FIG. 17F, the result is five leaves LRL, LRR, RR, LLL, LLR, and RL, one for each of the five points of star 860.

  FIG. 17G shows the illumination of individual points within each leaf. In FIG. 18A, each leaf is shaded to show its brightness. In FIG. 18B, the leaves are hierarchically ordered according to their brightness. From “darkest” to “lightest”, the leaves are arranged in the following order: LRL, LRR, LLR, RR, RL.

  Unfortunately, the quality of this order is strongly determined by the quality of the spatial hierarchy used for the simulation, which means that if the number of leaves in the hierarchy is much less than the number of chains N, this This is especially problematic because it causes the state to map to the same leaf.

4.3 Bucket sorting and space filling curves Bucket sorting can be used to guarantee linear time complexity. In the extension of the simple technique outlined in section 2.1 to the s dimension, multi-dimensional states are classified into buckets by a first dimension of states, after which the state of each bucket is converted to a second dimension. According to the above, it was classified into buckets and so on. This procedure works well, but has the problem that states that are close in the state space can be separated by the sorting procedure. In addition, the hierarchical structure of each dimension must be used, which causes a curse of dimension in the number of Markov chains that are simulated simultaneously.

  Thus, the state space becomes an equal volume cell or “voxel” that serves as a bucket. Each state bucket is found by truncating the state coordinates according to the resolution of the voxel grid. Note that this applies to both continuous and discrete state spaces. The enumeration of voxels by proximity brings the desired order to the state and can be performed in linear time for the number of voxels.

  The order in which voxels are listed by proximity is given by a space filling curve such as a Peano curve, a Hilbert curve, or H-indexing. These curves ensure that all voxels are visited exactly once and the overall path length is relatively short. For problems related to large shapes, which is the case in our own simulations, this can even be one of the few possibilities for generating fast, memory-efficient approximate solutions to the traveling salesman problem. However, these curves are rather expensive to evaluate and need to be tabulated to be efficient, or are not available for higher dimensions.

  Fortunately, the Z curve, also known as Lebesgue curve or Morton order, avoids these problems. Given the integer coordinates of a bucket in a multidimensional space, its one-dimensional Z-curve index is simply calculated by interleaving the coordinate values bit by bit, as shown in FIGS.

19A-C show a series of two-dimensional Z-curves 891, 893, 895 for 2 × 2 (FIG. 19A), 4 × 4 (FIG. 19B), and 16 × 16 (FIG. 19C) buckets. 890, 892, and 894 are shown. If the upper left is the origin (0, 0), the point marked by x has integer coordinates (3, 4), which corresponds to binary (011, 100) ₂ . Its binary Z curve index 100101 ₂ is calculated by interleaving the binary coordinates in bit units.

  This is easy and fast to compute for any dimension and problem size and does not require additional memory. The results were found not to be as good in the overall context as for example the Hilbert curve. However, the average split quality and the average / worst case log index range are equally good. The problem arises in highly symmetric scenes, such as Shirley's “Scene 6” shown in FIG. 20A, used for numerical experiments. The states on the wall parallel to the (x, y) plane are sorted very well, but the states placed on the other two walls parallel to the (y, z) plane are alternated by the curves. Visited, which may lead to correlation artifacts in some scenarios.

5. Numerical Evidence Following the discussion above, a Halton sequence, randomized by Cranley-Patterson rotation, replaced by Faure, was used to obtain an unbiased error estimate. For sorting, the Z curve order gave the best results for a particular setting and was used for the following experiments. It is further noted that the omission of randomization has been numerically verified that it has no significant effect on the accuracy of the results. In numerical experiments, four approaches for simulating Markov chains were compared.

  MC: Uniform random numbers generated by Mersenne Twister were used for classical Monte Carlo sampling.

  A high-dimensional Halton sequence, randomized by RQMC: Cranley-Patterson rotation, with substitution by Faure, was used, and pairs of elements were used to sample two-dimensional emission and scattering events.

  lo-dim RQMCS: A two-dimensional Halton sequence randomized by Cranley-Patterson rotation was used. The Z curve was used to enumerate the bucket sorted state.

  hi-dim RQMCS: A high-dimensional Halton sequence randomized by Cranley-Patterson rotation, with substitution by Faure, was used. The Z curve was used to enumerate the bucket sorted state.

  In the first experiment, a strong global illumination technique was used to calculate the path integral. FIG. 20A shows a test scene “Shirley 6”, which is a relatively simple light transport setting, and FIG. 21A shows a test scene “Invisible Date”, which is a more complicated light transport setting. 20B and 21B are graphs 901 and 911 showing the average RMS error for the master solution, and FIGS. 20C and 21C are graphs 902 and 912 showing the global variance (averaged over the entire image). 20D and FIG. 21D are graphs 903 and 913 showing pixel-based variance. The figure was obtained by averaging 10 independent runs with varying number of Markov chains. The measured numbers are convincing only for simple test scenes. Even in complex cases, the performance gain over Monte Carlo sampling cannot be measured because the number of independent runs is too small. More experiments were not possible due to excessive execution time.

  However, visual error tells a dramatically different story, as seen in FIGS. 22A-C, revealing the unmistakable superiority of the new technique, even in very difficult settings. 22A-C show a visual comparison for the test scene “Invisible Date” using 300 chains for simulation. 22A shows a test scene rendered using Monte Carlo, FIG. 22B shows a test scene rendered using randomized quasi-Monte Carlo, and FIG. 22C uses randomized quasi-Monte Carlo with sorting. Shows the rendered test scene.

22A-C, the only light source is not visible because it is located on the ceiling of the adjacent room. Because of the better distribution of photons, randomized quasi-Monte Carlo (RQMC) (FIG. 22B) performs better visually than Monte Carlo (MC) (FIG. 22A), as can be seen by the reduced shadow artifacts. RQMC (RQMCS) with sorting (using 256 ³ voxels for bucket sort) (FIG. 22C), more photons enter the second room and even the back of the door is very well illuminated So it is even better.

  This case is no exception and can be observed for many test cases. Standard error measurement is not an appropriate error measurement for visual quality, and this only highlights a known but unsolved problem in computer graphics.

  Figures 23A-D show measurements for a very difficult light transport problem, with photons traced directly from the light source and the final path segment connected to the camera (one technique of bidirectional path tracing). FIG. 23A shows a schematic diagram of a labyrinth test scene 940 with the floor and ceiling removed for illustration. The camera is placed in the lower room and the light source is placed at the other end of the connecting tunnel. Graph 950 shown in FIG. 23B is a “box-and-whisker” diagram, and is the average of the total radiation received by the camera of 1048576 simulated light paths for 50 independent implementations using each technique. The quantity (marked by x) is shown. The graph 960 shown in FIG. 23C is an enlarged view of the portion of interest in the graph shown in FIG. 23B. Table 970 shown in FIG. 23D used the L1 and L2 norms (dispersion) of the total radiation received by the camera and each pixel (256 × 256 pixel resolution), averaged over 50 independent implementations. The deviation from the master solution is finally displayed. In this difficult setting, the new method (RQMCS) is undoubtedly superior.

  Contrary to the above measurement, only one simultaneously simulated Markov chain is considered. Now, a sufficient number of experiments are computationally feasible, and the superiority of the new technique is clearly visible.

6). CONCLUSION It can be seen that the systems and techniques described above have many advantages over the prior art. These include the following: simplification of techniques that are intuitively focused on sequence characteristics, improved convergence, linear complexity, and application to light transport simulations. The systems and techniques described above are one step closer to formulating the technique as an integral, and in that respect it is relatively easy to show consistency.

The inventors have successfully simplified techniques for simultaneously simulating Markov chains and provided an intuition as to when and why state sorting can improve convergence. In addition, the technique is no longer limited by the curse of dimension, and the simulation can just use the transition probabilities P≡P _n that can change over time, so there is no limit to the homogenous Markov chain. Thus, the techniques described herein can be used to simulate an inhomogeneous Markov chain, i.e., where P≡P _n depends on time step n. .

  Experiments have shown that not all (t, s) -sequences or root inverse-based point sequences are equally good. This deserves further characterization.

  The techniques described herein are even simpler when rank 1 lattice sequences are applicable. However, previous constructions lack the (t, s) -sequence characteristics required for improved performance. According to a further aspect of the invention, it may be possible to construct a suitable rank-1 lattice sequence to achieve the desired result.

  It is also interesting to establish relevance to recent work in that direction, since the inventors use all directions, i.e., integrate over the product of spheres.

7). General Techniques FIGS. 24-30 are flowcharts of a number of techniques and sub-techniques according to the above aspects of the invention. It should be noted that some or all of the following techniques and sub-techniques, including individual boxes and groups of boxes, may be combined as desired to implement the invention described herein. .

FIG. 24 is a flowchart of a first general technique 1000 according to the present invention. Technique 1000 includes the following.
Box 1001: Simulate a Markov chain using a quasi-Monte Carlo method.
Box 1002: Sort status.
(Can include sorting by state, including proximity sorting)
Box 1003: Simulate multiple Markov chains simultaneously.
(A quasi-Monte Carlo point can be used. A quasi-Monte Carlo point can be selected using either a Halton sequence, a variant of the Halton sequence, a grid sequence, or a (t, s) -sequence. Any number of chains can be used. (This can include adding trajectories on the fly, such as by trajectory splitting. Russian roulette methods or other techniques to simulate particle absorption can be used)

FIG. 25 is a flowchart of a further technique 1100 according to the present invention for simultaneously simulating multiple Markov chains for an s-dimensional problem. Technique 1100 includes the following:
Box 1101: n: = 0
(Set n equal to 0)
Box 1102: X _{0, l} is initialized using the quasi-Monte Carlo point _xl .
Box 1103: Sort state vector using appropriate order.
Box 1104: n: = n + 1
(Increase n by 1)
Box 1105: Continue the chain by calculating X _{n, l} using the next quasi-Monte Carlo point _xl .
(Any selected point set or sequence with any chosen randomization and state ordering by proximity in state space can also be used)
Boxes 1103 to 1105 are executed in a loop.

FIG. 26 is a flowchart of an additional technique 1200, one or more of which may be performed in conjunction with technique 1100 shown in FIG. These techniques 1200 include:
Box 1201: Sort by brightness or state norm.
Box 1202: Sort using spatial hierarchy.
(May include using a binary hierarchy, or any other axis-parallel subset of a BSP tree, kD tree, or BSP tree structure, or a boundary volume hierarchy, or regular voxels. The binary hierarchy is selected. Can be constructed by recursively subdividing using the plane selected by the heuristic and then traversing the hierarchy in order to enumerate the leaves of the hierarchy in order of proximity. (It can also include left balancing the tree and saving the tree in an array data structure to allow for efficient traversal of the hierarchy)
Box 1203: Sort using buckets enumerated by at least one space filling curve.

FIG. 27 is a flowchart of a further technique 1300 according to the present invention. Technique 1300 includes the following.
Box 1301: Select n: = 0, point sequence _xl .
Box 1302: Instantiate randomization.
Box 1303: Initialize X _{0, l} using rand (x _l ).
Box 1304: Sort the state vector using some order.
Box 1305: n: = n + 1
Box 1306: Instantiate randomization.
Box 1307: Continue chaining by calculating X _{n, l} using rand (x _l ).
Boxes 1304-1307 are executed in a loop.

FIG. 28 is a flowchart of an additional technique 1400 that may be performed in conjunction with technique 1300 of FIG. The additional technique 1400 randomizes the quasi-Monte Carlo points at each time step n and includes:
Box 1401: Select a (t, s) -column of radix b = 2 from which subsequent samples are taken.
Box 1402: Sample is XORed with s-dimensional random vector, and a random vector is generated after each transition step.

FIG. 29 is a flowchart of a further technique 1500 according to the present invention. Technique 1500 includes the following.
Box 1501: An axis parallel bounding box delimits the object to be rendered.
Box 1502: recursively or iteratively splits the bounding box into successive left and right branches, each ending with a leaf, by applying a splitting plane that passes through the selected points of the object; To generate a plurality of leaves.
Box 1503: The leaves are ordered hierarchically according to their brightness.
Box 1504: Perform a bucket sort on a matrix of a selected size divided into a selected number of buckets.
Box 1505: Use the selected space filling curve to follow in the matrix, and within each bucket, a smaller space filling curve is used to follow within the individual cells in the matrix.

FIG. 30 is a flowchart of a further technique 1600 according to the present invention. The technique 1600 configures a simulation that enables a light transport simulation to synthesize a real-looking image and includes:
Box 1601: The integral equation underlying the light transport simulation is reformulated as a path integral so that sampling the path space corresponds to simulating a Markov chain, where the path is a ray tracing and Established by scattering events.
Box 1602: The initial distribution is determined by the emission characteristics of the light source or sensor, and the transition probability is determined by the bidirectional reflection distribution function and the bidirectional subsurface scattering distribution function for the surface in the image.
Box 1603: Process transmission effects by utilizing path integrals, initial distribution, and transition probabilities.
Box 1604: The path integral is solved by using high-dimensional quasi-Monte Carlo points, with or without any kind of randomization, or by padding low-dimensional quasi-Monte-Carlo points.
(The integral equation underlying the light transport simulation is regarded as a Fredholm integral equation or a Volterra integral equation)

  While the above description includes details that enable one skilled in the art to practice the invention, the description is intended to be exemplary in nature and to those skilled in the art who would benefit from these teachings. It should be understood that many modifications and variations are apparent. Accordingly, it is intended that the invention herein be determined only by the claims appended hereto, and that the claims be interpreted as broadly as permitted by the prior art.

Computer program list

CODE LISTING 2.2.1

void Triangle :: Transform ()
{
Point * p = (Point *) this;
Vector n3d;
Vector n_abs = n3d = (p [1] -p [0]) | (p [2] -p [0]);
// search largest component forprojection (0 = x, 1 = y, 2 = z)
uintCast (n_abs.dx) & = 0x7FFFFFFF;
uintCast (n_abs.dy) & = 0x7FFFFFFF;
uintCast (n_abs.dz) & = 0x7FFFFFFF;
// Degenerated Triangles must be handled (set edge-signs)
if (! ((n_abs.dx + n_abs.dy + n_abs.dz)> DEGEN_TRI_EPSILON))
// (! (...)> EPS) to handle NaN's
{
d = p [0] .x;
p0.u = -p [0] .y;
p0.v = -p [0] .z;
nu = nv = 0.0f;
e [0] .u = e [1] .u = e [0] .v = e [1] .v = 1.0f;
return;
}
U32 axis = 2;
if (n_abs.dx> n_abs.dy)
{
if (n_abs.dx> n_abs.dz)
axis = 0;
}
else if (n_abs.dy> n_abs.dz)
axis = 1;
Point p03d = p [0];
Point p13d = p [1];
Point p23d = p [2];
float t_inv = 2.0f / n3d [axis];
e [0] .u = (p23d [PlusOneMod3 [axis]]-p03d [PlusOneMod3 [axis]]) * t_inv;
e [0] .v = (p23d [PlusOneMod3 [axis + 1]]-p03d [PlusOneMod3 [axis + 1]]) * t_inv;
e [1] .u = (p13d [PlusOneMod3 [axis]]-p03d [PlusOneMod3 [axis]]) * t_inv;
e [1] .v = (p13d [PlusOneMod3 [axis + 1]]-p03d [PlusOneMod3 [axis + 1]]) * t_inv;
t_inv * = 0.5f;
nu = n3d [PlusOneMod3 [axis]] * t_inv;
nv = n3d [PlusOneMod3 [axis + 1]] * t_inv;
p0.u = -p03d [PlusOneMod3 [axis]];
p0.v = -p03d [PlusOneMod3 [axis + 1]];
d = p03d [axis] + nu * p03d [PlusOneMod3 [axis]] + nv * p03d [PlusOneMod3 [axis + 1]];
}


CODE LISTING 2.2.2

U32 * idx = pointer_to_face_indices;
U32 ofs = projection_case;
for (U32 ii = num_triData; ii; ii-, idx ++)
{
float t = (triData [* idx] .d -ray.from [ofs]
-triData [* idx] .nu * ray.from [PlusOneMod3 [ofs]]
-triData [* idx] .nv * ray.from [PlusOneMod3 [ofs + 1]])
/ (ray.d [ofs] + triData [* idx] .nu * ray.d [PlusOneMod3 [ofs]]
+ triData [* idx] .nv * ray.d [PlusOneMod3 [ofs + 1]]);
if (uintCast (t) -1> uintCast (result.tfar)) //-1 for + 0.0f
continue;
float h1 = t * ray.d [PlusOneMod3 [ofs]] + ray.from [PlusOneMod3 [ofs]]
+ triData [* idx] .p0.u;
float h2 = t * ray.d [PlusOneMod3 [ofs + 1]] + ray.from [PlusOneMod3 [ofs + 1]]
+ triData [* idx] .p0.v;
float u = h1 * triData [* idx] .e [0] .v -h2 * triData [* idx] .e [0] .u;
float v = h2 * triData [* idx] .e [1] .u -h1 * triData [* idx] .e [1] .v;
float uv = u + v;
if ((uintCast (u) | uintCast (v) | uintCast (uv))> 0x40000000)
continue;
result.tfar = t;
result.tri_index = * idx;
}


CODE LISTING 2.3.1

Point * p = (Point *) & triData [tri_index];
int boxMinIdx, boxMaxIdx;
// boxMinIdx and boxMaxIdx index thesmallest and largest vertex of the triangle
// in the component dir [0] of the splitplane
if (p [0] [dir [0]] <p [1] [dir [0]])
{
if (p [2] [dir [0]] <p [0] [dir [0]])
{
boxMinIdx = 2;
boxMaxIdx = 1;
}
else
{
boxMinIdx = 0;
boxMaxIdx = p [2] [dir [0]] <p [1] [dir [0]]? 1: 2;
}
}
else
{
if (p [2] [dir [0]] <p [1] [dir [0]])
{
boxMinIdx = 2;
boxMaxIdx = 0;
}
else
{
boxMinIdx = 1;
boxMaxIdx = p [2] [dir [0]] <p [0] [dir [0]]? 0: 2;
}
}
/ * If the triangle is in the split planeor completely on one side of the split plane
is decided without any numerical errors, ie at the precision the triangle is
entered to the rendering system.Usingepsilons here is wrong and not necessary.
* /
if ((p [boxMinIdx] [dir [0]] == split) && (p [boxMaxIdx] [dir [0]] == split)) // in split plane?
{
on_splitItems ++;
if (split <middle_split) // put tosmaller volume
left_num_divItems ++;
else
{
unsorted_border--;
U32 t = itemsList [unsorted_border];
right_num_divItems--;
itemsList [right_num_divItems] = itemsList [left_num_divItems];
itemsList [left_num_divItems] = t;
}
}
else if (p [boxMaxIdx] [dir [0]] <= split) // triangle completely left?
left_num_divItems ++;
else if (p [boxMinIdx] [dir [0]]> = split) // triangle completely right?
{
unsorted_border--;
U32 t = itemsList [unsorted_border];
right_num_divItems--;
itemsList [right_num_divItems] = itemsList [left_num_divItems];
itemsList [left_num_divItems] = t;
}
else
// and now detailed decision, trianglemust intersect split plane ...
{
/ * In the sequel we determine whether atriangle should go left and / or right, where we already know that it mustintersect the split plane in a line segment.
All computations are ordered so that themore precise computations are done first.Scalar products and cross productsare evaluated last.
In some situations it may be necessaryto expand the bounding box by
an epsilon. This, however, will blow up the required memory by large amounts.
If such a situation is encountered, itmay be better to analyze it numerically in order not to use any epsilons ...
Arriving here we know thatp [boxMaxIdx] [dir [0]] <split <p [boxMaxIdx] [dir [0]]
and that p [boxMidIdx] [dir [0]] \ in [p [boxMaxIdx] [dir [0]], p [boxMaxIdx] [dir [0]]].
We also know, that the triangle has anon-empty intersection with the current
voxel.The triangle also cannot lie inthe split plane, and its vertices cannot
lie on one side only.
* /
int boxMidIdx = 3-boxMaxIdx -boxMinIdx; // missing index, found by 3 = 0 + 1 + 2
/ * We now determine the vertex that isalone on one side of the split plane.
Depending on whether the lonely vertexis on the left or right side,
we have to later swap the decision, whether the
triangle should be going to the left orright.
* /
int Alone = (split <p [boxMidIdx] [dir [0]])? boxMinIdx: boxMaxIdx;
int NotAlone = 3-Alone-boxMidIdx;
// == (split <p [boxMidIdx] [dir [0]])? boxMaxIdx: boxMinIdx;
// since sum of idx = 3 = 0 + 1 + 2
float dist = split-p [Alone] [dir [0]];
U32 swapLR = uintCast (dist) >>31; // == p [Alone] [dir [0]]>split;
/ * Now the line segments connecting thelonely vertex with the remaining two verteces
are intersected with the split plane.a1and a2 are the intersection points.
The special case "if (p [boxMidIdx] [dir [0]] == split)" [yields ax / x, which could
be optimized] does not help at all sinceit only can happen as often as the highest
valence of a vertex of the mesh is ...
* /
float at = dist / (p [boxMidIdx] [dir [0]]-p [Alone] [dir [0]]);
float at2 = dist / (p [NotAlone] [dir [0]]-p [Alone] [dir [0]]);
float a1x = (p [boxMidIdx] [dir [1]] -p [Alone] [dir [1]]) * at;
float a1y = (p [boxMidIdx] [dir [2]] -p [Alone] [dir [2]]) * at;
float a2x = (p [NotAlone] [dir [1]] -p [Alone] [dir [1]]) * at2;
float a2y = (p [NotAlone] [dir [2]]-p [Alone] [dir [2]]) * at2;
// n is a vector normal to the line ofintersection a1a2 of the triangle
// and the split plane
float nx = a2y-a1y;
float ny = a2x-a1x;
// The signs indicate the quadrant of the vector normal to the intersection line
U32 nxs = uintCast (nx) >>31; // == (nx <0.0f)
U32 nys = uintCast (ny) >>31; // == (ny <0.0f)
/ * Numerical precision: Due tocancellation, floats of approximately same exponent
should be subtracted first, beforeadding something of a different order of
magnitude.All brackets in the sequelare ESSENTIAL for numerical precision.
Change them and you will see more errors in the BSP ...
pMin is the lonely point in thecoordinate system with the origin at
bBox.bMinMax [0]
pMax is the lonely point in thecoordinate system with the origin at
bBox.bMinMax [1]
* /
float pMinx = p [Alone] [dir [1]] -bBox.bMinMax [0] [dir [1]];
float pMiny = p [Alone] [dir [2]] -bBox.bMinMax [0] [dir [2]];
float pMaxx = p [Alone] [dir [1]] -bBox.bMinMax [1] [dir [1]];
float pMaxy = p [Alone] [dir [2]] -bBox.bMinMax [1] [dir [2]];
// Determine coordinates of the boundingbox, however, with respect to p + a1 being the origin.
float boxx [2];
float boxy [2];
boxx [0] = (pMinx + a1x) * nx;
boxy [0] = (pMiny + a1y) * ny;
boxx [1] = (pMaxx + a1x) * nx;
boxy [1] = (pMaxy + a1y) * ny;
/ * Test, whether line of intersection of the triangle and the split plane passes by the
this is done by indexing the coordinates of the bounding box by the
quadrant of the vector normal to theline of intersection.In fact this is
the nifty implementation of the 3d testintroduced by in the book with Haines:
"Real-Time Rendering"
By the indexing the vertices areselected, which are farthest from the line.
Note that the triangle CANNOT completelypass the current voxel, since it must have
a nonempty intersection with it.
* /
U32 resultS;
if (pMinx + MAX (a1x, a2x) <0.0f) // line segment of intersection a1a2 left of box
resultS = uintCast (pMinx) >>31;
else if (pMiny + MAX (a1y, a2y) <0.0f) // line segment of intersection a1a2 below box
resultS = uintCast (pMiny) >>31;
else if (pMaxx + MIN (a1x, a2x)> 0.0f) // line segment of intersection a1a2 right of box
resultS = (pMaxx>0.0f);
else if (pMaxy + MIN (a1y, a2y)> 0.0f) // line segment of intersection a1a2 above box
resultS = (pMaxy>0.0f);
else if (boxx [1 ^ nxs]> boxy [nys])
// line passes beyond bbox? => triangle can only be on one side
resultS = (a1y * a2x> a1x * a2y);
// sign of cross product a1 x a2 ischecked to determine side
else if (boxx [nxs] <boxy [1 ^ nys])
resultS = (a1y * a2x <a1x * a2y);
else
// Ok, now the triangle must be bothleft and right
{
stackList [currStackItems ++] = itemsList [left_num_divItems];
unsorted_border--;
itemsList [left_num_divItems] = itemsList [unsorted_border];
continue;
}
if (swapLR! = / * ^ * / resultS)
{
unsorted_border--;
U32 t = itemsList [unsorted_border];
right_num_divItems--;
itemsList [right_num_divItems] = itemsList [left_num_divItems];
itemsList [left_num_divItems] = t;
}
else
left_num_divItems ++;
}


CODE LISTING 2.4.1

Intersection Boundary :: Intersect (Ray & ray) //ray.tfar is changed!
{
// Optimized inverse calculation (saves2 of 3 divisions)
float inv_tmp = (ray.d.dx * ray.d.dy) * ray.d.dz;
if ((uintCast (inv_tmp) &0x7FFFFFFF)> uintCast (DIV_EPSILON))
{
inv_tmp = 1.0f / inv_tmp;
ray.inv_d.dx = (ray.d.dy * ray.d.dz) * inv_tmp;
ray.inv_d.dy = (ray.d.dx * ray.d.dz) * inv_tmp;
ray.inv_d.dz = (ray.d.dx * ray.d.dy) * inv_tmp;
}
else
{
ray.inv_d.dx = ((uintCast (ray.d.dx) &0x7FFFFFFF)> uintCast (DIV_EPSILON))?
(1.0f / ray.d.dx): INVDIR_LUT [uintCast (ray.d.dx) >>31];
ray.inv_d.dy = ((uintCast (ray.d.dy) &0x7FFFFFFF)> uintCast (DIV_EPSILON))?
(1.0f / ray.d.dy): INVDIR_LUT [uintCast (ray.d.dy) >>31];
ray.inv_d.dz = ((uintCast (ray.d.dz) &0x7FFFFFFF)> uintCast (DIV_EPSILON))?
(1.0f / ray.d.dz): INVDIR_LUT [uintCast (ray.d.dz) >>31];
}
Intersection result;
result.tfar = ray.tfar;
result.tri_index = -1;
//
// BBox-Check
//
float tnear = 0.0f;
worldBBox.Clip (ray, tnear);
if (uintCast (ray.tfar) == 0x7F7843B0) // ray.tfar == 3.3e38f // !!
return (result);
//
U32 current_bspStack = 1; // wegen emptystack case == 0
U32 node = 0;
//
// BSP-Traversal
//
const U32 whatnode [3] = {(uintCast (ray.inv_d.dx) >> 27) & sizeof (BSPNODELEAF),
(uintCast (ray.inv_d.dy) >> 27) & sizeof (BSPNODELEAF),
(uintCast (ray.inv_d.dz) >> 27) & sizeof (BSPNODELEAF)};
U32 bspStackNode [128];
float bspStackFar [128];
float bspStackNear [128];
bspStackNear [0] = -3.4e38f; // sentinel
do
{
// Ist Node ein Leaf (type <0) oder nur neVerzweigung (type> = 0)
while (((BSPNODELEAF &) bspNodes [node]). type> = 0)
{
// Split-Dimension (x | y | z)
U32 proj = ((BSPNODELEAF &) bspNodes [node]). Type &3;
float distl = ((((BSPNODELEAF &) bspNodes [node]). splitlr [whatnode [proj] >> 4]
-ray.from [proj]) * ray.inv_d [proj];
float distr = ((((BSPNODELEAF &) bspNodes [node]). splitlr [(whatnode [proj] >> 4) ^ 1]
-ray.from [proj]) * ray.inv_d [proj];
node = (((BSPNODELEAF &) bspNodes [node]). type-proj) | whatnode [proj];
// type & 0xFFFFFFF0
if (tnear <= distl)
{
if (ray.tfar> = distr)
{
bspStackNear [current_bspStack] = MAX (tnear, distr);
bspStackNode [current_bspStack] = node ^ sizeof (BSPNODELEAF);
bspStackFar [current_bspStack] = ray.tfar;
current_bspStack ++;
}
ray.tfar = MIN (ray.tfar, distl);
}
else
if (ray.tfar> = distr)
{
tnear = MAX (tnear, distr);
node ^ = sizeof (BSPNODELEAF);

}
else
goto stackPop;
}
//
// Faces-Intersect
... code omitted ...
//
//
// Hit gefunden?
//
do // !! NEEDS bspStackNear [0] = -3.4e38f; !!!!!
{
stackPop:
current_bspStack--;
tnear = bspStackNear [current_bspStack];
} while (result.tfar <tnear);
if (current_bspStack == 0)
return (result);
node = bspStackNode [current_bspStack];
ray.tfar = bspStackFar [current_bspStack];
} while (true);
}

Claims (37)

A computer graphics system for generating pixel values of pixels in an image;
The pixel values represent points in the scene as recorded on the image plane of the simulated camera;
The computer graphics system generates the pixel value of an image using a selected ray tracing method that includes simulating at least one ray from the pixel into the scene along a selected direction. Configured to
The ray-tracing method further includes the steps of calculating the intersection of a ray with a surface of an object in the scene, and simulating a trajectory of the ray that illuminates the object in the scene, simulating the trajectory The computer graphics system , wherein the step uses a Markov chain,
Simulating the Markov chain using a quasi-Monte Carlo method,
Said step of simulating a Markov chain comprises sorting states;
A computer graphics system wherein the step of sorting includes proximity sorting . The computer graphics system of claim 1, comprising simultaneously simulating a plurality of Markov chains . The computer graphics system of claim 2, wherein the step of simultaneously simulating a plurality of Markov chains comprises utilizing a quasi-Monte Carlo point . 4. The computer graphics system of claim 3, wherein the quasi-Monte Carlo points are selected using any of a Halton sequence, a Halton sequence variant, a grid sequence, or a (t, s) -sequence . The computer graphics system of claim 2, wherein the step of simultaneously simulating a plurality of Markov chains utilizes an arbitrary number of chains . The computer graphics system of claim 2, wherein the step of simultaneously simulating a plurality of Markov chains further comprises the step of adding additional trajectories on the fly . The computer graphics system of claim 2, wherein the step of simultaneously simulating a plurality of Markov chains further comprises trajectory splitting . The computer graphics system of claim 2, wherein the step of simultaneously simulating a plurality of Markov chains comprises utilizing a technique for simulating particle absorption . 9. The computer graphics system of claim 8, wherein the step of simultaneously simulating multiple Markov chains further comprises using a Russian roulette method . For the s-dimensional problem, said step of simultaneously simulating multiple Markov chains
n: = 0
Initialize X0, l using quasi-Monte Carlo point xl Loop:
-Sort the state vector using the appropriate order-n: = n + 1
The computer graphics system according to claim 2, comprising the step of continuing the chain by calculating Xn, l using the next quasi-Monte Carlo point xl . The computer graphics system of claim 10, comprising ordering states by proximity in the state space . The computer graphics system of claim 11, comprising randomization . 13. The computer graphics system of claim 12, comprising utilizing any selected point set or point sequence with any selected randomization . Randomization and the following techniques:
n: = 0, select a sequence of points xi Instantiate randomization Initialize X0, l using rand (xl) Loop sort the state vector using some order n: = N + 1
12. The computer graphics system of claim 11, comprising: instantiating randomization; applying continuation of chain by calculating Xn, l using rand (xl) . Randomizing the quasi-Monte Carlo points at each time step n;
Said step of randomizing comprises:
Selecting a (t, s) -sequence of radix b = 2 from which subsequent samples are taken;
The computer graphics system of claim 11, comprising performing an XOR operation on the sample with an s-dimensional random vector, wherein the random vector is generated after each transition step . The computer graphics system of claim 15, comprising adding a selected number of additional trajectories on the fly . The computer graphics system of claim 16, comprising trajectory splitting . The computer graphics system of claim 17, comprising utilizing a technique for simulating particle absorption . The computer graphics system of claim 18, comprising using a Russian roulette method . Sorting by brightness or state norm;
Sorting using a spatial hierarchy;
And sorting using buckets enumerated by at least one space-filling curve . 21. The computer graphics system of claim 20, wherein sorting using a spatial hierarchy includes utilizing a binary hierarchy . The computer graphics system of claim 21, wherein the spatial hierarchy includes any of a BSP tree, a kD tree, or any other axis-parallel subset of a BSP tree structure, or a boundary volume hierarchy, or regular voxels . The binary hierarchy recursively subdivides using the plane selected by the selected heuristic,
21. The computer graphics system of claim 20, wherein the computer graphics system is then constructed by traversing the hierarchy in order to enumerate the leaves of the hierarchy in proximity order . 24. The computer graphics system of claim 23, comprising left balancing the tree and storing the tree in an array data structure to allow efficient traversal of the hierarchy . The spatial hierarchy includes leaves;
25. The computer graphics system of claim 24, wherein sorting using the spatial hierarchy further comprises utilizing bucket sorting and a selected space filling curve . 25. The computer graphics system of claim 24, wherein sorting using the spatial hierarchy further comprises using regular voxels . Delimiting the object to be rendered by an axis parallel bounding box;
By recursively or iteratively dividing the bounding box into successive left and right branches, each ending with a leaf, by applying a dividing plane that passes through a selected point of the object Generating a plurality of leaves;
Ordering the leaves hierarchically according to their brightness;
Performing a bucket sort on a matrix of a selected size divided into a selected number of buckets;
Using a selected space filling curve to follow within the matrix, wherein within each bucket a smaller space filling curve is used to follow within individual cells within the matrix. including the door, the computer graphics system of claim 25. The computer graphics system of claim 2, wherein the simulation is operable for a configuration for solving an integral equation . 30. The computer graphics system of claim 28, wherein the simulation is operable to allow light transport simulation to synthesize a real-looking image . The integral equation underlying the light transport simulation is reformulated as a path integral so that sampling the path space corresponds to simulating a Markov chain, where the path is a ray tracing and scattering event. Established by
30. The computer of claim 29, wherein the initial distribution is determined by the emission characteristics of the light source or sensor, and the transition probability is determined by a bidirectional reflection distribution function and a bidirectional subsurface scattering distribution function for the surface in the image. Graphics system . 31. The computer graphics system of claim 30, comprising processing a transmission effect by utilizing the path integral, initial distribution, and transition probability . 31. The computer graphics system of claim 30, wherein the path integral is solved by utilizing a high-dimensional quasi-Monte Carlo point or by padding a low-dimensional quasi-Monte Carlo point . 33. The computer graphics system of claim 32, wherein the integral equation underlying a light transport simulation is considered a Fredholm integral equation or a Volterra integral equation . A computer graphics system for generating pixel values of pixels in an image;
The pixel values represent points in the scene as recorded on the image plane of the simulated camera;
The computer graphics system generates the pixel value of an image using a selected ray tracing method that includes simulating at least one ray from the pixel into the scene along a selected direction. Configured to
The ray-tracing method further includes calculating intersections of rays and surfaces of objects in the scene; and simulating trajectories of rays that illuminate objects in the scene;
The method in the computer graphics system, wherein the step of simulating a trajectory uses a Markov chain,
Simulating the Markov chain using a quasi-Monte Carlo method,
Said step of simulating a Markov chain comprises sorting states;
The method wherein the step of sorting includes proximity sorting. A computer graphics system for generating pixel values of pixels in an image;
The pixel values represent points in the scene as recorded on the image plane of the simulated camera;
The computer graphics system generates the pixel value of an image using a selected ray tracing method that includes simulating at least one ray from the pixel into the scene along a selected direction. Configured to
The ray-tracing method further includes calculating intersections of rays and surfaces of objects in the scene; and simulating trajectories of rays that illuminate objects in the scene;
The step of simulating a trajectory is a system in the computer graphics system using a Markov chain,
Means for simulating the Markov chain using a quasi-Monte Carlo method;
Said means for simulating a Markov chain comprises means for sorting states;
The system wherein the means for sorting comprises means for proximity sorting. A computer graphics system for generating pixel values of pixels in an image;
The pixel values represent points in the scene as recorded on the image plane of the simulated camera;

The computer graphics system generates the pixel value of an image using a selected ray tracing method that includes simulating at least one ray from the pixel into the scene along a selected direction. Configured to
The ray-tracing method further includes calculating intersections of rays and surfaces of objects in the scene; and simulating trajectories of rays that illuminate objects in the scene;
The step of simulating a trajectory is a computer executable program that causes a computer to function as the computer graphics system using a Markov chain ;
In the computer,
Computer code means for simulating the Markov chain using a quasi-Monte Carlo method ;
Said computer code means for simulating a Markov chain, computer code means for sorting states ;
A computer-executable program for causing the computer code means for sorting to implement computer code means for proximity sorting. The computer graphics system of claim 19 including simulating an inhomogeneous Markov chain .

Images