Signal reconstruction via operator guiding

TL;DR

This paper introduces a flexible operator-guided approach for signal reconstruction that enhances specific components, demonstrated through super-resolution examples using guiding images to improve resolution and detail.

Contribution

It proposes general guiding operators beyond orthogonal projectors, enabling practical and effective signal reconstruction, especially in super-resolution tasks.

Findings

Successful super-resolution of RGB and depth images using guiding operators

Guided reconstruction improves image resolution and detail

Method applicable to various signal enhancement scenarios

Abstract

Signal reconstruction from a sample using an orthogonal projector onto a guiding subspace is theoretically well justified, but may be difficult to practically implement. We propose more general guiding operators, which increase signal components in the guiding subspace relative to those in a complementary subspace, e.g., iterative low-pass edge-preserving filters for super-resolution of images. Two examples of super-resolution illustrate our technology: a no-flash RGB photo guided using a high resolution flash RGB photo, and a depth image guided using a high resolution RGB photo.

Equations

$$y=Ax+n,$$

$$\underset{x}{\min }\Vert Ax-y{\Vert }_{\omega }^2+\rho R(x),$$

$$A^T\omega (Ax-y)+\rho L(g)x=0,$$

$$(A^T[I+\beta L]A+\rho L)x=A^T[I+\beta L]y$$

$$A^{+}(Ax-y)+\rho L(g)x=0.$$

$$Sx-s+\rho L(g)x=0,$$

$$\underset{x\in {\mathbb{R}}^n}{\min }\Vert Sx-s{\Vert }_2^2+\rho x^{T}Lx,$$

$$\underset{x}{\min }{x}^{T}Lx\text{mbox},subjecttoSx-s=0,$$

$$(I-S)L(I-S)x=-(I-S)Ls,$$

$$(I-S)Lu=-(I-S)Ls,$$

$$x_0=s,x_{i+1}=s+(I-S)\text{mbox}BM3D(x_i),$$

$$(I-S)L(I-S)x=-(I-S)LSx,$$

$$S=\left[\begin{array}{cc}0& \\ & I\end{array}\right],L=\left[\begin{array}{cc}{L}_{11}& L_{12}\\ L_{21}& L_{22}\end{array}\right],x=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right],s=\left[\begin{array}{c}0\\ s_2\end{array}\right],$$

$$\left[\begin{array}{cc}\rho L_{11}& \rho L_{12}\\ \rho L_{21}& I+\rho L_{22}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}0\\ s_2\end{array}\right].$$

$$(I+\rho L/L_{11})x_2=s_2,\text{mbox}whereL/L_{11}=L_{22}-L_{21}{L}_{11}^{-1}{L}_{12},$$

Full text

I Introduction

Super-resolution (SR) refers to techniques that reconstruct high-resolution (HR) images from low-resolution (LR) images [1, 2, 3, 4, 5], determining high-frequency components and removing degradation caused by image acquisition in LR cameras. A single-image edge-preserving SR by interpolation, called the zooming problem in [6, 7], is not well-posed due to an ambiguity of determining high-frequency components, needed to preserve edges in the reconstructed high-resolution (HR) image. To remove the ambiguity and make the SR problem well-posed, one can introduce a guiding subspace, determined by frames, or via an action of an orthogonal projection on it; see, e.g., [8, 9].

Frame-less guided SR using an orthogonal projector on the guiding subspace is theoretically well justified, e.g., in [8, 9], but may be difficult in practice, where even the guiding subspace itself may be not explicitly known. We extend the approach of [8, 9] to guiding operators, which increase signal components in the guiding subspace, relative to those in a complementary subspace. As noticed in [9], the traditional Tikhonov’s regularization may be substituted with pre- and post-smoothing of sample-consistent reconstruction in case of noisy samples. We adopt a similar approach and provide alternative to [9] mathematical justification covering general guiding operators, not necessarily orthogonal projectors.

