Define $u := \frac{1}{d} (x_2 - x_1, y_2 - y_1, z_2 - z_1)^T$, in which $d$ is the distance between your two points. This means, $u$ is the unit vector in the direction of the line segment. You are looking for a rotation that maps $e_x = (1, 0, 0)$ to $u$.
The solution is not unique, as there are (infinitely) many rotations that achieve this, but you can single out the one that goes the most "direct" route - which is the one rotating about the axis perpendicular to both $u$ and $e_x$.
If you define $v := u \times e_x$ and $w := v \times u$ (using the cross product), the rotation matrix $R = (u, v, w)$ achieves exactly this (here, the vectors are the columns of $R$, which is a 3x3 matrix).
This matrix allows you to perform the desired rotation, and if all you need to do is rotate stuff, you should stick with it. However, if you really need to extract the rotation angles about the individual axes, it gets a bit ugly, bu see for example here for a way to do it.
As a first step, I'd move $P_1$ to the origin, so that the points become $P_1=(0,0,0)$, $P_2=(x_2-x_1, y_2-y_1, z_2-z_1)$, $P_3=(x_3-x_1, y_3-y_1, z_3-z_1)$. From there it's just a question of applying a rotation matrix.
A rotation matrix $\mathbf{R}$ which preserves all distances and shapes etc is
given by
$$
\mathbf{R}=\left(
\begin{array}
[c]{ccc}%
e_{1}^{x} & e_{2}^{x} & e_{3}^{x}\\
e_{1}^{y} & e_{2}^{y} & e_{3}^{y}\\
e_{1}^{z} & e_{2}^{z} & e_{3}^{z}%
\end{array}
\right)
$$
where the vectors $\left( e_{1}^{\alpha},e_{2}^{\alpha},e_{3}^{\alpha
}\right) $ for $\alpha\in\left\{ x,y,z\right\} $ are orthogonal and of unit
length, i.e.
$$
\sum_{i=1}^{3}e_{i}^{\alpha}e_{i}^{\beta}=\left\{
\begin{array}
[c]{ccc}%
1 & & \alpha=\beta\\
0 & & \alpha\neq\beta
\end{array}
\right. .%
$$
Applying $\mathbf{R}$, your rotated points $P_{i}^{\prime}=\left(
x_{i}^{\prime},y_{i}^{\prime},z_{i}^{\prime}\right) $ are given by
$$
\left(
\begin{array}
[c]{c}%
x_{i}^{\prime}\\
y_{i}^{\prime}\\
z_{i}^{\prime}%
\end{array}
\right) =\mathbf{R}\left(
\begin{array}
[c]{c}%
x_{i}\\
y_{i}\\
z_{i}%
\end{array}
\right)
$$
I suggest making the top row of $\mathbf{R}$ the unit
vector in the direction from $P_{1}$ to $P_{2}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{x}\\
e_{2}^{x}\\
e_{3}^{x}%
\end{array}
\right) =\frac{1}{\sqrt{\left( x_{1}-x_{2}\right) ^{2}+\left( y_{1}%
-y_{2}\right) ^{2}+\left( z_{1}-z_{2}\right) ^{2}}}\left(
\begin{array}
[c]{c}%
x_{1}-x_{2}\\
y_{1}-y_{2}\\
z_{1}-z_{2}%
\end{array}
\right)
$$
which would mean that the line from $P_{1}$ to $P_{2}$ gets mapped into the
$x$-axis.
Working out what you could use for the second and third rows of $\mathbf{R}$ is now just a question of solving some simple linear equations, to make sure that the line from $P_{1}$ to $P_{3}$ has no z component.
To solve for $\left( e_{1}^{z},e_{2}^{z},e_{3}^{z}\right)$, the vector
$e_{i}^{z}$ must have a zero dot-product with both $e_{i}^{x}$ and $\left(
x_{3}-x_{1},y_{3}-y_{1},z_{3}-z_{1}\right) $, so the equations are
$$
\begin{align*}
e_{1}^{x}e_{1}^{z}+e_{2}^{x}e_{2}^{z}+e_{3}^{x}e_{3}^{z} & =0\\
\left( x_{3}-x_{1}\right) e_{1}^{z}+\left( y_{3}-y_{1}\right) e_{2}%
^{z}+\left( z_{3}-z_{1}\right) e_{3}^{z} & =0
\end{align*}
$$
Hence
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{z}\\
e_{2}^{z}\\
e_{3}^{z}%
\end{array}
\right) =\lambda_{z}\left(
\begin{array}
[c]{c}%
\left( z_{3}-z_{1}\right) e_{2}^{x}-\left( y_{3}-y_{1}\right) e_{3}^{x}\\
\left( x_{3}-x_{1}\right) e_{3}^{x}-\left( z_{3}-z_{1}\right) e_{1}^{x}\\
\left( y_{3}-y_{1}\right) e_{1}^{x}-\left( x_{3}-x_{1}\right) e_{2}^{x}%
\end{array}
\right)
$$
for some $\lambda_z$, which should be determined so that $e_{i}^{z}$ is a vector
of unit length. Finally $e_{i}^{y}$ can be determined as the vector which is
perpendicular to both $e_{i}^{x}$ and $e_{i}^{z}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{y}\\
e_{2}^{y}\\
e_{3}^{y}%
\end{array}
\right) =\lambda_{y}\left(
\begin{array}
[c]{c}%
e_{3}^{z}e_{2}^{x}-e_{2}^{z}e_{3}^{x}\\
e_{1}^{z}e_{3}^{x}-e_{3}^{z}e_{1}^{x}\\
e_{2}^{z}e_{1}^{x}-e_{1}^{z}e_{2}^{x}%
\end{array}
\right)
$$
where again $\lambda_{y}$ is determined so that $e_{i}^{y}$ is a vector of
unit length.
Best Answer
Because OP is using OpenSCAD, the claims of the question aren't entirely correct: orientation can be changed not only by rotations about axes, but by other transformations as well.
I'm going to assume that the cylinder has been rotated/translated/whatever by the OP, who can surely do this, so that one end is at the origin, and the other end is at $(0,0,1)$. (Maybe OpenSCAD aligns it along $(0,1,0)$, but I don't think that adapting to this difference will present any challenges). I'm not going to write correct openSCAD code because I haven't used it in about a year, but I'll write something that should get you going in the right direction.
Here goes.
There's a
rotate(v, a)function in openSCAD that rotates around a vector $v$ by angle $a$. We'll use that.$$ \newcommand{zv}{\mathbf{\hat{z}}} \newcommand{yv}{\mathbf{\hat{y}}} \newcommand{sv}{\mathbf{{s}}} $$
Let's let
zvec(in code) or $\zv$ (in math) mean $[0,0,1]$, the initial axis of the cylinder. ThenCompute $v = Q - P$, the vector from the point $P$ to the point $Q$. If $v$ is zero, your input was bad (a length-zero cylinder makes no sense), and you'll have to handle this exceptional case.
Compute $u = \frac{v}{ \| v\|}$. In openSCAD, you do something like
v = (1/norm(v))*v, I believe. If $u$ turns out to be exactly $-\zv$, then this is another special case; because the cylinder is already aligned with the $z$-axis, all you need to do is rotate it about the $x$-axis (for instance) by 180 degrees to get it where we want; you can then skip to step 5.)Compute $w = \frac12 (\zv + u)$, or in code,
w = (0.5) * (u + zvec), and then $t = \frac{w}{\|w\|}$, ort = (1/norm(w)) * w. At this point, $t$ is a unit vector in the plane of $v$ and $\zv$, and is exactly halfway between those two vectors.Rotate about $t$ by 180 degrees to take the cylinder, which used to point along the $\zv$ direction, and make it point along the $v$ direction instead.
Scale everything by $\| v \| $, i.e.,
norm(v), which is the distance from $P$ to $Q$.Translate by $P$ to move the origin (one end of the cylinder) to $P$, and the other end (which is at the tip of $v = Q - P$) to the point $Q$.
All of that stuff works mathematically, but I should caution you that the special case in step 2 is problematic: testing equality of floating-point numbers is fraught with peril. If the vector $u$ turns out to be very near to $-\zv$, then the remaining computations may be numerically a bit unstable, and you really don't want that. On the other hand, if you're just making a model of your back porch to decide where to put the picnic table, this nasty case is unlikely to affect you, and you can stop reading.
To continue: to make this whole approach a little more numerically robust, I'll modify things from step 2 onwards: if the vector $u$ is too close to being $\pm \zv$, I'm going to rotate the cylinder to point along the $y$-axis first. I'll use $\yv$ to denote the vector $(0, 1, 0)$. Here's the revised version, math only. I leave it to you to translate to openSCAD.
Compute $v = Q - P$, the vector from the point $P$ to the point $Q$. If $v$ is zero, your input was bad (a length-zero cylinder makes no sense), and you'll have to handle this exceptional case. In fact, if $\| v \|$ is smaller than, say, $10^{-5}$, you should probably treat this as an error case. (I can't recall whether openSCAD uses IEEE floats or IEEE doubles or something else to represent things, but using $10^{-5}$ as a threshold should keep you completely safe regardless. You could probably get away with using $10^{-10}$.)
Compute $u = \frac{v}{ \| v\|}$. If the absolute value of the $z$-component of $u$ is no more than $0.8$, let $\sv = \zv$; if the absolute value of the $z$-component is more than $0.8$, then rotate the cylinder by $-90$ degrees about the $x$-axis, so that it now lies along the positive $y$ axis, and let $\sv = \yv$. In both cases, $\sv$ refers to "the starting orientation of the cylinder".
Compute $w = \frac12 (\sv + u)$, and $t = \frac{w}{\|w\|}$, At this point, $t$ is a unit vector in the plane of $v$ and $\sv$, and is exactly halfway between those two vectors.
Rotate about $t$ by 180 degrees to take the cylinder, which used to point along the $\sv$ direction, and make it point along the $v$ direction instead.
Scale everything by $d = \| v \| $, the distance from $P$ to $Q$. The cylinder now lies along the $PQ$ direction, with one end at the origin, and the other end $d$ units from the origin, i.e., it has the same length and orientation as the line segment $\overline{PQ}$.
Translate the cylinder by $P$ to move the origin (one end of the cylinder) to $P$, and the other end (which is at the tip of $v = Q - P$) to the point $Q$.
There's another method that involves constructing a $4 \times 4$ matrix directly, and using openSCAD's
multmatrixoperation, but it's probably a little messier to write and to understand, so unless you're unhappy with the solution above, I won't bother writing it out.