Template Matching via l 1 Minimization and Its Application to Hyperspectral Data

—Detecting and identifying targets or objects that are present in hyperspectral ground images are of great interest. Applications include land and environmental monitoring, mining, military, civil search-and-rescue operations, and so on. We propose and analyze an extremely simple and efﬁcient idea for template matching based on l 1 minimization. The designed algorithm can be applied in hyperspectral classiﬁcation and target detection. Synthetic image data and real hyperspectral image (HSI) data are used to assess the performance, with comparisons to other approaches, e.g. spectral angle map (SAM), adaptive coherence estimator (ACE), generalized-likelihood ratio test (GLRT) and matched ﬁlter. We demonstrate that this algorithm achieves excellent results with both high speed and accuracy by using Bregman iteration.


I. INTRODUCTION
I N recent decades, hyperspectral remote sensing technology has become highly developed.Remote sensing can be used to collect the reflective spectrum of the objects in a scene at specific wavelengths simultaneously, resulting in hundreds of digital images.The data collected from a hyperspectral sensor contains not only the visible spectrum, but also ultraviolet and infrared ranges as well.It is common to list the hyperspectral data in a three-dimensional array or "cube", with the first two dimensions corresponding to spatial dimensions and the third one corresponding to the spectrum.
Since the distinct materials leave different spectral signals, each of which is known as the spectral signature, there has been an increasing demand in using spectral imagery to detect specific objects or targets of interest.Template matching is one of the most important hyperspectral image (HSI) applications based on materials' spectral signatures.
In hyperspectral classification and especially target detection, the main task is to find the spatial pixels in threedimensional hyperspectral cube data for some given spectral signals of interest.However, this becomes difficult because of the uncertainty and variability of each material's spectral signature.The difficulties include the noise from atmospheric conditions, sensor influence, location, illumination and so on, all of which depend on when and where the image was taken.Another disadvantage is that the same materials will not generate identical spectral curves for physical reasons, especially between the laboratory data and actual image data.Occasionally, two different materials might unexpectedly reflect similar spectral signals, like two different types of vegetation.Additionally, the low resolution of the hyperspectral image make it reasonable to accept the fact that a single pixel may be composed of several different materials.These are usually called mixed pixels.The spectral signal of this mixed pixel is totally different from the spectrum of the pure materials.All of these facts increase the difficulties in processing HSI data, making classification and target detection highly challenging problems.
Previous work dealing with hyperspectral target detection includes both statistical and geometrical approaches.Statistical approaches are mostly based on the assumption that HSI data can be represented by a multivariate normal distribution.The statistical models employing two hypotheses H 0 and H 1 corresponding to given material absent and present are well known as the likelihood ratio (LR) test (optimum) and generalized-likelihood ratio test (GLRT) [1], [2] (adaptive) detectors.Some other related algorithms with details can be found in [3], [4], [5], [6], [7].Also, adaptive coherence estimator (ACE) [8] as an extension of generalized-likelihood ratio test (GLRT), employs a covariance matrix to identify the background, and achieves better result.However, in these models, prior knowledge including background information has to be obtained from training data, like the covariance matrix of the background.Where to set the threshold towards a reasonable compromise between probability of detection and probability of false alarms is also a question.In practice, there still exist serious challenges.
Many classical hyperspectral classifiers have been used in remote sensing, including minimum distance to mean (MD), spectral angle mapper (SAM), maximum likelihood (ML), Fisher linear likelihood and so on.In order to reduce the error, classification is usually performed in the Principal Component Analysis (PCA) and Generalized Eigenvalue (GEV) dimensionality reduced domains.
Recent methods have largely improved the efficiency of mixed pixel computation.The method presented in [16], [17], [18] using the l 1 unmixing model arose from the linear mixing model (LMM), which assumed that the observed pixel spectrum can be decomposed into a linear combination of the spectrum of the target and the background material.The distribution of each material decides the detection result.
Our method resembles SAM, but has an l 1 sparsity term which seems to improve the performance.We will discuss this further.The paper is organized as follows.A brief introduction of Bregman methods for l 1 minimization and the algorithm for template matching are presented in Section II.In Section III, we apply the proposed algorithm to a simulated data set and the real HSI image data sets.Comparisons with classical algorithms, SAM, ACE, GLRT and matched filter, are presented based on receiver operating characteristic (ROC) curves and detection results to show the obtained improvements.Finally, conclusions are given in Section IV.

