Clearly, these clouds have a very complicated, seemingly chaotic structure. They are made up of many small parts and have a fractal-like appearance: when you zoom in on a portion of the cloud near the edge, the zoomed in version of the image would look very similar to the original image.
Despite the chaos, we could still try to get some statistical information about the clouds in this image. We could for example try to determine what the typical size of a smaller cloud within the larger cloud is. To do this, we could e.g. decompose the image into different size scales. We can use some image manipulation tool to rescale the image to a very coarse version, like a 2x2 pixel image. This image will have pixels that contain the dominant colour within each large pixel in the original image. If we subtract this dominant colour from the original image, we are left with a version of the image without these largest structures. We can then repeat this procedure for increasingly more refined versions of the original image (always subtracting out the dominant components in the coarser levels), until we end up at the original image resolution. Every image in between will now contain a measurement of how strong features in the original image are at that specific resolution, and the size scale associated with that resolution tells us what the size scale of those features is.
A power spectrum is in essence a measure of the strength of the different features at the different resolutions. For each resolution, it is defined as the variance of the features on that resolution, i.e. how much total spread there is in values. A large spread means a strong signal and hence a lot of variation on this scale.
In order to do the proposed decomposition in resolutions, we require a mathematical technique called a Fourier transform. The signal on each resolution is measured by setting up a periodic function consisting of sines and cosines, with the period of these sines and cosines, called , linked to the resolution (a period corresponds to a sine wave that fits exactly once inside the frame of the image). For each value of , we get two Fourier amplitudes: the amplitudes of the corresponding sine and cosine waves that best fit the pixels in one image dimension. Since the image has two dimensions, itself has two components, so that we end up with four Fourier amplitudes per value of . The sine and cosine amplitudes are usually combined into a single complex number for convenience.
In principle, this Fourier decomposition can be done for an infinite number of values. In practice, there is however a clear limit: when the period of the sine waves is equal to the number of pixels in one dimension, then each oscillation of the sine wave exactly covers one pixel. If we where to increase the resolution even more, then multiple sine waves would cover a single pixel, and the resulting amplitude would no longer change by increasing the resolution even more. The maximum frequency that is useful in an image is hence , with the number of pixels in the image. In practice, this means we need exactly the same number of complex pixels in Fourier space to represent our image, as there are pixels in the original image.
For a multi-dimensional Fourier decomposition, there are many combinations of that represent the same resolution . However, they all encode the same information, so that we usually bin these together. The appropriate measure for the power for a specific value is than the average power within the corresponding bin, multiplied with the volume of that bin in space. For a two dimensional Fourier decomposition, this bin volume increases linearly with , in three dimensions it even increases as .
Note that a Fourier decomposition is not the only possible way of decomposing a data set. Other spatial decompositions are possible, and the power spectrum for those can also be called the power spectrum. When the data set is taken on the surface of the sphere (e.g. the cosmic microwave background data that is measured on the celestial sphere), then a decomposition in terms of spherical harmonics is used. Again, it is possible to calculate a power in spherical harmonics space that links to a specific angular scale, and the corresponding curve of power versus angular scale is also called the power spectrum. Hence the confusion I mentioned before.
In order to calculate the power spectrum for a data set, we have to do the following:
Below, I will illustrate these steps for the cloud image introduced above.
This first step depends a lot on the data set under study. For the cloud image, we want to read in the pixel data for the image. We can do this using
import matplotlib.image as mpimg image = mpimg.imread("clouds.png")
This will simply create a 1000x1000
numpy.array of floating point values that represents the grey scale level of each pixel. You can visualise this list using
import pylab as pl pl.imshow(image) pl.show()
Since this function will map the pixel values to the default colour scale, the result will look quite interesting.
To take the Fourier transform of our two dimensional image data array, we will use
numpy.fft. This library has a number of functions to handle one, two and multi dimensional data arrays. We will use the multi dimensional function here, as that makes it possible to generalise the technique to three dimensional data sets quite easily:
import numpy as np fourier_image = np.fft.fftn(image)
The Fourier image array now contains the complex valued amplitudes of all the Fourier components. We are only interested in the size of these amplitudes. We will further assume that the average amplitude is zero, so that we only require the square of the amplitudes to compute the variances. We can get these values like this
fourier_amplitudes = np.abs(fourier_image)**2
To bin the results found above in space, we need to know what the layout of the return value of
numpy.fft.fftn is: what is the wave vector corresponding to an element with indices
j in the return array? If we were to make things complicated for ourselves, we could try to figure this out from the online documentation.
Alternatively, we could not worry about this, and just use the utility function provided by
kfreq = np.fft.fftfreq(1000) * 1000
This will automatically return a one dimensional array containing the wave vectors for the
numpy.fft.fftn call, in the correct order. By default, the wave vectors are given as a fraction of 1, by multiplying with the total number of pixels, we convert them to a pixel frequency. To convert this to a two dimensional array matching the layout of the two dimensional Fourier image, we can use
kfreq2D = np.meshgrid(kfreq, kfreq)
Finally, we are not really interested in the actual wave vectors, but rather in their norm:
knrm = np.sqrt(kfreq2D**2 + kfreq2D**2)
For what follows, we no longer need the wave vector norms or Fourier image to be laid out as a two dimensional array, so we will flatten them:
knrm = knrm.flatten() fourier_amplitudes = fourier_amplitudes.flatten()
To bin the amplitudes in space, we need to set up wave number bins. We will create integer value bins, as is common:
kbins = np.arange(0.5, 501., 1.)
Note that the maximum wave number will equal half the pixel size of the image. This is because half of the Fourier frequencies can be mapped back to negative wave numbers that have the same norm as their positive counterpart. The
kbin array will contain the start and end points of all bins; the corresponding values are the midpoints of these bins:
kvals = 0.5 * (kbins[1:] + kbins[:-1])
To compute the average Fourier amplitude (squared) in each bin, we can use
import scipy.stats as stats Abins, _, _ = stats.binned_statistic(knrm, fourier_amplitudes, statistic = "mean", bins = kbins)
Remember that we want the total variance within each bin. Right now, we only have the average power. To get the total power, we need to multiply with the volume in each bin:
Abins *= 4. * np.pi / 3. * (kbins[1:]**3 - kbins[:-1]**3)
Finally, we can plot the resulting power spectrum as a function of wave number, (typically plotted on a double logarithmic scale):
This power spectrum has some interesting features. First of all, we can see that most of the power is located at large scales (small wave numbers). This makes sense, as the image is dominated by the large cloud structures. Towards smaller scales (larger wave numbers), the power drops of almost linearly. This is an expression of the fractal nature of the cloud: at lower resolutions the same types of patterns return, but at increasingly lower signal. There is a huge spike in the power spectrum at a wave number of , or a size scale of about 3 pixels. This is the scale at which the smallest cloud features at the edges of the clouds manifest.
Below is the full script to plot the power spectrum for the cloud image. I hope the part that calculates the power spectrum might be useful for other applications. An extension to higher dimensions is straightforward; simply replace the
knrm calculation with a higher order equivalent.
import matplotlib.image as mpimg import numpy as np image = mpimg.imread("clouds.png") fourier_image = np.fft.fftn(image) fourier_amplitudes = np.abs(fourier_image)**2 kfreq = np.fft.fftfreq(1000) * 1000 kfreq2D = np.meshgrid(kfreq, kfreq) knrm = np.sqrt(kfreq2D**2 + kfreq2D**2) knrm = knrm.flatten() fourier_amplitudes = fourier_amplitudes.flatten() kbins = np.arange(0.5, 501., 1.) kvals = 0.5 * (kbins[1:] + kbins[:-1]) Abins, _, _ = stats.binned_statistic(knrm, fourier_amplitudes, statistic = "mean", bins = kbins) Abins *= 4. * np.pi / 3. * (kbins[1:]**3 - kbins[:-1]**3) pl.loglog(kvals, Abins) pl.xlabel("$k$") pl.ylabel("$P(k)$") pl.tight_layout() pl.savefig("cloud_power_spectrum.png", dpi = 300, bbox_inches = "tight")