Progressive Sample Processing of Band Selection for Hyperspectral Image Transmission

Band selection (BS) is one of the important topics in hyperspectral image (HSI) processing. Many types of BS algorithms were proposed in the last decade. However, most of them were designed for off-line use. They can only be used with pre-collected data, and are sometimes ineffective for applications that require timeliness, such as disaster prevention or target detection. This paper proposes an online BS method that allows us obtain instant BS results in a progressive manner during HSI data transmission, which is carried out under band-interleaved-by-sample/pixel (BIS/BIP) format. Such a revolutionary method is called progressive sample processing of band selection (PSP-BS). In PSP-BS, BS can be done recursively pixel by pixel, so that the instantaneous BS can be achieved without waiting for all the pixels of an image. To develop a PSP-BS algorithm, we proposed PSP-OMPBS, which adopted the recursive version of a self-sparse regression BS method (OMPBS) as a native algorithm. The experiments conducted on two real hyperspectral images demonstrate that PSP-OMPBS can progressively output the BS with extremely low computing time. In addition, the convergence of BS results during transmission can be further accelerated by using a pre-defined pixel transmission sequence. Such a significant advantage not only allows BS to be done in a real-time manner for the future satellite data downlink, but also determines the BS results in advance, without waiting to receive every pixel of an image.


Introduction
Due to the use of hundreds of spectral bands, hyperspectral imaging (HSI) generally has enormous data volume and contains vast amount of information.This special characteristic results in several issues.First, the inter-band correlation of HSI is very high, and adjacent bands may contain redundant spectral information.This would lead to the well-known problem called the "curse of dimensionality" in data analysis.Second, the computational complexity of processing HSI data is very high.Third, storing HSI data usually requires large amount of disc space.Finally, transmitting HSI data requires higher bandwidths.Under such circumstances, removing partial data without significant loss of an image's spectral information is necessary.One commonly used approach is band selection (BS) [1].BS takes advantage of such high-band correlation to remove the redundant bands, in order to achieve a wide range of applications, such as dimensionality reduction, data storage, data transmission, target detection, and classification.
Many different types of BS algorithms [2][3][4][5][6][7][8][9][10][11][12][13][14][15] have been proposed in the past two decades.Most of them make an assumption that the BS problem is an optimization problem, which maximizes or minimizes a pre-defined objective function that can measure the amount of information or inter-band redundancy contained by the currently selected bands.For instance, Keshava et al. [2] adopts common distance measures, such as Euclidean distance and a spectral angle mapper, to measure the similarity of bands.Du et al. [3] uses linear regression and orthogonal subspace projection to sequentially select a new band based on the complementary vector space.Chang et al. [4] regards a spectral band as a signal, and designed a BS method based on a constrained energy minimization algorithm.Martínez-Usó et al. [5] introduces a hierarchical clustering-based method to select bands from the band groups that were pre-clustered by a specific method, with some informational measures.In addition, there has been a lot of research that adopts the perspective of information theory.The works in [6][7][8] adopt mutual information as the similarity index to find the maximum information spanned by the bands.In recent years, sparse regression models were used to perform BS.For instance, Sun et al. [9][10][11] uses a self-sparse regression (SSR) model to select bands.In an SSR problem, finding a new basis is equivalent to selecting the most representation bands in an HSI image.Lai et al. [12] proposed a SSR-based BS method, which adopts an orthogonal matching pursuit (OMP) algorithm to sequentially find the next band in an efficient manner.For the purpose of HSI data transmission, Chang et al. [13,14] proposed a progressive BS (PBS) method for spectral unmixing.In PBS, the bands are first prioritized by a certain criterion, and are transmitted by the order of prioritization.Du et al. [15] proposed a BS-based dimensional reduction method for change detection in multi-temporal hyperspectral images.Except for the above-mentioned work, there were still various kinds of BS research published in the literature [16][17][18][19][20][21][22][23][24][25].
Since many existing HSI algorithms were designed to deal with off-line problems, they usually require more computing time for pursuing extremely high accuracy.Undoubtedly, those sophisticated methods may not be appropriate for dealing with timeliness applications.Hence, developing real-time and efficient approaches is urgently needed.Based on that point, some research started to develop on-board computing approaches [26][27][28][29], where the data compression could be carried out on satellites to reduce the demand of downlink bandwidth, focusing additionally on real-time image classification with the aid of GPU acceleration [30,31].
On the other hand, processing HSI data in a real-time progressive fashion under the data transmission held from satellite to ground station has become another attractive topic.The term "progressive processing" for HSI data originated from [32].According to [32], the term refers to a type of signal processing method.It decomposes entire data processing into a finite number of stages, and processes data progressively stage by stage, in the sense that the results obtained by previous stages can be used to update or improve the results to be processed in subsequent stages.Recently, the concept of progressive processing was utilized to process HSI data.For example, in order to achieve real-time spectral unmixing during band transmission, a revolutionary concept called progressive band processing (PBP) [33], was developed.PBP assumes that the HSI data is transmitted by band sequential (BSQ) format, and processes the image immediately after receiving a new band.Under the PBP framework, many common topics, such as spectral unmixing, target detection, and anomaly detection could be realized [33][34][35].Unlike the traditional algorithms, which can only be implemented on the collected data, PBP makes use of recursive processing, to instantly process the current data segment.More specifically, PBP methods preserve the past information obtained at previous stage, and use it to accelerate the computation at the current stage.The most attractive feature of PBP is that the computing time will not significantly increase as the growth of number of received bands.Therefore, PBP algorithms have significant potential for the uses of applications in data transmission and communication.
Based on the importance of BS, and the urgent demand of real-time algorithms, this paper proposes a real-time progressive BS method for applications related to data transmission.To our best knowledge, there were no BS methods proposed for this purpose before.In our case, we assume that HSI data is transmitted in band-interleaved-by-pixel (BIP) or band-interleaved-by-sample (BIS) format, and is processed pixel by pixel (sample by sample).According to [32], such a process operated on BIP/BIS transmission is called progressive sample processing (PSP).To sum up, we developed a novel concept called progressive sample processing of band selection (PSP-BS), which combines BS with PSP. Figure 1 illustrates the difference of traditional BS and PSP-BS.
To realize PSP-BS, we had to choose an appropriate BS algorithm as the core of PSP-BS.In this paper, a self-sparse regression BS algorithm, called orthogonal matching pursuit-based BS (OMPBS) [12], has been adopted because (1) the OMPBS can sequentially find the most information-complementary bands from the image, without any prior knowledge or any complicated optimization scheme; and (2) The algorithm is easily re-formulated to fulfill recursive processing through mathematical decomposition of the optimization equations of OMPBS.In the later sections, we first introduce OMPBS, and then derive its recursive version, called recursive OMPBS (Re-OMPBS), as the key algorithm to realize PSP-BS.In addition, we make an assumption that the pixel transmission sequence (PTS) could be defined before transmission.We propose three PTSs to make BS results converge more quickly during the PSP-BS.We call Re-OMPBS combined with PTS PSP-OMPBS.
To evaluate the effectiveness of PSP-OMPBS, two real hyperspectral datasets were used in the experiments.We conducted a comparison of computing time and the corresponding BS accuracy at every time stage in a progressive manner.Land cover classification was also adopted to evaluate the quality of the selected bands.According to the experimental results, the computational efficiency of PSP-OMPBS is stably high so that the process could be carried out in real time.The BS accuracy can be improved early on by using particular PTSs.

