Sibson’s formula for higher order Voronoi diagrams

Sibson’s formula for higher order Voronoi diagrams Mercè Claverol1, Andrea de las Heras-Parrilla1, Clemens Huemer1, and Dolores Lara2 1 Universitat Politècnica de Catalunya [email protected], [email protected], [email protected] 2 Centro de Investigación y de Estudios Avanzados [email protected] Abstract Let S be a set of n points in general position in Rd . The orderk Voronoi diagram of S , Vk ( S ), is a subdivision of Rd into cells whose points have the same k nearest points of S . Sibson, in his seminal paper from 1980 (A vector identity for the Dirichlet tessellation), gives a formula to express a point Q of S as a convex combination of other points of S by using ratios of volumes of the intersection of cells of V2 ( S )and the cell of Q in V1 ( S ). The natural neighbour interpolation method is based on Sibson’s formula. We generalize his result to express Q as a convex combination of other points of S by using ratios of volumes from Voronoi diagrams of any given order. 1 Introduction Let S be a set of n points in general position in Rd , meaning no m of them lie in a ( m− 2)- dimensional flat for m = 2 , 3 , ..., d + 1 and no d + 2 of them lie in the same d -sphere, and let k be a natural number with 1 ≤k≤n− 1. Let σd denote the Lebesgue measure on Rd , to simplify we just write σ. The orderk Voronoi diagram of S , Vk ( S ), is a subdivision of Rd into cells such that points in the same cell have the same k nearest points of S . Thus, each cell f ( Pk )of Vk ( S )is defined by a subset Pk of S of k elements, where each point of f ( Pk )has Pk as its k closest points from S . Similarly, the ordered Voronoi diagram of order k of S , OVk ( S ), can be defined as a subdivision of Rd into cells such that points in the same cell have the same ordered k nearest points of S . Thus, each cell f ( ⟨Pk⟩ )of OVk ( S )is defined by an ordered subset ⟨Pk⟩ of size k of S , where the points are arranged in order of proximity starting from the closest to the farthest. Note that, by definition, the union of all the cells of OVk ( S )corresponding to the different permutations of a fixed subset of length k of S is the cell, f ( Pk ), associated to such subset in the (ordinary) order-kVoronoi diagram, Vk(S). See Figure 1. For the orderk Voronoi diagram of S , the region Rk ( ℓ )of Qℓ∈S is defined as the set of cells of Vk ( S )that have the point Qℓ as one of their k nearest neighbours. See Figure 2. For OVk ( S )we can define these regions in the same way. These regions are not necessarily convex but star-shaped, see [ 2 , 4 , 10 , 16 ], and it is known that R1 ( ℓ )is contained in the kernel of Rk ( ℓ ); see [ 3 ]. Also, these regions are related to Brillouin zones. For a given k , the region Rk ( ℓ ) \Rk−1 ( ℓ )is known as a Brillouin zone of Qℓ . Brillouin zones have been studied mainly for lattices but also for arbitrary discrete sets, see e.g. [6, 17]. Local coordinates based on Voronoi diagrams were introduced by Sibson [ 13 ]. He states that, given a set S of n points of Rd in general position, a point Qℓ∈S can be expressed as a convex combination of its nearest points of S . This is described next. Cells of V2 ( S )that intersect f ( {Qℓ} )in V1 ( S )are of the form f ( {Qℓ, Qj} ), i.e., cells defined by Qℓ and another point Qj , that we call its natural neighbour. These intersections give ratios of volumes which are the coefficients multiplying the corresponding natural neighbours in the convex 40th European Workshop on Computational Geometry, Ioannina, Greece, March 13–15, 2024. This is an extended abstract of a presentation given at EuroCG’24. It has been made public for the benefit of the community and should be considered a preprint rather than a formally reviewed paper. Thus, this work is expected to appear eventually in more final form at a conference with formal proceedings and/or in a journal. 21:2 Sibson’s formula for higher order Voronoi diagrams 15 12 13 31 34 43 41 14 51 25 52 21 45 54 23 32 Q3 Q4 Q5 Q1 Q2 135 153 134 314 341 431 132 123 125 152 215 213 312 251 521 512 514 154 541 451 415 145 143 413 231 321 513 Q3 Q4 Q5 Q1 Q2 Figure 1 For a set S = {Q1,··· , Q5} of five points in R2 . Each cell of OVk ( S )is labeled by the indices of its k nearest points of S . V1 ( S )is shown in black, V2 ( S )in green, and V3 ( S )in orange colour. Left: The cells of OV2 ( S )with the same nearest neighbour Qi from S form the cell f ( {Qi} ) in V1 ( S ). The cells of OV2 ( S )with the same subset P2 of two points of S (in any order) form the cell f(P2)of V2(S). Right: OV3(S)is shown together with V1(S),V2(S), and V3(S). Q2 Q3 Q4 Q5 Q6 Q8 Q9 Q0 V2(S) Q7 Q1 V1(S) R1(1) R2(1) Figure 2 R1 (1) is the cell f ( {Q1} )in V1 ( S ). R2 (1) is the union of cells of V2 ( S )that have Q1 as one of its two nearest neighbours. R1(1) ⊂R2(1). combination that expresses Qℓ . Volumes σ ( f ( {Qℓ, Qj} ) ∩f ( {Qℓ} )) are equal to the volumes given by the intersection of the cells of V1(S\{Qℓ})and f(Qℓ)in V1(S), see Figure 3. M. Claverol, A. de las Heras-Parrilla, C. Huemer and D. Lara 21:3 ▶Theorem 1.1. (Local coordinates property [13]). For a bounded cell f({Qℓ})of V1(S), Qℓ=X j=ℓ σ(f({Qℓ, Qj})∩f({Qℓ})) σ(f({Qℓ})) Qj(1) Sibson’s formula has been used to define the natural neighbour interpolation method [ 14 ]. Given a set of points and a function, this interpolation method provides a smooth approximation of new points to the function. Sibson’s algorithm uses the closest subset of the input set S\{Qℓ} to interpolate a query point, Qℓ , and applies weights based on the ratios of volumes provided by Theorem 1.1. Local coordinates and the natural neighbour interpolation method have been studied e.g. in [ 5 , 11 , 15 ], and they have many applications such as reconstruction of a surface from unstructured data or interpolation of rainfall data, see [9, 15]. Q3 Q2 Q` Q1 Q6 Q5 Q4 Figure 3 In R2 . Left: The initial Voronoi diagram V1 ( S\{Qℓ} )without query point Qℓ . Right: Colored areas given by the intersections of f ( {Qℓ} )and the cells of V1 ( S\{Qℓ} ), are the same as the ones given by the intersections of the cells of V2(S)(shown in dashed) with the cell f({Qℓ}). Aurenhammer gave a generalization of Sibson’s result to Voronoi diagrams of higher order, and more generally to power diagrams, see [ 1 ]. Aurenhammer’s formula allows to write a point Qℓ of S as a linear combination of other points of S . We state this in Theorem 2.1 and Corollary 2.2 below. The formula in Theorem 2.1 is defined in terms of OVk+1 ( S ). It is restated in Corollary 2.2 in terms of intersections of cells of Vk−1 ( S )and Vk+1 ( S )with a cell of Vk(S). This formula works for a bounded cell of Vk(S). Our main contribution is another generalization of Sibson’s result, stated in Theorem 2.3. In this theorem, we express a point Qℓ∈S as a convex combination of its neighbours of S using ratios of volumes in the region Rk ( ℓ ). Similar to Sibson’s formula that required the cell of the point Qℓ to be bounded, our formula requires its region Rk ( ℓ )to be bounded. For the case k= 1, Theorem 2.3 coincides with Theorem 1.1. This paper is organized as follows. Section 2 details the generalization of Sibson’s formula. In Section 3 we give a geometric interpretation of the formulas presented in Section 2 for point sets in the plane. Finally, Section 4 is on how the generalization of Sibson’s formula could be used for interpolation. Proofs are omitted in this abstract. EuroCG’24 21:4 Sibson’s formula for higher order Voronoi diagrams 2 Coordinates based on Voronoi diagrams In this section we present a generalization of Sibson’s formula that expresses a point using its neighbours of the Voronoi diagram of any given order. For this, we recall results from Aurenhammer [1] in Theorem 2.1 and Corollary 2.2, using a different notation. Let Fk+1 ( Pk ), to simplify F ( Pk ), be the set of cells of OVk+1 ( S ),1 ≤k≤n− 2, whose k nearest neighbours are the points of Pk⊂S in any order, and the ( k + 1)-th nearest neighbour is another point of S not in Pk . Let fi,j denote the union of cells of OVk+1 ( S ) whose k-th nearest neighbour is Qiand whose (k+ 1)-th nearest neighbour is Qj. ▶Theorem 2.1. ([1]) If all cells in F(Pk)are bounded in OVk+1(S), then X j fi,j ∈F(Pk) σ(fi,j)Qi=X i fi,j ∈F(Pk) σ(fi,j)Qj Q1 Q6 Q3 Q4 f5,1∪f5,3 f2,1∪f2,3∪f2,6 Q5 Q2 Q1 Q6 Q3 Q4 Q5 f2,1∪f5,1 f2,3∪f5,3 f2,6 Q2 Figure 4 Illustrating Theorem 2.1 for F ( {Q2, Q5} )in OV3 ( S ), where S is a set of six points in R2 . In this case the equation reduces to σ ( f5,1∪f5,3 ) Q5 + σ ( f2,1∪f2,3∪f2,6 ) Q2 = σ ( f2,1∪f5,1 ) Q1 + σ(f2,3∪f5,3)Q3+σ(f2,6)Q6. Left: cells grouped according to its k-nearest neighbour. Right: cells grouped according to its (k+ 1)-nearest neighbour. Note that, the subdivisions induced by Vk−1 ( S )at the interior of f ( Pk )correspond to grouping the cells of F ( Pk )in OVk+1 ( S )that have the same k -nearest neighbour. Also, the subdivisions induced by Vk+1 ( S )at the interior of f ( Pk )correspond to grouping the cells of F(Pk)in OVk+1(S)that have the same (k+ 1)-nearest neighbour. See Figure 4. By these observations, Theorem 2.1 can be stated as follows. ▶Corollary 2.2. ([1]) Let 2≤k≤n−2and let f(Pk)be a bounded cell of Vk(S). Then, X f(Pk−1)∈Vk−1(S) Qi∈Pk\Pk−1 σ(f(Pk−1)∩f(Pk))Qi=X f(Pk+1)∈Vk+1(S) Qj∈Pk+1\Pk σ(f(Pk+1)∩f(Pk))Qj Note that, by the relation between the Voronoi diagrams and the ordered Voronoi diagrams, Rk ( ℓ )is the set of cells of OVk+1 ( S )that have Qℓ as one of their k nearest neighbours from S , i.e., Rk ( ℓ ) = ∪Qℓ∈PkF ( Pk ). Based on this observation and Theorem 2.1 we can prove the following result. ▶Theorem 2.3. If Rk(ℓ)is a bounded region, then Qℓ=X i fi,j ∈Rk(ℓ) σ(fi,j) σ(Rk(ℓ))Qj. M. Claverol, A. de las Heras-Parrilla, C. Huemer and D. Lara 21:5 ▶Corollary 2.4. Let 1≤k≤n−2and let Rk(ℓ)be a bounded region. Then, Qℓ=X f(Pk)∈Rk(ℓ)X f(Pk+1)∈Vk+1(S) Qj∈Pk+1\Pk σ(f(Pk+1)∩f(Pk)) σ(Rk(ℓ)) Qj 3 A geometric interpretation In the following we examine the generalization of Sibson’s theorem to higher order Voronoi diagrams from Corollary 2.2 in more detail for cells f ( Pk )of Vk ( S ), when S is a point set in R2 . Divide both sides of the equation given in Corollary 2.2 by σ ( f ( Pk )); then, each side of the equation describes a point Hthat is a convex combination of points from S. We have H=X f(Pk−1)∈Vk−1(S) Qi∈Pk\Pk−1 σ(f(Pk−1)∩f(Pk)) σ(f(Pk)) Qi=X f(Pk+1)∈Vk+1(S) Qj∈Pk+1\Pk σ(f(Pk+1)∩f(Pk)) σ(f(Pk)) Qj(2) What can we say about this point H? Let f ( Pk )be an r -gon. Then S contains r points Q1, . . . , Qr , such that each edge of the r -gon lies on a perpendicular bisector between two of these r points, and each vertex, Cijℓ , of f ( Pk )is the center of a circle passing through three of them, Qi , Qj , and Qℓ ; see e.g. [ 3 ]. We denote with ∆( ABC )the triangle with vertices A , B , and C , and with □ ( ABCD ) the quadrilateral with vertices A, B, C and D, in cyclic order. Let us consider the case when f ( Pk )is a quadrilateral cell of Vk ( S )with vertices C123 , C124, C134 , and C234 , in cyclic order along the boundary of the quadrilateral cell f ( Pk ) = □ ( C123C124C134C234 ). One of the diagonals C123C134 and C124C234 is an edge of Vk−1 ( S ) and the other one of Vk+1 ( S ). Figure 5 shows an example. We refer to [ 3 , 7 ] for a more detailed discussion on the structure of cells of Vk(S). Corollary 2.2 states in this case that H=Q1·σ(∆(C123C134C234)) σ(□(C123C124C134C234)) +Q3·σ(∆(C123C124C134)) σ(□(C123C124C134C234)) =Q2·σ(∆(C124C134C234)) σ(□(C123C124C134C234)) +Q4·σ(∆(C124C234C123)) σ(□(C123C124C134C234)) (3) It follows that H is the intersection point of diagonals Q1Q3 and Q2Q4 of □ ( Q1Q2Q3Q4 ) . This implies that given a quadrilateral cell □ ( C123C124C134C234 )of Vk ( S ), the four corresponding points of S also form a convex quadrilateral, □ ( Q1Q2Q3Q4 ) . Moreover, we can show that areas of triangles with vertices in □ ( C123C124C134C234 )are proportional to areas of triangles with vertices in □(Q1Q2Q3Q4), also see [8, 12]. Let us then consider the case when f ( Pk )is a cell of Vk ( S )with more than four sides. Equation (2) gives a point H that can be expressed in two ways as convex combination of points of S . Let us look at a pentagonal cell f ( Pk ) = ( C123C134C145C245C125 )of Vk ( S ); See Figure 6. For r > 5the situation is similar. Corollary 2.2 here gives H=Q1·σ(□(C123C125C145C134)) σ( (C123C134C145C245C125)) +Q5·σ(∆(C125C245C145)) σ( (C123C134C145C245C125)) =Q2·σ(□(C245C125C123C234)) σ( (C123C134C145C245C125)) +Q4·σ(∆(C245C234C1345C145)) σ( (C123C134C145C245C125)) +Q3·σ(∆(C123C234C134)) σ( (C123C134C145C245C125)) EuroCG’24 21:6 Sibson’s formula for higher order Voronoi diagrams Q1 Q3 Q2 Q4 b23 b12 b34 b14 H C124 C123 C234 C134 Figure 5 The quadrilateral cell f ( Pk ) = □ ( C123C124C134C234 )of Vk ( S )is obtained by perpendicular bisector construction from {Q1, Q2, Q3, Q4} ⊂ S. Point H given by Equation (2) is the intersection point of diagonals Q1Q3and Q2Q4.Triangles with same color have proportional area. We get that H lies on the segment Q1Q5 and inside the triangle ∆( Q2Q3Q4 ). Furthermore, H divides the segment Q1Q5 in the same proportion as the edge C125C145 divides the pentagon ( C123C134C145C245C125 )into the the quadrilateral □ ( C125C145C134C123 )and the triangle ∆( C125C145C245 ). And H divides triangle ∆( Q1Q2Q3 )in the same proportion into triangles ∆( Q3HQ4 ),∆( Q2HQ3 ), and ∆( Q2HQ4 )as C234 divides ( C123C134C145C245C125 ) into □(C245C125C123C234),□(C245C234C134C145)and ∆(C134C234C123). 4 Towards higher order natural neighbour interpolation Sibson’s theorem (Theorem 1.1) gave rise to the natural neighbour interpolation method. Given a set of points S and known function values G ( Qj )for Qj∈S\{Qℓ} , the function value G ( Qℓ )of a point Qℓ is interpolated by G ( Qℓ ) = PjcjG ( Qj ), where the sum is over the natural neighbours Qj of Qℓ in V1 ( S ). The local coordinates cj are given by Theorem 1.1. Note that they satisfy Pjcj = 1 and cj≥ 0for all j . Then, Sibson’s natural neighour interpolation is given by G(Qℓ) = X j=ℓ σ(f({Qℓ, Qj})∩f({Qℓ})) σ(f({Qℓ})) G(Qj).(4) The generalization of Sibson’s formula given in Theorem 2.3 suggests to approximate the function value G ( Qi )by using the natural neighbours of higher order Voronoi diagrams. By using the region Rk(ℓ)for k > 1, we can estimate the function value of a point Qℓas G(Qℓ) = X i fi,j ∈Rk(ℓ) σ(fi,j) σ(Rk(ℓ))G(Qj).(5) Note that R1(ℓ) = f({Qℓ})in V1(S), and for k= 1 Equations (4) and (5) coincide. A better estimation can be obtained by using Theorem 2.3 in a combination of different values of k. We explore this for the 1-dimensional case. M. Claverol, A. de las Heras-Parrilla, C. Huemer and D. Lara 21:7 Q2 Q3 Q5 Q4 V1(S) V2(S) V3(S) OV3(S) F(P2) P2={Q1,Q5} Q1 C245 H C234 Q1 C245 C125 C123 C134 C145 Q2 Q3 Q5 Q4 Figure 6 (Left) OV3 ( S )for a set of five points S = {Q1, Q2, Q3, Q4, Q5}. For P2 = {Q1, Q5} , the grey region F ( P2 )of OV3 ( S )is the pentagonal cell f ( P2 )of V2 ( S ) . f ( P2 )is divided by an edge of V1 ( S )and is also divided by three edges of V3 ( S ). (Right) The point H lies on the segment Q1Q5 and inside the triangle ∆( Q2, Q3, Q4 ). Triangle areas of ∆( Q2HQ3 ) , ∆( Q3HQ4 )and ∆( Q2HQ4 ) are proportional to the areas of the three colored regions inside f ( P2 ), green, yellow, and pink, respectively. The lengths of segments HQ1 and HQ5 are proportional to the areas σ ( f ( P2 ) ∩f ( {Q1} )) and σ(f(P2)∩f({Q5})), respectively. Theorem 2.3, respectively Corollary 2.4, for dimension 1 reduces to the following statement. ▶Property 4.1. Let S={x0, x1,...x2ℓ}with x0< x1< . . . < x2ℓbe real numbers. Then, xℓ=1 x2ℓ−x0 ℓ−1 X i=0 xi(xℓ+1+i−xℓ+i)!+ 2ℓ X i=ℓ+1 xi(xi−ℓ−xi−ℓ−1)!!.(6) ▶ Remark. Property 4.1 has actually a more general statement. The assumption x0< x1< . . . < x2ℓis not needed. We denote points Qi of S as xi and their function values G ( Qi )as yi . When k = 1 we have Sibson’s classical nearest neighbour interpolation, which for dimension d = 1 is piecewise linear interpolation. Let x0, x1, . . . , x5 be six points on the real line in that order. And let x2< x < x3 be a query point whose function value G ( x )we want to interpolate. To avoid degenerate cases where bisectors between points coincide, we also assume that all midpoints ( xi + xj ) / 2with xi, xj∈ {S∪{x}} are different. Sibson’s classical formula, Equation (4), uses the two neighbours x2and x3of x, and gives the interpolation G1(x) = 1 x3−x2 (y2(x3−x) + y3(x−x2)),(7) i.e. point ( x, G1 ( x )) lies on the line segment connecting points ( x2, y2 )and ( x3, y3 ) . This can also be deduced from Property 4.1. Combining Equation (5) for k = 1 and k = 2, we obtain EuroCG’24 21:8 Sibson’s formula for higher order Voronoi diagrams G2(x) = 1 x4−x1+x3−x2 (y1(x3−x) + y2(x4−x) + y3(x−x1) + y4(x−x2)).(8) In the same way, combining Equation (5) for k= 1,k= 2, and k= 3, we obtain G3(x) = 1 x5−x0+x4−x1+x3−x2 (y0(x3−x) + y1(x4−x) + y2(x5−x) +y3(x−x0) + y4(x−x1) + y5(x−x2)).(9) Figure 7 shows an example of the interpolation formulas given in Equations (7), (8), and (9). 11 22 33 44 55 66 77 88 -2-2 -1-1 11 22 33 44 55 66 00 yy00 yy11 yy22 yy33 yy44 yy55 Figure 7 The generalized Sibson interpolation in R1 . In green: Sibson’s original interpolation, Equation (7), used only R1 ( x ). The blue segment shows the interpolation using R1 ( x )and R2 ( x ), given by Equation (8). Four points are used. The red segment shows the interpolation using R1 ( x ) , R2(x),and R3(x), given by Equation (9). Six points are used. We conclude with some comments on the proposed interpolation formulas. First, they appear in a natural way from the generalization of Sibson’s formula. This already makes it worth to study such generalized interpolation formulas. In Equations (8) and (9), the coefficients cj in Gi ( x ) = Pjcjyj , i = 2 , 3, satisfy Pjcj = 1 and cj≥ 0for every cj. We also mention that it can not be guaranteed that Gi ( x )coincides with Gi ( x2 )or with Gi ( x3 ), when x coincides with one of the endpoints of the interval, x2 or x3 , respectively. Though, we observe that in this case, the point farthest away from x on one side, drops from being used in the interpolation formula. This also holds for the classical case k= 1. Finally, we expect that the generalized interpolation formulas can have applications. For instance, when the used values for the interpolation are obtained by measurements and measurement inaccuracy can not be ruled out. Then reliability might be improved by using nearest neighbours from Vk(S)or by using Rk(x), instead of only V1(S). Acknowledgements This research has been supported by projects 2021SGR00266 and PID2019-104129GB-I00/ MCIN/ AEI/ 10.13039/501100011033. References 1 Franz Aurenhammer. Linear combinations from power domains. Geometriae dedicata, 28(1):45–52, 1988. M. Claverol, A. de las Heras-Parrilla, C. Huemer and D. Lara 21:9 2 Franz Aurenhammer and Otfried Schwarzkopf. A simple on-line randomized incremental algorithm for computing higher order voronoi diagrams. In Proceedings of the seventh annual symposium on Computational geometry, pages 142–151, 1991. 3 Mercè Claverol, Andrea de las Heras Parrilla, Clemens Huemer, and Alejandra MartínezMoraian. The edge labeling of higher order Voronoi diagrams. 2021. https://arxiv.org/ abs/2109.13002. 4 Herbert Edelsbrunner and Mabel Iglesias-Ham. Multiple covers with balls i: Inclusion– exclusion. Computational Geometry, 68:119–133, 2018. 5 Gerald Farin. Surfaces over Dirichlet tessellations. Computer aided geometric design, 7(1-4):281–292, 1990. 6 Gareth A Jones. Geometric and asymptotic properties of Brillouin zones in lattices. Bulletin of the London Mathematical Society, 16(3):241–263, 1984. 7 Der-Tsai Lee. On k-nearest neighbor Voronoi diagrams in the plane. IEEE Transactions on Computers, C-31(6):478–487, 1982. 8 Maria Flavia Mammana and Biagio Micale. Quadrilaterals of triangle centres. The Mathematical Gazette, 92:466–475, 2008. 9 Atsuyuki Okabe, Barry Boots, and Kokichi Sugihara. Nearest neighbourhood operations with generalized Voronoi diagrams: a review. International Journal of Geographical Information Systems, 8(1):43–71, 1994. 10 Atsuyuki Okabe, Barry Boots, Kokichi Sugihara, and Sung Nok Chiu. Spatial tessellations: concepts and applications of Voronoi diagrams. 2009. 11 Bruce R. Piper. Properties of local coordinates based on Dirichlet tessellations. In Geometric modelling, pages 227–239. Springer, 1993. 12 Olga Radko and Emmanuel Tsukerman. The perpendicular bisector construction, the isotopic point, and the Simson line of a quadrilateral. Forum Geometricorum, 12:161–189, 2012. 13 Robin Sibson. A vector identity for the Dirichlet tessellation. In Mathematical Proceedings of the Cambridge Philosophical Society, volume 87, pages 151–155. Cambridge University Press, 1980. 14 Robin Sibson. A brief description of natural neighbour interpolation. Interpreting multivariate data, pages 21–36, 1981. 15 Kokichi Sugihara. Surface interpolation based on new local coordinates. Computer-Aided Design, 31(1):51–58, 1999. 16 G Fejes Tóth. Multiple packing and covering of the plane with circles. Acta Math. Acad. Sci. Hungar, 27(1-2):135–140, 1976. 17 JJP Veerman, Mauricio M Peixoto, André C Rocha, and Scott Sutherland. On Brillouin zones. Communications in Mathematical Physics, 212:725–744, 2000. EuroCG’24