# Persistent Homology on Grassmann Manifolds for Analysis of Hyperspectral Movies

###### Abstract

The existence of characteristic structure, or shape, in complex data sets has been recognized as increasingly important for mathematical data analysis. This realization has motivated the development of new tools such as persistent homology for exploring topological invariants, or features, in large data sets. In this paper we apply persistent homology to the characterization of gas plumes in time dependent sequences of hyperspectral cubes, i.e. the analysis of 4-way arrays. We investigate hyperspectral movies of Long-Wavelength Infrared data monitoring an experimental release of chemical simulant into the air. Our approach models regions of interest within the hyperspectral data cubes as points on the real Grassmann manifold (whose points parameterize the -dimensional subspaces of ), contrasting our approach with the more standard framework in Euclidean space. An advantage of this approach is that it allows a sequence of time slices in a hyperspectral movie to be collapsed to a sequence of points in such a way that some of the key structure within and between the slices is encoded by the points on the Grassmann manifold. This motivates the search for topological features, associated with the evolution of the frames of a hyperspectral movie, within the corresponding points on the Grassmann manifold. The proposed mathematical model affords the processing of large data sets while retaining valuable discriminatory information. In this paper, we discuss how embedding our data in the Grassmann manifold, together with topological data analysis, captures dynamical events that occur as the chemical plume is released and evolves. Keywords: Grassmann manifold, persistent homology, hyperspectral imagery, signal detection, topological data analysis

## 1 Introduction

Hyperspectral imaging (HSI) technology allows the acquisition of information across the electromagnetic spectrum that is invisible to humans. In a very real sense, these cameras allow us to “see the unseen” by including wavelengths spanning ultraviolet and far infrared. In contrast, humans can observe a very limited range of the electromagnetic spectrum, i.e. wavelengths of approximately 400-700nm are visible to the human eye.

Multi- and hyper-spectral imaging technology has become widely available, and there is an increasing number of canonical data sets available for scientific analysis including, e.g. the AVIRIS Indian Pines^{1}^{1}1Available from https://engineering.purdue.edu/~biehl/MultiSpec. and the ROSIS University of Pavia^{2}^{2}2Available from http://www.ehu.es/ccwintco/index.php/ data sets. In addition, moving objects may be detected with devices such as the Fabry-Pérot Interferometer[10] which can capture hyperspectral movies at frame rates of approximately 5Hz. See Figure 1 for an illustration. The resulting 4-way arrays of spatial-spectral-temporal data provide a high fidelity view of our environment and may help in the monitoring of pollution in the air and water. An application that concerns us in this paper is the characterization of gaseous plumes as they are released into the environment.

Traditionally, one of the primary applications of hyperspectral image analysis consists of object detection and classification. The focus is generally on the identification of anomalous pixels in the image and the determination of the composition of the materials in the pixel. A range of mathematical tools have been developed for the analysis of hyperspectral images including, e.g. matched subspaces, the RX algorithm, and the adaptive cosine estimator[19]. More recently, manifold learning algorithms have been applied to hyperspectral images to exploit topology and geometry, i.e. mathematical shape, or signatures, in data at the pixel level [1, 18].

The subspace perspective is also taken in this paper, but in the direction of understanding the topology and geometry of the Grassmann manifold (Grassmannian) associated with hyperspectral images, i.e. the manifold parameterizing the -dimensional subspaces of -dimensional space. While we are motivated by ideas similar to those found in prior applications of manifold learning algorithms, e.g. [1, 18], our application data is not at the pixel level. By constructing subspaces of pixels we are able to exploit the rich metric structure of the Grassmannian based on measuring angles between subspaces. The advantage of this approach is that a set of pixels used to form a subspace is seen to capture the variability in the data missing in a single pixel observation.

An example that illustrates the power of this framework is the application to illumination spaces in the face recognition problem. The variation in illumination on an object may be approximated by a cone captured in a low-dimensional subspace. Subspace angles can be used to compute similarity of illumination spaces and the effect on classification accuracy was striking when applied to the CMU-PIE data set, even on ultra-low resolution images [4]. More recently, tools have been developed to represent points on Grassmannians via subspace means [20], or nested flags of subspaces [12]. In another application to video sequence data, we used the setting of the Grassmannian to extend an algorithm on vector spaces for detection of anomalous activities [25].

