I don't understand how you arrived at your rotation matrix, so I'll just go by your description of the transformed system and the image, which I think I understand.
If we choose $\phi=0$, $O'$ lies on the $x$-axis, so the $z'$-axis should be the negative $x$-axis, the $y'$-axis should be the $z$-axis, and the $x'$-axis should correspondingly be the negative $y$-axis. That corresponds to the transformation matrix
$$\left(\begin{array}{ccc}0&-1&0\\0&0&1\\-1&0&0\end{array}\right)\;.$$
This is not what your expression for the rotation matrix yields for $\phi=0$. Since the $y'$-axis always points along the $z$-axis and only the $x'$-axis and the $z'$-axis are rotated into each other if we change $\phi$, we can immediately infer the general form of the rotation matrix as
$$\left(\begin{array}{ccc}\sin\phi&-\cos\phi&0\\0&0&1\\-\cos\phi&-\sin\phi&0\end{array}\right)\;,$$
where only the sign of the sine was still to be determined and can be fixed by considering $\phi=\pi/2$.
Let us assume you have points $\vec{a}_1$ and $\vec{b}_1$ in the first coordinate system, and $\vec{a}_2$ and $\vec{b}_2$ in the second coordinate system. Let us also assume that the coordinate systems are not scaled; only rotated and translated. If the points define a line segment, it is from $\vec{a}_1$ to $\vec{b}_1$ in the first coordinate system, and from $\vec{a}_2$ to $\vec{b}_2$ in the second coordinate system.
Let us choose the rotation to center on points $\vec{a}$ in their respective coordinate systems; that is, we translate by $-\vec{a}_1$ before rotation, and by $\vec{a}_1$ after the rotation, and we only rotate $\vec{b}$ around $\vec{a}$.
(You could choose different centers for the rotations, of course. A typical case for many points is to rotate the point cloud around its centroid. For two points, it is much simpler to just pick one point as the rotation origin.)
We can find the axis of the rotation and the sine of the rotation angle using vector cross product. In this case,
$$\hat{r}\sin\varphi = \vec{c} = \frac{\left ( \vec{b}_1 - \vec{a}_1 \right ) \times \left ( \vec{b}_2 - \vec{a}_2 \right )}{\left\lVert \vec{b}_1 - \vec{a}_1 \right\rVert \; \left\lVert \vec{b}_2 - \vec{a}_2 \right\rVert}$$
where $\hat{r}$ is an unit vector, $r_x^2 + r_y^2 + r_z^2 = 1$. First, calculate $\vec{c}$, and then its norm:
$$\lVert \vec{c} \rVert = \sqrt{c_x^2 + c_y^2 + c_z^2}$$
We'll also need the sign of $\cos\varphi$ to determine the correct quadrant for $\varphi$. $\cos\varphi$ has the same sign (and zeroes) as
$$d = \left ( \vec{b}_1 - \vec{a}_1 \right ) \cdot \left ( \vec{b}_2 - \vec{a}_2 \right )$$
Then,
$$\begin{cases}
\vec{r} = \frac{\vec{c}}{\lVert \vec{c} \rVert}, & \varphi = \operatorname{asin}\left(\frac{1}{\lVert \vec{c} \rVert}\right), & d \ge 0 \\
\vec{r} = -\frac{\vec{c}}{\lVert \vec{c} \rVert}, & \varphi = -\operatorname{asin}\left(\frac{1}{\lVert \vec{c} \rVert}\right), & d \lt 0
\end{cases}$$
This axis-angle representation corresponds to a versor, or unit quaternion, $\mathbf{q}$, that represents the same rotation:
$$\mathbf{q} = \cos\left(\frac{\varphi}{2}\right) + \sin\left(\frac{\varphi}{2}\right) r_x \mathbf{i} + \sin\left(\frac{\varphi}{2}\right) r_y \mathbf{j} + + \sin\left(\frac{\varphi}{2}\right) r_z \mathbf{k}$$
Typically, the components are referred to as $(w,x,y,z)$ or $(r,x,y,z)$ or $(a,b,c,d)$, in that order (scalar component of quaternion first); in other words,
$$\begin{cases}
q_w = q_a = q_0 = \cos\left(\frac{\varphi}{2}\right) \\
q_x = q_b = q_1 = \sin\left(\frac{\varphi}{2}\right) r_x \\
q_y = q_c = q_2 = \sin\left(\frac{\varphi}{2}\right) r_y \\
q_z = q_d = q_3 = \sin\left(\frac{\varphi}{2}\right) r_z
\end{cases}$$
The Wikipedia page for Conversion between quaternions and Euler angles shows how to convert this to Euler angles, if we rotate around $x$ axis first, then around $y$ axis, and finally around $z$ axis:
$$\begin{cases}
\theta_X = \operatorname{atan2}\left( 2 ( q_w q_x + q_y q_z ) , 1 − 2 ( q_x^2 + q_y^2 ) \right) \\
\theta_Y = \operatorname{asin}\left( 2 ( q_w q_y − q_z q_x ) \right) \\
\theta_Z = \operatorname{atan2}\left( 2 ( q_w q_z + q_x q_y ) , 1 − 2 ( q_y^2 + q_z^2 ) \right)
\end{cases}$$
For completeness, the rotation matrix using this order is (also shown on the Wikipedia Conversion between quaternions and Euler angles page):
$$\mathbf{R} = \left[\begin{matrix}
C_Z C_Y & -S_Z C_X + C_Z S_Y S_X & S_Z S_X + C_Z S_Y C_X \\
S_Z C_Y & C_Z C_X + S_Z S_Y S_X & -C_Z S_X + S_Z S_Y C_X \\
-S_Y & C_Y S_X & C_Y C_X
\end{matrix}\right]$$
where
$$\begin{matrix}
C_X = \cos \theta_X & S_X = \sin \theta_X \\
C_Y = \cos \theta_Y & S_Y = \sin \theta_Y \\
C_Z = \cos \theta_Z & S_Z = \sin \theta_Z \\
\end{matrix}$$
to keep the notation simple.
The transformation as defined above involves a translation before and after rotation. We can combine them to the usual post-rotation translation by applying the inverse rotation to the pre-rotation translation, and adding the post-rotation translation.
We translate by $-\vec{a}_1$ before the rotation, and $\vec{a}_1$ after the rotation, so we need to apply the inverse rotation to the pre-rotation translation to get the combined post-rotation translation:
$$\vec{t} = \vec{a}_1 + \mathbf{R}^{-1}(-\vec{a}_1) = \vec{a}_1 - \mathbf{R}^T \vec{a}_1$$
Note that $\mathbf{R}^{-1} = \mathbf{R}^{T}$ because pure rotation matrices like $\mathbf{R}$ here are orthogonal.
If $\vec{a}_1 = (x, y, z)$ and $\vec{t} = (t_x, t_y, t_z)$, then
$$\begin{cases}
t_x = x (1 - C_Z C_Y) - y (S_Z C_Y) + z S_Y \\
t_y = x (S_Z C_X - C_Z S_Y S_X) + y (1 - C_Z C_X - S_Z S_Y S_X) - z (C_Y S_X) \\
t_z = -x (S_Z S_X + C_Z S_Y C_X) + y (C_Z S_X - S_Z S_Y C_X) + z (1 - C_Y C_X)
\end{cases}$$
Separating the translation to pre- and post-rotation parts, and then combining the two, as above, is very often useful: it tends to make the problems simpler to solve, as you can basically freely choose the center of the rotation in the two coordinate systems, work out the rotation that way, and finally use the transpose of the rotation matrix to convert the pre-rotation translation to post-rotation translation. (Just remember to multiply the transpose of the rotation matrix $\mathbf{R}$ and the pre-rotation translation vector.)
Best Answer
This question is somewhat related to this question, with the main difference that you also have an translation.
Lets first make sure that the notation is clear:
A vector $\vec{v}$ can be represented in the initial Cartesian coordinate system as,
$$ \vec{v} = \begin{bmatrix} \vec{e}_x \\ \vec{e}_y \\ \vec{e}_z \end{bmatrix} \cdot \begin{bmatrix} v_x \\ v_y \\ v_z \end{bmatrix} = v_x \vec{e}_x + v_y \vec{e}_y + v_z \vec{e}_z, $$
where $\vec{e}_x$, $\vec{e}_y$ and $\vec{e}_z$ are the orthonormal basis vectors of the initial Cartesian coordinate system.
Similar, the vector $\vec{v}$ can also be represented in the new Cartesian coordinate system as,
$$ \vec{v} = \begin{bmatrix} \vec{e}_u \\ \vec{e}_v \\ \vec{e}_w \end{bmatrix} \cdot \begin{bmatrix} v_u \\ v_v \\ v_w \end{bmatrix} + \vec{o} = (v_u + o_u) \vec{e}_u + (v_v + o_v) \vec{e}_v + (v_w + o_w) \vec{e}_w, $$
where $\vec{o}$ is the vector connecting the initial and new origin and $\vec{e}_u$, $\vec{e}_v$ and $\vec{e}_w$ are the orthonormal basis vectors of the new Cartesian coordinate system, however these can be written as a linear combination of the basis vectors of the initial Cartesian coordinate system,
$$ \vec{e}_u = e_{ux} \vec{e}_x + e_{uy} \vec{e}_y + e_{uz} \vec{e}_z, $$
$$ \vec{e}_v = e_{vx} \vec{e}_x + e_{vy} \vec{e}_y + e_{vz} \vec{e}_z, $$
$$ \vec{e}_w = e_{wx} \vec{e}_x + e_{wy} \vec{e}_y + e_{wz} \vec{e}_z. $$
If $e_{wx}$, $e_{wy}$ and $e_{wz}$ are unknown you can find it by taking the cross product. Assuming that the initial basis vectors are right-handed, then $\vec{e}_w$ can be found with,
$$ \vec{e}_w = \vec{e}_u \times \vec{e}_v = (e_{uy}e_{vz} - e_{uz}e_{vy}) \vec{e}_x + (e_{uz}e_{vx} - e_{ux}e_{vz}) \vec{e}_y + (e_{ux}e_{vy} - e_{uz}e_{vx}) \vec{e}_z, $$
The rotation of the basis vectors can be calculated with the help of the rotation matrix and can be found with the help of dyadic products,
$$ R = \vec{e}_u\vec{e}_x + \vec{e}_v\vec{e}_y + \vec{e}_w\vec{e}_z, $$
this yields a second order tensor. In order to calculate the rotation you would have to take the dot product, but if you instead of dyadic products use vector direct products you can use a normal matrix-vector multiplication. In MATLAB this would can be done with,
such that,
$$ \begin{bmatrix} v_u \\ v_v \\ v_w \end{bmatrix} = R \begin{bmatrix} v_x \\ v_y \\ v_z \end{bmatrix} + \begin{bmatrix} o_u \\ o_v \\ o_w \end{bmatrix}. $$