Image upsampling by learning pairs of low-resolution dictionaries using a structured subspace model

Abstract

A computational method is disclosed for producing a sequence of high-resolution (HR) images from an input sequence of low-resolution (LR) images. The method uses a structured subspace framework to learn pairs of LR dictionaries from the input LR sequence ‘and’ employ learned pairs of LR dictionaries into estimating HR images. The structured subspace framework itself is based on a pair of specially structured HR basis matrices, wherein a HR basis spans any HR image whose so-called polyphase components (PPCs) are spanned by the corresponding LR dictionary.

Claims

1. A computer implemented method for upsampling low resolution (LR) frames to create high resolution (HR) versions thereof, comprising the steps of: estimating a sequence of said HR frames, from an input sequence of said LR frames, using a structured subspace framework that leverages pairs of LR dictionaries, respectively, wherein one member of each of said pairs of LR dictionaries is directly created, and a second member of each of said pairs is indirectly created, from the input sequence of LR frames to create estimated HR versions of said frames.

2. The method of claim 1, wherein each pair of said LR dictionaries is comprised of a first LR dictionary that is learned from the input sequence of LR frames, and a second, different-resolution, LR dictionary that is learned from a second sequence of LR frames obtained by decimating an estimate of the HR sequence.

3. The method of claim 2, wherein the structured subspace framework is based upon estimating at least one of two representations (α and β) of a HR image in terms of two HR bases (V and W), by solving
Vα=Wβ where V is the HR basis matrix that spans the HR subspace of all HR images whose polyphase components (PPCs) are spanned by atoms of the first LR dictionary; and W is a second HR basis matrix spanning the HR subspace of all HR images whose PPCs are spanned by the atoms of the second LR dictionary.

4. The method of claim 3, wherein said HR frames can be estimated by estimating patches thereof using patch-based processing, where the structured subspace framework is applied locally per patch to leverage a pair of local dictionaries learned per patch.

5. The method of claim 4, wherein the first UR dictionary is learned from local LR image patches extracted from the input sequence of LR frames, whereas the second LR dictionary is learned from local LR patches extracted from the second different-resolution LR sequence obtained by decimating the previous estimate of the HR sequence.

6. The method of claim 2, wherein the estimation of the HR sequence is repeated.

7. A computer implemented computational method for estimating a sequence of high-resolution (HR) frames from an input sequence of low-resolution (LR) frames, using a structured subspace framework that leverages pairs of LR dictionaries, comprising the steps of: scanning into a computer memory the input sequence of LR frames; creating a sequence of HR frames using said structured subspace framework by estimating representations of HR images in terms of pairs of HR bases, each of which embeds a different LR dictionary, by solving, per image, the equation:
Vα=Wβ where: V is the first HR basis matrix constructed such that it spans any HR image whose polyphase components (PPCs) are spanned by a first LR dictionary Ψ; W is the second HR basis matrix constructed such that it spans any HR image whose PPCs are spanned by another different resolution second LR dictionary Φ; α is the representation of the HR image in terms of V; and, β is the representation of the HR image in terms of W; multiplying the basis matrix V with the representation α, computed from the previous step, to obtain the estimate of the corresponding HR image.

8. A method as in claim 7, further comprising: dividing a frame, to be estimated, into patches, and applying the estimation process for each individual patch, independently of other patches.

9. A method as in claim 8, wherein: said patches have size d×d pixels, and overlap, vertically and horizontally, by pq pixels, where: d=rpq, p is the desired upsampling factor, q=p+1, and r is a user-defined factor that determines the size of patches.

10. A method in accord with claim 9, wherein: per patch, a first LR dictionary is learned from a first set of local LR image patches that are extracted from the input sequence of LR images, whereas a second LR dictionary is learned from a second set of local LR image patches that are extracted from the second different-resolution LR sequence obtained by decimating the previous estimate of the HR sequence, in accord with the steps of: extracting a first set of LR patches from a spatiotemporal neighborhood within the input LR sequence, where said neighborhood is centered at the spatiotemporal location of the LR patch whose HR version is to be estimated, and LR patches extracted from a LR frame, that belongs to said neighborhood, can overlap vertically and horizontally by at most rq−1 pixels; wherein the M patches, from the first set of extracted LR patches, that are most similar to the LR patch whose HR version to be estimated, are used as atoms of the first LR dictionary; extracting a second set of LR patches from a spatiotemporal neighborhood within the second different-resolution LR sequence, where said neighborhood is centered at the spatiotemporal location of the LR patch, in the second LR sequence, that corresponds to the HR patch to be estimated, and LR patches extracted from a LR frame, that belongs to said neighborhood, can overlap vertically and horizontally by at most rp−1 pixels; wherein the N patches, from the second set of extracted LR patches, that are most similar to the LR patch in the second LR sequence that corresponds to the HR patch to be estimated, are used as atoms of the second LR dictionary.

11. A method as in claim 9, wherein, an estimate of a frame is obtained by: tiling its estimated patches in their respective locations, forming r.sup.2 different estimates of the same HR frame; and, combining all r.sup.2 estimates into one single estimate.

12. The method of claim 11, where the estimated HR frames are decimated by q to create an updated version of the second LR sequence, and the solution is repeated, resulting in updated estimates of the sequence of HR frames until the HR frames stop getting any better in accord with a predetermined criteria defined as when a user-defined number of iteration is exhausted, and where the very first estimate of the HR sequence is an initial guess.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 shows unknown polyphase components (PPCs) of an unknown (HR) frame vs. a set of captured (LR) frames.

(2) FIG. 2 is an illustration of the relationships between PPCs corresponding to two different subsampling rates, p and q=p+1.

(3) FIG. 3 is a Flow chart illustrating a basic solution baseline based on the proposed structured subspace model. Note: dashed boxes/lines indicate initial guess operations (which are performed only once, before starting the iterative solution).

(4) FIG. 4 is a series of LR to HR image upsampling experiment (p=3).

(5) FIG. 5 is a series of LR to HR image upsampling experiment (p=4).

(6) FIG. 6 is a series of sample frames of the primary (available) LR sequence X.

(7) FIG. 7 is a series of sample frames of the initial guess of the HR sequence .Math..sub.(0).

(8) FIG. 8 is a series of sample frames of the first estimate of the secondary (different/lower resolution) LR sequence Y.sub.(1).

(9) FIG. 9 is the 31.sup.st frame from the primary (originally available) sequence X vs. the 31.sup.st frame from the first estimate of the secondary LR sequence Y.sub.(1). The frames are displayed in this figure with their actual relative sizes.

(10) FIG. 10 is frame #31 from the primary sequence X, partitioned in S=r.sup.2=16 different patterns. The sub-images within (imaginary) black borders are the LR patches {x.sup.kl.sup.s}, where k=31, s=1, . . . , 16, and total number of patches (per pattern s) is 48, i.e. l.sub.s=1, . . . , 48.

(11) FIG. 11 shows multiple HR estimates of frame #31, obtained by tiling all estimated HR patches {û.sub.(1).sup.kl.sup.s}, k=31, s=1, . . . , 16, and l.sub.s=1, . . . , 48.

(12) FIG. 12 shows û.sub.(1).sup.k (the first iteration estimate) displayed against the initial estimate û.sub.(0).sup.k (k=31), where û.sub.(1).sup.k is obtained by computing the median of all estimates shown in FIG. 11.

(13) FIG. 13 highlights an estimated HR patch û.sub.(1).sup.kl.sup.s at k=31, s=4, and l.sub.s=20.

(14) FIG. 14 is an illustration of the extraction process of a pair of local LR training sets (needed for computing the HR patch highlighted in FIG. 13).

(15) FIG. 15 is the pair of LR dictionaries Ψ.sup.kl.sup.s and Φ.sub.(i).sup.kl.sup.s we used to compute the HR estimate û.sub.(1).sup.kl.sup.s highlighted in FIG. 13.

DETAILED DESCRIPTION OF THE INVENTION

(16) Estimating HR images based on LR dictionaries is a notion that is quite alien to conventional wisdom. Indeed, image upsampling would be an intrinsically simple problem if all it takes to solve it is to “learn” LR dictionaries from a LR sequence, but could it really be that simple? The answer by the authors of [1], [2] (above) was ‘no’, unless the issue of partial measurements for the PPCs is resolved via an imaging hardware modification.

(17) Adding the structured subspace perspective, however, reveals that no imaging hardware modification is necessary, and that the problem can indeed be solved solely by learning LR dictionaries. Before we proceed with solution details, we would like first to list a few assumptions, acronyms, and notational conventions that will be used in describing the details. 1—Without loss of generality, we assume that the HR image, we seek to estimate, is square, and with a dimension that is an integer multiple of pq, where p is the primary decimation factor (equivalently, p is also the desired upsampling factor), and q is the secondary decimation factor, with the added assumption that
q=p+1.  (1) Namely, if we let u denote the HR image, then u is a matrix of size d×d, where
d=rpq,  (2) and r is a positive integer. 2—The acronym PPCs stands for polyphase components. PPCs associated with the primary decimation factor p are called primary PPCs or PPPCs for short. SPPCs is the acronym for secondary PPCs, which are associated with the secondary decimation factor q. 3—An underlined symbol denotes a (column) vector. Underlined symbols are also used to denote the vector form of matrices. To elaborate, vec(.Math.) is conventionally used to denote the vectorization operator, which turns a matrix into its (column) vector form (by lexicographical reordering of its elements). Therefore, if symbol A is used to denote a matrix, then A is shorthand for vec(A). 4—The (vector) subspace that contains all (real) vectors (or signals) with dimension n is denoted by R.sup.n. For example, writing u∈R.sup.d.sup.2 is simply saying that u (the HR image, in vector form), is in R.sup.d.sup.2, and thus has dimension d.sup.2. Similarly, R.sup.m×n denotes the subspace of all (real) matrices of size m×n. 5—We denote the identity matrix with the letter I. If a subscript is added, then it signifies the dimension of the identity matrix. For example, I.sub.n is the identity matrix with dimension n. Also, for an all zero matrix, we use the letter O. 6—The transpose of a matrix (or a vector) is always denoted by adding the superscript ‘*’. For example, if A is a matrix, then its transpose is denoted by A*. 7—A binary matrix A∈R.sup.m×n is a matrix whose entries are either “0” or “1”. A binary matrix A that is tall or square (m≥n), with each of its columns containing a sole “1” can be encoded (parameterized) as follows. Let π.sub.m.sup.A(k) denote the row-index of the sole “1” in the k-th column of A, then the set of all n row-indices π.sub.m.sup.A(1:n), completely encodes A (1:n is shorthand for k=1, 2, . . . , n). Note, if m=n, we drop the subscript m. For example,

(18) ${\pi }_4^A\left(1\text{:}3\right)=\left[\begin{array}{ccc}2& 3& 1\end{array}\right]\mathrm{encodes}A=\left[\begin{array}{ccc}0& 0& 1\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right],\mathrm{and}$${\pi }^B\left(1\text{:}3\right)=\left[\begin{array}{ccc}3& 2& 3\end{array}\right]\mathrm{encodes}B=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 1& 0& 1\end{array}\right].$ 8—We denote the Kronecker product and the Hadamard (element-wise) product, using the symbols .Math. and ⊙, respectively. Also, the so-called ‘direct sum’ is denoted by .Math.. In particular, the notation .Math..sub.k=1.sup.sA is exactly equivalent to I.sub.s.Math.A. In other words, if A∈R.sup.m×n, then .Math..sub.k=1.sup.sA is simply a matrix of size sm×sn that has s ‘replicas’ of A along the main diagonal (and zeros everywhere else). 9—We use the expression mod(a,b) as a shorthand for ‘a modulo b’. Also, GCD(a,b) is the ‘greatest common divisor of a and b’. 10—Finally, if c is a real number, the smallest integer greater than or equal to c (the ‘ceiling’ of c) is denoted by ┌c┐. On the other hand, └c┘ denotes the largest integer smaller than or equal to c (i.e., └c┘ is the ‘floor’ of c).

1. A HR SUBSPACE BASED ON A LR DICTIONARY

(19) Laying the foundation for our proposed structured subspace-based framework begins with asking (and answering) the following question: suppose we have available a LR dictionary that spans the PPCs of the sought after HR image, how do we construct the ‘basis matrix’ of the subspace of all HR images whose PPCs are spanned by said dictionary?

(20) The answer to that question lies in the basic fact that, instead of simply interlacing its PPCs, a HR image, u∈R.sup.d×d can alternatively be constructed from its PPCs, by adding up zero-filled and shifted versions of the PPCs. Said differently, the answer to that question begins by analytically expressing the HR image in terms of its PPCs.

(21) Specifically, let {ƒ.sub.i}.sub.i=1.sup.p.sup.2∈R.sup.rq×rq denote the set of all p.sup.2 PPPCs of the HR image (corresponding to the primary decimation factor p), then

(22) $\begin{array}{cc}u=\underset{i=1}{\overset{p^2}{.Math.}}{S}_{m_p\left(i\right)}\left(Z_p{f}_i{Z}_p^{*}\right)S_{n_p\left(i\right)}^{*}=\underset{i=1}{\overset{p^2}{.Math.}}\left(S_{m_p\left(i\right)}{Z}_p\right)f_i\left(Z_p^{*}{S}_{n_p\left(i\right)}^{*}\right),& \left(3\right)\end{array}$

(23) where Z.sub.p∈R.sup.d×rq zero-fills the columns by a factor of p (zero-filling ‘upsamples’ a vector by inserting p−1 zeros between any two elements, and p−1 zeros after the last element), post-multiplying with Z.sub.p* zero-fills the rows, S.sub.m.sub.p.sub.(i)∈R.sup.d×d circularly shifts by m.sub.p(i) rows down, and post-multiplying by S.sub.n.sub.p.sub.(i)*∈R.sup.d×d circularly shifts by n.sub.p(i) columns to the right. Both S.sub.m.sub.p.sub.(i) and S.sub.n.sub.p.sub.(i)* are therefore (cyclic) permutation matrices, where m.sub.p(i) and n.sub.p(i) are in correspondence with the location of the i-th PPPC ƒ.sub.i on the HR grid, and are given by

(24) $\begin{array}{cc}{m}_p\left(i\right)=.Math.\frac{i}{p}.Math.-1,n_p\left(i\right)=\mathrm{mod}\left(i-1,p\right),& \left(4\right)\end{array}$

(25) Using the fact that zero-filling matrices and circular shifting matrices are binary matrices with the property that each column contains a sole “1”, the row-indices of the “1” elements give us a full description of such matrices. In particular, the expression
Π.sub.d.sup.Z.sup.p(k)=1+(k−1)p,
k=1,2, . . . ,rq  (5)

(26) which encodes the zero-filling matrix Z.sub.p, gives us all we need to know to construct Z.sub.p. Specifically, equation (5) literally says that Z.sub.p is of size d×rq, and that the k-th column has a “1” located at row Π.sub.d.sup.Z.sup.p(k).

