Networks for Nonlinear Diffusion

Networks for Nonlinear Diffusion Problems in Imaging

S. Arridge and A. Hauptmann
July 19, 2019
Abstract.

A multitude of imaging and vision tasks have seen recently a major transformation by deep learning methods and in particular by the application of convolutional neural networks. These methods achieve impressive results, even for applications where it is not apparent that convolutions are suited to capture the underlying physics.

In this work we develop a network architecture based on nonlinear diffusion processes, named DiffNet. By design, we obtain a nonlinear network architecture that is well suited for diffusion related problems in imaging. Furthermore, the performed updates are explicit, by which we obtain better interpretability and generalisability compared to classical convolutional neural network architectures. The performance of DiffNet tested on the inverse problem of nonlinear diffusion with the Perona-Malik filter on the STL-10 image dataset. We obtain competitive results to the established U-Net architecture, with a fraction of parameters and necessary training data.

S. Arridge and A. Hauptmann are with the Department of Computer Science; University College London, London, United Kingdom

1. Introduction

We are currently undergoing a paradigm shift in imaging and vision tasks from classical analytic to learning and data based methods. In particular by deep learning and the application of convolutional neural networks (CNN). Whereas highly superior results are obtained, interpretability and analysis of the involved processes is a challenging and ongoing task [18, 22].

In a general setting, Deep Learning proposes to develop a non-linear mapping between elements of two spaces , (which may be the same) parametrised by a finite set of parameters , which need to be learned:

(1.1)

This learning based approach is in marked contrast to classical methods with a physical interpretation of the process (1.1) for both computational modelling of data given a physical model (which we call a forward problem), and the estimation of parameters of a physical model from measured, usually noisy, data (which we call an inverse problem). Several recent papers have discussed the application of Deep Learning in both forward [21, 27, 32, 33, 36, 37] and inverse [19, 20, 25, 38, 39] problems. Several questions become essential when applying learned models to both forward and inverse problems including:

  • how and to what extent can learned models replace physical models?

  • how do learned models depend on training protocols and how well do they generalise?

  • what are appropriate architectures for the learned models, what is the size of the parameter set that needs to be learned, and how can these be interpreted?

In this paper we aim to answer some of these questions, by taking a few steps back and looking at an analytic motivation for network architectures. Here we consider in particular mappings between images in dimension , i.e. , for which there exist several widely used linear and nonlinear mappings defined by differential and/or integral operators. For example a general linear integral transform such as

(1.2)

includes stationary convolution as a special case if the kernel is translation invariant, i.e. . If the kernel depends on the image, i.e. then (1.2) becomes nonlinear. Alternatively, the forward model may be modelled as the end-point of an evolution process which becomes nonlinear if the kernel depends on the image state ; see eq.(2.7) below for a specific example.

Furthermore, a widely studied class of image mappings is characterised by defining the evolution through a partial differential equation (PDE) where flow is generated by the local structure of [23, 30, 34], such that

(1.3)

These problems in general do not admit an explicit integral transform representation and are solved instead by numerical techniques, but since they depend on a physical model, described by the underlying PDE, they can be analysed and understood thoroughly [35]. This also includes a statistical interpretation [17] as hyperpriors in a Bayesian setting [4]. Furthermore, models like (1.3) can be used to develop priors for nonlinear inverse problems, such as optical tomography [7, 13] and electrical impedance tomography [11].

Motivated by the success of these analytic methods to imaging problems in the past, we propose to combine physical models with data driven methods to formulate network architectures for solving both forward and inverse problems that take the underlying physics into account. We limit ourselves to the case where the physical model is of diffusion type, although more general models could be considered in the future. The leading incentive is given by the observation that the underlying processes in a neural network do not need to be limited to convolutions.

Similar ideas of combining partial differential equations with deep learning have been considered earlier. For instance learning coefficients of a PDE via optimal control [24], as well as deriving CNN architectures motivated by diffusion processes [5], deriving stable architectures by drawing connections to ordinary differential equations [9] and constraining CNNs [29] by the interpretation as a partial differential equation. Another interpretation of our approach can be seen as introducing the imaging model into the network architecture; such approaches have led to a major improvement in reconstruction quality for tomographic problems [2, 12, 15].

This paper is structured as following. In section 2 we review some theoretical aspects of diffusion processes for imaging and the inversion based on theory of partial differential equations and differential operators. We formulate the underlying conjecture for our network architecture, that the diffusion process can be inverted by a set of local non-stationary filters. In the following we introduce the notion of continuum networks in section 3 and formally define the underlying layer operator needed to formulate network architectures in a continuum setting. We draw connections to the established convolutional neural networks in our continuum setting. We then proceed to define the proposed layer operators for diffusion networks in section 4 and derive an implementable architecture by discretising the involved differential operator. In particular we derive a network architecture that is capable of reproducing inverse filtering with regularisation for the inversion of nonlinear diffusion processes. We examine the reconstruction quality of the proposed DiffNet in the following section 5 for an illustrative example of deconvolution and the challenging inverse problem of inverting nonlinear diffusion with the Perona-Malik filter. We achieve results that are competitive to popular CNN architectures with a fraction of the amount of parameters and training data. Furthermore, all computed components that are involved in the update process are interpretable and can be analysed empirically. In section 6 we examine the generalisablity of the proposed network with respect to necessary training data. Additionally, we empirically analyse the obtained filters and test our underlying conjecture. Section 7 presents some conclusions and further ideas.