II. TEMPLATE MATCHING APPROACH AND ALGORITHMS
In this section, we first briefly introduce the concept of Bregman iteration [19], and split Bregman [20], [21], then present our algorithm for template matching.This is done to give simple and effective methods for potential users to do the l 1 minimization.

A. Bregman Iteration and Algorithms
Bregman methods [22] applied to image processing [19] and compressive sensing (CS) [23], [24], [20], [25], [26] are an emerging field for solving sparse reconstruction and improving signal and image enhancement.Bregman iteration is based on the definition of Bregman distance.The Bregman distance for a convex function where p ∈ ∂J is a subgradient of J at the point v, and < •, • > is used as the inner product.In general Bregman distance shows the closeness between two points u and v, even though it doesn't satisfy the symmetry and triangle inequality properties of usual distance.
Consider a general constrained minimization problem: where H(u) is convex and differentiable energy function with zero as its minimum value.σ is the noise measure coming from the original problem.We transform it to an unconstrained minimization problem using: µ > 0 works as penalty weight.Bregman iteration solves a related problem iteratively by replacing J(u) by D J p k (u, u k ): where ∂H(u k+1 ) is a subgradient of H at u k+1 .The problems under consideration here use H(u) = 1 2 ||Au − f || 2 2 for some linear operator A and vector f .Convergence analysis of the iterative scheme is also provided in [19].It is also shown that under fairly weak conditions, H(u k ) → 0 as k → ∞.These can be easily solved by split Bregman.
Split Bregman was first introduced by Goldstein and Osher [20] for solving l 1 , TV, and related regularized problems and applied to various imaging problems [21], [26], [27].Consider of the problem where A ∈ R m×n , f ∈ R n be given, and m < n.The related unconstrained minimization problem is To solve this sparse reconstruction model, the iteration is generated by using an auxiliary variable d = u, and given by or where µ, λ > 0 work as penalties balancing the energy functions.Apply the Bregman formulation, we obtain the following algorithm (8) We"split" the l 1 and l 2 components of the first subproblem in (8) and perform optimization in an alternating fashion: (9) Set u 0 = b 0 = d 0 = 0, and define where the shrink + is designed for computing nonnegative solutions.
We have the iterations: Additionally, we can also enforce another Bregman iteration process, and obtain a two-layer inner-outer iteration style.The inner iteration is given as (9), using k as the inner iteration steps.And the outer iteration is given as , where u n+1 is generated by the inner iteration.n denotes the steps of outer iteration.
We know that ||u k − u k+1 || converges monotonically as k → ∞ in the inner iteration, so does ||u k − d k ||.It is also proved in [26] that the iteration converges even with only one inner step if 0 < λ < 1 ||A T A||2 .This is called Bregmanized operator splitting (BOS).In practice, we use a fixed number of inner iterations, usually 5-10 steps.µ is a parameter that balances the l 1 and l 2 components.It is determined by the actual data.This turns out to be a very efficient method for l 1 minimization, and we will use it below.We also found experimentally that λ = 100 ||A T A||2 is a good choice for the problem.
If A is fixed and saved as a prior knowledge, (λA T A+I) −1 in ( 12) could be precomputed, resulting in run-time savings.However, the size of square matrix (λA T A + I) is dependent on the number of columns of A. Note that if Au = f is underdetermined and the number of column is very large, it is impractical to compute the inverse.An alternative method could be used here.Instead of solving in ( 12), we insert two new variables y and g, Using (15), the computation of inverse matrix is only done for a much smaller size matrix of (λAA T + I).Equivalence between ( 14) and ( 15) can be proved easily through a few steps of simple linear transformation.

