Visualization of Multidimensional Shape and Texture Features in Laser Range Data Using Complex-Valued Gabor Wavelets.
AbstractThis paper describes a new method for visualization and analysis of multivariate laser range data using complex-valued non-orthogonal Gabor wavelets, principal component analysis and a topological mapping network. The initial data set that provides both shape and texture information is encoded in terms of both amplitude and phase of a complex valued 2D image function. A set of carefully designed oriented Gabor filters performs a decomposition of the data and allows for retrieving local shape and texture features. The feature vector obtained from this method is multidimensional and in order to evaluate similar data features, further subspace methods to transform the data onto visualizable attributes, such as R, G, B, have to be determined. For this purpose, a feature-based visualization pipeline is proposed consisting of principal component analysis, normalization and a topological mapping network. This process finally renders a R, G, B subspace representation of the multidimensional feature vector. Our method is primarily applied to the visual analysis of features in human faces_but is not restricted to that.
Introduction
One of the major challenges of vision research was, and still is, to develop methods for the automatic
modeling of complex geometric objects or scenes. In spite of countless efforts during the
last decades [2], [42], [32], [36] there is not yet a generic solution to this problem. Electronic
photogrammetry, however, has invented active vision methods, like laser range finders [40], that
are widely used to yield complex surface shape information elegantly and efficiently. Once only
available in military applications for robust terrain following navigation systems, different types
of laser scanners nowadays serve as helpful tools for automatic modeling tasks in many different
areas [57]. Recent advancements in the development of laser scanners allow the capture of shape
and color information of the object in question [59]. Hence, sophisticated methods to encode and
analyze range image maps have become essential requirements [5], [62].
Human faces in particular are attractive objects for laser range scanners [53], [45], since they contain
complex geometric features and intrinsic symmetries on the one hand [33], and are very appropriate
for human visual analysis on the other hand. Accurate geometric data from human faces
can be used to enhance the information in simple photographs and allow for the development of
higher quality techniques for facial surgery, facial recognition, facial reconstruction, simulation
of aging, etc. The question arises how to encode this shape and texture data into representations
that visualize important data features.
A lot of research has been done to find appropriate geometric model descriptions for regular [18]
and irregular [19], [47] surface data and nonuniform rational B-splines (NURBS) have been
shown to be a flexible scheme for controlling complex shapes. Yet, it turns out to be very difficult
to obtain robust features from such descriptions. Moreover, when considering range data, there
needs to be a way to also encode and analyze texture. Apart from that, wavelets and multiresolution
analysis, as proposed by [14], [29] or [41] have becomevery attractive for many applications.
In computer graphics mostly orthonormal wavelets and their separable 2D and 3D extensions
have been used for hierarchical data decomposition and approximation. Taking advantage of the
compact coding scheme and of the local support of the basis functions in [44], [26] and [60], different
approaches for volume rendering and isosurface reconstruction have been suggested. [51],
for instance, used wavelets for fast radiosity computations and [34] for control of 3D morphing.
Especially in image processing, the power of wavelets has been investigated for feature extraction
and analysis [7], [11], [55]. While in [11] orthonormal wavelets are stressed, [15] employes
nonorthogonal and nonseparable 2D-Gabor wavelets for image analysis. Other interesting re-search
is reported by [8], who combines graph matching with Gabor functions for face recogni-
tion. It is evident that in most cases the requirement for compact support and orthonormality
along with a smooth wavelet shape and a dense spatial orientation staggering for high performance
in both analysis and coding presents a problem. In those cases, where the analysis properties
of the basis function are superior to the approximation and coding behavior, non-orthogonal
wavelets have to be considered as well.
This paper addresses two fundamental areas of scientific visualization: First, it describes how to
extract multidimensional features from complex data sets, such as laser range images, using a
complex-valued decomposition scheme with Gabor wavelets. Secondly, it provides a generic
scheme to visualize these features in orthogonal subspaces. For this purpose we consider shape
and texture as amplitude and phase of a complex-valued 2D image function and perform a hierarchical
decomposition with a carefully designed set of Gabor wavelets. Although there is no
straightforward way to perfectly reconstruct the initial data set from the filter pyramid, the Gabor
wavelet has been shown to be a much more powerful feature descriptor than an orthonormal
wavelet [25]. Specifically in contrast to the non-separable 2D Gabor filters, the tensor product
extensions of orthogonal or biorthogonal wavelets [55] accomplish compact coding, but their
directional selectivity is much poorer with regard to the diagonal components. Moreover, Gabor
wavelets meet uniquely the lower bound of space-frequency resolution as it is stated by the Heisenberg
principle. Hence, Gabor wavelets provide advantages for any type of feature analysis,
that bases on local spectral estimates, such as the one we introduce in this paper. The complex-valued
interpretation of the image function allows to encode the range and texture information
elegantly and to compute the Gabor decomposition efficiently via FFT methods. In particular,
convolutions with complex-valued Gabor functions turn out to become simple multiplications
with Gaussians in Fourier space.
However, once the decomposition of the data is computed, the result is a multidimensional feature
vector at each surface point and the problem arises how to inspect it visually. To this end,
subspace projections have to be applied in terms of a principal component analysis of the wavelet
features. After normalization we use a modified topological mapping neural network, as proposed
in [23] or [28] to accomplish cluster analysis and mapping onto the R,G,B color space. The
generic feature-based visualization pipeline, introduced here, exemplifies to some degree the demand
for complex multidimensional feature-based visualization techniques, as stated in [50].
The organization of the paper is as follows: First of all we briefly elaborate the mathematical basics
of the Gabor function and its relationship to non-orthogonal wavelet decompositions. This
section also elucidates how to build 2DGabor wavelets in the frequency domain. Section 3 sheds
light on our encoding scheme for laser range data sets and illustrates how to get local data features
from this transform. Chapter 4 addresses the feature based visualization pipeline we employ with
special emphasis on principal component analysis and topological mapping networks. Finally,
in chapter 5 we present results on the performance of our method applied to range data of human
faces. R,G,B representations of multidimensional face features are depicted as derived from the
Gabor wavelets and how this method generalizes onto different data sets is outlined.
Mathematical Foundations
2.1 The Gabor Function
The Gabor function [21] used to decompose the range data in our approach is of fundamental
importance in signal processing [1] and has been widely used in a broad range of applications
[7], [15]. It provides an effective way to analyze images and has been elaborated as a framework
for understanding the orientation and spatial frequency selective properties of simple cortical
neurons [16]. The use of Gabor functions for data analysis therefore has its foundations both in
signal analysis and in biology. Gabor functions of different frequency range and orientation can
provide a non-orthogonal function basis for any finite energy function f # L 2 R n [3] and one way
to rapidly obtain data features is to compute a Gabor transform as an expansion of f. The 2D version
of the complex-valued Gabor function can be expressed as below. Its basis consists of a
Gaussian envelope of a harmonic oscillation term. Important inherent advantages are infinite
smoothness and exponential decay in frequency, but it has to be stated again that there is no
straightforward way to reconstruct the data from the expansion basis. The Gabor function can
be defined in both the spatial and frequency domain, although the latter one turns out to be easier.
In Cartesian coordinates, we obtain # G (x,y):
(1)
stand for the translation of its origin and (u 0 , specify the modulation coordinates in
the frequency domain. The effective width and length are given by (# x , # y ) that specify the elliptic
envelope.
Gabor functions are localized in both space and in frequency and they moreover uniquely achieve
the theoretical lower bound of uncertainty as it is dictated by the Heisenberg relation:
(2)
The resolutions in space # xy and in frequency # uv are computed by the 2nd order moments and
correspond to the variances and covariances of the Gabor functions.
dx dy
(v
Fourier transform of # G (x,y)
y, u, v: first order moments (centers)
does not yet account for spatial orientation, frequency # and orientation # of the
Gabor function are defined by
The coordinates x and y can be rotated by replacing them as
sin # 0
respectively. Note that rotation leads to cross terms xy and to a nonseparable expression.
The complex oscillation term represents the orthogonal quadrature components sine and cosine.
Fig. 1 illustrates the real (symmetric) and imaginary (antisymmetric) component of the Gabor
function from eq. 1. The shape of both functions strongly resembles an oriented wavelet function.
Fig. 1. a) Real (cosine) and
imaginary (sine) component of the complex valued 2D Gabor
function in the spatial domain for
2.
The effective number of oscillations under the Gaussian envelope can be expressed as # 0.
Fortunately, the Fourier transform G(u,v) of # G (x, y) has exactly the same functional form with
inverted or interchanged parameters. For functions located at the origin it reduces to a simple
Gaussian in the frequency plane. Introducing the spatial frequencies (u, v) we obtain the following
-#
with #
This relationship which is depicted in fig. 2, shows the Fourier transforms of the sine and of the
cosine wavelet respectively. Due to the symmetry and antisymmetry of the harmonics the negative
Gaussian peaks cancel each other and the remaining function reduces to one Gaussian controlling
7the wavelet's position, size and orientation. Evidently, implementations of Gabor filters
can be accomplished easily employing eq. 7.
Fig. 2. a) Fourier transform of the Gabor function of fig. 1.
It's decomposition into b) sine and c) cosine components.
From there, it is straightforward to analyze data using Gabor functions. The so-called Gabor-
transform [1] as it is often stressed in the signal processing literature turns out to be a short time
Fourier Transform (STFT) using a basis of Gabor functions. In the 2D-case, we get
-#
-#
where w is a Gaussian window of constant size # x , #
2.2 Definition of a Gabor Family
The constant size of the envelope that is assumed in eq. 8 leads to an uniform partitioning of the
frequency plane when expanding f. Thus the localization of the Gabor transform is equal through-out
the entire spectrum. However, sophisticated data analysis should provide a high spatial resolution
for high frequencies and a low resolution of lower frequencies. Therefore, we have to adapt
the size of the Gaussian window to its frequency position. Due to the scaling theorem of the Fourier
transform shifting and increasing the width and length of the Gaussian in frequency space go
in hand with different modulation and a decrease of width and length of the wavelet's envelope
in the spatial domain. In fig. 3a both frequency position and effective width of the Gaussian were
doubled. Fig. 3b illustrates the corresponding real part of the wavelet. Obviously, we can generate
a whole family of self-similar functions of this type, simply by shifting and scaling a Gaussian
in the frequency plane.
Fig. 3. The effect of scaling and shifting in frequency space:
a) Gabor function of different size, position and frequency in
Fourier domain.
Real part of the resulting wavelet in spatial domain
a # 2).
Hence, a scheme to construct a set of self-similar Gabor functions as an expansion basis for f can
be set up.
The idea of building an expansion basis from self-similar functions of different size and orientation
immediately leads to the generic wavelet transform, whose mathematical principles are described
in [13], [41], [37] or [14]. Basically, it can be considered as the inner product of any finite
energy function f with a set of self-similar basis functions # a,b , that are derived from each other
by scaling and shifting of one prototype using the parameters a, b # respectively. Its 2D representation
can be addressed as:
-#
with
# a,b (x, y) #|a x a y |
a x
a
The brackets <,> denote the inner product operation.
Its Fourier transform is given by
# a,b (u, v) # |a x a y |
Note that the normalization forces the bases to hold Parseval's energy equivalent.
In wavelet theory, the bases are usually assumed to satisfy the constraints of orthonormality
and band-pass behavior:
Note furthermore, that due to its definition the Gabor function has an intrinsic DC-fraction in
it's real (cosine) component. Hence, any wavelet derived from it will not meet eq. 13. More spe-
cifically, any two scaled and shifted versions of Gabor type will not satisfy eq. 12 as well.
2.3 Frequency Decomposition with Gabor Wavelets
Depending on the job for which they are tailored, additional constraints on the bases are often
demanded, such as compact support, smooth shape, fast decay, closed-form descriptions, etc. It
has turned out that the construction of bases satisfying these criteria is a non-trivial optimization
process and usually a compromise has to be found, as semi-orthonormality [56]. However, relaxing
the orthonormality and bandpass constraint, we can construct a wavelet basis using the Gabor
functions of different scale, shift and orientation [15]. In this case a complete set of self-similar
functions is defined by
Here, we assume a circular scaling with a a. The parameters represent transformations
of the mother wavelet in scale (a,m), translation (p,q) and orientation (#). The index m runs from
lower to higher frequencies in our paper, since the wavelets are constructed accordingly as in appendix
A.
As stated earlier, the most effective way to construct a multiresolution wavelet basis is to find
appropriate criteria for specifying the Gaussian in the frequency plane. Once constraints for setting
a and # are found, a complete set of Gaussians represents our filter pyramid and the decomposition
is obtained by subsequent filtering and inverse Fourier transform. In our implementa-
0tion, we define the decomposition using a lower frequency bound # 0 and an angular resolution

Figure

4 shows how the frequency plane is decomposed by superimposing Gaussian filters.
The mathematical details of the construction rules we propose here are described in Appendix
A at the end of this paper.
Fig. 4. Decomposition of the frequency plane with Gabor wavelets. The
parameters are #
# a # according to Appendix
A.
Let i(x,y) be a 2D image function and I(u,v) be it's Fourier transform, we obtain a set of convolution
products g m,i (x,y) by multiplying G m,i (u,v) with I(u,v) and subsequent inverse Fourier transform

-#
-#
G m,i (u, v) I(u, v) e
where G m,i (u,v) denotes the Fourier transform of # G
mpq# (x, y) according to Appendix A.
This scheme provides a more effective way of building up the wavelet pyramid - a way which
is preferable to computing all inner products of the scaled and shifted Gabor functions
explicitly. This is because the convolution products computed by eq.
the inner products of all shifted versions of a wavelet given at m and #. The implementation
scheme from above corresponds to a filter bank, indicated in fig. 5.
Fig. 5. Filter bank, that implements the Gabor decomposition.
The method introduced represents the feature extraction module in our visualization pipeline and
can be employed to decompose laser range data sets.
It should also be noted here that wrap-around problems and cut-off errors arising from the FFT
of finite data sets have to be avoided. Therefore, we strongly recommend reflecting the data at
their boundaries and fading them out exponentially. All computations in this paper are performed
that way.
3 Shape and Texture Coding using Gabor Wavelets
3.1 Characteristics of Laser Range Data
This section elaborates the decomposition of data sets using the Gabor wavelets introduced
above. Although our further investigations are focussed on laser range data, the decomposition
scheme as well as the feature based visualization pipeline is not restricted to that. We will show
that in particular when dealing with complex-valued data this method provides an elegant and
efficient method for visual data inspection. Consequently, any type of shape and texture informa-
tion, such as terrain data and aerial photographs, are well suited for this scheme.
We employ laser range data from a Cyberware laser scanner, which provides highly accurate
shape and color values of scanned objects. The data sets can be interpreted as R,G,B range images
defined on a cylindrical coordinate system, as depicted in fig. 6. The grid size is 512 x 480 at a
resolution in the range of about 0.3 mm. A typical scan is rendered in fig. 7, that shows the range
information (fig.7a) and the color mapped as a texture onto the shape (fig.7b).
The scanner moves a diffuse line light source around the object synchronously to the range detec-
tor. This is the reason why the illumination is preserved constant at each surface point and the
texture information can be interpreted in terms of local surface properties. It is not affected by
specific illumination conditions. It has to be mentioned, that the albedo of human facial skin is
very high and a respective reflection model has been developed by [31].
Due to the cylindrical coordinate system of the scanner, the data have to be projected into Cartesian
coordinates representing the 2D image function in (x,y). The projection can be considered
as taking the normalized z S -coordinate of the scanner system as the y-coordinate of the image
and the angular coordinate # S as it's x-coordinate. Consequently, this operation unrols the cylinder
onto a plane, as in fig. 6b.
where (z S , r S , # S ) denote the cylinder coordinates of the scanning system.
Fig. 6. a) Coordinate system employed for data acquisition with the laser
scanner.
b) Transformation of the range data onto a rectangular grid.
Fig. 7. Shape and texture obtained from a female face (Sylvia) and its
encoding scheme.
a) Range rendered with Gouraud-shading.
Texture-mapped color information.
c) Encoded range information.
d) Encoded grey-level texture information.
(raw laser range data: courtesy provided by Computer Graphics
Center, Darmstadt, Germany)
3.2 Complex-valued Coding of Shape and Texture Data
The question arises how to encode both color and shape in order to decompose them efficiently
with the Gabor pyramid. A straightforward way would be to treat R,G,B and range separately
in terms of four different decompositions. This is however very expensive and provides highly
correlated and high-dimensional features. Referring back to section 2, the Gabor function is composed
of a symmetric and an antisymmetric oscillation term, that account for the real and imaginary
components of the data, respectively. Moreover, sine and cosine wavelets are orthogonal to
each other. Hence, in our approach range and color are interpreted as amplitude and phase of a
complex-valued image function that is convoluted with a complex-valued wavelet.
Let r(x,y), g(x,y), b(x,y), s(x,y) be the color and range information at a particular point (x,y). In
order to compute the amplitude a(x,y) we extract luminance according to the rules of colorimetry
[61]:
In fact, the CCD sensor is not calibrated explicitly, but most systems lie near the CIE standard
of that equation. The range function is taken immediately as phase of the image function and the
description of i(x,y) is obtained as:
with its real part i R (x,y) and its imaginary part i I (x,y):
The range function is assumed to be normalized between 0 and 1.
Fig. 7c,d depicts also, how the data appears, as encoded by means of phase and amplitude. Although
information is lost by transforming color into the scalar-valued luminance, shape and
grey values are, however, still superior to color information when extracting features from objects
[32]. Furthermore the efficient encoding scheme provides computational advantages with
FFT methods.
Once the data is encoded as in eq. 19, we can compute it's Fourier transform I(u,v) and apply the
Gabor decomposition, as defined with eq. 16. This is illustrated in fig. 8, which presents a decomposition
of the data from fig. 7. The convolution products g m,i (x,y) are arranged according to the
orientation preference of their wavelet. The depth of the pyramid is M=4, at #/4 and the
DC parts are presented in the middle of each picture, respectively. Due to the complex-valued
results of g m,i (x,y), the inverse Fourier transform produces two resolution pyramids, one for the
real and one for the imaginary component, respectively. These pictures elucidate the orientation
and frequency selectivity of the wavelets. The low-pass functions at m=0 only account for rough
global features since their spatial location is rather bad. The high-pass wavelets at m=3, however,
reveal the fine-grain details in the data. Obviously, due to the non-orthogonality, the different
responses are correlated.
Fig. 8. Decomposition of Sylvia's image:
a) Real part of the Gabor pyramid.
Imaginary part of the Gabor pyramid.
Note, that the pyramid is not subsampled in this picture. A strict representation according to the
definitions of the WT, however would require further sampling of the g m,i (x,y) at the current rate
given by the choice of # and #.
3.3 Extracting Local Data Features
Evidently, this complex valued decomposition pyramid from fig. 8 gives a multiresolution view
onto the data. Due to the localization properties of the Gabor function a set of multidimensional
local data features can easily be derived from it just by evaluating the convolution products
m,i (x,y) at any point of interest. Fig. 9 illustrates the respective procedure of local feature extrac-
tion. Here, we assume the decomposition pyramid arranged into a multidimensional set, where
the dimension index k ranges from 1 (DC-part) to In the further sections
of the paper we refer to a feature vector as the set of scalar-valued coefficients
re
evaluated at any data point according to the illustrations of fig. 9.
This interpretation of the features implicitely transforms the former complex-valued vector into
a scalar-valued one of double length. Hence, the following steps in our pipeline rely on scalar-valued
computations although PCA, normalization and Kohonen map could be extended. This
however would restrict the method very much to the specific application of this paper without
any gain.
Fig. 9. Extraction of local data features from the wavelet decomposition

