Creation of Dihedral Escher-like Tilings Based on As-Rigid-As-Possible Deformation

2.1 Basis

Some of the shapes of tiles in a tiling may be congruent to each other. The smallest set of non-congruent tile shapes that make up a tiling is called the set of prototiles. A tiling consisting of k prototiles is called k-hedral. In particular, the cases of \(k=1\) and \(k=2\) are referred to as monohedral and dihedral, respectively.

Consider a repeating pattern in the plane. If a rigid motion of a figure still matches itself, then that motion is called a symmetry of the figure. For two congruent tiles A and B in a tiling, if a rigid motion that maps one onto the other is a symmetry of the tiling, A and B are called transitively equivalent. Transitive equivalence is an equivalence relation that partitions the tiles into equivalence classes called transitivity classes. When a tiling has k transitivity classes, it is called k-isohedral. In particular, the case \(k=1\) is referred to as isohedral.

2.2 Isohedral Tiling

Among monohedral tilings, the class of isohedral tilings is the simplest but flexible enough to represent most monohedral tilings we see in everyday life. According to the definition, every tile in an isohedral tiling is adjacent to its neighbors in the same way. This allows an isohedral tiling to be constructed by repeatedly translating one or a set of several adjacent tiles. For a more detailed description of isohedral tilings, see Kaplan and Salesin [2000] and Kaplan [2009].

Isohedral tilings can be classified into 93 different types [Grünbaum and Shephard 1987], referred to as IH1, IH2,..., IH93, based on the adjacency relationship between tiles. If the adjacency relationship is not considered (only possible tile shapes are considered), then the 93 isohedral types can be reduced to nine [Schattschneider 2004]. For each isohedral type, possible tile shapes can be characterized by a template. Figure 1 shows the templates for six of the nine isohedral types, while the remaining three were not used in this study (as explained in Section 4.4) and are omitted. For each isohedral type, the template is represented as a polygon, and the vertices and edges can be moved and deformed under the specified constraints to generate tile shapes, which can always tile the plane as long as there is no self-intersection. In the nine isohedral types, the edges of the templates are classified into two types depending on how they can be deformed, as follows:

Fig. 1.

View Figure

Figure 2 shows an example of a tile shape generated from the template of IH4 (ignore the points in the template here). For each isohedral type, the adjacency relationship between tiles can be represented by arranging the template so that edges of the same type overlap each other; Figure 3 shows such an example for IH4 and IH3 (ignore the dotted lines and symbols \(U_1\) and \(U_2\) here).

Fig. 2.

View Figure

Fig. 3.

View Figure

2.3 Dihedral Tiling

A dihedral tiling is at least 2-isohedral because differently shaped tiles necessarily belong to different transitivity classes. If a dihedral tiling is 2-isohedral, then, for each of the two prototiles, every tile is adjacent to its neighbors in the same way. In this sense, among dihedral tilings, 2-isohedral tilings have the simplest structure.

One straightforward way of generating dihedral tilings is to split every tile of an isohedral tiling into two shapes in the same way by adding a “splitting path” [Kaplan and Salesin 2004]. The resulting tiling, which is called a “split isohedral tiling,” is necessarily a dihedral tiling with the 2-isohedral property. A template of split isohedral tiling can be constructed from the template of any isohedral type by adding a line that starts and ends on the boundary of the template, where the added line can be freely deformed as long as no intersections occur. Figure 3 shows examples of the split isohedral tiling template based on IH4 and IH3.

The split isohedral tilings cover all dihedral tilings with the 2-isohedral property if the two prototiles occur in equal amounts. However, there also exist 2-isohedral dihedral tilings with different relative amounts of the two prototiles, which cannot be generated by the split isohedral method [Kaplan and Salesin 2004].

Among dihedral tilings with the 2-isohedral property, there are those with a special property where every tile is adjacent only to tiles of another shape, i.e., every A tile is completely surrounded by B tiles (and vice versa). This property allows A and B tiles to be painted in two different colors so that no two tiles of the same color are adjacent to each other. The class of dihedral tilings with this property is called Heaven and Hell patterns [Kaplan and Salesin 2004], named after Escher’s well-known dihedral tiling consisting of white angels and black devils.

Dress [1986] classified the Heaven and Hell patterns into 37 types, 29 of which can be represented as split isohedral tilings. These 29 types can be further summarized into twelve types by considering symmetric features as special cases of asymmetric features [Kaplan and Salesin 2004]. The templates of three of these types are constructed from the templates of IH1, IH2, and IH3 shown in Figure 1 by adding splitting paths. Figure 3 shows the template of Heaven and Hell patterns based on IH3, where the splitting path is placed so that the paired j edges are completely separated. Therefore, both ends of the splitting path are uniquely determined.

3.1 Formulations

Let the goal shape \(\tilde{W}\) be given as a triangular mesh. The number of points on the boundary of \(\tilde{W}\) is denoted as n, where these points should be arranged at approximately equal intervals and are numbered clockwise from 1 to n. The number of inner points of \(\tilde{W}\) is denoted as \(n^{\prime }\), where these points are arbitrarily numbered from \(n+1\) to \(n+n^{\prime }\). We denote the coordinates of the ith point of \(\tilde{W}\) as \(\mathbf {w}_i = (x^w_i, y^w_i)^{\top }\), and the coordinates of all points are represented as \(\mathbf {\tilde{w}} = {(x^w_1, y^w_1, \dots , x^w_{n+n^{\prime }}, y^w_{n+n^{\prime }})}^{\top }\). A tile shape \(\tilde{U}\) is also represented as a triangular mesh that has the same structure as \(\tilde{W}\). We denote the coordinates of the ith point of \(\tilde{U}\) as \(\mathbf {u}_i = (x_i, y_i)^{\top }\), and the coordinates of all points are represented as \(\mathbf {\tilde{u}} = {(x_1, y_1, \dots , x_{n+n^{\prime }}, y_{n+n^{\prime }})}^{\top }\). Figure 4 shows examples of \(\tilde{W}\) and \(\tilde{U}\). The vectors \(\mathbf {\tilde{w}}\) and \(\mathbf {\tilde{u}}\) are sometimes represented by (1) \(\begin{equation} \mathbf {\tilde{w}} = \begin{pmatrix} \mathbf {w} \\ \mathbf {w^{\prime }} \\ \end{pmatrix} \ \ \mbox{and} \ \ \mathbf {\tilde{u}} = \begin{pmatrix} \mathbf {u} \\ \mathbf {u^{\prime }} \\ \end{pmatrix}, \end{equation}\) where the vectors \(\mathbf {w}\) (\(\mathbf {u}\)) and \(\mathbf {w^{\prime }}\) (\(\mathbf {u^{\prime }}\)) are the coordinates of the boundary and inner points, respectively. The polygons consisting of the boundaries of \(\tilde{W}\) and \(\tilde{U}\) are denoted as W and U, respectively.