(27) Similarly, the shifting matrices S.sub.m.sub.p.sub.(i) and S.sub.n.sub.p.sub.(i) can be expressed as follows (without using the subscript d since S.sub.m.sub.p.sub.(i) and S.sub.n.sub.p.sub.(i) are square matrices)

(28) $\begin{array}{cc}{.Math.}^{S_{m_{\beta }\left(i\right)}}\left(k\right)=1+\mathrm{mod}\left(k+m_p\left(i\right)-1,d\right),k=1,2,\mathrm{.Math.},d.& \left(6\right)\\ {.Math.}^{S_{n_{\beta }\left(i\right)}}\left(k\right)=1+\mathrm{mod}\left(k+n_p\left(i\right)-1,d\right).k=1,2,\mathrm{.Math.},d.& \left(7\right)\end{array}$

(29) Using the fact that the matrix equation B=CAD* can be reshaped into the vector equation vec(B)=(D.Math.C)vec(A), we rewrite (3) to express u∈R.sup.d.sup.2 in terms of {ƒ.sub.i}.sub.i=1.sup.p.sup.2∈R.sup.r.sup.2.sup.q.sup.2 as follows.

(30) $\begin{array}{cc}\underset{\_}{u}=\underset{i=1}{\overset{p^2}{.Math.}}\left(\left(S_{n_p\left(i\right)}{Z}_p\right).Math.\left(S_{m_p\left(i\right)}{Z}_p\right)\right)\underset{\_}{f_i}.& \left(8\right)\end{array}$

(31) By the mixed-product property of the Kronecker product, we have
(S.sub.n.sub.p.sub.(i)Z.sub.p).Math.(S.sub.m.sub.p.sub.(i)Z.sub.p)=(S.sub.n.sub.p.sub.(i).Math.S.sub.m.sub.p.sub.(i))(Z.sub.p.Math.Z.sub.p).  (9)
Let
S.sub.p.sub.2.sup.i≙S.sub.n.sub.p.sub.(i).Math.S.sub.m.sub.p.sub.(i),  (10)
and
Z.sub.p.sub.2≙Z.sub.p.Math.Z.sub.p.  (11)
S.sub.p.sub.2.sup.i∈R.sup.d.sup.2.sup.×d.sup.2 therefore represents 2D (circular) shifting, by m.sub.p(i) rows down, and by n.sub.p(i) columns to the right, and Z.sub.p.sub.2∈R.sup.d.sup.2.sup.×r.sup.2.sup.q.sup.2 represents 2D zero-filling by a factor of p. Note that S.sub.p.sub.2.sup.i is a permutation matrix, since the Kronecker product of two permutation matrices (S.sub.m.sub.p.sub.(i) and S.sub.n.sub.p.sub.(i)) is a permutation matrix.

(32) Using (8)-(11), we have

(33) $\begin{array}{cc}\underset{\_}{u}=\underset{i=1}{\overset{p^2}{.Math.}}{S}_{p^2}^i{Z}_{p^2}\underset{\_}{f_i}.& \left(12\right)\end{array}$

(34) Note that (12) is the vector equation form of the matrix equation (3).

(35) Now, assume we have available a LR (primary) dictionary, Ψ∈R.sup.r.sup.2.sup.q.sup.2.sup.×M, with M linearly independent atoms as its columns, and that W spans the set of PPPCs, {ƒ.sub.i}.sub.i=1.sup.p.sup.2, and let {α.sub.i}.sub.i=1.sup.p.sup.2∈R.sup.M denote the corresponding representations in terms of Ψ, i.e.,
ƒ.sub.i=Ψα.sub.i for i=1,2, . . . ,p.sup.2.  (13)

(36) Equation (12) can therefore be rewritten as

(37) $\begin{array}{cc}\underset{\_}{u}=\underset{i=1}{\overset{p^2}{.Math.}}{S}_{p^2}^i{Z}_{p^2}\Psi {\underset{\_}{\alpha }}_i.& \left(14\right)\end{array}$

(38) Rewriting (14) in matrix-vector form, we get

(39) $\begin{array}{cc}\underset{\_}{u}=\left[\begin{array}{cccc}{S}_{p^2}^1{Z}_{p^2}\Psi & S_{p^2}^2{Z}_{p^2}\Psi & \mathrm{.Math.}& S_{p^2}^{p^2}{Z}_{p^2}\Psi \end{array}\right]\left[\begin{array}{c}{\underset{\_}{\alpha }}_1\\ {\underset{\_}{\alpha }}_2\\ \mathrm{.Math.}\\ {\underset{\_}{\alpha }}_{p^2}\end{array}\right].& \left(15\right)\end{array}$

(40) Or more concisely,

(41) $\begin{array}{cc}\underset{\_}{u}=V\underset{\_}{\alpha },& \left(16\right)\\ \mathrm{with}& \\ V\stackrel{\Delta }{=}\left[\begin{array}{cccc}{S}_{p^2}^1{Z}_{p^2}\Psi & S_{p^2}^2{Z}_{p^2}\Psi & \mathrm{.Math.}& S_{p^2}^{p^2}{Z}_{p^2}\Psi \end{array}\right],& \left(17\right)\\ \mathrm{and}& \\ \underset{\_}{\alpha }\stackrel{\Delta }{=}\left[\begin{array}{c}{\underset{\_}{\alpha }}_1\\ {\underset{\_}{\alpha }}_2\\ \mathrm{.Math.}\\ {\underset{\_}{\alpha }}_{p^2}\end{array}\right].& \left(18\right)\end{array}$

(42) Clearly, the matrix V∈R.sup.d.sup.2.sup.×p.sup.2.sup.M is the answer to the question we posed at the beginning of this section. In particular, if we pick any vector {acute over (α)}∈R.sup.p.sup.2.sup.M, then ú=V{acute over (α)} is a HR image with p.sup.2 PPPCs that correspond to p.sup.2 linear combinations {α.sub.i}.sub.i=1.sup.p.sup.2 of the M atoms of Ψ, where the coefficients of the linear combinations are stored in {acute over (α)}.

(43) In other words, if we let V denote the subspace of all HR images of dimension d.sup.2 whose PPPCs are in span(Ψ), then
V=span(V).  (19)

(44) Equation (17) gives us two options as to how to construct V. The first (needlessly computationally expensive) option is to construct the matrices, {S.sub.p.sub.2.sup.i}.sub.i=1.sup.p.sup.2 and Z.sub.p.sub.2 (using equations (4), (6), (7), (10) and equations (5), (11), respectively) and then perform p.sup.2 multiplications with Ψ.

(45) The second option is to simply perform the operations of zero-filling and shifting on the atoms of Ψ, without actually using any matrix operators. In other words, V can be constructed from Ψ without using a single floating point arithmetic. For example, the following Matlab code constructs V from Ψ without any calculations

(46) TABLE-US-00001 TABLE 1 Matlab code to construct V from Ψ % input: Psi (and r,p,q,M) d = r*p*q; Psi_vec2mat = reshape(Psi,r*q,r*q,M); V = zeros(d*d,p*p*M); zrs = zeros(d,d,M); c = 0; for i = 1:p  for j = 1:p   c = c+1;   atomZshift = zrs;   atomZshift(i:p:end,j:p:end,:) = Psi_vec2mat;   V(:, (c−1)*M+1:c*M) = reshape(atomZshift,d*d,M);  end end % output: V

(47) However, besides allowing for a formal description of V, the matrices {S.sub.p.sub.2.sup.i}.sub.i=1.sup.p.sup.2 and Z.sub.p.sub.2 pave the way for useful insights. To elaborate, note first that V can be factored as the product of two matrices as follows.
V=[S.sub.p.sub.2.sup.1Z.sub.p.sub.2S.sub.p.sub.2.sup.2Z.sub.p.sub.2 . . . S.sub.p.sub.2.sup.p.sup.2Z.sub.p.sub.2](⊕.sub.i=1.sup.p.sup.2Ψ).  (20)

(48) If we define
P≙[S.sub.p.sub.2.sup.1Z.sub.p.sub.2S.sub.p.sub.2.sup.2Z.sub.p.sub.2 . . . S.sub.p.sub.2.sup.p.sup.2Z.sub.p.sub.2],  (21)
and
V.sub.Ψ≙⊕.sub.i=1.sup.p.sup.2Ψ,  (22)
then
V=PV.sub.Ψ.  (23)

(49) Now, instead of wastefully computing the Kronecker product in (11), it can be verified (using the code (5) and equation (11)) that Z.sub.p.sub.2 is a binary matrix with the following code

(50) 0$\begin{array}{cc}{\Pi }_{d^2}^{Z_{p^2}}\left(k\right)=1+p\left(\left(.Math.\frac{k}{rq}.Math.-1\right)d+\mathrm{mod}\left(k-1,rq\right)\right),k=1,2,\mathrm{.Math.},r^2{q}^2.& \left(24\right)\end{array}$

(51) Similarly, the row-indices encoding the i-th (permutation) matrix, S.sub.p.sub.2.sup.i, are

(52) $\begin{array}{cc}{\Pi }^{S_{p^2}^i}\left(k\right)=1+\mathrm{mod}\left(\mathrm{mod}\left(n_p\left(i\right)d+m_p\left(i\right)+k-1,d^2\right),d\right)+\mathrm{mod}\left(d\left(.Math.\frac{k-1}{d}.Math.+n_p\left(i\right)\right),d^2\right),k=1,2,\mathrm{.Math.},d^2& \left(25\right)\end{array}$

(53) Now, let's turn our attention to the matrix P∈R.sup.d.sup.2.sup.×d.sup.2, defined in (21). P consists of p.sup.2 blocks, {S.sub.p.sub.2.sup.iZ.sub.p.sub.2}.sub.i=1.sup.p.sup.2, and can be encoded given (24) and (25). In particular, knowing the fact that a zero-filling (or ‘upsampling’) matrix acts as a ‘decimation’ matrix when we post-multiply with it, then all post-multiplying S.sub.p.sub.2.sup.i with Z.sub.p.sub.2 does is simply knock out (p.sup.2−1)r.sup.2q.sup.2 columns of S.sub.p.sub.2.sup.i. Said differently, Z.sub.p.sub.2, conveniently represented by its code (24), tells us which (of the remaining r.sup.2 q.sup.2) columns of S.sub.p.sub.2.sup.i are stored as the i-th block in P. Using all these facts together, one can deduce the following encoding of P.

(54) $\begin{array}{cc}{\Pi }^P\left(k\right)=1+\mathrm{mod}\left(\mathrm{mod}\left(n_p\left(\stackrel{\_}{i}\right)d+m_p\left(\stackrel{\_}{i}\right)+\stackrel{\_}{k}-1,d^2\right),d\right)+\mathrm{mod}\left(d\left(.Math.\frac{\stackrel{\_}{k}-1}{d}.Math.+n_p\left(\stackrel{\_}{i}\right)\right),d^2\right),\stackrel{\_}{k}=\mathrm{mod}\left(\ell -1,d^2\right)+1,\stackrel{\_}{i}=.Math.\frac{\ell }{d^2}.Math.,\ell =1+p\left(\left(.Math.\frac{\stackrel{\stackrel{\_}{\_}}{k}}{\mathrm{rq}}.Math.-1\right)d+\mathrm{mod}\left(\stackrel{\stackrel{\_}{\_}}{k}-1,\mathrm{rq}\right)\right)+\left(.Math.\frac{k}{r^2{q}^2}.Math.-1\right)d^2,\stackrel{\stackrel{\_}{\_}}{k}=\mathrm{mod}\left(k-1,r^2{q}^2\right)+1k=1,2,\mathrm{.Math.},d^2& \left(26\right)\end{array}$
where m.sub.p and n.sub.p are defined in (4).

(55) Equations (22), (23) and (26) offer an alternate formal and efficient description of how to construct V (in lieu of the computationally expensive option of using equations (4)-(7), (10), (11)). In particular, note that V.sub.Ψ (defined in (22)) is just a matrix that contains p.sup.2 replicas of Ψ on its diagonal (no calculations are required to construct V.sub.Ψ). Additionally, it can be proven that P is a permutation matrix and therefore V can be constructed simply by shuffling the rows of V.sub.Ψ according to Π.sup.P(k).

(56) By definition, each column and each row in a permutation matrix must have a single “1” element. In other words, if one can verify that the row-indices Π.sup.P (k) are all unique (no repetitions), then P is a permutation matrix. An easier approach to show that P is a permutation matrix, is to use the fact that a matrix A is a permutation matrix iff A is square, binary and orthogonal (A*A=I). Now, since P is square and binary, proving that P* P=I proves that P is a permutation matrix. In other words, in light of the block structure of P, we just need to show that

