Dimension of square rotated in 3D from projection on 2D

Question

I would like to track four points arranged as a square through 3D space. The square can rotate freely in the three dimensions. The example below left shows the square(blue) rotated 30 degrees about the x-axis, 60 degrees about the y-axis and 20-degrees about the z-axis. On the right is what the camera sees. The orange figure on left and right is the projection on xy plane.

To track the distance from camera, I need to know the dimensions of the blue square. We know everything about the orange figure, all angles and lengths. My problem is I can't make the connection to the blue square. I would appreciate some idea on how to calculate the blue square dimensions from knowing the orange angles and lengths.

Here is the Python 2.7 code to generate the figure:

(the code has been updated to implement the answer)

import numpy as np
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def show_the_rotated_square():
    #create the square in 2d with center at 0,0
    square_2D = [[-1, 1],[1, 1], [1, -1], [-1, -1]]

    #create a square in 3d so the square can be rotated in all three dimensions
    square_3D = [add_z(point) for point in square_2D]

    #create the three axis which the square will be rotated about
    axis = [[1, 0, 0],[0,1,0],[0,0,1]]

    X_ROTATE = 30
    Y_ROTATE = 60
    Z_ROTATE = 20

    #create a list of angles the square will be rotated for each axis
    theta = [np.radians(X_ROTATE),np.radians(Y_ROTATE),np.radians(Z_ROTATE)] 

    #rotate the 3d square about all three axis
    square_rotated_3D = rotate_point_list(square_3D,axis,theta)

    #create 2d version by removing z element
    square_rotated_2D = [remove_z(point) for point in square_rotated_3D]

#   print(angle_between_points(square_rotated_2D[0],square_rotated_2D[1],square_rotated_2D[2]))

    plot_3D(square_rotated_3D,121,"3D View, square rotated " + str(X_ROTATE) + " " + str(Y_ROTATE) + " " + str(Z_ROTATE),show_projection = True)
    plot_2D(square_rotated_2D,122,"Projection on XY plane (what camera sees)")

    solve(square_rotated_2D[0],square_rotated_2D[1],square_rotated_2D[3])
    solve(square_rotated_2D[1],square_rotated_2D[2],square_rotated_2D[0])
    solve(square_rotated_2D[2],square_rotated_2D[3],square_rotated_2D[1])
    solve(square_rotated_2D[3],square_rotated_2D[0],square_rotated_2D[2])

    plt.show()

#see math.stackexchange for discussion of solution
#https://math.stackexchange.com/questions/2879056/dimension-of-square-rotated-in-3d-from-projection-on-2d

def solve(pa,pb,pd):

    pb = pb - pa
    pd = pd - pa
    pa = pa - pa

    xb = pb[0]
    yb = pb[1]
    xd = pd[0]
    yd = pd[1]

    v = -(xb * xd + yb * yd)
    u = (xd**2 + yd**2) - (xb**2 + yb**2)

    a = 1
    b = -u
    c = -(v * v)

    #calculate discriminant
    d = b**2 - 4*a*c

    sol1 = (-b - math.sqrt(d)) / (2*a)
    sol2 = (-b + math.sqrt(d)) / (2*a)

    if (sol1>0):
        calc_square_side_length(sol1,v,xb,yb,xd,yd)
    if (sol2>0):
        calc_square_side_length(sol2,v,xb,yb,xd,yd)

def calc_square_side_length(s,v,xb,yb,xd,yd):
    b = math.sqrt(s)
    d = v / b

    ab = math.sqrt(xb**2 + yb**2)
    ad = math.sqrt(xd**2 + yd**2)

    ablen = math.sqrt(ab**2 + b**2)
    adlen = math.sqrt(ad**2 + d**2)

    print(ablen,adlen) # ablen and adlen should be same and be the side lengh of original 2d square, which is 2.0

#rotation_matrix code courtesy of:
#https://stackoverflow.com/questions/6802577/rotation-of-3d-vector/12261243


def rotation_matrix(axis, theta):
    """
    Return the rotation matrix associated with counterclockwise rotation about
    the given axis by theta radians.
    """
    axis = np.asarray(axis)
    axis = axis / math.sqrt(np.dot(axis, axis))
    a = math.cos(theta / 2.0)
    b, c, d = -axis * math.sin(theta / 2.0)
    aa, bb, cc, dd = a * a, b * b, c * c, d * d
    bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
    return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
                     [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                     [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])

def rotate_point_list(point_list,axis_list,theta_list):
    point_list_return = []
    for point in point_list:
    for axis,theta in zip(axis_list,theta_list):
        point = np.dot(rotation_matrix(axis,theta),point)
        point_list_return.append(point)

    return point_list_return