Fig. 4.

View Figure

Possible tile shapes \(\tilde{U}\) are parameterized using templates of isohedral tilings, where only the boundary points are constrained by the template, while the inner points are freely moved. Figure 2 shows an example of the template of IH4 represented as an n-point polygon; exactly one point must be placed on every corner of the template, and zero or more points are placed on each edge of the template, but each pair of J edges must be assigned the same number of points. We denote the lengths (numbers of segments divided by points) of the edges as \(k_1, k_2\), ....¹ The n points on the template are numbered clockwise from 1 to n, starting from one of the vertices of the template. We call this the original numbering.

The n points on the template can move freely under the constraints imposed by the template. For example, Figure 2 shows a tile shape U obtained by deforming the template. As in this example, the constraint conditions imposed by any template can be expressed as a homogeneous system of linear equations of the form \(A \mathbf {u} = \mathbf {0}\) with a matrix A [Koizumi and Sugihara 2011; Nagata and Imahori 2020]. Therefore, possible values of \(\mathbf {u}\) can be parameterized as the general solutions of this system in the form of (2) \(\begin{equation} \mathbf {u} = B \mathbf {\xi }, \end{equation}\) where B is a matrix and \(\mathbf {\xi }\) is a parameter vector.

Let I be a set of indices for the isohedral types and \(K_i \ (i \in I)\) be defined as the set of all possible configurations of \((k_1, k_2, \dots)\) for the template of IHi. For example, for IH4, template configuration is specified by \((k_1, k_2, k_3, k_4)\), from which \(k_5\) is automatically determined. As the matrix B depends on \(i \in I\) and \(k \in K_i\), we denote it as \(B_{ik}\). In addition, we consider n different numbering schemes for the n points of the template. Let \(J = \lbrace 1, 2, \dots , n\rbrace\) be a set of indices that specify the n different numbering schemes; if \(j \in J\) is selected, then the points of the template are renumbered such that the first point in the original numbering is assigned as j.² Indices k and j specify the correspondence between the edges of the template and the boundary of the goal shape W. The matrix B also depends on \(j \in J\) and is denoted as \(B_{ikj}\), which is obtained from \(B_{ik}\) simply by cyclically shifting up the row elements.

To measure the similarity between the tile shape \(\tilde{U}\) and goal shape \(\tilde{W}\), a distance function based on the ARAP deformation technique [Sorkine and Alexa 2007] was developed. Let \(\mathcal {N}(i)\) be the set of points of \(\tilde{W}\) (\(\tilde{U}\)) connected to point i. For each point i of \(\tilde{W}\) (\(\tilde{U}\)), cell \(\mathcal {W}_i\) (\(\mathcal {U}_i\)) is defined by the edges between point i and its neighboring points \(j \in \mathcal {N}(i)\). The distance function is then defined by (3) \(\begin{equation} d^2_{IR}(\tilde{U},\tilde{W}) = \frac{1}{2} \sum _{i=1}^{n+n^{\prime }} \alpha _i \sum _{j \in \mathcal {N}(i)} {\left\Vert ({\mathbf {u}}_i - {\mathbf {u}}_j) - R_i ({\mathbf {w}}_i - {\mathbf {w}}_j) \right\Vert }^2, \end{equation}\) where \(R_i\) denotes a rotation matrix. For each point i, \(R_i\) is determined to minimize (4) \(\begin{equation} \sum _{j \in \mathcal {N}(i)} {\left\Vert ({\mathbf {u}}_i - {\mathbf {u}}_j) - R_i ({\mathbf {w}}_i - {\mathbf {w}}_j) \right\Vert }^2, \end{equation}\) which measures the degree of local deformation between \(\mathcal {U}_i\) and \(\mathcal {W}_i\). Therefore, this distance function is defined as the weighted sum of the degrees of local deformation over all cells, where \(\alpha _i\) is the weight assigned to the cells. The value of \(\alpha _{i}\) was set to 1.0 (\(1 \le i \le n\)) and 0.5 (\(n+1 \le i \le n+n^{\prime }\)), and we used the same setting in the proposed method. This distance function is expected to be less sensitive to physically plausible deformations with well-preserved local structures between \(\tilde{U}\) and \(\tilde{W}\). This allows for the generation of satisfactory tilings even when large deformations of the goal shape are indispensable to form possible tile shapes.

For a given distance function, the exhaustive search of templates (EST) refers to minimizing the distance function between \(\tilde{U}\) and \(\tilde{W}\) under the constraint \({\mathbf {u}} = B_{ikj} {\mathbf {\xi }}\) for all possible configurations of \(i \in I\), \(k \in K_i\), and \(j \in J\). After performing the EST, the top \(n_\mathrm{top}\) (e.g., 10) tile shapes in terms of distance value are displayed to the user.

3.2 Optimization

