Directional Mean Curvature for Textured Image Demixing

by   Duy Hoang Thai, et al.
Duke University

Approximation theory plays an important role in image processing, especially image deconvolution and decomposition. For piecewise smooth images, there are many methods that have been developed over the past thirty years. The goal of this study is to devise similar and practical methodology for handling textured images. This problem is motivated by forensic imaging, since fingerprints, shoeprints and bullet ballistic evidence are textured images. In particular, it is known that texture information is almost destroyed by a blur operator, such as a blurred ballistic image captured from a low-cost microscope. The contribution of this work is twofold: first, we propose a mathematical model for textured image deconvolution and decomposition into four meaningful components, using a high-order partial differential equation approach based on the directional mean curvature. Second, we uncover a link between functional analysis and multiscale sampling theory, e.g., harmonic analysis and filter banks. Both theoretical results and examples with natural images are provided to illustrate the performance of the proposed model.



page 29

page 30

page 31

page 32

page 33

page 34

page 35

page 36


Multiphase Segmentation For Simultaneously Homogeneous and Textural Images

Segmentation remains an important problem in image processing. For homog...

A novel variational model for image registration using Gaussian curvature

Image registration is one important task in many image processing applic...

Improved Image Deblurring based on Salient-region Segmentation

Image deblurring techniques play important roles in many image processin...

Directional Global Three-part Image Decomposition

We consider the task of image decomposition and we introduce a new model...

Image Reconstruction via Discrete Curvatures

The curvature regularities are well-known for providing strong priors in...

Fast Multi-grid Methods for Minimizing Curvature Energy

The geometric high-order regularization methods such as mean curvature a...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.


Image Deblurring, Demixing Problem, Variational Calculus, Feature Extraction, Harmonic Analsysis, Sampling Theory, Multiscale Analysis, Mean Curvature, High-order PDE

1 Introduction

Processing textured images is difficult, but such images are common in many applications. For example, as part of the 2015-2016 forensic statistics program at the Statistical and Applied Mathematical Sciences Institute, much of the research dealt with images of fingerprints, balistic striations on bullets, and shoeprints, and all of these are typically textured.

