Image synthesis by rank-1 lattices

[0001] The present application claims the benefit of U.S. Provisional Patent Application Ser. No. 60/693,232, filed on Jun. 23, 2005, which is hereby incorporated by reference in its entirety.

[0002] Submitted herewith is a source code listing, which is incorporated herein in its entirety. The source code listing is referred to herein as the “Appendix,” and is organized into sections identified by a reference number in the format “1.1.1” or “1.1.1.1.”

[0003] 1. Field of the Invention

[0004] The present invention relates generally to methods and systems for image synthesis in and by digital computing systems, such as for motion pictures and other applications, and in particular, relates to methods, systems, devices, computer software, and display technology for image synthesis by rank-1 lattices.

[0005] 2. Description of Prior Art

[0006] The use of synthetic images has become increasingly important and widespread in motion pictures and other commercial and scientific applications. A synthetic image typically represents a two-dimensional array of digital values, called picture elements or pixels, and thus is commonly regarded as a two-dimensional function.

[0007] Image synthesis is the process of creating synthetic images from scenes. As a general matter, digital images are generated by rasterization or, in the case of photorealistic images of three-dimensional scenes, by ray tracing. Both approaches aim at determining the appropriate color for each pixel by projecting the original function into the pixel basis. Due to the discrete representation of the original function, the problem of aliasing arises, as described below.

[0008] In particular, computing the intensity of a single pixel requires an integration function over the pixel area. This integral is often highly complex and cannot be solved analytically, thus requiring numerical methods for solution, which may include Monte Carlo and quasi-Monte Carlo methods. However, typical numerical methods used in such applications have their own limitations and attendant problems.

[0009] Another computational approach used in computer graphics is lattice theory. For example, the accumulation buffer typically used in computer graphics exploits the structure of regular grids. In addition, efficient rasterization techniques have been designed by taking advantage of such grids. Certain known sampling techniques employ Fibonacci lattices, which are a type of rank-1 lattice.

[0010] It would be desirable to provide methods, systems, devices, and computer software using rank-1 lattices that can be generated rapidly and efficiently, and that are readily adapted to the requirements of image synthesis and other computer graphics applications.

[0011] It would also be desirable to provide such methods, systems, devices and software that enable the efficient application of rank-1 lattices to anti-aliasing, rasterization, and other computer graphics functions.

[0012] Aspects of present invention provide a number of techniques in which rank-1 lattices are used in computerized image processing, in particular in the context of image synthesis.

[0013] A first aspect of the invention is directed to techniques used to select rank-1 lattices for use in image synthesis. A number of lattice selection techniques are described, including techniques for searching Korobov lattices in two dimensions by maximized minimum distance, thereby achieving maximum uniformity in the sense of (M, μ)-uniform point sets, which is used in quasi-Monte-Carlo theory. Further described are techniques for searching the lattice basis, and analyzing the relation of shifted rank-1 lattices to (t, m, s)-nets.

[0014] Another aspect of the invention provides a technique for structuring an image pixel layout using rank-1 lattices, and using the pixel layout as the basis of a rasterization scheme. The described rasterization scheme may be implemented using traditional raster displays. In addition, there is described herein a new display technology based on the rasterization scheme that directly takes advantage of the structure of rank-1 lattices.

[0015] A further aspect of the invention provides a technique for anti-aliasing by rank-1 lattices, in which maximum minimum distance lattices are compared with respect to their convergence to other Monte Carlo and quasi-Monte Carlo point sets. The problem of aliasing varies from pixel to pixel. Thus, the described aspect of the invention includes a technique for adaptive refinement by rank-1 lattices. Further described is a scheme for depth adaptive filtering in a ray tracing environment.

[0016] FIG. 1 is a schematic diagram of a conventional digital processing system in which the present invention can be deployed.

[0017] FIG. 2 is a schematic diagram of a conventional personal computer, or like computing apparatus, in which the present invention can be deployed.

[0018] FIGS. 3A and 3B are a pair of graphs respectively illustrating a Fibonacci lattice and its dual lattice.

[0019] FIG. 4 is a flowchart illustrating an overall method in accordance with a first aspect of the present invention.

[0020] FIG. 5 is a diagram of a Fibonacci lattice L(34, 21).

[0021] FIGS. 6A and 6B show, respectively, Korobov lattices L(8, 1) and L(8, 5).

[0022] FIG. 7 is a graph illustrating the computation of the distance of the point (2, 4) to the origin for the lattice L(5, 2).

[0023] FIG. 8A shows a fundamental parallelepiped of a hexagonal lattice, and FIG. 8B shows an equilateral triangle representing half of the fundamental parallelepiped.

[0024] FIGS. 9A-C illustrate a lattice selected, respectively, in the unit square, in the unit square scaled to an actual region, and in the actual region.

[0025] FIGS. 10A and 10B illustrate a scaled sampling lattice before and after shifting to obtain exactly one lattice point per pixel.

[0026] FIG. 11 shows a table setting forth max-min-dist lattices with valid shifts.

[0027] FIG. 12 shows a table setting forth lattices produced according to a search technique under the condition of a valid shift.

[0028] FIG. 13 shows a graph of a periodically extended lattice covered with a family of equidistant parallel lines.

[0029] FIG. 14 shows a graph illustrating a lattice with its Minkowsky reduced basis.

[0030] FIG. 15 shows a graph illustrating a (0, 2, 2)-net in basis 5.

[0031] FIG. 16 shows a table setting forth lattice parameters for some (0, 2, 2)-nets in base b.

[0032] FIG. 17 shows a flowchart of a rasterization method using rank-1lattices in accordance with a further aspect of the invention.

[0033] FIGS. 18A and 18B are a pair of graphs respectively illustrating intensification of a pixel and a line on a raster grid.

[0034] FIGS. 19A and 19B are a pair of graphs respectively illustrating intensification of a pixel and a line on a lattice grid.

[0035] FIG. 20 is a table setting forth single steps performed in the rasterization of a line on a lattice grid.

[0036] FIG. 21 is a table setting forth variables used in Crow's rasterization algorithm for convex polygons.

[0037] FIG. 22 is a graph illustrating a polygon rasterized on lattice L(576, 155), in a resolution of 24×24.

[0038] FIGS. 23A and 23B are diagrams illustrating rasterization on lattice grids at different resolutions.

[0039] FIG. 24 is a diagram illustrating the smoothing of a polygon by splitting cells of the lattice into two triangles.

[0040] FIG. 25 shows a raster display built using a maximum minimum distance lattice.

[0041] FIG. 26A shows a layout of area light sources for monochrome displays, and FIG. 26B shows a layout for color displays, in a system for anti-aliasing by rank-1 lattices.

[0042] FIG. 27 shows a Voroni diagram of the maximum minimum rank-1lattice L(64, 19).

[0043] FIG. 28A shows a layout of area light sources for hexagonal pixels in a monochrome display. FIG. 28B shows a layout for a color display.

[0044] FIGS. 29A and 29B show diagrams of a display comprising four modules, each of which contains 256 cells. FIG. 31A shows a quadratic layout of the individual modules, and FIG. 31B shows a rectangular layout.

[0045] FIG. 30 shows a diagram of a system in which display modules are addressed by a demultiplexer.

[0046] FIG. 31 shows a table setting forth parameters for 2ⁿ displays and images.

[0047] FIG. 32 shows a flowchart of a method for anti-aliasing by rank-1 lattices according to a further aspect of the invention.

[0048] FIG. 33 shows an image of an infinite checkerboard used to illustrate anti-aliasing by rank-1 lattices.

[0049] FIG. 34 shows a diagram of a pixel samples with a 4×4 grid.

[0050] FIG. 35A shows a uniform 4×4 pattern and FIG. 35B shows a jittered grid with one random sample in each of the 4×4 strata.

[0051] FIG. 36A is a graph illustrating the reconstruction of an original frequency with additional noise, and FIG. 36B is a graph illustrating the recovery of only noise.

[0052] FIG. 37A shows a Larcher-Pillichshammer point set for n=2⁴ points, and FIG. 38B shows a randomized Larcher-Pillichshammer point set.

[0053] FIG. 38A shows a max-min-dist lattice for n=16 and a=3, and FIG. 38B shows a randomized lattice with the same parameters.

[0054] FIG. 39A is a graph showing the first 5×5 pixels of a screen with a resolution of 512×512 covered by the max-min-dist lattice L(252144, 549). FIG. 39B shows the shifted lattice with shift (259, 0).

[0055] FIG. 40 shows a method for adaptive refinement by rank-1 lattices according to a further aspect of the invention.

[0056] FIG. 41 shows a max-min-dist lattice L(8,5) determined in the sheared basis of the lattice L(1048576, 33919).

[0057] FIG. 42 shows a refinement geometry for a gradient criterion.

[0058] FIG. 43 shows a placement diagram of a refinement max-min-dist lattice for a gradient criterion.

[0059] FIG. 44 shows a graph of a decision function for the B-spline degree.

[0060] FIG. 45 shows a diagram illustrating the influence of a sampling point for both a support of 1 pixel radius and 1.5 pixel radius.

[0061] FIG. 46 shows a graph illustrating interpolation between successive B-spline filters of degree d and d+1 for d=0.

[0062] FIG. 47A shows a diagram of a rendering of a sloped diagonal line without knowledge of the sub-pixel structure. FIG. 47B shows the rendering after tripling the horizontal resolution.

[0063] FIG. 48 shows a diagram of an RGB decimation using one-pixel-wide box filters.

[0064] FIGS. 49-51 are a series of drawings illustrating lattices used to render an infinite checkerboard image.

[0065] The present invention relates to the generation and synthesis of images, such as for display in a motion picture or other dynamic display, using rank-1 lattices. The techniques described herein are practiced as part of a computer graphics system, in which a pixel value is generated for each pixel in an image. The pixel value is representative of a point in a scene as recorded on an image plane of a simulated camera. The computer graphics system is configured to generate the pixel value for an image using a selected methodology.

[0066] The following discussion describes methods, structures, systems, and display technology in accordance with these techniques. It will be understood by those skilled in the art that the described methods and systems can be implemented in software, hardware, or a combination of software and hardware, using conventional computer apparatus such as a personal computer (PC) or equivalent device operating in accordance with (or emulating) a conventional operating system such as Microsoft Windows, Linux, or Unix, either in a standalone configuration or across a network. The various processing means and computational means described below and recited in the claims may therefore be implemented in the software and/or hardware elements of a properly configured digital processing device or network of devices. Processing may be performed sequentially or in parallel, and may be implemented using special purpose or reconfigurable hardware.

[0067] Methods, devices or software products in accordance with the invention can operate on any of a wide range of conventional computing devices and systems, such as those depicted by way of example in FIG. 1 (e.g., network system 100), whether standalone, networked, portable or fixed, including conventional PCs 102, laptops 104, handheld or mobile computers 106, or across the Internet or other networks 108, which may in turn include servers 110 and storage 112.

[0068] In line with conventional computer software and hardware practice, a software application configured in accordance with the invention can operate within, e.g., a PC 102 like that shown in FIG. 2, in which program instructions can be read from CD ROM 116, magnetic disk or other storage 120 and loaded into RAM 114 for execution by CPU 118. Data can be input into the system via any known device or means, including a conventional keyboard, scanner, mouse or other elements 103.

[0069] The present discussion is divided into the following sections:

• • [0070] 1. Introduction
  • [0071] 2. Selection of Rank-1 Lattices
  • [0072] 3. Rasterization on Rank-1 Lattices
  • [0073] 4. Anti-Aliasing by Rank-1 Lattices
  • [0074] 5. Adaptive Refinement and Filtering by Rank-1 Lattices
  • [0075] 6. Conclusion

[0076] It is noted that these section headings are used in order to provide organization and clarity to the following discussion. It is not intended for the use of headings to define or otherwise limit the scope of the invention.

[0077] A synthetic image represents a two-dimensional array of digital values, called picture elements or “pixels” and thus can be regarded as a two-dimensional function. As used herein, the term “image synthesis” refers to the process of creating synthetic images from scenes.

[0078] In general, digital images are generated by a process known as “rasterization,” or in the case of photorealistic images of three dimensional scenes, by a process known as “ray tracing.” Both approaches aim at determining the appropriate color for each pixel by projecting the original function into the pixel basis. Due to the discrete representation of the original function the problem of aliasing arises, discussed below.

[0079] Computing the intensity I(x) of a single pixel P requires the following integration over the pixel area: $\mathrm{pixel}\mathrm{color}=I\left(x\right)=\frac{1}{P}{\int }_{P}f\left(x\right)dx.$

