Simple pupile function calculations

Note: This is a remake of my original post on pupil functions. I decided to write this remake to fix what I perceived to be some errors in thinking and to better address the problem of units in some of the physical quantitities. For the original, click here.

A pupil function is a theoretical tool for characterizing an imaging system. In simple terms, it is a mathematical model for any general arrangement of lenses and mirrors used to collect light from an object and form an image of that object in some other plane. A few of the reasons why pupil functions are useful include:

  1. they reduce a complicated optical system--such as a microscope--to a relatively simple, two-dimensional, and complex-valued function;
  2. they provide a convenient way to represent the aberrations present in the system;
  3. they are easy to simulate on a computer using fast Fourier transforms.

In this post I will show you how to write a simple program that computes the image of an isotropic point source of light from an optical system's pupil function.

Theoretical background

In scalar diffraction theory, light is represented as a three-dimensional function known as the scalar field. At every point in space $ \mathbf{r} $, the scalar field $ u \left( \mathbf{r} \right) $ is a single, complex value that represents the electric field at that point. The physical laws of diffraction require that the scalar field be described by two numbers, an amplitude and a phase, that are derived from the field's real and the imaginary parts, $ \text{Re} \left[ u \left( \mathbf{r} \right) \right] $ and $ \text{Im} \left[ u \left( \mathbf{r} \right) \right] $:

\begin{align*} A &= \sqrt{\text{Re} \left[ u \left( \mathbf{r} \right) \right]^2 + \text{Im} \left[ u \left( \mathbf{r} \right) \right]^2 } \\ \phi &= \arctan \left( \frac{\text{Im} \left[ u \left( \mathbf{r} \right) \right]}{\text{Re} \left[ u \left( \mathbf{r} \right) \right]} \right) \end{align*}

If we know the amplitude and phase at a given point, then we know the scalar field at that point. Despite the fact that scalar diffraction theory ignores the polarization of light, it does wonderfully well at describing a large range of optical phenomena.

For most problems in imaging, we don't really need to know the three-dimensional distribution of the field in all of space. Instead, we simplify the problem by asking how an optical system transforms the field in some two-dimensional object plane into a new field distribution in the image plane (see the figure below). Any changes in scale between these two planes are caused by the system's magnification; any blurring or distortion is caused by diffraction and possibly aberrations.

Optical system with pupils

The pupil function is the two-dimensional Fourier transform of the scalar field in the image plane when the object is a point source emitting light equally in all directions, i.e. an isotropic point source. Mathematically, the pupil function is written as

\begin{equation*} P \left(f_x, f_y \right) = \frac{1}{A}\iint_{-\infty}^{\infty} \text{PSF}_A \left( x, y \right) \exp \left[ -j 2 \pi \left( f_x x + f_y y\right) \right] \, dx \, dy \end{equation*}

where $ A $ is a normalizing constant, $ f_x $ and $ f_y $ represent spatial frequencies in the x- and y-directions, and $ j $ is the imaginary number. $ \text{PSF}_A \left( x, y \right) $ is known as the amplitude point spread function. Despite the intimidating name, $ \text{PSF}_A \left( x, y \right) $ is just the scalar field from the isotropic point source in the image plane. The pupil function and the amplitude point spread function form a Fourier transform pair, so we can also write

\begin{equation*} \text{PSF}_A \left(x, y \right) = A \iint_{-\infty}^{\infty} P \left( f_x, f_y \right) \exp \left[ j 2 \pi \left( f_x x + f_y y\right) \right] \, df_x \, df_y \end{equation*}

What all of this means is that we can compute the image of an on-axis, isotropic point source if we know the pupil function that describes the system: compute the two-dimensional Fourier transform of the pupil function and voilà, you have the image (or at least the field that will form the image).

The pupil punction is dimensionless; $ \text{PSF}_A $ is a field

By convention the pupil function is a dimensionless complex number and has a magnitude between 0 and 1. The amplitude PSF, however, has units of electric field, or Volts per distance, $ V / m $.

If you perform a dimensional analysis on the Fourier transform expressions above, you can see that the normalizing constant $ A $ has to have units of $ V \times m $. For example, if $ \text{PSF}_A $ has units of $ V / m $ and $ dx $ and $ dy $ both have units of $ m $, then $ A $ has units of $ V \times m $ and the pupil function is therefore dimensionless. Sometimes in the literature you will find that units or a normalizing constant are ignored or, worse, that the reader can sort of "fill them in" later. I prefer to be explicit and to define the $ \text{PSF}_A $ to be a field.

Pupil functions are not entrance or exit pupils

The pupil function and the pupil planes of a system are not the same thing. The entrance and exit pupils--which together are known as the pupil planes--are the planes in which the images of the system's aperture stop are located. The pupil function, however, is not an image of the aperture stop of the system; it's a 2D Fourier transform of a field.

There is never-the-less a relationship between the pupil function and the plane of the exit pupil. The pupil function represents the relative amplitude and phase of the field on the surface of a so-called reference sphere that intersects the optics axis in the plane of the system's exit pupil.

Pupil function simulations in python

The goal of this simulation will be simple: given a pupil function, a single wavelength, an optical system with a numerical aperture NA, and an amount of power that passes through the image plane, compute the image of an on-axis isotropic point source.

There are only a few steps needed to achieve our goal:

  1. define the simulation's input parameters;
  2. setup the image plane and pupil plane coordinate system;
  3. create the pupil plane and normalize it so that the field carries the desired amount of power;
  4. and compute the field in the image plane.

