An Entropy-Guided Residual Clustering Framework for Lossless Compression of Hyperspectral Images

1. Introduction

Hyperspectral remote sensing has become an essential technology for observing, measuring, and interpreting complex physical, chemical, and biological phenomena on the Earth’s surface. Unlike conventional RGB imaging, which records only three broad visible channels, hyperspectral sensors acquire hundreds or even thousands of narrow and contiguous spectral bands, typically spanning wavelengths from approximately 400 nm to 2500 nm. This rich spectral representation enables the identification of subtle material properties that are difficult or impossible to distinguish in standard images. As a result, hyperspectral imagery has been widely used in surveillance [1], mineral exploration [2], industrial inspection [3], precision agriculture [4], and environmental monitoring [5]. However, the same spectral richness that makes hyperspectral data valuable also produces extremely large data volumes. Efficient compression is therefore necessary to reduce storage requirements, accelerate transmission, and support practical deployment in airborne, satellite, and field-based remote sensing systems.

Compression methods for hyperspectral imagery are generally classified into lossy, near-lossless, and lossless approaches. Lossy compression can achieve high compression ratios by discarding information deemed visually or statistically less important, but it is irreversible. Although lossy methods may suppress sensor noise, they may also remove weak but scientifically meaningful spectral features or introduce compression artifacts. Such distortions can reduce the reliability of subsequent analysis, especially when the data are used for classification, target detection, anomaly identification, or quantitative retrieval. Near-lossless methods seek a balance between compression efficiency and bounded reconstruction error and have been investigated through techniques such as the Karhunen–Loève transform [6], principal component analysis, vector quantization, and JPEG2000-based coding [7]. These methods can improve rate–distortion performance and computational efficiency, but they still do not guarantee exact recovery of the original data.

Lossless compression, by contrast, preserves the original image without any numerical alteration. Every pixel value and every spectral band can be reconstructed exactly after decompression. This property is particularly important for scientific remote sensing, where raw observations may need to be reanalyzed using future algorithms, validated against field measurements, or archived as reference data. Lossless preservation is also critical in domains such as medical hyperspectral imaging [8] and material characterization [9], where small spectral variations may carry diagnostic or structural significance. Although lossless compression usually yields lower compression ratios than lossy alternatives, its ability to maintain complete data integrity makes it indispensable for applications where information loss is unacceptable.

Traditional lossless image compression is built on several major principles: entropy coding, dictionary-based coding, predictive coding, and transform-based decorrelation. Classical methods such as run-length encoding, Huffman coding, arithmetic coding, Lempel–Ziv family algorithms, and prediction by partial matching have strongly influenced many widely used compression standards and file formats [10]. Dictionary-based approaches, including LZ77, LZ78, and LZW, have been adopted in practical systems such as ZIP, PNG, and GZIP because of their simplicity and efficiency [11]. Entropy coding methods remain important because they assign shorter codewords to frequently occurring symbols, thereby reducing the average number of bits required for representation. However, hyperspectral images introduce additional challenges because redundancy exists not only within spatial neighborhoods but also across adjacent spectral bands.

In remote sensing, several lossless compression algorithms have been developed specifically for multispectral and hyperspectral imagery, including clustering vector quantization, predictive trees [12], differential pulse-code modulation methods [13], and other spatial–spectral prediction schemes [14]. Arnavut et al. [15] proposed a method in which one band is compressed using lossless JPEG. In contrast, subsequent bands are reordered based on the sorting permutation of previous bands and then encoded using move-to-front coding to reduce entropy. The CCSDS 123.0-B standard [16] is among the most widely recognized real-time lossless hyperspectral compression standards. It operates by first generating predictions of incoming samples and then encoding the differences between predicted and actual values. This approach has been reported to achieve compression ratios of approximately 50–60%. Álvarez-Cortés et al. [17,18] combined Regression Wavelet Analysis with the CCSDS-123.0-B-1 predictor and a lightweight contextual arithmetic coder, achieving notable coding gains over RWA combined with JPEG2000 [19]. Other studies have shown that integral wavelet transforms [20] and clustering-based methods can effectively reduce both spatial and structural redundancies, thereby improving the compression of multispectral images [21].

Predictive coding has become one of the most influential strategies for lossless hyperspectral image compression. In this approach, a pixel or spectral sample is estimated from previously encoded neighboring samples, and only the prediction residual is encoded. Since residuals usually have lower entropy than the original data, they can be represented more compactly. Context-dependent predictors further improve this process by adapting the prediction model to local image characteristics. Ulacha and Łazoryszczak [22] proposed a context-dependent linear prediction method based on mean absolute error minimization and reported better compression performance than conventional predictors. Rajput et al. [23] introduced content-aware prediction mechanisms that dynamically adjust to image characteristics, making predictive compression more suitable for efficient data sharing and distributed storage environments.

More recently, machine learning has opened new directions for lossless compression. Instead of relying solely on manually designed predictors or fixed statistical models, neural-network-based methods can learn complex dependencies directly from training data. Rhee et al. [24] proposed a duplex neural network architecture that predicts both pixel values and contextual information for lossless image compression, achieving performance improvements over traditional methods. Sun et al. [25,26] further surveyed and benchmarked neural universal compressors, showing that deep learning models are increasingly competitive for compressing diverse data sources. These approaches suggest a broader transition from conventional statistical compression toward adaptive and intelligent compression systems.

For hyperspectral and multispectral remote sensing data, compression performance depends heavily on how effectively spatial, spectral, and structural redundancies are removed. Recent studies have therefore focused on spectral decorrelation, predictive modeling, transform coding, and hybrid compression frameworks. Altamimi and Ben Youssef [27] reviewed lossless and near-lossless hyperspectral compression techniques and emphasized the importance of maintaining spectral fidelity. Nagendran et al. [28] proposed a lossless framework that combines spectral decorrelation with transform coding to reduce inter-band redundancy. These studies demonstrate that exploiting the strong correlation among neighboring spectral bands is essential for improving compression efficiency without sacrificing reconstruction accuracy.

The entropy-based lossless compression algorithms are based on the improved K-means algorithm (IKA), which efficiently improves the K-means clustering algorithm [21] by accounting for the spatial geometric relationships among remotely sensed pixels. The best class number can be predicted from entropy analysis of the probability distribution models for the clustering lossless compression algorithm (CLCA). The predicted best class number closely matches the experimental results, where the class number sweeps from 1 to a large value to locate the class number that has the highest lossless compression ratio (LCR).

This article is a revised and expanded version of a paper entitled “Entropy Analysis for Clustering Based Lossless Compression of Remotely Sensed Images” [29], which was presented virtually at 2021 IEEE Big Data Conference on Dec. 18, 2021.

