Image Blending

Image blending joins two images, or a region of one and the rest of another, into a single picture in which the join cannot be seen. It is the last step of panorama stitching, the core of exposure fusion and seamless cloning, and, as deepfake detection reports, a step shared by most face-manipulation methods, which paste an altered face into a background frame. The difficulty is that the two sources rarely agree: they differ in exposure and colour, and they carry different textures. A sharp cut shows every disagreement as a seam. A gradual cross-fade hides the seam but lays one texture over the other across the whole transition. The methods differ in how they divide this compromise between scales, or, in the gradient domain, in where they put the disagreement instead.

This article covers four methods: the hard cut, feathering (a linear cross-fade), multiband blending with Laplacian pyramids (Burt and Adelson, 1983), and Poisson or gradient-domain blending (Pérez, Gangnet and Blake, 2003). An example on synthetic images measures the compromise each makes. The image pyramid article builds the pyramids used here; image matting derives the compositing equation that every weighted blend uses.

A note on sources. Burt and Adelson’s multiresolution spline paper, Pérez et al.’s Poisson Image Editing, Porter and Duff’s compositing paper and Mertens et al.’s exposure fusion are on hosts this site’s research does not use, so they were not read. What is said about them comes from the published image pyramid, image matting and HDR imaging articles and from two articles in the journal Image Processing On Line: Di Martino, Facciolo and Meinhardt-Llopis (2016) on Poisson editing, and Hessel (2018) on exposure fusion. Each statement says which.

Purpose

Blending answers one question: given two images $A$ and $B$ that are each correct on their own side of a boundary, what should the output be near the boundary? Three things can go wrong, and the example below measures each.

• A visible seam. If $A$ and $B$ differ in brightness by an offset $\Delta$, a hard cut shows a step of height $\Delta$ along the boundary.
• Texture overlap. Wherever the output is a mixture of $A$ and $B$, both textures are present at reduced contrast, one over the other. When the two images are of different content, the mixture looks like a double exposure; we call it texture overlap, and leave “ghosts” to the doubled content of misaligned images, the sense in which Hessel (2018, §1) uses the word.
• A shift away from the sources. A method can hide the seam by changing the images themselves, so that one side, or both, no longer has its own brightness far from the boundary.

When to Use It

• Panoramas and mosaics. Overlapping photographs, aligned by a homography and warped onto a common surface, differ in exposure and vignetting. The image pyramid article lists pyramid blending among the compositing methods of panorama stitching, and computational photography treats stitching as registration, warping and blending.
• Exposure fusion. Several exposures of one scene are blended with per-pixel weights that favour well-exposed, saturated, contrasted pixels (Mertens et al., 2007, as HDR imaging and Hessel, 2018, §2, describe it). Hessel reports the authors’ observation that blending with the weights directly “would lead to strong seams”, which is why they blend in a pyramid.
• Inserting an object or region. Seamless cloning pastes a region of one photograph into another so that its colours adapt to the new surroundings (Di Martino et al., 2016, §2).
• When not to blend. If the region has a genuine soft boundary, such as hair or smoke, the problem is matting: the opacity must be estimated, not designed. If the images are misaligned, blending does not fix the misalignment; registration comes first (image registration).

How It Works

Weighted blending: hard cut and feathering

Every weighted blend is the compositing equation with a designed weight map $w(x) \in [0, 1]$:

$$I(x) = w(x)\, A(x) + \bigl(1 - w(x)\bigr)\, B(x).$$

This is the compositing equation of image matting, with $w$ in the role of the opacity $\alpha$; that article reports Chuang et al.’s description of it as the “over” operation and their credit to Porter and Duff (1984) for the alpha channel. In matting $\alpha$ is unknown and estimated; in blending $w$ is chosen. A hard cut takes $w = 1$ on $A$’s side and $0$ on $B$’s. Feathering replaces the step by a ramp of width $w_f$ across a straight seam at $x_s$,

$$w(x) = \operatorname{clip}\!\left(\tfrac{1}{2} - \frac{x - x_s}{w_f},\; 0,\; 1\right),$$

and for a seam of arbitrary shape the binary mask can be smoothed, for instance with a Gaussian filter (our construction; the example uses the straight ramp).

Feathering faces a trade-off that no choice of $w_f$ escapes. A brightness offset $\Delta$ is spread over the ramp, so the steepest brightness step falls as $\Delta / w_f$: a wide ramp hides the seam. But across the whole ramp both textures are present, so the width of the overlap grows as $w_f$. The example makes this exact: the product of the two measures is the same for every even width.

Multiband blending

Burt and Adelson’s multiresolution spline (1983) splits the trade-off by scale. As the image pyramid article describes it, each Laplacian pyramid band of the two images is blended separately, with a transition zone proportional to the wavelengths in that band, and the result is reconstructed: coarse differences such as overall brightness blend over a wide zone, while fine detail meets along a narrow seam. For a region of arbitrary shape, the binary mask $M$ is turned into a Gaussian pyramid whose levels serve as the per-band weights. With $L_l\{\cdot\}$ the Laplacian band at level $l$ and $G_l\{\cdot\}$ the Gaussian level,

