<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Euclidean geometry
Step1:  2D transformations
Step2:  Rotation around the origin
$$
R(\theta) =
\begin{pmatrix}
    \cos \theta & -\sin \theta \
    \sin \theta & \cos \theta
\end{pmatrix}
$$
<img src="https
Step3:  3D transformations
Step4:  Rotation around the x axis
$$
R_{\mathbf {x}}(\theta) =
\begin{pmatrix}
    1& 0           & 0 \
    0& \cos \theta & -\sin \theta \
    0& \sin \theta & \cos \theta
\end{pmatrix}
$$
See
Step5:  Rotation around a given axis
Rotation around the axis defined by the unit vector $\overrightarrow{u} \begin{pmatrix} u_{x} \ u_{y} \ u_{z} \end{pmatrix}$ (i.e. with $u_{x}^{2}+u_{y}^{2}+u_{z}^{2}=1$) by an angle $\theta$.
$$
R = 
\begin{pmatrix}
    u_{x}^{2}(1-c)+c       & u_{x}u_{y}(1-c)-u_{z}s & u_{x}u_{z}(1-c)+u_{y}s \
    u_{x}u_{y}(1-c)+u_{z}s & u_{y}^{2}(1-c)+c       & u_{y}u_{z}(1-c)-u_{x}s \
    u_{x}u_{z}(1-c)-u_{y}s & u_{y}u_{z}(1-c)+u_{x}s & u_{z}^{2}(1-c)+c
\end{pmatrix}
$$
where $c = \cos(\theta)$ and $s = \sin(\theta)$
See https
Step6:  Projections
Project 2D points on a line
Line equation
Step7:  Distance from a point to a line
Line defined by an equation
In the case of a line in the plane given by the equation $ax + by + c = 0$, where $a$, $b$ and $c$ are real constants with $a$ and $b$ not both zero, the distance from the line to a point $(x_0,y_0)$ is
$$\operatorname{distance}(ax+by+c=0, (x_0, y_0)) = \frac{|ax_0+by_0+c|}{\sqrt{a^2+b^2}}.$$
The point on this line which is closest to $(x_0,y_0)$ has coordinates
Step8:  Line defined by two points
If the line passes through two points $P_1 = (x_1, y_1)$ and $P_2 = (x_2, y_2)$ then the distance of $(x_0, y_0)$ from the line is
Step9:  Line defined by a point and an angle
Step10:  Project 3D points on a plane without perspective
TODO...
Project 3D points on a plane with perspective
The following variables are defined to describe this transformation
Step11:  Multiple points version
<ASSISTANT_TASK:>
Python Code:
# Import directives
#%pylab notebook
%pylab inline
pylab.rcParams['figure.figsize'] = (6, 6)
#import warnings
#warnings.filterwarnings('ignore')
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from ipywidgets import interact
Explanation: Euclidean geometry
End of explanation
def plot2d(x, fmt="ok"):
    plt.axis('equal')
    plt.axis([-5, 5, -5, 5])
    plt.xticks(np.arange(-5, 5, 1))
    plt.yticks(np.arange(-5, 5, 1))
    plt.axhline(y=0, color='k')
    plt.axvline(x=0, color='k')
    plt.plot(x[:,0], x[:,1], fmt)
    plt.grid()
Explanation: 2D transformations
End of explanation
# Define initial points
A = np.array([[0., 0.],
              [1., 0.],
              [1., 1.],
              [0., 0.]])
# Define the rotation angle
theta = np.radians(30)
# Define the rotation matrix
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])
# Rotate points
Aprime = np.dot(R, A.T).T
# Print and plot
print(A)
print(Aprime)
plot2d(A, fmt="-ok")
plot2d(Aprime, fmt="-or")
Explanation: Rotation around the origin
$$
R(\theta) =
\begin{pmatrix}
    \cos \theta & -\sin \theta \
    \sin \theta & \cos \theta