Although the reconstruction of the data from the Gabor decomposition requires the non-trivial
computation of a dual frame [3], it is still of interest to locally expand the data, as in fig. 10. In
this picture, the feature vectors are superimposed only within the area indicated around the right
eye and all others are neglected. Range and amplitude values are computed according to eq. 19
and mapped as displacements (fig. 10a) and texture (fig. 10b) onto a plain cylinder. Obviously,
they encode important data features, such as curvature or texture modulations around the center
of interest.
Fig. 10. Local reconstruction of the face shape and texture superimposing
the wavelet features within the demarcation of fig. 9:
a) Range around right eye.
b) Texture around right eye.
4 Visualization of Multidimensional Face Features
4.1

Overview

In the section above, we proposed a method for getting multidimensional data features by decomposing
them with multiresolution Gabor wavelets. Unfortunately, there is no immediate way to
visualize the topology of this feature space. Hence, we have to perform additional data analysis
steps to map the most important features down into a subspace that can be visualized.
The problem of visualizing multidimensional features in complex dataspaces is one of the key
issues in scientific visualization [43], [50] and is addressed in many applications, such as fluid
flow and tensor visualization [35], statistical data visualization [63] or multispectral imaging
[28]. The main point is to find orthogonal projection methods that preserve the most important
data features. One basic method that addresses this problem is the principal component analysis
(PCA). It optimizes the mapping procedure in a least-square error sense based on data statistics.
Consequently, we propose a visualization pipeline as figured out in fig. 11. This concept comprises
an embedded framework using different methods proposed by one of the authors in [25]
and [27]. Once the feature vectors g l are extracted from the decomposition, they are fed into a
subsequent processing pipeline of PCA-analysis, normalization and clustering into R,G,B color
space.
As mentioned earlier, the Gabor wavelets provide non-orthonormal expansions of the data and
the resulting feature vector is more or less highly correlated, depending on the choice of # and
#. The decorrelation can be accomplished by further expansion into its principal components.
Fig. 11. Feature-based visualization pipeline employed.
This first breakdown of the number of data dimensions is controlled by thresholding the eigenvalues
associated to each eigenvector. The coefficients of the feature vector transformed by the eigenvector
matrix are decorrelated and have to be normalized. Visualization of the data is figured
out by a self-organizing neural network, as proposed in [23]. This network clusters the data and
maps them automatically onto the R,G,B color space with the constraint of coding similar data
features in terms of similar colors. The method has the advantage of reducing similar features
to a limited set of clusters, rather than to visualize a whole subspace. Hence, the results obtained
are much more expressive than an immediate mapping of the first three eigenvectors into R,G,B.
4.2 Visualizing Principal Components
Generally, dimensionality reduction is strongly connected with finding subspaces satisfying respective
optimization criteria. In statistical data analysis [20], we can find several techniques for
this task. One of the most famous methods, taking the optimization as an eigenvalue problem,
is the principal component analysis or Karhunen-Lo-ve expansion. In general, this method aims
to find a subset of principal directions in a data set of any statistical distribution. The corresponding
basis vectors satisfy the constraint of orthonormality and the subspace defined by these vectors
diagonalizes the covariance matrix of the data. The transformation of an initial data set into
the space of principal components can be formulated as follows:
Let{ g l }: be a data set of the dimension K and g a mean vector:
The principal component analysis essentially looks for a set of orthonormal vectors with the
constraint of
and
eigenvectors and eigenvalues of the covariance matrix C, which can be
estimated as
Note, that its dimension is given by K # K.
Hence, the set of diagonalizes the covariance matrix, as
c 21 .
c K1
c K2
c 1K
c 2K .
k1 .
There are numerous methods for numerically handling this approach for huge covariance matrices
[20], such as singular value decomposition.
After solving this eigenvector problem, the feature vectors g l have to be projected into the eigen-space
spanned by the
and the coordinates in eigenspace are obtained by #
Conversely, any datum g l can be expanded as a linear combination of
The dimensionality reduction can now be achieved by ordering the eigenvectors according to
the absolute value of their corresponding eigenvalue # k and by taking only the most significant
subset for the data expansion in eq. 27.
The vector # can be interpreted as the decorrelated feature vector.
The remaining eigenvectors define a subspace that minimizes the average error of information
lost by the reduction of the number of dimensions.
This method is very popular for any kind of signal analysis and image coding, i.e. for removing
correlation from data like face images [23].
The principles of PCA force most of the data energy to be concentrated in a few significant eigen-
vectors. Therefore, once the PCA based feature vectors #
are computed for each initial g l they
have to be normalized for further processing. Since the a priori probability of these vectors is un-
known, we assume uniform distribution and normalize their components according to:
l
l
~
Fig. 12 shows the 6 most significant feature coordinates #
for l # {x, y} and the 3 least
significant ones as computed from the laser range data set from fig. 7. The respective eigenvalues
are indicated below the picture, as well.
Fig. 12. a) Coordinates of the feature vectors of Sylvia's image encoded
in eigenspace (see also fig. 7).
One straightforward way to visualize the multidimensional feature space of #
as it is derived
from the Gabor pyramid is to encode it into R,G,B. This is illustrated in fig. 13 where a smooth
color representation of the 3 most significant coordinates in eigenspace is depicted. Although this
presentation reveals a first sketch of similarities in the multidimensional feature space in terms
of similar colors, yet it does not provide a sophisticated feature extraction and the visual interpretation
of the resulting images remains difficult. Hence, additional clustering algorithms are
required to group similarities in the data. This is figured out by a topological mapping neural net-work
explained below.
Fig. 13. Distribution of the coordinates of the 3 most important eigen-vectors
of Sylvia's image encoded into R, G and B, respectively.
4.3 Topological Mapping
The 3D extension of the Kohonen map was introduced by the authors in [22] and employed in
different applications, such as [23] or [9]. It is basically a self-organizing network which is trained
with or without supervision. It aims at an organization of the input patterns to a topological structure
represented by its neurons, where the relations between different patterns are preserved as
much as possible.
The Kohonen map is a two-layered network. The first layer of neurons can be considered as similar
to a group of sensors picking up the data, entirely connected to a second, 3D-layer: the competitive
layer. Figure 14 shows the topology of the network. The weights associated with the connections
are adjusted during training where only one single neuron in the competitive layer can
be active at a time. This neuron represents the cluster to which the data set belongs in the spirit
of c-means clustering. Due to the training rules explained below, the Euclidian distance between
two neurons reacting on different input patterns can be a measure that determines the similarity
of the two patterns.
The training of the network is determined by presenting feature vectors #
randomly to the input
layer of the network whose connection weight vectors m h of all competitive neurons h are initialized
by random values. We choose K input neurons according to the data dimension and define
a Euclidean distance d h between #
and m h with
The neuron e with the minimum distance is then activated, where
The updating of the weights m hk associated with the neurons is only performed within a proximity
of e. This proximity N e (t) is reduced with increasing training time t. The updating
conforms to Eq. 32, where #(t) represents a time-dependent learning rate:
This rule refers directly to c-means clustering [17], where the m h represent the centroids. The
time-dependent neighborhood can be described for rectangular areas as follows:
where
Fig. 14. a) 3D Kohonen feature map and the representation of its neurons
as entities in the RGB color space.
b) Class assignment to single neurons.
The network has two essential properties related to our problem:
First, it separates clusters of the presented data by mean vectors m h that are associated as weights
to the neurons. Secondly, it performs a topological ordering of the competitive neurons in a sense
that neighboring neurons in the layer represent similar clusters in multidimensional space and
thus it achieves a further dimensionality reduction.
Since the neurons in the competitive layer are ordered topologically, neighboring neurons react
to similar data vectors in the input space. The mapping can be interpreted as reduction of any K-dimensional
input pattern into 3 dimensions, preserving the topology of the data as much as pos-
sible. The resolution of this discrete 3D space is given by the number of competitive neurons,
i.e., by the clusters.
Referencing the axes of the cube with the primaries R, G, and B, each neuron of the competitive
layer represents a discrete entity in R, G, B space, i.e., it corresponds to a particular color triplet.
Thus, the similarity of the color provided by the reacting neuron refers to a neighborhood in K-dimensional
space.
Hence, the dynamics of the network are strongly related to the PCA analysis. In this case, how-
ever, not only all feature vectors are projected into an orthogonal subspace, but a discrete subset
of centroids is computed, representing the features. This topological structure is further represented
in terms of R,G,B colors. Therefore, in our approach, the usage of the PCA for huge data
sets is limited to a rough scaling of the dimensionality and any further fine-tuning is accomplished
by the Kohonen map.

