\( \newcommand{\bm}{\mathbf} \) Home > Misc

Perspective drawing and Moiré patterns

For as long as I can remember, I've enjoyed seeing Moiré patterns in pedestal fans: the fan is covered front and back by a cage with radial spokes, and as the fan rotates or you move, the Moiré pattern moves also. I thought it'd be interesting to draw a fan in ggplot2 and enjoy the Moiré motion.

There are two interesting steps in creating the gif above. The first is making a perspective drawing in the first place, and the second is making the random-yet-smooth motion of the observer.

Perspective

What follows is my understanding – essentially derived from scratch – of how to take points in 3-space and plot them as seen by an observer. My brief perusal of the Wikipedia article revealed quite a few types of perspective drawing, and I haven't bothered to work out how this method fits in with them.

The basic idea is to imagine a small sphere with the observer at its centre. The observer is looking at some point; draw a line from that point to the observer; where that line intersects the surface of the sphere defines zero-latitude and zero-longitude. Choosing an orientation for the observer then fixes the equatorial plane and prime meridian. The position of a object in the observer's field of view is the longitude and latitude of where a line from that object to the osberver intersects the surface of the sphere.

If we imagine our observer at the origin, looking along the positive \(x\)-axis, and oriented so that "up" is the positive \(z\)-axis, then a latitude \(\phi\) and longitude \(\theta\) correspond to a point on the surface of the unit sphere as follows:

\begin{align*} x &= \cos(\phi) \cos(\theta), \\ y &= \cos(\phi) \sin(\theta), \\ z &= \sin(\phi). \end{align*}

We can therefore find the latitude and longitude of a point \(\bm{v}\) on the surface of the unit sphere by taking dot products with unit vectors pointing in the axis directions:

\begin{align} \nonumber \phi &= \arcsin(\bm{v \cdot k}), \\ \label{angle_dotproducts} \theta &= \arcsin\left(\frac{\bm{v \cdot j}}{\cos(\phi)}\right). \end{align}