$$L_l\{I\} = G_l\{M\}\, L_l\{A\} + \bigl(1 - G_l\{M\}\bigr)\, L_l\{B\}, \qquad l = 0, \dots, N - 1,$$

and $I$ is obtained by collapsing the pyramid, adding each band to the expanded level above it, from the top down. The top level $L_{N-1}$ is the coarsest Gaussian image, which Hessel (2018, §3) calls the residual. Hessel writes the same operation for $K$ images with normalized weights $w_k$ as $L_l\{v\} = \sum_k G_l\{w_k\}\, L_l\{u_k\}$ (eq. 6), and credits it to Ogden et al. (1985).

The mask’s Gaussian pyramid does the work. At the finest level the weight is the binary mask itself, so the finest band switches from $A$ to $B$ at the boundary. Each level up, the mask is smoothed once more and sampled half as densely, so (in our reading) the weights at level $l$ change over a zone roughly $2^l$ pixels wide; once expanded to full resolution, the coarse bands, which hold the brightness offset, cross over slowly.

Depth. How wide the coarsest transition is depends on how many levels the pyramid has. Hessel (2018, §8) traces a low-frequency halo in exposure fusion to the residual level: a deeper pyramid removes it but increases out-of-range values, while a shallower one leaves a thinner halo, which is more visible; the halo’s width is of the order of the width of the residual’s equivalent blur kernel. In our reading, the same holds for an offset between two images: the residual decides how far the offset is spread.

Range. The blended bands need not reconstruct to values inside the input range. Hessel quotes Mertens et al.’s 2009 journal version on this: pyramid reconstruction does not guarantee that the result lies in $[0, 1]$ even when the inputs do, and shifting and scaling the result risks reducing contrast (§8).

Poisson (gradient-domain) blending

Pérez et al. (2003) were the first to formulate gradient-domain editing as a variational problem, according to Di Martino et al. (2016, §1.1). Instead of mixing intensities, the output is asked to have prescribed gradients inside the region $\Omega$ and to agree with the destination image $f^*$ on its boundary. In Di Martino et al.’s notation (eqs. 1 and 2),

$$\min_{f} \int_{\Omega} \lVert \nabla f - \mathbf{v} \rVert^2 \, dx \quad \text{with} \quad f\big|_{\partial\Omega} = f^*\big|_{\partial\Omega}, \qquad\text{so that}\qquad \Delta f = \operatorname{div} \mathbf{v} \;\text{ in } \Omega,$$

the Poisson equation that gives the method its name. For seamless cloning the guidance field $\mathbf{v}$ is the gradient of the source; Di Martino et al. also show the average of the source and destination gradients and the larger of the two at each pixel (eq. 19). On a pixel grid, the Laplacian becomes the 5-point stencil and the problem a sparse linear system: for every pixel $p$ of $\Omega$,

$$\sum_{q \in N(p)} (f_p - f_q) = \sum_{q \in N(p)} v_{pq},$$

with $v_{pq}$ the source’s difference $B_p - B_q$ in cloning, and with $f_q$ taken from the destination when $q$ lies outside $\Omega$. Where $\Omega$ meets the edge of the image, Di Martino et al. replace the Dirichlet condition by a Neumann condition, zero derivative across the edge (§3.1). They also give a Fourier solution, which solves the equation on the whole rectangle with the source gradient inside $\Omega$ and the destination gradient outside (§3.2, eqs. 12 and 13).

Gradient-domain blending does not mix the two images at all inside $\Omega$: it keeps the source’s detail and adds to it the smoothest correction, a harmonic function, that makes it meet the boundary values. The Laplacian operator article describes the result as content that keeps its detail but blends in without a visible seam. What the correction does to the region’s brightness is a different matter, and the example measures it.

How to Apply It

The choices that matter in practice are the width of the transition for feathering, the number of levels for multiband blending, and the region and boundary for Poisson blending. The example measures all three on images whose ground truth is known exactly.

Example: seam, texture overlap and brightness shift

Images. Two synthetic grey images $A$ and $B$, $128 \times 256$ pixels, in floating point. $A$ has brightness $a = 0.35$ and $B$ brightness $b = 0.65$, an offset $\Delta = 0.30$. $A$ owns columns 0–127 and $B$ columns 128–255, so the seam lies between columns 127 and 128. Each image carries its own texture, a sum of sinusoids along $y$ with a whole number of cycles over the 128 rows, whose amplitudes vary smoothly and randomly along $x$ (seed 0). Two consequences are built in: every column of a texture has mean zero, and $A$’s and $B$’s textures use disjoint sets of frequencies, so they are exactly orthogonal column by column. There are two texture conditions:

• fine: periods 3.3–5.1 pixels, amplitudes smoothed along $x$ with $\sigma = 1$;
• coarse: periods 8.5–25.6 pixels, amplitudes smoothed with $\sigma = 4$.

