Frolian's blog standing on the shoulders of PyGiants

• About
• Archive

Building an Eiffel Tower made of Springs

Sat, 02 Nov 2019

Physics Numerical Physics Animations

As a bonus to these last days blogging about spring models, I've come to think of a fun little idea. Since it is pretty easy to define spring models, why not build one out of a picture, like that of the Eiffel tower? We could then have a bumpy Eiffel tower made from springs! Let's see if we can work this idea out.

Getting the Eiffel tower image and extracting a mask

To start, we need to download a picture of the Eiffel tower. We can use Wikipedia to get one!

import numpy as np
from PIL import Image
import requests

url = 'https://upload.wikimedia.org/wikipedia/commons/8/85/Tour_Eiffel_Wikimedia_Commons_%28cropped%29.jpg'
img = np.array(Image.open(requests.get(url, stream=True).raw))

Now let's make the image gray and do some image processing on it to extract just the shape of the Eiffel tower.

import holoviews as hv
hv.extension('matplotlib', 'bokeh', logo=False)

from skimage.transform import rescale
from skimage.color import rgb2gray

gray = rgb2gray(img)
gray = rescale(gray, scale=0.1)
hv.Raster(gray).opts(cmap='gray')

from skimage import exposure
from skimage.morphology import dilation

equalized = exposure.equalize_adapthist(gray, kernel_size=(5, 5))
dilated = dilation(equalized, selem=np.ones((3, 2)))

hv.Raster(dilated).opts(colorbar=True)

As we can see on this image, the Eiffel tower sticks out as mostly bright pixels. We can use that to extract just these points!

Extracting a triangular mesh

Let's select all the points in the image above a certain treshold.

points = np.nonzero(dilated > 0.7)
x_coords, y_coords = points[1], -points[0]

We will also just keep the pixels of interest, i.e. just the structure.

selection = (x_coords > 70) & (x_coords < 220) & (y_coords > -400)
x_coords, y_coords = x_coords[selection], y_coords[selection]

How mayn points do we end up with?

x_coords.shape

(11617,)

Let's now make and display a triangulation of the points we've ended up with.

from scipy.spatial import Delaunay

pts = np.c_[x_coords, y_coords].astype(float)
tris = Delaunay(pts) 

trimesh = hv.TriMesh((tris.simplices, pts))
trimesh.opts(hv.opts.TriMesh(node_size=0.1, edge_alpha=0.1))

Since some edges are pretty long, let's filter them out: we only want to keep the short edges.

def longest_edge(x0, x1, x2):
    """Returns length of longest edge in triangle (x0, x1, x2)."""
    x01 = x0 - x1
    x02 = x0 - x2
    return np.max([np.linalg.norm(x01), np.linalg.norm(x02)])

longests = []
for tri_index in range(tris.simplices.shape[0]):
    xi = tris.points[tris.simplices[tri_index]]
    longests.append(longest_edge(xi[0], xi[1], xi[2]))

Now, let's mask the triangulation.

mask = np.logical_not(np.where(np.array(longests) < 15, 0, 1))
simplices = tris.simplices[mask]
nodes = hv.Points(pts)
hv.TriMesh((simplices, nodes)).opts(hv.opts.TriMesh(node_size=0.1, edge_alpha=0.1))

This starts looking nice.

Making it move

Finally, we can make the tower move. To do that, we're just going to copy some of the routines I used in the previous notebook.

import numba 

@numba.jit()
def compute_spring_forces(x, spring_anchors, spring_lengths, spring_stiffnesses):
    """Computes the spring forces."""
    forces = np.zeros_like(x)
    for i in range(spring_anchors.shape[0]):
        a, b = spring_anchors[i]
        x_a, x_b = x[a], x[b]
        dist = x_a - x_b
        length = np.linalg.norm(dist) + 1e-4
        F = (length - spring_lengths[i]) * spring_stiffnesses[i] * dist / length
        forces[a] += -F
        forces[b] += F
    return forces

@numba.jit()
def time_step_integration(x, v, forces, masses, dt, damping):
    """Integrates the forces using semi-implicit Euler time integration with damping."""
    s = np.exp(-dt * damping) 
    v_new = s * v + dt * forces / masses.reshape(-1, 1)
    x_new = x + dt * v_new
    return x_new, v_new

@numba.jit()
def forward(x, spring_anchors, spring_lengths, spring_stiffnesses, v0=None, steps=100, dt=0.001, damping=20):
    new_x = x.copy()
    if v0 is None:
        new_v = np.zeros_like(x)
    else:
        new_v = v0.copy()
    snapshots = []
    snapshots.append(new_x.copy())
    for t in range(steps):
        forces = compute_spring_forces(new_x, spring_anchors, spring_lengths, spring_stiffnesses)
        new_x, new_v = time_step_integration(new_x, new_v, forces, masses, dt=dt, damping=damping)
        snapshots.append(new_x.copy())
    return snapshots

To use the above routines, we have to:

• define the point coordinates
• use the triangle information to connect points to each others with springs (that have the rest length determined by initial point positions)
• build a mass vector
• build a stiffness vector

x = pts.copy()
n_points = x.shape[0]
spring_anchors = []
spring_lengths = []
spring_stiffnesses = []
masses = []
for triangle in simplices:
    x0, x1, x2 = triangle
    p0, p1, p2 = pts[x0], pts[x1], pts[x2]
    l0, l1, l2 = np.linalg.norm(p1 - p0), np.linalg.norm(p2 - p1), np.linalg.norm(p0 - p2)
    spring_anchors.extend([[x0, x1], [x1, x2], [x2, x0]])
    spring_lengths.extend([l0, l1, l2])
    spring_stiffnesses.extend([1/l0, 1/l1, 1/l2])

spring_anchors = np.array(spring_anchors)
spring_lengths = np.array(spring_lengths)
spring_stiffnesses = np.array(spring_stiffnesses)
masses = np.ones((n_points,))

The last thing we need to define is a vector holding the initial velocities.

v0 = np.zeros_like(x)
v0[0:100][:, 1] = np.ones(v0[0:100][:, 0].shape)

Done, we can now run our simulation!

snapshots = forward(x, spring_anchors, spring_lengths, spring_stiffnesses, v0, dt=0.5, steps=1000, damping=0.01)

Let's look at some displacement plots from the simulation.

def make_time_plots(snapshots, node_number):
    return hv.Curve([snapshot[node_number, 0] for snapshot in snapshots], 
             label='node{}_x'.format(node_number), kdims='time', vdims='coordinate') * \
           hv.Curve([snapshot[node_number, 1] for snapshot in snapshots], 
             label='node{}_y'.format(node_number), kdims='time', vdims='coordinate')

make_time_plots(snapshots, 0) * make_time_plots(snapshots, 1110)

Let's look at the animation over time.

def plot_points(x, spring_anchors):
    """Plots all x as 2D points, connected by springs indicated by spring_anchors."""
    scatter = hv.Scatter((x[:, 0], x[:, 1])).opts(hv.opts.Scatter(padding=0.1, s=0.2))
    return scatter 

def plot_points2(x, spring_anchors):
    return hv.TriMesh((simplices, x)).opts(hv.opts.TriMesh(node_size=0.1, edge_alpha=0.1))

def make_animation(snapshots, n_frames):
    hmap = hv.HoloMap({i: plot_points2(snapshots[i], spring_anchors) for i in range(0, len(snapshots), n_frames)}, 
                      kdims=['snapshot'])
    return hv.output(hmap, holomap='scrubber')

make_animation(snapshots, n_frames=10)

– ▐◀◀ ▐◀ ◀ ▐ ▌ ▶ ▶▌ ▶▶▌ +

Once Loop Reflect