(57) $\begin{array}{cc}{\left(S_{p^2}^i{Z}_{p^2}\right)}^{*}{S}_{p^2}^j{Z}_{p^2}=\{\begin{array}{cc}I,& i=j\\ O,& i\ne j\end{array}& \left(27\right)\end{array}$

(58) The first part is trivial. In particular, S.sub.p.sub.2.sup.i is a permutation matrix, and therefore its columns are orthonormal. Obviously, deleting columns (by post-multiplying with Z.sub.p.sub.2), does not change the fact that the remaining columns are still orthonormal, i.e. (S.sub.p.sub.2.sup.iZ.sub.p.sub.2)*S.sub.p.sub.2.sup.iZ.sub.p.sub.2=1.

(59) As to the second part, recall that the function of zero-filling and shifting is to perform interlacing of PPCs using the sum in (12), which means that the non-zero elements of zero-filled and shifted PPCs cannot overlap. This implies that
(S.sub.p.sub.2.sup.iZ.sub.p.sub.2ƒ.sub.i)⊙(S.sub.p.sub.2.sup.jZ.sub.p.sub.2ƒ.sub.i)=O
.Math.(S.sub.p.sub.2.sup.iZ.sub.p.sub.2ƒ.sub.i)*S.sub.p.sub.2.sup.jZ.sub.p.sub.2ƒ.sub.i=0,   (28)
i.e. ƒ.sub.i*(Z.sub.p.sub.2*S.sub.p.sub.2.sup.i*S.sub.p.sub.2.sup.jZ.sub.p.sub.2)ƒ.sub.i=0. Now consider the special case where ƒ.sub.i=ƒ.sub.j=1 for some i≠j, where 1 is an all-one vector. In this case, we have 1*(Z.sub.p.sub.2*S.sub.p.sub.2.sup.i*S.sub.p.sub.2.sup.jZ.sub.p.sub.2)1=0, i.e. the sum of all elements of the matrix T=Z.sub.p.sub.2*S.sub.p.sub.2.sup.i*S.sub.p.sub.2.sup.jZ.sub.p.sub.2 must be zero, so what does that tell us about T? We know that the product of two permutation matrices is a permutation matrix and therefore S.sub.p.sub.2.sup.i*S.sub.p.sub.2.sup.j is a permutation matrix. Also, post-multiplying by Z.sub.p.sub.2 and pre-multiplying by Z.sub.p.sub.2* only deletes columns and rows, respectively. In particular, the product (S.sub.p.sub.2.sup.i*S.sub.p.sub.2.sup.j)Z.sub.p.sub.2 removes (p.sup.2−1)r.sup.2q.sup.2 columns of the permutation matrix S.sub.p.sub.2.sup.i*S.sub.p.sub.2.sup.j, which means that (S.sub.p.sub.2.sup.i*S.sub.p.sub.2.sup.j)Z.sub.p.sub.2 has only r.sup.2q.sup.2 binary rows, each of which containing a single “1”, with the remaining (p.sup.2−1)r.sup.2q.sup.2 rows being entirely zero-row vectors. Now, if we delete (p.sup.2−1)r.sup.2q.sup.2 rows of S.sub.p.sub.2.sup.i*S.sub.p.sub.2.sup.jZ.sub.p.sub.2, by pre-multiplying with Z.sub.p.sub.2*, we get T, and thus T can be either a binary matrix or an all zero matrix, depending on which rows are deleted from S.sub.p.sub.2.sup.i*S.sub.p.sub.2.sup.jZ.sub.p.sub.2 to get T. However, since the sum of all elements of T must be zero, then T cannot be a binary matrix, which leaves us with only one possibility, i.e. Z.sub.p.sub.2*S.sub.p.sub.2.sup.i*S.sub.p.sub.2.sup.jZ.sub.p.sub.2=O.

2. THE INTERSECTION SUBSPACE OF A PAIR OF HR SUBSPACES BASED ON A PAIR OF LR DICTIONARIES

(60) Let {g.sub.i}.sub.i=1∈R.sup.r.sup.2.sup.p.sup.2 denote the set of all q.sup.2 SPPCs of the HR image u (corresponding to the secondary decimation factor q), and assume we have available a (secondary) LR dictionary, Φ∈R.sup.r.sup.2.sup.p.sup.2.sup.×N, with N linearly independent atoms as its columns, and that Φ spans the set of SPPCs, i.e.
g.sub.i=Φβ.sub.i for i=1,2, . . . ,q.sup.2,  (29)
where {β.sub.i}.sub.i=1.sup.q.sup.2∈R.sup.N denote the corresponding representations of {g.sub.i}.sub.i=1.sup.q.sup.2 in terms of Φ.

(61) Also, let W denote the subspace of all HR images of dimension d.sup.2 whose SPPCs are in span(Φ), and let W∈R.sup.d.sup.2.sup.×q.sup.2.sup.N be the basis matrix of W, i.e.
W=span(W).  (30)

(62) Following the analysis presented in the previous section, we can construct W from Φ, simply by exchanging p and q. In particular, the following equations define W.
W=QW.sub.Φ,  (31)
where
W.sub.Φ=⊕.sub.i=1.sup.q.sup.2Φ,  (32)
and Q∈R.sup.d.sup.2.sup.×d.sup.2 is a permutation matrix with row-indices given by

(63) $\begin{array}{cc}{\Pi }^Q\left(k\right)=1+\mathrm{mod}\left(\mathrm{mod}\left(n_q\left(\stackrel{\_}{i}\right)d+m_q\left(\stackrel{\_}{i}\right)+\stackrel{\_}{k}-1,d^2\right),d\right)+\mathrm{mod}\left(d\left(.Math.\frac{\stackrel{\_}{k}-1}{d}.Math.+n_q\left(\stackrel{\_}{i}\right)\right),d^2\right),& \left(33\right)\\ \overline{k}=\mathrm{mod}\left(-1,d^2\right)+1,\stackrel{\_}{\mathrm{.Math.}}=.Math.\frac{d^2}{}.Math.,\ell =1+q\left(\left(.Math.\frac{\stackrel{\stackrel{\_}{\_}}{k}}{rp}.Math.-1\right)d+\mathrm{mod}\left(\stackrel{\stackrel{\_}{\_}}{k}=-1,rp\right)\right)+\left(.Math.\frac{k}{r^2{p}^2}.Math.-1\right)d^2,\stackrel{\stackrel{\_}{\_}}{k}=\mathrm{mod}\left(k-1,r^2{p}^2\right)+1k=1,2,\mathrm{.Math.},d^2\mathrm{where}& \\ m_q\left(\stackrel{\_}{i}\right)=.Math.\frac{\stackrel{\_}{i}}{q}.Math.-1,n_q\left(\stackrel{\_}{i}\right)=\mathrm{mod}\left(\stackrel{\_}{i}-1,q\right).& \left(34\right)\end{array}$

(64) Now, similarly to (16), the HR image can be expressed as

(65) $\begin{array}{cc}\underset{\_}{u}=W\underset{\_}{\beta }.\mathrm{where}& \left(35\right)\\ \underset{\_}{\beta }\stackrel{\Delta }{=}\left[\begin{array}{c}{\underset{\_}{\beta }}^1\\ {\underset{\_}{\beta }}^2\\ \mathrm{.Math.}\\ {\underset{\_}{\beta }}_{q^2}\end{array}\right].& \left(36\right)\end{array}$

(66) Now we are ready to ask the following question: Could it be possible that V and W are all we need to find u? If the answer is yes, then, effectively, the pair of LR dictionaries Ψ and Φ are all that is required to find the HR image.

(67) To answer that question, we begin by first combining equations (16) and (35) into one equation
Vα=Wβ,  (37)
which we rewrite as

(68) $\begin{array}{cc}\left[\begin{array}{cc}V& W\end{array}\right]\left[\begin{array}{c}\underset{\_}{\alpha }\\ -\underset{\_}{\beta }\end{array}\right]=\underset{\_}{O}.& \left(38\right)\end{array}$

(69) Equation (38) suggests that we need to study the nullspace of the augmented matrix [V W]∈R.sup.d.sup.2.sup.×(p.sup.2.sup.M+q.sup.2.sup.N). Before we explain, let Z denote the nullspace of [V W], and let Z∈R.sup.(p.sup.2.sup.M+q.sup.2.sup.N)×z be its basis, where z≙dim(Z), i.e.
Z≙null([V W]),  (39)
and
Z=span(Z),  (40)
and assume that the sought after HR image is large enough. Specifically, assume
d.sup.2=r.sup.2p.sup.2q.sup.2≥(p.sup.2M+q.sup.2N).  (41)

(70) Now, define
U≙V∩W,  (42)
and let U denote the basis for the intersection subspace U, i.e.
U=span(U).  (43)

(71) Since u∈V and u∈W, then u∈U, which means that if dim(U)=1, then finding U is tantamount to finding (a scaled version of) u. However, thus far, all we know about the dimensionality of the intersection subspace is
1≤dim(U)≤min(p.sup.2M,q.sup.2N),  (44)
where the lower bound is due to the fact that u∈U, while the upper bound is due to (42).

(72) Clearly, to gain additional insight, we need to turn our attention to the nullspace. In particular, let Z be partitioned into a column of two matrices, Z.sub.V∈R.sup.p.sup.2.sup.M×z and Z.sub.W∈R.sup.q.sup.2.sup.N×z, i.e.

(73) $\begin{array}{cc}Z=\left[\begin{array}{c}{Z}_V\\ Z_W\end{array}\right],\mathrm{then}& \left(45\right)\\ U=V{Z}_V,& \left(46\right)\end{array}$

(74) Equivalently,
U=WZ.sub.W.  (47)

(75) Therefore, if V has full rank p.sup.2M, or, equivalently, W has full rank q.sup.2N, then,

(76) $\begin{array}{cc}\begin{array}{c}z\stackrel{\Delta }{=}\dim \left(\mathrm{��}\right)=\mathrm{rank}\left(Z\right)=\mathrm{rank}\left(Z_V\right)=\mathrm{rank}\left(Z_W\right)\\ =\mathrm{rank}\left(U\right)=\dim \left(\mathrm{��}\right).\end{array}& \left(48\right)\end{array}$

(77) To verify that V is full rank, recall equations (22), (23) and that P*P=I, hence
V*V=⊕.sub.i=1.sup.p.sup.2(Ψ*Ψ)  (49)

(78) Therefore, rank(V)=rank(V*V)=p.sup.2rank(Ψ*Ψ)=p.sup.2rank(Ψ)=p.sup.2M. Similarly,
W*W=⊕.sub.i=1.sup.q.sup.2(Φ*Φ),  (50)
and thus rank(W)=q.sup.2N.

(79) Now equation (48) says that the dimensionality of the intersection subspace is equal to the nullity of the augmented matrix [V W]. Of course, the nullity of a matrix is simply the multiplicity of its zero singular values, so the question (which we posed before equation (37)) boils down to: what are the necessary and sufficient conditions for [V W] to have no more than 1 zero-singular value?

(80) The first three necessary conditions are obvious and are already satisfied. Namely, we obviously need the augmented matrix to have more rows than columns (41). Also, the columns of V must be linearly independent (rank(V)=p.sup.2M). Similarly, the columns of W must be linearly independent as well (rank(W)=q.sup.2N).

(81) The fourth necessary condition is also obvious and it pertains to the span of each of the pair of LR dictionaries Ψ and Φ. In particular, let {{tilde over (ƒ)}.sub.i}.sub.i=1.sup.p.sup.2 and {{tilde over (g)}.sub.i}.sub.i=1.sup.q.sup.2 denote the two sets of PPPCs and SPPCs, respectively, of a different HR image ũ≠cu, c∈R (namely, ũ is a different image and it is not simply a scaled version of u). If {{tilde over (ƒ)}.sub.i}.sub.i=1.sup.p.sup.2∈span(Ψ), then we know that ũ∈V. If, additionally, {{tilde over (g)}.sub.i}.sub.i=1.sup.q.sup.2∈span(Φ), then ũ∈W, and therefore, ũ∈W. But we already know that u ∈U, and therefore, in this case, dim(U)≥2. In other words, the fourth necessary condition necessitates that the pair of LR dictionaries cannot, simultaneously, accommodate the PPCs of any HR image other than u.

(82) Are there any more conditions? Ostensibly, the answer would require the daunting task of examining the very sparse structure of both V and W; particularly the interaction of both sparse structures within the augmented matrix [V W]. Instead, we derive equivalent forms of (38) to see if the associated matrices are more revealing.

(83) We proceed with the pre-multiplication of both sides of (38) with [V W]*, to obtain

(84) $\begin{array}{cc}\left[\begin{array}{cc}{V}^{*}V& V^{*}W\\ W^{*}V& W^{*}W\end{array}\right]\left[\begin{array}{c}\underset{\_}{\alpha }\\ -\underset{\_}{\beta }\end{array}\right]=\underset{\_}{O}.& \left(51\right)\end{array}$

(85) For convenience, and without loss of generality, we assume that atoms of a LR dictionary are orthonormal, i.e.
Ψ*Ψ=I,  (52)
and
Φ*Φ=I.  (53)

(86) Using (49) and (52) we get
V*V=I,  (54)

(87) Similarly (50), (53),
W*W=I.  (55)

(88) Consequently, equation (51) is rewritten as

(89) 0$\begin{array}{cc}\left[\begin{array}{cc}I& V^{*}W\\ W^{*}V& I\end{array}\right]\left[\begin{array}{c}\underset{\_}{\alpha }\\ -\underset{\_}{\beta }\end{array}\right]=\underset{\_}{O}.& \left(56\right)\end{array}$

(90) Now the top part of (56) gives us the equation
α=V*Wβ,  (57)
and the bottom part reveals that
β=W*Vα.  (58)

(91) Plugging (58) in (57), we get
(I−V*WW*V)α=O.  (59)

(92) Note that V*WW*V=(W*V)*W*V, and therefore, V*WW*V is a symmetric positive semidefinite (PSD) matrix with singular values between 0 and 1 (the upper limit on the largest singular value is due to (54) and (55)). Likewise, the matrix (I−V*WW*V) is also symmetric, PSD, with singular values between 0 and 1. Moreover, since (59) is derived from (51), which in turn is derived from (38), then Z.sub.V=null(I−V*WW*V). In other words, dim(U) is equal to the multiplicity of the smallest (“0”) singular value of (I−V*WW*V). Equivalently, dim(U) is equal to the multiplicity of largest (“1”) singular value of V*WW*V, which is the same multiplicity of the largest singular value of W*V, i.e.
dim(U)=multiplicity of σ.sub.1(W*V),  (60)
where σ.sub.1(.Math.) denotes the largest singular value of a matrix.

(93) So let's take a look at the structure of W*V and see if it could reveal anything about the multiplicity of its singular values (particularly, the largest one). Given equations (22), (23) and equations (31), (32) we can write
W*V=(⊕.sub.j=1.sup.q.sup.2Φ*)Q*P(⊕.sub.j=1.sup.p.sup.2Ψ).  (61)

(94) If we only consider the product (⊕.sub.j=1.sup.q.sup.2Φ*)(⊕.sub.j=1.sup.p.sup.2Ψ), then it can be proven (with relative ease) that GCD(p,q)=1 is a necessary condition for the largest singular value to be distinct. Unfortunately, the presence of the permutation matrix Q*P in the middle between (⊕.sub.j=1.sup.q.sup.2Φ*) and (⊕.sub.j=1.sup.p.sup.2Ψ) complicates the picture. Yet, the same necessary condition can still be observed to hold for σ.sub.i(W*V) to be distinct.

(95) For practical reasons that shall become apparent later, we are only interested in the case q=p+1, which automatically satisfies GCD(p,q)=1. Said differently, we are only interested in what the structure of W*V can reveal about the multiplicity of σ.sub.i (W*V) given the assumption (1). Before we proceed, let us first make the following definitions.

(96) Let A.sub.j.sup.i∈R.sup.r.sup.2.sup.×M and B.sub.j.sup.i∈R.sup.r.sup.2.sup.×N denote submatrices of Ψ, and Φ, respectively, such that
A.sub.j.sup.i=D.sub.q.sub.2R.sub.k.sub.q.sub.(i),k.sub.q.sub.(j)Ψ,  (62)
where D.sub.q.sub.2∈R.sup.r.sup.2.sup.×r.sup.2.sup.q.sup.2 represents 2D decimation by a factor of q, and R.sub.k.sub.q.sub.(i),k.sub.q.sub.(j)∈R.sup.r.sup.2.sup.q.sup.2.sup.×r.sup.2.sup.q.sup.2 represents 2D shifting by k.sub.g(i) rows up, and k.sub.g(j) columns to the left, where
k.sub.g(i)=mod(q−i,q).  (63)

(97) Said differently, if we obtain the 2D form of each atom (column) in W, and decimate (by q) each (2D) atom starting at the entry located in column k.sub.g(j)+1 and row k.sub.g(i)+1, and then vectorize the decimated atoms and stack them next to each other, we get

(98) Similarly,
B.sub.j.sup.i=D.sub.p.sub.2R.sub.k.sub.p.sub.(i),k.sub.p.sub.(j)Φ,  (64)
where D.sub.p.sub.2∈R.sup.r.sup.2.sup.×q.sup.2.sup.×r.sup.2.sup.q.sup.2 represents 2D decimation by a factor of p, and R.sub.k.sub.p.sub.(i),k.sub.p.sub.(j)∈R.sup.r.sup.2.sup.p.sup.2.sup.×r.sup.2.sup.p.sup.2 represents 2D shifting by k.sub.p(i) rows up, and k.sub.p(j) columns to the left, where
k.sub.p(i)=mod(p−i,p).  (65)

(99) Now, it can be verified that when (1) is satisfied, W*V∈R.sup.q.sup.2.sup.N×p.sup.2.sup.M block-Toeplitz matrix, i.e.

(100) $\begin{array}{cc}{W}^{*}V=\left[\begin{array}{ccccc}{T}_0& T_{-1}& T_{-2}& \mathrm{.Math.}& T_{1-p}\\ T_1& T_0& T_{-1}& \mathrm{.Math.}& T_{2-p}\\ T_2& T_1& T_0& \mathrm{.Math.}& T_{3-p}\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ T_{q-1}& T_{q-2}& T_{q-3}& \mathrm{.Math.}& T_0\end{array}\right],& \left(66\right)\end{array}$
with each block T.sub.i∈R.sup.qN×pM,1−p≤i≤q−1, being a q×p block-Toeplitz matrix itself,

(101) $\begin{array}{cc}{T}_i=\left[\begin{array}{ccccc}{C}_0^i& C_{-1}^i& C_{-2}^i& \mathrm{.Math.}& C_{1-p}^i\\ C_1^i& C_0^i& C_{-1}^i& \mathrm{.Math.}& C_{2-p}^i\\ C_2^i& C_1^i& C_0^i& \mathrm{.Math.}& C_{3-p}^i\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ C_{q-1}^i& C_{q-2}^i& C_{q-3}^i& \mathrm{.Math.}& C_0^i\end{array}\right],& \left(67\right)\end{array}$
where the j-th sub-block C.sub.j.sup.i∈R.sup.N×M,1−p≤j≤q−1, is given by
C.sub.j.sup.i=B.sub.j.sup.i*A.sub.j.sup.i,  (68)
and A.sub.j.sup.i and B.sub.j.sup.i are as defined in (62)-(65).

(102) In other words, when (1) is satisfied it can be verified that W*V is a Toeplitz-block-Toeplitz (a 2-level Toeplitz) matrix generated by the (p+q−1).sup.2=4p.sup.2 sub-blocks {C.sub.j.sup.i}.sub.i,j=1−p.sup.q−1 as detailed by equations (62)-(68). Moreover, equations (62)-(65), and (68) tell us that all these sub-blocks are unique (i.e. C.sub.j.sup.i≠C.sub.j.sup.ī if either ī≠i or j≠j).

(103) So now, instead of wondering about other conditions we might need for the (sparse) augmented matrix of (38) to have a nullity of exactly 1, the question becomes: given that q=p+1, are the aforementioned necessary conditions all we need for the (dense) Toeplitz-block-Toeplitz matrix W*V to have a distinct largest singular value?

(104) The relevant literature is rich with analysis and properties pertaining to Toeplitz matrices. However, to the best of our knowledge, when it comes to the question of the multiplicity of singular values of Toeplitz-structured matrices, the answer can only be given in the context of asymptotic distributions of singular values (and under assumptions that do not apply in our case). Therefore, we state the following empirically verifiable result.

(105) Empirical Result 1: Given that the basic assumption (1) is satisfied, i.e. q=p+1, the following conditions are both necessary and sufficient for σ.sub.i (W*V) to be distinct.

(106) $\text{1-}r\ge .Math.\frac{\sqrt{\left(p^{2}M+q^{2}N\right)}}{pq}.Math..$ This is a reformulation of (41) in light of the fact that r must be an integer. 2—rank(Ψ)=M (the atoms of Ψ must be linearly independent for columns of V to be linearly independent). 3—rank(Φ)=N (the atoms of Φ must be linearly independent for columns of W to be linearly independent). 4—At most, the pair (Ψ,Φ) simultaneously accommodate the PPCs of only one image.

(107) Note that if the pair (Ψ,Φ) simultaneously accommodate the PPCs of exactly one image, then σ.sub.i(W*V)=1, and, equivalently, dim(1)=1.

3. THE NEED FOR AN ARBITRATOR

(108) In the previous two sections, we have seen that the intersection subspace of two HR subspaces based on two (different-resolution) LR dictionaries is at most 1-dimensional (if the four conditions listed above are satisfied). Relying on this knowledge, we now turn our attention to finding the HR image.

(109) As is well-known if dim(U)=1, the easiest (and most computationally efficient) way to find a solution to (37) is to simply compute the right singular vector associated with the smallest (“0”) singular value of [V W]. Similarly, when (52) and (53) are satisfied, the solution to the alternative formulation (59) is the singular vector of (I−V*WW*V) associated with its smallest (“0”) singular value, which is also the right singular vector associated with the largest (“1”) singular value of W*V.

(110) In practice, however, the pair of LR dictionaries are ‘learned’ from training (example) LR images and thus they would never simultaneously exactly accommodate the PPCs of the sought after HR image ū (or any other image for that matter). In other words, the intersection subspace is always empty if the conditions of Empirical Result 1 are satisfied (since, in practice, the pair of LR dictionaries do not, simultaneously, exactly accommodate the PPCs of any image).

(111) Nevertheless, since σ.sub.i(W*V)<1 is distinct, the optimal (in the Frobenius norm sense) approximate solution to (59), remains the right singular vector of W*V associated with σ.sub.i(W*V). Alternatively, since the smallest singular value of [V W] is distinct iff σ.sub.i(W*V) is distinct, then the optimal solution of (37) remains the right singular vector of [V W] associated with its smallest singular value.

(112) Of course, our ultimate goal of solving (37) (or its equivalent formulation (59)) is to estimate u. However, an optimal solution to (37) does ‘not’ necessarily guarantee the optimality of the estimation of u. To explain with an (extreme) example, let u.sub.1 and u.sub.2 be two images such that

(113) $.Math.\frac{\underset{\_}{u}}{.Math.\underset{\_}{u}.Math.}-\frac{{\underset{\_}{u}}_1}{.Math.{\underset{¯}{u}}_1.Math.}.Math.⪡.Math.\frac{\underset{\_}{u}}{.Math.\underset{\_}{u}.Math.}-\frac{{\underset{\_}{u}}_2}{.Math.{\underset{\_}{u}}_2.Math.}.Math.,$
i.e. u.sub.1 is a much better approximation of u than is u.sub.2, and suppose that the available pair of LR dictionaries simultaneously accommodate the PPCs of u.sub.2 exactly, while the same pair of dictionaries only approximate the PPCs of u.sub.1 (albeit very well). In this case, the exact solution of (37) will lead to the inferior approximation u.sub.2.

(114) Again, in practice, the pair of (learned) dictionaries never do simultaneously accommodate the PPCs of any image, but the previous example helps to highlight the fact that solving (37) can be completely blind to whatever constitutes an optimal estimation of u. Put differently, if we are to rely entirely on (37) (to estimate the representation of the HR image in terms of the HR bases), then the learning process of the pair of LR dictionaries must be carefully designed to try to avoid picking inferior approximations of our sought after HR image.

(115) Before we suggest an alternative to requiring a painstaking learning process of the pair of LR dictionaries, we first note that the solution to the following optimization problem

(116) $\begin{array}{cc}\underset{\underset{\_}{\alpha },\underset{\_}{\beta }}{\min }{.Math.V\underset{\_}{\alpha }-W\underset{\_}{\beta }.Math.}^{2}s.t.{.Math.\underset{¯}{\alpha }.Math.}^2=1& \left(69\right)\end{array}$
is the same optimal solution to (59), i.e. the solution to (69) is also the right singular vector of W*V associated with σ.sub.i(W*V).

(117) Now consider the following modified optimization problem

(118) $\begin{array}{cc}\underset{\underset{\_}{\alpha },\underset{\_}{\beta }}{\min }\left({.Math.V\underset{\_}{\alpha }-W\underset{¯}{\beta }.Math.}^2+\mu {.Math.\frac{1}{p^2}\underset{i=1}{\overset{p^2}{.Math.}}\Psi {\underset{\_}{\alpha }}_i-\underset{\_}{\stackrel{\_}{f}}.Math.}^2\right)s.t.{.Math.\underset{¯}{\alpha }.Math.}^2=c>0,& \left(70\right)\\ \mathrm{where}\underset{\_}{\stackrel{\_}{f}}=\frac{1}{{\rho }^2}\underset{i=1}{\overset{p^2}{.Math.}}{\underset{\_}{f}}_i,& \end{array}$
where ƒ=1/p.sup.2∈Σ.sub.i=1.sup.p.sup.2ƒ.sub.i, the mean of the PPPCs of the HR image, is assumed to be known (for now), and μ is a control parameter. By adding the term

(119) ${.Math.\frac{1}{p^2}\underset{i=1}{\overset{p^2}{.Math.}}\Psi {\underset{\_}{\alpha }}_i-\underset{\_}{\stackrel{\_}{f}}.Math.}^2,$
we require the solution to be such that the mean of the estimated PPPCs be close to ƒ, which has the effect of solution ‘arbitration’ in favor of u.

(120) Of course, ƒ is unknown, but what we do have in practice is the LR image, whose HR version we seek, and which we denote by x. Now since the (unknown) PPPCs of the (unknown) HR image would be highly correlated with x, then their unknown mean, ƒ, would also be highly correlated with x. Moreover, if x is well approximated by Ψ, i.e. ∥x−Ψ(Ψ*x)∥≈0, then we can use the following expression for arbitration in favor of u,

(121) ${.Math.\frac{1}{p^2}\underset{i=1}{\overset{p^2}{.Math.}}\Psi {\underset{\_}{\alpha }}_i-\Psi \left({\Psi }^{*}\underset{\_}{x}\right).Math.}^2={.Math.\Psi \left(\frac{1}{p^2}\underset{i=1}{\overset{p^2}{.Math.}}{\underset{¯}{\alpha }}_i-{\Psi }^{*}\underset{¯}{x}\right).Math.}^2={.Math.\frac{1}{p^2}\underset{i=1}{\overset{p^2}{.Math.}}{\underset{¯}{\alpha }}_i-{\Psi }^{*}\underset{¯}{x}.Math.}^2={.Math.\underset{¯}{\alpha }-{\Psi }^{*}\underset{¯}{x}.Math.}^2,$
where J∈R.sup.M×p.sup.2.sup.M is a 1×p.sup.2 block matrix, scaled by the factor

(122) $\frac{1}{p^2},$
with all p.sup.2 blocks being equal to the identity matrix (of dimension M),

(123) 0$\begin{array}{cc}\stackrel{△}{=}\frac{1}{p^2}\underset{\underset{p^2\mathrm{blocks}}{︸}}{\left[\begin{array}{cccc}{I}_M& I_M& \mathrm{.Math.}& I_M\end{array}\right]},& \left(71\right)\end{array}$
and therefore,

(124) $\underset{\_}{\alpha }=\frac{1}{p^2}\underset{i=1}{\overset{p^2}{.Math.}}{\underset{\_}{\alpha }}_i.$

(125) We are now ready to write a practical version of optimization problem (70)

(126) $\begin{array}{cc}\underset{\underset{\_}{\alpha },\underset{\_}{\beta }}{\min }\left({.Math.V\underset{\_}{\alpha }-W\underset{¯}{\beta }.Math.}^2+\mu {.Math.\underset{\_}{\alpha }-{\Psi }^{*}\underset{\_}{x}.Math.}^2\right)s.t.{.Math.\underset{\_}{\alpha }.Math.}^2=c>0,& \left(72\right)\end{array}$
which has the closed-form solution (see the Equation Solution Section below)
{circumflex over (α)}=μ[μJ*J+(1−σ)I−V*WW*V].sup.−1(J*Ψ*x)  (73)
where σ is the smallest singular value of the matrix

(127) $\begin{array}{cc}\left[\begin{array}{cc}\mu {*}^{}+I-V^{*}{\mathrm{WW}}^{*}V& \frac{\mu }{.Math.{\Psi }^{*}\underset{\_}{x}.Math.}{*}^{}{\Psi }^{*}\underset{\_}{x}\\ \frac{\mu }{.Math.{\Psi }^{*}\underset{\_}{x}.Math.}{\left({*}^{}{\Psi }^{*}\underset{\_}{x}\right)}^{*}& \mu \end{array}\right]\in {ℝ}^{\left(p^{2}M+1\right)\times \left(p^{2}M+1\right)},& \left(74\right)\end{array}$
and J*Ψ*x∈R.sup.p.sup.2.sup.M is simply p.sup.2 replicas of Ψ*x∈R.sup.M, scaled by

(128) $\frac{1}{p^2}.$

(129) After estimating the representation of the HR image in terms of the HR basis V, all is left to get an estimate of the HR image is to use equation (16), i.e. û=V{circumflex over (α)}. However, recall (13) and (18) and partition {circumflex over (α)} into p.sup.2 vectors, each of length M, i.e.

(130) $\underset{¯}{\hat{\alpha }}=\left[\begin{array}{c}{\underset{\_}{\hat{\alpha }}}_1\\ {\underset{\_}{\hat{\alpha }}}_2\\ \mathrm{.Math.}\\ {\underset{\_}{\hat{\alpha }}}_{p^2}\end{array}\right].$
Therefore,
û=V{circumflex over (α)}=interlace{Ψ{circumflex over (α)}.sub.i}.sub.i=1.sup.p.sup.2.

(131) Consequently, the construction of neither V nor W is required here. Only the product W*V is needed, which can be very efficiently computed using equations (62)-(68).

(132) It must be re-emphasized, however, that problem formulation (72) is only one possible answer to the following question: How do we find a solution to (37), that is in favor of an optimal approximation of our sought after HR image? In other words, the HR basis matrices V and W only outline a structured subspace framework within which we have complete freedom to seek an estimate of u. Indeed, one can come up with a different objective function, combined with a different set of constraints (or lack thereof), in lieu of (72), but whatever optimization problem one chooses to solve, the matrices V and W will always be in the picture, one way or another.

(133) In that sense, the structured subspace framework (37) can be seen as a conduit through which the pair of ‘LR’ dictionaries (Ψ,Φ) can generate (an estimate of) our ‘HR’ image. Moreover, we shall see next that the same framework provides the foundation for training of said ‘pair’ of LR dictionaries from a ‘single’ LR sequence.

(134) Additionally, it goes without saying that any permutations of equation (37) are exactly equivalent to equation (37). Said differently, since a linear system of equations (e.g. (37)) remains exactly the same system of equations if you simply change the ordering of its equations, then any reformulation of V and W that amounts to permuting their rows does not change the solution space for which equation (37) was set up.

4. LEARNING PAIRS OF LR DICTIONARIES FROM A SINGLE LR SEQUENCE: THE BASELINE

(135) Thus far, we have introduced the (unorthodox) notion that a pair of LR dictionaries can indeed be all we need for estimating a HR image, using the structured subspace model. Moreover, we proposed infusing an element of arbitration to the solution framework, to avoid practical difficulties that might arise with learning a pair of LR dictionaries. In this section we show how the same subspace structure can additionally be used for training pairs of LR dictionaries in the context of estimating HR images from a single LR sequence of images.

(136) Before we proceed, we would like to highlight a notational change. Since we estimate the k-th HR frame, u.sup.k, in a sequence of K frames U=[u.sup.1 u.sup.2 . . . u.sup.K], by estimating patches of it (patch-based processing), we will add the superscript ‘kl.sub.s’ to our symbols to indicate that they are associated with estimating the l.sub.s-th HR patch, u.sup.kl.sup.s∈R.sup.d.sup.2, in the k-th frame u.sup.k, where d=rpq (recall (2)). Specifically, since there is more than one way to partition the same frame into patches, we end up with multiple estimates of the same frame by computing all sets of patches that correspond to all possible partitionings. Let S denote the number of ways a frame can be partitioned, with L patches per partitioning, we use the notation
û.sup.k=P({{û.sup.kl.sup.s}.sub.l.sub.s.sub.=1.sup.L}.sub.s=1.sup.S)  (75)
to indicate that the k-th HR frame is estimated by combining all SL estimated patches {{û.sup.kl.sup.s}.sub.l.sub.s.sub.=1.sup.L}.sub.s=1.sup.S, using some combination process P.

(137) Now, let L denote a learning function (or procedure) that learns a LR dictionary (with orthonormal atoms), using a LR training sample and a training set of LR image patches. In particular, if x.sup.k denotes the k-th LR frame in a sequence of LR frames, X=[x.sup.1 x.sup.2 . . . x.sup.K] that correspond to (2D) decimation of the HR sequence U by a factor of p, i.e. X=D.sub.p.sub.2U, where D.sub.p.sub.2 represents 2D decimation by a factor of p, then the LR patch x.sup.kl.sup.s∈R.sup.r.sup.2.sup.q.sup.2 (corresponding to the HR patch u.sup.kl.sup.s) is used as the training sample, around which a local set of M.sub.tr LR patches X.sup.kl.sup.s∈R.sup.r.sup.2.sup.q.sup.2.sup.×M.sup.tr is extracted from X. Specifically, we use the notation
Ψ.sup.kl.sup.s≡L(X.sup.kl.sup.s,x.sup.kl.sup.s)  (76)
to signify that the LR dictionary Ψ.sup.kl.sup.s∈R.sup.r.sup.2.sup.q.sup.2.sup.×M (expected to approximate the PPPCs of u.sup.kl.sup.s) is learned from the training set X.sup.kl.sup.s and the training sample x.sup.kl.sup.s, using the learning function L. Also, the expression
X.sup.kl.sup.s≡K(X,x.sup.kl.sup.s)  (77)
indicates that the training set of patches is extracted from X around x.sup.kl.sup.s using some strategy K.

(138) Similarly, the same learning function L can be used to learn the secondary dictionary Φ.sup.kl.sup.s∈R.sup.r.sup.2.sup.p.sup.2.sup.×N (to approximate the SPPCs of u.sup.kl.sup.s) if there exists a secondary sequence of LR images Y=[y.sup.1 y.sup.2 . . . y.sup.K], such that y.sup.k=D.sub.q.sup.2u.sup.k for k=1,2, . . . , K, where D.sub.q.sup.2 represents 2D decimation by a factor of q. Namely Φ.sup.kl.sup.s≡L(Y.sup.kl.sup.s,y.sup.kl.sup.s), where y.sup.kl.sup.s∈R.sup.r.sup.2.sup.p.sup.2 is the secondary training sample (corresponding to u.sup.kl.sup.s), and Y.sup.kl.sup.s∈R.sup.r.sup.2.sup.p.sup.2.sup.×Ntr is the secondary training set containing N.sub.tr LR patches, i.e. Y.sup.kl.sup.s∈K(Y,y.sup.kl.sup.s).

(139) Following the same notational change for estimates corresponding to patches, the symbol V.sup.kl.sup.s denotes the HR basis matrix constructed from Ψ.sup.kl.sup.s, i.e. V.sup.kl.sup.s=P(⊕.sub.i=1.sup.p.sup.2,Ψ.sup.kl.sup.s). Similarly, W.sup.kl.sup.s=Q(⊕.sub.i=1.sup.q.sup.2Φ.sup.kl.sup.s) is the secondary HR basis matrix. Now since we are going to use (72) to solve (37), we only need the construction of the product W*.sup.kl.sup.sV.sup.kl.sup.s, so to simplify presentation we use the following notation
W*.sup.kl.sup.sV.sup.kl.sup.s≡C(Ψ.sup.kl.sup.s,Φ.sup.kl.sup.s),  (78)
to indicate that the matrix product W*.sup.kl.sup.sV.sup.kl.sup.s is constructed from the learned LR dictionary pair (Ψ.sup.kl.sup.s,Φ.sup.kl.sup.s), using equations (62)-(68).

(140) Normally, however, an imaging system would only have one sensor, whose output is X, with p being simply the desired upsampling factor. In other words, a secondary sequence Y does not exist, and, therefore, neither does Y.sup.kl.sup.s, nor y.sup.kl.sup.s. This obstacle can nevertheless be overcome by starting the solution with an initial estimate of the HR sequence. Specifically, let
.Math..sub.(0)≡ε(X)  (79)
denote the initial estimate of the unknown HR sequence, based on the available LR sequence X, using some estimation process ε, then the first estimate of the secondary sequence is given by
Y.sub.(1)=D.sub.q.sub.2.Math..sub.(0),  (80)
and the first version of the secondary LR dictionary can be learned Φ.sub.(1).sup.kl.sup.s≡L(Y.sub.(1).sup.kl.sup.s,y.sub.(1).sup.kl.sup.s), where Y.sub.(1).sup.kl.sup.s and y.sub.(1).sup.kl.sup.s denote the training set of LR patches, and the training sample, respectively, extracted from Y.sub.(1).

(141) With both X and an initial guess .Math..sub.(0) of the HR sequence at hand, our structured subspace model can be used in an iterative joint estimation of the HR sequence ‘and’ training of secondary LR dictionaries, as follows.

(142) For a current estimate û.sup.kl.sup.s of the h HR patch in the k-th frame, first compute (recall (73))
{circumflex over (α)}.sup.kl.sup.s=μ[μJ*J+(1−σ.sub.(t).sup.kl.sup.s)I−V*.sup.kl.sup.sW.sub.(t).sup.kl.sup.sW.sub.(t)*.sup.kl.sup.sV.sup.kl.sup.s].sup.−1(J*Ψ*.sup.kl.sup.sx.sup.kl.sup.s)  (81)
with σ.sub.(t).sup.kl.sup.s being the smallest singular value of (recall (74))

(143) $\begin{array}{cc}\left[\begin{array}{cc}\mu {*}^{}+I-V^{*k{\ell }_s}{W}_{\left(t\right)}^{k{\ell }_s}{W}_{\left(t\right)}^{*k{\ell }_s}{V}^{k{\ell }_s}& \frac{\mu }{.Math.{\Psi }^{*k{\ell }_s}{\underset{\_}{x}}^{k{\ell }_s}.Math.}{*}^{}{\Psi }^{*k{\ell }_s}{\underset{\_}{x}}^{k{\ell }_s}\\ \frac{\mu }{.Math.{\Psi }^{*k{\ell }_s}{\underset{\_}{x}}^{k{\ell }_s}.Math.}{\left({*}^{}{\Psi }^{*k{\ell }_s}{\underset{\_}{x}}^{k{\ell }_s}\right)}^{*}& \mu \end{array}\right],& \left(82\right)\end{array}$
and obtain the current estimate of the HR patch via (recall (13), (16) and (18)),
û.sub.(t).sup.kl.sup.s=V.sup.kl.sup.s{circumflex over (α)}.sub.(t).sup.kl.sup.s=interlace{Ψ.sup.kl.sup.s{circumflex over (α)}.sub.i,(t).sup.kl.sup.s}.sub.i=1.sup.p.sup.2   (83)
where W*.sub.(t)*.sup.kl.sup.sV.sup.kl.sup.s≡C(Ψ.sup.kl.sup.s,Φ.sub.(t).sup.kl.sup.s), and Φ.sub.(t).sup.kl.sup.s≡L(Y.sub.(t).sup.kl.sup.s,y.sub.(t).sup.kl.sup.s), with Y.sub.(t).sup.kl.sup.s and y.sub.(t).sup.kl.sup.s being the training set, and training sample, respectively, extracted from the current secondary LR sequence Y.sub.(t)=D.sub.q.sub.2U.sub.(t−1), where .Math..sub.(t−1) is the previous estimate of the HR sequence.

(144) After computing all SL current estimates of patches, {{û.sub.(t).sup.kl.sup.s}.sub.l.sub.s.sub.=1.sup.L}.sub.s=1.sup.S, corresponding to all S partitionings of the k-th frame, find the current estimate of the k-th HR frame, û.sub.(t).sup.k≡P ({{û.sub.(t).sup.kl.sup.s}.sub.l.sub.s.sub.=1.sup.L}.sub.s=1.sup.S), for k=1,2, . . . , K, to get a current estimate .Math..sub.(t)=[û.sub.(t).sup.1 û.sub.(t).sup.2 . . . û.sub.(t).sup.K] of) the HR sequence.

