Since there is a $2$-$1$ homomorphism from $SU(2)\times SU(2)$ to $SO(4)$ there should be a way to recover the two $SU(2)$ matrices given an $SO(4)$ matrix.
I believe I could set this up as a system of equations using the map from above and solve for the coefficients of the $SU(2)$ matrices. However, I wonder if anyone has already done this or can point me to the formulas?
On
Here's a start to figuring out the formulas. The first thing is to understand the map $SU(2)\times SU(2)\to SO(4)$. Let $std$ denote the standard 2-dimensional complex representation of $SU(2)$. Then the external tensor product $V=std\otimes std$ is a natural 4-dimensional complex representation of $G=SU(2)\times SU(2)$. Concretely, it's given by the homomorphism $G\to GL(V)=GL_4(\mathbf{C})$ sending $(g,h)\in G$ to the "Kronecker product" $g\otimes h$, which is just the matrix of $g\otimes h$ with respect to the basis for $V=\mathbf{C}^2\otimes \mathbf{C}^2$ given by $e_1\otimes e_1, e_1\otimes e_2, e_2\otimes e_1, e_2\otimes e_2$, with $e_1,e_2$ the standard basis for $\mathbf{C}^2$. Now the standard representation $std$ actually lands in $SL_2(\mathbf{C})$, since $SU(2)\subset SL_2(\mathbf{C})$. But $SL_2(\mathbf{C})=Sp_2$, which means that the standard representation $std$ is self-dual, and thus $\mathbf{C}^2$ has an $SU(2)$-invariant, anti-symmetric bilinear form. Concretely this takes a pair of column vectors $v,w\in \mathbf{C}^2$ to the determinant of the matrix with columns $v,w$. Call this bilinear form $\omega:std\otimes std\to \mathbf{C}$. Then the external tensor product $\omega\otimes \omega:V\otimes V\to \mathbf{C}$ gives a $G$-invariant \textit{symmetric} bilinear form on $V$. In fact, this is an $SL_2(\mathbf{C})\times SL_2(\mathbf{C})$-invariant symmetric bilinear form on $V$. So we've shown that $SL_2(\mathbf{C})\times SL_2(\mathbf{C})$ actually lands in $SO(q)(\mathbf{C})$ for a certain 4-dimensional quadratic space $(V,q)$. (Here $q$ is the quadratic form associated to the sym. bilinear form $B=\omega\otimes\omega$.) By general Lie theory, this implies that the maximal compact $G$ of $SL_2(\mathbf{C})\times SL_2(\mathbf{C})$ lands in the maximal compact $SO(q)(\mathbf{R})$.
Of course, we really want a homomorphism $SU(2)\times SU(2)\to SO(4)$. To do this we must diagonalize $q$, over the complex number. It's not difficult to check that with respect to the tensor product basis for $V$ we were using above, the quadratic form $q$ is equal to (maybe a multiple of) $wz-xy$. That is, if $v\in V$ has coordinates $(w,x,y,z)$ with respect to the $e_i\otimes e_j$ basis, then $q(v)=wz-xy$. (Note that this is the "determinant" quadratic form!) You can diagonalize this form over $\mathbf{C}$ with the change of basis matrix
M=
[1 i 0 0]
[0 0 i -1]
[0 0 i 1]
[1 -i 0 0]'
unless I screwed up the arithmetic. So conjugating the representation $V$ by this matrix (thus not changing the isomorphism type of the repn), we obtain a homomorphism $SL_2(\mathbf{C})\times SL_2(\mathbf{C})\to SO(4,\mathbf{C})$ landing in the standard 4-dimensional special orthogonal group over the complex numbers. Now the maximal compact $G$ lands in the usual group $SO(4)$, and this gives the double cover you referred to in the question.
In summary, I think that the 2-1 map $SU(2)\times SU(2)\to SO(4)$ is given explicitly in matrices by taking the Kronecker product and then conjugating by $M$. If I didn't make an error, one should find in particular that the condition of being in $SU(2)$ forces a pair $(g,h)$ to have real matrix entries under this map! I haven't bothered to check this.
OK, so how does that actually let you find the 2 preimages of each rotation matrix in $SO(4)$? I guess it amounts to solving a whole bunch of nonlinear equations, so that's totally disgusting and impractical. To get a nice answer to your question, you need to analyze the structure of $SU(2)$ and of $SO(4)$. The latter is described thoroughly on the Wikipedia page for $SO(4)$. The former is quite simple and is well described in M. Artin's book Algebra. EDIT: I've deleted a bunch of stuff which Jason observed in a comment is totally bogus. Fortunately Robin gave another better answer which explains how to do this using quaternions.
On
In terms of SU(2) = unit quaternions, the answer can be found in the sections "Isoclinic decomposition" and "Relation to quaternions" of the Wikipedia article about SO(4). (This is equivalent to Robin Chapman's answer, I presume, but the formulas are a bit more explicit.)
One way to think of the 2-to-1 map from $SU(2)\times SU(2)$ is via quaternions. Regard $\mathbb{R}^4$ as $\mathbb{H}$, the set of quaternions, and $SU(2)$ as the group of unit (that is, norm $1$) quaternions. Then $SU(2)\times SU(2)$ acts on $\mathbf{H}$ by rotations in the following way: $$h\mapsto q_1 hq_2^{-1}$$ is a rotation of $\mathbb{R}^4=\mathbb{H}$ for $q_1$, $q_2\in SU(2)$. Then $(-q_1,-q_2)$ represents the same rotation as $(q_1,q_2)$.
Suppose we have a rotation $\phi$ of $\mathbb{H}$. If $\phi$ fixes $1$, it's effectively a rotation of the $\mathbb{R}^3$ spanned by $i$, $j$ and $k$ and is represented by the pair $(q,q)$ where $q=\cos(t/2)+h\sin(t/2)$ where $h=ai+bj+ck$ is a unit vector along the axis of rotation, and $t$ is the angle of rotation.
In general, let $q'$ be the unit quaternion defined by $\phi(1)$. Then $\psi:h\mapsto q'^{-1}\phi(h)$ is a rotation fixing $1$ so that one can find a quaternion $q$ such that $(q,q)$ represents $\psi$. Then $(q'q,q)$ represents $\phi$.