Deep Graph-Convolutional Image Denoising

TL;DR

This paper introduces a novel neural network architecture using graph convolution layers with dynamically computed graphs to incorporate non-local self-similarity for improved image denoising, achieving state-of-the-art results.

Contribution

It presents a new end-to-end trainable network with graph convolution layers that dynamically model non-local relationships, enhancing denoising performance over existing methods.

Findings

Achieves state-of-the-art denoising results on synthetic and real noise.

Demonstrates the effectiveness of dynamic graph convolution in capturing self-similarity.

Introduces a lightweight Edge-Conditioned Convolution to improve training stability.

Abstract

Non-local self-similarity is well-known to be an effective prior for the image denoising problem. However, little work has been done to incorporate it in convolutional neural networks, which surpass non-local model-based methods despite only exploiting local information. In this paper, we propose a novel end-to-end trainable neural network architecture employing layers based on graph convolution operations, thereby creating neurons with non-local receptive fields. The graph convolution operation generalizes the classic convolution to arbitrary graphs. In this work, the graph is dynamically computed from similarities among the hidden features of the network, so that the powerful representation learning capabilities of the network are exploited to uncover self-similar patterns. We introduce a lightweight Edge-Conditioned Convolution which addresses vanishing gradient and over-parameterization issues of this particular graph convolution. Extensive experiments show state-of-the-art performance with improved qualitative and quantitative results on both synthetic Gaussian noise and real noise.

Tables

Table 1. Table I: PSNR (dB) v. non-local neighborhood size (Set12)

$\sigma$	0-NN	4-NN	8-NN	12-NN	16-NN	20-NN
15	32.91	33.09	33.11	33.13	33.14	33.13
25	30.50	30.70	30.74	30.75	30.78	30.78
50	27.28	27.52	27.58	27.58	27.60	27.59

Table 2. Table II: Edge attention v. ECC + edge attention (8-NN). PSNR (dB).

$\sigma$	Edge attention only	Proposed
25	30.53	30.74

Table 3. Table III: Natural image denoising results. Metrics are PNSR (dB) and SSIM.

Dataset	Noise $\sigma$	BM3D	WNNM	TNRD	DnCNN	N3Net	NLRN	GCDN
Set12	15	32.37 / 0.8952	32.70 / 0.8982	32.50 / 0.8958	32.86 / 0.9031	- / -	33.16 / 0.9070	33.14 / 0.9072
	25	29.97 / 0.8504	30.28 / 0.8557	30.06 / 0.8512	30.44 / 0.8622	30.55 / -	30.80 / 0.8689	30.78 / 0.8687
	50	26.72 / 0.7676	27.05 / 0.7775	26.81 / 0.7680	27.18 / 0.7829	27.43 / -	27.64 / 0.7980	27.60 / 0.7957
BSD68	15	31.07 / 0.8717	31.37 / 0.8766	31.42 / 0.8769	31.73 / 0.8907	- / -	31.88 / 0.8932	31.83 / 0.8933
	25	28.57 / 0.8013	28.83 / 0.8087	28.92 / 0.8093	29.23 / 0.8278	29.30 / -	29.41 / 0.8331	29.35 / 0.8332
	50	25.62 / 0.6864	25.87 / 0.6982	25.97 / 0.6994	26.23 / 0.7189	26.39 / -	26.47 / 0.7298	26.38 / 0.7389
Urban100	15	32.35 / 0.9220	32.97 / 0.9271	31.86 / 0.9031	32.68 / 0.9255	- / -	33.42 / 0.9348	33.47 / 0.9358
	25	29.70 / 0.8777	30.39 / 0.8885	29.25 / 0.8473	29.97 / 0.8797	30.19 / -	30.88 / 0.9003	30.95 / 0.9020
	50	25.95 / 0.7791	26.83 / 0.8047	25.88 / 0.7563	26.28 / 0.7874	26.82 / -	27.40 / 0.8244	27.41 / 0.8160

Table 4. Table IV: Depth map denoising results. Metrics are PNSR (dB) and SSIM.

$\sigma$	Method	aloe	art	baby	cones	dolls	laundry	moebius	reindeer	Average
15	GCDN	40.74 / 0.9873	40.66 / 0.9886	41.64 / 0.9917	39.29 / 0.9832	40.70 / 0.9830	41.97 / 0.9842	42.07 / 0.9877	42.62 / 0.9915	41.21 / 0.9872
	NLRN	40.50 / 0.9844	40.48 / 0.9858	41.76 / 0.9899	39.50 / 0.9814	40.69 / 0.9800	41.96 / 0.9814	42.01 / 0.9848	42.44 / 0.9880	41.17 / 0.9845
	OGLR	40.82 / 0.9801	40.77 / 0.9821	40.90 / 0.9806	39.65 / 0.9774	40.41 / 0.9756	41.32 / 0.9764	41.48 / 0.9793	41.72 / 0.9823	40.88 / 0.9792
25	GCDN	37.12 / 0.9771	37.15 / 0.9788	37.50 / 0.9814	35.88 / 0.9697	37.05 / 0.9705	38.62 / 0.9730	38.39 / 0.9786	38.80 / 0.9836	37.56 / 0.9766
	NLRN	37.08 / 0.9720	37.01 / 0.9734	37.37 / 0.9797	36.09 / 0.9661	37.01 / 0.9646	38.42 / 0.9679	38.33 / 0.9723	38.65 / 0.9786	37.50 / 0.9718
	OGLR	36.67 / 0.9592	36.68 / 0.9649	36.29 / 0.9594	35.51 / 0.9545	36.41 / 0.9541	37.44 / 0.9541	37.17 / 0.9575	37.86 / 0.9655	36.75 / 0.9587
50	GCDN	33.37 / 0.9522	33.18 / 0.9536	32.23 / 0.9468	31.61 / 0.9379	32.37 / 0.9417	34.07 / 0.9526	33.73 / 0.9567	34.35 / 0.9672	33.11 / 0.9511
	NLRN	33.23 / 0.9444	32.86 / 0.9448	32.42 / 0.9534	31.53 / 0.9304	32.40 / 0.9347	34.15 / 0.9459	33.58 / 0.9475	34.37 / 0.9603	33.07 / 0.9452
	OGLR	32.24 / 0.9121	31.92 / 0.9129	31.23 / 0.9027	30.21 / 0.8926	31.44 / 0.8999	32.85 / 0.9051	32.46 / 0.9093	32.99 / 0.9191	31.92 / 0.9067

