Fast real-time onboard processing of hyperspectral imagery for ...
J Real-Time Image Proc (2009) 4:273–286 DOI 10.1007/s11554-008-0106-9
SPECIAL ISSUE
Fast real-time onboard processing of hyperspectral imagery for detection and classification Qian Du Æ Reza Nekovei
Received: 9 June 2008 / Accepted: 5 November 2008 / Published online: 10 December 2008 Springer-Verlag 2008
Abstract Remotely sensed hyperspectral imagery has many important applications since its high-spectral resolution enables more accurate object detection and classification. To support immediate decision-making in critical circumstances, real-time onboard implementation is greatly desired. This paper investigates real-time implementation of several popular detection and classification algorithms for image data with different formats. An effective approach to speeding up real-time implementation is proposed by using a small portion of pixels in the evaluation of data statistics. An empirical rule of an appropriate percentage of pixels to be used is investigated, which results in reduced computational complexity and simplified hardware implementation. An overall system architecture is also provided. Keywords Hyperspectral imagery Detection Classification Real-time processing Fast processing
1 Introduction A hyperspectral sensor collects hundreds of images with different wavelengths for the same area on the earth. Hyperspectral imaging is of great interest due to the fact that the high-spectral resolution (*10 nm) enables more
Q. Du (&) Department of Electrical and Computer Engineering, Mississippi State University, Starkville, MS, USA e-mail: [email protected] R. Nekovei Department of Electrical Engineering and Computer Science, Texas A&M University-Kingsville, Kingsville, TX, USA
powerful diagnostic capability in detection, classification, and quantification than the traditional multispectral imaging. However, its vast data volume and high data complexity make the data processing and analysis very challenging. In the traditional way, images are downlinked to ground stations, where data are processed and analyzed. This results in a long turn around time, which is a major bottleneck in remote sensing applications. Real-time processing and analysis are desired to provide immediate products to support fast decision-making [1–5]. It is highly desirable that onboard processing can be performed [6]. Thus, in our research, real-time data analysis is started right after data are received onboard. In other words, intermediate results are available before data collection is completed. In real-time data analysis, computation is timeconsuming because intermediate results are generated continuously. In this research, we will present a useful approach to reduce computation time such that hardware implementation can be significantly simplified, as preferred for onboard processing. It is worth mentioning that our research is focused on detection and classification, and the algorithms we implement in real-time include the computation of data correlation matrix R or covariance matrix R and their inverses, R-1 or R-1. Therefore, their real-time implementation becomes how to update the R-1 or R-1 as pixels being received. Although some studies have been presented [4, 5], we will focus on the reduction of computational complexity in real-time implementation in this paper. We will also provide different implementation strategies for data with different formats: band interleaved by pixels (BIP), band interleaved by lines (BIL), and band sequential (BSQ) [7].
123
274
J Real-Time Image Proc (2009) 4:273–286
The contributions of this paper include: 1.
2. 3.
the development of real-time implementation strategies for any detection or classification algorithms using data whitening and matched filters; the development of real-time implementation strategies for any remote sensing data formats; the development of an efficient strategy for reducing computational complexity.
This paper is organized as follows. Section 2 briefly introduces the popular detection and classification algorithms for hyperspectral imagery. Section 3 provides realtime implementation strategies for data in three different formats. Section 4 proposes an approach to save the computational complexity. Section 5 presents the system diagram. Section 6 draws the conclusions.
2 Detection and classification algorithms Assume a hyperspectral image has N pixels and L spectral bands. Let a pixel vector be denoted as x, and the entire data matrix X = [x1, x2,…, xN]. Let the sample correlation P matrix R be defined as R ¼ N1 Ni¼1 xi xTi ; which can be related to the sample covariance matrix R and sample mean l by R = R - l lT. If the target to be detected is known, it is denoted as d; if the p targets to be classified are known, they are denoted as D = [d1, d2,…, dp]. Seven algorithms are briefly summarized as follows. They have the common features: (1) using the term R-1 or R-1 for data whitening; (2) using the matched filter dT or DT for detection or classification. It is noteworthy that if the desired target is unknown, then the pixel itself xT acts as the term to be matched, which is the idea in the wellknown Reed–Xiaoli (RX) algorithm for anomaly detection. These algorithms are our research focus because: (1) they can outperform other existing algorithms due to their excellent performance in background suppression; (2) they are suitable to remote sensing images in an unknown circumstance since they require least prior information; (3) they can be easily implemented in real-time with simple and similar hardware architecture. •
RX algorithm [8, 9]:
The well-known RX algorithm [8] is a detector for anomalies, defined as pixels whose spectral signatures are different from their surroundings. The global detector can be written as wRX = xTR-1x [9]. •
Constrained energy minimization (CEM) [10]:
The CEM detector can be written as wCEM = (dTR-1d)-1R-1d. To detect if d is contained in a pixel x, we can simply apply wTCEMx, i.e., ðdT R1 dÞ1 dT R1 x:
123
•
Kelly’s generalized likelihood ratio test (KGLRT) [11]:
This generalized likelihood ratio test is given by T 1 2 d R x T 1 ; d R d 1 þ xT R1 x=N 0 where N’ is the number of samples used in the estimation of R. •
Adaptive matched filter (AMF) [12]:
When the number of samples N is a very large value, the KGLRT is reduced to a simple format: ðdT R1 dÞ1 ðdT R1 xÞ2 : We can see that it is close to the CEM except that the numerator has a power of two. •
Adaptive coherence estimator (ACE) [13]:
The estimator can be written as T 1 2 d R x T 1 : d R d xT R1 x It is similar to AMF except that a term similar to the RX algorithm is included in the denominator. •
Target-constrained (TCIMF) [14]:
interference-minimized
filter
The filter can be expressed as ðDT R1 DÞ1 ðDT R1 xÞ; where p dimensions correspond to p detection maps. It can not only detect targets but also separate them from each other. •
Constrained linear discriminant analysis (CLDA) [15]:
The classifier can be simply put as: ðDT R1 DÞ1 ðD R1 xÞ: It is similar to TCIMF except that R-1 is used rather than R-1. Its major advantage over the traditional classifiers, such as maximum likelihood classifier, is that it does not need to know the complete class information in the data [16]. In other words, it needs to know only the p targets to be classified, even when p is much less than the actual number of classes in the data. T
3 Real-time implementation In our research, we assume that an image is acquired from left to right and from top to bottom. Three real-time processing fashions will be discussed to fit the three remote sensing data formats: pixel-by-pixel processing for BIP format, line-by-line processing for BIL format, and band-by-band processing for BSQ formant. In the pixelby-pixel fashion, a pixel vector is processed right after it is received and the analysis result is generated within an acceptable delay; in the line-by-line fashion, a line of pixel vectors is processed after the entire line is received;
J Real-Time Image Proc (2009) 4:273–286
275
in the band-by-band fashion, a band is processed after it is received. To speed up the process and simplify hardware implementation, one strategy is to replace the data covariance matrix R with the data correlation matrix R as discussed in [4, 5], as both provide similar results. Thus, in order to implement the detection or classification algorithms in Sect. 2 in real-time, the key becomes the adaptation of R-1. In other words, R-1 at time t should be quickly calculated by updating the previous R-1 at time t - 1 using the data received at time t, without completely recalculating the R and R-1. As a result, the intermediate data analysis result (e.g., target detection or classification) is available in support of decision-making even when the entire data set is not received; and when the entire data set is received, the data analysis result is finalized.
~ 1 xt Þ; which In (3), L2 multiplications are needed for ðR t1 ~ 1 Þ; ðxT R ~ 1 xt Þ requires is the same as that for ðxTt R t t1 t1 additional L multiplications only; determining the product ~ 1 xt Þ1 ; and xT R ~ 1 ; ~ 1 xt ; ð1 þ xT R of the three terms: R
This format is easy to handle because a pixel vector of size L 9 1 is received continuously, which well fits the spectral-analysis based algorithms previously discussed in Sect. 2. ~ denotes NR, which is the correlation matrix Let R without being normalized by the number of pixels. Using ~ 1 in the algorithms in Sect. 2 does not change the result, R but simplifies the update process. Supposed that at time t a pixel vector, denoted as xt, is received. The data matrix Xt including all the pixels received up to the time t is Xt = [x1, x2,…, xt]. The sample correlation matrix at time t and t – 1 are denoted as Rt and Rt-1, respectively, and they are related as
and
The following Woodbury’s formula [17] can be used to ~ 1 : update R t 1 ðA þ BCDÞ1 ¼ A1 A1 B C1 þ DA1 B DA1 ð2Þ where A and C are two positive-definite matrices, and the sizes of matrices A, B, C and D allow the inverse operation of (A ? BCD). It should be noted that (2) is for the most general case. Actually, A, B, C, and D can be reduced to vector or scalar as long as (2) is applicable. ~ 1 in (1) can be By using the Woodbury’s formula, R t updated as: 1 ~ 1 ¼ R ~ 1 xt ~ 1 : ~ 1 R ~ 1 xt 1 þ xT R R xTt R ð3Þ t t1 t1 t1 t t1 ~ 1 xt Þ in (3) is a scalar. This means no Here ð1 þ xTt R t1 matrix inversion is involved in each adaptation. Equation (3) can be initiated with a small random matrix ~ 1 [4]. R 0
t1
If the data are in BIL format, we can simply wait for all the pixels in a line is received to form pixel vectors. Let M be the total number of pixels in each line. Let the line received at time t form a data matrix Yt = [xt1xt2,…, xtM]. Then (1) and (3) become ~t ¼ R ~ t1 þ Yt YT R t
ð1Þ
t
t1
3.2 BIL format
3.1 BIP format
~t ¼ R ~ t1 þ xt xT : R t
t
t1
requires L2 multiplications individually. Thus the total number of multiplications in the update for each pixel vector using (3) is 2L2 ? L. For the entire image, this number becomes N(2L2 ? L).
ð4Þ
1 ~ 1 ¼ R ~ 1 ~ 1 ~ 1 R ~ 1 Yt IMM þ YT R R YTt R t t1 t1 t1 t t1 Yt
ð5Þ
respectively, where IM9M is an M 9 M identity matrix. ~ 1 Yt Þ in (5) is a matrix. This means matrix ðIMM þ YTt R t1 inversion is involved in each adaptation. ~ 1 Yt Þ; In (5), ML2 multiplications are needed for ðR t1 T ~ 1 T ~ 1 which is the same as ðY R Þ; ðY R Yt Þ requires t
t1
t
t1
additional M2L multiplications only; determining the ~ 1 Yt Þ1 generally requires O(M3) inverse ðIMM þ YT R t
t1
multiplications [18]; determining the product of the three 1 ~ 1 Yt ; ðIMM þ YT R ~ 1 ~ 1 in (5), terms: R and YTt R t1 t1 t t1 Yt Þ requires (ML2 ? M2L) multiplications. Therefore, the total multiplications in the update for each line is 2(ML2 ? M2L) ? O(M3). For the entire image, this number becomes (2(ML2 ? M2L) ? O(M3))N/M = (2NL2 ? 2NML) ? O(M3)N/M. 3.3 BSQ format If the data format is BSQ, the sample correlation matrix R and its inverse R-1 have to be updated in a different way, because there is no completed pixel vector available before the entire data set is received. Let r1 denote the auto-correlation of pixels in Band 1, which is a scalar, calculated by the average of pixel squared values in Band 1. Let R2 be the correlation matrix between the first two bands. Then r1 can be related to R2 as: r1 r12 R2 ¼ ; r21 r2 where r2 is the auto-correlation of pixels in Band 2, r12 = r21 is the cross-correlation of corresponding pixel
123
276
J Real-Time Image Proc (2009) 4:273–286
values in Band 1 and 2. Therefore, Rt can be related to Rt-1 as Rt1 rt1;t Rt ¼ T ð6Þ rt1;t rt;t where rt is the auto-correlation of pixels in Band t, rt1;t ¼ ½r1;t ; ; rj;t ; ; rt1;t T is a (t - 1) 9 1 vector with rj,t being the cross-correlation of corresponding pixel values in Band j and t. Equation (6) shows that the dimension of Rt is increased as more bands are received. -1 When R-1 t-1 is available, it can be used to calculate Rt -1 by modifying Rt-1 with rt,t and rt-1,t. The following partitioned matrix inversion formula [17] can be used for this purpose. Let a matrix A be partitioned as " # A11 A12 A¼ : A21 A22 Then its inverse matrix A-1 can be calculated as
" A
1
¼
1 A11 A12 A1 22 A21 1 A22 A21 A1 A21 A1 11 A12 11
Let A11 = Rt-1, A22 = rt,t, A21 = rTt-1,t. According to (7)
A12 = rt-1,t,
1 Rt1 rt1;t rt;t rTt1;t 6 ¼4 1 rt;t rTt1;t R1 rTt1;t R1 t1 rt1;t t1
3 1 1 Rt1 rt1;t rt;t rTt1;t rt1;t rt;t 7 1 5 1 rt;t rTt1;t Rt1 rt1;t
1 T 1 ¼ R1 t1 þ aRt1 rt1;t rt1;t Rt1
ð9Þ
aR1 t1 rt1;t
#
a ð10Þ
Thus, all the elements in (10) can be generated by simple matrix multiplication. Actually, in this case, no
123
ð7Þ
and
1 In (8), a ¼ ðrt;t rTt1;t R1 is a scalar; as for t1 rt1;t Þ 1 T ðRt1 rt1;t rt;t rt1;t Þ ; the Woodbury’s formula in (2) can be used again, i.e., 1 Rt1 rt1;t rt;t rTt1;t 1 1 1 T 1 ¼ R1 rTt1;t R1 t1 þ Rt1 rt1;t rt;t rt1;t Rt1 rt1;t t1
and (8) becomes " 1 1 T Rt1 þ aR1 t1 rt1;t rt1;t Rt1 1 Rt ¼ arTt1;t R1 t1
1 T T 1 1 ðrt;t rTt1;t R1 t1 rt1;t Þ rt1;t Rt1 (since rt-1 Rt-1 does not need to be re-calculated); and additional (t - 1)2 multiplications for ðRt1 rt1;t rt;t rTt1;t Þ1 : Thus, the total number of multiplications in the update for the whole P image is N þ Lt¼2 ðN þ Nðt 1Þ þ 2ðt 1Þ2 þ2ðt 1ÞÞ NL þ ðN þ 2ÞLðL 1Þ=2 þðL 1ÞLð2L 1Þ=3 ¼ NL2 =2 þNL=2 þ 2ðL3 LÞ=3
# 1 A11 A12 A1 A12 A1 22 A21 22 1 A22 A21 A1 11 A12
2
R1 t
operation of matrix inversion is needed when reaching the final R-1. The intermediate result still can be generated by applying the R-1 to the first t bands. This means the t spectral features in these t bands are used for target detection and discrimination. This may help find targets at early processing stages. It is worth mentioning that R-1 is symmetric and three elements in the matrix in (8) need to be calculated. N and N(t - 1) multiplications are needed to calculate rt,t, and rt,t-1, respectively; (t - 1)2 ? (t - 1) multiplications for 1 a ¼ ðrt;t rTt1;t R1 t1 rt1;t Þ ; ðt 1Þ multiplications for
ð8Þ
4 Fast implementation 4.1 Computational complexity reduction for R-1 evaluation In this paper, we will propose another way to save the number of computations. In the offline version, the major computation time is spent in the calculation of the data correlation matrix R. For an image with N pixels and L bands, the number of multiplications for R is NL2. Its matrix inversion generally needs O(L3) multiplications, depending upon the methods employed [18]. For a remote sensing image, N L, then the key to reducing computational complexity is to use a part of pixels for the R calculation. Fortunately, remote sensing image usually has high-spatial correlation; using a small part of pixels can still capture the major data statistics.
J Real-Time Image Proc (2009) 4:273–286
277
Table 1 Computational Complexity in real-time and offline processing of different formats g % pixels are used
All the pixels are used Offline BIP
Real-time
NL2 ? o(L3) 2
3
Offline
N(2L2 ? L) 2
Real-time
g NL2 ? o(L3) 3
NL ? o(L )
(2NL ? 2NML) ? o(M )N/M
g NL ? o(L )
(2gNL2 ? 2gNML) ? o(M3)g N/M
BSQ
NL2 ? o(L3)
NL2/2 ? NL/2 ? 2(L3-L)/3
g NL2 ? o(L3)
g NL2/2 ? gNL/2 ? 2(L3-L)/3
4.2 Experiments 4.2.1 HYDICE experiment: target detection The hyperspectral digital imagery collection experiment (HYDICE) image scene shown in Fig. 1 was collected in Maryland in 1995 from the flight altitude of 10,000 ft with approximately 1.5 m spatial resolution in 0.4–2.5 lm spectral region. The atmospheric water bands with low signal-to-noise ratio were removed reducing the data dimensionality from 210 to 169. The image scene has 128 lines and the number of pixel in each line M is 64. This scene includes 30 panels arranged in a 10 9 3 matrix. Each panel is denoted by pij with row indexed by i = 1,…, 10 and column indexed by j = a, b, c. The three panels in the same row were made from the same material of size 3 m 9 3 m, 2 m 9 2 m, 1 m 9 1 m, respectively, which were considered to belong to the same class, Pi.
3
g N(2L2 ? L)
BIL
Table 1 summarizes the number of multiplications for each data format in real-time and offline processing, where g is the percentage of pixels participating in the R-1 evaluation. Here the offline processing means all the data are processed simultaneously; in other words, the processing is started after an entire data set is received. We can see that the computational complexity in offline processing of BIP and BIL format is much less than that in their realtime processing when all the pixels are used. However, when image size is small, it is possible for the real-time processing of BIP format to have less complexity if g is less than a certain value. When using g% of pixels, the computation reduction factor for all the formats is about g. In Sect. 4.2 we will show that when g is appropriately selected, the detection or classification performance will remain unchanged. When the image scene contains large homogeneous areas, then g can be small; otherwise, g should be greater such that the data statistics can be manifested by a portion of pixels. An empirical rule is to select five to ten times of the number of spectral dimensionality (i.e., the number of spectral bands). The low bound is defined such that the number of linearly independent pixels is equal to the number of spectral bands for the existence of R-1 (i.e., an R with full rank).
2
P1a, P1b, P1c P2a, P2b, P2c P3a, P3b, P3c P4a, P4b, P4c P5a, P5b, P5c P6a, P6b, P6c P7a, P7b, P7c P8a, P8b, P8c P9a, P9b, P9c P10a, P10b, P10c Fig. 1 HYDICE image scene with 30 panels
Figures 2 and 3 show the detection results using TCIMF when all the pixels and 10% pixels were involved in R-1 calculation, respectively. The difference is not perceivable. Table 2 lists the number of panel pixels correctly detected ND and the number of false alarm pixels NF when all the pixels, 10% pixels, and 5% pixels were used. The total number of panel pixels NP is 38. When all the pixels and 10% pixels were used, 28 pixels out of 38 could be detected without false alarm pixels. When 5% pixels were used, 26 pixels were detected with 4 false alarm pixels. It seems that only 10% pixels are needed for this image scene. Since this image scene has only 8,192 pixels. When 10% pixels are used, the spatial size is about 5 times of the spectral size. Real-time processing results are tabulated in Tables 3, 4 and 5. In Table 3 for pixel-by-pixel fashion, the result was the same as the one in offline processing when all the pixels were used for correlation matrix update; if only 10% pixels participated, then the result was slightly degraded with ND = 29 and NF = 13; if only 5% pixels were used, then the result was greatly degraded with NF = 303. In Table 4 for line-by-line processing, the result was the same as the one in offline processing when all the pixels were involved in the correlation matrix update; the performance was slightly better than the case of pixel-by-pixel since NF was reduced to 29 with 5% pixels. In Table 5 for band-by-band processing, the result was the same as the offline processing since the results can be compared only after all the bands are received. Based on this computer simulation, it seems that using 10% of pixels is a good choice for real-time processing.
123
278
J Real-Time Image Proc (2009) 4:273–286
Fig. 2 Detection result using TCIMF with all the pixels for R-1 calculation
Fig. 3 Detection result using TCIMF with 10% pixels for R-1 calculation
4.2.2 AVIRIS Cuprite experiment: classification The airborne visible/infrared imaging spectrometer (AVIRIS) Cuprite subimage scene of size 350 9 350
123
shown in Fig. 4 was collected in Nevada in 1997. The spatial resolution is 20 m. Originally it has 224 bands with 0.4–2.5 lm spectral range. After water absorption and lowSNR bands were removed, 189 bands were left. This image
J Real-Time Image Proc (2009) 4:273–286
279
Table 2 Performance comparison for HYDICE experiment 10%
Table 5 Real-time processing performance in HYDICE experiment (band-by-band)
Pixel percentage
100%
Panel
NP
ND
P1
3
2
0
2
0
2
0
P2
3
2
0
2
0
2
2
P3
4
3
0
4
0
3
0
NF
ND
5% NF
ND
NF
Pixel percentage
100%
10%
Panel
NP
ND
P1
3
2
0
2
0
2
0
P2
3
2
0
2
0
2
2
4 3
3 2
0 0
4 2
0 0
3 2
0 0
NF
ND
5% NF
ND
NF
P4
3
2
0
2
0
2
0
P5
6
5
0
6
0
4
0
P3 P4
P6
3
2
0
2
0
1
0
P5
6
5
0
6
0
4
0
P7 P8
4 4
3 3
0 0
3 3
0 0
3 3
2 0
P6
3
2
0
2
0
1
0
P7
4
3
0
3
0
3
2
0
P8
4
3
0
3
0
3
0
0
P9
4
3
0
3
0
3
0
4
P10
4
3
0
3
0
3
0
Total
38
28
0
28
0
26
4
P9 P10 Total
4 4 38
3 3 28
0 0 0
3 3 28
0
3
0
3
0
26
Table 3 Real-time processing performance in HYDICE experiment (pixel-by-pixel) Pixel percentage
100%
Panel
NP
ND
P1
3
2
P2
3
2
P3 P4
4 3
P5
10% NF
5%
ND
NF
ND
NF
0
2
0
2
0
0
2
0
3
41
3 2
0 0
3 2
0 0
2 3
201 1
6
5
0
5
0
5
0
P6
3
2
0
2
0
2
0
P7
4
3
0
4
13
4
59
P8
4
3
0
3
0
3
0
P9
4
3
0
3
0
3
1
P10
4
3
0
3
0
3
0
Total
38
28
0
29
13
30
303
Table 4 Real-time processing performance in HYDICE experiment (line-by-line) Pixel percentage
100%
Panel
NP
ND
P1
3
3
P2
3
P3
4
P4 P5
10% NF
5%
ND
NF
ND
NF
0
3
0
3
0
2
0
2
0
2
0
3
0
3
1
3
27
3
2
0
2
0
2
0
6
5
0
5
0
5
0
P6
3
2
0
2
0
2
0
P7
4
3
0
3
10
3
2
P8
4
3
0
3
0
3
0
P9
4
3
0
3
0
3
0
P10
4
3
0
3
0
3
0
Total
38
29
0
29
11
29
29
Fig. 4 AVIRIS Cuprite scene with five major minerals: alunite (A), buddingtonite (B), calcite (C), kaolinite (K), and muscovite (M)
scene is well understood mineralogically. At least five minerals were present: alunite (A), buddingtonite (B), calcite (C), kaolinite (K), and muscovite (M). Their approximate spatial locations of these minerals are marked in Fig. 4. Figures 5 and 6 show the results using CLDA when all the pixels and 1% pixels were involved in R-1 calculation, respectively. The difference is not perceivable. Since no pixel-level ground truth is available, the spatial correlation coefficient is employed to evaluate the similarity between the corresponding maps. A larger correlation coefficient means better performance. As listed in Table 6, the correlation coefficient is as high as 0.9 even when only 1% pixels was used. This image scene is large and contains 122,500 pixels. When 1% pixels are used, the spatial size is about 6.5 times of the spectral size. The ratio between the number of pixels
123
280
J Real-Time Image Proc (2009) 4:273–286
Fig. 5 Classification result using CLDA with all the pixels for R-1 calculation
Fig. 6 Classification result using CLDA with 1% pixels for R-1 calculation Table 6 Performance comparison using correlation coefficient for AVIRIS Cuprite experiment
Table 8 Real-time processing performance in AVIRIS Cuprite experiment (line-by-line)
Pixel percentage
10%
5%
1%
Pixel percentage
100%
10%
5%
1%
Alunite
0.9926
0.9827
0.8953
Alunite
0.8944
0.8802
0.8719
0.8306
Buddingtonite
0.9926
0.9847
0.8850
Buddingtonite
0.9484
0.9352
0.9274
0.8224
Calcite
0.9936
0.9861
0.9121
Calcite
0.9433
0.9174
0.8920
0.6930
Kaolinite
0.9903
0.9794
0.8711
Kaolinite
0.9430
0.9164
0.9190
0.8364
Muscovite
0.9905
0.9826
0.9362
Muscovite
0.9229
0.9179
0.8943
0.8156
Average
0.9919
0.9831
0.8999
Average
0.9304
0.9134
0.9009
0.7996
Table 7 Real-time processing performance in AVIRIS Cuprite experiment (pixel-by-pixel)
Table 9 Real-time processing performance in AVIRIS Cuprite experiment (band-by-band)
Pixel percentage
100%
10%
5%
1%
Pixel percentage
100%
10%
5%
1%
Alunite
0.8935
0.8843
0.8706
0.8082
Alunite
1.0000
0.9926
0.9827
0.8953
0.7321
Buddingtonite
1.0000
0.9926
0.9847
0.8850
1.0000
0.9936
0.9861
0.9121 0.8711
Buddingtonite
0.9477
0.9254
0.8655
Calcite
0.9422
0.8848
0.6456
0.6280
Calcite
Kaolinite
0.9424
0.9255
0.8836
0.8266
Kaolinite
1.0000
0.9903
0.9794
Muscovite
0.9218
0.8906
0.8831
0.6917
Muscovite
1.0000
0.9905
0.9826
0.9362
Average
0.9295
0.9021
0.8297
0.7373
Average
1.0000
0.9919
0.9831
0.8999
and the number of bands (189) may be an indictor for the appropriate percentage that should be chosen. Real-time processing results are tabulated in Tables 7, 8 and 9. We can see that the performance of line-by-line processing was better than the pixel-by-pixel processing; using 5% pixels the classification result was not
123
significantly changed in line-by-line processing and using 10% pixels the classification result was not significantly changed in pixel-by-pixel processing. Once again, the classification result for the BSQ format was the same as the offline processing. This means the pixel percentage for data statistics evaluation may need to be increased when an algorithm is implemented in real-time.
J Real-Time Image Proc (2009) 4:273–286
281
4.2.3 AVIRIS Lunar Lake experiment: anomaly detection
Table 10 Performance comparison using correlation coefficient for AVIRIS Lunar Lake experiment
The 224-band AVIRIS Lunar Lake image scene shown in Fig. 7 was used in the experiment. It is a subscene of 200 9 200 pixels cropped from the left upper corner of the Lunar Crater Volcanic Field in Northern Nye County, Nevada. Five classes are present: ‘‘red oxidized basaltic cinders’’, ‘‘rhyolite’’, ‘‘playa (dry lakebed)’’, ‘‘shade’’, and ‘‘vegetation’’. In addition, there is an anomaly class present. The number of bands was reduced to 158 after lowSNR and water absorption bands were removed. The RX algorithm was applied for anomaly detection. As shown in Fig. 8, there is no perceivable difference between the results using all the pixels and 2% pixels for R-1 calculation, where both the anomaly and some part of vegetation were found. Table 10 lists the correlation coefficients between the detection maps using all the pixels and a small portion of pixels. We can see that the correlation coefficient is as high as 0.8869 even when only 2% pixels were used. At this time, the number of pixels was
Pixel percentage
100%
10%
5%
2%
Offline
–
0.9805
0.9545
0.8869
Online (pixel-by-pixel)
1.0000
0.9805
0.9545
0.8869
Online (line-by-line)
1.0000
0.9805
0.9545
0.8869
Online (band-by-band)
1.0000
0.9805
0.9545
0.8869
reduced to 800, which is about 5 times of the number of bands. Real-time processing results are tabulated in Table 10. Since the offline result was the one using global statistics (global anomaly detector), the results were comparable only after all the data were received. We can see that in this case real-time processing for all the three formats provided the same result as in offline processing.
5 System architecture 5.1 System overview The system for fast real-time detection or classification can be summarized as follows. The block diagram is shown in Fig. 9. 1. 2.
Input the data format and choose the value of g. If the format is BIP, • •
• Fig. 7 The AVIRIS Lunar Lake scene of size 200 9 200 used in the experiment
3.
If the format is BIL, • • • •
4.
Initiate R-1 0 using a small random matrix of size L 9 L; Wait until the entire line of pixel vectors is received at time t; Use the pixels selected by the PRS module to update R-1 with (5); t Go to Step 5.
If the format is BSQ, • •
Fig. 8 Anomaly detection result using RX algorithm
Initiate R-1 0 using a small random matrix of size L 9 L; When xt is received at time t, check if it is selected by the pixel random selection (PRS) module; if yes, use (3) to update R-1 t ; if not, keep the previous -1 R-1 t-1 as Rt ; Go to Step 5.
•
Use the PRS module to select the pixels in each band that will participate in R-1 update; When the t-th band is received at time t, compute r1 and rj,t for j = 1,2,…, (t - 1); Update R-1 with (10). t
123
282
J Real-Time Image Proc (2009) 4:273–286
BIP
R t−1 update
usin g p ixels xt
Format Selection
Random Pixel Selection
BIL
R t−1 update
using data lines
Detector/ Classifier Construction
Detection/ Classification
R t−1 update
BS Q
using band correlation vectors
Data Storage Fig. 9 Real-time onboard detection or classification pipeline stages
5. 6.
Construct the detector or classifier in Sect. 2. Provide the detection or classification result for the pixels/lines/bands that have been received up to time t.
5.2 Hardware architecture The major system modules have been discussed in [4], including parallel matrix vector multiplier (PMVM), matrix–matrix multiplier (MMM), and matrix inversion (MI) models. Here, we only discuss the PRS module. We also provide hardware architecture for R-1 for three dift ferent formats, which are very similar. 5.2.1 PRS module Assume 10% pixels need to be randomly selected. The PRS module will provide a bit pattern with approximately 10% of pattern data being 1 s and 90% 0 s, where a value of 1 indicates the corresponding pixel is selected for calculations and a 0 means to ignore the pixel. As shown in Fig. 10, the module includes a random bit generator (RBG) [19] which employs Zener Diode as the source of its randomness to provide approximately 50% 1 s and 50% 0 s, followed by a counter-comparator circuit to select one out of every five 1 s. A simulation of the circuit generated a bit pattern with 10.6% 1 s and 89.4% 0 s is illustrated in
Fig. 10 A module for pixel random selection (PRS)
123
Fig. 10. By slightly modifying the comparator, different percentages can be easily simulated. 5.2.2 Hardware architecture for R-1 update for BIP t format Figure 11 illustrates a parallel pipeline hardware realization of (3) where data are arriving one pixel vector at the ~ 1 is updated each time a new time. In this module R t selected pixel vector is processed. Parallel Multiply Accumulators calculate result of the matrix vector multi~ 1 xt Þ; which is pipelined to calculate the plication, ðR t1 required scalar inner product and outer vector product ~ 1 ¼ R ~ 1 1 T ~ 1 ~ 1 xt ð1 þ xT R matrix: DR t t1 xt Þ xt Rt1 : t t1 5.2.3 Hardware architecture for R-1 update for BIL t format Figure 12 illustrates how BIP architecture can be modified for a real-time line by line approximation. The matrix vector manipulation has also been reordered to take advantage of data arriving one line of each band at a time. This will improve system throughput and the pipeline performance. Additionally, the architecture has been extended with vector multiply accumulators to be able to process the entire line.
J Real-Time Image Proc (2009) 4:273–286
283
~ ~ R t−−11 x t = x Tt R t−−11
Outer Product Select
MAC1
MV1 Reg.
*
MAC2
MV2 Reg.
*
Inner_en MUX
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
D E M U X
. . .
Inner_en
*
MVL Reg.
~ ∆R t−1
MUX MACL
Inner_en
~ R t−−11
xtt
DEMUX
Inner Product Register
MV1,L Reg.
Inverse Operation
* ~ (1 + xTt R t−−11x t ) −1
Fig. 11 Architecture for BIP iterative correlation matrix inversion
5.2.4 Hardware architecture for R-1 update for BSQ t format Figure 13 illustrates the hardware architecture for BSQ format, which is similar to the previous two architectures. The three major terms in (10) are generated, from which the inverse matrix R-1 t can be easily updated. In addition to the R-1 t-1, the input includes the correlation vector at time t which can be easily generated using MAC. It is noteworthy that the dimensionality of R-1 is increased as more t bands being received; as a result, some MACs and
multipliers in early stages are actually not used with inputs and outputs being zero.
6 Conclusion In this paper, we investigate real-time onboard processing of hyperspectral imagery for detection or classification. The detectors and classifiers we employ include the term of the inverse of data correlation or covariance matrix (R-1or R-1) for data whitening, followed by a matched filter for
123
284
J Real-Time Image Proc (2009) 4:273–286
Column Select
V_MAC1 c1 c2 : cM
+
*
~ YtT R t−−11
Outer Product Select
MAC
MV1 Reg.
Inner_en MAC
MUX
MV2 Reg.
V_MAC2
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
D E M U X
. . .
MAC
Inner_en MUX
MVL Reg.
V_MACL
Yt
DEMUX
Inner_en
~ R t−−11
~ ∆R t−1
Inner Product Register
*
MV1,L Reg.
Inverse Operation
(I + Y
T t
~ R t−−11 Yt
)
−1
Fig. 12 Architecture for BIL iterative correlation matrix inversion
the target(s) to be detected or classified. In real-time processing, data processing and analysis start immediately after pixels are received. Thus, the key is to update R-1 or R-1 in real-time, which has been investigated when data are formatted as BIP, BIL, or BSQ. To speed up the process and facilitate onboard hardware implementation, we propose the use of a small portion of pixels in the update of R-1. Experimental results show that with a small portion of pixels, comparable detection or classification results can be provided. For data with large
123
homogeneous areas, only a small number of pixels need to be used. To be conservative, in the offline processing, the percentage of selected pixels can be empirically chosen such that the number of pixels is about five to ten times the number of bands (which is coincident with the observation in [20] for training set size selection in speech signal processing); this percentage for real-time implementation may need to be increased to about 5 times the one in offline processing. We also provide the system architecture for the proposed fast real-time algorithms for different data formats.
J Real-Time Image Proc (2009) 4:273–286
285
(r
T t −1,t
R t−−11
)
Outer Product Select
MAC1
MV1 Reg.
*
MAC2
MV2 Reg.
*
Inner_en MUX
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
D E M U X
. . .
Inner_en
*
MVL Reg.
αR t−−11rt −1rtT−1R t−−11
MUX MACL
R
rt −1,t
−1 t −1
DEMUX
Inner_en
Inner Product Register
MV1,L Reg.
* αrtT−1,t R t−−11
Inverse Operation
α = (rt ,t − rtT−1,t R t−−11rt −1,t )
−1
Fig. 13 Architecture for BSQ iterative correlation matrix inversion
Similarities among BIP, BLP, and BSQ hardware architectures allow the use of a single parallel pipelined computation engine for all three data formats and reduce hardware complexity for a general purpose system.
References 1. Stellman, C.M., Hazel, G.G., Bucholtz, F., Michalowicz, J.V.: Real-time hyperspectral detection and cuing. In: Optical Engineering, vol. 39, no. 7, July (2000) 2. Chang, C.-I., Ren, H., Chiang, S.S.: Real-time processing algorithms for target detection and classification in hyperspectral imagery. In: IEEE Transactions on Geoscience and Remote Sensing, vol. 39, No. 4 (2001)
3. Du, Q., Ren, H.: Real-time constrained linear discriminant analysis to target detection and classification in hyperspectral imagery. In: Pattern Recognition, vol. 36, No. 1 (2003) 4. Du, Q., Nekovei, R.: Implementation of real-time constrained linear discriminant analysis to remote sensing image classification. In: Pattern Recognition, vol. 38, No. 4 (2005) 5. Du, Q.: Unsupervised real-time constrained linear discriminant analysis to hyperspectral image classification. In: Pattern Recognition, vol. 40, No. 5 (2007) 6. Subramanian, S., Gat, N., Ratcliff, A., Eismann, M.: Real-time hyperspectral data compression using principal components transformation. In: Proceedings of the AVIRIS Earth Science and Applications Workshop (2000) 7. Jensen, J.R.: Introductory Digital Image Processing: A Remote Sensing Perspective, 3rd edn. Prentice-Hall, Englewood Cliffs (2004) 8. Reed, I.S., Yu, X.: Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution. In: IEEE
123
286
9.
10.
11. 12.
13.
14.
15.
16.
17. 18.
J Real-Time Image Proc (2009) 4:273–286 Transactions on Acoustic, Speech and Signal Processing, vol. 38, No. 10 (1990) Chang, C.-I., Heinz, D.: Subpixel spectral detection for remotely sensed images. In: IEEE Transactions on Geoscience and Remote Sensing, vol. 38, No. 3 (2000) Farrand, W.H., Harsanyi, J.C.: Mapping the distribution of mine tailing in the coeur d’Alene river valley, Idaho through the use of constrained energy minimization technique. In: Remote Sensing of Environment, vol. 59, No. 1 (1997) Kelly, E.J.: An adaptive detection algorithm. In: IEEE Transactions on Aerospace and Electronic Systems, vol. 22, No. 1 (1986) Robey, F.C., Fuhrmann, D.R., Kelly, E.J., Nitzberg, R.: A CFAR adaptive matched filter detector. In: IEEE Transactions on Aerospace and Electronic Systems, vol. 28, No. 1 (1992) Kraut, S., Sharf, L.L.: The CFAR adaptive subspace detector is a scale-invariant GLRT. In: IEEE Transactions on Signal Processing, vol. 47, No. 9 (1999) Ren, H., Chang, C.-I.: A target-constrained interference-minimized approach to subpixel detection for hyperspectral images. In: Optical Engineering, vol. 39, No. 12 (2000) Du, Q., Chang, C.-I.: Linear constrained distance-based discriminant analysis for hyperspectral image classification. In: Pattern Recognition, vol. 34, No. 2 (2001) Du, Q.: Modified Fisher’s linear discriminant analysis for hyperspectral imagery. In: IEEE Geoscience and Remote Sensing Letters, vol. 4, No. 4 (2007) Golub, G.H., Van Loan, C.F.: Matrix Computations, 3rd edn. Johns Hopkins University Press (1996) Dremmel, J.E.: Applied Numerical Linear Algebra. Society for Industrial and Applied Mathematics, Philadelphia (1997)
123
19. Stipcevic, M.: Fast nondeterministic random bit generator based on weakly correlated physical events. In: Review of Scientific Instruments, vol. 75, No. 11 (2004) 20. Makhoul, J., Roucos, S., Gish, H.: Vector quantization in speech coding. In: Proceedings of the IEEE, vol. 73, No. 11 (1985)
Author Biographies Qian Du received her Ph.D. degree in electrical engineering from University of Maryland Baltimore County in 2000. She was with the Department of Electrical Engineering and Computer Science, Texas A&M University-Kingsville, from 2000 to 2004. She joined the Department of Electrical and Computer Engineering at Mississippi State University in Fall 2004, where she is currently an Associate Professor. Her research interests include remote sensing image analysis, pattern classification, data compression, and neural networks. Dr. Du is a senior member of IEEE, a member of SPIE, ASPRS, and ASEE. Reza Nekovei received his Ph. D. degree in Electrical Engineering from The University of Rhode Island in 1994. He is currently an Associate Professor at the Department of Electrical Engineering and Computer Science, Texas A&M University-Kingsville. He was previously with the Remote Sensing Laboratory, Graduate School of Oceanography at The University of Rhode Island. His research interests include VLSI Architectures, Distributed Systems, Remote Sensing, and Medical Imaging. He is a member of IEEE and Etta Kappa Nu honor society.
SPECIAL ISSUE
Fast real-time onboard processing of hyperspectral imagery for detection and classification Qian Du Æ Reza Nekovei
Received: 9 June 2008 / Accepted: 5 November 2008 / Published online: 10 December 2008 Springer-Verlag 2008
Abstract Remotely sensed hyperspectral imagery has many important applications since its high-spectral resolution enables more accurate object detection and classification. To support immediate decision-making in critical circumstances, real-time onboard implementation is greatly desired. This paper investigates real-time implementation of several popular detection and classification algorithms for image data with different formats. An effective approach to speeding up real-time implementation is proposed by using a small portion of pixels in the evaluation of data statistics. An empirical rule of an appropriate percentage of pixels to be used is investigated, which results in reduced computational complexity and simplified hardware implementation. An overall system architecture is also provided. Keywords Hyperspectral imagery Detection Classification Real-time processing Fast processing
1 Introduction A hyperspectral sensor collects hundreds of images with different wavelengths for the same area on the earth. Hyperspectral imaging is of great interest due to the fact that the high-spectral resolution (*10 nm) enables more
Q. Du (&) Department of Electrical and Computer Engineering, Mississippi State University, Starkville, MS, USA e-mail: [email protected] R. Nekovei Department of Electrical Engineering and Computer Science, Texas A&M University-Kingsville, Kingsville, TX, USA
powerful diagnostic capability in detection, classification, and quantification than the traditional multispectral imaging. However, its vast data volume and high data complexity make the data processing and analysis very challenging. In the traditional way, images are downlinked to ground stations, where data are processed and analyzed. This results in a long turn around time, which is a major bottleneck in remote sensing applications. Real-time processing and analysis are desired to provide immediate products to support fast decision-making [1–5]. It is highly desirable that onboard processing can be performed [6]. Thus, in our research, real-time data analysis is started right after data are received onboard. In other words, intermediate results are available before data collection is completed. In real-time data analysis, computation is timeconsuming because intermediate results are generated continuously. In this research, we will present a useful approach to reduce computation time such that hardware implementation can be significantly simplified, as preferred for onboard processing. It is worth mentioning that our research is focused on detection and classification, and the algorithms we implement in real-time include the computation of data correlation matrix R or covariance matrix R and their inverses, R-1 or R-1. Therefore, their real-time implementation becomes how to update the R-1 or R-1 as pixels being received. Although some studies have been presented [4, 5], we will focus on the reduction of computational complexity in real-time implementation in this paper. We will also provide different implementation strategies for data with different formats: band interleaved by pixels (BIP), band interleaved by lines (BIL), and band sequential (BSQ) [7].
123
274
J Real-Time Image Proc (2009) 4:273–286
The contributions of this paper include: 1.
2. 3.
the development of real-time implementation strategies for any detection or classification algorithms using data whitening and matched filters; the development of real-time implementation strategies for any remote sensing data formats; the development of an efficient strategy for reducing computational complexity.
This paper is organized as follows. Section 2 briefly introduces the popular detection and classification algorithms for hyperspectral imagery. Section 3 provides realtime implementation strategies for data in three different formats. Section 4 proposes an approach to save the computational complexity. Section 5 presents the system diagram. Section 6 draws the conclusions.
2 Detection and classification algorithms Assume a hyperspectral image has N pixels and L spectral bands. Let a pixel vector be denoted as x, and the entire data matrix X = [x1, x2,…, xN]. Let the sample correlation P matrix R be defined as R ¼ N1 Ni¼1 xi xTi ; which can be related to the sample covariance matrix R and sample mean l by R = R - l lT. If the target to be detected is known, it is denoted as d; if the p targets to be classified are known, they are denoted as D = [d1, d2,…, dp]. Seven algorithms are briefly summarized as follows. They have the common features: (1) using the term R-1 or R-1 for data whitening; (2) using the matched filter dT or DT for detection or classification. It is noteworthy that if the desired target is unknown, then the pixel itself xT acts as the term to be matched, which is the idea in the wellknown Reed–Xiaoli (RX) algorithm for anomaly detection. These algorithms are our research focus because: (1) they can outperform other existing algorithms due to their excellent performance in background suppression; (2) they are suitable to remote sensing images in an unknown circumstance since they require least prior information; (3) they can be easily implemented in real-time with simple and similar hardware architecture. •
RX algorithm [8, 9]:
The well-known RX algorithm [8] is a detector for anomalies, defined as pixels whose spectral signatures are different from their surroundings. The global detector can be written as wRX = xTR-1x [9]. •
Constrained energy minimization (CEM) [10]:
The CEM detector can be written as wCEM = (dTR-1d)-1R-1d. To detect if d is contained in a pixel x, we can simply apply wTCEMx, i.e., ðdT R1 dÞ1 dT R1 x:
123
•
Kelly’s generalized likelihood ratio test (KGLRT) [11]:
This generalized likelihood ratio test is given by T 1 2 d R x T 1 ; d R d 1 þ xT R1 x=N 0 where N’ is the number of samples used in the estimation of R. •
Adaptive matched filter (AMF) [12]:
When the number of samples N is a very large value, the KGLRT is reduced to a simple format: ðdT R1 dÞ1 ðdT R1 xÞ2 : We can see that it is close to the CEM except that the numerator has a power of two. •
Adaptive coherence estimator (ACE) [13]:
The estimator can be written as T 1 2 d R x T 1 : d R d xT R1 x It is similar to AMF except that a term similar to the RX algorithm is included in the denominator. •
Target-constrained (TCIMF) [14]:
interference-minimized
filter
The filter can be expressed as ðDT R1 DÞ1 ðDT R1 xÞ; where p dimensions correspond to p detection maps. It can not only detect targets but also separate them from each other. •
Constrained linear discriminant analysis (CLDA) [15]:
The classifier can be simply put as: ðDT R1 DÞ1 ðD R1 xÞ: It is similar to TCIMF except that R-1 is used rather than R-1. Its major advantage over the traditional classifiers, such as maximum likelihood classifier, is that it does not need to know the complete class information in the data [16]. In other words, it needs to know only the p targets to be classified, even when p is much less than the actual number of classes in the data. T
3 Real-time implementation In our research, we assume that an image is acquired from left to right and from top to bottom. Three real-time processing fashions will be discussed to fit the three remote sensing data formats: pixel-by-pixel processing for BIP format, line-by-line processing for BIL format, and band-by-band processing for BSQ formant. In the pixelby-pixel fashion, a pixel vector is processed right after it is received and the analysis result is generated within an acceptable delay; in the line-by-line fashion, a line of pixel vectors is processed after the entire line is received;
J Real-Time Image Proc (2009) 4:273–286
275
in the band-by-band fashion, a band is processed after it is received. To speed up the process and simplify hardware implementation, one strategy is to replace the data covariance matrix R with the data correlation matrix R as discussed in [4, 5], as both provide similar results. Thus, in order to implement the detection or classification algorithms in Sect. 2 in real-time, the key becomes the adaptation of R-1. In other words, R-1 at time t should be quickly calculated by updating the previous R-1 at time t - 1 using the data received at time t, without completely recalculating the R and R-1. As a result, the intermediate data analysis result (e.g., target detection or classification) is available in support of decision-making even when the entire data set is not received; and when the entire data set is received, the data analysis result is finalized.
~ 1 xt Þ; which In (3), L2 multiplications are needed for ðR t1 ~ 1 Þ; ðxT R ~ 1 xt Þ requires is the same as that for ðxTt R t t1 t1 additional L multiplications only; determining the product ~ 1 xt Þ1 ; and xT R ~ 1 ; ~ 1 xt ; ð1 þ xT R of the three terms: R
This format is easy to handle because a pixel vector of size L 9 1 is received continuously, which well fits the spectral-analysis based algorithms previously discussed in Sect. 2. ~ denotes NR, which is the correlation matrix Let R without being normalized by the number of pixels. Using ~ 1 in the algorithms in Sect. 2 does not change the result, R but simplifies the update process. Supposed that at time t a pixel vector, denoted as xt, is received. The data matrix Xt including all the pixels received up to the time t is Xt = [x1, x2,…, xt]. The sample correlation matrix at time t and t – 1 are denoted as Rt and Rt-1, respectively, and they are related as
and
The following Woodbury’s formula [17] can be used to ~ 1 : update R t 1 ðA þ BCDÞ1 ¼ A1 A1 B C1 þ DA1 B DA1 ð2Þ where A and C are two positive-definite matrices, and the sizes of matrices A, B, C and D allow the inverse operation of (A ? BCD). It should be noted that (2) is for the most general case. Actually, A, B, C, and D can be reduced to vector or scalar as long as (2) is applicable. ~ 1 in (1) can be By using the Woodbury’s formula, R t updated as: 1 ~ 1 ¼ R ~ 1 xt ~ 1 : ~ 1 R ~ 1 xt 1 þ xT R R xTt R ð3Þ t t1 t1 t1 t t1 ~ 1 xt Þ in (3) is a scalar. This means no Here ð1 þ xTt R t1 matrix inversion is involved in each adaptation. Equation (3) can be initiated with a small random matrix ~ 1 [4]. R 0
t1
If the data are in BIL format, we can simply wait for all the pixels in a line is received to form pixel vectors. Let M be the total number of pixels in each line. Let the line received at time t form a data matrix Yt = [xt1xt2,…, xtM]. Then (1) and (3) become ~t ¼ R ~ t1 þ Yt YT R t
ð1Þ
t
t1
3.2 BIL format
3.1 BIP format
~t ¼ R ~ t1 þ xt xT : R t
t
t1
requires L2 multiplications individually. Thus the total number of multiplications in the update for each pixel vector using (3) is 2L2 ? L. For the entire image, this number becomes N(2L2 ? L).
ð4Þ
1 ~ 1 ¼ R ~ 1 ~ 1 ~ 1 R ~ 1 Yt IMM þ YT R R YTt R t t1 t1 t1 t t1 Yt
ð5Þ
respectively, where IM9M is an M 9 M identity matrix. ~ 1 Yt Þ in (5) is a matrix. This means matrix ðIMM þ YTt R t1 inversion is involved in each adaptation. ~ 1 Yt Þ; In (5), ML2 multiplications are needed for ðR t1 T ~ 1 T ~ 1 which is the same as ðY R Þ; ðY R Yt Þ requires t
t1
t
t1
additional M2L multiplications only; determining the ~ 1 Yt Þ1 generally requires O(M3) inverse ðIMM þ YT R t
t1
multiplications [18]; determining the product of the three 1 ~ 1 Yt ; ðIMM þ YT R ~ 1 ~ 1 in (5), terms: R and YTt R t1 t1 t t1 Yt Þ requires (ML2 ? M2L) multiplications. Therefore, the total multiplications in the update for each line is 2(ML2 ? M2L) ? O(M3). For the entire image, this number becomes (2(ML2 ? M2L) ? O(M3))N/M = (2NL2 ? 2NML) ? O(M3)N/M. 3.3 BSQ format If the data format is BSQ, the sample correlation matrix R and its inverse R-1 have to be updated in a different way, because there is no completed pixel vector available before the entire data set is received. Let r1 denote the auto-correlation of pixels in Band 1, which is a scalar, calculated by the average of pixel squared values in Band 1. Let R2 be the correlation matrix between the first two bands. Then r1 can be related to R2 as: r1 r12 R2 ¼ ; r21 r2 where r2 is the auto-correlation of pixels in Band 2, r12 = r21 is the cross-correlation of corresponding pixel
123
276
J Real-Time Image Proc (2009) 4:273–286
values in Band 1 and 2. Therefore, Rt can be related to Rt-1 as Rt1 rt1;t Rt ¼ T ð6Þ rt1;t rt;t where rt is the auto-correlation of pixels in Band t, rt1;t ¼ ½r1;t ; ; rj;t ; ; rt1;t T is a (t - 1) 9 1 vector with rj,t being the cross-correlation of corresponding pixel values in Band j and t. Equation (6) shows that the dimension of Rt is increased as more bands are received. -1 When R-1 t-1 is available, it can be used to calculate Rt -1 by modifying Rt-1 with rt,t and rt-1,t. The following partitioned matrix inversion formula [17] can be used for this purpose. Let a matrix A be partitioned as " # A11 A12 A¼ : A21 A22 Then its inverse matrix A-1 can be calculated as
" A
1
¼
1 A11 A12 A1 22 A21 1 A22 A21 A1 A21 A1 11 A12 11
Let A11 = Rt-1, A22 = rt,t, A21 = rTt-1,t. According to (7)
A12 = rt-1,t,
1 Rt1 rt1;t rt;t rTt1;t 6 ¼4 1 rt;t rTt1;t R1 rTt1;t R1 t1 rt1;t t1
3 1 1 Rt1 rt1;t rt;t rTt1;t rt1;t rt;t 7 1 5 1 rt;t rTt1;t Rt1 rt1;t
1 T 1 ¼ R1 t1 þ aRt1 rt1;t rt1;t Rt1
ð9Þ
aR1 t1 rt1;t
#
a ð10Þ
Thus, all the elements in (10) can be generated by simple matrix multiplication. Actually, in this case, no
123
ð7Þ
and
1 In (8), a ¼ ðrt;t rTt1;t R1 is a scalar; as for t1 rt1;t Þ 1 T ðRt1 rt1;t rt;t rt1;t Þ ; the Woodbury’s formula in (2) can be used again, i.e., 1 Rt1 rt1;t rt;t rTt1;t 1 1 1 T 1 ¼ R1 rTt1;t R1 t1 þ Rt1 rt1;t rt;t rt1;t Rt1 rt1;t t1
and (8) becomes " 1 1 T Rt1 þ aR1 t1 rt1;t rt1;t Rt1 1 Rt ¼ arTt1;t R1 t1
1 T T 1 1 ðrt;t rTt1;t R1 t1 rt1;t Þ rt1;t Rt1 (since rt-1 Rt-1 does not need to be re-calculated); and additional (t - 1)2 multiplications for ðRt1 rt1;t rt;t rTt1;t Þ1 : Thus, the total number of multiplications in the update for the whole P image is N þ Lt¼2 ðN þ Nðt 1Þ þ 2ðt 1Þ2 þ2ðt 1ÞÞ NL þ ðN þ 2ÞLðL 1Þ=2 þðL 1ÞLð2L 1Þ=3 ¼ NL2 =2 þNL=2 þ 2ðL3 LÞ=3
# 1 A11 A12 A1 A12 A1 22 A21 22 1 A22 A21 A1 11 A12
2
R1 t
operation of matrix inversion is needed when reaching the final R-1. The intermediate result still can be generated by applying the R-1 to the first t bands. This means the t spectral features in these t bands are used for target detection and discrimination. This may help find targets at early processing stages. It is worth mentioning that R-1 is symmetric and three elements in the matrix in (8) need to be calculated. N and N(t - 1) multiplications are needed to calculate rt,t, and rt,t-1, respectively; (t - 1)2 ? (t - 1) multiplications for 1 a ¼ ðrt;t rTt1;t R1 t1 rt1;t Þ ; ðt 1Þ multiplications for
ð8Þ
4 Fast implementation 4.1 Computational complexity reduction for R-1 evaluation In this paper, we will propose another way to save the number of computations. In the offline version, the major computation time is spent in the calculation of the data correlation matrix R. For an image with N pixels and L bands, the number of multiplications for R is NL2. Its matrix inversion generally needs O(L3) multiplications, depending upon the methods employed [18]. For a remote sensing image, N L, then the key to reducing computational complexity is to use a part of pixels for the R calculation. Fortunately, remote sensing image usually has high-spatial correlation; using a small part of pixels can still capture the major data statistics.
J Real-Time Image Proc (2009) 4:273–286
277
Table 1 Computational Complexity in real-time and offline processing of different formats g % pixels are used
All the pixels are used Offline BIP
Real-time
NL2 ? o(L3) 2
3
Offline
N(2L2 ? L) 2
Real-time
g NL2 ? o(L3) 3
NL ? o(L )
(2NL ? 2NML) ? o(M )N/M
g NL ? o(L )
(2gNL2 ? 2gNML) ? o(M3)g N/M
BSQ
NL2 ? o(L3)
NL2/2 ? NL/2 ? 2(L3-L)/3
g NL2 ? o(L3)
g NL2/2 ? gNL/2 ? 2(L3-L)/3
4.2 Experiments 4.2.1 HYDICE experiment: target detection The hyperspectral digital imagery collection experiment (HYDICE) image scene shown in Fig. 1 was collected in Maryland in 1995 from the flight altitude of 10,000 ft with approximately 1.5 m spatial resolution in 0.4–2.5 lm spectral region. The atmospheric water bands with low signal-to-noise ratio were removed reducing the data dimensionality from 210 to 169. The image scene has 128 lines and the number of pixel in each line M is 64. This scene includes 30 panels arranged in a 10 9 3 matrix. Each panel is denoted by pij with row indexed by i = 1,…, 10 and column indexed by j = a, b, c. The three panels in the same row were made from the same material of size 3 m 9 3 m, 2 m 9 2 m, 1 m 9 1 m, respectively, which were considered to belong to the same class, Pi.
3
g N(2L2 ? L)
BIL
Table 1 summarizes the number of multiplications for each data format in real-time and offline processing, where g is the percentage of pixels participating in the R-1 evaluation. Here the offline processing means all the data are processed simultaneously; in other words, the processing is started after an entire data set is received. We can see that the computational complexity in offline processing of BIP and BIL format is much less than that in their realtime processing when all the pixels are used. However, when image size is small, it is possible for the real-time processing of BIP format to have less complexity if g is less than a certain value. When using g% of pixels, the computation reduction factor for all the formats is about g. In Sect. 4.2 we will show that when g is appropriately selected, the detection or classification performance will remain unchanged. When the image scene contains large homogeneous areas, then g can be small; otherwise, g should be greater such that the data statistics can be manifested by a portion of pixels. An empirical rule is to select five to ten times of the number of spectral dimensionality (i.e., the number of spectral bands). The low bound is defined such that the number of linearly independent pixels is equal to the number of spectral bands for the existence of R-1 (i.e., an R with full rank).
2
P1a, P1b, P1c P2a, P2b, P2c P3a, P3b, P3c P4a, P4b, P4c P5a, P5b, P5c P6a, P6b, P6c P7a, P7b, P7c P8a, P8b, P8c P9a, P9b, P9c P10a, P10b, P10c Fig. 1 HYDICE image scene with 30 panels
Figures 2 and 3 show the detection results using TCIMF when all the pixels and 10% pixels were involved in R-1 calculation, respectively. The difference is not perceivable. Table 2 lists the number of panel pixels correctly detected ND and the number of false alarm pixels NF when all the pixels, 10% pixels, and 5% pixels were used. The total number of panel pixels NP is 38. When all the pixels and 10% pixels were used, 28 pixels out of 38 could be detected without false alarm pixels. When 5% pixels were used, 26 pixels were detected with 4 false alarm pixels. It seems that only 10% pixels are needed for this image scene. Since this image scene has only 8,192 pixels. When 10% pixels are used, the spatial size is about 5 times of the spectral size. Real-time processing results are tabulated in Tables 3, 4 and 5. In Table 3 for pixel-by-pixel fashion, the result was the same as the one in offline processing when all the pixels were used for correlation matrix update; if only 10% pixels participated, then the result was slightly degraded with ND = 29 and NF = 13; if only 5% pixels were used, then the result was greatly degraded with NF = 303. In Table 4 for line-by-line processing, the result was the same as the one in offline processing when all the pixels were involved in the correlation matrix update; the performance was slightly better than the case of pixel-by-pixel since NF was reduced to 29 with 5% pixels. In Table 5 for band-by-band processing, the result was the same as the offline processing since the results can be compared only after all the bands are received. Based on this computer simulation, it seems that using 10% of pixels is a good choice for real-time processing.
123
278
J Real-Time Image Proc (2009) 4:273–286
Fig. 2 Detection result using TCIMF with all the pixels for R-1 calculation
Fig. 3 Detection result using TCIMF with 10% pixels for R-1 calculation
4.2.2 AVIRIS Cuprite experiment: classification The airborne visible/infrared imaging spectrometer (AVIRIS) Cuprite subimage scene of size 350 9 350
123
shown in Fig. 4 was collected in Nevada in 1997. The spatial resolution is 20 m. Originally it has 224 bands with 0.4–2.5 lm spectral range. After water absorption and lowSNR bands were removed, 189 bands were left. This image
J Real-Time Image Proc (2009) 4:273–286
279
Table 2 Performance comparison for HYDICE experiment 10%
Table 5 Real-time processing performance in HYDICE experiment (band-by-band)
Pixel percentage
100%
Panel
NP
ND
P1
3
2
0
2
0
2
0
P2
3
2
0
2
0
2
2
P3
4
3
0
4
0
3
0
NF
ND
5% NF
ND
NF
Pixel percentage
100%
10%
Panel
NP
ND
P1
3
2
0
2
0
2
0
P2
3
2
0
2
0
2
2
4 3
3 2
0 0
4 2
0 0
3 2
0 0
NF
ND
5% NF
ND
NF
P4
3
2
0
2
0
2
0
P5
6
5
0
6
0
4
0
P3 P4
P6
3
2
0
2
0
1
0
P5
6
5
0
6
0
4
0
P7 P8
4 4
3 3
0 0
3 3
0 0
3 3
2 0
P6
3
2
0
2
0
1
0
P7
4
3
0
3
0
3
2
0
P8
4
3
0
3
0
3
0
0
P9
4
3
0
3
0
3
0
4
P10
4
3
0
3
0
3
0
Total
38
28
0
28
0
26
4
P9 P10 Total
4 4 38
3 3 28
0 0 0
3 3 28
0
3
0
3
0
26
Table 3 Real-time processing performance in HYDICE experiment (pixel-by-pixel) Pixel percentage
100%
Panel
NP
ND
P1
3
2
P2
3
2
P3 P4
4 3
P5
10% NF
5%
ND
NF
ND
NF
0
2
0
2
0
0
2
0
3
41
3 2
0 0
3 2
0 0
2 3
201 1
6
5
0
5
0
5
0
P6
3
2
0
2
0
2
0
P7
4
3
0
4
13
4
59
P8
4
3
0
3
0
3
0
P9
4
3
0
3
0
3
1
P10
4
3
0
3
0
3
0
Total
38
28
0
29
13
30
303
Table 4 Real-time processing performance in HYDICE experiment (line-by-line) Pixel percentage
100%
Panel
NP
ND
P1
3
3
P2
3
P3
4
P4 P5
10% NF
5%
ND
NF
ND
NF
0
3
0
3
0
2
0
2
0
2
0
3
0
3
1
3
27
3
2
0
2
0
2
0
6
5
0
5
0
5
0
P6
3
2
0
2
0
2
0
P7
4
3
0
3
10
3
2
P8
4
3
0
3
0
3
0
P9
4
3
0
3
0
3
0
P10
4
3
0
3
0
3
0
Total
38
29
0
29
11
29
29
Fig. 4 AVIRIS Cuprite scene with five major minerals: alunite (A), buddingtonite (B), calcite (C), kaolinite (K), and muscovite (M)
scene is well understood mineralogically. At least five minerals were present: alunite (A), buddingtonite (B), calcite (C), kaolinite (K), and muscovite (M). Their approximate spatial locations of these minerals are marked in Fig. 4. Figures 5 and 6 show the results using CLDA when all the pixels and 1% pixels were involved in R-1 calculation, respectively. The difference is not perceivable. Since no pixel-level ground truth is available, the spatial correlation coefficient is employed to evaluate the similarity between the corresponding maps. A larger correlation coefficient means better performance. As listed in Table 6, the correlation coefficient is as high as 0.9 even when only 1% pixels was used. This image scene is large and contains 122,500 pixels. When 1% pixels are used, the spatial size is about 6.5 times of the spectral size. The ratio between the number of pixels
123
280
J Real-Time Image Proc (2009) 4:273–286
Fig. 5 Classification result using CLDA with all the pixels for R-1 calculation
Fig. 6 Classification result using CLDA with 1% pixels for R-1 calculation Table 6 Performance comparison using correlation coefficient for AVIRIS Cuprite experiment
Table 8 Real-time processing performance in AVIRIS Cuprite experiment (line-by-line)
Pixel percentage
10%
5%
1%
Pixel percentage
100%
10%
5%
1%
Alunite
0.9926
0.9827
0.8953
Alunite
0.8944
0.8802
0.8719
0.8306
Buddingtonite
0.9926
0.9847
0.8850
Buddingtonite
0.9484
0.9352
0.9274
0.8224
Calcite
0.9936
0.9861
0.9121
Calcite
0.9433
0.9174
0.8920
0.6930
Kaolinite
0.9903
0.9794
0.8711
Kaolinite
0.9430
0.9164
0.9190
0.8364
Muscovite
0.9905
0.9826
0.9362
Muscovite
0.9229
0.9179
0.8943
0.8156
Average
0.9919
0.9831
0.8999
Average
0.9304
0.9134
0.9009
0.7996
Table 7 Real-time processing performance in AVIRIS Cuprite experiment (pixel-by-pixel)
Table 9 Real-time processing performance in AVIRIS Cuprite experiment (band-by-band)
Pixel percentage
100%
10%
5%
1%
Pixel percentage
100%
10%
5%
1%
Alunite
0.8935
0.8843
0.8706
0.8082
Alunite
1.0000
0.9926
0.9827
0.8953
0.7321
Buddingtonite
1.0000
0.9926
0.9847
0.8850
1.0000
0.9936
0.9861
0.9121 0.8711
Buddingtonite
0.9477
0.9254
0.8655
Calcite
0.9422
0.8848
0.6456
0.6280
Calcite
Kaolinite
0.9424
0.9255
0.8836
0.8266
Kaolinite
1.0000
0.9903
0.9794
Muscovite
0.9218
0.8906
0.8831
0.6917
Muscovite
1.0000
0.9905
0.9826
0.9362
Average
0.9295
0.9021
0.8297
0.7373
Average
1.0000
0.9919
0.9831
0.8999
and the number of bands (189) may be an indictor for the appropriate percentage that should be chosen. Real-time processing results are tabulated in Tables 7, 8 and 9. We can see that the performance of line-by-line processing was better than the pixel-by-pixel processing; using 5% pixels the classification result was not
123
significantly changed in line-by-line processing and using 10% pixels the classification result was not significantly changed in pixel-by-pixel processing. Once again, the classification result for the BSQ format was the same as the offline processing. This means the pixel percentage for data statistics evaluation may need to be increased when an algorithm is implemented in real-time.
J Real-Time Image Proc (2009) 4:273–286
281
4.2.3 AVIRIS Lunar Lake experiment: anomaly detection
Table 10 Performance comparison using correlation coefficient for AVIRIS Lunar Lake experiment
The 224-band AVIRIS Lunar Lake image scene shown in Fig. 7 was used in the experiment. It is a subscene of 200 9 200 pixels cropped from the left upper corner of the Lunar Crater Volcanic Field in Northern Nye County, Nevada. Five classes are present: ‘‘red oxidized basaltic cinders’’, ‘‘rhyolite’’, ‘‘playa (dry lakebed)’’, ‘‘shade’’, and ‘‘vegetation’’. In addition, there is an anomaly class present. The number of bands was reduced to 158 after lowSNR and water absorption bands were removed. The RX algorithm was applied for anomaly detection. As shown in Fig. 8, there is no perceivable difference between the results using all the pixels and 2% pixels for R-1 calculation, where both the anomaly and some part of vegetation were found. Table 10 lists the correlation coefficients between the detection maps using all the pixels and a small portion of pixels. We can see that the correlation coefficient is as high as 0.8869 even when only 2% pixels were used. At this time, the number of pixels was
Pixel percentage
100%
10%
5%
2%
Offline
–
0.9805
0.9545
0.8869
Online (pixel-by-pixel)
1.0000
0.9805
0.9545
0.8869
Online (line-by-line)
1.0000
0.9805
0.9545
0.8869
Online (band-by-band)
1.0000
0.9805
0.9545
0.8869
reduced to 800, which is about 5 times of the number of bands. Real-time processing results are tabulated in Table 10. Since the offline result was the one using global statistics (global anomaly detector), the results were comparable only after all the data were received. We can see that in this case real-time processing for all the three formats provided the same result as in offline processing.
5 System architecture 5.1 System overview The system for fast real-time detection or classification can be summarized as follows. The block diagram is shown in Fig. 9. 1. 2.
Input the data format and choose the value of g. If the format is BIP, • •
• Fig. 7 The AVIRIS Lunar Lake scene of size 200 9 200 used in the experiment
3.
If the format is BIL, • • • •
4.
Initiate R-1 0 using a small random matrix of size L 9 L; Wait until the entire line of pixel vectors is received at time t; Use the pixels selected by the PRS module to update R-1 with (5); t Go to Step 5.
If the format is BSQ, • •
Fig. 8 Anomaly detection result using RX algorithm
Initiate R-1 0 using a small random matrix of size L 9 L; When xt is received at time t, check if it is selected by the pixel random selection (PRS) module; if yes, use (3) to update R-1 t ; if not, keep the previous -1 R-1 t-1 as Rt ; Go to Step 5.
•
Use the PRS module to select the pixels in each band that will participate in R-1 update; When the t-th band is received at time t, compute r1 and rj,t for j = 1,2,…, (t - 1); Update R-1 with (10). t
123
282
J Real-Time Image Proc (2009) 4:273–286
BIP
R t−1 update
usin g p ixels xt
Format Selection
Random Pixel Selection
BIL
R t−1 update
using data lines
Detector/ Classifier Construction
Detection/ Classification
R t−1 update
BS Q
using band correlation vectors
Data Storage Fig. 9 Real-time onboard detection or classification pipeline stages
5. 6.
Construct the detector or classifier in Sect. 2. Provide the detection or classification result for the pixels/lines/bands that have been received up to time t.
5.2 Hardware architecture The major system modules have been discussed in [4], including parallel matrix vector multiplier (PMVM), matrix–matrix multiplier (MMM), and matrix inversion (MI) models. Here, we only discuss the PRS module. We also provide hardware architecture for R-1 for three dift ferent formats, which are very similar. 5.2.1 PRS module Assume 10% pixels need to be randomly selected. The PRS module will provide a bit pattern with approximately 10% of pattern data being 1 s and 90% 0 s, where a value of 1 indicates the corresponding pixel is selected for calculations and a 0 means to ignore the pixel. As shown in Fig. 10, the module includes a random bit generator (RBG) [19] which employs Zener Diode as the source of its randomness to provide approximately 50% 1 s and 50% 0 s, followed by a counter-comparator circuit to select one out of every five 1 s. A simulation of the circuit generated a bit pattern with 10.6% 1 s and 89.4% 0 s is illustrated in
Fig. 10 A module for pixel random selection (PRS)
123
Fig. 10. By slightly modifying the comparator, different percentages can be easily simulated. 5.2.2 Hardware architecture for R-1 update for BIP t format Figure 11 illustrates a parallel pipeline hardware realization of (3) where data are arriving one pixel vector at the ~ 1 is updated each time a new time. In this module R t selected pixel vector is processed. Parallel Multiply Accumulators calculate result of the matrix vector multi~ 1 xt Þ; which is pipelined to calculate the plication, ðR t1 required scalar inner product and outer vector product ~ 1 ¼ R ~ 1 1 T ~ 1 ~ 1 xt ð1 þ xT R matrix: DR t t1 xt Þ xt Rt1 : t t1 5.2.3 Hardware architecture for R-1 update for BIL t format Figure 12 illustrates how BIP architecture can be modified for a real-time line by line approximation. The matrix vector manipulation has also been reordered to take advantage of data arriving one line of each band at a time. This will improve system throughput and the pipeline performance. Additionally, the architecture has been extended with vector multiply accumulators to be able to process the entire line.
J Real-Time Image Proc (2009) 4:273–286
283
~ ~ R t−−11 x t = x Tt R t−−11
Outer Product Select
MAC1
MV1 Reg.
*
MAC2
MV2 Reg.
*
Inner_en MUX
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
D E M U X
. . .
Inner_en
*
MVL Reg.
~ ∆R t−1
MUX MACL
Inner_en
~ R t−−11
xtt
DEMUX
Inner Product Register
MV1,L Reg.
Inverse Operation
* ~ (1 + xTt R t−−11x t ) −1
Fig. 11 Architecture for BIP iterative correlation matrix inversion
5.2.4 Hardware architecture for R-1 update for BSQ t format Figure 13 illustrates the hardware architecture for BSQ format, which is similar to the previous two architectures. The three major terms in (10) are generated, from which the inverse matrix R-1 t can be easily updated. In addition to the R-1 t-1, the input includes the correlation vector at time t which can be easily generated using MAC. It is noteworthy that the dimensionality of R-1 is increased as more t bands being received; as a result, some MACs and
multipliers in early stages are actually not used with inputs and outputs being zero.
6 Conclusion In this paper, we investigate real-time onboard processing of hyperspectral imagery for detection or classification. The detectors and classifiers we employ include the term of the inverse of data correlation or covariance matrix (R-1or R-1) for data whitening, followed by a matched filter for
123
284
J Real-Time Image Proc (2009) 4:273–286
Column Select
V_MAC1 c1 c2 : cM
+
*
~ YtT R t−−11
Outer Product Select
MAC
MV1 Reg.
Inner_en MAC
MUX
MV2 Reg.
V_MAC2
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
D E M U X
. . .
MAC
Inner_en MUX
MVL Reg.
V_MACL
Yt
DEMUX
Inner_en
~ R t−−11
~ ∆R t−1
Inner Product Register
*
MV1,L Reg.
Inverse Operation
(I + Y
T t
~ R t−−11 Yt
)
−1
Fig. 12 Architecture for BIL iterative correlation matrix inversion
the target(s) to be detected or classified. In real-time processing, data processing and analysis start immediately after pixels are received. Thus, the key is to update R-1 or R-1 in real-time, which has been investigated when data are formatted as BIP, BIL, or BSQ. To speed up the process and facilitate onboard hardware implementation, we propose the use of a small portion of pixels in the update of R-1. Experimental results show that with a small portion of pixels, comparable detection or classification results can be provided. For data with large
123
homogeneous areas, only a small number of pixels need to be used. To be conservative, in the offline processing, the percentage of selected pixels can be empirically chosen such that the number of pixels is about five to ten times the number of bands (which is coincident with the observation in [20] for training set size selection in speech signal processing); this percentage for real-time implementation may need to be increased to about 5 times the one in offline processing. We also provide the system architecture for the proposed fast real-time algorithms for different data formats.
J Real-Time Image Proc (2009) 4:273–286
285
(r
T t −1,t
R t−−11
)
Outer Product Select
MAC1
MV1 Reg.
*
MAC2
MV2 Reg.
*
Inner_en MUX
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
D E M U X
. . .
Inner_en
*
MVL Reg.
αR t−−11rt −1rtT−1R t−−11
MUX MACL
R
rt −1,t
−1 t −1
DEMUX
Inner_en
Inner Product Register
MV1,L Reg.
* αrtT−1,t R t−−11
Inverse Operation
α = (rt ,t − rtT−1,t R t−−11rt −1,t )
−1
Fig. 13 Architecture for BSQ iterative correlation matrix inversion
Similarities among BIP, BLP, and BSQ hardware architectures allow the use of a single parallel pipelined computation engine for all three data formats and reduce hardware complexity for a general purpose system.
References 1. Stellman, C.M., Hazel, G.G., Bucholtz, F., Michalowicz, J.V.: Real-time hyperspectral detection and cuing. In: Optical Engineering, vol. 39, no. 7, July (2000) 2. Chang, C.-I., Ren, H., Chiang, S.S.: Real-time processing algorithms for target detection and classification in hyperspectral imagery. In: IEEE Transactions on Geoscience and Remote Sensing, vol. 39, No. 4 (2001)
3. Du, Q., Ren, H.: Real-time constrained linear discriminant analysis to target detection and classification in hyperspectral imagery. In: Pattern Recognition, vol. 36, No. 1 (2003) 4. Du, Q., Nekovei, R.: Implementation of real-time constrained linear discriminant analysis to remote sensing image classification. In: Pattern Recognition, vol. 38, No. 4 (2005) 5. Du, Q.: Unsupervised real-time constrained linear discriminant analysis to hyperspectral image classification. In: Pattern Recognition, vol. 40, No. 5 (2007) 6. Subramanian, S., Gat, N., Ratcliff, A., Eismann, M.: Real-time hyperspectral data compression using principal components transformation. In: Proceedings of the AVIRIS Earth Science and Applications Workshop (2000) 7. Jensen, J.R.: Introductory Digital Image Processing: A Remote Sensing Perspective, 3rd edn. Prentice-Hall, Englewood Cliffs (2004) 8. Reed, I.S., Yu, X.: Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution. In: IEEE
123
286
9.
10.
11. 12.
13.
14.
15.
16.
17. 18.
J Real-Time Image Proc (2009) 4:273–286 Transactions on Acoustic, Speech and Signal Processing, vol. 38, No. 10 (1990) Chang, C.-I., Heinz, D.: Subpixel spectral detection for remotely sensed images. In: IEEE Transactions on Geoscience and Remote Sensing, vol. 38, No. 3 (2000) Farrand, W.H., Harsanyi, J.C.: Mapping the distribution of mine tailing in the coeur d’Alene river valley, Idaho through the use of constrained energy minimization technique. In: Remote Sensing of Environment, vol. 59, No. 1 (1997) Kelly, E.J.: An adaptive detection algorithm. In: IEEE Transactions on Aerospace and Electronic Systems, vol. 22, No. 1 (1986) Robey, F.C., Fuhrmann, D.R., Kelly, E.J., Nitzberg, R.: A CFAR adaptive matched filter detector. In: IEEE Transactions on Aerospace and Electronic Systems, vol. 28, No. 1 (1992) Kraut, S., Sharf, L.L.: The CFAR adaptive subspace detector is a scale-invariant GLRT. In: IEEE Transactions on Signal Processing, vol. 47, No. 9 (1999) Ren, H., Chang, C.-I.: A target-constrained interference-minimized approach to subpixel detection for hyperspectral images. In: Optical Engineering, vol. 39, No. 12 (2000) Du, Q., Chang, C.-I.: Linear constrained distance-based discriminant analysis for hyperspectral image classification. In: Pattern Recognition, vol. 34, No. 2 (2001) Du, Q.: Modified Fisher’s linear discriminant analysis for hyperspectral imagery. In: IEEE Geoscience and Remote Sensing Letters, vol. 4, No. 4 (2007) Golub, G.H., Van Loan, C.F.: Matrix Computations, 3rd edn. Johns Hopkins University Press (1996) Dremmel, J.E.: Applied Numerical Linear Algebra. Society for Industrial and Applied Mathematics, Philadelphia (1997)
123
19. Stipcevic, M.: Fast nondeterministic random bit generator based on weakly correlated physical events. In: Review of Scientific Instruments, vol. 75, No. 11 (2004) 20. Makhoul, J., Roucos, S., Gish, H.: Vector quantization in speech coding. In: Proceedings of the IEEE, vol. 73, No. 11 (1985)
Author Biographies Qian Du received her Ph.D. degree in electrical engineering from University of Maryland Baltimore County in 2000. She was with the Department of Electrical Engineering and Computer Science, Texas A&M University-Kingsville, from 2000 to 2004. She joined the Department of Electrical and Computer Engineering at Mississippi State University in Fall 2004, where she is currently an Associate Professor. Her research interests include remote sensing image analysis, pattern classification, data compression, and neural networks. Dr. Du is a senior member of IEEE, a member of SPIE, ASPRS, and ASEE. Reza Nekovei received his Ph. D. degree in Electrical Engineering from The University of Rhode Island in 1994. He is currently an Associate Professor at the Department of Electrical Engineering and Computer Science, Texas A&M University-Kingsville. He was previously with the Remote Sensing Laboratory, Graduate School of Oceanography at The University of Rhode Island. His research interests include VLSI Architectures, Distributed Systems, Remote Sensing, and Medical Imaging. He is a member of IEEE and Etta Kappa Nu honor society.