There is an extensive literature on methods for processing piecewise smooth images; e.g., image restoration by functional analysis [1, 2, 3, 4, 5, 6, 7, 8], image representation by harmonic analysis [9, 10, 11, 12, 13, 14, 15, 16] and image segmentation [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. Much less has been done with textured images.

Instead of being piecewise smooth, textured images have fine scale discontinuities, so that neighboring pixels in the image may take very different values in gray scale or color space. Examples are a shoeprint in sand, where the grainy structure of the sand provides the texture, or a fingerprint on leather, where the leather’s texture interacts with the signal, see [28].

To demix (i.e., to simultaneously deblur and decompose) textured images, we substantially extend the work in Thai and Gottschlich [29]. That paper provided a directional decomposition of gray-scale images into three parts, consisting of piecewise smooth structure , texture , and fine scale residual structure . The solution required minimization of the directional total variation norm, which entails the solution to a second-order partial differential equation.

The new approach is based upon [30], and also focuses upon gray-scale images. But now it decomposes the image into four parts, adding fine-scale noise structure to the previous three distinctions. It also models a blur operator (for simplicity, is assumed to be known, which is reasonable since in principle it can be measured for any camera in a specific application). Thus the demixing problem is formulated as

where is a convolution operator.

The solution minimizes a function of several well-chosen norms and uses directional mean curvature rather than a directional decomposition. This entails solution of a fourth-order partial differential equation (cf. Zhu and Chan [6, 22]) which addresses the “staircase effect” in the TV-L2 model [3]. Other high-order partial differential equation (PDE) approaches for image reconstruction are given in [31, 32, 33, 34, 35]. Importantly, this work finds a new connection between the minimization problem in demixing and multiscale harmonic analysis. This connection enables a more general theory and a new solution strategy.

Our model realistically reflects the mechanism of image capture. The true image is a combination of both smooth regions, texture, and fine scale structure. Blurring is inevitable, and the additive noise term accounts for other distortions (e.g., threshold variation in the CMOS or CCD light sensors, edge effects between pixels). The reconstruction error (the pixel-wise difference between the true image and the estimated image) is especially useful since it permits a quantitative basis for comparing the quality of demixing algorithms. When two algorithms make comparable assumptions about smoothness and the blur operator, then any structure that persists after demixing appears in the reconstruction error. Its mean squared error, or the eigenvalues of the estimated covariance matrix, enable one to determine which algorithm has successfully extracted more signal.

Despite the ability of the Euler-Lagrange equations associated with the variational model to achieve advanced performance in image analysis, numerical solution of this high order PDE is difficult. Following [22, 34], a numerical augmented Lagrangian method (ALM) is applied to split the directional mean curvature norm into several - and -norms, which are solvable by iterative shrinkage/thresholding (IST) algorithms [36, 37, 38]. Wu et al. [39] proved the equivalence between ALM, dual methods (e.g., Chambolle’s projection [40]), and the splitting Bregman method [8] for solving the convex optimization. (Note that ALM or the splitting Bregman method can be employed in the Douglas-Rachford splitting scheme [41].)

We focus on understanding the advantages of the fourth-order PDE approach to reconstruct a smooth image with sharp edges, and make two contributions:

  1. we provide a general solution to a mathematical model for textured image deconvolution and decomposition with four meaningful components, and

  2. we find a link between functional analysis and multiscale sampling theory in harmonic analysis and filter banks.

This employs a novel Directional Mean Curvature Demixing (DMCD) method for textured images corrupted by i.i.d. noise with a given blur kernel.

To develop these ideas, section 2 introduces the mathematical framework, and section 3 describes the use of directional mean curvature for demixing and the algorithm needed to fit the model. Section 4 makes a quantitative comparison between the demixing results from the proposed algorithm and results from a standard alternative algorithm. Section 5 describes the new model’s correspondence with multiscale harmonic analysis. Section 6 summarizes our conclusions. For readability, mathematical proofs are provided in the Appendix and we also refer the reader to [42, 29, 43, 44] for more detailed notation and literature review.

2 Notation and Mathematical Preliminaries

In this section we define the image and present the directional forward/backward difference operators. We also specify the directional gradient and divergence operators, the directional Laplacian operator, and the discrete directional -norm. Using the curvelet transform and pointwise shrinkage operators, these enable the mathematical derivation of the new demixing algorithm.

Given a discrete grayscale image of size , with the lattice

let be the Euclidean space whose dimension is given by the size of the lattice ; i.e.,

. The 2D discrete Fourier transform

acting on is

where the Fourier coordinates are defined on the lattice as

We let denote the discrete coordinates of the Fourier transform.

Given the matrix

the directional forward and backward difference operators (with periodic boundary condition and direction ) in a matrix form, i.e. , are

The transposed matrices of and are and , respectively. Their adjoint operators are .

The directional gradient and divergence operators are, respectively,

Note that the adjoint operator of is , i.e.,

The directional Laplacian operator is

Remind that given , the impulse responses of directional derivative and directional Laplacian operators in a continuous setting are

with are continuous version of which is numerically used instead.

Extending [45, 46, 47], and due to the discrete nature of images, Thai and Gottschlich [29] defines the discrete directional -norm with in the anisotropic version as


As stated in [45] and [46, p. 87], the space of bounded variation is suitable for the piecewise smooth image component , the G-space is suitable for the texture component , and the dual Besov space is suitable for the noise component , where

Since natural images are better described in terms of multi-scale and multi-direction, we apply the curvelet transform [10, 9, 12] instead of the wavelet transform in the dual Besov space, see [44, 29]. Motivated by the Dantzig selector [48], Aujol and Chambolle [46] and Thai and Gottschlich [29], find that the residual is better captured by (bounded by a constant ) in the curvelet domain, which can be represented by an indicator function on a feasible convex set as


This measure controls the maximum curvelet coefficient of the residual . Due to the curvlet transform, no assumption on the kind of noise is needed (e.g. Gaussian, Laplacian or weakly correlated noise), see [43, 42]. Following [49], Thai and Gottschlich [44] proposes a threshold based on extreme value theory for (2). Note that if

is normally distributed, the random variable

has the Gumbel distribution (since the curvelet of a Gaussian process is weakly correlated), see [49]. In general, for correlated noise, its distribution is a max-stable process [50, 51].

It is known that the solution of minimization is a shrinkage operator, see [36, 37, 52]. It can be defined in a matrix form as

with the point-wise operators. For example, a multiply pointwise operator of functions is . The discrete convolution of two functions is defined in a matrix form as

For simplicity, given , we denote

Given a curvelet transform and its inverse version [10] and a function , the curvelet shrinkage thresholding operator is defined as

Note that one can easily apply other kinds of harmonic analysis than the curvelet, e.g., the shearlet [53], steerable wavelet [15, 54, 55], contourlet [14, 56] and dual-tree complex wavelet [57]. Given a function , its time reversed function is defined as

For more mathematical background and notation, we refer the reader to [52, 44, 29].

3 Directional mean curvature for image demixing (DMCD)

We assume that the original image , which consists of piecewise smooth regions , texture and fine scale structure , is blurred by the operator and corrupted by noise as

Instead of the total variation norm [3], and following the higher order approaches [58, 31, 6, 22, 59], we propose a discrete directional mean curvature (DMC) norm to reconstruct the piecewise-smooth image as

with the matrix of ones (size ). According to [59], a discrete DMC is rewritten as

As in the DG3PD model [29], and following the image generation mechanism shown in Fig. 1, we define the discrete DMCD model for image deconvolution and decomposition problem as


Following [44, 47, 46, 29, 60, 45, 9], the discrete directional -norm measures texture in several directions, and the dual of a generalized Besov space in the curvelet domain captures the residual structure and noise . Note that the

-norm of the curvelet transform is a good measure for fine scale oscillating patterns (i.e., residual structure and noise), which can be either independent or “weakly” correlated and need not follow a Gaussian distribution. The

-norm for texture is handled by the approach of Vese and Osher [47]. According to Meyer [45] (or [61]), the oscillating components do not have small or -norm.

From the definition of the directional -norm (1) and the curvelet space for noise measurement (2), the constrained minimization (3) is rewritten as

Figure 1: Visualization of the innovation model for DMCD. The inverse operator does not exist in general, but the reconstructed image is obtained by the proposed minimization (3).

Similar to [29, 22, 59], in order to solve (3) we introduce five new variables where

with the indicator function on its feasible convex set

The augmented Lagrangian method (ALM) is applied to (3) by introducing Lagrange multipliers and positive parameters as


and the Lagrange function is

Due to multi-variable minimization, the numerical solution of (5) is obtained by applying the alternating directional method of multipliers through iteration to find

Given the initialization as , we solve the following ten subproblems and then update the seven Lagrange multipliers in each iteration.

The “-problem”: Fix and then solve


Given , we denote the discrete Fourier transforms of and by and , respectively. The minimizer of (6) is solved in the Fourier domain as



The “-problem”: Fix and then solve


Let denote the matrix-valued Dirac delta function evaluated at . Then a solution of (7) is


Using the inverse cumulative distribution function to get the quantile corresponding to probability

, one can choose an adaptive at each iteration as

The “-problem”: Fix and then solve


The solution of the minimization (8) is obtained by the shrinkage operator as


The “-problem”: Fix and solve


The minimizer of problem (10) is obtained by the shrinkage operator as

The “-problem”: Fix and solve


Similar to the “-problem”, we denote and as the discrete Fourier transforms of and with , respectively. The solution of this “-problem” (11) for each separable problem is


Note that is similar to the auto-correlation function in the Riesz basis [62].

The “-problem”: Fix and solve


Due to its separability, we consider the problem at and the solution of (12) is

The “