(145) TABLE-US-00002 TABLE 2 Pseudo code illustrating the basic solution baseline based on the proposed structured subspace framework. Input: A sequence     of K LR images, and the desired upsampling factor, p. Compute the initial estimate    .sub.(0) ≡ ε(   ) of the HR sequence. For t = 1:T  Obtain   .sub.(t) =    .sub.q.sub.2    .sub.(t−1).  For k = 1:K   For s = 1:S    For   .sub.s = 1:L         Learn the primary LR dictionary:     Extract       Extract     ≡    (   ).     Compute     ≡       (Note:     are iteration-independent).         Learn the secondary LR dictionary:     Extract       Extract     ≡         Compute     ≡             Estimate the HR image within the structured      subspace framework (37). We recommend (72)     Compute     ≡         Compute       Compute smallest singular value     of matrix (82).     Compute     (81).     Compute     (83).    End   End   Compute    End  Construct   End Output:     = [û.sub.(T).sup.1 û.sub.(T).sup.2 ... û.sub.(T).sup.K].

(146) Repeat until a prescribed number of iterations, T, has been reached to get the final estimate .Math. of the HR sequence,
.Math.≙.Math..sub.(T)=[û.sub.(t).sup.1û.sub.(t).sup.2 . . . û.sub.(t).sup.K].  (84)