2. Diffusion and flow processes for imaging

In the following we want to explore the possibility to include a nonlinear process as an underlying model for network architectures. Specifically, the motivation for this study is given by diffusion processes that have been widely used in imaging and vision. Let us consider here the general diffusion process in ; then on a fixed time interval with some diffusivity , to be defined later, we have

(2.1)
Remark 2.1.

When considering bounded domains we will augment (2.1) with boundary conditions on . We return to this point in section 4.

In the following we denote the spatial derivative by

(2.2)

Let us first consider the isotropic diffusion case, then the differential operator becomes the spatial Laplacian . In this case, the solution of (2.1) at time is given by convolution with a Green’s function

(2.3)

where in dimension and we recall that the convolution of two functions is defined by

(2.4)
Definition 2.2.

The Green’s operator is defined with the Green’s function as its kernel

(2.5)

by which we have

In the general case for an anisotropic diffusion flow (ADF) we are interested in a scalar diffusivity that depends on itself, i.e.

(2.6)

This is now an example of a non-linear evolution

(2.7)

where is now a non-stationary, non-linear and time-dependent kernel. In general there is no explicit expression for and numerical methods are required for the solution of (2.6).

Remark 2.3.

Not considered here, but a possible extension to (2.6) is where is a tensor, which, for takes the form

Furthermore, extensions exist for the case where is vector or tensor-valued. We do not consider these cases here, see [34] for an overview.

2.1. Forward Solvers

First of all, let us establish a process between two states of the function . Integrating over time from to yields

Note, that the left hand side can be expressed as and we denote the right hand side by an integral operator , such that

(2.8)

In the following we denote the solution of (2.1) at time instances as . Then we can establish a relation between two time instances of by

(2.9)

where denotes the identity and .

Since we cannot compute by (2.9) without the explicit knowledge of the (possibly time-dependent) diffusivity , it is helpful to rather consider a fixed diffusivity at each time instance , or in the nonlinear case; then by using the differential operator (2.2) we have an approximation of (2.8) by

We can now solve (2.6) approximately by iterating for time steps using either an explicit scheme

(2.10)

or an implicit scheme

(2.11)

Whereas (2.10) is stable only if CFL conditions are satisfied and (2.11) is unconditionally stable, they are both only accurate for sufficiently small steps . In fact, by the Neumann series, the schemes are equivalent in the limit

(2.12)

and coincide with the integral formulation of (2.9).

It’s also useful to look at the Green’s functions solutions.

Lemma 2.4.

Consider the isotropic case with . Then we may write, with time steps ,

(2.13)
Proof.

Take the Fourier Transform111Note the definition of Fourier Transform is chosen to give the correct normalisation so that and use the convolution theorem to give

Let us also note that in Fourier domain, by Taylor series expansion we have

and therefore in the spatial domain the finite difference step and the Gaussian convolution step are the same

(2.14)

2.2. Inverse Filtering

Let us now consider the inverse problem of reversing the diffusion process. That is we have and aim to recover the initial condition . This is a typical ill-posed problem as we discuss in the following.

2.2.1. Isotropic case

As the forward problem is represented as convolution in the spatial domain, the inverse mapping is a (stationary) deconvolution. We remind that , then the inversion is formally given by division in the Fourier domain as

(2.15)

However we note:

  • The factor is unbounded and hence the equivalent convolution kernel in the spatial domain does not exist.

  • (2.15) is unstable in the presence of even a small amount of additive noise and hence it has to be regularised in practice.

Nevertheless, let’s consider formally with that by Taylor series we get

Motivated by this we define an operator for the inversion process

(2.16)

Clearly coincides with the inverse of the implicit update in (2.11), and

(2.17)

is an estimate for the deconvolution problem which (in the absence of noise) is correct in the limit

(2.18)

2.2.2. Anisotropic case

In this case the diffused function is given by (2.7). Following Lemma 2.4 we may put

(2.19)

and we also have

(2.20)
Conjecture 2.5.

There exists a set of local (non-stationary) filters where

(2.21)

and where has only local support and such that

(2.22)
Remark 2.6 (Unsharp Masking).

We recall that a simple method for "deconvolution" is called Unsharp Masking which is usually considered as

for some blur value and sufficiently small . By similar methods as above, we find

We may choose to interpret the operators as a kind of "non-stationary unsharp masking".

2.3. Discretisation

We introduce the definition of a sparse matrix operator representing local non-stationary convolution

Definition 2.7.

is called a Sparse Sub-Diagonal (SSD) matrix if its non-zero entries are all on sub-diagonals corresponding to the local neighbourhood of pixels on its diagonal.

Furthermore we are going to consider that a class of SSD matrices with learned parameters can be decomposed as where is smoothing and is zero-mean; i.e. has one zero eigenvalue such that its application to a constant image gives a zero valued image

In the following we restrict ourselves to the typical 4-connected neighbourhood of pixels in dimension . For the numerical implementation we have the Laplacian stencil

