A Wavelet Multiresolution Technique for Texture Analysis and by yurtgc548


									  A Wavelet Multiresolution Technique for Texture Analysis and Segmentation of
                                 SAR Images

                          G. De Grandi (1), J. Fortuny(1), J.S. Lee(2) , D.L. Schuler(2)
                                     European Commission – Joint Research Centre
                                              21027 Ispra (VA), Italy
                                             Naval Research Laboratory
                                         Remote Sensing Division Code 7263
                                         Washington, DC 20375-5351, USA

A technique for multi-scale texture analysis and segmentation of SAR images is presented. Textural features are
extracted using a decomposition based on a wavelet frame. The feature vector is composed of local estimates of
normalized variances of the smooth image and of the wavelet coefficients. The technique is also extended to the
polarimetric case. To the purpose polarimetric synthesis is performed in the wavelet domain. In the synthesis,
selected antenna configurations that are known to be optimal with respect to texture diversity are used to construct the
feature vector. Controlled experiments, based on Monte Carlo simulations, are set up to assess the performance of the
technique. Finally, examples are given of the application of the technique to real SAR data (ERS-1 and SIR-C)
related to tropical swamp forests mapping and oil slick detection.

Keywords: texture, segmentation, wavelet frame, multi-scale representation, SAR, polarimetry.


The notion of texture is well defined at the level of our visual perception. The human eyes and brain are
able to delineate areas that exhibit common spatial patterns. However this notion which is fairly simple
fro m the intuitive point of view, eluded until now all efforts to cast it into a universal mathematical
definition. Despite this, several approaches in image processing have been devised that work reasonably
well in certain specific areas of application.
Narrowing our scope to hit a concrete target, we concern ourselves here with textures that are somehow
associated with the “roughness” of the – eventually multi-dimensional - surface that describes the
backscatter as a function of space in a radar image. From the statistical point of view this surface is
considered a random process in two dimensions and analyzed using local estimates of a second-order one-
point statistics (variance) and of a second-order two-points statistics, the structure function (field of
increments of a random process over a distance). The roughness associated with a feature of interest (e.g.
a natural target such as a forest) develops at some specific scale multiple of the resolution element, a fact
that suggests that measures must be taken at several scales to characterize and discriminate different
targets. Questions are then: how do we measure such textures, how do we take scale into account, how
do we detect structures or segment homogeneous regions provided by these measures? Additional
questions related in particular to Synthetic Aperture Radar (SAR) imagery are: how do we take speckle
into account and what is the role of polarimetry in texture segmentation? We will try here to give one
possible answer to those questions.
2.1. Texture Measures
Our approach is heavily inspired by by a texture segmentation and feature reduction technique presented
in [1][2]. The underpinning mechanism for the extraction of texture measures is a multi-resolution image
representation provided by a wavelet frame that acts as a differential operator [3]. Text ural features are
captured by spatially local estimates of: i) variance of the smooth signal, ii) variance of the wavelet
coefficients. Intuitively, these measures are sensitive to sharp transitions (contours of areas with different
reflectivity, point targets), and to edge density, smoothness and swing (at a given scale) within extended
In more detail a feature vector holding the texture measures is constructed as follows. We will consider at
first the case of scalar radar.
The wavelet frame decomposition is imp lemented numerically using an “á trous” algorithm and a filter
bank structure. The low-pass and high-pass filter coefficients are given in [3]. Let S j be the image at the

output of the low-pass filter (smooth image),      Wjx and Wjy the wavelet coefficients in the image row and
column directions (output of high pass filters) at scale s = 2 . We assume the convention that the j index

increases with scale. The following steps are performed: i) normalization of the wavelet coefficients by
the smooth signal at the same scale to cancel speckle influence; ii) introduction of non-linearity; iii)
smoothing by convolution with a separable discrete cubic B-spline kernel enlarged by a factor m . The
support of the discrete cubic B-spline is 4. Therefore averaging will take place in a neighborhood of
 m × 4 pixels.
