Let $\Omega_0,\Omega_1\subset\mathbb R^n$ be bounded uniformly convex domains with $C^\infty$ boundaries. Let $\rho_0\in C^\infty(\overline{\Omega_0})$ and $\rho_1\in C^\infty(\overline{\Omega_1})$ satisfy
\begin{align*}
0<\lambda\le \rho_0(x)\le \Lambda<\infty,\qquad
0<\lambda\le \rho_1(y)\le \Lambda<\infty
\end{align*}
on their respective closures, and assume the compatibility condition
\begin{align*}
\int_{\Omega_0}\rho_0\,d\mathcal L^n=\int_{\Omega_1}\rho_1\,d\mathcal L^n.
\end{align*}
Let $\phi:\Omega_0\to\mathbb R$ be the Brenier potential for the quadratic transport from $\rho_0\,d\mathcal L^n$ to $\rho_1\,d\mathcal L^n$, normalised by adding a constant. Then $\phi$ is smooth in $\Omega_0$, and the global second boundary value regularity theorem gives a representative with $\phi\in C^\infty(\overline{\Omega_0})$ and $\nabla\phi:\Omega_0\to\Omega_1$ a smooth diffeomorphism.