We illustrate our sample-consistent operator guided signal reconstruction with pre- and post-smoothing for imaging. An HR image is reconstructed from a noisy LR image using an HR image of the same scene but in different modality as a guidance in setting up an edge-preserving denoising filter. Several authors study a corresponding application, where the LR image is a depth image and the HR guidance image is an RGB image of the same scene; e.g., see [10, 11, 12, 13].

II Notation and prior work

A sampled (degraded) signal vector $y$ in signal processing is often represented by the linear model

[TABLE]

where the vector $x$ is the latent signal (image), $n$ is a noise vector, e.g., consisting of independent and identically distributed zero mean noise. The matrix $A$ is typically a product of down-sampling (e.g., decimation) and blurring operators.

The most frequently used restoration method is given by Tikhonov’s regularization model, see, e.g., [4, 7, 14],

[TABLE]

where the functional $R(x)$ is called a regularization term, and $\rho>0$ is a regularization parameter. The regularization term $R(x)$ aims at bounding and smoothing the solution $x$ of (2). The weighted norm $\|\cdot\|_{\omega}$ is defined by $\|x\|_{\omega}^{2}=x^{T}\omega x$, where $\omega$ is a positive definite matrix. If $A=I$ is the identity matrix, then (2) describes denoising.

Signal (image) processing often uses total variation due to its remarkable ability to preserve contours/edges of signals. The total variation term $R(x)$ approximates $\|\nabla x\|_{1}$; see, e.g., [7, 11, 14, 15]. Total variation filtering analogs can be set up in a framework of graph-based signal processing, e.g., by setting $R(x)=x^{T}L(g)x$, where $L(g)$ is a graph Laplacian matrix guided by a signal $g$, e.g., [16], that we define next.

A signal (image) is interpreted as an intensity function on $N$ vertices $V$ of a weighted graph $G=(V,E,W)$ consisting of a finite set $V$ of vertices (e.g., representing image pixels) and a finite set $E\subset V\times V$ of edges $(i,j)$ with typically nonnegative (cf., [17, 18]) weights $W(i,j)$, which measure similarity between vertices $i$ and $j$ in the graph. The intensity values of the signal form the vector $g=[g_{1},\ldots,g_{N}]^{T}$, where the vertices $V$ are arbitrarily numbered. The similarity weights form a symmetric $N\times N$ graph adjacency matrix $W=W(g)$.

Row-sums $d=W1_{N}$, where $1_{N}$ is the vector of ones, of $W$ define the diagonal $N\times N$ degree matrix $D={\rm diag}(d)$ of the graph. Desired for smoothing filters property $d=1_{N}$, i.e. $D=I$, of the graph adjacency matrix $W=W(g)$ can be ensured by scaling $W$ with positive diagonal matrix multipliers via Sinkhorn’s algorithm; e.g., [19]. Alternatively, $W$ can be substituted by the normalized graph adjacency matrix $D^{-1}W$, although the latter is technically non-symmetric.

The symmetric and positive semidefinite graph Laplacian matrix $L=L(g)$ is then defined as $L=D-W$. In graph-based signal processing, eigenvectors of $L$ serve as generalizations of basis functions of the Discrete Cosine Transform. If the guidance image $g$, used to determine $L(g)$, is aligned with an noisy image $x$ and shares the same edges, edge-preserving denoising of $x$ can be performed using filters based on $L(g)$. For example, the multiplication $W(g)x$ on $x$ amplifies spectral components in $x$ corresponding to large eigenvalues of $L$ and thus can be viewed as an approximate high-pass filter; cf. [20].

With $R(x)=x^{T}L(g)x$, the normal equations for the optimization problem (2) become

[TABLE]

where $L(g)$ is the graph Laplacian operator with a guidance signal $g$. A self-guided choice $g=x$ is nonlinear; e.g., [16].

