Pyxu_diffops#
Pyxu-diffops
Welcome to Pyxu-diffops, a Pyxu plugin implementing many diffusion operators which are popular in the PDE-based image processing community. Here, we provide a basic introduction to the PDE-based image processing techniques implemented by this plugin.
Diffusion operators#
This plugin provides an interface to deal with PDE-based image processing [Weickert], [Tschumperle]. For simplicity, throughout the documentation we consider a \(2\)-dimensional image \(\mathbf{f} \in \mathbb{R}^{N_{0} \times N_1}\), but higher dimensional signals could be considered. We denote by \(f_i\) the \(i\)-th pixel of \(\mathbf{f},\;i=0,\ldots,N_{tot}-1,\;N_{tot}=N_0N_1\). Furthermore, let \(\boldsymbol{\nabla} \mathbf{f} \in \mathbb{R}^{2 \times N_{0} \times N_1}\) be the image gradient; we denote by \((\boldsymbol{\nabla} \mathbf{f})_i \in \mathbb{R}^{2}\) the local gradient in correspondence of the \(i\)-th pixel.
Let us first consider the simple case of PDE-based image smoothing [Weickert]. Gaussian smoothing, or linear diffusion filtering, can be achieved computing the time evolution of the Partial Differential Equation (PDE)
where \(t\) represents an artificial time and the divergence operator acts on the local gradient. The above diffusion equation corresponds to the well-known heat equation. It achieves linear, isotropic, homogeneous smoothing. In practice, the time evolution can be computed discretizing the PDE in time; if we use the simple explicit Euler method with time-step \(\Delta t=\tau\) , we obtain the scheme
Letting the PDE evolve in time, we obtain increasingly smoother versions of the original image. Eventually, we would reach the steady state \(\frac{\partial\mathbf{f}}{\partial t}=0\), corresponding to a flat image. By early stopping, we can achieve arbitrarily smooth versions of \(\mathbf{f}\).
In Gaussian smoothing, the divergence is applied to the flux \(\boldsymbol{\nabla} \mathbf{f}\), which coincides with the image gradient. However, the flux can be locally modified to promote desired features by designing suitable diffusion tensors \(\mathbf{D}(\mathbf{f})\), leading to the diffusion equation
The diffusion tensor acts on the local gradient, biasing the diffusion process to enhance features of interest. Anisotropic, inhomogeneous smoothing can thus be achieved. Many diffusion tensors have been proposed. Furthermore, generalizations going beyond purely divergence-based diffusion terms have been investigated [Tschumperle].
The Pyxu-diffops module allows to consider diffusion processes that, in their most general form, can be written as
- where
\(\mathbf{D}_{out} = \mathbf{D}_{out}(\mathbf{f})\) is the outer diffusivity for the divergence term,
\(\mathbf{D}_{in} = \mathbf{D}_{in}(\mathbf{f})\) is the diffusion coefficient for the divergence term,
\(\mathbf{b} = \mathbf{b}(\mathbf{f})\) is the balloon force,
\(\mathbf{T}_{out} = \mathbf{T}_{out}(\mathbf{f})\) is the outer diffusivity for the trace term,
\(\mathbf{T}_{in} = \mathbf{T}_{in}(\mathbf{f})\) is the diffusion coefficient for the trace term,
\(\mathbf{H}(\mathbf{f})\) is the Hessian evaluated at \(\mathbf{f}\),
\(\mathbf{w}\) is a vector field assigning a \(2\)-dimensional vector to each pixel,
\(\mathbf{J}_\mathbf{w}\) is the Jacobian of the vector field \(\mathbf{w}\).
Pyxu-diffops offers a collection of diffusion operators, implemented as Pyxu differentiable functionals DiffFunc. These classes
exhibit an apply() and a grad() method.
The right-hand side of the above PDE represents the output of the grad() method applied to the image \(\mathbf{f}\).
To understand why the diffusion operators are implemented as functionals, i.e., why they exhibit an apply() method, we need to need to consider PDE-based restoration rather than simple smoothing.
PDE-based restoration#
In some cases, the right-hand side of the PDE (3) admits a potential, i.e., there exists a functional \(\phi(\mathbf{x})\) such that \(-\nabla\phi(\mathbf{x})\) is the right-hand side of Eq.:eq:grad_flow_generic. In this case,
the algorithm (1) can be seen as a gradient descent algorithm towards the minimizer of the functional \(\phi\). As already remarked, such minimizer will be the steady-state of the PDE (3). Typically, the steady-state of
diffusion operators is a flat solution.
However, non-trivial steady-states are achieved when the diffusion operators are used for PDE-based image restoration. Indeed, let us assume that we observe a noisy version \(\mathbf{y}\) of the true image \(\mathbf{x}\), so that \(\mathbf{y}=\mathbf{x}+\boldsymbol{\varepsilon}\), where \(\boldsymbol{\varepsilon}\) is the noise corrupting the image. If we know the noise distribution (Gaussian/Poisson/…), we can define a likelihood function \(l(\mathbf{y}\vert\mathbf{x})\) describing the likelihood of observing \(\mathbf{y}\) if the ground truth is \(\mathbf{x}\). Typically, the inverse problem of recovering \(\mathbf{x}\) from the noisy \(\mathbf{y}\) is an ill-posed problem. Regularization is a popular way of addressing such ill-posedness: the image restoration problem is formulated as the penalized optimization problem
where \(\phi\) is the so-called regularization functional and \(\lambda\) the regularization parameter determining the amount of imposed regularization. To minimize (4), we can consider the gradient flow
Now, any diffusion operator can be used to replace the term \(-\nabla \phi(\mathbf{x})\) in Eq.:eq:gradient_flow_penalized by the right-hand side of the PDE (3).
Therefore, diffusion operators can be used for PDE-based image restoration combining them with a suitable likelihood term. If the diffusion operators admit a potential \(\phi\), the penalized optimization problem (4) can be explicitly
defined. However, diffusion operators do not necessarily derive from a potential. As a matter of fact, many of them don’t. Nevertheless, they can still be used for PDE-based image smoothing based on Eq.:eq:gradient_flow_penalized, with the only caveat
that the resulting minimizer cannot be interpreted as the minimizer of an energy functional \(l(\mathbf{y}\vert\cdot)+\lambda\phi(\cdot)\).