Poisson Image Editing
By Sébastien Boisgérault, Mines ParisTech, under CC BY-NC-SA 4.0
April 25, 2017
Contents
Introduction
Poisson Image Editing refers to a family of methods introduced in (Pérez, Gangnet, and Blake 2003) that rely on the resolution of the Poisson equation to perform seamless editing of images. Typical use cases include: healing of images from damaged or missing data, image enhancements, e.g. the removal of skin blemishes in portraits, heavy image editing, for example the concealement of objects.
The GNU Image Manipulation Program (GIMP) “Heal tool” implements a Poisson Image Editing method1, described in the documentation as “a smart clone tool on steroids”, because:
Pixels are not simply copied from source to destination, but the area around the destination is taken into account before cloning is applied.

Let’s demonstrate Poisson image editing with GIMP. Say that we want to hide the watch in the image below, in a way one cannot guess that there was an item to begin with.

We first select a disk in a region of the image where there is no object, to capture the texture of the background behind the objects. GIMP displays this region with a dashed line (see “GIMP Heal tool in action” figure). We make sure that the disk is big enough to cover a significant part of the watch, and small enough to avoid the objects that are close to it. Then we graft this heal source on top of the lower part of the watch, then on its upper part. The final result is a seamless removal of the watch from the picture.



Modelling the Problem
Consider a circular image fragment ϕ that requires some editing. If the image resolution is good enough, it is sensible to model this fragment as a function ϕ of two continuous variables x and y, defined inside the unit disk D={(x,y)∈R2|x2+y2<1} and on its boundary ∂D.
If we deal with grayscale images, we can assume that ϕ is real-valued, even with values in [0,1] after a suitable normalization of the image intensity. Color images can easily be modelled too, for example as triples of real-valued functions, one for each channel in a RGB decomposition.
We search for a new image fragment ψ, that matches exactly the original fragment ϕ on the boundary ∂D and whose texture matches approximately the texture of a third image fragment χ in D. We assume that the texture of an image fragment is essentially captured by its gradient field; for the sake of simplicity, we measure the dissimilarity of two gradients with the quadratic mean of their difference.
To summarize, we are trying to solve: minψ∫D‖∇(ψ−χ)‖2 with ψ|∂D=ϕ|∂D.
If we set u=ψ−χ and f=(ϕ−χ)|∂D, we end up with a function u that solves the variational Dirichlet problem: the minimization of the Dirichlet energy of u, subject to the Dirichlet boundary condition f: minu∫D‖∇u‖2 with u|∂D=f. We may also search u as a solution of the Dirichlet problem: Δu=0 on D with u|∂D=f as both problems are – at least informally – equivalent. Indeed, given two candidate solutions u and u+ϵ to the variational problem, we have ∫D‖∇(u+ϵ)‖2=∫D‖∇u‖2+2∫D∇u⋅∇ϵ+∫D‖ϵ‖2, hence u is a minimizer if and only if ∀ϵ such that ϵ|∂D=0,∫D∇u⋅∇ϵ=0 Given Green’s first identity, this condition holds if and only if Δu=0 on D.
In the sequel we study the Dirichlet problem in the classic setting:
Definition – Dirichlet Problem.
Given a continuous function f:∂D→R, a classic solution of the Dirichlet problem is a function u:¯D→R which is continuous in ¯D, twice continously differentiable in D, and such that Δu=0 on D with u|∂D=f.Harmonic Functions
Definition – Harmonic functions.
A function u:Ω→R defined in an open subset Ω of C is harmonic if it is twice continuously differentiable and Δu=0 on Ω.There is a partial converse theorem:
Theorem – A harmonic function is locally the real part of a holomorphic function.
If u:Ω→R is harmonic, for any z∈Ω, there is an open neighbourhood V of z and a holomorphic function f:V→C such that u|V=Ref.Actually, we will show a stronger result: the global existence of a holomorphic function if the domain of definition of the harmonic function is simply connected. But first, we need the following lemma:
Lemma.
If u:Ω→R is harmonic, the function f:Ω→C defined by f(x+iy)=∂u∂x(x,y)−i∂u∂y(x,y) is holomorphic.Now, we may prove the converse theorem.
A consequence of the converse theorem:
Theorem – The Maximum Principle.
Let Ω be an open connected subset of C and u:Ω→R be an harmonic function. If u has a maximum or a minimum on Ω, then u is constant.The Dirichlet Problem
We will prove soon that the Dirichlet problem has a unique solution; what we can already prove is the:
Lemma – Uniqueness of the Solution.
There is at most one classic solution u to the Dirichlet problem with continuous boundary condition f.Harmonic/Fourier Analysis
Lemma – Elementary Solutions.
Let n∈N. If f(eiθ)=cos(nθ+ϕ), then u(reiθ)=rncos(nθ+ϕ) solves the Dirichlet problem with boundary condition f.