from which we have that is zero-mean. Similarly we will have for the numerical approximation of the matrix operator

Further we conjecture that in the numerical setting is approximated by the sum of identity plus a SSD matrix operator as

(2.23)

In the presence of noise, the inverse operator can be explicitly regularised by addition of a smoothing operation

(2.24)

Whereas in classical approaches to inverse filtering the regularisation operator would be defined a priori, the approach in this paper is to learn the operator and interpret it as the sum of a differentiating operator and a (learned) regulariser . This is discussed further in section 4 below

3. Continuum Networks

Motivated by the previous section, we aim to build network architectures based on diffusion processes. We first discuss the notion of (neural) networks in a continuum setting for which we introduce the concept of a continuum network as a mapping between function spaces. That is, given a function on a bounded domain with , we are interested in finding a non-linear parametrised operator acting on the function . We will consider in the following the case ; extensions to other spaces depend on the involved operations and will be the subject of future studies.

We will proceed by defining the essential building blocks of a continuum network and thence to discuss specific choices to obtain a continuum version of the most common convolutional neural networks. Based on this we will then introduce our proposed architecture as a diffusion network in the next chapter.

3.1. Formulating a continuum network

The essential building blocks of a deep neural network are obviously the several layers of neurons, but since these have a specific notion in classical neural networks, see for instance [31], we will not use the term of neurons to avoid confusion. We rather introduce the concept of layers and channels as the building blocks of a continuum network. In this construction, each layer consists of a set of functions on a product space and each function represents a channel.

Definition 3.1 (Layer and channels).

For , let be a set of functions for , . Then we call: the layer with channels and corresponding index set .

The continuum network is then built by defining a relation or operation between layers. In the most general sense we define the concept of a layer operator for this task.

Definition 3.2 (Layer operator).

Given two layers and , , with channel index set , respectively, we call the mapping with

a layer operator. If the layer operator depends on a set of parameters , then we write .

We note that for simplicity, we will not index the set of parameters, i.e generally stands for both all involved parameters of each layer separately, or the whole network. The classical structure for layer operators follows the principle of affine linear transformations followed by a nonlinear operation. Ideally the affine linear transformation should be parameterisable by a few parameters, whereas the nonlinear operation is often fixed and acts pointwise. A popular choice is the maximum operator also called the “Rectified Linear Unit”:

The continuum network is then given by the composition of all involved layer functions. For example in monochromatic imaging applications we typically have an input image and a desired output with several layer functions inbetween, that perform a specific task such as denoising or sharpening. In this case the input and output consists of one channel, i.e. ; consequently for colour images (in RGB) we have .

3.2. Continuum convolutional networks

Let us now proceed to discuss a specific choice for the layer operator, namely convolutions. Such that we will obtain a continuum version of the widely used convolutional neural networks, which we will call here a continuum convolutional network, to avoid confusion with the established convolutional neural networks (CNN). We note that similar ideas have been addressed as well in [2].

Let us further consider linearly ordered network architectures, that means each layer operator maps between consecutive layers. The essential layer operator for a continuum convolutional network is then given by the following definition.

Definition 3.3 (Convolutional layer operator).

Given two layers and with channel index set , respectively. Let and , with compact support in , be the layer operator’s parameters for all . We call the convolutional layer operator for layer , if for each output channel

(3.1)

with a point-wise nonlinear operator .

If the layer operator does not include a nonlinearity, we write . Now we can introduce the simplest convolutional network architecture by applying convolutional layer operators consecutively.

Definition 3.4 (-layer Continuum Convolutional Network).

Let , then we call the composition of convolutional layer operator, denoted by , a K-layer Continuum Convolutional Network, such that

(3.2)

In the following we will also refer to a -layer CNN as the practical implementation of a -layer continuum convolutional network. A popular network architecture that extends this simple idea is given by a resdiual network (ResNet) [16], that is based on the repetition of a 2-layer CNN with a residual connection, that consists of addition. That is, the network learns a series of additive updates to the input. The underlying structure in ResNet is the repeated application of the following residual block, given by

(3.3)

Note that the second convolutional layer does not include a nonlinearity. Furthermore, it is necessary that , but typically it is often chosen such that . The full continuum ResNet architecture can then be summarized as follows. Let , then the composition of residual blocks, denoted by , defines a -block continuum ResNet

(3.4)

Note that the two additional convolutional layers in the beginning and end are necessary to raise the cardinality of the input/output layer to the cardinality needed in the residual blocks. A complete K-block ResNet then consists of layers. Note that in the original work [16], the network was primarily designed for an image classification task rather than image-to-image mapping that we consider here.

4. DiffNet: discretisation and implementation

Next we want to establish a layer operator based on the diffusion processes discussed in chapter 2. This means that we now interpret the layers of the continuum network as time-states of the function , where is a solution of the diffusion equation (2.1). In the following we assume single channel networks, i.e. for all layers. Then we can associate each layer with the solution such that . To build a network architecture based on the continuum setting, we introduce the layer operator versions of (2.10), and (2.21):

Definition 4.1 (Diffusion and filtering layer operator).

