{'fraction_non_alphanumeric': 0.04002672903441363, 'fraction_numerical': 0.012094888072168393, 'mean_word_length': 4.594579439252336, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Super-resolving multiresolution images with band-independant geometry of multispectral pixels

A new resolution enhancement method is presented for multispectral and multi-resolution images, such as these provided by the Sentinel-2 satellites. Starting from the highest resolution bands, band-dependent information (reflectance) is separated from information that is common to all bands (geometry of scene elements). This model is then applied to unmix low-resolution bands, preserving their reflectance, while propagating band-independent information to preserve the sub-pixel details. A reference implementation is provided, with an application example for super-resolving Sentinel-2 data.


I. INTRODUCTION
A. Context and state of art Earth Observation missions typically operate at medium to low resolution ranges in order to favor both larger satellite swath and better temporal revisit of the same site (e.g. 3-4 days over Europe for the Sentinel-2 Satellite series). For each acquisition, optical constraints furthermore often restrict that only some spectral bands have maximal resolution. For example, a common case is to compensate the smaller pixel size of the higher spatial resolution bands (e.g. 10m/pixel) by capturing light over a larger spectrum range (e.g. 4 large bands in the red, green, blue and near infrared for Sentinel-2). The limit case being a single high-resolution panchromatic band (Pleiades, Spot, Landsat...). Narrower spectral bands, invaluable for specific measurements (e.g. chlorophyll or water vapor absorbtion wavelengths) are then provided at lower resolution (e.g. 20m/pixel and 60m/pixel for Sentinel-2). Yet, this trade-off on spectral and spatial resolutions may become a limiting factor for many Earth Observation applications, for example for getting accurate land cover classification at the highest resolution [1]. Some techniques have thus been devised in order to propagate the high-resolution spatial details to the lower-resolution dedicated bands while preserving their spectral content. These can be sorted along the following categories: -Probabilistic [7], [8]: The spectral information in each sub-pixel of an original low-resolution pixel is determined by maximizing a probabilistic model, constrained by the observed data at all bands and resolutions (possibly including the panchromatic band). A Bayesian formulation can be choosen to represent this constraint which, provided this does not becomes intractable, allows hyperparameters to be set according to prior knowledge.
-Sensor-based: If the sensor has a known point spread function (PSF), then deconvoluting it enhances the resolution of the acquired images [9]. But the PSF for many satellites can only be estimated empirically (Sentinel-2, Spot-5, Landsat-8 [1]). When that is the case, limits on sub-pixel detection can be established [1].
-Learning-based: These methods exploit local patterns in the low-resolution images, and propagate these features (e.g. edges) to infer the higher resolution image [10]. Many models may be used to "learn" the features: neural network [11], example-based [12] with kernel ridge regression [13], deep learning [14] and more references therein, including for cross-image learning. These methods can be applied for single image resolution enhancement, possibly with different channels (typically red, green, blue [14]). Filling details from learned (or duplicated) texture features might be very good to produce visually plausible results [12]. Their main problem, similar to in-painting with image-based examples [15], is that "hallucinated" [10] details do not necessarily correspond to true higher-resolution objects (esp. with non-local or crossimage features) and then become misleading pixels for land monitoring purposes.
-Scaling laws: Instead of learning local patterns, this method learns multi-scale relationships in the data such as local power laws between spatial extent and band values [18]. Such scailing laws are inferred from data above the aquisition resolution but, assuming the same laws remain valid below that resolution, these can then be used to infer sub-pixels of the original image. Very good results have been obtained with such methods for turbulent oceanic data [19], where energy cascades translates to power laws spanning multiple decades. However, for land monitoring purposes, usually no such physical interpretation can be found: for example a mixture of trees, houses and roads in a peri-urban environment is not locally scale-invariant.
-Frequency representations: Working in the frequency domain, whether with Fourier methods or using wavelet decompositions. An idea is to upsample the low-frequencies (e.g. with a bicubic filter) and preserve the high frequency components. Unfortunately, natural images are not statistically consistent: knowledge that there is a tree (or road, house...) a few hundred meters away (i.e. the wavelet support size) does not help subdivide a local pixel into its higher-resolution components.
-Panchromatic sharpening [2]: Using a very wide band with high-resolution in order to compensate the lower resolution of narrow bands. Multiple variants exist, from a simple renormalization of the multispectral bands [3] to more advanced unmixing techniques which estimate the contribution of each spectral band to the panchromatic one [4], possibly with preprocessing steps designed to uncorrelate each component [5], [2], or using angular spaces such as the hyperspherical color space [6].
The advantages of panchromatic sharpening are its efficiency, and its applicability even when only a single highresolution band is acquired. Many observation satellites thus include a panchromatic band, but some do not (e.g. ESA's Sentinel-2 series, specifications given in Appendix). In the absence of such a band, and given the inadequacies of the other methods presented above, another solution is needed.
B. Super-resolution of multispectral images in the absence of a panchromatic band In a multispectral measurement system, each pixel in a band B captures the light intensity over some part of the spectrum, according to some sensor sensitivity distribution s B (λ) for each wavelength λ: B(x, y) = s B (λ)I(λ, x, y)dλ. The light intensity I is reflected at wavelength λ by the pixel surface between x + r B and y + r B , where r B is the square pixel resolution of band B. Larger pixels thus collect more light and may be necessary for some bands with narrow s B spectral support. This very simplified model ignores many sources of optical distortions and the post-processing from satellite geometry to square pixels, but these effects are irrelevant in the context of this section.
When a panchromatic band P is available, light is collected over a wide spectrum support : s P (λ) = 0 for a large range of wavelengths [λ P min , λ P max ], usually covering all other bands λ P min ≤ λ Bmin ≤ λ Bmax ≤ λ P max . Collecting more light spectrally allows to reduce the pixel size r P < r B while still maintaining a minimal intensity to trigger the panchromatic sensor. Each band pixel B(x, y) thus covers an area with multiple panchromatic pixels P (x + i * r P , y + j * r p ) for i, j indices depending on the band resolution r B . The total light reaching the sensor for B(x, y) may thus be related to the values of each sub-pixels P (x+i * r P , y+j * r p ) by a function involving both spatial and spectral components. Panchromatic sharpening methods attempt to unmix this total light contribution by defining B(x, y) = f (H(x+i * r p , y+i * r p )), where H is a high-resolution version of band B. f is a coarse-graining function, for example f = i,j in the simplest case. More elaborated functions may be necessary, for example when combining both panchromatic sharpening and atmospheric corrections. The relations between s P (λ) and s B (λ) are used in order to further constrain H with values of P at the same pixel locations.
Unfortunately, for Sentinel-2 images, there is no panchromatic band. Methods using such bands, cited in the previous Section, are thus not directly useable. An idea would be to create a virtual P band by combining the four high-resolution bands B 10 at 10m/pixel, and then apply pansharpening methods for bands at 20m/pixel and 60m/pixel. But combining blue, green, red, and infrared bands into a virtual intensity would result in a virtual s P (λ) that only covers the original s B 10 (λ), and thus prone to spectral artifacts for super-resolving the values of the other bands B 20 and B 60 . Similarly, spatial details in that hypothetical P band would depend on how that virtual P is constructed from other B 10 bands (in particular, how infrared details are fused with visible light details).
In the absence of a real panchromatic band, a better idea is to design a new method that: -Explicitely encodes geometric details from available highresolution bands, as pixel properties independent from their reflectance; -Preserves the spectral content of each low-resolution band independently from the geometry.
The next sections introduce one method for acheiving this, with direct applicability to Sentinel-2 images. The method presented below works with only local information, hence it is not subject to the non-local effects mentionned in the previous section. The method relies on the observation that the proportion of objects of the same nature within a pixel area (e.g. 30% urban area, 70% trees), is a physical property of that pixel and therefore independent of the spectral band. Only the reflectance of these objects may change from band to band. Moreover, there is no reason why pixel boundaries would match natural object boundaries. The method thus identifies generic "shared" information between adjacent pixels, then commonalize the geometric aspects of these shared values across bands. High resolution bands are used to separate bandindependent information from band-dependant reflectance. The geometric information is then used to unmix the lowresolution pixels, while preserving their overall reflectance.
Section II presents the super-resolution problem and Section III details how that problem is addressed by the model introduced in this paper. Section IV indicates how to quantify the results quality, and section V shows super-resolution results for three different types of regions of interest (coastal, urban, agricultrural). Results are followed by a discussion in Section VI, which also demonstrates the influence of each step of the algorithm.

A. Super-resolution formulation
Let L be an observed low resolution image with N x /2 columns and N y /2 rows. We consider the problem of finding a high resolution image H with N x columns and N y rows. Each low resolution pixel L x/2,y/2 thus corresponds to 4 high resolution pixels, as depicted in Fig. 1. Averaging these pixels should give the original observed low-resolution pixel back: Images remotely sensed from satellites are subject to multiple transforms (including atmospheric corrections [24]) before being released as a useable product. These transforms are out of scope of the present document but may introduce correlations between high-resolution pixels (e.g. due to scattering), hence should be applied before (or integrated to) superresolution. Similarly, known PSF [9] should be deconvoluted in addition to the method presented below. In any case, Eq. 1 ensures that down-sampling by averaging the high-resolution solution will recover the observations. Eq. 1 is undetermined, with 3 free parameters per low resolution pixel. Some extra constraints are needed, which are extracted from available high-resolution data bands.

B. Shared information between neighbor pixels
Natural objects do not fall exactly on pixel boundaries. Therefore, some content is shared between nearby pixel values. This shared information is explicitly defined as in Fig. 2, left. For example, S x+1,y+1 is the reflectance corresponding to the shared part between high-resolution pixels H x,y , H x+1,y , H x,y+1 and H x+1,y+1 (compare Fig. 1 and Fig. 2). This particular element is fully within the observed low-resolution pixel L x/2,y/2 . Other shared values may span multiple lowresolution pixels. With this notation, there are (N x + 1) × (N y + 1) spatially shared values S. These located at the image boundaries, or which would span invalid pixels such as in the case of sensor failure, are simply not commonalized and effectively remain internal to the valid pixels. All shared values are expressed in reflectance units and constrained to the range of their respective band.
In remote sensing, the reflectance of each pixel is often considered to be a linear mix [20] of the reflectances of its constituents (e.g. a mix of vegetation and soil). Assuming the shared values correspond to unknown constituents spanning over pixel boundaries, and using this linear mixing model, the proportion of each shared value that is present in each pixel is thus determined by weights that are specific to that pixel (Fig. 2, right). This leads to the following mixing equation for the shared values and the weights: With the following constraints:

III. SOLVING THE SUPER-RESOLUTION PROBLEM
A. Separating band-specific information from information common to all bands The proportion of mixed elements within a pixel (e.g. 20% road / 80% vegetation) is a physical property of that pixel, but the reflectance of each element depends on the spectral band at which it is observed. Therefore, the weights are common to all bands, while shared values are band-dependant. Weights encode the geometric consistency of pixels accross bands. Shared values encode the spatial consistency of nearby pixels. The high-resolution data are used to fit the full mixing model, containing both weights and shared values. This step is presented below. The next section addresses how to un-mix low-resolution bands in order to produce the super-resolution result, reusing the weights fit from the high-resolution bands.
Starting from an observed high-resolution band H o , a downsampled version L d of the data is created with Eq. 1. The best mixing model is estimated by minimizing the difference between: a) the observed pixel values H o , and b) the resolution-enhanced values H r computed from the downsampled data L d . Let us subscript data specific to each band with an additional index β. Thus, L d , H o , H r and S are subscripted, but not the weights W . Solving this first problem is a constrained minimization, for k = 0 . . . 3, and β ∈ H the set of high-resolution bands: with each H r β,x,y term given by Eq.
2. An iterative solver [21] is used for constrained least squared error optimization 1 , allowing Eq. 3 to be enforced by a reparametrization and Eq. 4 by soft boundaries (a reference implementation is provided, link given in appendix). Initial weights for the iterative algorithm are set to 1 4 (i.e. equal influence to all shared values, see Fig. 2, left). The initial shared values S ini β are computed by averaging each high-resolution pixel H o β that 1 The Ceres solver [21] can be fine-tuned with many internal parameters. Extensive testing determined that conjugate gradients with a block Jacobi preconditionner give the best quality/processing time tradeoff for the superresolution problem presented in this article. These are set by default in the reference implementation.
partially covers S ini β in Fig. 2. For example, S ini β,x+2,y+1 = 1 2 H o β,x/2,y/2 + H o β,x/2+1,y/2 for x < N x − 1 and y < N y . At the end of this step, both shared values S opt β between high-resolution pixels, and weights W opt common to all bands, are available.