Figure

15 demonstrates the impact of clustering and visualizing the features from Sylvia's image.
In order to realize the geometric relationships, the output of the Kohonen map is rendered as a
texture onto the facial shape. The parameters were selected with 8x8x8 neurons and M=4,
In fig. 15a the inital 26-dimensional feature vector g l was feed immediately
into the network, which furnishes the R,G,B values at each surface point. The network clearly
distinguishes shape features, such as local curvature or textural consistency between the left and
right hemisphere of the face. The similarity of colors clearly appears in spatially coherent regions
unless in those of frequency components spread throughout the spectrum. Regions of high, low
and strong directional curvature are separated, such as around mouth, nose and eyes. Fig. 15b
contrasts to those results using a first rough breakdown of the dimensions from 26 to 7 through
the usage of the PCA. On the one side, the regions extracted appear even more homogeneous,
since some of the fine grain information was cut off by the PCA. On the other hand, however,
the topology preservation of the Kohonen map does not reveal as evidently as in the previous pic-
ture. This is because the PCA encodes the variance of the features which has not to correspond
to coherent facial regions.
Fig. 15. Cluster analysis and visualization of features derived from Syl-
via's image:
a) Immediate subspace mapping using the Kohonen map only.
Cascaded mapping with PCA and Kohonen map.
In order to stress the influence of the texture information and the importance of the PCA, fig.
depicts results derived from a segmentation of the texture data of Sylvia's image without including
range. Figure 16a shows the b/w texture image and fig.16b its segmentation using PCA pre-
processing. Due to local distortion in the texture, the image looks more noisy than those presented
in fig. 15. The smoothing properties of the range data cannot be harvested in this case. These results
strike even more in Fig. 16c where the features are segmented without PCA preprocessing.
One important aspect is related to the invariance against affine transforms. In principle, rotation
and scaling can be interpreted as a shifting of the respective coefficients in our feature vector. This
is because of the properties in Fourier space, where rotation is transformed into a rotated spectrum
and scaling results in a frequency modulation. Hence, the respective basis functions in the multi-resolution
pyramid respond in different frequency channels. Due to the finite dimensionality of
the feature vector, some coefficients are pushed out of the vector, whereas others will enter at the
lower and higher frequency bounds. Because of the Euclidian metrics employed for segmenta-
tion, we have to keep the most important coefficients in our vector in order to preserve the topology
and the distances in feature space. For this reason, a PCA based preprocessing helps to extract
those coefficients and to keep the segmentation tolerant against rotation and scaling. Translation
invariance of the segmentation is achieved by the localized features. An illustration of this important
subject is given in fig. 16d - 16f, where a scaled and rotated version of the texture data is
presented and the segmentations are generalized from, fig. 16b,c. Evidently, the color similarity
and thus the preservation of similar features is performed much better when using the PCA. For
instance, the high frequency regions around the eyes, nostrils and mouth appear in bluish colors
in contrast to cheeks and forehead, that come out in orange, red and green.
Fig. 16. Influence of range data and affine transforms on the segmentation