We outline how to minimize the distance function \(d^2_{IR}(\tilde{U},\tilde{W})\) under the constraint \(\mathbf {u} = B_{ikj} \mathbf {\xi }\). Although this optimization problem cannot be solved analytically because the matrices \(R_i \ (1 \le i \le n+n^{\prime })\) depend on \(\mathbf {\tilde{u}}\), a near-optimal solution can be obtained through a two-step iterative procedure described later.

If the rotation matrices \(R_i\) are fixed, then \(d^2_{IR}(\tilde{U},\tilde{W})\) can be written as a quadratic form of \(\tilde{\mathbf {u}}\): (5) \(\begin{equation} d^2_{IR}(\tilde{U},\tilde{W}) = {\tilde{\mathbf {u}}}^{\top } \tilde{G} {\tilde{\mathbf {u}}} - 2 {\tilde{\mathbf {u}}}^{\top } \tilde{T} {\tilde{\mathbf {w}}} + {\tilde{\mathbf {w}}}^{\top } \tilde{G} {\tilde{\mathbf {w}}}, \end{equation}\) where \(\tilde{G}\) and \(\tilde{T}\) are matrices. The matrix \(\tilde{T}\) is constructed as follows, where \(\tilde{T}.block(i,j)\) is the \(2 \times 2\) block matrix of \(\tilde{T}\) beginning at location \((i,j)\).

The matrix \(\tilde{G}\) is obtained by replacing \(R_i\) with the identity matrix in the above procedure.

Since the inner points of \(\tilde{U}\) can move freely, \(d^2_{IR}(\tilde{U},\tilde{W})\) can be written as a function of \({\mathbf {u}}\) by determining their positions to minimize this function. Let T and \(T^{\prime }\) be the first 2n rows and the last \(2n^{\prime }\) rows of \(\tilde{T}\), respectively. Then, \(d^2_{IR}(\tilde{U},\tilde{W})\) can be written as a quadratic form of \({\mathbf {u}^{\prime }}\): (6) \(\begin{align} d^2_{IR}(\tilde{U},\tilde{W}) =&\, {\mathbf {u}}^{\top } G {\mathbf {u}} + 2 {\mathbf {u}}^{\top } G^{\prime \prime } {\mathbf {u}^{\prime }} + {\mathbf {u}^{\prime }}^{\top } G^{\prime } {\mathbf {u}^{\prime }} \nonumber \nonumber\\ & -2 ({\mathbf {u}}^{\top } T {\tilde{\mathbf {w}}} + {\mathbf {u}^{\prime }}^{\top } T^{\prime } {\tilde{\mathbf {w}}}) + {\tilde{\mathbf {w}}}^{\top } \tilde{G} {\tilde{\mathbf {w}}}, \end{align}\) where \(G \ (\in \mathbb {R}^{2n \times 2n})\) is the top-left corner of \(\tilde{G}\), \(G^{\prime \prime } \ (\in \mathbb {R}^{2n \times 2n^{\prime }})\) is the top-right corner of \(\tilde{G}\), and \(G^{\prime } \ (\in \mathbb {R}^{2n^{\prime } \times 2n^{\prime }})\) is the bottom-right corner of \(\tilde{G}\). The matrices \(\tilde{G}\), G, and \(G^{\prime }\) are symmetric. By solving \(\frac{\partial d^2_{IR}(\tilde{U},\tilde{W})}{\partial {\mathbf {u}^{\prime }}} = 0\), we have the positions of the inner points: (7) \(\begin{equation} {\mathbf {u}^{\prime }} = {G^{\prime }}^{-1} (T^{\prime } {\tilde{\mathbf {w}}} - {G^{\prime \prime }}^{\top } {\mathbf {u}}). \end{equation}\) By substituting Equation (7) into Equation (6), we finally obtain the following energy function with respect to \({\mathbf {u}}\): (8) \(\begin{align} E_{IR}({\mathbf {u}}) =&\, \min _{\mathbf {u}^{\prime }} d^2_{IR}(\tilde{U},\tilde{W}) \nonumber \nonumber\\ =&\,\, {\mathbf {u}}^{\top } G_I {\mathbf {u}} -2 {\tilde{\mathbf {w}}}^{\top } H^{\top } {\mathbf {u}} + {\tilde{\mathbf {w}}}^{\top } (\tilde{G} - {T^{\prime }}^{\top } {G^{\prime }}^{-1} {T^{\prime }}) {\tilde{\mathbf {w}}}, \end{align}\) where \(G_I = G - G^{\prime \prime } {G^{\prime }}^{-1} {G^{\prime \prime }}^{\top }\) and \(H = T - G^{\prime \prime } {G^{\prime }}^{-1} T^{\prime }\). Because this is a quadratic form of \(\mathbf {u}\), the minimization of this function under the constraint \({\mathbf {u}} = B_{ikj} {\mathbf {\xi }}\) is straightforward.