(In the equation for \(\theta\), it is better to take the dot product with \(\bm{j}\) than with \(\bm{i}\), because in the latter case you then take an arccos, which doesn't distinguish between positive and negative values of \(\theta\).)

The fact that the observer will not in general be at the origin is not important, since we can just subtract the observer's position from the relevant position vector being plotted. The goal now is to rotate the axis vectors \(\bm{i}, \bm{j}, \bm{k}\) so that they match the orientation of the observer.

The observer in some direction; let \(\bm{r}\) be the vector from the observer to some point along this direction. In the earlier equations for latitude and longitude, the observer was looking along the positive \(x\)-axis. We'll rotate the axes so that \(\bm{i}\) points in the \(\bm{r}\)-direction in two steps.

First, we spin \(\bm{i}\) and \(\bm{j}\) axes around the \(\bm{k}\)-axis so that \(\bm{i}\) is parallel to \((r_x, r_y, 0)\). This is a rotation by an angle \(\alpha\) defined by

\[ \alpha = \arctan(r_y/r_x). \]

The second step is to tilt the \(xy\)-plane, so that \(\bm{i}\) is parallel to \((r_x, r_y, r_z)\). This tilt is of an angle \(\beta\) defined by

\[ \beta = \arcsin\left(r_z/|\bm{r}|\right). \]

The actual transformations of the axis vectors are

\begin{align*} \bm{i'} &= \cos(\alpha)\bm{i} + \sin(\alpha)\bm{j}, \\ \bm{j'} &= -\sin(\alpha)\bm{i} + \cos(\alpha)\bm{j}, \\ \bm{k'} &= \bm{k}, \end{align*}

and

\begin{align*} \bm{i''} &= \cos(\beta)\bm{i'} + \sin(\beta)\bm{k'}, \\ \bm{j''} &= \bm{j'}, \\ \bm{k''} &= -\sin(\beta)\bm{i'} + \cos(\beta)\bm{k'}. \end{align*}

It is the vectors \(\bm{k''}\) and \(\bm{j''}\) that are used in the dot products of equations (\ref{angle_dotproducts}) to decide where an object appears to the observer.

Smooth random motion

For the gif at the top of this page, the basic idea is to set random positions and velocities at regular intervals, and then smoothly interpolate the motion between the ends of these intervals. (And then to get the perfect gif loop, the last position and velocity is set equal to the first.) I used R's ggplot2 package to draw all the lines – not a choice I'd recommend for this exercise unless you want to practise ggplotting – and it won't plot a geom_segment if one of the edges is outside the plotting window. If the observer gets too close to the fan, then the edges of the fan may fall outside the plotting window, and some of the spokes won't be plotted at all. So (in my case!) the smooth interpolation of the motion shouldn't go near the fan itself.

(I haven't defined the fan, but it's not hard – just a collection of vertices that get joined together in a couple of different ways.)

The natural solution to the disappearing-spokes problem is to use polar co-ordinates, with the observer's random radial distances having a suitably large lower bound. We can use a smoothstep function to ensure that the time-derivative of the radial distance is continuous.

More interesting are the angular velocity and the z-component of the velocity; I'll treat these two as independent cylindrical co-ordinates, and the theory applies to each component. We have \(x_0 = x(0),\, x_1 = x(1),\, \dot{x}(0) = v_0,\, \dot{x}(1) = v_1\) all given (and set randomly). To make the motion a little smoother, I'll also require \(\ddot{x}(0) = \ddot{x}(1) = 0\). The idea is to write the acceleration as a power series in \(t\), truncate appropriately, and solve for the coefficients.

The acceleration being zero at \(t=0\) means that the power series has no constant term,

\begin{equation*} a(t) = \sum_{n=1}^{\infty} a_n t^n, \end{equation*}

and the acceleration being zero at \(t=1\) means that

\begin{equation} \label{powerseries1} \sum_{n=1}^{\infty} a_n = 0. \end{equation}

The velocity is

\[ v(t) = v_0 + \int_0^t a(t')dt', \]

and in particular,

\begin{align} \nonumber v_1 &= v(1) \\ \nonumber &= v_0 + \int_0^1 a(t)dt \\ \label{powerseries2} &= v_0 + \sum_{n=1}^{\infty} \frac{a_n}{n+1}. \end{align}

The position is

\[ x(t) = x_0 + \int_0^t v(t')dt', \]

and so

\begin{align} \nonumber x_1 &= x(1) \\ \nonumber &= x_0 + \int_0^1 v(t)dt \\ \nonumber &= x_0 + \int_0^1 \left(v_0 + \int_0^t a(t')dt'\right) dt \\ \nonumber &= x_0 + v_0 + \int_0^1 \sum_{n=1}^{\infty} \frac{a_n}{n+1}t^{n+1} dt \\ \label{powerseries3} &= x_0 + v_0 + \sum_{n=1}^{\infty} \frac{a_n}{(n+1)(n+2)}. \end{align}

Equations (\ref{powerseries1}), (\ref{powerseries2}), and (\ref{powerseries3}) are enough to satisfy the conditions placed on the motion. The easiest way to proceed is therefore to truncate the power series after the cubic term, so that there are three equations in three unknowns, \(a_1, a_2, a_3\). Empirically, it's possible to solve these equations up to a sign error with a pen and paper during a poorly attended early-year trivia night:

\begin{align*} a_1 &= -36v_0 - 24v_1 - 60(x_0 - x_1), \\ a_2 &= 96v_0 + 84v_1 + 180(x_0 - x_1), \\ a_3 &= -60v_0 + 60v_1 - 120(x_0 - x_1). \end{align*}

(I fixed the sign error.) R code follows. It takes about a second to produce each frame; it would be substantially quicker if smaller pictures were printed.

library(ggplot2)
library(EBImage)
library(grid)

tau = 2*pi

# Subscript 1 refers to start of edge;
# subscript 2 refers to end of edge.
list_pts_1 = list()
list_pts_2 = list()

# -----------------------------------
# First element of list: front spokes

num_pts_set1 = 120
angles = tau * (1:num_pts_set1)/num_pts_set1

radius_set1_1 = 0.25
radius_set1_2 = 0.05
protrusion1 = 0.1
centre_set1 = c(0, 0, 1)

pts_set1_1 = matrix(rep(0, 3*num_pts_set1), nrow = num_pts_set1, ncol = 3)
pts_set1_2 = matrix(rep(0, 3*num_pts_set1), nrow = num_pts_set1, ncol = 3)

pts_set1_1[, 1] = centre_set1[1]
pts_set1_1[, 2] = centre_set1[2] + radius_set1_1*cos(angles)
pts_set1_1[, 3] = centre_set1[3] + radius_set1_1*sin(angles)

pts_set1_2[, 1] = centre_set1[1] + protrusion1
pts_set1_2[, 2] = centre_set1[2] + radius_set1_2*cos(angles)
pts_set1_2[, 3] = centre_set1[3] + radius_set1_2*sin(angles)

# Need to make a full loop so that the geom_segment will be drawn:
pts_set1_1 = rbind(pts_set1_1, pts_set1_1[1, ])
pts_set1_2 = rbind(pts_set1_2, pts_set1_2[1, ])

list_pts_1[[1]] = pts_set1_1
list_pts_2[[1]] = pts_set1_2

# -----------------------------------
# Second element of list: back spokes

num_pts_set2 = 120
angles = tau * (1:num_pts_set1)/num_pts_set2

radius_set2_1 = 0.25
radius_set2_2 = 0.03
protrusion2 = -0.1
centre_set2 = c(0, 0, 1)

pts_set2_1 = matrix(rep(0, 3*num_pts_set2), nrow = num_pts_set2, ncol = 3)
pts_set2_2 = matrix(rep(0, 3*num_pts_set2), nrow = num_pts_set2, ncol = 3)

pts_set2_1[, 1] = centre_set2[1]
pts_set2_1[, 2] = centre_set2[2] + radius_set2_1*cos(angles)
pts_set2_1[, 3] = centre_set2[3] + radius_set2_1*sin(angles)

pts_set2_2[, 1] = centre_set2[1] + protrusion2
pts_set2_2[, 2] = centre_set2[2] + radius_set2_2*cos(angles)
pts_set2_2[, 3] = centre_set2[3] + radius_set2_2*sin(angles)

pts_set2_1 = rbind(pts_set2_1, pts_set2_1[1, ])
pts_set2_2 = rbind(pts_set2_2, pts_set2_2[1, ])

list_pts_1[[2]] = pts_set2_1
list_pts_2[[2]] = pts_set2_2

# -----------------------------------
# Third element of list: stand

num_pts_set3 = 4


pts_set3_1 = matrix(rep(0, 3*num_pts_set3), nrow = num_pts_set3, ncol = 3)
pts_set3_2 = matrix(rep(0, 3*num_pts_set3), nrow = num_pts_set3, ncol = 3)

pts_set3_1[1, ] = c(centre_set2[1] + protrusion2, -radius_set2_2, centre_set2[3])
pts_set3_1[2, ] = c(centre_set2[1] + protrusion2,  radius_set2_2, centre_set2[3])
pts_set3_1[4, ] = c(centre_set2[1] + protrusion2 - radius_set2_2, -radius_set2_2, centre_set2[3])
pts_set3_1[3, ] = c(centre_set2[1] + protrusion2 - radius_set2_2,  radius_set2_2, centre_set2[3])

pts_set3_2 = pts_set3_1
pts_set3_2[, 3] = 0

pts_set3_1 = rbind(pts_set3_1, pts_set3_1[1, ])
pts_set3_2 = rbind(pts_set3_2, pts_set3_2[1, ])

list_pts_1[[3]] = pts_set3_1
list_pts_2[[3]] = pts_set3_2

# -----------------------------------

pts_full_1 = list_pts_1[[1]]
pts_full_2 = list_pts_2[[1]]

num_sets = length(list_pts_1)
if (num_sets > 1) {
  for (i in 2:num_sets) {
    pts_full_1 = rbind(pts_full_1, list_pts_1[[i]])
    pts_full_2 = rbind(pts_full_2, list_pts_2[[i]])
  }
}

num_pts_total = length(pts_full_1[, 1])

# End of points definitions.
# ------------------------------------

# As well as drawing the edge from the start to the end of a 
# spoke, we also want to connect up the inside circle and the
# outside circle.  The draw_segment vector holds values used
# as the alpha for the geom_segment; 1 is opaque, 0 is invisibly
# transparent.  The only values we want zero are when we hop from
# one set of points (say, the front spokes) to the next (say, the
# stand).
draw_segment = rep(1, num_pts_total)
ct = 0
for (i in 1:length(list_pts_1)) {
  this_num_pts = length(list_pts_1[[i]][ , 1])
  ct = ct + this_num_pts
  draw_segment[ct] = 0
}


# Field-of-view limits; adjust as necessary.  Since the stand
# goes all the way to the floor, the fov_y_min is set quite low.
# To ensure that there's good definition in the Moire region,
# this means that a quite large image will be printed, and later
# on we'll use the EBImage package to crop it appropriately.
fov_x = 0.25
fov_y_min = -1
fov_y_max = 0.3


dot_prod = function(v1, v2) {
  return(sum(v1*v2))
}

smoothstep = function(x) {
  return(x*x*(3 - 2*x))
}


num_jumps = 9
frames_per_jump = 40

# Velocities at the "jump" points::
angular_vels = tau/2 - tau*runif(num_jumps)
z_vels = rnorm(num_jumps, 0, 0.3)

vels = unname(cbind(angular_vels, z_vels))
vels = rbind(vels, vels[1, ])

# The velocities are defined at the end of each "jump"; positions will be
# defined as we go; in between will be a smooth interpolation with
# the acceleration equal to zero at each "jump".

# Cylindrical coords (r, theta, z).

# Set the positions at the "jump" points.
# r shouldn't be too close to the fan;
# The angle is near 0 or tau/2 because that is where the best
# Moire patterns are;
# z is kept pretty close to the height of the fan, with more
# leeway at greater distances.

r_posns = 1.2 + 2*runif(num_jumps)
angle_posns = (1 + sign(rnorm(num_jumps,0.5,1)))*tau/2 + rnorm(num_jumps,0,0.15)
z_posns = 1 + r_posns * rnorm(num_jumps, 0, 0.1)

posns = unname(cbind(r_posns, angle_posns, z_posns))
posns = rbind(posns, posns[1, ])


full_frame_count = 1

for (jump_ct in 1:(num_jumps)) {
  r0 = posns[jump_ct, 1]
  r1 = posns[jump_ct+1, 1]
  
  x0 = posns[jump_ct, 2:3]
  x1 = posns[jump_ct + 1, 2:3]
  
  v0 = vels[jump_ct]
  v1 = vels[jump_ct + 1]
  
  # acceleration a(t) = a_1*t + a_2*t^2 + a_3*t^3
  a_1 = -36*v0 - 24*v1 - 60*(x0-x1)
  a_2 = 96*v0 + 84*v1 + 180*(x0-x1)
  a_3 = -60*(v0+v1) - 120*(x0-x1)
  
  
  for (frame_ct in 1:frames_per_jump) {
    t = frame_ct / frames_per_jump
    angle_z = x0 + v0*t + a_1*t^3/6 + a_2*t^4/12 + a_3*t^5/20
    
    r = r0 + smoothstep(t)*(r1 - r0)
    theta = angle_z[1]
    z = angle_z[2]
    
    viewer_posn = c(r*cos(theta), r*sin(theta), z)
    
    # I have the observer looking at the centre of the fan,
    # but you might want to change it:
    looking_at = c(0, 0, 1)
    
    # Rotation angles for the transformed axes:
    axis_r = looking_at - viewer_posn
    axis_dist = sqrt(sum(axis_r*axis_r))
    alpha = atan2(axis_r[2], axis_r[1])
    beta = asin(axis_r[3]/axis_dist)
    
    # Transform the axis vectors:
    i1 = c(1, 0, 0)
    j1 = c(0, 1, 0)
    k1 = c(0, 0, 1)
    
    i2 = i1*cos(alpha) + j1*sin(alpha)
    j2 = -i1*sin(alpha) + j1*cos(alpha)
    k2 = k1
    
    i3 = i2*cos(beta) + k2*sin(beta)
    j3 = j2
    k3 = -i2*sin(beta) + k2*cos(beta)
    
    # -------------------------------
    v_1 = pts_full_1
    v_2 = pts_full_2
    
    for (i in 1:3) {
      v_1[, i] = v_1[, i] - viewer_posn[i]
      v_2[, i] = v_2[, i] - viewer_posn[i]
    }
    
    # Normalise:
    v_1 = v_1 / apply(v_1, 1, function(x) sqrt(sum(x*x)))
    v_2 = v_2 / apply(v_2, 1, function(x) sqrt(sum(x*x)))
    
    phi_1 = apply(v_1, 1, function(x) asin(dot_prod(x, k3)))
    phi_2 = apply(v_2, 1, function(x) asin(dot_prod(x, k3)))
    
    theta_1_interm = apply(v_1, 1, function(x) dot_prod(x, j3))
    theta_1 = asin(theta_1_interm / cos(phi_1))
    
    theta_2_interm = apply(v_2, 1, function(x) dot_prod(x, j3))
    theta_2 = asin(theta_2_interm / cos(phi_2))
    
    
    persp_plot = ggplot() + coord_equal() +
      geom_segment(aes(x=theta_1, y=phi_1,
                       xend=theta_2, yend=phi_2)) +
      geom_segment(aes(x=theta_1, y=phi_1,
                       xend=c(theta_1[2:num_pts_total], theta_1[1]),
                       yend=c(phi_1[2:num_pts_total], phi_1[1])), alpha=draw_segment) +
      geom_segment(aes(x=theta_2, y=phi_2,
                       xend=c(theta_2[2:num_pts_total], theta_2[1]),
                       yend=c(phi_2[2:num_pts_total], phi_2[1])), alpha=draw_segment) +
      xlim(-fov_x, fov_x) + ylim(fov_y_min, fov_y_max) + 
      theme(line = element_blank(),
            text = element_blank(),
            line = element_blank(),
            title = element_blank(),
            rect = element_blank())
    
    # I wanted to extract just the plot region of the ggplot,
    # since there's no point having any margins.  Not sure if
    # it makes much difference though....
    gt <- ggplot_gtable(ggplot_build(persp_plot))
    ge <- subset(gt$layout, name == "panel")
    
    outfile <- sprintf("images/%03d.png", full_frame_count)
    png(filename=outfile, width=800, height=1200, units="px")
    grid.draw(gt[ge$t:ge$b, ge$l:ge$r])
    # print(persp_plot)
    dev.off()
    
    # Extract a 500x500-pixel region:
    temp_img = readImage(outfile)
    temp_img = temp_img[150:649, 60:559, ]
    writeImage(temp_img, outfile)
    
    print(sprintf("Done frame %d.", full_frame_count))
    full_frame_count = full_frame_count + 1
  }
}

Posted 2014-01-10.

Home > Misc