(147) Refer to Table 2 and FIG. 3 for a summary of the baseline.

(148) At this point, it is worthwhile to note that, all KSL estimates, {û.sup.kl.sup.s}, 1≤k≤K, 1≤s≤S, and 1←l.sub.s≤L, needed to create the HR sequence .Math..sub.(t), are computed completely independently from each other, and hence they can all be concurrently computed. Namely, if the computing hardware allows it, it is possible to parallel-process all KSL patches at once (per iteration).

(149) Moreover, it should be noted that the proposed baseline is the most straightforward, and least computationally expensive for implementing the structured subspace framework. However, more complex baselines can still benefit from the structured subspace framework by jointly estimating HR patches, for example. In other words, the solution baseline can be devised such that the estimation of any patch is made dependent on the estimation of other patches (e.g. neighboring patches within the same frame, or patches across multiple estimates of the same frame, or even patches across different frames). However, we do not believe that the improvement in results, if any, would justify the added complexity. Specifically, the real power of our solution lies in the unprecedented access it provides to very narrow training sets (that are extracted from a naturally highly correlated sequence of LR images).

5. A WORKING SOLUTION

(150) In the previous section we described a general baseline for iterative estimation of a HR sequence from a LR sequence, by learning pairs of LR dictionaries within a structured subspace framework. However, many details were (intentionally) left out. Specifically, which estimation process ε are we going to use to obtain an initial guess of the HR sequence? How about the learning function L? Regarding those local training sets, what are the specifics of their extraction 3C from a LR sequence? Which combination process P is going to be used to piece together all SL estimated HR patches per frame ? Indeed, what are the S different possible ways to partition a frame into patches?

(151) The discussion of such details was postponed so as to help appreciate the cornerstone role of the structured subspace model, and to emphasize the level of freedom in the design of the baseline's remaining components. Indeed, as we shall see in the experiments section, we can get impressive results despite using some of the simplest outlines for ε, K, L, and P. This, again, underscores the baseline's most valuable asset: the structured subspace model which, in effect, leverages a sequence of LR images as a very narrow training set for estimating HR images.

(152) A. The Initial Guess (ε)

(153) For an initial estimate (79) of the HR sequence from the available LR sequence, we choose ε to simply represent Lanczos interpolation (by a factor of p). One might try more advanced options, but besides added complexity, even advanced methods would not give appreciably better estimates compared to simple image interpolation methods (such as bicubic or Lanczos) when the LR sequence contains complex motion patterns, and the aliasing is relatively strong.

(154) B. Partitioning a Frame

(155) In its simplest form, patch-based processing works by dividing a frame we desire to estimate into patches (sub-regions), and then applying the proposed estimation process for each individual patch, independently from other patches. By simply putting each estimated patch in its place, we get an estimate of the HR frame (i.e. an estimate of a HR frame is obtained by ‘tiling’ its estimated patches).

(156) However, multiple estimates of the same frame can be obtained by dividing the frame into overlapping patches. To elaborate, since the size of a HR patch is d×d and d=rpq, if we use an overlap of pq pixels, vertically and/or horizontally, then a HR frame can be partitioned into patches in S=r.sup.2 different ways. In particular, since the location of the ‘leading’ patch u.sup.kl.sup.s.sup.=1 of the s-th partitioning determines the locations of the remaining patches (since patches are tiled next to each other, vertically and horizontally, per partitioning), we shall give a description of the location of the leading HR patch corresponding to partitioning s, as follows [2].

(157) $\begin{array}{cc}{u}^{k{\ell }_s=1}=u^k\left(v\left(s\right)pq+1\text{:}v\left(s\right)pq+d,h\left(s\right)pq+1\text{:}h\left(s\right)pq+d\right),\mathrm{where}& \left(85\right)\\ v\left(s\right)=.Math.\frac{s}{r}.Math.-1,h\left(s\right)=\mathrm{mod}\left(s-1,r\right),& \left(86\right)\end{array}$
and A(r.sub.1:r.sub.2, c.sub.1:c.sub.2) denotes the submatrix of A, from row r.sub.1 to row r.sub.2, and column c.sub.1 to column c.sub.2.

