Applying Modern PDE Techniques to Digital Image Restoration

Inpainting, or image interpolation, is a process used to reconstruct missing parts of images. Artists have long used manual inpainting to restore damaged paintings. Today, mathematicians apply partial differential equations (PDEs) to automate image interpolation. The PDEs operate in much the same way that trained restorers do: They propagate information from the structure around a hole into the hole to fill it in [1,2,3].

While artists can work directly on a painting, a PDE requires a mathematical representation of the subject matter, such as a digital image. A digital image is essentially a 2D matrix of integers, with each integer representing the color or grayscale value of an individual pixel. Holes in the image are represented by unknown values in the matrix. PDE-inpainting fills in these missing values based on the values of nearby pixels.

One of my first explorations of the use of PDEs for inpainting focused on the restoration of 14th-century frescoes. I have found that high-order, nonlinear PDEs enable a wide range of other applications for inpainting, in fields as diverse as medicine, astronomy, and cartography.

For me, MATLAB^® is a natural choice for inpainting with PDEs. MATLAB enables me to rapidly develop and test preliminary algorithms. The MATLAB operations in my scripts have a virtually one-to-one relationship with the steps in the algorithms I develop. This level of correspondence is all but impossible to achieve using a lower-level programming language. I can easily visualize the results in many different formats using graphs and plots, a type of graphical analysis that, with a lower-level language, often requires a second tool.

Simple Inpainting: Fourier's Heat Equation

Fourier's heat equation describes how temperature in a material changes. Given this PDE, knowledge of the thermal conductivity of the material, and knowledge of the initial temperature conditions u(t=0), it is possible to calculate the temperature at any point in the material for any given time by

\[u_t = \Delta u\]

The heat equation is also known as the diffusion equation because it describes how a value (in this case, temperature) diffuses across an area over a time interval.

With grayscale values replacing temperature, a slightly modified version of this equation can be used for basic inpainting. I implemented a simple MATLAB program that uses this technique. This program can restore an image with holes cut in it by propagating (or diffusing) the grayscale values surrounding each hole toward the hole's center (Figure 1).

This approach works well only when the areas surrounding each hole are homogenous. When the hole spans a sharp edge, the edge pixels are diffused, causing the edge to be lost in the restored section (Figure 2).

Handling Nonlinearity: The Perona-Malik Equation

Because the heat equation possesses only smooth solutions, it cannot preserve discontinuous image features such as edges. In particular, it cannot propagate edges across large holes because it relies on linear interpolation. In Figure 2, these edges (the dark parallel lines on either side of the hole) are not propagated into the hole by the inpainting process.

To address this shortcoming, I created a second MATLAB program that uses nonlinear diffusion. This program uses total variation flow [2], a technique in which the diffusivity constant depends upon the size of the image gradient, reducing the amount of diffusion that occurs near edges.

Starting with the initial guess \(u(t=0)\) for the image inside the hole, an inpainted image \(u(t=T)\) is constructed by evolving

\[u_t = \mathrm{div}\left(\frac{\nabla u}{|\nabla u|}\right)\]

This program, which uses a form of the Perona-Malik equation [4], preserves edges in an image and to some extent propagates them into missing regions (Figure 3).

Although total variation flow is a clear improvement over the heat equation approach, it still has disadvantages. Because the Perona-Malik equation is a second-order PDE, it does not perform well on edges that span large gaps or holes in an image (Figure 4). The equation cannot incorporate information about the direction in which to continue an edge. As a result, the edge fades away into the hole.

Factoring in Gradients: Fourth-Order Total Variation Equation

Higher-order PDEs address the deficiencies of the total variation flow approach. With a MATLAB program that uses fourth-order PDEs, we can take two boundary conditions into account instead of just one. Specifically, instead of factoring in only the grayscale value of the pixels near a hole, the program also factors in the gradient, which enables it to model the direction in which edges should be propagated into the hole.

Material science researchers often use the Cahn-Hilliard equation to model phase separation, a process in which the constituents of a fluid mixture spontaneously separate, coarsening the phases. The Cahn-Hilliard equation has been very successfully applied to the inpainting of binary structures with large gaps [5]. A generalization of this equation for grayvalue images [6] can be shown to be equivalent to a fourth-order version of total variation flow:

\[u_t = - \Delta \left(\mathrm{div}\left(\frac{\nabla u}{|\nabla u|}\right)\right)\]

I implemented this fourth-order total variation equation in MATLAB to reconstruct the edges even within large holes (Figure 5). The complete algorithm comprises about 200 lines of MATLAB code, and it uses only basic matrix operations and calls to fft() to perform fast Fourier transforms.

Figure 6 shows more elaborate examples using this inpainting equation. Although the number of missing pixels in the first example is only one tenth of the number of missing pixels in the second example, the correct inpainting is much harder to achieve for the former. This is because the diameter of the inpainting domain (the width of the holes) is the crucial factor rather than their area. Especially in the more difficult situation of the first example in Figure 6, fourth-order and nonlinear PDEs are of utmost importance.

The solution to the fourth-order total variation equation for inpainting is both space-dependent and time-dependent. Using an explicit method to solve this fourth-order PDE would be computationally prohibitive: Each time step involves complex matrix multiplication operations, and it would take thousands of time steps to ensure that the algorithm converged. The MATLAB algorithm that I developed uses a semi-implicit method for time discretization that linearizes the problem at every time step and then solves a linear system of equations [7]. For space discretization, the algorithm uses finite difference methods. These methods work well because images are essentially fixed, equidistant grids of pixels.

Applications of PDE-Based Inpainting

My main work with the Cahn-Hilliard inpainting began on a project to restore some 14th-century Viennese frescoes that were rediscovered in 1979 (Figure 7).

I recently worked on a project to map a network of roads from satellite images. Numerous objects, including cars and trees, obscured the roads in the images. Using MATLAB and Cahn-Hilliard inpainting, my colleagues and I developed an algorithm that automatically removed these obstructions from the satellite images. French astronomers are using a similar technique to remove stars in the Milky Way from astronomical images so as to obtain a clearer picture of the galaxies beyond. PDE-based inpainting algorithms are also applicable in the medical field to reconstruct images from highly undersampled MRI measurements.

New Trends in Inpainting

Not all inpainting methods rely on PDEs. For example, exemplar-based inpainting uses a copy-and-paste algorithm to fill in missing areas of an image. If an image has a hole where a missing shirt button should be, for example, an algorithm that uses exemplar-based inpainting could insert the button by copying one from another location in the image; a PDE-based algorithm could not.

Exemplar-based inpainting is a nonlocal technique–it uses information from the whole image, not just from the area surrounding the image. While exemplar-based inpainting provides impressive results in the restoration of textiles and repetitive structures [8], its ability to reconstruct structures without a given example is limited. For this reason, a synthesis of local (PDE) and non-local (exemplar-based) inpainting techniques is desirable, and current research is focusing on this approach.