A structure tensor for hyperspectral images

A Structure Tensor for Hyperspectral Images Maider Marin-McGee and Miguel Velez-Reyes Laboratory for Applied Remote Sensing and Image Processing, University of Puerto Rico-Mayaguez, PR, USA E-mail: {maider.marin, miguel.velez-reyes}@upr.edu ABSTRACT In this article, a structure tensor for hyperspectral images (HSI) is proposed. A weighted zero mean smoothed gradient to calculate the initial matrix field is used. The weights are constructed from the zero mean data by comparison with a normalized absolute value of the median. The problem with the classical definition is the assumption that all bands provide the same amount of edge information giving them the same weights. As a result many non-edge pixels are reinforced producing false edges. Therefore, other processes that depend on the structure tensor will be misguided. The proposed weights are selected from pixel’s spectral values that are greater than and closely around the absolute value of the median, reinforcing only the better candidates in the spectra to be edges in their respective spectral interval, therefore making the structure tensor a better edge discriminator. Comparisons of this method with the standard definition of the structure tensor and results of using the proposed method for nonlinear tensor anisotropic diffusion are presented. Index Terms— Hyperspectal Images, Structure Tensor 1. INTRODUCTION A structure tensor for vector valued images such as hyperspectral and multispectral images was proposed in [1]. The idea was to have one edge descriptor for all bands such that edge information along a given direction in one channel reinforces edge evidence in the other channels [1]. Therefore, any process using the edges found by the structure tensor will behave in the same way on each channel and have the information from all channels. This is in contrast with processing each band independently that can produce edges in different locations for different bands, generating edge dislocation and unwanted discontinuities [2]. The structure tensor is derived from the gradient. It provides the main directions of the gradient in a specified neighborhood of a point. Its eigen-decomposition provides information on how strongly the gradient is biased towards a particular direction which is known as coherence. Therefore, it can be used for both orientation estimation and image structure analysis [3]. It has proven its usefulness in many application fields such as corner detection [4], texture

трс΀к΀норр΀ллйк΀кͿккͿглпͨйй
Ылйкк
ΒΎΎΎ

analysis [5-7], diffusion filtering [8], [9] and optic flow estimation [7], [10], [11]. In the classical definition of the structure tensor, the fusion of all edge information in one band is done with a pixel-based averaging along the spectral values of its initial matrix field (see Section 2). It is based on the assumption that all bands provide the same amount of edge information and hence, all of them have the same weights. However, in a HSI depending on the spectrum and materials that compose the object, some bands may show all or part of an object while some of them may show no evidence of its existence. Correspondingly, some bands of the structure tensor will show a given edge and others will not. Therefore the fusion step should weigh the spectra depending on the amount of edge information they provide. The proposed method uses a pixel-based weighted zero mean gradient in which the weights are constructed from the data by comparison with a normalized median. This method is simple to understand, depends on the data, has a linear computational time with respect to number of pixels, and since the operations are done pixel-wise, it is totally parallelizable. Experimental results show that nonlinear tensor anisotropic diffusion is greatly improved by using the proposed method while also reducing the number of iterations need it to obtain a result. This paper is organized as follows: Section 2 is an overview of the structure tensor. Section 3 presents the proposed method and Section 4 shows the effectiveness of this method via practical examples with comparison to the standard approach. Finally, Section 5 presents the conclusions. 2. THE STRUCTURE TENSOR The local structure at a pixel can only be described considering the data of its neighborhood [12]. This fact introduces the concept of the structure tensor. Given a pixel, the task of knowing if it is part of an edge or a corner cannot be accomplished with only its gradient information. The incorporation of its neighbors’ gradient information is required. The local neighborhood information becomes more important when the data is corrupted by noise or other unwanted signals, since the structure has to be estimated in the presence of unreliable data. The edge descriptor known as the smoothed gradient I can be a candidate. This descriptor is the gradient of the smoothed image