[0080] Typically, this integral is quite complicated and cannot be solved analytically. Therefore numerical methods have to be used. As the integrands are very likely to be discontinuous, the standard tensor product quadrature does not work well, especially for higher dimensions. The square integrability of the underlying functions allows for computing these integrals by means of the so-called “Monte Carlo method,” independent of dimension. This approach uses probability theory for approximating the integral $\begin{array}{cc}{\int }_{I^s}f\left(x\right)dx,& \left(1.2\right)\end{array}$
where I^S=[0,1)^S represents the s-dimensional unit cube. The estimation of the above integral is based on constructing a random variable Y, a so called “estimator,” whose expectation EY corresponds to the desired integral. Choosing $Y:=\frac{f\left(X\right)}{p\left(X\right)},$
for some random variable X with probability density function $p\text{:}{I}^s\to {ℝ}_0^{+}$ $\mathrm{yields}$ $\mathrm{EY}={\int }_{I^s}\frac{f\left(x\right)}{p\left(x\right)}p\left(x\right)dx={\int }_{I^s}f\left(x\right)dx$
as the desired result. Considering n independent identically distributed (i.i.d.) variables Y_i, the central limit theorem states that the sum $\sum _{i=0}^{n-1}{Y}_i$
becomes asymptotically
N(n·EY,n·VY)
distributed for large n. Then the 3σ rule $P\left(\left\{n·\mathrm{EY}-3\mathrm{VY}\sqrt{n}<\sum _{i=0}^{n-1}{Y}_i<\mathrm{nEY}+3\mathrm{VY}\sqrt{n}\right\}\right)=P\left(\left\{\sum _{i=0}^{n-1}{Y}_i-\mathrm{EY}<\frac{3\mathrm{VY}}{\sqrt{n}}\right\}\right)\approx 0.997$
gives a method of computation for the integral 1.2 by $\mathrm{EY}={\int }_{I^s}f\left(x\right)\approx \sum _{i=0}^{n-1}\frac{f\left(x_i\right)}{p\left(x_i\right)}$
This means that the Monte Carlo method performs numerical integration by random sampling and averaging. According to the 3σ rule the approximation error is below $\frac{3\mathrm{VY}}{\sqrt{n}}$
with high probability for large n with the convergence rate being $O\left(\frac{1}{\sqrt{n}}\right).$
Let t_Y be the time spent for one realization y_i of Y. As the computational cost C_Y=t_YVY of the Monte Carlo algorithm depends on the variance of the estimator, there exist various variance reduction techniques, like separation of the main part or jittering to increase efficiency. Since there are no assumptions about smoothness or the dimension of the integrand, the Monte Carlo algorithm overcomes the problems of the standard tensor product quadrature. However, this method has only a slow convergence and there are no real random numbers on computers anyway. The quasi-Monte Carlo method exploits the fact that the independence of the sampling points being essential for the error bound of the Monte-Carlo method is not required for convergence. Instead the uniform distribution of the sampling pattern represents the crucial factor. Whereas independence, unpredictability and uniformity are characteristic for real random numbers on the unit interval, the so-called quasi-Monte Carlo points are deterministically generated by discarding the first two properties. These deterministic points are much more uniform than is possible for real random numbers. Examples for such constructions are the Hammersley point set or (t, m, s)-nets, which are based on radical inversion and rank-1 lattices, discussed below.

[0081] These point sets share the property of partitioning the unit cube in a uniform way. This kind of uniformity has been formalized as follows:

[0082] Definition 1 Let (X, B, μ) be an arbitrary probability space and let M be a nonempty subset of B. A point set P_n of n elements of X is called (M, μ)-uniform if $\sum _{i=0}^{n-1}{x}_M\left(x_i\right)=\mu \left(M\right)·n\mathrm{for}\mathrm{all}M\in M,$
where _XM(x_i)=1 if x_iεM, zero otherwise.
Known examples for (M, μ)-uniform point sets include radical inversion based points, points from the Cartesian product midpoint rule and rank-1 lattices. Using the probability space
([0, 1)^s,B,λ_s),
where B corresponds to the Borel-sets and λ_S to the s-dimensional Lebesgue-measure, this definition also applies in the context of computer graphics.

[0083] The quasi-Monte Carlo method principally looks like the Monte Carlo method, but uses deterministic points for numerical integration. This means that integral 1.2 is approximated as follows: ${\int }_{I^s}f\left(x\right)dx\approx \frac{1}{n}\sum _{i=0}^{n-1}f\left(x_i\right).$
By means of the (M, μ)-uniformity a deterministic error bound on L² can be obtained as follows:

[0084] Theorem 1 Let (X, B,μ) be an arbitrary probability space and let M={M₁ . . . M_k} be a partition of X with M_jεB for 1≦j≦k. Then for any (M, μ)-uniform point set P={X₀ . . . X_n-1} and any function ƒ on X we have ${\int }_{X}f\left(x\right)d\mu \left(x\right)-\frac{1}{n}\sum _{i=0}^{n-1}f\left(x_i\right)\le \sum _{j=1}^k\mu \left(M_j\right)\left({\sup }_{x\in M_j}f\left(x\right)-{\inf }_{x\in M_j}f\left(x\right)\right).$

[0085] The above discussion applies to the intensity of a single pixel. Simultaneously computing all pixels $g\left(y\right):={\int }_{I^s}f\left(x,y\right)dx$
is an integro-approximation problem, typically addressed by generalizing to a parametric integration. The generalization of the Monte Carlo method for random fields yields the so-called “method of dependent tests.” A whole function is estimated as follows: $g\left(y\right):={\int }_{I^s}f\left(x,y\right)dx\approx \frac{1}{n}\sum _{i=0}^{n-1}f\left(x_i,y\right),$
with the computational cost being $n·t_s·E{{\int }_{I^s}f\left(x,y\right)dx-\frac{1}{n}\sum _{i=0}^{n-1}f\left(x_i,y\right)}_2^2.$
Considering the quasi-Monte Carlo method, Theorem 1 has been generalized to parametric integration:

[0086] Theorem 2 Let (X, B, μ) be an arbitrary probability space and let M={M₁, . . . M_k} be a partition of X with M_jεB for 1≦j≦k. Then for any (M, μ)-uniform point set P={X₀ . . . X_n-1} and any function ƒ on X we have $\frac{1}{n}\sum _{i=0}^{n-1}f\left(x_i,y\right)-{\int }_{X}f\left(x_i,y\right)d\mu \left(x\right)\le \sum _{j=1}^k\mu \left(M_j\right)\left({\sup }_{x\in M_j}f\left(x,y\right)-{\inf }_{x\in M_j}f\left(x,y\right)\right).$
A lattice L is a discrete subset in ^s
L:=P_n+⊂^s
which is closed under addition and subtraction. For example, regular axis-aligned grids represent one type of simple lattices. Every lattice in ^s can be generated from a set of linear independent vectors {b₁, . . . , b_s} by considering all integer linear combinations of b₁, . . . , b_s, where b_jε^s for j=1, . . . , s $L\left(b_1,\dots ,b_s\right)=\left\{\sum _{j=1}^s{\lambda }_j{b}_j\text{:}{\lambda }_1,\dots ,{\lambda }_s\in \mathbb{Z}\right\}.$
The vectors {b₁, . . . , b_s} are one basis of the lattice. As they generate the lattice, they are also said to be generator vectors. The number s of basis vectors yields the dimension of L.

[0087] The absolute value of the determinant of the basis or rather “generator” matrix B whose row vectors are formed by b₁^τ, . . . , b_s^τ corresponds to the s-dimensional volume of the fundamental parallelepiped
Λ=Λ(b₁, . . . ,b_s)={γ₁b₁+ . . . +γ_sb_s:0≦γ_i≦1,1≦i≦s}
which is spanned by the basis vectors and the origin. By means of the fundamental parallelepiped ^s can be partitioned into lattice cells of the same volume and orientation:
{x+Λ}={yε^s |y=x+z, where zεΛ} for xεL

[0088] While the basis is not unique for a specific lattice, its determinant det L=|det B| represents an invariant quantity. This means that a set of linear independent vectors generates the lattice L if the determinant of the corresponding basis matrix equals the lattice determinant.

[0089] The row vectors of (B⁻¹)^τ build a basis of the dual, or polar, lattice of L, which is defined by
L^⊥:={hε^s:h·xε∀xεL}
The determinant of L^⊥ arises from the reciprocal of the determinant of L. The following interpretation of the dual lattice L^⊥ is useful for the spectral test, discussed below. Every point h in the dual lattice L^⊥ defines a set of (s−1) dimensional hyperplanes
{x:h·x=δ}, for δ=0, ±1, ±2, . . . ,
so that the following properties hold: Every point in L lies on one of the hyperplanes and every hyperplane contains at least one point of L. These hyperplanes are parallel and equidistant. Their distance equals
∥h∥₂⁻¹.
FIG. 3A shows a Fibonacci lattice 200 containing the family of hyperplanes of the dual lattice point h=(1, 2), and FIG. 3B shows its dual lattice 210.

[0090] Typical problems in lattice theory are the closest-point and the shortest vector problem, discussed below. The closest-point problem arises, for example, in communication theory where lattices are applied both for modulation and quantization. The shortest vector problem has been used, for example, in order to measure the quality of random number generators. Another application of this problem is found in cryptography. A further problem, discussed below, is that of lattice basis reduction and the related search for a lattice basis. Other lattice problems include determining the kissing number or finding the Voronoi-relevant vectors.

[0091] Aspects of the present invention provide techniques for using rank-1 lattices in image synthesis. FIG. 4 is a flowchart illustrating a method 220 according to the present invention. In step 221, a rank-1 lattice is selected in accordance with a maximized minimum distance function corresponding to pixels in a synthetic image to be generated. In step 222, pixel values are calculated in accordance with the selected rank-1 lattice. In step 223, a display-controlling electronic output signal is generated, based on the calculated pixel values.

[0092] Lattices are now discussed in greater detail, particularly techniques for selecting rank-1 lattices, including techniques for searching Korobov lattices in two dimensions by maximized minimum distance, thereby achieving maximum uniformity in the sense of (M, μ)-uniform point sets.

[0093] Some applications of rank-1 lattices in accordance with other aspects of the invention described herein require a “lattice basis,” i.e., a set of linearly independent vectors capable of generating the subspace defined by a selected rank-1 lattice. The following discussion therefore includes a technique for searching the lattice basis. In addition, the relation of shifted rank-1 lattices to (t, m, s)-nets is analyzed.

[0094] The “rank” of a lattice L in ^S is defined to be the minimum number of vectors which are necessary to generate L. Given only one s-dimensional generator vector gε^S, $x_i=\frac{1}{n}·g=\frac{1}{n}\left(g_1,\dots ,g_s\right)$
describes an s-dimensional lattice. Thus, such lattices are known as “rank-1-lattices.”

[0095] “Korobov lattices” provide an example of s-dimensional rank-1-lattices. Their generator vector has the form g=(1, a, a², . . . , a^s−1). n being the number of points, it follows that a Korobov lattice is uniquely determined by the tuple (n, a): $L\left(n,a\right)=\left\{\frac{i}{n}\left(1,a,a^2,\dots ,a^{s-1}\right)\mathrm{mod}1|0\le i<n\right\}$
Thereby, the mod 1 operation restricts the lattice to the unit square resulting in a one-periodic pattern. This means that rank-1 lattices are point symmetric to the origin. One advantage of such lattices is that they can be tiled seamlessly in order to fill the s-dimensional space and that they can be generated very quickly.

[0096] The “Fibonacci lattices,” which represent a two-dimensional Korobov lattice, are another example of rank-1-lattices. Based on the Fibonacci sequence F_k:=F_k-1+F_k-2 with F₂:=F₁:=1, the number of points is set to n=F_k, k≧2 and the generator vector is defined as
g=(1,F_k-1)

[0097] FIG. 5 shows a Fibonacci lattice 230 on the unit square with n=F₉=34 points and generator vector g=(1, F₈). The quality of rank-1 lattices is significantly influenced by the integer generator vector (g₁, . . . , g_s)ε^S. Lattices in Korobov form, which do not differ very much from other lattices in two dimensions, reduce the choice of the generator vector to one parameter. Therefore these lattices are of particular interest in the present discussion with respect to the search for two-dimensional lattices with predefined properties.

[0098] As described in detail below, rank-1-lattices are applied in the contexts of anti-aliasing and rasterization. The present discussion addresses the question as to what kind of rank-1-lattices to use.

[0099] The generator vector g of a rank-1-lattice can be chosen such that the resulting point set is of low discrepancy. In that case, the elements of this point set are called “good lattice points.” For example, Fibonacci lattices for s=2 and lattices whose generator vector is in Korobov form are of low discrepancy. However, few explicit constructions for good lattice points exist.

