Reproduce the color portrait and domain movie

I keep both panels under conserved Cahn–Hilliard dynamics and change the energy. r is the reference portrait and u the evolving RGB fields. The portrait energy is quadratic; the domain energy contains a double well and a gradient coefficient kappa=0.4. dt=0.1 and stabilization S=2. Channels are independent and all quantities are dimensionless. This reproduces the full movie source; the shorter introductory notebook remains available separately.

I develop this calculation in small steps. I run a cell, inspect its result and test it before adding the next part. Run the cells in order in Jupyter. A code cell defines a part of the calculation or carries it out; defining a function does not run it.

I permit free noncommercial teaching, learning and demonstrations with acknowledgment. Research and commercial use require my written permission. Earlier licensed releases retain their permissions. Read the code-use terms.

Copyright (c) 2026 Saswata Bhattacharyya.

Step 1: Prepare the Python environment

I install the packages listed in the next cell once in the same environment as the notebook. In Jupyter or Colab, remove the comment marker from the installation line below if a package is missing. NumPy stores arrays and evaluates numerical operations. Matplotlib displays plots. Pillow reads and writes images. This installation does not change the model parameters.

# %pip install numpy "pillow>=10.1" matplotlib
import numpy as np
assert np.isfinite(np.array([0.0, 1.0])).all()
print("NumPy is ready:", np.__version__)

NumPy is ready: 2.1.3

Step 2: Read the image and check its dimensions

I download my public portrait once, or replace this path with my own image. An RGB array has two position axes and one channel axis. The three channels are red, green and blue.

from pathlib import Path
from urllib.request import urlretrieve

work = Path("lesson-data")
work.mkdir(exist_ok=True)

def get_file(relative_path):
    # I download a public example only when its local copy is missing.
    destination = work / Path(relative_path).name
    if not destination.exists():
        urlretrieve("https://saswataiith.github.io/" + relative_path, destination)
    return destination

portrait = get_file("assets/images/saswata-landscape.png")
output = "portrait-abstraction-results"

Step 3: Prepare the next part of the calculation

I define the following operations. These are the same operations as in the original source, written on separate lines so their order is clear. Inspect their values before running the next cell.

from pathlib import Path
import sys
import numpy as np
from PIL import Image, ImageDraw, ImageFont

Step 4: A portable display font

I try the original font, then DejaVu Sans, then Pillow's default font. The font changes labels on the image only; it does not change a numerical field.

def portrait_font(size):
    fonts = [Path('/') / 'System' / 'Library' / 'Fonts' / 'Helvetica.ttc', 'DejaVuSans.ttf']
    for font in fonts:
        try:
            return ImageFont.truetype(str(font), size)
        except OSError:
            pass
    return ImageFont.load_default(size=size)

Step 5: Prepare inputs and output locations

I name the input data and output directory explicitly. Re-running these cells replaces files in that lesson's output directory. I keep the original image or composition field separate from the evolving array.

out = Path(output)
out.mkdir(parents=True, exist_ok=True)
r = np.asarray(Image.open(portrait).convert('RGB').resize((256, 312)), float) / 127.5 - 1
noise = np.random.default_rng(19).normal(0, 0.12, r.shape)
noise -= noise.mean((0, 1), keepdims=True)
u0 = r + noise
u = u0.copy()
mean0 = u.mean((0, 1))
axes = (0, 1)

Step 6: Set up the Fourier grid

I construct Fourier wave numbers from the grid spacing and domain size. A Laplacian multiplies each Fourier mode by minus its squared wave number. The zero Fourier mode represents the spatial mean.

kx = 2 * np.pi * np.fft.fftfreq(256)
ky = 2 * np.pi * np.fft.fftfreq(312)
k2 = (ky[:, None] ** 2 + kx[None, :] ** 2)[..., None]
fft = lambda x: np.fft.fft2(x, axes=axes)
ifft = lambda x: np.fft.ifft2(x, axes=axes).real
e0 = fft(noise)
N = 256 * 312
dt = 0.1
kap = 0.4

Step 7: Prepare the next part of the calculation

I define S, font, rows. These are the same operations as in the original source, written on separate lines so their order is clear. Inspect their values before running the next cell.

S = 2
font = portrait_font(23)
rows = []

Step 8: Define the energy diagnostic

I calculate the bulk and gradient contributions in the same convention as the source program. Fourier normalization is part of this definition. I compare this quantity through time within one model; energies from different functionals need not have the same meaning.

def energy(a):
    return np.mean((a * a - 1) ** 2 / 4) + kap / 2 * np.sum(k2 * abs(fft(a)) ** 2) / (3 * N * N)

Step 9: Advance the field and record the result

I apply the update in its original order. I record the mean, energy or other diagnostics at the specified intervals. I do not clip the numerical field. A saved picture is a display of the field, not a substitute for a conservation check.

for step in range(2001):
    if step % 25 == 0:
        t = step * dt
        left = r + ifft(e0 * np.exp(-k2 * t))
        energies = [0.5 * np.mean((left - r) ** 2), energy(u)]
        drift = np.max(abs(u.mean((0, 1)) - mean0))
        rows.append([t, *energies, drift, *u.mean((0, 1))])
        im = Image.new('RGB', (1120, 845), '#f5f7f2')
        d = ImageDraw.Draw(im)
        d.text((24, 16), 'CAHN-HILLIARD IN BOTH PANELS: THE ENERGY CHANGES', font=font, fill='#173d40')
        for x, a, title, e in [(24, left, 'Portrait energy: restores the face', energies[0]), (584, u, 'Double-well energy: makes domains', energies[1])]:
            d.text((x, 63), title, font=font, fill='#173d40')
            im.paste(Image.fromarray(np.uint8(np.clip((a + 1) * 127.5, 0, 255))).resize((512, 624)), (x, 108))
            d.text((x, 750), f'Mean energy density (dimensionless): {e:.6f}', font=font, fill='#173d40')
        d.text((24, 790), f'Time {t:.1f} | each channel mean conserved', font=font, fill='#173d40')
        im.save(out / f'frame-{step // 25:03d}.png')
    if step == 2000:
        break
    u = ifft((fft(u) - dt * k2 * fft(u ** 3 - u - S * u)) / (1 + dt * k2 * (S + kap * k2)))

Step 10: My check of this result

I check finite values, conservation or the known reference result as specified below. An assertion stops the notebook if the condition fails. Inspect the quantities and tolerance rather than suppressing a failure.

assert np.isfinite(u).all() and drift < 1e-10
v = np.array(rows)
assert np.max(np.diff(v[:, 2])) < 1e-09
np.savetxt(out / 'diagnostics.csv', rows, delimiter=',', header='time,portrait_energy,double_well_energy,max_channel_mean_drift,mean_R,mean_G,mean_B', comments='')
print('Final diagnostics:', rows[-1])

Final diagnostics: [200.0, np.float64(1.1488722250515884e-06), np.float64(0.054203544067834036), np.float64(2.736699755701011e-14), np.float64(0.25635801596276664), np.float64(0.23499795751631475), np.float64(0.09989267141149472)]

Step 11: Review and extend the calculation

I first restart the kernel and run all the required cells. Then I change one input, predict its effect and compare the result with that prediction. A successful check establishes the tested property; it does not establish convergence or equilibrium. I retain the complete script for command-line use. The notebook embeds the calculation directly and does not need that script to run.

Read the notebook execution checks → · More Python, Julia and C lessons →