I think you can do this in a more "quaternionic" way if you want:
Quick review of quaternion and $\mathrm{SO}_3(\mathbb R)$:
We write $\mathbb H$ for the quaternions, and $h \mapsto \bar{h}$ for the quaternionic conjugation map, that is, if $h = h_0+h_1 \mathbf i + h_2\mathbf j+h_3\mathbf k\in\mathbb H$ then $\bar{h} = h_0-h_1 \mathbf i - h_2\mathbf j-h_3\mathbf k$. Let us also set $\Re(h) = h_0$ and $\Im(h) = h-h_0$. It is well-known that $\mathbb H$ is a normed division algebra over $\mathbb R$ with norm
$$
N(h) = h.\bar{h} = h_0^2+h_1^2+h_2^2+h_3^2.
$$
Thus in particular $\mathbb H$ is a 4-dimensional $\mathbb R$-vector space with inner product given by $\langle h,h' \rangle = \Re(h\overline{h'})$ for which the basis $\{1,\mathbf i,\mathbf j,\mathbf k\}$ is orthonormal. Let $\mathbb U = \{h \in \mathbb H: N(h)=h\bar{h}=1\}$ be the group of unit quaternions and let $\mathbb I = \{h \in \mathbb H: \bar{h}=-h\} = \text{span}_{\mathbb R}\{\mathbf i,\mathbf j,\mathbf k\}$ be the
subspace of purely imaginary quaternions. We identify $\mathbb R^3$ with the standard "dot product" with $\mathbb I$.
It is easy to check that $\mathbb U$ acts on $\mathbb I$ by conjugation, that is, if $u \in \mathbb U, h\in \mathbb I$ then $uhu^{-1} = uh\bar{u} \in \mathbb I$ and that this yields a homomorphism $\pi\colon\mathbb U \to \mathrm{O}(\mathbb I)$. Since $\mathbb U$ is connected, $\pi(\mathbb U)\subseteq \mathrm{SO}(\mathbb I) \cong \mathrm{SO}_3(\mathbb R)$. Now if $u = \cos(\phi)+\sin(\phi)u_{im}$ where $u_{im}\in \mathbb U\cap \mathbb I$, then $\cos(\phi) \in Z(\mathbb H)$, hence $u_{im}$ and $u$ commute with each other, and hence $u.u_{im}u^{-1}=u_{im}$. Thus $\pi(u) \in \mathrm{SO}_3(\mathbb R)$ is a rotation which preserves the line $\mathbb R.u_{im}$. It is then easy to check that $\pi(u)$ is a rotation by $2\phi$ about $u_{im}$.
Finally, if $a,b \in \mathbb I$, then
$$
a.b = \Re(a.b) + \Im(a.b) = -\langle a,b\rangle + \Im(a.b).
$$
since $\langle h_1,h_2\rangle = \Re(h_1\overline{h_2}) = -\Re(h_1h_2)$ if $h_2 \in \mathbb I$. Moreover, $a^2 = -a(-a) = -a.\bar{a} = -N(a)$ is real, hence $\Im(a.a)=0$. Thus $\Im(a.b)$ is an alternating bilinear map $\mathbb I\times \mathbb I\to \mathbb I$, which becomes the vector cross product via the standard identification.
Obtaining the rotator:
If we are given $v_1,v_2$ linearly independent vectors together with $v_1' = R(v_1)$ and $v_2' = R(v_2)$, where we view all of these as elements of $\mathbb I$.
Let $w_1= v_1-v_1'=(I-R)(v_1)$ and $w_2 = (I-R)(v_2)$. If $w_1,w_2$ are linearly dependent, say $w_2=\lambda w_1$, then $\lambda v_1-v_2 = R(\lambda v_1-v_2)$ and hence $\mathbb R.(\lambda v_1 -v_2)$ is the axis of the rotation $R$. Otherwise $\text{Im}(I-R) = \text{ker}(I-R)^{\perp}$ is spanned by $\{w_1,w_2\}$, and $\Im(w_1w_2) = a$ spans the axis of rotation of $R$.
Once we have a basis $\{a\}$ of the axis of rotation, it is straight-forward to find the angle of rotaiton: set $u = a/\|a\|$ so that $\{u\}$ is an orthonormal basis of the axis of rotation and let $c_1 = \langle v_1 ,u\rangle$. Then if $p_1 = v_1 - c_1u$ and $p_1' = v_1' - c_1 u$ we have $p_1,p'_1=R(p_1) \in \mathbb R. u^{\perp}$ and if
$$
\cos(\alpha) = \langle p_1,p_1'\rangle/\|p_1\|^2
$$
then $R$ is the rotation by $\alpha$ about $u$ and so the quaternion rotator is then $\cos(\alpha/2)+\sin(\alpha/2).u$.