For a function $f \in W^{1,p}(\mathbb R^N)$, it is well-known that there exists a constant $C_N$ (dependent on $N$) such that $$ |f(x)-f(y)| \le C_N|x-y|(\mathcal M|\nabla f|(x) + \mathcal M|\nabla f|(y)), $$ where $\mathcal M$ is the Hardy-Littlewood maximal function.
How can this estimate be improved if we assume $f \in L^p(\mathbb R^N)$ and $\nabla f \in BMO(\mathbb R^N)$ instead of $f \in W^{1,p}(\mathbb R^N)$?