Table 5. Table V: Real image denoising (SIDD dataset)

	CBM3D	DnCNN	NLRN	GCDN
PSNR	38.73 dB	39.98 dB	41.24 dB	41.48 dB
SSIM	0.9587	0.9605	0.9652	0.9697

Equations

$${\mathbf{H}}_i^{l+1,\mathrm{NL}}$$

$$=\sum _{j\in {\mathcal{S}}_i^l}{\gamma }^{l,j\to i}\frac{{\mathbf{\Theta }}^{l,j\to i}{\mathbf{H}}_j^l}{|{\mathcal{S}}_i^l|},$$

$${\gamma }^{l,j\to i}=\exp \left(-\Vert {\mathbf{d}}^{l,j\to i}{\Vert }_2^2/\delta \right)$$

$${\mathbf{H}}_i^{l+1}=\frac{{\mathbf{H}}_i^{l+1,\mathrm{NL}}+{\mathbf{H}}_i^{l+1,\mathrm{L}}}{2}+{\mathbf{b}}^l,$$

$$f({\mathbf{H}}_i^l,{\mathbf{H}}_j^l)\propto \exp \left({\mathbf{u}}^T{\mathbf{H}}_i^l+{\mathbf{v}}^T{\mathbf{H}}_j^l+c\right),$$

$$f({\mathbf{H}}_i^l,{\mathbf{H}}_j^l)=\exp \left(-\frac{1}{2}({\mathbf{d}}^{l,j\to i}-\mathbf{\mu }{)}^T{\mathbf{\Sigma }}^{-1}({\mathbf{d}}^{l,j\to i}-\mathbf{\mu })\right),$$

$$f({\mathbf{H}}_i^l,{\mathbf{H}}_j^l)=\frac{\exp \left({\mathbf{H}}_i^{lT}{\mathbf{W}}_{\theta }^T{\mathbf{W}}_{\varphi }{\mathbf{H}}_j^l\right)}{{\sum }_{j\in {\mathcal{S}}_i}\exp \left({\mathbf{H}}_i^{lT}{\mathbf{W}}_{\theta }^T{\mathbf{W}}_{\varphi }{\mathbf{H}}_j^l\right)}{\mathbf{W}}_g,$$

$$\mathbf{A}=\mathbf{\Phi }\mathbf{\Lambda }{\mathbf{\Psi }}^T=\sum _s{\lambda }_s{\mathbf{\varphi }}_s{\mathbf{\psi }}_s^T,$$

$${\mathbf{\Theta }}^{l,j\to i}=\sum _{s=1}^r{\kappa }_s^{j\to i}{\mathbf{\theta }}_s^{j\to i,L}{\mathbf{\theta }}_s^{j\to i,R^T},$$

$${\mathbf{H}}_i^{l+1,\mathrm{NL}}$$

$$=\sum _{j\in {\mathcal{S}}_i^l}{\gamma }^{l,j\to i}\frac{{\sum }_{s=1}^r{\kappa }_s^{j\to i}{\mathbf{\theta }}_s^{j\to i,L}{\mathbf{\theta }}_s^{j\to i,R^T}{\mathbf{H}}_j^l}{|{\mathcal{S}}_i^l|},$$

$$\begin{array}{rl}{\mathbf{\theta }}_{s,u}^{j\to i,L}& \sim \mathcal{N}(0,F^l{\sigma }_L^2),\\ {\mathbf{\theta }}_{s,u}^{j\to i,R}& \sim \mathcal{N}(0,F^l{\sigma }_R^2),\\ {\kappa }_s^{j\to i}& \sim \mathcal{N}(0,F^l{\sigma }_{\kappa }^2),\end{array}$$

$${\mathbf{H}}_{i,u}^{l+1,\mathrm{NL}}\sim \mathcal{N}\left(0,\frac{1}{2}r{F}^{l^4}{\sigma }_L^2{\sigma }_R^2{\sigma }_{\kappa }^2\right),$$

$$\mathbf{y}=\mathbf{Ax}+\mathbf{n}$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}}{\mathrm{argmin}}\left[\frac{1}{2}\Vert \mathbf{y}-\mathbf{Ax}{\Vert }_2^2+\frac{\beta }{2}{\mathbf{x}}^T\mathbf{Lx}\right]$$

$${\mathbf{x}}^{(t+1)}$$

$$={\text{prox}}_g\left((\mathbf{I}-\alpha {\mathbf{A}}^T\mathbf{A}){\mathbf{x}}^{(t)}+\alpha {\mathbf{A}}^T\mathbf{y}\right)$$

$${\text{prox}}_g\left(\mathbf{\mu }\right)$$

$${\mathbf{x}}^{(t+1)}={\left(\mathbf{I}+\beta \mathbf{L}\right)}^{-1}\left[(\mathbf{I}-\alpha {\mathbf{A}}^T\mathbf{A}){\mathbf{x}}^{(t)}+\alpha {\mathbf{A}}^T\mathbf{y}\right].$$

$${\mathbf{x}}^n-{\mathbf{\nu }}^{(t+1)}=$$

$${\left(\mathbf{I}+\beta \mathbf{L}\right)}^{-1}\left[\left(\mathbf{I}-\alpha {\mathbf{A}}^T\mathbf{A}\right)\left({\mathbf{x}}^n-{\mathbf{\nu }}^{(t)}\right)+\alpha {\mathbf{A}}^T\mathbf{y}\right].$$

$${\mathbf{\nu }}^{(t+1)}={\left(\mathbf{I}+\beta \mathbf{L}\right)}^{-1}\left[\left(\mathbf{I}-\alpha {\mathbf{A}}^T\mathbf{A}\right){\mathbf{\nu }}^{(t)}+\beta \mathbf{L}{\mathbf{x}}^n\right].$$

$${\mathbf{H}}_i^{l+1,\mathrm{NL}}$$

Full text

I Introduction

