^{1}

^{1}

^{*}

^{2}

^{3}

^{1}

^{2}

^{3}

Edited by: Ke Shi, Old Dominion University, United States

Reviewed by: Jian Lu, Shenzhen University, China; Xin Guo, Hong Kong Polytechnic University, Hong Kong

This article was submitted to Mathematics of Computation and Data Science, a section of the journal Frontiers in Applied Mathematics and Statistics

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Image super-resolution is an image reconstruction technique which attempts to reconstruct a high resolution image from one or more under-sampled low-resolution images of the same scene. High resolution images aid in analysis and inference in a multitude of digital imaging applications. However, due to limited accessibility to high-resolution imaging systems, a need arises for alternative measures to obtain the desired results. We propose a three-dimensional single image model to improve image resolution by estimating the analog image intensity function. In recent literature, it has been shown that image patches can be represented by a linear combination of appropriately chosen basis functions. We assume that the underlying analog image consists of smooth and edge components that can be approximated using a reproducible kernel Hilbert space function and the Heaviside function, respectively. We also extend the proposed method to pansharpening, a technology to fuse a high resolution panchromatic image with a low resolution multi-spectral image for a high resolution multi-spectral image. Various numerical results of the proposed formulation indicate competitive performance when compared to some state-of-the-art algorithms.

