The invention relates generally to signal processing, and more particularly to estimating parameters of a graph spectral filter and processing a signal using such a graph-based filter,
Graph-based signal processing merges algebraic and spectral graph theoretic concepts with computational harmonic analysis to process signal using graphs.
An undirected graph G=(V, E) includes a set of nodes V={1, 2, . . . , N} connected by a set of edges E:{(i,j,wij)}, where (i, j, wij) the edge between nodes i and j having weights wij. An adjacency matrix W of the graph is an N×N matrix. A degree di of a node i is a sum of edge weights connected to node i.
A degree matrix D is a diagonal matrix: diag{d1, d2, . . . , dN}. A combinatorial Laplacian matrix is :=D−W. The normalized Laplacian matrix L:=D−1/2D<1/2 is a symmetric positive semi-definite matrix. The matrix L has an eigendecomposition L=U Λ Ut, where Ut is a graph Fourier transform {u1, u2, . . . , uN} that is composed of an orthogonal set of eigenvectors, and Λ=diag{λl, . . . λN} is a corresponding eigenvalue matrix. In the above notation Ut, the superscript t indicated a transpose of the matrix U.
The eigenvectors and eigenvalues of the Laplacian matrix provide a spectral interpretation of signals represented by the graph. The eigenvalues can be treated as graph frequencies that are always in the interval [0,2] on a real line.
Graph Fourier transform (GFT) and Graph Spectral Filter (GSF)
A signal {circumflex over (x)} is normalized as x=D1/2 {circumflex over (x)}. A graph Fourier transform (GFT) is defined as a projection of the normalized signal x onto the eigenvectors U of the graph:
{tilde over (x)}:=Utx, (1)
where {tilde over (x)} are the GFT coefficients. The inverse transform is
x:=U{tilde over (x)}, (2)
The signal x is normalized because the eigenvector U is defined, over the normalized Laplacian matrix L. By definition, W, D and can be regarded as operators over the unnormalized signal {circumflex over (x)}. The unnormalized signal {circumflex over (x)} can also be processed as
D½LD½{circumflex over (x)}=D½UΛUtD½{circumflex over (x)}=D½UΛ(Utx).
It is meaningful to analyze a transform by applying Ut over the normalized signal x, which is the GFT. The normalizer D1/2 is determined in the context of the graph G.
Normalization of an unnormalized signal is represented as
x:=D1/2{circumflex over (x)}, (3)
and the denormalization is represented as
{circumflex over (x)}:=D−1/2x, (4)
With the GSF and the inverse GSF, the signal x can be converted between a spatial domain and a graph spectral domain. The GSF is defined by the graph G, its adjacency matrix W, or the Laplacian matrix L.
As described later, any graph G has a corresponding filter design. Hence, a graph can be named by its corresponding filter name, e.g. bilateral graph filter refers to a graph with weights defined as bilateral coefficients.
Defining the graph does not change signal processing because it can be assumed as an equivalent representation of a conventional filter. However, with the GFT based on the graph G, a new spectral domain is defined. It is possible to perform the GST in the graph spectral domain as defined by
{tilde over (x)}out:=H1(Λ){tilde over (x)}in, (5)
where H1(Λ)=diag (h1(λi)). Or in another form:
{tilde over (x)}out(λi):=hi(λi){tilde over (x)}in(λi),
where h1(λi) can be a lowpass highpass or any other type of filter in the graph spectral domain.
Correspondingly, in the spatial domain, the GSF includes the GFT transform, the GSF filtering, and the inverse GFT transform as described below:
xout=UH1(Λ)Utxin:=H()xin, (7)
The filtering kernels in the spatial domain and spectral domain are linked as
H()=UH1(Λ)Ut (8)
or H1(Λ)=UtH()U (9)
In the above equations, the GSF and inverse GSF are defined over normalized signals. For unnormalized signals, the eigendecomposition may be performed over a combinatorial Laplacian matrix, which has different eigenvectors and eigenvalues.
Considerations for GSF Filter Design
The GSF can be used for image filtering and other applications. To apply graph based signal processing described above to design a graph based image filter, the following steps can be performed.
First, a graph is constructed. Then, a graph spectral filter H1 is designed in terms of the graph. The design of the graph can use any known image filter F. For an input image {circumflex over (x)}in, the output image {circumflex over (x)}out is a weighted average of pixels
where wij are the weighting coefficients between pixel i and j. For example, a bilateral graph G has bilateral weights w defined by:
The corresponding graph is defined as G=(V, E) where the nodes V represent the pixels of the input image and the edges E represent the weighting (filtering) coefficients between two pixels i and j. The adjacency matrix W of the graph is W=[wij]. With these definitions, the degree matrix D, the combinatorial Laplacian matrix , the normalized Laplacian matrix L, the eigenvalue matrix Λ, and the eigenvalue matrix U can be determined.
Identity Graph Spectral Filter
In an identity graph spectral filter (iGSF), the signal is unchanged in the graph spectral domain. That is, the signal filtered by iGSF is
Eqn. (14) can be rewritten as
xout=(I−L)xin (15)
where xin is the normalized input image, and xout is the normalized output image, and I is an identity matrix.
The normalization can be used to define the image filtering in terms of the non-negative definite matrix L, and thus has a spectral interpretation. The normalization also ensures that a constant signal, when normalized, is an eigenvector of L associated with zero eigenvalue.
By comparing Eqn. (7) and Eqn. (15) H for an identity GSF can be represented as:
H=I−L (16)
In the frequency domain, H1 for the identity GSF can be represented as:
As shown in
In an example application used to describe the embodiments, the signals are images, for example, noisy images to be enhanced by reducing noise. Specifically, the noisy images can be depth images.
To better understand the identity GSF, Ĥ and Ĥ1 are detailed in the unnormalized domain:
From Eqn. (27), the filtering kernel in frequency domain is composed of three parts.
In Part. A, the spectral signal is first inverse transformed to the spatial domain such that the normalizer is applied and then transformed to spectral domain.
The signal fed to part B is in the spectral domain, and Part B is given by H1=I−Λ in case of the iGSF. After filtering by H1 in the graph-spectral domain, the spectral signal is then fed to Part. C.
Part C inverse trans onus to the spatial domain such that the denormalizer is applied and then transformed hack to frequency domain.
The design of the graph G is an important task, because xout=UH1Utxin, U and Λ are all determined by the selected graph G.
Following the graph selection, another issue is to design the GSF filter H1 to satisfy requirements for particular applications. H1 can be a lowpass, highpass or other types of filter. This disclosure addresses the filter parameter estimation problem for the GSF filter.
High-Pass Graph Filter
In a high-pass graph Filter, the graph filter H is used to extract high frequency components from an input image xout=Hxin=UH1Utxin.
Overview of Graph Based Interpolation
Conventional methods for interpolating signals defined on graphs include a method based on sampling theorem for bandlimited graph signals, and a method based on regularization.
Interpolation Based on Graph Sampling Theorem
In the sampling theorem method, data interpolation is formulated as a signal reconstruction problem on a graph, where a graph signal is defined as the information associated with each node, e.g., scalar or vector values mapped to the set of vertices and edges of the graph. In the graph sampling method, w stands for the weight between nodei and j. The known set is S, and the unknown set is Sc=V−S. A signal is bandlimited to the graph frequency band [0, ω) on a graph G when its GFT Ut has support only at frequencies [0, ω). The space of ω-bandlimited signals is called the Paley-Wiener space and is denoted by PWω(G),
The partially known graph function can be formulated as a downsampled-upsampled (DU) signal. In a regular signal domain, the original signal is recovered from its DU signal by applying a low-pass filter. Similarly, in the graph domain, a low-pass filters can be designed to recover the original graph signal from the DU signal. An optimal cut-off frequency ω* can be determined as
ω*=σmin, (28)
where σ2min is a smallest singular value of (2)Sc.
Given the cut-off frequency ω*, the reconstruction of the graph signal is done by a least-square (LS) projection of the corresponding DU signal onto the PWω(G) space. Let k be the number of eigenvalues of the Laplacian matrix less than ω*.
The ω-bandlimited signal can be written, under appropriate permutation, as:
where Uk is the matrix containing the first k eigenvectors, and (Uk)s is the sub-matrix of Uk containing rows corresponding to the set S. Then, α can be obtained by determining the least square solution to:
x(S)=(Uk)sα (31)
That is,
α=(Uk)S†x(S) (32 )
where † indicates a pseudo-inverse of the matrix.
Hence, the unknown signal values are
x(Sc)=(Uk)S
x(Sc)=(Uk)S
Interpolation Based on Regularization
In the regularization method, the data interpolation is formulated as a regularization problem
denotes a downsampling operator. M is a size of an known subset of samples, and N is a size of a full set of samples. J can be represented by an appropriate permutation:
JM×N=(IM×M|0M×(N−M))M×N (36)
and xdu can be in the following form with appropriate permutation
xdu=[x(S)t, 0(Sc)t]t,
Then, H is a high pass regularization graph filter. Part A in equation 35 is a data-fitting term that determines the error between the reconstructed signal and the original signal at the known samples, and part B is the Euclidean norm of the output of a highpass graph filter H.
The solution is
Iterative methods can be used to solve the least square problem at a reduced complexity.
Given an underlying graph constructed for a signal as well as a data set of unenhanced and corresponding enhanced training signals, a graph spectral filter is estimated in the spectral domain. A key realization for designing an appropriate graph filter is to determine an optimal graph spectral filter (GSF).
Instead of using ad-hoc selection methods, the embodiments of the invention estimate the filter parameters by formulating the filter design as a solution to a least square (LS) problem. In an application where training data are available, the method provides a numerical way to estimate the optimal GSF filter.
The embodiments of the invention provide methods for processing a signal using a graph spectral filter (GSF).
In an example application used to describe the embodiments, the signals are images, for example, noisy images to be enhanced by reducing noise. Specifically, the noisy images can be depth images.
The invention can be used for applications such as image and video processing, and for signal processing in social, energy, transportation, sensor, or neural networks, where high-dimensional data naturally can be represented as vertices of weighted graphs. The invention can also be used to process epidemiological, census, logistics, and medical data.
As described above, the design of a filter H or a functional ƒ(H) for a signal represented as a graph, equivalently g(U, H1)) includes two phases. In a first phase, a graph transform domain is selected, based on the graph G for the signal to be processed, and a graph spectral filter (GSF) filter H1 is designed in the graph spectral domain using training data. We assume that the graph G is predetermined.
The embodiments of our invention estimate parameters H1 for the. GSF filter using training, data. The parameters are used to design the graph filter. Though there are many known filter designs that can be used, the selection is based on certain assumptions, e.g., the signal has certain properties or follows a particular distribution, which makes the selection difficult in practice.
Method Overview
Training
It is noted that training data are available for sonic applications. In this case, the selection, is simpler. For example, in a 3D video compression system, noise in depth images can come from down-upsampling due to a mixed resolution coding. If enhanced depth images, before encoding, are available along with the reconstructed depth images, then it is possible to form training data and use the images to train the parameters of the graph filter. Subsequently, the estimated graph filter can be used to denoise the depth images at a decoder without having access to the original depth images.
The embodiments of the invention provide a method for estimating the graph filter H or a GSF filter H1. It should be understood that the invention can be applied to other graph-based signal processing applications.
Training data includes M image pairs, each image has N=Width×Height pixels. Each image pair includes a noisy image 212 and a corresponding enhanced (denoised or clean) image 211. The filter design is formulated as determining an optimal graph filter:
represents the noisy training images in columnar form (after down and upsampling (du)), and E represents the errors between the noisy images and the enhanced images.
Given a graph 220, e.g., a bilateral filter or a joint bilateral filter, the graph transform Ut is known. Inserting Eqn. (8) into Eqn. (39) is equivalently to determining a best H1:
where H1 is the GSF filter in the graph spectral domain. The filter can be estimated from the training data set by solving the least square (LS) problem.
Because H1 is a diagonal matrix, we vectorize the martix as
h1={h(λ1), h(λ2), . . . , h(λN)
so that a conventional least square solution can be used. For this purpose, the columns of EtU are stacked, and a vector y is
y={y1t, y2t, . . . , yNt}t (43)
Similarly, the columns of XtduU are used to define a matrix A with (M×N) rows and N columns.
A=diag{α1, α2, . . . , αN} (44)
With this notation, Eqn. (42) can be written as
with a solution
{dot over (h)}1=A†y (46)
where † represents the pseudo-inverse. Thus, we formulate the filter estimation problem in Eqn. (39) and (42) with a solution provided by Eqn. (46),
Filtering Method
The method is shown in the
Training
During training, die method first estimates 210 the graph spectral filter (GSF) 203 as follows.
Step 1: A graph G is defined 220 heuristically for an input signal 204.
Step 2: Next, we derive 230 the graph Fourier transform (GFT) Ut, where Ut is an orthogonal set of eigenvectors of the corresponding normalized Laplacian matrix L. The GFT transform is an input for the estimation 210.
Step 3: A normalizer/denormalizer D1/2 and D−1/2 are derived 240.
Step 4: For each training pair of images, e.g., an enhanced image 205 and a noisy image 205′, an error signal E is obtained as a function of the image pairs, e.g. the difference image of the image pair, or simply just the enhanced images.
Step 5: A signal y 206 is obtained as a function of signal U 207 and F. 208.
Step 6: A signal A 209 is obtained as a function of signal U and the noisy training image 205′.
Step 7: Given signals y and A as inputs, the GSF filter parameters 211 are estimated by the GSF estimator 210 module using, e.g., a least square procedure.
Then, all the parameters for the modules used during the real time operation as shown on the left of
Operation
During operation 202 the GSF 203 is used as follows.
Step 1 : Provide an input image 204 to filtered or enhanced by the GSF filter 203.
Step 2: Perform normalization 250 on the input image using Eqn, (3).
Step 3: Perform a GFT transform 260 on the normalized input image using Eqn (1). The image obtained in the transform domain is provided as input for the next step.
Step 4: Conduct GSF filtering using the estimated parameters H1 of the filter 203.
Step 5: Perform an inverse GFT (iGSF) transform 270 on the filtered image in GFT transform domain using Eqn. (2).
Step 6: Perform denormalization 280 using Eqn. (4).
Step 7: Output 290 the processed image.
In the above, steps 3, 4, and 5 (260, 203, 270) constitute the graph filter H 265 as described above.
Normalization/Denormalization
In another implementation of graph based filtering, the GFT Ut is based on the combinatorial Laplacian matrix , instead, the filtering based on the normalized Laplacian L. In this case, the normalizes and denormalizer are an identity I, or equivalently, the not and denormalization modules are skipped in this case.
Iteration
The invention is not limited to use a closed solution to an LS solution. In a practical application, an iterative solution of the least squares problem may be preferred.
Image Block Classifying
It may be suboptimal to use a single filter parameter set for all images. Hence, we can classify image blocks into different defined categories, e.g., a block with a “flat” texture or an edge. Then, each category can be associated with one parameter set. Each category maintains one set of parameters for the GSF filter,
High-pass Bilateral Graph Filter
In one embodiment, the graph filter H is used to extract high frequency components from an input image xout=Hxin=UH1Utxin. As done in the prior art, the graph G is based on a bilateral filter (BF) with a “star” structure of the graph as shown in
The parameters of the high-pass GSF filter H1 can be estimated using the GSF design according to this invention. In the training data set, the enhanced images 205 are the enhanced versions without high frequency components. The error signal E is set equal to the difference image between the input image and enhanced image. Then Eqn. (46) can he used to solve the parameter estimation for the graph fiber,
Low-Pass Bilateral Graph Filter
Similarly, a low-pass bilateral graph filter can be estimated if the error signal E is the enhanced image rather the difference image.
High-Pass Joint Bilateral Graph Filter (JBF)
In this embodiment, errors in a noisy depth image, e.g. due to coding artifacts, are filtered. In contrast to the bilateral filter, the weights of the joint bilateral filter are determined from a guidance image Ĩ, e.g., a corresponding texture image. The weights can he represented by Eqn. (47), with xin replaced by
Similar to the high-pass bilateral graph filter, our GSF design is applied to estimate the parameters of GSF filter H1 to filter the noisy depth. The training data set is composed of noisy depth images and enhanced depth images. The error signal E is a difference depth image between the noisy depth and enhanced depth. Eqn. (46) is then used to estimate the parameter for the graph filter.
Filter for Depth Denoising Based on Regularization
In this embodiment, the graph filter is used to denoise an image. This problem can be formulated as a regularization problem:
where xdu is a noisy image. H is a regularization kernel that is designed as a high-pass graph filter with the underlying graph being a joint bilateral as in Eqn. (48). That is, the weights of the graph are determined from the corresponding texture image.
In Eqn. (49), part A determines the error between the noisy depth and the filtered depth image. Part B is the Euclidean norm of the output of a high pass filter H. Parameter ρ is predetermined, e.g. ρ=ρ0 to control the weights between the parts A and B. The solution, in closed form, is
{dot over (x)}=(I+ρHtH)−1xdu (50)
The design of H (or H1) is as for the high-pass joint bilateral filter described above, and the resulted filter (I+ρHtH)−1 is referred as a graph-based joint-bilateral filter (GB-JBF).
For an example of 3D video system, as shown in
Sealing Parameter ρ
It is noted that the scaling parameter ρ in Eqn. (49) may be difficult to determine because it does directly have a physical meaning. One way to avoid setting ρ is to solve the following similar problem:
With this formation, parameter δ has a clear physical meaning that represents the changes made on the signals that can be tolerated.
GB-JBU Upsampling Based on Regularization
In another embodiment as shown in
Pixel-level Method vs. Block-level Method
In this embodiment, two methods are described for different graphs, namely pixel-level and block-level methods. For the pixel-level method as shown in
For the block-level method, the entire image is partitioned into non-overlapping patches and all missing samples in a patch are linked to other known samples (pixels) within a predetermined distance, and a single graph is constructed for all missing samples (pixels) Within one graph,
Parameter ρ
Similar to Eqn. (51), one way to avoid setting δ is to solve the following similar problem:
With this formation, parameter ρ has a clear physical meaning that represents the changes made on the known samples that can be tolerated.
Although the invention has been described by way of examples of preferred embodiments, it is to be understood that various other adaptations and modifications can be made within the spirit and scope of the invention. Therefore, it is the object of the appended s to cover all such variations and modifications as come within the true spirit and scope of the invention.
Entry |
---|
Sandryhaila12, Discrete Signal Processing on Graphs [online], Copyright (c) 2012 IEEE [retrieved on Aug. 24, 2015]. Retrieved from the Internet:<URL:http://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=10&ved=0CGIQFjAJahUKEwi9yOmbl8HHAhXMVT4KHSD-APU&url=http%3A%2F%2Fusers.ece.cmu.edu%2F˜asandryh%2Fpapers%2 Ftsp13.pdf>. |
D.I. Shuman, S.K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst. The emerging field of signal processing on graphs: Extending highdimensional data analysis to networks and other irregular domains. Signal Processing Magazine, IEEE, 30(3):83-98, 2013. |
Sunil K. Narang, Akshay Gadde, and Antonio Ortega. Signal processing techniques for interpolation in graph structured data. In IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), May 2013. |
Sunil K. Narang, Akshay Gadde, Eduard Sanou, and Antonio Ortega. Localized iterative methods for interpolation in graph structured data. To appear in 1st IEEE Global Conference on Signal and Information Processing (GlobalSIP), Dec. 2013. |
Akshay Gadde, Sunil K. Narang, and Antonio Ortega. Bilateral filter: Graph spectral interpretation and extensions. ARVIX, Mar. 2013. |
Number | Date | Country | |
---|---|---|---|
20150112897 A1 | Apr 2015 | US |