Real-time constrained linear discriminant analysis ... - Semantic Scholar
Pattern Recognition 36 (2003) 1 – 12
www.elsevier.com/locate/patcog
Real-time constrained linear discriminant analysis to target detection and classi cation in hyperspectral imagery Qian Dua; ∗ , Hsuan Renb a Department
of Electrical Engineering and Computer Science, MSC 192, Texas A& M University-Kingsville, TX 78363, USA Biological and Chemical Center, US Army Aberdeen Proving Ground, MD 21010, USA
b Edgewood
Received 11 April 2001; received in revised form 2 November 2001; accepted 6 March 2002
Abstract In this paper, we present a constrained linear discriminant analysis (CLDA) approach to hyperspectral image detection and classi cation as well as its real-time implementation. The basic idea of CLDA is to design an optimal transformation matrix which can maximize the ratio of inter-class distance to intra-class distance while imposing the constraint that di2erent class centers after transformation are along di2erent directions such that di2erent classes can be better separated. The solution turns out to be a constrained version of orthogonal subspace projection (OSP) implemented with a data whitening process. The CLDA approach can be applied to solve both detection and classi cation problems. In particular, by introducing color for display the classi cation is achieved with a single classi ed image where a pre-assigned color is used to display a speci ed class. The real-time implementation is also developed to meet the requirement of on-line image analysis when the immediate data assessment is critical. The experiments using HYDICE data demonstrate the strength of CLDA approach in discriminating the targets with subtle spectral di2erence. ? 2002 Pattern Recognition Society. Published by Elsevier Science Ltd. All rights reserved. Keywords: Detection; Classi cation; Real-time processing; Constrained linear discriminant analysis (CLDA); Hyperspectral imagery
1. Introduction Hyperspectral imaging is a new technique in remote sensing, which generates hundreds of images with di2erent spectral channels for the same area on the earth. Because of very high spectral resolution, it is most likely to better discriminate di2erent materials using hyperspectral images than using multispectral images, which contains only tens of spectral channels. For instance, the airborne visible=infrared imaging spectrometer (AVIRIS) developed by NASA’s Jet Propulsion Laboratory has 224 spectral channels with 10 nm spectral resolution and 20 m spatial resolution while the hyperspectral digital image collection experiment (HYDICE) sensor developed by Naval Research Laboratory has 210 spectral channels with the same 10 nm spectral resolution and 1–4 m spatial resolution. AVIRIS is mainly used for ∗ Corresponding author. Tel.: +1-361-5932625; fax: +1-3615932110. E-mail address: [email protected] (Q. Du).
ecology, geology, global climate monitoring, land mapping, etc., while HYDICE can be used for surveillance and reconnaissance purposes because of its ner spatial resolution. One of the major challenging problems about hyperspectral image analysis is how to deal with the massive data volume while taking full advantage of plentiful spectral information [1]. Moreover, the size of targets generally is smaller than the spatial resolution or the area covered by a pixel, so image analysis in remote sensing actually copes with mixed pixel processing instead of pure pixel processing in standard digital images. In general, new algorithm development for hyperspectral imagery is necessary. Fisher’s linear discriminant analysis (LDA) is a standard technique in pattern recognition. It projects the original high dimensional data onto a low dimensional space, where all the classes are well separated by maximizing the Raleigh quotient, i.e., the ratio of between-class scatter matrix determinant to within-class scatter matrix determinant [2]. The application of Fisher’s LDA to hyperspectral images has been investigated for hyperspectral image classi cation
0031-3203/02/$22.00 ? 2002 Pattern Recognition Society. Published by Elsevier Science Ltd. All rights reserved. PII: S 0 0 3 1 - 3 2 0 3 ( 0 2 ) 0 0 0 6 5 - 1
2
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
by Du and Chang in Ref. [3]. The derived algorithm in Ref. [3] is referred to as constrained linear discriminant analysis (CLDA) algorithm in this paper. In CLDA the original high dimensional data is projected onto a low dimensional space as done by Fisher’s LDA. But di2erent classes are forced to be along di2erent directions in this low dimensional space. Thus all classes were expected to be better separated and the classi cation is achieved simultaneously with the CLDA transform. The transformation matrix in CLDA maximizes the ratio of inter-class distance to intra-class distance while satisfying the constraint that the means of di2erent classes are aligned with di2erent directions, which can be constructed by using an orthogonal subspace projection (OSP) method coupled with a data whitening process. The experimental results in Ref. [3] demonstrated that the CLDA algorithm could provide more accurate classi cation results than other popular methods in hyperspectral image processing, such as OSP classi er [4] and constrained energy minimization (CEM) operator [5]. Recently, a new approach, referred to as lter vector algorithm (FVA) was proposed by J. Bowles et al. and its advantage in hyperspectral image classi cation was validated in Ref. [6]. In our experiments we found out that the CLDA algorithm outperforms the FVA. Interestingly, the FVA can be viewed as a special case of the CLDA algorithm [7]. Although the CLDA algorithm showed its strength in hyperspectral image processing, further improvement can be made in two aspects. One is that it was only used for classi cation in Ref. [3], not detection. Another is that the p images are needed to produce by repeating the same procedure p times, one for each target classi cation. Therefore, in this paper we will extend and modify the CLDA algorithm to accomplish target detection and classify all the targets in a single image. Since recently the spectral and spatial resolutions of hyperspectral images are signi cantly improved, their applications have been expanded to battle eld and national defense area where the real-time data processing becomes inevitable so as to process data on-line and provide the information for immediate response [8,9]. In order to meet this demand, a real-time version of the CLDA algorithm is also developed. The remainder of this paper is organized as follows. Section 2 brieJy reviews the CLDA algorithm. Section 3 describes the target detector based on the CLDA algorithm. Section 4 focuses on the CLDA target classi er. In Section 5, the real-time implementation of the CLDA-based target detection and classi cation algorithm is presented. In Section 6, HYDICE data experiments are used to demonstrate the performance of the CLDA detector and classi er as well as their real-time implementation. Finally, Section 7 concludes this paper with several remarks. 2. Constrained linear discriminant analysis Assume that there are c classes and the kth class contains Nk patterns. Let N = N1 + · · · + Nc be the total number of
training patterns. The jth pattern in the kth class, denoted k k k T by xjk = [x1j ; x2j ; : : : ; xLj ] is an L-dimensional vector. Let c Nk k \ = 1=N k=1 j=1 xj be the global mean of training sam k k ples and \k = 1=Nk Nj=1 xj the mean of the kth class. We de ned J (F) to be the ratio of the inter-class distance to the intra-class distance after a linear transformation F, which is given by c 2 2=c(c − 1) c−1 i=1 j=i+1 F(\i ) − F(\j ) J (F) = (1) c Nk k 1=cN k=1 [ j=1 F(xj ) − F(\k )2 ] and F(x) = (WL×c )T x = [w1 ; w2 ; : : : ; wc ]T x;
(2)
where wi is an L × 1 vector. The optimal linear transformation F ∗ is the one that maximizes J (F) subject to tk = F(\k ) for all k; where
(3)
T
tk = 0 · · · 0 1 0 · · · 0 k
is a c × 1 unit column vector with one in the kth component and zeros elsewhere. F ∗ can be determined by wi∗ = \ˆ Ti P⊥ ˆi U
(4)
where ˆ ˆ T ˆ −1 U ˆ Ti P⊥ ˆ i = I − Ui (Ui Ui ) U
(5)
with Uˆ i = [ uˆ 1 · · · uˆ j · · · uˆ c ]j=i and I the identity matrix. The “∧” speci es the whitened data, i.e., xˆ = PwT x, where Pw is the data whitening operator and can be calculated by PwT = V2−1=2 V1T :
(6)
Here, V1 and V2 are eigenvector and eigenvalue matrices of the data sample covariance matrix = 1=N Ni=1 T (xi − \)(xi − \) , respectively. They can be determined by the eigendecomposition = V1 V2 V1T , where V1 = [v1 ; v2 ; : : : ; vL ] with vi being the ith eigenvector and V2 is a diagonal matrix with the ith diagonal item being the ith eigenvalue corresponding to the ith eigenvector, i.e., V2 = diag (1 ; 2 ; : : : ; L ).
3. Constrained linear discriminant analysis-based target detector Assume that an image scene contains p desired targets, say, man-made targets and q undesired targets, say, background objects. For derivation purpose, we treat the detection problem as a two-class classi cation problem: the p
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
desired targets as class !D with desired class mean matrix D=\1 \2 · · · \p and the q undesired targets as another class !U with undesired class mean matrix U = \p+1 \p+2 · · · \p+q . We have p + q = c, where c is the total number of di2erent target classes. The patterns in class !D are {xD } = {xk }pk=1 and the number of patterns ND in !D equals ND = N1 + N2 + · · · + Np . The patterns in class !U are {xU } = {xk }p+q k=p+1 and the number of patterns NU in !U equals NU = Np+1 + Np+2 + · · · + Np+q . The sample means \D and can be calculated by \D = \U of these two classes 1=ND pi=1 Ni \i and \U = 1=NU p+q i=p+1 Ni \i , respectively. For this two-class classi cation problem, Eq. (1) becomes J (F) = F(\D )−F(\U )2 ;
U D 1=2N Nj=1 F(xjD )−F(\D )2 + Nj=1 F(xjU )−F(\U )2 (7) where F(x) = wT x and w is an L × 1 transformation vector. Our goal is to minimize Eq. (7) while imposing the constraints w T D = 1T
(8)
and
3
The algorithm of CLDA-based target detection is presented as follows. (1) Identify the desired and undesired targets and construct D and U. (2) Calculate the data covariance matrix and generate the data whitening operator Pw via eigendecomposition of . (3) Whiten the data by applying Pw to the data matrix X as well as the desired and undesired target matrices, D ˆ Dˆ and U. ˆ and U, resulting in X, (4) The detection of desired targets can be achieved by PdT Xˆ with PdT constructed using Eq. (10). It should be noted that if D and U are determined by an unsupervised method in step (1) using the techniques in Ref. [10] or [11], this algorithm can be extended to an unsupervised version. In Appendix B, we prove that 1T (DT PU⊥ D)−1 DT PU⊥ = [1T
0T ](ST S)−1 ST ;
(11)
where S = [D U] is the entire target matrix. If the data is whitened, Eq. (11) becomes T ˆ −1 Dˆ T P⊥ˆ = [1T PdT = 1T (Dˆ PU⊥ˆ D) U
T ˆ −1 Sˆ T : 0T ](Sˆ S)
(12)
Since Sˆ = PwT S and Xˆ = PwT X,
wT U = 0T
T ˆ −1 Sˆ T Xˆ 0T ](Sˆ S)
(9)
PdT Xˆ = [1T
where 1 is a p ×1 column vector with all the elements equal to 1, i.e., 1=(1; 1; : : : ; 1)T and 0 is a q×1 column vector with
= [1T
0T ](ST Pw PwT S)−1 ST Pw PwT X
= [1T
P Td X; 0T ](ST −1 S)−1 ST −1 X = P
p
T
all the elements equal to 0, i.e., 0 = (0; 0; : : : ; 0) . If Eqs. (8)
where
and (9) are satis ed, it is easy to show that wT \D = 1 and wT \U = 0, which means the constraint in Eq. (3) is relaxed. The advantage of using Eqs. (8) and (9) as constraints is that \D and \U are unnecessary to be calculated. We notice that the numerator in Eq. (7) is a constant when the constraints are satis ed. Thus maximizing J (F) is equivalent to minimizing the denominator, which is equal to (1=2N )wT SW w with SW denoting the with-class scatter matrix de ned in Ref. [2]. Since SW + SB = ST = N · , minimizing (1=2N )wT SW w is equivalent to minimizing wT w, where SB and ST are the between-class scatter matrix and the total scatter matrix de ned in Ref. [2], respectively. After the data is whitened by Pw in Eq. (6), becomes the identity matrix. Hence the constrained detection problem can be described as: to nd a transformation vector w such that w2 is minimum with the constraints that wT D = 1 and wT U = 0. In Appendix A we derive the solution w∗ of this constrained detection problem, referred to as CLDA-based target detector Pd
P Td = [1T P
q
T
T
ˆ −1 Dˆ Pˆ⊥ : PdT = 1T (Dˆ PUˆ⊥ D) U
(10)
0T ](ST −1 S)−1 ST −1
(13)
(14)
which means applying Pd to the whitened data is equivalent P d to the original data. Therefore we can avoid to applying P P Td speci ed in Eq. (14). the data whitening process by using P So a modi ed algorithm for CLDA-based target detection can be described as follows. (1) Identify the desired and undesired targets and construct D, U and S. (2) Calculate the data covariance matrix and its inverse. P d using Eq. (14) and apply it to the original (3) Construct P P Td X. data matrix X by P
4. Constrained linear discriminant analysis-based target classi&er In Section 3 the CLDA-based target detector can detect the p targets of interests but not necessarily distinguish one from another. Instead of using the technique in Ref. [3] for
4
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
classi cation, where the same algorithm should be repeated p times and generate p classi ed gray-scaled images, one for each target, we will classify these targets in a single image and colors are used to separate them. By simply replacing the p × 1 vector 1 in Eq. (10) with a p × 3 color matrix C, the CLDA-based target classi er Pc can be expressed as T
T
ˆ −1 Dˆ P⊥ˆ : PcT = CT (Dˆ PU⊥ˆ D) U
(15)
The ith row of C speci es the color for the ith target display in the classi ed image. For instance, (1; 0; 0)T is used to represent red color, (0; 1; 0)T green and (0; 0; 1)T blue, etc. If T 1 0 0 1 1 0 C=0 1 0 1 0 1 ; 0 0 1 0 1 1 it means that the six targets will be shown in red, green, blue, yellow, magenta and cyan, respectively. The resulting image contains three components for the three primary colors, red (R), green (G) and blue (B). As a consequence, the same algorithm is needed to perform only once for classifying all the targets. The algorithm of CLDA-based target classi cation is described as follows: (1) Identify the desired and undesired targets and construct D and U. (2) Predetermine the colors used for each target display in D and construct the color matrix C. (3) Calculate the data covariance matrix and generate the data whitening operator Pw via eigendecomposition of . (4) Whiten the data by applying Pw to the data matrix X as well as the desired and undesired target matrices, D ˆ Dˆ and U. ˆ and U, resulting in X, (5) The detection of desired targets can be achieved by PcT Xˆ with PcT constructed using Eq. (15). Using a similar derivation in Appendix B, the equivalent P Tc can be written as CLDA-based target classi er P P Tc P
= [C
T
0Tq×3 ](ST −1 S)−1 ST −1
(16)
which can be directly applied to the original data. Here, 0q×3 is a q × 3 matrix with all the elements equal to 0. So a modi ed algorithm for CLDA-based target classi cation can be described as follows. (1) Identify the desired and undesired targets and construct D, U and S. (2) Predetermine the colors used for each target display in D and construct the color matrix C. (3) Calculate the data covariance matrix and its inverse. P c using Eq. (16) and apply it to the original (4) Construct P P Tc X. data matrix X by P
The major advantage of using Eq. (15) or (16) is that the classi cation can be nished with a single color image. If we use the technique in Ref. [3] to classify the p desired targets with p gray-scaled images, we have to adopt a similar classi ers p times. For the purpose Pc of convenient comparison, we apply a corresponding P for such multiple-pass processing to the original image. It can be easily derived that for the ith target classi caP Tc = (0 · · · 0 1 0 · · · 0)(ST −1 S)−1 ST −1 , where tion P i
(0 · · · 0 1 0 · · · 0) is a 1 × (p + q) vector with a single 1 at i
the ith element, 1 6 i 6 p. Comparing it to Eq. (16) where [CT 0Tq×3 ] is a 3 × (p + q) matrix, it can be determined that the ratio of computation times between multiple-pass processing without color and one-pass with color is p : 3, which means when the number of desired targets p is greater than 3, the proposed one-pass approach can always reduce the computational complexity and fasten the processing speed.
5. Real-time implementation On-line image analysis has received considerable attention recently [8,9]. It can provide immediate data processing and assessment, which is particularly useful in military applications where the immediate data assessment is crucial. In this section we develop the real-time implementation of the CLDA algorithm. Since a data whitening process is carried out based on the entire data matrix in the CLDA algorithm, it seems diScult to implement it directly in real-time. So we use modi ed versions to avoid the data whitening process. Assumed that the target matrix S in Eqs. (14) and (16) is known a priori. The task is to derive an algorithm for their real-time implementation. The following matrix inversion lemma in Ref. [12] is used. Lemma 1. Let A and B be two positive-de9nite M × M matrices related by A = B−1 + CD−1 CT ;
(17)
where D is another positive-de9nite N × N matrix; and C is an M × N matrix. A−1 can be calculated by A−1 = B − BC (D + CT BC)−1 CT B:
(18)
The sample correlation matrix R is de ned as R = 1=N Ni=1 xi · xiT , which can be related to and sample mean \ by = R − \ · \T :
(19)
Using the data matrix X, Eq. (19) can be written as N · = X · XT − N · \ · \T . If ˜ denotes N · and R˜ denotes N · R,
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
then V˜ = R˜ − N · \ · \T :
(20)
Suppose that at time t we receive the pixel vector xt . The data matrix Xt at time t is Xt = [x1 ; x2 ; : : : ; xt ] with N vectors. The sample mean, sample correlation and covariance matrices at time t are denoted as \t , Rt and t , respectively. Then Eq. (20) becomes V˜ t = R˜ t − N · \t · \tT :
(21)
At time t + 1 we receive a new pixel vector xt+1 , then the data matrix Xt+1 is Xt+1 = [Xt xt+1 ]. Eq. (21) is updated by T V˜ t+1 = R˜ t+1 − (N + 1) · \t+1 · \t+1
√ √ T = Xt+1 · Xt+1 − (( N + 1)\t+1 ) · (( N + 1)\t+1 )T T = Xt · XtT + xt+1 · xt+1
√ √ − (( N + 1)\t+1 ) · (( N + 1)\t+1 )T T = R˜ t + Pt+1 · Q · Pt+1 ;
(22)
where √ Pt+1 = xt+1 ; ( N + 1)\t+1
Q=
1 0 0 −1
(24)
If we compare Eq. (22) to Eq. (18) and realize Q−1 = Q, −1 V˜ t+1 , R˜ t , Pt+1 , Q in Eq. (22) correspond to A, B, C, D in Eq. (18), respectively, using Eq. (19) we have −1 −1 −1 T ˜ −1 T ˜ −1 V˜ t+1 = R˜ t −R˜ t Pt+1 (Q+Pt+1 Rt Pt+1 )−1 Pt+1 Rt :
(25) The \t+1 in Pt+1 can be updated by \t+1 =
N · \t + xt+1 : N +1 −1 R˜ t
In order to update inversion lemma again:
(26) in Eq. (25), we use the matrix
R˜ t = R˜ t−1 + xt · xtT
(27)
and R˜ t , R˜ t−1 , xt and 1 correspond to A, B, C and D in Eq. (18), respectively, we have −1 −1 −1 −1 −1 R˜ t = R˜ t−1 − R˜ t−1 xt (1 + xtT R˜ t−1 xt )−1 xtT R˜ t−1 :
(1) Using Eq. (28) to update the inverse of the sample cor−1 relation matrix R˜ t at time t. (2) Using Eq. (26) to update the sample mean \t+1 at time t + 1. (3) Using Eq. (25) to update the inverse of the sample co−1 variance matrix V˜ t+1 at time t + 1. (4) Using Eq. (14) to calculate the CLDA-based target deP Td )t+1 or using Eq. (16) for target classi er tector (P T P c )t+1 at time t + 1. (P P Tc )t+1 to xt+1 at time t + 1. P Td )t+1 or (P (5) Apply (P The above algorithm is based on pixel-by-pixel processing. If the data processing is carried out line-by-line from left to right and top to bottom in an image, the line received at time t forms a data matrix Yt = [xt1 xt2 · · · xtM ] where M is the total number of pixels in each line. Assume that the number of pixels received up to time t − 1 is N , then Eq. (25) remains the same except that √ Pt+1 = Yt+1 ; ( N + M )\t+1 (29) and Q=
:
Hence a real-time CLDA-based target detection= classi cation algorithm is described as follows.
(23)
and
(28)
5
IM ×M 0M ×1 01×M −1
;
(30)
where IM ×M is an M × M identity matrix and 0M ×1 is an M × 1 vector with all the elements equal to zero. Eqs. (26) and (28) become N · \t + M i=1 x(t+1)i \t+1 = (31) N +M and −1 −1 −1 −1 −1 R˜ t =R˜ t−1 −R˜ t−1 Yt (IM ×M +YtT R˜ t−1 Yt )−1 YtT R˜ t−1 ;
(32) respectively. The line-by-line real-time implementation algorithm is almost the same as the pixel-by-pixel processing −1 except that the updated V˜ t+1 is applied to the line data matrix Yt+1 received at time t + 1. 6. Experiments The image scene of size 128 × 64 shown in Fig. 1(a) was collected in Maryland in 1995 from the Jight altitude of 10,000 ft with approximately 1:5 m spatial resolution. 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, so the total number of pixel vector N is 128 × 64 = 8192. This scene includes 30 panels arranged in a 10 × 3 matrix. Each element in this matrix
6
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
Fig. 1. (a) A HYDICE image scene which contains 30 panels; (b) Spatial locations of 30 panels provided by ground truth; (c) Spectra of P1; P2; P3; P4; P5; P6; P7; P8; P9; P10.
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 pia ; pib ; pic were made from the same material of size 3 m × 3 m, 2 m × 2 m, 1 m × 1 m, respectively, which could be considered as one class, Pi. The ground truth map in Fig. 1(b) shows the precise locations of the panel center
pixels which were used to generate the target signatures. Fig. 1(c) plots the spectra of these 10 panel classes. We can see that they are very close. In addition, there were two nature objects: grass eld as the background and tree line to the left. Therefore the entire target matrix S contains 12 signatures.
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
7
Fig. 2. The detection results when the 10 panel classes are of interest: 2(a) CLDA detection result; 2(b) CEM detection result; 2(c) FVA detection result; 2(d) OSP detection result.
Fig. 3. The detection results when the panel classes in the odd rows are of interest: 3(a) CLDA detection result; 3(b) CEM detection result; 3(c) FVA detection result; 3(d) OSP detection result.
6.1. Experiment I 6.1.1. Target detection In the rst example about target detection, we assumed that the 10 panels were the desired targets constructing the designed target matrix D = [P1; P2; : : : ; P10], and the undersigned target matrix consisted of two background signatures, i.e., U = [Grass; Tree]. The detection result of the CLDA is shown in Fig. 2(a) and Figs. 2(b–d) are the detection results using CEM, FVA and OSP. We found that the results using FVA and OSP are very close, and both could detect the man-made targets, but also extracted many background signatures. The results from CLDA and CEM-based detectors were signi cantly better because they correctly detected the panels in the rst and second column while rejecting all the background pixels. Small abundance fractions of the panels in the third column were detected because they are of size 1 m × 1 m, smaller than the area 1:5 m × 1:5 m covered by a pixel. In order to demonstrate the capability of grouping di2erent panels, we randomly selected part of panels as desired in the second example. We assumed that the panels in the odd rows belonged to the class of interest, i.e., the desired target matrix D = [P1; P3; P5; P7; P9] and the undesired target matrix U =
[P2; P4; P6; P8; P10; Grass; Tree]. The detection results are shown in Fig. 3(a–d). Now CEM had problem with detecting these ve rows of panels from the others since their spectra are very similar. However, the CLDA still worked very well and provided an accurate detection result. The provided ground truth was used for quantitative study. First we converted the gray-level detected image into a binary image with varied thresholds. Then we compared the results to the ground truth and calculated the correctly detected panel pixels and false alarms. Table 1 lists the number of false alarms required for all the desired panel pixel to be detected in above two examples. It further validates that the CLDA algorithm signi cantly outperforms the FVA and OSP algorithms and slightly outperforms the CEM detector.
Table 1 False alarms required for panel detection No. of panel pixels Example 1 38 Example 2 21
CLDA
CEM
OSP
FVA
84 52
88 58
2268 2218
2446 2369
8
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 background (a) CLDA
(b) FVA
Fig. 4. The classi cation results with color: 4(a) CLDA classi cation result; 4(b) FVA classi cation result.
up because it has diSculty in di2erentiating the panels from each other, which resulted in inappropriate color assignment.
Fig. 5. The classi cation result shown in 3D using CLDA.
6.1.2. Target classi9cation The strength of the CLDA-based classi er was investigated and the comparative study has been conducted in Ref. [2]. Here, we demonstrate how to classify di2erent targets in a single image with color. Fig. 4(a) shows the classi cation results, where the 10 panel classes were represented by di2erent colors and background classes (grass eld and tree line) are in black. Fig. 5 is a 3-D classi ed image with the third dimension speci ed by target abundance fraction. The magnitudes of the left-most panels are the greatest because of their largest size, while the sizes of the right-most panels are very small so that their magnitudes are the lowest. The FVA was also applied for the classi cation, whose result is shown in Fig. 4(b). The colors were totally messed
6.1.3. Real-time implementation Figs. 6(a–d) depicts a sequence of line-by-line real-time processed images that discriminate these 10 panel classes in di2erent stages with the CLDA algorithm. It rst classi ed P1 in red and all the background objects in black, then it classi ed P2 and displayed it in yellow, after that P3 was classi ed in green, and so on. If we compare the nal on-line classi cation in Fig. 6(d) with the o2-line classi cation in Fig. 4(a), we nd that there are some tiny di2erences between them. In particular, the di2erence in P1 and P2 classi cation is perceivable, where the panel sizes of the classi ed P1 and P2 are smaller in Fig. 6(d). The reason is that the on-line classi er (here it is about the −1 calculation) is generated and updated based on the received pixels instead of all the pixels in the entire scene as in the o2-line case. However, the discrepancies brought about by the on-line processing are tolerable. The entire real-time simulation was nished very quickly. It took a Sun Ultra10 440 MHz workstation only about 50 s. 6.2. Experiment II In the second experiment with HYDICE data, the scene ◦ was rotated by 90 counterclockwise as shown in Fig. 7. The 10 panel classes now were ordered from left to right. P1 was in the rst column, P2 the second and so on. 1 m×1 m panels are in the rst row, 2 m × 2 m the second and 3 m × 3 m
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
9
Fig. 6. Real-time classi cation I.
The rst row of panels were classi ed rst in pre-assigned di2erent colors. Since they are of size 1 m ×1 m and smaller than a pixel area, the color shades are light. p1b , p2b , p3b , p4b and p5b in the second row are among the next group of classi ed panels, which were displayed in red, yellow, green, magenta and cyan as pre-assigned, respectively. After all the 10 panels in the second row were classi ed, those in the third row were to be classi ed. p1a , p2a , p3a , p4a and p5a were displayed in the same color as p1b , p2b , p3b , p4b and p5b as the same classes respectively. If we compared the nal image in Fig. 8(j) with Fig. 6(h), we found that these two results were very close. This experiment demonstrates that the performance of the proposed real-time algorithm is irrelevant to the data processing order. 7. Conclusion
Fig. 7. (a) The image scene used in Experiment II after rotating ◦ the image in Fig. 1(a) 90 counterclockwise; (b) Spatial locations of 30 panels.
the third. Only real-time implementation may be a2ected since the data are processed from left to right and top to bottom. The classi ed image sequence is shown in Fig. 8.
This paper presents CLDA-based target detection and classi cation algorithms as well as their real-time implementation for hyperspectral imagery. The CLDA-based target detector is derived by treating desired and undesired classes as two groups. As a result, no new sample means are needed to recalculate. The CLDA-based target classi er introduces the color for display such that great convenience is o2ered by visual discrimination of di2erent targets in a single image. All the CLDA-based detection and classi cation
10
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
Fig. 8. Real-time classi cation II.
algorithms have great potential to be implemented in real-time, which is very important in defense and intelligence related applications. Both pixel-by-pixel and line-by-line real-time processing algorithms are developed. The experiment results with HYDICE data demonstrated its excellent performance in hyperspectral image detection, classi cation and real-time implementation.
Acknowledgements The authors would like to thank Professor Chein-I Chang at the University of Maryland Baltimore County for providing the image data used in the experiments. This work was performed while the second author held a National Research Council Research Associateship Award at Edgewood Chemical and Biological Center, US Army.
Appendix A The following theorem for orthogonal decomposition of vectors is being used [13]: Theorem. Given ; a vector subspace of n-dimensional Euclidean space Rn ; every n × 1 vector y can be expressed uniquely in the form y = u + v; where u ∈ and v ∈ ⊥ . Here, the L × 1 vector w can be represented as w = wpfs + wofs = wpus + wous + wofs ;
(A.1)
where wpfs and wofs are the projections of w onto the subspace spanned by the signature vectors (feature subspace) and its orthogonal complement, respectively. wpfs can be further decomposed as the sum of its projections onto the undesired subspace wpus and its orthogonal complement wous .
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
Since the undesired signatures U are orthogonal to wous and wofs , the constraint in Eq. (9) becomes wT U = (wpus )T U = 0T :
(A.2)
11
If the partition is S = [S1 S2 ], then T S1 S1 ST1 S2 : ST S = ST2 S1 ST2 S2
(B.2)
Assumed that the signatures in U are linear independent. wpus in Eq. (34) must be a zero vector with all the elements equal to 0 such that Eq. (A.2) can hold. The next task is to determine wous and wofs in Eq. (A.1). The projector of the orthogonal complement of the undesired signature subspace can be constructed by PU⊥ =I−UU# , where U# = (UT U)−1 UT is the pseudo-inverse of U. Therefore, wous is proportional to PU⊥ and can be represented as ⊥ wous = PU a, where a is a vector to be determined. Since the desired signatures D are orthogonal to wofs , the constraint in Eq. (8) becomes
Comparing Eq. (B.2) with Eq. (B.1), A, B, D, E, F in Eq. (B.1) are
⊥ T ⊥ wT D = (wous )T D = (PU a) D = aT PU D = 1T :
F = (ST1 S1 )−1 ST1 S2 :
(A.3)
A = ST1 S1
(B.3)
B = ST1 S2
(B.4)
D = ST2 S2
(B.5)
E = ST2 S2 − ST2 S1 (ST1 S1 )−1 ST1 S2 = ST2 (I − S1 (ST1 S1 )−1 ST1 )S2 = ST2 PS⊥1 S2
(B.6) (B.7)
Substituting Eqs. (B.3) – (B.7) into Eq. (B.2), (ST S)−1 =
(ST1 S1 )−1 + (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ST2 S1 (ST1 S1 )−1
−(ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1
−(ST2 PS⊥1 S2 )−1 ST2 S1 (ST1 S1 )−1
(ST2 PS⊥1 S2 )−1
(B.8)
and T
(S S)
−1 T
T
S = (S S) =
ST1
ST2
(ST1 S1 )−1 ST1 + (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ST2 S1 (ST1 S1 )−1 ST1 − (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ST2
Sx1 + Sx2 − Sx3 (ST2 PS⊥1 S2 )−1 ST2 PS⊥1
;
(B.9)
⊥ Therefore, a can be estimated by aT = 1T (PU D)# = 1T (DT PU⊥ D)−1 DT PU⊥ . As for wofs , it is orthogonal to wous . In order to make w2 minimal, wofs should also be a zero vector. Hence,
w=w
ous
=
⊥ PU a
=
−(ST2 PS⊥1 S2 )−1 ST2 S1 (ST1 S1 )−1 ST1 + (ST2 PS⊥1 S2 )−1 ST2
=
−1
PU⊥ D(DT PU⊥ D)−1 1:
(A.4)
where Sx1 =(ST1 S1 )−1 ST1 , Sx2 =(ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ×ST2 S1 (ST1 S1 )−1 ST1 and Sx3 =(ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ST2 . Using the matrix inversion lemma and matrix operations: (ST1 PS⊥2 S1 )−1 ST1 PS⊥2 =(ST1 S1 − ST1 S2 (ST2 S2 )−1 ST2 S1 )−1 ST1 PS⊥2
Appendix B The following partitioned matrix inversion formula is used [13]: −1 −1 A B A + FE−1 FT −FE−1 ; (B.1) = BT D E−1 −E−1 FT where E = D − BT A−1 B, F = A−1 B.
=[(ST1 S1 )−1 + (ST1 S1 )−1 ST1 S2 (ST2 S2 − ST2 S1 ×(ST1 S1 )−1 ST1 S2 )−1 ST2 S1 (ST1 S1 )−1 ]ST1 PS⊥2 =[(ST1 S1 )−1 + (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ×ST2 S1 (ST1 S1 )−1 ]ST1 PS⊥2
12
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
=[(ST1 S1 )−1 ST1 + (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1
column vector with all the elements equal to 0, i.e., 0 = (0; 0; : : : ; 0)T .
×ST2 S1 (ST1 S1 )−1 ST1 ](I − S2 (ST2 S2 )−1 ST2 ) =[Sx1 + Sx2 ](I −
q
S2 (ST2 S2 )−1 ST2 )
=Sx1 + Sx2 − Sx1 S2 (ST2 S2 )−1 ST2 − Sx2 S2 (ST2 S2 )−1 ST2 =Sx1 + Sx2 − (ST1 S1 )−1 ST1 S2 [(ST2 S2 )−1 +(ST2 PS⊥1 S2 )−1 ST2 S1 (ST1 S1 )−1 ST1 S2 (ST2 S2 )−1 ]ST2 =Sx1 + Sx2 − (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 [ST2 PS⊥1 S2 (ST2 S2 )−1 + ST2 S1 (ST1 S1 )−1 ST1 S2 (ST2 S2 )−1 ]ST2 =Sx1 + Sx2 − (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 [ST2 PS⊥1 S2 (ST2 S2 )−1 + ST2 PS1 S2 (ST2 S2 )−1 ]ST2 =Sx1 + Sx2 − (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 [ST2 S2 (ST2 S2 )−1 ]ST2 =Sx1 + Sx2 − (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ST2 =Sx1 + Sx2 − Sx3 which means the rst element in Eq. (B.9) equals (ST1 PS⊥2 S1 )−1 ST1 PS⊥2 , i.e., T ⊥ −1 T ⊥ (S1 PS2 S1 ) S1 PS2 T −1 T (S S) S = : (B.10) (ST2 PS⊥1 S2 )−1 ST2 PS⊥1 If D = S1 and U = S2 are L × p and L × q matrices, respectively, T ⊥ −1 T ⊥ (D PU D) D PU T T T −1 T T T [1 0 ](S S) S = [1 0 ] · (UT PD⊥ U)−1 UT PD⊥ = 1T (DT PU⊥ D)−1 DT PU⊥ ;
(B.11)
where 1 is a p × 1 column vector with all the elements equal to 1, i.e., 1 = (1; 1; : : : ; 1)T and 0 is a q × 1 p
References [1] A.F.H. Goetz, G. Vane, J.E. Solomon, B.N. Rock, Imaging spectrometry for earth remote sensing, Science 228 (4704) (1985) 1147–1153. [2] R.O. Duda, P.E. Hart, Pattern Classi cation Scene Anal, Wiley, New York, 1973. [3] Q. Du, C.-I Chang, Linear constrained distance-based discriminant analysis for hyperspectral image classi cation, Pattern Recognition 34 (2) (2001) 361–373. [4] J.C. Harsanyi, C.-I Chang, Hyperspectral image classi cation and dimensionality reduction: an orthogonal subspace projection IEEE Trans. Geosci. Remote Sensing 32 (4) (1994) 779–785. [5] W.H. Farrand, J.C. Harsanyi, Mapping the distribution of mine tailing in the coeur d’Alene river valley, Idaho through the use of constrained energy minimization technique, Remote Sensing Environ. 59 (1997) 64–76. [6] J. Bowles, P. Palmadesso, J. Antoniades, M. Baumback, Use of lter vectors in hyperspectral data analysis, Proc. SPIE 2553 (1995) 148–152. [7] C.-I Chang, Hyperspectral Imaging: Subpixel Detection and Classi cation,KluwerAcademicPublishers, Dordrecht, 2002. [8] C.M. Stellman, G.G. Hazel, F. Bucholtz, J.V. Michalowicz, Real-time hyperspectral detection and cuing, Opt. Eng. (2000) 1928–1935. [9] C.-I Chang, H. Ren, S.S. Chiang, Real-time processing algorithms for target detection and classi cation in hyperspectral imagery, IEEE Trans. Geosci. Remote Sensing 39 (4) (2001) 760–768. [10] H. Ren, C.-I Chang, An unsupervised orthogonal subspace projection to target detection and classi cation in unknown environment, Spectroradiometric Symposium, San Diego, November 2–7, 1997. [11] D.C. Heinz, C.-I. Chang, Fully constrained least squares linear mixture analysis for material quanti cation in hyperspectral imagery, IEEE Trans. Geosci. Remote Sensing 39 (3) (2001) 529–545. [12] S. Haykin, Adaptive Filter Theory, 3rd Edition, Prentice-Hall, Englewood Cli2s, NJ, 1996. [13] R.J. Muirhead, Aspects of Multivariate Statistical Theory, Wiley, New York, 1982.
About the Author—QIAN DU received the Ph.D. degree in Electrical Engineering from University of Maryland Baltimore County in 2000. She is currently an assistant professor at the Department of Electrical Engineering and Computer Science, Texas A & M University-Kingsville. Her research interests include signal and image processing, pattern recognition, neural networks and remote sensing. She is a member of IEEE, SPIE and the honor society of Phi Kappa Phi. About the Author—HSUAN REN received the B.S. degree in Electrical Engineering from the National Taiwan University, Taipei, Taiwan, R.O.C., in 1994 and the M.S. and Ph.D. degrees from University of Maryland Baltimore County in 1998 and 2000, respectively, all in Electrical Engineering. He is currently on National Research Council (NRC) Associateship and supported by U.S. Army Edgewood Chemical Biological Center as a postdoctoral fellow. He has two patents pending on hyperspectral target detection and image classi cation. Dr. Ren is also a part-time instructor at the Department of Computer Science and Electrical Engineering, University of Maryland Baltimore County. His research interests include data compression, signal and image processing and pattern recognition. He is a member of IEEE and SPIE. He is also a member of the honor society of Phi Kappa Phi.
www.elsevier.com/locate/patcog
Real-time constrained linear discriminant analysis to target detection and classi cation in hyperspectral imagery Qian Dua; ∗ , Hsuan Renb a Department
of Electrical Engineering and Computer Science, MSC 192, Texas A& M University-Kingsville, TX 78363, USA Biological and Chemical Center, US Army Aberdeen Proving Ground, MD 21010, USA
b Edgewood
Received 11 April 2001; received in revised form 2 November 2001; accepted 6 March 2002
Abstract In this paper, we present a constrained linear discriminant analysis (CLDA) approach to hyperspectral image detection and classi cation as well as its real-time implementation. The basic idea of CLDA is to design an optimal transformation matrix which can maximize the ratio of inter-class distance to intra-class distance while imposing the constraint that di2erent class centers after transformation are along di2erent directions such that di2erent classes can be better separated. The solution turns out to be a constrained version of orthogonal subspace projection (OSP) implemented with a data whitening process. The CLDA approach can be applied to solve both detection and classi cation problems. In particular, by introducing color for display the classi cation is achieved with a single classi ed image where a pre-assigned color is used to display a speci ed class. The real-time implementation is also developed to meet the requirement of on-line image analysis when the immediate data assessment is critical. The experiments using HYDICE data demonstrate the strength of CLDA approach in discriminating the targets with subtle spectral di2erence. ? 2002 Pattern Recognition Society. Published by Elsevier Science Ltd. All rights reserved. Keywords: Detection; Classi cation; Real-time processing; Constrained linear discriminant analysis (CLDA); Hyperspectral imagery
1. Introduction Hyperspectral imaging is a new technique in remote sensing, which generates hundreds of images with di2erent spectral channels for the same area on the earth. Because of very high spectral resolution, it is most likely to better discriminate di2erent materials using hyperspectral images than using multispectral images, which contains only tens of spectral channels. For instance, the airborne visible=infrared imaging spectrometer (AVIRIS) developed by NASA’s Jet Propulsion Laboratory has 224 spectral channels with 10 nm spectral resolution and 20 m spatial resolution while the hyperspectral digital image collection experiment (HYDICE) sensor developed by Naval Research Laboratory has 210 spectral channels with the same 10 nm spectral resolution and 1–4 m spatial resolution. AVIRIS is mainly used for ∗ Corresponding author. Tel.: +1-361-5932625; fax: +1-3615932110. E-mail address: [email protected] (Q. Du).
ecology, geology, global climate monitoring, land mapping, etc., while HYDICE can be used for surveillance and reconnaissance purposes because of its ner spatial resolution. One of the major challenging problems about hyperspectral image analysis is how to deal with the massive data volume while taking full advantage of plentiful spectral information [1]. Moreover, the size of targets generally is smaller than the spatial resolution or the area covered by a pixel, so image analysis in remote sensing actually copes with mixed pixel processing instead of pure pixel processing in standard digital images. In general, new algorithm development for hyperspectral imagery is necessary. Fisher’s linear discriminant analysis (LDA) is a standard technique in pattern recognition. It projects the original high dimensional data onto a low dimensional space, where all the classes are well separated by maximizing the Raleigh quotient, i.e., the ratio of between-class scatter matrix determinant to within-class scatter matrix determinant [2]. The application of Fisher’s LDA to hyperspectral images has been investigated for hyperspectral image classi cation
0031-3203/02/$22.00 ? 2002 Pattern Recognition Society. Published by Elsevier Science Ltd. All rights reserved. PII: S 0 0 3 1 - 3 2 0 3 ( 0 2 ) 0 0 0 6 5 - 1
2
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
by Du and Chang in Ref. [3]. The derived algorithm in Ref. [3] is referred to as constrained linear discriminant analysis (CLDA) algorithm in this paper. In CLDA the original high dimensional data is projected onto a low dimensional space as done by Fisher’s LDA. But di2erent classes are forced to be along di2erent directions in this low dimensional space. Thus all classes were expected to be better separated and the classi cation is achieved simultaneously with the CLDA transform. The transformation matrix in CLDA maximizes the ratio of inter-class distance to intra-class distance while satisfying the constraint that the means of di2erent classes are aligned with di2erent directions, which can be constructed by using an orthogonal subspace projection (OSP) method coupled with a data whitening process. The experimental results in Ref. [3] demonstrated that the CLDA algorithm could provide more accurate classi cation results than other popular methods in hyperspectral image processing, such as OSP classi er [4] and constrained energy minimization (CEM) operator [5]. Recently, a new approach, referred to as lter vector algorithm (FVA) was proposed by J. Bowles et al. and its advantage in hyperspectral image classi cation was validated in Ref. [6]. In our experiments we found out that the CLDA algorithm outperforms the FVA. Interestingly, the FVA can be viewed as a special case of the CLDA algorithm [7]. Although the CLDA algorithm showed its strength in hyperspectral image processing, further improvement can be made in two aspects. One is that it was only used for classi cation in Ref. [3], not detection. Another is that the p images are needed to produce by repeating the same procedure p times, one for each target classi cation. Therefore, in this paper we will extend and modify the CLDA algorithm to accomplish target detection and classify all the targets in a single image. Since recently the spectral and spatial resolutions of hyperspectral images are signi cantly improved, their applications have been expanded to battle eld and national defense area where the real-time data processing becomes inevitable so as to process data on-line and provide the information for immediate response [8,9]. In order to meet this demand, a real-time version of the CLDA algorithm is also developed. The remainder of this paper is organized as follows. Section 2 brieJy reviews the CLDA algorithm. Section 3 describes the target detector based on the CLDA algorithm. Section 4 focuses on the CLDA target classi er. In Section 5, the real-time implementation of the CLDA-based target detection and classi cation algorithm is presented. In Section 6, HYDICE data experiments are used to demonstrate the performance of the CLDA detector and classi er as well as their real-time implementation. Finally, Section 7 concludes this paper with several remarks. 2. Constrained linear discriminant analysis Assume that there are c classes and the kth class contains Nk patterns. Let N = N1 + · · · + Nc be the total number of
training patterns. The jth pattern in the kth class, denoted k k k T by xjk = [x1j ; x2j ; : : : ; xLj ] is an L-dimensional vector. Let c Nk k \ = 1=N k=1 j=1 xj be the global mean of training sam k k ples and \k = 1=Nk Nj=1 xj the mean of the kth class. We de ned J (F) to be the ratio of the inter-class distance to the intra-class distance after a linear transformation F, which is given by c 2 2=c(c − 1) c−1 i=1 j=i+1 F(\i ) − F(\j ) J (F) = (1) c Nk k 1=cN k=1 [ j=1 F(xj ) − F(\k )2 ] and F(x) = (WL×c )T x = [w1 ; w2 ; : : : ; wc ]T x;
(2)
where wi is an L × 1 vector. The optimal linear transformation F ∗ is the one that maximizes J (F) subject to tk = F(\k ) for all k; where
(3)
T
tk = 0 · · · 0 1 0 · · · 0 k
is a c × 1 unit column vector with one in the kth component and zeros elsewhere. F ∗ can be determined by wi∗ = \ˆ Ti P⊥ ˆi U
(4)
where ˆ ˆ T ˆ −1 U ˆ Ti P⊥ ˆ i = I − Ui (Ui Ui ) U
(5)
with Uˆ i = [ uˆ 1 · · · uˆ j · · · uˆ c ]j=i and I the identity matrix. The “∧” speci es the whitened data, i.e., xˆ = PwT x, where Pw is the data whitening operator and can be calculated by PwT = V2−1=2 V1T :
(6)
Here, V1 and V2 are eigenvector and eigenvalue matrices of the data sample covariance matrix = 1=N Ni=1 T (xi − \)(xi − \) , respectively. They can be determined by the eigendecomposition = V1 V2 V1T , where V1 = [v1 ; v2 ; : : : ; vL ] with vi being the ith eigenvector and V2 is a diagonal matrix with the ith diagonal item being the ith eigenvalue corresponding to the ith eigenvector, i.e., V2 = diag (1 ; 2 ; : : : ; L ).
3. Constrained linear discriminant analysis-based target detector Assume that an image scene contains p desired targets, say, man-made targets and q undesired targets, say, background objects. For derivation purpose, we treat the detection problem as a two-class classi cation problem: the p
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
desired targets as class !D with desired class mean matrix D=\1 \2 · · · \p and the q undesired targets as another class !U with undesired class mean matrix U = \p+1 \p+2 · · · \p+q . We have p + q = c, where c is the total number of di2erent target classes. The patterns in class !D are {xD } = {xk }pk=1 and the number of patterns ND in !D equals ND = N1 + N2 + · · · + Np . The patterns in class !U are {xU } = {xk }p+q k=p+1 and the number of patterns NU in !U equals NU = Np+1 + Np+2 + · · · + Np+q . The sample means \D and can be calculated by \D = \U of these two classes 1=ND pi=1 Ni \i and \U = 1=NU p+q i=p+1 Ni \i , respectively. For this two-class classi cation problem, Eq. (1) becomes J (F) = F(\D )−F(\U )2 ;
U D 1=2N Nj=1 F(xjD )−F(\D )2 + Nj=1 F(xjU )−F(\U )2 (7) where F(x) = wT x and w is an L × 1 transformation vector. Our goal is to minimize Eq. (7) while imposing the constraints w T D = 1T
(8)
and
3
The algorithm of CLDA-based target detection is presented as follows. (1) Identify the desired and undesired targets and construct D and U. (2) Calculate the data covariance matrix and generate the data whitening operator Pw via eigendecomposition of . (3) Whiten the data by applying Pw to the data matrix X as well as the desired and undesired target matrices, D ˆ Dˆ and U. ˆ and U, resulting in X, (4) The detection of desired targets can be achieved by PdT Xˆ with PdT constructed using Eq. (10). It should be noted that if D and U are determined by an unsupervised method in step (1) using the techniques in Ref. [10] or [11], this algorithm can be extended to an unsupervised version. In Appendix B, we prove that 1T (DT PU⊥ D)−1 DT PU⊥ = [1T
0T ](ST S)−1 ST ;
(11)
where S = [D U] is the entire target matrix. If the data is whitened, Eq. (11) becomes T ˆ −1 Dˆ T P⊥ˆ = [1T PdT = 1T (Dˆ PU⊥ˆ D) U
T ˆ −1 Sˆ T : 0T ](Sˆ S)
(12)
Since Sˆ = PwT S and Xˆ = PwT X,
wT U = 0T
T ˆ −1 Sˆ T Xˆ 0T ](Sˆ S)
(9)
PdT Xˆ = [1T
where 1 is a p ×1 column vector with all the elements equal to 1, i.e., 1=(1; 1; : : : ; 1)T and 0 is a q×1 column vector with
= [1T
0T ](ST Pw PwT S)−1 ST Pw PwT X
= [1T
P Td X; 0T ](ST −1 S)−1 ST −1 X = P
p
T
all the elements equal to 0, i.e., 0 = (0; 0; : : : ; 0) . If Eqs. (8)
where
and (9) are satis ed, it is easy to show that wT \D = 1 and wT \U = 0, which means the constraint in Eq. (3) is relaxed. The advantage of using Eqs. (8) and (9) as constraints is that \D and \U are unnecessary to be calculated. We notice that the numerator in Eq. (7) is a constant when the constraints are satis ed. Thus maximizing J (F) is equivalent to minimizing the denominator, which is equal to (1=2N )wT SW w with SW denoting the with-class scatter matrix de ned in Ref. [2]. Since SW + SB = ST = N · , minimizing (1=2N )wT SW w is equivalent to minimizing wT w, where SB and ST are the between-class scatter matrix and the total scatter matrix de ned in Ref. [2], respectively. After the data is whitened by Pw in Eq. (6), becomes the identity matrix. Hence the constrained detection problem can be described as: to nd a transformation vector w such that w2 is minimum with the constraints that wT D = 1 and wT U = 0. In Appendix A we derive the solution w∗ of this constrained detection problem, referred to as CLDA-based target detector Pd
P Td = [1T P
q
T
T
ˆ −1 Dˆ Pˆ⊥ : PdT = 1T (Dˆ PUˆ⊥ D) U
(10)
0T ](ST −1 S)−1 ST −1
(13)
(14)
which means applying Pd to the whitened data is equivalent P d to the original data. Therefore we can avoid to applying P P Td speci ed in Eq. (14). the data whitening process by using P So a modi ed algorithm for CLDA-based target detection can be described as follows. (1) Identify the desired and undesired targets and construct D, U and S. (2) Calculate the data covariance matrix and its inverse. P d using Eq. (14) and apply it to the original (3) Construct P P Td X. data matrix X by P
4. Constrained linear discriminant analysis-based target classi&er In Section 3 the CLDA-based target detector can detect the p targets of interests but not necessarily distinguish one from another. Instead of using the technique in Ref. [3] for
4
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
classi cation, where the same algorithm should be repeated p times and generate p classi ed gray-scaled images, one for each target, we will classify these targets in a single image and colors are used to separate them. By simply replacing the p × 1 vector 1 in Eq. (10) with a p × 3 color matrix C, the CLDA-based target classi er Pc can be expressed as T
T
ˆ −1 Dˆ P⊥ˆ : PcT = CT (Dˆ PU⊥ˆ D) U
(15)
The ith row of C speci es the color for the ith target display in the classi ed image. For instance, (1; 0; 0)T is used to represent red color, (0; 1; 0)T green and (0; 0; 1)T blue, etc. If T 1 0 0 1 1 0 C=0 1 0 1 0 1 ; 0 0 1 0 1 1 it means that the six targets will be shown in red, green, blue, yellow, magenta and cyan, respectively. The resulting image contains three components for the three primary colors, red (R), green (G) and blue (B). As a consequence, the same algorithm is needed to perform only once for classifying all the targets. The algorithm of CLDA-based target classi cation is described as follows: (1) Identify the desired and undesired targets and construct D and U. (2) Predetermine the colors used for each target display in D and construct the color matrix C. (3) Calculate the data covariance matrix and generate the data whitening operator Pw via eigendecomposition of . (4) Whiten the data by applying Pw to the data matrix X as well as the desired and undesired target matrices, D ˆ Dˆ and U. ˆ and U, resulting in X, (5) The detection of desired targets can be achieved by PcT Xˆ with PcT constructed using Eq. (15). Using a similar derivation in Appendix B, the equivalent P Tc can be written as CLDA-based target classi er P P Tc P
= [C
T
0Tq×3 ](ST −1 S)−1 ST −1
(16)
which can be directly applied to the original data. Here, 0q×3 is a q × 3 matrix with all the elements equal to 0. So a modi ed algorithm for CLDA-based target classi cation can be described as follows. (1) Identify the desired and undesired targets and construct D, U and S. (2) Predetermine the colors used for each target display in D and construct the color matrix C. (3) Calculate the data covariance matrix and its inverse. P c using Eq. (16) and apply it to the original (4) Construct P P Tc X. data matrix X by P
The major advantage of using Eq. (15) or (16) is that the classi cation can be nished with a single color image. If we use the technique in Ref. [3] to classify the p desired targets with p gray-scaled images, we have to adopt a similar classi ers p times. For the purpose Pc of convenient comparison, we apply a corresponding P for such multiple-pass processing to the original image. It can be easily derived that for the ith target classi caP Tc = (0 · · · 0 1 0 · · · 0)(ST −1 S)−1 ST −1 , where tion P i
(0 · · · 0 1 0 · · · 0) is a 1 × (p + q) vector with a single 1 at i
the ith element, 1 6 i 6 p. Comparing it to Eq. (16) where [CT 0Tq×3 ] is a 3 × (p + q) matrix, it can be determined that the ratio of computation times between multiple-pass processing without color and one-pass with color is p : 3, which means when the number of desired targets p is greater than 3, the proposed one-pass approach can always reduce the computational complexity and fasten the processing speed.
5. Real-time implementation On-line image analysis has received considerable attention recently [8,9]. It can provide immediate data processing and assessment, which is particularly useful in military applications where the immediate data assessment is crucial. In this section we develop the real-time implementation of the CLDA algorithm. Since a data whitening process is carried out based on the entire data matrix in the CLDA algorithm, it seems diScult to implement it directly in real-time. So we use modi ed versions to avoid the data whitening process. Assumed that the target matrix S in Eqs. (14) and (16) is known a priori. The task is to derive an algorithm for their real-time implementation. The following matrix inversion lemma in Ref. [12] is used. Lemma 1. Let A and B be two positive-de9nite M × M matrices related by A = B−1 + CD−1 CT ;
(17)
where D is another positive-de9nite N × N matrix; and C is an M × N matrix. A−1 can be calculated by A−1 = B − BC (D + CT BC)−1 CT B:
(18)
The sample correlation matrix R is de ned as R = 1=N Ni=1 xi · xiT , which can be related to and sample mean \ by = R − \ · \T :
(19)
Using the data matrix X, Eq. (19) can be written as N · = X · XT − N · \ · \T . If ˜ denotes N · and R˜ denotes N · R,
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
then V˜ = R˜ − N · \ · \T :
(20)
Suppose that at time t we receive the pixel vector xt . The data matrix Xt at time t is Xt = [x1 ; x2 ; : : : ; xt ] with N vectors. The sample mean, sample correlation and covariance matrices at time t are denoted as \t , Rt and t , respectively. Then Eq. (20) becomes V˜ t = R˜ t − N · \t · \tT :
(21)
At time t + 1 we receive a new pixel vector xt+1 , then the data matrix Xt+1 is Xt+1 = [Xt xt+1 ]. Eq. (21) is updated by T V˜ t+1 = R˜ t+1 − (N + 1) · \t+1 · \t+1
√ √ T = Xt+1 · Xt+1 − (( N + 1)\t+1 ) · (( N + 1)\t+1 )T T = Xt · XtT + xt+1 · xt+1
√ √ − (( N + 1)\t+1 ) · (( N + 1)\t+1 )T T = R˜ t + Pt+1 · Q · Pt+1 ;
(22)
where √ Pt+1 = xt+1 ; ( N + 1)\t+1
Q=
1 0 0 −1
(24)
If we compare Eq. (22) to Eq. (18) and realize Q−1 = Q, −1 V˜ t+1 , R˜ t , Pt+1 , Q in Eq. (22) correspond to A, B, C, D in Eq. (18), respectively, using Eq. (19) we have −1 −1 −1 T ˜ −1 T ˜ −1 V˜ t+1 = R˜ t −R˜ t Pt+1 (Q+Pt+1 Rt Pt+1 )−1 Pt+1 Rt :
(25) The \t+1 in Pt+1 can be updated by \t+1 =
N · \t + xt+1 : N +1 −1 R˜ t
In order to update inversion lemma again:
(26) in Eq. (25), we use the matrix
R˜ t = R˜ t−1 + xt · xtT
(27)
and R˜ t , R˜ t−1 , xt and 1 correspond to A, B, C and D in Eq. (18), respectively, we have −1 −1 −1 −1 −1 R˜ t = R˜ t−1 − R˜ t−1 xt (1 + xtT R˜ t−1 xt )−1 xtT R˜ t−1 :
(1) Using Eq. (28) to update the inverse of the sample cor−1 relation matrix R˜ t at time t. (2) Using Eq. (26) to update the sample mean \t+1 at time t + 1. (3) Using Eq. (25) to update the inverse of the sample co−1 variance matrix V˜ t+1 at time t + 1. (4) Using Eq. (14) to calculate the CLDA-based target deP Td )t+1 or using Eq. (16) for target classi er tector (P T P c )t+1 at time t + 1. (P P Tc )t+1 to xt+1 at time t + 1. P Td )t+1 or (P (5) Apply (P The above algorithm is based on pixel-by-pixel processing. If the data processing is carried out line-by-line from left to right and top to bottom in an image, the line received at time t forms a data matrix Yt = [xt1 xt2 · · · xtM ] where M is the total number of pixels in each line. Assume that the number of pixels received up to time t − 1 is N , then Eq. (25) remains the same except that √ Pt+1 = Yt+1 ; ( N + M )\t+1 (29) and Q=
:
Hence a real-time CLDA-based target detection= classi cation algorithm is described as follows.
(23)
and
(28)
5
IM ×M 0M ×1 01×M −1
;
(30)
where IM ×M is an M × M identity matrix and 0M ×1 is an M × 1 vector with all the elements equal to zero. Eqs. (26) and (28) become N · \t + M i=1 x(t+1)i \t+1 = (31) N +M and −1 −1 −1 −1 −1 R˜ t =R˜ t−1 −R˜ t−1 Yt (IM ×M +YtT R˜ t−1 Yt )−1 YtT R˜ t−1 ;
(32) respectively. The line-by-line real-time implementation algorithm is almost the same as the pixel-by-pixel processing −1 except that the updated V˜ t+1 is applied to the line data matrix Yt+1 received at time t + 1. 6. Experiments The image scene of size 128 × 64 shown in Fig. 1(a) was collected in Maryland in 1995 from the Jight altitude of 10,000 ft with approximately 1:5 m spatial resolution. 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, so the total number of pixel vector N is 128 × 64 = 8192. This scene includes 30 panels arranged in a 10 × 3 matrix. Each element in this matrix
6
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
Fig. 1. (a) A HYDICE image scene which contains 30 panels; (b) Spatial locations of 30 panels provided by ground truth; (c) Spectra of P1; P2; P3; P4; P5; P6; P7; P8; P9; P10.
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 pia ; pib ; pic were made from the same material of size 3 m × 3 m, 2 m × 2 m, 1 m × 1 m, respectively, which could be considered as one class, Pi. The ground truth map in Fig. 1(b) shows the precise locations of the panel center
pixels which were used to generate the target signatures. Fig. 1(c) plots the spectra of these 10 panel classes. We can see that they are very close. In addition, there were two nature objects: grass eld as the background and tree line to the left. Therefore the entire target matrix S contains 12 signatures.
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
7
Fig. 2. The detection results when the 10 panel classes are of interest: 2(a) CLDA detection result; 2(b) CEM detection result; 2(c) FVA detection result; 2(d) OSP detection result.
Fig. 3. The detection results when the panel classes in the odd rows are of interest: 3(a) CLDA detection result; 3(b) CEM detection result; 3(c) FVA detection result; 3(d) OSP detection result.
6.1. Experiment I 6.1.1. Target detection In the rst example about target detection, we assumed that the 10 panels were the desired targets constructing the designed target matrix D = [P1; P2; : : : ; P10], and the undersigned target matrix consisted of two background signatures, i.e., U = [Grass; Tree]. The detection result of the CLDA is shown in Fig. 2(a) and Figs. 2(b–d) are the detection results using CEM, FVA and OSP. We found that the results using FVA and OSP are very close, and both could detect the man-made targets, but also extracted many background signatures. The results from CLDA and CEM-based detectors were signi cantly better because they correctly detected the panels in the rst and second column while rejecting all the background pixels. Small abundance fractions of the panels in the third column were detected because they are of size 1 m × 1 m, smaller than the area 1:5 m × 1:5 m covered by a pixel. In order to demonstrate the capability of grouping di2erent panels, we randomly selected part of panels as desired in the second example. We assumed that the panels in the odd rows belonged to the class of interest, i.e., the desired target matrix D = [P1; P3; P5; P7; P9] and the undesired target matrix U =
[P2; P4; P6; P8; P10; Grass; Tree]. The detection results are shown in Fig. 3(a–d). Now CEM had problem with detecting these ve rows of panels from the others since their spectra are very similar. However, the CLDA still worked very well and provided an accurate detection result. The provided ground truth was used for quantitative study. First we converted the gray-level detected image into a binary image with varied thresholds. Then we compared the results to the ground truth and calculated the correctly detected panel pixels and false alarms. Table 1 lists the number of false alarms required for all the desired panel pixel to be detected in above two examples. It further validates that the CLDA algorithm signi cantly outperforms the FVA and OSP algorithms and slightly outperforms the CEM detector.
Table 1 False alarms required for panel detection No. of panel pixels Example 1 38 Example 2 21
CLDA
CEM
OSP
FVA
84 52
88 58
2268 2218
2446 2369
8
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 background (a) CLDA
(b) FVA
Fig. 4. The classi cation results with color: 4(a) CLDA classi cation result; 4(b) FVA classi cation result.
up because it has diSculty in di2erentiating the panels from each other, which resulted in inappropriate color assignment.
Fig. 5. The classi cation result shown in 3D using CLDA.
6.1.2. Target classi9cation The strength of the CLDA-based classi er was investigated and the comparative study has been conducted in Ref. [2]. Here, we demonstrate how to classify di2erent targets in a single image with color. Fig. 4(a) shows the classi cation results, where the 10 panel classes were represented by di2erent colors and background classes (grass eld and tree line) are in black. Fig. 5 is a 3-D classi ed image with the third dimension speci ed by target abundance fraction. The magnitudes of the left-most panels are the greatest because of their largest size, while the sizes of the right-most panels are very small so that their magnitudes are the lowest. The FVA was also applied for the classi cation, whose result is shown in Fig. 4(b). The colors were totally messed
6.1.3. Real-time implementation Figs. 6(a–d) depicts a sequence of line-by-line real-time processed images that discriminate these 10 panel classes in di2erent stages with the CLDA algorithm. It rst classi ed P1 in red and all the background objects in black, then it classi ed P2 and displayed it in yellow, after that P3 was classi ed in green, and so on. If we compare the nal on-line classi cation in Fig. 6(d) with the o2-line classi cation in Fig. 4(a), we nd that there are some tiny di2erences between them. In particular, the di2erence in P1 and P2 classi cation is perceivable, where the panel sizes of the classi ed P1 and P2 are smaller in Fig. 6(d). The reason is that the on-line classi er (here it is about the −1 calculation) is generated and updated based on the received pixels instead of all the pixels in the entire scene as in the o2-line case. However, the discrepancies brought about by the on-line processing are tolerable. The entire real-time simulation was nished very quickly. It took a Sun Ultra10 440 MHz workstation only about 50 s. 6.2. Experiment II In the second experiment with HYDICE data, the scene ◦ was rotated by 90 counterclockwise as shown in Fig. 7. The 10 panel classes now were ordered from left to right. P1 was in the rst column, P2 the second and so on. 1 m×1 m panels are in the rst row, 2 m × 2 m the second and 3 m × 3 m
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
9
Fig. 6. Real-time classi cation I.
The rst row of panels were classi ed rst in pre-assigned di2erent colors. Since they are of size 1 m ×1 m and smaller than a pixel area, the color shades are light. p1b , p2b , p3b , p4b and p5b in the second row are among the next group of classi ed panels, which were displayed in red, yellow, green, magenta and cyan as pre-assigned, respectively. After all the 10 panels in the second row were classi ed, those in the third row were to be classi ed. p1a , p2a , p3a , p4a and p5a were displayed in the same color as p1b , p2b , p3b , p4b and p5b as the same classes respectively. If we compared the nal image in Fig. 8(j) with Fig. 6(h), we found that these two results were very close. This experiment demonstrates that the performance of the proposed real-time algorithm is irrelevant to the data processing order. 7. Conclusion
Fig. 7. (a) The image scene used in Experiment II after rotating ◦ the image in Fig. 1(a) 90 counterclockwise; (b) Spatial locations of 30 panels.
the third. Only real-time implementation may be a2ected since the data are processed from left to right and top to bottom. The classi ed image sequence is shown in Fig. 8.
This paper presents CLDA-based target detection and classi cation algorithms as well as their real-time implementation for hyperspectral imagery. The CLDA-based target detector is derived by treating desired and undesired classes as two groups. As a result, no new sample means are needed to recalculate. The CLDA-based target classi er introduces the color for display such that great convenience is o2ered by visual discrimination of di2erent targets in a single image. All the CLDA-based detection and classi cation
10
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
Fig. 8. Real-time classi cation II.
algorithms have great potential to be implemented in real-time, which is very important in defense and intelligence related applications. Both pixel-by-pixel and line-by-line real-time processing algorithms are developed. The experiment results with HYDICE data demonstrated its excellent performance in hyperspectral image detection, classi cation and real-time implementation.
Acknowledgements The authors would like to thank Professor Chein-I Chang at the University of Maryland Baltimore County for providing the image data used in the experiments. This work was performed while the second author held a National Research Council Research Associateship Award at Edgewood Chemical and Biological Center, US Army.
Appendix A The following theorem for orthogonal decomposition of vectors is being used [13]: Theorem. Given ; a vector subspace of n-dimensional Euclidean space Rn ; every n × 1 vector y can be expressed uniquely in the form y = u + v; where u ∈ and v ∈ ⊥ . Here, the L × 1 vector w can be represented as w = wpfs + wofs = wpus + wous + wofs ;
(A.1)
where wpfs and wofs are the projections of w onto the subspace spanned by the signature vectors (feature subspace) and its orthogonal complement, respectively. wpfs can be further decomposed as the sum of its projections onto the undesired subspace wpus and its orthogonal complement wous .
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
Since the undesired signatures U are orthogonal to wous and wofs , the constraint in Eq. (9) becomes wT U = (wpus )T U = 0T :
(A.2)
11
If the partition is S = [S1 S2 ], then T S1 S1 ST1 S2 : ST S = ST2 S1 ST2 S2
(B.2)
Assumed that the signatures in U are linear independent. wpus in Eq. (34) must be a zero vector with all the elements equal to 0 such that Eq. (A.2) can hold. The next task is to determine wous and wofs in Eq. (A.1). The projector of the orthogonal complement of the undesired signature subspace can be constructed by PU⊥ =I−UU# , where U# = (UT U)−1 UT is the pseudo-inverse of U. Therefore, wous is proportional to PU⊥ and can be represented as ⊥ wous = PU a, where a is a vector to be determined. Since the desired signatures D are orthogonal to wofs , the constraint in Eq. (8) becomes
Comparing Eq. (B.2) with Eq. (B.1), A, B, D, E, F in Eq. (B.1) are
⊥ T ⊥ wT D = (wous )T D = (PU a) D = aT PU D = 1T :
F = (ST1 S1 )−1 ST1 S2 :
(A.3)
A = ST1 S1
(B.3)
B = ST1 S2
(B.4)
D = ST2 S2
(B.5)
E = ST2 S2 − ST2 S1 (ST1 S1 )−1 ST1 S2 = ST2 (I − S1 (ST1 S1 )−1 ST1 )S2 = ST2 PS⊥1 S2
(B.6) (B.7)
Substituting Eqs. (B.3) – (B.7) into Eq. (B.2), (ST S)−1 =
(ST1 S1 )−1 + (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ST2 S1 (ST1 S1 )−1
−(ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1
−(ST2 PS⊥1 S2 )−1 ST2 S1 (ST1 S1 )−1
(ST2 PS⊥1 S2 )−1
(B.8)
and T
(S S)
−1 T
T
S = (S S) =
ST1
ST2
(ST1 S1 )−1 ST1 + (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ST2 S1 (ST1 S1 )−1 ST1 − (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ST2
Sx1 + Sx2 − Sx3 (ST2 PS⊥1 S2 )−1 ST2 PS⊥1
;
(B.9)
⊥ Therefore, a can be estimated by aT = 1T (PU D)# = 1T (DT PU⊥ D)−1 DT PU⊥ . As for wofs , it is orthogonal to wous . In order to make w2 minimal, wofs should also be a zero vector. Hence,
w=w
ous
=
⊥ PU a
=
−(ST2 PS⊥1 S2 )−1 ST2 S1 (ST1 S1 )−1 ST1 + (ST2 PS⊥1 S2 )−1 ST2
=
−1
PU⊥ D(DT PU⊥ D)−1 1:
(A.4)
where Sx1 =(ST1 S1 )−1 ST1 , Sx2 =(ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ×ST2 S1 (ST1 S1 )−1 ST1 and Sx3 =(ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ST2 . Using the matrix inversion lemma and matrix operations: (ST1 PS⊥2 S1 )−1 ST1 PS⊥2 =(ST1 S1 − ST1 S2 (ST2 S2 )−1 ST2 S1 )−1 ST1 PS⊥2
Appendix B The following partitioned matrix inversion formula is used [13]: −1 −1 A B A + FE−1 FT −FE−1 ; (B.1) = BT D E−1 −E−1 FT where E = D − BT A−1 B, F = A−1 B.
=[(ST1 S1 )−1 + (ST1 S1 )−1 ST1 S2 (ST2 S2 − ST2 S1 ×(ST1 S1 )−1 ST1 S2 )−1 ST2 S1 (ST1 S1 )−1 ]ST1 PS⊥2 =[(ST1 S1 )−1 + (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ×ST2 S1 (ST1 S1 )−1 ]ST1 PS⊥2
12
Q. Du, H. Ren / Pattern Recognition 36 (2003) 1 – 12
=[(ST1 S1 )−1 ST1 + (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1
column vector with all the elements equal to 0, i.e., 0 = (0; 0; : : : ; 0)T .
×ST2 S1 (ST1 S1 )−1 ST1 ](I − S2 (ST2 S2 )−1 ST2 ) =[Sx1 + Sx2 ](I −
q
S2 (ST2 S2 )−1 ST2 )
=Sx1 + Sx2 − Sx1 S2 (ST2 S2 )−1 ST2 − Sx2 S2 (ST2 S2 )−1 ST2 =Sx1 + Sx2 − (ST1 S1 )−1 ST1 S2 [(ST2 S2 )−1 +(ST2 PS⊥1 S2 )−1 ST2 S1 (ST1 S1 )−1 ST1 S2 (ST2 S2 )−1 ]ST2 =Sx1 + Sx2 − (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 [ST2 PS⊥1 S2 (ST2 S2 )−1 + ST2 S1 (ST1 S1 )−1 ST1 S2 (ST2 S2 )−1 ]ST2 =Sx1 + Sx2 − (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 [ST2 PS⊥1 S2 (ST2 S2 )−1 + ST2 PS1 S2 (ST2 S2 )−1 ]ST2 =Sx1 + Sx2 − (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 [ST2 S2 (ST2 S2 )−1 ]ST2 =Sx1 + Sx2 − (ST1 S1 )−1 ST1 S2 (ST2 PS⊥1 S2 )−1 ST2 =Sx1 + Sx2 − Sx3 which means the rst element in Eq. (B.9) equals (ST1 PS⊥2 S1 )−1 ST1 PS⊥2 , i.e., T ⊥ −1 T ⊥ (S1 PS2 S1 ) S1 PS2 T −1 T (S S) S = : (B.10) (ST2 PS⊥1 S2 )−1 ST2 PS⊥1 If D = S1 and U = S2 are L × p and L × q matrices, respectively, T ⊥ −1 T ⊥ (D PU D) D PU T T T −1 T T T [1 0 ](S S) S = [1 0 ] · (UT PD⊥ U)−1 UT PD⊥ = 1T (DT PU⊥ D)−1 DT PU⊥ ;
(B.11)
where 1 is a p × 1 column vector with all the elements equal to 1, i.e., 1 = (1; 1; : : : ; 1)T and 0 is a q × 1 p
References [1] A.F.H. Goetz, G. Vane, J.E. Solomon, B.N. Rock, Imaging spectrometry for earth remote sensing, Science 228 (4704) (1985) 1147–1153. [2] R.O. Duda, P.E. Hart, Pattern Classi cation Scene Anal, Wiley, New York, 1973. [3] Q. Du, C.-I Chang, Linear constrained distance-based discriminant analysis for hyperspectral image classi cation, Pattern Recognition 34 (2) (2001) 361–373. [4] J.C. Harsanyi, C.-I Chang, Hyperspectral image classi cation and dimensionality reduction: an orthogonal subspace projection IEEE Trans. Geosci. Remote Sensing 32 (4) (1994) 779–785. [5] W.H. Farrand, J.C. Harsanyi, Mapping the distribution of mine tailing in the coeur d’Alene river valley, Idaho through the use of constrained energy minimization technique, Remote Sensing Environ. 59 (1997) 64–76. [6] J. Bowles, P. Palmadesso, J. Antoniades, M. Baumback, Use of lter vectors in hyperspectral data analysis, Proc. SPIE 2553 (1995) 148–152. [7] C.-I Chang, Hyperspectral Imaging: Subpixel Detection and Classi cation,KluwerAcademicPublishers, Dordrecht, 2002. [8] C.M. Stellman, G.G. Hazel, F. Bucholtz, J.V. Michalowicz, Real-time hyperspectral detection and cuing, Opt. Eng. (2000) 1928–1935. [9] C.-I Chang, H. Ren, S.S. Chiang, Real-time processing algorithms for target detection and classi cation in hyperspectral imagery, IEEE Trans. Geosci. Remote Sensing 39 (4) (2001) 760–768. [10] H. Ren, C.-I Chang, An unsupervised orthogonal subspace projection to target detection and classi cation in unknown environment, Spectroradiometric Symposium, San Diego, November 2–7, 1997. [11] D.C. Heinz, C.-I. Chang, Fully constrained least squares linear mixture analysis for material quanti cation in hyperspectral imagery, IEEE Trans. Geosci. Remote Sensing 39 (3) (2001) 529–545. [12] S. Haykin, Adaptive Filter Theory, 3rd Edition, Prentice-Hall, Englewood Cli2s, NJ, 1996. [13] R.J. Muirhead, Aspects of Multivariate Statistical Theory, Wiley, New York, 1982.
About the Author—QIAN DU received the Ph.D. degree in Electrical Engineering from University of Maryland Baltimore County in 2000. She is currently an assistant professor at the Department of Electrical Engineering and Computer Science, Texas A & M University-Kingsville. Her research interests include signal and image processing, pattern recognition, neural networks and remote sensing. She is a member of IEEE, SPIE and the honor society of Phi Kappa Phi. About the Author—HSUAN REN received the B.S. degree in Electrical Engineering from the National Taiwan University, Taipei, Taiwan, R.O.C., in 1994 and the M.S. and Ph.D. degrees from University of Maryland Baltimore County in 1998 and 2000, respectively, all in Electrical Engineering. He is currently on National Research Council (NRC) Associateship and supported by U.S. Army Edgewood Chemical Biological Center as a postdoctoral fellow. He has two patents pending on hyperspectral target detection and image classi cation. Dr. Ren is also a part-time instructor at the Department of Computer Science and Electrical Engineering, University of Maryland Baltimore County. His research interests include data compression, signal and image processing and pattern recognition. He is a member of IEEE and SPIE. He is also a member of the honor society of Phi Kappa Phi.