Given two layers and , such that and , then a diffusion layer operator , with parameters , is given by

(4.1)

Similarly, an inverse filtering layer operator, with parameters is given by

(4.2)

Note that this formulation includes a learnable time-step and hence the time instances that each layer represents changes. That also means that a stable step size is implicitly learned, if there are enough layers. In the following we discuss a few options on the implementation of the above layer operator, depending on the type of diffusivity.

Remark 4.2.

The assumption of a single channel network, i.e. for all , can be relaxed easily. Either by assuming for some and all layers, or by introducing a channel mixing as in the convolutional operator (3.1).

4.1. Discretisation of a continuum network

Let us briefly discuss some aspects on the discretisation of a continuum network; let us first start with affine linear networks, such as the convolutional networks discussed in section 3.2. Rather than discussing the computational implementation of a CNN, (see e.g. the comprehensive description in [8]), we concentrate instead on an algebraic matrix-vector formulation that serves our purposes.

For simplicity we concentrate on the two-dimensional case here. Let us then assume that the functions in each layer are represented as a square -by- image and we denote the vectorised form as . Then any linear operation on can be represented by some matrix ; in particular we can represent convolutions as a matrix.

We now proceed by rewriting the convolutional operation (3.1) in matrix-vector notation. Given two layers and with channel index set , respectively, then we can represent the input layer as vector and similarly the output layer . Let represent the sum of convolutions in (3.1), then we can write the layer operator in the discrete setting as matrix-vector operation by

where denotes the identity matrix of suitable size. Furthermore, following the above notation, we can introduce a stacked matrix consisting of all and a matrix by stacking the diagonal matrices . Then we can represent the whole convolutional layer in the discrete setting as

(4.3)

Now the parameters of each layer are contained in the two matrices and .

4.2. Learned Forward and Inverse Operators

Let us now discuss a similar construction for the diffusion layers. For the implementation of the diffusion network, we consider the explicit formulation in (2.10) with the differential operator approximated by the stencil (including the time-step )

(4.4)

We use zero Neumann boundary conditions on the domain boundary

(4.5)

Then we can represent (4.1) by

(4.6)

The basis of learning a diffusion network is now given as estimating the diagonals of and the time-step . This can be done either explicitly as for the CNN or indirectly by an estimator network, as we will discuss next.

4.3. Formulating DiffNet

Let us first note, that if we construct our network by strictly following the update in (4.6), we restrict ourselves to the forward diffusion. To generalise to the inverse problem we consider

(4.7)

Additionally, there are two fundamentally different cases for the diffusivity we need to consider before formulating a network architecture to capture the underlying behaviour. These two cases are

  • Linear diffusion; spatially varying and possible time dependence, .

  • Nonlinear diffusion; diffusivity depending on the solution , .

In the first case, we could simply to try learn the diffusivity explicitly, to reproduce the diffusion process. In the second case, this is not possible and hence an estimation of the diffusivity needs to be performed separately in each time step from the image itself, before the diffusion step can be performed. This leads to two conceptually different network architectures.

Linear Diffusion Network
Figure 1. Illustration for a linear 3-layer Diffusion network. In this case we learn the filters as the diffusivity for each layer explicitly.

The linear case i.) corresponds to the diffusion layer operator (4.1) and is aimed to reproduce a linear diffusion process with fixed diffusivity. Thus, learning the mean-free filter suffices to capture the physics. The resulting network architecture is outlined in Figure 1. Here, the learned filters can be directly interpreted as the diffusivity of layer and are then applied to to produce .

Nonlinear Diffusion Network (DiffNet) k-layer CNNk-layer CNNk-layer CNN
Figure 2. Illustration for a nonlinear 3-layer Diffusion network. Here the filters are implicitly estimated by a small -layer CNN and then applied to the image in the filtering layer.

In the nonlinear case ii.) we follow the same update structure, but now the filters are not learned explicitly, they are rather estimated from the input layer itself, as illustrated in Figure 2. Furthermore, since this architecture is designed for inversion of the nonlinear diffusion process, we employ the generalised stencil . Then, given layer , the filters are estimated by a small CNN from , which are then applied following an explicit update as in (4.6) to produce . Note that the diagonals in the update matrix are produced by the estimation CNN. We will refer to this nonlinear filtering architecture as the DiffNet under consideration in the rest of this study.

4.3.1. Implementation

The essential part for the implementation of a diffusion network is to perform the update (4.6) with either or . For computational reasons it is not practical to build the sparse diagonal matrix and evaluate (4.6), we rather represent the filters and as an -image and apply the filters as pointwise matrix-matrix multiplication to a shifted and cropped image, according to the position in the stencil. This way, the zero Neumann boundary condition (4.5) is also automatically incorporated.

For the linear diffusion network, we would need to learn the parameter set , consisting of filters and time steps, explicitly. This has the advantage of learning a global operation on the image where all parameters are interpretable, but it comes with a few disadvantages. First of all, in this form we are limited to linear diffusion processes and a fixed image size. Furthermore, the parameters grow with the image size, i.e. for an image of size we need parameters per layer. Thus, applications may be limited.