B. Estimating shared values from low-resolution data
Shared values S opt are found by optimization on highresolution data, so they cannot be estimated directly on the low-resolution bands with the above procedure. Instead, the relation between S opt and nearby low-resolution pixels can be learned from downsampled high-resolution bands L d . That relation is also expressed as a geometric property common to all bands, so it can be used in order to produce a first estimate S f it for the low-resolution bands. More specifically, a second set of mixing coefficients V is built in order to fit S opt β,x,y from low-resolution pixels L d β,n at nearby locations n ∈ N (x, y). See Fig. 3, with N being either the corner, middle or inner variant depending on the position of (x, y) with respect to the low-resolution reference pixel. These neighborhoods hopefully capture local features (e.g. edges), in the form of up to 9 coefficents v x,y,n for each highresolution pixel (x, y) in the image: With this global optimization, the set of coefficients V encodes how the shared values are related to their low-resolution neighborhood, independently of the spectral band. They are fit from the high-resolution bands, and then propagated to the low-resolution bands b ∈ L in order to provide a first estimate S f it for the shared values in each band b: The fit from Eq. 6 is rarely perfect and values S f it β,x,y = n∈N (x,y) v n L d β,n can also be computed for the original highresolution bands. The ratios q β,x,y = S opt β,x,y /S f it β,x,y are then exploited in order to mimick the panchromatic sharpening method in [3], but using the multiple high-resolution bands instead. For Sentinel-2, no panchromatic band is available to encompass the spectrum of all low-resolution bands. An idea is to empirically replace the panchromatic band by a combination of high-resolution bands. Problems mentionned in Section I-B are adressed by weighting bands that yield close spectral responses for the shared values. For each low resolution band b ∈ L, and for each shared value, a normalized proximity measure is defined as x,y |, where the normalization is performed by using the maximum discrepancy over all high resolution bands α ∈ H. Combining the high-resolution sharpening ratios is then performed by geometric averaging, using these proximity measures as weighting factors: This overall average sharpening ratio is used as a prefactor for setting corrected shared values S cor b,x,y =q b,x,y S f it b,x,y for the low resolution bands. Having now estimated high-resolution shared values for the low-resolution bands, it is a simple matter of combining these S cor b,x,y with the weights W by applying Eq. 2, in order to produce super-resolved pixels. A final rescaling of these H r b,x,y is performed so as to ensure reflectance is preserved (Eq. 1).

