Multichannel reconstruction from nonuniform samples with application to   image recovery

Dong Cheng; Kit Ian Kou

arXiv:1901.01102·math.CA·December 20, 2024

Multichannel reconstruction from nonuniform samples with application to image recovery

Dong Cheng, Kit Ian Kou

PDF

TL;DR

This paper develops a theory for multichannel interpolation from nonuniform samples, enabling improved image recovery and interpolation accuracy, especially in highly corrupted images, by leveraging multiple sampling models and derivatives.

Contribution

It introduces a novel multichannel interpolation framework for nonuniform samples, extending previous uniform sampling methods and applying it to enhance image reconstruction quality.

Findings

01

Effective interpolation from nonuniform samples demonstrated

02

Superior image recovery compared to median filter methods

03

Lower errors in highly corrupted image reconstruction

Abstract

The multichannel trigonometric reconstruction from uniform samples was proposed recently. It not only makes use of multichannel information about the signal but is also capable to generate various kinds of interpolation formulas according to the types and amounts of the collected samples. The paper presents the theory of multichannel interpolation from nonuniform samples. Two distinct models of nonuniform sampling patterns are considered, namely recurrent and generic nonuniform sampling. Each model involves two types of samples: nonuniform samples of the observed signal and its derivatives. Numerical examples and quantitative error analysis are provided to demonstrate the effectiveness of the proposed algorithms. Additionally, the proposed algorithm for recovering highly corrupted images is also investigated. In comparison with the median filter and correction operation treatment, our…

Figures9

Click any figure to enlarge with its caption.

Tables3

Table 1. Table 1: Reconstruction results using the different interpolation formulas

