---
title: Session.  3D Objects in Python
categories: session
---

+ [Previous: [3D Modelling]()]-[Up: [Overview]()]-[Next: [3D Modelling Part II]()]

# Learning Objectives

**Key learning outcome**
Improve the understanding of the mathematical descriptions of
3D Motion, by testing implementations in Python.

**Secondary outcomes**
if you have time, it is worth browsing different libraries and frameworks
to build, visualise, and animate scenes using 3D objects.

# Reading
 
+ Ma (2004) Chapter 2.4-2.5
+ (Szeliski Chapter 2 until Section 2.1.3 inclusive)

# Briefing

**New Material** 

1. [Representations of 3D Motion]()
2. [Homogeneous Coordinates]()

## Recap: a simple object

Remember, last week we worked with this data structure in Python:

```
In [1]: print(vertices)
[[-1.   0.5  0.5]
 [ 1.   0.5  0.5]
 [ 0.  -0.5  0.5]
 [-1.   0.5  0.5]
 [ 1.   0.5  0.5]
 [ 0.   0.5 -0.5]
 [-1.   0.5  0.5]
 [ 0.  -0.5  0.5]
 [ 0.   0.5 -0.5]
 [ 1.   0.5  0.5]
 [ 0.  -0.5  0.5]
 [ 0.   0.5 -0.5]]
```

The rows are points in 3D. Note that there are only four distinct points.
If we divide the matrix into sets of three rows, each triplet defines a
triangle.  These four triangles form the faces of an irregular tetrahedron.

This is a standard way to define a 3D object.  More complex objects need more
triangles.

Note that Ma (2004) have vertices as *column* vectors, while
we here use row vectors.  This means that we need to transpose
matrices when we translate textbook formulæ to python formulæ.

Consider, for instance a rotation matrix

$$
\begin{align}
R =
\begin{bmatrix}
  \cos(\pi/6) & -\sin(\pi/6) & 0 \\
  \sin(\pi/6) & \cos(\pi/6) & 0 \\
  0 & 0 &  1
\end{bmatrix}
\end{align}
$$

```python
from numpy import *
R = array( [[ cos(pi/6), -sin(pi/6), 0], [ sin(pi/6), cos(pi/6), 0], [0,0,1]] )
```

This matrix is orthonormal, which we can check:

```python
matmul( R, transpose(R) )
```
To rotate a vector $\vec{v}$, we would calculate $\vec{u}= R\cdot \vec{v}$, where
$\vec{u}$ and $\vec{v}$ are column vectors.  If you have a shape $V$ where the columns
are points, it could be rotated as $U=R\cdot V$. In python this has to become

```python
newvertices = vertices @ R.T
```

(Note, `@` is shorthand for the `np.matmul()` function.)

Translations are straight forward; adding a row vector to a matrix, python
makes the same addition to every line.

```python
newvertices = vertices + array([1,1,0])
```

(Note that if you want to make both the translation and the rotation, the 
translation must be added to the result of the rotation and not to the original 
`vertices` prior to rotation.)


Recall, to visualise,

```python
ob1 = Poly3DCollection(vertices, linewidths=1, alpha=0.2)
ob2 = Poly3DCollection(newvertices, linewidths=1, alpha=0.2)
ob1.set_facecolor( [1, 0.5, 0.5] )
ob2.set_facecolor( [0.5, 1, 0.5] )
ob1.set_edgecolor([0,0,0])
ob2.set_edgecolor([0,0,0])
ax.add_collection3d(ob1)
ax.add_collection3d(ob2)
plt.show()
```

## Homogeneous co-ordinates

As we know, a lot of operations are simpler in homogeneous co-ordinates.

We can add a one to the vertices as follows:

```python
In [2]: vertices.shape
Out[2]: (12, 3)

In [3]: hvertices = hstack( [vertices, ones([12,1])] )

In [4]: print( hvertices)
[[-1.   0.5  0.5  1. ]
 [ 1.   0.5  0.5  1. ]
 [ 0.  -0.5  0.5  1. ]
 [-1.   0.5  0.5  1. ]
 [ 1.   0.5  0.5  1. ]
 [ 0.   0.5 -0.5  1. ]
 [-1.   0.5  0.5  1. ]
 [ 0.  -0.5  0.5  1. ]
 [ 0.   0.5 -0.5  1. ]
 [ 1.   0.5  0.5  1. ]
 [ 0.  -0.5  0.5  1. ]
 [ 0.   0.5 -0.5  1. ]]

In [5]: 
```