[0100] As mentioned above, rank-1 lattices are (M, μ)-uniform point sets. This property of rank-1 lattices can be seen by considering the partitions {x+Λ} of I^S which are induced by the fundamental parallelepiped. Maximum uniformity in this sense can be achieved by maximizing the mutual minimum distance
d_min(P_n):=min_0≦i<nmin_0≦j<n∥x_j−x_i∥2
among all sample points P_n. Thereby, the Euclidean norm is calculated on the unit torus.

[0101] Hence, rank-1-lattices are used herein in the Korobov form: $\begin{array}{cc}{x}_i=\left(\frac{i}{n},\frac{a·i}{n}\mathrm{mod}1\right)\in {\left[0,1\right)}^2,\mathrm{for}0\le i<n& \left(2.3\right)\end{array}$
featuring a so-called “max-min-dist-property” and identify them by search. Since a lattice in Korobov form is uniquely determined by the tuple (n, a), this means that for a fixed number of points n the search is conducted for such an aε{1, . . . , n−1} that the minimum distance between two lattice points achieves the maximum on the unit torus. The more the minimum distance between two lattice points grows, the more the resulting lattice approximates a hexagonal lattice, thereby increasing the uniform distribution of the point set.

[0102] FIG. 6A shows a Korobov lattice 240 for L(8, 1), in which a=8 and n=1, and FIG. 6B shows a Korobov lattice 250 for L(8, 5), in which a=8 and n=5. For the FIG. 6A lattice 240, the points lie on the diagonal of the unit square, whereas the points of the max-min-dist lattice shown in the FIG. 6B lattice 240 are more uniformly distributed.

[0103] The task of calculating the minimum distance between two lattice points is related to a well known problem in lattice theory. The “closest-point problem” consists in finding a vector c for a given lattice L and a given input point xε^S such that
∥x−c∥≦∥x−y∥∀yεL
Calculating the minimum distance in a lattice L corresponds to determining the shortest vector in L which can be accomplished by a simple modification of the closest-point problem. For that purpose, only vectors v≠0 are considered. Then the closest point to the origin yields the shortest vector, that is, the minimum distance in L. Closest-point search methods have been described in the literature for lattices without a regular structure. But the methods presented there are not appropriate to rank-1 lattices, as they provide a special structure and are only considered for low dimensions. In the following discussion, several approaches for lattice searching are described that exploit both the periodic and regular structure and the application of these lattices.

[0104] There is now described a technique for searching for rank-1 lattices on the unit square. Section 2.1.1 of the Appendix, sets forth a code listing for a basic search on the unit square which, for a given nε, considers all possible lattices (n, a), 2≦a<n−1. Values for a>n do not need to be taken into account because they produce the lattices (n, a), 1≦a<n due to the modulo operation. Furthermore, assuming n>3, a=1 and a=n−1 are excluded from search, since the lattice points lie on the two diagonals of the unit square. a=0 would cause the lattice to collapse to the origin and therefore is omitted as well.

[0105] The choice of the parameter a for a max-min-dist lattice is not unique, that is, there exist several possibilities that result in the same maximum minimum distance. Therefore, a set a_opt is used in order to store the set
a_opt:={a|max_a{min_∥x∥{∥y∥,yεL(n,a)}}, 2≦a<n−1}
found so far.

[0106] For each lattice, the algorithm determines the length of the shortest vector by looping over all lattice points x_i, 1≦i<n. Although according to equation 2.3 the lattice is defined on the unit torus, all computations in the § 2.1.1 code listing are performed in integer arithmetic in order to avoid rounding errors due to floating point imprecision. This means that the present discussion focuses on integer lattices on [0, n)² instead of [0, 1)²:
x_i=(i,a·i mod n)  (2.4)