In this paper, we address the question of the existence of topological signatures in the setting of hyperspectral movies mapped to the Grassmannian. Our approach builds on applying the Grassmannian architecture to hyperspectral movies that has shown promise in preliminary work [7, 6]. Here, our focus is on application of persistent homology (PH) to the characterization of the evolution of chemical plumes as acquired by hyperspectral movie data sets. As in the application to face recognition, we encode a single frame of a hyperspectral movie (or a collection of pixels of a single frame in the movie) as a point on the Grassmann manifold. We speculate that this manifold representation affords a form of compression of information while capturing essential topological structure. We consider the application of this approach to the characterization of chemical signals as measured by the Long-Wavelength Infrared (LWIR) data set [10]. Our goal is to establish the existence of topological signatures that can provide insight into the evolution of complex 4-way data arrays.

## 2 Persistent Homology

Homology is an invariant that measures features of a topological space and can be used to distinguish distinct spaces from one another [16]. Persistent homology encodes a parameterized family of these homological features. It is a computational approach to topology that allows one to answer basic questions about the structure of point clouds in data sets at multiple scales [3]. This procedure involves (1) interpreting a point cloud as a noisy sampling of a topological space, (2) creating a global object by forming connections between proximate points based on a scale parameter, (3) determining the topological structure made by these connections, and (4) looking for structures that persist across different scales. PH has been used to understand the topological structure of data arising from applications including [17, 21, 11, 8, 22, 24]. For a detailed discussion of homology, see [16], and for further discussions of persistent homology, see [15, 14, 3].

One way to associate a family of topological objects with a point cloud is to use the points to construct a family of nested simplicial complexes. The Vietoris-Rips complex builds a simplicial complex with vertices as the data points and higher dimensional -simplices formed whenever points have pairwise distances less than . The -dimensional holes of this simplicial complex generate a homology group whose rank, known as the -th Betti number, counts the number of -dimensional holes. For instance, measures the number of connected components (clusters) of the point cloud, while indicates the existence of topological circles (loops), or periodic phenomenon. To avoid picking a specific scale , persistent homology seeks structures that persist over a range of scales, exploiting the fact that as grows, so do the simplicial complexes indexed by the parameters . Thus, PH tracks homology classes of the point cloud along the scale parameter, indicating at which a th order hole appears and for which range of values it persists. The Betti numbers, as functions of the scale , can be visualized in a distinct barcode for each dimension [15].

Figure 2 is an example of the and barcodes generated for a point cloud sampled from a circle. Each horizontal bar begins at the scale where a topological feature first appears and ends at the scale where the feature no longer remains. The th Betti number at any given parameter value is the number of bars that intersect the vertical line through . Short-lived features are often considered as noise while those features persisting over a large range of scale represent true topological characteristics. In the case of , at small values of there will be a distinct bar for each point, as the simplicial complex consists of isolated vertices. At large values of , only one bar remains, as all data will eventually connect into a single component. For the circle, which correspond to the number of connected components and number of loops, respectively, shown by the longest (persistent) horizontal bars in each plot. We use JavaPlex, a library for computing PH and TDA in this paper[23].

## 3 The Geometry of the Grassmann Manifold

The (real) Grassmann manifold is a parameterization of all -dimensional subspaces of -dimensional space [13]. A point on can be represented by a tall matrix with the property that where is an element of the equivalence class consisting of all matrices of the form with , the orthogonal group that consists of orthogonal matrices [13].

Hyperspectral data is a 3-way cube that can be mapped to points in a Grassmannian in a variety of ways. Here, we select a subset of frequencies . For each of the frequencies we propose to “vectorize” the spatial components to form an matrix . It is assumed that the construction is such that so subspaces don’t overlap trivially. To map to a matrix representing a point on the Grassmann we compute any orthogonal basis for the column space of . For instance, the matrix in the thin singular value decomposition provides one option as a representation of a point on the Grassmanian .