During image acquisition, we often loose resolution due to the limited density of the imaging sensors and the blurring of the acquisition lens. Image super resolution, a post process to increase image resolution, is an efficient way to achieve a resolution beyond what the hardware offers. A detailed overview of this area is outlined in Farsiu et al. [

In recent literature, image super-resolution (SR) methods typically fall into three categories: interpolation based methods, learning based methods, and constrained reconstruction methods. Interpolation based methods [

Increased accessibility to satellite data has led to a wide range of applications in multispectral (MS) SR. These applications include target detection, scene classification, and change detection. A multitude of algorithms have been presented in recent literature to address this task. However, the performance of many of these algorithms are data dependent with respect to existing quantitative assessment measures. With the increasing advancements in technology, there have been varying modes of image acquisition systems. In remote sensing, satellites carry sensors that acquire images of a scene in multiple channels across the electromagnetic spectrum. This gives rise to what is referred to as

Most satellites also acquire a high-spatial resolution image known as a

The main assumption of CS methods is that the spatial and spectral information of the high spectral resolution image can be separated in some space [

In this paper, we extend the key ideas in the SR model in Deng et al. [

In Deng et al. [

Using the same basis as in Deng et al. [

In this paper, we use the RKHS approximations to model the smooth components of the image whiles preserving the non-smooth components by the approximated Heaviside function. For the pansharpening experiments, We also preprocess the data using fast state-of-the-art MRA techniques [

The experiments that follow show that these contributions yield better enhanced resolution results with faster convergence speeds at all scales.

The paper is organized as follows. In section 2, we review the regular functions used in decomposing the image. The proposed model is described in section 3 with numerical results outlined in section 4.

In this section, we review the mathematical functions used in the formulation of our model. Similar as in Deng et al. [

RKHS can be used to define feature spaces to help compare objects that have complex structures and are hard to distinguish in the original feature space. Wahba [

Let

By Taylor's theorem with remainder for

for some _{+} = _{+} = 0 otherwise. Let

and ^{m} denotes the

therefore for any

Let

The space

with norm

so that for any

The reproducing kernel for

with norm

therefore for any _{0} + _{1} with

Let _{jν} = ϕ_{ν}(_{j}), _{ji} = ξ_{i}(_{j}), then

where

where the term μ^{T}_{1}(_{2}(

Extending the one-dimensional spline model to two-dimensional data, we consider a similar discretization and set up a model similar to that defined in (11). Let _{i} ∈ [0, 1] × [0, 1] be a discretization for an observed noisy image. From Wahba [

where

The null space of _{m} is the _{m} is spanned by monomials ϕ_{1}, ϕ_{2}, …, ϕ_{6} given by

Duchon [_{μ} exists for (12) with representation

where _{m} is defined as

and

Based on the work of Duchon [

where _{i, ν} = ϕ_{ν}(_{i}) and _{i, j} = _{m}(_{i}, _{j}). _{m}(_{i}(

The one-dimensional Heaviside step function is defined as

Due to the singularity at

We refer to this as the approximated Heaviside function (AHF), an approximation to ϕ(

with variables θ and

Illustrating the surfaces corresponding to the approximated Heaviside functions for varying pairs of θ and

We define

Thus, given the low resolution input ^{M} and

where _{1} to enforces the sparsity of as all the orientations may not be present for any given image. Once the coefficients have been obtained, we have

where

In this section, we present our proposed models for SISR and pansharpening based on the functions defined in section 2. We propose a 3D patch-based approach that infers the HR patches from LR patches that has been grouped into classes based on their structural similarity. As a result, we impose similarity constraints within the classes so that the coefficients for neighboring patches in a group do not differ substantially. In addition, we use sparse constrained optimization techniques that simplify the minimization of the resulting energy functional by solving a series of subproblems, each with a closed form solution. Two algorithms will be discussed. The first is a 3D SR model for true color images and MS images with at most four bands. The second algorithm is an extension of our SISR model to pansharpening for MS images with at least four spectral bands. This is a two-step hybrid approach that incorporates CS and MRA methods into the first algorithm to obtain an enhanced resolution result.

Given a LR image _{i} = ω_{i}, where ω is a scale factor and λ is the number of bands as follows. We begin by decomposing the LR image into a set of overlapping patches _{p}. Next, for a suitably chosen

Clustering the set of patches

Following the clustering, we consider each set of patches _{i} with 1 ≤

where _{p} = ω_{p} for some scale factor ω ∈ ℕ, _{i,pj} and _{i}. Here, _{Ii} represents the total number of patches in the set

where ^{ℓ} and ^{h}. The coefficient matrices are obtained by solving the following minimization problem for each patch indexed by

with tr(·) as the trace and adaptive weights

where _{1}, …_{λNIi}] and μ_{1}, μ_{2}, μ_{3} ≥ 0. Note that _{i} are likely to vary smoothly within the same class of images. When _{i} and _{j} fall into the same class, _{ij} tends to be larger, forcing the next iterate of _{i} and _{j} to be close as well. The second regularization term is a structure guided regularity following from (16) to guarantee that the coefficients of the residual components do not grow too large. Finally, we impose sparsity constraints using the ||·||_{1, 1} norm [

since we assume that edges are sparse in the image patches. The resulting minimization can be solved using

Due to the non-differentiability of the ||·||_{1, 1} norm, we introduce a new variable

By introducing the augmented Lagrangian and completing the squares, we change the constrained optimization problem (29) into the following non-constrained optimization where

with step length γ ≥ 0. ADMM reduces (30) to solving the following separable subproblems:

Given that _{1}, …_{λNIi}], _{1}, …_{λNIi}], _{1}, …_{λNIi}] and we can solve for each column of

where _{(λ,Ii)} matching the column index of _{j}. This gives a solution

where ^{ℓ}.

From which we get

Similarly, this gives the solution

A solution to this problem arising from the group sparsity constraint can be obtained by applying ℓ_{1}

where

for any ^{n} and γ > 0.

The update for

The estimated set of patches

Analog three-dimensional single image SR model

Given a LR MS image

where

Pansharpening model.

To begin with the simulation procedure, we test the performance of the proposed algorithm by generating the reduced resolution image from a given image for some chosen downscaling factors. We then compute the enhanced resolution result using the proposed algorithm and compare its performance with other methods. The experiments are implemented using MATLAB(R2018a).

For the SISR experiments we compare our three-dimensional single image super-resolution algorithm to bicubic interpolation and the RKHS algorithm by Deng et al. [

We undertake the pansharpening experiments using semi-synthetic test data. This may be attributed to the fact that existing sensors cannot readily provide images at the increased spatial resolution that is desired. A generally acceptable procedure attributed to Ranchin and Wald [

The enhanced result should be as identical as possible to the original multispectral image once it has been degraded.

The enhanced result should be identical to the corresponding high resolution image that a capable sensor would observe at the increased resolution.

Each band of the enhance result should be as identical as possible to the corresponding bands of the image that would have been observed by a sensor with high resolution.

Flow diagram of the experimental procedure for synthetic datasets.

We apply Algorithm 2 to the observed multispectral data and compare it with existing state-of-the-art approaches.

Parameter selection is essential to obtain good results. However, to show the stability of our proposed method, we fix the parameters for all the experiments. Fine tuning of the parameters might lead to slightly better results. Parameter choices for the optimization of our proposed algorithms are outlined as follows:

Scale factor (ω): In the experiments shown above, we set ω = 2 for the SISR problem, and set ω = 4 for the pansharpening problem. Quantitative measures obtained suggest that for ω > 4 both problems yield unsatisfactory results as the spatial and spectral measures are less competitive. However, the choice of ω depends on the dimensions of the observed data. For large ω, we can select larger patch sizes to improve upon the target.

Patch size (_{p}): Patch sizes can vary according to the size of the given image. We use a default patch size, _{p} × _{p} × λ in our experiments. Reasonable quantitative measures were obtained for 6 ≤ _{p} ≤ 8 for the two problems considered. Outside the given bounds, the algorithms either take time to converge or yield unsatisfactory results.

Overlap size: The limits of the overlapping region range from 2 to 4 are limited by the patch size. We set the overlap size to 2.

Number of classes (

Construction of approximation matrices ^{−3}, and set 20 evenly distributed angles on [0, 2π] in computing the fine and coarse scale matrices for

Regularization penalties μ_{1}, μ_{2}, μ_{3} and γ: Without knowledge of optimal values for the penalty, we sweep through a range of values to determine the optimal result. For the experiments considered, μ_{1}, μ_{2}, μ_{3}, γ ∈ {10^{−8}, ⋯ , 10^{2}} and choose the one that leads to the best results evaluated by quantitative and qualitative measures. The step length γ for updating of

Preprocessing for Algorithm 2: Our proposed pansharpening method provides a flexible way to improve the results for both CS and MRA algorithms. In the experiments shown above, we preprocess the images using the GS algorithm. Future work will detail which preprocessing step works best for multispectral images.

To measure the similarity of the enhanced resolution image (

Correlation coefficient (CC): The correlation coefficient measures the similarity of spectral features between the reference and pansharpened image. It is defined as

where _{i} and _{Xi} and _{i} are _{Xi}, respectively.

Root mean square error (RMSE):

where

Peak signal-to-noise ratio (PSNR): It is an expression of the ratio between the maximum possible value of the signal and the distorting noise that affects the quality of the representation.

Structural Similarity (SSIM): The SSIM index is computed by examining various windows of the reference and target image. It is a full reference metric designed to improve upon the traditional methods such as PNSR and MSE.

Spectral angle mapper (SAM): This measure denotes the absolute value of the angle between the reference and estimated spectral vectors.

Let

Erreur relative globale adimensionnelle de synthése (ERGAS): The relative dimensionless global error in the synthesis reflects the overall quality of the pansharpened image.

where _{i} and

In this section, we compare the visual results and quantitative measures of the proposed approaches with some existing super-resolution algorithms. The results of Algorithm 1 will be compared to bicubic interpolation and the RKHS model by Deng et al. [

In

Decomposition of components for the grayscale test image. From top to bottom and left to right: the original, ^{h}^{h}^{h}

Results for

Quantitative Results with for single image super-resolution with ω = 2.

Bicubic | 22.6453 | 0.0737 | 0.9700 | 4.9851 | |

RKHS [ |
22.8789 | 0.0718 | 0.9770 | 5.1105 | |

Proposed | 23.1566 | 0.0695 | 0.9780 | 4.9757 | |

Bicubic | 29.7291 | 0.0326 | 0.9453 | 1.3056 | |

RKHS [ |
29.9922 | 0.0317 | 0.9564 | 1.3559 | |

Proposed | 30.7616 | 0.0290 | 0.9621 | 1.2173 | |

Bicubic | 32.3586 | 0.0241 | 0.9740 | 1.5280 | |

RKHS [ |
32.6024 | 0.0234 | 0.9810 | 1.5230 | |

Proposed | 33.2650 | 0.0217 | 0.9890 | 1.4663 | |

Bicubic | 29.6166 | 0.0330 | 0.9765 | 1.3992 | |

RKHS [ |
30.1480 | 0.0311 | 0.9788 | 1.6617 | |

Proposed | 31.2364 | 0.0274 | 0.9801 | 1.3770 | |

Bicubic | 31.4943 | 0.0266 | 0.9807 | 2.2896 | |

RKHS [ |
31.6163 | 0.0263 | 0.9820 | 2.6242 | |

Proposed | 31.9927 | 0.0251 | 0.9829 | 2.4506 | |

∞ | 0 | 1 | 0 |

For pansharpening, we compare with 10 other methods in the literature using five different quantitative metrics (see

Quantitative Results with for pansharpening ω = 4.

PCA | 5.01 | 0.93 | 38.19 | 4.40 | 30.14 | |

IHS | 5.29 | 0.93 | 35.72 | 4.17 | 30.26 | |

Brovey | 4.66 | 0.94 | 35.35 | 4.07 | 30.55 | |

GS | 4.75 | 0.96 | 34.48 | 4.00 | 30.86 | |

Indusion | 4.52 | 0.92 | 37.47 | 4.36 | 30.11 | |

HPF | 3.64 | 0.96 | 27.48 | 3.23 | 32.52 | |

SFIM | 3.61 | 0.96 | 27.18 | 3.17 | 32.75 | |

AWLP | 3.29 | 0.96 | 23.63 | 2.84 | 33.61 | |

GLP | 3.23 | 0.97 | 2.61 | |||

RKHS [ |
23.32 | 34.15 | ||||

Proposed | ||||||

PCA | 3.52 | 0.96 | 20.34 | 2.64 | 33.25 | |

IHS | 5.04 | 0.89 | 27.10 | 3.56 | 29.72 | |

Brovey | 4.49 | 0.90 | 24.12 | 3.15 | 30.63 | |

GS | 3.51 | 0.96 | 20.82 | 2.76 | 32.70 | |

Indusion | 3.89 | 0.94 | 25.28 | 3.25 | 31.54 | |

HPF | 3.32 | 0.96 | 20.84 | 2.69 | 33.16 | |

SFIM | 3.27 | 0.96 | 20.60 | 2.65 | 33.27 | |

AWLP | 0.96 | 33.61 | ||||

GLP | 3.27 | 0.96 | 20.46 | 2.71 | 32.78 | |

RKHS [ |
3.06 | 19.52 | 2.56 | |||

Proposed | ||||||

PCA | 4.83 | 0.90 | 28.39 | 4.75 | 36.82 | |

IHS | 4.96 | 0.87 | 28.73 | 4.78 | 36.36 | |

Brovey | 5.23 | 0.85 | 29.62 | 4.91 | 36.00 | |

GS | 4.76 | 0.90 | 28.12 | 4.71 | 36.90 | |

Indusion | 4.54 | 0.83 | 29.41 | 4.88 | 36.33 | |

HPF | 3.99 | 0.91 | 23.13 | 3.89 | 38.65 | |

SFIM | 3.95 | 0.92 | 22.76 | 3.83 | 38.78 | |

AWLP | 4.31 | 0.92 | 22.73 | 3.85 | 39.03 | |

GLP | 3.76 | |||||

RKHS [ |
0.92 | 21.47 | 3.58 | |||

Proposed | 39.28 | |||||

PCA | 11.66 | 0.85 | 0.06 | 14.05 | 28.34 | |

IHS | 6.30 | 0.69 | 0.04 | 18.67 | 25.49 | |

Brovey | 0.77 | 0.05 | 16.46 | 25.88 | ||

GS | 5.77 | 0.81 | 0.05 | 15.30 | 26.52 | |

Indusion | 4.64 | 0.88 | 0.04 | 12.70 | 28.23 | |

HPF | 4.53 | 0.89 | 0.04 | 12.30 | 28.52 | |

SFIM | 4.68 | 0.69 | 0.19 | 102.38 | 22.38 | |

AWLP | 8.39 | 0.86 | 0.06 | 15.04 | 27.57 | |

GLP | 4.73 | 0.88 | 0.04 | 12.74 | 28.21 | |

RKHS [ |
3.94 | |||||

Proposed | ||||||

PCA | 4.50 | 0.96 | 62.23 | 4.31 | 30.98 | |

IHS | 5.04 | 0.95 | 65.75 | 4.48 | 30.60 | |

Brovey | 4.86 | 0.95 | 65.01 | 4.49 | 30.43 | |

GS | 4.52 | 0.96 | 62.52 | 4.33 | 30.94 | |

Indusion | 4.56 | 0.91 | 74.77 | 5.14 | 29.38 | |

HPF | 4.10 | 0.95 | 57.14 | 3.96 | 31.72 | |

SFIM | 4.14 | 0.95 | 56.80 | 3.94 | 31.78 | |

AWLP | 4.62 | 0.95 | 52.82 | 3.88 | 32.26 | |

GLP | 3.90 | 50.92 | 3.55 | 32.70 | ||

RKHS [ |
0.96 | |||||

Proposed | ||||||

0 | 1 | 0 | 0 | ∞ |

Note that we did not conduct comparison with deep neural network based methods. The proposed approach is based on single image modeling while the deep neural network approach relies on a lot of external data. We don't think it is a fair comparison.

In this paper, we have proposed a technique for single image SR by modeling the image as a linear combination of regular functions in tandem with sparse coding. We show that the proposed scheme is an improvement upon similar existing approaches as it outperforms these algorithms. Besides this advantage, the proposed methods can also benefit image decomposition.

In future work, we will apply the model to multispectral and hyperspectral image SR where sparse coding can be useful due to the redundant nature of the input.

Some of the datasets for this study can be found in

WG contributes on providing the main ideas including the mathematical models and algorithms proposed in the paper. RL was a former Ph.D. student of WG. Under WG's guidance, RL did computer programming to implement the ideas and conducted numerical experiments to test the effective of the proposed method compared to the state of the art. XZ and CG are domain scientists who help us evaluate the numerical results and provide comments.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.