Before we go further, it's worth pointing out that the pupil function and $ \text{PSF}_A $ are obtained by calculating the continuous Fourier transform of one another. On a computer, however, it's often easiest to compute what's called a discrete Fourier tranform via the fast Fourier transform (FFT) algorithm. The continuous Fourier transform and the FFT are not, strictly speaking, the same thing. Therefore, we should expect from the start that there may be small differences between the computed $ \text{PSF}_A $ and the analytical calculation.

With this in mind, we'll start by importing a few scientific libraries like Numpy and Scipy.

In [1]:
%pylab inline
import sys
from numpy.fft import fft2, fftshift
import scipy
from scipy.integrate import simps
import seaborn as sns # Used only to set up plots
sns.set_context(context = 'talk')'dark_background')
plt.rcParams['figure.facecolor'] = '#272b30'
plt.rcParams['image.cmap'] = 'viridis' 

print('Python version:\n{}\n'.format(sys.version))
print('Numpy version:\t\t{}'.format(np.__version__))
print('matplotlib version:\t{}'.format(matplotlib.__version__))
print('Scipy version:\t\t{}'.format(scipy.__version__))
print('Seaborn version:\t{}'.format(sns.__version__))
Populating the interactive namespace from numpy and matplotlib
Python version:
3.5.2 |Continuum Analytics, Inc.| (default, Jul  2 2016, 17:53:06) 
[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]

Numpy version:		1.11.2
matplotlib version:	1.5.3
Scipy version:		0.18.1
Seaborn version:	0.7.1

Step 1: Define the input parameters

Next, we need to define a few parameters that will determine the output of the simulations. These parameters are:

  1. wavelength Units are $ \mu m $.
  2. NA Numerical aperture of the system. No units.
  3. pixelSize The length of a square pixel in the object space. Units are $ \mu m $.
  4. numPixels The number of pixels in your camera. This will be assumed to be even.
  5. power The total power carried by the field in Watts, $ W $.

Note that pixel size is defined as the size of a pixel in the object space. Defining it like this is often more intuitive than defining the pixel size in the image space. For example, when you see a scale bar on a microscopy image, the distances correspond to object space distances, not the actual distance that the image spans on the piece of paper or computer screen. Furthermore, there is no problem with working with an object space pixel size in the image plane since the coordinate systems in the object and image planes of our perfect optical system can be easily mapped onto one another by a linear scaling by the system's magnification.

We don't neccessarily need to use a camera as the detector. Since we are limited to working with discrete arrays in the computer, though, it's convenient to say that we have a camera as a detector since each pixel is a discrete sample of the field.

In addition to the above parameters, we'll assume that the object, the imaging system, and the image plane are all in air. We'll define a constant $ Z_0 = 376.73 \, \Omega $ which is known as the impedance of free space or the vacuum impedance. This is the constant of proportionality between the power carried by the scalar field and the integral of its absolute square in the image plane:

\begin{equation*} P_0 = \frac{1}{Z_0} \iint_{-\infty}^{\infty} \left| \text{PSF}_A \left( x, y \right) \right|^2 \, dx \, dy \end{equation*}

Of course, air does not really have the same impedance as vacuum, but the two values are close enough. I have also used $ P_0 $ to denote the power because I already used $ P \left( f_x, f_y \right) $ to denote the pupil function.

In [2]:
# Setup the simulation parameters
wavelength = 0.68   # microns
NA         = 0.7    # Numerical aperture of the objective
pixelSize  = 0.1    # microns
numPixels  = 2048   # Number of pixels in the camera; keep this even
power      = 0.1    # Watts
Z0         = 376.73 # Ohms; impedance of free space

Step 2: Create the image and pupil plane coordinate systems

Our simulation will transform the values in a two-dimensional square array of complex numbers (the pupil function) into a new two-dimensional array of complex numbers of the same size. Before we do this, however, let's first determine the coordinates of each pixel.

Since we specified the pixel size and the number of pixels, it's easiest to start with the image plane coordinates. We will define the origin of our coordinate system to lie at the center of the array, which, for an even number of pixels and array indexes that start at zero, lies halfway between the pixels $ \left( \frac{\text{numPixels}}{2} \right) - 1 $ and $ \left( \frac{\text{numPixels}}{2} \right) $ in both the horizontal and vertical directions.

In [3]:
# Create the image plane coordinates
x = np.linspace(-pixelSize * numPixels / 2, pixelSize * numPixels / 2, num = numPixels, endpoint = True)

We only need to create a single, one-dimensional array to represent the coordinates because our image and pupil plane arrays will be square; we can use the same array to represent the coordinates in both the horizontal and vertical directions.

With the image plane coordinates taken care of, the next question is: what are the pupil plane coordinates? This question is a frequent source of frustration for students (and full-time scientists). I won't go into the details in this post, but instead will just tell you the two rules you need to remember for Fourier optics simulations

  1. The number of elements in the pupil function array is the same as the number of elements in the image plane array.
  2. For an even number of array elements, the frequency values along either principle direction in the pupil function run from $ -\frac{f_S}{2} $ to $ f_S \left( \frac{1}{2} - \frac{1}{\text{numPixels}} \right) $ with the spacing between discrete coordinate values equal to $ \frac{f_S}{\text{numPixels}} $.

$ f_S $ is called the sampling frequency and is equal to one divided by the spacing between image space coordinates. We can go ahead now and compute the frequency-space coordinate values.