For some isohedral types (IH2, IH3, IH5, and IH6), however, the energy function \(E_{IR}({\mathbf {u}})\) should be rotation invariant because the parameterized tile shape U can only appear in a specific orientation. Let \(R_n(\theta) \ (\in \mathbb {R}^{2n \times 2n})\) be the matrix that rotates n points represented by \({\mathbf {u}}\) by angle \(\theta\). A rotation-invariant version of \(E_{IR}({\mathbf {u}})\) is defined and calculated as follows: (9) \(\begin{align} E_{IRP}({\mathbf {u}}) =&\, \min _{\theta } E_{IR}(R_n(\theta) {\mathbf {u}}) \nonumber \nonumber\\ =&\,\, {\mathbf {u}}^{\top } G_I {\mathbf {u}} - 2 \sqrt {{\mathbf {u}}^{\top } H \tilde{V} H^{\top } {\mathbf {u}}} + {\tilde{\mathbf {w}}}^{\top } (\tilde{G} - {T^{\prime }}^{\top } {G^{\prime }}^{-1} {T^{\prime }}) {\tilde{\mathbf {w}}}, \end{align}\) where \(\tilde{V} = \tilde{\mathbf {w}} \tilde{\mathbf {w}}^{\top } + \tilde{\mathbf {w_c}} \tilde{\mathbf {w_c}}^{\top }\) and \(\mathbf {\tilde{w_c}} = (y^w_1, -x^w_1, \ldots , y^w_{n+n^{\prime }}, -x^w_{n+n^{\prime }})^{\top }\). The minimization of \(E_{IRP}({\mathbf {u}})\) under the constraint \({\mathbf {u}} = B_{ikj} {\mathbf {\xi }}\) can be reduced to finding the maximum eigenvalue for a generalized eigenvalue problem given by \({B_{ikj}}^{\top } H \tilde{V} H^{\top } B_{ikj} {\mathbf {\xi }} = \lambda {B_{ikj}}^{\top } G_I B_{ikj} {\mathbf {\xi }}\).

The two-step iterative procedure is then constructed as follows:

In the above procedure, the matrix \(\tilde{T}\) is initialized by setting the rotation matrices \(R_i\) to the identity matrix, and the iteration is repeated until the value of the energy function converges.

Note that the energy function \(E_{IR}({\mathbf {u}})\) is rotation invariant in the sense that the rotation matrices \(R_i\) are optimized separately during the two-step iterative procedure. Therefore, the same result can be obtained using either energy function \(E_{IR}({\mathbf {u}})\) or \(E_{IRP}({\mathbf {u}})\). However, for some isohedral types (IH2, IH3, IH5, and IH6), the value of the energy function converges faster with the latter, especially when the global rotation angle \(\theta\) is large.

4.1 Parametrization of Tile Shapes

Given two goal shapes, corresponding tile shapes are defined for each. The two goal shapes and two tile shapes are represented in the manner described in Section 3.1. Let the two goal shapes represented as triangular meshes be denoted as \(\tilde{W}_1\) and \(\tilde{W}_2\). The other variables, constants, and symbols listed in Table 1 are also denoted in the same way.

For each isohedral type IHi, possible tile shapes \(U_1\) and \(U_2\) are parameterized by a template constructed by adding a splitting path to the template of IHi. For example, Figure 5(a) shows the template of split isohedral tilings constructed from the template of IH4 shown in Figure 2, where the templates of the tile shapes \(U_1\) and \(U_2\) are indicated by the same symbols. Recall that in the template of IHi, the lengths of the edges are denoted as \((k_1, k_2\), ...), which are abbreviated as k. In the split isohedral case, two parameters, s and \(k_0\), are added. The parameter s specifies the position of the end of the splitting path in the original numbering of the template of IHi (see Figure 2), while the position of another end is denoted as e, which is determined automatically. The parameter \(k_0\) specifies the length of the splitting path. The points of \(U_1\) (\(U_2\)) are indexed clockwise from 1 to \(n_1\) (\(n_2\)) starting from the point s (e). Let n denote the number of points on the perimeter of the template, which is determined by \(n = n_1 + n_2 - 2k_0\). The value of e is given by \(e = s+n_1-k_0\). Note that if \(n \lt e\), then n is subtracted from e, and the same applies to other variables thereafter as needed.

Fig. 5.

View Figure

Given that the positions of points on the template of IHi with the original numbering are parameterized by \(\mathbf {u} = B_{ik} \mathbf {\xi }\) and points on the splitting path (not including both ends) can be freely moved, possible tile shapes \(U_1\) and \(U_2\) for a given template configuration \((i,k_0,k,s)\) can be parameterized as follows: (10) \(\begin{equation} \mathbf {u}_1 = B_1 \mathbf {\xi } \ \ \mbox{and} \ \mathbf {u}_2 = B_2 \mathbf {\xi }, \end{equation}\) where the matrices \(B_1\) and \(B_2\) are given by

(11)

As in the Escherization case, we need to consider \(n_1\) (\(n_2\)) different numbering schemes for the points of \(U_1\) (\(U_2\)). Let \(J_1 = \lbrace 1, 2, \dots , n_1\rbrace\) be the set of indices that specify the \(n_1\) different numbering schemes; if \(j_1 \in J_1\) is selected, then the points of \(U_1\) are renumbered such that the first point (point s) is assigned as \(j_1\). For each value of \(j_1 \in J_1\), the matrix \(B^{(1)}_{ik_0ks}\) depends on \(j_1\) and is denoted as \(B^{(1)}_{ik_0ksj_1}\), which is obtained from \(B^{(1)}_{ik_0ks}\) simply by cyclically shifting up the row elements. For the tile shape \(U_2\), \(J_2\), and \(B^{(2)}_{ik_0ksj_2} \ (j_2 \in J_2)\) are defined in the same manner.

4.2 Energy Function

Recall that the energy function \(E_{IR}({\mathbf {u}})\) (Equation (8)) can be regarded as a metric of the similarity between the tile shape U (represented as a polygon) and goal shape \(\tilde{W}\) (represented as a mesh), where the positions of the inner points of the tile shape \(\tilde{U}\) (represented as a mesh) are determined automatically. For the parameterized tile shapes \(U_1\) and \(U_2\), we evaluate the similarity between \(U_1\) and \(\tilde{W}_1\) (\(U_2\) and \(\tilde{W}_2\)) using a similar function.

Recall that the energy function \(E_{IRP}({\mathbf {u}})\) (Equation (9)) is a rotation-invariant version of \(E_{IR}({\mathbf {u}})\) and this property is always required in the dihedral case because the parameterized tile shapes \(U_1\) and \(U_2\) cannot be rotated independently. In addition, scale invariance is also required because \(U_1\) and \(U_2\) cannot be scaled independently. Therefore, we need to modify the energy function \(E_{IR}({\mathbf {u}})\) so that it is scale and rotation invariant.

