SSymmetry, Integrability and Geometry: Methods and Applications SIGMA (2009), 075, 23 pages Image Sampling with Quasicrystals
Mark GRUNDLAND † , Jiˇr´ı PATERA ‡ , Zuzana MAS ´AKOV ´A § and Neil A. DODGSON †† Computer Laboratory, University of Cambridge, UK
E-mail:
[email protected], [email protected] ‡ Centre de Recherches Math´ematiques, Universit´e de Montr´eal, Canada
E-mail: [email protected] § Department of Mathematics FNSPE, Czech Technical University in Prague, Czech Republic
E-mail: zuzana.masakova@fjfi.cvut.cz
Received December 15, 2008, in final form July 06, 2009; Published online July 20, 2009doi:10.3842/SIGMA.2009.075
Abstract.
We investigate the use of quasicrystals in image sampling. Quasicrystals pro-duce space-filling, non-periodic point sets that are uniformly discrete and relatively dense,thereby ensuring the sample sites are evenly spread out throughout the sampled image.Their self-similar structure can be attractive for creating sampling patterns endowed witha decorative symmetry. We present a brief general overview of the algebraic theory of cut-and-project quasicrystals based on the geometry of the golden ratio. To assess the practicalutility of quasicrystal sampling, we evaluate the visual effects of a variety of non-adaptiveimage sampling strategies on photorealistic image reconstruction and non-photorealistic ima-ge rendering used in multiresolution image representations. For computer visualization ofpoint sets used in image sampling, we introduce a mosaic rendering technique.
Key words: computer graphics; image sampling; image representation; cut-and-project qua-sicrystal; non-periodic tiling; golden ratio; mosaic rendering
Non-periodic tilings have emerged as an important mathematical tool in a variety of computergraphics applications [27]. They have proven especially useful in the design of sampling algo-rithms, where they serve to direct the spatial distribution of rendering primitives by enforcingspatial uniformity while precluding regular repetition. Recently, Wang tilings [3, 19, 23, 24, 25],Penrose tilings [37], Socolar tilings [38] and polyominoes [39] have been used to generate pointsets for non-periodic sampling. In one of the earliest applications of non-periodic tilings in com-puter graphics, Penrose tilings [12, 14, 44] were employed by Rangel-Mondragon and Abas [42]in the design of decorative patterns inspired by Islamic art. They had effectively reinvented themedieval trade secrets of the craftsmen of fifteenth century Islamic mosques [28] who created byhand highly intricate mosaics closely resembling quasicrystal tilings only discovered by modernscience in the late twentieth century. Wang tilings [12, 14] were first introduced by Stam [48]in order to enable wave texture patches to cover water surfaces of arbitrary size without the ap-pearance of regularly repeating artifacts. Further computer graphics applications of non-periodictilings include texture mapping and synthesis [3, 23, 24, 25, 48, 49], photorealistic rendering us-ing environmental maps [24, 37], and non-photorealistic rendering using stippling [23, 24]. Forcomputer graphics applications, non-periodic tilings have usually been generated by geometricconstructs, such as matching rules and hierarchical substitution [11, 14]. In this work, we presentthe cut-and-project method of generating quasicrystals as an alternative algebraic approach tothe production of non-periodic tilings and point sets (Fig. 1). This algebraic approach has the a r X i v : . [ c s . C V ] J u l M. Grundland, J. Patera, Z. Mas´akov´a and N.A. Dodgson
Point set Voronoi diagramDelaunay graph Delaunay triangulation
Figure 1.
Quasicrystal tilings produced using spatial proximity graphs. In these visualizations, a non-periodic, rotationally symmetric point set (top left) is depicted as a planar tiling induced by a Voronoidiagram (top right), a Delaunay graph (bottom left), and a Delaunay triangulation (bottom right). Thisset of 1035 points comprises a cut-and-project quasicrystal derived from the standard root lattice ofthe non-crystallographic Coxeter group H . Its viewing window is a square centered at the origin withradius 1, while its acceptance window is a decagon centered at the origin with radius τ + τ , where τ is the golden ratio. The visualization of the quasicrystal tilings reveals some remarkable properties. It iswell known that quasicrystals can exhibit five and ten fold rotational symmetry, an impossible feat forany periodic tiling. Recently, it has been shown analytically that a quasicrystal Delaunay graph can yielda non-periodic tiling with four distinct tile shapes [29]. As illustrated by this visualization, a Delaunaytriangulation of a cut-and-project quasicrystal can yield a non-periodic tiling with just three distinct tileshapes. mage Sampling with Quasicrystals 3advantages of being straightforward to implement, easy to calculate, and readily amenable torigorous mathematical analysis. Moreover, it may be directly extended to higher dimensions aswell as adaptive sampling applications, although this is outside the scope of our present work.We choose to base our method on the algebra of the golden ratio τ = (1 + √ While a periodic point set is characterized by its translational symmetries, a non-periodic pointset admits no translational symmetries. For use in image sampling, we focus on non-periodicpoint sets that are determined by their inflation symmetries. In such a non-periodic pointset, a fixed configuration of sample sites can be repeated at different scales to generate a self-similar pattern. The simplest way of producing non-periodic point sets is to use hierarchicalsubstitution tilings [13, 14]. For instance, hierarchical substitution can be readily applied to thefamous Penrose tiling [12, 14, 37, 42, 44]. The strategy starts with a small set of polygonal tiles.The tile set is carefully designed so that each tile can be decomposed by geometric subdivisioninto smaller instances of itself and the other tiles. A hierarchy is formed whereby an existingtile becomes the parent of new child tiles. Starting with an initial configuration of the tilesscaled to cover the image plane, the tiling is refined through an iterative process of deflationand substitution. The tile vertices or centroids are used to derive a point set from the tiling.The choice of the initial configuration appears mirrored in the global structure and symmetryproperties of the tiling and the resulting point set.We focus on a more general class of non-periodic point sets corresponding to cut-and-projectquasicrystals [2, 40, 44]. The cut-and-project method was originally introduced by Meyer [31] inthe context of harmonic analysis and it was later adapted for generating quasicrystals by Moodyand Patera [34]. The Fibonacci chain and Penrose tilings can be regarded as special cases ofsuch quasicrystals. In our work, we employ the standard root lattice of the non-crystallographicCoxeter group H , a group of reflections taken from Lie algebra theory. This approach to qua-sicrystals can be used to relate discrete, non-periodic point sets and tilings (Fig. 1) with thelevel sets of continuous, non-periodic functions (Fig. 2). To produce a 2D cut-and-project qua-sicrystal, a 4D periodic lattice is projected on a suitable 2D plane that is irrationally orientedwith respect to the lattice. Using this method, we obtain a dense subspace consisting of inte-ger coefficient linear combinations of the vertices of a regular decagon centered on the origin,which are the roots of the non-crystallographic Coxeter group H . This construction ensuresthat the coordinates of all its elements can be expressed using only integers and an irrationalnumber, the golden ratio. To obtain a finite 2D quasicrystal, we select only those elements ofthe dense subspace contained in a specified bounded region, called the viewing window, that aremapped by an everywhere discontinuous algebraic transformation, called the star map, to an-other specified bounded region, called the acceptance window. Through the gradual expansionof a rotationally symmetric acceptance window, a quasicrystal sequence can be uniquely orderedby radial distance and angle from then center of the acceptance window in order to produce auniformly space-filling point set in the viewing window. As an important practical consequence,this property directly enables progressive sampling. It could also potentially enable adaptivesampling by varying the radius of the acceptance window according to an application depen-dent importance map defined for the viewing window. A 2D quasicrystal can also be expressedmage Sampling with Quasicrystals 5 Phase function contour map Phase function level set
Quasicrystal phase function: f ( z ) = (cid:80) j =0 exp (cid:0) πi ( ζ j ) · (2 τ z ) (cid:1) definedfor regular decagon vertices ζ j = exp (cid:0) πi j (cid:1) and golden ratio τ = (1 + √ u · v = ( uv + uv ) and complex conjugate u + u i = u − u i . Figure 2.
Quasicrystal construction using a continuous quasicrystal phase function. Observe thata quasicrystal point set is contained within a level set of a continuous phase function f : C → R definedin the complex plane [35]. This continuous phase function is formulated using the roots of the non-crystallographic Coxeter group H , which comprise the vertices ζ j of a regular decagon. To generatea quasicrystal point set, start by placing the first point at the origin. For each newly placed point x ∈ C ,consider the candidate points z = x + ζ j with j ∈ { , . . . , } , which are the vertices of a regular decagoncentered at x , and only accept the candidates for which the phase function f ( z ) ≥ T exceeds a desiredthreshold T that controls the density of the resulting point set. as a subset of a Cartesian product of two 1D quasicrystals, in accordance with the fact thatthe points that lie on any straight line through a 2D quasicrystal correspond to some linearlytransformed 1D quasicrystal. Therefore, in practice, a 2D quasicrystal (Fig. 4) can be generatedfrom a 2D lattice of 1D quasicrystals. Meanwhile, a 1D quasicrystal (Fig. 3) is produced bytaking a strip of a 2D periodic lattice, having finite width and irrational slope, and orthogonallyprojecting its points onto a line of the same slope. The resulting 1D quasicrystal is composedof at most three distinct tiles. It is easy to generate 1D quasicrystal points using an iterativenumerical algorithm. Alternatively, it is possible to exploit the self-similar structure of a 1Dquasicrystal, viewing it as the fixed point of a set of substitution rules that act recursively on afinite alphabet of possible tile arrangements.Based on the geometry and algebra of the golden ratio, these quasicrystal point sets exhibitsome remarkable properties. They display pentagonal and decagonal rotational symmetries,which cannot occur in any periodic point set. Originally, the theory of quasicrystals was moti-vated by solid state physics as a model of the non-periodic geometric structures that describe thesymmetries of a new kind of long-range atomic order discovered in certain metallic alloys [46].While translational symmetries define periodic crystals, inflation symmetries can be used to de-scribe quasicrystals based on algebraic irrational numbers, such as the golden ratio. When the M. Grundland, J. Patera, Z. Mas´akov´a and N.A. Dodgson Cut-and-project principle:
To generate a non-periodic point set in an n -dimensionalspace, take a region of a periodic point set in a2 n -dimensional space and orthogonally project it onto anirrationally oriented n -dimensional subspace. For the re-sulting point set in n -dimensional space to be discreterather than dense, the region of 2 n -dimensional space un-dergoing projection, typically a cylinder, must be boundedalong the directions of projection. Ruler and graph paper construction:
Start with the 2D integer lattice. Consider three parallellines with slope τ , the inner one passing through the origin(0 , c, d ),while the visible length of the strip defines the viewingwindow V = [ l, r ]. Observe that when the lattice pointscontained in the strip are orthogonally projected onto thecentral axis, the resulting 1D sequence Σ Ω is non-periodicif and only if the slope τ of the strip is irrational. Golden ratio: τ = (cid:1) √ (cid:2) ≈ .
618 and its conjugate τ (cid:1) = (cid:1) − √ (cid:2) ≈ − .
618 are the solutions of x = x + 1. Golden integers: Z [ τ ] = { a + bτ | a, b ∈ Z } is an Euclidean domain that is dense in R . Conjugate golden integers: Z [ τ (cid:1) ] = { a + bτ (cid:1) | a, b ∈ Z } is the set of conjugates ( a + bτ ) (cid:1) = a + bτ (cid:1) = a − bτ − , where τ (cid:1) + τ = 1 and τ (cid:1) τ = − Cut-and-project 1D quasicrystals: Σ Ω = (cid:3) a + bτ ∈ Z [ τ ] (cid:4)(cid:4)(cid:4) a + bτ (cid:1) ∈ Ω ∩ Z [ τ (cid:1) ] (cid:5) quasicrystal is specified by a bounded acceptance window Ω = [ c, d ). Duality of 1D quasicrystals: x k ∈ (Σ Ω ∩ V ) ⊂ Z [ τ ] restricted to the bounded viewing interval V = [ l, r ] implies a dual quasicrystal x (cid:1) k ∈ (Σ (cid:1) V ∩ Ω) ⊂ Z [ τ (cid:1) ] contained in the bounded acceptance interval Ω = [ c, d ). Translation and scaling of 1D quasicrystals: Σ [ c,d ) + λ = Σ [ c + λ (cid:1) ,d + λ (cid:1) ) for λ ∈ Z [ τ ] and ξ (cid:1) Σ [ c,d ) = Σ [ ξc,ξd ) for ξ = τ k and k ∈ Z . Stepping function for 1D quasicrystals:
Assume a standard acceptance window Ω = [ c, d ) such that 0 ∈ Ω and d − c ∈ [1 , τ )and observe that a 1D quasicrystal is an arrangement of just three possible tiles. x (cid:1) k +1 = x (cid:1) k + 1 if x (cid:1) k ∈ [ c, d − ⇒ x k +1 − x k = 1 (small tile) x (cid:1) k + 1 + τ (cid:1) if x (cid:1) k ∈ [ d − , c − τ (cid:1) ) ⇒ x k +1 − x k = 1 + τ (large tile) x (cid:1) k + τ (cid:1) if x (cid:1) k ∈ [ c − τ (cid:1) , d ) ⇒ x k +1 − x k = τ (medium tile) (cid:1)(cid:2) xσ ( x ) c d − c − τ (cid:1) dd + τ (cid:1) c +1 d Iterative construction of 1D quasicrystals:
For a standard acceptance window Ω = [ c, d ), containing the origin c ≤ < d inside an interval of suitablewidth 1 ≤ d − c < τ , a quasicrystal sequence x k ∈ Σ Ω is generated by successively applying the stepping function σ : Ω → Ω to obtain the conjugate of the right neighbor x (cid:1) k +1 = σ ( x (cid:1) k ) or the inverse stepping function σ − : Ω → Ωto obtain the conjugate of the left neighbor x (cid:1) k − = σ − ( x (cid:1) k ). Stepping algorithm for 1D quasicrystals:
1. Translate and scale the desired acceptance Ω and viewing V intervals to make the acceptance window Ωinto a standard acceptance window.2. Starting from x = 0, apply the stepping function x (cid:1) k +1 = σ ( x (cid:1) k ) and its inverse x (cid:1) k − = σ − ( x (cid:1) k ) to generatethe consecutive quasicrystal points until both edges of the translated and scaled viewing window are reached.3. Reverse the translation and scaling to obtain the desired quasicrystal Σ Ω . Figure 3.
Algorithm for 1D cut-and-project quasicrystals. mage Sampling with Quasicrystals 7
Figure 4.
Algorithm for 2D cut-and-project quasicrystals.
M. Grundland, J. Patera, Z. Mas´akov´a and N.A. Dodgson
Evolution of a 1D quasicrystal point set (cid:54) from [0,1) (cid:58)(cid:32) to [0, ) (cid:58) (cid:87)(cid:32) . Evolution of a 1D quasicrystal point set Σ Ω from Ω = [0 ,
1) to Ω = [0 , τ ).Evolution of a Delaunay graph of a 2D quasicrystal, where 10 new sites are added at each generation.Evolution of a Delaunay triangulation of a 2D quasicrystal, where the acceptance window is a decagon centeredat the origin that is expanded by a factor of τ , from radius τ + τ on the left and to radius τ + τ on the right. Figure 5.
Evolution of cut-and-project quasicrystals reveals their self-similar structure. mage Sampling with Quasicrystals 9acceptance window is a convex region, every point of a quasicrystal can be viewed as a centerof inflation symmetry. A quasicrystal can have no translational symmetries and no periodicsubsets. Moreover, it can be partitioned into subsets such that each subset forms a valid qua-sicrystal. Consider the local configurations of tiles in an infinite quasicrystal mosaic formed bya Voronoi diagram or a Delaunay triangulation. According to the repetitivity property impliedby suitable regularity conditions, each fragment is repeated infinitely many times in the mosaic.Yet no single fragment is ever sufficient to determine the structure of the whole mosaic becauseevery finite fragment, no matter its size, occurs in an uncountable infinity of nonequivalentmosaics. Furthermore, according to the Delaunay point set property, quasicrystals are bothuniformly discrete and relatively dense, creating space-filling tilings.As Delaunay point sets, quasicrystals are particularly well suited to image sampling. Theyenforce both a minimal and a maximal distance between each sample site and its closest neigh-boring site. In quasicrystal sampling, we rely on the golden ratio to ensure symmetry and selfsimilarity, which are generally absent when other than algebraic irrational numbers are used withthe cut-and-project method to produce non-periodic point sets. By taking this approach, we canendow a rendition with a decorative symmetry, which viewers may find attractive in the contextof non-photorealistic rendering. Compared with the regular grids of periodic point sets, theself-similar, space-filling structure of non-periodic point sets (Fig. 5) appears less monotonous,especially during progressive image rendering. In effect, the geometric structure of quasicrys-tal sampling eliminates the possibility of aliasing artifacts regularly repeating in the rendition.However, in quasicrystal sampling, fixed local configurations of sample sites can be repeated atmultiple places and orientations, albeit not at regular intervals, with the potential to yield somerecurring, anisotropic aliasing artifacts. Although we did do so in this work, for photorealisticimage reconstruction, it is preferable to avoid inducing global rotational symmetry in the sam-pling pattern, which is done by ensuring the viewing window does not contain the origin whenthe acceptance window is symmetric with respect to the origin.
We now present a qualitative evaluation of quasicrystal sampling in the context of non-adaptivesampling strategies for use in image representation. In this application, non-adaptive samplingstrategies serve as building blocks for interactive sampling, adaptive sampling, and importancesampling techniques. In effect, they dictate the placement of sample sites in image regions sam-pled at a uniform resolution. For this purpose, a non-adaptive sampling strategy should satisfyseveral image representation objectives. The sample sites should be distributed in a manner thatfairly and accurately represents the image. The sampling pattern should be evenly space-fillingin order to enable progressive image rendering. Without any preconceptions about the distribu-tion of visually salient features in the image, the same amount of information should be devotedto capturing each part of the image. Hence, the number of sample sites placed in any region ofthe image should be proportional to its area, so the sampling density remains the same through-out the image. Both globally and locally, the placement of sample sites should be uniform andisotropic while still allowing for a variety of different sample site configurations. To preventaliasing, the sample sites should not be arranged into fixed configurations that visibly repeatlocally or globally. To prevent clustering, a minimum distance between sample sites should bemaintained throughout the image. Assuming that the correlation between pixels decreases withdistance, a sample site should be placed close to the centroid of its Voronoi polygon to enable itssampled color to be most representative of its region of influence when the image is reconstructedusing a local interpolation technique. Naturally, the relative importance of these considerationsdepends on the requirements of a particular application. We have formulated these prioritiesbased on our previous work on a multiresolution image representation [16] designed to simultane-0 M. Grundland, J. Patera, Z. Mas´akov´a and N.A. Dodgsonously support both photorealistic image reconstruction and non-photorealistic image rendering.While previous qualitative surveys [10] and quantitative comparisons [47] of image samplinghave focused on applications in photorealistic image reconstruction, they did not cover qua-sicrystal sampling and farthest point sampling. Furthermore, their tests were not carried out onan image representation of digital photographs and they did not specifically address the needsof non-photorealistic image rendering.We compared quasicrystal sampling to a number of standard non-adaptive image samplingstrategies [10], reflecting different approaches to the inherent trade-off between aliasing andnoise. For our evaluation, we chose approaches that exemplify divergent aims in sampling.For each approach, we selected a representative implementation. As noted below, alternativeimplementations are certainly possible but they are likely to produce similar qualitative results.From the deterministic to the stochastic, we tested a range of sampling strategies (Fig. 7):1.
Periodic sampling [36] aims for global regularity. Our implementation relies on a squarelattice refined in scan line order. An alternative implementation could use a hexagonallattice, the densest periodic lattice in the plane.2.
Quasicrystal sampling [44] aims for local regularity. Our implementation relies on thecut-and-project method applied using the golden ratio. An alternative implementationcould use a Penrose tiling produced using a hierarchical substitution algorithm.3.
Farthest point sampling [9] aims for spatial uniformity. Our implementation relies onthe principle of progressively sampling at the point of least information, placing each newsample site at the point farthest from any preceding sample site, which is necessarily avertex of the Voronoi diagram of the preceding sample sites. An alternative implementationcould position sample sites to conform to a centroidal Voronoi diagram so that each samplesite is placed at the centroid of its Voronoi polygon.4.
Jittered sampling [22] aims for local variability. Our implementation relies on a fullrandom displacement of a square lattice refined in scan line order. An alternative imple-mentation could use a partial random displacement of a hexagonal lattice.5.
Quasirandom sampling [41] aims for low discrepancy. Our implementation relies on theHalton sequence. An alternative implementation could use a Sobol sequence.6.
Random sampling [10] aims for global variability. Our implementation relies on a uni-form distribution. An alternative implementation could use a random walk on a unitsquare with toroidal boundary conditions.To perform a qualitative evaluation of the image sampling strategies, we applied a numberof computer visualization techniques. For each non-adaptive sampling strategy, we visualize itssample sites (Fig. 7) in the spatial domain using a Voronoi diagram (Fig. 8) and in the frequencydomain using a Fourier power spectrum (Fig. 9). We examined the visual effects of applyingthe image sampling strategies in the context of various image rendering techniques that are usedin multiresolution image representations. For photorealistic image reconstruction [1], we testedthe accuracy of Shepard’s interpolation (Fig. 13), an inverse distance weighted interpolationmethod that applies the Voronoi diagram to determine the color of each pixel based on itsfour nearest sample sites, as well as Gouraud shading (Fig. 14), a piecewise linear interpolationmethod that applies the Delaunay triangulation to determine the color of each pixel based onits three surrounding sample sites. To quantitatively assess the results of these widely used localinterpolation techniques, we relied on the peak signal-to-noise ratio (PSNR). This standardimage fidelity metric [45] estimates the accuracy of the rendition according to the negativelogarithm of the mean squared error between the rendered and actual RGB color values ofmage Sampling with Quasicrystals 11
Figure 6.
Techniques for visualizing the spatial configuration of sample sites: Voronoi diagram (farleft), Delaunay triangulation (center left), Delaunay star shapes (center right), and our mosaic renderingtechnique (far right). the pixels. Hence, higher peak signal-to-noise ratio scores are considered better. For non-photorealistic image rendering [16], we experimented with a simple “paint strokes” renderingstyle (Fig. 12), which applies geometric subdivision to the Delaunay triangulation of the samplesites and then uses linear and nonlinear interpolation to emphasize the transitions between thesampled colors.Our mosaic rendering style (Fig. 6) offers a new computer visualization tool for evaluating thespatial properties of point sets, such as the sample sites produced by the various image samplingstrategies. To produce a mosaic rendering, we apply geometric subdivision to the Delaunaytriangulation of the sample sites. The midpoints of the edges of each Delaunay triangle are joinedto form three outer triangles and one inner triangle. Each outer triangle is rendered with thecolor sampled at its vertex of the original Delaunay triangle, while the central triangle is coloredblack. As each sample site is represented by a star-shaped polygonal tile, the resulting mosaicappears packed as tightly as possible, with the black central triangles serving as grout betweenthe tiles. As a sample site’s local neighborhood (Fig. 6, far left) comprises the surrounding samplesites connected to it by edges in the Delaunay triangulation (Fig. 6, center left), the sample site’smosaic tile (Fig. 6, far right) is shaped to reflect the star of its surrounding Delaunay triangles(Fig. 6, center right). For instance, a sample site’s mosaic tile is a convex or concave polygonaccording to whether its neighboring sites are arranged in a convex or concave configuration.As a visualization tool, the advantage of mosaic rendering is that the layout of the star-shapedmosaic tiles makes the spatial properties of a sampling strategy easier to see at a glance than thetriangles of a Delaunay triangulation or the convex polygons of a Voronoi diagram. The sizesof the mosaic tiles are indicative of the uniformity of sampling, as coarsely sampled regions giverise to large tiles and finely sampled regions give rise to small tiles, making common defects suchas clustering, undersampling, and oversampling easy to detect. The orientations of the mosaictiles are indicative of the isotropy of sampling, as the preferred directions of the sampling arerevealed in the preferred rotations of the tiles, making global or local grid structures easy todetect. The shapes of the mosaic tiles are indicative of the heterogeneity of sampling, as thelocal configurations of neighboring sites uniquely determine the tile polygons, making repetitivepatterns easy to detect. For instance, farthest point sampling produces tiles of uniform size andsimilar shape to create the appearance of a pebble mosaic, while quasicrystal sampling yieldsa decorative tiling with just a small set of possible tile shapes.Our qualitative evaluation of image sampling strategies uses seven criteria (Fig. 10) knownto affect the visual quality of photorealistic image reconstruction and non-photorealistic imagerendering [16]:1.
Accurate reconstruction requires the rendition to faithfully represent the likeness of theoriginal image. It is a necessary but not sufficient condition of success in both photorealis-tic and non-photorealistic image rendering. This objective appears to be closely associated2 M. Grundland, J. Patera, Z. Mas´akov´a and N.A. Dodgson
Figure 7.
Non-adaptive sampling strategies: periodic sampling (top left), quasicrystal sampling (topcenter), farthest point sampling (top right), jittered sampling (bottom left), quasirandom sampling (bot-tom center), and random sampling (bottom right). Sampling starts with the dark blue sites and finisheswith the light green sites. with uniform coverage and centroidal regions. It is assessed by measuring the peak signal-to-noise ratio for the results of photorealistic image reconstruction (Figs. 13 and 14). Itsmost pronounced effects can also be observed in the results of non-photorealistic imagerendering (Fig. 12). When the resolution of sampling is uniform, periodic sampling yieldsthe most accurate image interpolation (in Fig. 14, for the regular square grid, this takesplace when there are 33 = 1089, 65 = 4225, and 129 = 16641 samples). However,when regions of varying resolution arise during progressive refinement, the accuracy ofperiodic sampling can substantially deteriorate. Similar behavior is observed in jitteredsampling since it applies random perturbations to a periodic point set. Farthest pointsampling produces image reconstructions that are nearly as accurate as periodic sampling,but its performance does not diminish during progressive refinement. Intermediate ac-curacy is offered by quasicrystal sampling, which appears to be slightly more accuratethan quasirandom sampling. The least accurate reconstructions are produced by jitteredand random sampling. Given its popularity in computer graphics implementations, thepoor performance of jittered sampling is rather disappointing. Of course, the accuracy ofjittered sampling can always be made closer to that of periodic sampling by reducing theamount of random displacement, which risks reintroducing the aliasing artifacts of periodicsampling. In general, accuracy is improved by uniformity and reduced by randomness, aneffect that can be readily seen as producing tight or loose image stylization.2. Progressive ref inement requires the sample sites to smoothly fill the available space,avoiding abrupt changes in appearance as new sample sites are sequentially added tothe rendition. This objective serves to enable a multiresolution image representation tosupport progressive rendering of compressed images based on an incremental sampling ofmage Sampling with Quasicrystals 13the image data. It is assessed by examining the spatial layout of the sequence of samplesites (Fig. 7). Under ideal circumstances, progressive refinement should yield a smoothcurve for the peak signal-to-noise ratio (Fig. 14). The best progressive refinement resultsare produced by farthest point sampling and quasicrystal sampling, as these methodsmaintain a uniform sampling density by ensuring that new sample sites are placed in thelargest empty spaces between the existing sample sites. Quasirandom sampling provesslightly less proficient, as it places some sample sites very close together while keepingothers far apart. Random sampling is even less effective due to its tendency to locallycluster sample sites. The regular grids used in periodic and jittered sampling are notsuitable for smooth progressive refinement, especially when they are refined in scan lineorder. While other refinement schemes can be applied to regular grids, such as refinementin random order, their intrinsic symmetry makes it difficult to smoothly increase thesampling density throughout the image.3.
Uniform coverage requires the sample sites to be evenly distributed regardless of posi-tion, avoiding configurations that place sample sites too close or too far from their nearestneighbors. This objective is assessed using Voronoi diagrams (Fig. 8) as well as mosaicrenderings (Fig. 11). Its effect determines the sizes of brush strokes in non-photorealisticimage rendering (Fig. 12). Uniform coverage is associated with a Fourier power spectrum(Fig. 9) that displays an empty ring around the central spike, as low frequencies are at-tenuated in favor of a threshold frequency corresponding to the most commonly observednearest neighbor distance between sample sites. Although a blue noise spectrum can en-sure uniform coverage, it is not a necessary condition. When the resolution of samplingis uniform, periodic sampling generates uniform coverage, as its mosaic tiles are all ex-actly the same size. However, periodic sampling cannot sustain uniform coverage duringprogressive refinement. By design, farthest point sampling maintains uniform coverageat all times, as its mosaic tiles are all approximately the same size. Quasicrystal sam-pling maintains nearly as uniform coverage, as its mosaic tiles are limited to just a fewcomparable sizes. While quasirandom and jittered sampling strive to uphold a uniformdensity of sampling, they nevertheless are less effective at providing uniform coverage, astheir mosaic tiles come in many sizes. In the case of jittered sampling, uniform coveragecan be improved by reducing the amount of random displacement. Random samplingdoes not give uniform coverage, as its mosaic tiles exhibit the greatest range of differentsizes.4.
Isotropic distribution requires the sample sites to be evenly distributed regardless oforientation, avoiding configurations that align sample sites along globally or locally pre-ferred directions. This objective is assessed using Voronoi diagrams (Fig. 8) as well asmosaic renderings (Fig. 11). Its effect is to determine the orientations of brush strokesin non-photorealistic image rendering (Fig. 12). An isotropic distribution produces aFourier power spectrum (Fig. 9) that displays a rotational symmetry around the cen-tral spike, as the power at each frequency does not depend on its orientation. Al-though a blue noise spectrum can ensure isotropic distribution, it is not a necessarycondition. Random sampling is the most isotropic, as its sample sites are both locallyand globally uncorrelated. Farthest point and jittered sampling are nearly as isotropic,as their sample sites can exhibit slight local alignment. Farthest point samples canappear to be placed in roughly hexagonal local configurations. Jittered samples canappear to retain some of the structure of the underlying square grid, as the isotropyof jittered sampling reflects the amount of random displacement used to generate thesampling. Quasirandom sampling has intermediate isotropy, as its sample sites can ex-hibit slight global alignment, which can be verified in the lack of radial symmetry in4 M. Grundland, J. Patera, Z. Mas´akov´a and N.A. Dodgson
Figure 8.
Voronoi diagrams of image sampling strategies applied to a color spiral test image: peri-odic sampling (top left), quasicrystal sampling (top center), farthest point sampling (top right), jitteredsampling (bottom left), quasirandom sampling (bottom center), and random sampling (bottom right). its Fourier power spectrum. Periodic sampling and quasicrystal sampling do not haveisotropic distribution since their sample sites are globally aligned along predeterminedaxes.5.
Blue noise spectrum requires the sample sites to be distributed similarly to a Poissondisk distribution, a random point field conditioned on a minimum distance between thepoints. According to this objective, for an image sampling strategy to provide effectiveantialiasing for image rendering, it should attempt to mimic the idealized distribution ofphotoreceptors in the human eye. Usually implying both uniform coverage and isotropicdistribution, a blue noise spectrum is highly desirable in many computer graphics ap-plications, particularly photorealistic image reconstruction. It is assessed by examiningthe Fourier power spectrums of the sampling strategies (Fig. 9) for a radially symmetricprofile that concentrates noise in the high frequencies while attenuating the power of thelow frequencies, thereby eliminating the aliasing artifacts associated with low frequencypatterns that can appear distracting to the eye. In effect, a blue noise spectrum exhibits adisk of low power around the origin, surrounded by roughly constant power at the higherfrequencies. Its effects can be judged according to the amount of aliasing present in photo-realistic image reconstruction (Fig. 13) and non-photorealistic image rendering (Fig. 12).Farthest point sampling has a Fourier power spectrum that is closest to a blue noise spec-trum. Jittered sampling attempts to replicate the blue noise spectrum, but it does notclearly exhibit the threshold frequency ripple around the central spike. Quasirandom iseven less successful because its Fourier power spectrum lacks radial symmetry. Periodic,quasicrystal, and random sampling have Fourier power spectrums that do not resemble theblue noise spectrum. The Fourier power spectrums of periodic and quasicrystal samplingmage Sampling with Quasicrystals 15
Figure 9.
Fourier power spectrums of image sampling strategies: periodic sampling (top left), qua-sicrystal sampling (top center), farthest point sampling (top right), jittered sampling (bottom left),quasirandom sampling (bottom center), and random sampling (bottom right). reflect the spatial structures and directional symmetries of the lattices used to place thesample sites. On the other hand, the white noise spectrum of random sampling assignsroughly the same power to all frequencies.6.
Centroidal regions require sample sites to be well centered with respect to their Voronoipolygons, approximating a centroidal Voronoi diagram. Typically associated with uniformcoverage and accurate reconstruction, this objective is popular in non-photorealistic imagerendering. Sampling strategies, such as periodic sampling, that produce centroidal regionscan still be prone to aliasing artifacts since centroidal regions do not guarantee a blue noisespectrum. Centroidal regions can be readily assessed using the Voronoi diagrams (Fig. 8).The effects can also be observed in the shapes of tiles in mosaic rendering (Fig. 11) andbrush strokes in non-photorealistic image rendering (Fig. 12). Periodic sampling, placingeach sample site at the same distance from all of its nearest neighbors, generates exact cen-troidal Voronoi regions. Quasicrystal and farthest point sampling produce approximatelycentroidal Voronoi regions. To place sample sites close to the center of their Voronoi poly-gons, quasicrystal sampling relies on local symmetries while farthest point sampling relieson nearest neighbor distance. Jittered sampling and quasirandom sampling have difficultyensuring centroidal regions because they are less effective at maintaining a minimal near-est neighbor distance. Finally, random sampling can only produce centroidal regions bychance.7.
Heterogeneous conf igurations require sample sites to be placed in a variety of differentlocal arrangements, avoiding regularly or randomly repeating the same sampling patterns.While this objective is not traditionally a concern in photorealistic image reconstruc-tion, it helps to prevent non-photorealistic image rendering from appearing too perfect,seemingly mechanical and artificial. For instance, it helps to give a vibrant appearanceto brush stroke rendering. Typically, when sampling strategies yield centroidal regions,6 M. Grundland, J. Patera, Z. Mas´akov´a and N.A. Dodgson
Sampling Accurate Progressive Uniform Isotropic Blue Noise Centroidal HeteroeneousStrategies Reconstruction Refinement Coverage Distribution Spectrum Regions ConfigurationsPeriodic (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70)
Quasicrystal (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70)
Farthest Point (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70)
Jittered (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70)
Quasirandom (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70)
Random (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70) (cid:70)(cid:70) (cid:70) (cid:70) (cid:70)
Superior (cid:70) (cid:70) (cid:70)
Good (cid:70) (cid:70)
Fair (cid:70)
Poor
Figure 10.
Qualitative evaluation of image sampling strategies. they also tend to produce homogeneous configurations and vice versa, illustrating an ap-parent trade-off between these competing objectives. Heterogeneous configurations areassessed using the Voronoi diagrams (Fig. 8) and mosaic renderings (Fig. 11). Their effectis also visible in the arrangement of brush strokes in non-photorealistic image rendering(Fig. 12). Random sampling produces the most heterogeneous local configurations, as itsmosaic tiles exhibit the greatest variety of shapes. Jittered and quasirandom sampling arenearly as heterogeneous, since their mosaic tiles are almost as widely varied, though fewof them are exceptionally large in size. Quasicrystal and farthest point sampling are farless heterogeneous. By upholding local symmetries, quasicrystal sampling causes samplesites to have only a few possible local configurations, resulting in mosaic tiles that haveonly a few possible shapes. By upholding nearest neighbor distance, farthest point sam-pling causes sample sites to have similar local configurations, resulting in mosaic tiles thatlook very much alike, mostly convex and rounded. Based on repetitions of a single localconfiguration, periodic sampling is entirely homogeneous.Our qualitative analysis (Fig. 10) indicates that quasicrystal sampling offers a useful com-promise between the ordered behavior of standard periodic sampling using a regular squarelattice and the disordered behavior of standard Monte Carlo sampling using jittered, quasiran-dom, or random sampling. Compared with periodic sampling, quasicrystal sampling displaysa greater variety of local sample site configurations resulting in smoother progressive refinement,although its sampling patterns are somewhat less uniform, leading to lower accuracy of imagereconstruction. Compared with jittered, quasirandom, and random sampling, quasicrystal sam-pling displays more uniform coverage resulting in better accuracy of image reconstruction, al-though its sampling patterns are anisotropic, exhibiting significantly less variety of local samplesite configurations. By virtue of its deterministic construction, quasicrystal sampling does notsuffer from the variability that can affect the results of random sampling, jittered sampling,and, to a much lesser degree, farthest point sampling. Nevertheless, its lack of a blue noisepower spectrum renders it rather susceptible to aliasing artifacts. Research on quasicrystalsampling based on the Penrose tiling [37] suggests that it may be possible to partially alleviatethis problem by taking advantage of the symmetries and the repetitions of the local sample siteconfigurations in order to systematically displace the sample sites in a manner that improvesthe spectral properties of the sampling pattern.Based on our qualitative evaluation of the various non-adaptive sampling strategies (Fig. 10),we recommend a blue noise sampling strategy, such as farthest point sampling, for general usein image representation. In particular, farthest point sampling does not perform poorly onany of our seven evaluation criteria. Overall, our qualitative evaluation of non-adaptive sam-pling strategies is in broad agreement with previous studies, which did not consider quasicrystalsampling. They emphasized the importance of Poisson disk distributions [10] and low discrep-mage Sampling with Quasicrystals 17ancy distributions [47], which are exemplified in our evaluation by farthest point sampling andquasirandom sampling respectively. Hence, the good overall performance of these two tech-niques should come as no surprise. Farthest point sampling performs better than quasirandomsampling on six out of the seven evaluation criteria. For the majority of our evaluation crite-ria, quasicrystal sampling performs no better than farthest point sampling and no worse thanquasirandom sampling. Nevertheless, from a practical point of view, quasicrystal samplingis significantly simpler to implement and calculate than farthest point sampling, which relieson maintaining complex geometric data structures to keep track of the vertices of a Voronoidiagram. This could be an important consideration for imaging applications on mobile de-vices that have limited processing and storage capabilities. From a theoretical point of view,the deterministic algebraic construction of quasicrystals renders their sampling patterns par-ticularly well suited to mathematical analysis. Presenting possibilities for future research, thecut-and-project method could be adapted for higher dimensional sampling or adaptive samplingapplications.In future work, it would also be interesting to explore the relationship between local sym-metry and sampling quality. The cut-and-project method can be used to generate non-periodicpoint sets with different symmetries, not just the pentagonal and decagonal symmetries asso-ciated with the golden ratio, as shown in this work. Just as for periodic sampling it wouldbe interesting to compare the image reconstruction accuracy of square and hexagonal grids,for non-periodic sampling it would be interesting to compare our decagonal quasicrystal tilingwith the dodecagonal Socolar tiling, which was recently proposed for use in sampling applica-tions [38].
Cut-and-project quasicrystals present new possibilities for image sampling in computer gra-phics. This non-periodic sampling approach deterministically generates uniformly space-fillingpoint sets, ensuring that sample sites are evenly distributed throughout the image. It offersa useful compromise between predictability and randomness, between the standard periodicsampling and the standard Monte Carlo sampling methods. Although blue noise sampling cangenerate higher quality sampling patterns for photorealistic image reconstruction, quasicrys-tal sampling may prove much simpler to implement and calculate. In the context of non-photorealistic image rendering, quasicrystal sampling may prove attractive for its symmetryproperties.
Acknowledgements
We are grateful to the referees for valuable remarks that led to improvements of the originalmanuscript. We also acknowledge the financial support of the Natural Sciences and EngineeringResearch Council of Canada, le Fonds Qu´eb´ecois de la Recherche sur la Nature et les Technolo-gies, as well as the grants MSM6840770039 and LC06002 of the Ministry of Education of theCzech Republic, and GA201/09/0584 of the Czech Science Foundation. We are also grateful forthe support of the MIND Research Institute and the Merck Frosst Company. Mark Grundlandfurther acknowledges the financial assistance of the Celanese Canada Internationalist Fellow-ship, the British Council, the Cambridge Commonwealth Trust, and the Overseas ResearchStudent Award Scheme. The images were provided by FreeFoto.com and the Waterloo BragZone.8 M. Grundland, J. Patera, Z. Mas´akov´a and N.A. Dodgson
Periodic sampling Quasicrystal samplingFarthest point sampling Jittered samplingQuasirandom sampling Random sampling
Figure 11.
Image sampling strategies rendered using the mosaic style (4225 samples ≈ . mage Sampling with Quasicrystals 19 Periodic sampling Quasicrystal samplingFarthest point sampling Jittered samplingQuasirandom sampling Random sampling
Figure 12.
Image sampling strategies rendered using the “paint strokes” style (4225 samples ≈ . Periodic sampling (PSNR = 18 .
58) Quasicrystal sampling (PSNR = 18 . .
55) Jittered sampling (PSNR = 18 . .
24) Random sampling (PSNR = 17 . Figure 13.
Image sampling strategies rendered using Shepard interpolation (4225 samples ≈ . mage Sampling with Quasicrystals 21 Sample Sites (N) R e c o n s t r u c t i o n A cc u r a c y ( P S N R ) Periodic samplingQuasicrystal samplingFarthest point samplingJittered samplingQuasirandom samplingRandom sampling
Fidelity of Image Reconstruction: Lena
Sample Sites (N) R e c o n s t r u c t i o n A cc u r a c y ( P S N R ) Periodic samplingQuasicrystal samplingFarthest point samplingJittered samplingQuasirandom samplingRandom sampling
Fidelity of Image Reconstruction: Park
Figure 14.
Quantitative evaluation of image sampling strategies rendered using Gouraud shading.
Supplementary materials
For a full resolution version of this paper, along with supplementary materials, please visit: . References [1] Amidror I., Scattered data interpolation methods for electronic imaging systems: a survey,
J. ElectronicImaging (2002), 157–176.[2] Chen L., Moody R.V., Patera J., Non-crystallographic root systems, in Quasicrystals and Discrete Geometry(Toronto, ON, 1995), Fields Inst. Monogr. , Vol. 10, Amer. Math. Soc., Providence, RI, 1998, 135–178.[3] Cohen M.F., Shade J., Hiller S., Deussen O., Wang tiles for image and texture generation, in Proceedingsof SIGGRAPH, 2003, 287–294.[4] Cook R.L., Stochastic sampling in computer graphics,
ACM Trans. Graphics (1986), 51–72.[5] Deussen O., Hiller S., van Overveld C., Strothotte T., Floating points: a method for computing stippledrawings, in Proceedings of EUROGRAPHICS, 2000, 41–50.[6] Dippe M.A.Z., Wold E.H., Antialiasing through stochastic sampling, in Proceedings of SIGGRAPH, 1985,69–78.[7] Du Q., Faber V., Gunzburger M., Centroidal Voronoi tessellations: applications and algorithms, SIAM Rev. (1999), 637–676.[8] Dunbar D., Humphreys G., A spatial data structure for fast Poisson-disk sample generation, in Proceedingsof SIGGRAPH, 2006, 503–508.[9] Eldar Y., Lindenbaum M., Porat M., Zeevi Y.Y., The farthest point strategy for progressive image sampling, IEEE Trans. Image Process. (1997), 1305–1315.[10] Glassner A., Principles of digital image synthesis, Vol. 1, Morgan Kaufmann, 1995.[11] Glassner A., Aperiodic tiling, IEEE Computer Graphics and Applications (1998), no. 3, 83–90.[12] Glassner A., Penrose tiling, IEEE Computer Graphics and Applications (1998), no. 4, 78–86.[13] Goodman-Strauss C., Aperiodic hierarchical tilings, in Foams, Emulsions, and Cellular Materials (Carg`ese,1997), NATO Adv. Sci. Inst. Ser. E Appl. Sci. , Vol. 354, Kluwer Acad. Publ., Dordrecht, 1999, 481–496.[14] Grunbaum B., Shephard G.C., Tilings and patterns, WH Freeman, 1987.[15] Grundland M., Style and content in digital imaging: Reconciling aesthetics with efficiency in image repre-sentation, VDM, 2008.[16] Grundland M., Gibbs C., Dodgson N.A., Stylized multiresolution image representation,
J. Electronic Imaging (2008), 013009, 1–17.[17] Hausner A., Simulating decorative mosaics, in Proceedings of SIGGRAPH, 2001, 573–580.[18] Hausner A., Pointillist halftoning, in Proceedings of the International Conference on Computer Graphicsand Imaging, 2005, 134–139.[19] Hiller S., Deussen O., Keller A., Tiled blue noise samples, in Proceedings of Vision, Modeling and Visuali-zation, 2001, 265–271.[20] Hiller S., Hellwig H., Deussen O., Beyond stippling – methods for distributing objects on the plane, inProceedings of EUROGRAPHICS, 2003, 515–522.[21] Jones T.R., Efficient generation of Poisson-disk sampling patterns, J. Graphics Tools (2006), no. 2,27–36.[22] Klassen R.V., Filtered jitter, Computer Graphics Forum (2000), no. 4, 223–230.[23] Kopf J., Cohen-Or D., Deussen O., Lischinski D., Recursive Wang tiles for real-time blue noise, in Pro-ceedings of SIGGRAPH, 2006, 509–518.[24] Lagae A., Dutre P., A procedural object distribution function, ACM Trans. Graphics (2005), 1442–1461.[25] Lagae A., Dutre P., An alternative for Wang tiles: colored edges versus colored corners, ACM Trans.Graphics (2006), 1442–1459.[26] Lagae A., Dutre P., A comparison of methods for generating Poisson disk distributions, Computer GraphicsForum (2008), no. 1, 114–129. mage Sampling with Quasicrystals 23 [27] Lagae A., Kaplan C.S., Fu C.-W., Ostromoukhov V., Deussen O., Tile-based methods for interactive appli-cations, SIGGRAPH 2008 Class Notes, ACM, 2008.[28] Lu P.J., Steinhardt P.J., Decagonal and quasi-crystalline tilings in medieval Islamic architecture, Science (2007), no. 5815, 1106–1110.[29] Mas´akov´a Z., Patera J., Zich J., Classification of Voronoi and Delone tiles of quasicrystals. III. Decagonalacceptance window of any size,
J. Phys. A: Math. Gen. (2005), 1947–1960.[30] McCool M., Fiume E., Hierarchical Poisson disk sampling distributions, in Proceedings of Graphics Interface,1992, 94–105.[31] Meyer Y., Algebraic numbers and harmonic analysis, North-Holland, 1972.[32] Mitchell D.P., Spectrally optimal sampling for distribution ray tracing, in Proceedings of SIGGRAPH, 1991,157–164.[33] Mojsilovic A., Soljanin E., Color quantization and processing by Fibonacci lattices, IEEE Trans. ImageProcess. (2001), 1712–1725.[34] Moody R.V., Patera J., Quasicrystals and icosians, J. Phys. A: Math. Gen. (1993), 2829–2853.[35] Moody R.V., Patera J., Dynamical generation of quasicrystals, Lett. Math. Phys. (1996), 291–300.[36] Ostromoukhov V., Mathematical tools for computer-generated ornamental patterns, Lecture Notes in Com-puter Science , Vol. 1375, 1998, 193–223.[37] Ostromoukhov V., Donohue C., Jodoin P.M., Fast hierarchical importance sampling with blue noise proper-ties, in Proceedings of SIGGRAPH, 2004, 488–495.[38] Ostromoukhov V., Building 2D low-discrepancy sequences for hierarchical importance sampling using dode-cagonal aperiodic tiling, in Proceedings of GRAPHICON, 2007, 139–142.[39] Ostromoukhov V., Sampling with polyominoes, in Proceedings of SIGGRAPH, 2007, 078, 1–6.[40] Patera J., Non-crystallographic root systems and quasicrystals. in The Mathematics of Long-Range Aperio-dic Order (Waterloo, ON, 1995),
NATO Adv. Sci. Inst. Ser. C Math. Phys. Sci. , Vol. 489, Kluwer Acad.Publ., Dordrecht, 1997, 443–465.[41] Press W., Teukolsky S.A., Vetterling W.T., Flannery B.P., Numerical recipes in C, 2nd ed., CambridgeUniversity Press, 1992.[42] Rangel-Mondragon J., Abas S.J., Computer generation of Penrose tilings,
Computer Graphics Forum (1988), no. 1, 29–37.[43] Secord A., Weighted Voronoi stippling, in Proceedings of the Second International Symposium on Non-Photorealistic Animation and Rendering, 2002, 37–43.[44] Senechal M., Quasicrystals and geometry, Cambridge University Press, Cambridge, 1995.[45] Sharma G., Digital color imaging handbook, CRC Press, 2003.[46] Shechtman D., Blech I., Gratias D., Cahn J.W., Metallic phase with long-range orientational order and notranslational symmetry, Phys. Rev. Lett.53