Methods. A hard cut; feathering with ramps 4, 16 and 64 pixels wide; multiband blending with $N = 3$, 5 and 8 pyramid levels (8 reaches a $1 \times 2$ top level), using the hard-cut mask; and Poisson cloning of $B$ into columns 128–255, with $B$’s differences as guidance, $A$’s values in column 127 as the boundary and zero derivative across the image’s edges. The pyramid is built as in the image pyramid article, with Burt and Adelson’s five-tap kernel ($a = 0.4$), except at the borders: ours extends the image by reflection about the centre of the edge pixel (SciPy’s mode="mirror", the whole-sample symmetric extension among the boundary rules image interpolation lists) in both REDUCE and EXPAND; with this extension EXPAND returns a constant image unchanged, and a check line prints that the Laplacian bands of a constant image are zero.

Measures. Let $m(x)$ be the mean of column $x$ of a result. Because the textures have zero column means, $m$ is the brightness the blend gives each column.

• Seam visibility $S$: the largest jump $|m(x+1) - m(x)|$ anywhere in the image, divided by $\Delta$. A hard cut has $S = 1$.
• Texture overlap $G$: in each column, the least-squares coefficients $\alpha_A(x)$ and $\alpha_B(x)$ of the two textures in the result (exact projections, since the textures are orthogonal to each other and to constants); $G = \sum_x \max\bigl(0, \min(\alpha_A, \alpha_B)\bigr)$, in columns. It counts how many columns’ worth of both textures appear at once.
• Far-field shift $F_\text{left}$, $F_\text{right}$: the largest deviation of $m(x)$ from the side’s own brightness ($a$ or $b$) over the 32 outermost columns of each side, divided by $\Delta$.

Checks. Before the comparison the script checks the data (zero column means, orthogonality, non-degenerate columns, values inside $[0, 1]$, and that the projection recovers $\alpha_A = 1$, $\alpha_B = 0$ on $A$ alone) and every component: the bands of a constant image are zero, the pyramid reconstructs its input, multiband blending with a constant mask returns $A$ or $B$, the Poisson system is solved to rounding error, and cloning $A + \Delta$ onto $A$ returns $A$ (the gradients are the same, so the boundary forces the original). Feathering is also checked against Pillow 12.0.0’s Image.composite on 8-bit copies. Image.composite pastes $A$ through the mask onto a copy of $B$, and for an 8-bit mask Pillow computes $\operatorname{round}\bigl((B \cdot (255 - m) + A \cdot m)/255\bigr)$ for a mask value $m$ (the BLEND macro of ImagingUtils.h, applied by paste_mask_L in Paste.c); the NumPy expression is the same rounding.

What to expect. Two results follow from the construction (our derivations).

• Feathering. With zero-mean texture columns and a weight that depends on $x$ only, $m(x) = w(x)\,a + (1 - w(x))\,b$, so each step of the brightness is $\Delta$ times the step of the weight. The hard cut has one step of $\Delta$: $S = 1$. A ramp of width $w_f$ has steps of $1/w_f$: $S = 1/w_f$. Orthogonality makes the coefficients exactly $\alpha_A = w(x)$ and $\alpha_B = 1 - w(x)$; for an even width the ramp takes the values $\tfrac12 \pm (j + \tfrac12)/w_f$, $j = 0, \dots, w_f/2 - 1$, so $G = 2\sum_j \bigl(\tfrac12 - (j + \tfrac12)/w_f\bigr) = w_f/4$. Hence $S \cdot G = 1/4$ for every even width: no ramp has both $S \le 0.1$ and $G \le 1$. Even the 64-pixel ramp covers only columns 96–159, clear of the 32 outermost columns of each side, so $F = 0$.
• Poisson. Write the solution on $B$’s side as $f = B + h$. The discrete equations then say that $h$ is harmonic, with the values $A - B$ in column 127 and zero derivative at the image’s edges. Summing the equations over a column, the vertical terms cancel (the top and bottom edges have zero derivative), so the column mean of $h$ satisfies the one-dimensional Laplace equation: it is linear in $x$, starts at the mean of $A - B$ in column 127, which is $a - b = -\Delta$, and has zero slope at the right edge. It is therefore $-\Delta$ in every column. Every column of $B$’s side has mean $b - \Delta = a$: there is no brightness step anywhere ($S = 0$), $A$’s side is untouched ($F_\text{left} = 0$), and $B$’s far side has moved by the whole offset ($F_\text{right} = 1$).

The other hypotheses were written from the sources before any run:

• H1: the feathering results above, in both textures.
• H2: for the fine texture, multiband blending with $N = 8$ has $S < 1/64$ (below the widest ramp’s) and $G < 1$ (below the narrowest ramp’s) at once, as transition zones proportional to wavelength, in the image pyramid article’s description of Burt and Adelson’s spline, imply.
• H3: for each $N$, multiband $G$ for the coarse texture is at least twice $G$ for the fine texture, since the coarse wavelengths are about four times longer.
• H4: following Hessel (2018, §8), in both textures multiband $S$ falls strictly as $N$ grows, and the far-field shift $\max(F_\text{left}, F_\text{right})$ is larger at $N = 8$ than at $N = 3$.
• H5a: the Poisson results above, to $10^{-9}$.
• H5b: Poisson $G$ is below 1 for the fine texture and larger for the coarse texture than for the fine (our reasoning: a mismatch at the boundary in a sinusoid of $k$ cycles over the height decays into the region roughly like $e^{-2\pi k d / H}$ at distance $d$, faster for finer texture).