The scale- and rotation-invariant version of the energy function \(E_{IR}({\mathbf {u}})\) is defined and calculated as follows: (12) \(\begin{align} E({\mathbf {u}}) & = \min _{s, \theta } E_{IR}(sR_n(\theta) {\mathbf {u}}) \nonumber \nonumber\\ & = {\tilde{\mathbf {w}}}^{\top } (\tilde{G} - {T^{\prime }}^{\top } {G^{\prime }}^{-1} {T^{\prime }}) {\tilde{\mathbf {w}}} - \frac{{\mathbf {u}}^{\top } H \tilde{V} H^{\top } {\mathbf {u}}}{{\mathbf {u}}^{\top } G_I \mathbf {u}}. \end{align}\) Let the optimal values of s and \(\theta\) in Equation (12) be denoted as \(s^*\) and \(\theta ^*\). These values are given by (13) \(\begin{equation} s^* \cos \theta ^* = \frac{{\tilde{\mathbf {w}}}^{\top } H^{\top } {\mathbf {u}}}{{\mathbf {u}}^{\top } G_I \mathbf {u}} \quad \mbox{and} \quad s^* \sin \theta ^* = \frac{{\tilde{\mathbf {w_c}}}^{\top } H^{\top } {\mathbf {u}}}{{\mathbf {u}}^{\top } G_I \mathbf {u}}, \end{equation}\) from which \(s^*\), \(\cos \theta ^*\), and \(\sin \theta ^*\) are easily obtained. The derivation of Equations (12) and (13) is presented in Appendix A.

Let \(E_1({\mathbf {u}}_1)\) and \(E_2({\mathbf {u}}_2)\) be the energy functions that measure the similarity between \(U_1\) and \(\tilde{W}_1\) and between \(U_2\) and \(\tilde{W}_2\), respectively. That is, \(E_1({\mathbf {u}}_1)\) is defined by (14) \(\begin{equation} E_1({\mathbf {u}_1}) = {\tilde{\mathbf {w}}_1}^{\top } (\tilde{G}_1 - {T^{\prime }_1}^{\top } {G^{\prime }_1}^{-1} {T^{\prime }_1}) {\tilde{\mathbf {w}}_1} - \frac{{\mathbf {u}}_1^{\top } H_1 \tilde{V}_1 H_1^{\top } {\mathbf {u}}_1}{\mathbf {u}_1^{\top } G_{1I} \mathbf {u}_1}, \end{equation}\) and \(E_2({\mathbf {u}}_2)\) is defined in the same manner. We define an objective function for the dihedral Escherization problem as follows: (15) \(\begin{equation} E({\mathbf {u}_1}, {\mathbf {u}_2}) = \frac{n_1}{a_1} E_1({\mathbf {u}}_1) + \frac{n_2}{a_2} E_2({\mathbf {u}}_2), \end{equation}\) where \(a_1 = {\mathbf {w}_1}^{\top } G_{1I} \mathbf {w}_1\) and \(a_2 = {\mathbf {w}_2}^{\top } G_{2I} \mathbf {w}_2\). The coefficients \(\frac{1}{a_1}\) and \(\frac{1}{a_2}\) normalize the two energy functions such that \(\frac{1}{a_1} E_1({\mathbf {u}}_1)\) and \(\frac{1}{a_2} E_2({\mathbf {u}}_2)\) take values between 0 and 1 when the rotation matrices \(R_i\) are all set to the identity matrix. The coefficients \(n_1\) and \(n_2\) balance the magnitude of change in the two energy functions. Having \(n_1 \lt n_2\) without this adjustment is likely to minimize \(E({\mathbf {u}_1}, {\mathbf {u}_2})\) in favor of \(\frac{1}{a_1} E_1({\mathbf {u}}_1)\) over \(\frac{1}{a_2} E_2({\mathbf {u}}_2)\) because local changes at the boundary of two tile shapes affect \(\frac{1}{a_1} E_1({\mathbf {u}}_1)\) more significantly than \(\frac{1}{a_2} E_2({\mathbf {u}}_2)\).

After minimizing the energy function \(E({\mathbf {u}_1}, {\mathbf {u}_2})\) in the manner described later, tile shapes \(U_1\) and \(U_2\) are optimized to be similar to \(\tilde{W}_1\) and \(\tilde{W}_2\), respectively, in the sense of the ARAP deformation energy. Depending on the template configuration, various tile shapes are obtained. Figure 6 shows three pairs of tile shapes \(U_1\) and \(U_2\); the two on the left, which are included in the top-10 solutions of the EST (see Section 4.4), have small energy function values while the one on the right has a much larger value than them.

Fig. 6.

View Figure

When their mesh representations \(\tilde{U}_1\) and \(\tilde{U}_2\) are required, we need to compute \(\mathbf {u^{\prime }}_1\) and \(\mathbf {u^{\prime }}_2\) from \({\mathbf {u}}_1\) and \({\mathbf {u}}_2\), which can be achieved as follows (see Figure 7 for illustration). First, \(U_1\) and \(U_2\) are scaled and rotated according to \((s_1^*, \theta _1^*)\) and \((s_2^*, \theta _2^*)\). Then, the coordinates of the inner points are computed from the resulting tile shapes using Equation (7). Finally, the obtained tile meshes are re-scaled and re-rotated according to \((\frac{1}{s_1^*}, -\theta _1^*)\) and (\(\frac{1}{s_2^*}, -\theta _2^*\)) to obtain \(\tilde{U}_1\) and \(\tilde{U}_2\).

Fig. 7.

View Figure

