10.8 Poisson image editing🔗⧉

Cut a toucan out of one photo and drop it into another, onto a branch in a different forest, say. Paste the bird's pixels onto the new background and the composite fails in one very specific way: a crisp halo runs all the way around the toucan, exactly where its pixels abut the new scene's. The bird was lit by a different sun, white-balanced by a different camera, and exposed a touch differently, and every one of those mismatches becomes a sudden jump in brightness and color right at the boundary. The interior of the toucan looks perfect. It is the seam that betrays the forgery.

The obvious patch is feathering, where you blur the boundary so the two images cross-fade, but it does not really work. A blur has a single scale, so it either leaves the broad color mismatch plainly visible or it smears the toucan's beak into the leaves behind it. What we actually want is for the pasted bird to flow out of its new surroundings: to keep its own texture and shading but adopt the new scene's local color and light at the rim, so that there is simply no jump for the eye to catch. That, and many other seamless edits, follows from one idea: do not edit the pixels; edit their gradients, then reconstruct the image whose gradients best match.

Poisson editing fits the source gradients while matching destination boundary values. Those requirements can conflict, so the reconstructed gradients and colors can differ from the source. The method reduces some seams; it cannot guarantee an invisible composite.

10.8.1 Edit gradients, not pixels🔗⧉

Here is the single observation the whole chapter stands on. The human visual system is far more sensitive to local contrast, meaning how much a pixel differs from its neighbors, than to any pixel's absolute value. What makes the toucan look like a toucan is not that one patch of its beak happens to be RGB $(0.91, 0.43, 0.08)$; it is that there is a sharp transition from a dark feather to a bright beak right here. The absolute level can drift. A whole region can grow lighter or pick up a color cast, and as long as the differences between neighbors survive, the eye reads the same scene. This is precisely why lightness constancy works: we judge a surface by its contrast against its surroundings, not by the raw quantity of light reaching us. A sheet of paper looks white in deep shade and in full sun even though the shadowed paper sends our eye far less light than a sunlit lump of coal sends from the next room.

So we move our object of manipulation from the image to its gradient field. Write $\mathbf{v}$ for the gradient field we want the answer to have (the guidance field) and $f$ for the image we solve for. Ideally $\nabla f = \mathbf{v}$ everywhere. The catch, explained in the sidebar just below, is that this is almost never exactly solvable: a field of two numbers per pixel is generally not the gradient of any single image. So we settle for the next best thing and ask for the image whose gradients come closest, in the least-squares sense:

Read back in words: find the image $f$ whose gradients deviate as little as possible, summed over the whole image, from the gradient field we prescribed. This is exactly an inverse problem of the kind in Chapter 6.7, Image formation as a linear system (a quadratic energy in the unknown image) and minimizing a quadratic means setting its derivative to zero, which gives a linear system. That system is the Poisson equation,

where $\nabla^2$ is the Laplacian (the second-derivative operator) and $\operatorname{div}\mathbf{v}$ is the divergence of the guidance field (the net outflow of a vector field); if divergence is rusty, see Refreshers → Calculus. We derive the discrete form in the next section; for now keep its shape in view: we ask that the Laplacian of the output match the divergence of the field we prescribed. When $\mathbf{v}$ really is the gradient of some source image $g$, this just says the output's Laplacian should equal the source's, since $\operatorname{div}\nabla g = \nabla^2 g$.

10.8.2 Seamless cloning: the headline application🔗⧉

The toucan-into-the-forest problem now has a clean recipe, due to Pérez, Gangnet and Blake (2003). State it exactly, because everything else in the chapter varies it. We have a destination image $f^{*}$, a source image $g$, and a region $\Omega$ (the mask) where we want to paste. Rather than copy the source's pixels into $\Omega$, we copy its gradients, set the guidance field to $\mathbf{v} = \nabla g$ inside $\Omega$, and solve for the interior values $f$ that best match those gradients, subject to one constraint: at the boundary $\partial\Omega$ of the region, $f$ must equal the destination. That boundary constraint is a Dirichlet condition,