corner regions) [12]. Figure 1(d) shows the coherence orientation for the fingerprint image and also how well the fingerprint singularity (minutiae) is described. Jρ = J0 with ρ=0 is the structure tensor used for edge enhancing diffusion. Jρ with ρ > 0 is the structure tensor used for coherence enhancing diffusion, optical flow, etc. [8]. 3. PROPOSED METHOD

Figure 1. Smoothed Gradient vs. Structure Tensor in a fingerprint image of size 227× 227. (a) Original image. (b) Gradient orientation σ = 0.5. (c) Gradient orientation σ = 2.5. (d) Structure Tensor σ = 0.5, ρ = 4.

I  G * I , where the * is a convolution and G is a Gaussian of variance   0 . However, the smoothed gradient can suffer from cancellation effects as shown in Figure 1 (b)-(c) which illustrate gradient orientation using grey values with vertical gradients depicted in black and horizontal ones in white. If σ is small as in (b), then high fluctuation of noise remains; as σ gets larger as in (c), then it is useless, since neighboring gradients will have the same orientation but opposite signs (direction) thus cancelling out one another. To avoid that effect, the outer product of the smoothed gradient, which is a symmetric, positive semidefinite matrix, is considered in the structure tensor, and known as the initial matrix field [12]

 I2 J 0  I I T   x

Ix I y

Ix I y  , I y2 

(1)

subscripts denote partial derivatives and Iσ=[Ix Iy]. For the vector-valued image I = (I1,…,Im) with m bands, the structural information from all channels is coupled by taking the average over all matrices, in this case the initial matrix field is given by

 m 2  I x ,i m i 1 J 0  J 0i   m  i 1  I x ,i I y ,i

i 1

 I y ,i i 1 . m 2 I y ,i i 1  m

I

x ,i

(2)