The smoothing window is fixed, and therefore in this respect the approach is not multi-scale. However,
within the averaging window measures of textures that occur at different dyadic scales are taken into
account by the frame decomposition. In this sense the approach is multi-scale.
In summary, for each of the two scales in the decomposition the feature vector will be composed of:

                               (S )
                                                      Wjx   
                                                                                 W jy   

                       v1j =                v2 j   =                   v3 j =                          (1)
                                        2             Sj                       Sj     
                                S p−1
                                  j                                                   
The same pixel size (equal to the original image pixel size) is kept for all scales and features, resulting in
2.2 Feature Reduction
Feature reduction is obtained by an approximate solution of the multiple discriminant transform (MDT)
as proposed in [2]. This estimate of the MDT solution does not require a-priori knowledge about the
regions (e.g. the within-region and between region scatter matrices) and relies on the availability of local
variance estimates at two scales. On the other hand, for this reason, the reduced features are only
available at scale 2m. For trials where the characteristic scale of texture is too small to be handled at
resolution 2m using the approximate MDT, a Karhunen-Loève (KLT) transform is used instead.
2.3 Extension to the Polarimetric Case
We are also interested in dealing with texture that arises from a mixture of polarimetrically different
scattering mechanisms [4]. The polarimetric information is included by augmenting the feature vector
with local estimates derived from the wavelet decomposition of the backscattered intensity synthesized at
a number of antenna configurations. An analytical tool presented in [4] - the polarimetric texture
signature – indicates that optimal polarization states for enhancing differences in polarimetric text ure are
the linear cross-polarized states with orientation angle ψ = 0 (HV) and ψ = 45 . The latter
                                                                             o                   o

configuration will henceforth be referred to as xpol45. Accordingly the feature vector is constructed by
stacking the local estimates at two scales as in (1) derived from the synthesized power at the two optimal
polarization states.
We have also experimented with a second alternative, where texture measures derived from the Cloude-
Pottier scattering mechanisms decomposition - in particular entropy - are used as entries in the feature
The frame decomposition of the synthesized power images at different polarization states is performed
using the following technique. The polarimetric power synthesis and the wavelet fra me operators
            I = T(ψ T , χ T , ψ R , χ R ) ⋅ M PW ⋅ I = PW ⋅ T(ψ , χ) ⋅ M = T( ψ, χ) ⋅ PW ⋅ M           (2)

where I is the synthesized backscattered power,       PW is the frame operator, M is the Stokes operator and
 T is the linear transformation holding information on the transmit and receive antenna configurations
(Stokes vectors). Therefore the wavelet transform of the backscattered intensity at any polarization state
can be derived from the wavelet transform of the Stokes operator, or, in other words, polarimetric
synthesis is performed in the wavelet domain. The advantage of this approach is that the wavelet
transform, which is the computationally heavy part, needs to be performed only once. Once the wavelet
coefficients for the Stokes operator have been computed then the coefficients of the backscattered
intensity for all desired polarization states can be derived by applying the synthesis operator, which is
computationally quite lighter.
A number of controlled experiments based on simulated polarimetric SAR imagery were performed.
These experiments allow us to explore the potential of the technique with respect to different detection
and segmentation problems and to assess the performance of alternatives in the processing method. In this
short communication, we report next about one case concerning the simulation of a fragmented forest.
3.1 Fragmented Forest
In this experiment we simulate the polarimetric scattering from a natural extended target – a fragmented
forest. In this case we assume that spatial modulation of the radar return is caused by openings in a dense
canopy where the electromagnetic wave can penetrate down to a bushy understory. The model consists
therefore of two vegetation classes – referenced as forest and understory - with different structural
parameters. Scattering at P-band for each class is characterized by a polarimetric covariance matrix that
is computed using the University of Texas at Arlington (UTA) wave scattering model. Details on the
model parameters can be found in [4]. We expect to generate polarimetric texture by a mixture of the
returns from the two classes because different scattering mechanisms come into play. For instance, for the
forest class the VV backscatter is the result of both direct scattering from the canopy and indirect
scattering branch-soil, while at HH polarization the dominating mechanism is indirect branch-soil
scattering. This creates polarimertric diversity between the two channels. For the understory class, direct
and indirect scattering from the branches and direct scattering from the surface affect the HH and VV
channel equally.

  a)                                       b)                                    c)

   d)                                      e)