Denoising is a staple among image processing problems and its importance cannot be overstated. Despite decades of work and countless methods, it still remains an active research topic because its purpose goes far beyond generating visually pleasing pictures. Denoising is fundamental to enhance the performance of higher-level computer vision tasks such as classification, segmentation or object recognition, and is a building block in the solution to various problems [1, 2, 3, 4]. The recent successes achieved by convolutional neural networks (CNNs) extended to this problem as well and have brought a new generation of learning-based methods that is redefining the state of the art. However, it is important to learn the lessons of past research on the topic and integrate them with the new deep learning techniques. In particular, classic denoising methods such as BM3D [5] showed the importance of exploiting non-local self-similar patterns. However, the convolution operation underpinning all CNNs architectures [6, 7, 8, 9] is unable to capture such patterns because of the locality of the convolution kernels. Only very recently, some works started addressing the integration of non-local information into CNNs [10, 11, 12, 13].

This paper presents a denoising neural network, called GCDN, where the convolution operation is generalized by means of graph convolution, which is used to create layers with hidden neurons having non-local receptive fields that successfully capture self-similar information. Graph convolution is a generalization of the traditional convolution operation when the data are represented as sitting over the vertices of a graph. In this work, every pixel is a vertex and the edges in the graph are dynamically computed from the similarities in the feature space of the hidden layers of the network. This allows us to exploit the powerful representational features of neural networks to discover and use latent self-similarities. With respect to other CNNs integrating non-local information for the denoising task, the proposed approach has several advantages: i) it creates an adaptive receptive field for the pixels in the hidden layers by dynamically computing a nearest-neighbor graph from the latent features; ii) it creates dynamic non-local filters where feature vectors that may be spatially distant but close in a latent vector space are aggregated with weights that depend on the features themselves; iii) the aggregation weights are estimated by a fully-learned operation, implemented as a subnetwork, instead of a predefined parameterized operation, allowing more generality and adaptability. Starting from the Edge-Conditioned Convolution (ECC) definition of graph convolution, we propose several improvements to address stability, over-parameterization and vanishing gradient issues. Finally, we also propose a novel neural network architecture which draws from an analogy with an unrolled regularized optimization method.

A preliminary version of this work appeared in [14]. There are several differences with the work in this paper. The architecture of the network is improved by drawing an analogy with proximal gradient descent methods, and it is significantly deeper. Moreover, we propose several solutions to address the ECC overparameterization and computational issues. Finally, we also present an in-depth analysis of the network behavior and greatly extended experimental results.

This paper is structured as follows. Sec. II provides some background material on graph-convolutional neural networks and state-of-the-art denoising approaches. Sec. III describes the proposed method. Sec. IV analyzes the characteristics of the proposed method and experimentally compares it with state-of-the-art approaches. Finally, Sec. V draws some conclusions.

II Related work

II-A Graph neural networks

Inspired by the overwhelming success of deep neural networks in computer vision, a significant research effort has recently been made in order to develop deep learning methods for data that naturally lie on irregular domains. One case is when the data domain can be structured as a graph and the data are defined as vectors on the nodes of this graph. Extending CNNs from signals with a regular structure, such as images and video, to graph-structured signals is not straightforward, since even simple operations such as shifts are undefined over graphs.

One of the major challenges in this field is defining a convolution-like operation for this kind of data. Convolution has a key role in classical CNNs, thanks to its properties of locality, stationarity, compositionality, which well match prior knowledge on many kinds of data and thus allow effective weight reuse. For this reason, defining an operation with similar characteristics for graph-structured data is of primary importance in order to obtain effective graph neural networks. The literature has identified two main classes of approaches to tackle this problem, namely spectral or spatial. In the former case [15, 16, 17], the convolution is defined in the spectral domain through the graph Fourier transform [18]. Fast polynomial approximations [16] have been proposed in order to obtain an efficient implementation of this operation. Graph-convolutional neural networks (GCNN) with this convolution operator have been successfully applied in problems of semi-supervised node classification and link prediction [17, 19]. The main drawback of these methods is that the graph is supposed to be fixed and it is not clear how to handle the cases where the structure varies. The latter class of approaches overcomes this issue by defining the convolution operator in the spatial domain [20, 21, 22, 23, 24, 25]. In this case, the convolution is performed by local aggregations, i.e. a weighted combination of the signal values over neighboring nodes. Since in this case the operation is defined at a neighborhood level, the convolution remains well-defined even when the graph structure varies. Many of the spatial approaches present in the literature [22, 23, 24] perform local aggregations with scalar weights. Instead, [20] proposes to weight the contributions of the neighbors using edge-dependent matrices. This makes the convolution a more general function, increasing its descriptive power. For this reason, in this paper we employ the convolution operator proposed in [20]. However, in order to obtain an efficient operation, we introduce several approximations that reduce its computation complexity, memory occupation, and mitigate vanishing gradient issues that arise when trying to build very deep architectures.

II-B Image denoising

The literature on image denoising is vast, as it is one of most classic problems in image processing. Focusing on the recent developments, we can broadly define two categories of methods: model-based approaches and learning-based approaches.

Model-based approaches traditionally focused on defining hand-crafted priors to carefully capture the salient features of natural images. Early works in this category include total variation minimization [26], and bilateral filtering [27]. Non-local means [28] introduced the idea of non-local averaging according to the similarity of local neighborhood. The popular BM3D [5] expanded on the idea by collaborative filtering of the matched patches. WNNM [29] used nuclear norm minimization to enforce a low-rank prior. Finally, some works recently introduced graph-based regularizers [30] to enforce a measure of smoothness of the signal across the edges of a graph of patch or pixel similiarities. Many of the most successful model-based approaches are non-local, i.e., they exploit the concept of self-similarity among structures in the image beyond the local neighborhood.