a) Texure data of Sylvia's image.
b) Segmentation with PCA.
c) Segmentation without PCA.
d) Scaled and rotated version of the image.
Generalization with PCA.
f) Generalization without PCA.
Besides non-supervised clustering, classification is an essential part in feature-based visualiza-
tion. In this case, the method has to match previously trained patterns. For supervised classifica-
tion, however, each neuron - and also each cluster centroid - has to be assigned to a certain class,
depending on the definition of the user. This can be done by interactive selection of training areas
and by a majority voting of each neuron stimulated by the training set (see fig. 14b). After this,
each neuron has an associated class and the network is able to classify. However, during the organization
process the goal was to find a limited set of centroids representing the data in a c-means
sense rather than to find optimal placements of the decision boundaries in a minimum error
(Bayes) sense. For this reason, the network can be once again trained with a supervised postprocessing
in order to move corresponding centroids towards the Bayes decision boundary and to
improve the classification result.
This postprocessing is well known as learning vector quantization (LVQ 1, 2, and the self-stabilizing
type 3 can be described as follows:
For a given input pattern #
be the closest centroids to #
. Wemodify the centroids
according to
The window is defined as a symmetric area around the midplane of m h and m j . Then #
falls into
the window if
Where d h and d j are the two distances of #
# to m h and m j . The threshold # is calculated according
to eq. 40 and the relative window size # is chosen to be about 20%.
A detailed study of LVQ and of related methods can be found in [38].
Applications
5.1 Face Data Base
The goal of the following investigation reports how the previously introduced method performs
and generalizes on multiple range image data sets. For this purpose, a data base of 10 entities was
built up, as shown in fig. 17. It consists of both male and female range images that were clipped
to the depicted masks. Both range and color information is presented. The left 7 images were
employed to train the method, whereas the right 3 images were used to investigate the behavior
of the generalization. The Gabor decomposition was applied separately on each of the 7 training
samples and the Kohonen map was trained by randomly selecting 1 mil. feature samples from
the resulting pyramids. We used a constant map size of 8x8x8 neurons, i.e. clusters. The parameters
5for the decomposition were selected to # and # and an additional overlap
of 2 was added to tune the localization of the wavelets.
Fig. 17. Initial face data base a) with and b) without texture information.
Upper rows and left: Range images used for training,
lower rows right: Range images used for generalization.
5.2 Results
In order to emphasize the influence of the PCA figure 18a illustrates the results obtained by the
method first without PCA. In this case, the features were mapped immediately from 34 dimensions
to R,G,B using the Kohonen map. Obviously, the wavelets extract similar facial features
that appear as similar colors on each face, such as tip of the nose in red-orange colors, or the
cheeks that come up either in blue-green or violet. The mouth, mostly revealed in red-pink and
moustache boundaries are extracted, as are eye brows (light orange). It is interesting to verify,
that the left and right hemispheres are clearly distinguished and because of the different colors
associated, belong to different clusters. Some of the results reveal a more or less symmetric or
antisymmetric color map, depending on the person analyzed. At this point, our method can be
contrasted with curvature analysis, or slicing methods. In particular, when the results on the training
data are compared with those obtained on the generalization data set (lower row), this technique
evidently outperforms robust results, since color similarities referring to anatomic regions
are rather well preserved. In general, due to the influence of the high-frequency components,
within a single face there is hardly any large coherently clustered region. The similaritiy of colors
in neighbored clusters is, however, preserved in regions of low spatial frequencies.

