# Compression of Multispectral Images with Comparatively Few Bands Using Posttransform Tucker Decomposition

Authors: Jin Li; Fei Xing; Zheng You
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Up to now, data compression for the multispectral charge-coupled device (CCD) images with comparatively few bands (MSCFBs) is done independently on each multispectral channel. This compression codec is called a “monospectral compressor.” The monospectral compressor does not have a removing spectral redundancy stage. To fill this gap, we propose an efficient compression approach for MSCFBs. In our approach, the one dimensional discrete cosine transform (1D-DCT) is performed on spectral dimension to exploit the spectral information, and the posttransform (PT) in 2D-DWT domain is performed on each spectral band to exploit the spatial information. A deep coupling approach between the PT and Tucker decomposition (TD) is proposed to remove residual spectral redundancy between bands and residual spatial redundancy of each band. Experimental results on multispectral CCD camera data set show that the proposed compression algorithm can obtain a better compression performance and significantly outperforms the traditional compression algorithm-based TD in 2D-DWT and 3D-DCT domain.


---

## Body

## 1. Introduction

Multispectral CCD images can be used to identify the targets using their geometric feature and spectrum information simultaneously. Multispectral CCD imaging technology has extensive application in ground target detection, geology monitoring, climate monitoring, military reconnaissance, scientific research, and so on [1, 2]. With the increase of application requirements, the performance indexes of multispectral CCD camera, such as field of view, coverage width, spatial resolution, spectral resolution, time resolution, and radiation resolution, continue improving. These lead to rapid increase of amount of multispectral image data. However, the amounts of memory and the channel bandwidth used on the satellite are limited. Therefore, the suitable compression approaches play an important role in the successful applications of multispectral CCD camera used on the satellite [3–5].Multispectral CCD camera has two types of instruments: (1) comparatively many spectral channels and (2) comparatively few spectral channels [6, 7]. Therefore, multispectral CCD images also have types: (1) multispectral CCD images with comparatively few bands (MSCFBs) and (2) multispectral CCD images with comparatively large number of bands (MSCLBs). MSCFBs and MSCLBs both have spectral redundancy and spatial redundancy. For MSCLBs, the correlation between the different bands is usually stronger than spatial correlation among pixels in the same band [8]. There have existed several well-known MSCLBs compression algorithms, like principal components analysis (PCA) [9]. The PCA approach combines all bands into one single spectral matrix. The covariance matrix of the spectral matrix can be obtained. The inverse of eigenvector matrix of covariance matrix multiplies by the spectral matrix to obtain new bands called principal components. The number of principal components is equal to the number of bands. There is no correlation between different components. Therefore, the PCA is mainly used to remove spectral redundancy. In order to remove the spectral and spatial correlation simultaneously, the MSCLBs compression algorithms use three dimensions (3D) transform approach, like 3D-DCT [10], 3D-DWT [11, 12], and Karhunen-Loeve transform (KLT) + 2D transform [13]. It is generally accepted that the 3D-DWT-based approach is one of best compression algorithms for MSCLBs. After each dimension of MSCLBs performs 1D-DWT, the MSCLBs are decomposed into several different resolution 3D subimages. The 3D subimages have pyramid and multiresolution characteristics. These characteristics are very suitable for the latter entropy encoding and can complete progressive encoding.The compression algorithms based on 3D-DWT, such as 3D-EZW, 3D-SPIHT [14], 3D-SPECK [15, 16], and 3D-EBCOT, have been widely used to compress the MSCLBs. In [17], Tang and Pearlman show that the 3D-SPECK approach can improve 0.061~0.173 dB more than 3D-SPIHT at 0.2~2.0 bpppb. It is generally accepted that JPEG2000 algorithm [18] is the best algorithm for 2D still images. Recently, JPEG2000 algorithm has extended to MSCLBs. JPEG2000 algorithm in 3D model [19] uses 3D-DWT combined with 3D-EBCOT to compress MSCLBs. JPEG2000 algorithm is usually considered as a reference standard to evaluate other algorithms. In [20], Du and Fowler proposed a compression approach based on the JPEG2000 combined with PCA. The proposed algorithm can get better rate distortion (RD) performance than JPEG2000 combined with DWT for removing spectral correlation. Up to now, the only compression standard for MSCFBs and MSCLBs is CCSDS-MHDC. The CCSDS-MHDC [21] algorithm mainly targets space multispectral images data compression. The CCSDS-MHDC algorithm uses FL-based approach to trade off between computation complexity and compression performance. The CCSDS-MHDC algorithm is suitable for the application of on-board multispectral camera. However, the CCSDS-MHDC is only suitable for lossless compression. The CCSDS-MHDC can only provide limited compression ratios which are not a reasonable value in high-resolution remote sensing multispectral camera.Recently, a tensor decomposition (TD) [22] method is applied to the hyperspectral images coding. Hyperspectral images are considered as a third-order nonnegative tensor. And then the nonnegative tensor is decomposed by Tucker method to reach the compression results. In [23], Karami et al. proposed a new hyperspectral image compression method based on 3D-DCT combined with TD. Their core idea is that the TD is used to decompose 3D-DCT coefficients tensor to remove spectral and spatial redundancies. In [24], Karami et al. proposed a new hyperspectral image compression method based on 2D-DWT combined with TD. 2D-DWT performs on each band to remove spatial correlation. The TD is used to decompose wavelet subband tensor to remove spectral correlation. These methods are also suitable for MSCLBs. However, they are too complex to realize in on-board multispectral camera.For MSCFBs, the comparatively few bands, like three or four bands, are captured by on-board CCD camera. These proposed MSCLBs compression algorithms are not suitable for MSCFBs because they have comparatively few bands. Up to now, compression methods for MSCFBs are done independently on each multispectral channel. This compression codec is called a “monospectral compressor.” The monospectral compressor does not remove spectral redundancy stage. In order to reach the excellent compression performance, MSCFBs compression algorithm must remove both spatial and spectral redundancies. In this paper, we propose an efficient compression approach for MSCFBs. In our approach, the 1D-DCT is performed on spectral dimension to exploit the spectral information, and the posttransform (PT) in 2D-DWT domain is performed on each spectral band to remove the spatial redundancies. MSCFBs are decomposed into several subtensors by PT. Then, these subtensors are decomposed by multilevel Tucker decomposition (TD) to reach compression results.


## 2. Proposed Algorithm

### 2.1. The Imaging Principle of MSCFBs

Figure1 shows the imaging principle of MSCFBs. The long linear multispectral CCD has four parallel CCD linear arrays. Each array is one imaging spectral band. Four imaging spectral bands are red, blue, green, and near-infrared, respectively. The light radiated or reflected by the hundreds of kilometers of linear arrays of ground pixels is concentrated onto optical thin film of CCD detector through an optical system. The space and spectrum distributing of ground targets radiation, acquired by CCD detector, can be expressed as
(1)
B
(
x
,
y
,
θ
,
λ
)
=
1
π
[
E
0
(
θ
,
λ
)
ρ
(
x
,
y
,
λ
)
τ
α
(
λ
)
*
h
α
(
x
,
y
,
λ
)
]
+
N
(
θ
,
λ
)
,
where θ is the solar height angle, λ is wavelength, E
0 is ground luminometer, ρ is spectral reflectance, τ
α is atmospheric transmittances, h
α is the point spread function, and N is radiant flux. The optical thin film of each linear array CCD allows the light of corresponding wavelength through. B (x
,
y
,
θ
,
λ) of each band is captured by a linear CCD array, and one dimensional spatial and spectral information is gained. When the CCD camera scans the ground target, spatial dimension information is gained. Therefore, MSCFBs are considered as 3D images, Y
∈
R
I
1
×
I
2
×
I
3, which is shown in Figure 2.Figure 1
Imaging principle of MSCFBs.Figure 2
MSCFBs properties.Based on the imaging principle of MSCFBs, MSCFBs have the spectral redundancies between the adjacent spectral bands and spatial redundancies between the adjacent pixels in each spectral band. Therefore, the compression algorithm must remove the spectral redundancies and spatial redundancies simultaneously.


### 2.2. Spectral Decorrelation