Learning-based approaches use training data to learn a model for natural images. The popular K-SVD algorithm [31] learns a dictionary in which natural patches have a sparse representation, and therefore casts image denoising as a sparse coding problem on this learned dictionary. The TNRD method [32] uses a nonlinear reaction diffusion model with trainable filters. An early work with neural networks [33] used a multilayer perceptron discriminatively trained on synthetic Gaussian noise and showed significant improvements over model-based methods. More recently, CNNs have achieved remarkable performance. Zhang et al. [6] showed that the residual structure and the use of batch normalization [34] in their DnCNN greatly helps the denoising task. Following the DnCNN, many other architectures have been proposed, such as RED [7], MemNet [8] and a CNN working on wavelet coefficients [9]. However, those CNN-based methods are limited by the local nature of the convolution operation, which is unable to increase the receptive field of a neuron-pixel to model non-local image features. This means that CNNs are unable to exploit the self-similar patterns that were proven to be highly successful in model-based methods. Very recently, a few works started addressing this issue by trying to incorporate non-local information in a CNN. NN3D [10] uses a global post-processing stage based on a non-local filter after the output of a denoising CNN. This stage performs block matching and filtering over the whole image denoised by the CNN. This is clearly suboptimal as the non-local information does not contribute to the training of the CNN. UNLNet [11] introduces a trainable non-local layer which collaboratively filters image blocks. However, performance is limited by the selection of matching blocks from the noisy input image instead of the feature space, and ultimately UNLNet does not improve over the performance of the simpler DnCNN. N3Net [12] introduces a continuous nearest-neighbor relaxation to create a non-local layer. Finally, NLRN [13] proposes a non-local module that uses the distances among hidden feature vectors of a search window around the pixel of interest to aggregate such vectors and return the output features of the pixel. However, there are significant differences with respect to the work in this paper. First, they use all the pixels in the search window instead of only a number of nearest neighbors, which means that their receptive field cannot dynamically adapt to the content of the image. Then, while in both works the feature aggregation weights are dynamically computed from the features themselves, NLRN uses an explicitly-parameterized function with learnable parameters, in contrast to this work where the function is fully learned as a dedicated sub-network. These choices increase the adaptivity of the proposed non-local operations, which result in better performance around edges.

III Proposed denoiser

III-A Overview

An overview of the proposed graph-convolutional denoiser network (GCDN) can be seen in Fig. 1. The structure will be explained more in detail in Sec. III-D where an analogy is drawn between unrolled proximal gradient descent with a graph total variation regularizer and the proposed network architecture. At a first glance, the network has a global input-output residual connection whereby the network learns to estimate the noise rather than successively clean the image. This has been shown [6] to improve training convergence for the denoising problem.

The main feature of the proposed network is the use of graph-convolutional layers where the graphs are dynamically computed from the feature space. The graph-convolutional layer, described in Sec. III-B, creates a non-local receptive field for each pixel-neuron, so that pixels that are spatially distant but similar in the feature space created by the network can be merged.

An important block of the proposed network is the preprocessing stage at the input. It can be noticed that the first layers of the network are classic 2D convolutions rather than graph convolutions. This is done to create an embedding over a receptive field larger than a single pixel and stabilize the graph construction operation, which would otherwise be affected by the input noise. The preprocessing stage has three parallel branches that operate on multiple scales, in a fashion similar to the architectures in [35] and [36]. The multiscale features are extracted by a sequence of three convolutional layers with filters of size $3\times 3$, $5\times 5$, and $7\times 7$, depending on the branch. After a final graph-convolutional layer, the features are concatenated.

The remaining network layers are grouped into an HPF block and multiple LPF blocks, named after the analogy with highpass and lowpass graph filters described in Sec. III-D. These blocks have an initial $3\times 3$ convolutional layer followed by three graph-convolutional layers sharing the same graph constructed from the output of the convolutional layer. All layers are interleaved by Batch Normalization operations [34] and leaky ReLU nonlinearities. Notice that the LPF blocks have themselves a residual connection to help backpropagation, as in ResNet architectures [37]. The final layer is a graph-convolutional layer mapping from feature space to the image space.

III-B Graph-convolutional layer

The operation performed by the graph-convolutional layer is summarized in Fig. 2. The two inputs to the graph-convolutional layer are the feature vectors $\mathbf{H}^{l}\in\mathbb{R}^{F^{l}\times N}$ associated to the $N$ image pixels at layer $l$ and the adjacency matrix of a graph connecting image pixels. In this work, the graph is constructed as a $K$-nearest neighbor graph in the feature space. For each pixel, the Euclidean distances between its feature vector and the feature vectors of pixels inside a search window are computed and an edge is drawn between the pixel and the $K$ pixels with smallest distance. Using this method, we obtain a $K$-regular graph $\mathcal{G}^{l}(\mathcal{V},\mathcal{E}^{l})$, where $\mathcal{V}$ is the set of vertices with $|\mathcal{V}|=N$ and $\mathcal{E}^{l}\subseteq\mathcal{V}\times\mathcal{V}$ is the set of edges. We also assume that the edges of $\mathcal{G}^{l}$ are labeled, i.e. there exists a function $\mathcal{L}:\mathcal{E}^{l}\to\mathbb{R}^{F^{l}}$ that assigns a label to each edge. In this work, we define the edge labeling function as the difference between the two feature vectors, i.e. $\mathcal{L}(i,j)=\mathbf{H}_{j}^{l}-\mathbf{H}_{i}^{l}=\mathbf{d}^{l,j\to i}$. A classic $3\times 3$ local convolution processes the local neighborhood to provide its estimate of the output feature vector for the current pixel, while the feature vectors of the non-local pixels connected by the graph are aggregated by means of the edge-conditioned convolution (ECC) [20]. Notice that the 8 local neighbors of the pixel are excluded from graph construction as they are already used by the local convolution. The non-local aggregation is computed as:

[TABLE]

where $\mathcal{F}_{\mathbf{w}^{l}}^{l}:\mathbb{R}^{F^{l}}\to\mathbb{R}^{F^{l+1}\times F^{l}}$ is a fully-connected network that takes as input the edge labels and outputs the corresponding weight matrix $\bm{\Theta}^{l,j\to i}=\mathcal{F}_{\mathbf{w}^{l}}^{l}\left(\mathcal{L}(i,j)\right)\in\mathbb{R}^{F^{l+1}\times F^{l}}$, $\mathbf{w}^{l}$ are the weights parameterizing network $\mathcal{F}^{l}$, and $\mathcal{S}_{i}^{l}$ is the set of neighbors of node $i$ in the graph $\mathcal{G}^{l}$. The scalar $\gamma^{j\to i}$ is an edge-attention term computed as:

[TABLE]

where $\delta$ is a cross-validated hyper-parameter. This term is reminiscent of the edge attention mechanism from the graph neural network literature [38] and it serves the purpose of stabilizing training by underweighting the edges that connect nodes with distant feature vectors. Note that this term could, in principle, be learned by the $\mathcal{F}$ network but we found that decoupling it and making it explicitly dependent on feature distances in an exponential way, accelerated and stabilized training. Also notice that in Sec. IV we show that this term alone, i.e. without weight matrices $\bm{\Theta}$, is not powerful enough to reach good performance. Moreover, it is worth mentioning that the edge weights $\bm{\Theta}$ and the edge-attention term $\gamma$ depend only on the edge labels. This means that two pairs of nodes with the same edge labels will have the same weights, resulting in a behaviour similar to weight sharing in classical CNNs.