Graph-based interpretations of denoising filters are common. Classical bilateral and total variation filters for image denoising are usually constructed via explicit formulas for the weights $W(i,j)$ on a priori determined graph edges $E$; e.g., [12, 16]. More recent guided filter [20] is defined directly via explicit description of its action $W(g)x$ on a given image $x$, while its edges $E$ and formulas for the weights $W(i,j)$ are then derived for theoretical purposes only. Sophisticated filters, e.g., BM3D [21], are determined exclusively by functions that implement their action, so it may be difficult to obtain explicit formulas for their weights $W(i,j)$. Exact low-pass filters, where orthogonal projectors represent $W$ and $W1_{N}=1_{N}$, also fit the graph-based framework, although some weights $W(i,j)$ are negative.

System (3) with the weight matrix $\omega=I+\beta L(g)$, where $\beta\geq-1$ is a parameter, is used in [19]. The resulting system

[TABLE]

is solved by the conjugate gradient (CG) method; see [19].

Connecting Tikhonov’s regularization equation (2) to frame-based reconstruction, one can set $\omega=\left(AA^{T}\right)^{+}$, where the operation + denotes the Moore-Penrose pseudoinverse, which gives $A^{+}=A^{T}\left(AA^{T}\right)^{+}$ and turns (3) into

[TABLE]

Denoting $s=A^{+}y$ and introducing the orthogonal projector $S=A^{+}A$, we equivalently rewrite equation (4) as

[TABLE]

which is the normal equation for Tikhonov’s regularization

[TABLE]

—a particular case of optimization problem (2).

The authors of [8, 9] investigate the case of problem (6), where $L=I-W$ is an orthogonal projector. In particular, they prove that when the regularization parameter $\rho$ vanishes unconstrained minimization (6) reduces to

[TABLE]

which can be interpreted as graph-harmonic sample-consistent signal extension, since the quadratic form $x^{T}Lx$ can be viewed as energy of the signal $x$, defined by the graph Laplacian $L$.

Since $S$ is an orthogonal projector, we have $Sx=s=Ss$ and the orthogonal decomposition $x=s+(I-S)x$ helps to show that the minimizer $x$ in (7) solves the system

[TABLE]

where the matrix $(I-S)L(I-S)$ is symmetric positive semidefinite. The special structure of equation (8) allows applying CG method to the linear system

[TABLE]

with an initial approximation from the null-space $\mbox{Null}(S)$ of $S$ to iteratively approximate the solution $u$ within the subspace $\mbox{Null}(S)$. The sample-consistent reconstruction $x$ in (7) is then given by the orthogonal sum $x=s+u$; see [8, 9] for details.

The authors of [21] propose solving a self-guided nonlinear version of (7) via simple iteration

[TABLE]

where $\mbox{BM3D}(x_{i})$ is an application of the BM3D filter to $x_{i}$.

Post-processing by $\alpha x+(1-\alpha)Wx$ with $\alpha=1/(\rho+1)$ is proved in [9] to solve Tikhonov’s regularization problem (6). Since $\alpha\approx 1-\rho$ for small $\rho$, the post-processing is approximated by $x-\rho Lx$ when $\rho\to 0$. The corresponding arguments in [9] rely on the assumptions that $L=I-W$ and that $W$ is an orthogonal projector. We provide an alternative analysis, dropping these assumptions, in the next section.

III Tikhonov’s regularization demystified

Equation (5) implies $(I-S)Lx=0$, i.e.

[TABLE]

which differs from (8) only in the right-hand side, since the sample consistency $Sx=s$ is not enforced in (5).

In an orthonormal basis of $\mathbb{R}^{N}$, where

[TABLE]

equation (5) takes the following block form,

[TABLE]

Assuming that the block $L_{11}$ is invertible, we solve the equivalent to (11) equation $L_{11}\,x_{1}=-L_{12}\,x_{2}$ for $x_{1}=-L_{11}^{-1}L_{12}\,x_{2}$ and then for $x_{2}$ we obtain the following equation

[TABLE]

is the Schur complement of the block $L_{11}$ in the matrix $L$.