The mapping described above allows us to construct a sequence of points on , each one taken from the same spatial location in the 3-way array of hyperspectral pixels or from the same frame of a hyperspectral movie. The pairwise distances between the points in this sequence are computed in terms of the principal angles between the subspaces. The implementation of the Grassmannian framework is, in part, motivated by the rich metric structure of a variety of distance measures including the chordal, geodesic, and Fubini-Study distances, which are all functions of the principal angles between the subspaces [2, 9].

The experiments in this paper use the (pseudo)distance between two subspaces measured by the smallest principal angle. This (pseudo)distance has been effective in other numerical experiments [4, 6], and in fact, we observed, in the experiments in this paper, that using it resulted in stronger topological signals than other distance measures. Once a distance matrix for the points on the Grassmannian is computed, we apply PH to determine topological structure. In particular, we explore barcodes to estimate the number of connected components and barcodes to detect topological circles. The goal is to associate physical properties in the HSI image that relate to these structures.

## 4 Experimental Results

In this section, we apply PH to Long-Wavelength Infrared (LWIR) multispectral movies, each of an explosive release of a chemical and resulting toxic plume which travels across the horizon of the scene [10]. The simulants released included Triethyl Phosphate (TEP) and Methyl Salicylate (MeS) in quantities of 75kg and 150kg, respectively. The LWIR data sets are captured using an interferometer in the 8-11 range of the electromagnetic spectrum. A single frame, or data cube, of this movie consists of pixels collected at IR bands. A given movie is a sequence of data cubes consisting of pre-burst and post-burst frames.

The purpose of this paper is not to propose a new algorithm for detecting chemical plumes but rather to investigate the topological features associated with a known plume. The data processing workflow consists of the following steps: (1) band selection, (2) identification of the region containing the chemical plume, (3) mapping data to the Grassmannian, (4) computing (pseudo)distances on the Grassmannian using the smallest principal angle, (5) determination of PH and barcodes, and, finally, (6) interpretation of the structure in the data as encoded by the topological invariants. We describe more detail of steps (1) and (2) below.

Band Selection. We applied the sparse support vector machine (SSVM) algorithm for optimal in situ band selection, i.e. the SSVM identifies wavelengths that best discriminate the plume from the natural background [5]. In another approach, we visually choose bands which have the strongest plume signal in data cubes which have had the background removed and thus, have visible plume.

Plume detection. The location of the chemical plume in the post-burst cubes is determined using the well-known adaptive-cosine-estimator (ACE) [19]. The ACE detector is one of the benchmark hyperspectral detection algorithms which computes the squared cosine of the angle between the whitened test pixel and the whitened target’s spectral signature. Based on a chosen threshold, an ACE score indicates if the chemical is present in the test pixel. Figure (a)a depicts an image corresponding to a cube without a plume, and Figure (b)b depicts a cube with a chemical plume detected by the ACE.

### 4.1 Experiment on Triethyl Phosphate Movie

We first consider the 561 frame multispectral movie of the data collection event of chemical Triethyl Phosphate (TEP) being released into the air. The data consists of the raw, unpreprocessed data including background clutter. It was determined that the wavelengths (nm) were optimal for discriminating TEP from background using the SSVM band selection algorithm. In this experiment we determine barcodes using all 561 TEP cubes, where subcubes have been extracted from regions of each data cube along the plume location region.

The barcode in Figure (a)a arises from the 561 Grassmanian points corresponding to the left horizon region in each data cube, limited by pixel rows 124 to 127 and pixel columns 34 to 41. This region belongs to the area when a plume forms and first becomes visible at frame 112 as detected by the ACE. At scale , there are 31 bars corresponding to 31 connected components on , with 28 isolated points from frames 111 to 142, one cluster containing frames {134, 135, 137}, one cluster containing frame 519, and another containing all other frames. At scale , we have 19 bars corresponding to 19 connected components on , with 18 isolated frames from 112 to 129, and one cluster containing all the rest. These bars persist for a large range of parameter value (to just beyond ), indicating a large degree of separation. At , we have 13 clusters with 11 isolated frames 112, 114 to 118, 120 to 123 and 125, one cluster of frames {119, 124}, and another containing everything else; see also [7]. Cubes following frame 111 are where the plume first occurs with the highest concentration of chemical, and the composition of the plume changes quickly as time progresses. PH detects separation of these cubes from pre-plume cubes and those cubes where plume no longer remains at multiple scales.