Orthogonal Matching Pursuit-Based BS (OMPBS)
The OMPBS [12] is a sequential BS method based on a self-sparse regression (SSR) model [36].It aims to find a set of representative bands that can represent all bands based on the minimization of reconstructed errors in the SSR model, by using an orthogonal matching pursuit (OMP) algorithm [37].Suppose L is the number of bands and N is number of pixels in a hyperspectral image.Let b i ∈ R N×1 presents the ith band vector.In an SSR model, both the observation matrix as well as the dictionary matrix are the set of band vectors denoted by B = [b 1 , b 2 , . . . ,b L ] and the coefficient matrix is denoted by C. Obviously, the goal is to find the optimal coefficient matrix C that minimizes the reconstructed error ||B − BC|| 2  F .Since only k bands need to be selected, we impose the sparsity constraint ||C|| 0,2 p where l 0 /l 2 norm ||C|| 0,2 counts the number of the non-zero rows in C. Thus, our sparse-BS optimization problem is where ||B − BC|| 2 F is the reconstructed error, ||C|| 0,2 p is the sparsity constraint, and p is the sparsity level.The sub-optimization problem (find C) is solved by least square estimator Q = P T P −1 P T B, where P is temporary basis matrix composed by the candidate bands.The original OMPBS algorithm is composed of a doubly-nested loop.So, it is required to solve L!/(L − p)! times of (1), so the overall computation is time-consuming.To improve its efficiency, we revised the original algorithm by residual perspective.The revised algorithm (Algorithm 1) is shown as follows.
2. At j-th iteration, calculate the row norm vector k ∈ L × 1 of matrix R T j R j .3. Find the row index b j whose value is maximum at k. Then set

Solve the sub optimization problem
F by Q j = (P T j P j ) −1 P T j B. 5. Update the residual matrix R j+1 = B − P j Q j .6.If j = p, break; otherwise, set j ← j + 1 and go to Step 2.

Progressive Sample Processing-Based OMPBS (PSP-OMPBS)
The proposed PSP-OMPBS method is fully introduced in this section.It consists of two parts, as shown in Figure 2. The first part controls the data transmission order, in which the pixel transmission sequence will be formed.The second part is the real-time BS computing core, recursive OMPBS (Re-OMPBS).Section 3.1 introduces the derivation of Re-OMPBS, where we convert some steps in OMPBS to their recursive forms.Section 3.2 describes the condition for implementing recursive equations introduced in Section 3.1.Section 3.3 summaries the details of the Re-OMPBS algorithm.Finally, the formation of the pixel transmission sequence is introduced in Section 3.4.The center part is the proposed PSP-OMPBS method which includes two blocks: formation of the pixel transmission sequence (PTS), and implementation of the recursive OMPBS (Re-OMPBS) algorithm.The transmission is carried out from the PTS block to the Re-OMPBS block.