This result can be readily extended to finite trigonometric series: if for some finite sequence of real-valued coefficients (an) and (ϕn), f:∂D→R can be decomposed as f(θ)=∑nancos(nθ+ϕn) then by linearity, the function u(reiθ)=∑nanrncos(nθ+ϕn) is a solution of the corresponding Dirichlet problem. Now, if f∗ is merely continuous but can be uniformly approximated to the precision ϵ by the finite trigonometric series f: supθ|f∗(eiθ)−f(eiθ)|≤ϵ and if u∗ is a solution to the Dirichlet problem with boundary condition f∗, then by linearity, u∗−u is a solution to the Dirichlet problem with boundary condition f∗−f; the maximum principle then provides sup|z|≤1|u∗(z)−u(z)|≤ϵ. Now, Fejér’s theorem ensures that the approximation of f∗ by finite trigonometric series can be achieved with an arbitrary small precision ϵ>0. Hence, we can build an abritrarily precise uniform approximation of the solution to the Dirichlet problem.

The Poisson Kernel
The main result of this section:
Theorem – Solution of the Dirichlet Problem.
The Dirichlet problem with a continuous boundary condition f has a unique classic solution u, given in the unit disk by the Poisson integral u(reiθ)=12π∫π−π1−r21−2rcos(θ−α)+r2f(eiα)dα.Let’s start with some definitions:
Definition – Poisson Kernel.
The Poisson kernel is the function P defined in the unit disk by P(reiθ)=1−r21−2rcosθ+r2.Definition – Poisson Integral Operator.
The Poisson integral operator P maps a continuous function f defined on the boundary of the unit disk to the the function P[f] defined inside the unit disk by P[f](reiθ)=12π∫π−πP(rei(θ−α))f(eiα)dα.Accordingly, the main result of this section may be restated as:
Theorem – Solution of the Dirichlet Problem.
Let f:∂D→R be continuous. The function u:¯D↦R defined as u|D=P[f] and u|∂D=f is the unique classic solution of the Dirichlet problem with boundary condition f.The proof of this fundamental result requires several lemmas.
Lemma – Poisson Kernel, Alternate Representations.
The Poisson kernel satisfies: P(z)=11−z−11−1/¯z=Re[1+z1−z]. Hence, it is harmonic.Lemma – Poisson Kernel, Fourier Series.
The Poisson kernel has the following (locally uniformly convergent) Fourier expansion P(reiθ)=+∞∑n=−∞r|n|einθ.Lemma – Poisson Integral, Harmonicity.
Let f:∂D→R be continuous. The function P[f] is harmonic in D.Lemma – Poisson Integral, Boundary Values.
For any continuous function f:∂D→R, the function u:¯D↦R defined as u|D=P[f] and u|∂D=f is continous on ¯D.Proof (see e.g. Rudin (1987)).
Appendix
Linear Gradient
What happens if you heal a region in the left of the image below with a source taken from the right of the image?

Analytic Boundary Condition
We search for approximate solutions to the Dirichlet problem associated with f(eiθ)=14[1+35−4cosθ].
Check that f is defined and continuous in ∂D; compute its range.
Find a holomorphic function, defined in C∗, that extends f and show that it is unique; we also denote f this extension.
Show that f(z)=12[12−z+12−z−1]; determine the Laurent series expansion of f in a neighbourghood of ∂D.
Find a finite trigonometric series that approximates f with the precision ϵ=10−2; provide an approximate solution of the Dirichlet problem associated to f with the same precision.
References
Notes
the GIMP developers state that they use the method described in (Georgiev 2005), a variant of the original design that is invariant under relighting, but a closer examination of the source code shows that they actually use the original method.↩
for example, consider a polyline γ=[z0→z1→⋯→zn−1→zn] of Ω that joins z0 and zn=z; we may assume that it is simple (that the function γ is injective). If r>0 is small enough, the open connected set V=γ([0,1])+D(0,r) is included in Ω, connected and simply connected.↩