\end{pmatrix}
$$
<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/6/64/Mog_rotacion_vector.jpg/263px-Mog_rotacion_vector.jpg"></img>
$$
\begin{pmatrix}
    A'{x} \
    A'{y}
\end{pmatrix}
=
\begin{pmatrix}
    \cos \theta & -\sin \theta \
    \sin \theta &\cos \theta \
\end{pmatrix}
\begin{pmatrix}
    A_{x} \
    A_{y}
\end{pmatrix}
$$
See: https://fr.wikipedia.org/wiki/Matrice_de_rotation
End of explanation
def plot3d(x, axis=None, fmt="ok"):
    if axis is None:
        fig = plt.figure()
        axis = axes3d.Axes3D(fig)
        
    axis.scatter(x[:,0], x[:,1], x[:,2], fmt)
    axis.plot(x[:,0], x[:,1], x[:,2], fmt)
Explanation: 3D transformations
End of explanation
# Define initial points
A = np.array([[0., 0., 0.],
              [1., 0., 0.],
              [0., 0., 0.],
              [0., 1., 0.],
              [0., 0., 0.],
              [0., 0., 1.]])
# Define the rotation angle
theta = np.radians(90)
# Define the rotation matrices
Rx = np.array([[1., 0.,            0.],
               [0., np.cos(theta), -np.sin(theta)],
               [0., np.sin(theta), np.cos(theta)]])
Ry = np.array([[np.cos(theta),  0., np.sin(theta)],
               [0.,             1., 0.           ],
               [-np.sin(theta), 0., np.cos(theta)]])
Rz = np.array([[np.cos(theta), -np.sin(theta), 0.],
               [np.sin(theta),  np.cos(theta), 0.],
               [0.,             0.,            1.]])
# Rotate points
Ax = np.dot(Rx, A.T).T
Ay = np.dot(Ry, A.T).T
Az = np.dot(Rz, A.T).T
# Plot
fig = plt.figure()
ax = axes3d.Axes3D(fig)
plot3d(A, axis=ax, fmt="-ok")
plot3d(Ax, axis=ax, fmt=":or")
plot3d(Ay, axis=ax, fmt=":og")
plot3d(Az, axis=ax, fmt=":ob")
ax.text(1, 0, 0, "x", color="r")
ax.text(0, 1, 0, "y", color="g")
ax.text(0, 0, 1, "z", color="b")
Explanation: Rotation around the x axis
$$
R_{\mathbf {x}}(\theta) =
\begin{pmatrix}
    1& 0           & 0 \
    0& \cos \theta & -\sin \theta \
    0& \sin \theta & \cos \theta
\end{pmatrix}
$$
See: https://fr.wikipedia.org/wiki/Matrice_de_rotation
Rotation around the y axis
$$
R_{\mathbf {y}}(\theta )
=
\begin{pmatrix}
    \cos \theta  & 0 & \sin \theta \
    0            & 1 & 0 \
    -\sin \theta & 0 & \cos \theta
\end{pmatrix}
$$
See: https://fr.wikipedia.org/wiki/Matrice_de_rotation
Rotation around the z axis
$$
R_{\mathbf {z}}(\theta )
=
\begin{pmatrix}
    \cos \theta & -\sin \theta & 0 \
    \sin \theta & \cos \theta  & 0 \
    0           & 0            & 1
\end{pmatrix}
$$
See: https://fr.wikipedia.org/wiki/Matrice_de_rotation
Example
End of explanation
# Define initial points
A = np.array([[0., 0., 0.],
              [1., 0., 0.],
              [0., 0., 0.],
              [0., 1., 0.],
              [0., 0., 0.],
              [0., 0., 1.]])
# Define the rotation angle
theta = np.radians(10)
u = np.array([1., 1., 0.])
ux, uy, uz = u[0], u[1], u[2]
c = np.cos(theta)
s = np.sin(theta)
# Define the rotation matrices
R = np.array([[ux**2 * (1-c) + c,     ux*uy * (1-c) - uz*s,  ux*uz * (1-c) + uy*s],
              [ux*uy * (1-c) + uz*s,  ux**2 * (1-c) + c,     uy*uz * (1-c) - ux*s],
              [ux*uz * (1-c) - uy*s,  uy*uz * (1-c) + ux*s,  uz**2 * (1-c) + c]])