In [4]:
# Create the Fourier plane
dx = x[1] - x[0]    # Sampling period, microns
fS = 1 / dx         # Spatial sampling frequency, inverse microns
df = fS / numPixels # Spacing between discrete frequency coordinates, inverse microns
fx = np.arange(-fS / 2, fS / 2, step = df) # Spatial frequency, inverse microns

Step 3: Create the pupil function

In nearly all imaging systems, the pupil is circular because its optical elements are circular. The radius of the pupil function is the ratio of the system's numerical aperture to the wavelength of the light, $ \frac{\text{NA}}{\lambda} $ (Hanser, 2004). Perfect systems like the one we are modeling here have a pupil with a constant value everywhere inside this circle and zero outside of it.

We can simulate such a pupil by making a circular mask with a radius of $ \frac{\text{NA}}{\lambda} $. The mask is one inside the circle and zero outside of it.

In [5]:
# Create the pupil, which is defined by the numerical aperture
fNA             = NA / wavelength # radius of the pupil, inverse microns
pupilRadius     = fNA / df        # pixels
pupilCenter     = numPixels / 2   # assumes numPixels is even
W, H            = np.meshgrid(np.arange(0, numPixels), np.arange(0, numPixels)) # coordinates of the array indexes
pupilMask       = np.sqrt((W - pupilCenter)**2 + (H - pupilCenter)**2) <= pupilRadius

Define the power carried by the scalar field

I mentioned in the theoretical background above that the total optical power carried by the field is the two dimensional integral of the absolute square of the field divided by the impedance. If we want to set the power as an input of the simulation, we need to first normalize our pupil values by this integral.

Parseval's theorem tells us that we can integrate over $ \text{PSF}_A \left( x, y \right) $ and $ A \times P \left( f_x, f_y \right) $, i.e. the image plane field or the normalized pupil, respectively, and get the same number:

\begin{equation*} \iint_{-\infty}^{\infty} \left| \text{PSF}_A \left( x, y \right) \right|^2 \, dx \, dy = A^2 \iint_{-\infty}^{\infty} \left| P \left( f_x, f_y \right) \right|^2 \, df_x \, df_y \end{equation*}

Now that we have the pupil, we can perform a numerical integration over it using Simpson's rule to find the normalizing constant. We then multiply the pupil by the square root of this constant times our desired value for the power to set the total power carried by the field.

In [6]:
# Compute normalizing constant
norm_factor = simps(simps(np.abs(pupilMask)**2, dx = df), dx = df) / Z0 # A^2
print('Normalization constant:\t\t{:.4f} W'.format(np.sqrt(norm_factor)))

# Renormalize the pupil values
pupil = pupilMask * (1 + 0j)
normedPupil = pupil * np.sqrt(power / norm_factor)

new_power = simps(simps(np.abs(normedPupil)**2, dx = df), dx = df) / Z0
print('User-defined power:\t\t{:.4f} W'.format(power))
print('Power now carried by field:\t{:.4f} W'.format(new_power))
Normalization constant:		0.0940 W
User-defined power:		0.1000 W
Power now carried by field:	0.1000 W
In [7]:
# Show the pupil
ax = plt.imshow(np.abs(pupil), extent = [fx[0], fx[-1], fx[0], fx[-1]])
cb = plt.colorbar(ax)
cb.set_label('$P \, ( f_x, f_y ) $')
plt.xlabel('$f_x$, $\mu m^{-1}$')
plt.ylabel('$f_y$, $\mu m^{-1}$')

# Compute the power
power_pupil = simps(simps(np.abs(normedPupil)**2, dx = df), dx = df) / Z0
print('Power in pupil plane: {:.4f} W'.format(power_pupil))
Power in pupil plane: 0.1000 W

Step 4: Compute the image of the point source

With the pupil and coordinate systems established, we are now ready to compute the image of the isotropic point source that this system produces.

To do this, we need to perform a few easy but important steps. In the first step, we will shift the origin of the pupil from the center of the array to the array indexes $ \left( 0, 0 \right) $ using the ifftshift() function. The reason we do this is that fft2() expects that the zero frequency value lies at the origin of the array. The two-dimensional FFT of the shifted pupil is then computed, producing a new array with the zero frequency at array indexes $ \left( 0, 0 \right) $ and the Nyquist frequency $ f_S / 2 $ in the middle of the array's axes (numpy.fft.fft2 - NumPy v1.11 Manual, accessed on 2016-11-15). We finish by shifting the origin back to the center of the array using fftshift() so that it makes sense when we visualize the results.

The final step is to multiply the result by the square of the spacing between frequency coordinates. This step ensures that the power is preserved during the FFT operation (Schmidt, 2010, pp. 15-18).

Chaining the functions together, these steps look like:

psf_a = fftshift(fft2(ifftshift(normedPupil))) * df**2

The image is the irradiance in the image plane, which is the absolute square of the field, divided by the vacuum impedance. Its units are power per area, or in this case $ W / \mu m^2 $. In optics, the irradiance is what is measured by a camera, not the field.

image = np.abs(psf_a)**2 / Z0
In [8]:
psf_a = fftshift(fft2(ifftshift(normedPupil))) * df**2
image = np.abs(psf_a)**2 / Z0

# Show the image plane
img = plt.imshow(image, interpolation='nearest', extent = [x[0], x[-1], x[0], x[-1]])
cb = plt.colorbar(img)
plt.gca().set_xlim((-2, 2))
plt.gca().set_ylim((-2, 2))

plt.xlabel('x, $\mu m$')
plt.ylabel('y, $\mu m$')
cb.set_label('Irradiance, $W / \mu m^2$')

