Unsupervised Linear Unmixing for Hyperspectral Data Exploitation Mingkai Hsueh

1. Unsupervised Linear Unmixing for Hyperspectral Data Exploitation Mingkai Hsueh Remote Sensing Signal and Image Processing Laboratory Department of Computer Science and Electrical Engineering University of Maryland Baltimore County 1000 Hilltop Circle, Baltimore, MD 21250

2. Outline • Introduction to Hyperspectral Image Processing and its Applications • Endmember Extraction • Pixel Purity Index Algorithm (PPI) • Block of Skewers (BOS) based PPI • Anomaly Detection • Anomaly Detection Algorithms and its real-time implementation • Speed-up of Adaptive Causal Anomaly Detection • Conclusions

3. 4000 3000 Reflectance 2000 1000 0 300 600 900 1200 1500 1800 2100 2400 Wavelength (nm) 4000 3000 Reflectance 2000 1000 0 300 600 900 1200 1500 1800 2100 2400 Wavelength (nm) 5000 4000 3000 Reflectance 2000 1000 0 300 600 900 1200 1500 1800 2100 2400 Wavelength (nm) Hyperspectral Image Mixed pixel (soil + mineral) Water Mixed pixel (trees + soil)

4. Applications of Hyperspectral Image Processing • Applications • Man-made objects: canvas, camouflage, military vehicles in defense applications • Toxic waste, oil spills in environmental monitoring • Landmines • Trafficking in law enforcement • Chemical/biological agent detection • Special species in agriculture, ecology

5. Types of Signatures • Endmembers: Pure signatures for a spectral class used for spectral unmixing • Anomalies: Signals/signatures spectrally distinct from their surroundings, i.e., abnormality. • rare minerals in geology • abnormal activities in military applications.

6. Part I : Endmember Extraction

7. Endmember Extraction • An endmember pixel is defined as a pixel with idealized, pure spectral signature for a class.

8. Pixel Purity Index (PPI) • The idea of PPI was first proposed by Boardman and has been one of most popular endmember extraction algorithms (EEAs) due to its publicity and availability in ENVI software. • For the PPI to work effectively, a large number of dot products between skewers (random vectors) and data sample vectors are required.

9. PPI Algorithm NPPI(e2)=0 NPPI(e2)=1 e2 Maximum Projection skewer2 Maximum Projection skewer3 skewer1 Minimum Projection e3 NPPI(e3)=1 NPPI(e3)=3 NPPI(e3)=0 NPPI(e3)=2 e1 Minimum Projection NPPI(e1)=2 NPPI(e1)=0 NPPI(e1)=1 Minimum Projection

10. Block of Skewer – An Example • Given K1 & K2 are skewers, K3 ~ K6 are linear combinations of K1 & K2 and “r” is a pixel vector from hyperspectral image cube. K3= – K1 – K2 K4= – K1 + K2 K5= + K1 – K2 K6= + K1 + K2 • In stead of generating more skewers for projection, we perform linear combination between the projection results of the generated skewers.

11. C-BOS Dskewer = a1Iskewer1 + a2Iskewer2 + a3Iskewer3 • Ideally, a1, a2 and a3 can be any real numbers. However, given limited hardware resource, fixed-point implementation are usually preferred. • Here we constrain the coefficients to +1 or -1 to form 8 combinations as a cube shown on the left. (1, -1, 1) (-1, -1, 1) (-1, 1, 1) (1, 1, 1) z (0, 0, 0) x y (1, -1, -1) (-1, -1, -1) (-1, 1, -1) (1, 1, -1)

12. P-BOS (0, 1, 0) (1, 0, 0) (0, 0, -1) (0, 0, 1) (-1, 0, 0) (0, -1, 0)

13. Skewer Redundancy • In the previous example, there exist redundancy among skewers. K3= – K1 – K2 K6= K1 + K2 K4= – K1 + K2 K5= K1 – K2 • Same thing happens to the C-BOS. D1= + I1 + I2 +I3 D5 = – I1 – I2 –I3 D2= + I1 + I2 –I3 D6 = – I1 – I2 +I3 D3= + I1 – I2 +I3 D7 = – I1 + I2 –I3 D4= + I1 – I2 –I3 D8 = – I1 + I2 +I3