Finally, we combine the feature vector estimated by the non-local aggregation with the one produced by the local convolution to provide the output features as follows

[TABLE]

where $\mathbf{H}_{i}^{l+1,\mathrm{L}}$ is the output of the $3\times 3$ local convolution for the node $i$ and $\mathbf{b}^{l}\in\mathbb{R}^{F^{l}}$ is the bias.

The advantages of the ECC with respect to other definitions of graph convolution are trifold: i) the edge weights depend on the edge label, ii) it allows to compute an affine transformation along every edge, and iii) the edge weight function is highly general since it does not have a predefined structure. By making the edge weights depend on the input features, the ECC implements an adaptive filter which can be more complex than the non-adaptive local filters. Moreover, the second advantage is due to the fact that $\bm{\Theta}^{l,j\to i}$ is an edge-dependent matrix, making the convolution operation more general than other non-local aggregation methods using scalar edge weights. Among such methods we can find GCN [17], GIN [22], MoNet [23], and FeastNet [24]. Finally, the $\mathcal{F}$ function is a general function which can be learned to be the optimal one for the denoising task by the function approximation capability of the subnetwork implementing it. This is in contrast with other methods where the function predicting the edge weights is fixed with some learnable parameters. For example, FeastNet [24] employs scalar edge weights computed using the following function

[TABLE]

where $\mathbf{u},\mathbf{v}\in\mathbb{R}^{F^{l}}$ and $c\in\mathbb{R}$ are learnable parameters. Instead, MoNet [23] employs a Gaussian kernel as follows

[TABLE]

where $\bm{\Sigma}\in\mathbb{R}^{F^{l}\times F^{l}}$ and $\bm{\mu}\in\mathbb{R}^{F^{l}}$ are learnable parameters. Also NLRN [13] uses a Gaussian kernel to perform non-local aggregations. We can consider this operation as a graph convolution where each pixel is connected to all the other pixels in its search window and the edge weights are defined as follows

[TABLE]

where $\mathbf{W}_{\theta},\mathbf{W}_{\phi}\in\mathbb{R}^{t\times F^{l}}$ and $\mathbf{W}_{g}\in\mathbb{R}^{F^{l+1}\times F^{l}}$ are learnable parameters.

III-C Lightweight Edge-Conditioned Convolution

As seen in the previous section, the function $\mathcal{F}$ has a key role in the ECC because it defines the weights for the neighborhood aggregation. In the original definition of ECC [20], the function $\mathcal{F}$ is implemented as a two-layer fully connected network. This definition raises some relevant issues. In the following, we will describe in detail these issues and present two possible solutions.

III-C1 Circulant approximation of dense layer

The first issue is related to the risk of over-parameterization. The dimension of the input of the $\mathcal{F}$ network is $F^{l}$, while the dimension of its output is $F^{l+1}\times F^{l}$. This means that the number of weights of the network depends cubically on the number of features. Therefore, the number of parameters quickly becomes excessively large, resulting in vanishing gradients or overfitting.

To address the over-parameterization problem we propose to use a partially-structured matrix for the last layer, instead of an unstructured one. We impose that this matrix is composed of multiple stacked partial circulant matrices, i.e., matrices where only a few shifted versions of the first row are used instead of all the possible ones of the full square matrix. Fig. 3 shows the structure of the approximated matrix. Using this approximation, the only free parameters are in the first row of each partial circulant matrix. If only $m$ shifts per partial circulant matrix are allowed, we reduce the number of parameters by a factor $m$. Thus, if the unstructured dense matrix has $F^{l}F^{l+1}\times F^{l}$ parameters, with the proposed approximation the number of parameters drops to $\frac{F^{l}F^{l+1}}{m}\times F^{l}$. Similar approaches to approximate fully connected layers have already been studied in the literature [39, 40]. In particular, [39] shows that imposing a partial circulant structure does not significantly impact the final performance in a classification problem. Indeed, there are connections with results stating that random partial circulant matrices implement stable embeddings almost as well as fully random matrices [41, 42, 43].

III-C2 Low-rank node aggregation

The second issue related to the $\mathcal{F}$ network regards memory occupation and computations. In order to perform the ECC operation, we have to compute a weight matrix $\bm{\Theta}^{l,j\to i}$ for each edge $j$ of every neighborhood $\mathcal{N}_{i}$ of every image in the batch. If we consider a $K$-regular graph and a batch of $B$ images with $N$ pixels each, the memory occupation needed to store all the matrices $\bm{\Theta}^{l,j\to i}$ as single-precision floating point tensors is equal to $B\times N\times K\times F^{l+1}\times F^{l}\times 4$ bytes and this quantity can easily become unmanageable. To give an idea of the required amount of memory, let us consider an example with $B=16$, $N=1024$, $K=8$, $F^{l}=F^{l+1}=66$, then the memory required to store all the matrices $\bm{\Theta}^{l,j\to i}$ for only one graph-convolutional layer is around 2 GB.

In order to solve this issue, we propose to impose a low-rank approximation for $\bm{\Theta}^{l,j\to i}$. Let us consider the singular value decomposition of a matrix

[TABLE]

where $\bm{\phi}_{s}$ and $\bm{\psi}_{s}$ are the left and right singular vectors and $\lambda_{s}$ the singular values. We can obtain a low-rank approximation of rank $r$ by keeping only the $r$ largest singular values and setting the others to zero. Therefore, the approximation is reduced to a sum of $r$ outer products. Inspired by this fact, we define $\bm{\Theta}^{l,j\to i}$ as follows

[TABLE]

