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.