When two projections in a C*-algebra are "almost" Murray-von Neumann equivalent, they are equivalent
One can prove that $p,q$ are Murray -von-Neumann equivalent if they fulfill $$ \|aa^*-p\|<\frac 12,\quad \|a^*a-q\|<\frac 12. $$ The proof is standard (for example it is outlined in the comment of @ougao), I was just lucky to find $\|(1/2-p)^{-1}\|=2$.
We first assume that the $C^*$-algebra $A$ is unital.
First note that $(1-2p)^2=1$, hence $(1-2p)$ is its own inverse and $\|1-2p\|=1$. It follows that the inverse of $\frac 12-p$ is $2-4p$ and that $\|2-4p\|=2$. If $\|aa^*-p\|<\frac 12$, then $\|aa^*\|< \frac 32$ and $Spec(aa^*)\subset [0,\frac 32[\subset\Bbb{R}$. Moreover, since $\frac 12-p$ is invertible and $\|p-aa^*\|<\frac{1}{\|(\frac 12-p)^{-1}\|}$, we have that $\frac 12-p+p-aa^*=\frac 12-aa^*$ is invertible. This means that the spectrum of $aa^*$ has a gap at $\frac 12$, i.e., an open interval around $\frac 12$ which does not intersect $Spec(aa^*)$, since $Spec(aa^*)$ is closed.
Hence there exists an $r\in]0,\frac 12[$ such that $$ Spec(aa^*)\subset B_r(0)\cup B_r(1). $$ Similarly we find that $$ Spec(a^*a)\subset B_{r_1}(0)\cup B_{r_1}(1), $$ and we can assume that $r=r_1$.
We define on $D:=(B_{r}(0)\cup B_{r}(1))\cap\Bbb{R}$ the continuous function $$ f(x)=\left\{ \begin{array}{lr} 0& \text{if $x\in B_r(0)$}\\ 1/\sqrt{x}& \text{if $x\in B_r(1)$.} \end{array} \right. $$ Note that on $D$ the function $xf^2(x)$ is a projection with a distance less than $r$ to the identity function, i.e.,
$x^2f^4(x)=xf^2(x)$.
$|xf^2(x)-x|\le r$.
Then we apply functional calculus and set $v:= f(aa^*)a$. We have that $vv^*=aa^*f^2(aa^*)$ and moreover, by 1., $$ (vv^*)^2=(aa^*)^2f^4(aa^*)=aa^*f^2(aa^*). $$ Note that for any continuous function $f$ on $D$ we have $a^*f(aa^*)a=a^*af(a^*a)$ (this holds for polynomials and $f$ is the limit of polynomials).
Therefore we have $v^*v=a^*f^2(aa^*)a=a^*af^2(a^*a)$ and so, by 1., $$ (v^*v)^2=(a^*a)^2f^4(a^*a)=a^*af^2(a^*a). $$ So $vv^*$ and $v^*v$ are equivalent projections (they are the range projections of $aa^*$ and $a^*a$ respectively) and we compute $$ \|vv^*-p\|\le \|vv^*-aa^*\|+\|aa^*-p\|< r+\frac 12<1, $$ where $\|vv^*-aa^*\|\le r$ follows from 2.
Similarly $$ \|v^*v-q\|\le \|v^*v-a^*a\|+\|a^*a-q\|< r+\frac 12<1, $$ where $\|v^*v-a^*a\|\le r$ follows from 2.
Since $A$ is unital, there are unitary elements $u,s$ with $uvv^*u^*=p$ and $v^*v=sqs^*$. Then $(uvs)^*uvs=q$ and $uvs(uvs)^*=p$.
Finally we prove the general case. If $A$ is non-unital, embed it into its unitalization $A^\sim$, and there $p$ and $q$ are Murray-von-Neumann equivalent, i.e., there exists $b\in A^\sim$, such that $bb^*=p$ and $b^*b=q$. But then $b=pb$, since $\|(b-pb)(b-pb)^*\|=0$, and so, since $A$ is an ideal in $A^\sim$ we have that $b\in A$, hence $p,q$ are Murray-von-Neumann equivalent in $A$.
Here is my own attempt on the problem, which turns out to be essentially San’s argument above (I only realized this when writing my response). I would therefore like to focus on demystifying the continuous functional calculus that takes place in the argument and on making portions of it more transparent. All credit is due to San, of course, for posting the first solution. My only goal here is to explain my thought processes to show that the problem yields to a systematic attack and is therefore not exceedingly difficult.
We will prove the stronger version of the problem as presented in San’s answer. We will also relax the condition that $ A $ is unital.
As all projections have norm $ 1 $, playing around with the Triangle Inequality yields $$ \| a a^{*} \| = \| a^{*} a \| < \frac{3}{2}. $$ As $ a a^{*},a^{*} a \in A_{\geq 0} $, we thus have $$ \sigma(a a^{*}),\sigma(a^{*} a) \subseteq \left[ 0,\dfrac{3}{2} \right). $$ Sam has already shown that $ \dfrac{1}{2} \notin \sigma(a a^{*}) $ and $ \dfrac{1}{2} \notin \sigma(a^{*} a) $, under the extra assumption that $ A $ is unital. However, we do not have to require that $ A $ is unital, for once we have argued in $ A^{\sim} $ to get $ \dfrac{1}{2} \notin \sigma(a a^{*}) $ and $ \dfrac{1}{2} \notin \sigma(a^{*} a) $, we can switch back to $ A $. The remainder of the argument does not rely on the existence of a unit. However, when applying the continuous functional calculus in the absence of a unit, any function introduced must satisfy the property that its value at $ 0 $ is $ 0 $.
By a standard Banach-algebra result, we get $$ \sigma(a a^{*}) \cup \{ 0 \} = \sigma(a^{*} a) \cup \{ 0 \}. $$ We can thus find an $ r \in \left[ 0,\dfrac{1}{2} \right) $ and an $ s \in \left( \dfrac{1}{2},\dfrac{3}{2} \right) $ such that $$ \sigma(a a^{*}),\sigma(a^{*} a) \subseteq D \stackrel{\text{df}}{=} [0,r] \cup \left[ s,\dfrac{3}{2} \right]. $$ Then the function $ f: D \to \Bbb{R} $ defined by $$ \forall x \in D: \quad f(x) \stackrel{\text{df}}{=} \begin{cases} 0 & \text{if $ x \in [0,r] $}, \\ 1 & \text{if $ x \in \left[ s,\dfrac{3}{2} \right] $}, \end{cases} $$ is continuous on both $ \sigma(a a^{*}) $ and $ \sigma(a^{*} a) $. Applying the continuous functional calculus, we have well-defined projections $ f(a a^{*}) $ and $ f(a^{*} a) $, because $ f = f^{*} = f^{2} $. These are called spectral projections, and they exist for self-adjoint elements that have a gap in their spectrum.
Now, $ f(a a^{*}) $ and $ f(a^{*} a) $ are Murray-von Neumann equivalent. To prove this, we want to exploit the fact that $ f(a a^{*}) a = a f(a^{*} a) $, as already explained in this other post. Suppose that there is a positive continuous function $ g: D \to \Bbb{R} $ such that $ f(x) = x [g(x)]^{2} $ for all $ x \in D $. If we let $$ b \stackrel{\text{df}}{=} g(a a^{*}) a, $$ then \begin{align} b b^{*} & = [g(a a^{*}) a] [a^{*} \overline{g}(a a^{*})] \\ & = [g(a a^{*}) a] [a^{*} g(a a^{*})] \quad (\text{As $ g $ is positive.}) \\ & = g(a a^{*}) (a a^{*}) g(a a^{*}) \\ & = f(a a^{*}). \quad (\text{As $ g(x) x g(x) = f(x) $ for all $ x \in D $.}) \end{align} Similarly, \begin{align} b^{*} b & = [a^{*} \overline{g}(a a^{*})] [g(a a^{*}) a] \\ & = [a^{*} g(a a^{*})] [g(a a^{*}) a] \quad (\text{As $ g $ is positive.}) \\ & = [g(a^{*} a) a^{*}] [a g(a^{*} a)] \\ & = g(a^{*} a) (a^{*} a) g(a^{*} a) \\ & = f(a^{*} a). \quad (\text{As $ g(x) x g(x) = f(x) $ for all $ x \in D $.}) \end{align} Of course, we still have to find $ g $, and the most natural (in fact, only) choice is $$ \forall x \in D: \quad g(x) \stackrel{\text{df}}{=} \begin{cases} 0 & \text{if $ x \in [0,r] $}, \\ \dfrac{1}{\sqrt{x}} & \text{if $ x \in \left[ s,\dfrac{3}{2} \right] $}, \end{cases} $$ which is the same function exhibited in San’s answer.
Observe next that $ \| p - f(a a^{*}) \| < 1 $. This is because \begin{align} \| p - f(a a^{*}) \| & \leq \| p - a a^{*} \| + \| a a^{*} - f(a a^{*}) \| \\ & \leq \| p - a a^{*} \| + \| \text{id}|_{D} - f \|_{\infty} \quad (\text{By the continuous functional calculus.}) \\ & < \frac{1}{2} + \frac{1}{2} \quad (\text{Compare the graphs of $ \text{id} $ and $ f $.}) \\ & = 1. \end{align} Similarly, $ \| q - f(a^{*} a) \| < 1 $.
Projections that are within distance $ 1 $ of each other are homotopy-equivalent (this is folklore). It is also not hard to show that homotopy equivalence implies Murray-von Neumann equivalence. Therefore, $$ p \sim f(a a^{*}) \sim f(a^{*} a) \sim q. $$