14. S-BOS • Decompose cube into 6 squares (-1, -1, 1) (1, -1, 1) (-1, -1, 1) (1, -1, 1) (1, 1, 1) (1, -1, 1) Right Back Top (-1, -1, -1) (1, -1, -1) (-1, 1, 1) (1, 1, 1) (1, 1, -1) (1, -1, -1) (-1, 1, 1) (1, 1, 1) (-1, -1, -1) (1, -1, -1) (-1, 1, 1) (-1, -1, 1) Front Bottom Left (-1, 1, -1) (1, 1, -1) (-1, 1, -1) (1, 1, -1) (-1, 1, -1) (-1, -1, -1)

15. T-BOS • A tetrahedron shape can be considered as half of a pyramid. • A cube can be shifted so that coordinates of 8 vertices are holding values 0 or 1. z (1, 0, 1) (0, 0, 1) (0, 1, 1) (1, 1, 1) (1, 0, 0) (0, 0, 0) x (0, 1, 0) (1, 1, 0) y

16. HYDICE Data • HYDICE (Hyperspectral Digital Imagery Collection Experiment) • 15 panels of five types with three different materials. • They are arranged into a matrix in such a way that each row represents 3 panels of the same type with three different sizes, 3m3m, 2m2m, 1m1m. Each column represents 5 panels of different types with the same size. Anomaly Original image Target masked image

17. Experiments with Real Image (Cont’d) • Experimental results of HYDICE real image (a) C-BOS (b) S-BOS (c) P-BOS (d) T-BOS

18. A large amount of independent dot products make it particularly suitable for FPGA implementation due to the readily parallel design architecture. BOS based PPIHardware module Dot product module MINMAX MINMAX Dskewer Generator . . . . . . MINMAX

19. Dot-Product Module 2nd pixel vector 1st band 2nd band 3rd band 4th band 1st pixel vector 1st band 2nd band 3rd band 4th band Iskewer 1 Dskewer Generator PE unit PE PE PE Iskewer 2 PE unit PE PE PE Iskewer 3 PE unit PE PE PE

20. S-BOS Dskewer Generators P1 P1 P2 P2 P3 P3 –P1 – P2 + P3 P1 –P2 + P3 –P1 + P2 – P3 –P1 – P2 – P3

21. FPGA Implementation • Four different Dskewer generators are implemented in XESS XSB-300E board which carries a Spartan II E (XC2S300E) FPGA.

22. Computational Complexity • The computational complexity is calculated based on the number of multiplications and additions performed with different BOS design. • K is the total number of skewers, L isthe number of spectral bands and N is the total number of pixel vectors.

23. Part II : Anomaly Detection

24. RX Algorithm • RX algorithm basically performs the Mahalanobis distance that is specified by (ri-)T× (K)-1 × (ri -) • The required mean vector μhinder the possibility of implementing the algorithm in real-time fashion.

25. Causal RX Filter (CRXF) • By replacing the covariance matrix by correlation matrix, we can achieve the real-time processing. • The functional form of CRXF riT× (Ri)-1 × ri • The major drawback is that if a detected anomaly remains on the image to be processed, it may decrease the detectability of the following anomalies.

26. Adaptive Causal Anomaly Detector (ACAD) • ACAD has the same functional form as does CRXF, except the sample correlation matrix R’ is formed by all the arrived pixel vectors except the detected anomalous target pixel vectors that have been removed. riT× (R’i)-1 × ri • An anomalous target map is generated at the same time as the detection process takes place.

27. HYDICE Data • HYDICE (Hyperspectral Digital Imagery Collection Experiment) • 15 panels of five types with three different materials. • They are arranged into a matrix in such a way that each row represents 3 panels of the same type with three different sizes, 3m3m, 2m2m, 1m1m. Each column represents 5 panels of different types with the same size. Anomaly Original image Target masked image