Fig. 1. Simulated fragmented forest data set. Mask used to define regions holding a radiometrically “pure” class (closed
canopy forest) and a mixture of two classes (fragmented forest) (a). HV amplitude image (b). Feature vector constructed from
the synthesized xpol45 image: Variance of smooth image at scale 22 (c). Textural segmentation using K -means clustering
with 5 categories. Results when clustering is applied to the first KLT transform in the AMDT procedure: Clusters belonging
to high textures related to fragmented forest (grey to black) (d). Clusters belonging to low textures related to closed-canopy
forest (grey to white) (e). Segmentation when clustering is applied using the full AMDT procedure (f).

A simulated one-look scattering matrix image for this trial is constructed in the following way. Random
deviates of scattering matrix samples, corresponding to the covariance matrices computed by the UTA
model for the two vegetation classes, are generated by a Monte Carlo technique. A square array of these
samples is partitioned into two sub-sets by a curvilinear mask. The first subset is filled with scattering
matrix samples belonging to a pure class (either forest or understory). The second set is filled with a
mixture of samples belonging to the two classes. The mixture is laid out spatially simply alternating
pixels of each class. For the sake of brevity we report here only the case where the pure class is forest.
Thus the segmentation problem concerns the ability to discriminate pristine forest from fragmented forest.
Since analysis using the polarimetric texture signature indicates that the xpol45 configuration should give
maximum of the normalized variance of backscatter for this type of polarimetric mixture [4], while the
HV configuration should give a minimum, we build feature vectors corresponding to these polarization
The mask used for mixing the two classes in the simulated image is shown in Fig. 1 and the HV
backscatter image is reported in Fig. 1(b). This image suggests the difficulty of the segmentation problem,
in particular if only one polarization were available and only backscatter were used in the classification.
In this experiment we use 2 scales and a smoothing kernel with m = 4 . The feature vector component
related to the variance of the smooth image at scale 2 is shown in Fig. 1(c).
Both the wavelet coefficients and the variance of the smooth image feature a significant between-regions
variance. At the coarser scale the smooth image gives mostly information on the edge between regions, a
fact which points at the role of scale in texture analysis.
A segmentation trial has been set up, using the approximated multiple discriminant transform (AMDT)
for feature reduction (see section 2.2) and a K-means clustering procedure with 5 classes. At first we
analyze the results when clustering is applied to the output of the first KLT transform (see Fig. 1(de)).
The fragmented forest region is well segmented by the categories corresponding to high textures. Pixels
corresponding to these clusters are shown in Fig. 1(d ), while pixels belonging to the low-texture
categories are shown in Fig. 1(e). There are some commissions that are tantamount to 4% of the class
population. Omissions are 2.8 %.
Segmentation results for the same 5 categories problem but when the full AMDT procedure is used are
shown in Fig. 1(f). Here the segmentation accuracy is much higher (no omission or commis sion errors) at
the expense of loss in resolution due to the additional smoothing at a scale factor 2 m which is a requisite
for obtaining an approximate MDT solution.
The question whether texture can play a significant role in thematic information extraction from SAR data
is on the table of the SAR community since a long time. Indeed most of the applications reported in the
literature rely on classifications algorithms that assume a piece-wise continuous model of the underlying
radar reflectivity. However some natural targets do not exhibit radiometrically “pure” backscatter, and
therefore difficulties arise in mapping these targets into a pure category unless texture is taken specifically
into account. A case in point is provided by the closed canopy dense tropical forests in Central Africa. A
second case, in an entirely different application area, is provided by surfactant slicks in the ocean. In the
following we illustrate some results related to those two cases, thus hoping to contribute an answer to the
question about the relevance of texture.
4.1 Swamp Forest Mapping in the Congo River Floodplain
In tropical Africa and specifically in the Congo River floodplain, the swamp forest upper canopy appears
as a smooth surface when imaged by a radar with a resolution of some 30 m. On the other hand the
primary lowland rain forest’s canopy is much more inhomogeneous, being composed by many and
diverse competing species. Therefore it appears as a more textured surface at the resolution of the radar.
This effect is enhanced in imagery acquired by a C-band radar, which, due to lower penetration, senses
only a shallow slice of the upper canopy. This difference in backscatter statistical properties can be
exploited to discriminate in a very effective way the two forest classes based on texture measures. The
physics behind this effect is explained in [5].
In Fig. 2(a) an ESA ERS backscatter image of an area in the Congo floodplain is shown. The image was
resized by pixel aggregation to 100 m pixel spacing.
The approximate multi discriminant transform reduces the feature vector to one component, which is
shown in Fig. 2(b). Notice how in the backscatter image the differences in smoothness between the two
types of forest are barely vis ible, whereas the texture image captures neatly these differences (see dark
ribbon along the river, and the fern-like patterns spreading from the river bed into the forest domain).The
ratio of between- to within-regions variances in the texture image is high enough to allow us to use this
image as input product for the segmentation stage. In our trial a K-Means algorithm with four clusters
produces the result shown in Fig. 2    (c). Notice how the swamp forest fine structures along the water
drainage patterns are well preserved in the map.
The swamp forests of the Congo basin were previously mapped in connection with the Global Land
Cover Map 2000, using texture features based on block window estimates of the K distribution order
parameter. Comparison of the relative accuracy of the two approaches is under way.
4.2 Oil Slicks Detection
It is well known that detection of oil slicks based on one-point statistics of SAR backscatter values only is
marred by the fact that several other features, such as waves, ship wakes, floor topography, present
similar range of values. To achieve reliable detection the classification algorithm must usually be
compounded with contextual information (e.g. two-point statistics, morphological operators) and/or take
into account ancillary data, such as ship routes and wind data.
                    a)                                       b)                                       c)

