Euclidean geometry



In [ ]:

    
# 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

2D transformations



In [ ]:

    
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()

Rotation around the origin

$$ R(\theta) = \begin{pmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{pmatrix} $$

$$ \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



In [ ]:

    
# 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")

3D transformations



In [ ]:

    
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)

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



In [ ]:

    
# 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")

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



In [ ]:

    
# 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")

Projections

Project 2D points on a line

Line equation



In [ ]:

    
@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()

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.



In [ ]:

    
# 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)

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.



In [ ]:

    
# TODO...

Line defined by a point and an angle



In [ ]:

    
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)

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 homogeneous coordinates, 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



In [ ]:

    
# 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")

Multiple points version



In [ ]:

    
@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")