def angle_between_points(p1,p2,p3):
    v0 = p1 - p2
    v1 = p3 - p2

    angle = np.math.atan2(np.linalg.det([v0,v1]),np.dot(v0,v1))
    return np.degrees(angle)

def plot_3D(square,arrangement,title='',show_projection = False):
    x = [point[0] for point in square]
    y = [point[1] for point in square]
    z = [point[2] for point in square]
    x.append(x[0])
    y.append(y[0])
    z.append(z[0])

    sa = fig.add_subplot(arrangement, projection='3d')
    sa.plot(x,y,z,color="darkblue")

    if show_projection is True:
        colors = ["#FF0000", "#00FF00","#0000FF","#000000"]
        for xx,yy,zz,cc in zip(x,y,z,colors):
            sa.scatter(xx,yy,zz,color=cc)
            xl = []
            xl.append(xx)
            xl.append(xx)
            yl = []
            yl.append(yy)
            yl.append(yy)
            zl = []
            zl.append(zz)
            zl.append(0)
            sa.plot(xl,yl,zl,color = cc)

        z = [0 for point in x]
        sa.plot(x,y,z,color="orange")

    AXIS_LIM = 2

    sa.set_xlabel('X')
    sa.set_ylabel('Y')
    sa.set_zlabel('Z')
    sa.set_xlim(-AXIS_LIM,AXIS_LIM)
    sa.set_ylim(-AXIS_LIM,AXIS_LIM)
    sa.set_zlim(-AXIS_LIM,AXIS_LIM)
    sa.set_title(title)


def plot_2D(square,arrangement,title='',color = 'darkorange'):
    x = [point[0] for point in square]
    y = [point[1] for point in square]
    x.append(x[0])
    y.append(y[0])
    sa = fig.add_subplot(arrangement)
    sa.plot(x,y,color = color)
    colors = ["#FF0000", "#00FF00","#0000FF","#000000"]
    for xx,yy,cc in zip(x,y,colors):
        sa.scatter(xx,yy,color=cc)

    AXIS_LIM = 2


    sa.set_xlabel('X')
    sa.set_ylabel('Y')
    sa.set_xlim(-AXIS_LIM,AXIS_LIM)
    sa.set_ylim(-AXIS_LIM,AXIS_LIM)
    sa.set_title(title)

def remove_z(point):
    return np.delete(point,2)

def add_z(point):
    return np.append(point,0)

fig = plt.figure()


show_the_rotated_square()

Answer 1

Call $ABCD$ (named counterclockwise) the ''orange'' parallelogram (see figure below). Consider the prismatic volume $V$ (with four vertical sides) intersecting the horizontal plane along this parallelogram. Up to a translation, one can assume that the coordinates of our points are : $$A(0,0,0),B(x_B,y_B,0),D(x_D,y_D,0),B'(x_B,y_B,b),D'(x_D,y_D,d)$$ with unknowns $b$ and $d$ (positive or negative).

The constraints are

$$(1) \ \ \|\overrightarrow{AB'}\|^2=\|\overrightarrow{AD'}\|^2 \ \ \text{(same - squared - lengths) and}$$

$$(2) \ \ \overrightarrow{AB'}.\overrightarrow{AD'}=0 \ \ \text{(orthogonality).}$$

(dot product = 0 means orthogonality)

(1) becomes : $(1') \ \ (x_B^2+y_B^2)+b^2=(x_D^2+y_D^2)+d^2.$

and (2) becomes $ (2') \ \ x_Bx_D+y_By_D+bd=0.$

(classical formula for dot product).

System (1')+(2') of two equations with two unknowns is rather easy to solve. Here is how. Let us write it under the form :

$$b^2-d^2=u \ \ \text{and} \ \ bd=v$$

$u$ and $v$ being constants.

Let us plug $d=\dfrac{v}{b}$ (we can assume WLOG that $b \neq 0$) into the first equation, giving $b^2-\dfrac{v^2}{b^2}=u$ ; set $\beta=b^2$ yielding $\beta-\dfrac{v^2}{\beta}=u$, i.e., a quadratic equation that we solve for $\beta$ ; we then take (plus or minus) its square root to get $b$, out of which we deduce the other unknown $d=\dfrac{v}{b}$.

One can wonder what happens if $\beta$ is negative (square root of a negative number...). This case should never happen, first of all because parallelogram $ABCD$ comes from a ''real'' (with quotes !) projection of a ''real'' square that our equations must find out. But also for a theoretical reason that I have not the time now to write down in detail.