The division by m has been dropped for simplicity. The structure tensor for a neighborhood of scale ρ is computed by Gaussian convolution of the components of J0 with a Gaussian kernel Gρ in any of the cases (3) J  J 0 * G This definition has two advantages: (i) Robustness under noise as result from smoothing the resulting matrix field; (ii) Additional information is obtained by integrating local orientation so it becomes possible to distinguish areas where structures are oriented uniformly (like in regions with edges) from areas where structures have different orientations (like

The initial matrix field in Equation (2) describes the fusion step in which our method is concentrated. For symmetrically distributed data, the mean is a good measure but as in the case of hyperspectral data this assumption is not necessarily true for the majority of the pixels. On the other hand, the median can handle both cases, and it gives a better idea of any general tendency in the data, since 50% of values are above it, and 50% below it no matter how the data is distributed. We propose to use a weighted zero mean gradient with components defined by Ix = Ix ‒ I̅ x and Iy = Iy ‒ I̅ y with the upper bar denoting the mean, and the under bar zero mean variable, defined as:  wx 0   I x  wI       wx I x wy I y  ,

0 wy  I y  to calculate the initial matrix field and

 m 2 2  wx ,i I x ,i m i 1 J 0  wi J 0i wi T   m  i 1  wxy ,i I x ,i I y ,i

i 1

 I x ,i I y ,i i 1 , (4) m 2 w2y ,i I y ,i i 1  m

w

xy ,i

with the weights wx, wy calculated using the zero mean gradient components Ix̱, Iy respectively and wxy= wxwy. Clearly J0, as defined in Equation (4), is still a symmetric, positive semi-definite matrix. The weight for a pixel p in a gradient component, denoted by wp, is calculated as follows:  p median( p )  1  Step 1. w p   (5) p  median( p ) median( p )  1  In Equation (5), a threshold of 1 was chosen to avoid division by small numbers. Next wp was normalized  wp max w p  1   (6) Step 2. w p   wp max w p  1  max w p  

 

 

 

  0 wp ,i  a Step 3. wp ,i   for a  0. (7)   wp ,i wp ,i  a Absolute values were used since the gradient has negative values. After Step 2, the range of values of wp includes negative values such that min(gradient component) ≤

min(wp) < 0 and max(wp)=1. The threshold a is chosen by finding the histogram of wp and clipping its minimum value at the first bin in the tail which height is  1% of the total of pixels. 4. NONLINEAR TENSOR ANISOTROPIC DIFFUSION In this section, the nonlinear tensor anisotropic diffusion (NTAD) is described by the partial differential equation: u (8)    Du   0, in Ω t where u, the filtered version of image f (x) with domain Ω, is the solution of the initial boundary value problem for the diffusion equation with f as the initial condition with the Neumann boundary condition, i.e. flux is zero outside the boundary Ω. D is a positive semi definite symmetric matrix known as diffusion tensor defined as: 0 v  D   v w  1 (9)  .

0  2  w [v w] are the orthogonal eigenvectors of the structure tensor, Jρ, defined in Equation (3), and κ1 and κ2 are functions of its respective eigenvalues µ1, µ2. Edge enhancing diffusion uses the structure tensor, J0, defined in Equation (2). In the experimental results, a finite volume discretization [13] was used and eigenvalues defined by [8] 1 1  0    3.31488 1   1  exp !! (  /  ) 4 "" 1  0 # 1 $   2  1 / 3 5. EXPERIMENTS The initial matrix field of a region of Forrest Radiance I (FRI) run 5 is calculated for the standard formulation and using proposed method. FRI is then processed with two iterations of NTAD. 5.2. Forrest Radiance I A portion of FRI run 5 is used. This image was collected by the airborne imaging spectrometer HYDICE of a forest region at the U.S. Army Aberdeen Proving Ground in 1995 by the HYMSMO (Hyperspectral MASINT Support to Military Operations) program. It contains 120×120 pixels with 145 spectral bands in the 400-2500 nm range. Bands with low signal to noise ratio were rejected as in [14]. This region of interest was used because targets are smaller and discriminating edges are difficult because the image has a strong vertical texture. Figure 2 Shows initial matrix field components after two iterations of NTAD with σ = 0.5, step size = 0.5,  = 0.3. All the images where scaled to 1 for qualitative comparison. The left column of Figure 2-shows the initial matrix field components related to Ix2, Ixy, Iy2

Figure 2. Initial matrix field components Ix2, Ixy, Iy2 calculated using the standard formulation (left column) vs. proposed method (right column). Results after two iterations of NTAD. FIRST ROW: average of Ix2 bands. MIDDLE ROW: average of Ixy bands. LAST ROW: average of Iy2 bands. (Zoom in to see details).

calculated using Equation (2). The right column of Figure 2 shows the same components calculated using Equation (4) with weights defined in Equation (7) and a = −2.5. Results of the standard method for components Ix2, Iy2, enhance nonedge pixels so that false edges appear in unwanted signals and actual edges are hardly distinguishable. The proposed approach for Ix2 produces better results than the standard method. Target edges can also be distinguished and the vertical texture is less preserved. The results for Iy2 are much better than the standard method. Targets, including the smallest size target, can be completely distinguished and unwanted signals decrease dramatically. For component Ixy, the contrast is enhanced and all the targets can be distinguished since the weighting process rescales this component from negative values where the standard processing technique produces positive values.

This definition is based on a weighted zero mean smoothed gradient with weights that depend on a normalized median. Results in real HSI showed better discrimination of edges in all components of the initial matrix field. Results show that fewer unwanted signals are enhanced and more actual edges are enhanced compared with the standard method. The proposed structure tensor was used to improve NTAD producing a strong diffusion around edges and preserving the majority of them. It was also shown that this process accelerates the diffusion and the weighting process produces more uniform diffusion . 7. ACKNOWLEDGMENTS This work is supported by the U.S. Dept of Homeland Security under Coop. Agreement 2008-ST-061-ML0002. 8. REFERENCES

Figure 3. Results of nonlinear tensor anisotropic diffusion (NTAD) applied to FRI-run 5. RGB composite with bands 49-3718. (a) Original image (b) NTAD using standard definition. (c): NTAD using zero mean gradient without weights. (d) NTAD using weighted zero mean.

Figure 3 shows results for NTAD after two iterations with the same parameters as Figure 2. Figure 3(b) shows NTAD using the classical structure tensor defined in Equation (2) while Figure 3(d) shows NTAD using the weighted zero mean method, see Equations (4)-(7). Clearly both methods preserve the edges of almost all targets. The proposed algorithm produce a stronger diffusion to the nontarget texture but it smooths the smallest target. It will take approximately 5 more iterations for the classical approach to obtain a result that are visually comparable to our proposed algorithm. ∝=0.3 was chosen to fit the classical approach however it is not necessarily the best value for the proposed method. Figure 3(c) results from NTAD using a zero mean smoothed gradient but skipping the weighting process described in Equations (5)-(7). Figure 3(c) was done to check if there was any improvement using the weighting process. Comparing Figures 3(c) and (d) we can see that, under the same circumstances, the weighting process produces a more uniform smoothing and covers more regions than using a non-weighted zero mean smoothed gradient. 6. CONCLUSIONS A definition of the initial matrix field for HSI has been proposed to incorporate edge information from all bands and obtain a better edge discriminating structure tensor for HSI.

[1] S. Di Zenzo (1986), “A note on the gradient of a multiimage,” Computer Vision, Graphics, and Image Processing, vol. 33, no. 1, pp. 116 - 125. [2] J. Weickert (1994), Scale-Space Properties of Nonlinear Diffusion Filtering with a Diffusion Tensor. University of Kaiserslautern, P.O. [3] T. Brox and J. Weickert (2002), “Nonlinear Matrix Diffusion for Optic Flow Estimation,” in Proc 24th DAGM-Symposium on Pattern Recognition, London, UK, 446–453. [4] W. Fröstner and E. Gülch (1987), “A Fast Operator for Detection and Precise Location of Distinct Points, Corners and Centers of Circular Features,” in Proc. ISPRS Intercommission Workshop on Fast Processing of Photogrammetric Data, 281-305. [5] T. Lindeberg (1993), Scale-Space Theory in Computer Vision, Springer. [6] A. R. Rao and B. G. Schunck (1989), “Computing oriented texture fields,” in Proc. CVPR '89, San Diego, CA, USA, 61-68. [7] J. Bigün, G. H. Granlund, and J. Wiklund (1991), “Multidimensional Orientation Estimation with Applications to Texture Analysis and Optical Flow,” IEEE Trans. Pattern Anal. Mach. Intell. , 13, 775–790. [8] J. Weickert (1998), Anisotropic Diffusion in Image Processing. Teubner-Verlag. [9] D. Tschumperlé (2006), “Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE’s,” IJCV, 68, 65-82. [10] H. Scharr and H. Spies (2005), “Accurate optical flow in noisy image sequences using flow adapted anisotropic diffusion,” Signal Processing: Image Communication, 20(6), 537-553. [11] T. Preusser and M. Rumpf (1999), “Anisotropic nonlinear diffusion in flow visualization,” in Proc Visualization,. 325 -539. [12] T. Brox et al. (2006), “Adaptive Structure Tensors and their Applications,” in Visualization and Processing of Tensor Fields, J. Weickert and H. Hagen, Eds. Springer Berlin Heidelberg, 17-47. [13] O. Drblíková and K. Mikula (2007), “Convergence Analysis of Finite Volume Scheme for Nonlinear Tensor Anisotropic Diffusion in Image Processing,” SIAM J. Numer. Anal., 46 (1), 37– 60. [14] J. Kerekes, A. Cisz, and R. Simmons (2005), “A comparative evaluation of spectral quality metrics for hyperspectral imagery,” Proc. of SPIE, 5806, 469–480.