where $\bm{\theta}_{s}^{j\to i,L}\in\mathbb{R}^{F^{l}}$, $\bm{\theta}_{s}^{j\to i,R}\in\mathbb{R}^{F^{l+1}}$, $\kappa_{s}^{j\to i}\in\mathbb{R}$ and $1\leq r\leq F^{l}$. Notice that the approximation in (3) ensures that the rank is at most $r$ rather than exactly enforcing a rank-$r$ structure, because we do not impose orthogonality between $\bm{\theta}_{s}^{j\to i,L}$ and $\bm{\theta}_{s}^{j\to i,R}$, even though random initialization makes them quasi-orthogonal. Using this approximation, we can redefine the $\mathcal{F}$ network in such a way that it outputs $\bm{\theta}_{s}^{j\to i,L},\bm{\theta}_{s}^{j\to i,R},\kappa_{s}^{j\to i}$ for $s=1,2,\dots,r$. In particular, we redefine the second layer of the $\mathcal{F}$ network: instead of having a single fully connected layer that outputs the entire matrix $\bm{\Theta}^{l,j\to i}$, we have three parallel fully connected layers that separately output $\bm{\theta}_{s}^{j\to i,L}$, $\bm{\theta}_{s}^{j\to i,R}$ and $\kappa_{s}^{j\to i}$, as shown in Fig. 4. The advantage of this approximation is that we only need to store $\bm{\theta}_{s}^{j\to i,L}$, $\bm{\theta}_{s}^{j\to i,R}$ and $\kappa_{s}^{j\to i}$ instead of the entire matrix $\bm{\Theta}^{l,j\to i}$, drastically reducing the memory occupation to $B\times N\times K\times r(2F^{l}+1)\times 4$ bytes. If we consider the example presented above and set $r=10$, the memory requirement drops from 2 GB to 700 MB. Another advantage of this approximation is that it also leads to a significant reduction of the computation burden, because we never have to actually compute all the matrices $\bm{\Theta}^{l,j\to i}$. In fact, the neighborhood aggregation can be reduced as follows

[TABLE]

where the computational cost of the full operation on the first line is $O(F^{l}F^{l+1})$, instead the cost of the decoupled operation on the second line is $O(r(F^{l}+F^{l+1}))$. Finally, this approximation also helps to reduce the number of parameters of the last layer of the $\mathcal{F}$ network since the output has size $r(F^{l}+F^{l+1}+1)$ instead of $F^{l+1}F^{l}$.

When we employ the new structure of the $\mathcal{F}$ network, we need to pay special attention to the weight initialization. In particular, we have to carefully define the variance of the random weight initialization of the three parallel layers to avoid scaling problems. We define $\mathbf{W}_{0}$ as the weight matrix of the first layer of the $\mathcal{F}$ network, and $\mathbf{W}^{L}$, $\mathbf{W}^{R}$ and $\mathbf{W}^{\kappa}$ as the weight matrices of the three parallel fully connected layers. Let us suppose that $\mathbf{d}^{j\to i}_{t}$ has been normalized to be approximately a standard Gaussian, i.e., $\mathbf{d}^{j\to i}_{t}\sim\mathcal{N}(0,1)$ for $t=1,\dots,F^{l}$, and that $\mathbf{W}^{0}$ has been initialized using Glorot initialization [44], i.e., $\mathbf{W}^{0}_{uv}\sim\mathcal{N}\left(0,\frac{1}{F^{l}}\right)$ with $u,v=1,\dots,F^{l}$. Let us also assume that $\mathbf{W}^{L}_{uv}\sim\mathcal{N}(0,\sigma^{2}_{L})$, $\mathbf{W}^{R}_{uv}\sim\mathcal{N}(0,\sigma^{2}_{R})$, and $\mathbf{W}^{\kappa}_{u}\sim\mathcal{N}(0,\sigma^{2}_{\kappa})$. Then, we obtain

[TABLE]

where $s=1,\dots,r$. Finally, considering the aggregation formula in Eq. (III-C2) leads to the following result:

[TABLE]

with $u=1,\dots,F^{l+1}$. In Eq. (5), we can observe that the variance of $\mathbf{H}_{i,u}^{l+1,\mathrm{NL}}$ depends on the fourth power of the number of features. This term can easily become extremely large, therefore it is important to set $\sigma^{2}_{L}$, $\sigma^{2}_{R}$ and $\sigma^{2}_{\kappa}$ in such a way that they can balance it. In this work, we set $\sigma^{2}_{L}=\sigma^{2}_{R}=\frac{1}{F^{l^{2}}}$ and $\sigma^{2}_{\kappa}=\frac{2}{r}$. This allows us to obtain $\mathbf{H}_{i,u}^{l+1,\mathrm{NL}}\sim\mathcal{N}(0,1)$ with $u=1,\dots,F^{l+1}$.

III-D Analogy with unrolled graph smoothness optimization

The neural network architecture presented in Sec. III-A can be seen as a generalization of few iterations of an unrolled proximal gradient descent optimization method, which is widely used to solve linear inverse problems in the form of

[TABLE]

being $\mathbf{x}$ the clean image, $\mathbf{A}$ a forward model (e.g., a degradation such as blurring, downsampling, compressed sensing, etc.) and $\mathbf{n}$ a noise term. A well-known technique to recover $\mathbf{x}$ from $\mathbf{y}$ is to cast the problem as a least-squares minimization problem with a regularization term that models some prior knowledge about the image. One such regularizer is graph smoothness. Considering a graph with Laplacian matrix $\mathbf{L}$ where edges connect pixels that are deemed correlated according to some criterion, the graph smoothness $\mathbf{x}^{T}\mathbf{L}\mathbf{x}$ is the graph equivalent of the total variation measure, indicating how much $\mathbf{x}$ varies across the edges of the graph. Natural images where the graph connects the local neighborhood typically have lowpass behavior, resulting in a low graph smoothness value. Reconstruction is therefore cast as:

[TABLE]

The functional in Eq. (7) is in the form of a sum of two terms ($f(\mathbf{x})+g(\mathbf{x})$) and can be minimized by means of proximal gradient descent [45] which alternates a gradient descent step over $f$ and a proximal mapping over $g$:

[TABLE]

Solving for the proximal mapping operator results in the following update equation:

[TABLE]

In order to match the framework of residual networks, let us define the least-squares solution $\mathbf{x}^{n}=\mathbf{A}^{+}\mathbf{y}=\left(\mathbf{A}^{T}\mathbf{A}\right)^{-1}\mathbf{A}^{T}\mathbf{y}$ and perform a change of variable whereby the optimization estimates the residual of the least squares solution, i.e., $\bm{\nu}^{(t)}=\mathbf{x}^{n}-\mathbf{x}^{(t)}$. Hence, we can rewrite Eq. (8) as:

[TABLE]

Finally, the following update equation can be derived:

[TABLE]

This update can be visualized as in Fig. 5a and is composed of two major operations involving the signal prior:

$\mathbf{L}\mathbf{x}^{n}$: the graph Laplacian can be seen as a graph highpass filter applied to $\mathbf{x}^{n}$; 2. 2.

$\left(\mathbf{I}+\beta\mathbf{L}\right)^{-1}$: this term can be seen as a graph lowpass filter. In order to see this, let us use the matrix inversion lemma as $\left(\mathbf{I}+\beta\mathbf{L}\right)^{-1}=\left(\mathbf{I}+\beta\mathbf{U}\bm{\Lambda}\mathbf{U}^{H}\right)^{-1}=\mathbf{I}-\mathbf{U}\left(\beta^{-1}\bm{\Lambda}^{-1}+\mathbf{I}\right)^{-1}\mathbf{U}^{H}=\mathbf{U}\left[\mathbf{I}-\left(\beta^{-1}\bm{\Lambda}^{-1}+\mathbf{I}\right)^{-1}\right]\mathbf{U}^{H}$, where $\mathbf{U}$ is the graph Fourier transform. The term $\mathbf{I}-\left(\beta^{-1}\bm{\Lambda}^{-1}+\mathbf{I}\right)^{-1}$ is a diagonal matrix whose entries are equal to $\frac{1}{\beta\lambda_{i}+1}$ where $\lambda_{i}$ are the eigenvalues of the graph Laplacian, and the lowpass behavior is due to decreasing value of such entries for increasing $\lambda$.

For the denoising problem, we can set $\mathbf{A}=\mathbf{I}$ and obtain the update shown in Fig.5b. The network architecture proposed in Sec. III-A draws from this derivation by unrolling a finite number of Eq. (9) iterations and generalizing the lowpass and highpass filters with learned graph filters interleaved by nonlinearities. In Sec. IV we experimentally show that the learned filters actually show an approximate highpass and lowpass behavior.

IV Experimental Results

IV-A Training details

The training protocol follows the one used in [6]. The network is trained with patches of size $42\times 42$ randomly extracted from 400 images from the train and test partitions of the Berkeley Segmentation Dataset (BSD) [46], withholding the 68 images in the validation set for testing purposes (BSD68 dataset). The loss function is the mean squared error (MSE) between the denoised patch output by the network and the ground truth. Each model is trained for approximately 800000 iterations with a batch size of 8. The Adam optimizer [47] has been used with an exponentially decaying learning rate between $10^{-4}$ and $10^{-5}$. The behavior of the graph-convolutional layer is slightly different between training and testing for efficiency reasons. During training all pairwise distances are computed among the feature vectors corresponding to the pixels in the patch. On the other hand, testing is “fully convolutional”, as every pixel has a search window centered around it and neighbors are identified as the closest pixels in such search window. The search window size is $43\times 43$, roughly comparable to the patch size used in training. This procedure is slightly suboptimal as some pixels might suffer from border effects during training (their search windows are not centered around them) but it is advantageous in terms of speed and memory requirements. Reflection padding is used for all 2D convolutions to avoid border effects. The $\delta$ parameter in the edge attention term in Eq. (2) is set to a value equal to 10. The number of features used in all convolutional layers is 132, except for the three parallel branches of the preprocessing stage which have 44 features. The number of circulant rows in the circulant approximation of dense layers in the $\mathcal{F}$ network is $m=3$. The low-rank approximation uses $r=11$ terms. During training, we noticed that the proposed lightweight ECC presented in Sec. III-C is extremely useful. In fact, without it, the network suffered from vanishing gradient problems even with a significantly lower number of layers.

IV-B Feature analysis

In this section we study the properties of the features in the hidden layers of the network.

IV-B1 Adaptive receptive field

We first analyze the characteristics of the receptive field of a single pixel. Since the proposed network employs graph-convolutional layers, the shape of the receptive field is not fixed as in classical CNNs, but it depends on the structure of the graph. In Fig. 6 we show two examples of the receptive field of a single pixel for the graph-convolutional layers in an LPF block with respect to the input of the block. Instead, in Fig. 7 we show the receptive field of a single pixel for the layers in the HPF and in the first LPF blocks with respect to the output of the preprocessing block. We can clearly see that the receptive field is adapted to the characteristics of the image: if we consider a pixel in a uniform area, its receptive field will mostly contain pixels that belong to similar regions; instead if we consider a pixel on an edge, its receptive field will be mainly composed of other edge pixels. This is beneficial to the denoising task as it allows to exploit self-similarity and it descends from the use of a nearest neighbor graph, connecting each pixel to other pixels with similar features. Notice that differently from algorithms performing block matching in the pixel space, we compute distances between feature vectors which can capture more complex image characteristics. This can be seen in Fig. 8 where we compute the Euclidean distances between the feature vector of the central pixel and the feature vectors of the other pixels in the search window. We notice that the distances reflect the type of edge that includes the central pixel, e.g., a pixel sitting on a horizontal edge will detect as closest other pixels sitting on horizontal edges. This is due to the visual features learned by the network and would not happen in pixel-space matching. Thanks to the adaptability of the receptive field, graph convolution can be interpreted as a generalization of the block matching operation performed in other non-local denoising methods, such as BM3D [5].

IV-B2 Filter analysis

We also study the behavior of the LPF and HPF operators. In particular, we are interested in validating the analogy made in Sec. III-D. We compute the discrete Fourier transform (DFT) of the feature maps at the output of these operators. As an example, Fig. 9 shows the log-magnitude of the coefficients of three feature maps at the output of the HPF block and of the first LPF block. The energy of the DFT coefficients of the LPF feature maps is concentrated in the low frequencies, thus showing a lowpass behavior. Instead, the coefficients of the HPF feature maps show a typical highpass behavior, having the energy concentrated along few directions. This substantiates our claim that the learned convolutional layers actually approximate nonlinear highpass and lowpass operators.

IV-B3 Edge prediction