MSCFBs have two types of correlations: spatial and spectral correlations. The spatial correlation is defined as the correlation between the adjacent pixels in each spectral band. The spatial line correlation of MSCFBs can be expressed as(2)
R
(
p
,
0
)
=
∑
x
=
1
M
∑
y
=
1
N
(
f
(
x
,
y
)
-
f
¯
)
×
(
f
(
x
+
p
,
y
)
-
f
¯
)
∑
x
=
1
M
∑
y
=
1
N
(
f
(
x
,
y
)
-
f
¯
)
2
.The spatial column correlation of MSCFBs can be expressed as(3)
R
(
0
,
q
)
=
∑
x
=
1
M
∑
y
=
1
N
(
f
(
x
,
y
)
-
f
¯
)
×
(
f
(
x
,
y
+
q
)
-
f
¯
)
∑
x
=
1
M
∑
y
=
1
N
(
f
(
x
,
y
)
-
f
¯
)
2
,
where M is the line number of each band, N is the pixel number of each line, f
(
x
,
y
) is the function of the band gray, and f
¯ is the average value of the band gray. The spectral correlation is defined as the correlation between the pixels located in the same spatial position in the different spectral band, which can be expressed as(4)
h
(
p
,
q
)
=
(
∑
x
=
1
M
∑
y
=
1
N
(
f
(
x
+
p
,
y
+
q
)
-
f
¯
)
00000000
×
(
g
(
x
,
y
)
-
g
¯
)
∑
y
=
1
N
)
×
(
(
∑
x
=
1
M
∑
y
=
1
N
[
f
(
x
,
y
)
-
f
¯
]
2
)
0000
×
(
∑
x
=
1
M
∑
y
=
1
N
[
g
(
x
,
y
)
-
g
¯
]
2
)
)
-
1
/
2
,
where f
(
x
,
y
) and g
(
x
,
y
) are the value of pixel located in (
x
,
y
) in two adjacent bands and h
(
0,0
) is the cross-correlation function of f
(
x
,
y
) and g
(
x
,
y
). Based on (2)~(4), the Landsat-7 Thematic Mapper (TM) images are chosen to perform an experiment of the spatial and spectral correlation. Figure 3 shows the test results.The spatial and spectral correlation test result.
(a)
Spatial correlation test result
(b)
Spectral correlation test resultThe correlations between the adjacent bands for an MSCFB are in the range[
0.318
–
0.926
], and the line and column correlations between the adjacent pixels for an MSCFB are in the ranges [
0.685
–
0.955
] and [
0.609
–
0.963
], respectively. Therefore, the MSCFBs have higher spatial correlation than spectral correlation. Because the spectral correlations are comparatively weak, the traditional 3D-DWT-based compression algorithms are not suitable for MSCFBs. When the 3D-DWT-based approaches compress Landsat-7 images, the good compression results are not obtained. Simultaneously removing the spectral redundancies and spatial redundancies is primary target for the successful implementation of a high-performance compression method.Because the MSCFBs have a smaller number of bands, 1D-DCT is performed on the spectral dimension to decorrelate the spectral bands. The number of bands is denoted asP; the size of each band is denoted as M
×
N. An MSCFB is composed of M
×
N 1D spectral vector (Length P). The definition of the 1D-DCT for an input 1D spectral vector A
i (i
=
1,2
,
…
,
M
×
N) and output image B
i  (i
=
1,2
,
…
,
M
×
N) is
(5)
B
i
r
=
α
r
∑
k
=
0
P
-
1
A
i
k
cos
⁡
(
π
(
2
k
+
1
)
r
2
P
)
,
α
r
=
{
1
P
,
r
=
0
,
2
P
,
1
≤
r
≤
P
-
1
.
The B
i
r values are called the DCT coefficients of A
i. Equation (5) is performed M
×
N times to complete the spectral transform. In order to speed up calculation, we proposed an iterative method. The iterative equation can be expressed as
(6)
B
i
+
1
r
=
∑
k
=
0
P
-
1
A
i
+
1
k
∑
k
=
0
P
-
1
A
i
k
×
B
i
r
.
The DCT coefficients of the current spectral vector can be computed by the intermediate and output data obtained from the previous spectral vector.


### 2.3. Spatial Decorrelation

The 2D-DWT can decompose the image into multilevel wavelet to obtain lower resolution subband and detail subbands. The multilevel decomposition can be regarded as the process of low pass and high pass filtering. At each level, the low subband produced by previous level performs the high pass filtering to obtain detail information called wavelet coefficients, while performing the low pass filtering linked with scaling function to obtain the approximate information called scaling coefficients. The DWT has been widely employed exploiting the spatial correlations for remote sensing image, such as JPEG2000 and CCSDS-IDC. In this paper, we apply a 2D-DWT coupled with a posttransform to each band of MSCFBs. In our approach, the 2D-DWT is performed on each image band to reduce spatial correlations, and then the transformed wavelet coefficients are transformed by the posttransform approach to reduce remaining correlations in wavelet subbands.After remote sensing MSCFBs are transformed by 2D-DWT, there exits residual correlation between coefficients in a small neighborhood (see Figure4). Statistical dependencies between DWT coefficients have been studied for many years. In [25], correlations between nearby wavelet coefficients are reported to be in the range [
0.01
–
0.54
] at a distance less than or equal to 3 pixels. We found even wider range, and here we provide a more detailed discussion regarding this topic. We use a Pearson correlation coefficient ρ to analyze statistical dependency of DWT coefficients, which can be expressed as
(7)
ρ
x
,
y
=
δ
x
y
(
δ
x
δ
y
)
,
where δ
x
,
y denotes the covariance between the variables X and Y and δ
x and δ
y denote the standard deviation of X and Y, respectively. According to project experience, three-level 2D-DWT is appropriate for on-board compressor and we used three-level 2D-DWT in this paper.Figure 4
The wavelet transform of high-resolution image with a zoom on some coefficients.Three-level 2D-DWT is performed on each image band to produce one low-frequency subband (denoted as LL) and nine high-frequency subbands (denoted asLH
⁡
L, HL
L, HH
L, L
=
1,2
,
3). The test is performed from ten MSCFBs. The residual correlations between wavelet coefficients in a small neighborhood of wavelet subbands are shown as in Figure 5. In Level L (L
=
1), the residual directional correlation within 16-connected region is ρ (HL1, HH1, LH1, |
Δ
x
|
≤
2, |
Δ
y
|
≤
2) > 0.3. In Level L (L
=
2), the residual directional correlation within 16-connected region is ρ (HL2, HH2, LH2, |
Δ
x
|
≤
2, |
Δ
y
|
≤
2) > 0.4. In Level L (L
=
3), the residual directional correlation within 16-connected region is ρ (HL3, HH3, LH3, |
Δ
x
|
≤
2, |
Δ
y
|
≤
2) > 0.3. These imply a strong redundancy among neighboring 4 coefficients. The results in Figure 5 also indicate that the residual correlations between nearby wavelet coefficients are much weaker at a distance equal to 5 pixels.The correlation between the adjacent wavelet coefficients.
(a)
L
=
1
(b)
L
=
2
(c)
L
=
3To obtain the optimum compression performance for MSCFBs, the compression algorithms must fully consider the abovementioned statistical properties. In [26, 27], EBCOT has been reported to be very efficient method to remove these residual redundancy values. However, its implementation complexity is too high. In [28], Delaunay has proposed a posttransform to remove remaining redundancies between the adjacent wavelet coefficients within a small region. After the wavelet transform of an image, wavelet coefficients or scaling coefficients subbands are divided into several coefficients blocks. Each block has 4 × 4 wavelet coefficients. Posttransform performs on each coefficients block to remove the residual redundancies. From testing a large number of remote sensing images, the selected size of wavelet coefficients block is the best to a low complexity and efficient compression method. There are two major reasons: (1) the correlation between the adjacent wavelet coefficients is very weak when a distance is greater than or equal to 5 pixels; (2) the larger size of wavelet coefficients block leads to the higher computational complexity, and the smaller size of block results in increasing the number of coefficients blocks and the size information. In addition, since the posttransform is performed on coefficients blocks in wavelet transform domain, not on image domain, the compression algorithm has no additional overhead time to perform deblocking filter.The core idea behind posttransform compression is that wavelet coefficients blocks further are transformed using one group particular direction basis (such as Bandelet, DWT, DCT, and PCA) in a dictionary. First, a 2D-DWT is applied to an image. Next, blocks of 4 × 4 DWT coefficients are projected onN
B orthonormal bases B
b of the dictionary D. Then, the minimizing of the Lagrangian cost function is calculated to select the best posttransformed block. Finally, the posttransformed coefficients are encoded by entropy coding method. Each 4 × 4 DWT coefficients block f is considered as an input vector f
=
[
f
1
,
f
2
,
…
,
f
M
]
T
∈
R
M with M
=
16. The M vectors of the basis B
b are noted ϕ
m
b with m
∈
[
1
,
M
-
1
]. The posttransformed block f
b can be expressed as follows:
(8)
f
b
=
∑
m
=
1
16
〈
f
,
ϕ
m
b
〉
·
ϕ
m
b
.
Since a dictionary has N
B bases, N
B
+
1 (including one original block) posttransformed blocks f
b can be obtained. Among all the posttransformed blocks f
b, the best posttransformed block f
b
* is selected according to calculating the minimizing of the Lagrangian rate-distortion cost. The minimizing of the Lagrangian rate-distortion cost can be expressed as
(9)
L
(
f
q
b
)
=
D
(
f
q
b
)
+
λ
·
R
(
f
q
b
)
,
where f
q
b denotes the quantized posttransformed coefficients, q is the quantization step, and D
(
f
q
b
) denotes the square error between the posttransformed coefficients f
b and the quantized posttransformed coefficients f
q
b. λ is a Lagrangian multiplier and R
(
f
q
b
) denotes the required bit rates for encoding f
q
b and the associated side information b
*.In the posttransform, the dictionary has multiple bases. The better the compression performance is, the more the number of bases is and the higher the computational complexity is. However, on-board compression requires a low computational complexity. The space MSCFBs compressor allows the dictionaryD of posttransform to have only one basis. In [25], Delaunay et al. proposed a simple posttransform-based compression method. They only perform one-level wavelet transform and only use a Hadamard basis. Then, the posttransform only performs on the coefficients subbands. Compared with the compression method only using the DWT, the PSNR can only improve 0.4 dB~0.6 dB. In this paper, to obtain a posttransform with the low computational complexity and still have high compression performance, we thus consider a very simple dictionary containing only one dynamic base, which is Hadamard basis at the low bit rates, and DCT basis at the high bit rates.In [29], Delaunay has shown that the Lagrangian approach for selecting the best posttransformed block has two main drawbacks. First, the estimation of the required bit rate for R
(
f
q
b
) is computationally intensive and not always accurate. Second, calculating the Lagrangian function depends on the quantization step q. However, the coder does not define the quantization step when coding. Therefore, Delaunay proposed an l
1-norm minimization approach to select the best posttransformed block. In this paper, we propose a new method based on l
p-norm minimization approach to evaluate best posttransform, which replaces an l
1-norm minimization. It is very suitable for the space MSCFBs compressor which has the low complexity constraints.At the low compression bit rate, the bit rateR can be expressed as
(10)
R
≈
γ
0
M
with
γ
0
=
7
,
where M is the number of nonzero posttransformed coefficients. We propose an l
0-norm (p
=
0) minimization approach to select the best posttransformed block. The selected best posttransformed block has the minimal sum of coefficients magnitude, which can be expressed as
(11)
f
b
*
=
arg
min
⁡
f
b
,
b
∈
[
0
,
N
B
]
∥
f
b
∥
0
,
∥
f
b
∥
0
=
∑
m
=
1
M
I
m
with
I
m
=
{
1
,
⌊
|
a
b
[
m
]
|
⌋
≠
0
,
0
,
⌊
|
a
b
[
m
]
|
⌋
=
0
,
where a
b
[
m
] is a posttransformed coefficient. At the high compression bit rate, the bit rate R can be expressed as
(12)
R
≈
-
∑
i
=
1
N
p
i
lo
g
2
p
i
-
lo
g
2
q
.
We use the l
1-norm (p
=
1) minimization approach to select the best posttransformed block. The selected best posttransformed block has the minimal sum of coefficients magnitude, which can be expressed as
(13)
f
b
*
=
arg
min
⁡
f
b
,
b
∈
[
0
,
N
b
]
∥
f
b
∥
1
s
.
t
.
∥
f
b
∥
1
=
(
∑
m
=
0
15
|
a
b
[
m
]
|
)
.Figure6 shows the proposed posttransform architecture. Each high-frequency subband is performed by posttransform. The bit-rate comparator decides the type of coding bit rate. We define the coding bit rate greater than or equal to 0.5 bpp as the type of high bit rate (type 1) and that less than 0.5 bpp as the type of low bit rate (type 0). The posttransform is performed using DCT basis and   l
-
1-norm minimization when coding type 1 and Hadamard basis and l
-
0-norm minimization when coding type 0.Figure 6
The proposed posttransform architecture.


## 2.1. The Imaging Principle of MSCFBs

Figure1 shows the imaging principle of MSCFBs. The long linear multispectral CCD has four parallel CCD linear arrays. Each array is one imaging spectral band. Four imaging spectral bands are red, blue, green, and near-infrared, respectively. The light radiated or reflected by the hundreds of kilometers of linear arrays of ground pixels is concentrated onto optical thin film of CCD detector through an optical system. The space and spectrum distributing of ground targets radiation, acquired by CCD detector, can be expressed as
(1)
B
(
x
,
y
,
θ
,
λ
)
=
1
π
[
E
0
(
θ
,
λ
)
ρ
(
x
,
y
,
λ
)
τ
α
(
λ
)
*
h
α
(
x
,
y
,
λ
)
]
+
N
(
θ
,
λ
)
,
where θ is the solar height angle, λ is wavelength, E
0 is ground luminometer, ρ is spectral reflectance, τ
α is atmospheric transmittances, h
α is the point spread function, and N is radiant flux. The optical thin film of each linear array CCD allows the light of corresponding wavelength through. B (x
,
y
,
θ
,
λ) of each band is captured by a linear CCD array, and one dimensional spatial and spectral information is gained. When the CCD camera scans the ground target, spatial dimension information is gained. Therefore, MSCFBs are considered as 3D images, Y
∈
R
I
1
×
I
2
×
I
3, which is shown in Figure 2.Figure 1
Imaging principle of MSCFBs.Figure 2
MSCFBs properties.Based on the imaging principle of MSCFBs, MSCFBs have the spectral redundancies between the adjacent spectral bands and spatial redundancies between the adjacent pixels in each spectral band. Therefore, the compression algorithm must remove the spectral redundancies and spatial redundancies simultaneously.


## 2.2. Spectral Decorrelation

MSCFBs have two types of correlations: spatial and spectral correlations. The spatial correlation is defined as the correlation between the adjacent pixels in each spectral band. The spatial line correlation of MSCFBs can be expressed as(2)
R
(
p
,
0
)
=
∑
x
=
1
M
∑
y
=
1
N
(
f
(
x
,
y
)
-
f
¯
)
×
(
f
(
x
+
p
,
y
)
-
f
¯
)
∑
x
=
1
M
∑
y
=
1
N
(
f
(
x
,
y
)
-
f
¯
)
2
.The spatial column correlation of MSCFBs can be expressed as(3)
R
(
0
,
q
)
=
∑
x
=
1
M
∑
y
=
1
N
(
f
(
x
,
y
)
-
f
¯
)
×
(
f
(
x
,
y
+
q
)
-
f
¯
)
∑
x
=
1
M
∑
y
=
1
N
(
f
(
x
,
y
)
-
f
¯
)
2
,
where M is the line number of each band, N is the pixel number of each line, f
(
x
,
y
) is the function of the band gray, and f
¯ is the average value of the band gray. The spectral correlation is defined as the correlation between the pixels located in the same spatial position in the different spectral band, which can be expressed as(4)
h
(
p
,
q
)
=
(
∑
x
=
1
M
∑
y
=
1
N
(
f
(
x
+
p
,
y
+
q
)
-
f
¯
)
00000000
×
(
g
(
x
,
y
)
-
g
¯
)
∑
y
=
1
N
)
×
(
(
∑
x
=
1
M
∑
y
=
1
N
[
f
(
x
,
y
)
-
f
¯
]
2
)
0000
×
(
∑
x
=
1
M
∑
y
=
1
N
[
g
(
x
,
y
)
-
g
¯
]
2
)
)
-
1
/
2
,
where f
(
x
,
y
) and g
(
x
,
y
) are the value of pixel located in (
x
,
y
) in two adjacent bands and h
(
0,0
) is the cross-correlation function of f
(
x
,
y
) and g
(
x
,
y
). Based on (2)~(4), the Landsat-7 Thematic Mapper (TM) images are chosen to perform an experiment of the spatial and spectral correlation. Figure 3 shows the test results.The spatial and spectral correlation test result.
(a)
Spatial correlation test result
(b)
Spectral correlation test resultThe correlations between the adjacent bands for an MSCFB are in the range[
0.318
–
0.926
], and the line and column correlations between the adjacent pixels for an MSCFB are in the ranges [
0.685
–
0.955
] and [
0.609
–
0.963
], respectively. Therefore, the MSCFBs have higher spatial correlation than spectral correlation. Because the spectral correlations are comparatively weak, the traditional 3D-DWT-based compression algorithms are not suitable for MSCFBs. When the 3D-DWT-based approaches compress Landsat-7 images, the good compression results are not obtained. Simultaneously removing the spectral redundancies and spatial redundancies is primary target for the successful implementation of a high-performance compression method.Because the MSCFBs have a smaller number of bands, 1D-DCT is performed on the spectral dimension to decorrelate the spectral bands. The number of bands is denoted asP; the size of each band is denoted as M
×
N. An MSCFB is composed of M
×
N 1D spectral vector (Length P). The definition of the 1D-DCT for an input 1D spectral vector A
i (i
=
1,2
,
…
,
M
×
N) and output image B
i  (i
=
1,2
,
…
,
M
×
N) is
(5)
B
i
r
=
α
r
∑
k
=
0
P
-
1
A
i
k
cos
⁡
(
π
(
2
k
+
1
)
r
2
P
)
,
α
r
=
{
1
P
,
r
=
0
,
2
P
,
1
≤
r
≤
P
-
1
.
The B
i
r values are called the DCT coefficients of A
i. Equation (5) is performed M
×
N times to complete the spectral transform. In order to speed up calculation, we proposed an iterative method. The iterative equation can be expressed as
(6)
B
i
+
1
r
=
∑
k
=
0
P
-
1
A
i
+
1
k
∑
k
=
0
P
-
1
A
i
k
×
B
i
r
.
The DCT coefficients of the current spectral vector can be computed by the intermediate and output data obtained from the previous spectral vector.


## 2.3. Spatial Decorrelation

The 2D-DWT can decompose the image into multilevel wavelet to obtain lower resolution subband and detail subbands. The multilevel decomposition can be regarded as the process of low pass and high pass filtering. At each level, the low subband produced by previous level performs the high pass filtering to obtain detail information called wavelet coefficients, while performing the low pass filtering linked with scaling function to obtain the approximate information called scaling coefficients. The DWT has been widely employed exploiting the spatial correlations for remote sensing image, such as JPEG2000 and CCSDS-IDC. In this paper, we apply a 2D-DWT coupled with a posttransform to each band of MSCFBs. In our approach, the 2D-DWT is performed on each image band to reduce spatial correlations, and then the transformed wavelet coefficients are transformed by the posttransform approach to reduce remaining correlations in wavelet subbands.After remote sensing MSCFBs are transformed by 2D-DWT, there exits residual correlation between coefficients in a small neighborhood (see Figure4). Statistical dependencies between DWT coefficients have been studied for many years. In [25], correlations between nearby wavelet coefficients are reported to be in the range [
0.01
–
0.54
] at a distance less than or equal to 3 pixels. We found even wider range, and here we provide a more detailed discussion regarding this topic. We use a Pearson correlation coefficient ρ to analyze statistical dependency of DWT coefficients, which can be expressed as
(7)
ρ
x
,
y
=
δ
x
y
(
δ
x
δ
y
)
,
where δ
x
,
y denotes the covariance between the variables X and Y and δ
x and δ
y denote the standard deviation of X and Y, respectively. According to project experience, three-level 2D-DWT is appropriate for on-board compressor and we used three-level 2D-DWT in this paper.Figure 4
The wavelet transform of high-resolution image with a zoom on some coefficients.Three-level 2D-DWT is performed on each image band to produce one low-frequency subband (denoted as LL) and nine high-frequency subbands (denoted asLH
⁡
L, HL
L, HH
L, L
=
1,2
,
3). The test is performed from ten MSCFBs. The residual correlations between wavelet coefficients in a small neighborhood of wavelet subbands are shown as in Figure 5. In Level L (L
=
1), the residual directional correlation within 16-connected region is ρ (HL1, HH1, LH1, |
Δ
x
|
≤
2, |
Δ
y
|
≤
2) > 0.3. In Level L (L
=
2), the residual directional correlation within 16-connected region is ρ (HL2, HH2, LH2, |
Δ
x
|
≤
2, |
Δ
y
|
≤
2) > 0.4. In Level L (L
=
3), the residual directional correlation within 16-connected region is ρ (HL3, HH3, LH3, |
Δ
x
|
≤
2, |
Δ
y
|
≤
2) > 0.3. These imply a strong redundancy among neighboring 4 coefficients. The results in Figure 5 also indicate that the residual correlations between nearby wavelet coefficients are much weaker at a distance equal to 5 pixels.The correlation between the adjacent wavelet coefficients.
(a)
L
=
1
(b)
L
=
2
(c)
L
=
3To obtain the optimum compression performance for MSCFBs, the compression algorithms must fully consider the abovementioned statistical properties. In [26, 27], EBCOT has been reported to be very efficient method to remove these residual redundancy values. However, its implementation complexity is too high. In [28], Delaunay has proposed a posttransform to remove remaining redundancies between the adjacent wavelet coefficients within a small region. After the wavelet transform of an image, wavelet coefficients or scaling coefficients subbands are divided into several coefficients blocks. Each block has 4 × 4 wavelet coefficients. Posttransform performs on each coefficients block to remove the residual redundancies. From testing a large number of remote sensing images, the selected size of wavelet coefficients block is the best to a low complexity and efficient compression method. There are two major reasons: (1) the correlation between the adjacent wavelet coefficients is very weak when a distance is greater than or equal to 5 pixels; (2) the larger size of wavelet coefficients block leads to the higher computational complexity, and the smaller size of block results in increasing the number of coefficients blocks and the size information. In addition, since the posttransform is performed on coefficients blocks in wavelet transform domain, not on image domain, the compression algorithm has no additional overhead time to perform deblocking filter.The core idea behind posttransform compression is that wavelet coefficients blocks further are transformed using one group particular direction basis (such as Bandelet, DWT, DCT, and PCA) in a dictionary. First, a 2D-DWT is applied to an image. Next, blocks of 4 × 4 DWT coefficients are projected onN
B orthonormal bases B
b of the dictionary D. Then, the minimizing of the Lagrangian cost function is calculated to select the best posttransformed block. Finally, the posttransformed coefficients are encoded by entropy coding method. Each 4 × 4 DWT coefficients block f is considered as an input vector f
=
[
f
1
,
f
2
,
…
,
f
M
]
T
∈
R
M with M
=
16. The M vectors of the basis B
b are noted ϕ
m
b with m
∈
[
1
,
M
-
1
]. The posttransformed block f
b can be expressed as follows:
(8)
f
b
=
∑
m
=
1
16
〈
f
,
ϕ
m
b
〉
·
ϕ
m
b
.
Since a dictionary has N
B bases, N
B
+
1 (including one original block) posttransformed blocks f
b can be obtained. Among all the posttransformed blocks f
b, the best posttransformed block f
b
* is selected according to calculating the minimizing of the Lagrangian rate-distortion cost. The minimizing of the Lagrangian rate-distortion cost can be expressed as
(9)
L
(
f
q
b
)
=
D
(
f
q
b
)
+
λ
·
R
(
f
q
b
)
,
where f
q
b denotes the quantized posttransformed coefficients, q is the quantization step, and D
(
f
q
b
) denotes the square error between the posttransformed coefficients f
b and the quantized posttransformed coefficients f
q
b. λ is a Lagrangian multiplier and R
(
f
q
b
) denotes the required bit rates for encoding f
q
b and the associated side information b
*.In the posttransform, the dictionary has multiple bases. The better the compression performance is, the more the number of bases is and the higher the computational complexity is. However, on-board compression requires a low computational complexity. The space MSCFBs compressor allows the dictionaryD of posttransform to have only one basis. In [25], Delaunay et al. proposed a simple posttransform-based compression method. They only perform one-level wavelet transform and only use a Hadamard basis. Then, the posttransform only performs on the coefficients subbands. Compared with the compression method only using the DWT, the PSNR can only improve 0.4 dB~0.6 dB. In this paper, to obtain a posttransform with the low computational complexity and still have high compression performance, we thus consider a very simple dictionary containing only one dynamic base, which is Hadamard basis at the low bit rates, and DCT basis at the high bit rates.In [29], Delaunay has shown that the Lagrangian approach for selecting the best posttransformed block has two main drawbacks. First, the estimation of the required bit rate for R
(
f
q
b
) is computationally intensive and not always accurate. Second, calculating the Lagrangian function depends on the quantization step q. However, the coder does not define the quantization step when coding. Therefore, Delaunay proposed an l
1-norm minimization approach to select the best posttransformed block. In this paper, we propose a new method based on l
p-norm minimization approach to evaluate best posttransform, which replaces an l
1-norm minimization. It is very suitable for the space MSCFBs compressor which has the low complexity constraints.At the low compression bit rate, the bit rateR can be expressed as
(10)
R
≈
γ
0
M
with
γ
0
=
7
,
where M is the number of nonzero posttransformed coefficients. We propose an l
0-norm (p
=
0) minimization approach to select the best posttransformed block. The selected best posttransformed block has the minimal sum of coefficients magnitude, which can be expressed as
(11)
f
b
*
=
arg
min
⁡
f
b
,
b
∈
[
0
,
N
B
]
∥
f
b
∥
0
,
∥
f
b
∥
0
=
∑
m
=
1
M
I
m
with
I
m
=
{
1
,
⌊
|
a
b
[
m
]
|
⌋
≠
0
,
0
,
⌊
|
a
b
[
m
]
|
⌋
=
0
,
where a
b
[
m
] is a posttransformed coefficient. At the high compression bit rate, the bit rate R can be expressed as
(12)
R
≈
-
∑
i
=
1
N
p
i
lo
g
2
p
i
-
lo
g
2
q
.
We use the l
1-norm (p
=
1) minimization approach to select the best posttransformed block. The selected best posttransformed block has the minimal sum of coefficients magnitude, which can be expressed as
(13)
f
b
*
=
arg
min
⁡
f
b
,
b
∈
[
0
,
N
b
]
∥
f
b
∥
1
s
.
t
.
∥
f
b
∥
1
=
(
∑
m
=
0
15
|
a
b
[
m
]
|
)
.Figure6 shows the proposed posttransform architecture. Each high-frequency subband is performed by posttransform. The bit-rate comparator decides the type of coding bit rate. We define the coding bit rate greater than or equal to 0.5 bpp as the type of high bit rate (type 1) and that less than 0.5 bpp as the type of low bit rate (type 0). The posttransform is performed using DCT basis and   l
-
1-norm minimization when coding type 1 and Hadamard basis and l
-
0-norm minimization when coding type 0.Figure 6
The proposed posttransform architecture.


## 3. Tucker Decomposition in Posttransform Domain

In [30], the posttransformed blocked coefficients are encoded by an adaptive arithmetic coding method. In [29], the posttransformed blocked coefficients are encoded by the bit-plane encoding (BPE) method. However, the posttransforms destroy the zero-tree structure. The PSNR using BPE is less than 0.2 dB when compared with the method using DWT and even worse at high bit rates. A basis vector ordering approach has been used. This ordering is obtained by training or learning with thousands of wavelet coefficients blocks from a large number of remote sensing images. The ordering-based strategy has two major shortcomings for space MSCFBs compression. For the first weakness, the ordering of basis vector has huge calculations and is not always accurate. For the second weakness, the ordering highly depends on a training set with thousands of wavelet coefficients blocks of remote sensing images. However, the training set is not captured when coding.Indeed, after 1D-DCT, 2D-DWT, and posttransform, the MSCFBs still have residual spatial and spectral correlation. In order to obtain farther excellent compression performance, we consider the posttransform coefficients of MSCFBs as a three-order 3D tensor. The tensor performs Tucker decomposition (TD) to remove residual spatial and spectral redundancies simultaneously.A third-order tensor is denoted asY
∈
R
I
1
×
I
2
×
I
3. The third-order tensor is decomposed into a low dimension core tensor multiplied by three factor matrices using TD (see Figure 7). The core tensor is denoted as G
∈
R
J
1
×
J
2
×
J
3. The factor matrices are denoted as A
(
n
)
=
[
a
1
(
n
)
,
a
2
(
n
)
,
…
,
a
J
n
(
n
)
]
∈
R
I
n
×
J
n (n
=
1,2
,
3).Figure 7
Third-order tensor Tucker decomposition.The third-order TD can be expressed as follows [31]:
(14)
A
=
∑
j
1
=
1
J
1
∑
j
2
=
1
J
2
∑
j
3
=
1
J
3
g
j
1
j
2
j
3
a
j
1
(
1
)
∘
a
j
2
(
2
)
∘
a
j
3
(
3
)
+
E
=
G
×
1
A
(
1
)
×
2
A
(
2
)
×
3
A
(
3
)
+
E
=
G
×
{
A
}
+
E
=
Y
^
+
E
≈
Y
^
,
where J
1, J
2, J
3 are dimensions of the core tensor G, tensor Y
^ is an approximation of tensor Y, and tensor E is an error tensor. The tensor Y
^ depends on the (J
1
,
J
2
,
J
3) values. Equation (14) can be solved by the following optimal problem:
(15)
min
⁡
1
2
∥
Y
-
G
×
1
A
(
1
)
×
2
A
(
2
)
×
3
A
(
3
)
∥
F
2
s
.
t
.
.
G
∈
R
J
1
×
J
2
×
J
3
,
A
(
n
)
∈
R
I
n
×
J
n
,
n
=
1,2
,
3
.The size of each band is denoted asI
1
×
I
2, and the number of band is denoted as I
3. Each band performs three-level 2D-DWT to produce nine high-frequency wavelet subbands (HL
L, LH
⁡
L, HH
L, L
=
1,2
,
3) and one low-frequency wavelet subband (LL). All high subbands perform the proposed posttransform. The size of each subband is I
1
/
2
L
×
I
2
/
2
L. The size of posttransformed block is 4
×
4. All bands HL
L, LH
⁡
L, and HH
L (L
=
1,2
,
3), which have been posttransformed, are considered as one tensor, respectively. Here, 9 tensors are denoted as LH
⁡
L
′, HL
L
′, and HH
L
′ (L
=
1,2
,
3), respectively (see Figure 8). The size of each tensor is I
1
/
2
L
×
I
2
/
2
L
×
I
3. Note that the elements of each tensor are the posttransformed coefficients.Figure 8
High-frequency subband as a third-order tensor after transform and posttransformed block as a subtensor.Since the spatial dimension of each tensor isI
1
/
2
L
×
I
2
/
2
L,  the tensor decomposition computing cost is very high. To reduce the computational complexity, in this paper, we propose a multilevel TD approach. First, each tensor is divided into several subtensors (see Figures 8 and 9). In each tensor, each subband is one posttransformed block. Therefore, the size of each subtensor is 4
×
4
×
I
3. Each tensor contains I
1
/
2
L
+
2
×
I
2
/
2
L
+
2 subtensors, and the 9 tensors have ∑
L
=
1
3
(
I
1
/
2
L
+
2
×
I
2
/
2
L
+
2
) three-order subtensors. In the first level, all subtensors for each tensor are performed by TD to produce the I
1
/
2
L
+
2
×
I
2
/
2
L
+
2 core tensor and 3
×
I
1
/
2
L
+
2
×
I
2
/
2
L
+
2 factor matrices. The produced core tensors are organized into a new tensor, which is performed by TD in the next level. In the proposed codec we are using one level of TD. Since the spatial dimension of each subtensor is 4 × 4, TD of each subtensor has low complexity.Figure 9
One-level TD.To find optimal component matricesA and G, we apply HALS-NTD algorithm [31] to each subtensor in this paper. The HALS-NTD algorithm can be shown as Algorithm 1.Algorithm 1
Input: (1) tensor Y; (2) tensor size I
1
×
I
2
×
I
3, where I
1
=
4
,
I
2
=
4; (3) core tensor size (
J
1
,
J
2
,
J
3
).
Output:  (1) 3 factors A
(
n
)
∈
R
I
n
×
J
n (n
=
1,2
,
3); (2) a core tensor G
∈
R
J
1
×
J
2
×
J
3.
Begin
Initialization A
(
n
)
G for n
=
1,2
,
3;
Normalize all a
j
n
(
n
) for n
=
1,2
,
3;
E
_
=
Y
-
Y
^
Do {
For n
=
1 to 3 do
For j
n
=
1 to J
n  do
Y
(
n
)
(
j
n
)
=
E
(
n
)
+
a
j
n
(
n
)
[
G
(
n
)
]
j
n
U
⊗
-
n
T
a
j
n
(
n
)
⟵
[
Y
(
n
)
(
j
n
)
[
(
G
×
-
n
{
A
}
)
(
n
)
]
j
n
T
]
+
a
j
n
(
n
)
⟵
a
j
n
(
n
)
∥
a
j
n
(
n
)
∥
1
E
n
⟵
Y
(
n
)
(
j
n
)
-
u
j
n
(
n
)
[
G
(
n
)
]
j
n
A
⊗
-
n
T
End
End
G
⟵
[
G
⊛
(
Y
⊘
Y
^
)
×
{
A
T
}
]
+
For each j
(
n
)
=
1,2
,
…
,
J
(
n
), For n
=
1,2
,
3 do
g
j
1
j
2
j
3
⟵
g
j
1
j
2
j
3
+
E
×
1
a
j
1
(
1
)
×
2
a
j
2
(
2
)
×
3
a
j
3
(
3
)
E
⟵
E
+
Δ
g
j
1
j
2
j
3
a
j
1
(
1
)
∘
a
j
2
(
2
)
∘
a
j
3
(
3
)
End
}  While  (the stopping condition is reached)
end


### 3.1. Deep Coupling between Posttransform and Tucker Decomposition

The dimensions of each tensor or each subtensor (J
1
,
J
2
,
J
3) are determined by compression ratio (CR). The lower values of J
1
,
J
2
,
J
3 make compression algorithm have a higher CR. The dimensions of 9 tensors, that is, LH
⁡
L
′, HL
L
′, and HH
L
′ (L
=
1
,
2,3), are denoted as J
1
LH
⁡
L
′, J
2
HL
L
′, and J
3
HH
L
′ (L
=
1
,
2
,
3). The CR can be defined as a ratio value between the required bit numbers for representing original image data and the required bit numbers for representing compressed bit-streams, which can be expressed as
(16)
CR
=
I
1
×
I
2
×
I
3
B
,
B
=
∑
L
=
1
3
(
I
1
2
×
J
1
LH
⁡
L
′
+
I
2
2
×
J
2
LH
⁡
L
′
+
I
3
×
J
3
LH
⁡
L
′
0000
+
J
1
LH
⁡
L
′
×
J
2
LH
⁡
L
′
×
J
3
LH
⁡
L
′
0000
+
I
1
2
×
J
1
HL
L
′
+
I
2
2
×
J
2
HL
L
′
+
I
3
×
J
3
HL
L
′
0000
+
J
1
HL
L
′
×
J
2
HL
L
′
×
J
3
HL
L
′
0000
+
I
1
2
×
J
1
HH
L
′
+
I
2
2
×
J
2
HH
L
′
+
I
3
×
J
3
HH
L
′
0000
+
J
1
HH
L
′
×
J
2
HH
L
′
×
J
3
HH
L
′
I
2
2
)
.In order to efficiently determine the size of all core tensors of subtensor when each tensor performs TD, we proposed a deep coupling between posttransform and TD not only to determine the size of the core tensor and reduce the side information of posttransform but also to code these core tensors and complete the bit-rate control. Figure10 shows the proposed deep coupling between the posttransform and TD.Figure 10
Proposed deep coupling algorithm.The bit-rate allocation module allocates the bit rates for each tensor. The information evaluation module evaluates the information of each tensor. Now the target bit rates for different tensors can be allocated based on their information contents. Then, the allocated bit rates of each tensor are allocated to each subtensor, and the dimensions of each subtensor can be determined according to its allocated bit rates.In the first place, the information content inLL subband I
LL is calculated through an l
-
1-norm approach. Let f
(
x
,
y
) denote the coefficients in LL subband. The information content in LL subband can be calculated as
(17)
I
LL
=
∥
f
LL
∥
1
=
∑
x
=
1
I
1
/
8
∑
y
=
1
I
2
/
8
|
f
(
x
,
y
)
|
.In the second place, the information content of 9 tensors (LH
⁡
L
′, HL
L
′, HH
L
′, L
=
1,2
,
3) is calculated. Let I
LH
⁡
L
′, I
HL
L
′, and I
HH
L
′, L
=
1,2
,
3, denote the information content of 9 tensors, respectively. Since the selected representation f
b
* reflects the image information, the information content of 9 tensors can be evaluated through f
b
*. The dimensions of each tensor are I
1
×
2
-
L
×
I
2
×
2
-
L
×
I
3. Each tensor contains I
1
×
2
-
L
-
2
×
I
2
×
2
-
L
-
2 posttransformed blocks. Let f
i
b
* denote the selected representation of ith posttransform block. The information content of the tensor LH
⁡
L
′ can be calculated as
(18)
I
LH
⁡
L
′
=
∑
i
=
1
I
1
×
I
2
×
2
-
2
(
L
-
2
)
f
i
b
*
.
Other tensors can be likewise calculated.In the third place, the weight of bit-rate allocation for each tensor is acquired through(19)
r
LH
⁡
L
′
=
I
LH
⁡
L
′
I
LL
+
∑
L
=
1
3
LH
⁡
L
′
+
HL
L
′
+
HH
L
′
,
r
HL
L
′
=
I
HL
L
′
I
LL
+
∑
L
=
1
3
LH
⁡
L
′
+
HL
L
′
+
HH
L
′
,
r
HH
L
′
=
I
H
H
L
′
I
LL
+
∑
L
=
1
3
LH
⁡
L
′
+
HL
L
′
+
HH
L
′
.Finally, the allocated bit rates of each tensor are allocated to each subtensor. The bit rate assigned to theith subtensor can be calculated as
(20)
r
i
LH
⁡
L
′
=
f
i
b
*
∑
j
=
1
I
1
×
I
2
×
2
-
2
(
L
-
2
)
f
j
b
*
×
r
LH
⁡
L
′
,
r
i
HL
L
′
=
f
i
b
*
∑
j
=
1
I
1
×
I
2
×
2
-
2
(
L
-
2
)
f
j
b
*
×
r
HL
L
′
,
r
i
HH
L
′
=
f
i
b
*
∑
j
=
1
I
1
×
I
2
×
2
-
2
(
L
-
2
)
f
j
b
*
×
r
HH
L
′
.The side informationb
* denotes the posttransform used on this block. Since the posttransform has only one basis at low and high bit rates, the side information is one bit per block and is not entropy coded. Table 1 shows the side information cost. The side information always has 1 bpb of bit-rate overhead in whole coding process. The side information occupies 12.5% of the total bit rate at 0.5 bpp, 4.16% at 1.5 bpp, and only 2.08% at 3 bpp, respectively. Therefore, the side information has more effect on compression performance at low bit rates than at high bit rates. Indeed, these nonsignificant blocks do not require the side information because their coefficients are all zeros. Coincidentally, there are more nonsignificant blocks at low bit rate than high bit rate. Therefore, the allocated bit rate for side information can be reduced. As the number of nonsignificant blocks decreases, the required bit rate of side information will certainly increase.Table 1
Side information cost associated to the lp norm approach.
Target bit rate
0.5 bpp
1.5 bpp
3.0 bpp
Side information cost (bit per block)
1.0 bit
1.0 bit
1.0 bit
Side information budget (% of the total bit rate)
12.5%
4.16%
2.08%In this paper, we use an improved BPE method to encode the core tensor. And then the side information is integrated into the bit-streams. We modify the SPITH algorithm [32] to complete the bit-plane coding. The bit-plane coding processes each core tensor of each subtensor at a time. After one core tensor is processed, the next core tensor is processed. For each time coding, the improved BPE defines two encoding paths: (1) significant path (SP) and (2) the refining path (RP). We define a significant state threshold T. The element of core tensor is denoted as g
j
1
j
2
j
3 at the location (j
1
,
j
2
,
j
3). We use S
n
(
T
) to denote an index of significant state when the threshold is T
=
2
n in the nth bit plane. The S
n
(
T
) can be expressed as
(21)
S
T
(
n
)
=
{
1
,
T
≤
|
g
j
1
j
2
j
3
|
≤
2
T
,
0
,
otherwise
.
In the nth bit plane, when S
T
(
n
)
=
1, the g
j
1
j
2
j
3 is a significant element, which passes into the SP. When S
T
(
n
)
=
0, the g
j
1
j
2
j
3 is insignificant and passes into the RP. The elements in RP perform BPE in the next bit plane. In the next bit plane, n
=
n
-
1, T
=
2
n
-
1. The coding method is the same with the previous bit plane. When traversing all the bit planes, the elements in SP are moved from the core tensor and encoded by entropy coding approach. If the SP for a core tensor has any elements, the side information is embedded into compressed bit-streams. Otherwise, the side information is not embedded into compressed bit-streams.


### 3.2. Our Proposed Architecture

Figure11 shows the whole architecture of the proposed posttransform Tucker decomposition- (PTTD-) based compression method in this paper. The whole process of PTTD-based approach can be divided into four steps as follows.Figure 11
Compression algorithm architecture for MSCFBs using the proposed PTTD method.In the first step, MSCFBs perform 1D-DCT in the spectral dimension to remove the spectral correlation.In the second step, each spectral band of the MSCFBs performs 3-level 2D-DWT to remove the spatial correlation. And 10 wavelet subbands (one low-frequency and nine high-frequency subbands) for each band can be obtained.In the third step, all high-frequency subbands are divided into several blocks. Each block has 4 × 4 DWT coefficients. Each block performs the proposed posttransform.In the fourth step, each posttransformed subtensor performs TD to remove the residual spectral and spatial redundancies.In the fifth step, the posttransformed coefficients in low-frequency subband are encoded by DPCM. The factor matrices perform quantization. The core tensors perform the bit-plane coding and entropy coding using adaptive arithmetic coding.In the sixth step, the bit streams are packed and transmitted to the decoder via signal channels.


## 3.1. Deep Coupling between Posttransform and Tucker Decomposition

The dimensions of each tensor or each subtensor (J
1
,
J
2
,
J
3) are determined by compression ratio (CR). The lower values of J
1
,
J
2
,
J
3 make compression algorithm have a higher CR. The dimensions of 9 tensors, that is, LH
⁡
L
′, HL
L
′, and HH
L
′ (L
=
1
,
2,3), are denoted as J
1
LH
⁡
L
′, J
2
HL
L
′, and J
3
HH
L
′ (L
=
1
,
2
,
3). The CR can be defined as a ratio value between the required bit numbers for representing original image data and the required bit numbers for representing compressed bit-streams, which can be expressed as
(16)
CR
=
I
1
×
I
2
×
I
3
B
,
B
=
∑
L
=
1
3
(
I
1
2
×
J
1
LH
⁡
L
′
+
I
2
2
×
J
2
LH
⁡
L
′
+
I
3
×
J
3
LH
⁡
L
′
0000
+
J
1
LH
⁡
L
′
×
J
2
LH
⁡
L
′
×
J
3
LH
⁡
L
′
0000
+
I
1
2
×
J
1
HL
L
′
+
I
2
2
×
J
2
HL
L
′
+
I
3
×
J
3
HL
L
′
0000
+
J
1
HL
L
′
×
J
2
HL
L
′
×
J
3
HL
L
′
0000
+
I
1
2
×
J
1
HH
L
′
+
I
2
2
×
J
2
HH
L
′
+
I
3
×
J
3
HH
L
′
0000
+
J
1
HH
L
′
×
J
2
HH
L
′
×
J
3
HH
L
′
I
2
2
)
.In order to efficiently determine the size of all core tensors of subtensor when each tensor performs TD, we proposed a deep coupling between posttransform and TD not only to determine the size of the core tensor and reduce the side information of posttransform but also to code these core tensors and complete the bit-rate control. Figure10 shows the proposed deep coupling between the posttransform and TD.Figure 10
Proposed deep coupling algorithm.The bit-rate allocation module allocates the bit rates for each tensor. The information evaluation module evaluates the information of each tensor. Now the target bit rates for different tensors can be allocated based on their information contents. Then, the allocated bit rates of each tensor are allocated to each subtensor, and the dimensions of each subtensor can be determined according to its allocated bit rates.In the first place, the information content inLL subband I
LL is calculated through an l
-
1-norm approach. Let f
(
x
,
y
) denote the coefficients in LL subband. The information content in LL subband can be calculated as
(17)
I
LL
=
∥
f
LL
∥
1
=
∑
x
=
1
I
1
/
8
∑
y
=
1
I
2
/
8
|
f
(
x
,
y
)
|
.In the second place, the information content of 9 tensors (LH
⁡
L
′, HL
L
′, HH
L
′, L
=
1,2
,
3) is calculated. Let I
LH
⁡
L
′, I
HL
L
′, and I
HH
L
′, L
=
1,2
,
3, denote the information content of 9 tensors, respectively. Since the selected representation f
b
* reflects the image information, the information content of 9 tensors can be evaluated through f
b
*. The dimensions of each tensor are I
1
×
2
-
L
×
I
2
×
2
-
L
×
I
3. Each tensor contains I
1
×
2
-
L
-
2
×
I
2
×
2
-
L
-
2 posttransformed blocks. Let f
i
b
* denote the selected representation of ith posttransform block. The information content of the tensor LH
⁡
L
′ can be calculated as
(18)
I
LH
⁡
L
′
=
∑
i
=
1
I
1
×
I
2
×
2
-
2
(
L
-
2
)
f
i
b
*
.
Other tensors can be likewise calculated.In the third place, the weight of bit-rate allocation for each tensor is acquired through(19)
r
LH
⁡
L
′
=
I
LH
⁡
L
′
I
LL
+
∑
L
=
1
3
LH
⁡
L
′
+
HL
L
′
+
HH
L
′
,
r
HL
L
′
=
I
HL
L
′
I
LL
+
∑
L
=
1
3
LH
⁡
L
′
+
HL
L
′
+
HH
L
′
,
r
HH
L
′
=
I
H
H
L
′
I
LL
+
∑
L
=
1
3
LH
⁡
L
′
+
HL
L
′
+
HH
L
′
.Finally, the allocated bit rates of each tensor are allocated to each subtensor. The bit rate assigned to theith subtensor can be calculated as
(20)
r
i
LH
⁡
L
′
=
f
i
b
*
∑
j
=
1
I
1
×
I
2
×
2
-
2
(
L
-
2
)
f
j
b
*
×
r
LH
⁡
L
′
,
r
i
HL
L
′
=
f
i
b
*
∑
j
=
1
I
1
×
I
2
×
2
-
2
(
L
-
2
)
f
j
b
*
×
r
HL
L
′
,
r
i
HH
L
′
=
f
i
b
*
∑
j
=
1
I
1
×
I
2
×
2
-
2
(
L
-
2
)
f
j
b
*
×
r
HH
L
′
.The side informationb
* denotes the posttransform used on this block. Since the posttransform has only one basis at low and high bit rates, the side information is one bit per block and is not entropy coded. Table 1 shows the side information cost. The side information always has 1 bpb of bit-rate overhead in whole coding process. The side information occupies 12.5% of the total bit rate at 0.5 bpp, 4.16% at 1.5 bpp, and only 2.08% at 3 bpp, respectively. Therefore, the side information has more effect on compression performance at low bit rates than at high bit rates. Indeed, these nonsignificant blocks do not require the side information because their coefficients are all zeros. Coincidentally, there are more nonsignificant blocks at low bit rate than high bit rate. Therefore, the allocated bit rate for side information can be reduced. As the number of nonsignificant blocks decreases, the required bit rate of side information will certainly increase.Table 1
Side information cost associated to the lp norm approach.
Target bit rate
0.5 bpp
1.5 bpp
3.0 bpp
Side information cost (bit per block)
1.0 bit
1.0 bit
1.0 bit
Side information budget (% of the total bit rate)
12.5%
4.16%
2.08%In this paper, we use an improved BPE method to encode the core tensor. And then the side information is integrated into the bit-streams. We modify the SPITH algorithm [32] to complete the bit-plane coding. The bit-plane coding processes each core tensor of each subtensor at a time. After one core tensor is processed, the next core tensor is processed. For each time coding, the improved BPE defines two encoding paths: (1) significant path (SP) and (2) the refining path (RP). We define a significant state threshold T. The element of core tensor is denoted as g
j
1
j
2
j
3 at the location (j
1
,
j
2
,
j
3). We use S
n
(
T
) to denote an index of significant state when the threshold is T
=
2
n in the nth bit plane. The S
n
(
T
) can be expressed as
(21)
S
T
(
n
)
=
{
1
,
T
≤
|
g
j
1
j
2
j
3
|
≤
2
T
,
0
,
otherwise
.
In the nth bit plane, when S
T
(
n
)
=
1, the g
j
1
j
2
j
3 is a significant element, which passes into the SP. When S
T
(
n
)
=
0, the g
j
1
j
2
j
3 is insignificant and passes into the RP. The elements in RP perform BPE in the next bit plane. In the next bit plane, n
=
n
-
1, T
=
2
n
-
1. The coding method is the same with the previous bit plane. When traversing all the bit planes, the elements in SP are moved from the core tensor and encoded by entropy coding approach. If the SP for a core tensor has any elements, the side information is embedded into compressed bit-streams. Otherwise, the side information is not embedded into compressed bit-streams.


## 3.2. Our Proposed Architecture

Figure11 shows the whole architecture of the proposed posttransform Tucker decomposition- (PTTD-) based compression method in this paper. The whole process of PTTD-based approach can be divided into four steps as follows.Figure 11
Compression algorithm architecture for MSCFBs using the proposed PTTD method.In the first step, MSCFBs perform 1D-DCT in the spectral dimension to remove the spectral correlation.In the second step, each spectral band of the MSCFBs performs 3-level 2D-DWT to remove the spatial correlation. And 10 wavelet subbands (one low-frequency and nine high-frequency subbands) for each band can be obtained.In the third step, all high-frequency subbands are divided into several blocks. Each block has 4 × 4 DWT coefficients. Each block performs the proposed posttransform.In the fourth step, each posttransformed subtensor performs TD to remove the residual spectral and spatial redundancies.In the fifth step, the posttransformed coefficients in low-frequency subband are encoded by DPCM. The factor matrices perform quantization. The core tensors perform the bit-plane coding and entropy coding using adaptive arithmetic coding.In the sixth step, the bit streams are packed and transmitted to the decoder via signal channels.


## 4. Experimental Results

### 4.1. Experimental Scheme

In order to test the performance of the proposed algorithm in this paper, we use the research and development ground test equipment to build experimental system. Figure12 shows the structure of experimental system to test the performance of the proposed algorithm. The experimental system is composed of image simulation resource, multispectral compression system, ground test equipment, DVI monitor, and server. The server injects the remote sensing image into the image simulation resource. The image simulation resource adjusts the size of image and line frequency of output to simulate the output of CCD. The multispectral compression system receives the test remote images from the image simulation resource and then compresses them to verify the proposed compression algorithm. The ground test equipment completes image decompression to gain the reconstructed images. The reconstructed images are transferred to the server through camera link bus. Finally, the compression performance is analyzed in the server.Figure 12
The structure of experiment system.The server is a high performance computer using 4.00 GHz and eight-core AMD processor and 28.0 GB of RAM. The image simulation resource has 64 GB NAND flash arrays to store remote sensing images. It outputs image as line-to-line way. The maximum number of pixels in each line is 4096, and the line frequency is 0~8 KHz. In the multispectral compression system, the bit rate can be set to 2.0~0.25 bpp. The working frequency of system is 88 MHz.


### 4.2. Multispectral Image Compression Validation and Analysis

In order to verify the validation of the proposed algorithm, the server injects a series of the AVIRIS multiband remote sensing images into the image simulation resource. The group of multispectral images have four bands; the size of each band is 512 × 512. The bit depth of pixels is 8 bpp (bits per pixel). The bit rate is set to 1 bpp. Figure13 demonstrated the four original and reconstructed bands. From the displayed images, the original image and reconstructed image are almost not different because the proposed compression algorithm has a high signal-to-noise ratio. The PSNR reaches 48.61 dB, 47.19 dB, 47.07 dB, and 46.38 dB, respectively, in Figures 13(e)–13(h). After each band is encoded by the proposed algorithm, the dynamic range of pixels is 0~2 bits. And most pixels are about 1 bit. So, the proposed compression algorithm is feasible for multispectral CCD image.Encoding multiband images at 1.0 bpp.
(a)
Original band 1
(b)
Original band 2
(c)
Original band 3
(d)
Original band 4
(e)
Reconstructed band 1
(f)
Reconstructed band 2
(g)
Reconstructed band 3
(h)
Reconstructed band 4


### 4.3. Multispectral Image Compression Algorithm Performance Analysis

To objectively evaluate the performance of the proposed deep coupling-based compression scheme, the server injects multiple groups of multispectral images having different texture characteristics into the image simulation resource. The compression rate is set to 2.0~0.25 bpp. We use the PSNR formula to evaluate the rate distortion performance of the proposed deep coupling-based compression approaches objectively. The PSNR can be expressed as(22)
PSNR
=
10
lo
g
10
(
P
MSE
)
(
dB
)
,
MSE
=
1
R
×
C
∑
i
=
0
R
-
1
∑
j
=
0
C
-
1
(
x
i
,
j
-
x
i
,
j
′
)
2
,
where P is equal to 2
B
-
1, B is the quantization bits of original images, x
i
,
j denotes the value of the pixel of original images at (
i
,
j
), x
i
,
j
′ denotes the value of the pixel of reconstructed images at (
i
,
j
), and R and C denote the number of rows and columns of multispectral images, respectively.In the first experimental part, in order to test the compressed performance of the proposed approach, we use three groups of SPOT-1 multispectral images having different texture characteristics. Each group is composed of red, green, blue, and near-infrared spectral bands. The size of each band is 512 × 512. The bit depth of pixels is 8 bpp (bits per pixel). Table2 demonstrates the tested PSNR results of our approach at different bit rates. According to Table 2, the average PSNR reaches upwards of 40 dB. Therefore, the proposed algorithm based on deep coupling approach has the good compression performance and can satisfy the requirements of design index.Table 2
Tested PSNR results of our compression.
Images
PSNR (dB)
Average PSNR (dB)
Bit rate (bits/pixel)
SPOT-1
42.99
42.53
0.25
SPOT-2
43.70
SPOT-3
40.90
SPOT-1
47.07
47.06
0.5
SPOT-2
48.86
SPOT-3
45.26
SPOT-1
50.95
51.27
1
SPOT-2
51.93
SPOT-3
50.93
SPOT-1
53.64
53.74
2
SPOT-2
54.75
SPOT-3
52.82In the second experimental, in order to compare compressed performance of our approach, we use several proposed compression algorithms, AT-3DSPIHT [33], 3D-DWT [34], 3D-PCA [35], 3D-SPECK [36], and SA-DCT [37], to compress remote sensing MSCFBs. We use the AVIRIS multiband remote sensing images in JPL laboratory. Each group is 256
pixels
×
288
pixels
×
4; the depth of every pixel is 8 bits. The images are Low Altitude, Lunar Lake, Jasper Ridge, and Cuprite. The compression bit rate is set to 2.0~0.25 bpp. Figure 14 shows the comparison results of average PSNR for different compression methods. From Figure 14, due to the full usage of posttransform and TD, the proposed PTTD-based compression algorithm has the best compression performance. The PTTD-based method can improve 0.3~1.3 dB of average PSNR compared with other compression methods at 2.0~0.25 bpp.The test and comparison result of multiband images: (a) is PSNR difference between 3D-DWT and other three compression codecs, and (b) is PSNR difference between 3D-PCA and other two compression codecs.
(a)
(b)Overall, the PTTD-based compression method has excellent compression performance. It can fully meet the requirement of space multispectral CCD camera.


### 4.4. The Spectral Distortion Analysis

We use the max similarity of spectrum (MSS) to evaluate the distortion of multispectral images for our algorithm. LetX be original remote sensing MSCFBs and let Y be reconstructed remote sensing MSCFBs. The MSS can be expressed as
(23)
h
MSS
=
max
⁡
i
,
j
f
(
i
,
j
)
2
+
(
1
-
c
(
i
,
j
)
2
)
2
,
f
(
i
,
j
)
=
1
n
λ
∑
λ
(
X
(
i
,
j
,
λ
)
-
Y
(
i
,
j
,
λ
)
)
2
,
c
(
i
,
j
)
=
∑
λ
(
X
(
i
,
j
,
λ
)
-
μ
i
,
j
)
(
Y
(
i
,
j
,
λ
)
-
μ
i
,
j
′
)
(
n
λ
-
1
)
δ
i
,
j
δ
i
,
j
′
,
where i, j, and λ are line, column, and band number of multispectral image, respectively, n
λ is the total number of bands, and μ
i
,
j and δ
i
,
j are the average value and mean square of the corresponding pixel in the spectral dimension, respectively.The greater the value of MSS is, the greater the distortion of spectra is. The MSS between original and reconstructed multispectral images for the image group using different approach is shown as in Figure15.Figure 15
Comparison of MSS performance.From Figure15, our approach has the lowest MSS, so the proposed algorithm can protect spectral information well.Figure16 demonstrates the reconstructed third band of two remote sensing images, which is 1024 × 1024; the depth of every pixel is 8 bits. The bit rate is set to 0.5 bpp.Figure 16
Comparison of reconstruction performance.In order to view the detailed information of reconstructed image, the part of image is100
×
100. From Figure 14, under the low bit rate, 3DSPIHT-based approach causes the visible ringing effects and the blurred edge. However, our approach does not cause the blurred edge. So, the proposed algorithm can protect spectral information well.


### 4.5. Proposed Algorithm Complexity Analysis and Compression Time

In our method, 1D-DCT is performed on spectral dimension to exploit the spectral information. It is reported that DCT has the lower amount of computations than DWT. When one of the most efficient methods is used to perform DCT, 8 × 8 points of 2D-DCT use only 462 time adding operations and 54 time multiplication operators. This means that each point of 2D-DCT needs only 0.84 time multiplication operators and 7.22 adding operations. Therefore,N bands of MSCFBs performing 1D-DCT need N
×
0.84 multiplications and N
×
7.22 additions.After applying 1D-DCT, three-level 2D-DWT is applied to the spatial bands. LetL be decomposition level of 2D-DWT. For the I
1
×
I
2
×
I
3 MSCFBs performing 2D-DWT, the calculation complexity is O
(
8
N
I
1
I
2
I
3
(
1
-
2
-
2
L
)
/
6
). In the PTTD-based method, a 9/7 2D-DWT is used to transform each band of MSCFBs. And we perform three levels of decomposition. Therefore, the calculation complexity of the PTTD-based method is O
(
9
I
1
I
2
I
3
/
7
).After MSCFBs performing 2D-DWT, we use the Hadamard posttransform in low bit rate. Let each DWT coefficients block haveM coefficients. The Hadamard-based posttransform needs the O (M
log
⁡
M) operations. In high bit rate, we use the DCT posttransform. 2D 4 × 4 DCT requires only 27 multiplications and 231 additions. In addition, the best posttransform selection uses only 2
×
(
M
-
1
) time adding operations and one comparison operator. Therefore, the complexity of posttransform stage is M
(
log
⁡
2
M
+
3
)
+
2 operations for each block with M coefficients.After 2D-DWT coefficients performing posttransform, the dimension of each subtensor is4
×
4
×
I
3. Let I
A be the average dimensions of core tensor G. The complexity of the calculation of tensor decomposition is O
(
J
A
I
A
3
), where
(24)
I
A
=
1
3
(
4
+
4
+
I
3
)
,
J
A
=
1
3
(
J
1
A
+
J
2
A
+
J
3
A
)
.
Therefore, the calculation complexity of TD in the PTTD-based method is O
(
J
A
I
A
3
)
×
∑
L
=
1
3
(
I
1
/
2
L
+
2
×
I
2
/
2
L
+
2
).In the following, we test and analyze the calculation times of the proposed PTTD-based compression method. Note that, in this test experiment, they only evaluate the calculation time of the PTTD-based compression because it is not the best optimized implementation on FPGA-based compression system. The evaluations testing performs the lossy compression for an image of size 3072 × 128 at 1.0 bpp on FPGA EVM board with system clock frequency at 88 MHz. Table3 demonstrates the comparison of complexity between the proposed multispectral compression algorithm and others.Table 3
The results of complexity comparison.
Methods
Times
Ours
0.042 us/sample
KTL [38]
0.102 us/sample
3DSPIHT [3]
0.062 us/sample
JPEG2000 [39]
0.181 us/sampleFrom Table3, the processing time of our algorithm reaches 0.042 us/sample; the data throughput is 23.81MSPS, which is higher than JPEG2000, KLT, and 3DSPIHT, so our approach has the lowest complexity. In our project, the multispectral CCD camera works at an orbits altitude of 500 km, the scroll angle is -
40
°~+
40
°, the latitudes are -
70
°~+
70
°, and the line working frequency is 1.8094 kHz~0.85946 kHz. In the line working frequency, capturing the image of 128 × 3072 required 70.74 ms. Using our approach, compressing the image of 128
×
3072 requires 16.51 ms. So, our approach can process the four bands images simultaneously. This meets the requirement of the project.In addition, we use the XC2V6000-6FF1152 FPGA to implement the proposed algorithm. The design language is VerilogHDL, the development platform is ISE8.2, and synthesis tool is XST. Table4 demonstrates the occupancy of resources of our approach.Table 4
The occupancy of resources.
Resources
Utilization rate
Slices
23655/33792 (70%)
Four input LUTs
45281/67584 (67%)
BRAM
116/144 (80%)From Table4, the LUTs occupy 67%, slices occupy 70%, and BRAM occupy 80%. Various indicators are lower than 95%, which meet the requirement of our project.


## 4.1. Experimental Scheme

In order to test the performance of the proposed algorithm in this paper, we use the research and development ground test equipment to build experimental system. Figure12 shows the structure of experimental system to test the performance of the proposed algorithm. The experimental system is composed of image simulation resource, multispectral compression system, ground test equipment, DVI monitor, and server. The server injects the remote sensing image into the image simulation resource. The image simulation resource adjusts the size of image and line frequency of output to simulate the output of CCD. The multispectral compression system receives the test remote images from the image simulation resource and then compresses them to verify the proposed compression algorithm. The ground test equipment completes image decompression to gain the reconstructed images. The reconstructed images are transferred to the server through camera link bus. Finally, the compression performance is analyzed in the server.Figure 12
The structure of experiment system.The server is a high performance computer using 4.00 GHz and eight-core AMD processor and 28.0 GB of RAM. The image simulation resource has 64 GB NAND flash arrays to store remote sensing images. It outputs image as line-to-line way. The maximum number of pixels in each line is 4096, and the line frequency is 0~8 KHz. In the multispectral compression system, the bit rate can be set to 2.0~0.25 bpp. The working frequency of system is 88 MHz.


## 4.2. Multispectral Image Compression Validation and Analysis

In order to verify the validation of the proposed algorithm, the server injects a series of the AVIRIS multiband remote sensing images into the image simulation resource. The group of multispectral images have four bands; the size of each band is 512 × 512. The bit depth of pixels is 8 bpp (bits per pixel). The bit rate is set to 1 bpp. Figure13 demonstrated the four original and reconstructed bands. From the displayed images, the original image and reconstructed image are almost not different because the proposed compression algorithm has a high signal-to-noise ratio. The PSNR reaches 48.61 dB, 47.19 dB, 47.07 dB, and 46.38 dB, respectively, in Figures 13(e)–13(h). After each band is encoded by the proposed algorithm, the dynamic range of pixels is 0~2 bits. And most pixels are about 1 bit. So, the proposed compression algorithm is feasible for multispectral CCD image.Encoding multiband images at 1.0 bpp.
(a)
Original band 1
(b)
Original band 2
(c)
Original band 3
(d)
Original band 4
(e)
Reconstructed band 1
(f)
Reconstructed band 2
(g)
Reconstructed band 3
(h)
Reconstructed band 4


## 4.3. Multispectral Image Compression Algorithm Performance Analysis

To objectively evaluate the performance of the proposed deep coupling-based compression scheme, the server injects multiple groups of multispectral images having different texture characteristics into the image simulation resource. The compression rate is set to 2.0~0.25 bpp. We use the PSNR formula to evaluate the rate distortion performance of the proposed deep coupling-based compression approaches objectively. The PSNR can be expressed as(22)
PSNR
=
10
lo
g
10
(
P
MSE
)
(
dB
)
,
MSE
=
1
R
×
C
∑
i
=
0
R
-
1
∑
j
=
0
C
-
1
(
x
i
,
j
-
x
i
,
j
′
)
2
,
where P is equal to 2
B
-
1, B is the quantization bits of original images, x
i
,
j denotes the value of the pixel of original images at (
i
,
j
), x
i
,
j
′ denotes the value of the pixel of reconstructed images at (
i
,
j
), and R and C denote the number of rows and columns of multispectral images, respectively.In the first experimental part, in order to test the compressed performance of the proposed approach, we use three groups of SPOT-1 multispectral images having different texture characteristics. Each group is composed of red, green, blue, and near-infrared spectral bands. The size of each band is 512 × 512. The bit depth of pixels is 8 bpp (bits per pixel). Table2 demonstrates the tested PSNR results of our approach at different bit rates. According to Table 2, the average PSNR reaches upwards of 40 dB. Therefore, the proposed algorithm based on deep coupling approach has the good compression performance and can satisfy the requirements of design index.Table 2
Tested PSNR results of our compression.
Images
PSNR (dB)
Average PSNR (dB)
Bit rate (bits/pixel)
SPOT-1
42.99
42.53
0.25
SPOT-2
43.70
SPOT-3
40.90
SPOT-1
47.07
47.06
0.5
SPOT-2
48.86
SPOT-3
45.26
SPOT-1
50.95
51.27
1
SPOT-2
51.93
SPOT-3
50.93
SPOT-1
53.64
53.74
2
SPOT-2
54.75
SPOT-3
52.82In the second experimental, in order to compare compressed performance of our approach, we use several proposed compression algorithms, AT-3DSPIHT [33], 3D-DWT [34], 3D-PCA [35], 3D-SPECK [36], and SA-DCT [37], to compress remote sensing MSCFBs. We use the AVIRIS multiband remote sensing images in JPL laboratory. Each group is 256
pixels
×
288
pixels
×
4; the depth of every pixel is 8 bits. The images are Low Altitude, Lunar Lake, Jasper Ridge, and Cuprite. The compression bit rate is set to 2.0~0.25 bpp. Figure 14 shows the comparison results of average PSNR for different compression methods. From Figure 14, due to the full usage of posttransform and TD, the proposed PTTD-based compression algorithm has the best compression performance. The PTTD-based method can improve 0.3~1.3 dB of average PSNR compared with other compression methods at 2.0~0.25 bpp.The test and comparison result of multiband images: (a) is PSNR difference between 3D-DWT and other three compression codecs, and (b) is PSNR difference between 3D-PCA and other two compression codecs.
(a)
(b)Overall, the PTTD-based compression method has excellent compression performance. It can fully meet the requirement of space multispectral CCD camera.


## 4.4. The Spectral Distortion Analysis

We use the max similarity of spectrum (MSS) to evaluate the distortion of multispectral images for our algorithm. LetX be original remote sensing MSCFBs and let Y be reconstructed remote sensing MSCFBs. The MSS can be expressed as
(23)
h
MSS
=
max
⁡
i
,
j
f
(
i
,
j
)
2
+
(
1
-
c
(
i
,
j
)
2
)
2
,
f
(
i
,
j
)
=
1
n
λ
∑
λ
(
X
(
i
,
j
,
λ
)
-
Y
(
i
,
j
,
λ
)
)
2
,
c
(
i
,
j
)
=
∑
λ
(
X
(
i
,
j
,
λ
)
-
μ
i
,
j
)
(
Y
(
i
,
j
,
λ
)
-
μ
i
,
j
′
)
(
n
λ
-
1
)
δ
i
,
j
δ
i
,
j
′
,
where i, j, and λ are line, column, and band number of multispectral image, respectively, n
λ is the total number of bands, and μ
i
,
j and δ
i
,
j are the average value and mean square of the corresponding pixel in the spectral dimension, respectively.The greater the value of MSS is, the greater the distortion of spectra is. The MSS between original and reconstructed multispectral images for the image group using different approach is shown as in Figure15.Figure 15
Comparison of MSS performance.From Figure15, our approach has the lowest MSS, so the proposed algorithm can protect spectral information well.Figure16 demonstrates the reconstructed third band of two remote sensing images, which is 1024 × 1024; the depth of every pixel is 8 bits. The bit rate is set to 0.5 bpp.Figure 16
Comparison of reconstruction performance.In order to view the detailed information of reconstructed image, the part of image is100
×
100. From Figure 14, under the low bit rate, 3DSPIHT-based approach causes the visible ringing effects and the blurred edge. However, our approach does not cause the blurred edge. So, the proposed algorithm can protect spectral information well.


## 4.5. Proposed Algorithm Complexity Analysis and Compression Time

In our method, 1D-DCT is performed on spectral dimension to exploit the spectral information. It is reported that DCT has the lower amount of computations than DWT. When one of the most efficient methods is used to perform DCT, 8 × 8 points of 2D-DCT use only 462 time adding operations and 54 time multiplication operators. This means that each point of 2D-DCT needs only 0.84 time multiplication operators and 7.22 adding operations. Therefore,N bands of MSCFBs performing 1D-DCT need N
×
0.84 multiplications and N
×
7.22 additions.After applying 1D-DCT, three-level 2D-DWT is applied to the spatial bands. LetL be decomposition level of 2D-DWT. For the I
1
×
I
2
×
I
3 MSCFBs performing 2D-DWT, the calculation complexity is O
(
8
N
I
1
I
2
I
3
(
1
-
2
-
2
L
)
/
6
). In the PTTD-based method, a 9/7 2D-DWT is used to transform each band of MSCFBs. And we perform three levels of decomposition. Therefore, the calculation complexity of the PTTD-based method is O
(
9
I
1
I
2
I
3
/
7
).After MSCFBs performing 2D-DWT, we use the Hadamard posttransform in low bit rate. Let each DWT coefficients block haveM coefficients. The Hadamard-based posttransform needs the O (M
log
⁡
M) operations. In high bit rate, we use the DCT posttransform. 2D 4 × 4 DCT requires only 27 multiplications and 231 additions. In addition, the best posttransform selection uses only 2
×
(
M
-
1
) time adding operations and one comparison operator. Therefore, the complexity of posttransform stage is M
(
log
⁡
2
M
+
3
)
+
2 operations for each block with M coefficients.After 2D-DWT coefficients performing posttransform, the dimension of each subtensor is4
×
4
×
I
3. Let I
A be the average dimensions of core tensor G. The complexity of the calculation of tensor decomposition is O
(
J
A
I
A
3
), where
(24)
I
A
=
1
3
(
4
+
4
+
I
3
)
,
J
A
=
1
3
(
J
1
A
+
J
2
A
+
J
3
A
)
.
Therefore, the calculation complexity of TD in the PTTD-based method is O
(
J
A
I
A
3
)
×
∑
L
=
1
3
(
I
1
/
2
L
+
2
×
I
2
/
2
L
+
2
).In the following, we test and analyze the calculation times of the proposed PTTD-based compression method. Note that, in this test experiment, they only evaluate the calculation time of the PTTD-based compression because it is not the best optimized implementation on FPGA-based compression system. The evaluations testing performs the lossy compression for an image of size 3072 × 128 at 1.0 bpp on FPGA EVM board with system clock frequency at 88 MHz. Table3 demonstrates the comparison of complexity between the proposed multispectral compression algorithm and others.Table 3
The results of complexity comparison.
Methods
Times
Ours
0.042 us/sample
KTL [38]
0.102 us/sample
3DSPIHT [3]
0.062 us/sample
JPEG2000 [39]
0.181 us/sampleFrom Table3, the processing time of our algorithm reaches 0.042 us/sample; the data throughput is 23.81MSPS, which is higher than JPEG2000, KLT, and 3DSPIHT, so our approach has the lowest complexity. In our project, the multispectral CCD camera works at an orbits altitude of 500 km, the scroll angle is -
40
°~+
40
°, the latitudes are -
70
°~+
70
°, and the line working frequency is 1.8094 kHz~0.85946 kHz. In the line working frequency, capturing the image of 128 × 3072 required 70.74 ms. Using our approach, compressing the image of 128
×
3072 requires 16.51 ms. So, our approach can process the four bands images simultaneously. This meets the requirement of the project.In addition, we use the XC2V6000-6FF1152 FPGA to implement the proposed algorithm. The design language is VerilogHDL, the development platform is ISE8.2, and synthesis tool is XST. Table4 demonstrates the occupancy of resources of our approach.Table 4
The occupancy of resources.
Resources
Utilization rate
Slices
23655/33792 (70%)
Four input LUTs
45281/67584 (67%)
BRAM
116/144 (80%)From Table4, the LUTs occupy 67%, slices occupy 70%, and BRAM occupy 80%. Various indicators are lower than 95%, which meet the requirement of our project.


## 5. Conclusion

In this paper, we propose an efficient compression approach for MSCFBs. In our approach, the 1D-DCT is performed on spectral dimension to exploit the spectral information, and the PT in 2D-DWT domain is performed on each spectral band to exploit the spatial information. A deep coupling approach between the PT and TD is proposed to remove residual spectral redundancy between bands and residual spatial redundancy of each band. Experimental results show that the PTTD-based method has the average PSNR reaching upwards of 40 dB. And the PTTD-based method can improve 0.3~1.3 dB of average PSNR compared to other compression methods at 2.0~0.25 bpp. The processing time of our algorithm reaches 0.042 us/sample; the data throughput is 23.81MSPS, which is higher than JPEG2000, KLT, and 3DSPIHT. Therefore, the proposed PTTD-based compression algorithm has the excellent compression performance and the low calculation complexity. It is very suitable for the space MSCFBs compression.



---

*Source: 296474-2014-08-20.xml*