# Rotate points
Ar = np.dot(R, A.T).T
# Plot
fig = plt.figure()
ax = axes3d.Axes3D(fig)
plot3d(A,  axis=ax, fmt="-ok")
plot3d(np.array([np.zeros(3), u]),  axis=ax, fmt="--ok")
plot3d(Ar, axis=ax, fmt=":or")
ax.text(1, 0, 0, "x", color="k")
ax.text(0, 1, 0, "y", color="k")
ax.text(0, 0, 1, "z", color="k")
Explanation: Rotation around a given axis
Rotation around the axis defined by the unit vector $\overrightarrow{u} \begin{pmatrix} u_{x} \ u_{y} \ u_{z} \end{pmatrix}$ (i.e. with $u_{x}^{2}+u_{y}^{2}+u_{z}^{2}=1$) by an angle $\theta$.
$$
R = 
\begin{pmatrix}
    u_{x}^{2}(1-c)+c       & u_{x}u_{y}(1-c)-u_{z}s & u_{x}u_{z}(1-c)+u_{y}s \
    u_{x}u_{y}(1-c)+u_{z}s & u_{y}^{2}(1-c)+c       & u_{y}u_{z}(1-c)-u_{x}s \
    u_{x}u_{z}(1-c)-u_{y}s & u_{y}u_{z}(1-c)+u_{x}s & u_{z}^{2}(1-c)+c
\end{pmatrix}
$$
where $c = \cos(\theta)$ and $s = \sin(\theta)$
See https://fr.wikipedia.org/wiki/Matrice_de_rotation
End of explanation
@interact(a=(-5., 5., 0.1), b=(-5., 5., 0.1), c=(-5., 5., 0.1))
def plot(a, b, c):
    plt.axis('equal')
    plt.axis([-5, 5, -5, 5])
    plt.xticks(np.arange(-5,5,1))
    plt.yticks(np.arange(-5,5,1))
    plt.axhline(y=0, color='k')
    plt.axvline(x=0, color='k')
    x = np.array([-10., 10.])
    f = lambda x: a/(-b) * x + c/(-b)
    
    try:
        plt.plot(x, f(x))
    except ZeroDivisionError:
        print("b should not be equal to 0")
    plt.grid()
Explanation: Projections
Project 2D points on a line
Line equation
End of explanation
# Setup the plot
def plot(a, b, c, p, p2):
    plt.axis('equal')
    plt.axis([-5, 5, -5, 5])
    plt.xticks(np.arange(-5,5,1))
    plt.yticks(np.arange(-5,5,1))
    plt.axhline(y=0, color='k')
    plt.axvline(x=0, color='k')
    x = np.array([-10., 10.])
    f = lambda x: a/(-b) * x + c/(-b)
    
    plt.plot(x, f(x))
    plt.scatter(*p)
    # Plot the projection point
    plt.scatter(*p2)
    plt.plot((p2[0], p[0]), (p2[1], p[1]))
    #plt.arrow(*p2, *p) # TODO: doesn't work...
    plt.grid()
# Define the distance and projection functions
def distance(a, b, c, p):
    d1 = (a*p[0] + b*p[1] + c)
    d2 = math.sqrt(math.pow(a, 2.) + math.pow(b, 2.))
    d = abs(d1)/d2
    return d
def projection(a, b, c, p):
    p2 = ((b*(b*p[0] - a*p[1]) - a*c)/(math.pow(a,2.)+math.pow(b,2.)),
          (a*(-b*p[0] + a*p[1]) - b*c)/(math.pow(a,2.)+math.pow(b,2.)))
    return p2
