I'm interested in a mechanical procedure (i.e., algorithm) for sampling uniformly from $\{y \in \mathcal{M} : d_\mathcal{M}(x, y) \le 1\}$ for some fixed $x \in \mathcal{M}$, using well-known tools. The measure is the one induced by the corresponding Riemannian metric $\{g_x\}_{x \in \mathcal{M}}$.
When does such an algorithm exist and when is it non-trivial to find one?
We could start from the following schema:
Sample a tangent vector uniformly from $\{ u \in \mathcal{T}_x \mathcal{M} : \lVert u \rVert_x \le 1 \}$.
$y = \textrm{exp}_x(u)$, where $\textrm{exp}(\cdot)$ is the exponential map.
Clearly the first step is non-trivial, but that aside, would this approach yield uniform samples in the unit ball around $x$?
And then, is it possible to specify a class of manifolds for which it is possible to automate the first step? Perhaps for embedded submanifolds where we could sample uniformly from the ambient space and then somehow restrict to $\mathcal{T}_x \mathcal{M}$ (and apply the inverse Riemannian metric if different than the inherited one)?
Take for instance the Lorentz model of hyperbolic geometry: $$ \mathcal{H}^n = \{ x \in \mathbb{R}^{n+1} : \langle x, x \rangle_\mathcal{L} = -1, x_0 > 0 \},\quad\textrm{with}\ \langle x, y \rangle_\mathcal{L} = -x_0 y_0 + \sum_{i=1}^n x_i y_i. $$ The tangent space is $\mathcal{T}_x \mathcal{H}^n = \{ v \in \mathbb{R}^{n+1} : \langle x, v \rangle_\mathcal{L} = 0 \}$. How could one sample from the unit ball of tangent vectors uniformly?
Your proposed scheme will be ruined by curvature. For example, consider a 2-sphere of radius $1/\pi,$ so that the disc $$\exp(B_1) = \{\exp(u) : \|u\| \le 1\}$$ perfectly covers it. The outer region $$0.5 \le \|u\| \le 1$$ of the unit disc in the tangent space contains $3/4$ of the measure, so samples drawn using the uniform distribution on $B_1$ would hit this region $3/4$ of the time; but this corresponds to only half of the sphere in the Riemannian measure. (In your hyperbolic example you would find the opposite - the outer regions have much more Riemannian measure than "exponential" measure.)
If we use the exponential map on $B_1$ as a coordinate system $x$, this is reflected in the equation $$d\mu_g = \sqrt{\det g}\ d^nx. \tag 1$$ Here $d\mu_g$ is the Riemannian measure associated to $g,$ while $d^nx$ is the standard/Euclidean/uniform measure on the tangent/coordinate space. Thus, unless you happen to have $\det g=1$ in these coordinates, the measures are not the same.
However, equation $(1)$ should tell you how you can sample from $d\mu_g:$ if you have an expression for $g$ in normal coordinates, then $w=\sqrt{\det g}$ is a scalar function you can compute, so we have reduced the problem to that of sampling a weighted measure $d\mu_g = w\ d^nx,$ which you should be able to do programmatically (e.g. via the rejection method).
(One caveat: this only works if the injectivity radius at $x$ is at least $1$; i.e. if the exponential map is injective on $B_1.$ Otherwise, different points in $B_1$ will map to overlapping points on the manifold, so you'd have to restrict to a subset of $B_1$ making the exponential map cover each point in the unit ball exactly once.)