For the nonlinear architecture of DiffNet, where the filters depend on the image at each time step, we introduced an additional estimator consisting of a -layer CNN. This CNN gets an image, given as layer , as input and estimates the filters . The architecture for this -layer CNN as estimator is chosen to be rather simplistic, as illustrated in Figure 3. The input consists of one channel, which is processed by convolutional layers with 32 channels and a ReLU nonlinearity, followed by the last layer without nonlinearity and 5 channels for each filter, which are represented as matrices of the same size as the input . In particular for the chosen filter size of we have exactly convolutional parameters and biases per diffusion layer. That is for a 5-layer CNN 29.509 parameters independent of image size.

K-layer CNN for filter estimation-layers32323232convReLU(conv)
Figure 3. Architecture of the -layer CNN used as diffusivity estimator in the nonlinear diffusion network (DiffNet).

5. Computational experiments

In the following we will examine the reconstruction capabilities of the proposed DiffNet. The experiments are divided into a simple case of deconvolution, where we can examine the learned features, and a more challenging problem of recovering an image from its nonlinear diffused and noise corrupted version.

5.1. Deconvolution with DiffNet

We first examine a simple deconvolution experiment to determine what features the DiffNet learns in an inverse problem. For this task we will only consider deconvolution without noise.

Ground truthKernel Convolved imageNetwork output
Figure 4. Illustration of the deconvolution problem for a simple ball. Top row shows the image space and the bottom row shows the corresponding absolute value of the Fourier coefficients. All images are plotted on their own scale.

The forward problem is given by (2.1) with zero Neumann boundary condition (4.5) and constant diffusivity . For the experiment we choose , which results in a small uniform blurring, as shown in Figure 4. We remind that for the isotropic diffusion, the forward model is equivalent to convolution in space with the kernel , see (2.3). As it is also illustrated in Figure 4, convolution in the spatial domain is equivalent to multiplication in Fourier domain. In particular, high frequencies get damped and the convolved image is dominated by low frequencies. Hence, for the reconstruction task without noise, we essentially need to recover the high frequencies.

-layer CNN-layer CNN-layer CNNReLUImage spaceFourier space
Figure 5. Illustration of the deconvolution process with three layers of DiffNet. The left column shows the processed image and intermediate steps. The right column shows the corresponding absolute value of Fourier coefficients. All images are plotted on their own scale.

The training and test data for DiffNet consists of simple disks of varying radius and contrast. The training set consists of 1024 samples and the test set of an additional 128, each of size . The network architecture is chosen following the schematic in Figure 2, with 3 diffusion layers and a final projection to the positive numbers by a ReLU layer. The filter estimator is given by a 4-layer CNN, as described in section 4.3.1. All networks were implemented in Python with TensorFlow [1].

The input to the network is given by the convolved image without noise and we have minimised the -loss of the output to the ground-truth image. The optimisation is performed for about 1000 epochs in batches of 16 with the Adam algorithm and initial learning rate of and a gradual decrease to . Training on a single Nvidia Titan Xp GPU takes about 24 minutes. The final training and test error is both at a PSNR of 86.24, which corresponds to a relative -error of . We remind that this experiment was performed without noise.

The result of the network and intermediate updates for one example from the test data are illustrated in Figure 5. We also show the filters computed as the output of the trained CNN in each layer, . The output of the last diffusion layer is additionally processed by a ReLU layer to enforce positivity in the final result. It can be seen that the network gradually reintroduces the high frequencies in the Fourier domain; especially the last layer mainly reintroduces the high frequencies to the reconstruction. It is interesting to see, that the learned filters follow indeed the convention that the central filter is of different sign than the directional filters. This enforces the assumption that the filter consists of a mean free part and a regularising part, which should be small in this case, since we do not have any noise in the data. Lastly, we note that the final layer, before projection to the positive numbers, has a clearly negative part around the target, which will be cut off resulting in a sharp reconstruction of the ball.

5.2. Nonlinear diffusion

Let us now consider the nonlinear diffusion process with the Perona-Malik filter function [26] for (2.1) with zero Neumann boundary condition (4.5). In this model the diffusivity is given as a function of the gradient

(5.1)

with contrast parameter . We mainly concentrate here on the inverse problem of restoring an image that has been diffused with the Perona-Malik filter and contaminated by noise.

For the experiments we have used the test data from the STL-10 database [6], which consists of 100,000 RGB images with resolution . These images have been converted to grayscale and divided to 90,000 for training and 10,000 for testing. The obtained images were then diffused for 4 time steps with and . A few sample images from the test data with the result of the diffusion are displayed in Figure 6. The task is then to revert the diffusion process with additional regularisation to deal with noise in the data.

Figure 6. Samples from the test data for learning the inversion of nonlinear diffusion without noise. Mean PSNR for reconstructed test data with DiffNet is: 63.72
Figure 7. Samples from the test data for learning the inversion of nonlinear diffusion with noise. Mean PSNR for reconstructed test data with DiffNet is: 34.21