# Define the line and the point
a = 2.
b = 1.
c = -2.
p = (-4., 2.)
# Compute the distance and the projection point on the line
d = distance(a, b, c, p)
p2 = projection(a, b, c, p)
print("Distance:", d)
print("Projection point:", p2)
# Plot the line and the point
plot(a, b, c, p, p2)
Explanation: Distance from a point to a line
Line defined by an equation
In the case of a line in the plane given by the equation $ax + by + c = 0$, where $a$, $b$ and $c$ are real constants with $a$ and $b$ not both zero, the distance from the line to a point $(x_0,y_0)$ is
$$\operatorname{distance}(ax+by+c=0, (x_0, y_0)) = \frac{|ax_0+by_0+c|}{\sqrt{a^2+b^2}}.$$
The point on this line which is closest to $(x_0,y_0)$ has coordinates:
$$x = \frac{b(bx_0 - ay_0)-ac}{a^2 + b^2}$$
and
$$y = \frac{a(-bx_0 + ay_0) - bc}{a^2+b^2}.$$
For mor information, see https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line.
End of explanation
# TODO...
Explanation: Line defined by two points
If the line passes through two points $P_1 = (x_1, y_1)$ and $P_2 = (x_2, y_2)$ then the distance of $(x_0, y_0)$ from the line is:
$$\operatorname{distance}(P_1, P_2, (x_0, y_0)) = \frac{|(y_2-y_1)x_0-(x_2-x_1)y_0+x_2 y_1-y_2 x_1|}{\sqrt{(y_2-y_1)^2+(x_2-x_1)^2}}.$$
The denominator of this expression is the distance between $P_1$ and $P_2$.
The numerator is twice the area of the triangle with its vertices at the three points, $(x_0, y_0)$, $P_1$ and $P_2$.
The expression is equivalent to $h=\frac{2A}{b}$, which can be obtained by rearranging the standard formula for the area of a triangle: $A = \frac{1}{2} bh$, where $b$ is the length of a side, and $h$ is the perpendicular height from the opposite vertex.
End of explanation
def angle_point_to_equation(angle_degree, p):
    angle_radian = math.radians(angle_degree)
    
    a = math.tan(angle_radian)
    b = -1
    c = -math.tan(angle_radian) * p[0] + p[1]
    
    return a, b, c
angle_degree = 30
p0 = (3, 2)
a, b, c = angle_point_to_equation(angle_degree, p0)
p = (-4., 2.)
# Compute the distance and the projection point on the line
d = distance(a, b, c, p)
p2 = projection(a, b, c, p)
print("Distance:", d)
print("Projection point:", p2)
# Plot the line and the point
plot(a, b, c, p, p2)
plt.scatter(*p0)
Explanation: Line defined by a point and an angle
End of explanation
# Define initial points to project
a = np.array([0., 1., 2.])
# Define camera's position
c = np.array([0., 0., 0.])
# Define viewer's position
e = np.array([0., 0., -1.])
# Define the orientation of the camera
theta = np.array([np.radians(0),
                  np.radians(0),
                  np.radians(0)])
theta_x, theta_y, theta_z = theta[0], theta[1], theta[2]
# Define the rotation matrices
Rx = np.array([[1., 0.,               0.],
               [0., np.cos(theta_x),  np.sin(theta_x)],
               [0., -np.sin(theta_x), np.cos(theta_x)]])
Ry = np.array([[np.cos(theta_y),  0., -np.sin(theta_y)],
               [0.,               1., 0.           ],
               [np.sin(theta_y),  0., np.cos(theta_y)]])
Rz = np.array([[np.cos(theta_z),  np.sin(theta_z), 0.],
               [-np.sin(theta_z), np.cos(theta_z), 0.],
               [0.,               0.,              1.]])
d = np.dot(Rx, Ry)
d = np.dot(d, Rz)
d = np.dot(d, a-c)
## TODO: which version is correct ? The one above or the one below ?
#d = a - c
#d = np.dot(Rz, d)
#d = np.dot(Ry, d)
#d = np.dot(Rx, d)
print("d:", d)
b = np.array([e[2]/d[2] * d[0] - e[0],
              e[2]/d[2] * d[1] - e[1]])