After the plume is released, the plume drifts to the right in the multispectral movie as time progresses. We now consider a plume patch corresponding to a horizon region located to the right of the original plume location discussed above. That is, a patch is drawn from pixel rows 124 to 127 and pixel columns 75 to 82 for each of the 561 data cubes in the TEP movie. This data is embedded in , and PH is implemented to uncover the structure of the data. Figure (b)b contains the -dimensional barcode. Analyzing connected components as varies, we observe that they differ from those found in the previous experiment, see Figure (a)a. At scale , we have 52 connected components on corresponding to 47 isolated points from 119 to 141, 145 to 165, and 170 to 172. The other points are connected into four smaller clusters {142,143,144}, {166,167}, {168,169}, and {173,174}, and one cluster containing all the other points. At scale , there are 30 connected components on the Grassmannian, including 25 isolated points from 119 to 127, 129 to 140, 149, and 151 to 156; four clusters each containing {128,136-138}, {141-150}, {157,158}, {162-164}; and one cluster containing all the rest. Further, at scale , the barcode plot has 5 bars that persist over a large range of values, namely, up to a little beyond : 4 isolated points from frame 121 to 124 and one cluster containing all the rest.

We observe that for this region, PH separates points from frame 119 and later, in contrast to the frames separated from frame 112 in the previous experiment (Figure (a)a). Note that points corresponding to frames 112 to 118 are “plume-free” as the plume does not reach this region until frame 119. It is also interesting to note that points corresponding to frames 121 to 124 are kept isolated for a large range of scales, i.e. they are far away from each other and the rest of the points. PH, under the Grassmannian framework, treats these frames as being the most distinct in this region.

### 4.2 Experiments Detecting A Loop in Methyl Salicylate Movie

The next two experiments consider the multispectral movie of the data collection event of chemical Methyl Salicylate (MeS) being released into the air, consisting of 829 frames. Here we use 3 out of 20 wavelength bands {10.57,10.68,10.94} (nm) that were determined by visual inspection of a background-removed data cube where plume was present. These bands, in particular, were selected as strong plume signal was visible at these corresponding wavelengths. In this movie, the plume first becomes visible at frame 32.

In the first experiment, we construct a sliding window along the horizon, where the plume is released, in both a frame with and without a plume present (frames 32 and 1, respectively) to compare the topological structure of each. This sliding window is constructed by selecting patches of each frame limited by rows 125-128 and columns 190-245 where each new point samples 8 columns, incrementing by one. Each patch is then embedded into and the topological structure is analyzed with PH. In this experiment and the next, our focus is on the information which measures the number of loops present in the data.

Observe in Figure 5 that no persistent topological circle is present in the barcode of frame 1, while a persistent loop is present in the barcode of frame 32. This is interpreted as follows. In frame 32, where a plume is present, the sliding window first constructs points in of the natural background, then traverses through points that contain plume, finally returning to points of the natural background. This creates a closed loop in . This behavior is captured in the topological structure of the plume cube. On the other hand, the sliding window in frame 1 only has points in of the natural background, and thus, no persistent loop is formed in this space.

We mention that this experiment was done on background removed frames. In analysis with raw data, loops were not as prevalent with this framework. However, the next experiment does in fact use raw data in our analysis.

In the second experiment, we consider the first one hundred frames of the MeS cubes and focus on a “plume location” patch of size , limited by pixel rows 125 to 128 and pixel columns 217 to 224, embedded into for each cube.

Figure (a)a displays the and barcodes from applying PH to this Grassmannian data. A fairly persistent bar appears in the barcode that begins at and ends at . This represents a loop through the data in . All other bars are considered as noise. Let us inspect this loop further. It begins once all of the data has been connected into a single component (refer to in the barcode). The maximum pairwise distance–measured by the smallest angle between subspaces–for this data is 0.0308. This loop persists until just under half this distance.