Derivation of Recursive OMPBS (Re-OMPBS)
In this section, we aim to integrate OMPBS into PSP.According to the definition, PSP-BS can instantly acquire the current BS result during the data transmission by using two types of information: the processed information provided in previous stage (i.e., the (n − 1)th stage), and the new information provided in current stage (i.e., nth stage).More specifically, PSP-BS can immediately update the BS result when a new pixel (the nth pixel) is received.Assume that n presents the number of currently received pixels.We would like to utilize the past information obtained in the (n − 1)th stage to reduce the computational burden of the nth stage.In doing so, the optimization process of OMPBS must be broken down.

Decomposition of the Least Square Estimator
In the following derivation, we assume that j is a fixed number and do not annotate it for each term.In the OMPBS algorithm, there are two calculations that can be decomposed.The first one is the least square estimator in Step 4: ( As we know, computing Equation ( 2) is required for each iteration j in the inner loop of OMPBS.This computation would be time-consuming when n is large.This would be unfavorable to progressive processing.To resolve this issue, we must convert Equation (2) to its recursive version.To do that, we define two matrices: B(n) represents the total pixel data after the nth pixel is received, and P(n) presents the data's temporary basis matrix, composed by selected band vectors.Let r n = [b n1 , b n2 , . . . ,b nL ] be the new pixel and p n = p n1 , p n2 , . . ., p np .Next, we define auxiliary matrix H(n) Substituting ( 4) into (2), we get

Now use Woodberry's Identity [38]
A to decompose Equation (5).Let A = P(n − 1) T P(n − 1), c = p T n , v = A −1 and ρ = c T A −1 c, and Equation (6) becomes where . We expand Equation (7), and further simplify it by 1+ρ n|(n−1) 1+ρ n|(n−1) Multiply the 3rd term of Equation ( 8) by H(n − 1)H(n − 1) −1 , and then re-organize it.Finally we obtain 1+ρ n|(n−1) 1+ρ n|(n−1) Equation ( 9) is the recursive version of Equation ( 2).In Equations ( 7)-( 9), Q(n − 1) and H(n − 1) are the past information obtained at the (n − 1)th stage; v n|(n−1) and ρ n|(n−1) are the so-called innovation information [32], whose calculation is involved with both previous information and current information; and p n and r n are provided by the new incoming pixel.Once Q(n) is obtained, we calculate H(n) = H(n − 1) + p T n p n for the use in the next stage.It should be noted that Equation ( 9) can only be used when the updating condition is met.This detail will be presented in Section 3.2.

Decomposition of the Residual Multiplication Term
On the other hand, in OMPBS, the other calculation that can be simplified is the residual Step 5 of OMPBS.Calculating this term would be time-consuming when n becomes large.According to OMPBS, the residual term is formed by We expand R(n) T j+1 R(n) j+1 in terms of Equation ( 10): There are four terms required to be computed in Equation (11).The first term is B(n) T B(n), which can be obtained by the following recursive equation Similarly, the second term can be represented by where Q j (n) T is supposed to be known by Equation ( 9).The third term B(n is the transposition of the second term, so it is unnecessary to re-compute it.Finally, the fourth term can be calculated by In conclusion, term R(n) T j+1 R(n) j+1 can be quickly calculated by using Equations ( 11)-( 14).

Condition for Applying Recursive Equations
Unlike related work about PBP [33][34][35], where the recursive calculation can be applied to each iteration of PBP, Equation ( 9) of Re-OMPBS can only be implemented under particular conditions.To understand it, we assume that Ω j (n) presents the BS result at the jth iteration in the nth stage, and the result was just obtained.Now we move to calculate Q j (n) for finding the residual.If the condition Ω j (n) = Ω j (n − 1) holds, we can use Equation ( 9) to calculate Q j (n).More specifically, if the selected bands' indices of the jth iteration of the previous stage are the same as the band indices of the jth iteration in the current stage, we can use the recursive equation to save computing time.This is simply because the sparse coefficient matrix Q j (n) can only be "updated" when the currently-selected bases (i.e., band indices) are exactly the same as the previous ones.More specifically, the information in the horizontal dimension of temporary matrices P j (n − 1) and P j (n) are aligned so that the information of previous stage can be shared with the current stage.In the following, we call it the "updating condition".
Satisfying the updating condition not only reduces computing time for Q j (n), but also slightly accelerates Equations ( 12)-( 14).One key of implementing Re-OMPBS is that we have to store the computing results obtained at the (n − 1)th stage, such as terms , for applying in Equation ( 9) and Equations ( 11)-( 14), when the updating condition is held at the nth stage.One should note that meeting this condition is not necessary for Re-OMPBS.For those cases that do not meet the conditions, we directly use the One might think that the probability of reaching these conditions would be low.However, this is not the case.In fact, when the amount of received pixels n reaches a certain amount, a single new incoming pixel will not obviously influence the BS result.In other words, the BS result is expected to be consistent.This phenomenon will be demonstrated in our experiment section.Moreover, Re-OMPBS is regarded as the "approximate" version of OMPBS.Theoretically, there would be some tiny numerical errors accumulated in Q j (n) if the updating condition has been reached continuously for a long period.In this case, an occasional lack of meeting the updating condition exactly gives the chance to reset the numerical error to zero.

Re-OMPBS Algorithm
Based on the above-mentioned content, the Re-OMPBS algorithm (Algorithm 2) is listed below.

Design of Pixel Transmission Sequence (PTS)
The original pixel transmission order in BIP/BIS format is based on a horizontal line scan [35].If we regard a hyperspectral image as a 2D matrix, the pixel transmission starts from the upper left corner.Once the first row has been transmitted, then it transmits the second row, and so on.However, using this sequence for Re-OMPBS might mean that the BS cannot converge until the almost all pixels are transmitted.Ideally, in a progressive processing, we might expect to see the final results as soon as possible.In this case, adopting an appropriate sampling technique in the spatial domain to select the representative pixels for early transmission is necessary.There are three pixel transmission sequences we considered in this paper.

Original Band-Interleaved-by-Sample/Pixel (BIS/BIP) Sequence
As mentioned above, the original pixel transmission sequence is the same as the storage sequence in the BIS/BIP format.Suppose the size of a hyperspectral image cube I is x × y × L. The number of total pixels is N = x × y.We reshape I to an L × N matrix form in row-major order, [r 1 r 2 . . .r N ] L×N , where the column vectors are spectral pixels of I.The BIS/BIP transmission sequence is defined by {r 1 , r 2 , . . . ,r N }.A simple illustration is shown in Figure 3a.

Step Sequence
The step method refers to uniformly selecting pixels in an original sequence with interval k.
To do this, we uniformly segment pixel set [r 1 r 2 . . .r N ] into k groups.For each transmission, one pixel is selected from each group in turn.If k = 5, the resulting transmission pixel sequence would be {r 1 , r 6 , r 11, . . .r 2 , r 7 , r 12... }, as shown in Figure 3b.As a result, the transmitted pixels are the samples uniformly picked from the pixel set.Theoretically, using such sequence would make the PSP-BS converge more quickly, because in this case fewer received pixels can represent the image better.

Block Sequence
Since the step method only considers image sampling in the horizontal direction, we further consider another method, which can consider both vertical and horizontal directions simultaneously, called the block method.The block method first splits the image into several b × b square blocks, as shown in Figure 3c.Then the new pixel is picked from each block in turn.In such a design, the sampling process is performed on the full spatial domain.Theoretically, it would shorten convergence time of PSP-BS, too.In the later sections, we call the PSP-OMPBS implemented with original BIS/BIP sequence, step sequence, and block sequence PSP-OMPBS, S-PSP-OMPBS, and B-PSP-OMPBS, respectively.
Finally, we summarize the distinct features of OMPBS and the proposed PSP-OMPBS in Table 1.

Table 1.
The basic properties of orthogonal matching pursuit-based BS (OMPBS) and the proposed PSP-OMPBS.

Description
A traditional BS method which is implemented on a pre-obtained image data.
A PSP-BS method.It is used in the moment of spectral sample transmission or collecting.

Use of past information No Yes
Use of PTS No Yes (Original/Step/Block) Least square estimator used in the optimization Residual term R T j R j calculation Direct computation Recursive Equations ( 11)-( 14)

Capable for progressive processing
No, the computing time will increase with the data volume.Yes

Experiments
In this section, the proposed PSP-OMPBS is tested on two publicly-available real hyperspectral images.To evaluate the performance of PSP-OMPBS, we conduct three different studies, including BS accuracy analysis, land cover classification, and computational efficiency.According to the introduction, the development of PSP-BS is not for achieving outstanding BS for a particular data analysis, such as image classification or spectral unmixing, but for the real-time BS monitoring during data transmission.In this case, the experiment is conducted based on self-comparison, instead of comparing with state-of-the-art methods.To evaluate the BS accuracy and efficiency of PSP-OMPBS, we still can adopt OMPBS as the compared method (baseline).
In the following progressive experiments, we set n = 100-207,400 for the Pavia data, and n = 100-21,025 for the Purdue data.For the p value, we set 16 for the Pavia data and 34 for the Purdue data, according to the virtual dimensionality (VD) algorithm [39] with false alarm 0.01.Three PTS methods mentioned in Section 3.4 were adopted for PSP-OMPBS methods.We empirically set k as 200 and 100 and b as 50 and 25 for the Pavia and Purdue data, respectively.Table 2 lists the parameters used in the experiment.

BS Accuracy Analysis
The accuracy of the instant BS result during data transmission is an index to evaluate the effectiveness of PSP-OMPBS.Theoretically, the BS would gradually converge to the final BS result (i.e., the result of OMPBS performed on the complete image cube) over time.To observe this phenomenon, we use the results of OMPBS implemented on both images, as the ground truths.For instance, Pavia's ground truth is Ω Pavia OMPBS (N) = {91, 62, 16, 1, 34, 3, 73, 105, 5, 46, 85, 8, 83, 2, 11, 78}.To evaluate BS correctness, the accuracy index (ACC) is defined by Base on Equation (15), ACC(n) indicates the ratio of the target bands in overall p-selected bands, when the number of received pixels is n.Higher ACC values mean higher BS correctness.Figure 6a-d plots the ACC curves of OMPBS, PSP-OMPBS, S-PSP-OMPBS, and B-PSP-OMPBS, all implemented on the Pavia data.Several observations can be seen: 1.
Comparing Figure 6a with Figure 6b, the tendency of the ACC curves of OMPBS and PSP-OMPBS are nearly the same.This implies that the derivation of Re-OMPBS is correct.In fact, these two curves are still slightly different in some regions.Based on our study, the difference was caused by the numerical errors produced by using the recursive equations derived from Woodbury's Identity.2.
In the OMPBS and PSP-OMPBS curves, the ACC values kept at 40~60% when n is less than 1.5 × 10 5 .This is mainly due to the fact that these values were generated by using the original BIS/BIP PTS.The BS results could not stabilize quickly.

3.
The ACC results of using step sequence are shown in Figure 6c.We find that the ACC values located between 60~80% at n ∈ [2000,125,000].After receiving 1.25 × 10 5 pixels, the ACCs increased to 90%.In other words, using the PTS formed by uniform sampling drastically accelerated the speed of convergence.It suggests that using full pixels is not necessary to obtain the correct BS result using PSP-OMPBS.The nearly-correct BS result could be obtained during transmission.4.
Using a block pixel transmission sequence also accelerated the BS convergence.Figure 6d shows the corresponding ACC results of B-PSP-OMPBS.Similar conclusions can be drawn.It only requires 40,000 pixels to reach 65-70% accuracy.After receiving 60,000 pixels, the ACC grew to over 80%.The overall ACC performance is undoubtedly better than OMPBS, PSP-OMPBS, and even S-PSP-OMPBS.5.
In Figure 6a-d , it can be seen that the ACC curves remain flat in some segments.In those regions, the BS results kept the same.That is, the new incoming pixels did not change the BS results.6.
All the curves reached 100% at n = N = 207,400.
In the remote sensing community, the Purdue image was considered as a tough image for algorithm evaluation because of its noise and heavily-mixed properties.Figure 7a-d plots the ACC curves of OMPBS, PSP-OMPBS, S-PSP-OMPBS, and B-PSP-OMPBS, implemented on the Purdue data, respectively.Theoretically, the larger p and the special properties of the Purdue image, may lead to different consequences.Compared with the Pavia case, we have several interesting findings.Firstly, the oscillation of the four ACC curves is more notable.This is probably due to the noise property of the Purdue image.Secondly, the tendency of OMPBS and PSP-OMPBS curves are the same.This is in accordance with our expectations, because both methods are essentially the same.However, their ACC values are slightly different at some places.The PSP-OMPBS curve of the Purdue image is more unstable and lower at middle region.We think such a phenomenon is caused by the larger setting of p.The greater p is set, the more numerical errors will be accumulated while selecting later bands at each n.Thirdly, and most importantly, we found that the advantages of using sampled sequences (i.e., step and block) is not obvious in the Purdue experiment.We think the strange phenomenon is caused by the homogeneity of spectral profiles of the 16 ground classes.In this case, using any kind of PTS may obtain analogous pixel sets, so that the BS results are similar.From this point of view, in the case of the Purdue data, the OMPBS selects bands based more on spectral variation, instead of spatial/geographical variation.Despite this issue, it is still observed that S-PSP-OMPBS and B-PSP-OMPBS outperformed OMPBS a little bit with regard to ACC stability.The ACC trends of S-PSP-OMPBS and B-PSP-OMPBS did not drop after they reached over 80%, particularly at n ∈ [6000,13,000].In conclusion, the ACC performance of PSP-OMPBS methods varies with the properties of images.Using PSP-OMPBS on the images with lower noise seemed to produce more stable ACC curves.The sampled pixel sequences were more suitable for the images with larger heterogeneity in the spatial domain.The spectral similarity of ground classes was another important factor.Finally, we can find the curve of B-PSP-OMPBS did not end at 100% when n = N = 21,025.This was caused by the numerical error.Tables 3 and 4 list the corresponding BS results of the Pavia and Purdue experiments, where n is selected by 50% of N for each case.The last row shows the BS ground truth.It can be seen that the BS accuracies of PSP-OMPBSs are better than OMPBS at the 50% middle of transmission.

Land Cover Classification
Image classification is a common procedure to evaluate the effectiveness of a BS approach.For each dataset, we implemented supervised classification using a PSP-OMPBS-selected band Ω P (n), with a support vector machine (SVM) classifier [40].The radial basis function (RBF) kernel was adopted for SVM with the selected parameter [σ/4, σ/2, σ, 2σ, 4σ] where r is calculated by the average pairwise distance among training data σ = E x i − x j .The training samples were obtained by randomly selecting 10% of date samples in each class, according to the ground truth maps shown in Figures 4b and 5b.The other 90% of samples were used as the test samples.
For measuring the classification performance, average accuracy (AA), overall accuracy (OA), and Cohen's kappa coefficient were used as the performance metrics.Since PSP-OMPBS provided 207,400 and 21,025 BS results for two datasets respectively, we simply chose the BS results using 20%, 40%, 60%, 80%, and 100% of N for the experiment.Three PSP-OMPBS algorithms were considered.The classification result of using full bands was also used for the comparison.
Table 5 lists all the AA, OA, and kappa coefficient values of the SVM results performed on PSP-OMPBS selected bands for both the Pavia and Purdue datasets.We observed that using the selected bands of PSP-OMPBS preserved sufficient spectral information to achieve comparative classification performance, compared to the results of using full bands.The highest accuracies of the AA, OA, and kappa coefficient values reached were 0.907, 0.868, and 0.874, respectively, which is close to the full bands results of 0.914, 0.880, and 0.884, respectively.This simply implies that the bands selected by PSP-OMPBSs are informatively complementary.
On the other hand, in the Purdue case, the most accurate results for the AA, OA, and kappa coefficient values were 0.736, 0.600, and 0.693, respectively, and were slightly worse than the full bands results of 0.756, 0.586, and 0.716.This might be simply because only using 34 bands is not enough for the Purdue image classification.This issue can be resolved by increasing p, fine-tuning kernel parameters, or use other types of SVM classifiers that utilize spectral-spatial information to improve overall performance.Such a study is beyond the scope of this paper, and thus is not included.
To evaluate classification performance in visual assessment, Figures 8 and 9 plot the classification maps of SVM using the bands selected by PSP-OMPBS, S-PSP-OMPBS, and B-PSP-OMPBS at n = 60% of N for the Pavia and Purdue data, respectively.The maps using full bands are included for comparison.In the Pavia case, it can be seen that using PSP-OMPBS-selected bands could generate nearly the same classification maps compared to those obtained by using full bands.In the Purdue case, the quality of the PSP-OMPBS-generated maps are a little worse than the full bands map, due to the insufficient number of selected bands.To analyze the classification performance in the progressive manner, Figure 10 plots the kappa OA/AA/kappa coefficient curves of the SVM classification using uniformly selected BS results, performed by PSP-OMPBSs on the n-axis with interval 25 for the Pavia dataset.In those figures, the x-axis (n) indicates band set Ω p (n), and the y-axis denotes the corresponding accuracy metrics of the classification.Several observations can be found.First, all the metric curves of PSP-OMPBS and S-PSP-OMPBS are not stable in the beginning n ∈ [1,200] because of the poor BS quality that occurs with a number of received pixels.With too few pixels, PSP-OMPBS could not select correct bands to fulfill the complementary spectral information.This resulted in unstable and lower classification performance.In contrast, the curves of B-PSP-OMPBS could be consistent because the transmitted pixels were uniformly sampled from the image.Second, the overall averaged performances of S-PSP-OMPBS and B-PSP-OMPBS are slightly better than OMP-BS for these three criteria.Third, when n is greater than 250, the curves of all three methods tended to be consistent.According to our extended study, all the curves will stay roughly at the same level, with extremely low deviation in the future time n ∈ [600,N].Similarly, Figure 11 shows the OA, AA, and kappa coefficient curves for the Purdue experiment.Similarly, there are some interesting findings.First, the classification performance is low at the beginning n ∈ [1,100].Second, the Purdue image seemed to require more pixels to stabilize the BS quality.We can find that the curves of PSP-OMPBS, S-PSP-OMPBS, and B-PSP-OMPBS could not converge at n ∈ [100,600] in each figure, particular Figure 10b.This is probably because of the heavily mixed and noisy properties of the Purdue image.Each newly incoming pixel may easily disturb the result of band selection.According to our investigation, this phenomenon will last up to n = 4000.Finally, the S-PSP-OMPBS and B-PSP-OMPBS did not outperform PSP-OMPBS in the Purdue case.We think this is caused by the homogeneity of spectral profiles of the 16 ground classes.In this case, using a sampled transmission sequence could not significantly improve the collecting of spectral information.Figure 12a shows the required computing time of progressive OMPBS, PSP-OMPBS, S-PSP-OMPBS, and B-PSP-OMPBS, implemented on the Pavia dataset, where the x-axis presents the n value and the y-axis presents the corresponding computing time.It can be observed that the computing time of OMPBS significantly increases when n increases.This is due to the increase in size of B, which results in greater computational complexity for OMPBS.The variation of computing time increases with n, too.On the contrary, by virtue of recursive processing, the PSP-BSs produced almost flat curves.The computing time is stably under 0.1 s in most n indices.This implies that PSP-BS has the potential to be run in a real-time manner during transmission.Figure 12b-d further shows the individual time curve for PSP-OMPBS, S-PSP-OMPBS, and B-PSP-OMPBS, respectively.It could be seen that the computing times of most of the n regions stick on the straight trend line.In those cases, the BS results remained the same, so that the recursive equations could be applied continuously.On the other hand, we can find some of the "peak phenomenon" shown in the curves.Those peaks occurred when the updating condition was not reached.Thus, the optimization problems were solved by non-recursive equations instead.The extra computing time is required.Similarly, Figure 13a shows the computing time of progressive OMPBS, PSP-OMPBS, S-PSP-OMPBS, and B-PSP-OMPBS implemented on the Purdue dataset, and Figure 13b-d further shows the individual curves of the three PSP-OMPBSs.Similarly, the three PSP-OMPBS methods required significantly less computing time than OMPBS.The required time for PSP-OMPBSs at is less than 0.07 s from n = 1-21,025, while OMPBS needs over 1 s when n is larger than 14,000.Again, this implies the superior computational efficiency provided by PSP-OMPBS.Comparing Figure 13 with Figure 12, we can interestingly find that the peak phenomenon occurred more frequently in the Purdue case.This is probably due to the noisy property of the Purdue image.Under the circumstances, the BS result varies easily if the new incoming pixel is a noisy pixel.This results in the instability of BS on the n-axis, and thus reduces the opportunity of using recursive equations.Fortunately, the peak phenomenon disappeared gradually after sufficient pixels were received, since in later periods of transmission, the received pixels represented the image well.In other words, the BS results tended to be more stable.It is worth mentioning that PSP-OMPBSs required a little more computing time compared to OMPBS at the beginning of transmission.Figure 14a,b shows the zoom-in computing time plots of the Pavia and Purdue datasets at n = 100-1000 and n = 100-600.There is a noticeable intersection between the OMPBS and PSP-BS curves at roughly n = 300.When n is less than 300, PSP-OMPBS requires more computing time because it needs to perform an additional logic operation (i.e., Step 4 in PSP-OMPBS algorithm), which leads to an increase of overall computing time.As n grows, the proportion of computing logic operations decreases, and this additional burden is relatively diminished.Finally, Table 6 lists the overall accumulative computing time of the progressive experiments of Figures 12 and 13.It was found that using PSP-BS could significantly reduce the amount of calculation.In addition, there was no significance difference between three PSP-OMPBS methods.

Graphical User Intervace Design
In order to analyze the relationship between the ground location of the received pixels, the produced values of quantitative indexes, and the BS results of progressive processing, a Matlab graphical user interface (GUI) was developed, as shown in Figure 15.It allows users to load different image data, input p values, and choose different PTS methods (original, step, block) with the corresponding parameters.Once all the inputs are loaded and the stage button is pressed, the PSP-OMPBS starts to do real-time BS simulation.The red square in the left will show the image scene, where the red dots present the locations the received pixels.In the top-right corner, the green square will record the statistics, including the number of received pixels (n), the processing time, and ACC value for the current n.The time curve of a short period is drawn.In the bottom-right corner, the yellow square shows the BS results of the previous stage and the current stage, in which the red bins denote the instant results of PSP-OMPBS and the black bins denote the BS ground truth.
For an easy illustration, the experiment shown in Figure 15 was performed by B-PSP-OMPBS on the Pavia dataset.The parameters were set by p = 10 and b = 10.Over time, we can observe that the image gradually filled with red dots (transmitted sample pixels), and the ACC increased until it reached 100%.As a result, the BS can be fully monitored in the whole progressive process.

Discussion
With the perspective of HSI application, there are several advantages of using PSP-BS.First, we can observe all the BS results in the whole transmission, where the local BS results (i.e., the BS results in a particular time segment) would not be missed.This is one of the most attractive features in progressive processing [34,35].For instance, some "weaker bands" would appear during the transmission and disappear later on.Those bands might be very important for analyzing some specific types of ground materials, and should be considered to be selected.Besides, by virtue of progressive processing, the PSP-BS can easily be implemented on specific spatial regions (i.e., sub-image) within the image scene.Analyzing the partial data is also important.The second advantage of PSP-BS is the saving of storage volume.After PSP-BS is done, the redundant bands can be immediately removed, without waiting to apply the BS algorithm again, which saves time and storage space.The third advantage of PSP-BS is the saving of transmission bandwidth.Once the BS results converge in the early stage of transmission, we can ask the transmitter to re-transmit the rest of the image data only the selected spectral dimensions, in order to save bandwidth.This is particularly important if the transmission bandwidth is limited, or the volume of original data cube is extremely large.
Similar to PBP, that PSP-BS obeys the principle that the time of processing one band must be less than the time of transmitting one band; in PSP-BS, the processing time at each stage must be less than the time to receive a new spectral pixel.In fact, both PBP and PSP methods suffer from the computation complexity issue when the image size or spectral dimension increases.So a balance between computation capability and transmission bandwidth is a prerequisite for the real-time process.This issue has not been emphasized in the related literature.Unlike the PBP works [33][34][35], the PSP-OMPBS algorithm is further involved with the iterations related to sequential search.The computing time may increase with the drastic increase of p. Thus, the computation time must be further reduced to satisfy the prerequisite.Fortunately, the calculation of PSP-OMPBS can be accelerated by parallel processing.For instance, the task of Step 2, and any terms about matrix multiplication in other steps, can be allocated to different cores of the central processing unit (CPU) to compute.The calculation of the matrix inverse can also be reduced by using a graphics processing unit (GPU).If the processing time cannot be reduced to fulfill the real-time process, we still can perform BS with larger n interval.In this case, deriving new recursive equations is necessary.We leave that to our future work.

Conclusions
This paper presents an instant BS method based on progressive sample processing (PSP), called PSP-OMPBS.PSP-OMPBS uses a recursive algorithm to efficiently accelerate computation.It processes BS pixel-by-pixel according to the BIS/BIP format, or by re-arranged pixel transmission sequences.Unlike traditional BS methods, which must re-implement the total received data, PSP-OMPBS can immediately obtain BS results when receiving a new pixel by referring to past information.The experiments conducted on two hyperspectral datasets show that PSP-OMPBS can instantly output the BS of each stage with very low computing time.Besides, by adopting different types of pixel transmission sequences, it is proved that using sampled sequences can significantly accelerate the BS convergence speed.Those advantages have allowed PSP-OMPBS to be applied effectively in a real-time manner during pixel transmission, and have potential for real-time BS monitoring.Our future work includes reducing the computational complexity by parallel processing, minimizing the numerical errors, considering other pixel transmission strategies, fusing different BS algorithms with PSP, and extending PSP to other topics in the field of hyperspectral image processing.

Figure 1 .
Figure 1.An illustration of the difference of traditional band selection (BS) and progressive sample processing of band selection (PSP-BS).(a) The traditional BS algorithm is implemented on the collected hyperspectral image (HSI) data.(b) The PSP-BS is implemented in data sample transmission.In PSP-BS, the real-time BS monitoring can be achieved.

Figure 2 .
Figure 2. Conceptual flowchart of progressive sample processing-based OMPBS (PSP-OMPBS).The center part is the proposed PSP-OMPBS method which includes two blocks: formation of the pixel transmission sequence (PTS), and implementation of the recursive OMPBS (Re-OMPBS) algorithm.The transmission is carried out from the PTS block to the Re-OMPBS block.

Figure 3 .
Figure 3. Illustration of the pixel transmission sequence of a 9 × 9 hyperspectral image formed by different methods.(a) A band-interleaved-by-sample/pixel (BIS/BIP) sequence.(b) Step method with k = 5.(c) The block method with b = 3.

4. 1 .
Hyperspectral Dataset and Experimental SettingThe first dataset used in the experiments was a real hyperspectral image, which was collected by the reflective optics system imaging spectrometer (ROSIS) optical sensor over an urban area of the University of Pavia.The Pavia image measures 610 × 340, with very high spatial resolution of about 1.3 m per ground pixel.The original data contains 115 spectral bands, with a spectral range from 0.43 to 0.86 µm.After removing 12 noisy bands, the remaining 103 bands were used for the experiments.
Figure 4a-c respectively show the image scene of the 50th band, the geometric locations of all target classes, and the corresponding spectral signatures.According to the ground truth in Figure 5b, there are 9 classes in this image scene, consisting of several urban targets, such as vegetation, soil, and roads.The second dataset used in our experiments is another real hyperspectral image data, Purdue's Indiana Indian Pine test site, which was collected by the airborne visible-infrared imaging spectrometer (AVIRIS) system.It has been extensively studied in the literature and provides a good candidate for those who are interested in algorithm design and analysis.It has a 20 m spatial resolution and a 10 nm spectral resolution in the range of 0.4-2.5 µm with size 145 × 145 pixel vectors, taken from an area of mixed agriculture and forestry in Northwestern Indiana, U.S. It was recorded in June 1992 with 220 bands, among which bands 104-108 and 150-162 were removed, whereas the remaining 202 bands were retained.Figure 5a,b shows the image of band 20 and the ground truth map, respectively.The ground truth map contains 16 crops classes and one background class.

Figure 4 .
Figure 4. ROSIS image scene: University of Pavia.(a) Band 80.(b) Ground truth map for nine classes.(c) Spectral signatures of nine classes.

4. 4 .
Computing Time To validate the real-time capability in BIS/BIP transmission, we first compare three PSP-OMPBSs with the OMPBS performed in a progressive manner.The number of transmitted pixels, n, was set from 1 to N for image.The experiments were conducted on a computer with Intel i7-4790 3.6 GHz CPU, 16 GB RAM, Windows 7 and Matlab 2015.The values of computing time are reported by the average of ten random runs.

Figure 14 .
Figure 14.Zoom-in plots of computing time at the beginning of progressive process: (a) the Pavia data, and (b) the Purdue data.

Table 2 .
The parameters used in the experiments.

Table 3 .
List of the 16 bands, selected by four different methods, implemented on the Pavia data at n = 103,700.

Table 4 .
List of the 34 bands, selected by four different methods, implemented on the Purdue data at n = 10,513.

Table 5 .
Classification performance, in the form of the overall accuracy (OA)/average accuracy (AA)/kappa coefficient, of support vector machine (SVM) performed on the PSP-OMPBS-selected bands Ω p (n) in different n stages.

Table 6 .
The overall required computing time (in seconds) in the experiments.