read back as: pin the edge of the pasted region to the surrounding image's actual pixel values. Inside, solving $\nabla^2 f = \operatorname{div}\mathbf{v}$ fills in values that carry the source's contrasts (the toucan's feathers, the gloss of its beak) yet flow smoothly out of the destination at the seam.

Matching boundary values removes the explicit value discontinuity there. The adjustment inside the mask is generally spatially varying, not just one additive offset, and can alter both color and gradients. Boundary texture, mask placement and guide mismatch still matter.

There is a plain way to see what the solver is doing, and it matters later. Suppose for a moment the destination's boundary values happened to coincide with the source's. Then copying the gradients would simply reproduce the source: no jump, done. In general they do not coincide, so think about the discrepancy: at each boundary point there is a difference between what the destination wants ($f^{*}$) and what the source provides ($g$). The Poisson solve quietly spreads that boundary discrepancy smoothly across the interior (it adds a gentle correction surface, determined by the mismatch around the whole boundary) so the patch keeps much of the source structure while changing its gradients by the gradient of that correction. We will call that correction surface the membrane and use it for speed near the end of the chapter.

A note on encoding, following our standing rule of declaring the space every edit lives in: for an ordinary clone where source and destination are similarly exposed, the usual gamma values are fine. But when tones or exposures differ substantially, especially for gradient-domain HDR, we integrate in log, for reasons made precise in the next section.

10.8.3 Reconstruction = solving the Poisson equation🔗⧉

To turn $\nabla^2 f = \operatorname{div}\mathbf{v}$ into something a computer can solve, we discretize it, and the discretization is satisfyingly plain. Approximate the gradient with finite differences, meaning the difference between a pixel and its left neighbor, and between a pixel and its top neighbor. We deliberately use the simplest kernels, $[-1,\,1]$ and its vertical transpose, not the wider Sobel kernels: Sobel is tuned for analyzing gradients and averages over too large a support, whereas here we want to stay as close as possible to the raw pixel differences.

Use a discrete gradient and its transpose from the same boundary convention. The positive-semidefinite normal operator $D^\top D$ has the interior stencil $[-1,2,-1]$ in one dimension and is the negative of the usual discrete Laplacian. With edge guidance defined as $v_{pq}=g_p-g_q$, the four-neighbor interior equation is

where $v_{pq}$ is the prescribed difference across the edge from $p$ to its neighbor $q$ (for a clone, $v_{pq}=g_p-g_q$, the source's own difference). A note on the sign convention. We use the positive discrete Laplacian, $|N_p|\,f_p-\sum_{q} f_q$ (the standard $[-1,2,-1]$ stencil) which approximates $-\partial^2/\partial x^2$, the negative of the continuous operator. Under this convention both sides of the continuous identity $\nabla^2 f=\operatorname{div}\mathbf v$ flip sign together, so the linear system above is unchanged; this matches Pérez. Read the equation back: four times each interior pixel's value, minus the sum of its four neighbors, should equal the prescribed difference sum there. Stack one such equation per pixel and you have a large, sparse, structured linear system, sparse because each equation touches only five pixels and structured because every row is the same Laplacian stencil, just shifted (a convolution). Only the right-hand side depends on the guidance field; the matrix does not.

The system needs boundary conditions to be well-posed, and which one you use depends on the task. For cloning we already have it: the Dirichlet condition $f|_{\partial\Omega}=f^{*}|_{\partial\Omega}$ pins the region's boundary to the destination. When instead we integrate a guidance field over a whole image with no enclosing destination, as in gradient-domain HDR, there is no boundary to pin, and we use a Neumann condition (zero normal derivative at the image border), which fixes the overall brightness only up to a constant.

The single equation per pixel already hands us the most transparent way to solve the system. Solve the per-pixel relation $4 f_p - \sum_{q\in N(p)} f_q = b_p$ for $f_p$ (where $b_p = \sum_{q\in N(p)} v_{pq}$ is the divergence right-hand side, computed once) and each pixel is just the average of its four neighbors plus a quarter of $b_p$. Sweep that update over every interior pixel, holding the boundary fixed, and repeat until the image stops changing:

b[p] = ∑ over q in N(p) of v[p, q]      # divergence, computed once
initialize f (e.g. f = destination, or zero in the interior)
while not converged:
    for each interior pixel p:
        f[p] = ¼ · (∑ over q in N(p) of f[q] + b[p])
    keep f fixed on the boundary
return f

Read it back: each pass nudges every pixel toward the average of its neighbors, carrying the prescribed divergence with it. The in-place update shown is Gauss–Seidel. Jacobi reads every neighbor from a separate previous-iteration buffer. It is correct but slow because it propagates information inward from the boundary one pixel-ring per sweep, so a large region takes many passes to settle. That sluggishness is exactly what the better solvers below fix. As a sanity check on the whole construction, set the guidance field to zero ($\mathbf{v}=0$): the Poisson equation collapses to Laplace's equation $\nabla^2 f = 0$, whose solution inside a region, with the boundary pinned, is the smoothest possible interpolation of the boundary values, a membrane stretched across the hole. In 1-D that membrane is just a straight line between the two endpoints; in 2-D it is the soap-film surface those boundary values support. Hold onto that membrane; it returns.

Stop using a residual criterion and an iteration cap, not only a fixed visual guess. For the constrained system $Af=b$, report $\|Af-b\|/\max(\|b\|,\epsilon)$ together with an absolute tolerance. Keep fixed boundary values out of the unknown vector or account for their contributions on the right-hand side. Use the actual neighbor count where the computational domain truncates the stencil.

The same solver as the regression chapter🔗⧉

Never form a dense matrix. A one-megapixel region has a million unknowns, and a dense Laplacian would hold $10^{12}$ entries, almost all of them zero. There are two workable representations. A sparse matrix, with five nonzeros per row, is what a direct factorization needs, and it suits the case where the same system is solved many times with different right-hand sides. For iterative solvers you can go further and never build a matrix at all: keep the image as an image, and use the linear algebra only to say which operations to perform. Each of them turns out to be an ordinary image operation. Applying the Laplacian to the current guess is a convolution with the 5-point stencil, restricted to the mask. The right-hand side is the divergence of the guidance field, computed once. The inner products that conjugate gradient needs are sums of pixelwise products over the region. The whole solve reduces to repeated stencil convolutions, a few image-wide sums, and scalar arithmetic. Whichever representation you choose, use the same boundary convention in every operator: the stencil, the right-hand side and the residual must agree about which pixels are unknown and how many neighbors each has.

This is exactly the machinery of Chapter 6.7, Image formation as a linear system, and the same three solver choices apply:

• Conjugate gradient (CG), optionally with a preconditioner, is the default solver. It only ever needs to apply the Laplacian to a vector (a convolution), it exploits the system's sparsity and symmetry, and it converges far faster than plain gradient descent, which crawls because it propagates information inward from the boundary one pixel-ring at a time and tends to zigzag.
• Multigrid: coarse-to-fine, solving on a small image first and refining down the pyramid, which fixes precisely gradient descent's slow long-range propagation. This is the pyramid idea of Chapter 6.6, Pyramids and wavelets put to work as a solver.
• An FFT solver on a rectangular periodic or Neumann domain. The Laplacian is a convolution, and convolution diagonalizes in Fourier (the sine waves are its eigenvectors), so the Poisson solve becomes a per-frequency division: one FFT, a divide, and one inverse FFT. This ties straight back to Part 6, Linear methods and Fourier analysis.

Colorization, matting and regularized deblurring also lead to sparse quadratic systems. Their operators differ, but understanding symmetry, nullspaces and preconditioning here prepares us to solve those systems too.

On a periodic rectangle the discrete Laplacian is diagonalized by the discrete Fourier transform. Rectangular Neumann or Dirichlet problems use compatible cosine or sine transforms and different boundary spectra. An arbitrary cloning mask does not become diagonal under one ordinary image FFT.

For the periodic unit-grid five-point Laplacian, the spectrum is $2\cos\omega_x+2\cos\omega_y-4$. The continuous expression $-(\omega_x^2+\omega_y^2)$ belongs to a different discretization and should not be substituted silently. At zero frequency, enforce the compatibility condition and choose the mean separately. An inconsistent nonzero mean on the right-hand side cannot be solved by simply dividing by an epsilon. Figure 10.8.4⛶ fullscreen should use this same discrete convention.

Which space to integrate in: linear vs. gamma vs. log🔗⧉

The interactive playground makes all of this concrete (Figure 10.8.5⛶ fullscreen): it lets you paste one photo into another, paint the paste region, and flip the solve between linear, gamma, and Chong's illumination-invariant log color space, and toggle mixed gradients, watching the seam dissolve and the re-leveling change as you switch. Keep it open as you read the next points.

State the processing domain and its interpretation before fitting gradients. Linear, encoded and log-domain differences define different objectives.

Linear-light gradient fitting penalizes additive radiometric differences. This can be appropriate, but it gives large radiances more influence than a relative-contrast objective. The domain is a modeling choice, not a universally right or wrong encoding.

For positive radiances, a log-domain fit treats ratios as differences: $\log 200-\log100=\log2000-\log1000$. Define the log base and a floor for zero or negative values. That floor changes dark-region behavior, and exponentiation can amplify errors. Compare the result after the intended output transform; log fitting does not guarantee exposure-independent perception.

A perception-based color space🔗⧉

The linear-versus-log point is really a special case of a broader question: in what representation does gradient magnitude best track perceived contrast? Hamilton Chong, Steven Gortler and Todd Zickler (2008) answer it for color, building a color space designed so that the gradient's magnitude is covariant with perceived contrast, a perception-based, illumination-invariant space in which to take and integrate gradients. Doing the Poisson solve there is the best-behaved choice when color, not just luminance, is at stake; it generalizes "use log" from a single channel to the full color manifold. We will not need its machinery, but it is the principled endpoint of the encoding discussion: choose the space so that the thing you are matching, the gradient, means the same to the math as it does to the eye.

One consequence is immediate, and it is the property that makes such a space "best-behaved" (Figure 10.8.6⛶ fullscreen). Clone the same source into the same destination, but vary the source's exposure by $\pm 2$ EV first. In linear light the cloned gradients scale with the exposure, so the re-leveled paste shifts and clips from one exposure to the next; in gamma it shifts too, differently. But in a log (and, for color, Chong) space a $\times 2^{\text{EV}}$ exposure is a constant additive offset, and the Poisson solve sees only gradients, which an offset leaves untouched, so the cloned result is identical at every exposure under that idealized model, provided scaling causes no clipping and no value crosses a log floor. The encoding has bought exposure-invariance for free: the edit no longer depends on how bright the source happened to be.

The chosen intensity domain changes which differences are preserved. A change from 100 to 200 and one from 1000 to 2000 have different additive differences but the same ratio. Their log differences are equal. Integrating gradients of log intensity can therefore be useful when the intended structure is relative contrast, as in gradient-domain tone mapping. Integrating linear values answers a different question. Zeros, negative values, and the absolute-level constraint require explicit treatment; changing the domain does not by itself supply the missing boundary conditions.

10.8.4 Poisson vs. pyramid blending🔗⧉

That membrane is also the cleanest lens on how Poisson blending differs from the other classic seam-hider, Laplacian-pyramid blending (Chapter 6.6, Pyramids and wavelets, Burt & Adelson 1983 (spline blending)). Both make a paste seamless by treating slow variation differently from fine detail. They differ on one point: what becomes of the source's absolute level, its DC.

Poisson is gradient-domain with a pinned boundary. It reproduces the source's gradients and lets the membrane absorb whatever offset the boundary demands, so the source's DC is discarded and reset by the destination. Adding a constant to the source changes nothing. Laplacian-pyramid blending instead cross-fades the actual pixel values, band by band, with the mask feathered wider at coarser scales; the coarsest band contains the DC and any remaining low-frequency variation, so the source's brightness is kept, merely feathered against the destination.

Paste a flat, constant patch and the two diverge completely (Figure 10.8.7⛶ fullscreen). For Poisson the source gradient is zero, so the solve collapses to the bare membrane, meaning Laplace's equation with the destination on the boundary, and the patch's value is irrelevant: the region is filled by smoothly interpolating its surroundings (a straight line in 1-D, a soap film in 2-D), and the constant vanishes. The pyramid has no boundary to satisfy; its coarsest band simply copies the source DC, so a softened version of the constant survives, feathered into the destination over the wide coarse mask. It is this chapter's big lesson L10.4show (the eye cares about gradients, not absolute values) seen from both sides: Poisson keeps only the gradients and throws the level away, while the pyramid keeps the level and blends it gently. The practical consequence follows directly: Poisson re-lights the insert to match its new surroundings (seamless, but it washes out or bleeds when the boundary varies wildly, and it erases flat regions), whereas the pyramid preserves the insert's appearance and only dissolves the seam (robust and solver-free, but a genuine brightness mismatch lingers as a soft halo). Cloning wants the first; panorama stitching (where both halves are equally real) wants the second.

10.8.5 Advanced techniques🔗⧉

Everything so far changed the image by changing the guidance field $\mathbf{v}$ and re-using one solver. The gradient-domain view is that a long list of seemingly different edits are all just different choices of $\mathbf{v}$.

Mixing gradients🔗⧉

In a plain clone we set the guidance field to the source's gradient everywhere inside the region. But sometimes the destination's own texture should show through, as when we paste a holey or lacy object, a transparent overlay, or text onto a textured surface, we want the background to fill the gaps rather than be flattened by the source's near-zero gradients there. Pérez's mixing-gradients variant handles this with a one-line change: at each pixel, keep whichever gradient is stronger, the source's or the destination's,

then solve the same Poisson system. Note carefully what is being selected: the comparison of magnitudes decides which vector gradient to copy, and the whole vector (both the $x$- and $y$-component) comes along. It is a vector selection, not a per-component scalar maximum: we never mix the $x$ of one source with the $y$ of the other. Read back: at each pixel, take the gradient of whichever image has more going on there. Where the source has strong structure (the strokes of the text), its gradients win and the source shows; where the source is flat but the destination has texture (the gaps between strokes), the destination's gradients win and the background shows through. The seam still vanishes, because we still solve Poisson with the destination boundary, we have only changed which gradients we ask the result to match (Figure 10.8.8⛶ fullscreen).

Screened Poisson: pinning the low frequencies🔗⧉

There is a blind spot built into everything above. Solving the Poisson equation matches the image's gradients, which is exactly why it excels at edges and fine detail. But gradients say nothing about the absolute level. Only the constant mode lies in the gradient operator's nullspace on a connected domain. Other low frequencies have small responses, so a slowly varying error can be penalized weakly even though it is not invisible to the operator. That is the same fact behind "Poisson discards the source's DC" above. It means a pure gradient-domain solve leaves the low frequencies underdetermined, pinned only indirectly through the boundary. When there is no trustworthy boundary (large regions, whole-image edits, gradient fields that don't integrate cleanly), the result can drift, smear, or low-frequency-bleed.

The fix is the screened Poisson equation: add a data (fidelity) term that asks the result to also stay close to a reference image $f^{*}$, not just match the target gradients. Minimize

whose Euler–Lagrange equation is $(\lambda - \nabla^2)\,f = \lambda f^{*} - \operatorname{div}\mathbf{v}$. This is the Laplacian "screened" by the $\lambda$ term (the same form as a screened/Yukawa potential in physics). The data term is a high-pass-vs-low-pass knob: where the gradient term dominates it reproduces edges and detail as before, but the $\lambda(f-f^{*})$ term anchors the low frequencies and the DC to the reference, so the solution can no longer float away. The weight $\lambda$ trades the two off: small $\lambda$ recovers ordinary Poisson (gradients only), large $\lambda$ pulls back toward the original, and because the screened operator is better-conditioned than the bare Laplacian (no zero eigenvalue at DC), it is also numerically friendlier and admits fast Fourier / multigrid / convolution-pyramid solvers. This data-plus-gradient energy is the workhorse of gradient-domain editing frameworks ([@bhat-etal-2008-screened; @bhat-etal-2010-gradientshop]) and continues to be refined for blending and tone tasks (Screened-Poisson gradient-domain method (BMVC 2017)). The lesson dovetails with L10.4show: gradients carry the structure the eye cares about, but you still need one low-frequency anchor to know where to put it.

The healing brush is a covariant derivative🔗⧉

Photoshop 7's healing brush preceded the Poisson-editing paper. Georgiev's account Georgiev 2004/2005 explains both a perceptual motivation and an implementation distinction: matching values at a boundary need not match the slopes there. His historical implementation uses a biharmonic correction to improve derivative continuity, initialized from a harmonic solution. This is a useful example of a product requirement changing the differential equation, rather than merely changing its solver. See Georgiev's account.

The covariant-derivative argument asks how to compare spatial changes after accounting for a reference image's variation. For a positive reference $g$, its connection uses the relative derivative $\nabla g/g$. That motivates transporting texture in a way that respects changes in illumination. It does not make the resulting energy identical to unweighted least squares on log gradients, and the paper's practical biharmonic algorithm is distinct from that theoretical derivation. These are historical descriptions, not a specification of current Photoshop.

Conservative-gradient simplification: the membrane shortcut🔗⧉

A general Poisson problem needs a linear solve, which is what kept early gradient-domain editing off the interactive path. But recall the sidebar: when the guidance field is conservative, meaning it is the exact gradient of some image $g$, as in a plain clone, the problem simplifies dramatically. We saw the reason already, in the boundary-discrepancy picture: pasting $\nabla g$ and solving Poisson gives exactly the same answer as adding a correction surface to the naive paste, where that correction surface is the smooth membrane that interpolates the boundary mismatch $f^{*}-g$ into the interior. And the membrane solves Laplace's equation,

Then return $g+u$. The harmonic correction has the same interior linear operator as the original Poisson problem, so changing the right-hand side does not automatically improve its condition number. The useful opportunity is to approximate a smooth correction on a coarser representation or reuse a fixed-domain solve Farbman et al. 2009. Check the approximation near complicated boundaries.

Other gradient-domain edits🔗⧉

Once you see editing as change the guidance field, reuse the solver, the family keeps growing, and several of these get their own treatment later in the book:

• HDR tone compression (Fattal, Lischinski and Werman 2002): a high-dynamic-range scene has a few enormous gradients (the window, the lamp) swamping the rest. Attenuate the large gradients (leave small ones alone, shrink big ones) then integrate the modified field. The huge range collapses while local contrast survives, because we manipulated the gradients and let the absolute level fall out of the solve. Integrate in log, where attenuation is well-behaved.
• Photomontage (Agarwala et al. 2004): stitch many photographs into one by choosing, per region, which source contributes its gradients, then solving Poisson over the whole composite. This is the explicit meeting point of the three legs of this part: a graph-cut seam Boykov et al. 2001 (the seam-optimization chapter) decides where to switch between sources, and a gradient-domain blend (this chapter) hides whatever residual mismatch the seam leaves: cut along the edges, then reconstruct across them.
• Gradient-domain panorama stitching (Levin, Zomet, Peleg & Weiss 2004, Seamless Image Stitching in the Gradient Domain): the same idea aimed squarely at panoramas. Rather than blend the overlapping photographs' values, copy each photo's gradients into its region and reconstruct one image whose gradients best match the patchwork: the exposure and white-balance jumps between shots dissolve into an invisible smooth shift, because (again L10.4show) only the gradients are kept. Levin et al. compare several gradient-domain cost functions and optimization schemes and show the seam vanishing where a plain feathered blend would leave a visible band; it is the panorama-scale instance of everything in this chapter.
• And further afield, all built the same way: shadow and reflection removal (zero out the gradients that belong to the shadow edge, re-integrate) and flash / no-flash fusion (borrow gradients from the sharp flash image while keeping the ambient color). Every one of them changes only $\mathbf{v}$ and reuses the same sparse solver.

Decolorization (Color2Gray): keeping isoluminant contrast🔗⧉

One more gradient-domain edit, and a less obvious one: a color-to-grayscale conversion built from this chapter's machinery rather than an edit of a photo. Two colors can be isoluminant, meaning different hues at the same luminance, like a red and a green tuned to match in brightness, and the naive conversion (take the luminance, or desaturate) maps both to the identical gray, erasing a contrast the viewer plainly saw: a red shape on an equal-luminance green field simply vanishes into a flat gray (think of a pie chart, a subway map, or red berries among green leaves). Gooch, Olsen, Tumblin & Gooch 2005 preserve such differences through a least-squares objective over selected pixel pairs:

Here $K$ is the set of compared pairs and $\delta_{ij}$ is a desired signed difference derived from lightness and chrominance in CIELAB. The user can control which chromatic differences become lighter or darker; that ordering is not determined by luminance alone. The additive offset remains free and is chosen to keep the result close to the source lightness.

This shares gradient-domain editing's idea of reconstructing values from desired differences. It is not generally the same nearest-neighbor Poisson problem: the original method can compare all pixel pairs. Long-range comparisons distinguish disconnected isoluminant regions that a purely local objective can miss, at substantial computational cost. The choice of pair set changes the operator and suitable solver. See the original Color2Gray formulation (Figure 10.8.9⛶ fullscreen).