For all experiments we have used the same network architecture of DiffNet using the architecture illustrated in Figure 2. By performing initial tests on the inversion without noise, we have found that 5 diffusion layers with a 4-layer CNN, following the architecture in 3, gave the best trade-off between reconstruction quality and network size. Increasing the amount of either layers, led to minimal increase in performance. Additionally, we have used a ReLU layer at the end to enforce non-negativity of the output, similarly to the last experiment. We emphasise that this architecture was used for all experiments and hence some improvements for the high noise cases might be expected with more layers. All networks were trained for 18 epochs, with a batch size of 16, and -loss. For the optimisation we have used the Adam algorithm with initial learning rate of and a gradual decrease to . Training on a single Nvidia Titan Xp GPU takes about 75 minutes.

As benchmark, we have performed the same experiments with a widely used network architecture known as U-Net [28]. This architecture has been widely applied in inverse problems [3, 10, 19, 20] and is mainly used to post-process initial directly reconstructed images from undersampled or noisy data. For instance by filtered back-projection in X-ray tomography or the inverse Fourier transform in magnetic resonance imaging [14]. The network architecture we are using follows the publication [19] and differs from the original mainly by a residual connection at the end. That means the network is trained to remove noise and undersampling artefacts from the initial reconstruction. In our context, the network needs to learn how to remove noise and reintroduce edges. For training we have followed a similar protocol as for DiffNet. The only difference is that we started with an initial learning rate of with a gradual decrease to . Training of U-Net takes about 3 hours.

OriginalDiffused (No Noise)ReconstructedResidual
Figure 8. Comparison of reconstruction quality for reconstruction from nonlinear diffused image without noise. Both networks are trained on the full set of 90,000 images. PNSR: DiffNet 65.34, U-Net 61.08
OriginalDiffused ( Noise)ReconstructedResidual
Figure 9. Comparison of reconstruction quality for reconstruction from 1% noise contaminated nonlinear diffused image. Both networks are trained on the full set of 90,000 images. PNSR: DiffNet 34.96, U-Net 35.27

The reconstruction results, for some samples of the test data, with DiffNet can be seen in Figure 6 for the case without noise and in Figure 7 for 1% noise on the diffused images. A comparison of the reconstruction results with U-Net and DiffNet are shown in Figure 8 for the test without noise and in Figure 8 for 1% noise. Qualitatively, the reconstructed images are very similar, as can be seen in the residual images in the last column. The leftover noise pattern for both networks is concentrated on the fine structures of the ship. Quantitatively, for the noise free experiment DiffNet has an increase of 4dB in PSNR compared to the result of U-Net, 65.34 (DiffNet) compared to 61.08 (U-Net). For the case with 1% noise, the quantitative measures are very similar. Here U-Net has a slightly higher PSNR with 35.27 compared to DiffNet with 34.96. A thorough study of reconstruction quality of both networks follows in the next section as well as some interpretation of the learned features in DiffNet.

6. Discussion

First of all we note that the updates in DiffNet are performed explicitly and that the CNN in the architecture is only used to produce the filters . This means that DiffNet needs to learn a problem specific processing, in contrast to a purely data driven processing in a CNN. Consequently, the amount of necessary learnable parameters is much lower. For instance the 5-layer DiffNet used for inversion of the nonlinear diffusion in section 5.2 has 101,310 parameters, whereas the used U-Net with a filter size of has a total of 34,512,705 parameters; i.e DiffNet uses only of parameters compared to U-Net and hence the learned features can be seen to be much more explicit. In the following we discuss some aspects that arise from this observation, such as generalisability and interpretability.

Inverse Problem (1% noise)Test error vs. training size
Figure 10. Generalisation plot for the forward and inverse problem of nonlinear diffusion and varying noise levels. Test error depending on the amount of training data, for both DiffNet and U-Net.

6.1. Generalisability

To test the generalisation properties of the proposed DiffNet we have performed similar experiments as in section 5.2 for nonlinear diffusion, but with increasing amounts of training data. Under the assumption that DiffNet learns a more explicit update than a classic CNN, we would expect also to require less training data to achieve a good test error. To certify this assumption, we have examined 4 settings of nonlinear diffusion with the Perona-Malik filter: learning the forward model, learning to reconstruct from the diffused image without noise, as well as with 0.1% and 1% noise. We then created training data sets of increasing size from just 10 samples up the full size of 90,000. For all scenarios we have trained DiffNet and U-Net following the training protocol described in 5.2. Additionally, we have aborted the procedure when the networks started to clearly overfit the training data.

Results for the four scenarios are shown in Figure 10. Most notably DiffNet outperforms U-Net clearly for the forward problem and the noise free inversion, by 4dB and 3dB, respectively. For the noisy cases both networks perform very similar for the full training data size of 90,000. The biggest difference overall is that DiffNet achieves its maximum test error already with 500-1,000 samples independent of the noise case, whereas the U-Net test error saturates earlier with higher noise. In conclusion we can say, that for the noisy cases both networks are very comparable in reconstruction quality, but for small amounts of data the explicit nature of DiffNet is clearly superior.

6.2. Interpretation of learned filters

Since all updates are performed explicitly with the output from the filter estimation CNN, we can interpret some of the learned features. For this purpose we show the filters for the ship image from section 5.2 for the three inversion scenarios under consideration. In Figure 11 we show the sum of all learned filters, i.e. . If the network would only learn the mean-free differentiating part, then these images should be zero. That implies that the illustrated filters in Figure 11 can be related to the learned regularisation . Additionally, we also show the diagonal filters in Figure 12.