(158) Therefore, the corresponding leading LR patch x.sup.kl.sup.s.sup.=1 of the k-th LR frame x.sup.k is
x.sup.kl.sup.s.sup.=1=x.sup.k(v(s)q+1:v(s)q+rq,h(s)q+1:h(s)q+rq).  (87)

(159) Similarly, the corresponding leading LR patch y.sub.(t).sup.kl.sup.s.sup.=1 of the k-th (secondary) LR frame y.sub.(t).sup.k is s
y.sub.(t).sup.kl.sup.s.sup.=1=y.sub.(t).sup.k(v(s)p+1:v(s)p+rp,h(s)p+1:h(s)+rp).  (88)
C. Extracting Local Training Sets (K)

(160) We follow the simple procedure K described in [2], wherein (overlapping) LR patches are extracted from within the spatiotemporal neighborhood of the patch x.sup.kl.sup.s, as follows. Pick the κ LR frames in X that are (temporally) closest to the k-th frame. Pick the (spatial) location of the patch x.sup.kl.sup.s as the origin. Extract, from each one of the κ frames, (2b+1).sup.2 spatially overlapping patches, starting at the origin, and moving (vertically and horizontally, and in opposite directions) with a step of only one pixel, and for as far as b pixels from the origin.

(161) The result is a (local) training set X.sup.kl.sup.s containing M.sub.tr=κ(2b+1).sup.2 LR patches that are extracted from within a spatiotemporal neighborhood around x.sup.kl.sup.s, with temporal width of κ frames, and vertical (horizontal) spatial width of (2b+1) pixels. The following pseudocode describes K more explicitly.

(162) Suppose x.sup.kl.sup.s=x.sup.k(r.sub.0+1:r.sub.0+rq, c.sub.0+1:c.sub.0+rq), then a local training set can be extracted from X, around x.sup.kl.sup.s, as follows.

(163) TABLE-US-00003 (89) For i.sub.k = −0.5κ + k : 0.5κ + k − 1  For i.sub.r = −b:b   For i.sub.c = −b:b    A = x.sup.i.sup.k (r.sub.0 + i.sub.r + 1 : r.sub.0 + i.sub.r + rq, c.sub.0 + i.sub.c + 1 : c.sub.0 + i.sub.c + rq).    Put vec(A) as a column in       End  End End

(164) Similarly, Y.sub.(t).sup.kl.sup.s is extracted, around y.sub.(t).sup.kl.sup.s, from Y.sub.(t), using the same exact procedure (with N.sub.tr=M.sub.tr=κ(2b+1).sup.2).

(165) D. Learning the Dictionaries (L)

(166) Dictionary learning (L) can be based on feature selection (FS), feature extraction (FE), or a combination of both. In our case, FS far outweighs FE. In particular, PPCs are highly correlated signals, so learning a dictionary that represents them well would be heavily dependent on choosing those features (LR patches) in the training set that would capture the subtle differences between the PPCs [2].

(167) To elaborate, since the PPCs are expected to be highly correlated with x.sup.kl.sup.s, we use it as a training sample for selecting the M«M.sub.tr‘best’ LR patches in X.sup.kl.sup.s. To that end, we use best individual feature (BIF) selection based on L1 distance between

(168) $\frac{{\underset{¯}{x}}^{k{\ell }_s}}{.Math.{\underset{¯}{x}}^{k{\ell }_s}.Math.}$
and normalized LR patches in X.sup.kl.sup.s. In other words, we select a subset X.sub.m.sup.kl.sup.s with only M LR patches out of the M.sub.tr=κ (2b+1).sup.2 LR patches in X.sup.kl.sup.s, where the elements in the selected subset are simply the normalized LR patches in the training set that have the smallest L1 distance from the normalized training sample.

(169) Now, since some of our equations are based on the (simplifying) assumption that the dictionary atoms are orthonormal, all is left to do after selecting the best M LR patches is to simply orthonormalize them. This can be done, for example, using the Gram Schmidt method (QR factorization). However, we shall be using the singular value decomposition (SVD), since it keeps the door open for further dimensionality reduction if so desired. In other words, our local dictionary Ψ.sup.kl.sup.s is simply the first M left singular vectors of the selected subset X.sub.M.sup.kl.sup.s. If {tilde over (M)}<M is desired, then we simply keep the first {tilde over (M)} left singular vectors instead.

(170) Likewise, Φ.sub.(t).sup.kl.sup.s is the first N left singular vectors of the N-element subset Y.sub.N,(t).sup.kl.sup.s selected from the (normalized) training set Y.sub.(t).sup.kl.sup.s, using the BIF method based on the smallest L1 distance from

(171) $\frac{{\underset{\_}{y}}_{\left(t\right)}^{k{\ell }_s}}{.Math.{\underset{¯}{y}}_{\left(t\right)}^{k{\ell }_s}.Math.}.$

(172) Clearly, the dictionary learning process described here is exceedingly simple, yet it works surprisingly well. It must be noted, however, that this simplicity of dictionary learning is afforded to us by the ‘narrowness’ of the training sets, which, in turn, is attributed to the ‘high correlation’ between images in a LR sequence and the PPCs of a HR frame (FIG. 1). Yet, without the structured subspace framework as a tool to harness the power of narrow training LR sets for estimating HR frames, we would have been confined to the imaging hardware requirement of [1], [2].

(173) E. Combining Multiple Estimates of a HR Frame (P)

(174) Given Ψ.sup.kl.sup.s and Φ.sub.(t).sup.kl.sup.s (and x.sup.kl.sup.s) for 1≤l.sub.s≤L, we can estimate û.sub.(t).sup.kl.sup.s for 1≤l.sub.s≤L, and tile all estimated patches into the s-th estimate û.sub.s,(t).sup.k of the k-th HR frame. Repeating this for all s (1≤s≤S=r.sup.2), we get S=r.sup.2 estimates {û.sub.s,(t).sup.k}.sub.s=1.sup.r.sup.2 of the k-th HR frame.

(175) So now, what do we do with all these estimates of the same frame ? One of the simplest options is to average all r.sup.2 estimates. We, nevertheless, choose to follow the solution proposed in [2] to get the final (per iteration t) estimate û.sub.(t).sup.k of the k-th HR frame, as follows.

(176) 0$\begin{array}{cc}{\underset{¯}{\hat{u}}}_{\left(t\right)}^k=\underset{{\underset{¯}{u}}_{\left(t\right)}^k}{\min }\underset{s=1}{\overset{r^2}{.Math.}}\underset{w_1}{.Math.}\underset{w_2}{.Math.}{.Math.H{\underset{¯}{u}}_{\left(t\right)}^k-R_{w_1,w_2}{\underset{\_}{\hat{u}}}_{s,\left(t\right)}^k.Math.}_1,& \left(90\right)\end{array}$

(177) where R.sub.w.sub.1.sub.w.sub.2 shifts the s-th estimate by w.sub.1 and w.sub.2 pixels in the horizontal and vertical direction, respectively, −(ω−1)/2≤w.sub.1,w.sub.2≤(ω−1)/2, and H represents a blurring kernel, to be defined by the user, to counteract the softening effect of median filtering (across all r.sup.2 estimates, and within a spatial window of size ω×ω, per pixel).

(178) F. Solution Parameters

(179) Now we turn to the solution parameters, and start with the integer scalar r. Given that, in practice, p is the user-defined upsampling factor, and given that q=p+1 (1), then r is the only parameter left to determine the size of the HR patch (2) or, in particular, r determines how small (41) the HR patch can be. Before we continue, now is the time to address the question: why would we want to work on patches of a frame (patch-based processing), instead of upsampling an entire frame as a whole (i.e. with r being large enough to engulf the entire frame)?

(180) The answer to this question is in [2], but we reproduce it here for convenience. First, the short version: the more ‘dynamic’ the LR sequence, the greater the need for patch-based processing, with particularly small patches. The longer version of the answer would require revisiting FIG. 1 to notice that the inter-correlations between the PPCs are much greater than the correlations that exist between the LR frames (of a dynamic scene) and the PPCs (of a frame). In other words, a ‘training’ set of LR frames would not capture well the subtle differences between the PPCs of a HR frame as a whole. Note that in addition to a basic common component the PPCs might share, these subtle differences are, ultimately, what gives us the HR frame. Working on smaller sub-regions (patches), however, allows us to harness, from the LR frames, local training sets (as described above) from which local dictionaries can be learned for estimating HR patches. In other words, the smaller the patch, the more ‘local’ the training sets (and the dictionaries) get.

(181) On the other hand, larger patches (larger r) allow for the advantage of a more overdetermined (38). So to keep things simple, when the desired upscaling factor p=3, we fix the dimensionality of all LR dictionaries at 3q.sup.2, i.e. M=N=3q.sup.2, and let r=4. For p=4, we use r=3 and M=N=2q.sup.2.

(182) Now, how large a local training set can be (recall N.sub.tr=M.sub.tr=κ(2b+1).sup.2)? Ideally, we would want to search for ‘best’ LR patches in as large a spatiotemporal neighborhood as possible, i.e. κ and b both can be as large as possible. However, and even though using BIF selection based on L1 distance from a training sample is a very low complexity FS policy, with more and more LR patches to choose from, the computational cost of total FS (performed for all patches) can get quite significant if K and b are too large. For that reason, we pick κ=30 and b=5.

(183) Finally, for combining multiple estimates of the same frame (90), we choose H to be a Gaussian kernel with size 5 and standard deviation of 1, and we pick ω=3.

(184) G. The Control Parameter μ

(185) Recall that equation (81) and (82) are simply the patch-based, iteration-dependent versions of equations (73) and (74), respectively. Yet, the control parameter μ in these equations appears patch-independent as well as iteration-independent. Can this indeed be the best strategy for selecting μ? To use a fixed value across all patches and throughout all iterations?

(186) First, recall that μ determines how close to the reference LR image should the mean of the PPCs of the estimated image be. Therefore, if μ is too large, we run the risk of an estimated HR image that has PPCs that are too close to each other (a solution with some residual aliasing). On the other hand, in practice, if μ is not high enough, the end result can be too smooth, with low contrast, especially for those parts of the scene that are too dynamic (large scene changes across captured frames). In other words, instead of using the symbol μ in equations (81) and (82), we could use μ.sub.(t).sup.kl.sup.s to account for the patch and iteration dependency in selecting the control parameter.

(187) Nevertheless, choosing the “right” value for μ is not easy. According to experiments, the range of good values can be quite large (0.1≤μ≤5), and a good choice can be dependent on many factors including high spatial frequency content per patch, patch-size, and selected dictionary atom diversity. Currently, we use the following empirical formula for determining good values for μ (per-patch and per-iteration).

(188) $\begin{array}{cc}{\mu }_{\left(t\right)}^{k{\ell }_s}=\min \left(\max \left({\mu }_0\frac{d^2}{p^{2}M+q^{2}N}\frac{\gamma }{\tau },0.1\right),5\right),\mathrm{where}& \left(91\right)\\ \tau =\frac{{.Math.y_{\left(t\right)}^{k{\ell }_s}**\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right].Math.}_F^2}{r^2{p}^2},& \left(92\right)\end{array}$
** denotes 2D convolution, ∥.Math.∥.sub.F is the Frobenius norm, μ.sub.0 is a user-defined scalar that is fixed across all iterations and all patches (we pick μ.sub.0=20), and γ is a measure of the diversity of selected atoms computed as follows. 1—Find the correlation coefficient between x.sup.kl.sup.s and each atom in X.sub.M.sup.kl.sup.s. 2—Assign a value of 1 for correlation coefficients higher than 0.99, a value of 3 for correlation coefficients between 0.95 and 0.99, and a value of 10 for correlation factors below 0.95. 3—γ is the average of all assigned values.

(189) No matter how good an empirical formula for selecting μ.sub.(t).sup.kl.sup.s might be, if the value for μ.sub.0 is too high for some frames (or sub-regions of a frame), we would end up with estimated HR frames that contain regions with some residual aliasing. If that is the case, we propose the following remedy.

(190) After obtaining the final estimate .Math.≙.Math..sub.(T)=[û.sub.(T).sup.1 û.sub.(T).sup.2 . . . û.sub.(T).sup.K] of the HR sequence (we pick T=4), we start the solution over, but this time we use .Math. as our input sequence of LR images, and we keep all the parameters the same except for μ and r. Specifically, we keep μ fixed at the low value of 0.2, and use r=7 for p=3, and r=5 for p=4. This should take care of residual aliasing. The final sequence can then be resized down to the original desired frame size (i.e. each frame can be decimated by a factor of p).

6. RESULTS

(191) In this section we test the proposed solution using the ‘Suzie’ sequence of HR frames of size 486×720 each. First, we decimate each frame by a factor of 6 (vertically and horizontally). The top row of FIG. 4 shows some of the LR frames magnified by a factor of 3 using Lanczos interpolation. Corresponding upsampling results of the proposed solution using an upsampling factor of p=3 are shown in the bottom row.

(192) We repeat the same experiment using a decimation factor of 8, and an upsampling factor of p=4. FIG. 5 shows Lanczos-magnified LR images vs. upsampling.

(193) If we compare both FIG. 4 and FIG. 5, it becomes immediately evident that the upsampling quality is higher for p=3 compared to p=4. Of course, this can be justified by the fact that the LR sequence of FIG. 4 is less degraded than FIG. 5 LR sequence. But why, specifically in our case, is it easier to get better results for lower upsampling factors? It all goes back to (local) LR dictionaries' ability to capture the subtle differences between PPCs. Specifically, the higher the upsampling factor needs to be, the higher the number of PPCs, and the more difficult it becomes to represent the subtle details differentiating between PPCs, especially when the scene is highly dynamic. Moreover, with larger upsampling factors, the patch size has to be larger (recall (2) and (41)), and the larger the patches the more difficult it becomes to find good example images for learning local dictionaries.

(194) However, this gap in upsampling quality (e.g. p=4 vs. p=3) would become smaller when the LR sequence is not highly dynamic (specifically, when scene changes, across a set of captured frames, are not too large because the frame rate is sufficiently high).

7. CONCLUSION

(195) The basic principle of learning-based image upsampling is fairly simple. In general, learning-based methods use the available LR image as ‘partial measurements’ (of the unknown HR image) to estimate the representation of the HR image in terms of a HR dictionary learned from a training set of HR example images. While these methods certainly differ in the details, the two most important features that account for most of the variation in results are the training sets they use (generic vs. narrow—special—sets) and the learning process of the HR dictionaries.