2. Methods

2.1. Improve K-Means Algorithm (IKA)

In the K-means algorithm [21], every pixel calculates the Euclidean distance with all clustering centers to determine its assignment, so clustering a large number of classes is a relatively time-consuming task. Actually, classes far from the pixel do not need to be analyzed, and stable classes do not need re-clustering. To create one more class, the class with the most irregular spatial distribution or the largest standard deviation can be divided according to its spreading direction. These strategies can help the clustering converge faster.

In the improved K-means algorithm (IKA) (shown in Figure 1), it is supposed that N_s is the total number of pixels attending clustering; b is the number of bands; n_c is the number of classes; C_i is the set of pixels in the class i; N(C_i) is the number of pixels in class i; c_i is the center of class i; c_i^l is the value of lth band for c_i; c_ilast^l is the lth band value of c_i in the previous iteration; e_i^jis the jth pixel in class i (e_i^j∈C_i); e_i^jl is the value of lth band for pixel e_i^j; C_ij is the set of classes that need to make re-clustering judgment with pixel e_i^j; d_i^j is the distance between pixel e_i^j∈C_i and clustering center c_i; d_ij is the distance from class i to class j; d_i^jk is the distance from the pixel e_i^jto the center of class k; array Pd[i] =max({d_i^j|e_i^j∈C_i}) (i∈[0,n_c]) is a collection of the maximum deviate distances for the n_c classes, and P_i={e_i^m |d_i^m=max({d_i^j|e_i^j∈C_i})) corresponds to the pixel e_i^m with the maximum of d_i^j in the class i; E(d_i) is the mean of {d_i^j|e_i^j∈C_i}; lastsign[i], used to judge inter-class stability, is the modification mark of class i in the previous iteration; sign[i], used to judge clustering convergence, is the modification mark of class i in the current iteration; ε_c is the maximal iterative error of clustering centers.

1

N_c initial clustering centers c_i (i∈[0, N_c-1]) are created according to the initial spatial division.

(a)

If the previous clustering record with the class number N_c′ exists, the previous clustering environment, including class centers, pixel status, and other important information, is rebuilt, setting n_c=N_c′. Otherwise, all pixels are grouped into a single class C₀, with n_c=1, meaning the one-class space has N_s pixels. The clustering center c₀ is calculated as (1) with i=0, and the initial clustering status, lastsign, is assigned with lastsign [0]=1. The distances d₀^j(j∈[0, N_s-1]) (from (2) with i=0) between the intra-class pixels e₀^jand their center c₀ are sorted in descending order, and the maximal value of d₀^j(j∈[0, N_s-1]) is saved to Pd [0] corresponding to the farthest pixel P₀ to the center c₀.

$${\mathit{c}}_i^l={\displaystyle \frac{1}{N\left({\mathit{C}}_i\right)}}{\sum }_{j=1}^{N\left({\mathit{C}}_i\right)}{\mathit{e}}_i^{jl},1<l<b,{\mathit{e}}_i^j\in {\mathit{C}}_i$$

(1)

$$d_i^j={\left[\sum _{l=1}^b\left({\mathit{e}}_i^{jl}-{\mathit{c}}_i^l\right)\right]}^{1/2}$$

(2)

(b)

The pixel e_i^m in the class C_i that is the farthest from the center c_i is noted as P_i, and this maximal distance for class i is recorded in Pd[i] (i∈ [ 0, n_c-1]). Along the direction $\stackrel{⃑}{P_i{C}_i}$ = c_i -e_i^m in (3), this class is divided according to the vertical plane passing through the middle point between the farthest pixel e_i^m and the point that has a distance of E(d_i) with c_i along the inverse direction of $\stackrel{⃑}{P_i{C}_i}$, creating two new classes C_nc and C_i′ as (4) (see Figure 2a). C_i′ will take the place of C_i, letting C_i=C_i′. Their sizes are N(C_nc) and N(C_i), respectively. The two new classes are updated, including the following values: for an instance of the new class C_i, the class center c_i (1), the distance d_i^j between the intra-class pixel e_i^j ∈ C_i and its center (2), E(d_i), Pd[i], and P_i.

$$\mathit{{\rm I}}={\displaystyle \frac{\overrightarrow{P_i{C}_i}}{\left|\overrightarrow{P_i{C}_i}\right|}}$$

(3)

$$\begin{gathered} \forall {\mathit{e}}_i^j\in {\mathit{C}}_i,j\ne m,\mathrm{if}{d}_i^j<{dm}_i-{\displaystyle \frac{d_i^m+E\left(d_i\right)}{2}}\Rightarrow {{\mathit{C}}_i}^{\text{'}}={{\mathit{C}}_i}^{\text{'}}\cup \left\{{\mathit{e}}_i^j\right\} \\ \mathrm{else}\mathrm{if}\left|\left({\mathit{e}}_i^j-{\mathit{e}}_i^m\right)•\mathit{I}\right|<{\displaystyle \frac{d_i^m+E\left(d_i\right)}{2}}\Rightarrow {\mathit{C}}_{n_c}={\mathit{C}}_{n_c}\cup \left\{{\mathit{e}}_i^j\right\}\mathrm{else}{{\mathbf{C}}_i}^{\text{'}}={{\mathit{C}}_i}^{\text{'}}\cup \left\{{\mathit{e}}_i^j\right\} \end{gathered}$$

(4)

Let lastsign[i]=lastsign[n_c]=1 and n_c=n_c+1. When the number of classes is reached, i.e., n_c==N_c, terminate the class division stage and proceed to step 2; otherwise, proceed to step 1b and continue the division process.

2.

The clustering centers before division or in the previous iteration are recorded in c_ilast, letting c_{ilast =}c_i (i∈[0, N_c - 1]), and the class modification mark in the current iteration, sign, is set as sign[i]=0 (i∈[0, N_c-1]).

3.

The distance matrix D={d_ij, i,j∈[0, N_c - 1]} is created by calculating the Euclidean distance as (5) between any two class centers.

$$d_{ij}={\left[{\sum }_{l=1}^b\left({\mathit{c}}_i^l-{\mathit{c}}_j^l\right)\right]}^{1/2},i\ne j,d_{ij}\in \mathit{D}$$

(5)

4.

For each pixel e_i^j∈C_iin class i∈[0, N_c - 1], d_ij(j≠i, j∈[0, N_c - 1]) are sorted in ascending sequence to create one array C_ij (for k∈C_ij, exists k≠i, k∈[0, N_c - 1], d_i^j>0.5d_ik), defining the search order for pixel e_i^j, where the class with smaller distance to the class i is in the preceding position (see Figure 2b). The class centers with a distance to the class center c_i larger than 2d_i^j are not under consideration.

