choice of data structure for domino tilings
There is an alternative, more compact kind of array representation which might be better if you are working with diamonds as defined in the paper you linked (for axis aligned rectangles you need to find the smallest covering diamond).

The idea is to checkerboard the plane and assign a 2-bit number to each "black" cell, such that the number represents the orientation of the domino relative to the cell. Consider for example:


where $0$ means "horizontal domino pointing east", $1$ means "vertical domino pointing north", etc. Tilting the numbers clockwise yields the array representation:

$$A = \begin{pmatrix} 0 &amp; 0 &amp; 0 &amp; 0 \\ 2 &amp; 3 &amp; 3 &amp; 0 \\ 2 &amp; 3 &amp; 0 &amp; 2 \\ 2 &amp; 1 &amp; 0 &amp; 1 \\ 2 &amp; 2 &amp; 2 &amp; 1 \\ \end{pmatrix}$$

To find two dominos forming a square search for a $0$ adjacent to a $2$ or a $1$ adjacent  to a $3$. Now, if you deal with big tesselations and are into implemtations sped up by bit operations, you can do this:

Let $A$ be a $\{0,1,2,3\}^{m \times n}$ matrix, $b$ be the bit mask $\dots 01010101$, then we can represent the $i^{th}$ row of $A$ as a single integer by using bit-concatenation:

$$r_i = \bigvee_{j=1}^n (a_{ij} &lt;&lt; 2(j-1))$$

where $\vee$ denotes a bitwise OR and $&lt;&lt;$ a leftshift. Denoting $\oplus$ as bitwise XOR we can calculate

$$s_i=(r_i \text{ mod } 2^{2n}) \oplus (r_i &gt;&gt; 2) \oplus b \\
t_i=r_i \oplus r_{i+1}  \oplus b, ~~i &lt; m$$

Then domino $a_{ij}$ and domino $a_{ij+1}$ form a square iff bits $(s_i)_{2j-2}$ and $(s_i)_{2j-1}$ are both $1$. Likewise domino $a_{ij}$ and domino $a_{i+1j}$ do so iff bits $(t_i)_{2j-2}$ and $(t_i)_{2j-1}$ are both $1$.

Example calculation for rows 1,2 (little-endian!):

r_1    00 00 00 00
r_2    10 11 11 00
b      01 01 01 01
xor    -----------
       11 10 10 01


Therefore $a_{11}$ and $a_{21}$ form a square.