import numpy as np
import numpy.linalg as la
# for in-line plots
%matplotlib inline
# for plots in a window
# %matplotlib qt
import matplotlib.pyplot as pt
from mpl_toolkits.mplot3d import Axes3D
Make two random 3D vectors:
np.random.seed(13)
x = np.random.randn(3)
y = np.random.randn(3)
Make them orthonormal:
y = y - y.dot(x)/x.dot(x)*x
x = x/la.norm(x)
y = y/la.norm(y)
Check:
print(y.dot(x))
print(la.norm(x))
print(la.norm(y))
Plot the two vectors:
fig = pt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlim3d([-3, 3])
ax.set_ylim3d([-3, 3])
ax.set_zlim3d([-3, 3])
xy = np.array([x, y]).T
ax.quiver(
0, 0, 0,
xy[0], xy[1], xy[2],)
Make an array with the cornerpoints of a cube:
points = np.array([
[-1,-1,-1],
[-1,-1,1],
[-1,1,-1],
[-1,1,1],
[1,-1,-1],
[1,-1,1],
[1,1,-1],
[1,1,1],
])
Plot them:
fig = pt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:,0], points[:, 1], points[:, 2])
Construct the projection matrix:
Q = np.array([
x,y,np.zeros(3)
]).T
print(Q)
P = Q.dot(Q.T)
print(P)
Check that $P^2=P$:
la.norm(P.dot(P)-P)
Project the points, assign to proj_points
:
proj_points = np.einsum("ij,nj->ni", P, points)
fig = pt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:,0], points[:, 1], points[:, 2])
ax.scatter(proj_points[:,0], proj_points[:, 1], proj_points[:, 2], color="red")
xy = np.array([x, y]).T
ax.quiver(
0, 0, 0,
xy[0], xy[1], xy[2],)