Lebesgue Measure of Image of Unit Square under Continuous Map 
Problem. Let $h\in C(\mathbb{R})$ be a continuous function, and let
  $\Phi:\Omega:=[0,1]^{2}\rightarrow\mathbb{R}^{2}$ be the map defined
  by \begin{align*} 
\Phi(x_{1},x_{2}):=\left(x_{1}+h(x_{1}+x_{2}),x_{2}-h(x_{1}+x_{2})\right) \tag{1}
\end{align*} 
  What is the (Lebesgue) measure (denoted by $\left|\cdot\right|$) of the set $\Phi(\Omega)$?

Implicit in the problem statement is that $\Phi(\Omega)$ is Lebesgue measurable, but this is obvious as $\Phi$ is the composition of continuous functions, so $\Phi(\Omega)$ is a compact subset of $\mathbb{R}^{2}$.
If we knew that $h$ was $C^{1}$, then the Jacobian of $\Phi$ is
\begin{align*}
J\Phi(x_{1},x_{2})=\begin{bmatrix} {1+h'(x_{1}+x_{2})} & {h'(x_{1}+x_{2})} \\ {-h'(x_{1}+x_{2})} & {1-h'(x_{1}+x_{2})} \end{bmatrix}
\end{align*}
which has determinant $1$ everywhere. Whence,
\begin{align*}
\left|\Phi(\Omega)\right|=\int_{\mathbb{R}^{2}}\chi_{\Phi(\Omega)}(x_{1},x_{2})\mathrm{d}x_{1}\mathrm{d}x_{2}=\int_{\mathbb{R}^{2}}\chi_{\Omega}(x_{1},x_{2})\left|J\Phi(x_{1},x_{2})\right|\mathrm{d}x_{1}\mathrm{d}x_{2}=1
\end{align*}
My thought was to approximate $h$ uniformly in a neighborhood of $[-2,2]$ by a smooth function $g$, so that $\left|h(x_{1}+x_{2)}-g(x_{1}+x_{2})\right|<\epsilon$ for all $(x_{1},x_{2})\in\Omega$, given $\epsilon>0$. Denote the analogue of $\Phi$ with $g$ instead of $h$ by $\Phi_{\epsilon}$. One can verify that
$$\forall x=(x_{1},x_{2})\in\Omega,\quad \left|\Phi(x)-\Phi_{\epsilon}(x)\right|=\sqrt{2}\left|h(x_{1}+x_{2})-g(x_{1}+x_{2})\right|<\sqrt{2}\epsilon$$
So $\Phi(\Omega)\subset U_{\epsilon}:=\left\{y\in\mathbb{R}^{2} : d(y,\Phi_{\epsilon}(\Omega))<\sqrt{2}\epsilon\right\}$ and $\Phi_{\epsilon}(\Omega)\subset V_{\epsilon}:=\left\{y\in\mathbb{R}^{2} : d(y,\Phi(\Omega))<\sqrt{2}\epsilon\right\}$.
We obtain that
\begin{align*}
1=\left|\Phi_{\epsilon}(\Omega)\right|\leq\int_{\mathbb{R}^{2}}\chi_{V_{\epsilon}}(x)\mathrm{d}x
\end{align*}
Since $\Phi(\Omega)$ is compact, in particular closed, $\chi_{V_{\epsilon}}\rightarrow\chi_{\Phi(\Omega)}$ a.e. as $\epsilon\downarrow 0$. From monotone convergence, we obtain that
\begin{align*}
1\leq\left|\Phi(\Omega)\right| \tag{2}
\end{align*}
My issue is with establishing the reverse inequality. I know that
\begin{align*}
\left|\Phi(\Omega)\right|\leq\left|U_{\epsilon}\right|=1+\left|\left\{x\in\mathbb{R}^{2}:0<d(x,\Phi_{\epsilon}(\Omega))<\sqrt{2}\epsilon\right\}\right|
\end{align*}
But I do not know how to control the measure of set $\left\{x\in\mathbb{R}^{2}: 0<d(x,\Phi_{\epsilon}(\Omega))<\sqrt{2}\epsilon\right\}$ as $\epsilon\downarrow 0$, as $\Phi(\Omega)\setminus\Phi_{\epsilon}(\Omega)$ is contained in this set for all $\epsilon>0$.
 A: HINT:
