Let's say I have two boxes in 3D space. I know the Transformation of Box1 = $T^{world}_{box_0}$ and Box4 = $T^{box_0}_{box_4}$.
Our lovely designer want's me to write a tool that fills the the space between them for given amount of boxes (let's say 4 for demonstration) and each of them should transform such that $$T^{box_0}_{box_1} = T^{box_1}_{box_2} = T^{box_2}_{box_3} = T^{box_3}_{box_4}$$ and etc... So basically this should be my equation :
$$ T^{world}_{box_0} . T^{box_0}_{box_4} = T^{world}_{box_4} $$
$$ T^{world}_{box_0} . T^{box_0}_{box_1}. T^{box_1}_{box_2}. T^{box_2}_{box_3}. T^{box_3}_{box_4} = T^{world}_{box_4} $$
since all transformation between boxes should be equal : $$ T^{box_0}_{box_1} = T^{box_1}_{box_2} = T^{box_2}_{box_3} = T^{box_3}_{box_4} = T_R $$
$$ T^{world}_{box_0} . T_R. T_R. T_R. T_R = T^{world}_{box_4} $$
$$ T_R. T_R. T_R. T_R = T^{box_0}_{box_4} $$ $$ (T^{box_0}_{box_4})^{1/4} = T_R $$
So how do I calculate $(T^{box_0}_{box_4})^{1/n}$ ?
If your transformations $\mathbf{T}$ consist only of rotations and translations (and optionally global scaling), then the translation part is linear, and the rotation part easily solved via quaternions or bivectors.
See the Wikipedia Rotation matrix article for how to convert between a quaternion and a rotation matrix. (In a computer program, converting from a matrix to a quaternion or a bivector is best implemented in three ways, chosen at runtime depending which of the three diagonal entries has the largest magnitude.)
If $\mathbf{q}_0$ describes the initial orientation (rotation), and $\mathbf{q}_1$ the final orientation, then you can linearly interpolate between the two ($\lambda = 0 \dots 1$) to get an unnormalized quaternion, $$\mathbf{q}_u = (1 - \lambda) \mathbf{q}_0 + \lambda \mathbf{q}_1$$ which you can directly convert to a rotation matrix.
More generally, a $4 \times 4$ 3D transformation matrix $\mathbf{T}$ can be described as $$\mathbf{T} = \left[ \begin{matrix} \vec{u} & \vec{v} & \vec{w} & \vec{t} \\ 0 & 0 & 0 & 1 \\ \end{matrix} \right ]$$ where $\vec{u}$ is the $x$ axis unit vector after the transformation, $\vec{v}$ is the $y$ axis unit vector after the transformation, $\vec{w}$ is the $z$ axis unit vector after the transformation, and $\vec{t}$ is the translation.
You can find the corresponding orthonormal basis $[ \hat{u} ~ \hat{v} ~ \hat{w} ]$ from $\vec{u}$, $\vec{v}$, $\vec{w}$: $$\begin{aligned} \hat{u} &= \frac{\vec{u}}{\left\lVert \vec{u} \right\rVert} \\ \hat{v} &= \frac{\vec{\nu}}{\left\lVert \vec{\nu} \right\rVert}, \quad \vec{\nu} = \vec{v} - \left( \vec{v} \cdot \hat{u} \right) \hat{u} \\ \hat{w} &= \hat{u} \times \hat{v} \\ \end{aligned}$$ Then, $$\left\lbrace ~ \begin{aligned} \lambda_{uu} &= \vec{u} \cdot \hat{u} \\ \lambda_{vv} &= \vec{v} \cdot \hat{v} \\ \lambda_{vu} &= \vec{v} \cdot \hat{u} \\ \lambda_{ww} &= \vec{w} \cdot \hat{w} \\ \lambda_{wv} &= \vec{w} \cdot \hat{v} \\ \lambda_{wu} &= \vec{w} \cdot \hat{u} \\ \end{aligned} \right. \quad \Rightarrow \quad \left\lbrace ~ \begin{aligned} \vec{u} &= \lambda_{uu} \hat{u} \\ \vec{v} &= \lambda_{vv} \hat{v} + \lambda_{vu} \hat{u} \\ \vec{w} &= \lambda_{ww} \hat{w} + \lambda_{wv} \hat{v} + \lambda_{wu} \hat{u} \\ \end{aligned} \right. $$ This way, you can do the pure rotation as described earlier, linearly interpolating both the translation and the six non-orthogonality parameters $\lambda$, between the two endpoints.
This does have a preferred order of axes (which in above is $x$, $y$, and finally $z$). There are six possible orders – $xyz$, $xzy$, $yxz$, $yzx$, $zxy$, $zyx$ – and using the one that minimizes $\max\left( \lvert \lambda_{vu} \rvert, \lvert \lambda_{wv} \rvert, \lvert \lambda_{wu} \rvert \right)$ (i.e., minimizes the deviation from the basis set, choosing "the best" basis set as a result) might be a good idea. (I suspect it yields "the best" results, unless additional information is available on how the initial and final transformations differ from an orthonormal matrix; but I have no proof of this.)