The checks and their tolerances were fixed, with the values a correct implementation reaches on these images, before the pilot run (the added check described next, before the repeated pilot). One change was made before the first run, after the pilot run. The EXPAND step as first written blurred the zero-inserted image with SciPy’s mode="reflect", which repeats the edge sample; with zero insertion this puts a sample where a zero belongs, so EXPAND did not return a constant image near the border. The defect was found while the pilot’s results were examined: it changed the pilot’s multiband rows, and some multiband hypotheses failed in that pilot. None of the planned checks could see it, because reconstruction is exact for any EXPAND and a constant mask hides the bands. On the editor’s decision, as the repair of a broken component, the extension was changed to mode="mirror" in REDUCE and EXPAND and the check on the bands of a constant image was added; the repaired script was piloted once more, that pilot was fixed in advance as final whatever it printed, and its script, unchanged, is the first run shown here. The hypotheses and criteria were not changed. The script uses NumPy 2.3.5, SciPy 1.16.3 and Pillow 12.0.0 and runs in under a second; two runs gave identical output.

import os

os.environ["OMP_NUM_THREADS"] = "1"  # one BLAS thread, set before NumPy is imported
os.environ["OPENBLAS_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"

import numpy as np
from scipy import ndimage, sparse
from scipy.sparse.linalg import spsolve
from PIL import Image

H, W = 128, 256
SEAM = W // 2                  # A owns columns 0..127, B owns columns 128..255
a, b = 0.35, 0.65              # the two images' brightness levels
DELTA = b - a
rng = np.random.default_rng(0)
y, x = np.arange(H)[:, None], np.arange(W)


def texture(ks, sigma_x, rms=0.08):
    """Sum of whole-cycle sinusoids along y with amplitudes that vary smoothly along x.
    Every column has zero mean, and textures built from disjoint k sets are orthogonal column by column."""
    t = np.zeros((H, W))
    for k in ks:
        c = ndimage.gaussian_filter1d(rng.standard_normal(W), sigma_x, mode="reflect")
        s = ndimage.gaussian_filter1d(rng.standard_normal(W), sigma_x, mode="reflect")
        t += c * np.cos(2 * np.pi * k * y / H) + s * np.sin(2 * np.pi * k * y / H)
    return t * rms / np.sqrt(np.mean(t**2))


TEX = {  # (texture of A, texture of B)
    "fine": (texture([25, 29, 33, 37], 1.0), texture([27, 31, 35, 39], 1.0)),    # periods 3.3-5.1 px
    "coarse": (texture([5, 9, 13], 4.0), texture([7, 11, 15], 4.0)),            # periods 8.5-25.6 px
}

# ---- Laplacian pyramid with Burt and Adelson's 5-tap kernel (a = 0.4), borders reflected about the edge pixel ----
KERNEL = np.array([0.05, 0.25, 0.4, 0.25, 0.05])


def blur(img):
    # mode="mirror" (d c b | a b c d | c b a) keeps zero-inserted samples on the even grid, so EXPAND returns constants
    out = ndimage.convolve1d(img, KERNEL, axis=0, mode="mirror")
    return ndimage.convolve1d(out, KERNEL, axis=1, mode="mirror")


def reduce(img):
    return blur(img)[::2, ::2]


def expand(img, shape):
    up = np.zeros(shape)
    up[::2, ::2] = img
    return 4.0 * blur(up)


def gauss_pyr(img, n):
    p = [img]
    for _ in range(n - 1):
        p.append(reduce(p[-1]))
    return p


def lap_pyr(img, n):
    g = gauss_pyr(img, n)
    return [g[k] - expand(g[k + 1], g[k].shape) for k in range(n - 1)] + [g[-1]]


def collapse(lap):
    img = lap[-1]
    for band in reversed(lap[:-1]):
        img = band + expand(img, band.shape)
    return img


# ---- The blending methods ----
def feather(w):
    """Weight of A per pixel: a hard step (w = 0) or a linear ramp w pixels wide centred on the seam."""
    if w == 0:
        return np.broadcast_to((x < SEAM).astype(float), (H, W))
    return np.broadcast_to(np.clip(0.5 - (x - (SEAM - 0.5)) / w, 0, 1), (H, W))


def multiband(A, B, wA, n):
    """Blend each Laplacian band with the matching level of the mask's Gaussian pyramid."""
    la, lb, gm = lap_pyr(A, n), lap_pyr(B, n), gauss_pyr(wA, n)
    return collapse([g * p + (1 - g) * q for g, p, q in zip(gm, la, lb)])


def poisson_clone(A, B, seam):
    """Columns >= seam: solve the discrete Poisson equation with B's gradients as guidance,
    A's values on the column left of the seam (Dirichlet) and zero normal derivative on the image border."""
    h, w = A.shape
    nw = w - seam
    idx = np.arange(h * nw).reshape(h, nw)
    yy, xx = np.mgrid[0:h, seam:w]
    rows, cols, vals = [], [], []
    rhs, diag = np.zeros(h * nw), np.zeros(h * nw)
    for dy, dx in [(-1, 0), (1, 0), (0, -1), (0, 1)]:
        qy, qx = yy + dy, xx + dx
        ok = (qy >= 0) & (qy < h) & (qx < w)          # neighbours outside the image are dropped (Neumann)
        p = idx[ok]
        np.add.at(diag, p, 1.0)
        np.add.at(rhs, p, B[yy[ok], xx[ok]] - B[qy[ok], qx[ok]])   # guidance: B's differences
        inner = qx[ok] >= seam
        rows.append(p[inner])
        cols.append(idx[qy[ok][inner], qx[ok][inner] - seam])
        vals.append(-np.ones(inner.sum()))
        np.add.at(rhs, p[~inner], A[qy[ok][~inner], qx[ok][~inner]])  # known boundary values from A
    rows.append(np.arange(h * nw)); cols.append(np.arange(h * nw)); vals.append(diag)
    M = sparse.csr_matrix((np.concatenate(vals), (np.concatenate(rows), np.concatenate(cols))),
                          shape=(h * nw, h * nw))
    f = spsolve(M.tocsc(), rhs)
    out = A.copy()
    out[:, seam:] = f.reshape(h, nw)
    return out, np.abs(M @ f - rhs).max()


# ---- Measures ----
def measures(r, TA, TB):
    m = r.mean(axis=0)                                   # column-mean brightness
    S = np.abs(np.diff(m)).max() / DELTA                 # steepest brightness step, in units of the offset
    aA = (r * TA).sum(0) / (TA**2).sum(0)                # least-squares weight of A's texture per column
    aB = (r * TB).sum(0) / (TB**2).sum(0)
    G = np.clip(np.minimum(aA, aB), 0, None).sum()      # columns' worth of both textures at once
    FL = np.abs(m[:32] - a).max() / DELTA                # far-field brightness shift, left and right
    FR = np.abs(m[-32:] - b).max() / DELTA
    return S, G, FL, FR


# ---- Checks ----
print("data: max |column mean| %.1e, max |cos(T_A, T_B)| per column %.1e" % (
    max(np.abs(t.mean(0)).max() for pair in TEX.values() for t in pair),
    max(np.abs((TA * TB).sum(0) / np.sqrt((TA**2).sum(0) * (TB**2).sum(0))).max() for TA, TB in TEX.values())))
for name, (TA, TB) in TEX.items():
    A, B = a + TA, b + TB
    aA = (A * TA).sum(0) / (TA**2).sum(0)
    aB = (A * TB).sum(0) / (TB**2).sum(0)
    print(f"data {name}: min column RMS / RMS {min((np.sqrt((t**2).mean(0)) / np.sqrt((t**2).mean())).min() for t in (TA, TB)):.3f}, "
          f"values in [{min(A.min(), B.min()):.3f}, {max(A.max(), B.max()):.3f}], "
          f"fit on A: max|alpha_A - 1| {np.abs(aA - 1).max():.1e}, max|alpha_B| {np.abs(aB).max():.1e}")
    err_rec = max(np.abs(collapse(lap_pyr(A, n)) - A).max() for n in (3, 5, 8))
    err_ones = max(max(np.abs(multiband(A, B, np.ones((H, W)), n) - A).max(),
                       np.abs(multiband(A, B, np.zeros((H, W)), n) - B).max()) for n in (3, 5, 8))
    c9 = max(np.abs(band).max() for c in (a, b) for n in (3, 5, 8) for band in lap_pyr(np.full((H, W), c), n)[:-1])
    print(f"check {name}: Laplacian bands of a constant image (levels below the top, N = 3, 5, 8) {c9:.1e}")
    _, res = poisson_clone(A, B, SEAM)
    same, res2 = poisson_clone(A, A + DELTA, SEAM)
    print(f"check {name}: pyramid reconstruction {err_rec:.1e}, multiband with constant mask {err_ones:.1e}, "
          f"Poisson residual {max(res, res2):.1e}, Poisson clone of A + offset vs A {np.abs(same - A).max():.1e}")

TA, TB = TEX["fine"]
A8 = np.rint(np.clip(a + TA, 0, 1) * 255).astype(np.uint8)
B8 = np.rint(np.clip(b + TB, 0, 1) * 255).astype(np.uint8)
pil_diff = 0
for w in (0, 4, 16, 64):
    m8 = np.ascontiguousarray(np.rint(255 * feather(w)).astype(np.uint8))
    pil = np.asarray(Image.composite(Image.fromarray(A8), Image.fromarray(B8), Image.fromarray(m8))).astype(int)
    ref = np.rint((A8 * m8.astype(float) + B8 * (255.0 - m8)) / 255).astype(int)
    pil_diff = max(pil_diff, np.abs(pil - ref).max())
print(f"check: Pillow Image.composite vs NumPy feathering on 8-bit copies, max difference {pil_diff} codes")

# ---- The comparison ----
methods = [("hard cut", lambda A, B: feather(0) * A + (1 - feather(0)) * B)]
methods += [(f"feather w={w}", lambda A, B, w=w: feather(w) * A + (1 - feather(w)) * B) for w in (4, 16, 64)]
methods += [(f"multiband N={n}", lambda A, B, n=n: multiband(A, B, feather(0), n)) for n in (3, 5, 8)]
methods += [("Poisson", lambda A, B: poisson_clone(A, B, SEAM)[0])]

results = {}
for name, (TA, TB) in TEX.items():
    A, B = a + TA, b + TB
    print(f"\n{name} texture        S/Delta      G   F_left  F_right")
    for label, blend in methods:
        S, G, FL, FR = measures(blend(A, B), TA, TB)
        results[name, label] = (S, G, FL, FR)
        print(f"  {label:16s} {S:9.3e} {G:7.3f} {FL:8.4f} {FR:8.4f}")

def expected(w):  # derived: hard cut S = 1, G = 0; a ramp of width w has S = 1/w, G = w/4
    return (1.0, 0.0) if w == 0 else (1 / w, w / 4)


cut = [(w, results[name, "hard cut" if w == 0 else f"feather w={w}"]) for name in TEX for w in (0, 4, 16, 64)]
print("\nhard cut and feathers: max |S - expected| %.1e, max |G - expected| %.1e, max F %.1e" % (
    max(abs(r[0] - expected(w)[0]) for w, r in cut), max(abs(r[1] - expected(w)[1]) for w, r in cut),
    max(max(r[2], r[3]) for w, r in cut)))
print("multiband G(coarse) / G(fine): " + ", ".join(
    f"N={n} {results['coarse', f'multiband N={n}'][1] / results['fine', f'multiband N={n}'][1]:.3f}" for n in (3, 5, 8)))
print("Poisson: S/Delta %.1e (fine), %.1e (coarse); F_left %.1e (fine), %.1e (coarse); "
      "|F_right - 1| %.1e (fine), %.1e (coarse)" % (
          results["fine", "Poisson"][0], results["coarse", "Poisson"][0],
          results["fine", "Poisson"][2], results["coarse", "Poisson"][2],
          abs(results["fine", "Poisson"][3] - 1), abs(results["coarse", "Poisson"][3] - 1)))

Output:

data: max |column mean| 2.0e-16, max |cos(T_A, T_B)| per column 2.6e-15
data fine: min column RMS / RMS 0.383, values in [0.069, 0.940], fit on A: max|alpha_A - 1| 1.0e-14, max|alpha_B| 1.4e-14
check fine: Laplacian bands of a constant image (levels below the top, N = 3, 5, 8) 0.0e+00
check fine: pyramid reconstruction 0.0e+00, multiband with constant mask 0.0e+00, Poisson residual 1.7e-15, Poisson clone of A + offset vs A 2.0e-13
data coarse: min column RMS / RMS 0.375, values in [0.053, 0.927], fit on A: max|alpha_A - 1| 3.8e-15, max|alpha_B| 2.4e-15
check coarse: Laplacian bands of a constant image (levels below the top, N = 3, 5, 8) 0.0e+00
check coarse: pyramid reconstruction 0.0e+00, multiband with constant mask 0.0e+00, Poisson residual 2.1e-15, Poisson clone of A + offset vs A 2.0e-13
check: Pillow Image.composite vs NumPy feathering on 8-bit copies, max difference 0 codes

fine texture        S/Delta      G   F_left  F_right
  hard cut         1.000e+00   0.000   0.0000   0.0000
  feather w=4      2.500e-01   1.000   0.0000   0.0000
  feather w=16     6.250e-02   4.000   0.0000   0.0000
  feather w=64     1.563e-02  16.000   0.0000   0.0000
  multiband N=3    1.250e-01   0.100   0.0000   0.0000
  multiband N=5    3.023e-02   0.100   0.0000   0.0000
  multiband N=8    2.559e-03   0.100   0.2139   0.5226
  Poisson          1.018e-14   0.353   0.0000   1.0000

coarse texture        S/Delta      G   F_left  F_right
  hard cut         1.000e+00   0.000   0.0000   0.0000
  feather w=4      2.500e-01   1.000   0.0000   0.0000
  feather w=16     6.250e-02   4.000   0.0000   0.0000
  feather w=64     1.563e-02  16.000   0.0000   0.0000
  multiband N=3    1.253e-01   1.206   0.0000   0.0000
  multiband N=5    2.971e-02   1.232   0.0000   0.0000
  multiband N=8    2.373e-03   1.235   0.2142   0.5228
  Poisson          1.092e-14   1.256   0.0000   1.0000

hard cut and feathers: max |S - expected| 2.4e-15, max |G - expected| 5.0e-13, max F 2.2e-15
multiband G(coarse) / G(fine): N=3 12.080, N=5 12.322, N=8 12.350
Poisson: S/Delta 1.0e-14 (fine), 1.1e-14 (coarse); F_left 9.3e-16 (fine), 7.4e-16 (coarse); |F_right - 1| 6.5e-13 (fine), 6.8e-13 (coarse)

Every check held: the data properties hold to rounding error (largest column mean $2.0 \times 10^{-16}$, largest cosine between the textures $2.6 \times 10^{-15}$, smallest column RMS 0.375 of the overall RMS against a floor of 0.2, values inside $[0, 1]$, projection error at most $1.4 \times 10^{-14}$); the bands of a constant image, the pyramid reconstruction and the constant-mask blends are exact ($0.0$); the Poisson residual is at most $2.1 \times 10^{-15}$ and the clone of $A + \Delta$ returns $A$ to $2.0 \times 10^{-13}$; and Pillow’s composite matches the NumPy feathering exactly (0 codes). The tolerances were fixed in advance: $10^{-12}$ for the exact quantities, $10^{-10}$ for the Poisson solve, $10^{-9}$ for the derived feathering and Poisson values, a floor of 0.2 for the column RMS, and zero codes for Pillow.

Feathering (H1: held, in both textures). The hard cut and the three ramps reproduce the derivation to $2.4 \times 10^{-15}$ in $S$ and $5.0 \times 10^{-13}$ in $G$, with $F$ at $2.2 \times 10^{-15}$. Reading down the rows, each fourfold widening of the ramp divides $S$ by four and multiplies $G$ by four: making the seam less steep costs texture overlap in exact proportion. The rows are the same for both textures, since feathering treats every frequency alike.

Multiband against feathering (H2: held). For the fine texture, the 8-level blend has $S = 2.56 \times 10^{-3}$, about a sixth of the widest ramp’s $1/64$, and $G = 0.100$, a tenth of the narrowest ramp’s. No ramp comes close to both at once; their product is about $2.6 \times 10^{-4}$, against $1/4$ for every ramp. This is, in numbers, the behaviour the image pyramid article describes for Burt and Adelson’s spline: the offset, which lives in the coarse bands, crosses over slowly, while the fine texture, which lives in the fine bands, switches almost as sharply as in a hard cut.

Wavelength (H3: held at every $N$). The coarse texture’s overlap is about twelve times the fine texture’s at every depth (12.080, 12.322, 12.350), well above the factor of two the hypothesis asked for. The coarse texture sits in coarser bands, whose weights cross over across wider zones, so its two textures overlap over more columns. Multiband blending does not escape the trade-off for texture that is itself coarse: its $G$ of about 1.2 is slightly larger than the 4-pixel ramp’s 1.0, although its $S$ is smaller (half the 4-pixel ramp’s at three levels, under a hundredth at eight).

The fine texture’s $G$ prints 0.100 at all three depths. By construction, the two finest bands and the three finest mask levels are the same in the 3-, 5- and 8-level pyramids, and only the coarser part differs; the printed values show that the extra coarse levels leave the fine texture’s overlap unchanged at three decimals. The example does not show why; in our reading, little of a texture with periods of 3–5 pixels survives the two smoothing steps that lead to the third level.

Depth (H4: held, both parts, in both textures). $S$ falls from 0.125 at three levels to about 0.030 at five and 0.0026 at eight for the fine texture, and almost identically for the coarse one. At three levels the seam is as steep as an 8-pixel ramp’s ($1/8 = 0.125$), but with a twentieth of that ramp’s texture overlap for the fine texture (0.100 against $8/4 = 2$); at five it is as steep as a ramp about 33 pixels wide ($1/0.0302$). The price of depth appears in $F$: up to five levels both far sides keep their own brightness ($F = 0$), but at eight levels the offset is spread across the whole image, and the far columns have moved by 0.21 $\Delta$ on $A$’s side and 0.52 $\Delta$ on $B$’s. This is Hessel’s trade-off in another form: the residual level decides how far the offset spreads. The two sides do not move equally. In our reading, the cause is the sampling grid: each REDUCE keeps the even columns, so the $1 \times 2$ top level holds samples taken at columns 0 and 128, and column 128 is the first column of $B$, next to the seam; the coarsest weights are therefore not symmetric about the seam.

Poisson (H5a: held; H5b: held). Poisson cloning removes the brightness step completely ($S$ about $10^{-14}$, the solver’s rounding error), but by moving the whole offset into $B$’s side: $F_\text{right} = 1$ to $7 \times 10^{-13}$, while $A$’s side is untouched. On $B$’s side the picture keeps $B$’s texture and takes $A$’s brightness. Its texture overlap is 0.353 for the fine texture and 1.256 for the coarse, larger for the coarse texture as H5b expected. By construction the harmonic correction starts from $A - B$ in the boundary column, which contains both textures there; in our reading, that mismatch carried into $B$’s side is the overlap, and the coarse mismatch reaches further. Both are larger than the multiband overlaps (0.100 and about 1.2), and the fine one is about a third of the 4-pixel ramp’s.

The table also says what each method is for. If the two sides must keep their own brightness, as when both are trusted measurements, feathering or a shallow pyramid keeps $F = 0$, and at the same seam steepness the pyramid does so with less texture overlap (0.100 and 1.206 at three levels against 2 for the 8-pixel ramp). If the inserted region should take on its surroundings’ brightness, as in cloning an object into a new photograph, Poisson blending does exactly that, without any brightness step.

Variants

• Alpha compositing. When the weight is a genuine opacity rather than a design choice, blending becomes compositing with a matte, the forward model of image matting. Matting estimates the opacity and the foreground colour; blending chooses the weights.
• Exposure fusion. Mertens et al. (2007) blend a bracketed set of exposures with per-pixel weights, the product of contrast, saturation and well-exposedness measures, normalized to sum to one, through the Laplacian pyramid (Hessel, 2018, §2, eqs. 4–6). HDR imaging describes the measures and where exposure fusion fits against radiance merging and tone mapping. Hessel describes how many pyramid levels the authors’ implementation uses (§3, eq. 10) and analyses its halo and range artefacts (§8).
• Choice of guidance field. For gradient-domain blending, Di Martino et al. (2016, eq. 19) compare replacing the destination’s gradients with the source’s, averaging the two, and keeping the larger at each pixel.
• Gradient-domain stitching. Levin, Zomet, Peleg and Weiss (2004), as Di Martino et al. (2016, §1.1) describe their work, stitch overlapping images with cost functions defined on gradients, combining gradients rather than intensities to reduce global inconsistencies and the effects of illumination and camera response.
• Fourier solvers. Morel, Petro and Sbert’s Fourier implementation, as Di Martino et al. describe it (§1.1, §3.2), solves the Poisson equation on the whole image with Neumann conditions at its edges, which allows an exact and fast solution even when the region has a complex shape; Morel et al. showed empirically that it is perceptually equivalent to the restricted solution outside the region.
• Learned blending. Face-swapping systems may train a blending network: deepfakes describes FSGAN, whose blending network is trained with a Poisson blending loss.

Best Practices

• Match the transition to the content’s scale. In our reading, the example’s rows are the design rule in miniature: a transition narrower than the brightness difference’s scale leaves a seam, and one wider than the texture’s scale leaves overlap. Multiband blending is the method that sets the width per scale.
• Choose the pyramid depth for the offset, not for the image size. A deeper pyramid spreads a brightness offset further. If the far parts of both images must keep their brightness, stop the pyramid before the residual spans the whole image (in our reading of the example and of Hessel, 2018, §8).
• Decide where the disagreement should go. Weighted methods share it between the two sides; Poisson cloning hands all of it to the inserted region. Neither is wrong; the application decides.
• Check the pyramid on a constant. A pyramid reconstructs its input exactly whatever its border handling, so reconstruction tests do not detect a faulty EXPAND. The Laplacian bands of a constant image should be zero; the example prints this check.
• Blend in linear light when brightness must be physically right. Image matting reports Szeliski’s advice to linearize colours before matting; in our reading the same applies to any weighted average of two images, since averaging gamma-encoded values does not average light.
• Register first. Every method here assumes the two images agree where they are joined; alignment errors belong to image registration, not to blending.

Pitfalls

• Feathering too wide. A wide ramp hides the seam but shows both textures across the whole ramp, as the example shows. With images of the same scene that are slightly misaligned, the same mixing produces doubled edges (our reasoning; the example uses unrelated textures).
• Out-of-range values. Pyramid reconstruction can leave the input range even when the weights sum to one; clipping saturates, and rescaling loses contrast (Hessel, 2018, §8, quoting Mertens et al.).
• Halos from the residual level. A shallow pyramid leaves a thin, visible low-frequency transition; a deep one spreads the offset across the image (Hessel, 2018, §8; the example’s $N = 3$ and $N = 8$ rows).
• Border handling in EXPAND. EXPAND inserts zeros between samples and smooths. The image’s extension at its edges must keep the inserted samples and zeros on their own grid, or EXPAND will not return a constant image near the border. The example uses reflection about the edge pixel for this reason.
• Sampling-grid asymmetry. With power-of-two sizes and subsampling that keeps even indices, the coarse levels do not sit symmetrically about a seam in the middle (our reading of the example’s $N = 8$ row). The image pyramid article’s pitfalls on odd sizes and level alignment apply.
• Poisson colour shifts and bleeding. Poisson cloning changes the inserted region’s brightness and colour to match its boundary (the example’s $F_\text{right} = 1$). Where the boundary crosses content that does not match, the mismatch spreads into the region; Di Martino et al. (2016, §1.1) describe later work by Tao et al. that reduces such bleeding artefacts by placing the residual errors where they are least visible.

Evidence and Impact

Multiband blending is one of the compositing methods of panorama stitching (image pyramid, image warping), and exposure fusion is built on it. Hessel (2018, §2) reports that in his experiments, after an affine rescaling of the result to $[0, 1]$, the fused images generally presented no artefacts, and he analyses the cases with a low-frequency halo and with out-of-range values (§8). Poisson editing has grown into a family of gradient-domain methods: Di Martino et al. (2016, §1.1) review applications to stitching, matting, image fusion, completion and video editing that followed Pérez et al. (2003).

Blending also leaves traces that detectors use. The face X-ray of Li, Bao et al. (2020), as deepfake detection describes it, rests on the observation that most face-manipulation methods blend an altered face into a background through a soft mask $M$, and marks the blending boundary with $4M(1 - M)$, a map that is zero inside and outside the mask and peaks where both sources contribute. Like the texture overlap $G$ of the example, it is large exactly where the two sources are mixed.

The example shows the limits of what any single blend can do. Feathering trades a seam for overlap at a fixed rate; multiband blending beats that rate by handling each scale separately but spreads a brightness difference further the deeper the pyramid goes; and Poisson blending removes the seam entirely by changing the inserted image. The choice is about where the two images’ disagreement should end up.