I don't have a good reference at hand but I can explain the procedure. I'll use a quadruple to denote a Wang tile $T = (N,E,S,W)$ referring to the North, East, South and West edges of a tile respectively.
Easy direction: Given a set of Wang tiles $\{T_1, \ldots, T_k\}$, let $B_h$ be the set of pairs $(T_i,T_j)$ such that $T_i(E) \neq T_j(W)$ and let $B_v$ be the set of pairs $(T_i,T_j)$ such that $T_i(S) \neq T_j(N)$. Now let $\mathcal{F}_h$ be the patch of tiles with $T_j$ to the right of $T_i$, where $(T_i, T_j) \in B_h$ and let $\mathcal{F}_v$ be the patch of tiles with $T_j$ below $T_i$, where $(T_i,T_j) \in B_v$. We may now take our set of forbidden patches $\mathcal{F}$ to be given by $\mathcal{F}_h \cup \mathcal{F}_v$.
Harder direction: Let $\mathcal{A}$ be our alphabet. Given a set of forbidden patches $\mathcal{F} = \{P_i \mid 1 \leq i \leq k\}$, without loss of generality, we may assume that each of the patches is supported on an $n \times n$ square (as otherwise, just take all possible extensions of the patches to an $n \times n$ square such that $n$ is big enough to contain all the patches). Now, let $\mathcal{L}^{n,n}$ be the set of legal $n \times n$ patches. That is, $\mathcal{L}^{n,n} = \mathcal{A}^{n,n}\setminus \mathcal{F}$. Let $P_N$ be the patch one gets by deleting the northern-most row of $P$. Similarly for $P_E$, $P_S$ and $P_W$ for their respective rows/columns.
For every $P \in \mathcal{L}^{n,n}$, we define a Wang tile $T_P$ by
$$T_P(N) = (P_S, P_W, P_N, P_E).$$
So, the northern edge of $T_P$ contains all of the information of the patch $P$ except for the southern-most row, etc.
It's not too hard to see that two Wang tiles $T_P$ and $T_Q$ can only be placed next to each other if their corresponding North-South or East-West pair are identical, which is the same as saying that the patches $P$ and $Q$ overlap in a legal way.
There is now a simple local rule that takes one from a tiling with the Wang tiles $T_P$, $P \in \mathcal{L}^{n,n}$ to an element of the original SFT $X_{\mathcal{F}}$ whereby one just needs to rebuild the original patch $P$ from $T_P$ (which one can do using the labels of each of the edges) and then outputs the central element of $P$ (or round down if $n$ is even).
No, you cannot three-color that tiling.
Here's a finite part of the tiling from page 10 top left of the article, the tile numbered 2 here is the one darkened on that figure. This part cannot be three-colored.
To prove this, start from the three adjacent tiles marked 0, 1, 2, and call their colors yellow, blue and red respectively. Tile 0 is surrounded by a sequence of adjacent tiles 1, 2, 3, 4, 5, 6, so these must be colored blue and red alternatingly, 5 being blue, 6 being red. Tile 7 is adjacent to tile 5 and 6 and so must be yellow; while tile 8 is adjacent to tile 1 and 6 and so also must be yellow. Now tile 9 is adjacent to tiles 6 and 8 and so must be blue, but tile 10 is adjacent to 6 and 7 and so also must be blue, but 9 and 10 are adjacent so this coloring is invalid.
Best Answer
I wrote this in a hurry and am sure it's hard to read and has dozens of typos, feel free to ask for clarifications.
As I mentioned, the periodic tiling problem is undecidable whether or not you require that all tiles are used. I guess what you wanted to know was whether the periodic tiling problem is undecidable if we do not necessarily require that all tiles are used in the periodic tiling. I know that this is undecidable, and I have only seen a proof of this variant. I don't actually know what's in the paper the linked MO answer cites, but I assume it's this variant.
Now let me clarify my comment. I claimed that if we add the assumption that all tiles are used in the tiling, the problem stays undecidable. This comes from the following theorem.
Properties $P$ invariant under scaling include more or less anything that deals with things like tileability, periodicity or computability of configurations.
What I mean by the promise stuff is that if $P$ is undecidable, then it stays undecidable even if you are only required to answer correctly assuming that all tiles are used in every valid tiling, and for tile sets without this property you are allowed to answer whatever. It is trivial that if the promise variant is undecidable then the original is as well, because if you know that you cannot even solve $P$ when you are only given tile sets where all tilings use all tiles, then of course you cannot solve $P$ for general tile sets. The other direction requires a proof, and I give one below.
This theorem in particular implies that the periodic tiling problem stays undecidable if we require that all tiles are used, because if you could solve the problem of existence of a periodic tiling using all tiles, then you could solve the promise version of the periodic tiling problem (as under the promise it does not matter whether we require all tiles in the solutions, since they are going to have all tiles anyway).
Note that I do not claim in general that if you have an undecidable problem $P$ that asks about the existence of a single tiling (with some additional property like periodicity) that uses all tiles, then also the variant of $P$ where we ask for a general tiling with the same property is undecidable**.
For example for the usual tiling problem the variant that requires all tiles are used is somehow easier to prove undecidable: the seeded tiling problem is roughly equivalent to tiling under this assumption, and Wang correctly solved the seeded tiling problem and incorrectly conjectured the answer to the general problem...
**: This is also true in high generality since most things about tilings are undecidable in both variants anyway, but one would have to make it more precise what kind of properties we are looking at.
Now let me explain the proof of the theorem.
For a set of Wang tiles $T$, write $X_T$ for the subshift of finite type they define, i.e. the set of all valid tilings $X_T \subset T^{\mathbb{Z}^2}$. Under the shift action $\sigma_{(a,b)}(x)_{(c,d)} = x_{a+c,b+d}$, $(\mathbb{Z}^2, X_T)$ is a topological dynamical system, meaning $\mathbb{Z}^2$ acts continuously and $X_T$ is a compact metrizable space. For two topological dynamical systems $(\mathbb{Z}^2, X)$, $(\mathbb{Z}^2, Y)$, write $(\mathbb{Z}^2, X) \cong (\mathbb{Z}^2, Y)$ if these systems are topologically conjugate, meaning there is a homeomorphism $\phi : X \to Y$ which conjugates the actions like $\phi(\sigma_{(a,b)}(x)) = \sigma_{(a,b)}(\phi(x))$ for all $x \in X$.
For a topological dynamical system $(\mathbb{Z}^2, X)$, write $\sqrt[(m,n)]{X}$ for the scaling (of shape $(m,n)$), namely the system whose points are equivalence classes $X \times \mathbb{Z} \times \mathbb{Z} / \equiv$ where $(x,a,b) \equiv (x,a',b') \iff (a \equiv a' \bmod m \wedge b \equiv b' \bmod n)$, the action of shift $(1, 0)$ is $(x,a,b) \mapsto (x,a+1,b)$ and the action of $(0,1)$ is $(x,a,b) \mapsto (x,a,b+1)$. If $X$ is a subshift, intuitively the scaling $\sqrt[(m,n)]{X}$ is just the system where we have replaced symbols in configurations of $X$ with $m \times n$-blocks.
The theorem follows immediately from this lemma. Namely, for the nontrivial direction of the theorem, if you can solve the promise version of a problem $P$, given an instance $T$ for the usual version of $P$, construct $U$ and solve the promise version of $P$ for it. Assuming $P$ is preserved under scaling $X \mapsto \sqrt[(m,n)]{X}$, this indeed solves the original problem.
Possible there is some cheap trick that gives this lemma, but I don't see it, so my construction is a little complicated. The general form of the macrotile is inspired by [1], though the proof is not really related to that, rather we just need various very basic tricks from self-assembly tilings, tilings and cellular automata theory.
Proof sketch. The idea is to simulate the original tiles $T$ with large macrotiles that mostly look the same but carry different data (coding a tile from $T$). Intuitively, we simulate the original tiles in software. We just need to ensure that the hardware uses all of its basic computational resources (tile types) in the computation, no matter what tiles $T$ we are simulating. I will ensure that the computation uses the same tiles always, by performing the computation that depends on the particular tile $t \in T$ being simulated with a simple unary counter (the set $T$ will be coded into a look-up table basically).
First, look at this picture, which represents a typical macrotile that simulates a tile from $T$:
(This is drawn freehand with touchpad in MS Paint in about 10 seconds, I may draw a proper picture later, or not.) This represents a very large $m \times n$ area (it's not actually a rectangle due to the complex boundary, but it is a contiguous fundamental domain for the $m \times n$ grid). The square(ish) gray area in the middle is the computation zone where we will perform computation. The pink lines are wires that will route information between adjacent macrotiles. They are only one cell thick, and I've drawn just two on each side but in general you'll need $k = \lceil \log_2 c \rceil$ of them, where $c$ is the number of colors in $T$-tiles (we assume a Wang tile is determined by its colors).
Note that if we ensure that all tiles belong to such a macrotile, then already just due to the complex shape of the boundary, all tilings will consist of a regular grid of macrotiles. If we then make each macrotile correspond to a unique $T$-tile and synchronize its colors with its neighbors (by sending bits along the wires), then clearly the resulting system is topologically conjugate to $\sqrt[(m,n)]{X_T}$ as required (where again $m \times n$ are the dimensions of the macrotile). Of course we then need to ensure that all tiles are used in all tilings.
The main body of the macrotile (non-wire, non-computation zone) is hard-coded: In each position we use a unique tile, and we synch all information between two such tiles. So if a tile from one of the connected components of the macrotile (separated by the wires and computation zone, drawn in green in the picture) appears in a tiling, it builds the entire component around it. The Wang tile colors used at the boundary of the macrocell don't actually matter as long as they do not prevent putting two macrotiles next to each other.
Now let's discuss the wires. For the purpose of discussion, orient the left and top wires from the boundary to the computation zone, and the west and bottom wires from the computation zone to the boundary. The colors to the left of a wire (w.r.t. this orientation) are of a different color than the ones below, but (importantly) we use the same color throughout the left (resp. right) side of the wire slot. We will also have only a small number of wire tiles that form the actual wire, so that we can easily ensure all of them are used in every tiling. I show a couple of examples of wire tiles:
Here, red is the left color and blue is the right color. We send information by alternation of a black and gray color. Here we see a wire tiles for when the wire goes straight to the right, and then the wire tiles for when a wire is going to the right and turns upward. So the point is that we send a bit of information by alternating 0 and 1 (represented by black and gray in the picture), which means that we will see both types of wire tiles no matter which bit we are sending. Note that some of the wire tiles turn, so we need to make the wire squiggly enough that also all turns that happen happen at both parities. For different wires (of which we have $4k$ many in total), we may use different colors on their sides. Then it is easy to see that if even one wire tile or tile belonging to the body appears in a valid configuration, then a complete macrotile must appear containing it.
We have also positioned the wires on the west and east so that they match exactly, so that for two neighboring macrotiles, the wires send and receive information between neighboring macrotiles. If total length of the wire between computation zones of neighboring macrotiles (containing one wire from each of the two neighboring macrotiles) is even, then you see the same bit in the wire tile color at both ends.
What is left is to figure out what we should put in the hole that is the computation zone (the gray rectangle in the first picture). The idea is that we have to some bits on the bottom coming from wires, $4k$ bits in total (not necessarily in consecutive positions), and some of the resulting binary strings should allow us to tile computation zone (precisely those combinations that code a tile of $T$), others not.
We will have two things going on in the computation zone. Vertically, we propagate a counter that starts at the bottom (so just one color per height). The point of this is just so that the computation zone tiles alone do not allow any infinite tilings, since going upward you run out of counter values and have to write in some body tile. We have this running on a separate layer of the tile set. A computation zone tile should not know its horizontal position, however, because then there would be a unique tile type in each position, and since we want all tiles to appear in each macrotile, this is problematic. Note that we can pick more or less freely where the wire bits come in at the bottom of the computation zone (by choosing how we route the wires), and in other positions we can feed in arbitrary constant symbols.
So all in all, the current situation is that on the bottom of the computation zone, we have a sequence of symbols [ w_0 b_1 w_1 b_2 ... b_{4k} w_{4k} ] over some finite alphabet $S \subset \mathbb{N}$ with $S \ni 0,1$ in the bottom colors, where $b_i \in \{0, 1\}$ is the $i$ bit, and $w_i \in S^*$ are some constant words. We will now start simulating a nearest-neighbor cellular automaton, i.e. in the sequence of colors we apply in each step a transformation where each cell's new value only depends on its own value and the values of neighboring cells. The reason we do computation with a cellular automaton rather than say a Turing machine is that it's easy to guarantee that we see all tile types on each row, by making sure each symbol triple appears somewhere (which is trivial to guarantee by picking the words $w_i$). I won't explain how to simulate a cellular automaton with a tile set, but basically just use the side colors to share values of neighboring tiles.
All we want to achieve with our cellular automaton is that the effect of each $b_i$ is amplified. More precisely, in the first $C$ steps (steps refer to the height at which we are in the computation zone) the $b_i$s do not interact (say, the distance $|w_i|+1$ between $b_i$ and $b_{i+1}$ is larger than $10 C$), and the effect of flipping $b_i = 0$ to $b_i = 1$ is that we add some $m_i$ to the total sum of the symbols at step $C$. We also want that for $i < j$, the effect of $i$ is amplified much more than the effect of $j$. In fact we want $m_j > \sum_{i < j} m_i$ (so that from the final sum we can deduce the values $b_i$). To get such behavior from a cellular automaton, we can take $S = \{0,1,2\}$, and use a simple rule where each cell with a bit (i.e. value in $\{0,1\}$) just copies the value from the right neighbor if the right neighbor also contains a bit, and fixes its value if the right neighbor contains $2$. And $2$s always just stay put. Now we can put each $b_i$ at the right end of an interval between two $2$s, and pick the lengths of intervals suitably, and indeed the effect of $b_i$ will be amplified by the length of the interval. To make sure that all tiles appear on each row (which we need to do, since the tiles know their height), for each zone containing a $b_i$ we can have another equivalent zone containing $0$ and $1$ in the corresponding slot.
At the end of step $C$, we turn the $2$s into $1$s (i.e. we use a different rule on this step; this is fine since we know the vertical height due to the height counter). So now we just have a sequence of bits again. Because $m_j > \sum_{i < j} m_i$, the total sum of these symbols determines the original symbols $b_i$. Now, we start using a different, even simpler, cellular automaton after the $C$th step. Namely we can simply use the traffic jam CA that moves the $1$-bits greedily to the left, i.e. the elementary CA rule $$ \begin{array}{lll}1 & 1 & 1 \\ & 1 & \end{array}, \; \; \begin{array}{lll}1 & 1 & 0 \\ & 1 & \end{array}, \; \; \begin{array}{lll}1 & 0 & 1 \\ & 1 & \end{array}, \; \; \begin{array}{lll}1 & 0 & 0 \\ & 0 & \end{array}, \; \; \begin{array}{lll}0 & 1 & 1 \\ & 0 & \end{array}, \; \; \begin{array}{lll}0 & 1 & 0 \\ & 0 & \end{array}, \; \; \begin{array}{lll}0 & 0 & 1 \\ & 1 & \end{array}, \; \; \begin{array}{lll}0 & 0 & 0 \\ & 0 & \end{array} $$
It is easy to ensure that the original values of $b_i$ do not affect what symbols are seen on each step, as they make up a small fraction of the total mass. Finally, at a suitable step $C'$ where we are sure all the relevant $1$s (dealing with the mass coming from the $b_i$s) have converged to the left of the computation zone, we end the computation, and simply send a signal upward from the rightmost $1$ (again by using a different rule depending on the height).
There are a couple of details here: First, at the step $C'$, the mass coming from the $b_i$'s has to have converged, so typically (if some $b_i = 0$) it has already stayed put for some time when step $C'$ is reached. If there is nothing else going on, then this will definitely lead into some rows not using all the tiles. So we need to have a lot of $1$s still coming in from the right at this point (which will not have time to join the mass on the left by step $C'$ no matter what the values of the $b_i$ are).
Second, the upward signal sent from the rightmost $1$ of the leftmost run of $1$s on step $C'$ needs to be sent immediately on this step, i.e. we cannot use a cellular automaton that sends a signal from the leftmost row to the right, as then the length of the run would determine when the signal is sent and now it would depend on the initial values of the $b_i$ which row sends this signal. Using Wang tiles, we can luckily send also horizontal signals, so this is indeed possible.
On step $C'$, we can safely erase all computation other than the upward signal. Finally, note that on the north border of the computation zone the tiles are part of the body. We can choose whether a particular combination of the $b_i$ leads to a valid tiling by picking whether the tile that the signal reaches with that combination allows a bottom color that receives the signal.
I have explained during the construction why indeed the tilings of with the resulting tile set are in one-to-one correspondence with tilings of $T$, and at each step we were careful (hopefully careful enough) to ensure that all tiles are used in every valid tiling. Square.
References
[1] Durand, Bruno; Romashchenko, Andrei; Shen, Alexander, Fixed-point tile sets and their applications, J. Comput. Syst. Sci. 78, No. 3, 731-764 (2012). ZBL1244.05049.