Some readers may think that the energy function \(E_1({\mathbf {u}}_1)\) (and \(E_2({\mathbf {u}}_2)\)) need not to be rotation invariant because the rotation matrices \(R_i\) are optimized. As in the Escherization problem (see the end of Section 3.2), rotation invariance is preferred because this reduces the number of iterations required in the two-step iterative procedure. Furthermore, we observed that without this property, the solution sometimes got stuck in a bad local optimum during the two-step iterative procedure, although this was a rare occurrence.

As shown in Figures 6 and 7, good tile shapes \(U_1\) and \(U_2\) tend to be optimized so that the lengths of all edges of \(U_1\) and \(U_2\) are roughly the same (even if \(\tilde{W}_1\) and \(\tilde{W}_2\) are given in arbitrary size) because some of the edges of \(U_1\) and \(U_2\) are shared at the boundary. Therefore, on condition that the goal shape \(\tilde{W}_2\) is scaled so that the average length of the edges of \(W_2\) is equal to that of \(W_1\), we would have roughly the same result even if scale invariance is dropped from the energy functions \(E({\mathbf {u}_1})\) and \(E({\mathbf {u}_2})\).³ Nevertheless, scale invariance is preferred because there also exist the tile boundary where \(U_1\) (\(U_2\)) is adjacent to itself.

4.3 Optimization

For each configuration \((i,k_0,k,s,j_1,j_2)\), we want to solve the following optimization problem: (16) \(\begin{align} \begin{split}\mbox{min} & \quad E({\mathbf {u}_1}, {\mathbf {u}_2}) \\ \mbox{s.t.} & \quad \mathbf {u}_1 = B^{(1)}_{ik_0ksj_1}\mathbf {\xi }, \ \mathbf {u}_2 = B^{(2)}_{ik_0ksj_2}\mathbf {\xi }.\end{split} \end{align}\) As in the Escherization problem (see Section 3.2), we use the two-step iterative procedure to solve this optimization problem.

Let the matrices \(\tilde{T}_1\) and \(\tilde{T}_2\) be initialized by setting the rotation matrices \(R_i\) to the identity matrix for both \(\tilde{W}_1\) and \(\tilde{W}_2\). In the first step, the matrices \(\tilde{T}_1\) and \(\tilde{T}_2\) (and in turn, \(T^{\prime }_1\), \(T^{\prime }_2\), \(H_1\), and \(H_2\)) are fixed in the energy function \(E({\mathbf {u}_1}, {\mathbf {u}_2})\). By substituting the constraint equations into the objective function, the optimization problem to be solved is equivalent to the following problem: (17) \(\begin{equation} \text{argmax}_{\mathbf {\xi }} \left\lbrace \frac{n_1}{a_1} \frac{{\mathbf {\xi }}^{\top } B_1^{\top } H_1 \tilde{V}_1 H_1^{\top } B_1 \mathbf {\xi } }{{\mathbf {\xi }}^{\top } B_1^{\top } G_{1I} B_1 \mathbf {\xi }} + \frac{n_2}{a_2} \frac{{\mathbf {\xi }}^{\top } {B_2}^{\top } H_2 \tilde{V}_2 H_2^{\top } B_2 \mathbf {\xi } }{{\mathbf {\xi }}^{\top } B_2^{\top } G_{2I} B_2 \mathbf {\xi }} \right\rbrace , \end{equation}\) where \(B_1 = B^{(1)}_{ik_0ksj_1}\) and \(B_2 = B^{(2)}_{ik_0ksj_2}\). By solving this optimization problem, we obtain the tile shapes \(U_1\) and \(U_2\). The solution of this optimization problem will be explained later.

In the second step, the rotation matrices \(R_i\) are determined to minimize Equation (4) for the mesh representations of the tile shapes \(U_1\) and \(U_2\) obtained in the first step. Note that the mesh representations used here are those obtained after scaling and rotating \(U_1\) and \(U_2\) according to \((s_1^*, \theta _1^*)\) and (\(s_2^*, \theta _2^*)\) in the procedure shown in Figure 7. Then, the matrices \(\tilde{T}_1\) and \(\tilde{T}_2\) are updated.

In our setting, the first and second steps are repeated three times. Tile shapes \(U_1\) and \(U_2\) are finally determined at the first step of the final iteration. If \(U_1\) and \(U_2\) contain self-intersections, then they are discarded. After the matrices \(\tilde{T}_1\) and \(\tilde{T}_2\) are updated at the second step of the final iteration, the value of the energy function \(E({\mathbf {u}_1}, {\mathbf {u}_2})\) is re-computed to evaluate the obtained tile shapes.

Now, we describe how to solve Equation (17). For symmetric matrices B and D and positive definite matrices W and V, a function (18) \(\begin{equation} \frac{{\mathbf {x}}^{\top } B \mathbf {x}}{{\mathbf {x}}^{\top } W \mathbf {x}} + \frac{{\mathbf {x}}^{\top } D \mathbf {x}}{{\mathbf {x}}^{\top } V \mathbf {x}} \end{equation}\) is called the sum of a generalized Rayleigh quotient, and the function to be maximized in Equation (17) is of this type. In general, this type of function has many local optima [Nguyen et al. 2016; Zhang 2013], and it is difficult to maximize it.

To maximize the objective function in Equation (17), we employed the conjugate gradient method with the Fletcher–Reeves rule. Utilizing the structure of the objective function, we constructed an efficient algorithm to perform the line search performed in the process of the conjugate gradient method; if some matrix computations are completed in \(O(n^2)\) time at the beginning of the linear search, then each of the function values of the points to be sampled can be computed in a constant time. Because the objective function is scale invariant, \(\mathbf {\xi }\) is normalized each time the line search is completed to avoid the divergence of \(\mathbf {\xi }\).

