I posted this answer to a similar question on sci.math. I will transcribe the question and the summary of the solution below. For this problem, we don't need to compute $r$, just set it to $1$.
Least-Squares Conformal Multilinear Regression
Given $\{ P_j : 1 \le j \le m \}$ and $\{ Q_j : 1 \le j \le m \}$, two sets of
points, we want to find a conformal map, defined by a linear map, $M$,
and a vector, $R$, which maps one set of points to the other via
$$
Q = P M + R\tag{1}
$$
where we require that $M M^T = r^2 I$ and that the square residue
$$
\sum_{j=1}^m\left|P_jM+R-Q_j\right|^2\tag{2}
$$
is minimal. Note that $(1)$ requires that $P$ and $Q$ are row vectors.
Summary of the Method
To find the least squares solution to $P M + R = Q$ for a given set of
$\{ P_j \}$ and $\{ Q_j \}$, under the restriction that the map be conformal,
we first compute the centroids
$$
\overline{P}=\frac1m\sum_{j=1}^mP_j\qquad\text{and}\qquad
\overline{Q}=\frac1m\sum_{j=1}^mQ_j
$$
Next, compute the matrix
$$
\begin{align}
S
&=\sum_{j=1}^m\left(Q_j-\overline{Q}\right)^T\left(P_j-\overline{P}\right)\\
&=\sum_{j=1}^mQ_j^TP_j-m\overline{Q}^T\overline{P}
\end{align}
$$
Let the Singular Value Decomposition of $S$ be
$$
S=UDV^T
$$
Next compute $\{ c_k \}$ with
$$
\begin{align}
c_k
&=\sum_{j=1}^m\left[\left(P_j-\overline{P}\right)V\right]_k\left[\left(Q_j-\overline{Q}\right)U\right]_k\\
&=\sum_{j=1}^m\left[P_jV\right]_k\left[Q_jU\right]_k-m\left[\overline{P}V\right]_k\left[\overline{Q}U\right]_k
\end{align}
$$
and define
$$
a_k = \mathrm{sgn}( c_k )
$$
Let $I_k$ be the matrix with the $(k,k)$ element set to $1$ and all the
other elements set to $0$. Then calculate
$$
E=\sum_{k=1}^na_kI_k
$$
Compute the orthogonal matrix
$$
W=VEU^T
$$
If $\det(W) < 0$ but $\det(W) > 0$ is required, change the sign of the $a_k$
associated with the $c_k$ with the smallest absolute value.
If required, compute $r$ by
$$
r\sum_{j=1}^m\left|P_j-\overline{P}\right|^2=\sum_{j=1}^m\left\langle\left(P_j-\overline{P}\right)W,Q_j-\overline{Q}\right\rangle
$$
or equivalently
$$
r\left(\sum_{j=1}^m\left|P_j\right|^2-m\left|\overline{P}\right|^2\right)
=\sum_{j=1}^m\left\langle P_jW,Q_j\right\rangle-m\left\langle\overline{P}W,\overline{Q}\right\rangle
$$
Finally, we have the desired conformal map $Q = P M + R$ where
$$
M = r W
$$
and
$$
R = \overline{Q} - \overline{P} M
$$
More information, easier computation
Suppose you want to map $\{P_i\}_{i=1}^3$ to $\{Q_i\}_{i=1}^3$, and the distances between the $P_i$'s and $Q_i$'s are the same. Compute a fourth point by
$$
P_4=P_1+(P_2-P_1)\times(P_3-P_1)
$$
and
$$
Q_4=Q_1+(Q_2-Q_1)\times(Q_3-Q_1)
$$
Then create the matrix $P$ whose columns are $P_2-P_1$, $P_3-P_1$, and $P_4-P_1$.
Also create the matrix $Q$ whose columns are $Q_2-Q_1$, $Q_3-Q_1$, and $Q_4-Q_1$.
Then $x\mapsto QP^{-1}x+(Q_1-QP^{-1}P_1)$ maps the source points to the destination points.
Given the axes of the geometric body (hereafter: "body axes"), described with respect to some fixed global axes, you should be able to derive the corresponding rotation matrix that would rotate the global axes onto the body axes and vice versa.
In other words, if the body axes are $b_1, b_2, b_3$ and the global axes are $g_1, g_2, g_3$, then you want a rotation $R$ that does
$$R(g_1) = b_1, R(g_2) = b_2, R(g_3) = b_3$$
Or perhaps vice versa. Since inversion is the same as transposition for rotations, vice versa isn't hard to do either.
Now, again, let the body axes be described in terms of the global axes like so:
$$b_1 = B^{11} g_1 + B^{21} g_2 + B^{31} g_3$$
And similarly for $b_2, b_3$. Then, the matrix is basically computed for you. If you wrote out the above relations in matrix language, using the $g$ basis, you would have
$$R \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} B^{11} \\ B^{21} \\ B^{31}\end{bmatrix}$$
You should recognize now that this is just the first column of the matrix $R$. The columns of $R$ are just the body axes' unit vectors, expressed in terms of the global basis vectors.
Now that the matrix is known, you can use well-established algorithms to translate that into quaternions, if you should need to.
Best Answer
Suppose you want to find a rotation matrix $R$ that rotates unit vector $a$ onto unit vector $b$.
Proceed as follows:
Let $v = a \times b$
Let $s = \|v\|$ (sine of angle)
Let $c = a \cdot b$ (cosine of angle)
Then the rotation matrix R is given by: $$R = I + [v]_{\times} + [v]_{\times}^2\frac{1-c}{s^2},$$
where $[v]_{\times}$ is the skew-symmetric cross-product matrix of $v$, $$[v]_{\times} \stackrel{\rm def}{=} \begin{bmatrix} \,\,0 & \!-v_3 & \,\,\,v_2\\ \,\,\,v_3 & 0 & \!-v_1\\ \!-v_2 & \,\,v_1 &\,\,0 \end{bmatrix}.$$
The last part of the formula can be simplified to $$ \frac{1-c}{s^2} = \frac{1-c}{1-c^2} = \frac{1}{1+c}, $$ revealing that it is not applicable only for $\cos(\angle(a, b)) = -1$, i.e., if $a$ and $b$ point into exactly opposite directions.