[0107] Due to the modulo operation the lattice L(n, a)⊂[0, n)² is periodic. In order to find the length of the shortest vector, there is now computed the distance from the lattice point x_i to the origin. Calculating on the unit torus, it is necessary to consider the point x_i from a different quadrant if i>n−i or a·i>n−a·i mod n: $\begin{array}{cc}\mathrm{mindist}\left(x_i\right)={\tilde{x}}_i=\{\begin{array}{cc}\left(i,a·i\right)\mathrm{mod}n& \mathrm{if}i\le n-i\mathrm{and}a·i\le n-a·i\mathrm{mod}n\\ \left(i,a·i-n\right)\mathrm{mod}n& \mathrm{if}i\le n-i\mathrm{and}a·i>n-a·i\mathrm{mod}n\\ \left(i-n,a·i\right)\mathrm{mod}n& \mathrm{if}i>n-i\mathrm{and}a·i\le n-a·i\mathrm{mod}n\\ \left(i-n,a·i-n\right)\mathrm{mod}n& \mathrm{if}i>n-i\mathrm{and}a·i>n-a·i\mathrm{mod}n\end{array}& \left(2.5\right)\end{array}$

[0108] FIG. 7 shows a graph 260 illustrating the above technique for the lattice L(5, 2). In order to compute the distance of the point x₂=(2, 4) to the origin, due to the periodicity, the point (0, 5) is taken as the new origin. Using the new origin corresponds to calculating the Euclidean norm of the point x₂=(2, 4−5) relative to the origin of the first quadrant at the point (0, 0).

[0109] The distance of the lattice point {tilde over (x)}_i to the origin is computed by means of the squared Euclidean norm ∥x∥₂²=x₁²+x₂². Both the squares and their sum might produce an integer overflow during computation which is handled in the following way:

• • [0110] The data type for the variables is either chosen to be unsigned int (UINT), especially for the tuple (n,a), or unsigned long long int (ULLINT) for the entities involved in the Euclidean norm calculation.
  • [0111] An integer overflow for x_j² will occur if x_j>√{square root over (MAX_ULLINT)}=0xFFFFFFFF.
  • [0112] As x+y≧x and x+y≧y holds true for the sum of two elements and
    MAX_—ULLINT+1=0, z+MAX_—ULLINT=z−1
  • [0113] an overflow is evident if x+y<x or x+y<y.
    In the last two cases, the § 2.2.1 code proceeds to the next lattice point. The multiplication a·i mod n in equation 2.4 might also be problematic with respect to overflow. Therefore it is resolved into the incremental sum $a·i=\sum _{i=0}^{n-1}a$
    being computed iteratively and the modulo operation is applied in each iteration of the inner for-loop by $\mathrm{var}\mathrm{mod}n=\{\begin{array}{cc}\mathrm{var}& \mathrm{if}\mathrm{var}<n\\ \mathrm{var}-n& \mathrm{else}\end{array}$
    Once the shortest vector is found, its length is compared to the length of the vector having maximum minimum distance maxmindist so far and possibly remember the related parameter a. But this is not done until the Boolean flag noOverflow is set true the first time in order to assure the new maxmindist not to become MAX_ULLINT which is the initial value for mindist in the § 2.2.1 code listing. Depending on whether the new maxmindist is greater or equal to the old one the vector a_opt is cleared before storing the new value.

[0114] The complexity of the § 2.2.1 code listing is O(n²) which turns out to be far too slow when the lattice search is executed for large n. Hence, it is desirable to restrain the search domain. In the first place, it is possible to further reduce the number of a's being examined. If a=â is a valid solution to the search problem, so is a=n−â having the same max-min-dist because of symmetry. Thus it is sufficient to range a from 2 to n/2.

[0115] In addition, it is not necessary to loop over all lattice points in order to determine the length of the shortest vector because of the “sphere packing” argument. “Sphere packing” means to find an arrangement of spheres, or circles in two dimensions, which do not overlap each other with some or all of them being mutually tangent. In the present case, the circles are centered at the lattice points, and their radius is chosen such that they just touch each other. This is accomplished by r=mindist/2. Hence, finding the maximum minimum distance is a sphere packing problem. In terms of sphere packing, it is necessary to maximize the ratio of the volume of the sphere to the volume of the fundamental parallelepiped (in two dimensions) $\frac{r^2\pi }{\det B}=\frac{{\mathrm{mindist}}^2\pi }{4\det B},$
where B is the basis matrix. As the hexagonal lattice solves the sphere packing problem in ², the largest maximum minimum distance is achieved in this type of lattice. Therefore it is not necessary to examine lattice points which are outside the circle around the origin with radius l, where l is the distance between two points in the hexagonal lattice. This quantity can be evaluated in the following way:

[0116] Each lattice cell, i.e., the fundamental parallelepiped, is divided into two equilateral triangles. The lattice covers an area of n². Hence each of the n lattice cells (one for each point) has an area of $\det \left(B\right)==\frac{\mathrm{area}}{\mathrm{number}\mathrm{of}\mathrm{lattice}\mathrm{points}}=\frac{n^2}{n}=n.$
So the area A of an equilateral triangle is n/2. Now it is possible to calculate the length l of the triangle sides which is the sought quantity: $A=\frac{h·l}{2}=\frac{\sqrt{l^2-{\left(\frac{l}{2}\right)}^2}·l}{2}=\frac{n}{2}⟹l=\sqrt{\sqrt{\frac{4}{3}}·n}$

[0117] FIGS. 8A and 8B are a pair of diagrams summarizing the above-described process. FIG. 8A shows a fundamental parallelepiped 270 of a hexagonal lattice. As shown in FIG. 8B, one half of the fundamental parallelepiped represents an equilateral triangle 280. h indicates the height of the triangle. l represents the length of one triangle side.

[0118] Thus, the outer for-loop in the § 2.2.1 code listing can be stopped at $\sqrt{\sqrt{\frac{4}{3}}·n},i.e1\le i\le \sqrt{\sqrt{\frac{4}{3}}·n}.$
In combination, complexity of the basic search algorithm has been reduced from O(n²) to $O\left(\frac{n}{2}·\sqrt{\sqrt{\frac{4}{3}}·n}\right)=O\left(n^{1.5}\right).$

[0119] Thus far, no conditions have been imposed on the lattice parameter a apart from maximizing the minimum distance in L(n, a). If we demand that n and a are relatively prime, which means that their greatest common divisor equals 1, that is
gcd(n,a)=1,
the resulting lattice points will be stratified equidistantly on each coordinate axis. Hence the minimum distance can be bounded to $\mathrm{mindist}\ge \frac{1}{n}.$
However, this condition represents a restriction in regard to finding the maximum minimum distance lattice for a given nε.

[0120] For some n the basic search algorithm yields quite a huge number of max-min-dist lattices. In order to reduce the number of possible lattices, the § 2.2.1 code listing is extended such that also the distance of the second shortest vector is maximized. In the present case, considering the third shortest vector brings no further improvement.

[0121] There is now discussed the use of weighted norms. The discussion up to this point has focused on the search for max-min-dist lattices on the unit square. The resulting lattice points can be applied for numerical integration on a quadratic domain, for example. In the technique described further below, lattices are used covering an entire computer screen, which is usually not scaled quadratically. But simply scaling a lattice ε[0, 1)² to a rectangular region may not result in a uniform distribution of the points in the sense of (M, μ)-uniformity. Non-uniform distribution is illustrated in FIGS. 9A-C. FIG. 9A shows a graph 290 illustrating the lattice with n=786432 points, selected in the unit square. FIG. 9B shows a graph 300 illustrating the lattice in the unit square, scaled to the actual region. FIG. 9C shows a graph 310 illustrating the lattice in the actual region.

[0122] It is thus necessary to account for the scale of the actual region when calculating the distance of the lattice points to the origin. For that purpose, the Euclidean norm is replaced by a weighted norm. Therefore, the described technique is conducted as if on the unit square, but with the squared Euclidean norm
∥z∥₂²=z_x²+z_y²
modified by two scaling factors in order to make up for the region size. These scaling factors are derived by dividing both the x- and y-Resolution through their greatest common divisor (gcd) ${z}_2^2={\lambda }_x^2·z_x^2+{\lambda }_y^2·z_y^2,\mathrm{where}$ ${\lambda }_x=\frac{\mathrm{xRes}}{\gcd \left(\mathrm{xRes},\mathrm{yRes}\right)}$ $\mathrm{and}$ ${\lambda }_y=\frac{\mathrm{yRes}}{\gcd \left(\mathrm{xRes},\mathrm{yRes}\right)}$
The introduction of scaling factors brings about only a small modification of the basic search algorithm. Section 2.1.2 of the Appendix shows a code listing including an extended search considering the actual region size. The § 2.1.2 code listing incorporates the search on an arbitrary region size as well as the improvements discussed above. The dots in the § 2.1.2 code listing represent the unchanged part of the § 2.1.1 code listing.

[0123] As mentioned above, one use of the described lattices is to sample an entire display screen. Assuming a screen resolution of xRes×yRes and a number of lattice points n=xRes·yRes, it is possible to encounter the situation of not sampling every pixel, that is, the situation that there are pixels left without any sample. FIG. 10A shows a graph 320 illustrating this situation. In FIG. 10A, the lattice L(25, 7)⊂[0,1)² is scaled to [0,5)² and samples the screen with a resolution of 5×5. Pixel (0, 1) is missed, while pixel (0, 0) contains two samples.

[0124] Similarly, the same number of samples per pixel will not always be achieved, when the number of lattice points is set to n=xRes·yRes·α with α≧1. More generally, it is desirable to obtain the same number of points εL(n, a)⊂[0,1)² per interval $E:=\prod _{j=0}^1\left[\frac{a_j}{\mathrm{xRes}},\frac{a_j+1}{\mathrm{yRes}}\right).\mathrm{for}\mathrm{integers}0\le a_0<\mathrm{xRes},0\le a_1<\mathrm{yRes},$
where xRes and yRes refer to the screen resolution.

[0125] There are a number of solutions to this problem. The solutions have in common the search for a valid shift of the whole lattice. In this context, “valid” means that for the shifted lattice, there is an equal number of lattice points per interval being described in the above equations. In order to prevent rounding errors, calculations will again be performed on [0, n)² instead of on [0, 1)².

[0126] Due to the periodicity of the lattice it is sufficient to consider shifts of maximally one pixel width which is easier than searching in the fundamental parallelepiped. The search domain, i.e., one pixel, is discretized by looking at the low dimensional projections of the lattice points. If these are perfect in the sense of equally stratifying each coordinate axis, the distance between two lattice points is xRes/n on the x-axis and yRes/n on the y-axis, respectively. This means that there are $\mathrm{numX}={\left(\frac{\mathrm{xRes}}{n}\right)}^{-1}=\frac{n}{\mathrm{xRes}}$
projected lattice points on the x-axis and $\mathrm{numY}={\left(\frac{\mathrm{yRes}}{n}\right)}^{-1}=\frac{n}{\mathrm{yRes}}$
projected lattice points on the y-axis within one pixel. The discretization provided by those projections is highlighted by the small grid at the lower left corner of the FIG. 10A graph 320, which represents the shifts to be examined. In the following we discuss three potential search strategies.

[0127] Ideally for a given max-min-dist lattice L(n, a) we find a valid shift s=(s_x, s_y). For that purpose all possible shifts are looped over and the number of lattice points per pixel are counted. After that, the number of points in pixel (0, 0) is compared to the other count buffer entries. As soon as one of them does not equal the number of lattice points in the first pixel, comparing can be aborted, and a switch is made to the next lattice shift. The procedure stops as soon as a valid shift is found and returns the shift on [0, n)². Section 2.1.3 of the Appendix shows a code listing of a search shift for a given lattice.

[0128] FIG. 10B shows a graph 330, in which a shift of (10/25,10/25) has been added to the lattice L(25, 7) before scaling it from [0, 1)² to [0, 5)². After the shift is added, there is exactly one lattice point per pixel. However, lattice shifts can be found only for very few n, and the bigger n grows, the less probable it is that there will be a valid shift. FIG. 11 shows a table 340 listing some maximum minimum distance lattices, for which valid shifts have been found, ad well as the underlying screen resolution.

[0129] Thus, a different way to search is needed. One possibility is to first search for a set of a's for which the resulting lattices contain at least one valid shift and then perform the max-min-dist search among these candidates. But since there is sought a lattice
{tilde over (L)}({circumflex over (n)},a)
which approximates the max-min-dist lattice
L({circumflex over (n)},a)
best for a given
n={circumflex over (n)}ε
the following strategy turns out to be significantly better. In a first step the list of a's
[a|2≦a<n−1]
is sorted with respect to the minimum distance in L(n, a) using heap sort. After that, the sorted list is passed through in descending order and the shift search is performed for each a. Thereby, the first a with a valid shift yields the desired result. The implementation of this search strategy is incorporated in the code listing set forth in section 2.1.4 of the Appendix.

[0130] In the § 2.1.4 code listing, the number of lattice points per pixel is counted again. If one pixel holds more than n/(xRes·yRes) points, the process can switch over to the next lattice shift at once because in this case one pixel will surely contain less than n/(xRes·yRes) samples. The dots in the § 2.1.4 code listing represent unchanged parts of the previous code listings.

[0131] Due to the three nested for-loops, the § 2.1.4 code listing has a complexity of O(n³) and is very slow for large n. But the search for shifts can be shortened considerably by exchanging loop (1) and (2). That means that the search domain is passed through in a different order. This is based on the observation that the y-component of at least one shift is zero in all examined cases. FIG. 12 sets forth a table 350 listing lattices that result from the § 2.1.4 code listing. In the FIG. 12 table, shifts which have a non-zero y-component were obtained by not stopping the shift search after the first shift was found. It took approximately five days to compute the lattice L(4194304, 182183) on a Intel Pentium 4, 2.8 GHz, despite switching loops (1) and (2).

[0132] Though achieving the same number of lattice points per pixel, it should be borne in mind that the resulting lattices mostly will not satisfy the maximum minimum distance property, but rather will approximate it as well as possible.

[0133] There is now discussed the use of the so-called “spectral test” in the selection of rank-1 lattices. The “spectral test” measures the quality of linear congruential random number generators which are defined by a tuple (a, c, m, X₀) and a recurrence formula
X_n+1=(aX_n+c)mod m, n≧0
where mε is the modulo, aε{1, . . . , m−1} is the multiplier, cε{0, . . . m−1} is the increment and X₀ε{0, . . . , m−1} is the starting value. Numbers between 0 and 1 are given by
u_n=X_n/m.

[0134] Based on a sequence U_n of period m, for all t>0 the spectral test analyzes the set of all m points
{(u_n,u_n+1, . . . ,u_n+t−1)|0≦n<m}
which is formed by successive overlapping t-tuples of {u₀, u₁, . . . } thereby considering all possible starting values. To simplify theory, the spectral test makes the assumption that there is either a linear congruential sequence of full period m with c≠0 or that m is prime and c=0 and the period length is m−1. In order to assure that there are always m points, the origin is added to the point set in the latter case because X₀=0 must be excluded as starting value. Now the above point set can be written as $\begin{array}{cc}{L}_t=\left\{\frac{1}{m}\left(x,s\left(x\right),s\left(s\left(x\right)\right),\dots ,s^{t-1}\left(x\right)\right)|0\le x<m\right\},\mathrm{where}& \left(2\text{.}17\right)\\ s\left(x\right)=\left(ax+c\right)\mathrm{mod}m& \left(2\text{.}18\right)\end{array}$
It will be seen that the pattern generated by the point set L_t has lattice structure. The set for t=2 and c=0 yields the two-dimensional Korobov lattice. In general L_t for c=0 corresponds to the point set which is produced by a rank-1 lattice with generator vector (1, a, a², . . . , a^t−1) mod m. In two dimensions the lattice points can be covered by a small number of parallel lines which can be chosen from many different families. In higher dimensions the families of lines generalize to families of hyperplanes. The goal of the spectral test is to determine the t-dimensional accuracy v_t which is defined as the reciprocal of the maximum distance between hyperplanes, taken over all families of parallel (t−1)-dimensional hyperplanes covering all points. In the following discussion, it will be derived that v_t corresponds to the length of the shortest vector in the lattice which is dual to the lattice of equation 2.17. This means that the spectral test can be applied for the search for max-min-dist lattices.

[0135] The parameter c of a linear congruential random number generator does not have any effect on the spectral test, since it only brings about a shift of all elements of L_t so that their relative distances are kept. This can be seen when resolving the recurrence $\begin{array}{c}{s}^t\left(x\right)=s\left(s(\dots s\left(x\right)\dots )\right)=\left(a^{t}x+\left(1+a+\dots +a^{t-1}\right)c\right)\mathrm{mod}m\\ =a^{t}x\mathrm{mod}m+\left(1+a+\dots +a^{t-1}\right)c\mathrm{mod}m\\ =a^{t}x\mathrm{mod}m+\gamma ,\end{array}$ $\mathrm{with}$ $\gamma =\left(1+a+\dots +a^{t-1}\right)c$
being constant. Therefore we can set c=0 which means that we directly deal with lattices in Korobov form.

[0136] Since it is two-dimensional max-min-dist lattices that are being searched for, it is sufficient to examine the two-dimensional point set $\begin{array}{cc}{L}_2=\left\{\frac{1}{m}\left(x,s\left(x\right)\right)|0\le x<m\right\},c=0& \left(2\text{.}19\right)\end{array}$
and to calculate v₂, i.e., the reciprocal of the maximum distance between lines, taken over all families of parallel straight lines that cover the points of L₂. At first, the mod m-operation which is hidden in equation 2.18 is exchanged with mod 1 by writing $L_2=\left\{\left(\frac{x}{m},\frac{ax}{m}\right)\mathrm{mod}1|0\le x<m\right\}$
Then, the mod 1-operation is disposed of by extending the set L₂ periodically, i.e., by covering the whole ² with copies of the two dimensional hypercube: $\begin{array}{cc}{\tilde{L}}_2=\left\{\left(\frac{x}{m}+k_1,\frac{ax}{m}+k_2\right)|\mathrm{interger}x,k_1,k_2\right\}& \left(2\text{.}21\right)\end{array}$
This is illustrated in the graph 300 shown in FIG. 13, in which the point set $\left\{\left(\frac{i}{5},\frac{ai}{5}\right)\mathrm{mod}1|0\le i<5\right\}$
is extended periodically, and in which the periodically extended lattice is covered with one family of equidistant parallel lines. The set 2.21 can further be simplified by changing (x, k₁, k₂) to (x+k₁m, 0, k₂−ak₁) since $\frac{x+k_{1}m}{m}=\frac{x}{m}+k_1$ $\mathrm{and}$ $\frac{a\left(x+k_{1}m\right)}{m}+k_2-a{k}_1=\frac{ax}{m}+k_2$
thereby setting the variable k₁ to zero. Doing so allows the above point set to be described by $\begin{array}{cc}{\tilde{L}}_2=\left\{y_1{v}_1+y_2{v}_2|\mathrm{interger}{y}_1·y_2\right\},\mathrm{where}{v}_1=\frac{1}{m}\left(1,a\right)\mathrm{and}{v}_2=\left(0,1\right).& \left(2\text{.}22\right)\end{array}$
It will be seen that set 2.22 contains the m points of the original set 2.19 which satisfy
L₂={tilde over (L)}₂∩[0,1)².
It is now desired to find the maximum distance between adjacent parallel lines among those families of lines that cover all the points of {tilde over (L)}₂. Such a family can be defined by a non-zero vector u=(u₁, u₂) which is perpendicular to all contained lines. Then the set of points on a particular line can be described by
{(x₁,x₂)|x·u=x₁u₁+x₂u₂=q},
where q is a unique constant for each line. Since the lines are equidistant with one of them passing through the origin it is possible to adjust the length of u so that now qε defines all the lines in the family. This means that we can calculate the distance between adjacent lines as the minimum distance from the origin to the line for q=1:
min_{real x1,x2}{√{square root over (x₁²+x₂²)}|x₁u₁+x₂u₂1}.  (2.25)
According to Cauchy's inequality,
(x₁u₁+x₂u₂)²≦(x₁²+x₂²)(u₁²+u₂²)
x_j=u_j/(u₁²+u₂²) for j=1, 2 yields the minimum in equation 2.25, since $x_1^2+x_2^2\ge \frac{{\left(x_1{u}_1+x_2{u}_2\right)}^2}{u_1^2+u_2^2}\ge \frac{1^2}{u_1^2+u_2^2}.$
Hence, the distance between neighboring lines is $\frac{1}{\sqrt{u_1^2+u_2^2}}=\frac{1}{\mathrm{length}\left(u\right)}{=}_{\mathrm{def}}\frac{1}{v_2}.$
The distance will achieve its maximum if
√{square root over (u₁²+u₂²)}
is minimal. For this reason the two-dimensional accuracy v₂ corresponds to the length of the shortest vector u which defines a family of hyperplanes {x·u|integer q} covering all elements of the point set 2.22. Such a vector u=(u₁, u₂) must satisfy the following conditions:

[0137] 1. u≠0,

[0138] 2. v·u=integer ∀vε{tilde over (L)}₂,

[0139] 3. in particular u₁ and u₂ must be integers, because (0,1)ε{tilde over (L)}₂ and (1,0)ε{tilde over (L)}₂ $4.v_1\in {\tilde{L}}_2⟹v_1·u=\frac{1}{m}\left(u_1+a{u}_2\right)=\mathrm{integer},i.e.u_1+a{u}_2\equiv 0\left(\mathrm{mod}m\right).$
Conversely, given any non-zero integer vector u=(u₁,u₂) which fulfills condition 4, the corresponding family of lines will cover all elements of {tilde over (L)}₂, since the dot product (y₁v₁+y₂v₂)·u will result in an integer for all integers y₁,y₂. The spectral test thus needs to find the minimal value $\begin{array}{cc}\begin{array}{c}{v}_2^2={\min }_{\left(u_1,u2\right)\ne \left(0,0\right)}\left\{u_1^2=u_2^2|u_1+a{u}_2\equiv 0\right\}\\ ={\min }_{\left(x_1,x2\right)\ne \left(0,0\right)}\left\{{\left(m{x}_1-a{x}_2\right)}^2+x_2^2\right\},\end{array}& \left(2\text{.}28\right)\end{array}$
where u₁=mx₁−ax₂ and u₂=x₂ is a parametric and general solution of the basic congruence, since mx₁−ax₂+ax₂≡0 (mod m).

[0140] In order to find the minimal value v₂=min_{(x1,x2)≠(0,0)}{(mx₁−ax₂)²+x₂²} this problem is formulated as the more general problem of minimizing the quadratic form
ƒ(x₁,x₂)=(w₁₁x₁+w₂₁x₂)²+(w₁₂x₁+w₂₂x₂)²=(x₁w₁+x₂w₂)·(x₁w₁+x₂w₂)  (2.29)
over all non zero integer vectors (x₁, x₂), where $W=\left(\begin{array}{c}{w}_1^{\top }\\ w_2^{\top }\end{array}\right),w_1=\left(\begin{array}{c}{w}_{11}\\ w_{12}\end{array}\right)\mathrm{and}{w}_2=\left(\begin{array}{c}{w}_{21}\\ w_{22}\end{array}\right)$
is any non-singular matrix. w₁₁=m, w₁₂=0, w₂₁=−a and w₁₂=1 yield the special form of equation 2.28.

[0141] Section 2.1.5 of the Appendix sets forth a code listing of a technique for determining the value of the two-dimensional accuracy v₂. It incorporates a “Euclidean step,” which represents a simplification of the two-dimensional spectral test with respect to higher dimensions and is based on Gauss's notion of a reduced binary quadratic form. It can be proved that all integer solutions to the congruence u₁+au₂≡0 have the form u₁=x₁w₁₁+x₂w₂₁ and u₂=x₁w₁₂+x₂w₂₂ for integers x₁, x₂, given the matrix W such that w₁₁+aw₁₂=w₂₁+aw₂₂≡0 and det(W)=m. If additionally W satisfies the condition
2|w₁₁w₂₁+w₁₂w₂₂|≦w₁₁²+w₁₂²≦w₂₁²+w₂₂²,
(u₁, u₂)=(w₁₁, w₁₂) minimizes u₂₁+u₂₂ over all non-zero integer solutions.

[0142] The technique in the § 2.1.5 code listing takes advantage of these considerations by reducing the quadratic form shown in equation 2.29 until (w₁₁,w₁₂) represents the minimal solution of u₂₁+u₂₂ as stated above.

[0143] As v₂ is the length of the shortest vector which defines a family of lines covering all elements of {tilde over (L)}₂ and due to condition 2, the two-dimensional accuracy can also be regarded as the length of the shortest vector in
L^⊥={u|u·vε∀vε{tilde over (L)}},
i.e., in the lattice which is dual to {tilde over (L)}:
v₂=min_0≠uεL⊥∥u∥₂
Since every point u in L^⊥ defines a set of parallel and equidistant lines such that every point in {tilde over (L)} lies on one of the lines, the shortest vector uεL^⊥ defines the family of lines with maximum distance. Every other choice of u brings the lines closer together. There are now selected two neighboring lines of a family and one lattice point on each of them such that for one fixed point p₁ on the first line the other point p₂ on the second line is closest to p₁. Then the distance between p₁ and p₂ decreases as the distance between two adjacent lines increases. Therefore the shortest vector in L^⊥ also indicates the shortest vector in {tilde over (L)}₂. In fact, v₂, i.e., the length of the shortest vector in L^⊥, equals the length of the shortest vector in {tilde over (L)} up to some scaling factor. This can be proved by considering the basis B of {tilde over (L)} and the basis B^⊥ of L^⊥. Let $B=\frac{1}{m}\left(\begin{array}{cc}a& c\\ b& d\end{array}\right),\mathrm{with}{b}_1=\frac{1}{m}\left(\begin{array}{c}a\\ c\end{array}\right)\mathrm{and}{b}_2=\frac{1}{m}\left(\begin{array}{c}b\\ d\end{array}\right)$
be the Minkowsky reduced basis matrix of {tilde over (L)}. Then the column vectors of (B⁻¹)^⊥ represent the basis of the dual lattice L^⊥. With $\frac{ad}{m^2}-\frac{cd}{m^2}=\frac{1}{m}=\det \left(B\right)=\mathrm{volume}\left(A\right)$
being the determinant of B, the inverse of the basis matrix of {tilde over (L)} can be calculated by $B^{-1}=\frac{1}{m}\left(\begin{array}{cc}d& -c\\ -b& a\end{array}\right){\left(\frac{1}{m}\right)}^{-1}=\left(\begin{array}{cc}d& -c\\ -b& a\end{array}\right)$
Therefore the basis matrix of L^⊥ is $B^{\perp }=\left(\begin{array}{cc}d& -b\\ -c& a\end{array}\right),\mathrm{with}{b}_1^{\perp }=\left(\begin{array}{c}d\\ -c\end{array}\right)\mathrm{and}{b}_2^{\perp }=\left(\begin{array}{c}-b\\ a\end{array}\right)$
Computing the length of the basis vectors of the dual lattice, we notice that they differ from b₁ and b₂ just by the factor 1/m: $b_1^{\perp }·\frac{1}{m}=\sqrt{b^2+d^2}=b_2,b_2^{\perp }·\frac{1}{m}=\sqrt{c^2+a^2}=b_1.$
The basis vectors of the dual lattice result from both rotating b₁·n and b₂·n by ±90° $b_2^{\perp }=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)·\left(\begin{array}{c}a\\ c\end{array}\right)=\left(\begin{array}{c}-c\\ a\end{array}\right)\mathrm{and}{b}_1^{\perp }=\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)·\left(\begin{array}{c}b\\ d\end{array}\right)=\left(\begin{array}{c}d\\ -b\end{array}\right)$
which means that the relative distances between the lattice points in the dual lattice correspond to those in {tilde over (L)}. Since b₁ is the shortest vector in {tilde over (L)}, so is b₂^⊥ in {tilde over (L)}^⊥. Conversely, if we find the shortest vector in the dual lattice L^⊥, its length yields the length of the shortest vector in {tilde over (L)}multiplied by m. This fact is exploited in the search for max-min-dist lattices on the unit square, as set forth in the § 2.1.1 code listing. Instead of calculating the minimum distance by looping over all lattice points, the two-dimensional spectral test is applied. As [0, m)² is computed, the two-dimensional accuracy delivers the sought quantity. The rest of the basic search technique remains unchanged, as shown in the code listing set forth in Section 2.1.6 of the Appendix.

