I am trying to construct the Rauzy Fractal (http://en.wikipedia.org/wiki/Rauzy_fractal), I have a Tribonacci word generator and have the stairs constructed but I can't seem to get the projection onto a 2D plane correct.
On the 4th slide of this: http://kmwww.fjfi.cvut.cz/jn08/slides/Thuswaldner.pdf they depict a bounding cube, my assumption was that I would project each point (x,y,z) onto (x',y') along the line (xMax, yMax, zMax)->(0,0,0). But they mention a Contracting Plane what does this terminology mean? Are all points projected along the same line or each to their own?
For this I am using the unit cube orthogonal projection matrix shown on the wiki: http://en.wikipedia.org/wiki/Orthographic_projection with Ortho(-1, 1, -1, 1, 1, -1) being my cube and (x/xMax, y/yMax, z/zMax) being my normalised point to be projected. I think I've got this part completely wrong. The result would agree...
Code: https://bitbucket.org/snippets/NeuralOutlet/8qEa/rauzy-fractal-problem

Should I even be using the projection matrix?
In this answer, I'm going to outline two approaches to generate the Rauzy fractal via a projection. The first, I think, answers your question more directly describing both the direction along which to project and the oblique plane to projection onto. There's also a related orthogonal projection which, I think, is simpler and easier to code.
Speaking of code, all the code to generate the pictures in this post, along with several other worked examples, can be found in this Observable notebook.
The direct approach
First off, I guess you're dealing with the following substitution: $$1\to 12, \: 2\to 13, \: 3\to 1.$$ This extends to words via concatenation. For example: $$ \sigma(1213) = 12\:13\:12\:1 = 1213121. $$
We can iterate this operation from an initial seed, say
1to generate a sequence of words:This converges to a fixed word of infinite length, call it $W$. We can take the 1s, 2s, and 3s in this sequence as a list of forwards, lefts, and ups to create an infinite staircase through 3D that looks something like so:
The color of each point indicates the direction of the last step. Note that the same symbol 1, 2, or 3 never appears too many times in a row in the word $W$. As a result, the staircase never meanders too far away from an average line, as drawn in the figure. Using this, we can now describe how the Rauzy fractal is generated at an intuitive level at least. The idea is to project the points in the staircase along the average line and onto a plane. When we do so, we get the following image of the Rauzy fractal in 3D space:
In order to describe the projection in more detail, you need to understand how the substitution is related to a so-called incidence matrix, which for your example is $$ M = \left( \begin{array}{ccc} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{array} \right). $$
More generally, the entry $m_{ij}$ in row $i$ and column $j$ of $M$ is $$ m_{ij} = \# \text{ occurrences of } i \text{ in } \sigma(j). $$ Now, the eigenvalues of $M$ are $$ \begin{aligned} \lambda_0 &\approx 1.839 \text{ and} \\ \lambda_{\pm} &\approx -0.42 \pm 0.606,i. \end{aligned} $$
The one real eigenvalue tells us there is a one-dimensional subspace that is invariant under the action of $M$ and that, furthermore, $M$ is expansive on that subspace, since the eigenvalue is larger than $1$. The pair of complex conjugate eigenvalues tell us there is a two dimensional subspace that is also invariant that involves some kind of rotational behavior. The fact that $|\lambda_{\pm}|<1$ tells us that $M$ is eventually contractive on that subspace, i.e. $M^k$ is contractive on the subspace for sufficiently large $k$.
We can get more information on the geometric action of $M$ from the eigenvectors associated with the eigenvalues: $$ \vec{v}_0 \approx \begin{pmatrix}1\\0.544\\0.296\end{pmatrix} \: \text{ and } \: \vec{v}_{\pm} \approx \begin{pmatrix}1\\-0.772\\-0.648\end{pmatrix} \pm \begin{pmatrix}0\\-1.115\\1.721\end{pmatrix}i $$
Specifically, the vector $\vec{v}_0$ spans the invariant one-dimensional subspace of $\mathbb R^3$ on which $M$ is expansive. This one-dimensional subspace also happens to be exactly the line drawn in the illustration of the staircase above; that is the direction along which we should project the staircase to obtain the 3D image of the Rauzy fractal. The real and imaginary parts of the complex eigenvector $\vec{v}_{\pm}$ span the invariant two dimensional subspace on which $M$ is eventually contractive. That is the subspace onto which we should project the staircase to obtain the 3D image of the Rauzy fractal.
To actually perform the projection, we'll use multiplication by a projection matrix. To form that matrix, we can use the column vectors just described as a change of basis matrix applied to a simple projection onto the $yz$-plane. Thus, our projection matrix $P$ is
$$ \begin{aligned} P &\approx \begin{pmatrix}1&1&0\\0.544&-0.772&-1.115\\0.296&-0.648&1.721\end{pmatrix} \begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{pmatrix} \begin{pmatrix}1&1&0\\0.544&-0.772&-1.115\\0.296&-0.648&1.721\end{pmatrix}^{-1} \\ &\approx \begin{pmatrix}0.382&-0.519&-0.336\\-0.336&0.718&-0.183\\-0.183&-0.153&0.901\end{pmatrix}. \end{aligned} $$
The picture above that shows the Rauzy fractal in 3D was produced by multiplying each of the points in the staircase by $P$.
An alternative projection
Since the Rauzy fractal lives in a plane, it's natural to try to generate a two-dimensional picture of it. One simple way to do so would be to rotate the 3D image we already have so that it is parallel to the $xy$-plane and then drop the third coordinate. There is a much more direct approach, though, that takes advantage of the complex structure of the eigenvectors and leads to a much simpler computation and code.
As it turns out, if we take the transpose of the matrix $M$ and compute its complex eigenvectors, then they are orthogonal to the one-dimensional real eigenspace of $M.$ This suggests that we might project the staircase onto one of those, grab the scalar component, and plot the real and imaginary components in the plane.
Now, there's a very simple formula for the projection of a vector $\vec{v}$ onto the vector $\vec{u}$ that takes on an even simpler form when $\vec{u}$ is normalized so that its magnitude is one, namely: $$ \text{proj}_{\vec{u}}\vec{v} = \frac{\vec{v}\cdot\vec{u}}{\vec{u}\cdot{\vec{u}}}\vec{u} = (\vec{v}\cdot\vec{u})\vec{u}. $$ We just want the scalar component, namely the $\vec{v}\cdot\vec{u}.$ Thus, once we have the stair case and matrix $M,$ all we need to do is compute a complex eigenvector of $M^T,$ normalize it to obtain $\vec{u}$ and then compute the dot product $\vec{u}\cdot\vec{v}$ where $\vec{v}$ ranges all the points in the staircase.
Here's the result:
Finally, it's worth mentioning that this second approach is equivalent to a process called valuation in this paper by Sirvent and Wang. They, too, work with the transposed incidence matrix $$ M^T = \left( \begin{array}{ccc} 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 0 \\ \end{array} \right) $$ one of whose eigenvectors is approximately $$\omega = \langle \omega_1, \omega_2, \omega_3 \rangle =\langle -0.412 + 0.61 i, 0.223 - 1.12i, 1\rangle.$$
Given a finite word $U$ they define $|U|_i$ to denote the number of occurrences of the symbol $i$ in $U$. Then, according to equation 1.3 in the paper, the valuation $E$ is defined by $$E(U) = \sum_{i=1}^3 |U|_i \omega_i.$$ Note that the $\omega_i$s are simply the components of the vector $\omega$ and have complex values. So clearly, $E(U)$ is a complex number.
Now, if we start with a long (but finite) approximation to the infinite fixed word $W$ and apply valuation all the left shifts that word, we get the exact same picture.