Periodic Spline-Based Frames for Image Restoration

We present a design scheme to generate tight and semi-tight frames in the space of discrete-time periodic signals, which are originated from four-channel perfect reconstruction periodic filter banks. The filter banks are derived from interpolating and quasi-interpolating polynomial and discrete splines. Each filter bank comprises one linear phase low-pass filter (in most cases interpolating) and one high-pass filter, whose magnitude’s response mirrors that of a low-pass filter. In addition, these filter banks comprise two band-pass filters. We introduce the notion of local discrete vanishing moments (LDVM). In the tight frame case, analysis framelets coincide with their synthesis counterparts. However, in the semitight frames, we have the option to swap LDVM between synthesis and analysis framelets. The design scheme is generic and it enables us to design framelets with any number of LDVM. The computational complexity of the framelet transforms, which consists of calculating the forward and the inverse fast Fourier transforms, practically does not depend on the number of LDVM and does depend on the size of the impulse response filters. The designed frames are used for image restoration tasks, which were degraded by blurring, random noise and missing pixels. The images were restored by the application of the Split Bregman Iterations method. The frames performances are evaluated. A potential application of this methodology is the design of a snapshot hyperspectral imager that is based on a regular digital camera. All these imaging applications are described.

1. Introduction.Restoration of multidimensional signals that were corrupted and/or damaged and/or noised is a major challenge that the signal/image processing community faces nowadays when rich multimedia content is the most popular data being transmitted over diverse networks types including mobile.Quality degradation in multidimensional signals can come from sampling, acquisition, transmission through noisy channels, to name some.Restoration of multidimensional signals includes denoising, deblurring, recovering missing or damaged samples or fragments (inpainting in images), resolution enhancement and super resolution.The processing goals are to improve the visual perception of still and video signals.
Hyperspectral imaging is a relatively new field of investigation.Its imager has a large number of wavebands in a given wavelength range.For example, the range can be 10, hundreds or even thousands of wavebands.Spectral imaging has a multitude of applications in many fields including scientific research, engineering and medical equipment ranging from a color camera to sophisticated astronomical equipment.Currently, known spectral imagers usually use dynamic spectral filters, arrays, scanning procedures and multiple-lens optical schemes that are quite complicated and expensive.However, a combination of optics with proper digital processing described in this paper is able to provide "snapshot" spectral imagers based on regular digital cameras [16].
Recent developments in wavelet frames (framelets) analysis provide innovative and powerful tools to face faithfully and robustly the above challenges.Framelets produce redundant expansions for multidimensional signals that, in particular, provide an additional sparsity to the signals representation due to better adaptation abilities and due to redundant representations.A valuable advantage, which redundant representations hold, is their abilities to restore missing and incomplete information and to represent efficiently and compactly the data.
Frame expansions of signals demonstrate resilience to the coefficients disturbances and losses [18,17,21,1].Thus, frames can serve as a tool for error correction for signals transmitted through lossy channels.
Implicitly, this resilience is utilized in signal/image restoration, which is based on a prior assumption that a frame expansion of a given signal/image is sparse.In principle, only part of the samples/pixels is needed to (near) perfect object restoration.This approach, which is a variation of Compressive Sensing methodology ( [25]), proved to be extremely efficient for image restoration.Practically, compressed sensing approach is implemented via minimization of a parameterized functional where the sparse representation is reflected in the l 1 norm of the transform coefficients.• 1 minimization does not provide an explicit solution and can be resolved only by iterative methods.The split Bregman iteration (SBI) scheme, which was recently introduced in [15], provides for that a fast and stable algorithm.Variations of this scheme and its application to image restoration using wavelet frames are described in [29,23], to mention a few.A variety of impressive results on image restoration were reported in the last couple of years.A survey is given in [28] while a recent development is described in [23].
The SBI methodology proved to be useful in the design of "snapshot" spectral imagers based on regular digital cameras.It enables to reconstruct a spectral cube from a dispersed-diffused sensors' array.
Due to applications diversity, it is important to have a library of wavelet frames in order to select a frame that fits best a specific task.Forward and inverse transforms in iterative algorithms are repeated many times, therefore, members in this library must have fast and stable transforms implementation.Waveforms symmetry with the availability of vanishing moments are also important in order to avoid distortions when thresholding is used.To satisfy these requirements, most of the framelet systems in practical use operate with the compactly supported framelets and the transforms are implemented by finite (and short) impulse response (FIR) oversampled filter banks [12,9,10,13].A common requirement is to have tight frames.Increase in the number of vanishing moments in framelets requires an increase in the filters' length, which results in higher computational cost.It also can produce boundary artifacts in the processed images.
The oversampled perfect reconstruction (PR) filter banks generate wavelet-type frames in the space of discrete-time signals [11,8].Utilizing infinite impulse response (IIR) filter banks with a relaxation of the tightness requirement provides a number of additional opportunities.Properties such as symmetry, interpolation, flat spectra, combined with fine time-domain localization of framelets can be easily achieved as well as a high number of vanishing moments [6,7].In these papers, the key point is the low-pass filters design, which, being applied to the even subarray of a signal, well approximate the odd subarray (the prediction filters).A natural source for such filters are the discrete and polynomial splines.A number of 3-framelet systems was derived from the discrete splines in [6,7].The transforms are implemented in a fast way using recursive filtering.Non-compactness of the waveforms supports is compensated by their exponential decay as time goes to infinity.In principle, any number of vanishing moments can be derived but the implementation computational cost grows fast.
This drawback can be overcome by switching to a periodic setting, which is the subject of this paper.A variety of four-channel PR filter banks, where the lowpass filters are derived from interpolating and quasi-interpolating polynomial and discrete splines, are designed.These filter banks generate a library of 4-framelet tight and semi-tight frames with diverse properties.
The transforms implementation is reduced to application of the direct and the inverse fast Fourier transforms (FFT).Actually, the implementation cost does not depend on the number of vanishing moments in the framelets.Therefore, framelets utilization with any number of vanishing moments and filter banks with IIR filters become practical.
A number of framelets systems is explicitly designed and described in this paper.However, the design scheme is generic and this framelets library can be easily expanded.Preliminary results in the periodic discrete-time frame design with three framelets are reported in [32].Recently, framelets in the space of periodic continuous functions are studied in [19].
The designed framelets libraries were tested for image restoration and demonstrated a high quality performance.The framelet libraries enable us to select a frame that best fits each specific application.For example, in the restoration experiments with relatively smooth images such as "Lena" and "Window", the best results were achieved with compactly supported framelets derived from pseudo-splines ( [13]) and from the quasi-interpolating quadratic splines.However, in the experiments with images comprising a fine texture, such as "Barbara" and "Fingerprint", the framelets with high number of vanishing moments, which are derived from the discrete interpolating splines, significantly outperformed all others.Note that in many experiments, the semi-tight frames outperformed tight frames.For example, in the experiments with "Lena" and "Window" images, the best performance was demonstrated by the semi-tight frame derived from pseudo-splines, while in numerical simulation of the snapshot hyperspectral imaging, the semi-tight frame derived from the quasi-interpolating quadratic splines produces the best quality output.
The paper is organized as follows: The introductory Section 2 recalls the notion of periodic discrete-time signals and of periodic filter (p-filter).The PR periodic filter banks (p-filter banks) are presented and characterized via the polyphase representation of their discrete Fourier transform (DFT).Frames in the discrete-time periodic signals space, which originated from PR oversampled p-filter banks, are also introduced in Section 2. Design of four-channel p-filter banks, which generate tight and semi-tight frames, is presented in Section 3. Section 4 discusses the restoration of sampled polynomials by low-pass p-filters and introduces the notion of local discrete vanishing moments (LDVM).A collection of FIR and IIR p-filters from the interpolating and quasi-interpolating polynomial and discrete splines is derived in Section 5. A number of four-channel p-filter banks with spline p-filters that generate tight and semi-tight frames, is described in section 6.The designed tight and semi-tight frames are used in Section 7 for restoration of images, which were degraded by blurring, random noise and missing pixels.The frames performances are evaluated.The section is concluded by an extensive discussion on the restoration experiments.Section 8 describes experiments on numerical simulation of the snapshot hyperspectral imaging.
2. Preliminaries: Periodic discrete-time signals and filters.In this section, we present notions and some known facts.
The We use the notation for the DFT of signals belonging to the space Relations between signals and their DFT are: Circular convolution: Finite differences: The forward finite differences of a signal x is Inverse Problems and Imaging Volume 9, No. 3 (2015), 661-707 The even-order central finite differences are Polyphase representation of the DFT The signals are called the even and the odd polyphase components of the signal x ∈ Π[N ], respectively.

