IEEE Geoscience and Remote Sensing Magazine - March 2021 - 59

to restore the unbiased sharpened image. Thus, a hazecorrected version of the BT in (4), called BT-H, has recently been proposed [59]. Equation (4) is modif ied
as follows:

	

R
MS k - H
\
MS k = R
MS k + I - H k · ^P - I Lh
L
p

R
MS k - H k
^
= I - H · P - H ph + H k,
L

p

(5)

where H k denotes the constant haze for the kth MS band
and H p is the haze both of the intensity component and
of PAN [59]. For further details about this approach and
the estimation of the haze values, interested readers can
refer to the original article [59] and the related literature
[56], [60].
GRAM-SCHMIDT
The GS orthogonalization procedure has been used
for pansharpening [24], taking the name GS spectral
sharpening. This procedure exploits the intensity component, I L , as the first vector of the new orthogonal
basis. The orthogonalization processes one MS vector
at a time, finding its projection on the (hyper)plane
defined by the previously found orthogonal vectors
and its orthogonal component such that the sum of
the orthogonal and projection components is equal to
the zero-mean version of the original vectorized band.
Pansharpening is completed by substituting I L with the
histogram-matched P before performing the inverse
transformation.
The procedure can be recast into the flowchart in
Figure 1. Thus, the fusion process is described by (1) with
the following injection gains for k = 1, f, N [30]:
	

Gk = gk 1 =

cov (R
MS k, I L)
var (I L) 1, (6)

in which 1 is an all-one matrix; cov ^ X, Y h indicates the covariance between two images, X and Y; and var ^ X h is the
variance of X .
Several versions of GS are achieved by changing the
way to generate I L : by averaging the MS components
(i.e., w i = 1/N , for all i = 1, f, N) in the original GS [24];
by a weighted average of the MS bands minimizing the
MSE with respect to a low-pass-filtered version of a PAN
image in the GSA [30] and by applying GSA separately
to each cluster obtained through k -means clustering applied to the MS image in the context-adaptive GSA (CGSA) [32].
BAND-DEPENDENT SPATIAL-DETAIL-BASED
APPROACHES
The BDSD-based approaches [33] start from an extended
version of the generic formulation (1) for all k = 1, f, N .
Thus, we have
MARCH 2021

IEEE GEOSCIENCE AND REMOTE SENSING MAGAZINE

	

\
MS k = R
MS k + g k d P - / w k, i R
MS i n . (7)
N

i=1

The optimal minimum MSE joint estimation of the
weights, w k, i , and scalar gains, g k , in (7) would encompass
the use of the unknown target image, R
MS k . Thus, the problem is solved at RR by exploiting the original MS image as
a reference image [33]. A context-adaptive BDSD (C-BDSD)
has also been proposed in [34]. A physically constrained
optimization (BDSD-PC), improving the quality of the fusion product, particularly for eight-band MS data sets, has
been recently proposed in [52].
PARTIAL REPLACEMENT ADAPTIVE COMPONENT
SUBSTITUTION
The concept of partial replacement of the intensity component is introduced in [29]. The PAN image is not directly
used for CS. Instead, the algorithm uses a weighted sum of
the PAN and the kth MS band to calculate the kth sharpened band in (1). For this reason, this method is referred to
as PRACS.
MRA
MRA methods extract the PAN details exploiting the difference between P and its low-pass spatial version, PL .
Namely, the fused image is obtained for each k = 1, f, N
as follows:
	

\
MS k = R
MS k + G k·^P - PLh . (8)

The different approaches belonging to this class are characterized by the way to calculate PL and by the employed
injection gains " G k ,k = 1, f, N . In a very general setting, PL
is achieved through an iterative decomposition scheme
(called MRA) that aims at constructing a sequence of 2D
signals with progressively reduced information through the
repeated application of some analysis operators.
The general fusion scheme for MRA approaches is displayed in Figure 2. We can note blocks devoted to the interpolation of the MS image to reach the PAN scale; the
calculation of the low-pass version PL of the PAN image
depending on the resolution ratio, R; the computation of
the band dependent injection gains, " G k ,k = 1, f, N; and the
fusion according to (8). It is worth noting that there is a
difference if PL is decimated-interpolated (dashed box) or
not. The first solution is usually much more appreciated
due to the possibility of compensating for the aliasing of
the MS image through the fusion process [127]. Finally,
MRA-based fused products reduce the spectral distortion
but often compromise the spatial rendering. They are also
robust to temporal misalignments, which is a valuable feature when multiplatform fusion is considered [128].
ADDITIVE WAVELET LUMINANCE PROPORTIONAL
The undecimated " à trous " method [44] has recently
emerged as a very effective one for extracting PAN details
59
IEEE Geoscience and Remote Sensing Magazine - March 2021

Table of Contents for the Digital Edition of IEEE Geoscience and Remote Sensing Magazine - March 2021