We would expect that with increasing noise level, the filters will incorporate more smoothing to deal with the noise; this implies that the edges get wider with increasing noise level. This can be nicely observed for the diagonal filters in Figure 12. For the smoothing in Figure 11, we see that the first layer consists of broader details and edges that are refined in the noise free case for increasing layers. In the noisy case the later layers include some smooth features that might depict the regularisation necessary in the inversion procedure. It is generally interesting to observe, that the final layer shows very fine local details, necessary to restore fine details for the final output.

Layer 1Layer 2Layer 3Layer 4Layer 5
Figure 11. Filter updates for different noise levels. Each image is displayed on its own scale.
Layer 1Layer 2Layer 3Layer 4Layer 5
Figure 12. Obtained diagonal filters for different noise levels. Each filter is displayed on its own scale.

Finally, we have computed training data of a wider noise range to examine the regularisation properties of the learned network. For this we have taken the full 90,000 training samples and contaminated the diffused image with noise ranging from to 10% noise. As we conjectured in section 2.3, the learned update filters can be decomposed to a mean-free operation and a smoothing part . This implies that the magnitude of has to increase with higher noise. To examine this conjecture, we have taken (fixed) 32 samples from the reconstructed test data for each noise level and computed an estimate of as the sum , i.e. the deviation from the mean-free part. Furthermore, we interpret the smoothing level as the mean of the absolute value of . The resulting graph of smoothing versus noise level is shown in Figure 13. As we have conjectured, the estimate of increases clearly with the noise level and hence we believe our interpretation of the learned filters as the composition of a mean-free part and a smoothing necessary for ill-posed inverse problems is valid.

Increasing smoothing with noise level
Figure 13. Estimate of the smoothing level for increasing noise in the inverse problem. Computed over a sample of 32 images from the test data.

7. Conclusions

In this paper we have explored the possibility to establish novel network architectures based on physical models other than convolutions; in particular we concentrated here on diffusion processes. As main contributions, we have introduced some non-linear forward mappings, modelled through learning rather than just through PDEs or integral transforms. We have reviewed (regularised) inverse diffusion processes for inverting such maps. In particular, we have conjectured that these inverse diffusion processes can be represented by local non-stationary filters, which can be learned in a network architecture. More specific, these local filters can be represented by a sparse sub-diagonal (SSD) matrix and hence efficiently used in the discrete setting of a neural network. We emphasise that even though we have concentrated this study on a specific structure for these SSD matrices based on diffusion, other (higher order) models can be considered.

We obtain higher interpretability of the network architecture, since the image processing is explicitly performed by the application of the SSD matrices. Consequently, this means that only a fraction of parameters is needed in comparison to classical CNN architectures to obtain similar reconstruction results. We believe that the presented framework and the proposed network architectures can be useful for learning physical models in the context of imaging and inverse problems, especially, where a physical interpretation of the learned features is crucial to establish confidence in the imaging task.

Acknowledgements

We thank Jonas Adler and Sebastian Lunz for valubale discussions and comments. We acknowledge the support of NVIDIA Corporation with one Titan Xp GPU.