B. Template Matching Approach
The template matching approach basically employs the l 1 sparsity idea arising in compressed sensing [23], [28].Suppose the hyperspectral image is given as N by N pixels, each containing an M channel spectrum (a spectral band), as a three-dimensional cube.We aim to locate the pixels of interest in the image corresponding to a given spectral signal f , which is a specific object containing the same M channel spectrum.
We arrange A as an M by N 2 matrix, with generally M < N 2 .The signal a i,1≤i≤N 2 is the column of the spectrum matrix A, and corresponds to each pixel in the hyperspectral image.f is a given material's spectral signal.The aim is to find the nonzero component u i in the solution u ∈ R N 2 of the constrained minimization problem: The u i,1≤i≤N 2 corresponds to the same pixel as a i .By searching for the nonzero component u i , we can locate the pixels of interest.With the introduction of l 1 minimization, the algorithm becomes much more robust.This is different from the Basis Pursuit problem (4), which is used in compressed sensing.We are not sloving the exact solution of Au = f , but approaching the solution of min u ||Au − f || 2 2 , u ≥ 0 until the tolerance is reached.There are many methods proposed to do this numerical computation.In this paper, we apply split Bregman to the following unconstrained minimization problem.
To get a nonnegative solution, shrink + should be employed in the algorithm.In practice, the solution for hyperspectral data turns out to be nonnegative even without the nonnegative constraint employed here because of the nonnegativity of the data.
Surprisingly this simple idea works well because the l 1 norm term forces the energy to find a sparse solution u for µ large enough, which will enable the solution to pick up the matching pixels compared with the desired spectral signal f , and ignore nonmatching pixels.In detection applications, it is highly possible that the matching pixels are sparsely distributed, or the number of matching pixels is limited, which makes the problem difficult and interesting.The sparseness of the matching pixels means the sparseness of u.We say that given a proper value of µ, the solution results in u with a limited number of positive components, which just correspond to the matching pixels in the image with a given spectral signal f .However, even in the nonsparse (numerous) case, this also works well.We can easily illustrate this by a simple clean case of classification.We randomly choose n pixels in the image and assign them to be the material of interest.We assume the spectrum of these pixels identically equal to f , or proportional to f , which means kf , 0 < k < 1 works as a simulated shadow effect.After solving the l 1 minimization problem, we get a sparse solution u with each nonzero component u i corresponding to exact the same n pixel we chosen.The values of these u i are, close to 1 n after normalization.Of course, it does make sense since the algorithm aims to find the proper solution u to make ||Au − f || 2 as small as possible, while ignoring the spectral signatures distinct from f .
In the actual computation, the only significant parameter that should be considered is the quantity µ to control the sparsity of the solution, while λ could be chosen fixed.There exists a suitable range of µ for each set of HSI data.In detection applications, the number of matching pixels is much smaller than the total number of image pixels, µ could vary in a much large range.The result wouldn't be affected so much.But in the classification case, sometimes the number of matching pixels is numerous, the result depends largely on µ.If µ is chosen to be too large, the algorithm will only identify a few correct matches.To deal with this problem, we can remove the detected matches and do the algorithm again on the remaining pixels using a new matrix with fewer columns.
Simultaneously, if the background spectral signatures are quite distinct from f , we get no false alarms despite the change of µ.Otherwise, it will depend on µ.However, in practice the number of incorrect matches turns out to be much fewer than the number of correct matches, and the coefficients u i of most incorrect matches are much smaller than those of the correct matches, usually by 1 or 2 orders of magnitude.Thus we can set up a simple threshold to identify the correct result.
We notice the similarity of the template matching model proposed here and the unmixing model in [16], [17], [18].However, the connotation of parameters in the two models are totally different: A in the template matching model is listed as matrix of HSI data, while A in the unmixing model is arranged as endmember matrix; f in the template matching model denotes the spectrum of target, but in the unmixing model is the spectrum of one pixel in HSI data.In addition, the former model is typically underdetermined, but the later is generally overdetermined.Moreover, the indications of two solutions are not the same: one measures the similarities between spectrum, and the other illustrates the fractional values of endmembers occupying the pixel.Therefore, these are two different methods and two different applications.
In contrast to the other methods, the template matching algorithm involves neither any requirement of background information nor the assumption of multivariate normal distribution of the background and designed material during the whole computation process.Only the spectrum of material of interest is needed in the algorithm.The method can deal with the case when the matches are sparse or numerous and when background models are unavailable and unreliable.
Moreover, the template matching algorithm only matches the spectrum of image pixels with the spectrum of material interested, and doesn't involve spatial information.So the method is parallelizable and can work locally under the condition of inhomogeneous background.Thus the l 1 minimization based algorithm is realistically applicable.