Suppose vi have the rotational matrix $R$ as above, and want
to translate by a vector $\vec{a} = [1,1,0]^T$.
We can build the transform matrix as

```python
a = array( [[1,1,0]] )
H1 = hstack( [R, transpose(a)] )
H = vstack( [ H1, array([0,0,0,1] ) ] )
```

Using this transform, we can transform the vertices
in homogeneous co-ordinates, and also retrieve an
inhomogeneous form.

```python
newhvertices = transpose( matmul( H, transpose(hvertices) ) )
newvertices = newhvertices[:,:3]
```

## STL files and the STL library

```sh
pip3 install python-utils
pip3 install numpy-stl
```


STL is a standard format for writing surface meshes.
The following example defines the same tetrahedron as we have
been working with, but it creates a mesh object using the
numpy-stl library.

```python
import numpy as np
from stl import mesh

corners = [ [-1,0.5,0.5], [+1,0.5,0.5], [0,-0.5,0.5], [0,0.5,-0.5] ]

vertices = np.array( corners )

face1 = [ 0, 1, 2 ]
face2 = [ 0, 1, 3 ]
face3 = [ 0, 2, 3 ]
face4 = [ 1, 2, 3 ]

faces = np.array( [ face1, face2, face3, face4 ] )

th = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
for i, f in enumerate(faces):
    for j in range(3):
        print( (i,j), vertices[f[j],:])
        th.vectors[i][j] = vertices[f[j]]

th.save('tetrahedron.stl')
```

Here, we save it to file.  We can also show it.

```python
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

ob = Poly3DCollection(th.vectors, linewidths=1, alpha=0.2)

ob.set_facecolor( [0.5, 0.5, 1] )
ob.set_edgecolor([0,0,0])


plt.ion()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.show()

ax.add_collection3d(ob)
```

This is essentially the same code as we have used before,
just note how we access the vectors from  the `th` object.
Note that the `th.vectors` array is three-dimensional:

```
In [2]: cube.vectors
Out[2]: 
array([[[-1. ,  0.5,  0.5],
        [ 1. ,  0.5,  0.5],
        [ 0. , -0.5,  0.5]],

       [[-1. ,  0.5,  0.5],
        [ 1. ,  0.5,  0.5],
        [ 0. ,  0.5, -0.5]],

       [[-1. ,  0.5,  0.5],
        [ 0. , -0.5,  0.5],
        [ 0. ,  0.5, -0.5]],

       [[ 1. ,  0.5,  0.5],
        [ 0. , -0.5,  0.5],
        [ 0. ,  0.5, -0.5]]], dtype=float32)
```

This makes it a lot more difficult to use low-level matrix
operations to rotate and translate the object.

### Reshaping

```python
In [25]: th.vectors.shape
Out[25]: (4, 3, 3)

In [26]: th.vectors.reshape(12,3)
Out[26]: 
array([[-1. ,  0.5,  0.5],
       [ 1. ,  0.5,  0.5],
       [ 0. , -0.5,  0.5],
       [-1. ,  0.5,  0.5],
       [ 1. ,  0.5,  0.5],
       [ 0. ,  0.5, -0.5],
       [-1. ,  0.5,  0.5],
       [ 0. , -0.5,  0.5],
       [ 0. ,  0.5, -0.5],
       [ 1. ,  0.5,  0.5],
       [ 0. , -0.5,  0.5],
       [ 0. ,  0.5, -0.5]], dtype=float32)

In [27]: 
```

To be certain of what actually happens, it may be easiest
to reshape the matrix before rotation and translation, and
then reshape it again (`.reshape(4,3,3)` in this case)

### Loading STL files

```python
teapot = mesh.Mesh.from_file("teapot.stl")
```

Note that large meshes are computationally demanding and one
is likely to need a visualisation engine which is designed
for meshes, rather than `matplotlib`.

### More reading