Above you can see the image of an isotropic point source. The image is not a point but rather a blurred spot in the center of the image plane due to diffraction at the pupil.

Verify the results

The first thing we can do to check whether the above result is correct is compute the power over the image above. By Parseval's theorem, it should be the same as the integral of the normalized pupil function divided by the vacuum impedance.

In [9]:
power_image = simps(simps(image, dx = dx), dx = dx)
print('Power in pupil plane: {:.4f} W'.format(power_pupil))
print('Power in object plane: {:.4f} W'.format(power_image))
Power in pupil plane: 0.1000 W
Power in object plane: 0.1000 W

So far so good. The next thing that we can do to verify these results is to calculate the sampled values of the analytical solution to this problem. Scalar diffraction theory predicts that the solution for the irradiance of the field diffracted by a circular aperture is an Airy disk:

\begin{equation*} I \left( r \right) = I_0 \left[ \frac{2 J_1 \left( X \right)}{X} \right]^2 \end{equation*}

where $X = \frac{2 \pi r \, \text{NA}}{\lambda} $, $ r = \sqrt{x^2 + y^2} $ is the radial coordinate, and $ J_1 \left( r \right) $ is called the first-order Bessel function of the first kind (Weisstein, Mathworld, accessed on 2016-11-16). This function does not exist inside Python's scientific libraries, so we will need to create it.