Because the sum of a generalized Rayleigh quotient generally has many local maxima, a conjugate gradient method does not yield a global optimal solution. However, we infer that the objective function in Equation (17) appearing in the two-step iterative procedure is convex for almost all template configurations. This is because we have confirmed that the conjugate gradient method always converges to the same point, even when starting from multiple randomly generated initial points. This is not theoretically guaranteed, but at least we can expect that this objective function will not have many local optima because \(\tilde{V}_1\) and \(\tilde{V}_2\) are matrices of rank 2 (see Equation (9)).

4.4 EST and HH-EST

For the dihedral Escherization problem, we define the EST as solving the optimization problem (16) for all possible configurations of \((i,k_0,k,s,j_1,j_2)\). If the template configuration is restricted to Heaven and Hell patterns, then we call this HH-EST.

In the EST of the Escherization problem [Nagata and Imahori 2021], nine isohedral types were considered, but most of the top-ranked solutions were generated from the three isohedral types IH4, IH5, and IH6 (see Figure 1). This is because the templates of these isohedral types, especially IH4, can be deformed with the highest degrees of freedom. Therefore, we consider only these isohedral types for the EST of the dihedral Escherization problem.

However, performing the EST is computationally impractical because the number of all possible configurations of \((i,k_0,k,s,j_1,j_2)\) is huge. If we consider IH4, then k is the abbreviation of \((k_1, k_2, k_3, k_4)\), and, therefore, the order of the number of all possible template configurations is \(O(n^8)\), where n is defined here as \(n_1+n_2\). In the same way, this order is \(O(n^7)\) for IH5 and IH6. For example, if \(n_1=n_2=60\), then the number of template configurations that must be examined in the EST is approximately \(1.7 \times 10^{11}\). We estimate that it would take approximately 4 years to perform the EST using the modern PC used in our experiments.

We developed an alternative approach to obtain results in a reasonable computation time. The basic idea is to solve the optimization problem (16) only for promising template configurations, which are selected from all possible configurations based on a low-cost evaluation function. This process is described in the next section.

The overall procedure of the EST is summarized as follows.

One thing needs to be added to this procedure. If one of the goal shapes is reflected, then different results are obtained. Therefore, reflection of \(\tilde{W}_2\) is also considered.

Of course, there is no guarantee that the top-\(n_\mathrm{top}\) solutions obtained with this approach match the true top-\(n_\mathrm{top}\) solutions. This concern is verified in Section 6.4.

When Heaven and Hell patterns are considered, we use only the three isohedral types IH1, IH2, and IH3 among the 12 possible (see Section 2.3) because the templates of split isohedral tilings constructed from other isohedral types have lower degrees of freedom in possible deformations. For the HH-EST, the value of s is fixed so that the end of the splitting path is at the specified position (see Figure 3). Because \(k = (k_1, k_2)\) for IH1, IH2, and IH3, and the value of s is fixed, the order of the number of all possible template configurations is \(O(n^5)\). For example, if \(n_1=n_2=60\), then the number of template configurations that must be examined is approximately \(3.0 \times 10^{8}\). Although this number is much smaller than in the case of the EST, it is still too large to perform the HH-EST in a realistic computation time. As in the case of the EST, we developed a method to select promising template configurations. Unlike the case of the EST, \(n_1\) must be equal to \(n_2\), and reflection of \(\tilde{W}_2\) is not necessary.

6.1 Experimental Settings

We developed two algorithms to perform the EST and HH-EST. Note that these algorithms solve the optimization problem (16) only for the promising \(n_\mathrm{conf}\) template configurations selected through the preliminary estimate described in Section 5. In our setting, \(n_\mathrm{conf}\) was set to \(10^6\). These algorithms were applied to all 15 pairs that could be composed of the six goal figures shown in Figure 8. In this figure, the corresponding mesh representations are also displayed, where the numbers of points on the boundaries are all 60. After performing the EST or HH-EST, the top-10 solutions (tile shapes and tilings) were displayed, where the textures of the tile figures were generated by deforming the textures of the goal figures according to the mesh deformation. The results presented here are those that the authors thought were the best among the top-10 solutions, unless otherwise stated.

Fig. 8.

View Figure

The algorithms were implemented in C++ with the Eigen library on Ubuntu 20.04 and executed on PCs with Intel Core i7-12700 CPU (12 Cores/20 Threads/2.1 GHz) and 32 GB of DDR4-3200 memory. The program can be executed in parallel using OpenMP. The source code is available at https://github.com/nagata-yuichi/Dihedral-Escher-Tiling.

6.2 Main Results

First, the results of the EST are described. Figure 9 shows the tile figures obtained by the EST for the 15 pairs of goal figures. The comparison here focuses on the differences between the tile figures and goal figures rather than the tiling results, so individual tile figures are shown scaled and rotated to match the goal figures. Although the obtained tile figures are deformed from the goal figures to a greater or lesser extent, they retain the local structures of the goal figures even when the deformations are somewhat large. As a result, all tile figures are generally obtained by natural deformations of the goal figures. Figure 10 shows the tiling results created from the tile figures presented in Figure 9, where congruent tile figures are painted in two different tones so that tiles of the same tone are adjacent to each other as little as possible. The resulting tilings are generally satisfactory because the tile figures retain the features of the goal figures well. More detailed results, including all of the top-10 solutions, are provided in the online supplementary file 1. In addition, results for goal shapes other than the six used in the experiment are also included in the supplementary file 2 to verify the robustness of the EST.

Fig. 9.

View Figure

Fig. 10.

View Figure

View Figure

The artistry of the tiling may also depend on how the tile figures are arranged. For example, the displayed tiling for the goal figures Pegasus and camel (fifth row of Figure 10) may seem uninteresting because the two types of tile figures are so tidily arranged, and it might be more interesting to have two types of tile figures arranged in an intricate manner. For example, interesting tilings are often obtained when two adjacent A tiles are completely surrounded by B tiles, e.g., see the tiling for the goal figures camel and squirrel (third row of Figure 10). Figure 11 shows tilings of this type and individual tile figures for three pairs of goal figures, for which tilings with relatively regular tile layout are presented in the last row of Figure 10. Note that these tilings were selected from the top-100 solutions obtained by the EST, but it is also possible to display only this type of tiling. While the individual tile figures appear somewhat unnatural compared to those presented in the last row of Figure 9, the artistry of the tilings is enhanced.