print("b:", b)
# Alternative to compute b
Rf = np.array([[1., 0., -e[0]/e[2], 0.],
               [0., 1., -e[1]/e[2], 0.],
               [0., 0., 1.,         0.],
               [0., 0., 1./e[2],    0.]])
f = np.dot(Rf, np.concatenate([d, np.ones(1)]))
b = np.array([f[0]/f[3],
              f[1]/f[3]])
print("b:", b)
plot2d(np.array([b, b]), "ok")
Explanation: Project 3D points on a plane without perspective
TODO...
Project 3D points on a plane with perspective
The following variables are defined to describe this transformation:
$\mathbf{a}_{x,y,z}$: the 3D position of a point A that is to be projected.
$\mathbf{c}_{x,y,z}$: the 3D position of a point C representing the camera.
$\mathbf{\theta}_{x,y,z}$: the orientation of the camera (represented by Tait–Bryan angles).
$\mathbf{e}_{x,y,z}$: the viewer's position relative to the display surface which goes through point C representing the camera (i.e. $\mathbf{e}$ doesn't change with the position of the camera C and typically $\mathbf{e} = (0, 0, -1)^T$ for a view angle of $2*\cos^{-1}(1/\sqrt{2})$ radians when the projection is made on the x/y plane).
Which results in:
$\mathbf {b} _{x,y}$: the 2D projection of a $\mathbf{a}$
$$
\begin{bmatrix}
    \mathbf {d} {x} \
    \mathbf {d} {y} \
    \mathbf {d} {z} \
\end{bmatrix}
=
\begin{bmatrix}
    1 & 0                          & 0 \
    0 & \cos(\mathbf{\theta}{x})  & \sin(\mathbf{\theta}{x}) \
    0 & -\sin(\mathbf{\theta}{x}) & \cos(\mathbf{\theta}{x}) \
\end{bmatrix}
~
\begin{bmatrix}
    \cos(\mathbf{\theta}{y}) & 0 & -\sin(\mathbf{\theta}{y}) \
    0                         & 1 & 0 \
    \sin(\mathbf{\theta}{y}) & 0 & \cos(\mathbf{\theta}{y}) \
\end{bmatrix}
~
\begin{bmatrix}
    \cos(\mathbf{\theta}{z})  & \sin(\mathbf{\theta}{z}) & 0 \
    -\sin(\mathbf{\theta}{z}) & \cos(\mathbf{\theta}{z}) & 0 \
    0                          & 0                         & 1 \
\end{bmatrix}
~
\left(
\begin{bmatrix}
    \mathbf{a}{x} \
    \mathbf{a}{y} \
    \mathbf{a}{z} \
\end{bmatrix}
-
\begin{bmatrix}
    \mathbf{c}{x} \
    \mathbf{c}{y} \
    \mathbf{c}_{z} \
\end{bmatrix}
\right)
$$
The transformed point $\boldsymbol d$ of $\boldsymbol d$ can then be projected onto the 2D plane using the formula (here, x/y is used as the projection plane; literature also may use x/z):
$$
\begin{array}{lcl}
    \mathbf{b}_x & = & \frac{\mathbf{e}_z}{\mathbf{d}_z} \mathbf{d}_x - \mathbf{e}_x \
    \mathbf{b}_y & = & \frac{\mathbf{e}_z}{\mathbf{d}_z} \mathbf{d}_y - \mathbf{e}_y \
\end{array}
$$
Or, in matrix form using <a href="https://en.wikipedia.org/wiki/Homogeneous_coordinates">homogeneous coordinates</a>, the system
$$
\begin{bmatrix}
    \mathbf{f}_x \
    \mathbf{f}_y \
    \mathbf{f}_z \
    \mathbf{f}_w