Figure

18b depicts the same investigation, but in this case with a PCA downmapping of the dimensions
from an initial 34 to 7, preserving 99 % of the signal energy encoded in the eigenvec-
tors. The topology preservation of the Kohonen map is not as obvious as in the previous figure,
since the PCA accounts for the maximum variance in the features. One of the advantages of the
PCA is that the regions clustered by the neural network appear in most cases more homogeneous
and coherent because some fine grain detail information gets lost with the downscaling proce-
dure. In both computations, the most discontinuous regions appear around the eyes. This is due
to the discontinuouities in range, which spread the local spatial frequencies throughout the spec-
trum. The pupils of most of the test persons, however, are detected in similar colors.
Fig. 18. Face features as similarities of colors mapped as textures to the
range data.
a) without PCA and
including PCA.
Other interesting investigations can be accomplished by trying to retrieve similar regions and
landmarks in the face data base. For this purpose the Kohonen map was trained with a supervised
LVQ postprocessing on the results of fig. 18b using 100000 additional cycles. During this proce-
dure, only features within the training areas demarcated onto the face in fig. 19a were fed randomly
into the network. Four different feature types were selected, one each for tip of the nose
(blue), nostrils (yellow), pupils (green) and corner of the mouth (red). Figure 19b depicts the re-
sults, where the distribution of the trained features throughout each face is presented. It can be
stated, that all regions are properly retrieved, such as the tip of the nose or pupils. Considering
the especially small training areas, features similar to the corner of the mouths are detected mostly
around the eyes, eye brows, and moustache, where discontinouities along with high local curvature
in range can also be located. Furthermore, green regions appear systematically around the
throat of each person. This is due to an artifact arising from the linear interpolation of the scanned
data, which usually lacks in the green region. Although both nostrils and pupils are dark regions
in the texture map, they are clearly separated in all images. The diagram in fig. 20 compares again
the results within the training areas of fig. 19. Besides the two peaks within the nostrils the overall
error is relatively small.
Fig. 19. Detection of similar features extracted from a supervised classification
postprocess.
a) Training areas demarcated on faces.
b) Detection of the similar regions.
Fig. 20. Error rate for matching the significant regions within the training
areas.
Conversely, if a robust detection of previously learned training areas is required for matching or
recognition, a procedure can be set up as indicated in fig. 21. In this simulation, the method was
trained exclusively on Sylvia's image and the DC components were neglected in the feature vec-
tor. Hence, the segmentation looks noisy, no color similarities can be recognized, and the result
is badly shaped for visual inspection. The LVQ applied subsequently, was trained to distinguish,
left and right eye (green and blue), tip of the nose (red), and left and right corner of the mouth
(pink and yellow). Fig. 21c shows the performance of the detection. The landmarks are retrieved
almost perfectly without further errors because of the highly localized features.
In order to illustrate the importance of the range and texutre information for the outperformance
of the method the results are contrasted to those achieved by using range and texture data only,
as depicted in fig. 21d and fig. 21e. Although both computations are performed with the same
parameters, yet the retrieval of the landmarks is figured out much worse than in fig. 21c.
Note again, that the texture data is recorded under normalized conditions and hence represents
facial surface properties and is not affected by specific illumination.
Fig. 21. Detection of facial feature points in Sylvia's face:
a) RGB representation of the segmentation.
Training areas.
c) Performance of detection of the landmarks.
d) Performance of detection of the landmarks using range data
only.
e) Performance of detection of the landmarks using texture data
only.
All computations of this paper had been performed on a SGI-Indigo 2 workstation with a R4400
processor and 64 MB main memory. The file-oriented implementation of the method requires
to generate the appropriate filters, which took about 205 CPU seconds. Generation and normalization
of the feature vectors was figured out in 61 CPU seconds and the PCA took 48 s. The self-organization
of the Kohonen map was computed in 157 s with 200000 cycles and 7 coefficients
whereas the same procedure needed 472 s for the full feature vector without PCA.Hence the PCA
increases the computational performance of the method. The segmentation and classification was
done in 42 s respectively in 83 s. Note, that the figures are valid for one image.
6 Conclusion
We introduced a generic and robust method for feature-based visualization of multidimensional
data sets by decompositions with Gabor wavelets. Results have been presented, that show the capabilities
of the method, even in very complex data sets. Although it was primarily applied to
encode and analyze laser range data from human faces it is not restricted to and can be extended
to any kind of complex data, such as 3D- or 4D-velocity fields. However, one of the main aspects
when encoding data in the way we propose is, that there is no immediate way to perfectly reconstruct
8them in the general case. Furthermore, the current scheme does not accomplish any sub-sampling
of the pyramid and requires extensive memory. This is specifically a serious criteria
in multidimensional applications. [3] for instance, shows that wavelet frame theory, however,
defines bounds for the design of dual function spaces for a perfect reconstruction of Gabor
frames. Since Gabor wavelets are powerful feature extractors, further research has to be conducted
to find perfect reconstruction frames within these strict limitations and to achieve more
compact coding schemes. Other alternatives have to be considered, such as nonseparable extensions
of biorthonormal wavelets preserving directional selectivity.
Another important aspect is how the method performs under noise. Due to the spectral estima-
tion-like method, we expect muchmore robust results than those, for instance, given by curvature
analysis techniques for range data. In this context, a common strategy is to weight the wavelet
coefficients according to their frequency localization when using them for coding. For data anal-
ysis, however, the weights have to be selected much more carefully and strongly depend on the
signal characteristics. This leads immediately to higher-order spectral estimation methods,
which should be in the scope of future activities as well.