Lastly, we measure how much the true graph constructed by pixel or patch similarities on the noiseless image is successfully predicted by the graph constructed from the feature vectors in the hidden layers. In order to construct the true graph of the image, we first compute the average pixel value of a 5$\times$5 window centered at the considered pixel, for every pixel in the image, and then we use the obtained values to compute a nearest neighbor graph with Euclidean distances. We then compare the true graph with the graph computed in the hidden layers of the network. Fig. 10 shows the percentage of edges correctly identified as function of the number of neighbors considered for the true graph. We can notice that the accuracy of the prediction decreases in layers closer to the output. This is due to the fact that we use a residual network that estimates the noise instead of approximating the clean image. In fact, the network learns to successively remove the latent correlations in the feature space, and as a consequence, the graph becomes more random in the later layers.

IV-C Ablation studies

We study the impact of various design parameters on denoising performance. First, Table I shows the PSNR on the Set12 testing set as function of the number of neighbors used by the graph convolution operation for several values of the noise standard deviation $\sigma$. Each model has been independently trained for the specified number of neighbors. It can be noticed that increasing the number of neighbors improves the denoising performance up to a saturation point, and then the performance slightly decreases. This shows that there an optimal neighborhood size and that it is important to employ only a small number of neighbors, in order to select only pixels with similar characteristics. This is in contrast with the NLRN method which uses all the pixels in the search window.

Then, we study the relative impact on performance of the edge aggregation matrices $\bm{\Theta}$ in Eq. (III-B) with respect to using only the edge attention scalar $\gamma$. Table II reports the PSNR achieved on Set12 by the proposed method with the non-local aggregation performed as in Eq. (III-B) and a variant where the aggregation is computed as:

[TABLE]

Both methods use a non-local graph with 8 nearest neighbors. We can notice that the edge attention term alone achieves a worse PSNR with respect to GCDN by approximately 0.2 dB, even though it improves over a model without non-local neighbors (see Table I for the corresponding 0-NN value). This shows the advantage of using a trainable affine transformation, such as $\bm{\Theta}$ in Eq. (III-B), instead of a scalar weight function with a predefined structure.

Finally, we remark that we do not compare with respect to the full ECC without the approximations introduced in Sec. III-C because it suffers from vanishing gradient problems, rendering training unstable even for a much smaller number of layers, and it would be computationally prohibitive.

IV-D Comparison with state of the art

In this section we compare the proposed network with state-of-the-art models for the Gaussian denoising task of grayscale images. We train an independent model for each noise standard deviation, which is assumed to be known a priori for all methods. We fix the number of neighbors for the proposed method to 16. The reference methods can be classified into model-based algorithms such as BM3D [5], WNNM [29], TNRD [32] and recent deep-learning methods such as DnCNN [6], N3Net [12] and NLRN [13]. In particular, among the deep-learning methods, N3Net and NLRN propose non-local approaches. All results have been obtained running the pretrained models provided by the authors, except for N3Net at $\sigma=15$ which is unavailable. Table III reports the PSNR and SSIM values obtained for the Set12, BSD68 and Urban100 standard test sets. It can be seen that the proposed method achieves state-of-the art performance and works especially well at low to medium levels of noise. This can be explained by a higher difficulty in constructing a meaningful graph from the noisy image at higher noise levels. We also notice that the proposed method achieves strong results on the Urban dataset. This dataset contains higher resolution images with respect to the other two and is mainly composed of photos of buildings and other regular structures where exploiting self-similarity is very important. In addition, it is also worth mentioning that the proposed method provides a better visual quality. In many cases, the proposed method has a higher SSIM score, even if NRLN has better performance in terms of PSNR. This can also be noticed in Fig. 11, which shows a visual comparison on an image from the Urban100 dataset. In general, the images produced by the proposed algorithm present sharper edges and smoother content in uniform areas. We can notice that many areas in the photos from Urban100 have approximately piecewise smooth characteristics. It is well known that image processing algorithms based on graphs are well suited for piecewise smooth content (see, e.g., [48, 49] in the context of compression and [30] for denoising). To further show this point, we study the performance of the proposed method for denoising of depth maps, e.g., generated by time-of-flight cameras. The OGLR algorithm [30] based on a graph smoothness regularizer achieved state-of-the-art results among model-based algorithms for this specific task where it is essential to preserve edge sharpness while simultaneously smoothing the flat areas. Table IV reports the PSNR and SSIM results achieved on a standard set of depth maps111http://vision.middlebury.edu/stereo/data/.. It can be seen that the proposed method outperforms both NLRN and OGLR, even at high levels of noise. Also, we can notice that OGLR displays competitive performance at low noise levels, but its visual quality significantly degrades when in presence of stronger noise. Fig. 12 shows a visual comparison where it can be seen that GCDN produces sharper edges while also providing a very smooth background.

IV-E Real image denoising

Real image noise is generally more challenging than synthetic Gaussian noise. There are multiple contributions such as quantization noise, shot noise, fixed-pattern noise [50, 1], dark current, etc. that make it overall signal-dependent. It has been observed [51, 52] that deep learning methods trained on synthetic Gaussian noise perform poorly in presence of real noise. However, suitable retraining with real data generally improves their performance. In this section, we study the behavior of the proposed network in a blind denoising setting with real noisy images acquired by smartphones. We retrain the proposed method, NLRN and DnCNN on the SIDD dataset [52] composed of 30000 high-resolution images acquired by smartphone cameras at varying illumination and ISO levels. The authors provide clean and carefully registered ground truths for all the available scenes, so that it is possible to perform a supervised training. We create training and testing subsets from the sRGB images in the SIDD dataset by selecting a range of noise levels. Our training set is composed of 3500 crops of size $512\times 512$ whose RMSE with respect to the ground truth is below 15. The testing set is composed of 25 random crops of size $512\times 512$ with noise in the same range as the training set. Table V reports the results for CBM3D [53], DnCNN, NLRN and the proposed GCDN. Notice that CBM3D is not a blind method, so we provide an estimate of the noise standard deviation, as computed by a noise estimation algorithm [54]. We can notice that the proposed method achieves better results and this is confirmed by the visual comparison in Fig. 13.

V Conclusions

In this paper, we presented a graph-convolutional neural network targeted for image denoising. The proposed graph-convolutional layer allows to exploit both local and non-local similarities, resulting in an adaptive receptive field. We showed that the proposed architecture can outperform state-of-the-art denoising methods, achieving very strong results on piecewise smooth images. Finally, we have also considered a real image denoising setting, showing that the proposed method can provide a significant performance gain. Future work will focus on extending the proposed architecture to other inverse problems, such as super-resolution [55, 56].