We conclude the following from this experiment. The first few frames start with a fixed background. Then, the plume begins to form, spreading through the plume patch until the plume no longer remains in the sampled region. The remaining cubes then return to a fixed background, reflecting periodic behavior in the data. This collection of cubes traces out a closed loop, encoded in the Grassmann manifold . PH captures this loop in the persistent bar. Figure (b)b displays a schematic of one possibility in the equivalence class of the edges that form this loop. While not all data cubes are present, we notice that those cubes immediately following 31 connect to one another sequentially. This is when the chemical is first released and begins to evolve. Cubes before this frame (where no plume is present) do not follow sequentially and connect with later cubes which no longer contain plume in the sampled plume patch. That is, the time dependent information of ‘pre-plume’ and ‘post-plume’ cubes–which simply contain information about the natural background and not the evolving plume–is not as important as ‘plume’ cubes.

## 5 Conclusion

We propose a geometric and topological model for capturing dynamical changes in hyperspectral movies. The HSI data cubes (or a sequence of pixel patches) are viewed as a sequence of points on the Grassmann manifold. The tools of persistent homology are then applied to capture topological novelty in the setting of the Grassmann manifold. This approach models cubes as points, a technique that permits the processing of potentially large amounts of data while retaining basic dynamical structure.

The dynamic structure recorded by the multispectral movie of the gas plume consisting of the simulant Triethyl Phosphate was illuminated in the barcodes. Frames containing the plume were identified as topological singletons, i.e. isolated points on the manifold for large ranges of scale. Grassmann points before the release, as well as long after the release, appeared as clusters of points. At a location to the right of this region, we see that later frames had a similar behavior, indicating that the geometric model of the Grassmannian allows the dynamics of the scene to be effectively characterized in a topological sense.

In the next two experiments, we use the barcode on the movie of the release of Methyl Salicylate mapped to the Grassmannian to reveal that a closed loop is present on the manifold, again reflecting the evolution of the plume. First, we consider a sliding window of pixels along the plume location region and observe that a loop is present in a frame with a plume unlike a frame without a plume. Second, we consider a patch of pixels in each of the first one hundred frames and observe a closed loop that encompasses frames immediately following the release of the chemical in a sequential manner. We mention that in other HSI movies in the LWIR data set, when the amount of chemical released was not as much as in the MeS cubes, the signal of this loop was not as strong. These experiments illustrate that the use of the Grassmann manifold together with PH provide insight into the presence and concentration of chemical contamination in a hyperspectral movie.

## Acknowledgments

This paper is based on research partially supported by the National Science Foundation grants DMS-1228308, DMS-1322508. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

## References