Acknowledgement

The authors wish to thank the Computer Graphics Center in Darmstadt, from where the raw laser
range data set was provided. Thanks also to S. Spanier-Mason for proofreading the English
manuscript. The authors are also grateful to the referees for providing very useful criticism.



--R


Academic Press


"A reference model for the visualization of multi-dimensional data,"
of the EUROGRAPHICS
"Invariant surface characteristics for 3D object recognition in range images,"
Computer Vision
"A survey of curve and surface methods in CAGD,"
Geometric Design
"Multichannel texture analysis using localized spatial filters,"


shington: IEEE

in medical imaging


"Texture analysis and classification with tree-structured wavelet transform,"



"The wavelet transform, time-frequency localization and signal analysis,"

"High confidence visual recognition of persons by a test of statistical independence,"
IEEE Trans.
"Two-dimensional spectral analysis of cortical receptive field profiles,"



"Scattered data interpolation and applications: A tutorial and survey,"
Hagen and D.

"Theory of communication,"

"Subspace methods for the visualization of multidimensional data sets"
shop
"Integrated volume rendering and data analysis in wavelet space,"

"Multiscale image texture analysis in wavelet space,"
First IEEE Conf.



of human faces - A case study for man-machine-communication
Singapore: World Scientific
"Visualization of multidimensional data sets using a neural network,"



