Inductive Proof
We shall show that
$$f_n=(1+x)\,(xy)^{-(n-1)}\,(1+xy)^{2n-1}$$
for all positive integers $n$. We start with
$$\begin{align}f_1&=1+x+y+xy=(1+xy)+(x+y)
\\&=(1+xy)+x\,(1+xy)=(1+x)\,(1+xy)
\\&=(1+x)\,(xy)^{-(1-1)}\,(1+xy)^{2\cdot 1-1}\,.\end{align}$$
Thus, the claim is true for $n=1$. We now suppose that the claim holds for $n=k$ for some positive integer $k$.
By induction hypothesis,
$$\begin{align}f_{k+1}&=(1+x+y+xy)\,f_k\\&=(1+x+y+xy)\,(1+x)\,(xy)^{-(k-1)}\,(1+xy)^{2k-1}\,.\end{align}$$
We can see that
$$(1+x+y+xy)\,(1+x)\,(xy)^{-(k-1)}=(1+x)(1+y)(1+x)\,(xy)^{-(k-1)}\,.$$
Now,
$$(1+y)(1+x)=1+x+y+yx=(1+yx)+y(1+yx)=(1+y)(1+yx)\,,$$
so
$$\begin{align}(1+y)(1+x)\,(xy)^{-(k-1)}&=(1+y)(1+yx)\,(yx)^{k-1}\\&=(1+y)\,(yx)^{k-1}\,(1+yx)\\&=(1+y)\,(yx)^k\,(1+xy)\,.\end{align}$$
Therefore,
$$\begin{align}f_{k+1}&=(1+x)\,\big((1+y)\,(yx)^k\,(1+xy)\big)\,(1+xy)^{2k-1}\\&=(1+x)(1+y)\,(yx)^k\,(1+xy)^{2k}\,.\end{align}$$
However, as $(1+x)(1+y)=(1+x)(1+xy)$, we get
$$(1+x)(1+y)\,(yx)^k=(1+x)(1+xy)\,(yx)^k=(1+x)\,(yx)^k\,(1+xy)\,.$$
This shows that
$$\begin{align}
f_{k+1} &=\big((1+x)\,(yx)^k\,(1+xy)\big)\,(1+xy)^{2k}
\\
&=(1+x)\,(xy)^{-\big((k+1)-1\big)}\,(1+xy)^{2(k+1)-1}\,.
\end{align}$$
The proof is now complete.
Now, we have
$$\begin{align}f_n&=(1+x)\,(xy)^{-(n-1)}\,(1+xy)^{2n-1}\\&=(xy)^{-(n-1)}\,(1+xy)^{2n-1}+x\,(xy)^{-(n-1)}\,(1+xy)^{2n-1}\,.\end{align}$$ The constant term of $f_n$ can only come from the constant term of $(xy)^{-(n-1)}(1+xy)^{2n-1}$. Therefore, the constant term of $f_n$ is
$$\binom{2n-1}{n-1}=\frac{1}{2}\,\binom{2n}{n}\,.$$
In fact, the coefficients of $(xy)^r$ and $x\,(xy)^r$ in $f_n$, where $r\in\mathbb{Z}$, are both equal to $$\displaystyle\binom{2n-1}{n-1+r}=\frac{n+r}{2n}\,\binom{2n}{n+r}\,.$$
Geometric Proof
For $\theta\in\mathbb{R}$, let $\rho_\theta$ denote the counterclockwise rotation of the Euclidean plane $E:=\mathbb{R}^2$ by the angle $\theta$. We also write $\sigma_\theta$ for the reflection about the line $\ell_\theta$ through the origin that makes the angle $\theta$ (measured in the counterclockwise direction) with the horizontal axis. Define $\varpi_\theta$ to be the orthogonal projection onto $\ell_\theta$.
Note that
$$\varpi_\theta=\frac{1+\sigma_\theta}{2}\,,$$ where $1$ also denotes the identity map on $E$. We have some basic identities:
$$\rho_{\theta_1}\rho_{\theta_2}=\rho_{\theta_1+\theta_2}\,,$$
$$\sigma_{\theta_1}\sigma_{\theta_2}=\rho_{2\theta_1-2\theta_2}\,,$$
$$\sigma_{\theta_1}\rho_{\theta_2}=\sigma_{\theta_1-\frac{\theta_2}{2}}\,,$$
and
$$\rho_{\theta_1}\sigma_{\theta_2}=\sigma_{\frac{\theta_1}{2}+\theta_2}\,,$$
for all $\theta_1,\theta_2\in\mathbb{R}$.
Let $\alpha$ and $\beta$ be variable angles. We consider $x:=\sigma_\alpha$ and $y:=\sigma_\beta$. Then, $x^2=1$, $y^2=1$,
$$xy=\sigma_\alpha\sigma_\beta=\rho_{2\alpha-2\beta}\,,\text{ and }yx=\sigma_\beta\sigma_\alpha=\rho_{2\beta-2\alpha}\,.$$
It is easily verified that no two terms of the form $(xy)^r$ or $x\,(xy)^r$, where $r\in\mathbb{Z}$, are equal (as functions of $\alpha$ and $\beta$).
For a positive integer $n$, we perform an $n$-step transformation on the Euclidean plane $E$ such that, at each step, we can leave $E$ alone (this action, or rather, nonaction, corresponds to the identity $1$), reflect $E$ via the map $x$ or $y$, or rotate it via $xy$. Let $T_n$ be the set of all possible $n$-step transformations.
We want to count all elements of $T_n$ that ending up leaving $E$ unchanged (that is, the composition of the successive maps performed at each step is the identity $1$). Let this number be $q(n)$. For each $\mathbf{t}\in T_n$, write $\tau_\mathbf{t}$ be the final mapping resulting from the $n$-step transformation $\mathbf{t}$.
Observe that, for all $v\in E$, we see that
$$\sum_{\mathbf{t}\in T_n}\,\tau_\mathbf{t}(v)=(1+x+y+xy)^n(v)=f_n(v)\,.$$
Now observe that
$$1+x+y+xy=(1+x)(1+y)=2^2\,\left(\frac{1+x}{2}\right)\,\left(\frac{1+y}{2}\right)\,.$$
Because
$$\varpi_\alpha=\frac{1+x}{2}\text{ and }\varpi_\beta=\frac{1+y}{2}\,,$$
we conclude that $$f_n(v)=2^2\,(\varpi_\alpha\varpi_\beta)^n(v)\,.$$
Observe that
$$\varpi_\alpha\varpi_\beta\varpi_\alpha=\big(\cos(\alpha-\beta)\big)^2\,\varpi_\alpha\,.$$
By identifying $E$ with $\mathbb{C}^2$, we can write each $v$ as $r\,\text{e}^{\text{i}\phi}$ for some $r\geq 0$ and $\phi\in\mathbb{R}$. Observe that
$$\varpi_\alpha\varpi_\beta(r\,\text{e}^{i\phi})=\cos(\alpha-\beta)\,\cos(\phi-\beta)\,\text{e}^{\text{i}\alpha}\,r\,.$$ This means
$$f_n(r\,\text{e}^{i\phi})=2^{2n}\,\big(\cos(\alpha-\beta)\big)^{2n-1}\,\cos(\phi-\beta)\,\text{e}^{\text{i}\alpha}\,r\,.$$
By writing
$$\begin{align}\cos(\phi-\beta)&=\cos\big((\phi-\alpha)+(\alpha-\beta)\big)\\&=\cos(\phi-\alpha)\,\cos(\alpha-\beta)-\sin(\phi-\alpha)\,\sin(\alpha-\beta)\,,\end{align}$$
we obtain
$$\int_{0}^{2\pi}\,f_n(r\,\text{e}^{i\phi})\,\text{d}\beta=2^{2n}\,\left(\int_0^{2\pi}\,\cos^{2n}(\alpha-\beta)\,\text{d}\beta\right)\,\cos(\phi-\alpha)\,\text{e}^{\text{i}\alpha}\,r\,.$$
From this result, we have
$$\begin{align}\int_{0}^{2\pi}\,f_n(r\,\text{e}^{i\phi})\,\text{d}\beta&=2^{2n}\,\frac{\pi}{2^{2n-1}}\,\binom{2n}{n}\,\cos(\phi-\alpha)\,\text{e}^{\text{i}\alpha}\,r\\&=2\pi\,\binom{2n}{n}\,\cos(\alpha-\phi)\,\text{e}^{\text{i}(\alpha-\phi)}\,r\,\text{e}^{\text{i}\phi}\,.\end{align}$$
Hence,
$$\begin{align}\int_0^{2\pi}\,\int_0^{2\pi}\,f_n(r\,\text{e}^{i\phi})\,\text{d}\beta\,\text{d}\alpha&=2\pi\,\binom{2n}{n}\,\left(\int_0^{2\pi}\,\cos(\alpha-\phi)\,\text{e}^{\text{i}(\alpha-\phi)}\,\text{d}\alpha\right)\,r\,\text{e}^{\text{i}\phi}\\&=2\pi^2\,\binom{2n}{n}\,r\,\text{e}^{\text{i}\phi}\,.\end{align}$$
That is,
$$\frac{1}{(2\pi)^2}\,\int_0^{2\pi}\,\int_0^{2\pi}\,f_n(r\,\text{e}^{i\phi})\,\text{d}\beta\,\text{d}\alpha=\frac{1}{2}\,\binom{2n}{n}\,r\,\text{e}^{\text{i}\phi}\,.$$
However, for each $\mathbf{t}\in T_n$,
$$\frac{1}{(2\pi)^2}\,\int_0^{2\pi}\,\int_0^{2\pi}\,\tau_\mathbf{t}(r\,\text{e}^{i\phi})\,\text{d}\beta\,\text{d}\alpha=\left\{\begin{array}{ll}
r\,\text{e}^{i\phi}&\text{if }\tau_\mathbf{t}=1\,,\\
0&\text{if }\tau_\mathbf{t}\neq 1\,.
\end{array}\right.$$
Therefore,
$$\begin{align}\frac{1}{(2\pi)^2}\,\int_0^{2\pi}\,\int_0^{2\pi}\,f_n(r\,\text{e}^{i\phi})\,\text{d}\beta\,\text{d}\alpha&=\sum_{\mathbf{t}\in T_n}\,\frac{1}{(2\pi)^2}\,\int_0^{2\pi}\,\int_0^{2\pi}\,\tau_\mathbf{t}(r\,\text{e}^{i\phi})\,\text{d}\beta\,\text{d}\alpha\\&=q(n)\,r\,\text{e}^{\text{i}\phi}\,,\end{align}$$
Hence,
$$q(n)=\frac{1}{2}\,\binom{2n}{n}\,.$$
Remark. Let $G$ be the group $\big\langle x,y\,\big|\,x^2=1\text{ and }y^2=1\big\rangle$. Denote by $\mathcal{E}$ the $\mathbb{R}$-vector space of continuous maps $F:(\mathbb{R}/2\pi\mathbb{Z})^2\to \text{End}_\mathbb{R}(E)$. Basically, this approach uses the fact that the map $\Phi:G\to \text{GL}(\mathcal{E})$ given by
$$\Big(\big(\Phi(x)\,F\big)(\alpha,\beta)\Big)\,u:=\sigma_\alpha \Big(F(\alpha,\beta)\,u\Big)$$
and
$$\Big(\big(\Phi(y)\,F\,\big)(\alpha,\beta)\Big)\,u:=\sigma_\beta \Big(F(\alpha,\beta)\,u\Big)\,,$$
where $\alpha,\beta\in \mathbb{R}/2\pi\mathbb{Z}$, $F\in\mathcal{E}$, and $u\in E$,
is a faithful representation of $G$.