References

  • [1] M. Abadi et al., TensorFlow: Large-scale machine learning on heterogeneous systems, 2015. Software available from https://www.tensorflow.org/.
  • [2] J. Adler and O. Öktem, Solving ill-posed inverse problems using iterative deep neural networks, Inverse Problems, 33 (2017), p. 124007.
  • [3] S. Antholzer, M. Haltmeier, and J. Schwab, Deep learning for photoacoustic tomography from sparse data, Inverse Problems in Science and Engineering, (2018), pp. 1–19.
  • [4] D. Calvetti and E. Somersalo, Hypermodels in the Bayesian imaging framework, Inverse Problems, 24 (2008), p. 034013.
  • [5] Y. Chen, W. Yu, and T. Pock, On learning optimized reaction diffusion processes for effective image restoration, in Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 5261–5269.
  • [6] A. Coates, A. Ng, and H. Lee, An analysis of single-layer networks in unsupervised feature learning, in Proceedings of the fourteenth international conference on artificial intelligence and statistics, 2011, pp. 215–223.
  • [7] A. Douiri, M. Schweiger, J. Riley, and S. Arridge, Local diffusion regularization method for optical tomography reconstruction by using robust statistics, Opt. Lett., 30 (2005), pp. 2439–2441.
  • [8] I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning, MIT Press, 2016. http://www.deeplearningbook.org.
  • [9] E. Haber and L. Ruthotto, Stable architectures for deep neural networks, Inverse Problems, 34 (2017), p. 014004.
  • [10] S. J. Hamilton and A. Hauptmann, Deep d-bar: Real time electrical impedance tomography imaging with deep neural networks, IEEE Transactions on Medical Imaging, (2018).
  • [11] S. J. Hamilton, A. Hauptmann, and S. Siltanen, A data-driven edge-preserving D-bar method for electrical impedance tomography, Inverse Problems and Imaging, 8 (2014), pp. 1053–1072.
  • [12] K. Hammernik, T. Klatzer, E. Kobler, M. P. Recht, D. K. Sodickson, T. Pock, and F. Knoll, Learning a variational network for reconstruction of accelerated MRI data, Magnetic resonance in medicine, 79 (2018), pp. 3055–3071.
  • [13] A. Hannukainen, L. Harhanen, N. Hyvönen, and H. Majander, Edge-promoting reconstruction of absorption and diffusivity in optical tomography, Inverse Problems, 32 (2015), p. 015008.
  • [14] A. Hauptmann, S. Arridge, F. Lucka, V. Muthurangu, and J. Steeden, Real-time cardiovascular MR with spatio-temporal artifact suppression using deep learning-proof of concept in congenital heart disease., Magnetic resonance in medicine, (2018).
  • [15] A. Hauptmann, F. Lucka, M. Betcke, N. Huynh, J. Adler, B. Cox, P. Beard, S. Ourselin, and S. Arridge, Model-based learning for accelerated, limited-view 3-d photoacoustic tomography, IEEE transactions on medical imaging, 37 (2018), pp. 1382–1393.
  • [16] K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for image recognition, in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  • [17] T. Helin and M. Lassas, Hierarchical models in statistical inverse problems and the Mumford–Shah functional, Inverse problems, 27 (2010), p. 015008.
  • [18] C. Hofer, R. Kwitt, M. Niethammer, and A. Uhl, Deep learning with topological signatures, in Advances in Neural Information Processing Systems, 2017, pp. 1634–1644.
  • [19] K. H. Jin, M. T. McCann, E. Froustey, and M. Unser, Deep convolutional neural network for inverse problems in imaging, IEEE Transactions on Image Processing, 26 (2017), pp. 4509–4522.
  • [20] E. Kang, J. Min, and J. C. Ye, A deep convolutional neural network using directional wavelets for low-dose X-ray CT reconstruction, Medical Physics, 44 (2017).
  • [21] Y. Khoo, J. Lu, and L. Ying, Solving parametric PDE problems with artificial neural networks, ArXiv, 1707.03351v2 (2017).
  • [22] B. Kim, M. Wattenberg, J. Gilmer, C. Cai, J. Wexler, F. Viegas, et al., Interpretability beyond feature attribution: Quantitative testing with concept activation vectors (tcav), in International Conference on Machine Learning, 2018, pp. 2673–2682.
  • [23] R. Kimmel, Numerical geometry of images: Theory, algorithms, and applications, Springer Science & Business Media, 2003.
  • [24] R. Liu, Z. Lin, W. Zhang, and Z. Su, Learning pdes for image restoration via optimal control, in European Conference on Computer Vision, Springer, 2010, pp. 115–128.
  • [25] T. Meinhardt, M. Moeller, C. Hazirbad, and D. Cremers, Learning proximal operators: Using denoising networks for regularizing inverse imaging problems, in International Conference on Computer Vision, 2017, pp. 1781–1790.
  • [26] P. Perona and J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence, 12 (1990), pp. 629–639.
  • [27] M. Raissi and G. E. Karniadakis, Hidden physics models: Machine learning of nonlinear partial differential equations, ArXiv, 1708.00588v2 (2017).
  • [28] O. Ronneberger, P. Fischer, and T. Brox, U-net: Convolutional networks for biomedical image segmentation, in International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, 2015, pp. 234–241.
  • [29] L. Ruthotto and E. Haber, Deep neural networks motivated by partial differential equations, arXiv preprint arXiv:1804.04272, (2018).
  • [30] G. Sapiro, Geometric partial differential equations and image analysis, Cambridge university press, 2006.
  • [31] S. Shalev-Shwartz and S. Ben-David, Understanding machine learning: From theory to algorithms, Cambridge university press, 2014.
  • [32] J. Sirignano and K. Spiliopoulos, DGM: A deep learning algorithm for solving partial differential equations, ArXiv, 1708.07469v1 (2017).
  • [33] J. Tompson, K. Schlachter, P. Sprechmann, and K. Perlin, Accelerating Eulerian fluid simulation with convolutional networks, ArXiv, 1607.03597v6 (2017).
  • [34] J. Weickert, Anisotropic Diffusion in Image Processing, Teubner Stuttgart, 1998.
  • [35] J. Weickert, B. T. H. Romeny, and M. A. Viergever, Efficient and reliable schemes for nonlinear diffusion filtering, IEEE transactions on image processing, 7 (1998), pp. 398–410.
  • [36] E. Weinan, H. Jiequn, and J. Arnulf, Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations, ArXiv, 1706.04702v1 (2017).
  • [37] Y. Wu, P. Zhang, H. Shen, , and H. Zhai, Visualizing neural network developing perturbation theory, ArXiv, 1802.03930v2 (2018).
  • [38] J. C. Ye, Y. Han, and E. Cha, Deep convolutional framelets: A general deep learning framework for inverse problems, SIAM Journal on Imaging Sciences, 11 (2018), pp. 991–1048.
  • [39] B. Zhu, J. Z. Liu, S. F. Cauley, B. R. Rosen, and M. S. Rosen, Image reconstruction by domain-transform manifold learning, Nature, 555 (2018), pp. 487–489.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
321510
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description