"Visualization of large data sets"

"Reflection from layered surfaces due to subsurface scattering"
Computer Graphics

"Machine identification of human faces,"

"Wavelet-based volume morphing,"
"Visualization of vector and tensor data sets,"
Visualization, in Rosenblum

"An overview of wavelet based multiresolution analyses,"
Department of Mathematics
"The self-organizing map,"
"Texture classification by wavelet packet signatures,"


"A theory for multiresolution signal decomposition: The wavelet representation,"



Computer Graphics (special issue)
"Volumetric shape description of range data using 'blobby model',"
"Multiscale 3D edge representation of volume data by a DOGwavelet,"
on Volume
"Learning object models from appearance,"
"Scattered data modelling,"
"Visualization in scientific and engineering computation,"
"Equilibrium and interpolation solutions using wavelet bases,"


Scientific Visualization.
"Wavelet Radiosity"


"Recordering and visualizing complex shape from range data"
Tokyo: Springer
"Frequency domain volume rendering,"
"Multiresolution feature extraction and selection for texture segmentation,"
and
"A family of polynomial spline wavelet transforms,"

"Facial surface scanner,"
"Wavelets and filter banks: Theory and design,"


"A multiresolution framework for volume rendering,"

New York: John Wiley


"Visualizing structure in high-dimensional multivariate data,"
of Research and Development
--TR