Total samples	$f$	$f^{'}$	Pattern	$δ_{1}$ (Variance)	$δ_{2}$ (Variance)
$36$	$36$	$0$	RN1	$0.8560$	$0.8358$
$36$	$36$	$0$	GN1	$0.5548$ ( $0.0026$ )	$0.5535$ ( $0.0026$ )
$36$	$36$	$0$	U1	$0.5120$	$0.5116$
$36$	$18$	$18$	RN2	$0.6163$	$0.6159$
$36$	$18$	$18$	GN2	$1.0752$ ( $0.0346$ )	$1.0550$ ( $0.0361$ )
$36$	$18$	$18$	U2	$0.9241$	$0.8381$
$54$	$54$	$0$	RN1	$0.1955$	$0.1922$
$54$	$54$	$0$	GN1	$0.1501$ ( $1.11 \times 10^{- 4}$ )	$0.1498$ ( $1.09 \times 10^{- 4}$ )
$54$	$54$	$0$	U1	$0.1376$	$0.1376$
$54$	$27$	$27$	RN2	$0.1830$	$0.1830$
$54$	$27$	$27$	GN2	$0.2829$ ( $3.4 \times 10^{- 3}$ )	$0.2794$ ( $3.4 \times 10^{- 3}$ )
$54$	$27$	$27$	U2	$0.2582$	$0.2483$
$72$	$72$	$0$	RN1	$0.0437$	$0.0426$
$72$	$72$	$0$	GN1	$0.0321$ ( $3.77 \times 10^{- 6}$ )	$0.0320$ ( $3.73 \times 10^{- 6}$ )
$72$	$72$	$0$	U1	$0.0291$	$0.0291$
$72$	$36$	$36$	RN2	$0.0355$	$0.0355$
$72$	$36$	$36$	GN2	$0.0647$ ( $1.18 \times 10^{- 4}$ )	$0.0636$ ( $1.22 \times 10^{- 4}$ )
$72$	$36$	$36$	U2	$0.0557$	$0.0520$
$108$	$108$	$0$	RN1	$0.0018$	$0.0017$
$108$	$108$	$0$	GN1	$0.0013$ ( $4.56 \times 10^{- 9}$ )	$0.0013$ ( $4.50 \times 10^{- 9}$ )
$108$	$108$	$0$	U1	$0.0012$	$0.0012$
$108$	$54$	$54$	RN2	$0.0014$	$0.0014$
$108$	$54$	$54$	GN2	$0.0028$ ( $1.81 \times 10^{- 7}$ )	$0.0028$ ( $1.83 \times 10^{- 7}$ )
$108$	$54$	$54$	U2	$0.0023$	$0.0021$

Table 2. Table 2: Comparison of several existing interpolation methods.

Different methods	Untruncated implementation	Applicable to nonuniform samples	Applicable to multichannel samples	Closed form of interpolating functions
Proposed method	yes	yes	yes	yes
Single-channel interpolation by FFT [25]	yes	yes¹	N/A	N/A
GSE [30, 14]	N/A	yes²	yes	yes
Classical nonuniform sampling on real line [1, 9]	N/A	yes	N/A	N/A
Single-channel nonuniform trigonometric interpolation [16]	yes	yes	N/A	yes

Table 3. Table 3: Quantitative measurements for quality of image recovery results.

	GN1 + CRT	GN2 + CRT	MED + CRT	CRT +MED
RMSE: $δ$	0.0570	0.0488	0.0884	0.0875
PSNR: $ρ$	30.56	31.91	26.76	26.84
CC : $γ$	0.9874	0.9908	0.9719	0.9726

Equations226

f (t) = n = - \infty \sum \infty f (t_{n}) S_{n} (t)

f (t) = n = - \infty \sum \infty f (t_{n}) S_{n} (t)

S_{n} (t) = \frac{G ( t )}{G ^{'} ( t _{n} ) ( t - t _{n} )}, G (t) = (t - t_{0}) k = 1 \prod \infty (1 - \frac{t}{t _{k}}) (1 - \frac{t}{t _{- k}})

S_{n} (t) = \frac{G ( t )}{G ^{'} ( t _{n} ) ( t - t _{n} )}, G (t) = (t - t_{0}) k = 1 \prod \infty (1 - \frac{t}{t _{k}}) (1 - \frac{t}{t _{- k}})

f (t) = \frac{1}{2 M + 1} k = 0 \sum 2 M f (\frac{k}{2 M + 1}) \frac{( - 1 ) ^{k} sin π ( 2 M + 1 ) t}{sin π ( t - \frac{k}{2 M + 1} )} .

f (t) = \frac{1}{2 M + 1} k = 0 \sum 2 M f (\frac{k}{2 M + 1}) \frac{( - 1 ) ^{k} sin π ( 2 M + 1 ) t}{sin π ( t - \frac{k}{2 M + 1} )} .

(f, h) := \frac{1}{2 π} \int_{T} f (t) \overline{h (t)} d t, \forall f, h \in L^{2} (T) .

(f, h) := \frac{1}{2 π} \int_{T} f (t) \overline{h (t)} d t, \forall f, h \in L^{2} (T) .

f (t) = n \in Z \sum a (n) e^{i n t}

f (t) = n \in Z \sum a (n) e^{i n t}

(f, h) = n \in Z \sum a (n) \overline{b (n)},

(f, h) = n \in Z \sum a (n) \overline{b (n)},

(f * h) (t) := \frac{1}{2 π} \int_{T} f (s) h (t - s) d s = n \in Z \sum a (n) b (n) e^{i n t} .

(f * h) (t) := \frac{1}{2 π} \int_{T} f (s) h (t - s) d s = n \in Z \sum a (n) b (n) e^{i n t} .

H f (t) := \frac{1}{2 π} p.v. \int_{T} f (s) cot (\frac{t - s}{2}) d s = n \in Z \sum (- i \mbox s g n (n)) a (n) e^{i n t}

H f (t) := \frac{1}{2 π} p.v. \int_{T} f (s) cot (\frac{t - s}{2}) d s = n \in Z \sum (- i \mbox s g n (n)) a (n) e^{i n t}

f (t) = n \in I^{N} \sum a (n) e^{i n t}, I^{N} = {n : N_{1} \leq n \leq N_{2}} .

f (t) = n \in I^{N} \sum a (n) e^{i n t}, I^{N} = {n : N_{1} \leq n \leq N_{2}} .

h_{m} (t)

h_{m} (t)

g_{m} (t)

I_{j} = {n : N_{1} + (j - 1) K \leq n \leq N_{1} + j K - 1} .

I_{j} = {n : N_{1} + (j - 1) K \leq n \leq N_{1} + j K - 1} .

f (t) = \frac{1}{K} m = 1 \sum M p = 0 \sum K - 1 g_{m} (\frac{2 π p}{K}) y_{m} (t - \frac{2 π p}{K})

f (t) = \frac{1}{K} m = 1 \sum M p = 0 \sum K - 1 g_{m} (\frac{2 π p}{K}) y_{m} (t - \frac{2 π p}{K})

H_{n}^{- 1} = q_{11} (n) q_{21} (n) ⋮ q_{M 1} (n) q_{12} (n) q_{22} (n) ⋮ q_{M 2} (n) \dots \dots \dots q_{1 M} (n) q_{2 M} (n) ⋮ q_{M M} (n) .

H_{n}^{- 1} = q_{11} (n) q_{21} (n) ⋮ q_{M 1} (n) q_{12} (n) q_{22} (n) ⋮ q_{M 2} (n) \dots \dots \dots q_{1 M} (n) q_{2 M} (n) ⋮ q_{M M} (n) .

y_{m} (t) = n \in I^{N} \sum r_{m} (n) e^{i n t}

y_{m} (t) = n \in I^{N} \sum r_{m} (n) e^{i n t}

r_{m} (n) = {q_{mj} (n + K - j K), 0 if n \in I_{j}, j = 1, 2, \dots, M, if n \in / I^{N} .

r_{m} (n) = {q_{mj} (n + K - j K), 0 if n \in I_{j}, j = 1, 2, \dots, M, if n \in / I^{N} .

T_{1} (f, 2 m_{0}, α, t) = p = 0 \sum m_{0} - 1 f (t_{p}) y_{1, α} (t - t_{p}) + f (α + t_{p}) y_{2, α} (t - t_{p}) .

T_{1} (f, 2 m_{0}, α, t) = p = 0 \sum m_{0} - 1 f (t_{p}) y_{1, α} (t - t_{p}) + f (α + t_{p}) y_{2, α} (t - t_{p}) .

y_{1, α} (t)

y_{1, α} (t)

y_{2, α} (t)

y_{1, α} (t_{q} - t_{p})

y_{1, α} (t_{q} - t_{p})

y_{2, α} (t_{q} - t_{p})

f (t) = 0.05 t (t - 2 π) (0.04 t^{2} + 0.02 t^{3} + cos (3 sin t)), t \in [0, 2 π) .

f (t) = 0.05 t (t - 2 π) (0.04 t^{2} + 0.02 t^{3} + cos (3 sin t)), t \in [0, 2 π) .

H_{n} = [e^{i n α} e^{i (n + m_{0}) α} i n i (n + m_{0})] for n \in I_{1} .

H_{n} = [e^{i n α} e^{i (n + m_{0}) α} i n i (n + m_{0})] for n \in I_{1} .

H_{n}^{- 1} = [\frac{e ^{- i n α} ( m _{0} + n )}{m _{0} + n - n e ^{i m_{0} α}} \frac{i e ^{i m_{0} α}}{m _{0} + n - n e ^{i m_{0} α}} - \frac{n e ^{- i n α}}{m _{0} + n - n e ^{i m_{0} α}} - \frac{i}{m _{0} + n - n e ^{i m_{0} α}}] .

H_{n}^{- 1} = [\frac{e ^{- i n α} ( m _{0} + n )}{m _{0} + n - n e ^{i m_{0} α}} \frac{i e ^{i m_{0} α}}{m _{0} + n - n e ^{i m_{0} α}} - \frac{n e ^{- i n α}}{m _{0} + n - n e ^{i m_{0} α}} - \frac{i}{m _{0} + n - n e ^{i m_{0} α}}] .

v_{1, n, α} (t)

v_{1, n, α} (t)

v_{2, n, α} (t)

y_{k, α} (t - t_{p}) = \frac{1}{m _{0}} n = N_{1} \sum N_{1} + m_{0} - 1 v_{k, n, α} (t) e^{- i n \frac{2 π p}{m _{0}}}, k = 1, 2.

y_{k, α} (t - t_{p}) = \frac{1}{m _{0}} n = N_{1} \sum N_{1} + m_{0} - 1 v_{k, n, α} (t) e^{- i n \frac{2 π p}{m _{0}}}, k = 1, 2.

T_{2} (f, 2 m_{0}, α, t) := p = 0 \sum m_{0} - 1 f (α + t_{p}) y_{1, α} (t - t_{p}) + f^{'} (t_{p}) y_{2, α} (t - t_{p}),

T_{2} (f, 2 m_{0}, α, t) := p = 0 \sum m_{0} - 1 f (α + t_{p}) y_{1, α} (t - t_{p}) + f^{'} (t_{p}) y_{2, α} (t - t_{p}),

T_{2} (f, 2 m_{0}, 0, t) p = 0 \sum m_{0} - 1 f (t_{p}) y_{1, 0} (t - t_{p}) + f^{'} (t_{p}) y_{2, 0} (t - t_{p}),

T_{2} (f, 2 m_{0}, 0, t) p = 0 \sum m_{0} - 1 f (t_{p}) y_{1, 0} (t - t_{p}) + f^{'} (t_{p}) y_{2, 0} (t - t_{p}),

y_{1, 0} (t)

y_{1, 0} (t)

y_{2, 0} (t)

H = b_{1} (N_{1}) b_{1} (N_{1} + 1) ⋮ b_{1} (N_{1} + M - 1) b_{2} (N_{1}) b_{2} (N_{1} + 1) ⋮ b_{2} (N_{1} + M - 1) \dots \dots ⋱ \dots b_{M} (N_{1}) b_{M} (N_{1} + 1) ⋮ b_{M} (N_{1} + M - 1) .

H = b_{1} (N_{1}) b_{1} (N_{1} + 1) ⋮ b_{1} (N_{1} + M - 1) b_{2} (N_{1}) b_{2} (N_{1} + 1) ⋮ b_{2} (N_{1} + M - 1) \dots \dots ⋱ \dots b_{M} (N_{1}) b_{M} (N_{1} + 1) ⋮ b_{M} (N_{1} + M - 1) .

H = e^{i N_{1} t_{1}} e^{i (N_{1} + 1) t_{1}} ⋮ e^{i (N_{1} + M - 1) t_{1}} e^{i N_{1} t_{2}} e^{i (N_{1} + 1) t_{2}} ⋮ e^{i (N_{1} + M - 1) t_{2}} \dots \dots ⋱ \dots e^{i N_{1} t_{M}} e^{i (N_{1} + 1) t_{M}} ⋮ e^{i (N_{1} + M - 1) t_{M}} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Multichannel reconstruction from nonuniform samples with application to image recovery

Dong Cheng [email protected] Department of Mathematics, Faculty of Science and Technology, University of Macau, Macao, China

Kit Ian Kou [email protected] Department of Mathematics, Faculty of Science and Technology, University of Macau, Macao, China

Abstract

The multichannel trigonometric reconstruction from uniform samples was proposed recently. It not only makes use of multichannel information about the signal but is also capable to generate various kinds of interpolation formulas according to the types and amounts of the collected samples. The paper presents the theory of multichannel interpolation from nonuniform samples. Two distinct models of nonuniform sampling patterns are considered, namely recurrent and generic nonuniform sampling. Each model involves two types of samples: nonuniform samples of the observed signal and its derivatives. Numerical examples and quantitative error analysis are provided to demonstrate the effectiveness of the proposed algorithms. Additionally, the proposed algorithm for recovering highly corrupted images is also investigated. In comparison with the median filter and correction operation treatment, our approach produces superior results with lower errors.

Keywords: Interpolation, nonuniform sampling, FFT, trigonometric polynomial, error analysis, derivative, image recovery.

Mathematics Subject Classification (2010): 42A15, 94A12, 65T50, 94A08.

1 Introduction

Sampling and reconstruction are used as fundamental tools in data processing and communication systems. Classical uniform sampling [1] is an effective tool to recover signals and has been applied in many applications. However there are various instances that reconstruction of signals from their nonuniform samples are required, such as computed tomography [2], magnetic resonance [3] and radio astronomy [4]. Numerous approaches have been proposed in the literature to reconstruct bandlimited signals from nonuniform samples [5, 6, 7, 8]. A widely known nonuniform sampling theorem [1, 9] may be stated as follows. Let $\{t_{n}\}_{n\in\mathbb{Z}}$ be a sequence of real numbers such that $\left|t_{n}-\frac{n\pi}{\sigma}\right|<\frac{\pi}{4\sigma}$ , then a $\sigma$ -bandlimited function can be reconstructed by

[TABLE]

where

[TABLE]

and (1.1) converges uniformly on any compact subset of $\mathbb{R}$ . The series (1.1) is an extension of Lagrange interpolation. Unlike the uniform sampling case, (1.1) contains infinitely many terms and the interpolating functions $S_{n}(t)$ involve complicated components. These factors bring difficulties for exact reconstruction of a bandlimited function from nonuniform samples. To give a simpler approximation for reconstructing a bandlimited function from nonuniform samples, the authors in [10] proposed a new kind of sinc interpolation method and they restricted $S_{n}(t)$ in (1.1) to be of the form $\mathrm{sinc}[\sigma(t-\tilde{t}_{n})]$ , where $\tilde{t}_{n}=nT+\zeta_{n}$ and $\zeta_{n}$ is a sequence of random variables independent of $G(t)$ . This restriction guarantees that the interpolating functions only consist of translation of sinc function, just like most cases of uniform interpolation [11, 12, 13, 14]. However, the restriction strategy simplifies reconstruction problem but introduces error inevitably. To overcome the error, several methods for determining suitable $\tilde{t}_{n}$ were analyzed [10]. By the similar idea, the sinc interpolation for nonuniform samples in fractional Fourier domain was studied in [15].

In a real application, there are only finitely many samples, albeit with large amount, are given in a bounded region. Interpolating by sinc functions or any other bandlimited functions has some limitations. On the one hand, the bandlimited interpolating functions cannot be time limited by the uncertainty principle , thereby the approximation error is introduced. On the other hand, the interpolating functions $S_{n}(t)$ for nonuniform samples are complicated and there is no closed form in general. Therefore, treating a finite amount of data as samples of a periodic function is a convenient and feasible choice [16]. As we know, sines and cosines are classical periodic functions, they have wide applications in modern science. It is no exaggeration to say that trigonometry pervades the area of signal processing. We know that $\{e^{\bm{i}nt}:n\in\mathbb{Z}\}$ is an orthogonal system and is complete in square integrable functions space on unit circle $\mathbb{T}$ , i.e., ${L}^{2}(\mathbb{T})$ . Besides, these functions possess elegant symmetries and concise frequency meanings. These desirable properties of trigonometric functions could make interpolation much simpler and more effective [16, 17]. In fact, in the early 1841, Cauchy first proved a sampling interpolation theorem on trigonometric polynomials [18]. It may be recognized as the headstream of sampling theory [19]. Cauchy’s result states that if $f(t)=\sum_{\left|n\right|\leq M}c_{n}e^{2\pi\bm{i}tn}$ , then it can be written as a sum of its sampled values $f(\frac{k}{2M+1})$ , $0\leq k\leq 2M$ , each multiplied by a interpolating function. That is,

[TABLE]

Certain studies have been given to the problem of interpolating finite length samples by trigonometric functions or discrete Fourier transform. In a series of papers [20, 21, 22], the sinc interpolation of discrete periodic signals were extensively discussed. Although referred to as sinc interpolation, the resulting interpolating functions are trigonometric. In [23], the authors decomposed a periodic signal in a basis of shifted and scaled versions of a generating function. Moreover, an error analysis for the approximation method was also addressed. A generalized trigonometric interpolation was considered in [17] to make a good approximation for non-smooth functions. Recently, the nonuniform sampling theorems for trigonometric polynomials were presented [16, 24]. Selva [25] proposed a FFT-based interpolation of nonuniform samples. However, this method is valid only for nonuniform samples lying in a regular grid rather than for non-uniformly distributed data in the general sense.

In all of above mentioned interpolation methods for finite length discrete points, only the samples of original function are processed. As an extension of trigonometric interpolation, a multichannel interpolation of finite length samples was suggested in [26]. This novel method makes good use of multifaceted information (such as derivatives, Hilbert transform) of function and is capable of generating various useful interpolation formulas by selecting suitable parameters according to the types and amount of collected data. In addition, it can be used to approximate some integral transformations (such as Hilbert transform). A fast algorithm based on FFT makes multichannel interpolation more effective and stable. However, only the cases of uniform sampling were considered in [26]. There is a need to extend multichannel interpolation such that non-uniformly distributed data can be processed.

The purpose of this paper is to establish the theory of multichannel interpolation for non-uniformly distributed data. We will consider two kinds of nonuniform sampling patterns: recurrent and generic nonuniform sampling. Meanwhile, each kind of nonuniform sampling involves two types of samples: its own nonuniform samples and derivative’s samples. There are four nonuniform interpolation formulas will be analyzed. All closed-form expressions of interpolating functions are derived. Some examples are also demonstrated. We show that the trigonometric polynomial (also called periodic bandlimited function) of finite order can be exactly reconstructed by the proposed interpolation formulas provided that the total number of samples is enough. Error analysis of the reconstruction for non-bandlimited square integrable functions are analyzed. Concretely, the contributions of this paper may be summarized as follows:

We propose four types of interpolation formulas for non-uniformly distributed data. The proposed formulas involves not only samples of $f$ but also samples of $f^{\prime}$ , where $f$ is the function to be reconstructed. If the given data is sampled from a periodic bandlimited function, then we arrive at a perfect reconstruction provided that the amount of data is larger than the bandwidth. 2. 2.

We analyze the error that arise in reconstructing a non-bandlimited function by the proposed formulas. In particular, a comparison of performance on reconstructing square integrable functions (not necessarily to be bandlimited) by these formulas is made. 3. 3.

Applying the proposed interpolation formulas, we develop algorithms for the recovery of damaged pixels which are non-uniformly located in a degraded image. The algorithms perform well and can be efficiently implemented. Thus they could be good pre-processing methods for some more sophisticated approaches (such as deep learning) in the image recovery problem.

This paper is organized as follows. In Section 2 some preparatory knowledge of Fourier series and multichannel interpolation are reviewed. Section 3 and 4 formulate four types of interpolation formulas for non-uniformly distributed data. The numerical examples and error analysis are presented in Section 5. The application of proposed interpolation method to image recovery is shown in Section 6. Finally, conclusion will be drawn in Section 7.

2 Preliminaries

2.1 Fourier series

Without loss of generality, we will consider the functions defined on unit circle $\mathbb{T}$ . Let $L^{2}(\mathbb{T})$ be the totality of square integral functions defined on $\mathbb{T}$ . It is known that $L^{2}(\mathbb{T})$ is a Hilbert space embedded with the inner product

[TABLE]

For $f\in L^{2}(\mathbb{T})$ , it can be expanded as

[TABLE]

where the Fourier series is convergent to $f$ in $L^{2}$ norm. The general version of Parseval’s identity is of the form

[TABLE]

where $\{a(n)\}$ and $\{b(n)\}$ are Fourier coefficients of $f$ and $h$ respectively. The convolution theorem manifests as

[TABLE]

The circular Hilbert transform [27, 28] is an useful tool in harmonic analysis and signal processing. It is defined by the singular integral

[TABLE]

where sgn is the signum function taking values $1$ , $-1$ or [math] for $n>0$ , $n<0$ or $n=0$ respectively. From the definition, we see that it is simple and straightforward to compute Hilbert transform for trigonometric functions. Thus trigonometry-based interpolation can be availably used to approximate Hilbert transform as well.

2.2 Multichannel interpolation

Multichannel interpolation proposed in [26] is about the reconstruction problem of finite order trigonometric polynomials. To maintain consistent terminology with the classical case, in what follows, a finite order trigonometric polynomial is called a periodic bandlimited function, or briefly a bandlimited function. Let $\mathbf{N}=(N_{1},N_{2})\in\mathbb{Z}^{2}$ , in the sequel we denote by $B_{\mathbf{N}}$ the totality of bandlimited functions with the following form:

[TABLE]

The bandwidth of $f$ is defined by the cardinality of $I^{\mathbf{N}}$ , denoted by $\mu(I^{\mathbf{N}})$ .

For $1\leq m\leq M$ , let

[TABLE]

Suppose that $\frac{N_{2}-N_{1}+1}{M}=K\in\mathbb{N}^{+}$ , we cut $I^{\mathbf{N}}$ into pieces as $I^{\mathbf{N}}=\bigcup_{j=1}^{M}I_{j}$ , where

[TABLE]

The multichannel interpolation indicates that a bandlimited function $f\in B_{\mathbf{N}}$ can be reconstructed by samples of $g_{m}$ , namely,

[TABLE]

provided that $M\times M$ matrix $\mathbf{H}_{n}=\left[b_{m}(n+jK-K)\right]_{jm}$ is invertible for every $n\in I_{1}$ . Here, the interpolating functions are constructed by the elements of $\mathbf{H}_{n}^{-1}$ . We denote the inverse matrix as

[TABLE]

The interpolating function $y_{m}$ for $1\leq m\leq M$ is given by

[TABLE]

where

[TABLE]

In the following sections, using the powerful technique of multichannel interpolation, we present four types of nonuniform interpolation formulas. Since there are some similar concepts involved in the following parts, several notations may appear repeatedly with minor difference. We particularly remark that a notation could have different meanings across different parts.

3 Multichannel interpolation of recurrent non-uniformly distributed data

The recurrent nonuniform sampling often arises in time-interleaved analog-to digital converting process [29, 10]. As for recurrent nonuniform sampling, a classical result that has to be mentioned is the Papoulis’ generalized sampling expansion (GSE) [30]. The differences between GSE and the multichannel interpolation are mainly as follows:

•

The GSE involves infinite summation and is applied to recovering functions defined on whole real line. Therefore the truncation is inevitable in practice. The multichannel interpolation is about reconstructing a finite length function from a finite number of samples.

•

There is a FFT-based fast algorithm to implement the multichannel interpolation. The implementation of GSE is more complicated.

•

The multichannel interpolation can be extended to the generic nonuniform sampling case (see Section 4). However, to the authors’ knowledge, there is no generic nonuniform sampling formula based on GSE.

In the following two subsections, based on multichannel interpolation technique, we derive two interpolation formulas associated with recurrent nonuniform samples: one concerns derivative of function and the other does not. Throughout Section 3, let $m_{0}\in\mathbb{N}^{+}$ and $t_{p}=\frac{2\pi p}{m_{0}}$ for $p=0,1,\dots,m_{0}-1$ .

3.1 Recurrent nonuniform samples

By setting $b_{1}(n)=1$ , $b_{2}(n)=e^{\bm{i}n\alpha}$ with $0<\alpha<\frac{2\pi}{m_{0}}$ in (2.1) and applying multichannel interpolation, it is easy to have the interpolation formula for recurrent non-uniformly distributed data:

[TABLE]

The resulting interpolating functions are:

[TABLE]

It is noted that for the case $\alpha=\frac{\pi}{m_{0}}$ , the formula (3.1) reduces to the uniform sampling interpolation. Another fact is that if $m_{0}$ is larger than the half bandwidth of $f$ , then the reconstruction is exact. Most often, one may have no need to compute the interpolating functions, since $y_{i,\alpha}(t)$ in (3.2) and (3.3) can be implemented by FFT efficiently [26].

Importantly, the interpolation consistency holds for the formula (3.1). Namely, the following identities hold:

[TABLE]

where $p,q=1,2,\dots,m_{0}$ and $\delta_{pq}=1$ for $p=q$ and $\delta_{pq}=0$ otherwise. By direct computation from interpolating functions (3.2) and (3.3), formulas (3.4) and (3.5) hold. A concrete example is depicted in Figure 1. Here, $m_{0}=4,~{}\alpha=\frac{\pi}{2m_{0}}$ and the original function is given by

[TABLE]

We see that the red dash-dot line passes through all the red circles.

3.2 Recurrent nonuniform samples and derivatives

In this part we consider a kind of recurrent multichannel interpolation which involves nonuniform samples and derivatives. Let $b_{1}(n)=e^{\bm{i}n\alpha}$ and $b_{2}(n)=\bm{i}n$ , then we have a matrix defined by

[TABLE]

It is easy to get its inverse as

[TABLE]

Set

[TABLE]

It should be noted that if $\alpha=\frac{\pi}{m_{0}}$ , then the denominator in $v_{i,n,\alpha}(t)$ , i.e., $m_{0}+n-ne^{\bm{i}m_{0}\alpha}=m_{0}+2n$ would become [math] for $n=-\frac{m_{0}}{2}$ . Otherwise, we have the interpolating functions

[TABLE]

Moreover, a bandlimited function $f\in B_{\mathbf{N}}$ can be exactly reconstructed by

[TABLE]

provided that $m_{0}\geq\frac{\mu(I^{\mathbf{N}})}{2}$ .

As mentioned in [26], the interpolating functions $y_{k,\alpha}(t-t_{p})$ for $k=1,2$ can be calculated by taking FFT for $v_{k,n,\alpha}(t)$ with respect to $n$ . When $\alpha=0$ , the formula (3.7) reduces to a kind of multichannel interpolation for uniformly distributed data $\{f(t_{p})\},\{f^{\prime}(t_{p})\}$ :

[TABLE]

where

[TABLE]

We illustrate the interpolation formula (3.7) in Figure 2 for recurrent non-uniformly distributed data of $f(t)$ given by (3.6). The red circles represent the samples of $f(t)$ . The reconstructed function (in red dash-dot line) passes through all the red circles. Besides, the blue line and red dash-dot line have the same slope at the particular positions (shown by black asterisks), and the $t$ -coordinates of red circles and black asterisks are interlaced and bunched.

4 Multichannel interpolation of generic non-uniformly distributed data

Although referred to as nonuniform, there are restrictions on location of samples for recurrent nonuniform sampling case. The distribution of samples, to some extent, is still regular. Moreover, as mentioned in [31], recurrent nonuniform samples can be regarded as a combination of several mutual delayed sequences of uniform samples. In this part, we consider a more general interpolation formula which is applicable to generic non-uniformly distributed data. Thanks to the finite summation in (2.2), it is possible to consider a specific case. Let $M=\mu(I^{\mathbf{N}})$ , $I_{1}=\{N_{1}\}$ , $K=1$ in (2.2), then we construct a matrix

[TABLE]

Under this setting, one may drive various nonuniform sampling interpolation formulas provided that $\mathbf{H}$ is invertible. The key points are how to determine whether $\mathbf{H}$ is invertible and how to calculate the inverse. Unlike $\mathbf{H}_{n}$ in the normal case, $\mathbf{H}$ is a large complex-valued matrix with high condition number in general. Therefore, in order to achieve a stable reconstruction, it is not feasible to compute the inverse of $\mathbf{H}$ by numerical methods.

4.1 Generic nonuniform samples

Let $0\leq t_{1}<t_{2}<\dots<t_{M}<2\pi$ and $b_{p}(n)=e^{\bm{i}n{t_{p}}}$ for $1\leq p\leq M$ . We have the following matrix:

[TABLE]

It is easy to show that the determinant of $\mathbf{H}$ is

[TABLE]

That means that $\mathbf{H}$ is invertible. Denote by $z_{p}(k)$ the $(p,k)$ -th element of $\mathbf{H}^{-1}$ . It can be shown that

[TABLE]

Therefore we conclude that a periodic bandlimited function $f\in B_{\mathbf{N}}$ can be exactly reconstructed from its $M\geq\mu(I^{\mathbf{N}})$ non-uniformly distributed samples. The interpolation formula is given by

[TABLE]

where

[TABLE]

We can compute $N(>M)$ function values of $h_{p}(t)$ by taking $N$ fast Fourier transform for $\{z_{p}(k)\}_{k}$ through zero padding. By applying some trigonometric identities, the interpolating function $h_{p}(t)$ can be simplified into a simper form. By a few basic calculations,

[TABLE]

Therefore

[TABLE]

Note that the left hand side of (4.2) is a trigonometric polynomial with respect to $t$ . Expanding the product, we have

[TABLE]

where

[TABLE]

By similar arguments to (4.2), we have

[TABLE]

It follows that

[TABLE]

This formula is consistent with the result presented in [16] by selecting specific values for parameters $N_{1}$ and $M$ . In comparison to the proof of this result in [16], the proposed derivation is simpler and more understandable.

We illustrate the interpolation formula (4.1) in Figure 3 for nonuniform samples of $f(t)$ given by (3.6). The red circles represent the randomly selected nonuniform samples of $f(t)$ . The reconstructed function (in red dash-dot line) passes through all the red circles. For the case $t_{p}=\frac{2\pi(p-1)}{M}$ , $0\leq p\leq M-1$ , the formula (4.1) reduces to the uniform sampling interpolation given in [26].

4.2 Generic nonuniform samples and derivatives

The fact that a bandlimited function could be reconstructed from the values of the function and its derivative is well known [30]. However, the samples involved in such a theorem are uniformly distributed. Let $t_{1},t_{2},\dots,t_{m_{0}}$ be arbitrary $m_{0}$ non-uniformly spaced points on $[0,2\pi)$ . Suppose that $f\in B_{\mathbf{N}}$ with $\mu(I^{\mathbf{N}})\leq M=2m_{0}$ . There is a question of whether $f$ can be perfectly reconstructed from the samples of itself and its first derivative (i.e., $\{f(t_{p}),f^{\prime}(t_{p})\}_{p=1}^{m_{0}}$ ). It is tantamount to asking whether $\mathbf{H}$ is invertible. Here $\mathbf{H}=[v_{kj}]$ is a $M$ -th order square matrix with

[TABLE]

The answer is affirmative. In this subsection, we derive the main result of the current paper: interpolation for non-uniformly distributed samples of a function and its derivative. The interpolating functions are presented in closed-form and the error of reconstructing a non-bandlimited function by proposed formula will be discussed in the next section.

Let $\widetilde{\mathbf{H}}=[\widetilde{v}_{kj}]$ with

[TABLE]

It is easy to see that

[TABLE]

Note that $\det\widetilde{\mathbf{H}}$ is a function of $t_{1},t_{2},\dots,t_{m_{0}}$ . The following lemma gives a recursive relation of $\det\widetilde{\mathbf{H}}$ .

Lemma 4.1

Let $\widetilde{\mathbf{H}}$ be given above. Then its determinant satisfies the following recursive relation

[TABLE]

Proof. Applying some column operations to $\widetilde{\mathbf{H}}$ , it follows that $\det{\widetilde{\mathbf{H}}}$ is equal to the determinant of following matrix

[TABLE]

Subtracting the multiple $e^{\bm{i}t_{1}}$ of row $(k-1)$ from row $k$ for $k=2m_{0},2m_{0}-1,\cdots,2$ successively, we remove first column without changing the determinant:

[TABLE]

where $x_{k}=e^{\bm{i}t_{k}}-e^{\bm{i}t_{1}}$ for $k=2,3,\cdots,m_{0}$ . Subtracting the multiple $\frac{e^{\bm{i}t_{k}}}{x_{k}}$ of column $(2k-2)$ from column $(2k-1)$ and extracting $x_{k}$ from column $(2k-2)$ and $(2k-1)$ for $k=2,3,\cdots,m_{0}$ successively, we reach the result of

[TABLE]

where $\widetilde{\mathbf{H}}^{(1)}$ equals

[TABLE]

For $\widetilde{\mathbf{H}}^{(1)}$ , extracting $e^{\bm{i}t_{1}}$ from the first column and subtracting the multiple $e^{\bm{i}t_{1}}$ of row $(k-1)$ from row $k$ for $k=2m_{0}-1,2m_{0}-2,\cdots,2$ successively, we remove first column of $\widetilde{\mathbf{H}}^{(1)}$ and reach the result of

[TABLE]

where $\widetilde{\mathbf{H}}^{(2)}$ equals

[TABLE]

Subtracting the multiple $\frac{e^{\bm{i}t_{k}}}{x_{k}}$ of column $(2k-3)$ from column $(2k-2)$ and extracting $x_{k}$ from column $(2k-3)$ and $(2k-2)$ for $k=2,3,\cdots,m_{0}$ successively, we get

[TABLE]

Then the recursive relation (4.3) follows from (4.5), (4.6) and (4.7). The proof is complete. $\Box$

Since $\det{\widetilde{\mathbf{H}}}(t_{m_{0}})=\det\begin{bmatrix}1&N_{1}\\ e^{\bm{i}t_{m_{0}}}&(N_{1}+1)e^{\bm{i}t_{m_{0}}}\end{bmatrix}=e^{\bm{i}t_{m_{0}}}$ . By induction, we conclude that

[TABLE]

It follows that

[TABLE]

Therefore $\mathbf{H}$ is invertible. Let $w_{j}(k)$ denote the $(j,k)$ element of $\mathbf{H}^{-1}$ . We define the interpolating functions $\phi_{p}(t)$ and $\psi_{p}(t)$ as follows:

[TABLE]

Then we have a theorem about nonuniform multichannel interpolation as follows.

Theorem 4.2

Let $0\leq t_{1}<t_{2}<\dots<t_{m_{0}}<2\pi$ be non-uniformly distributed points. Suppose that $f\in B_{\mathbf{N}}$ with $\mu(I^{\mathbf{N}})\leq M=2m_{0}$ . Then it can be exactly recovered by the following interpolation formula

[TABLE]

To derive the closed form expressions of $\psi_{p}$ and $\phi_{p}$ , the direct approach is to compute the inverse of $\mathbf{H}$ . This is, as discussed earlier, not a feasible approach. Fortunately, the interpolating functions can be computed tactfully by introducing some auxiliary matrices. Firstly, we need to compute cofactor matrix of $\mathbf{C}$ defined by (4.4). Constructing an auxiliary matrix $\mathbf{A}$ by substituting the second column of $\mathbf{C}$ with $[1,e^{\bm{i}t},\cdots,e^{\bm{i}(2m_{0}-1)t}]$ will bring convenience to the computation:

[TABLE]

On the one hand, by similar arguments to the computation of $\det\widetilde{\mathbf{H}}$ , we get that

[TABLE]

On the other hand, the cofactor expansion of $\det\mathbf{A}$ along the second column gives:

[TABLE]

where $C_{k2}$ is the $(k,2)$ cofactor of $\mathbf{C}$ . By comparing the coefficients of $e^{\bm{i}(k-1)t}$ in (4.2) and (4.12), we obtain the expression of $C_{k2}$ for $k=1,2,\cdots,2m_{0}$ . For example,

[TABLE]

Let $H_{kj}$ denote the $(k,j)$ cofactor of $\mathbf{H}$ . Note that $\mathbf{C}$ can be constructed from $\mathbf{H}$ by applying some column operations. We immediately have the following relations:

[TABLE]

It is well known that the elements of $\mathbf{H}^{-1}$ can be expressed by cofactors of $\mathbf{H}$ , namely

[TABLE]

By similar arguments to (4.2) and (4.12), we have that

[TABLE]

Plugging (4.13) and (4.15) into (4.8) and applying (4.16), it follows that

[TABLE]

More efforts are needed to compute $\psi_{p}(t)$ due to the partial derivative operation in (4.14). For simplicity, we denote Eq.(4.16) and $\frac{\partial}{\partial t_{p}}\prod_{q=1,q\neq p}^{m_{0}}(e^{\bm{i}t_{q}}-e^{\bm{i}t_{p}})$ by $\xi_{p}(t)$ and $\gamma_{p}$ respectively. By some direct computations, we have that

[TABLE]

and

[TABLE]

From (4.13) and (4.14), it follows that

[TABLE]

Plugging (4.18) and (4.15) into (4.9), and applying (4.17) and (4.16), we get that

[TABLE]

Next we shall verify that $\phi_{p}$ and $\psi_{p}$ satisfy the following interpolation consistency:

[TABLE]

for $1\leq p\neq q\leq m_{0}$ . This consistency guarantees that

[TABLE]

even if the reconstruction is not exact. We only give the validation of (4.20) and omit the proof of (4.21) for the sake of brevity. It is obvious that $\phi_{p}(t_{p})=\phi_{p}(t_{q})=0$ . Let

[TABLE]

It is easy to check that $x_{p}^{\prime}(t_{p})=z_{p}(t_{p})=1$ and $x_{p}(t_{p})=z_{p}(t_{q})=z_{p}^{\prime}(t_{q})=0$ . Therefore

[TABLE]

Figure 4 illustrates multichannel interpolation of non-uniformly distributed data and its interpolation consistency. The blue line displays function given by (3.6). The nonuniform grid points are randomly selected as $(t_{1},t_{2},t_{3},t_{4})=(0.2998,1.5866,3.4062,5.0281)$ . The red dash-dot line presents the interpolated result for $(t_{p},f(t_{p}))$ , $p=1,2,3,4$ . We can see that not only the red dash-dot line pass through all the data points but also it is tangent to the blue line at each point.

5 Numerical examples and error analysis

5.1 Numerical examples

According to the types of samples, we abbreviate the interpolation formulas (3.1), (3.7), (4.1) and (4.10) as RN1, RN2, GN1 and GN2 respectively for simplicity. Specially, (4.1) and (4.10) are respectively abbreviated as U1 and U2 if the samples are uniformly spaced. Figure 5 illustrates the inclusion relations of these formulas.

We use the aforementioned formulas to reconstruct non-bandlimited functions. As in [26], we select

[TABLE]

as the test function. Let $f(t)=\Re[\Phi(e^{\bm{i}t})]$ , then its Hilbert transform is $\mathcal{H}f(t)=\Im[\Phi(e^{\bm{i}t})]$ by the theory of Hardy space. In the following, we compare the performance of proposed several formulas for reconstructing $f$ and $\mathcal{H}f$ . The results are listed in Table 1. Denote by $\hat{f}(t)$ the reconstructed result, then the relative mean square error (RMSE) is given by

[TABLE]

Similarly, we denote by $\delta_{2}$ the RMSE for reconstructing the Hilbert transform $\mathcal{H}f$ . In the experiments, $\alpha$ is selected as $\frac{\pi}{N}$ for RN1 and RN2, if the total number of samples is $N$ . For GN1 and GN2, the nonuniform grids are randomly generated by

[TABLE]

where $\zeta_{n}$ and $\eta_{n}$ are i.i.d. sequences of random variables with uniform distribution on $(0,\frac{2\pi}{3N})$ and $(0,\frac{4\pi}{3N})$ respectively. To give a more comprehensive presentation for GN1 and GN2, we repeat each experiment of generic nonuniform sampling for 100 times. Accordingly, $\delta_{1}$ and $\delta_{2}$ of GN1 and GN2 are averaged over these 100 times experiments, and the corresponding variances are also provided.

Some results for reconstructing $f$ and $\mathcal{H}f$ are depicted in Figure 6 and 7. Visually there is no much difference among these reconstructed results by different formulas provided that the same number of samples are used. Roughly, some conclusions could be drawn from the numerical results as follows.

If the same amount of data is employed to reconstruct $f$ (or $\mathcal{H}f$ ), the fluctuations of RMSE caused by the different data types and data distribution patterns are not significant. In other words, the amount of data is the chief factor that affects performance of the reconstruction. 2. 2.

The more grid points the data is distributed on, the better performance of the reconstruction behave. We can see this by comparing the reconstructed results of RN2 and GN2. In addition, the more even the data distribution is, the better performance of the reconstruction behave. We can see this by comparing the reconstructed results of GN1 and U1, or GN2 and U2. 3. 3.

In general, reconstructing a function from its own samples performs slightly better than the reconstruction that involves other types of data. This can be seen from the reconstructed results of GN1 and GN2.

The last two conclusions are certainly based on the premise that the same amount of data is used for reconstruction. And an additional observation is that $\delta_{1}$ and $\delta_{2}$ are nearly equal in each experiment， since the Fourier coefficients of $f$ and $\mathcal{H}f$ have the same absolute value for all $n\in\mathbb{Z}\setminus\{0\}$ .

5.2 Error analysis

In the previous subsection, we presented the reconstruction errors for the proposed interpolation formulas experimentally. In this part, we will give the error estimations analytically which are very important to the reliability of the reconstruction methods.

We denote by $f_{\tau}(t)=f(t-\tau)$ the shifted function of $f$ . Let $\mathcal{T}_{\mathbf{N}}$ be a reconstruction operator corresponding to any one of the aforementioned interpolation formulas. Here $\mathbf{N}$ represents the location of Fourier coefficients for reconstructed function $\mathcal{T}_{\mathbf{N}}f$ . It is easy to see that

[TABLE]

Thus $\mathcal{T}_{\mathbf{N}}f(t-\tau)\neq\mathcal{T}_{\mathbf{N}}f_{\tau}(t)$ . Not just for GN2, most of the other interpolation formulas are not shift-invariant in general. Therefore the MSE defined by

[TABLE]

is not independent on $\tau$ . There is no doubt that $\varsigma(f,\mathbf{N},\tau)$ is ${2\pi}$ periodic in $\tau$ . Note that the time shift $\tau$ could be viewed as the phase difference of $f$ and $f_{\tau}$ . And the exact phase of a function or a signal is generally unknown in most practical applications [23]. Hence, we need to compute the averaged error

[TABLE]

As can be seen from the previous section that the derivation of GN2 is more arduous than the others. In the following, we derive the expression of averaged error for GN2. From (5.3), (4.8) and (4.9), we rewrite $\mathcal{T}_{\mathbf{N}}f_{\tau}(t)$ as

[TABLE]

It is noted that $\mathbf{N}$ is equal to $\{N_{1},N_{1}+2m_{0}-1\}$ in the above formula. Applying the Parseval’s identity, we have that

[TABLE]

Similarly,

[TABLE]

where

[TABLE]

To simplify (5.4), (5.5) and (5.6), we need to introduce some identities:

[TABLE]

Integrating the both sides of (5.4) and (5.6) on $\mathbb{T}$ with respect to $\tau$ and making use of the above identities, we get that

[TABLE]

and

[TABLE]

From the definition of $w_{j}(k)$ , for any $n\in I^{\mathbf{N}}$

[TABLE]

It follows that

[TABLE]

and

[TABLE]

Therefore the square of the averaged error for GN2 is given by

[TABLE]

where

[TABLE]

Similarly, we can get the averaged errors for GN1 and RN2 respectively as

[TABLE]

where

[TABLE]

with $k_{n}=\text{fix}\left(\frac{n-N_{1}}{m_{0}}\right)+1$ and $\text{fix}(x)$ rounds $x$ to the nearest integer toward zero. Note that the other interpolation formulas can be subsumed in the above three cases, therefore we obtain all the averaged errors of six aforementioned formulas. The sequences $Er(\text{U1},\mathbf{N},n)$ , $Er(\text{GN1},\mathbf{N},n)$ , …, $Er(\text{RN2},\mathbf{N},n)$ are depicted graphically in Figure 8. Here $\mathbf{N}=(N_{1},N_{2})=(-31,32)$ , thereby $m_{0}=32,M=64$ . For GN1 and GN2, the nonuniform grids are randomly generated by (5.1) and (5.2) respectively. The domain for each $Er(\mathbf{N},n)$ plotted in Figure 8 is set as $\{N_{2}+1\leq n\leq N_{2}+3\mu(I^{\mathbf{N}})\}$ . The theoretical analysis of error is in accord with the result of numerical examples and therefore the conclusions made in the previous subsection are underpinned.

The proposed interpolation method involves non-uniformly spaced multichannel samples. There are notable existing interpolation methods involving nonuniform or multichannel samples. We provide the Table 2 to compare these existing results. Among the numerous sampling or interpolation methods , we only present several typical types in Table 2. It is noted that the representative references listed here are far from complete.

6 Application to image recovery

In the previous sections, we dealt with techniques for reconstructing a continuous function from different types of discrete samples. In this section, we introduce a simple application of the proposed interpolation formulas to image recovery. To begin, consider Figure 9 (b), which is severely degraded because of the damaged pixels. Suppose that the damaged pixels are non-uniformly located. The goal of this part is to recover the missing pixels via interpolation.

Note that the proposed formulas are one-dimensional, we have to compute interpolation result for each row of image first, and then apply interpolation for each column by using the same operations. As the distantly separated image regions are irrelevant virtually, we should treat the reconstruction problem locally. In the following, the test image is set to be Lena ( $256\times 256$ ), and it is degraded by wiping out $43.5\%$ randomly selected pixels, see Figure 9 (b). Each row of image is divided into $32$ equal parts, namely $8$ pixels per part. Repeating interpolation process through the image pieces produced by dividing, we obtain values for all the missing pixels. Applying the same operations to each column, we have another reconstructed result. It is noted that the dividing treatment has an additional benefit that it makes computation complexity linear in the size of image.

It is natural to average two reconstructed results. Besides, we need to convert interpolation result into unsigned 8-bit integer type. A direct way for such a conversion is based on

[TABLE]

where $\mathcal{I}_{xy}$ is the intensity value at location $(x,y)$ . For a more elaborate conversion, we introduce a correction for the values produced by interpolation. Let $\Lambda_{xy}$ be the $3\times 3$ neighborhood centered on $(x,y)$ , the correction is defined as

[TABLE]

where $\mathcal{I}_{xy}$ and $\hat{\mathcal{I}}_{xy}$ are the intensity values at location $(x,y)$ before and after correction respectively. From the definition, this correction is certain to be convergent after finite iterations. In practice, more fortunately, it can be convergent generally by $3$ or $4$ iterations.

Note that the damaged pixels can be also viewed as impulse noise (also called salt-and-pepper noise) in an image. It is known that the median filter, which is a very useful order-statistic filter in image processing, is particularly effective in the reduction of impulse noise [32, 33]. Basically, to perform median filtering at $(x,y)$ is to determine the median for values of the pixel in $\Lambda_{xy}$ and assign that median to $(x,y)$ in the filtered image.

We are in position to compare the performance of interpolation method and median filtering in the problem of restoring damaged pixels. Specifically, we consider three methods: GN1, GN2 and median filter (MED for short). In general, GN2 requires a prerequisite condition of differentiable since it involves derivative. It would be stretching a point to describe a digital image as a set of samples of a smooth (differentiable) function. Nevertheless, the introduction of difference (also called derivative in some literature without ambiguity) for the original digital image could help to preserve more useful information in the reconstructed image. The experimental results are shown in Figure 9 and Table 3.

Here CRT represents the correction operation. We use relative mean square error (RMSE) $\delta$ , peak signal to noise ratio (PSNR) $\rho$ and correlation coefficient (CC) $\gamma$ , to measure the quality of reconstructed images. They are defined respectively as:

[TABLE]

where $\mathcal{I}$ , $\mathcal{I}_{r}$ denote original and reconstructed image respectively, $\mathcal{I}^{0}$ , $\mathcal{I}_{r}^{0}$ denote their averaged pixel values, and $\left\lVert\cdot\right\rVert_{F}$ denotes Frobenius norm, and $L_{1}$ and $L_{2}$ are the number of rows and columns of $\mathcal{I}$ .

From Table 3, we conclude that the interpolation-based method GN1 performs significantly better than the median filtering method. If there is some information about gradient of original image available to be utilized, the performance of image recovery can be improved further by GN2. These conclusions are also reflected in Figure 9 visually.

It is noted that we consider the image recovery problem only from the point where a digital image is degraded by simply wiping out some pixel values. Besides, the material about recovery methods developed in this section is far from exhaustive. Even so, the nonuniform-interpolation-based image recovery methods perform well and are easily implemented. It is conceivable that these methods could be integrated into some more comprehensive image recovery approaches. These further explorations, although of importance in image processing, are beyond the scope of this paper.

7 Conclusion

Several interpolation formulas associated with non-uniformly distributed data are presented. If the signal to be reconstructed is bandlimited, then it is possible to reconstruct the entire signal by sampling it with the total number of samples larger than the corresponding bandwidth. For the case of non-bandlimited signal, quantitative error analysis for reconstructing is also analyzed. It has been shown that the introducing derivative samples of function can improve reconstruction result significantly. As an application, several nonuniform-interpolation-based algorithms for recovering a certain kind of corrupted images are demonstrated. The performance is satisfactory.

Bibliography33

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] A. I. Zayed, Advances in Shannon’s sampling theory . CRC press, 1993.
2[2] E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Physics in Medicine & Biology , vol. 53, no. 17, p. 4777, 2008.
3[3] Z.-P. Liang and P. C. Lauterbur, Principles of magnetic resonance imaging: a signal processing perspective . SPIE Optical Engineering Press, 2000.
4[4] B. F. Burke and F. Graham-Smith, An introduction to radio astronomy . Cambridge University Press, 2009.
5[5] J. Yen, “On nonuniform sampling of bandwidth-limited signals,” IRE Transactions on circuit theory , vol. 3, no. 4, pp. 251–257, 1956.
6[6] K. Yao and J. Thomas, “On some stability and interpolatory properties of nonuniform sampling expansions,” IEEE Transactions on Circuit Theory , vol. 14, no. 4, pp. 404–408, 1967.
7[7] A. J. Jerri, “The Shannon sampling theorem - its various extensions and applications: a tutorial review,” Proceedings of the IEEE , vol. 65, no. 11, pp. 1565–1596, 1977.
8[8] H. G. Feichtinger and K. Gröchenig, “Irregular sampling theorems and series expansions of band-limited functions,” Journal of Mathematical Analysis and Applications , vol. 167, no. 2, pp. 530–556, 1992.