If you are given a metric $d_g$ induced by a Riemannian metric $g$, then you can reobtain $g$ in the following way: I will wlog assume that $(M,g)$ is embedded in Euclidean space with $g$ the induced metric (this could be avoided, but the calculations might be messier).
Fix $p\in M$ and let $v\in T_pM$ be arbitrary. Then there exists a smooth path $\gamma: \mathbb R \to M$ such that $\gamma(0) = p$ and $\dot \gamma(0) = v$. Now we use the fact that
$\lim_{q\to p} \frac{d_g(p,q)}{|p-q|} = 1$
in the embedded setting (here $|\,\cdot\,|$ denotes Euclidean distance). This implies that
$\lim_{t\to 0} \frac{d_g(p,\gamma(t))}{t} = \lim_{t\to 0} \frac{d_g(p,\gamma(t))}{|p-\gamma(t)|} \frac{|p-\gamma(t)|}{t} = \lim_{t\to 0} \frac{|tv + O(t^2)|}{t} = |v| = \sqrt{g_p(v,v)}$
So if $d_g$ comes from a Riemannian metric, then we can recover the length of vectors via
$\Vert v \Vert_{g(p)} = \sqrt{g_p(v,v)} = \lim_{t\to 0} \frac{d_g(p,\gamma(t))}{t}$
Note that although the proof used a (hypothetical) isometric embedding of $M$ into Euclidean space, the resulting formula is purely intrinsic. Furthermore, it follows from the polarization identity that this determines $g_p$ uniquely (if it exists).
Thus, given an arbitrary metric $d$ on your manifold, you can in principle use the above expression (together with the polarization identity) and check whether it determines a smooth Riemannian metric $g$ on your manifold as you vary $p$. If it does determine a Riemannian metric $g$, then it would at least seem plausible if the associated metric $d_g$ where $=d$.
Add.: Indeed, suppose this is the case. I.e. let us assume that we can define $\Vert v\Vert_{g(p)}$ as indicated above, that this norm varies continuously on $M$ and that it is induced by an inner product. Then for an arbitrary curve $\gamma: [a,b]\to M$ and a partition $a=t_0 of $[a,b]$, we get
$\sum_{i=1}^n d(\gamma(t_i),\gamma(t_{i-1})) = \sum_{i=1}^n \frac{d(\gamma(t_i),\gamma(t_{i-1}))}{t_i - t_{i-1}} (t_i - t_{i-1})\to \int_a^b \Vert \dot \gamma \Vert_{(g\circ \gamma)(t)} \, dt$
as the partition becomes finer and finer. So the RHS approaches the length of the curve in the newly obtained Riemannian manifold $(M,g)$ in this case. On the other hand the expression on the LHS approaches the length of $\gamma$ in the metric space $(M,d)$ if we choose appropriate partitions! Hence $L^{(M,d)}(\gamma) = L^{(M,g)}(\gamma)$ for all smooth curves $\gamma$.
Let me summarize:
Let $M$ be a manifold and suppose $d$ is a metric on $M$ such that $d(p,q) = \inf_{\gamma} L(\gamma)$ for all $p,q \in M$, where the infimum is taken over all smooth paths running from $p$ to $q$. Then there is a Riemannian metric $g$ on $M$ such that $d_g = d$ if and only if
- $d^2$ is smooth, $\lim_{t\to 0} \frac{d(\gamma(t),p)}{t}$ exists for all curves running through $p$ at $t=0$ and
- the following "polarization identity" holds at every point $p\in M$: Given $v,w \in T_pM$ and any four curves $\gamma_v, \gamma_w, \gamma_{v+w}$ and $ \gamma_{v-w}$ in $M$ running through $p$ at $t=0$ with velocities $v,w,v+w$ and $ v-w$, respectively. Then \begin{align} \left[\frac{d}{dt}\Big|_{t=0} d(\gamma_{v+w}, p)\right]^2 &+ \left[\frac{d}{dt}\Big|_{t=0} d(\gamma_{v-w}, p)\right]^2 \\ & \qquad = 2\left[\frac{d}{dt}\Big|_{t=0} d(\gamma_{v}, p)\right]^2 +2\left[ \frac{d}{dt}\Big|_{t=0} d(\gamma_{w}, p)\right]^2 \end{align}