Generation of Simulated Images

Build patterns, understand their rules, and explore how changing those rules changes the image. I start with a single circle, then move through grids, stripes, checkerboards, criss-cross bands, random particles, and periodic fields. Additional examples vary checkerboard aspect ratio and orientation, curve the stripes, and change their white-to-black width ratio.

These controlled synthetic images help us learn how shape, spacing, orientation, arrangement, and phase fraction are represented by pixels. They can be useful test inputs for image analysis and microstructure statistics. No physical length scale or material evolution model is assigned here.

See the patterns first

Click a labelled example to jump to its geometric rule, Python code, generated output, and a change to try.

1. Single circle 2. Horizontal row of circles 3. Square grid of circles 4. Horizontal and vertical stripes 5. Checkerboard 6. Diagonal stripes 7. Criss-cross diagonal bands 8. Staggered rows of circles 9. Randomly placed circles 10. Smooth periodic field and a two-phase pattern 11. Rectangular checkerboard 12. Rotated checkerboard 13. Wavy stripes 14. Stripes with unequal widths

Start with NumPy and Matplotlib

Use a Jupyter notebook or Google Colab. If needed, install NumPy and Matplotlib first. Run the setup once, then run any example below.

import numpy as np
import matplotlib.pyplot as plt

N = 400
# First axis = row (y); second axis = column (x).
coordinates = np.arange(N)
x, y = np.meshgrid(coordinates, coordinates)

I use meshgrid throughout this tutorial. It gives every pixel an x-coordinate and a y-coordinate in two full 400 × 400 arrays. We apply a geometric rule directly to these arrays. In each binary image, black is 0 (False) and white is 1 (True). Lengths are in pixels.

The examples use the usual image convention: the first array index is the row, and rows increase downward. All outputs below were generated by the displayed numerical rules.

A small example before the full image

np.arange lists equally spaced integers; np.meshgrid makes full coordinate arrays from that list. Let us try just three coordinates so that we can see every value:

coordinates = np.arange(3)
print(coordinates)
# [0 1 2]

x, y = np.meshgrid(coordinates, coordinates)
print(x)
# [[0 1 2]
#  [0 1 2]
#  [0 1 2]]
print(y)
# [[0 0 0]
#  [1 1 1]
#  [2 2 2]]

At row 1, column 2, x is 2 and y is 1. The same idea extends to our full image. The lowercase x and y used below have the same roles as uppercase X and Y in the Fourier examples. When a new function appears, I will explain it with a small example.

1. Single circle

Select pixels whose distance from the center is no greater than the radius.

Target pattern

Python code