C. Super-resolving 60m/pixel bands
In this setup, each low-resolution pixel corresponds to 36 values at 10m/pixel. Solving this directly is not tractable, but an indirect solution with an intermediate step is feasible: -In a first pass, a 60m/pixel band is super-resolved to 20m/pixel. There are then 9 sub-pixels to infer for each lowresolution pixel, with 8 free parameters ( Fig. 4 and Eq. 9 below). However, there are also 10 bands at 20m/pixel (including the 4 bands at 10m/pixel, downsampled). These provide enough constraints for the inference of weight values, common to all these bands, computed with a modified version of the method presented below.
-In a second pass, the 20m/pixel solution from the first pass is super-resolved down to 10m/pixel with the method described in the previous sections.
Adapting the above notation to the 9 sub-pixels problem, in this section the low resolution bands are b ∈ L at 60m/pixel, while the high resolution bands β ∈ H consist of the 20m/pixel bands (either original or down-sampled from 10m/pixel). Preserving the reflectance imposes (Fig. 4): H r b is the intermediate super-resolution solution at 20m/pixel for band b. Fig. 4 also shows the relation between L b , the high-resolution pixels H r b and the shared values. With this convention, the weights W follow exactly the same (x, y) pattern with respect to S and H as in Fig. 2 Figure 4. Introducing the indexing and the relations between the low and high resolution pixels, for the case of tripling the resolution. The shared values between pixels are also indicated. Weights are still internal to each high-resolution pixel with exactly the same structure as in Fig. 2, right. available high-resolution bands, is thus also the same as above. All 10 bands β ∈ H are used as constraints for Eq. 5. A difference lies in estimating V from nearby pixels. There are still four neighborhood patterns of the "corner" type (see Fig. 3), for shared valuesS x,y , S x+3,y , S x,y+3 and S x+3,y+3 . But there are now two "middle" neighborhood patterns for each side of the lower-resolution pixel (e.g. S x+1,y , S x+2,y ), as well as four "inner" neighborhood patterns instead of one (see Fig. 4, right). With these definitions, Eq. 6 is solved as before. Averaging the sharpening ratios now involves all 10 bands β ∈ H (instead of 4 in the previous section), but this does not change Eq. 8. Thus, solving the 60m→20m superresolution problem is performed with very little adaptation.
Once computed, the intermediate solutions at 20m/pixel are further processed by a final 20m→10m super-resolution step, as described in the above sections.
-The normalized mean squared error ERGAS(x, y) = 100R 1 x is the reference image, and y is the image to be tested.
-The spectral angle SAM = arccos x·y x y , considering the x and y images as vectors. SAM is given in degrees in the next Section.
-When a panchromatic band P is available, the quality with no reference where D λ is a spectral distortion index and D s is a spatial distortion index. But D s depends on P so, in the present case with no panchromatic band available, QN R is not useable. In the following sections, Q, ERGAS and SAM are computed. The typical methodology for using these quantifiers is to down-sample an image, super-resolve this down-sampled image back to the original resolution, and compare with the original data. Typically, the highest-resolution bands are used for downsampling/super-resolution. But, due to the way superresolution is performed in this paper, downsampling 10m bands to 20m and super-resolving them back for comparison is not acceptable. Indeed, that very step is already included as part of the optimization in Eq. 5. Moreover, shared details coming from all original 10m bands are taken into account in Eqs. 6-8. Hence, a test that uses these same 10m bands as a basis for comparison wouldn't be fair. The method is thus adapted as follows : -The 20m bands are downsampled to 40m and the four 10m bands are downsampled to 20m.
-Each downsampled 40m band is super-resolved back to 20m using the four original 10m bands that were downsampled at 20m.
-Each original 20m band x is compared with its downsampled/super-resolved version y This way, the original 10m bands are only used in order to build y, but are not themselves the basis of comparison x. The downside is that this method does not quantify directly the quality of final product (the super-resolution of all bands to 10m), but a proxy for it.
Statistics are given for each 20m band, then globally averaged over all bands : geometric average for Q, using the given formula for ERGAS, and arithmetic angle average for SAM .