+ [Documentation for numpy-stl](https://pypi.org/project/numpy-stl/)
+ [A simple tutorial](https://spltech.co.uk/how-to-create-a-3d-model-of-a-photo-using-python-numpy-and-google-colab-part-i/)

# Exercise

## Rotations and translations

### Exercise 1.

Consider the translation vector $\vec{v}=[0,0,5]$
and the rotation matrix
$$R =
\begin{bmatrix}
-1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & -1 
\end{bmatrix}
$$
which rotates 180⁰ around the $y$-axis.

What is the difference between

1. first translating by $\vec{v}$ and then rotating by $R$, and
1. first rotating by $R$ and then translating by $\vec{v}$?

Try to reason first, and then check it in python.

### Exercise 2.

Make a 3D object `ob` and display it.

For each of the following matrices $R_i$, do the following:

1.  Define the 3D object `ob1` by multiplying `ob` by $R_i$ and
    translating by a vector $\vec{v}=[2,0,0]^T$.
2.  Display `ob1` together and `ob` and describe in your own
    words what the transformation has done.
    Is it a rigid body motion?
3.  Calculate the determinant $|R_i|$.  You can use
    `numpy.linalg.det` in python.
4.  Check if $R_i$ is orthonormal, i.e. whether
    $R_iR_i^T=I$.

When you have done a couple of example:

+ reflect on what the determinant tells you about the transformation, and
+ identify which matrices $R_i$ are rotations.


$$R_1 =
\begin{bmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 
\end{bmatrix}
$$

$$R_2 =
\begin{bmatrix}
\frac{\sqrt2}2 & -\frac{\sqrt2}2  & 0 \\
\frac{\sqrt2}2 & \frac{\sqrt2}2  & 0 \\
0 & 0 & 1 
\end{bmatrix}
$$

$$R_3 =
\begin{bmatrix}
\frac{\sqrt2}2 & \frac{\sqrt2}2  & 0 \\
\frac{\sqrt2}2 & -\frac{\sqrt2}2  & 0 \\
0 & 0 & 1 
\end{bmatrix}
$$

$$R_4 =
\begin{bmatrix}
\frac{\sqrt2}2 & 0 & -\frac{\sqrt2}2 \\
0 & 1 & 0 \\
\frac{\sqrt2}2 & 0 & \frac{\sqrt2}2  & 0 
\end{bmatrix}
$$

$$R_5 =
\begin{bmatrix}
1 & 0 & 0 \\
0 & 0 & 1 \\
0 & 1 & 0 \\
\end{bmatrix}
$$

$$R_6 =
\begin{bmatrix}
0 & 1 & 0 \\
1 & 0 & 0 \\
0 & 0 & 1 \\
\end{bmatrix}
$$

$$R_7 =
\begin{bmatrix}
2 & 0 & 0 \\
0 & 2 & 0 \\
0 & 0 & 2 \\
\end{bmatrix}
$$

$$R_8 = R_6\cdot R_7 =
\begin{bmatrix}
0 & 2 & 0 \\
2 & 0 & 0 \\
0 & 0 & 2 \\
\end{bmatrix}
$$

$$R_9 =
\begin{bmatrix}
\frac13 & 0 & 0 \\
0 & \frac13 & 0 \\
0 & 0 & \frac13 \\
\end{bmatrix}
$$

$$R_{10}=R_2\cdot R_4$$

## Homogenous Coordinates 

### Exercise 3.

Let $R_i$ be defined as in Exercise 1.
And let $\vec{b}=[1,1,0]^T$ and $\vec{c}=[0,1,1]^T$.

Define the homogeneous matrices $H_1$ and $H_2$ to represent the
following operations

1. rotate by $R_2$ and translate by $\vec{b}$
2. rotate by $R_4$ and translate by $\vec{c}$

Visualise again the object `ob` as in Exercise 1,
as well as `ob1` and `ob2` resulting from the
transforms $H_1$ and $H_2$ respectively.

What happens if you transform `ob` by the matrix
$H=H_1H_2$?  
Try to imagine the transformation first.
Then run it in python and see if you are correct or not.

Finally, try the transformation $H'=H_2H_1$.
What do you think will happen?  What actually happens?

### Exercise 4.

Define the homogeneous matrices

1.  $H_1$ which rotates by $30^\circ$
    around the $z$-axis and translates by $(2,0,0)$.
1.  $H_2$ which rotates by $60^\circ$
    around the $x$-axis and translates by $(-2,0,0)$.

Using the existing 3D object `ob`, test the two 
transformations $H_1H_2$ and $H_2H_1$.
Do they give the same result, or different results?  Why?

# Additional Libraries and Documentation

+ [Tiny 3D Engine](https://tiny-3d-engine.readthedocs.io/en/latest/)
  allows you to experiment with Scene Graphs without any fluff.
+ [PyVista](https://docs.pyvista.org/) is able to visualise large and
  complex surface meshes (STL files)
+ [plotly](https://plotly.com/python/getting-started/)
  is a popular alternative to `matplotlib`.
  It is more modern and efficient, but a little harder to get started
  with.

# Debrief

+ We use experiments in python to check our reasoning.
+ If you are stuck, please show us the experiment,
   and we discuss what sense we can make thereof.