If $\rho\to 0$, we conclude that $x_{2}\to s_{2}$ and $x_{1}\to-L_{11}^{-1}L_{12}s_{2}$, so that the limit vector $x$ solves (7) and (8), where $x_{2}=s_{2}$, i.e. $Sx=s$. For small $\rho$, we have $x_{2}\approx s_{2}-\rho L/L_{11}\,s_{2}$, i.e. the solution $x$ of Tikhonov’s regularization equation (5) depends linearly on $\rho$. The minuend $L_{22}\,s_{2}$ in $L/L_{11}\,s_{2}$ corresponds to the vector $SLs$. Thus, rather than solving (5) directly, for small $\rho$ one can try sample-consistent reconstruction via solving (7) or (8), but with an a priori denoised sampled signal $s-\rho SLs$.

IV The problem and proposed algorithms

We are given two images of the same scene, but in different modalities. The first image $y$ has low resolution and may be noisy. The second image $g$ is of high resolution and is noise-free. The images are registered (or aligned) by means of a linear downsampling transform, that is, there is a well-conditioned matrix $A$ such that the images $Ag$ and $y$ are aligned. We want to reconstruct an image $x$ of the same resolution as $g$ from the image $y$ so that the image $g$ serves as a guidance image in the regularization term of the Tikhonov’s regularization reconstruction model (4).

In contrast to the traditional approach directly solving (4), we use our arguments above to take the sample-consistent reconstruction approach (7) from [8] as a starting point, but for a more general case where the Laplacian $L$ is an arbitrary positive semi-definite operator, not necessarily an orthogonal projector. We notice that the proposed in [8, 9] reduction of (7) to (9) also works for this general case, so we solve (9).

Noisy LR samples $y$ should evidently be smoothed prior to computing the sample-consistent HR reconstruction $x$ by (9). Dealing with noise, we always pre-process $y$ via denoising, by analogy with $y-\rho ALA^{+}y$ proposed in Section III. We also find it helpful in some tests to post-process our sample-consistent reconstruction $x$, similar to subtracting $x-\rho Lx$ highly-oscillatory contributions described in Section II.

Let $CG_{m}(\mathcal{A},v)$ denote a function, which implements $m$ iterations of CG to solve the equation $\mathcal{A}(u)=v$ with a given linear operator $\mathcal{A}(u)$. The main cost per iteration is the cost of evaluation of $\mathcal{A}(u)$. Our reconstruction method is as follows.

[TABLE]

The level $x^{T}1_{N}$ of the DC-component $1_{N}$ in the reconstructed HR signal $x$ can be adjusted during post-processing to match that in the LR sample $y$, if necessary.

Compared to [8], where $W$ and $L=I-W$ are assumed to be orthogonal projectors, our approach allows choosing more general smoothing filters $W$, but there are still limitations:

• •

Negative entries in the filter matrix $W$ are allowed, but the resulting Laplacian $L=D-W$ needs to be symmetric positive semidefinite, in order for the minimization of its quadratic form $x^{T}Lx$ in (7) to make sense, and to satisfy assumptions for CG convergence in Algorithm 1.

• •

DC-invariant filters, i.e. satisfying $W1_{N}=1_{N}$, lead to $D=I$ and thus to the trivial construction of $L=I-W$, such as, e.g., the image guided filter of [20], implemented in the MATLAB Image Processing Toolbox, and the total variation filter presented in [16]. Otherwise, the filter DC-component evaluation $d=W1_{N}$ is needed to determine $D=\mbox{diag}(d)$ in $L=D-W$.

• •

Iterative filters are allowed, e.g., a DC-invariant smoothing filter $W$, satisfying $W1_{N}=1_{N}$, generates the DC-invariant polynomial filter $P_{n}(I-W)$ where $P_{n}(\cdot)$ is a polynomial of degree $n$, satisfying $P_{n}(1)=1$, for example, $P_{2}(W)=W^{2}$. However, the polynomial $P_{n}(\cdot)$ needs to remain fixed during the CG iterations.