[0144] As the spectral test only needs a few cycles to determine v₂, the technique set forth in the § 2.1.6 code listing is much faster than the basic search technique set forth in § 2.1.1, even when the § 2.1.1 technique incorporates the optimizations discussed above. However, by using the § 2.1.5 technique, the § 2.1.6 technique is restricted to the unit square.

[0145] There is now described a technique for searching the basis for a lattice. So far, the description of the invention has focused on techniques for lattice search. As described in further detail below, aspects of the invention include rasterization on rank-1 lattices and adaptive sampling. Thus, it is necessary to select a good basis for a given lattice. The process of determining a lattice basis according to a given criterion is called “reduction.” For the present purposes, it is advantageous if the basis vectors are as short as possible, resulting in reasonably hexagonal lattice cells. However, as a general matter, it is possible to simply choose the two shortest vectors because they do not necessarily form a basis.

[0146] This problem has been formalized in the following way. A basis G=(g₁, . . . , g_s)^⊥ with g₁, . . . , g_s being row vectors is said to be “Hermite reduced” if L(g₁, . . . , g_s)=L(g₁′, . . . , g_s′) and ∥g_j′∥=∥g_j∥ for jε{1, . . . , i−1}∥g_i∥≦∥g_i′∥ holds true for all iε{1, . . . , s} and for all G′=(g₁′, . . . , g₂′)^⊥. This means that a generator matrix G is reduced if the sequence ∥g₁∥, . . . , ∥g_s∥ is the first element in a list where the corresponding sequences of all generator matrices of the same lattice are lexicographically ordered. In such a reduced basis the first basis vector conforms to the shortest non-zero lattice vector. Although a Hermite reduced basis exists for every lattice, no explicit construction is known.

[0147] A less strict reduction criterion is used for the presently described basis search. A basis G=(g₁, . . . , g_s)^⊥ is reduced according to this criterion if

[0148] 1. g₁ is the shortest non-zero vector in L and

[0149] 2. for 2≦i≦s g_i is the shortest vector in L such that (g₁, . . . , g_i) may be extended to a basis of L.

This so-called “Minkowsky reduced basis” always contains the shortest non-zero lattice vector.

[0150] There are still other reduction techniques that may be used, including the “Korkine-Zolotareff reduction” and the “LLL-reduction.” The choice of reduction technique depends on the underlying problem. Although the problem of finding a Minkowsky reduced basis is NP-hard, for the present purposes it is sufficient to use this reduction, since it is restricted to the two dimensional case and can exploit the special structure of rank-1 lattices.

[0151] To find a Minkowsky reduced basis in two dimensions, a search is first conducted for the shortest non-zero lattice vector g₁. Then, the subsequent basis vector g₂ is determined by searching the second shortest lattice vector in L which is linearly independent of g₁. This can be done by modifying the basic search technique. For a given lattice L(n, a) all lattice points x_i are looped over in order to calculate mindist(x_i) like in equation 2.5, discussed above. The length of both the shortest and the second shortest vector are checked for, updating the basis vectors if necessary. But before updating the second basis vector, it is necessary to make sure that it is linearly independent to the first basis vector so that the conditions of a Minkowsky reduced basis hold. Special care has to be taken for this test because the simple computation of the determinant of the two two-dimensional basis vectors might produce an integer overflow. By definition two vectors x=(x₁, x₂) and y=(y₁, y₂) are linearly dependent if there exists a scalar λ≠0 such that
x=λ·y.
In other words, that means that they are linearly dependent if one vector is a multiple of the other one. This is the case if for $\tilde{x}=\left(\frac{x_1}{\gcd \left(x_1,y_1\right)},\frac{x_2}{\gcd \left(x_2,y_2\right)}\right)\mathrm{and}$ $\tilde{y}=\left(\frac{y_1}{\gcd \left(x_1,y_1\right)},\frac{y_2}{\gcd \left(x_2,y_2\right)}\right).{\tilde{x}}_1={\tilde{x}}_2\mathrm{and}{\tilde{y}}_1={\tilde{y}}_2$
the two resulting vectors {tilde over (x)} and {tilde over (y)} are equal and if in addition sgn(x₁)=sgn(y₁) and sgn(x₂)=sgn(y₂) or sgn(x₁)≠sgn(y₁) and sgn(x₂)≠sgn(y₂).

[0152] Section 2.2.1 of the Appendix sets forth a code listing that checks for the linear independence of two basis vectors. If the signs of the vector components were not involved, the § 2.2.1 code listing would take the vectors (1, 1) and (1,−1) for linearly dependent, for example. FIG. 14 shows a graph 370 illustrating the lattice L(25, 7) and the two basis vectors that result from the basis search algorithm. As shown in FIG. 15, b₁=(3, −4) and b₂=(4, 3).

[0153] Section 2.2.2 of the Appendix shows a code listing in which a search is conducted for a Minkowsky reduced basis. It should be remembered that computations are always performed on [0, n)² and therefore also get the basis vectors on [0, n)². Concerning shifted lattices, the basis search algorithm can be applied in the same way because a shift just represents a translation of the whole lattice and does not influence the basis.