5.

For each neighbor class k in C_ij, if lastsign[i]==0 and lastsign[k]==0 (k∈C_ij), meaning the relationship between classes j and k is stable, k is removed from C_ij because it is not necessary to make Euclidean analysis between them. Otherwise, the pixel e_i^j must be re-clustered: the distance d_i^jk between e_i^j and the center c_k is calculated. For the class l≠i with the minimal value {l|d_i^jl=min{{d_i^jk|k∈C_ij},d_i^j}} (6), pixel e_i^j is combined with the class l whose center is nearest to it, modifying class sets: C_i=C_i-{e_i^j}, C_l=C_l+{e_i^j}, and setting the class modification marks as sign[i]=1 and sign[l]=1.

$$d_i^{jk}={\left[{\sum }_{l=1}^b{\left({\mathit{e}}_i^{jl}-{\mathit{c}}_k^l\right)}^2\right]}^{1/2},k\in {\mathit{C}}_{ij}$$

(6)

6.

For each class i∈[0, N_c-1] whose sign[i]==1, c_i and d_i^j (e_i^j∈C_i) are updated according to (5) and (6). If |c^l_ilast-c^l_i|< ε_c (l∈[1, b]), class i is assumed stable, setting sign[i]=0.

7.

If none of the classes are changed, sign[i]==0 (i∈[0, N_c-1]), or the iteration number reaches the limit, the clustering process ends and proceeds to step 8; otherwise, it proceeds to step 2.

8.

The previous clustering status lastsign is updated to lastsign=sign to record class stability in the current iteration.

9.

Current clustering results are saved, including the class structure and all the variables.

Initial clustering continues working on the previous clustering results and divides the class with the maximum deviation value in P[]. This way not only saves initial clustering time but also preserves the stable relationships between classes formed by the clustering process.

In this algorithm, only classes whose correlations show changes will be re-clustered; intra-class pixels merely compute distances to neighboring class centers. So, during the iterative process, the statuses of more classes are stable, leading to higher clustering velocity.

2.2. Process of Clustering Lossless Compression Algorithm

Based on the fast IKA clustering for high class numbers, the clustering lossless compression algorithm (CLCA) could ultimately reduce intra-class redundancy. However, statistical Huffman coding could utilize this redundancy to improve the lossless compression ratio (LCA). Since the relationship between neighboring bands is strongest, clustering centers could eliminate spectral redundancy between them.

Supposed that m, n, b are respectively the length, width, and band number of remotely sensed images, X_{j, k}ⁱ is the kth band data of the jth sample in the ith class, I_kⁱ is the kth band data of the ith class center, N(i) is pixel number of the ith class, s_ij is the class serial number of pixel (i, j) (spatial coordinates). Then the detailed process of CLCA is provided as follows:

1

For remotely sensed images, apply IKA clustering, creating N_c class centers I_i, i∈[0, N_c-1] (choose the intra-class pixel nearest to the original computed clustering center as the center of this class, this has less error compared with directly being rounded to integer), and classification map S={s_ij∈[0, N_c - 1],i∈[0, m-1],j∈[0, n-1]}. And apply statistical Huffman coding to S at the same time.

2

Code clustering centers:

(a)

Compute spectral residue of each center:

$$d{I}_k^i=I_k^i-I_{\mathrm{k}-1}^i,i\in [0,N_c-1],k\in [1,b-1]$$

(7)

(b)

Compute the residue of the first band data between all classes:

$$d{I}_0^i=I_0^i-I_0^{\mathrm{i}-1},i\in [1,N_c-1]$$

(8)

(c)

Apply Huffman statistical coding to the residue data $dI$, and $I_0^0$ is saved separately.

3.

Code intra-class pixels:

(a)

Compute residue between each intra-class pixel and its center:

$$d{X}_{\mathrm{j},\mathrm{k}}^i=X_{\mathrm{j},\mathrm{k}}^i-I_k^i,i\in [0,N_c-1],j\in [0,N\left(i\right)-1],k\in [0,b-1]$$

(9)

(b)

Apply statistical Huffman coding to the residue $dX$

The preceding two steps of this algorithm, 1) and 2), could realize VQ compression, whose compression ratio (CR) could reach hundreds of times. Using the preceding clustering results, we could quickly determine CR for different class numbers to identify the best class number.

2.3. Probability Model Analysis of CLCA

Assuming that the number of pixels is N_s and the number of classes is N_c, the probability distribution model of the residue data for CLCA is analyzed as follows.

• Intra-class residue:

When the number of classes reaches a pre-defined value N_c, the intra-class residue is a set X={0, ±1, ±2, …}, the element of which is a variable x∈X approximately taking on a Gaussian distribution $N(0,{\sigma }_1^2)$. Maximum entropy with a limited mean and variance is provided in (10) [6].

$$H_1={\mathit{log}}_2\left({\sigma }_1\sqrt{2\pi e}\right)$$

(10)

Then, the total number of bits for the intra-class residue is given by (11).

$$B_1=N_{s}b{H}_1=N_{s}b{\mathit{log}}_2\left({\sigma }_1\sqrt{2\pi e}\right)$$

(11)

2.

Class center residue:

Supposed that A=(f_max - f_min) (f_max and f_min respectively represent the maximal and minimal values of all bands for all class centers), if the number of bands is high enough, the class center residue data x∈[-A, A] approximately take on a Gauss distribution $N(0,{\sigma }_2^2)$ (Figure 3a), whose maximum entropy with a limited mean and deviance is (12).

$$H_2={\mathit{log}}_2\left({\sigma }_2\sqrt{2\pi e}\right).$$

(12)

Then the total number of bits for the class center residue is (13).

$$B_2=N_{c}b{H}_2=N_{c}b{\mathit{log}}_2\left({\sigma }_2\sqrt{2\pi e}\right)$$

(13)

3.

Classification map:

The classification map data X={0, 1, 2, …, N_c} take on a uniform distribution with probability density function p₃(x)=1/N_c, whose entropy is H₃=log₂N_c, so the total number of coded bits is

B₃=N_sH₃=N_slog₂N_c

(14)

4.

The total number of bits in CLCA is given by (14).

$$B=B_1+B_2+B_3=N_{s}b{log}_2\left({\sigma }_1\sqrt{2\pi e}\right)+N_{c}b{log}_2\left({\sigma }_2\sqrt{2\pi e}\right)+N_s{log}_2\left(N_c\right)$$