III. NUMERICAL EXPERIMENTS AND RESULTS
In the following section, the split Bregman algorithm indicated in Section II is implemented in Matlab to assess this applicability of l 1 minimization for hyperspectral template matching.Simulated data and real hyperspectral images data will be tested.
The code was written in MATLAB using 3 GHz, Intel Core 2 Duo CPU.The most time-consuming part is a matrix inverse computation.But the matrix A is dependent on the HSI data, which could be used as prior knowledge.Depending on the size of data to be tested, one can choose to use the alternative method (15).Furthermore, because the problem is parallelizable, one can also split the image into smaller patches, and the original A into smaller size matrices.
Concerning the computational complexity, we note that, given A ∈ R m,n , m < n.The complexity of the algorithm per iteration using the alternative method ( 15) is thus O(m 2 n).

A. Simulated Data Set
The simulated data set we tested here is based on the a HYDICE hyperspectral image called Urban [29] which is publicly available.Firstly we performed a data correction by deleting the water absorption noise bands from the original 210 bands, leaving a 163-band image for actual computation.We chose a 50 by 50 pixels patch image in Urban to perform the background, basically including dirt road, tree and grass materials, as shown in the left of Fig 1 .So the size of matrix A to be computed was 163 by 2500, and the solution u ∈ R 2500 .
We manually chose another material to perform the target or material of interest, and here we used a material of metal rooftop appeared in another part of Urban to be distinct from the background, as shown in Fig 1 .Ten pixels in the patch were randomly selected as the locations of matching pixels, and were assigned to be metal rooftop.After that, we finished the test data simulation.
Solving the constrained minimization problem via split Bregman algorithm generates a sparse signal u.This enables us to locate the matched pixels by just checking the nonzero components u i .In the following, we list the different kinds of examples we tested, and show the performance of the algorithm under different cases.The result is also compared with other two methods.We used fixed value of parameter λ and varying µ in the experiments.In practice, we found choosing λ larger will speed up the algorithm, here, we use λ = 100 ||A T A||2 , where works well even without the proof of convergence.
Besides the time used to compute the matrix inversion, it took only 1-2 seconds to compute the sparse u of the simulated data to determine the targets using split Bregman.
1) Robustness to Noise: To simulate the effect of the environmental variations, illumination and so on, we added noise n ∈ N (0, σ) to the spectrum of 10 simulated pixels to make this simulation more physically meaningful.In order to characterize the noise level, we use SNR (signal to noise ratio) defined as where u is the mean of the signal, and σ is the standard deviation of the noise.With µ = 0.01, and SNR=20.3,we got the the solution u with the exact 10 nonzero components u i to be [0.1029,0.0832, 0.1107, 0.1119, 0.0906, 0.1143, 0.0759, 0.1294, 0.1065, 0.0849], the other u i = 0, and the l 1 norm of u = 1.0101.And each of the nonzero u i corresponded to one correct simulated pixel.
As the noise increased, we can still locate the correct pixels, but a few incorrect positive components u i may arise.We introduce some measure quantification to evaluate the performance of algorithm, true positive rate (TPR) and false positive rate (FPR), which are defined as where T P is the pixels detected correctly (true positive), F P is the false alarm pixels (false positive), T N is background pixels with zero components (true negative) and F N is the matching pixels with zero components (false negative).In the test SNR=10, we obtained the average T P R = 98.6% and F P R = 0.004% based on 100 computations.Increasing the noise to SNR=5, we found the algorithm gave the average T P R = 98.1% and F P R = 0.87%.The number of these incorrect pixels detected is limited since they just work to reduce the residual.Moreover, the magnitudes of the incorrect positive components is relatively small.A simple threshold can be employed to reduce the number of F P to zero.
In addition, we can also add a small perturbation to f , or change the magnitude of f .The algorithm still works well for moderate changes.
2) The Choice of µ: µ is a parameter that balances the sparsity of the solution and how closely Au describes the spectrum signal f.A good choice for µ is determined by the actual data.If the number of matching pixels is very large or the background materials have similar spectrum with f , variation of µ will affect the result.Larger values of µ may only detect a subset of of these correct pixels.This encourages us to process the algorithm in the following iterative way: i Choose µ large, get the solution u and identify the correct pixels.ii Remove correct pixels identified in i from A, get a matrix A 2 with fewer columns.Repeat i until you get as many good pixels as possible.
Experimentally, in a few iterations (2 or 3), we obtain all the correct pixels.The last iteration usually contains both correct and incorrect pixels.
On the other hand, smaller µ can solve the case that matching pixels are relatively numerous in the image.Under that situation, we are not supposed to look for a sparse solution, but a quite non-sparse solution.
In the follow contrast examples, we first chose 10 pixels out of 2500 to be the assigned matching pixels, with SNR=20.3.The result of computation found all the 10 correct pixels with 0 false alarms, even when µ changed from 1e-2 to 1e-6.Secondly we chose 2000 pixels out of 2500 to be the assigned matching pixels, with SNR=20.3.The difference was when µ changed from 1e-2 to 1e-4, the number of nonzero u i changed from 100 to 2003 (also containing a few incorrect pixels), as shown in Table I.The different results shows the power of µ in controlling the sparsity of solution.
3) Mixed Pixels Test: In many situations, the designed material only partially occupies the mixed pixels.We generated synthetic mixed pixels by adding one or two background materials proportionately to the pure pixels according to a linear mixing model (LMM).We tested the algorithm under different proportions.
We have observed that the algorithm can obtain the most correct matching pixels if the amount of pure material's  contribution to the mixed pixels is higher than the others, ranging from 95% to 50%.The result is shown in Table II (using SNR=20.3 and µ = 1e − 2 ).