[0154] There is now discussed the use of (0, m, 2)-nets in conjunction with rank-1 lattices. The present description of the invention has already set forth one motivation for shifted lattices with regard to sampling a whole screen in section. A closer look at the lattice in the FIG. 10B graph 330. FIG. 15 shows a graph 380 illustrating a refinement of the grid structure of FIG. 10B and shows the lattice L(25, 7)⊂[0, 1)² shifted by (10/25, 10/25). The lattice points perfectly fill the fine grid structure in a low dimensional projection. This is indicated by the hollow lattice point which is projected onto the x- and y-axis with no other lattice point lying on these lines. Furthermore, each interval $\begin{array}{cc}\left[\frac{i}{5},\frac{i+1}{5}\right)\times \left[\frac{j}{5},\frac{j+1}{5}\right),& 0\le i,j<5\end{array}$
also contains exactly one lattice point. Hence the lattice has the characteristics of both stratified and Latin hypercube sampling points.

[0155] This can be formalized by (0, m, 2)- or more generally (t, m, s). An elementary interval is defined by $E:=\prod _{j=1}^s\left[\frac{a_j}{b^{n_j}},\frac{a_j+1}{b^{n_j}}\right),$
where s≧1 is the dimension, b≧2 represents the base and n_j≧0 and 0≦a_j<b^nj are integers. The volume of E is given by ${\lambda }_s\left(E\right)=\frac{1}{b^{\sum _{j=1}^s{n}_j}}.$
For two integers 0≦t≦m, a finite point set of b^m points in s dimensions is called a (t, m, s)-net in base b, if every elementary interval of volume
λ_s(E)=b^t−m
contains exactly b^t points. As smaller values of t assure a better uniform distribution, t is said to be the quality parameter.

[0156] Keeping these definitions in mind, the lattice shown in the FIG. 15 graph 380 represents a (0, 2, 2)-net in base 5, since every elementary of volume $b^{t-m}=5^{0-2}=\frac{1}{25}$
contains exactly one lattice point. This observation led to a search of all possible shifts for the lattices
L(n,a),n=2^m,mε{4,6,8,10}
and test them for (0, m, 2)-nets in base 2, i.e., check if every elementary interval of volume $2^{-m}=\frac{1}{n}$
contains exactly one point. It will be seen that the examined lattices need to fulfill the condition gcd(n, a)=1, since otherwise it would not be possible to obtain perfect one dimensional projections. But none of the above lattices turned out to be a (0, m, 2)-net in base 2, which means that (0, m, 2)-nets cannot simply be found by shifting rank-1 lattices. Nevertheless (0, m, 2)-nets for m=2 are readily obtained if the base b is prime and gcd(b², a)=1 for the shifted lattices L(b², a), as it is the case for the lattice shown in the FIG. 15 graph 380. This is because then there are only three different elementary intervals of volume l/b² $\begin{array}{ccc}1.& \left[\frac{i}{b},\frac{i+1}{b}\right)\times \left[\frac{j}{b},\frac{j+1}{b}\right),& 0\le i,j<b\\ 2.& \left[\frac{i}{b^2},\frac{i+1}{b^2}\right)\times \left[0,1\right),& 0\le i<b\\ 3.& \left[0,1\right)\times \left[\frac{i}{b^2},\frac{i+1}{b^2}\right),& 0\le i<b\end{array}$
which naturally contain just one lattice point. FIG. 16 shows a table 390 listing some further (0, 2, 2)-nets in base b.

[0157] A raster display is typically formed from a matrix of picture elements (pixels) representing the whole screen. In order to render an image correctly, a rasterizer approximates the mathematical ideal primitives, such as lines, triangles or circles described in terms of vertices, on a Cartesian grid by assigning the appropriate colors to sets of pixels. The rasterizer converts the two-dimensional vertices in screen space into pixels on the screen. In current approaches, a pixel is typically represented as a square on the integral raster [0, . . . , ResX)×[0, . . . , ResY).

[0158] An aspect of the present is directed to a rasterization technique based on rank-1 lattices. According to this approach, rasterization is performed on the lattice cells of a rank-1 lattice instead of on a rectangular grid; i.e., the pixel layout is structured by rank-1 lattices. FIG. 17 is a flowchart of a method 400 according to this aspect of the invention.

[0159] In step 401, a rank-1 lattice is selected in accordance with a maximized minimum distance function corresponding to points in a synthetic image to be generated. In step 402, a point p in Cartesian coordinates is transformed into the basis of the rank-1 lattice. In step 403, it is determined the lattice cell into which the transformed point falls. In step 404, the Cartesian coordinates are converted into rank-1 lattice coordinates. In step 405, the rank-1 lattice coordinates are used to rasterize the image.

[0160] In the following discussion, the rasterization technique is described in conjunction with a Digital Differential Analyzer (DDA), which is a scan conversion technique for lines, and in conjunction with Crow's algorithm, which is a rasterization technique for convex polygons. The described technique can be performed using traditional raster displays. In addition, as discussed below, the technique may also form the basis of a new concept of display technology.

[0161] A scan-conversion technique for lines is based on computing the coordinates of those pixels that best approximate an ideal, straight line drawn on the two-dimensional raster. Typically, a line is characterized by a start point, an end point, a width, and by and additional style attributes. The present discussion focuses on lines having a width of one pixel, without additional attributes and without anti-aliasing.

[0162] The simplest technique for scan-converting lines is a Digital Differential Analyzer (DDA), which is based on the Taylor series expansion of a function $\begin{array}{cc}P\left(t+\tau \right)=P\left(t\right)+\tau \frac{dP\left(t\right)}{dt}+O\left({\tau }^2\right).& \left(3.1\right)\end{array}$

[0163] Since rasterization is performed on an integer grid, the endpoints of the lines are assumed to be integer coordinates. For a given pair of (integer) points P₀=(x₀, y₀) and P₁=(x₁, y₁), a straight line can explicitly described by $P\left(t\right)=\left(\begin{array}{c}{x}_0+t\left(x_1-x_0\right)\\ y_0+t\left(y_1-y_0\right)\end{array}\right),\mathrm{for}t\in \left[0,1\right].$
Using equation 3.1, a point (t+τ, P(t+τ)) on the line can incrementally be calculated from its preceding point (t, P(t)): $P\left(t+\tau \right)=P\left(t\right)+\tau \frac{dP\left(t\right)}{dt}=P\left(t\right)+\tau \left(\begin{array}{c}{x}_1-x_0\\ y_1-y_0\end{array}\right)=P\left(t\right)+\tau \left(\begin{array}{c}\Delta x\\ \Delta y\end{array}\right)$
In order to draw the line incrementally from its start to end point, the rasterization process has to be discretized inn steps, where n=max{|x₁−x₀|, |y₁−y₀|}. This means that the pixel coordinates approximating the ideal line are incrementally computed by $x_i+1=x_i\frac{\Delta x}{n}\mathrm{and}$ $y_i+1=y_i\frac{\Delta y}{n},\mathrm{for}0\le i<n.$
If for example |Δx|>|Δy|, n=|x₁−x₀|. Hence x_i is incremented by ±1 in each step, whereas y_i is incremented by a fraction $0\le \frac{\Delta y}{n}\le 1.$
After calculating the coordinates, the pixel (└x_i┘, └y_i┘), which is closest to the ideal line, is colored. For visualization, a pixel is interpreted as a sample in the center of the cell formed by two adjacent grid lines. FIG. 18A shows a graph 410, in which the pixel (2, 2) is intensified on the raster grid. FIG. 18B shows a graph 420 illustrating the rasterization of the line with integer endpoints (0, 0) and (3, 1). The rasterization process consecutively intensifies the pixels $\left(0,0\right),\left(1,\lfloor \frac{1}{3}\rfloor \right),\left(2,\lfloor \frac{2}{3}\rfloor \right),$
and (3,1), since n=3, dx=3 and dy=1. Section 3.1.1 of the Appendix sets forth a source code listing for the DDA.

[0164] In the above discussion, the DDA has been described with respect to traditional raster displays. According to an aspect of the invention, the rasterization technique is converted to rank-1 lattices by changing the basis in which the rasterization is performed. This change means that the rasterizer switches from a Cartesian coordinate system to a coordinate system formed by the basis of a corresponding rank-1 lattice. Hence, it will be seen that, according to the present aspect of the invention, the rasterization process does not operate on pixels any more, but rather on lattice cells.

[0165] For visualization, each lattice cell is interpreted as a sample in the center of the fundamental parallelepiped spanned by the Minkowsky reduced lattice basis. Each cell (x, y) is addressed in lattice coordinates, i.e., in the coordinate system formed by the lattice basis. The use of lattice coordinates is illustrated in FIGS. 19A and 19B. FIG. 19A shows a graph 430 illustrating the intensification of lattice cell (0, 2) on the lattice grid. FIG. 19B shows a graph 440 illustrating rasterization of the line extending from (1, 3) to (7, 14) on the max-min-dist lattice L(576, 155). The intensified (shaded) lattice cells in FIG. 19B represent the scan-converted line.

[0166] In order to specify a line as simple as possible and independent of the lattice on which the rasterization is performed, the scan-conversion technique assumes that its endpoints are given in Cartesian coordinates on the unit square. The lattice coordinates are obtained by transforming these points into the lattice basis and by determining the lattice cells into which the transformed points fall. The coordinates of these two lattice cells then define the start and endpoint of the line. Section 3.1.2 of the Appendix sets forth a source code listing for DDA rasterization on rank-1 lattices. FIG. 20 shows a table 450 setting forth the single steps performed by the § 3.1.2 code listing for the rasterization of the line in FIG. 19B from (1, 3) to (7, 14). The last line of the FIG. 20 table 450 contains the lattice cells approximating the ideal line.

[0167] In FIG. 19B, the ideal line is originally defined by the points $P_0=\left(\frac{2.5}{24},\frac{1.5}{24}\right)$
and $P_1\left(\frac{12.5}{24},\frac{5.5}{24}\right).$
The transformation in lattice coordinates yields {tilde over (P)}₀=(1.25, 3.0833) and {tilde over (P)}₁=(7.4999, 14.5). As {tilde over (P)}₀ lies within the lattice cell (1, 3) and {tilde over (P)}₁ within (7, 14), the line is identified by the start point (1, 3) and the end point (7, 14).