V. RESULTS
Three use cases were selected for testing the algorithm: -A coastal environment, the delta of the Eyre river (France); -A urban area, the city of Bordeaux; -Fields, in the Bordeaux peri-urban area. All these regions of interest are tested on an image acquired by the Sentinel 2A satellite on 2016/08/22 and processed with the "sen2cor" atmospheric correction utility [24]. Figs. 5,8 and 11 show the three selected regions as a composite images from the 10m/pixel visible bands, where each blue, green, red component was scaled between 1% and 99% of the original reflectance.
As an indication of computational performance, processing all the bands in either of these 408x300 pixel areas takes on average 1min17s on a machine with twelve 1.9GHz cores.

A. Coastal Area
Results for the coastal environment are displayed in Figs

VI. DISCUSSION
Details not present in the original bands are immediately visible in all images, especially at the largest super-resolution 60m->10m (e.g. Fig 5, middle-right). These details correspond to the band-independant information that was extracted from the other bands, and propagated to these images. Although each band presents different reflectance properties (in particular, B1 and B9), the exact same weights and geometric information extracted from Eqs. 5,6,8 were applied to all super-resolved images. This example demonstrates how the method correctly extracts band-independant information that encodes image details, while preserving the reflectance of each band (Eqs. 1,9).
Quantitative indicators are given for three typical land cover types. The method works best for agricultural environments, with large uniform areas, and worst in urban environments. Even then, Fig. 10, right, shows details that are very well recovered. Comparison with quantitative indicators from other works using panchromatic sharpening should be taken with caution. We estimate only the 40m->20m super-resolution proxy for reasons mentionned in Section IV. Nevertheless, compared to the litterature [22], [23], [6], our global averages of 0.941 ≤ Q ≤ 0.994, 2.08 ≤ ERGAS ≤ 4.49 and 2.21 ≤ SAM ≤ 4.87 in the different regions of interest for Sentinel-2 data are consistent with the state of art for panchromatic sharpening, albeit without a panchromatic band.
The method presented here includes multiple steps that are not trivial, but they are all necessary. In order to demonstrate this, a final experiment is performed on the coastal area, enlarged to encompass nearby cities and fields (Fig. 14). The 60m band 1 is used in order to best visualize the effect of each step. A fist idea would be to apply ratio sharpening (Eq. 8) directly on data values, instead of spatially shared values, so as to simulate panchromatic sharpening. This would simplify the method drastically. The result of this experiment is shown in Fig. 14, top-right, to be compared with the correct superresolved result in Fig. 5, middle-right. The role of spatial consistenty is immediately highlighted: without the shared values trick, unacceptable pixel blocks are clearly visible in the result image. Conversely, why isq ratio sharpening useful? Fig. 14, bottom-left shows that it is in fact quite important for recovering the fine structures. Given that importance, one may then question the usefulness of extracting weights W as bandindependant information, especially since we also compute reverse weights V in a second step. Why would these W encode image details? Fig. 14, bottom-right shows the result of simply setting these weights to 1 4 and S opt to the average values, as described in Section III-A, while maintainingq ratio sharpening. As expected, details are also smoothed out. Weights W are defined between high-resolution pixels, hence encode high-resolution details. The reverse weights V encode larger range patterns present in surrounding low-resolution pixels. Results presented in this section use the 60m/pixel band 1, for which two super-resolution steps are applied. This choice was made so as to enhance and clearly highlight the influence of shared values S, ofq ratio sharpening, and of weights W , V , on the final result. For 20m/pixel bands only one super-resolution step is applied, but all parts of the method are still needed for good results.