In [10]:
from scipy.special import j1 as bessel1
def airyDisk(x,y, NA = 0.5, wavelength = 0.532):
    """Computes a 2D airy disk pattern.
    x, y : array of int, array of int
        Coordinates where the function will be evaluated.
    NA   : float
        The system's numerical aperture.
    wavelength: float
        The wavelength of the light; same units as x and y.
    result : array of float
    r      = np.sqrt(x**2 + y**2)
    X      = 2 * np.pi * r * NA / wavelength
    result = (2 * bessel1(X) / X)**2
        # Replace value where divide-by-zero occurred with 1
        result[np.logical_or(np.isinf(result), np.isnan(result))] = 1
    except TypeError:
        # TypeError is thrown when single integers--not arrays--are passed into the function
        result = np.array([result])
        result[np.logical_or(np.isinf(result), np.isnan(result))] = 1
    return result

Now we can go ahead and visually compare our image plane calculations with the airy disk. If we subtract one from the other, we should get all zeros.

In [11]:
from mpl_toolkits.axes_grid1 import make_axes_locatable
In [12]:
# Subtraction by dx/2 places the origin at the edge of a pixel, not a center
X, Y = np.meshgrid(x - dx/2, x - dx / 2, indexing = 'xy')

fig, ax = plt.subplots(nrows = 1, ncols = 2, figsize = (12, 8))

img = ax[0].imshow(image, interpolation='nearest', extent = [x[0], x[-1], x[0], x[-1]])

divider = make_axes_locatable(ax[0])
cax = divider.append_axes("right", size="5%", pad=0.05)
cb0 = plt.colorbar(img, cax = cax)

ax[0].set_xlim((-2, 2))
ax[0].set_ylim((-2, 2))

ax[0].set_xlabel('x, $\mu m$')
ax[0].set_ylabel('y, $\mu m$')

I0  = np.max(image)
img = ax[1].imshow(I0 * airyDisk(X,Y, NA = NA, wavelength = wavelength), interpolation = 'nearest', extent = [x[0], x[-1], x[0], x[-1]])

divider = make_axes_locatable(ax[1])
cax  = divider.append_axes("right", size="5%", pad=0.05)
cb1  = plt.colorbar(img, cax = cax)

ax[1].set_xlim((-2, 2))
ax[1].set_ylim((-2, 2))

ax[1].set_xlabel('x, $\mu m$')

cb1.set_label('Irradiance, $W / \mu m^2$')

/home/kmdouglass/anaconda3/envs/nikola/lib/python3.5/site-packages/ipykernel/ RuntimeWarning: invalid value encountered in true_divide
In [13]:
# Find and plot the difference between the simulated image and the theoretical one
I0 = np.max(image)
diffImg = image - I0 * airyDisk(X,Y, NA = NA, wavelength = wavelength)

plt.imshow(diffImg, interpolation = 'nearest', extent = [x[0], x[-1], x[0], x[-1]])
cb = plt.colorbar()
cb.set_label('Irradiance, $W / \mu m^2$')
plt.xlabel('x $\mu m$')
plt.xlim((-2, 2))
plt.ylabel('y $\mu m$')
plt.ylim((-2, 2))
plt.title('Difference between simulation and theory')
/home/kmdouglass/anaconda3/envs/nikola/lib/python3.5/site-packages/ipykernel/ RuntimeWarning: invalid value encountered in true_divide

From the plots above you can see that the simulations perform pretty well at finding the image of the point source. In fact, the differences in the final irradiance values differ by at most a small fraction of a percent.

However, the pattern in the difference image is not random, so these differences are probably not round-off errors. In fact, they come from the minor but important detail that we are calculating a discrete Fourier transform with the FFT command, whereas the theory predicts that the point spread function is a continuous Fourier transform of the pupil.

Final remarks

The calculations I discussed here rely on quite a few simplifying assumptions. The first is that they are based on scalar diffraction theory. The real electromagnetic field is polarized and is described not by a scalar but by a three dimensional, complex-valued vector. The theory of vector diffraction describes how such a field propagates through an optical system; naturally, it is more complicated than scalar diffraction theory.

Isotropic point sources are also an idealization of reality. Single fluorescent molecules, for instance, do not emit light equally in all directions but rather in a pattern known as dipole radiation.

Finally, the optical system as modeled here was a linear, shift-invariant system. This means that the same pupil function can actually be used to compute the image of any off-axis point source as well; we would only need to apply a linear phase ramp to the pupil to account for this. Real optical systems are only approximatley linear and shift-invariant, which means that a more complete model would assign a different pupil function to each object point. This then allows for aberrations that vary as a function of object position, but requires more information to fully describe the system.

Field-dependent photophysics - video

Below you can see a video of a raw STORM movie that I acquired in the lab a few months ago on a standard, inverted laser fluorescence microscope. The sample consists of AlexaFluor 647-labeled tubulin in Cos7 cells at a frame rate of 100 fps. (I have slowed things down a bit to more easily see what's going on.) Go ahead and watch the video now, I'll discuss a little bit about what you're seeing aftwards.

In the beginning of the movie you can see the moment when the laser beam first hits the sample. Since the illumination geometry is just widefield epi-illumination, the laser beam has a Gaussian-like irradiance profile. (As explained in the link, irradiance is the amount of energy the laser beam carries per differential unit area. It is often referred to as intensity, though intensity actually refers to a different quantity in the field of radiometry.) This means that the fluorophores in the middle of the beam are exposed to a higher irradiance than those at the edges.

The effect of the laser beam profile for STORM is quite clear when watching the video. The photoswitching rates of the fluorophores to the long-lived dark state are much higher in the center of the beam than at the edges. This leads to a field-dependent density of molecules in the fluorescence emitting state, which will ultimately affect one's ability to precisely and accurately localize the single molecule emissions. If you look closely, you can also see that the brightness of the fluorophores is also higher near the beam center, meaning that these molecules will be localized more precisely because we recorded more photons from them.

All of these features mean that one must take care when performing quantitative analyses on STORM datasets, because the density, precision, and accuracy of the single molecule position estimates will vary across the field of view.

Perlin Noise

Since my last post, I have been experimenting with the noise package from Casey Duncan. This package is the easy solution for generating Perlin noise that I had hoped for. It seems simple and fast, though one downside at least to me is that it doesn't play so nicely with numpy arrays. This is not a big deal, though. It was probably not developed for scientists but rather for game design and/or computer graphics.

So, let's use this library to continue exploring the question of why correlated noise often makes for a better model of the real world.

In [1]:
%pylab inline
import seaborn as sns
from matplotlib import animation
from noise import pnoise1, pnoise3
from mpl_toolkits.axes_grid1 import make_axes_locatable
plt.rcParams['savefig.dpi'] = 120
Populating the interactive namespace from numpy and matplotlib

Generating Perlin noise with the noise library

The basic function pnoise1() takes a single float argument x. It won't work with numpy arrays. For this reason, I'm generating the noise signal in a for loop.

In [3]:
numPoints = 100
time      = np.linspace(0, 10, num = numPoints)

# Generate the signal one element at a time
signal    = np.zeros(numPoints)

for ctr, t in enumerate(time):
    signal[ctr] = pnoise1(t)

plt.plot(time, signal)

The octaves and persistence parameters determine how much energy is contained at frequencies in each successive octave of the noise's power spectrum. These parameters determine the presence and strength, respectively, of the finer details in the noise signal. Before going into them with more detail, we can vary them and observe the outcome.

Here's the result of varying the number of octaves.

In [4]:
signal2 = np.zeros(numPoints)
signal3 = np.zeros(numPoints)
for ctr, t in enumerate(time):
    signal[ctr]  = pnoise1(t, octaves = 1) + 1
    signal2[ctr] = pnoise1(t, octaves = 2)
    signal3[ctr] = pnoise1(t, octaves = 4) - 1

plt.plot(time, signal,  label = '1 octave',  linewidth = 2)
plt.plot(time, signal2, label = '2 octaves', linewidth = 1.5)
plt.plot(time, signal3, label = '4 octaves', linewidth = 1.0)

In the above plot, and besides constant offsets to aid in visualization, I added one and three additional octaves to the noise, keeping the persistence fixed. Increasing the octaves makes the curve appear more 'jagged' and noisy on the smaller scales, much like extending the frequency range of the filter applied to the random signal from my last post. Over longer time scales, the curve more or less follows the original one.

Below you can see what happens when we use four octaves but increase the persistence gradually. The amplitude of the spikes at small time scales increases with increasing persistence.

In [5]:
for ctr, t in enumerate(time):
    signal[ctr]  = pnoise1(t, octaves = 4, persistence = 0.25) +1
    signal2[ctr] = pnoise1(t, octaves = 4, persistence = 0.5)
    signal3[ctr] = pnoise1(t, octaves = 4, persistence = 1) - 1

plt.plot(time, signal,  label = 'persistence = 0.25',  linewidth = 2)
plt.plot(time, signal2, label = 'persistence = 0.50', linewidth = 1.5)
plt.plot(time, signal3, label = 'persistence = 1.00', linewidth = 1.0)

We're now in a good position to summarize a few of the main characteristics about Perlin noise. I'll use the term 'feature' as a loose definition for the hills and valleys seen in the curves.

  1. Perlin noise is a random signal that is correlated over short time scales.
  2. The number of octaves determines the time/length scales at which features in the signal are present.
  3. The persistence determines the amplitude of the features.

The existence of features at smaller-and-smaller scales and how strongly these features decrease in amplitude with scale are the fractal-like properties of the noise signal. So Perlin noise is essentially a way of generating random fractals on a computer. Unlike the deterministic fractals that we normally think of, these are statistical.

In the final parts of this post, let's see how we can use correlated noise to model realistic motions of particles.

Random waveforms in time

Because Perlin noise is continuous, we can easily generate a nice-looking, traveling Perlin noise waveform and animate it. To do this, I use FuncAnimation from matplotlib's animation API. Every time FuncAnimation is called, I add a constant to the time array, which shifts the curve to the left. I'm also plotting a small particle undergoing a 1D walk and whose position is equivalent to the waveform value at the left-most side of the axis.

In [6]:
# Generate initial axis for the waveform
numPoints = 1000
time      = np.linspace(0, 10, num = numPoints)
signal    = np.zeros(numPoints)
fig, ax = plt.subplots(nrows = 1, ncols = 1)
line,   = ax.plot([], [])
ax.set_xlim((np.min(time), np.max(time)))
ax.set_ylim((-1, 1))

# Generate initial axis for the particle
divider   = make_axes_locatable(ax)
pax       = divider.append_axes('left', size = '5%', pad = 0.05)
particle, = pax.plot([], [], 'o')
pax.set_xlim((-0.5, 0.5))
pax.set_ylim((-1, 1))

# Initialize the figure with an empty, single tuple of line
def init():
    return line, particle

# This function is called for each frame of the animation
def animate(frame):
    dt = 0.01 * frame
    for ctr, t in enumerate(time):
        signal[ctr] = pnoise1(t + dt, octaves = 4, persistence = 0.5)
    # Update waveform and particle
    line.set_data(time, signal)
    particle.set_data(np.array([0]), signal[0])
    return line, particle

# Create the animation
anim = animation.FuncAnimation(fig,
                               frames = 1000,
                               init_func = init,
                               interval = 20,
                               blit = True)

# Save the animation
myWriter = animation.FFMpegWriter(fps = 30, extra_args=['-vcodec', 'libx264'])'traveling_Perlin_wave.mp4', writer = myWriter)

Although random, we can see that the particle moves smoothly and continiously. If I were add more octaves or increase the persistence of the walk, the particle would move more erraticly because I would have introduced more high frequency noise like the pink noise example in my previous post.

Random Perlin and Brownian Walks

Finally, let's compare a random walk created using Perlin noise to one created with white noise, which is completely uncorrelated. A walk generated with white noise is the same as good ol' Brownian motion. By comparing the two, you should be able to see why Perlin noise often better approximates natural phenomena: things just don't look quite natural when modeled with white noise.

In [7]:
# Initialize the axis and two different walk trajectories
fig, ax = plt.subplots(nrows = 1, ncols = 1)
pWalk, = ax.plot([], [], 'o') # Perlin walk
bWalk, = ax.plot([], [], 'ro') # Brownian walk
ax.set_xlim((-100, 100))
ax.set_ylim((-100, 100))
ax.legend([pWalk, bWalk], ['Perlin walker', 'Brownian walker'])

start, stop, step = 0, 100, 0.01
time = np.arange(start, stop, step)

# Create arrays for the x and y positions of the Perlin walker
px   = np.zeros(len(time))
py   = np.zeros(len(time))
for ctr, t in enumerate(time):
    # offset of 1000 ensures x and y are uncorrelated
    px[ctr] = pnoise1(t,        octaves = 4, persistence = 0.5)
    py[ctr] = pnoise1(t + 1000, octaves = 4, persistence = 0.5)
px = np.cumsum(px)
py = np.cumsum(py)

# The Brownian walker
bx = np.random.randn(len(time))
by = np.random.randn(len(time))
bx = np.cumsum(bx)
by = np.cumsum(by)

def init():
    return line, particle

# This function is called for each frame of the animation
def animate(frame):
    # Update the walker positions
    pWalk.set_data(px[frame], py[frame])
    bWalk.set_data(bx[frame], by[frame])
    return pWalk, bWalk

# Create the animation
anim = animation.FuncAnimation(fig,
                               frames = 1000,
                               init_func = init,
                               interval = 20,
                               blit = True)

# Save the animation
myWriter = animation.FFMpegWriter(fps = 30, extra_args=['-vcodec', 'libx264'])'walkers.mp4', writer = myWriter)


The Perlin walker is just a calm, casual guy out for a walk, whereas the Brownian walker is your officemate after five cups of coffee.

So there you have it! Perlin noise is an interesting tool for generating correlated noise forms on a computer. It's used in a number of applications in computer graphics, animation, and games to produce realistic-looking patterns.

I only mentioned it briefly, but one of the reasons its patterns look more realistic than white noise forms is that it has fractal properties. In a later post, I'll explore what I mean by this by using the FFT tools and the pink noise example I generated in my last post.

Correlated noise and the FFT

I've been perusing the online coding book The Nature of Code recently. So far I really like it, especially its approach to demonstrating how one may use mathematics and code to model natural phenomena in a fun a way. I thought I would take a brief break from optics-based posts to blog about some thoughts I've had recently about one of the concepts in the book.

In the introduction I came across the concept of Perlin noise, a concept which is named after its inventor Ken Perlin who, at the time of its creation, was working on computer-generated graphics for the movie Tron. Perlin developed his namesake algorithm for generating artificial yet natural-looking textures, something that previously was difficult to do with purely deterministic mathematical tools. He later won an Academy Award for this work because it greatly enabled computer artists to model realistic textures and motion.

In [1]:
%pylab inline
from numpy.fft import fft, fftshift, ifft, ifftshift'dark_background')
Populating the interactive namespace from numpy and matplotlib

Uncorrelated noise

One reason that Perlin noise models many natural phenomena so well is that its randomness is correlated. To see what I mean by this, I'll first plot a single realization of a purely random signal composed of random integers between 0 and 1. The samples from this signal are not correlated.

In [2]:
numSamples = 100

x       = np.arange(numSamples)
samples = np.random.rand(numSamples)

plt.plot(x, samples, '.', markersize = 10)
plt.xlabel('Sample number')

In the plot above, you can see that the value of each sample from our random signal is uncorrelated with all the other samples. This situation is similar to flipping a coin; the result of a coin flip (heads or tails) doesn't depend on the previous or next flip of the coin.

If we want to use this random signal to model something natural, such as a landscape, we would be disappointed:

In [3]:
plt.plot(x, samples, linewidth = 2)
plt.xlabel('Sample number')

Our purely random landscape is much too jagged to even come close to approximating a real-looking landscape.

Correlated noise

Perlin noise is a way of generating a random signal whose samples are correlated over short distances but uncorrelated over longer distances. I haven't yet found any quick and easy implementations for generating Perlin noise in Python, but you can find many examples of on the web, such as 1D Perlin noise images from Google Image search. (By "quick and easy", I mean something I can install from Anaconda and use with only a couple minutes of reading the docs.)

A long time ago I wrote a computer simulation that modeled the motion of a small dipole in a three-dimensional, random electromagnetic field. This involved moving the dipole in a ballistic manner over short times, where the particle's mean-squared-displacement (MSD) evolved quadratically with time, and smoothly rounding off to a random, diffusive motion characterized by a linear time dependence of the MSD at long times. I have also done some partial coherence modeling of illumination systems, which also involves generating correlated noise.

The point is that I have often needed to generate correlated noise on a computer, but had never come across the concept of Perlin noise until now. Traditionally, I have used fast Fourier transform (FFT)-based approaches, so I am now wondering whether the generation of correlated noise with the FFT and Perlin noise is more-or-less equivalent.

Generating correlated noise using FFT's

Below is a relatively simple code snippet that generates a correlated noise signal. It was inspired by work publised by Xiao and Voelz in 2006 for modeling partially coherent laser beam propagation through turbulence.

In [4]:
M  = 2**10  # Size of the 1D grid
L  = 10      # Physical size of the grid
dx = L / M  # Sampling period
fs = 1 / dx # Sampling frequency
df = 1 / L  # Spacing between frequency components
x  = np.linspace(-L/2, L/2,   num = M, endpoint = False)
f  = np.linspace(-fs/2, fs/2, num = M, endpoint = False)

sigma_f = 0.1
sigma_r = 1

# Define a Gaussian function and an uncorrelated random function in frequency space
F = fftshift(np.exp(-np.pi**2 * sigma_f**2 * f**2))
R = np.random.randn(f.size) + 1j * np.random.randn(f.size)

# Convert the product of the two functions back to real space
noise = 2 * np.pi * ifftshift(ifft(F * R)) * sigma_r / dx / np.sqrt(df)

# Plot the real and imaginary parts of the noise signal
plt.plot(x, np.real(noise), linewidth = 2, label = 'real')
#plt.plot(x, np.imag(noise), linewidth = 2, label = 'imaginary')
plt.xlabel('Spatial position')
plt.legend(loc = 'best')

This algorithm actually produces two independent and random signals that are correlated over distances comparable to sigma_f; one signal comes from the real part and the other from the imaginary part. I only plotted the real part because we only need one. At distances larger than sigma_f the signal is uncorrelated. The strength of the fluctuations is determined by sigma_r, which is the square root of the variance of the underlying, uncorrelated random signal.

To better understand how this algorithm it works, I think it's illustrative to plot the product of F and R from the code above before they are inverse Fourier transformed.

In [5]:
plt.plot(f, np.real(ifftshift(R)),   label = 'real part of R')
plt.plot(f, np.real(ifftshift(F*R)), label = r'real part of F $\times$ R')
plt.plot(f, ifftshift(F),            linewidth = 4, label = 'F')
plt.xlim((-20, 20))

As you can see, the frequency space multiplication is essentially modulating a zero-mean, Gaussian random signal with a Gaussian envelope centered at the origin. We then compute the inverse Fourier transform of this product, scale it by a few factors to ensure that energy is conserved, and take the real or imaginary parts to get two different, random, yet smooth signals. If you're familiar with Fourier transforms, you will probably recognize that the frequency space multiplication of F and R is equivalent to the real space convolution of their individual Fourier transform pairs. Thus, we've essentially smoothed a random signal with a Gaussian blurring kernel.

Though we could easily have done this with a convolution and avoided the need to think in frequency space, the visualization above is illustrative for this reason: we have essentially modified the amplitudes of the various spectral components of the uncorrelated random signal, producing correlated noise in the process.

What happens if we shape the random signal's spectrum differently, for example by using a top hat function instead of a Gaussian? Below is a plot the real part of the random signal generated from the top hat filtering compared to the previous one.

In [6]:
# Define a Gaussian function and an uncorrelated random function in frequency space
F_tophat = np.ones(f.size)

# Set values outside of the width of the tophat to zero.
F_tophat[np.logical_or(f < -0.5 / sigma_f, f > 0.5 / sigma_f)] = 0
F_tophat = fftshift(F_tophat)

# Convert the product of the two functions back to real space
noise_tophat = 2 * np.pi * ifftshift(ifft(F_tophat * R)) * sigma_r / dx / np.sqrt(df)

# Plot the real and imaginary parts of the noise signal
plt.plot(x, np.real(noise_tophat), linewidth = 1, label = 'Tophat, real')
plt.plot(x, np.real(noise),        linewidth = 2, label = 'Gaussian, real')
plt.xlabel('Spatial position')

We still have smooth, large-scale fluctuations that are about the same size as the inverse width of the top hat, but now there are additional high-frequency fluctuations which lead to a more "jagged" appearance in the blue, top hat curve above.

For completeness, here's how the real part of the spectrum is reshaped by multiplying by a top hat.

In [7]:
plt.plot(f, np.real(ifftshift(R)),          label = 'real part of R')
plt.plot(f, np.real(ifftshift(F_tophat*R)), label = r'real part of F $\times$ R')
plt.plot(f, ifftshift(F_tophat),            linewidth = 4, label = 'F')
plt.xlim((-40, 40))

Reasoning about the differences

So far we've learned that we could create correlated noise by filtering a zero-mean, Gaussian random signal in frequency space. We also learned that the nature of the resulting signals changes when the shape of the frequency space filter is modified. In the case of the top hat filter, the resulting random signal became a bit more jagged. This makes sense because the amplitudes of the mid-range spectral components were larger for the top hat filtered signal; there was no smooth Gaussian rolloff but instead a hard cutoff at +/- 0.5 / sigma_f. Even though the Gaussian-filtered signal had a larger bandwidth, the high frequency amplitudes were diminished exponentially, so they essentially didn't matter.

The connection with Perlin noise

Conceptually, Perlin noise is generated by adding together curves with different amplitudes and frequencies as explained here. This is similar to what we did in the examples above: we took sine waves with random phases and different frequencies, modulated their amplitudes, and then added them together to get a realization of a smooth random signal.

So what are the differences between the two approaches then? Here's what I can think of for now.

  1. FFT's use sine waves for their basis functions; Perlin noise uses curves that have been interpolated from random numbers.
  2. The frequencies of the sine waves in the FFT approach have a fixed spacing; the frequencies in Perlin noise are on an octave scale (each successive frequency is twice as much as the previous one).

This last point probably means that we can't apply the FFT approach to generate Perlin noise without some modifications because the FFT expects frequency samples at fixed, linear intervals. Once I find a Python implementation of Perlin noise, I'll try to better understand how it's generated and see whether it's fundamentally different than the FFT approach above.

One last example: 1/f noise

I realized the day after writing the draft of this post that the above discussion is useful for illustrating some of the "colored" noises. In this case, we use the ideas presented above to generate a realization of 1/f noise, also known as pink noise or fractal noise.

Can you guess how the amplitude of the spectrum of 1/f noise varies?

Surprise! It actually varies as $ \frac{1}{f^{1/2}} $. This is because the power spectrum of pink noise varies as 1/f; the amplitude of the spectral components is the square root of the average power delivered by each component. The proper Wikipedia sections on energy and power spectral densities explain this. The difficulty in actually implementing this filter, though, is what happens as we approach the DC value $ f = 0 $ in frequency space. At this value the filter becomes infinite.

Instead, we can easily approximate a realization of pink noise with a filter that is one near the origin and quickly rolls off into a $ f^{-1/2} $ behavior. A function that does this is

$$H \left( f \right) = \frac{1}{\sqrt{1 + af}}$$

The value $ a $ determines how quickly the function moves from a weak dependence on frequency to a $ f^{-1/2} $ dependence. For us, we should set it to a value that is relatively large compared to the frequency spacing in Fourier space. This type of filter is just a special case of the Chebyshev filter that is well-known in signal processing circles.

In [8]:
# Define the pink noise frequency-space filter for our random uncorrelated signals
def pinkNoiseFilter(freq):
    freq = np.abs(freq)
    return 1 / np.sqrt(1 + 10 * freq)
pinkNoiseFilter = np.vectorize(pinkNoiseFilter)
In [9]:
F_pinkNoise = fftshift(pinkNoiseFilter(f))

# Convert the product of the two functions back to real space
noise_pink = 2 * np.pi * ifftshift(ifft(F_pinkNoise * R)) * sigma_r / dx / np.sqrt(df)

# Plot the real and imaginary parts of the noise signal
plt.plot(x, np.real(noise_pink), linewidth = 1, label = '1/f, real')
plt.plot(x, np.real(noise),      linewidth = 2, label = 'Gaussian, real')
plt.xlabel('Spatial position')

Now we have the same general trends in the 1/f noise signal as in the realization of the Gaussian signal, but with a lot of high frequency noise added.

We can see that the high frequency noise is coming from poor suppression of high frequencies in the filtered spectrum. Remember that this is the spectrum of just one realization of the noise, this time without the real part of the random signal for clarity.

In [10]:
plt.plot(f, np.real(ifftshift(F_pinkNoise*R)), label = r'real part of F $\times$ R')
plt.plot(f, ifftshift(F_pinkNoise),            linewidth = 4, label = 'F')
plt.xlim((-40, 40))

Pretty cool! Though I still haven't made a good connection with Perlin noise, I now have a better understanding of how to generate many different kinds of realizations of noise signals on a computer. The general recipe is to

  1. generate a completely uncorrelated, zero-mean random signal in frequency space,
  2. multiply this signal by a filter that shapes the amplitudes of the signal,
  3. and perform an inverse FFT to bring the signal back into real space. We get two random signals for the price of one; one from the real part and one from the imaginary part.