\end{bmatrix}
=
\begin{bmatrix}
    1 & 0 & -\frac{\mathbf{e}_x}{\mathbf{e}_z} & 0 \
    0 & 1 & -\frac{\mathbf{e}_y}{\mathbf{e}_z} & 0 \
    0 & 0 & 1                                  & 0 \
    0 & 0 & \frac{1}{\mathbf{e}_z}             & 0
\end{bmatrix}
\begin{bmatrix}
    \mathbf{d}_x \
    \mathbf{d}_y \
    \mathbf{d}_z \
    1
\end{bmatrix}
$$
in conjunction with an argument using similar triangles, leads to division by the homogeneous coordinate, giving
$$
\begin{array}{lcl}
    \mathbf{b}_x & = & \mathbf{f}_x / \mathbf{f}_w \
    \mathbf{b}_y & = & \mathbf{f}_y / \mathbf{f}_w \
\end{array}
$$
The distance of the viewer from the display surface, $\mathbf{e}{z}$, directly relates to the field of view, where $\alpha = 2 \cdot \tan^{-1}(\frac{1}{\mathbf{e}{z}})$ is the viewed angle.
Note: here we assume that you map the points $(-1,-1)$ and $(1,1)$ to the corners of your viewing surface.
See: https://en.wikipedia.org/wiki/3D_projection#Perspective_projection
End of explanation
@interact(theta_x=(-90., 90., 1.), theta_y=(-90., 90., 1.), theta_z=(-90., 90., 1.))
def projection(theta_x, theta_y, theta_z):
    
    # Define initial points to project
    A = np.array([[-1., 0., 1.],
                  [ 1., 0., 1.],
                  [-1., 0., 2.],
                  [ 1., 0., 2.],
                  [-1., 0., 5.],
                  [ 1., 0., 5.],
                  [-1., 0., 15.],
                  [ 1., 0., 15.]])
    # Define camera's position
    c = np.array([0., -2., 0.])
    C = np.tile(c, (A.shape[0], 1))
    # Define viewer's position
    e = np.array([0., 0., -1.])
    # Define the orientation of the camera
    theta = np.radians(np.array([theta_x, theta_y, theta_z]))
    theta_x, theta_y, theta_z = theta[0], theta[1], theta[2]
    # Define the rotation matrices
    Rx = np.array([[1., 0.,               0.],
                   [0., np.cos(theta_x),  np.sin(theta_x)],
                   [0., -np.sin(theta_x), np.cos(theta_x)]])
    Ry = np.array([[np.cos(theta_y),  0., -np.sin(theta_y)],
                   [0.,               1., 0.           ],
                   [np.sin(theta_y),  0., np.cos(theta_y)]])
    Rz = np.array([[np.cos(theta_z),  np.sin(theta_z), 0.],
                   [-np.sin(theta_z), np.cos(theta_z), 0.],
                   [0.,               0.,              1.]])
    d = np.dot(Rx, Ry)
    d = np.dot(d, Rz)
    d = np.dot(d, (A-C).T)
    ## TODO: which version is correct ? The one above or the one below ?
    #d = a - c
    #d = np.dot(Rz, d)
    #d = np.dot(Ry, d)
    #d = np.dot(Rx, d)
    print("d:", d)
    b = np.array([e[2]/d[2] * d[0] - e[0],
                  e[2]/d[2] * d[1] - e[1]])
    print("b:", b)
    # Alternative to compute b
    Rf = np.array([[1., 0., -e[0]/e[2], 0.],
                   [0., 1., -e[1]/e[2], 0.],
                   [0., 0., 1.,         0.],
                   [0., 0., 1./e[2],    0.]])
    # Add a line of ones
    d = np.vstack([d, np.ones(d.shape[1])])
    f = np.dot(Rf, d)
    b = np.array([f[0]/f[3],
                  f[1]/f[3]])
    print("b:", b)
    plot2d(b.T, "ok")
    plot2d(b.T, "-k")
Explanation: Multiple points version
End of explanation
<END_TASK>