Periodic filters. A linear operator
We use the notation h for both a p-filter and its IR {h[k]}.The DFT of the IR of the p-filter h is ĥ and it is called the frequency response (FR).
Periodic filtering of a signal reduces to multiplication in the frequency domain: The FR of a p-filter can be represented in a polar form ĥ[n] = ĥ[n] e i arg( ĥ[n]) , where the positive N -periodic sequence ĥ[n] is called the magnitude response (MR) of h and the real-valued 2π-periodic sequence arg( ĥ[n]) is called the phase response of h.A p-filter is referred to as a linear phase if its phase response is linear in n.If the IR of a filter h is symmetric or antisymmetric within the interval k = −N/2, ..., N/2 − 1, then h is a linear phase filter.

Multirate
is called upsampling the signal x by factor M .In the rest of the paper, we apply down(up) sampling by factor 2. Then, the downsampled signal (↓ 2)x = x 0 is the even polyphase component of the signal x, while the upsampled signal (↑ , where Application of the p-filter h to a signal y where Interpolating p-filters: If the DFT of the even polyphase component of a p-filter h is constant, i.e. ĥ0 [n] 1 ≡ C, then the p-filter is called interpolating.In this case, Eq. (2.4) implies that the DFT of the zero polyphase component of the output  (2.5), are used as an input to the synthesis p-filter bank and the output signal is x = x, then the pair of analysis-synthesis p-filter banks form a perfect reconstruction (PR) p-filter bank.
If the number of channels S equals to the downsampling factor M then the p-filter bank is said to be critically sampled.If S > M then the p-filter bank is oversampled.Critically sampled PR p-filter banks are used in wavelet analysis, while oversampled PR p-filter banks serve as a source for discrete-time wavelet frames design.
In the rest of the paper, we deal with p-filter banks, whose downsampling factor is M = 2 and h0 and h 0 are the low-pass filters.

Characterization of p-filter banks. Assume that
H def = hs , s = 0, ..., S − 1, is an analysis p-filter bank with the downsampling factor 2. Then, its application to a signal x ∈ Π[N ] produces S signals from Π[N/2] in the following way: Assume that produces a signal from Π[N ] in the following way: Equations (2.6) and (2.7) can be rewritten in a matrix form by where the S ×2 analysis and the 2×S synthesis polyphase matrices are, respectively, If the relations (2.8) hold for all n ∈ Z when I 2 is the 2 × 2 identity matrix, then Thus, Eq. (2.8) is the condition for the pair H, H of p-filter banks to form a PR p-filter bank.
If the frame bounds A and B are equal to each other then the frame is said to be tight.
If the system Φ is a frame then there exists another frame Φ def = {φ l } L−1 l=0 in the space Π[N ] such that any signal x ∈ Π[N ] can be expanded into the sum x = L−1 l=0 x, φl φ l .The analysis Φ and the synthesis Φ frames can be interchanged.Together they form the so-called bi-frame.If the frame is tight then Φ can be chosen as Φ = c Φ.
If the elements { φl } of the analysis frame Φ are not linearly independent (L > N ) then many synthesis frames can be associated with a given analysis frame.In this case, the expansions x = L−1 l=0 x, φl φ l provide redundant representations of the signal x.
It was established in [11] that the PR filter banks operating in the space l 1 of decaying discrete-time signals generate frames for this space.A similar fact was proved in [32] for the p-filter banks operating in Π[N ]. (2.9) ψs where hs [l] and {h s [l]} are the impulse responses of the p-filters hs and h s , respectively.Then, Eqs.(2.6) and (2.7) imply that where (2.10) The expansion coefficients are the inner products x, ψs The notation • T means matrix transposition.If the condition in Eq. (2.12) is satisfied, then the synthesis filter bank can be chosen to be equal to the analysis filter bank (up to a constant factor).
If S > 2 then the representation in Eq. (2.11) of a signal from Π[N ] is redundant.The redundancy ratio for one-level frame transform is ρ = S/2.
The signals ψs [1] and ψ s [1] , s = 0, ..., S − 1, are called the analysis and synthesis discrete-time framelets of the first decomposition level, respectively.To increase the redundancy of a signal representation, the frame transform is applied to the low frequency signal y 0 where The signals ψs [2] and ψ s are called the analysis and the synthesis discrete time periodic framelets of the second decomposition level, respectively.
The iterated transform with the polyphase matrices P[µ] [n] and , where (2.13) The synthesis and analysis framelets are derived iteratively 3. Design of four-channel p-filter banks for frames generation.In this section, we discuss the design of four-channel p-filter banks that generate tight and semi-tight frames in the space Π[N ] of periodic signals.
3.1.Four-channel perfect reconstruction p-filter banks.The PR condition for a pair of the analysis H = h0 , h1 , h2 , h3 and the synthesis H = h 0 , h 1 , h 2 , h 3 p-filter banks is expressed via their polyphase matrices (3.1) for all n ∈ Z.The matrix product in Eq. (3.1) can be split into two products: According to Theorem 2.
Design of a four-channel filter bank begins from a linear phase low-pass filter 1 is assumed to be a rational function of ω n = e 2πin/N with real coefficients that has no poles for n ∈ Z.In addition, assume that ĥ0 [n] is symmetric about the swap n → −n that implies that ĥ0 T and the product The simplest way to satisfy this condition is to define Due to symmetry of ĥ0 [n], the FR becomes ˆh1 It follows from the assumption in Eq. (3.3) that the product . Thus, design of the PR p-filter bank is reduced to factorization of the matrix Q[n].

Diagonal factorization of the matrix Q[n]
. There are many ways to factorize the matrix Q[n].One way is to define the matrices P 23 [n] and P23 [n] to be diagonal: , which means that the odd polyphase components of the p-filters h2 and h 2 as well as the even polyphase components of the p-filters h3 and h 3 vanish.Consequently, we have to find four sequences ĥ2 Tight frame p-filter banks: If the following inequality holds then, due to symmetry of the rational functions ĥ0 0 [n] 1 and ĥ0 0 [n] 1 , we have 1 − a[n] 1 = P (cos 2πn/N )/R(cos 2πn/N ), where R is strictly positive and P is nonnegative polynomials.Due to Riesz's Lemma [24], the polynomials can be factorized as , where p and q are polynomials with real coefficients and q(ω n ) does not have roots for n ∈ Z.
Thus, we can define ĥ2 Thus, we obtain the PR p-filter bank, whose FRs are where P 23 [n] where Then, the product A consequence of this choice is the fact that the p-filters h 2 and h 3 have a linear phase: ĥ2 The following symmetry properties hold.
Proposition 3.2.The impulse response (IR) of the p-filter h 2 is symmetric about 1/2, while the IR of h 3 is antisymmetric about 1/2.
Inverse Problems and Imaging Volume 9, No. 3 (2015), 661-707 Proof.The IR of the p-filter The claim about the p-filter h 3 is proved similarly.
Since P23 [n] = P 23 [n] T , the four-channel p-filter bank generates a tight frame in the space Π[N ].
The above symmetric factorization scheme is, in essence, similar to the scheme presented in [13].Certainly, other schemes are possible.
3.1.3.Interpolating p-filter banks.Assume that the low-pass p-filter h 0 is interpolating and its frequency response is ĥ0 where the sequence f [n] 1 is a rational function of ω n = e 2πin/N that has no poles for n ∈ Z and Then, the sequence t[n] 1 is factorized to become When the diagonal factorization of the matrix Q[n] is applied, the polyphase submatrices are In the case when the inequality Eq. (3.7) holds and the sequence t (3.11), while the matrices P 23 [n] = P23 [n] T are given in Eq. (3.9).
The key point in the design of a four-channel PR p-filter bank that generates a (semi-)tight frame in the signal space Π[N ] is in the definition of a relevant lowpass p-filter h 0 = h0 .Once we have a low-pass filter, the rest of the p-filters is designed in a way described above.Prior to the presentation of a family of the p-filters originating from polynomial and discrete splines, we briefly discuss how the concepts of the polynomials restoration and vanishing moments can be adapted to the periodic discrete-time setting.

Restoration of sampled polynomials and discrete vanishing moments.
Restoration of polynomials by low-pass filtering coupled with their elimination by respective high-and band-pass filters, which constitute filter banks (vanishing moment property), provide a sparse representation of signals/images in wavelet and frame bases.This is important, for example, for data compression and signals/images restoration.

Inverse Problems and Imaging
Volume 9, No. 3 (2015), 661-707 Restoration of sampled polynomials.Certainly, sampled polynomials do not belong to Π[N ] and p-filters cannot be applied to them.However, the following fact holds, which, in a sense, is a discrete periodic counterpart of the classical Fix-Strang condition [31].
Proposition 4.1.Assume that the FR of a low-pass filter h can be represented as ĥ where m is some natural number and α[n] is a rational function of ω n , which has no poles for n ∈ Z and α[0] = 0. Assume p is a signal from Π[N ], which coincides with a sampled polynomial Proof.Equation (2.1) implies that the sequence (ω n − 1) m is the DFT of the finite difference and ( . Under conditions of the proposition, the sequence {α[n]} can be regarded as the FR of a low-(all-)pass p-filter a.Then, application of the p-filter h to a signal x ∈ Π[N ] can be represented as The finite difference of the signal p is where m is a natural number and α[n] is a rational function of ω n , which has no poles for n ∈ Z and α[0] = 0. Assume p is a signal from Π[N ], which coincides with a sampled polynomial P m−1 of degree m − 1 at the interval p Proof. is similar to the proof of Proposition 4.1.
where m is a natural number and ᾱ[n] is a rational function of ω n that has no poles for n ∈ Z and ᾱ[0] = 0.

Spline-based low-and high-pass p-filters.
The idea behind the design.Assume that the FR of a low-pass interpolating p-filter Denote by In order for the high-pass filter to eliminate smooth signals (for example, fragments of polynomials), the p-filter f should be "predictive" in a sense that its application to the even polyphase component of the signal should "predict" the odd polyphase component and vice versa.
Polynomial and discrete splines are natural sources to derive prediction p-filters from.The idea is to construct either a polynomial or a discrete spline, which (quasi-)interpolates the even samples of a signal, and to predict the odd samples of the signal by spline's values at the midpoints between the (quasi-)interpolation points.Such an approach was explored for the design of the biorthogonal wavelet transforms in [2,5,4].The prediction filters derived from the discrete splines were used in [6,7] for the design of non-periodic interpolating wavelet frames with threechannel filter banks.5.1.Prediction p-filters derived from polynomial splines.A polynomial spline of order p on the grid {t k } is a function that has p − 2 continuous derivatives.The spline of order p coincides at the intervals (t k , t k+1 ) between the grid points with polynomials P p−1 k of degree p − 1.
5.1.1.Periodic interpolating splines.Denote by S p K the space of splines of order p ∈ N defined on the uniform grid {k}, which are periodic with the period K = N/2 = 2 j−1 .A basis in this space is constituted from shifts of B-splines.
The first-order non-periodic B-spline β 1 (t) on the grid {k} is the indicator function on the interval (−1/2, 1/2).The B-spline of order p is derived by the iterated It is supported on the interval (−p/2, p/2), is strictly positive inside this interval and is symmetric about zero where it has its single maximum.There exists an explicit representation of the B-spline of order p: is the periodization of the compactly supported function β p (t).
Any spline S p (t) from S p K is represented by Inverse Problems and Imaging Volume 9, No. 3 (2015), 661-707 If the spline S p (t) interpolates the even polyphase component x 0 = {x[2l]} of a signal x ∈ Π[N ] then its coefficients q[k] can be explicitly calculated via the application of DFT: Approximation properties of interpolating splines are well investigated.In particular, the non-periodic interpolating spline of order p, which consists of piecewise polynomials of degree p − 1, restores these polynomials.It means that splines of order p, which interpolate a polynomial of degree p − 1, coincides with this polynomial.For periodic splines, this property holds locally.This observation justifies the choice of interpolating splines as a source for the design of prediction filters.
To be specific, if the spline S p (t) ∈ S p K , given in Eq. (5.2), interpolates the even polyphase component is the DFT of the B-spline {B p (k + 1/2)} sampled at midpoints between grid points.
The continuous counterparts of the sequences u p [n] 1 and v p [n] 1 , which are the discrete-time Fourier transforms of the non-periodic B-splines, were introduced and studied in [27].Some additional properties of the sequences u p [n] 1 and v p [n] 1 are established in [2,34].Denote Obviously, the sequence f p c [n] 1 satisfies the conditions in Eq. (3.10).Proposition 5.1 ( [2]).If the spline order is either p = 2r or p = 2r − 1 then where r is a natural number and χ p [n] and γ p [n] are rational functions of ω n , which have no poles for n ∈ Z and no root at n = 0.This proposition coupled with Propositions 4.1, 4.3 and Remark 4.1 imply the following Corollary.
Corollary 5.2.The low-pass h 0 and the high-pass h 1 p-filters, whose FR are given in Eqs.(5.5) and (5.6), respectively, locally restore and eliminate sampled polynomials of degree 2r − 1, respectively.Therefore, the p-filter f p c , whose FR f p c [n] 1 is defined by Eq. (5.4), is a proper candidate to be utilized as a prediction p-filter.

Inverse Problems and Imaging
Volume 9, No. 3 (2015), 661-707 Remark 5.1.It is emphasized that the p-filters derived from the odd order 2r − 1 splines restore (eliminate) sampled polynomials of the same 2r − 1 degree as the the p-filters derived from the even order 2r splines.This is a consequence of the socalled super-convergence property of the odd order interpolating splines [34], which claims that the approximation order of such splines at the midpoints between the interpolation points is higher than those at the remaining points of the intervals between interpolation points.Linear spline p = 2:: Quadratic interpolating spline p = 3:: Cubic interpolating spline p = 4:: Comment: We observe that the high-pass p-filters h 1 derived from either quadratic or cubic splines locally eliminate sampled cubic polynomials.However, the structure of the quadratic FR is much simpler than the that of cubic FR.Therefore, in a number of applications, the p-filters derived from quadratic spline are advantageous over the cubic spline p-filters.Interpolating spline of fourth degree p = 5:: It is readily verified that the low-and high-pass p-filters locally restore and eliminate sampled polynomials of fifth degree, respectively.The p-filters, originated from higher order splines, are designed by the application of the DFT to the sampled B-splines.5.2.Prediction p-filters derived from discrete splines.We can see from the above examples that the structure of the spline-based p-filters becomes more complicated as far as the generating spline's order increases, which is necessary for restoration(elimination) of sampled polynomials of higher degrees.However, there is a way to explicitly design a family of prediction p-filters, which provide low-and high-pass p-filters, which restore (eliminate) sampled polynomials of higher degrees.Thus, the corresponding framelets can have any number of LDVMs.
The design scheme is, in essence, similar to the scheme given in Section 5.1.1.The difference is that it is based on the discrete rather than on the polynomial interpolating splines.The discrete splines are the signals from Π[N ], whose properties mimic the propertes of the polynomial splines.The discrete splines and the p-filters' design are described in full details in [4] while in the current paper we summarize the result.
Define the prediction p-filters f 2r d via the frequency responses: (5.11) sin 2r πn/N cos 2r πn/N + sin 2r πn/N .(5.12) If the p-filter f 2r d , whose frequency response is given in Eq. (5.11), is used as a prediction filter in the lifting scheme then the PR p-filter bank is The following proposition is the consequence of Eq. (5.12).
Proposition 5.3.The low-pass p-filters h 0 derived from the prediction p-filter f 2r d by Eqs.(5.13) locally restore sampled polynomials of degree 2r − 1.The highpass p-filters h 1 locally eliminate such polynomials.Consequently, the framelets We emphasize that Eq. ( 5.13) provides an explicit expression for p-filters with an arbitrary approximation accuracy.
Remark 5.3.One-pass non-periodic Butterworth filters were used in [20] for the design of orthogonal non-symmetric wavelets.The computations in [20] were conducted in time domain using recursive filtering.Biorthogonal periodic wavelets derived from discrete splines are presented in [4], while non-periodic biorthogonal wavelets are introduced in [5].

Examples of non-interpolating p-FIR p-filters.
The IR of all the p-filters introduced in Sections 5.1 and 5.2, except for the linear spline p-filters, are infinite, i.e. they occupy the whole Z.In this section, we introduce a few p-filters, whose IR are finite up to periodization (p-FIR p-filters).One such p-filters is originated from a linear interpolating spline.Another way to design the p-FIR p-filters is to use the splines, which are quasi-interpolating quasi-interpolate rather than interpolating the even polyphase component of a signal x and to predict the odd polyphase component by the spline's values at midpoints between grid nodes.The quasi-interpolating splines are studied in [3,33].We describe a design of the p-FIR p-filters originated from quadratic splines.Design of p-filters using quasi-interpolating splines of higher order is similar but their IR is much longer.5.3.1.Quadratic quasi-interpolating spline.The quadratic spline, which interpolates the even polyphase component x 0 of a signal x, is represented as Due to Eqs. (5.3) and (5.8), the DFT of the coefficients is Two initial terms from the series in Eq. (5.14) are taken and denoted as (5.15) The spline S 3 (t) = , where the DFT of the coefficients is q We define two types of low-and high-pass p-filters.
1. Non-interpolating low-pass p-filter:: Define the p-filters by their FR: (5.17) The IR of the p-filters h 0 and h 1 comprise nine terms (up to periodization).The low-pass h 0 and the high-pass h 1 p-filters locally restore and eliminate sampled cubic polynomials, respectively.2. Interpolating low-pass p-filter:: The same approximation order can be achieved by p-filters that have shorter IR.For this, the low-pass p-filter h 0 is defined to be interpolating: (5.18) Once the factorization is accomplished with T [n] 1 = T [n] 1 , the FRs of the band-pass p-filters hs and h s , s = 2, 3, are defined as The PR pair H, H of the p-filter banks, where H = h 0 , h 1 , h2 , h3 and we use the "Symmetric" design of the p-filter banks H that generates tight frames: where Recall that impulse responses of the p-filters hs and h s , s = 0, 1, 2, 3, are the framelets ψs [1] and ψ s [1] of the first level, respectively.If the low-pass p-filter is interpolating and the prediction p-filter f is available then the design of the p-filter bank is reduced to factorization of the sequence The impulse and the magnitude responses of PR p-filter banks H = {h s } , s = 0, 1, 2, 3, which generate tight frames, are displayed in Fig. 6.8.
The impulse and the magnitude responses of p-filters hs and h s , s = 2, 3, that together with the p-filters h s , s = 0, 1, generate semi-tight frames, are displayed in Figs.6.1-6.7.All the Figs.6.1-6.7,except Fig. 6.5, are structured identically.Left pictures display the IR of the p-filters in the following order: Left to right h 2 −→ h2 −→ h 3 −→ h3 .Right pictures display the MR of h 2 and h3 (which coincide with each other) (dashed lines), and the MR of h2 and h 3 (which coincide with each other) (solid line).The top pictures pair illustrates the p-filters arising from symmetric factorization, while the bottom pair does the same when the anti symmetric factorization is used.Figure 6.5 comprises only one pair of pictures, which correspond to the symmetric factorization.

Four-channel p-filter banks with p-FIR p-filters.
A few examples of p-filter banks and p-filter banks with FIR p-filters are given in this section.Linear interpolating spline: The prediction p-filter is f The p-filter h 0 locally restores first-degree sampled polynomials, while the p-filters h 1 and h 2 locally eliminate them.The p-filter h 3 locally eliminates only constants.The framelets ψ 1 [1] and ψ 2 [1] , which are the IRs of the p-filters h 1 and h 2 , respectively, have two LDVMs.Either of them is symmetric.The framelet ψ 3 [1] is antisymmetric and has one LDVM.
The IR of the p-filters h s , s = 0, 1, 2, 3, and their magnitude responses are displayed in Fig. 6.8 (top line).

Inverse Problems and Imaging
Volume 9, No. 3 (2015), 661-707 Quadratic quasi-interpolating spline (interpolating low-pass p-filter): The FR of the prediction p-filter f low-and high-pass p-filters h 0 and h 1 are given in Eqs.(5.16) and (5.18), respectively.The p-filter h 0 locally restores sampled polynomials of third degree while the p-filter h 1 locally eliminates them.Thus, the frameletψ 1 [1] has four LDVMs.The sequence W [n] 1 to be factorized is The following factorization modes of FRs of the p-filters h 2 and h 3 are defined in Eq. (6.2) with A The p-filter bank H = {h s } , s = 0, 1, 2, 3, generates a tight frame.The p-filters h 3 and h 2 locally eliminate sampled polynomials of first and second degrees, respectively.The framelets ψ 3 [1] and ψ 2 [1] , which are impulse responses of the p-filters h 3 and h 2 , have three and two LDVMs, respectively.IRs of the p-filters h s , s = 0, 1, 2, 3, and their MRs are displayed in Fig. 6.4 (second from top).MRs of the p-filters h 2 and h 3 mirror each other.

Symmetric factorization of W
, which provides equal number (two) of LDVMs to the analysis and to the synthesis framelets, are: .
FRs of the p-filters hs and h s , s = 2, 3, are given in Eq. (6.1), where T
FRs of the p-filters hs and h s , s = 2, 3, are given in Eq. (6.1).The pairs of the p-filter banks H = {h s } , s = 0, 1, 2, 3, and H = h 0 , h 1 , h2 , h3 generate semi-tight frames.Figure 6.1 displays the impulse responses of the p-filters hs and h s , s = 2, 3, which are the discrete time framelets of the first level and their MR.

Quadratic quasi-interpolating spline (non-interpolating low-pass p-filter):
Frequency responses of the p-filter h 0 and h 1 are given in Eq. (5.17).The p-filter h 0 locally restores sampled cubic polynomials while the p-filter h 1 locally eliminates them.Thus, the framelet ψ 1 [1] has four LDVMs.

Antisymmetric factorization (semi-tight frame) of t[n]
, which assigns three LDVMs to the analysis framelet ψ2 [1] leaving only one LDVM to the synthesis framelet ψ 2 [1] and vice versa for the framelets ψ3 [1] and FRs of the p-filters hs and h s , s = 2, 3, are defined in Eq. (6.1). Figure 6.2 displays the IRs of the p-filters hs and h s , s = 2, 3, which are the discrete-time framelets of the first level with their MRs.Pseudo-spline: FRs of the non-interpolating low-pass p-filter h 0 and the corresponding high-pass p-filter h 1 are given by Eq. (5.22).The p-filter h 0 locally restores sampled cubic polynomials, while the p-filter h 1 locally eliminates sampled fifth-degree polynomials.Thus, the framelet ψ 1 [1] has six LDVMs.We have The following factorizations modes are possible: FRs of the p-filters h 2 and h 3 are defined in Eq. (6.2) with The p-filters h 2 and h 3 locally eliminate sampled polynomials of first and second degrees, respectively, thus, the framelets ψ 2 [1] and ψ 3 [1] have two and three LDVMs, respectively.

Antisymmetric factorization (semi-tight frame) t[n]
assigns three LDVMs to the analysis framelet ψ2 [1] leaving only one LDVM to : FRs of the p-filters hs and h s , s = 2, 3, are defined in Eq. (6.1).[k], s = 0, 1, 2, are symmetric, while the framelets ψ 3 [1] [k] are antisymmetric.The magnitude responses of the p-filters h 0 and h 1 mirror each other and the same is true for the band-pass pair h 2 and h 3 .In the caption of the figure, abbreviation QIS means quasi-interpolating spline.6.2.Four-channel p-filter banks with p-IIR p-filters.Unlike non-periodic setting, the implementation cost of IIR p-filters is no higher than the implementation cost of p-FIR p-filters.However, giving up the requirement of having finite impulse response provides additional flexibility in the design of p-filter banks with needed properties.
For the design of IIR filter, we use the interpolating low-pass p-filters h 0 .The frequency responses of the p-filters are structured as ĥ0 where f [n] 1 is the frequency response of some prediction p-filter.Prediction p-filters were derived from interpolating polynomial and discrete splines, were described in Sections 5.1 and 5.2, respectively.= ω 2n + 6 + ω −2n .The prediction p-filter f 3 c and the p-filters h 0 and h 1 are defined in Eq. (5.8).The sequence The p-filter h 0 locally restores sampled cubic polynomials while the p-filter h 1 locally eliminates them.Thus, the framelet FRs of the p-filters h 2 and h 3 are defined in Eq. (6.2) where The p-filter h 3 locally eliminates sampled quadratic polynomials, thus the framelet ψ 3 [1] , has three LDVMs, while the framelet ψ 2 [1] , has two LDVMs.The IRs of the p-filters h s , s = 0, 1, 2, 3, and their MRs are displayed in Fig.

(top). 2. Antisymmetric factorization (semi-tight frame
which assigns three LDVMs to the analysis framelet ψ2 [1] leaving only one  LDVM to the synthesis framelet ψ 2 [1] and vice versa for the framelets ψ3 [1] and ψ 3 [1] : FRs of the p-filters hs and h s , s = 2, 3, are defined in Eq. (6.1). Figure 6.5 displays IRs of the p-filters hs and h s , s = 2, 3, which are discrete-time framelets of the first level and their MRs.
Cubic interpolating spline, p = 4: The prediction p-filter f 4 c and the p-filters h 0 and h 1 are defined in Eq. (5.9).The p-filter h 0 locally restores the sampled cubic polynomials while the p-filter h 1 locally eliminates them.Thus, the framelet ψ 1 [1] has four LDVM.We have Comparing Eq. (6.4) with Eq. ( 6.3), we observe that the numerator of the sequences W [n] 1 in both equations is the same.Therefore, factorization of W [n] 1 of the cubic interpolating spline is similar to the factorization of the quadratic quasiinterpolating spline.Here are some factorization examples: FRs of the p-filters h 2 and h 3 are defined in Eq. ( 6.2) where The framelet ψ 3 [1] has three LDVMs, while the framelet ψ 2 [1] has two LDVMs.IRs of the p-filters h s , s = 0, 1, 2, 3, and their MRs are displayed in Fig. 6.8 (second from top).

Antisymmetric factorization (semi-tight frame) W
assigns three LDVMs to the analysis ψ2 [1] and one LDVM to the synthesis framelet ψ 2 [1] and vice versa for ψ3 [1] and ψ 3  [1] : The p-filter h 0 locally restores the sampled polynomials of the fifth degree while the p-filter h 1 locally eliminates them.Thus, the framelet ψ 1 FRs of the p-filters h 2 and h 3 are defined in Eq. (6.2) where [1] has four LDVMs and ψ 3 [1] has three LDVMs.IRs of the p-filters h s , s = 0, 1, 2, 3, and their MRs are displayed in Fig. 6.8 (third from top).

Symmetric factorization (semi-tight frame) W
assigns four LDVMs to the analysis ψ2 [1] and two LDVMs to the synthesis framelet ψ 2 [1] and vice versa for ψ3 [1] and ψ 3 [1] : , where all the framelets ψs [1] and ψ s [1] , s = 2, 3, have three LDVMs: In both the symmetric and the antisymmetric cases, FRs of the p-filters hs and h s , s = 2, 3, are defined in Eq. (6.1). Figure 6.7 displays IRs of the p-filters hs and h s , s = 2, 3, which are the discrete time framelets of the first level and their MR. where where The interpolating low-pass p-filter h 0 locally restores sampled polynomials of degree 2r − 1, while the high-pass p-filter h 1 eliminates them (the framelet ψ 1 [1] has 2r LDVMs).The p-filter bank H = {h s } , s = 0, 1, 2, 3, generates a tight frame.If r is even then the framelet ψ 2 [1] has r + 1 LDVMs, while the framelet ψ 3 [1] has r LDVMs, and vice versa for r odd.Display of tight frame IIR p-filters: Figure 6.8 displays the impulse and the magnitude responses of IIR four-channel p-filter banks that generate tight frames.The p-filters are derived from the polynomial interpolating splines of orders 3,4,5 and from the discrete splines of orders 6, 8, 10 and 12. Impulse responses of the p-filters h s are the corresponding discrete-time framelets ψ s [k] is antisymmetric.MRs of the p-filters h 0 and h 1 mirror each other and the same is true for the band-pass pair h 2 and h 3 .Note that, as the spline order is growing, the MRs shapes of the p-filters h 0 and h 1 are approaching rectangles while the MRs of the p-filters h 2 and h 3 are shrinking.In the caption of Fig. 6.8, IS means interpolating spline and DS means discrete spline.7. Experimental results.In this section, we present experimental results where the developed framelets are applied to solve two image processing problems.
1. Image restoration from an input, which was degraded by blurring aggravated by random noise and by random loss of significant number of pixels.The images are transformed by periodic frames designed in Section 6, which are extended to the 2D setting in a standard tensor product way.The goal of our experiments is to compare between the performances of different tight and semi-tight frames in identical conditions.2. Hyperspectral snapshot imaging.
7.1.Experiments on image restoration.In this section, we compare between the performance of different wavelet frames applied to a restoration task of degraded images.
7.1.1.Outline of the restoration scheme.The images are restored by the application of the split Bregman iteration (SBI) scheme [15] that uses the analysis based approach (see for example [30]).
Denote by u = {u[κ, ν]} the original image array to be restored from the degraded array f = K u + ε, where K denotes the 2D discrete convolution operator of the array u with a kernel k = {k[κ, ν]} and ε = {e k,n } is a random error array.K denotes the conjugate operator of K that implements the discrete convolution with Right pictures display the magnitude responses of these p-filters: h 0 and h 1 (dashed lines), h 3 (dash-dot line) and h 2 (solid line).Top: quadratic IS.Second from top: cubic IS.Third from top: four degree IS.Center: DS of order 6.Third from bottom: DS of order 8. Second from bottom: DS of order 10.Bottom: discrete spline of order 12 the transposed kernel k T .If some number of pixels are missing then the image u should be restored from the available data (7.1) where P Λ denotes the projection on the remaining pixels set.The solution scheme is based on the assumption that the original image u has a sparse representation in a frame domain.Denote by F the frame expansion operator of the image u where C An approximated solution to Eq. (7.1) is derived via the unconstrained minimization of the functional where • 1 and • 2 are the l 1 and l 2 norms of the sequences, respectively.If Denote by T ϑ the soft thresholding operator: Following [30], we solve the minimization problem in Eq. ( 7.2) by an iterative algorithm.We begin with the initialization u 0 = 0, d 0 = b 0 = 0.Then, ( The linear system in the first line of Eq. ( 7.3) is solved by the application of the conjugate gradient algorithm.The operations in the second and third lines are straightforward.The choice of the parameters λ and µ depends on experimental conditions.
7.1.2.Results.The restoration algorithms were applied to the "Window", "Barbara", "Lena" and "Fingerprint" images.These images were blurred by the application of convolution that uses either the motion kernel from Matlab or a Gaussian kernel.Then, the images are further degraded by random removal of a large number of pixels.In some experiments, the degradation was aggravated by adding zero mean random noise.
In the experiments, we compare between the performances of a number of tight frames (TF) and semi-tight frames (STF) using four-channel filter banks, which are designed in Section 6.For comparison, we also included the frame, which uses the three-channel filter banks originated from interpolating second-order spline ( [26])see Fig. 7.1.
The framelets in Table 1 are denoted as follows: The proximity between an image ũ and the original image u is evaluated visually and by the Peak-Signal-to-Noise ratio (PSNR) Restoration experiments for the "Window" image: This image was taken from [30].The colored image is of size of size 512 × 512 × 3 given in Red, Green and Blue wavebands.3), which were applied separately to each of the RGB components.The conjugate gradient solver used 100 iterations.The tight and semi-tight frames listed in Table 1 were tested.The decomposition was implemented down to the fifth level.The restored PSNR results are given in Table 2.The best PSNR results (42.16 dB) was achieved by using the semi-tight frames S 4 8,2 that were derived from the pseudo-spline using the antisymmetric factorization.The pseudo-spline tight frame T 4  8,0 produced almost the same result and results produced by the frames S 4  6,2 and S 4 6,0 , which were derived from the quadratic QIS, are very close.The performance of the tight frames T 4 1,0 and T 3  1,0 , which were derived from the piece-wise linear splines, were significantly worse as well as the performance of the tight frames T 4 100,0 and T 4 120,0 that have many LDVMs.Figure 7.2 displays the restoration result.Visually, the restored image hardly can be distinguished from the original one.2. PSNR results after the restoration of the "Window" image from a blurred input where of 30% of its pixels were randomly removed   4. PSNR results after the restoration of the "Barbara" image from a blurred input where of 50% of its pixels were randomly removed.
gradient solver used 10 iterations.Tight and semi-tight frames listed in Table 1 were tested.The restored PSNR results are given in Table 3. Decomposition is implemented down to the second level.As in the previous experiment, the best PSNR 27.674 dB result was achieved by the application of the four-channel semitight frame S 4  8,2 derived from the pseudo-spline using anti-symmetric factorization.Again, the pseudo-spline tight frame T 4  8,0 produced almost the same result and results produced by the frames S 4  6,2 and S 4 6,0 derived from the quadratic QIS are very close.However, in this experiment, performance of the majority of the remaining frames does not differ much from the performance of the above mentioned ones.It is not true for the tight frames T 4 100,0 and T 4 120,0 , which have many LDVMs.The restoration result, which uses the frame S 4  8,2 , is displayed in Fig. 7.3.We observe that the noise is completely removed and the image is deblurred.The colors in the image are restored satisfactorily.
Results from the restoration of the "Barbara" image.Experiment 2.1, Image blurred and pixels removal: In these experiment, the grayscale "Barbara" image was restored after it was blurred by a convolution with the Gaussian kernel (MATLAB function fspecial('gaussian',[5 5],5)) and its PSNR became 23.32 dB.Then, 50% of its pixels were randomly removed.This reduces the PSNR to 7.55 dB.Random noise was not added.The image was restored by 50 SBI using the parameters λ = 0.001, µ = 0.005 in Eq. (7.3).The conjugate gradient solver used 100 iterations.Tight and semi-tight frames listed in Table 1 were tested.The decomposition is implemented down to the second level.The restored PSNR results are given in Table 4. Unlike the experiments with the "Window" image, the best performance was achieved by using the tight frames T 4 100,0 and T 4 120,0 , which are derived from the interpolating discrete splines of tenth and twelfth orders, respectively.Recall that, for the frame T 4 100,0 , the framelet ψ  spline.The restoration result from the application of the frame T 4 100,0 is displayed in Fig. 7.4.We observe that fine texture of the image, which is undistinguishable in the blurred image, is restored accurately.Results from the restoration of the "Lena" image.Experiment 3.1, Image strongly blurred and distorted by curves: The source image was strongly burred by convolution with a Gaussian kernel (MATLAB function fspecial('gaussian', [12,12],12)) to get PSNR=23.64dB.The blurred image was distorted by randomly drawn curves to get PSNR=16.44 dB.Noise was not added.The image was restored by 70 SBI with the parameters λ = 0.03 and µ = 0.17 by the application of the four-channel frames listed in Table 1.The decomposition was implemented down to the fifth level.The conjugate gradient solver used 15 iterations.The PSNR results are given in Table 5.The best results PSNR=29.227 dB and PSNR=29.222dB were achieved by the application of the four-channel semitight frames S 4  8,2 and S 4 6,2 derived from the the pseudo-spline and from the quadratic quasi-interpolating spline, respectively.The respective tight frames T 4  8,0 and T 4 6,0 produced a little bit of a lower PSNR.However, the rest of the tested frames, except of T 4 100,0 and T 4 120,0 , performed similarly.5. PSNR results after the restoration of the "Lena" image from a blurred input distorted by randomly drawn curves result.We observe that the corrupting curves are completely removed and the image becomes deblurred.
Results from the restoration of the "Fingerprint" image.Experiment 4.1, Image blurred, strongly noised and pixels missing: The "Fingerprint" image was affected by a strong zero-mean white noise with STD σ = 20 after being blurred by convolution with the Gaussian kernel (MATLAB function fspecial('gaussian',[5 5],5)) to get PSNR=19.65 dB.Then, 50% of its pixels were randomly removed and this produced PSNR=9.05dB.The image was restored by 50 SBI with the parameters λ = 2.3 and µ = 0.45 using the tight and semi-tight  6. PSNR results after the restoration of the "Fingerprint" image from a blurred noised input where 50% of its pixels were randomly removed frames listed in Table 1.The decomposition is implemented down to the sixth level.The conjugate gradient solver used 11 iterations.The PSNR results are given in Table 6.Here, the distribution of the PSNR results is very similar to the distribution from the "Barbara" Experiment 2.1.The best PSNR=23.80 dB and PSNR=23.78 dB results were achieved by the application of the tight frames T 4 100,0 and T 4  4,0 derived from the tenth-order interpolating discrete spline and interpolating polynomial spline of fifth order.Close results were produced by the semi-tight frame S 4  4,2 and the tight frame T 4 120,0 .Figure 7.6 displays the restoration result "Window" and "Lena" images: For these images (Experiments 1.1, 1.2 and 3.1), which are relatively smooth, the best results were produced by using the semi-tight frames that were derived from the quasi-interpolating quadratic spline and the pseudo-spline.The low-pass p-filters in generating p-filter banks are non-interpolating and have a long IR (9 taps).Note that in all the three experiments, the semi-tight frames outperformed their tight counterparts.In these experiments, the tight frames T 4 100,0 and T 4 120,0 with high number of LDVMs, which are derived from interpolating discrete splines, produced significantly inferior results."Barbara" and "Fingerprint" images: Both images have fine texture, which are satisfactorily restored after strong blurring, corruption by noise ("Fingerprint") and after random removal of 50% of their pixels.For this pair of images, the distribution of the restoration results between frames are different.The best results were produced by using the tight frames T 4 100,0 , T 4 120,0 and T 4 4,0 .All these frames outperformed the frames that were derived from the pseudo-spline and from the quadratic QIS, which dominated the achievement of the best results in the Experiments 1.1,1.2 and 3.1.The frames T 4 1,0 and T 3  1,0 , which were derived from the piece-wise linear splines, demonstrated the worst performance.
All the three tight frames T 4 100,0 , T 4 120,0 and T 4 4,0 , originated from interpolating splines and the low-pass p-filters in generating filter banks, are interpolating.All the filter banks use p-filters whose IR are infinite.The framelets in all these frames have many LDVMs.In non-periodic setting, filtering with such filters is implemented in a recursive mode (see for example [5]).The transfer functions of non-periodic counterparts of the corresponding filter banks are rational functions, where numerators and denominators are the Laurent polynomials of high degrees.Therefore, for the transform's implementation, many recursive passes are needed that results in a high computational cost.
However, the computational cost of the recursive implementation of IIR filtering is comparable and, sometimes even lower than the computational cost of FIR filtering, where filters have similar properties.Some results of such a comparison are given in [35].
There is a well known belief that infinite IR cause the corresponding wavelets to have large support in comparison to the wavelets generated by FIR filters.However, the impulse responses of IIR filters with rational transfer functions and the corresponding wavelets decay exponentially as time grows ( [35]).Therefore, the effective support of these wavelets is compact.Figure 7.7 displays the impulse responses of the low-pass filters h 0 4,0 (IIR) and h 0 8,0 (FIR), which generate the frames T 4 4,0 and T 4 8,0 , respectively.Recall that the IIR filter h 0 4,0 locally restores fifth-degree polynomials, while the FIR filter h 0 8,0 locally restores only cubic polynomials.A great advantage of the periodic setting is that the computational cost of the transform's implementation does not depend on the spline order and on the number of LDVMs.There is no difference between the application of FIR and IIR pfilters.For example, the CPU time for implementing 50 Bregman iterations on the 'Fingerprint' image with simple FIR filter bank T 4  1,0 takes 15.71 seconds, while the same procedure with the complicated IIR filter bank T 4 120,0 takes 15.70 seconds.This fact proves the flexibility in choosing the proper frame.
There are some reservations in the image processing community about using periodic setting.There are claims that, while time-domain implementation of wavelet and frame transforms require some extension of signal/image beyond the boundaries, the periodic extension may produce undesirable artifacts near the boundaries.But this is not the case in our scheme.Actually, we do not extend the processed images.Filtering is implemented via the discrete Fourier transform that is perfectly invertible.No boundary artifacts are produced.This statement is illustrated in Fig. 7.8.This figure consists of three sections separated by vertical lines.In the top left section, column # 15 from the original "Barbara" image is displayed versus the column taken from the image restored by 50 SBI in Experiment 2.1.The two images in the bottom display fragments adjacent to the upper and lower edges of the column.In the central section, the same exposition was done for column # 370 from "Barbara".The right section displays column # 470 from the "Fingerprint" image, which was restored from blurred and severely noised input in Experiment 4.1.We see that in the "Barbara" columns, the restored lines coincide with the original lines at the edges.In the "Fingerprint" column, the discrepancy between the restored and the original curves at the edges is small and does not exceed the discrepancy at the internal points.No boundary effects persist.Both images were restored using the tight frame T 4 100,0 that was originated from the tenth-order discrete interpolating spline.
Increase above 50 of the number of SBIs does not contribute to the restoration quality.In some cases, as in the "Window" and "Barbara" experiments, a big number of iterations of the conjugate gradient solver increased the restoration quality.However, in the experiments with the strongly degraded "Lena" and "Fingerprint" images, increase in the number of iterations above 15 depleted the quality.
The choice of the regularization parameters λ and µ is of critical importance.The restoration quality is sensitive to them.Even a small change in these parameters significantly affects the restoration quality.The number of decomposition levels is also important.8. Snapshot spectral imaging.In this section, we briefly outline a potential application for spline-based wavelet frames: developing of a snapshot hyperspectral imaging (HSI) camera.HSI is commonly known as imaging with a large number of wavelengths in a given wavelength range.For example, the resolution range can be between 10 and hundreds or even thousands of wavelengths."Snapshot" HSI imagers perform simultaneous (instantaneous) acquisition of spatial and spectral data in a single snapshot.The acquired data is processed that results in a spectral cube?of a source object.The spectral cube includes light intensity data in two spatial dimensions and one spectral dimension and it is expressed as a three-dimensional (3D) matrix.

Inverse Problems and Imaging
Volume 9, No. 3 (2015), 661-707 An algorithm, which transforms a snapshot image taken by a regular digital camera into a spectral cube, is presented in [16] and it will be described in full details in one of our next papers.Here, we present the experimental results from a mathematical modeling of the imaging procedure.The snapshot spectral imager comprises of an imaging lens, a dispersed image sensor and a diffuser inserted in the optical path between a source image and a pixelated (as.e.g. in a digital camera) image sensor.The diffuser provides a dispersed-diffused image at the dispersed image sensor.A plurality of spectral images of the source object in different spectral bands is obtained in a single shot by reconstruction involving the dispersed-diffused image.The diffuser may be included internally in a digital camera.
Each pixel in a dispersed-diffused image includes a linear mixture of spectral and spatial information from surrounding pixels.A reconstruction of the spectral cube is performed using an optimization process aimed to compensate for the underdetermined nature of the problem.8.1.Mathematical model of the HSI snapshot.As source images we used the hyperspectral (HS) digital images from [14] that we downloaded from the cite The physical diffuser is a diffraction grating with a specific design of grooves.Its action on the optical flow is modeled by multiplication of the source matrix Z 0 with a block-wise Toepliz matrix A = (A 1 , A 2 , ..., A 33 ), where each sub-matrix A l is of size 1024×256.The produced matrix Y def = A•Z 0 = 33 l=1 A l •Z l 0 , of size 1024×256 is regarded as a dispersed-diffused image (DDI), which serves as an input to the restoration processing.The DDI is displayed in Fig. 8.2.8.2.Outline of the restoration scheme.Unlike the unconstrained minimization, which was utilized in Section 7 for restoration of corrupted images, we reconstruct the HS image from the DDI Y by solving the following constrained minimization problem: where F X denotes a 2D direct frame transform applied to blocks of a block-wise matrix X = X 1 , X 2 , ..., X 33 T : (8.2) F X = D = D 1 , D 2 , ..., D 33 T , D l def = F X l , l = 1, ..., 33.
The inverse frame transform applied to blocks of a block-wise matrix D is  (8.4).Different tight and semi-tight frames were tested.The best restoration results were achieved by the semi-tight S 4  6,2 originated from the quadratic quasi-interpolating spline.The frame transform was deconmposed down into the fourth level.Figure 8.3 displays the restored spectral cube converted to RGB color image versus the similarly converted original spectral cube.We observe that both images are very similar to each other and the colors are restored.The PSNR computed for all the pixels of the spectral cubes is 28.18 dB.Restoration of the spectra at four control points, which are marked on Fig. 8.3, is displayed in Fig. 8.4.We observe that the spectra at control points are satisfactorily restored.The diagram in Fig. 8.5 illustrates the restoration quality of different wavebands from the spectral cube.It shows the distribution of the PSNR values between the 33 wavebands.All the bands except band #1 (400÷410 nm), which is strongly noised, are satisfactorily restored.The highest PSNR=32.59 dB is observed at the restored band #19 (580 ÷ 590 nm).

2. 1 .
Signals and transforms.Notations: Π[N ] is the N -dimensional vector space consisting of N -periodic realvalued sequences x def = {x[k]}, k ∈ Z that are called signals.Throughout the paper, N = 2 j , j ∈ N and ω def = e 2πi/N .The inner product and the norm in Π[N ] are defined by x, h def signal δ[k], such that δ[0] = 1, l ∈ Z and δ[k] = 0 for k = 0, is called the Kronecker delta.Note that any discrete-time signal y of limited duration L can be embedded into a space Π[N ], N ≥ L, by the periodization ỹ def = ỹ[k] = l∈Z y[k + lN ] .The discrete Fourier transform (DFT) of a signal x and its inverse (IDFT) are x[n] = N −1 k=0 x[k] ω −nk and x[k] = 1/N N −1 n=0 x[n] ω nk , respectively, and they are calculated by the fast Fourier transform (FFT) algorithm.Since the signals are real-valued, the complex conjugate DFT becomes x[n] * = x[−n].

1 .Remark 3 . 1 .
The PR p-filter bank defined by Eq. (3.8) generates a semi-tight frame in the space Π[N ].The rational function t[n] 1 of ω n is symmetric about the change n → −n.Therefore, it can be factorized into product of two symmetric rational functions T [n] 1 and T [−n] 1 .An additional advantage of the semi-tight design is the option to swap approximation properties between analysis and synthesis framelets.3.1.2.Symmetric factorization of the matrix Q[n].If the inequality in Eq. (3.7) holds then the sequence t[n] 1 can be factorized as t[n] 1 = T [n] 1 T [−n] 1 and the following factorization of the matrix Q[n] is possible.Define, (3.9)

Definition 4 . 4 .
If a high(band)-pass p-filter g satisfies the conditions of Proposition 4.3, we say that the p-filter g locally eliminates sampled polynomials of degreem − 1.If a framelet ψ [1] [k] def = g[k], k ∈ Z, we say that the framelet ψ[1] has m local discrete vanishing moments (LDVM).

Remark 4 . 1 .
The statements of Propositions 4.1 and 4.3 concerning the low-pass h and high-pass g filters, respectively, remain true if the frequency responses are represented as ĥ which is the DFT of the sampled Bspline {B p (k)}, is strictly positive.Due to symmetry and positiveness of the Bspline, u p [n] 1 = u p [−n] 1 , thus, it is a cosine polynomial with positive coefficients.It is symmetric about K/2, where it reaches its single minimum.The maximal value of u p [n] 1 is 1, being reached when n = 0.

5. 1 . 2 .
Examples of spline based p-filters.The sequences {v p [n] 1 } and {u p [n] 1 }, whose ratio constitute the FR f p c [n] 1 of a prediction p-filter f p c , are calculated explicitly via the application of DFT to the B-splines sampled at the points {k +1/2} and {k}, respectively.The samples of the B-splines of any order are readily calculated from Eq. (5.1).

3 c
sin 2r πn/N cos 2r πn/N +sin 2r πn/N .All the p-filters have a linear phase.Except for the simplest case r = 1, the impulse responses of the p-filters are infinite.Nevertheless, they are well localized in time domain.Note, that the prediction p-filter f 4 d coincides with the p-filter f originating from the quadratic polynomial spline (see Eq. (5.8)).
is called the quadratic spline quasi-interpolating the signal x 0 ∈ Π[N/2].Denote the samples of the spline s 0 [k] def = S 3 (k) and s 1 [k]

Figure 6 .
Figure 6.3 displays the IRs of the p-filters h s , which are discrete-time framelets of the first level with their MR.Display of tight frame FIR p-filters: Figure 6.4 displays the impulse and the magnitude responses of the p-FIR four-channel p-filter banks that generate tight frames.The impulse responses of the p-filters h s are the corresponding discretetime framelets ψ s [1] [k], s = 0, 1, 2, 3, of the first decomposition level.The framelets ψ s [1][k], s = 0, 1, 2, are symmetric, while the framelets ψ 3[1] [k] are antisymmetric.The magnitude responses of the p-filters h 0 and h 1 mirror each other and the same is true for the band-pass pair h 2 and h 3 .In the caption of the figure, abbreviation QIS means quasi-interpolating spline.

Figure 6 . 4 . 2 . 1 .
Figure 6.4.Impulse and magnitude responses of the p-FIR fourchannel p-filter banks that generate tight frames.Left pictures display the impulse responses of p-filters where, from left to right we have: h 0 −→ h 1 −→ h 2 −→ h 3 .Right: pictures display the magnitude responses of these p-filters: h 0 and h 1 (dashed lines), h 3 (dash-dot line) and h 2 (solid line).Top: linear spline.Second from top: quadratic QIS (interpolating low-pass p-filter).Second from bottom: quadratic QIS (non-interpolating low-pass p-filter).Bottom: pseudo-spline

Figure 6 . 6 displays
Figure 6.6 displays IRs of the p-filters hs and h s , s = 2, 3, which are discrete-time framelets of the first level, and their MRs.Interpolating spline of fourth degree, p = 5: Denote Ω 5 c [n] 1

Figure 6 . 8 .
Figure 6.8.Impulse and magnitude responses of the IIR p-filter banks that generate tight frames.IS means interpolating spline and DS means discrete spline.Left pictures display the impulse responses of the p-filters.Left to right in the left picture:h 0 −→ h 1 −→ h 2 −→ h 3 .Right pictures display the magnitude responses of these p-filters: h 0 and h 1 (dashed lines), h 3 (dash-dot line) and h 2 (solid line).Top: quadratic IS.Second from top: cubic IS.Third from top: four degree IS.Center: DS of order 6.Third from bottom: DS of order 8. Second from bottom: DS of order 10.Bottom: discrete spline of order 12

def=
F u, C = {c[κ, ν]}, is the set of frame coefficients.Denote Inverse Problems and Imaging Volume 9, No. 3 (2015), 661-707 by F the operator of reconstruction of the image u from the set C. We get F C = u, and F F = I, where I is the identity operator.
k − xk ) 2 dB where M is the number of pixels in the image (in our experiments, M = 512 2 ), {x k } M k=1 are the original pixels of the image u and {x k } M k=1 are the pixels of the image ũ.

Figure 7 . 8 .
Figure 7.8.Left section: Top -original column # 15 from the"Barbara" image (dashed line) versus the column from the restored image.Bottom: left -left fragments of the lines; rightright fragments of the lines.Central section: The same for column # 370 from "Barbara".Right section: The same for column # 470 from the "Fingerprint" image

Figure 8 . 3 .
Figure 8.3.Left: Fragment of the original HS image converted to RGB of size 256 × 256 pixels; four control points are marked.Right: Restored fragment, PSNR=28.18 dB

Table 3 .
= 5, was added to each of the RGB matrices separately.The PSNR was 22.96 dB.Then, 30% of the pixels from each matrix were randomly removed from each color component.This reduces the PSNR to 9.42 dB.The image was restored by 40 SBI, using the parameters λ = 0.3, µ = 0.07 in Eq.(7.3).The conjugate PSNR results after the restoration of the "Window" image from a blurred noised input where of 30% of its pixels in each color component were randomly removed.PSNR 29.71 29.91 29.76 30.50 30.48 29.21 Table