$$\Phi^{-1}(y_1, y_2) = (y_1- h(y_1, y_2), y_2 + h(y_1, y_2)\,)
$$
Adding some details: 
If $\Phi_n \overset{\text{u}}{\to} \Phi$ and $K$ compact then $\mu( \Phi(K)) \ge \limsup_n \mu(\Phi_n(K))$ like you showed with $\epsilon$-neighborhoods. Now take $\Phi_n$ smooth, measure preserving like you showed, and get $\mu(\Phi(K)) \ge \mu(K)$.
Conversely, take $h_n$ smooth converging to $h$. Then $\Phi^{-1}_n \to \Phi^{-1}$. Therefore, we have 
$$\mu(K) = \mu (\Phi^{-1}(\Phi(K))) \ge \limsup_n \mu (\Phi_n^{-1}(\Phi(K))$$
However, the smooth $\Phi_n$ preserve the measure so $\mu (\Phi_n^{-1}(\Phi(K)))= \mu( \phi(K))$. 
$\bf{Added:}$ What OP used was to approximate $\Phi$ with smooth maps that preserve the area, and then show that $\mu(\Phi(\Omega))\ge \limsup_n \mu(\Phi_n(\Omega))$. This follows from the fact that $\Phi_n(\Omega)\to \Phi(\Omega)$ in the Hausdorff metric and then using the fact that the volume is upper continuous for that metric. However, we do have in general 
$$\mu(\Phi(\Omega))=\lim_{n\to \infty} \mu(\Phi_n(\Omega))$$
whenever $\Phi_n$ converges to $\Phi$ uniformly on the square $\Omega$. Indeed, for any continous map $\Psi$ defined on the square $\Omega$ the image $\Psi(\Omega)$ contains all the points $y \in \mathbb{R}^2$ whose index relative to the image of the boundary $\Psi(\partial\Omega)$ is not zero.
A: Here is another proof which is not so much in the spirit of the OP
(approximating $h$ by smooth functions), but which I nevertheless
find interesting, in particular since we can now allow $h$ to be
an arbitrary (Borel) measurable function. Furthermore, we show that
$\Phi$ preserves the measure of arbitrary Borel measurable sets.
Let us first factor the map $\Phi$. In fact, if we set
$$
\Psi:\mathbb{R}^{2}\to\mathbb{R}^{2},\left(x,y\right)\mapsto\left(\begin{matrix}x+h\left(y\right)\\
y
\end{matrix}\right),
$$
then we have
\begin{eqnarray*}
\left(\begin{matrix}1 & 0\\
-1 & 1
\end{matrix}\right)\Psi\left(\left(\begin{matrix}1 & 0\\
1 & 1
\end{matrix}\right)\left(\begin{matrix}x\\
y
\end{matrix}\right)\right) & = & \left(\begin{matrix}1 & 0\\
-1 & 1
\end{matrix}\right)\Psi\left(\begin{matrix}x\\
x+y
\end{matrix}\right)\\
 & = & \left(\begin{matrix}1 & 0\\
-1 & 1
\end{matrix}\right)\left(\begin{matrix}x+h\left(x+y\right)\\
x+y
\end{matrix}\right)\\
 & = & \left(\begin{matrix}x+h\left(x+y\right)\\
x+y-\left[x+h\left(x+y\right)\right]
\end{matrix}\right)\\
 & = & \left(\begin{matrix}x+h\left(x+y\right)\\
y-h\left(x+y\right)
\end{matrix}\right)=\Phi\left(\begin{matrix}x\\
y
\end{matrix}\right).
\end{eqnarray*}
Since the linear maps induced by the matrices $\left(\begin{smallmatrix}1 & 0\\
-1 & 1
\end{smallmatrix}\right),\left(\begin{smallmatrix}1 & 0\\
1 & 1
\end{smallmatrix}\right)$ are measure preserving, it suffices to show that $\Psi$ is measure
preserving. Note that $\Psi$ is Borel measurable with Borel measurable
inverse
$$
\Psi^{-1}:\mathbb{R}^{2}\to\mathbb{R}^{2},\left(x,y\right)\mapsto\left(\begin{matrix}x-h\left(y\right)\\
y
\end{matrix}\right).
$$
Now, for a Borel measurable set $M\subset\mathbb{R}^{2}$ and $y\in\mathbb{R}$,
let
$$
M_{y}:=\left\{ x\in\mathbb{R}\,\mid\,\left(x,y\right)\in M\right\} 
$$
be the section of $M$ at height $y$. By Cavallieri's principle (which
just means "by applying Fubini's theorem to the indicator function
$\chi_{M}$"), we then have
$$
\lambda_{2}\left(M\right)=\int_{\mathbb{R}}\lambda_{1}\left(M_{y}\right)\,{\rm d}y,
$$
where $\lambda_{1},\lambda_{2}$ denote the $1$- and $2$-dimensional
Lebesgue measure.
Now, note
\begin{eqnarray*}
x\in\left[\Psi\left(M\right)\right]_{y} & \Longleftrightarrow & \left(x,y\right)\in\Psi\left(M\right)\\
 & \Longleftrightarrow & \left(\begin{matrix}x-h\left(y\right)\\
y
\end{matrix}\right)=\Psi^{-1}\left(x,y\right)\in M\\
 & \Longleftrightarrow & x-h\left(y\right)\in M_{y}\\
 & \Longleftrightarrow & x\in M_{y}+h\left(y\right),
\end{eqnarray*}
so that translation invariance of the Lebesgue measure yields
$$
\lambda_{1}\left(\left[\Psi\left(M\right)\right]_{y}\right)=\lambda_{1}\left(M_{y}+h\left(y\right)\right)=\lambda_{1}\left(M_{y}\right).
$$
By Cavallieri's principle as above, we conclude 
$$
\lambda_{2}\left(\Psi\left(M\right)\right)=\int\lambda_{1}\left(\left[\Psi\left(M\right)\right]_{y}\right)\,{\rm d}y=\int\lambda_{1}\left(M_{y}\right)\,{\rm d}y=\lambda_{2}\left(M\right).
$$
So essentially, the idea was that since all $\Psi$ does is a measurable
shift (depending on $y$) in $x$-direction, all sections of $\Psi\left(M\right)$
are just translates of the sections of $M$ and hence $\lambda_{2}\left(\Psi\left(M\right)\right)=\lambda_{2}\left(M\right)$
by Cavallieri.
Concluding remark: Since $\Psi$ is measure preserving on Borel sets,
it is not hard to show that $\Psi$ indeed maps Lebesgue measurable
sets to Lebesgue measurable sets and is also measure preserving on
these sets.