Fig. 2. Swamp forest mapping in the Congo basin. The original ERS -1 backscatter image (down -sampled to 100 m) (a).
The texture image at the output of the approximate multi discriminant transform (b). Segmentation into 4 clusters using
K-Means clustering (c).
We want here to explore the potential of a simpler approach based on the enrichment of the feature vector
by texture measures and polarimetric target decomposition.
For the experiment, a fully polarimetric SIR-C image of the English Channel, C-band, is used (courtesy
JRC-EMSL). The synthesized HH amplitude image is shown in Fig. 3             (a). We can identify by visual
inspection a fairly big oil slick (highlighted area A), features associated with waves (highlighted area B),
and several strong point targets related to ships (highlighted area C).
In a first trial the feature vector was constructed using the xpol45 configuration. More precisely, the
original covariance matrix data (multi-look complex) of the JPL distribution kit were averaged by 4x4
aggregation, and converted to the Stokes operator representation. The wavelet frame operator was then
applied to the Stokes data (see section 2.3). Finally by polarimetric synthesis the xpol45 amplitude image
was generated. The pixel spacing is 50 m. In the feature vector construction a smoothing kernel dilated by
m = 4 is used. The variance of the wavelet coefficients along the image row direction W2x at scale
2 2 is shown in Fig. 3(b) and the component along the column direction W2y in Fig. 3(c). The point
targets are very well captured, and the contour of the slick is clearly identified. However if we look at
W2y we see that other linear structures, like the wave front corresponding to area B, give similar
response. Therefore identification of the slick would require a second level of pattern recognition.
As a second trial we added to the feature vector texture meausures derived from the Cloude-Pottier
decomposition based on the coherency matrix. In particular we us the entropy image. It turns out that
entropy increases significantly in the presence of surfactant slicks due to dampening of resonant Bragg
scattering that dominates in clear waters [6]. The wavelet coefficients related to entropy allow now for the
slick detection with a good signal to nois e ratio (SNR = 78.4 dB) and can be captured unequivocally by
simple image thresholding (Fig. 4(a )). The other side of the coin is that entropy gives a worse signal to
noise for the point scatterers. The principal component of the composite feature vector (entropy and the
xpol 45 intensity) achieves optimal detection of both features (see Fig. 4(b)) with SNR=31.1 dB for the
edges related to the oil slicks.
Results of this experiment seem to indicate that the combination of polarimetry, coherent decompositions
and texture measures provide a robust vehicle for the oil-slick detection problem. This qualitative
statement will have to be of course quantified by a detection performance analysis based on a rigorous
validation exercise.
We have presented a simple but effective technique for textural segmentation of SAR imagery. The
technique exploits the ability of a wavelet frame to capture “correlation“ properties of pixels in
neighborhoods of different sizes, the aggregation of features extracted from SAR backscatter at different
polarizations, and a mechanism to represent the feature vector in a space that approximately optimizes the
between-regions variance to the within-regions variance. The merit of the approach lies on: i) simplicity
of the basic mechanism for constructing texture measures, which does not rest on a specific image model
for the underlying radar reflectivity; ii) efficient algorithms for the numerical imple mentation; iii) ease of
use (only a few free parameters to be set). The combination of these aspects make the approach flexible
enough to cope with different thematic applications and computationally suitable also for large scale real-
life radar mapping applications.



      Fig. 3. SIR -C C-band multi-look HH amplitude image acquired over the English Channel (a) (courtesy JRC-EMSL).
      Variance of the wavelet coefficients along the image row direction, estimated using the xpol45 and a smoothing
      factor m=4 (b). Variance of the wavelet coefficients in the image column direction (c).



       Fig. 4. Segmentation by image thresholding of the feature vector component holding the normalized variance of
       entropy (a). First component if the KLT transform of the feature vector constructed using both the xpol45 synthesized
       amplitude and the entropy (b).


     [1]     M. Unser, “Texture classification and segmentation using wavelet frames,” IEEE Trans. Image
             Processing, vol. 4, no. 11, pp. 1549-1560, Nov. 1995.
     [2]     M. Unser and M. Eden, “Multiresolution feature extraction and selection for texture
             segmentation”, IEEE Trans. Image Processing, vol. 11, no. 7, pp. 717-728, July 1989.
     [3]     S. Mallat and S. Zhong, “Characterization of signals from multiscale edges,” IEEE Trans.
             Pattern Anal. Machine Intell., vol. 14, no. 7, pp. 710-732, 1992.
     [4]     G.F. De Grandi, J-S. Lee, D. Schuler, and E. Nezry, “Texture and speckle statistics in
             polarimetric SAR synthesized images,” IEEE Trans. Geosci. And Remote Sensing, vol. 41, no. 9,
             pp. 2070-2088, Sept. 2003.
     [5]     M.L. Williams, S. Quegan, D. Blacknell, “Intensity statistics in the distorted Born
             approximation: application to C-band images of woodland”, Waves in Random Media, vol. 7, no.
             4, 1997.
     [6]     D.L. Schuler, J.S. Lee, and G. De Grandi, “Spiral Eddy Detection using Surfactant Slick Patterns
             and Polarimetric SAR Image Decomposition Techniques”, Proc. Geoscience and Remote
             Sensing Symposium ,IGARSS'04, IEEE International, vol 1, 20-24 Sept. 2004.

To top