--CTR
M. H. Gross , T. C. Sprenger , J. Finger, Visualizing Informationon a Sphere, Proceedings of the 1997 IEEE Symposium on Information Visualization (InfoVis '97), p.11, October 18-25, 1997
Ming Xi Tang, Visualization and Genetic Algorithms in Minimax Theory for Nonlinear Functionals, Journal of Scientific Computing, v.18 n.1, p.49-68, February
Rolf M. Koch , Markus H. Gross , Friedrich R. Carls , Daniel F. von Bren , George Fankhauser , Yoav I. H. Parish, Simulating facial surgery using finite element models, Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, p.421-428, August 1996
P. Cignoni , C. Montani , R. Scopigno , C. Rocchini, A general method for preserving attribute values on simplified meshes, Proceedings of the conference on Visualization '98, p.59-66, October 18-23, 1998, Research Triangle Park, North Carolina, United States
Markus H. Gross , Roger Gatti , Oliver Staadt, Fast Multiresolution Surface Meshing, Proceedings of the 6th conference on Visualization '95, p.135, October 29-November 03, 1995
Markus H. Gross , Oliver G. Staadt , Roger Gatti, Efficient Triangular Surface Approximations Using Wavelets and Quadtree Data Structures, IEEE Transactions on Visualization and Computer Graphics, v.2 n.2, p.130-143, June 1996
Sylvain Fischer , Filip roubek , Laurent Perrinet , Rafael Redondo , Gabriel Cristbal, Self-Invertible 2D Log-Gabor Wavelets, International Journal of Computer Vision, v.75 n.2, p.231-246, November  2007
gabor wavelets;wavelet transform;shape;feature extraction;subspace mapping;principal components;and texture analysis;multidimensional visualization;feature-based visualization