- [1] Bachmann, C.M., Ainsworth, T.L., Fusina, R.A.: Exploiting manifold geometry in hyperspectral imagery. Geoscience and Remote Sensing, IEEE Transactions on 43(3), 441–454 (2005)
- [2] Björck, Å., Golub, G.H.: Numerical methods for computing angles between linear subspaces. Mathematics of Computation 27(123), 579–594 (1973)
- [3] Carlsson, G.: Topology and data. Bulletin of the American Mathematical Society 46(2), 255–308 (2009)
- [4] Chang, J.M., Kirby, M., Kley, H., Peterson, C., Draper, B.A., Beveridge, J.: Recognition of digital images of the human face at ultra low resolution via illumination spaces. ACCV II, 733–743 (2007)
- [5] Chepushtanova, S., Gittins, C., Kirby, M.: Band selection in hyperspectral imagery using sparse support vector machines. In: Proc. SPIE. vol. 9088, pp. 90881F–90881F–15 (2014)
- [6] Chepushtanova, S., Kirby, M.: Classification of hyperspectral imagery on embedded Grassmannians. In: Proc. of the 2014 IEEE WHISPERS Workshop. Lausanne, Switzerland (June 2014)
- [7] Chepushtanova, S., Kirby, M., Peterson, C., Ziegelmeier, L.: An application of persistent homology on Grassmann manifolds for the detection of signals in hyperspectral imagery. In: Proceedings of the IEEE International Geoscience and Remote Sensing Symposium (IGARSS). Milan, Italy (July 2015)
- [8] Chung, M.K., Bubenik, P., Kim, P.T.: Persistence diagrams of cortical surface data. In: Information Processing in Medical Imaging. pp. 386–397. Springer (2009)
- [9] Conway, J.H., Hardin, R.H., Sloane, N.J.A.: Packing lines, planes, etc.: packings in Grassmannian spaces. Experimental Mathematics 5, 139–159 (1996)
- [10] Cosofret, B.R., Konno, D., Faghfouri, A., Kindle, H.S., Gittins, C.M., Finson, M.L., Janov, T.E., Levreault, M.J., Miyashiro, R.K., Marinelli, W.J.: Imaging sensor constellation for tomographic chemical cloud mapping. Applied Optics 48, 1837–1852 (2009)
- [11] Dabaghian, Y., Memoli, F., Frank, L., Carlsson, G.: A topological paradigm for hippocampal spatial map formation using persistent homology. PLoS computational biology 8(8), e1002581 (2012)
- [12] Draper, B., Kirby, M., Marks, J., Marrinan, T., Peterson, C.: A flag representation for finite collections of subspaces of mixed dimensions. Linear Algebra and its Applications 451, 15–32 (2014)
- [13] Edelman, A., Arias, T.A., Smith, S.T.: The geometry of algorithms with orthogonality constraints. SIAM J. Matrix Anal. Appl. 20(2), 303–353 (1998)
- [14] Edelsbrunner, H., Harer, J.: Persistent homology - a survey. Contemporary mathematics 453, 257–282 (2008)
- [15] Ghrist, R.: Barcodes: the persistent topology of data. Bulletin of the American Mathematical Society 45(1), 61–75 (2008)
- [16] Hatcher, A.: Algebraic topology. Cambridge University Press (2002)
- [17] Heath, K., Gelfand, N., Ovsjanikov, M., Aanjaneya, M., Guibas, L.J.: Image webs: Computing and exploiting connectivity in image collections. In: Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on. pp. 3432–3439. IEEE (2010)
- [18] Ma, L., Crawford, M.M., Tian, J.: Local manifold learning-based-nearest-neighbor for hyperspectral image classification. Geoscience and Remote Sensing, IEEE Transactions on 48(11), 4099–4109 (2010)
- [19] Manolakis, D.: Signal processing algorithms for hyperspectral remote sensing of chemical plumes. In: Acoustics, Speech and Signal Processing, 2008. ICASSP 2008. IEEE International Conference on. pp. 1857–1860 (March 2008)
- [20] Marrinan, T., Draper, B., Beveridge, J.R., Kirby, M., Peterson, C.: Finding the subspace mean or median to fit your need. In: Computer Vision and Pattern Recognition (CVPR), 2014 IEEE Conference on. pp. 1082–1089. IEEE (2014)
- [21] Perea, J.A., Harer, J.: Sliding windows and persistence: An application of topological methods to signal analysis. Foundations of Computational Mathematics pp. 1–40 (2013)
- [22] Singh, G., Memoli, F., Ishkhanov, T., Sapiro, G., Carlsson, G., Ringach, D.L.: Topological analysis of population activity in visual cortex. Journal of vision 8(8), 11 (2008)
- [23] Tausz, A., Vejdemo-Johansson, M., Adams, H.: JavaPlex: A research software package for persistent cohomology. In: Hong, H., Yap, C. (eds.) Proceedings of ICMS 2014. pp. 129–136. Lecture Notes in Computer Science 8592 (2014)
- [24] Topaz, C.M., Ziegelmeier, L., Halverson, T.: Topological data analysis of biological aggregation models. PLoS ONE 10(5), e0126383 (05 2015), http://dx.doi.org/10.1371%2Fjournal.pone.0126383
- [25] Wang, K., Thompson, J., Peterson, C., Kirby, M.: Identity maps and their extensions on parameter spaces: Applications to anomaly detection in video. In: Proceedings Science and Information Conference. pp. 345–351 (2015)