The formula is premised on the wavelet being analytic, or being nonzero only over non-negative frequencies: ${\hat\psi} (\omega < 0) = 0$. (Note all wavelets also have ${\hat \psi (0)}=0$ per the admissibility criterion). From Daubechies et al,
$$
\begin{align}
\int_0^\infty W_f(a, b) a^{-1} da & = \frac{1}{2\pi} \int_{-\infty}^{\infty} \int_{0}^{\infty}{\hat f}(\xi) \overline{ {\hat \psi}(a\xi)} a^{-1} da d\xi \tag{1a} \\
&= \frac{1}{2\pi} \int_0^\infty \int_0^\infty {\hat f}(\xi) \overline{ {\hat \psi} (a\xi)} e^{jb\xi} a^{-1} da d\xi \tag{1b} \\
&= \int_0^\infty \overline{ {\hat \psi}(z)} \frac{dz}{z} \cdot \frac{1}{2\pi} \int_0^\infty {\hat f}(\xi) e^{jb\xi} d\xi \tag{1c} \\
\end{align}
$$
where a rather confusing change of variables was replaced with a clearer one, $z=a\xi,\ dz=\xi da$. Note that this c.o.v. is only valid if integrals share limits, which is enabled by $\psi$'s analyticity; otherwise, $z$ would integrate $-\infty $ to $\infty$, but $z$ includes $a$, which only spans $0$ to $\infty$. An identity was also utilized in splitting integrals.
Assuming $f$ is real, we have ${\hat f(\xi)} = \overline{\hat f (-\xi)}$, thus from $0$ to $\infty$,
$$
\frac{1}{4\pi} \Re e \left\{ \int_0^\infty {\hat f(\xi)} e^{jb\xi} d\xi \right\} = f(b),
$$
and setting $C_\psi = \int_0^\infty \overline{ {\hat \psi (z)}} z^{-1} dz$ thusly derives
$$
\begin{align}
f(b) &= \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty W_f(a, b) a^{-1} da \right\} \tag{2a} \\
&= \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty f \star \overline{\psi_a(b)} \frac{da}{a} \right\} \tag{2b} \\
&= \boxed{ \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty \left< f(t), \psi(t) \right> \frac{da}{a} \right\} } \tag{2c} \\
\end{align}
$$
Since $W_f$ sweeps $b$ over all $t$, above equivalently holds for $f(t)$.
The paper you linked uses $a^{-3/2}$ -- Indeed, as it assumes an L2 norm on the wavelet, whereas MATLAB uses L1; further info can be found here.
What about the $\delta$? -- The $\delta$ is from MATLAB writing the expression as a special case of the more general inverse, as it shows above near it, with $\left<f, g\right>$, which is derived here.
and the $2 \Re e$? -- Good question - perhaps MATLAB's $C_{\psi, \delta}$ is a bit different, but per above definition we clearly divide by two, not multiply. If in doubt, try both.
Interpretation: it helps to observe that $b$ in all derivations can be replaced with a constant - so the relation holds pointwise; if we integrate along $a$, i.e. a column of the 2D time-frequency representation, we recover the function at that point, $f(b)$.
Notice in the double-integral iCWT we use a dual wavelet to "undo" the inner product with the CWT wavelet, analogous to $e^{+j\omega t}$ in $\mathcal{F}^{-1}$; here, we don't use any inverting kernel - the transform undoes itself when integrated over all scales. This implies some symmetry property on the wavelet $\psi$ in encoding information across scales - though I'm unsure how this exactly works out; the most important is analyticity which permits this in the first place, so one may search there for an answer.
Discretized case with exponentially distributed $a$
From same page of linked paper, the discretized linear case reads
$$
f(b) \approx \frac{1}{2} \Re e \left\{ C_\psi^{-1} \sum_k W_f (a_k, b) \frac{(\Delta a)_k}{a_k} \right\}, \tag{3}
$$
with $(\Delta a)_k$ taking role of $da$. For the linear case, $(\Delta a)_k = \Delta a$, but not for the log case, as the increment size increases for greater $k$. That's the only distinction between log and linear.
In practice, however, it's better to directly use the discretized log-scaled transform - which derives as follows:
$$
\begin{align}
& z = e^a, \ dz = a da \\
& \frac{da}{a^2} = \frac{d \ln (z)}{z} \tag{4}
\end{align}
$$
Suppose $z(a) = 2^{a/32}$. Note that the log-plot of an exponential scale is linear, so the difference $\Delta(\ln (z))$ or $d(\ln (z))$ is constant, and equal to $\ln (z(1))da$. To better see this, let's derive directly for this example, using $n_v$ for "number of voices":
$$
z(a) = 2^{a/n_v} \rightarrow dz(a) = z \frac{\ln(2)}{n_v} da = \ln(a_0)da \tag{5}
$$
When discretizing, $da=1$, and we end up with
$$
\boxed{f[] \approx \frac{\ln (a_0)}{2 C_\psi} \Re e \left\{\sum_{k=1}^{K} \frac{1}{a_0^k} W_f (a_0^k, .) \right\} } \tag{6}
$$
same as Mallat Eq 4.67 (note his derivation is for double-integral, but same scaling argument applies here).