Check out my answer to a somewhat similar problem here. You will need to generalize
$$ x = A_x m^2 + B_x mn + C_x n^2 + D_x m + E_x n + F_x $$
to
$$ x_{new} = A_x x^2 + B_x y^2 + C_x z^2 + D_x x y + E_x x z + F_x y z + G_x x + H_x y + I_x z + J_x $$
Also the same for $y$ and $z$. This will lead to needing to solve a 10x10 inverse. If it is linear, which will be confirmed if the second order coefficients are all zero or near zero, then you can put your results in matrix form like this:
$$
\left[
\begin{array}{c}
x_{new} \\
y_{new} \\
z_{new} \\
\end{array}
\right]
=
\left[
\begin{array}{cccc}
G_x & H_x & I_x & J_x \\
G_y & H_y & I_y & J_y \\
G_z & H_z & I_z & J_z \\
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
x \\
y \\
z \\
1 \\
\end{array}
\right]
$$
If you prefer, this will also work:
$$
\left[
\begin{array}{c}
x_{new} \\
y_{new} \\
z_{new} \\
\end{array}
\right]
=
\left[
\begin{array}{ccc}
G_x & H_x & I_x \\
G_y & H_y & I_y \\
G_z & H_z & I_z \\
\end{array}
\right]
\cdot
\left[
\begin{array}{c}
x \\
y \\
z \\
\end{array}
\right]
+
\left[
\begin{array}{c}
J_x \\
J_y \\
J_z \\
\end{array}
\right]
$$
If it is linear, recalculate with this equation and its $y$ and $z$ counterparts:
$$ x_{new} = G_x x + H_x y + I_x z + J_x $$
Then you will only have to solve a 4x4 inverse and will get a better fit.
In the 10x10 case, you will need at least 10 samples. In the 4x4 case, at least 4. The more you have, the better you can tell if your results are good. This technique produces a best fit solution.
It is also possible to generalize to higher powers, but of course this will lead to much larger inverses needing to be solved and more samples. You can think of this approach as a generalized multidimensional Taylor series best fit.
Hope this helps,
Ced