28. CRXF Results row 8 row 16 row 24 row 32 row 40 row 48 row 56 row 64

29. ACAD Results row 8 row 16 row 24 row 32 row 40 row 48 row 56 row 64

30. ACAD Target Map row 8 row 16 row 24 row 32 row 64 row 40 row 48 row 56

31. ACAD Hardware Design Ri= Ri-1+ri× riT Auto Correlator (Ri)-1= (Qi × Riupper )-1 = ( Riupper)-1× QiT QR Matrix Inverse Abundance Calculation δACAD(ri)= riT×(RiT)-1×ri Anomalous Target Discriminator tK ≤τ

32. Matrix Inversion Lemma • Let A be the current correlation matrix and r be the incoming pixel vector. (A+BCD)-1 = A-1 – A-1B(C-1+DA-1B)-1 DA-1 By Woodbury’s identity, set B a column vector, C a scalar of unity, and D a row vector  (A+rrT)-1 = A-1 – (A-1rrT A-1) / (1+rTA-1r)

33. Matrix Inversion Lemma (Cont’d) • With Matrix Inversion Lemma (MIL), we only need to compute • Using MIL the matrix inversion is reduced to matrix multiplications. • Simulation is provided to evaluate the performance of MIL.

34. ACAD Hardware Design Ri= Ri-1+ri× riT Auto Correlator (Ri)-1= (Qi × Riupper )-1 = ( Riupper)-1× QiT QR Matrix Inverse Abundance Calculation δACAD(ri)= riT×(RiT)-1×ri Anomalous Target Discriminator tK ≤τ Matrix Inversion Lemma

35. Speed-up of MIL • We use two versions of the MATLAB program to perform the ACAD on the same image cube. One uses the MATLAB inv() function and another one uses the MIL. • As we can see, the speed-up is about “2” times faster for the 64x64 HYDICE image than the one without MIL.

36. Conclusions • The Matrix Inversion Lemma has been successfully applied to reduce the matrix inversion performed by Adaptive Causal Anomaly Detection (ACAD) into matrix multiplications. • Since the Causal RX Filter (CRXF) and Real-time CEM (Constrained Energy Minimization) previously proposed in Wang [2003] also involve inverse matrix computation, the same MIL-based approach can be also applied to reduce the computational load.

37. Conclusions (Cont’d) • New block design of BOS including Square based and Tetrahedron based BOS have been introduced to improve the drawback of the Pyramid-based and Cube-based BOS design. • The FPGA design and implementation of the four BOS design has also been evaluated and analyzed.

38. Future Work • An effective Dimensionality Reduction (DR) or Band Selection (BS) may need to reduce the number of bands to an acceptable range so that we can further reduce the computation cost in both applications. • Heterogeneous platform may be also considered to reduce the design time and possibly achieve better performance.

39. Projects Conducted in RSSIPL • Joint Service Agent Water Monitor • Mission • Develop GUI image analysis software for detecting Biological Threat Agent on Handheld Assays • Ported developed algorithms onto embedded system, Stargate Gateway (SPB400, Linux single board computer) with external hand held scanner device. • Sponsor • US Army Edgewood Chemical and Biological Center (ECBC) • ANP Technologies, Inc.

40. Software for Detecting Agents

41. Projects Conducted in RSSIPL (Cont’d) • Multi-band Multi-threat warning sensor • Mission • Developed detection algorithms for missile and grenade images captured from real-time Multispectral imaging system. • Developed MATLAB based GUI for image analysis. • Sponsor • Surface Optics Corporation (SOC)

42. Projects Conducted in RSSIPL (Cont’d)

43. Projects Conducted in RSSIPL (Cont’d) • Multi-band Multi-threat warning sensor • Mission • Developed detection algorithms for missile and grenade images captured from real-time Multispectral imaging system. • Developed MATLAB based GUI for image analysis. • Sponsor • Surface Optics Corporation (SOC)