4) Comparisons:
The spectral angle map (SAM) involves computing the normalized inner product the cosine of the angle between two vectors, the spectrum of target f and pixels a i , which is always positive because spectral vectors contain only nonnegative components.The adaptive coherence estimator (ACE) algorithm uses a covariance matrix to qualify the background.It computes a ratio at each pixel where f is the target of interest, a i is a pixel in the test image, Σ is the covariance matrix of background with zero mean.Thresholding is done after T (a i ) is generated for each pixel.Then SAM and ACE denote those pixels as positive for which the value of T (a i ) exceeds a threshold τ , a training data, and denote the others as negative.The performances of SAM and ACE are highly dependent on the training data.  in subimage of SNR=5 are all located on the y axis, which are probably invisible.We can see with small enough noise, all three algorithms can detect the correct matching pixels well, while ACE needs a good threshold τ to decrease the false alarms.Increasing the noise to SNR=5, we found that l 1 template matching essentially shows zero probability of false positives and worked better than SAM and ACE, which intended to obtain low detection rate when the noise increased.

B. HSI Data Set
1) Fake Leaves: We performed l 1 template matching on the Fake Leaves indoor image provided by Surface Optics Company [30] using the SOC-700 hyperspectral imaging system as shown in Fig 3 .The HSI data includes 640x640 pixels and 120 bands in the 400-900 nm range.Under RGB model, fake leaves look exactly same as true leaves, however, true leaves have a strong feature in the near infrared.
To perform the algorithm, we selected 15 patches that contain most of the objects present in the image for simplicity.This is shown in Fig 4 .Using the template matching algorithm, we detected almost all the correct pixels.We used a red color to show the pixels detected to be fake leaf.As shown in the image, we detected 678 correct and missed around 16 fake leaf pixels among 6000 pixels, the true positive rate equals was 97.69%.All of the missed pixels were located between the fake leaf and true leaves.
2) NGA HSI Data Set: We worked with the hyperspectral image data provided by NGA, as an experimental data for classification given the materials' signatures.It contains 3349   pixels with 106 bands.A total of 9 endmembers are extracted from the data, including Coniferous, Deciduous, Grass, Lake1, Lake2, Crop, Road (Asphalt), Concrete, Gravel.To make matters worse, most of them look quite similar to each other as seen in Fig 5 .We performed l 1 minimization and split Bregman to all 9 endmember targets and achieved satisfactory results, which were shown in Table III.
To get a better result, we used µ = 1e − 4, λ = 100 ||A T A|| 2 .Based on the evidence that the coefficients of incorrect pixels were relatively small, a simple threshold can be set up to decrease the number of incorrect pixels as in Table IV.In order to see the effect of the algorithm, we compared the result of the spectral angle map (SAM) and adaptive coherence estimator (ACE) using ROC curves while changing the threshold τ .for target detection scoring [35], which is available at http://dirs.cis.rit.edu/blindtest/.The self test data used in the experiment was collected in the small town of Cooke City, Montana, USA as in Fig 9 .Multiple vehicles and fabric panels were provided as the target.
We choose a 60x60 spatial ground truth with 126 spectral bands, and material "3m Red Cotton" as the target to evaluate the results between five algorithm: l 1 template matching, spectral angle map (SAM), adaptive coherence estimator (ACE), generalized-likelihood ratio test (GLRT) and matched filter.The parameters of all the algorithms have been tuned to be the best.IV.CONCLUSION In this paper, we propose an l 1 minimization based approach for template matching and apply it to hyperspectral classification and target detection.This algorithm is simple to understand and implement, and only involves a few lines of code.Numerical results show that the algorithm is efficient in finding the designed pixels in a hyperspectral image data for a specific material, without using any background information.There are many other applications of this algorithm worth exploring.We will extend the applicability of this algorithm in the future.The limitation of the algorithm is that it requires the materials' spectral signatures to do the computation which is dependent on the existence of endmember extraction methods.

Fig. 1 .
Fig. 1.Simulated data image.Left: RGB view of patch chosen in the Urban.Right: Spectrum of materials in the patch.Red: target-metal rooftop.Yellow: dirt road.Blue: trees.Green: grass.

Fig. 2 .
Fig. 2. ROC Curves of SAM in blue and ACE in green vs l 1 template matching in red for simulated dataset.From left to right, SNR=20, SNR=5.

Fig. 4 .
Fig. 4. Fake leave detection result.Left: Patches chosen as test pixels.Right: Red pixels are detected to be fake leaves.

Fig. 6 .
Fig. 6.Smith Island: Up: RGB view of Smith Island data.Low: Ground truth of 22 materials.
Fig 8 shows the comparison between the three experiment's result for all 22 materials.As shown in Fig 8, almost all the red lines stay to the upper left position of green lines and blue lines, which means the experiments get better result than SAM and ACE.The reason for this is because our algorithm is much more robust when the noise increases.4) RIT self test Data and Comparisons: The data set selected for this comparative evaluation is one that is used
Fig 10 shows l 1 template matching outperforms the other matching algorithms.

TABLE III CLASSIFICATION
RESULTS FOR HSI DATA WITHOUT THRESHOLDING