CONCLUSION
This article presents a super-resolution method based on exploiting both the local consistency between neighborhood pixels and the geometric consistency of sub-pixel constituents across multispectral bands. Figs. 5-13 show the result of applying this method to Sentinel-2 images, in order to bring all bands from 20m/pixel and 60m/pixel down to the highest resolution at 10m/pixel. The algorithm is however generic and could be applied to other multi-resolution and multispectral satellite images. Further work could include the usage of secondary images, taken from a satellite with low temporal resolution but with a higher spatial resolution, in order to extract the geometric information used for the super-resolution. Assuming the pixel geometry does not change much between these acquisitions, then Sentinel-2 images could be enhanced below 10m/pixel. This form of multi-satellite temporal superresolution would combine high temporal frequency with high spatial resolution. Another trail of research would be to incorporate the super-resolution algorithm directly within the atmospheric correction step [24], rather than applying it as a separate stage. Indeed, using higher-resolution pixels instead of low-resolution bands for calibrating the atmospheric correction may lead to better accuracy.

ACKNOWLEDGMENTS
The author thanks Hussein Yahia and Dharmendra Singh for useful discussions and feedback while carrying this work.

SOURCE CODE
The source code for this super-resolution algorithm is provided as Free/Libre software, under either (your choice) the lesser GNU public licence v2.1 or more recent, or the CeCILL-C licence. The library is written and usable directly in C++ and it is also wrapped in a Python package. A ready to use Python script for super-resolving Sentinel-2 images is provided. Band 10 is dedicated for cloud detection and removed by the bottom-of-atmosphere correction utility [24]. The algorithm presented in this paper nevertheless takes it into account when applied to the unprecessed, top-of-atmosphere images.          are not used, that demonstrates their role in maintaining spatial consistency. Bottom-left: Incorrect image obtained when only mean square optimisation is applied, that shows the importance of theq ratio sharpening step described in Section III-B. Bottom-right: Incorrect image obtained when the estimation of weights W is omitted while maintaning theq ratio sharpening step.