44. Publication • Book Chapter • J. Wang, M. Hsueh and C.-I Chang, “FPGA Design for Second-order Statistics Based Target Detection Algorithm for Hyperspectral Imagery Applications,” High Performance Computing in Remote Sensing, Chapman & Hall/CR, Oct 2007. • J. Wang, M. Hsueh and C.-I Chang, “FPGA Implementation for Real-time Orthogonal Subspace Projection for Hyperspectral Imagery Applications,” High Performance Computing in Remote Sensing, Chapman & Hall/CR, Oct 2007.

45. Publication (cont’d) • Journal • C.-I Chang and M. Hsueh, “Characterization of Anomaly Detection in Hyperspectral Imagery,” Sensor Review, Volume 26, Issue 2, pp. 137-146, 2006. • M. Hsueh and C.-I Chang, “Field Programmable Gate Arrays for Pixel Purity Index Using Blocks of Skewers for Endmember Extraction in Hyperspectral Imagery,” International Journal of High Performance Computing Applications, Dec 2007. (to appear) • C.-I Chang, M. Hsueh, F. Chaudhry, W. Liu, C.-C. Wu, G. Solyar, “A pyramid-based block of skewers for pixel purity index for endmember Extraction in hyperspectral imagery,” International Journal of High Speed Electronics and Systems. (to appear) • M. Hsueh and C.-I Chang, “Adaptive Causal Anomaly Detection on Reconfigurable Computing,” IEEE Transaction on Industrial Electronics. (To be submitted)

46. Publication (cont’d) • Conference • M. Hsueh and C.-I Chang, “FPGA implementation of Adaptive Causal Anomaly Detection,” 2006 CIE Annual Convention, Newark, NJ, Sep 16, 2006. • C.-I Chang, M. Hsueh, F. Chaudhry, W. Liu, C. C. Wu, A. Plaza and G. Solyar, “A Pyramid-based Block of Skewers for Pixel Purity Index for Endmember Extraction in Hyperspectral Imagery,” 2006 International Symposium on Spectral Sensing Research, Bar Harbor, ME, May 29 to Jun 2, 2006. • D. Valencia, A. Plaza, M. A. Vega-Rodriguez, R. M. Perez and M. Hsueh, “FPGA Design and Implementation of a Fast Pixel Purity Index Algorithm for Endmember Extraction in Hyperspectral Imagery,” SPIE Optics East, Boston, MA, Oct 23-26 2005. • L. Wu, J. Wang, B. Ramakrishna, M. Hsueh, J. Liu, Q. Wu, C. Wu, M. Cao, C. Chang, J. L. Jensen, J. O. Jensen, H. Knapp, R. Daniel, R. Yin, “An embedded system developed for hand held assay used in water monitoring,” SPIE Optics East, Boston, MA, Oct 23-26, 2005.

47. Publication (cont’d) • Conference • M. Hsueh and C.-I Chang, “Adaptive Causal Anomaly Detection for Hyperspectral Imagery”, IEEE International Geoscience and Remote Sensing Symposium, Alaska, Sep 19-26, 2004. • M. Hseuh, A. Plaza, J. Wang, S. Wang, W. Liu, C.-I Chang, J. L. Jensen and J. O. Jensen, “Morphological algorithms for processing tickets by hand held assay,” OpticsEast, Chemical and Biological Standoff Detection II (OE120), Vol. 5584, Philadelphia, PA, Oct 25-28, 2004. • C.-I Chang, H. Ren, M. Hsueh, F. D’Amico and J.O. Jensen, “A Revisit to Target-Constrained Interference-Minimized Filter,” 48th Annual Meeting, SPIE International Symposium on Optical science and Technology, Imaging Spectrometry IX ( AM110), San Diego, CA, Aug 3-8, 2003. • S. T. Sheu, M. Hsueh, “An Intelligent Cell Checking Policy for Promoting Data Transfer Performance in Wireless ATM Networks,” IEEE ATM Workshop '99, Kochi City, Kochi, Japan, May 24-27, 1999.

48. Thank you!!