radius = 50
pattern = (x - N//2)**2 + (y - N//2)**2 <= radius**2

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. The inequality creates a Boolean mask: True inside the circle and False outside.

Create an empty image

np.zeros creates an array filled with zeros. dtype=bool stores False initially; assigning True marks a white region.

small = np.zeros((2, 3), dtype=bool)
print(small)
# [[False False False]
#  [False False False]]

Try it. Change radius from 50 to 80. How does the white area change? Predict the result before running the code.

2. Horizontal row of circles

Place equal circles along the middle row, with a gap between neighboring circles.

Target pattern

Python code

radius, spacing = 12, 40
pattern = np.zeros((N, N), dtype=bool)
for cx in range(spacing//2, N, spacing):
    circle = (x - cx)**2 + (y - N//2)**2 <= radius**2
    pattern |= circle

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. The loop moves the center horizontally. |= adds each circle by taking the logical union of masks.

Try it. Set radius to 20: neighboring circles now touch because the diameter equals the spacing. Predict the result before running the code.

3. Square grid of circles

Repeat circles at equal intervals along both image axes.

Target pattern

Python code

radius, spacing = 12, 50
pattern = np.zeros((N, N), dtype=bool)
for cy in range(spacing//2, N, spacing):
    for cx in range(spacing//2, N, spacing):
        pattern |= (x - cx)**2 + (y - cy)**2 <= radius**2

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. Nested loops visit every combination of row and column centers.

Try it. Change spacing to 40 while keeping the radius fixed. Compare particle count and white area. Predict the result before running the code.

4. Horizontal and vertical stripes

Alternate white and black bands. Show both orientations explicitly.

Target pattern

Python code

width = 20
horizontal = (y//width) % 2 == 0
vertical = (x//width) % 2 == 0
pattern = horizontal

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. // gives the band index; % 2 alternates between even and odd bands. Both coordinate arrays already have shape (N, N), so each condition directly produces a complete image.

Turn the bands vertically

pattern = vertical

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Try it. Change width to 40. The full repeating period is twice the band width. Predict the result before running the code.

5. Checkerboard

Divide the image into square tiles. A tile is white when its row index plus column index is even.

Target pattern

Python code

tile = 40
pattern = ((x//tile + y//tile) % 2) == 0

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. Each neighboring tile changes the sum by one, so its color switches. Diagonal neighbors retain the same color.

Try it. Use different tile widths along x and y to produce rectangular checks. Predict the result before running the code.

6. Diagonal stripes

Use the sum of row and column coordinates to draw repeating diagonal bands.

Target pattern

Python code

width = 20
pattern = ((x + y)//width) % 2 == 0

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. Constant x+y defines a diagonal line. width is measured in the coordinate sum; the perpendicular band width is width/√2 pixels.

Try it. Replace x+y with x-y. What happens to the direction? Predict the result before running the code.

7. Criss-cross diagonal bands

Combine two families of diagonal bands to make a diamond network.

Target pattern

Python code

period, band = 80, 10
diagonal_1 = (x + y) % period < band
diagonal_2 = (x - y) % period < band
pattern = diagonal_1 | diagonal_2

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. | means logical OR: a pixel is white if it belongs to either family. At crossings the bands join.

Try it. Replace | with &: only the intersections remain. Compare the two patterns. Predict the result before running the code.

8. Staggered rows of circles

Offset every alternate row by half the horizontal spacing.

Target pattern

Python code

radius, spacing = 12, 50
row_spacing = 43
pattern = np.zeros((N, N), dtype=bool)
for row, cy in enumerate(range(25, N, row_spacing)):
    offset = (row % 2) * (spacing//2)
    for cx in range(25 + offset, N, spacing):
        pattern |= (x - cx)**2 + (y - cy)**2 <= radius**2

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. enumerate supplies a row counter as well as the center coordinate. Alternate offsets change the arrangement without changing circle size. The row spacing approximates a triangular arrangement.

Count rows while placing particles

enumerate is a Python function. It supplies a counter starting at zero together with each value; the counter tells us whether a row is even or odd.

for row, value in enumerate([25, 68, 111]):
    print(row, value)
# 0 25
# 1 68
# 2 111

Try it. Set offset to zero to recover aligned rows. Predict the result before running the code.

9. Randomly placed circles

Use a reproducible random seed to choose particle centers.

Target pattern

Python code

radius, count = 15, 45
rng = np.random.default_rng(7)
centers = rng.integers(radius, N-radius, size=(count, 2))
pattern = np.zeros((N, N), dtype=bool)
for cy, cx in centers:
    pattern |= (x - cx)**2 + (y - cy)**2 <= radius**2

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. Each pair contains a row and column center. Circles may overlap; their union is white, so the white area is not generally count × π × radius².

Generate reproducible random integers

default_rng creates a random-number generator. The seed 7 reproduces the sequence when a new generator is created. integers(0, 5, ...) chooses integers from 0 through 4; the upper limit 5 is excluded. size=(2, 2) asks for two pairs.

rng = np.random.default_rng(7)
centers = rng.integers(0, 5, size=(2, 2))
print(centers)
# [[4 3]
#  [3 4]]

Try it. Change the seed to 8, then increase count. Distinguish a new realization from a change in particle density. Predict the result before running the code.

10. Smooth periodic field and a two-phase pattern

Build a smooth field, then turn it into two regions by applying a threshold.

Target pattern

Python code

u = np.linspace(0, 1, N, endpoint=False)
X, Y = np.meshgrid(u, u)
field = np.sin(2*np.pi*4*X) + np.cos(2*np.pi*6*Y)
threshold = 0.3
pattern = field > threshold

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. There are four cycles along x and six along y. endpoint=False avoids duplicating the periodic endpoint. The threshold assigns each pixel to one of two regions.

Sample a smooth function

linspace makes a specified number of equally spaced coordinates. endpoint=False leaves out the final endpoint. sin and cos evaluate trigonometric functions using radians; np.pi is the number π. meshgrid works exactly as in the setup, but now uses fractional coordinates.

u = np.linspace(0, 1, 4, endpoint=False)
print(u)
# [0.   0.25 0.5  0.75]
print(np.sin(0), np.cos(0))
# 0.0 1.0

View the continuous field before thresholding

plt.imshow(field, cmap='viridis', origin='upper')
plt.colorbar(label='Field value')
plt.show()

View the same field as a surface

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
surface = ax.plot_surface(X, Y, field, cmap='viridis')
ax.set(xlabel='x', ylabel='y', zlabel='Field value')
fig.colorbar(surface, ax=ax, shrink=0.7)
plt.show()

Try it. Raise the threshold to 0.8. Predict whether the white fraction increases or decreases. Predict the result before running the code.

11. Rectangular checkerboard

Keep the alternating tile rule, but give tiles different widths and heights.

Target pattern

Python code

tile_width, tile_height = 50, 25
pattern = ((x//tile_width + y//tile_height) % 2) == 0

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. Separate coordinate scales stretch the square checkerboard into rectangular checks. Alternation is preserved even though the pattern now has different length scales along its two axes.

Try it. Swap 50 and 25. Compare the directions of the long tile edges. Predict the result before running the code.

12. Rotated checkerboard

Rotate the coordinate system before applying the checkerboard rule. This produces alternating diamond-shaped tiles.

Target pattern

Python code

angle = np.deg2rad(45)
xc, yc = x - N/2, y - N/2
u_rot = xc*np.cos(angle) + yc*np.sin(angle)
v_rot = -xc*np.sin(angle) + yc*np.cos(angle)
tile = 40
pattern = ((np.floor(u_rot/tile) + np.floor(v_rot/tile)) % 2) == 0

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. The two rotated coordinates remain perpendicular. floor gives the tile index on both positive and negative coordinates. This is a filled checkerboard; the criss-cross example instead contains narrow crossing bands.

Convert an angle and identify a tile

deg2rad converts degrees to radians for sin and cos. floor rounds downward to the nearest integer, including for negative coordinates. It gives the tile index after rotation.

print(np.deg2rad(90))
# 1.5707963267948966 (π/2 radians)
print(np.floor([-1.2, 0.3, 2.8]))
# [-2.  0.  2.]

Try it. Change the angle to 30 degrees. Rotation changes orientation while keeping the nominal tile size. Predict the result before running the code.

13. Wavy stripes

Shift horizontal bands by a sinusoidal amount that varies across the columns.

Target pattern

Python code

period, white_width = 40, 20
amplitude, wavelength = 12, 160
shift = amplitude*np.sin(2*np.pi*x/wavelength)
pattern = (y - shift) % period < white_width

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. Every column shifts the same stripe rule up or down. amplitude controls the displacement; wavelength controls the distance between successive crests. The stated band width is measured vertically, not normal to the curved boundary.

Try it. Set amplitude to zero to recover straight stripes. Then halve wavelength to make the bands bend more frequently. Predict the result before running the code.

14. Stripes with unequal widths

Use a repeating period with a white band narrower than the black band.

Target pattern

Python code

period, white_width = 50, 15
pattern = y % period < white_width

plt.figure(figsize=(5, 5))
plt.imshow(pattern, cmap='gray', vmin=0, vmax=1,
           origin='upper', interpolation='nearest')
plt.xlabel('Column (x)')
plt.ylabel('Row (y)')
plt.tight_layout()
plt.show()
print('White area fraction:', pattern.mean())

Generated output

How it works. Modulo locates each row within one period. White occupies the first 15 rows out of 50, giving a white area fraction of 0.30 here because the image contains eight complete periods.

Try it. Set white_width to 35 without changing period. The repetition spacing stays fixed while the white fraction changes. Predict the result before running the code.

Compare the patterns in Fourier space

I will now show how these familiar patterns appear in Fourier space. Let us begin with one image at a time. You can copy each complete example below into a notebook cell and run it. Once you understand the first example, you will notice that the next two use the same steps; only the line that creates the pattern changes.

Before we run the code, look at the stripes. In which direction would you move to cross from a white band into a black band? Keep your answer in mind—we will soon connect it to the positions of the Fourier spots.

Horizontal stripes: complete simple code

import numpy as np
import matplotlib.pyplot as plt

N = 400
coordinates = np.arange(N)
X, Y = np.meshgrid(coordinates, coordinates)

pattern = (Y // 20) % 2 == 0

# Display the original pattern.
plt.figure(figsize=(6, 6))
plt.imshow(pattern, cmap='gray')
plt.axis('off')
plt.show()

# Calculate the Fourier transform.
centered = pattern.astype(float) - pattern.mean()
F = np.fft.fft2(centered)
F = np.fft.fftshift(F)
amplitude = np.abs(F)

# Brighten weak values without changing the Fourier data.
relative = amplitude / amplitude.max()
bright = relative ** 0.3

# For this 400-pixel image, the center index is 200.
view = bright[160:241, 160:241]
plt.figure(figsize=(6, 6))
plt.imshow(view, cmap='gray', vmin=0, vmax=1,
           interpolation='nearest', origin='lower')
plt.axis('off')
plt.colorbar(label='Display brightness (not raw amplitude)')
plt.show()

Checkerboard: complete simple code

import numpy as np
import matplotlib.pyplot as plt

N = 400
coordinates = np.arange(N)
X, Y = np.meshgrid(coordinates, coordinates)

pattern = ((X // 40 + Y // 40) % 2) == 0

# Display the original pattern.
plt.figure(figsize=(6, 6))
plt.imshow(pattern, cmap='gray')
plt.axis('off')
plt.show()

# Calculate the Fourier transform.
centered = pattern.astype(float) - pattern.mean()
F = np.fft.fft2(centered)
F = np.fft.fftshift(F)
amplitude = np.abs(F)

# Brighten weak values without changing the Fourier data.
relative = amplitude / amplitude.max()
bright = relative ** 0.3

# For this 400-pixel image, the center index is 200.
view = bright[160:241, 160:241]
plt.figure(figsize=(6, 6))
plt.imshow(view, cmap='gray', vmin=0, vmax=1,
           interpolation='nearest', origin='lower')
plt.axis('off')
plt.colorbar(label='Display brightness (not raw amplitude)')
plt.show()

Criss-cross: complete simple code

import numpy as np
import matplotlib.pyplot as plt

N = 400
coordinates = np.arange(N)
X, Y = np.meshgrid(coordinates, coordinates)

pattern = ((X + Y) % 80 < 10) | ((X - Y) % 80 < 10)

# Display the original pattern.
plt.figure(figsize=(6, 6))
plt.imshow(pattern, cmap='gray')
plt.axis('off')
plt.show()

# Calculate the Fourier transform.
centered = pattern.astype(float) - pattern.mean()
F = np.fft.fft2(centered)
F = np.fft.fftshift(F)
amplitude = np.abs(F)

# Brighten weak values without changing the Fourier data.
relative = amplitude / amplitude.max()
bright = relative ** 0.3

# For this 400-pixel image, the center index is 200.
view = bright[160:241, 160:241]
plt.figure(figsize=(6, 6))
plt.imshow(view, cmap='gray', vmin=0, vmax=1,
           interpolation='nearest', origin='lower')
plt.axis('off')
plt.colorbar(label='Display brightness (not raw amplitude)')
plt.show()

Let us understand the code, one step at a time

np.arange(400) makes the coordinate list 0, 1, 2, …, 399. np.meshgrid turns it into two 400 × 400 coordinate arrays: X gives each pixel’s column and Y gives its row. For a tiny 3 × 3 example:

X = [[0, 1, 2],    Y = [[0, 0, 0],
     [0, 1, 2],         [1, 1, 1],
     [0, 1, 2]]         [2, 2, 2]]

// is floor division: it supplies the band or tile index. % is the remainder: remainder 0 or 1 alternates the colors. For criss-cross, | means logical OR: make a pixel white when either band condition is true. The comparisons produce True and False arrays.

pattern.astype(float) converts True and False to the numbers 1.0 and 0.0. pattern.mean() computes their average. Subtracting this average removes the uniform background and suppresses the central zero-frequency spot.

np.fft.fft2(centered) calculates the two-dimensional Fourier transform. Its output contains complex numbers encoding strength and phase. np.fft.fftshift(F) rearranges that array so zero frequency is at the center; it does not change the transform coefficients.

np.abs(F) takes the magnitude of each complex value. amplitude.max() finds the strongest magnitude. Dividing by it gives relative values between 0 and 1. This assumes a nonconstant pattern, as in our three examples. An entirely uniform image has a zero transform after mean subtraction, so this division would need a separate zero check.

Try the new array operations on tiny inputs

values = np.array([0, 1, 1, 0])
print(values.mean())  # 0.5: add the values and divide by 4
print(values.max())   # 1: the largest value

flags = np.array([True, False])
print(flags.astype(float))  # [1. 0.]

print(np.abs(3 + 4j))  # 5.0: √(3² + 4²)

small = np.ones((2, 2))  # a 2 × 2 array of ones
F_small = np.fft.fft2(small)
print(F_small)
# [[4.+0.j 0.+0.j]
#  [0.+0.j 0.+0.j]]
print(np.fft.fftshift(F_small))
# [[0.+0.j 0.+0.j]
#  [0.+0.j 4.+0.j]]

np.array converts a Python list into a NumPy array; np.ones creates an array filled with ones. This constant 2 × 2 image has only a zero-frequency contribution: its four values add to 4. fftshift moves that contribution from the first position to the central index. The notation j denotes the imaginary unit; abs converts a complex coefficient to its magnitude.

How can we make the weak spots brighter?

bright = relative ** 0.3 applies the same fractional power to every pixel. Values between 0 and 1 become brighter: 0 and 1 stay unchanged, but weak intermediate values increase.

Try power 0.2 to brighten weak spots further, 0.5 for a milder enhancement, or 1.0 to display the original relative magnitudes. Zero remains black. This creates a display array; amplitude and F retain the original data. A log or decibel scale is another useful option, but is unnecessary for this first exercise.

Let us zoom in to see the spots clearly

An isolated Fourier peak can occupy one pixel. Showing all 400 × 400 bins in a small figure makes that pixel difficult to see. The slice [160:241, 160:241] selects 81 rows and 81 columns around index 200. Python excludes the stop index 241, so the last included index is 240. The displayed offsets run from −40 to +40 cycles per image.

These are the actual bin values, displayed larger by zooming. interpolation='nearest' keeps their edges crisp. Try [180:221, 180:221] for a tighter ±20 view, or [100:301, 100:301] for a wider ±100 view. The crop omits frequencies outside its window; it does not remove them from F.

plt.figure(figsize=(6,6)) starts a six-inch square figure. imshow displays the array. cmap='gray' maps low values to black and high values to white. vmin=0, vmax=1 fixes those endpoints. origin='lower' puts the lowest displayed row at the bottom. axis('off') hides image indices, colorbar shows the brightness scale, and show displays the finished plot.

Each pattern is normalized independently: compare frequency locations and relative strengths, not absolute amplitudes across plots. Multiplying F by a constant will not help here because dividing by the maximum cancels that constant.

What are these Fourier spots telling us?

I find it helpful to think of the Fourier transform as a way of asking an image: “Which repeating patterns are present in you?” It checks waves that repeat at different rates and in different directions. A strong match gives a bright spot. Let us see how we can predict the positions of these spots using a little counting and simple mathematics.

1. What is spatial frequency?

Look at one white stripe followed by one black stripe. Together they form one complete repeating unit. The distance covered by that unit is called the period. In our example, each band is 20 pixels wide. So the complete period is 20 + 20 = 40 pixels. Notice that we must count both colors, not just one band.

How many such units fit into our 400-pixel-high image? Let us divide the image height by the period:

repetitions = image length / period
            = 400 / 40
            = 10 cycles per image

This is why the first stripe peaks occur 10 frequency bins from the center. If we double each band width to 40 pixels, the complete period becomes 80 pixels and the image contains only 400/80 = 5 cycles. The peaks move from ±10 to ±5, closer to the center. In physical units the same rule is frequency = 1/period.

2. Why do horizontal stripes give vertical Fourier peaks?

Did you notice that the Fourier spots form a vertical line even though the stripes are horizontal? Try tracing a horizontal path across one stripe. The color stays the same. Now trace a vertical path: you encounter white, black, white, black. The repetition occurs as you move along y. This is the direction that the Fourier spots indicate.

To put this observation into a simple equation, I will first use a smooth cosine wave with the same repetition rate:

g(x, y) = cos(2π × 10 × y / 400)
         = cos(2π × 10 × y / N)

Here x does not appear at all. Its horizontal frequency is zero; its vertical frequency is 10 cycles per image. Its Fourier dots therefore lie at (horizontal frequency, vertical frequency) = (0, +10) and (0, −10). The binary stripes contain this fundamental frequency plus higher harmonics.

Why a pair of dots? A real cosine contains two opposite complex-wave directions with equal weight. In mathematical notation, cos θ = (e^iθ + e^−iθ)/2. You need not use complex numbers to read the plot: for any real image, opposite Fourier positions have equal magnitudes. The two dots describe one real oscillation.

3. Why does the transform detect that repetition rate?

How does the transform recognize a repetition? Imagine placing a trial wave over our image. We multiply their values at each position and add the results. When the trial wave follows the actual pattern, the contributions add up strongly. When it repeats at a different rate, positive and negative contributions cancel over complete periods. Let us check this idea with our 10-cycle cosine wave. The transform also tests sine waves, so it can recognize a pattern even when that pattern is shifted.

For a pure wave g(y) = cos(2π × 10 × y / N):

test frequency 10:
sum g(y) × cos(2π × 10 × y / N)
  = sum cos²(2π × 10 × y / N)  → large positive result

test a different integer frequency, such as 7:
sum cos(2π × 10 × y / N) × cos(2π × 7 × y / N)
  = 0 over the N equally spaced samples

The cancellation follows the identity cos A cos B = [cos(A−B) + cos(A+B)]/2: nonzero integer-frequency waves sum to zero over the full sampled image. A finite image that does not contain whole periods can spread the signal into neighboring bins; this is called spectral leakage.

4. Why does a checkerboard give a grid of dots?

With 40-pixel tiles, each horizontal and vertical black–white sequence has period 80 pixels: five cycles across a 400-pixel image. A checkerboard combines alternation along x with alternation along y. Its centered pattern is proportional to the product of two square waves, one in each direction.

The sharp checkerboard may look complicated at first. I will use two smooth cosine waves to show the main idea. The following identity tells us what happens when we multiply a wave along x by a wave along y:

cos(2π × 5x/N) × cos(2π × 5y/N)
 = ½ cos[2π × (5x + 5y)/N]
 + ½ cos[2π × (5x − 5y)/N]

Thus the first Fourier positions are (±5, ±5), not simply dots on the two axes. Changing a cosine to a sine changes phases but retains these frequency positions. The sharp tile edges add odd harmonics in both directions, giving combinations such as (±15, ±5), (±5, ±15), and (±15, ±15). This produces the observed grid. A smooth single product would have only the four fundamental dots.

5. Why do criss-cross bands give diagonal families?

One band family uses x+y. Along a line x+y = constant, the value stays unchanged. Moving in the (1,1) direction changes x+y, so it crosses the bands. The other family uses x−y and varies in the (1,−1) direction.

First family:  cos[2π(x+y)/80]
Second family: cos[2π(x−y)/80]

Since N = 400, N/80 = 5:
First fundamental pair:  (+5,+5), (−5,−5)
Second fundamental pair: (+5,−5), (−5,+5)

These frequency directions are normal (perpendicular) to their bands. Sharp bands create further multiples along the diagonals. There is an additional detail in our code: logical OR combines the two binary families. For binary masks A and B, OR = A+B−AB. The product AB adds mixed-frequency combinations, so the complete criss-cross spectrum contains more dots than two diagonal lines alone.

6. Why do sharp boundaries create higher harmonics?

Now look at a black–white boundary. Its color changes abruptly, whereas a sine wave changes smoothly. One smooth wave cannot reproduce that sharp boundary. We need to add faster waves to build the edges. These extra repetition rates are called harmonics. For equal-width black and white bands, the familiar square-wave expansion is:

square wave ≈ mean + C × [sin θ + (1/3)sin 3θ
                              + (1/5)sin 5θ + …]

The constant C depends on the chosen amplitude convention. For our stripes, θ = 2π × 10y/N, so harmonics occur at ±10, ±30, ±50, and so on, subject to the sampled image's frequency limit. Higher harmonics are weaker: their amplitudes decrease approximately as 1, 1/3, 1/5, … before discretization effects. Unequal white and black widths generally introduce even harmonics too.

7. What is the central DC spot, and why remove it?

At zero frequency, the test wave is constant: it never oscillates. The transform therefore adds all the pixel values. For a binary image, this measures how much white area is present:

F(0,0) = sum of all image pixels
       = N² × mean pixel value

Subtracting the mean makes that sum zero. This removes the uniform background contribution and lets us focus on variations. It does not change the nonzero Fourier coefficients, apart from numerical roundoff.

8. How does our simple display brighten the spots?

You may still find some spots difficult to see. I first divide all magnitudes by the largest one, giving values between 0 and 1. Then I raise them to power 0.3. A weak value such as 0.01 becomes approximately 0.25, so it appears gray instead of nearly black. We have changed how the result is displayed; the Fourier transform itself remains unchanged. The central ±40 crop enlarges frequency bins; it includes stripe harmonics at ±10 and ±30 but excludes ±50 and above. See the operation-by-operation explanation and brightness table above for the exact code.

Read more in my Fourier transform post. For further exploration, try reconstructing colored patterns with PCA, or browse the ImageSimulation examples.

What these examples do—and what they leave open

Each pattern is a controlled test case, not a claim that all materials look this way. The main limitation is idealization: most examples repeat exactly, and random circles can overlap. Rasterized boundaries also depend on pixel resolution. Before using a generated image for quantitative work, check area fraction, feature dimensions, connectivity, overlaps, and sensitivity to resolution. Opposite image edges do not automatically match; periodic boundary conditions require an explicit check, especially for rotated or curved patterns.

Connect these patterns to microstructures

Circles can represent idealized particle sections; stripes can represent lamellar arrangements; criss-cross bands can represent connected networks. These are geometric analogues. Real microstructures additionally have distributions of sizes, shapes, crystallographic orientations, and interfaces. A 2D white area fraction should not automatically be called a 3D volume fraction without appropriate stereological assumptions.