(196) In contrast, the work in [1], [2], uses the same learning-based upsampling principle, although for different target signals: the PPCs of the HR image, which are LR signals themselves. This shift of focus allowed the authors of [1], [2] to exploit a sequence of LR images as a source of very narrow training sets, from which highly efficient LR dictionaries can be easily learned to represent the PPCs (of sought after HR images). The major obstacle to their idea was the lack of ‘partial measurements’ required to find the representations of the PPCs in terms of the LR dictionaries. To overcome this obstacle, they proposed a hardware modification such that the imaging device is equipped with another (different resolution) LR sensor (and a beam splitter) to provide the much needed partial measurements for their target images (the PPCs).

(197) Compared to [1], [2], in current work, we shift the focus back to estimating the representation of the HR image itself, and in terms of a pair of HR bases, each of which is embedding a (different resolution) LR dictionary. The end result is a structured subspace framework through which pairs of LR dictionaries are used to estimate HR images directly, thus circumventing the very restrictive hardware assumption required in [1],[2], while still benefiting from the huge advantage of harvesting narrow training sets from a LR sequence of images.

Applied Example of Working Solution

(198) In this section, we give a visualization of the basic solution baseline whose details were laid out above in “A Working Solution”. To help with understanding, it is important to keep the big picture in mind, and which we summarize as follows:

(199) We estimate the entire HR sequence by estimating patches (pieces) of each frame.

(200) Because we estimate each HR patch based on a pair of (local, different resolution) LR dictionaries, we need a pair of LR sequences (of different resolution) from which to extract the pair of local dictionaries (the algorithm takes in a pair of LR dictionaries, per patch, as input, and spits out an estimate of the HR patch).

(201) The reason we need to iterate the solution is because we do NOT have the secondary (different resolution) LR sequence (we keep estimating a better version of the missing secondary LR sequence).

(202) The Initial Guess (ε)

(203) Let's first have a look at sample LR frames of an available LR sequence (we call it the primary LR sequence) in FIG. 6.

(204) To get the first estimate of the secondary (missing) sequence Y.sub.(1), we need to first get an initial guess of the HR sequence .Math..sub.(0). This can be done, for example, by using 2D interpolation (by a factor of p, we choose p=3 here) of each frame x.sup.k of the primary (available) LR sequence X. Using Matlab, this can be done as follows:

(205) TABLE-US-00004 (93) for k = 1:K  û.sub.(0).sup.k = imresize (x.sup.k,p,‘lanczos3’); end

(206) Now, let's take a look at the sample frames of the initial estimate of the HR sequence, shown in FIG. 7. Comparing FIG. 7 to FIG. 6, it becomes clear that the initial guess of the HR frames don't look appreciably better than the available (primary) LR frames, but this shall suffice for an initial guess. Note: each frame shown in FIG. 7 is actually p.sup.2 times bigger than the corresponding frame shown in FIG. 6.

First Iteration (t=1)

(207) To obtain the first estimate of the secondary (different/lower resolution) LR sequence, we simply decimate each frame of the HR sequence by a factor of q (q=p+1=4 here). Using Matlab, this can be done using the following code:

(208) TABLE-US-00005 (94) for k = 1:K  y.sub.(1).sup.k = û.sub.(0).sup.k (1:q:end,1:q:end); end

