How to construct a projection matrix to the midpoint of two points with rotation Picture of Scenario
I have an initial and terminal point.  I have a unit cube located at the origin.  I want to construct a projection matrix to move the cube to the origin, rotate, and scale it so that it becomes a rectangular prism that touches the initial and terminal point.  I can use the following functions to manipulate a matrix as well as anything else that could be useful:


*

*translate(mat: matrix, vec: vector) => matrix

*rotate(mat: matrix, angle: radians, axis: vector) => matrix

*scale(mat: mat, ratio: vector) => matrix

*radians(degrees: real) => real

*degrees(radians: real) => real
I can translate the cube to the midpoint(I, T), and scale it to the proper length easily with the help of distance(I, T).  However I am struggling to figure out how to do the rotation.  In my picture I drew yellow dots that are projected from the cube at the origin to the prism between the point I and T.
You could imagine a terminal point coming out of the paper.  How would you determine the axis of rotation?  How would you determine how many radians to rotate around the axis of rotation?
These posts seem to be the most relevant that I have found so far:


*

*Find rotation matrix to match points in parallel projection

*Finding an appropriate axis of rotation for two points such that they can be rotated and translated to overlay a given line
EDIT: Another Picture of the Scenario
Solution: https://stackoverflow.com/a/61404282/12750911
 A: Source
Spektre has answered this question in full detail here.
Cube to be Transformed

Summary of Solution
We want to compute transform matrix that will convert A,B into C,D.  While were at it we can also scale down the axes that are perpendicular to the vector AB to get a thickness.  We assume AB is on the x-axis.  The y-axis is vertical, and z-axis comes out of the screen.
We can then perform the transformation with the following sequence of functions.  The following functions are based on OpenGL Mathematics notation.
Define the points
A = vec3(-0.5, 0.0, 0.0);
B = vec3(+0.5, 0.0, 0.0);
C = p1;
D = p2;

Set M equal to the identity matrix
M = mat4(1.0);

Translate the unit cube to the midpoint of the vector CD.
center = 0.5 * (C + D);
M = translate(M, center);

Cross the rejections of AB and CD to get a normal.
p = B - A;
q = D - C;
n = cross(p, q);

We only rotate around the normal if it is not zero.
if (n != vec3()) {
    a = angle(normalize(p), normalize(q));
    M = rotate(M, a, n);
}

Finally, we scale the cube to be a rectangular prism that touchs the two points.
auto constexpr thickness = 0.05;
M = scale(M, vec3(0.5 * distance(C, D), thickness, thickness));

Apply the Matrix M
Now M * A = C.  This means we can multiple all the vertices on the cube by M to get the transformed cube.
End Result