• •

Self-guided filters lead to nonlinear operators $\mathcal{A}(\cdot)$, requiring special care, such as proposed in [16].

V Numerical experiments

Super-resolution for flash-no flash images

We have carried out numerical tests in MATLAB with two registered color images from the MATLAB own directory, toysflash.png and toysnoflash.png of the size $684\times 912$, taken respectively with and without flash. The first image is our HR guidance image. We downsample the second image by factor $4$ in both dimensions choosing every $4$-th pixel with the command image(1:4:end,1:4:end). The result of downsampling is our LR image. Both images are scaled so that their intensities lie in the range $[0,1]$. Figures 1 and 2 display the guidance and the ground truth images. The two images of the same scene have different modalities owing to big differences in color representations of the objects. For example, the yellow ball looks white in the flash image.

To produce a noisy downsampled LR image, a Gaussian noise with the default parameters, zero mean and variance $0.01$, has been added to the LR image. The degraded image is displayed in Figure 3. We pre-smooth the noisy LR image by the image guided filter with the default parameters and then apply Algorithm 1 to the pre-smoothed image.

We use a single application of the guided filter function imguidedfilter from the MATLAB Image Processing Toolbox as the smoothing filter $W(g)$ for color images. The function parameters are as follows: width of neighborhoods is $7$, the smoothing value is $10^{-6}$.

The relative residual is $10^{-3}$ after $m=20$ CG iterations, which translates into $20$ applications of imguidedfilter. Figure 4 shows the reconstructed image, which has the same desired colors as the ground truth image. The PSNR referred to the ground truth HR image equals $24.32$.

Super-resolution of a depth image

Algorithm 1 has been also evaluated on the Art image set from the simulated Middlebury 2007 data sets extensively used in [11]. The guidance image is the gray component of the HR RGB image Art in Figure 5. The LR noisy depth image is shown in Figure 6. The upsampling shown in Figure 7 has been computed by the code used in [11], which implements Tikhonov’s regularization by the generalized total variation.

To produce Figure 8 we use the guided Total Variation (TV) filter described in [16] as the guided smoothing filter in Algorithm 1 as well as during LR pre-smoothing and HR post-smoothing. The smoothing parameter used is $10^{-4}$ at all three stages, i.e. inside Algorithm 1, as well as for the pre- (and post-) smoothing, with 1/480/850 function evaluations, correspondingly. The number $m$ of CG iterations is $2900$ to achieve the relative residual $2\cdot 10^{-14}$.

Comparison Figure 7 has a smaller PSNR, but a bit less noisy and has sharper edges, than our Figure 8, as the generalized TV filter used in [11] is more powerful smoother compared to the TV filter from [16] that we use in these tests. However, Figure 7 has noticeable artifacts, clearly coming from sharp edges in RGB Figure 5, which are absent in LR Figure 6. In our approach, sharp edges in the guidance image $g$ only affect the filter weights and thus cannot possibly pollute our reconstruction, as confirmed in Figure 8.

VI Conclusion

We propose algorithms for guided iterative reconstruction of signals, illustrated with reconstruction of a higher resolution image from a lower resolution noisy sample. Guiding is given by a smoothing filter and plays the role of regularization. Least squares minimization of the data term $\|Ax-y\|_{2}$ is substituted by the equality constraint $Ax=y$, allowing one to control sample consistency and eliminating the problem of choosing the regularization parameter in Tikhonov’s regularization.

Noisy lower resolution samples are pre-smoothed. The obtained sample-consistent high resolution reconstruction can be post-smoothed, if needed. Iterations are based on conjugate gradients, giving the optimal performance. General smoothing guidance filters allow flexibility in designing reconstructions with desirable properties. Initial numerical experiments demonstrate feasibility of the proposed technology for super-resolution for flash-no flash images and for depth reconstruction guided by RGB images.