The link that you refer to does not describe the boundary as the attractor of a simple IFS, rather it describes a collection of portions of the boundary as the invariant list of a digraph IFS - or directed graph iterated function system. I am not an expert on sofic systems, but I believe that an analysis of that digraph IFS yields the type of description that you want.

Also, I can't agree with your dimension computation as your polynomial factors
$$2 x^3-x+1=(x+1) \left(2 x^2-2 x+1\right)$$
revealing only one real root, which is negative. The analysis below computes the dimension using something that looks like an entropy computation.

**Digraph self-similarity**

Here's an image of the twin dragon surrounded by 6 copies of itself.

The boundary consists of 6 parts, each of which is an intersection between the original twin dragon and one of the translated copies. That collection of 6 sets is exactly the invariant list of the following digraph IFS:

Thus, for example, the second piece consists of one copy of itself together with two copies of the first. To each edge corresponds one of two possible functions, namely
$$f_0(z)=\frac{1+i}{2}z \: \text{ or } \: f_1(z)=\frac{1+i}{2}(z + 1).$$

To compute the dimension of the boundary, we count the number $N_n$ of walks of length $n$ through the graph. The graph is strongly connected so it's adjacency matrix is irreducible. The Perron-Frobenius theorem guarantees that there is a positive eigenvalue $\lambda$ that is strictly larger than the absolute values of all the other eigenvalues. Furthermore, $N_n$ grows like $\lambda^n$. The size of a neighborhood generated by a path of length $n$ is $2^{-n/2}$. The dimension of the boundary is then given by
$$\lim_{n\rightarrow\infty}\frac{\log(N_n)}{\log(2^{n/2})} = 2\frac{\log(\lambda)}{\log(2)}.$$
The characteristic polynomial of the adjacency matrix is
$$\lambda ^6-2 \lambda ^5+\lambda ^4-4=(\lambda +1)
\left(\lambda ^2-2 \lambda +2\right) \left(\lambda
^3-\lambda ^2-2\right).$$
The largest root of the characteristic polynomial coincides with the largest root of $\lambda^3-\lambda ^2-2$ and is approximately $\lambda=1.69$. The dimension of the boundary is
$$2\frac{\log(\lambda)}{\log(2)} \approx 1.52363.$$

As a sanity check, we should be able to generate a portion of the boundary using your formulation. It takes just a few lines of Mathematica code to do so and it seems worth a look.

Here's the whole set:

```
n = 15;
TPic = ListPlot[{Re[#], Im[#]} & /@
(Tuples[{0, 1}, n].Table[((1 + I)/2)^k, {k, 0, n - 1}]),
PlotStyle -> RGBColor[1/5, 1/3, 4/5]]
```

Note that `Tuples[{0, 1}, n]`

generated a list of all tuples of zeros and ones of length $n$. We then took a dot product of that with a list of powers of $(1+i)/2$. Let's now try to do the same thing but, rather than using *all* tuples, we'll use only those tuples generated by a walk through the directed graph above starting at position 6. The code is a bit more involved than I should post here but that's the basic idea behind how I generated the the following pic:

Looks promising.

The description of the boundary of a self-similar tile as a collection of digraph self-similar sets was formulated the paper of Strichartz and Wang below. I wrote an exposition with Mathematica implementation here.

Strichartz,
R. and
Wang,
Y. Geometry
of
self-affine tiles I.
Indiana University
Mathematics
Journal.
1999
.
7:1-2
3.