Fig. 11.

View Figure

Next, we present the results of the HH-EST. Figure 12 shows tile figures and the corresponding tilings obtained by the HH-EST. To save space, only the results for the six pairs of goal figures that produced relatively good tilings are presented, where the six pairs were chosen so that each goal figure was selected twice. See the online supplementary file 1 for the results for all 15 pairs of goal figures. Compared to the results of the EST (mostly fourth and fifth rows of Figures 9 and 10), the obtained tile figures seem to have unnatural deformations, but this is expected because the templates have fewer degrees of freedom in possible deformations. Intuitively, it is difficult to construct satisfactory tilings of Heaven and Hell patterns unless the two goal shapes are suitably compatible because the two types of tiles are always adjacent at the boundaries. Nevertheless, the obtained tile figures are sufficiently reminiscent of the goal figures, while the resulting tilings also have the artistry inherent in the Heaven and Hell patterns.

Fig. 12.

View Figure

6.3 Influence of the Values of \(n_1\) and \(n_2\)

For the main results, the numbers of boundary points of the goal figures were all set to 60 (\(n_1=n_2=60\)), but \(n_1\) and \(n_2\) can have different values for the EST. To verify the effect of different \(n_1\) and \(n_2\), mesh representations with \(n_1=90\) and 120 were also created.

Figure 14 shows the tile meshes \(\tilde{U}_1\) and \(\tilde{U}_2\) obtained for a pair of the goal figures elephant \((n_1=60, 90,\) and \(120)\) and seahorse \((n_2=60)\), where \(\tilde{U}_1\) and \(\tilde{U}_2\) are drawn so that the area of \(\tilde{U}_2\) is the same. As shown in this figure, the size of \(\tilde{U}_1\) increases as the value of \(n_1\) increases; the ratios of the area of \(\tilde{U}_1\) to the area of \(\tilde{U}_2\) are \(0.81 \ (n_1=60)\), \(2.52 \ (n_1=90)\), and \(3.82 \ (n_1=120)\). This is a natural consequence because each tile mesh is generated so that points of the boundary have an approximately equal interval and parts of the boundaries are shared by the two tile meshes. Figure 13 shows the tiling results obtained from the tile meshes presented in Figure 14.

Fig. 13.

View Figure

Fig. 14.

View Figure

6.4 Analysis of Search Restrictions on the EST

In the proposed algorithm for the EST, the optimization problem (16) is solved only for \(n_\mathrm{conf} \ (=10^6)\) template configurations selected in advance. Therefore, the true \(n_\mathrm{top} \ (=10)\) solutions that would be obtained by the EST without restrictions may not have been obtained. We want to verify how many of the true top-10 solutions were actually obtained. However, obtaining the true top-10 solutions was impractical, so we regard the top-10 solutions obtained when \(n_\mathrm{conf}\) was set to \(5\times 10^8\) as the pseudo true top-10 solutions. This assumption seems reasonable because no change was observed in the top-10 solutions when \(n_\mathrm{conf}\) was increased beyond \(10^8\) for the pairs of goal figures we tested below.

Table 2 shows how many template configurations that yielded the pseudo true top-10 solutions were included in the promising template configurations obtained for various values of \(n_\mathrm{conf}\). The results for the pseudo true top-100 solutions are also listed for reference. We can see that more than half of the pseudo true top-10 solutions were found with the default setting of \(n_\mathrm{conf}=10^6\). As shown in the table, increasing \(n_\mathrm{conf}\) reduces the number of pseudo true top-10 solutions that are overlooked. However, the computation time required to solve the optimization problem (16) increases in proportion to the value of \(n_\mathrm{conf}\), while the computation time for obtaining the \(n_\mathrm{conf}\) template configurations did not change significantly if \(n_\mathrm{conf} \lt 10^7\).

Finally, we mention the computation time for the EST with the default setting (\(n_\mathrm{conf}=10^6\)). Table 3 shows the computation times for (i) selecting template configurations and (ii) the optimization process for the selected template configurations, where the average computation times for the goal figure pairs in Table 2 are presented. Results are shown without and with parallelization, respectively. The parallel version was executed with 20 threads, although the effect of parallelization was almost saturated at 10 threads. If the tree search was not introduced in the algorithm of selecting template configurations, then the computation time without parallelization was 18.6 hours (\(n_1=60, n_2=60\)), which indicates that the pruning strategies are working effectively. As for the HH-EST, the computation time without parallelization for selecting template configurations was 8.1 seconds.

6.5 Comparison with Other Methods

As described in Section 1, at least four methods [Kaplan and Salesin 2004; Kawade and Imahori 2015; Lin et al. 2017; Liu et al. 2020] that can tackle the dihedral Escherization problem have been proposed. We compared the EST and HH-EST with Liu et al. [2020] because this method seems to produce the highest-quality dihedral Escher-like tilings among these four methods.

Liu et al. [2020] can only generate dihedral tilings of Heaven and Hell patterns. As described in their paper, the compatibility of two goal shapes is important to obtain satisfactory dihedral tilings of Heaven and Hell patterns, and one example that they considered difficult to obtain good results for was presented. Figure 15 shows this example, where two goal shapes and two tile shapes obtained are displayed. We applied the EST and HH-EST to the mesh representations (\(n_1 = n_2 = 75\)) created for this pair of goal shapes. The obtained tile shapes are also shown in Figure 15. As a result of a questionnaire of 12 people, 10 answered that the result of the EST, HH-EST, and Liu et al. [2020] were better, in that order.

Fig. 15.

View Figure