lets assume your 2 points are defining the cube's internal diagonal u (blue).
From these 2 points we already know 2 vertexes of cube and their average is center of the cube p.
Now we need to select plane in which both u,v diagonals will be by using any direction not (anti)parallel to diagonal u and use it as normal of the plane n (green). In such cases I usually use either (1.0,0.0,0.0) or (0.0,1.0,0.0) depending on if the first is parallel or not (that can be decided with dot product). You can also use camera forward and right directions so the cube aligns the same way regardless on camera orientation.
Now to obtain the red diagonal we just rotate the blue one by ~70.529deg = 2*atan(1/sqrt(2)) around align vector n and center p.
Finally to obtain last missing 2 diagonals (U,V in code) we just rotate the first two (u,v) by 90deg around n' and center p. The n' is just direction in the midle u and -v so n'=u-v.
Here small C++ code for this using my GLSL style vector math (you can use any other like GLM or even own you just need v+v,v-v,v*c,v/c,dot,cross,length) also note angles a0,a1 are in radians:
//---------------------------------------------------------------------------
vec3 A=vec3(-0.7,-0.8,+0.5), // input points
B=vec3(+0.5,+0.6,-0.4);
vec3 pnt[8]; // output cube points
//---------------------------------------------------------------------------
// return rotated p around axis (center is (0,0,0))
vec3 rotate(float ang,vec3 axis,vec3 p)
{
float c=cos(ang),s=sin(ang);// precompute goniometrics
vec3 u,v,w; // bazis vectors
float x,y,z,xx,yy;
w=normalize(axis); // w is rotation axis
u=vec3(1.0,0.0,0.0); // u = any unit align vector not parallel to w
if (fabs(dot(u,w))>0.75) u=vec3(0.0,1.0,0.0);
u=cross(u,w); // make u paralel to w and align vector
v=cross(u,w); // make v paralel to u and w
// convert to basis vectors local coordinates
xx=dot(p,u);
yy=dot(p,v);
z =dot(p,w);
// rotate XY plane
x=xx*c-yy*s;
y=xx*s+yy*c;
// convert back to global coordinates
p=(x*u)+(y*v)+(z*w);
return p;
}
//---------------------------------------------------------------------------
// compute pnt[] from A,B
void diagonal2cube()
{
const float a0=1.230959417340774682134929178248; // 70.5288 deg
const float a1=1.5707963267948966192313216916398; // 90.0000 deg
float a;
vec3 p,u,v,U,V,n;
p=0.5*(A+B); // center of cube
u=A-p; // half diagonal from center to pnt[0]
a=length(u); // half size of cube's diagonal
u=normalize(u); // normalize u
n=vec3(1.0,0.0,0.0);// n = any unit align vector not parallel to u
if (fabs(dot(u,n))>0.75) n=vec3(0.0,1.0,0.0);
n=cross(u,n); // make n paralel to u and align vector
u*=a; // rescale u back to original size
v=rotate(a0,n,u); // v is u rotated by 70.5288 deg around n
n=u-v; // n is axis perpendicular to n and in the same plane as u,v
U=rotate(a1,n,u); // U is u rotated by 90.0 deg around n
V=rotate(a1,n,v); // V is v rotated by 90.0 deg around n
pnt[0]=p+u; // construct the cube vertexes
pnt[6]=p-u;
pnt[1]=p+v;
pnt[7]=p-v;
pnt[3]=p+U;
pnt[5]=p-U;
pnt[2]=p+V;
pnt[4]=p-V;
}
//---------------------------------------------------------------------------
Here preview of generated cube using OpenGL:
