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.
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!
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.
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:
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)
As a last step, we can amplify the displacement and animate that.
snapshots_amplified = [10 * (snapshot - snapshots[0]) + snapshots[0] for snapshot in snapshots]
make_time_plots(snapshots_amplified, 0) * make_time_plots(snapshots_amplified, 1110)
make_animation(snapshots_amplified, n_frames=10)
(apologies for the long time it takes to generate these animations: the TriMesh function is pretty slow...)
Nice! A wiggly Eiffel tower, just as I hoped to see it. Objective accomplished.
This post was written using the Jupyter Notebook. Its content is BSD-licensed. You can see a static view or download this notebook with the help of nbviewer at 20191031_EiffelTowerTriangulation.ipynb.