(15)

B₃ has a logarithmic relationship with the number of classes N_c, so the change of N_c has little influence on B₃. When the class number N_c is small, σ1 changes noticeably, decreasing as Nc increases. So, B depends on B₁ and decreases as Nc increases. When the class number N_c is high enough, B increases linearly with the increase in Nc.

2.4. Least Entropy Lossless Compression Algorithm

The Least Entropy Lossless Compression Algorithm (LELCA) calculates the class number with the fewest compression bits for CLCA. Suppose that if N_c increases m times, (N_c, σ₁) turns into (N_c^′, σ₁^′) described in equations (16), where k and x are to be determined.

$$\left\{\begin{array}{c}{{\sigma }_1}^{\text{'}}={\displaystyle \frac{1}{k{m}^x}}{\sigma }_1\\ {N_c}^{\text{'}}=m{N}_c\end{array}\right.$$

(16)

From (16), the class number N_c and the intra-class pixel standard deviation σ₁ are related as described in Equation (17).

$$x\left({\mathit{ln}{N}_c}^{\text{'}}-\mathit{ln}{N}_c\right)=\mathit{ln}{\sigma }_1-{\mathit{ln}{\sigma }_1}^{\text{'}}-\mathit{ln}k$$

(17)

Equation (17) can also be expressed as Equation (18).

$$xd\left(\mathit{ln}{N}_c\right)=-d\left(\mathit{ln}{\sigma }_1\right)-\mathit{ln}k$$

(18)

N_c can be derived as Equation (19).

$$N_c={\displaystyle \frac{k_2}{{\left({\sigma }_1{k}^{{\sigma }_1}\right)}^{1/x}}}$$

(19)

There are three unknown parameters, x, k, and k₂, which can be calculated using three initial conditions, as described by Equation (20).

$$\left\{\begin{array}{c}{N}_c={\displaystyle \frac{k_2}{{\left({\sigma }_1{k}^{{\sigma }_1}\right)}^{1/x}}}\\ N_c{|}_{{\sigma }_1={\sigma }_{1n_1}}=n_1,N_c{|}_{{\sigma }_1={\sigma }_{1n_2}}=n_2,N_c{|}_{{\sigma }_1={\sigma }_{1n_3}}=n_3\end{array}\right.$$

(20)

Solving the equations (20), we have equations (21) and (22).

$$\mathit{ln}{N}_c-\mathit{ln}{n}_1=\left(\mathit{ln}{n}_2-\mathit{ln}{n}_1\right){\displaystyle \frac{\mathit{ln}{\sigma }_{1n_1}-\mathit{ln}{\sigma }_1+\left({\sigma }_{1n_1}-{\sigma }_1\right)\mathit{ln}k}{\mathit{ln}{\sigma }_{1n_1}-\mathit{ln}{\sigma }_{1n_2}+\left({\sigma }_{1n_1}-{\sigma }_{1n_2}\right)\mathit{ln}k}}$$

(21)

and,

$$\mathit{ln}k={\displaystyle \frac{\left(\mathit{ln}{n}_3-\mathit{ln}{n}_1\right)\left(\mathit{ln}{\sigma }_{1n_1}-\mathit{ln}{\sigma }_{1n_2}\right)-\left(\mathit{ln}{n}_2-\mathit{ln}{n}_1\right)\left(\mathit{ln}{\sigma }_{1n_1}-\mathit{ln}{\sigma }_{1n_3}\right)}{-\left(\mathit{ln}{n}_3-\mathit{ln}{n}_1\right)\left({\sigma }_{1n_1}-{\sigma }_{1n_2}\right)+\left(\mathit{ln}{n}_2-\mathit{ln}{n}_1\right)\left({\sigma }_{1n_1}-{\sigma }_{1n_3}\right)}}$$

(22)

Then we can derive Equation (23).

$${\displaystyle \frac{d{\sigma }_1}{d{N}_c}}=-{\displaystyle \frac{\mathit{ln}{\sigma }_{1n_1}-\mathit{ln}{\sigma }_{1n_2}+\left({\sigma }_{1n_1}-{\sigma }_{1n_2}\right)\mathit{ln}k}{N_c\left(\mathit{ln}{n}_2-\mathit{ln}{n}_1\right)}}\times {\displaystyle \frac{{\sigma }_1}{1+{\sigma }_1\mathit{ln}k}}$$

(23)

We can set dB/dN_c = 0 to obtain the optimal class number that yields the highest compression ratio. From Equation (15),

$${\displaystyle \frac{dB}{d{N}_c}}={\displaystyle \frac{N_{s}b}{{\sigma }_1\mathit{ln}2}}{\displaystyle \frac{d{\sigma }_1}{d{N}_c}}+b{\mathit{log}}_2\left({\sigma }_2\sqrt{2\mathsf{\pi }\mathrm{e}}\right)+{\displaystyle \frac{N_s}{N_c\mathit{ln}2}}=0$$

(23)

The best number of classes can be predicted by Equation (24), where σ₂ is the standard deviation of clustering centers’ residue.

$$N_{\mathrm{cbest}}={\displaystyle \frac{N_{s}b\left[\mathit{ln}{\sigma }_{1n1}-\mathit{ln}{\sigma }_{1n2}+\left({\sigma }_{1n1}-{\sigma }_{1n2}\right)\times \mathit{ln}k\right]-N_s\left(\mathit{ln}{n}_2-\mathit{ln}{n}_1\right)\times \left(1+{\sigma }_1\times \mathit{ln}k\right)}{b\mathit{ln}\left({\sigma }_2\sqrt{2\mathsf{\pi }\mathrm{e}}\right)\times \left(\mathit{ln}{n}_2-\mathit{ln}{n}_1\right)\times \left(1+{\sigma }_1\times \mathit{ln}k\right)}}$$

(24)

The LELCA algorithm is proposed as follows based on Equation (24) to obtain the best class number and the corresponding best LCR.

• Choose initial three class numbers: n₁, n₂, n₃.
• Cluster the space into n₁ classes and use CLCA to get σ_1n1, σ_2n1, and Huffman coding length L₁; cluster the space into n₂ classes and use CLCA to get σ_1n2, σ_2n2, and Huffman coding length L₂; cluster the space into n₃ classes and use CLCA to get σ_1n3, σ_2n3, and Huffman coding length L₃.
• According to the prediction Equation (24), assuming σ₁ = σ_1n3 and σ₂ = σ_2n3, N_cbest is calculated. If N_cbest is in the list of the initial and predicted class numbers, go to step 7 to terminate.
• Cluster the space into N_cbest classes and use CLCA to get its corresponding σ_1best, σ_2best, and code length L_best.
• Update n₁, n₂, n₃: (1) sort all the initial and predicted class numbers from low value to high value as an array; (2) search the array from left to right; the class number with the lowest code length is set as n₃; (3) the previous two class numbers on the left in the array are set as n₁ and n₂.
• Go to step 2.
• In the list of the initial and predicted class numbers, class number n with the smallest Huffman coding length L is the best.

2.5. Conditions of Ultimate Division of Classes

When continuing dividing the intra-class pixels which shows a Gaussian distribution when class number is small, more intra-class pixels take on a spatially uniform distribution (see Figure 3b), noting that intra-class residue x∈[0, R] does not take on a uniform distribution). According to the Equation (54) in Appendix A, in b-dimensional space V_b, the probability distribution is p₁(x) in (25), with entropy H₁ in (26).

$$p_1\left(x\right)={\displaystyle \frac{b{x}^{b-1}}{R^b}}$$

(25)

$$H_1=-{\displaystyle \underset{0}{\overset{R}{\int }}}{\displaystyle \frac{b{x}^{b-1}}{R^b}}\log \left({\displaystyle \frac{b{x}^{b-1}}{R^b}}\right)dx={\displaystyle \frac{1}{\ln 2}}\left(\ln R+{\displaystyle \frac{b-1}{b}}-\ln b\right)$$

(26)

If R≈σ₁, H₁ is similar in form to the H₁ in (10). When the class number reaches N_c0 (here, σ₁=σ₁₀), if σ₁ is reduced m times, then σ₁^′=σ₁/m, N_c0^′=m^bN_c0. The class number N_c and the standard deviation of intra-class pixels have the relation of Equation (27).

$${\left({\displaystyle \frac{{N_{c_0}}^{\text{'}}}{N_{c_0}}}\right)}^{1/b}=m={\displaystyle \frac{{\sigma }_{10}}{{{\sigma }_{10}}^{\text{'}}}}$$

(27)

The (27) can be further expressed as (28) along with an initial condition.

$$\left\{\begin{array}{c}{\displaystyle \frac{1}{b}}\mathrm{d}\mathrm{l}\mathrm{n}{\mathrm{N}}_c=-d\mathit{ln}{\sigma }_1\\ N_c{|}_{{\sigma }_1={\sigma }_{10}}=N_{c_0}\end{array}\right.$$

(28)

Then

$$N_c={\displaystyle \frac{k}{{\sigma }_1^b}}$$

(29)

According to initial condition, (29) can be rewritten as

$$k=N_{c_0}{{\sigma }_{1_0}}^b$$

(30)

Then

$$N_c={\displaystyle \frac{N_{c_0}{{\sigma }_{1_0}}^b}{{\sigma }_1^b}}$$

(31)

Then

$$B_1+B_3=\left[N_s{\mathit{log}}_2{\sigma }_1^b+N_s\left(b-1-b\mathit{ln}b\right)/\mathit{ln}2\right]+N_s{\mathit{log}}_2{N}_c=N_s{\mathit{log}}_2\left(N_{c_0}{{\sigma }_{1_0}}^b\right)+N_s\left(b-1-b\mathit{ln}b\right)/\mathit{ln}2$$

(32)

which is a constant. Then

$$B=\left[N_s{\mathit{log}}_2\left(N_{c_0}{{\sigma }_{1_0}}^b\right)+N_s\left(b-1-b\mathit{ln}b\right)/\mathit{ln}2\right]+N_{c}b{\mathit{log}}_2\left({\sigma }_2\sqrt{2\pi e}\right)$$

(33)

Here, total entropy B is determined by B₂. When one class is added, B increases by $b{\mathit{log}}_2\left({\sigma }_2\sqrt{2\mathsf{\pi }\mathrm{e}}\right)$ (bits).

3. Results

Sook Lake 220-band AVIRIS images (256×256, 16-bit, obtained from [29], Figure 4a) were selected for algorithm analysis. The lossless compression ratio (LCR) is defined as (34).

$$LCR{\displaystyle \frac{=Originalsize}{/Losslesscompressedsize}}$$

(34)

3.1. IKA

In the IKA clustering algorithm, with ε_c=0.01 and N_c=100, it converges after 214 iterations; LCR is 2.3075, and the average iteration time is 8 seconds. Figure 4b is the gray balanced result of the rebuilt 60th band image using the IKA clustering results of 100 class centers, where the intra-class pixels are present uniformly using their class center. This is a near-lossless compression method.

Suppose N_c=200. Then the initial clustering further divides the space based on the clustering results of N_c=100, taking only 3s, converging after 210 iterations. Each iteration takes an average of 9s, and LCR reaches 2.3607. So, an increased class number has little influence on clustering iteration speed; however, traditional K-means needs 6 minutes (100 classes) for one iteration, and clustering time increases abruptly with the increased class number.

3.2. CLCA

The clustering a class number makes full use of historical clustering results, enabling lossless compression of each class number at high velocity. Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 shows the clustering results for class numbers from 1 to 1780, with a maximum of 20 iterations, and ε_c is set to 0.02. Figure 5 shows that intra-class residue B₁ has an inverse tendency to the class number, because σ_{1 decreases} with the increase of class number. The more classes in space, the smaller the occupied area of a class. Therefore, Figure 5 meets the supposition of the probability model in Equation (11).

In Figure 6, the clustering centers’ residue B₂ shows an approximately linear relationship with N_c, consistent with the assumptions of the probability model in Equation (13).

In Figure 7, the classification map coding B₃ has a logarithmic relationship with the number of classes, N_c, as given by the Equation (14), so when N_c is a large number, the change of N_c has little influence on B₃.

When the number of classes N_c increases, σ₁ and B₁ decrease (Figure 5), but B₂ and B₃ increase (Figure 6 and Figure 7). From Figure 8, the total length of the code and codebook reaches a minimum at N_c=1092, where the LCR reaches a maximum of 2.4152 (Figure 9). When the class number is less than 1092, the reduced amount of B₁ is larger than the increased amount of B₂+B₃, leading to an increased LCR; when the class number is higher than 1092, the reduced amount of B₁ is less than the increased amount of B₂+B₃.

3.3. Estimation of the Best Compression of CLCA

The least-entropy lossless compression algorithm (LELCA) is used to estimate the optimal class number N_c that yields the best compression with CLCA. In Table 1, the left table shows the clustering results for different class numbers, from the input and the predicted N_cbest using (24); the right table shows the execution of LELCA, where the number of iterations for clustering a class number has no limit, and ε_c is set to 0.0001. The execution process for LELCA is provided below and in Table 1.

• Setting n₁=200, n₂=300, n₃=400 as input parameters, their σ_1n1=35.407818, σ_1n2=33.676117, and σ_1n1=32.581664. According to Equation (24), N_cbest=1023, whose σ₁= 29.458689, σ₂ =232.382148, L_best=11959639 and LCR is 2.411096.
• Updating n₁=300, n₂=400, n₃=1023, according to Equation (24), N_cbest=934, whose σ₁= 29.763009, σ₂ =233.284102, L_best=11966192 and LCR is 2.409776.
• Updating n₁=400, n₂=934, n₃=1023, according to Equation (24), N_cbest=1040, whose σ₁= 29.402593, σ₂ =232.143746, L_best=11958672 and LCR is 2.411291.
• Updating n₁=934, n₂=1023, n₃=1040, according to Equation (24), N_cbest=1064, whose σ₁= 29.360868, σ₂ =231.907104, L_best=11962256 and LCR is 2.410569.
• Updating n₁=934, n₂=1023, n₃=1040, according to Equation (24), N_cbest=1064, which has been clustered in step 4. The LELCA converges and stops the iteration.

Therefore, the predicted best class number is 1040, with an LCR of 2.411291. The predicted best class number, 1040, is different from the best class number, 1092, obtained from sweeping the range of class numbers (Figure 8 and Figure 9). The reasons include (1) the prediction only considers σ₁ and σ₂ for calculating N_cbest, (2) the maximal iteration number in CLCA introduces noise in the sweep analysis of the range of class numbers. More information from the code length will enhance the accuracy of LELCA.

The iteration number decreases with the increased class numbers until reaching the best class number, 1040, whose iteration number is only 35; while the iteration number for the class number 200 is 199. When the class number is higher than the best class number, 1040, the iteration number increases.

3.4. Comparison with DPCM/D2PCM Algorithms

Two lossless compression algorithms, CLCA and DPCM, have a coding discrepancy of ΔB (equation 35), where B’ is from (58) in the Appendix B and B is from (15). Under normal conditions, σ₁<σ₁^′, i.e., ΔB>0; when the number of bands b is extremely high, likely, σ₁>σ₁^′, but current remote sensing equipment cannot meet this requirement. So, the CLCA algorithm is more efficient than the DPCM algorithm.

$$∆B=B^{\text{'}}-B\approx N_s\left(b-1\right){\mathit{log}}_2\left({{\sigma }_1}^{\text{'}}\sqrt{2\mathsf{\pi }\mathrm{e}}\right)-N_{s}b{\mathit{log}}_2\left({\sigma }_1\sqrt{2\mathsf{\pi }\mathrm{e}}\right)\approx N_{s}b{\mathit{log}}_2\left({\displaystyle \frac{{{\sigma }_1}^{\text{'}}}{{\sigma }_1}}\right)$$

(35)

When using DPCM, σ₁^′=177, B^′=14,048,912 bytes, and the lossless compression ratio is 2.0525, so CLCA is more efficient than DPCM. However, for D²PCM, σ₁^′ is 190.957, and LCR is 1.9963. Therefore, spectral second-order residue cannot increase LCR. Table 2 shows that LELCA yields the best class number and highest LCR but is more time-consuming because it uses an iterative method to determine the optimal class number.

5. Discussion

The experimental results demonstrate that entropy-guided clustering provides an effective mechanism for removing both spectral and spatial redundancy in hyperspectral images while maintaining exact reconstruction. Unlike conventional predictive coding approaches, the proposed framework exploits pixel similarities via clustering and represents both class centers and intra-class deviations using residual coding. This strategy transforms highly correlated hyperspectral data into lower-entropy representations that can be efficiently encoded using Huffman coding. The results obtained from the AVIRIS Sook Lake dataset indicate that clustering-based residual representation provides a practical alternative to traditional predictive lossless compression techniques.

The improved K-means algorithm (IKA) plays a crucial role in enabling efficient clustering with a large number of classes. Conventional K-means requires distance calculations between each pixel and all cluster centers, leading to rapidly increasing computational complexity as the number of clusters grows. By leveraging prior clustering results, preserving stable inter-class relationships, and restricting pixel reassignment to neighboring classes, the proposed IKA substantially reduces unnecessary computations while maintaining convergence speed, even for large class counts. The results for 100- and 200-class problems demonstrate that increasing the number of clusters does not significantly increase the iteration time, making the algorithm suitable for high-dimensional hyperspectral datasets.

The experimentally observed behavior of the compression components agrees well with the theoretical probability models derived in Section 2. As the number of classes increases, the intra-class standard deviation decreases, resulting in a reduction of the intra-class residual coding length. At the same time, the coding lengths associated with class-center residuals and classification maps increase. These opposing trends explain the existence of an optimal class number. The minimum total coding length observed near 1092 classes confirms the entropy analysis and validates the assumptions underlying the CLCA probability model. Therefore, the proposed entropy framework provides not only a compression method but also a theoretical explanation of the relationship between the number of clusters and coding efficiency.

The least-entropy lossless compression algorithm (LELCA) further demonstrates that exhaustive searches over all class numbers are unnecessary. Based on the entropy model and a limited number of initial clustering results, LELCA predicted an optimal class number of 1040 and achieved a compression ratio of 2.4113, which is very close to the maximum compression ratio of 2.4152 obtained by exhaustive analysis. The small discrepancy between the predicted and experimental optimum can be attributed to simplifying assumptions in the entropy model and the limited number of iterations used during sweeping analysis. Nevertheless, the prediction accuracy indicates that entropy-based parameter estimation provides a practical means of determining near-optimal clustering configurations with significantly reduced computational cost.

Compared with DPCM and D2PCM, the proposed clustering-based framework achieved noticeably higher compression ratios. DPCM and D2PCM primarily exploit spectral correlations between neighboring bands, whereas CLCA simultaneously utilizes spatial and spectral redundancy through clustering and residual coding. Experimental results show that DPCM achieved a lossless compression ratio of 2.0525, D2PCM obtained 1.9963, and the proposed method achieved compression ratios exceeding 2.41. These results suggest that clustering-based residual representation is more effective than first- and second-order spectral differencing for the AVIRIS dataset considered in this study. Furthermore, the degradation observed in D2PCM indicates that higher-order spectral differencing does not necessarily improve compression performance, as it may increase entropy rather than reduce it.

Although the proposed approach provides improved compression performance, several limitations should be noted. First, the entropy model assumes approximately Gaussian distributions for residual data, which may not accurately describe all remotely sensed images. Second, Huffman coding is used as the entropy coder; more sophisticated coders, such as arithmetic coding or context-adaptive coding, may provide additional coding gains. Third, the experimental evaluation was performed on a single AVIRIS dataset, and broader validation using images from different sensors and resolutions would strengthen the generality of the conclusions. In addition, the computational complexity associated with clustering large-scale hyperspectral datasets may become significant for onboard or real-time applications.

Overall, the results demonstrate that residual clustering and entropy analysis provide a unified framework for lossless hyperspectral image compression. The agreement between theoretical predictions and experimental observations confirms the validity of the proposed probability models and the effectiveness of entropy-guided cluster optimization. Future research may focus on adaptive clustering strategies, context-based entropy coding, parallel implementations, and machine-learning-assisted prediction models. Combining entropy-guided clustering with deep neural compression techniques may further improve compression efficiency and enable practical deployment in next-generation hyperspectral remote sensing systems.

6. Conclusions

Hyperspectral remote sensing images contain substantial spatial and spectral redundancy, making efficient lossless compression essential for long-term storage, transmission, and scientific analysis. This study presented a residual clustering-based lossless compression framework that combines an improved K-means clustering algorithm (IKA), a clustering-based lossless compression algorithm (CLCA), and a least-entropy lossless compression algorithm (LELCA) to determine the optimal number of clusters. By exploiting prior clustering results, preserving stable inter-class relationships, and limiting distance calculations to neighboring classes, the proposed IKA significantly accelerates convergence and enables efficient clustering for large numbers of classes.

Based on the clustering results, CLCA enhances intra-class redundancy and removes spectral correlation through residual coding and Huffman entropy coding. Probability distribution models were developed for the intra-class residuals, class-center residuals, and classification map, providing theoretical insight into the relationship between entropy and the number of clusters. Building on these models, LELCA predicts the class number that minimizes the total coding length without requiring exhaustive searches.

Experiments using the 220-band AVIRIS Sook Lake hyperspectral dataset demonstrated that the proposed approach effectively balances compression efficiency and computational complexity. The maximum lossless compression ratio obtained by CLCA was 2.4152 at 1092 classes, while LELCA predicted an optimal class number of 1040 and achieved a compression ratio of 2.4113, showing close agreement with exhaustive analysis. Compared with DPCM and D2PCM, whose compression ratios were 2.0525 and 1.9963, respectively, the proposed clustering-based method provided superior compression performance while preserving exact reconstruction of the original data.

The results indicate that entropy-guided clustering is an effective mechanism for reducing spatial and spectral redundancy in hyperspectral images. In addition to improving compression efficiency, the proposed framework provides a practical means of estimating optimal clustering parameters from entropy characteristics. Future work will investigate adaptive clustering strategies, more advanced entropy coders such as arithmetic coding and context-based coding, and deep-learning-assisted prediction models to improve compression efficiency further and extend the framework to larger hyperspectral datasets and real-time onboard remote sensing applications.

Appendix A.1. The Density Function of Equally Distributed Space

Pixel is

• 2D space

For an equally distributed two-dimensional circular space with maximal radius R at pixel density $\sigma$ (Figure 3b), for the annulus between radius $x\left(x\in \left[0,R\right]\right)$ and $x+dx$, the pixel number is $\left[\pi {\left(x+dx\right)}^2-\pi x^2\right]\sigma =\left[{\left(dx\right)}^2+2xdx\right]\pi \sigma$. The total number of pixels in that circular field is $\pi R^2\sigma$. Then we have one Equation (36).

$${\int }_0^R\pi \sigma \left[{\left(dx\right)}^2+2xdx\right]=\pi R^2\sigma$$

(36)

Or

$${\int }_0^R\left(\pi \sigma dx\right)dx+{\int }_0^{R}2\pi \sigma xdx=\pi R^2\sigma$$

(37)

As dx is very small, we can assume (38).

$${\int }_0^R\left(\pi \sigma dx\right)dx\approx 0$$

(38)

So, the number of pixels at position x is $2x\pi \sigma$, and the density function is given by (39).

$$f\left(x\right)={\displaystyle \frac{2x\pi \sigma }{\pi R^2\sigma }}={\displaystyle \frac{2x}{R^2}},\left(x\in \left[0,R\right]\right)$$

(39)

2.

3D space

For an equally distributed sphere space with maximal radius R at pixel density $\sigma$ (Figure 3b), for the annulus between radius $x\left(x\in \left[0,R\right]\right)$ and $x+dx$, the pixel number is $\left[{\displaystyle \frac{4}{3}}\pi {\left(x+dx\right)}^3-{\displaystyle \frac{4}{3}}\pi x^3\right]\sigma =\left[{\left(dx\right)}^3+3x{\left(dx\right)}^2+3x^{2}dx\right]{\displaystyle \frac{4}{3}}\pi \sigma$. The total number of pixels in that circular field is ${\displaystyle \frac{4}{3}}\pi R^3\sigma$. Then we have Equation (40).

$${\int }_0^R\left[{\left(dx\right)}^3+3x{\left(dx\right)}^2+3x^{2}dx\right]{\displaystyle \frac{4}{3}}\pi \sigma ={\displaystyle \frac{4}{3}}\pi R^3\sigma$$

(40)

Or,

$${\int }_0^R{\displaystyle \frac{4}{3}}\pi \sigma {\left(dx\right)}^{2}dx+{\int }_0^R\left(4\pi \sigma xdx\right)dx+{\int }_0^{R}4\pi \sigma x^{2}dx={\displaystyle \frac{4}{3}}\pi R^3\sigma$$

(41)

As dx is very small, we can assume (38).

$${\int }_0^R{\displaystyle \frac{4}{3}}\pi \sigma {\left(dx\right)}^{2}dx+{\int }_0^R\left(4\pi \sigma xdx\right)dx\approx 0$$

(42)

So, the number of pixels at position x is $4x^2\pi \sigma$, and the density function is (43).

$$f\left(x\right)={\displaystyle \frac{4x^2\pi \sigma }{\frac{4}{3}\pi R^3\sigma }}={\displaystyle \frac{3x^2}{R^3}},\left(x\in \left[0,R\right]\right)$$

(43)

3.

The n-dimensional Space

Following the same method, the density function in n-dimensional space is given by (44).

$$f_n\left(x\right)={\displaystyle \frac{n{x}^{n-1}}{R^n}}$$

(44)

Assume that in n-dimensional space, there are n independent variables, x₁, …, x_n, described by (45) and (46).

$$\mathfrak{R}:x_1^2+x_1^2+\dots +x_n^2\le R^2$$

(45)

$$\left\{\begin{array}{c}{x}_1=r\mathit{cos}{\varphi }_1\\ x_2=r\mathit{sin}{\varphi }_1\mathit{cos}{\varphi }_2\\ \dots \dots \dots \dots \\ x_{n-1}=r\mathit{sin}{\varphi }_1\mathit{sin}{\varphi }_2\dots \mathit{sin}{\varphi }_{n-2}\mathit{cos}{\varphi }_{n-1}\\ x_n=r\mathit{sin}{\varphi }_1\mathit{sin}{\varphi }_2\dots \mathit{sin}{\varphi }_{n-2}\mathit{sin}{\varphi }_{n-1}\end{array}\right.$$

(46)

We can get (47).

$$J={\displaystyle \frac{\partial \left(x_1,x_2,\dots ,x_n\right)}{\partial \left(r,{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_n\right)}}=r^{n-1}{\mathit{sin}}^{n-2}{\varphi }_1{\mathit{sin}}^{n-3}{\varphi }_2\dots \mathit{sin}{\varphi }_{n-2}$$

(47)

Then we can get the volume in (48).

$$V=\int \dots \int d{x}_{1}d{x}_2\dots d{x}_n=\int \dots \int r^{n-1}{\mathit{sin}}^{n-2}{\varphi }_1{\mathit{sin}}^{n-3}{\varphi }_2\dots \mathit{sin}{\varphi }_{n-2}d{\varphi }_{1}d{\varphi }_2\dots d{\varphi }_{n-2}dr={\int }_0^R{r}^{n-1}dr\left[\int \dots \int {\mathit{sin}}^{n-2}{\varphi }_1{\mathit{sin}}^{n-3}{\varphi }_2\dots \mathit{sin}{\varphi }_{n-2}d{\varphi }_{1}d{\varphi }_2\dots d{\varphi }_{n-2}\right]={\displaystyle \frac{R^n}{n}}A$$

(48)

where

$$A=\int \dots \int {\mathit{sin}}^{n-2}{\varphi }_1{\mathit{sin}}^{n-3}{\varphi }_2\dots \mathit{sin}{\varphi }_{n-2}d{\varphi }_{1}d{\varphi }_2\dots d{\varphi }_{n-2}$$

(49)

For an equally distributed n-sphere space with maximal radius R and pixel density $\sigma$, in the annulus between radius $x\left(x\in \left[0,R\right]\right)$ and $x+dx$, the pixel number is (50).

$$\left[{\displaystyle \frac{{\left(x+dx\right)}^n}{n}}A-{\displaystyle \frac{x^n}{n}}A\right]\sigma =\left[n{x}^{n-1}dx+C_n^2{x}^{n-2}{\left(dx\right)}^2+\dots +{\left(dx\right)}^n\right]{\displaystyle \frac{A\sigma }{n}}$$

(50)

According to (48), the total number of pixels in that circular field is ${\displaystyle \frac{R^n}{n}}A\sigma$. Then we have Equation (51).

$${\int }_0^R\left[n{x}^{n-1}dx+C_n^2{x}^{n-2}{\left(dx\right)}^2+\dots +{\left(dx\right)}^n\right]{\displaystyle \frac{A\sigma }{n}}={\displaystyle \frac{R^n}{n}}A\sigma$$

(51)

Or,

$${\int }_0^{R}n{x}^{n-1}dx+{\int }_0^R{C}_n^2{x}^{n-2}{\left(dx\right)}^2+\dots +{\int }_0^R{\left(dx\right)}^n=R^n$$

(52)

As dx is very small, we can assume (53).

$${\int }_0^R{C}_n^2{x}^{n-2}{\left(dx\right)}^2+\dots +{\int }_0^R{\left(dx\right)}^n\approx 0$$

(53)

So, the pixel number at position x is $A\sigma x^{n-1}$, and the density function is (54).

$$f\left(x\right)={\displaystyle \frac{A\sigma x^{n-1}}{\frac{R^n}{n}A\sigma }}={\displaystyle \frac{n{x}^{n-1}}{R^n}},\left(x\in \left[0,R\right]\right)$$

(54)

Appendix A.2. Probability Distribution Model of DPCM/D2PCM Algorithms

According to [4], the detailed process of the DPCM lossless compression algorithm for hyperspectral images is as follows.

• Calculate spectral residue of each pixel according to (55).

$$d{I}_k=I_k-I_{\mathrm{k}-1},k\in [1,b-1]$$

(55)

2.

Apply S+P integral wavelet transform [30] to the 1st band data, deleting spatial redundancy.

3.

Apply Huffman statistical coding to residue dI_k (k∈[1, b-1]) and solely save I(0,0).

D²PCM calculates the second-time residue of spectral data and makes amendments to the preceding steps: step 1 is added by $d^2{I}_k=d{I}_k-d{I}_{\mathrm{k}-1},k\in [2,b-1]$, and in step 3, the residue is updated as dI₁, d²I_k (k∈[2, b-1]).

The probability distribution model for DPCM/D2PCM can be derived as follows.

1.

Spectral residue x∈[-A, A] of each pixel takes on a Gaussian distribution N(0, σ₁^′2) [5] (see Figure 2c), whose entropy is ${H_1}^{\text{'}}={\mathit{log}}_2\left({{\sigma }_1}^{\text{'}}\sqrt{2\mathsf{\pi }\mathrm{e}}\right)$ , so the total bit number is (56).

$${B_1}^{\text{'}}=N_{s}b{H_1}^{\text{'}}=N_{s}b{\mathit{log}}_2\left({{\sigma }_1}^{\text{'}}\sqrt{2\mathsf{\pi }\mathrm{e}}\right)$$

(56)

2.

The 1st band residue isX={0, ±1, ±2, …}; the variable x∈Xtakes on a Gaussian distribution N(0, σ₂^′2); then the total number of bits for the intra-class pixels’ residue is (57).

$${B_2}^{\text{'}}=N_{s}b{H_2}^{\text{'}}=N_{s}b{\mathit{log}}_2\left({{\sigma }_2}^{\text{'}}\sqrt{2\mathsf{\pi }\mathrm{e}}\right)$$

(57)

3.

The total bit number of DPCM/ D²PCM algorithm is (58).

$$B^{\text{'}}={B_1}^{\text{'}}+{B_2}^{\text{'}}=N_{s}b{\mathit{log}}_2\left({{\sigma }_1}^{\text{'}}\sqrt{2\mathsf{\pi }\mathrm{e}}\right)+N_{s}b{\mathit{log}}_2\left({{\sigma }_2}^{\text{'}}\sqrt{2\mathsf{\pi }\mathrm{e}}\right)$$

(58)