[0168] Given a point (x, y) in Cartesian coordinates and the lattice basis (b1 b2) with $b1=\left(\begin{array}{c}{b}_{11}\\ b_{12}\end{array}\right)\mathrm{and}b2=\left(\begin{array}{c}{b}_{21}\\ b_{22}\end{array}\right),\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{cc}{b}_{11}& b_{21}\\ b_{12}& b_{22}\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right),$
where (u, v) represents the coordinates of the point (x, y) in the lattice basis. As the basis matrix (b1b2) is not singular, the coordinates (u, v) are the result of $\left(\begin{array}{c}u\\ v\end{array}\right)={\left(\begin{array}{cc}{b}_{11}& b_{21}\\ b_{12}& b_{22}\end{array}\right)}^{-1}\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{c}{\mathrm{xb}}_{22}-{\mathrm{yb}}_{21}\\ {\mathrm{yb}}_{11}-{\mathrm{xb}}_{12}\end{array}\right)·\frac{1}{b_{11}{b}_{22}-b_{21}{b}_{b12}}.$
This is implemented in the code listing set forth in section 3.1.3 of the Appendix, which additionally handles the transformation into the coordinate system of a shifted lattice. As a lattice shift only translates the whole lattice without affecting the basis, the shift vector (sx, sy)ε[0, n)² is also transformed into the lattice basis. Then the transformed shift is applied to (u, v). The lattice cell (c_x, c_y) which a point (u, v) in lattice coordinates belongs to is computed by rounding down its components, i.e,
c_x=└u┘ and c_y=└v┘.
This technique is as a method called determineLatticeCell that is implemented in the code listing set forth in section 3.1.4 of the Appendix. Once both the start and end point of the line are converted into lattice coordinates, the rasterization is realized the same way as in the § 3.1.1 code listing, discussed above. The only difference consists in the computation of the next point (x_i+1, y_i+1) on the line. Instead of $x_{i+1}=\mathrm{xi}+\frac{\Delta x}{n}$
and $y_{i+1}=y_i+\frac{\Delta y}{n},$
there are calculated $\begin{array}{c}{x}_{i+1}=\frac{x_{i·n}+\Delta x}{n}\mathrm{and}\\ y_{i+1}=\frac{y_{i·n}+\Delta y}{n}.\end{array}$
As the denominator never changes, only the nominator has to be accumulated in every step. The division is only performed to determine the lattice cell which has to be intensified. That way the accumulation of rounding errors can be avoided. Assuming that |Δy|>|Δx|, it follows that n=|Δy|. $\begin{array}{c}{y}_{i+1}=\frac{y_i·n+\Delta y}{n}\\ =\frac{y_i·n+\mathrm{sgn}\left(\Delta y\right)·n}{n}\\ =y_i+\mathrm{sgn}\left(\Delta y\right).\mathrm{where}\mathrm{sgn}\left(x\right)=\{\begin{array}{cc}1& \mathrm{if}x\ge 0\\ -1& \mathrm{else}\end{array}\end{array}$
In that case y_i just has to be incremented by sgn(Δy) in each step, wherefore y_nom is initialized by y₀ and y_denom by 1 in the above algorithm. The same applies if |Δx|≧|Δy|. Applying the § 3.1.2 code to the line shown in the FIG. 19B graph 390, yields:
Δx=7−1=6;Δy=14−3=11; n=max {6,11}=11, x_nom=1·11=11, x_denom=11, y_nom=3 and y_denom=1.

[0169] Due to the representation of a lattice cell as the center of the fundamental parallelepiped 0.5 is added to both x_i+1 and y_i+1 in the § 3.1.2 code, i.e., $\begin{array}{c}x=x_{i+1}\\ =\frac{x_{\mathrm{nom}}}{x_{\mathrm{denom}}}+0.5\mathrm{and}\\ y=\mathrm{yi}+1\\ =\frac{y_{\mathrm{nom}}}{y_{\mathrm{denom}}}+0.5.\end{array}$
The coordinates of the lattice cells which approximate best the ideal line are identified analog to the code set forth in § 3.1.1, namely, by (└x┘,└y┘).

[0170] The lattice-based rasterization technique is now described with respect to Crow's algorithm. A polygon is rasterized by finding the pixels which fill the polygon given by a set of vertices and edges. Crow's algorithm assumes that the vertices are ordered in a consistent way and that the polygons to be converted are convex. However, because every polygon can be decomposed into convex polygons, this requirement does not present a limitation.

[0171] Crow's algorithm takes advantage of the coherence in both the horizontal and the vertical scanline. Looping over a horizontal scanline from bottom to top, the algorithm determines in each scanline which pixels are inside the polygon by means of an incremental Jordan test for horizontal coherence. The Jordan test is based on the Jordan curve theorem, which states that any ray from inside a polygon will intersect the edges of a polygon an odd number of times. As only convex polygons are considered, there is just one span per scanline which is bordered by points V scanLeft and V scanRight on the current left and right edges. These points are found by exploiting the vertical coherence by means of a DDA approach.

[0172] This technique starts by finding the vertex bot with the smallest y-value. As long as there are vertices not having been considered yet the algorithm proceeds by determining an appropriate left edge V left→V nextLeft and an appropriate right edge V right→V nextRight, where V left=V right=bot at the beginning. While scanning upward, it is checked whether these edges span the current scanline and the appropriate span is drawn. Crow's algorithm for a standard pixel display is set forth in the code listing shown in section 3.2.1 of the Appendix, which assumes the vertices to be ordered counterclockwise. FIG. 21 shows a table 420 setting forth a description of the variables used in the § 3.2.1 code listing.

[0173] In order to rasterize convex polygons on rank-1 lattices the same approach is used as above. This means that the algorithm is performed in the coordinate system spanned by the basis of the particular rank-1 lattice L(n, a) instead of the Cartesian coordinate system. Thereby the first basis vector b₁ of the Minkowsky reduced basis of L(n, a) corresponds to the horizontal and the second basis vector b₂ to the vertical coordinate axis. The rasterization algorithm takes as arguments an array of vertices and the size of this array.

[0174] According to the code listing set forth in § 3.1.2, the vertices which have to be ordered counterclockwise are assumed to be given in Cartesian coordinates on the unit square. After transforming these vertices into the lattice basis, the rasterization can be implemented the same way as in the code listing of § 3.2.1 up to one exception, i.e., the basic rasterization algorithm remains the same, as shown in the code listing in § 3.2.4. After finding the vertex bot with the smallest y-value, the algorithm scans upwards by maintaining the active edges V left→V nextLeft and V right→V nextRight. Due to the assumption that the vertices are ordered counterclockwise, the indices of V nextLeft and V nextRight are calculated by li=(li+n−1) mod n and ri=(ri+1) mod n. But as the vertices are originally specified on the unit square of the Cartesian coordinate system, this order may change in the lattice basis. Although the vertices will still be ordered in a consistent way after the coordinate system transformation, it might be necessary to exchange V scanLeft and V scanRight in order to assure that lx≦rx when drawing the current polygon span. This means that
lx=ceil(l.m_—X);
rx=ceil(r.m_—X);
has to be changed to
lx=l.m_—X<r.m_—X?ceil(l.m_—X):ceil(r.m_—X);
rx=l.m_—X<r.m_—X?ceil(r.m_—X):ceil(l.m_—X);
in the code listing set forth in section 3.2.2 of the Appendix. Furthermore, it is necessary to call drawLatticeCell (x,y) instead of drawPixel (x,y).

[0175] FIG. 22 shows a graph 470 of a polygon defined by the points $\begin{array}{c}{P}_0=\left(\frac{10}{48},\frac{4}{48}\right),\\ P_1=\left(\frac{40}{48},\frac{8}{48}\right),\\ P_2=\left(\frac{20}{48},\frac{34}{48}\right)\mathrm{and}\\ P_3=\left(\frac{1}{48},\frac{24}{48}\right)\end{array}$
on the unit square. After transformation into the basis of the lattice L(576, 155) the polygon is specified by {tilde over (P)}₀=(3.125, 5.708), {tilde over (P)}₁=(15, 21), {tilde over (P)}₂=(−1.875, 17.375) and {tilde over (P)}₃=(−7.062, 5.979). FIGS. 23A and 23B show graphs 480 and 490 displaying the same polygon rendered in higher resolutions.

[0176] The principle for the rasterization on rank-1 lattices, which has been described in this chapter considering as example the DDA and the algorithm of Crow, can also be applied for other rasterization algorithms, like for example for the Bresenham algorithm for lines and the y-x rasterization algorithm for polygons. One disadvantage of this rasterization approach is that edges being almost vertical or rather horizontal appear very jagged due to the shape of the lattice cells. One possible solution to that problem consists in anti-aliasing, as it is done for traditional raster displays. Taking a closer look at the fundamental parallelepiped reveals another possibility. Using Delaunay triangulation, by splitting the lattice cells along their shorter diagonals, i.e., into two triangles, the overall resolution can be doubled which results in considerably smoother edges. This is shown in FIG. 24, which shows a graph 500, where this approach has been applied to the polygon shown in FIG. 23A.

[0177] As mentioned above, the described techniques provide the basis for a new display technology. The first computer monitors created images by moving an electron beam directly along the paths between specified points on the display. Thus these vector displays were able to generate sharp and smooth lines. However they only drew wire frames and could only display scenes of limited complexity. Even commonly used for scientific instruments like oscilloscopes, they have been replaced by raster displays. Such a display usually is constituted by a rectangular grid of pixels which are represented as squares on this raster. In contrast to vector displays images are created by scan line algorithms. Moreover the pixels can be addressed individually and they can be colored by means of individual red, green and blue phosphors, in the case of a CRT display.

[0178] Instead of rasterizing on a rectangular grid, the rasterization can be performed on lattice points, as discussed above. Whereas this method can be simulated on traditional raster displays by means of a software solution, it can even be performed on current graphics hardware in the following way. Since the rasterizer is only capable of operating on rectangular grids, in a first step the scene has to be transformed into the lattice basis, which in fact corresponds to a shear of the rectangular grid. After this change of frame the rasterization can be performed on the graphics hardware as usual. In order to display the rasterized scene, the resulting image has to be transformed back into the pixel basis. This procedure yields the same results as if the rasterization would have directly been performed on a rank-1 lattice.

[0179] On the basis of these considerations we introduce a new concept of display technology. Instead of a matrix of square pixels the raster display is conceived as a matrix of lattice cells, i.e., the single picture elements are represented as the cells which are induced by the Delaunay triangulation of the rank-1 lattice points, i.e., the Voronoi diagram. This means that the raster display is conceptualized as a rank-1 lattice which is selected by maximized minimum distance. Thus the picture elements approximate the hexagonal grid for an increasing number of pixels. By structuring the pixel layout by rank-1 lattices, the display takes advantage of the filter properties of rank-1 lattices, discussed below, featuring smaller wave vectors as the traditional rectangular grids. Hence functions whose Fourier coefficients decrease with increasing frequency can be approximated in an especially good way. Since rank-1 lattices are available for any number n of points or rather lattice cells, the number of pixels can be chosen freely. FIG. 25 illustrates a layout 510 of a display module consisting of 256 lattice cells.

[0180] This concept can technically be implemented in a number of ways. A first implementation comprises a display of point light sources, for example RGB LEDs, which are arranged in the center of each lattice cell, indicated by the small black dots in FIG. 25. Since the rank-1 lattice is chosen by maximum minimum distance, the RGB-LEDs are distributed over the screen in an optimal way. A second implementation composes the display of area light sources covering the single picture elements. This yields a technique for thin-film transistor (TFT) displays and liquid crystal displays (LCDs), respectively, which may be realized by means of organic light-emitting diodes (OLEDs) for instance.

[0181] As mentioned above, the overall resolution of a rank-1 lattice display can be doubled by splitting the lattice cells along their shorter diagonal, thus reducing aliasing artifacts, as shown in the FIG. 24 graph 500. This principle can directly be realized in hardware by adapting the layout of the area light sources. This is shown in FIGS. 26A and 26B for a monochrome display 520 and a color display 530.

[0182] So far the pixel layout has been determined by the cells induced by the Delaunay triangulation of the rank-1 lattice points. By means of the Voronoi diagram, which is dual to the Delaunay triangulation in ², the ² can be partitioned into almost hexagonal cells centered around each lattice point. FIG. 27 shows a Voronoi diagram 540 of the maximum minimum distance rank-1 lattice L(64, 19). Structuring the pixel layout in this sense leads to substantially hexagonal pixel lattices. As a principle of nature, ommatidia in compound eyes or the visual sensors in the retina are arranged in this manner, for example. The concept of hexagonal pixels has already been studied in the context of hexagonal image processing and is used, for example, in SmartSlab LED panels.

[0183] One advantage of the described pixel layout is the smoother representation of curved objects. This is due to the fact that all neighbors have a uniform distance from each other, i.e., they are connected in a consistent way. By means of this pixel layout, color displays can simply be realized by adapting the number of pixels to the number of basis colors used by the display. This means that the total number of pixels is computed by
number of pixels·color cells.

[0184] FIG. 28A shows a layout 550 of area light sources for hexagonal pixels for a monochrome display. FIG. 28B shows a layout 560 of area light sources for hexagonal pixels for a color display.

[0185] In a third implementation, the layout of the sensors, i.e., the photosites, of a charge-coupled device (CCD) camera take the form of a rank-1 lattice cell which is either induced by the Delaunay triangulation or the Voronoi diagram.

[0186] In a fourth implementation, rank-1 lattice patterns are used in projector technology. As an example, the hinge-mounted microscopic mirrors of digital light processing (DLP) integrated circuits could be arranged in this manner.

[0187] In a fifth implementation, the described patterns can be applied in the context of 3D displays, such as the Sharp LL-151-3D monitor, and other monitors disclosed at www.sharpsystems.com, in order to realize both the display layout and the parallax barrier, by which the 3D mode is enabled.

[0188] Since rank-1 lattices can be tiled seamlessly, it is possible to compose a display of k modules, each of which having the same number of lattice cells. FIGS. 29A and 29B show displays 570 and 580 according to this aspect of the invention. In a “2ⁿ display,” the number of cells for a display module is selected to be a power of 2. This type of display module has the advantage that the single cells easily can be addressed by means of a demultiplexer. If the number k of modules is set to a power of 2 as well, the single modules can be controlled the same way. FIG. 30 is a diagram of a system 590 according to this aspect of the invention. Such displays can be produced economically by fabricating small modules of the same layout and the same number of lattice cells which can easily be assembled to a display of desired resolution.

[0189] More generally, the concept of 2ⁿ displays perfectly fits all aspects of computer technology, taking advantage of memory layout, cache lines, addressing, and the like. Storing images according to this scheme leads to the concept of “2ⁿ images,” which equally benefit from the advantages of 2ⁿ displays. So the O(2ⁿ) memory requirements ideally fit paging (memory alignment). As a further example, storing a sequence of textures as 2⁰ . . . 2ⁿ images naturally supports MipMapping which thus can be performed in a very efficient way.

[0190] FIG. 31 shows a table 600 that lists the parameters for some 2ⁿ display modules and images, respectively. The listed sizes represent the proposed resolution for the single modules according to which the maximum minimum distance rank-1 lattices have been selected.

[0191] FIG. 32 shows a flowchart of a method 610 according to a further aspect of the invention for anti-aliasing by rank-1 lattices. In step 611, rank-1 lattices are selected in accordance with a maximized minimum distance function corresponding to points in a synthetic image to be generated. In step 612, the selected rank-1 lattices are compared with respect to their convergence to either Monte Carlo or quasi-Monte Carlo point sets. In step 613, an anti-aliasing function is applied, based on the result of this comparison. The following discussion develops the theoretical basis for the method, and describes further aspects of the present invention relating to anti-aliasing.

[0192] FIG. 33 shows a synthesized image 620 of a checkerboard of infinite size that has been used to develop and illustrate techniques for anti-aliasing by rank-1 lattices. Given a “real” checkerboard of infinite size, ideally, one would perceive a smooth and continuous image averaged to half grey at the horizon. However, a synthesized image generated by a computer is a digital signal. Since a computer is only capable of storing discrete values of finite precision, a computer-generated image must perforce also be discrete. Thus, in generating an image, a signal that is typically aperiodic and discontinuous must be converted into a digital signal.

[0193] In order to generate an image of the checkerboard it is necessary to sample the original signal. This transition from continuous to discrete may induce artifacts, e.g., in form of jaggies or Moiré patterns. These effects are known as “aliases.” The general process of displaying a continuous signal a begins by gathering samples of a at specific positions. That is, a sample pattern s is multiplied with a. The resulting signal b consists of a zero-valued function with spikes at the sampling positions. In order to recover the original signal, b is convolved with a reconstruction filter r yielding the signal c=b*r. In order to reproduce the signal on a display, the resulting signal c is resampled with a set of samples reflecting the frame buffer memory locations. Finally, the convolution of r with the display reconstruction filter gives the actually displayed signal. In practice, the sampled signal b is directly reconstructed for each single pixel, thereby omitting the resampling step.

[0194] The following discussion begins with an overview of the signal processing theory relating to aliasing. The overview followed by a description of anti-aliasing by rank-1 lattices, including a comparison of the convergence of test-patterns. A method is then described for adaptive sampling by rank-1 lattices. In addition, the filter properties of rank-1 lattices are described, and a method is presented for depth adaptive filtering in a ray tracing environment. As described below, this method provides a technique for automatically choosing the appropriate reconstruction filter depending on the distance between camera and objects. Finally, anti-aliasing by rank-1 lattices is examined on patterned displays by means of the Clear-Type technology, which incorporates a filtering technique called RGB-decimation.

[0195] The ideal sampling and reconstruction of a one-dimensional continuous signal is now discussed. The two-dimensional case is very similar to the one-dimensional case, from which it differs only in its complexity. By “ideal sampling and reconstruction,” it is meant that the continuous, band-limited signal is uniformly sampled and reconstructed in such a way that the original signal is re-obtained in the end.

[0196] The process of ideal sampling and reconstruction can be summarized by the following equation:
ƒ(x)·s(x)*r(x)=ƒ(x),
where ƒ is the original signal or function, respectively, s is the uniform sampling function and r is the reconstruction filter. As already mentioned, the signal ƒ needs to be band-limited. The reason for this requirement is that the Fourier transform of a uniformly sampled signal turns out to be periodic in frequency space. A function ƒ is said to be band-limited if
F(ω)=0 for |ω|>ω₀,
where F(ω) is the Fourier transform of ƒ(x). This means that the spectrum
F(ω)=(ƒ(x))
of ƒ has finite support. The frequency ω₀ is called the cutoff frequency for ƒ. If ω₀ =∞, the function ƒ is not band-limited.

[0197] In the first step, the function ƒ is uniformly sampled by the Shah-function, as given in the next equation: $s\left(x\right)=\sum _{k=-\infty }^{\infty }\delta \left(x-k\Delta x\right).$
The sampling is done by multiplying ƒ with the Shah-function: $\begin{array}{c}{f}_s\left(x\right)=f\left(x\right)·s\left(x\right)\\ =f\left(x\right)·\sum _{k=-\infty }^{\infty }\delta \left(x-k\Delta x\right).\end{array}$
To calculate the Fourier transform of ƒ_s(x), according to the convolution theorem the Fourier transforms of ƒ(x) and s(x) have to be convolved in frequency space. $\begin{array}{c}ℱ\left(f_s\left(x\right)\right)=F\left(\omega \right)*S\left(\omega \right)\\ =F\left(\omega \right)*\frac{1}{\Delta x}\sum _{k=-\infty }^{\infty }\delta \left(\omega -\frac{k}{\Delta x}\right)\\ =F_s\left(\omega \right)\end{array}$
Calculating the convolution integral we get $\begin{array}{c}F\left(\omega \right)*S\left(\omega \right)={\int }_{-\infty }^{\infty }F\left(\hat{\omega }\right)S\left(\omega -\hat{\omega }\right)d\hat{\omega }\\ =\frac{1}{\Delta x}\sum _{k=-\infty }^{\infty }{\int }_{-\infty }^{\infty }F\left(\hat{\omega }\right)\delta \left(\omega -\hat{\omega }-\frac{k}{\Delta x}\right)d\hat{\omega }\\ =\frac{1}{\Delta x}\sum _{k=-\infty }^{\infty }F\left(\omega +\frac{k}{\Delta x}\right)\end{array}$
So the Fourier transform of ƒ_s(x) is a periodic repetition of the spectrum of ƒ(x) with repetition period 1/Δx. If $\frac{1}{\Delta x}\ge 2{\omega }_0,$
the copies of F(ω) will not overlap. Therefore F(ω)=ΔxF_s(ω) holds true for ωε[−ω₀, ω₀]. But if $\frac{1}{\Delta x}<2{\omega }_0$
the periodic repetition of F(ω) will overlap and accumulate. In order to reconstruct ƒ(x) from its discrete samples, it does not suffice to perform the inverse Fourier transform ⁻¹ due to the periodic repetition of F(ω). ⁻¹(F_s(ω)) would result in a signal being zero everywhere but at λΔx, λε. Therefore the center copy of F(ω) needs to be isolated. Since F(ω)=ΔxF_s(ω) for ωε[−ω₀, ω₀],if $\frac{1}{\Delta x}\ge 2{\omega }_0,$
F(ω) can be recovered by multiplying F_s(ω) with a box spectrum $R\left(\omega \right)=\{\begin{array}{cc}\Delta x,& \mathrm{for}\omega \in \left[-{\omega }_0,{\omega }_0\right]\\ 0,& \mathrm{otherwise}.\end{array}$
In frequency space the reconstruction is performed by $\begin{array}{c}{F}_s\left(\omega \right)·R\left(\omega \right)=\frac{1}{\Delta x}\sum _{k=-\infty }^{\infty }F\left(\omega +\frac{k}{\Delta x}\right)·R\left(\omega \right)\\ =\{\begin{array}{ccc}0,& & \omega >{\omega }_0\\ \frac{1}{\Delta x}{\sum _k^{\infty }}_{\forall \omega +\frac{k}{\Delta x}\le {\omega }_0}& F\left(\omega +\frac{k}{\Delta x}\right)\Delta x,& \mathrm{otherwise}\end{array}\\ =\{\begin{array}{cc}0,& \omega >{\omega }_0\\ F\left(\omega \right),& \mathrm{otherwise}\end{array}=F\left(\omega \right)\end{array}$
As the inverse Fourier transform of the box spectrum results in the sinc-function, we reconstruct ƒ(x) in time space by convolving ƒ_s(x) with the sinc-function ${ℱ}^{-1}\left(F_s\left(\omega \right)·R\left(\omega \right)\right)=f_s\left(x\right)*r\left(x\right)=2{\omega }_0\Delta x\sum _{k=-\infty }^{\infty .}f\left(k\Delta x\right)\mathrm{sinc}2{\mathrm{\pi \omega }}_0\left(x-k\Delta x\right)=f\left(x\right).$
However, this only works if $\frac{1}{\Delta x}\ge 2{\omega }_0.$
When the sampling frequency is too low, i.e. $\frac{1}{\Delta x}<2{\omega }_0,$
the box filtering will cut out F(ω)+F(ω′), for ω′≠0 instead of F(ω). This means that some additional energy, a so called “alias,” is added to F(ω), which makes it impossible to recover F(ω) from its samples. These considerations are summarized in Shannon's sampling theorem:

[0198] Theorem 3 (One dimensional uniform Sampling Theorem) A continuous function ƒ may be perfectly reconstructed from its uniform samples with sampling frequency 1/Δx, if ƒ is bandlimited with cut-off frequency ω₀ and $\frac{1}{\Delta x}\ge 2{\omega }_0.$
Simplified, the theorem states that a band-limited function ƒ may be perfectly reconstructed from its uniform samples if the sampling frequency is at least twice as high as the maximum frequency ω₀ in the spectrum of ω₀. ω₀ is also knows as the “Nyquist frequency.” As used herein, the term “undersampling” refers to sampling below the Nyquist frequency. If a signal is sampled above the Nyquist frequency, it is said to be “oversampled.” The reconstruction filter depends on whether we calculate in time or in frequency space. In time space it is necessary to convolve the correctly sampled signal with the sinc-function, whereas in frequency space, it is necessary to multiply F_s(ω) with a box spectrum.

[0199] Even if the sampling theorem is fulfilled, artifacts may occur in the reconstructed signal. Let for example $R\left(\omega \right)=\{\begin{array}{cc}2\Delta x,& \mathrm{for}\omega \in \left[-{\omega }_0,{\omega }_0\right]\\ 0,& \mathrm{otherwise}\end{array}$
be the reconstruction filter used in frequency space. Then obviously, F_s(ω)·R(ω)≠F(ω), since energy from above ω₀ adds up to F(ω). In such a case, the term used is “reconstruction errors” because of poor reconstruction filters, as opposed to aliasing due to an improper sampling rate.

[0200] However, it is necessary to consider that the sampling theorem deals with continuous functions, which usually is not the case for the functions used in computer graphics. This is demonstrated by the simple example of rendering a checkerboard.

[0201] The checkerboard image 620 displayed in FIG. 33 was uniformly sampled with one sample placed in the center of every pixel and a resolution of 512×512. So the sample pattern represents a regular grid, such as the grid 670 shown in FIG. 35A. In terms of the sampling theorem, this means that we have a sampling frequency of 1/p, where p stands for the space between two pixels. In order to meet the sampling theorem 1/p should be twice as high as the cutoff frequency in the signal. But since the checkerboard provides many sharp edges, the problem arises of having a signal with compact support in signal space which cannot simultaneously have a compact support in frequency space. This means that the signal cannot be band-limited, thereby featuring increasingly high frequencies. Hence there are frequencies >1/p which cannot be represented directly and which turn into aliases, as shown in FIG. 33. In the following discussion, there are described various techniques for anti-aliasing by comparing different sampling and filtering schemes.

[0202] As discussed above, the majority of cases feature $\frac{1}{\Delta x}\le {\omega }_0,$
because typically signals of finite width need to be handled in computer graphics. Generating a synthetic image, it is necessary to somehow determine one single color for each pixel. One approach for reconstructing ƒ(x) as exactly as possible is to compute the area of all the polygon fragments visible in the pixel x_p analytically and to average their colors according to the resulting area proportion. However, this approach has proven to be expensive to calculate or, in some cases, impossible due to a missing analytic expression of the visible area. Hence, an alternative approach is to use a numerical solution. This can be done by choosing a number of sampling points N per pixel and estimating the area of the polygon fragments in x_p from these samples. Since the samples are taken at a higher frequency than once per pixel, the described sampling is referred to herein as “supersampling.” One common supersampling scheme is to use an n×n grid as sampling pattern within each pixel and average the samples to one value by weighting them appropriately. This can by described by the “pixel integral” $\begin{array}{c}\left(f·\tilde{s}*\tilde{r}\right)\left(x_p\right)={\int }_{-\infty }^{\infty }f\left(x\right)\left(\frac{1}{N}\sum _{i=0}^{N-1}\delta \left(x-x_i\right)\right)dx\\ =\frac{1}{N}\sum _{i=0}^{N-1}{\int }_{-\infty }^{\infty }f\left(x\right)\delta \left(x-x_i\right)dx\\ =\frac{1}{N}\sum _{i=0}^{N-1}f\left(x_i\right),\end{array}$
where in this case the sampling function {tilde over (s)} represents the fine uniform grid and the reconstruction filter {tilde over (r)} weights each sample equally with 1/N.