(209) Because of the size of the figures, comparing FIG. 8 to FIG. 6, it might not be easy to see that the first estimate of the secondary LR sequence is lower resolution than the primary LR sequence, so in FIG. 9, we compare bigger versions of one of the frames (frame #31) from the primary and secondary sequences. It might still not be easy to see, but careful examination reveals the resolution is different.

(210) Partitioning a Frame

(211) So now that we have two (different resolution) LR sequences as input, X and Y.sub.(1), how do we proceed? Like (almost) everyone else, we do not estimate an entire HR frame as a whole, instead we estimate patches (pieces) of each frame, so the next question is: how do we partition a frame into patches to be estimated?

(212) In Section 5 (A Working Solution), we provided a description of the simplest partitioning scheme (equations (85)-(88)). For a pictorial illustration, we pick frame #31 of the primary sequence X and display, in FIG. 10, the S=r.sup.2 possible ways of partitioning the same frame (we pick r=4 here, so each primary LR patch x.sup.kl.sup.s has size rq×rq=16×16).

(213) The so-called ‘leading’ patch per partitioning pattern s is shown as a brightened square in FIG. 10. We use the description ‘leading’ patch because it simply tells us how to cut up a frame (its location defines the partitioning pattern). So a frame is estimated in patches (pieces), and there are multiple ways to cut up a frame into patches.

(214) Now, what happens when we compute all the HR patches corresponding to all the squares (patches) shown in FIG. 10? We get S=r.sup.2=16 estimates of the same HR frame.

(215) But why don't we simply use just one partitioning pattern? There's nothing wrong with that except that getting more estimates of the same frame (by partitioning it using S=r.sup.2 patterns instead of just one pattern), allows us to reduce estimation errors since different patterns would give different estimation errors and therefore having multiple estimates of the same frame allows us to reduce estimation errors. There are many different ways one could use multiple estimates of the same frame to reduce estimation errors (equation (90), for example, is only one of them), and most people simply take the mean—or the median—of all estimates.

(216) To illustrate the benefits of having multiple estimates of the same frame, we show in FIG. 11 the S=r.sup.2=16 HR estimates of frame #31 after tiling (in place) all estimated HR patches {û.sub.(1).sup.kl.sup.s}, corresponding to all LR patches {x.sup.kl.sup.s}, shown in FIG. 10. Some of the most obvious estimation errors are present in some of the ‘edge’ patches, but beyond these obvious errors, subtle errors are present (and different) per the same image part. For example, if we concentrate on the mouth region, we′d see that it looks a little different across all 16 estimates.

(217) FIG. 12 shows what we get if we compute the median of all 16 estimates of frame #31 shown in FIG. 11 vs. the initial estimate of the same frame (obtained using Lanczos interpolation).

(218) But how did we compute each of the HR patches shown in FIG. 11? For example, how did we compute the HR patch û.sub.(1).sup.kl.sup.s at k=31, s=4, and l.sub.s=20? Let's highlight the patch in question in FIG. 13.

(219) Extracting Local Training Sets (K)

(220) Recall that, in our solution, estimating a HR patch requires a local pair of (different resolution) dictionaries to be extracted from a pair of LR sequences (of different resolution). We already had LR sequence X (FIG. 6) and from which we created Y.sub.(i) (FIG. 8). So now, we visualize how we extract the training sets X.sup.kl.sup.s and Y.sub.(1).sup.kl.sup.s from the LR sequences X and Y.sub.(1), respectively, in preparation for computing the highlighted patch shown in FIG. 13.

(221) First of all, recall parameters K (the width of the temporal neighborhood of frames around frame k from which training LR patches are to be extracted) and b (2b+1 is the width and height, in pixels, of the spatial window whose center is defined by the location of the LR patch corresponding to the HR patch we need to estimate). Here, we choose K=30 frames, and b=5 pixels. FIG. 14 illustrates the extraction process (X) of example LR patches to be used as elements of the pair of training sets.

(222) Learning the Dictionaries

(223) With (2b+1).sup.2=121 LR patches extracted from each frame in a predefined temporal neighborhood of size K=30 frames, the size of the entire training set X.sup.kl.sup.s is 30×121=3,630 LR patches (similarly, training set Y.sub.(1).sup.kl.sup.s has 3,630 LR patches).

(224) Clearly, however, we do not need all patches in a training set. Indeed, just by looking at FIG. 14, it becomes evident that most of the extracted patches are irrelevant, anyway. But then why do we extract so many LR patches? The point is to find LR patches that are similar to the LR patch we want to upsample, so the extraction process described here is our way for ‘fishing’ for good example LR patches. And although most of the ‘catch’ is useless, it doesn't matter as long as we have enough ‘good’ LR patches in our training set.

(225) Before we proceed, here we use M=N=2q.sup.2=32 as the number of (hopefully) ‘good’ enough LR patches to be chosen from each training set. All is left to do is define a measure of ‘goodness’. The following steps summarize the dictionary “learning” process we described in Section 5 (A Working Solution).

(226) First, execute the following pseudocode:

(227) TABLE-US-00006 For i = 1: κ(2b + 1).sup.2 $d_x\left[i\right]={.Math.\frac{X^{k{\ell }_s}\left[i\right]}{.Math.X^{k{\ell }_s}\left[i\right].Math.}-\frac{x^{k{\ell }_s}}{.Math.x^{k{\ell }_s}.Math.}.Math.}_1$ (95) $d_y\left[i\right]={.Math.\frac{Y_{\left(1\right)}^{k{\ell }_s}\left[i\right]}{.Math.Y_{\left(1\right)}^{k{\ell }_s}\left[i\right].Math.}-\frac{y_{\left(1\right)}^{k{\ell }_s}}{.Math.y_{\left(1\right)}^{k{\ell }_s}.Math.}.Math.}_1$ End

(228) where the LR patches x.sup.kl.sup.s and y.sub.(1).sup.kl.sup.s correspond to the HR patch you want to estimate (these are the highlighted areas in the left and right LR frames, respectively, shown in the middle row of FIG. 14).

(229) The i-th entry in d.sub.x is thus the value of the L1 distance between the i-th normalized LR patch in training set X.sup.kl.sup.s and normalized (reference) LR patch x.sub.(1).sup.kl.sup.s. Similarly, the i-th entry in d.sub.y is the value of the L1 distance between the i-th normalized LR patch in training set Y.sub.(1).sup.kl.sup.s and normalized LR patch y.sub.(1).sup.kl.sup.s.

(230) Second, pick the M LR patches, in the training set X.sup.kl.sup.s, that correspond to the smallest M entries in d.sub.x. These are the atoms of your local primary dictionary Ψ.sup.kl.sup.s. Similarly, the atoms of the local secondary dictionary Φ.sub.(1).sup.kl.sup.s are the N LR patches, in the training set Y.sub.(1).sup.kl.sup.s, that correspond to the smallest N entries in d.sub.y.

(231) Third, since it's easier to work with orthonormal dictionaries, use SVD decomposition to orthonormalize the atoms of each dictionary. After converting atoms, of each dictionary, from their 2D form to their vector form, and stacking these (per dictionary) as column-vectors next to each other, we can use the following Matlab code to orthonormalize our dictionaries:
[Ψ.sup.kl.sup.s,˜,˜]=svd(Ψ.sup.kl.sup.s,0);
[Φ.sub.(1).sup.kl.sup.s˜,˜]=svd(Φ.sub.(1).sup.kl.sup.s,0);  (96)

(232) FIG. 15 visualizes the pair of dictionaries we used to compute the highlighted patch in FIG. 13.

(233) Computing the HR Patch

(234) All the previous steps were taken to set the stage for the proposed structured subspace framework to do its work: converting a pair of (different resolution) LR dictionaries into a HR image. In the following, we describe the most straightforward application of the structured subspace solution.

(235) First, construct the HR basis matrix V.sup.kl.sup.s from the LR dictionary Ψ.sup.kl.sup.s, and the HR basis matrix W.sub.(1).sup.kl.sup.s from the LR dictionary Φ.sub.(1).sup.kl.sup.s, by running the following Matlab code:
V.sup.kl.sup.s=LR2HRbasis(Ψ.sup.kl.sup.s,r,p,q,M);
Ψ.sub.(1).sup.kl.sup.s=LR2HRbasis(Φ.sub.(1).sup.kl.sup.s,r,q,p,N);  (97)
Recall that we chose p=3 (and thus q=p+1=4), r=4, and M=N=2q.sup.2=32. LR2HRbasis is a predefined function written using Matlab as follows:

(236) TABLE-US-00007 (98) % function to construct HR basis from LR dictionary function V = LR2HRbasis(PSI,r,p,q,M) d = r*p*q; PSI = reshape(PSI,r*q,r*q,M); % this line converts % the vectorized % atoms (columns) % of PSI back to % their 2D form V = zeros(d*d,p*p*M); zrs = zeros(d,d,M); c = 0; for i = 1:p  for j = 1:p   c = c+1;   atomZshift = zrs;   atomZshift(i:p:end,j:p:end,:) = PSI;   V(:,(c−1)*M+1:c*M) = reshape(atomZshift,d*d,M);  end end

(237) Second, solve the following equation whichever meaningful way possible: V.sup.kl.sup.sα.sub.(1).sup.kl.sup.s=W.sub.(1).sup.kl.sup.sβ.sub.(1).sup.kl.sup.s. Recall: V.sup.kl.sup.s is the HR basis matrix that spans any HR image (patch) whose PPCs are spanned by the LR dictionary Ψ.sup.kl.sup.s. W.sub.(1).sup.kl.sup.s is the HR basis matrix that spans any HR image (patch) whose PPCs are spanned by another (different resolution) LR dictionary Φ.sup.kl.sup.s. α.sub.(1).sup.kl.sup.s is the representation of (first estimate of) the HR image (patch) in terms of V.sup.kl.sup.s. β.sub.(1).sup.kl.sup.s is the representation of the same HR image in terms of W.sub.(1).sup.kl.sup.s. Solving this equation is the equivalent of searching for the pair of HR images (one image in span(V.sup.kl.sup.s), and the other image in span(W.sub.(1).sup.kl.sup.s)) that are ‘closest’ to each other. So what remains is to define the “distance measure”, or the ‘yardstick’ with which distance is measured. We believe a good yardstick for measuring distance, in this context, is defined by optimization problem (72). Using this yardstick (72), it just so happens that we do not need to construct V nor W, only the multiplication W*V is needed which, due to the particular structure of V and W, can be very efficiently computed using equations (62)-(68). However, generally speaking, other yardsticks one might choose to use (other than (72)) might require the explicit construction of V and W. For this reason, we shall ignore, in this demonstration, the special computational efficiency advantage imparted by (72). Since we're working on patches, x.sup.kl.sup.s is the reference LR image here (the highlighted square in the LR frame shown on the left of the middle row of FIG. 14) whose HR version (the highlighted square in FIG. 13) we seek to compute.

(238) The following is the Matlab code we used for computing the HR image patch highlighted in FIG. 13.
û.sub.(1).sup.kl.sup.s=computeHRpatch(V.sup.kl.sup.s,W.sub.(1).sup.kl.sup.s,Ψ.sup.kl.sup.s,x.sup.kl.sup.s,μ.sub.(1).sup.kl.sup.s,r,p,q,M);  (99)
where the control parameter μ.sub.(1).sup.kl.sup.s determined according to (91), and computeHRpatch is a predefined function as follows:

(239) TABLE-US-00008 % function to solve Vα = Wβ, according to % optimization problem (72). function u = computeHRpatch(V,W,PSI,x,mu,p,M) A = W′*V; A = A′*A; x = x(:); % get the vector form of % the reference LR image. b = PSI′*x; % i.e. b = Ψ*x. brep = repmat(mu/p/p/norm(b)*b,p*p,1);                $\%\mathrm{this}\mathrm{is}\frac{\mu }{.Math.{\Psi }^{*}\underset{\_}{x}.Math.}{𝒥}^{*}{\Psi }^{*}\underset{\_}{x}.$ J = 1/p/p*repmat(eye(M),1,p*p); % this is eqn (71). (100) J = J′*J; E = zeros(p*p*M+1); % construct the matrix in (74). E(1:p*p*M,1:p*p*M) = mu*J +eye(p*p*M) −A; E(1:p*p*M,p*p*M+1) = brep; E(p*p*M+1,1:p*p*M) = brep′; E(p*p*M+1,p*p*M+1) = mu; sigma = svd(E); sigma = sigma(end); A = mu*J + (1−sigma)*eye(p*p*M) −A; brep = repmat(mu/p/p*b,p*p,1); alpha = A \ brep; % this is eqn (73). u = V * alpha; u = reshape(u,r*p*q,r*p*q); % get the 2D form of % the estimated HR image.
Parallel Processing:
Thus far, we explained how we got the HR patch highlighted in FIG. 13, which is patch û.sub.(1).sup.kl.sup.s at k=31, s=4, and l.sub.s=20, but that's just one patch among 48 patches (l.sub.s=1, . . . , 48) per partitioning pattern s=4, and we have S=16 partitioning patterns with 48 patches per partitioning (FIG. 11). In our simple baseline, all of these HR patches are computed independently from each other, so if the computing hardware allows it, they can all be computed simultaneously, and for all K frames too. So if we have K=100 frames, we can compute 100*16*48=76,800 patches simultaneously, and then arrange them as 16 HR estimates per frame, and then compute the median of all 16 estimates per frame to get the first estimate of the our HR sequence .Math..sub.(1). Of course, this can also be done in series, computing one HR patches at a time. In any case, parallel computing is beyond the scope of this invention.

2nd and Subsequent Iterations

(240) Compared to the first iteration, the second iteration (and subsequent iterations) are no different except for working with an updated estimate of the secondary LR sequence, which is obtained by simply decimating the previous estimate of the HR sequence (e.g. we get the updated sequence Y.sub.(2) by decimating each frame in .Math..sub.(1) by a factor of q).

(241) In particular, the primary LR sequence X, remains the same for all iterations, so unless we adopt a different training set extraction strategy (X) or a different dictionary learning procedure (L) per iteration, there is no need to repeat the training/learning process of local primary LR dictionaries. For example, if we stored Ψ.sup.kl.sup.s shown in FIG. 15, we do not need to repeat its extraction process for any iteration beyond the first one. This is why the symbol Ψ.sup.kl.sup.s is missing the subscript (t).

Equation Solution Section

(242) To solve the optimization problem (72) using Lagrange multipliers, we first rewrite (72) as

(243) $\begin{array}{cc}\underset{\underset{¯}{\alpha },\underset{¯}{\beta }}{\min }\left({.Math.V\underset{¯}{\alpha }-W\underset{¯}{\beta }.Math.}^2+\mu {.Math.\left[\begin{array}{cc}& {\Psi }^{*}\underset{¯}{x}\end{array}\right]\left[\begin{array}{c}\underset{\_}{\alpha }\\ -1\end{array}\right].Math.}^2\right)s.t.{.Math.\underset{¯}{\alpha }.Math.}^2=c>0.\mathrm{Now}\mathrm{let}& \left(A\mathrm{.1}\right)\\ \underset{¯}{\tilde{\alpha }}=\left[\begin{array}{c}\underset{\_}{\alpha }\\ -1\end{array}\right].Math.\underset{¯}{\alpha }=E\underset{¯}{\tilde{\alpha }},\mathrm{where}& \left(A\mathrm{.2}\right)\\ E=\left[\begin{array}{cc}{I}_{p^{2}M}& \underset{\_}{O}\end{array}\right]\in {ℝ}^{p^{2}M\times \left(p^{2}M+1\right)}.\mathrm{Also},\mathrm{let}& \left(A\mathrm{.3}\right)\\ R=\left[\begin{array}{cc}& {\Psi }^{*}\underset{¯}{x}\end{array}\right],& \left(A\mathrm{.4}\right)\end{array}$
and rewrite (A.1) as

(244) $\begin{array}{cc}\underset{\widetilde{\underset{¯}{\alpha }},\underset{¯}{\beta }}{\min }\left({.Math.VE\underset{¯}{\tilde{\alpha }}-W\underset{¯}{\beta }.Math.}^2+\mu {.Math.R\widetilde{\underset{¯}{\alpha }}.Math.}^2\right)s.t.{.Math.\widetilde{\underset{¯}{\alpha }}.Math.}^2=\tilde{c}>0.& \left(A\mathrm{.5}\right)\end{array}$
Let ƒ({tilde over (α)}β)=∥VE{tilde over (α)}−Wβ∥.sup.2+μ∥R{tilde over (α)}∥.sup.2, and recall (54), (55). Therefore, the objective function ƒ ({tilde over (α)},β) can be rewritten as
ƒ({tilde over (α)},β)={tilde over (α)}*(E*E+μR*R){tilde over (α)}−2{tilde over (α)}*E*V*Wβ+∥β∥.sup.2.   (A.6)

(245) Also, let
g({tilde over (α)},β)=∥{tilde over (α)}∥.sup.2−{tilde over (c)}  (A.7)

(246) Since we only have one constraint (A.7), there is only one Lagrange multiplier, which we denote with λ, and the Lagrangian is therefore
L({tilde over (α)},β,λ)=ƒ({tilde over (α)},β)−λg({tilde over (α)},β).  (A.8)

(247) Now the gradient of ƒ({tilde over (α)},β), with respect to {tilde over (α)} and β is

(248) $\begin{array}{cc}{\nabla }_{\widetilde{\underset{¯}{\alpha }},\underset{¯}{\beta }}f=\left[\begin{array}{c}2\left(E^{*}E+\mu R^{*}R\right)\widetilde{\underset{¯}{\alpha }}-2E^{*}{V}^{*}W\underset{¯}{\beta }\\ 2\left(\underset{¯}{\beta }-W^{*}VE\underset{¯}{\tilde{\alpha }}\right)\end{array}\right]\mathrm{Also},& \left(A\mathrm{.9}\right)\\ {\nabla }_{\underset{¯}{\tilde{\alpha }},\underset{¯}{\beta }}g=\left[\begin{array}{c}2\underset{\_}{\tilde{\alpha }}\\ \underset{\_}{O}\end{array}\right]\mathrm{Therefore},& \left(A\mathrm{.10}\right)\\ {\nabla }_{\widetilde{\underset{¯}{\alpha }},\underset{¯}{\beta ,}\lambda }ℒ=2\left[\begin{array}{c}\left(E^{*}E+\mu R^{*}R\right)\widetilde{\underset{¯}{\alpha }}-E^{*}{V}^{*}W\underset{¯}{\beta }-\lambda \widetilde{\underset{¯}{\alpha }}\\ \underset{¯}{\beta }-W^{*}VE\widetilde{\underset{¯}{\alpha }}\\ 0.5\left(\tilde{c}-{.Math.\underset{¯}{\tilde{\alpha }}.Math.}^2\right)\end{array}\right].& \left(A\mathrm{.11}\right)\\ {\nabla }_{\widetilde{\underset{¯}{\alpha }},\underset{¯}{\beta },\lambda }ℒ=\underset{¯}{O}.Math.\left(E^{*}E+\mu R^{*}R\right)\underset{¯}{\tilde{\alpha }}-E^{*}{V}^{*}W\underset{¯}{\beta }=\lambda \widetilde{\underset{¯}{\alpha }}& \left(A\mathrm{.12}\right)\\ \underset{¯}{\beta }=W^{*}VE\widetilde{\underset{¯}{\alpha }},& \left(A\mathrm{.13}\right)\\ {.Math.\widetilde{\underset{¯}{\alpha }}.Math.}^2=\tilde{c}.& \left(A\mathrm{.14}\right)\end{array}$

(249) Combining (A.12), (A.13), we get
(E*E+μR*R−E*V*WW*VE){tilde over (α)}=λ{tilde over (α)}.  (A.15)

(250) We can now obtain the value of the objective function, by plugging (A.13) in (A.6), and then using (A.15) and (A.14)

(251) $\begin{array}{cc}{\widetilde{\underset{¯}{\alpha }}}^{*}\underset{\underset{\lambda \widetilde{\underset{¯}{\alpha }}}{︸}}{\left(E^{*}E+\mu R^{*}R-E^{*}{V}^{*}W{W}^{*}VE\right)\widetilde{\underset{¯}{\alpha }}}=\lambda \underset{\underset{\tilde{c}}{︸}}{{.Math.\widetilde{\underset{¯}{\alpha }}.Math.}^2}=\lambda \tilde{c}.& \left(A\mathrm{.16}\right)\end{array}$

(252) Now note that

(253) $\begin{array}{cc}{E}^{*}E=\left[\begin{array}{cc}{I}_{p^{2}M}& \underset{\_}{O}\\ {\underset{\_}{O}}^{*}& 0\end{array}\right],E^{*}{V}^{*}W{W}^{*}\mathrm{VE}=\left[\begin{array}{cc}{V}^{*}W{W}^{*}V& \underset{\_}{O}\\ {\underset{\_}{O}}^{*}& 0\end{array}\right],\mathrm{and}& \left(A\mathrm{.3}\right)\\ R^{*}R=\left[\begin{array}{cc}{*}^{}& {*}^{}{\Psi }^{*}\underset{¯}{x}\\ {\left({*}^{}{\Psi }^{*}\underset{¯}{x}\right)}^{*}& {.Math.{\Psi }^{*}\underset{¯}{x}.Math.}^2\end{array}\right],& \left(A\mathrm{.4}\right)\end{array}$
and therefore equation (A.15) is rewritten as

(254) 0$\begin{array}{cc}\left[\begin{array}{cc}\mu {*}^{}+I-V^{*}W{W}^{*}V& \mu {*}^{}{\Psi }^{*}\underset{¯}{x}\\ {\mu \left({*}^{}{\Psi }^{*}\underset{¯}{x}\right)}^{*}& \mu {.Math.{\Psi }^{*}\underset{¯}{x}.Math.}^2\end{array}\right]\widetilde{\underset{¯}{\alpha }}=\lambda \widetilde{\underset{¯}{\alpha }}.& \left(A\mathrm{.17}\right)\end{array}$

(255) Since any singular pair of the (symmetric) matrix in (A.17) is a solution, and since we seek to minimize the error (A.16), we pick the solution

(256) $\widehat{\tilde{\underset{¯}{\alpha }}}=\sqrt{\tilde{c}}\underset{¯}{v},$
where ν is the last singular vector of the matrix in (A.17), with associated (smallest) singular value σ.sub.ν (and thus

(257) $\sqrt{\tilde{c}}=-\frac{1}{e},$
where e is the last element in v).

(258) Now recall (A.2) and rewrite (A.17) as

(259) $\begin{array}{cc}\left[\begin{array}{cc}\mu {*}^{}+I-V^{*}W{W}^{*}V& \mu {*}^{}{\Psi }^{*}\underset{¯}{x}\\ {\mu \left({*}^{}{\Psi }^{*}\underset{¯}{x}\right)}^{*}& \mu {.Math.{\Psi }^{*}\underset{¯}{x}.Math.}^2\end{array}\right]\left[\begin{array}{c}\underset{\_}{\hat{\alpha }}\\ -1\end{array}\right]={\sigma }_v\left[\begin{array}{c}\underset{\_}{\hat{\alpha }}\\ -1\end{array}\right].& \left(A\mathrm{.18}\right)\end{array}$

(260) The top part of equation (A.18) is
(μJ*J+I−V*WW*V){circumflex over (α)}−μJ*Ψ*x=σ.sub.ν{circumflex over (α)}(A.19)
and therefore
{circumflex over (α)}=μ[μJ*J+(1−σ.sub.ν)I−V*WW*V].sup.−1(J*Ψ*x).  (A.20)

(261) Now, for numerical stability reasons, the computation of the smallest singular value of the matrix in A.17 might be inaccurate (because the last row and the last column in the matrix in A.17 have much higher values compared to μJ*J+I−V*WW*V), so instead of using Ψ*x in optimization problem A.1 we use its normalized version

(262) $\frac{{\Psi }^{*}\underset{¯}{x}}{.Math.{\Psi }^{*}\underset{¯}{x}.Math.}.$
So equation A.18 becomes

(263) $\begin{array}{cc}[\begin{array}{cc}\mu {*}^{}+I-V^{*}W{W}^{*}V& \frac{\mu }{.Math.{\Psi }^{*}\underset{¯}{x}.Math.}{*}^{}{\Psi }^{*}\underset{¯}{x}\\ \frac{\mu }{.Math.{\Psi }^{*}\underset{¯}{x}.Math.}{\left({*}^{}{\Psi }^{*}\underset{¯}{x}\right)}^{*}& \mu \end{array}]\left[\begin{array}{c}\underset{\_}{\hat{\alpha }}\\ -1\end{array}\right]=\sigma \left[\begin{array}{c}\underset{\_}{\hat{\alpha }}\\ -1\end{array}\right].& \left(A\mathrm{.21}\right)\end{array}$
where σ is the smallest singular value the matrix in A.21. Therefore, the corresponding final solution (after scaling back by ∥Ψ*x∥) is
{circumflex over (α)}=μ[μJ*J+(1−σ)I−V*WW*V].sup.−1(J*Ψ*x).  (A.22)

System Requirements for Implementing Algorithm on Laptop

(264) Laptop information:

(265) OS Name Microsoft Windows 10 Home Version 10.0.17134 Build 17134

(266) System Manufacturer Dell Inc.

(267) System Model G3 3579

(268) System Type x64-based PC

(269) Processor Intel® Core™ i7-8750H CPU @ 2.20 GHz, 2208 Mhz, 6 Core(s).

(270) Installed Physical Memory (RAM) 16.0 GB

(271) Matlab information:

(272) MATLAB Version 9.4 (R2018a)

(273) Add-on toolbox: Parallel Computing Toolbox Version 6.12 (R2018a)

(274) This embodiment of the invention uses the very basic Matlab package+the parallel computing toolbox to take advantage of the laptop's 6-core CPU. Matlab is known as a scientific computing programming language. The computations in the algorithm are basic, so it can be done using other (free) programming platforms such as Anaconda (the Python distribution) or C or Java, etc. The 100 frame Suzie sequence in the experiment section can be upsampled using Matlab+the listed laptop in a couple of hours (or even half an hour if we change some parameters and sacrifice a bit of quality for speed). Nevertheless, the algorithm's complexity is considered comparatively low (compared to the norm in the field of upsampling), but it is still too high for real-time/low-power applications, which makes optimization for faster computations potentially very useful. One might also want to consider other solutions (beyond computational optimization). For instance, images captured by a low power device (e.g. an iPhone or a satellite) can be dealt with using cloud computing. As much as ×10-100 speed up can be obtained by converting Matlab code to C code. The algorithm is parallelizable, so besides writing the algorithm in C, using a GPU (or the cloud) can considerably speed up calculations.