Incorporating Orientation Selectivity in Wavelet
Transforms for Multiresolution Fourier Analysis of
Images
Andrew Calway
Computer Vision Laboratory,
Linkoping University, S-581 83 Linkoping, Sweden
September 26, 1995
Abstract
The problem of incorporating orientation selectivity into transforms which provide local frequency representation of image regions over a range of spatial scales is investigated. It is shown that this can be achieved if the local spectra are de ned on a log-polar coordinate lattice and that by appropriate choice of window functions, the spectra can be designed to be steerable in arbitrary orientations. In addition, the resulting class of transforms can be de ned to be invertible, be based on window functions having good localisation in both the spatial and spatial frequency domains, and be eciently implemented us-ing FFT techniques. Results of usus-ing one such transform for linear feature extraction demonstrate its e ectiveness when dealing with oriented features.
1 Introduction
For several years now image representations which combine spatial and spatial fre-quency information have been the subject of considerable interest in image processing. This interest has developed along essentially two paths: the use of single resolution forms such as 2-d versions of the short-time Fourier transform (STFT) and Gabor representation [1]-[3]; and the development of multiresolution techniques, culminat-ing in the recent interest in the wavelet transform (WT) [4]-[8]. In both instances, there is a recognition of the advantages provided by local frequency methods in the analysis of images, and, in the case of the latter, a realization that such methods are best employed over a range of spatial scales.
More recently, these approaches have been taken a step further by the introduction of the multiresolution Fourier transform (MFT) [9]-[13]. This is a generalised form of WT in which the properties of the two above mentioned approaches are combined to yield a general multiresolution framework. It achieves this by combining a number of single resolution representations de ned at di erent scales into one hierarchical entity; yielding a structure which contains local spectra referring to di erent sized spatial regions. As a result, the transform can then be seen either as a `stack' of STFTs or a set of `embedded' WTs each de ned over a range of spatial and spatial frequency resolutions. The advantage of this approach is that it combines the completeness of the local spectra provided by a STFT, with the multiresolution properties underlying the traditional form of WT. As was demonstrated in [13], this generalisation enables considerable extension of the applicability of such representations and provides a potentially more versatile framework in which to base suitable image models.
An additional characteristic which is of interest in image processing is that of ori-entation selectivity. The importance of oriented features in analysing images is well known and has considerable support from evidence concerning the workings of the human visual system [14]. Consequently, it has been suggested that as well as em-ploying local frequency methods over di erent scales, suitable image representations should also re ect the orientation properties of local regions. This has led to several investigations into such representations, particularly within the traditional form of the WT [15][16]. In a similar manner, given that the previous implementation of a 2-d MFT was based on a cartesian separable scheme [11], it is pertinent to inquire whether orientation selectivity might also be incorporated within an MFT framework, while still retaining the advantages of generalisation mentioned above.
The purpose of this report is to describe one such class of MFT. Orientation selectivity is incorporated by de ning the local spectra on a log-polar coordinate lattice and the resulting structure resembles a set of WTs each with a given radial and angular resolution. It is consequently more suited to applications involving oriented features and has the added advantage of being steerable in arbitrary orientations if appropriate forms of window functions are adopted. This latter property results in an ability to provide invariant orientation estimation in the case of linear features such as lines and edges. The transforms also retain the invertibilityproperty of their cartesian separable relations, and can be based upon window functions having good localisation in both the spatial and spatial frequency domains. Moreover, by appropriate choice of these functions and spatial sampling parameters, they can be eciently implemented using FFT techniques.
The following section gives a brief overview of the MFT and its extension to two or more dimensions. The de nition of a class of orientation selective MFTs is then presented and examples given for both synthetic and natural images. The report concludes with the results of using such transforms for linear feature extraction.
2 Multiresolution Fourier Analysis
2.1 The Continuous MFT
It is instructive to begin by reviewing the form of the MFT for the case of 1-d continuous signals. Denoting the latter byx(), this is given by [13]
^ x(;!;) = 1 2 Z 1 1 w(( ))x( )expf |! gd (1)
where w() is a window function typically chosen to be real, even, and normalised
such that Z
1
1
w()d > 0 (2)
andw(0) > 0. Readers familiar with the STFT will recognise the similarity with (1); the only di erence being the introduction of the scale parameter, ie for a given value of, ^x(;!;) is a conjoint representation in and !. In fact, closer examination also reveals characteristics of a WT. To see this, note that the MFT coecients correspond to projections of the signal onto a family offrame vectors[7] which transform amongst themselves under translation, dilatation and modulation, ie
^ x(;!;) =hw!;xi (3) where w!( ) = 1 2w(( ))exp f|! g (4)
and h;i denotes the inner product. The above two equations indicate clearly the
relationship between the MFT and its precursors the STFT and WT: the former is derived using frame vectors which transform under a group of symmetries which include as subgroups those associated with the latter, ie the Weyl-Heisenberg and ane groups respectively [7].
It is useful to consider this relationship further in terms of the notion of phase space, the conjoint space de ned by the coordinates (;!). In the case of the STFT, this space can be seen as being linearly represented by translated and modulated versions of some baseband window function; for the WT, the representation is derived from a family of `wavelets'which correspond to translated and dilatated forms of some basic wavelet, leading to a representation which is linear in and logarithmic in ! (also known as thescale space). However, from above, it can be seen that that de ned by the MFT combines both of these characteristics, simultaneously representing the phase space over a range of scales inboth and !. It can therefore be viewed either as a `stack' of STFTs with varying size of baseband window, or `embedded' WTs derived from di ering (ie modulated) basic wavelets. Moreover, the resulting transform space now has an additional dimension, yielding the 3-d space (;!;). Readers should note here that although related, the signi cance of the scale parameter di ers from
that normally associated with the WT; it is now no longer `tied' to the frequency coordinate !. The implications of this amalgamated and consequently overcomplete representation are discussed more fully in [13] and considered here in the context of image analysis in section 2.3.
In common with the STFT and other WTs, the MFT is an invertible transform. The signal x() can be recovered via the synthesis equation
x() = 12Z 1 0 p()d1 2 Z 1 1 Z 1 1 ^ x( ;!;)v(( ))expf|!gd d! (5)
wherev() is a synthesis window satisfying
Z 1
1
w()v( )d = 1 (6)
and p() is a scale density function with Z
1
1
p()d = 1 1 1 (7)
Note in particular that the presence of this density function further emphasizes the `freeing' of the scale parameter from the frequency coordinate; it means that the recovery of the signal from its MFT is based upon coecients de ned over the full range of scales in the phase space.
2.2 Discrete Signals, Windows and Sampling
In practice it is necessary to deal with `real' signals in the form of sequences of samples and a corresponding discrete representation of the phase space. The discrete MFT of the sequencex(k) is then de ned as
^
x(i(n);!j(n);(n)) = X
k wn(k i(n))x(k)exp
f |k!j(n)g (8)
wherek is summed over the set of integers or a nite set thereof (in which case modulo N arithmetic is adopted, 0k < N). The parameters i(n), !j(n), and (n) de ne
points in the 3-d space (;!;) and are each dependent upon the scale index n. The window function wn(k) is similarly scale dependent and is assumed to approximate
scaled versions of the continuous function w(), ie wn(k)(n)
1
2w((n)k) (9)
The main implication of sampling the phase space concerns the invertibility of the resulting transform: it imposes strict conditions upon both the form of the window
function wn(k) and the sampling intervals de ned by i(n), !j(n) and (n). For
example, although in general the inverse is de ned over scale as in (5), in practice it often suces to consider the recovery of x(k) for each value of (n), ie from each
`level' of the transform. In this case, the conditions for exact inversion are those corresponding to that for the discrete STFT [17]. Assuming a regular sampling of the and ! axes, these centre around the limit of the product between the sampling intervals in each domain, ie
(n)!(n)2 (10)
where
(n) = i+1(n) i(n) !(n) = !j+1(n) !j(n) (11)
The required window functions then depend upon the value of the above product: when the equality holds (corresponding to the Nyquist limit) the discrete versions of the frame vectors in (4) for some value of(n) are required to be linearly independent and thus form a complete basis set; as the product is reduced, there is less constraint on the windows (since they move closer together in phase space), resulting in atightor
snug frame [7]. Thus, as would be expected, there exists the usual trade-o between the `optimality'of the representation and the limitations imposed upon the underlying frame vectors.
Up to now the properties of the window functions have been discussed in terms of the invertibility of the transform. However, it is often the case that other desirable qualities are required and the choice of window is then based upon a further trade-o between optimality, invertibility, and the desired form of function. A good example of this concerns the simultaneous locality of the window in both the and ! domains, since the achieved resolution of a discrete representation depends upon the `spread' of the window in the phase space. As is well known, this is limited by the uncertainty principle in signal processing which de nes the so-called unit cells in phase space [7]. Hence, the trade-o then takes the form of nding a window function suciently localised and giving an invertible transform. For example, in [13], the windows were based upon versions of the prolate spheroidal wavefunctions which maximise locality in terms of energy concentration and possess well behaved synthesis counterparts. A related property is the need for a window with nite support, either in the-coordinate (index limited) or !-coordinate (bandlimited), to enable ecient implementation of the transform, either via spatial convolution (eg using a QMF based approach [5][8]) or using FFT methods [13].
2.3 The MFT in Multiple Dimensions
The most straightforward extension of the MFT to multiple dimensions is to use a direct analogy with the 1-d form, ie for the m-d signal x(~) this is given by [13]
^ x(~;~!;) = m2 Z 1 1 w((~ ~))x(~ )expf |~ :~!gd~ (12) 5
where ~, ~ and ~! are now vectors in m-d space and `:' indicates the scalar product. The transform is therefore de ned within a (2m+1)-d continuum which represents the m-d signal x(~) in terms of the coordinates (~;~!;). This is then readily de ned for the discrete case by employing regularly spaced sampling along each of the cartesian coordinate axes as follows
^
x(~i(n);~!j(n);(n)) =X
k wn(~k ~i(n))x(~k)exp
f |~k:~!j(n)g (13)
where ~i(n) and ~!j(n) are now m-d discrete cartesian vectors and the sampling along
each of the axes is dependent on the scale indexn, as are the window functions wn(~k),
which in turn are just the cartesian (tensor) products of the 1-d windows in (9). However, although the simplicity of the above approach is appealing, the increase in dimensionality does allow for considerable scope for adopting alternative structures which, for instance, maybe more suited to some underlying symmetries of a particular application. A good example of this is the importance of oriented features in image analysis and the consequent implicationthat image representations should incorporate some form of orientation selectivity; an example which is of course the motivation behind the present work.
It is appropriate to conclude this section by considering the methodology behind the use of the MFT when applied to image analysis. As discussed in the introduc-tion, multiresolution techniques are now generally accepted as being important for solving analysis problems, albeit a rather loose consensus based on a number of dif-ferent variants. Moreover, the MFT clearly represents a general form of these various approaches. There is however an important distinction to be appreciated between how the latter is used in practice and the previous utilisation of such representations. This derives from the other precursor to the MFT, namely the use of local frequency domain methods in feature extraction. The signi cance of this work, see eg [18][19], was the combination of feature models based on frequency domain characteristics and the assumption that an image consists of local regions which contain single events (such as lines, edges, or texture). Prompted by the success of this work, the MFT was introduced as a means of tackling the problem of scale, ie it is impossible a priori to select a size of region which will t all single event regions. Instead, it is necessary to employ a structure which incorporates local frequency information over a range of scales. The use of such a representation then becomes based on essentially two related operations: the estimation of parameters for some appropriate single event model (using the local spectra of the MFT); and the subsequent selection of the `op-timal' scale (eg (n) in (13)) for each of the assumed nite number of single event regions in the image [13].
Viewed in this way, the requirement for the generality provided by the MFT is unambiguous. Indeed, only in limited cases can other conjoint schemes be said to provide the required representation. Of course, this conclusion depends critically
upon the acceptance that the above constitutes an e ective approach to image anal-ysis, although as demonstrated by the results presented in [13], there is considerable evidence to support such an assertion.
3 An Orientation Selective Transform
The remainder of this report concerns the de nition and application of an orientation selective version of the MFT. By this is meant a class of transform in which the local spectra are de ned on a log-polar coordinate lattice. Such transforms will, for example, explicitly re ect the presence of oriented features in image regions over a range of scales and enable invariant estimation of these features. At the same time, they retain the properties relating to invertibility and good locality possessed by the cartesian separable form.
3.1 The Continuous Case
To get a feeling for how an orientation selective MFT can be constructed, it is useful to begin by considering the continuous case. Given the requirement for a log-polar representation of the frequency domain, it is perhaps not surprising that this is best expressed in a form similar to that for a traditional WT, ie for the image x(~) it is de ned as ^ x(~;;;s) = 1 2 Z 1 1 ws((~ ~))x(~ )d~ (14) where the parameters , , and s de ne the orientation, radial scale and polar fre-quency resolution respectively. The rst thing to note from this equation is that, as previously, the coecients are derived from the inner product of the image with a set of window functions, denoted by ws(~) in this case. In order to incorporate orientation selectivity, these windows are de ned such that they are con ned to some radial and angular portion of the frequency domain, de ned bys and respectively, thus giving a set of coecients which are orientation sensitive (see Fig. 1a). More-over, scaling these functions by then gives a representation of the frequency domain in terms of log-polar coordinates. The relationship with a traditional WT is thus clear: for a given value of s, (14) is simply a WT de ned in each of the orientations [8]. However, the presence of the parameter s enables a variation in the frequency resolution which is characteristic of an MFT. To see this, Fig. 1b shows examples of the frequency response of the window functionsws(~) for some value of and various values ofs. As can be seen, the purpose of the latter is to `scale' the base wavelets on each level of the transform in terms of their radial and angular bandwidths, leading to a consequent increase or decrease in polar resolution. In other words, for each value of s, the set of coecients form a 2-d WT with a speci c resolution in radial and
s (a) (b) ws 1 w s 2 ws 3 i1 i2 wO2 wsO1 i1 i2 O O O1 1 1
Figure 1: Examples of frequency domain support for the window functions ws(~). For some value ofs, the windows have a given radial and angular bandwidth and are oriented in angles as shown in (a). The e ect of s is to scale these bandwidths as shown in (b), leading to a corresponding change in the polar resolution.
angular frequency, a property which is clearly analogous to the cartesian separable case, in which each level corresponds to a STFT with a given resolution in frequency (cf (13)).
3.2 Window Functions and Sampling
Of course, the rst matter to attend to in practice is to de ne a discrete version of (14). This can be de ned as follows
^
x(~i(n;m);k(m);n(m);s(m)) =X
j wknm(~j ~i(n;m))x(~j) (15)
where the discrete windowswknm(~i) are approximations to the set of scaled continuous
functions wknm(~i)n(m)wkm(n(m)~i) (16) for wkm(n(m)~i)w k(m) s(m) (n(m)~i) (17)
and the indices k, n, and m de ne the orientation, radial scaling and polar resolu-tion respectively. Note in particular the dependence of the spatial sampling vectors ~i(n;m) on both n and m. This derives from the fact that the scaling of the windows
m
i1 i2
Figure 2: Frequency domain decomposition corresponding to the discrete orientation selective MFT for various values of m. Each segment represents the ideal region of support for the scaled windowswknm.
by both n(m) and s(m) enable a variation in the spatial sampling rate that can be
employed for di erent values of these two indices. This is covered in greater detail later, although for now it can be appreciated from Fig. 2, which shows examples of the frequency domain decomposition by the scaled windows for di erent values ofm, and by noting the analogy with traditional forms of pyramid representation.
The next matter to consider concerns the form of the window functions wkm(~) in (16). There are of course a numberof choices. For example, in termsof maximisingthe localisation in both domains, adopting suitable versions of either the Gabor functions [1] or the prolate spheroidal wavefunctions [20] would be possible options. Indeed, the latter have been used previously for linear feature estimation [21]. Alternatively, one could seek some form of QMF approach along the lines of recent WTs such as adopted in [15] or [8]. However, although these examples have a number of properties to recommend them, the functions adopted here are based on alternative criteria: namely that they should be polar separable in the frequency domain and de ned such that the set of scaled functions for a given value ofm sum to a constant value in the radial and angular directions. The advantages of this analytic approach are twofold: rst, it enables a class of invertible transforms to be readily de ned; and secondly, it allows the interpolation of intermediate coecients in arbitrary orientations with negligible error (steerability), leading to an ability to provide invariant estimation of linear feature orientation. Moreover, although maximum simultaneous localisation is no longer a prerequisite, it is possible to choose such functions so that they possess good localisation in both the spatial and spatial frequency domains.
The functions used here are related to those used successfully for orientation estimation by Knutsson [19]. They are de ned in the continuous frequency domain as
^
wkm(~!) = rm()qkm() m > 0 (18)
where the coordinates ~! are given by
~! = (cos; sin) (19)
for > 0 and 0 < 2. The function rm() determines the radial response of the
window rm() = ( cos2 ln(=(m)) 2Arlns(m) jln(=(m))j< Arlns(m) 0 else s(m) > 0 (20)
where(m) > 0 de nes the inner radial frequency on level m, and qkm() is the angular response qkm() = ( cos2 ( k(m)) 2Aq(m) j k(m)j< Aq(m) 0 else (21) where k(m) = (k + 1=2)(m) 0k < K(m) (22)
and(m) = 2=K(m). The parameters Ar > 0 and Aq> 0 are integers and determine
the radial logarithmic and angular bandwidths of the windows in terms of multiples of s(m) and (m) respectively. Thus, in the latter case, Aqde nes the angular overlap of
the windows as shown in Fig. 3. In a similar manner, if the discrete scaling parameter n(m) in (15) is de ned as
n(m) = s(m)n 1n L(m) (23)
then Ar determines the overlap in the radial direction as shown in Fig. 4.
Recall that a criterion for the windows was that for a given value ofm they should sum to a constant value in both the angular and radial directions. That this is so for the above functions can be readily derived from (20) and (21), and is clearly apparent from Figs. 3 and 4. Speci cally,
L(m) X n=1 rm(=n(m)) = Ar (m)L (m)(m)(m) (24) and K (m) 1 X k=0 qkm() = Aq (m)=2 i 2 (m)=2 (25) 10
0 0.5 1 0 0.5 1 1.5 2 2.5 3 (a) 0 0.5 1 0 0.5 1 1.5 2 2.5 3 (b)
Figure 3: Examples of the angular functionqkm() for (m) = =8, 1k < 8 and (a)
Aq= 1 (b) Aq = 2. 0 0.5 1 0 0.5 1 1.5 2 2.5 3 (a) 0 0.5 1 0 0.5 1 1.5 2 2.5 3 (b)
Figure 4: Examples of the radial functionrm() for s(m) = 1:2, (m) = =8 and (a)
Ar= 1 (b) Ar = 2.
Moreover, from these equations, it is then clear that the set of 2-d functions ^wkm(~!=n(m))
also sum to a constant within some radial frequency band
L(m) X n=1 K(m) 1 X k=0 ^ wkm(~!=n(m)) = ArAq (m) L (m)(m)(m) (26)
Examples of 2-d windows functions based on the above equations are presented in section 3.3.
The nal matter to consider concerns the form of the spatial sampling vectors ~i(n;m) in (15). The most straightforward approach is to de ne these on a standard
rectangular lattice, ie using integer arithmetic and assuming a nite lattice of size N1(n;m) N 2(n;m) ~i(n;m) = (i0 1(n;m); (i=N1(n;m))2(n;m)) 0 i < N 12(n;m) (27) where i0 =i mod N 1(n;m), N12(l;n) = N1(l;n)N2(l;n) and 1(n;m) and 2(n;m)
de ne the sampling intervals along each coordinate axis. In order to decide upon suitable values for the latter, it is useful to consider the e ect of the spatial sampling in the frequency domain. Towards this, recall from (15) that the subset of transform coecients for given values of k, n, and m correspond to the sampled output of a lter with impulse responsewknm(n(m)~j), ie the transform as a whole can be viewed
as a number of lterbanks de ned at di erent scales. Noting that sampling in the spatial domain implies a `folding' of the 2-d spectrum [22], the e ect of the former can be viewed as in Fig. 5, which shows the bounded region of sizeN1(n;m)
N
2(n;m)
into which the spectrum folds as a result of sampling at intervals of 1(n;m) and
2(n;m) along each axis. For simplicity, and without loss of generality, this is shown
centred about the respective lter's centre frequency.
From the above, an initial condition for determining suitable values for the sam-pling intervals can then be based on the requirement for an invertible transform. If the rectangular region in Fig. 5 is de ned by the set of discrete coordinates f~!j jj 2 R(k;n;m)g and the window functions in (16) are chosen such that
^
wknm(~!j) = 0 j 62 R(k;n;m) (28)
then from the 2-d sampling theorem it can be shown that the image x(~j) can be
recovered for some value ofm using x(~j) = xLP(~j;m) + K (m) 1 X k=0 L(m) X n=1 xkn(~j;m) + xHP(~j;m) (29) where xkn(~j;m) =N 12(n;m) 1 X i=0 ^ x(~i(n;m);k(m);n(m);s(m))vnm(~j ~0 i(n;m)) (30) 12
1 2 ( , ) i1 i2 x( , ) i2 2 i1 1 x( , ) i1 i2 i1 i2 N N 1 2
Figure 5: Subsampling the 2-d function x(i1;i2) at intervals of 1 and 2 along
each axis results in a `folding' of the spectrum into a bounded region of sizeN1 N 2, whereN11 =N22 =N and 0 i 1;i2 < N. for ~0 i(n;m) = (1(n;m)i1(n;m); 2(n;m)i2(n;m)) (31)
and synthesis windowsvnm(~j) given by
vnm(~j) = sin(j1N1(n;m)) j1N1(n;m) ! sin(j2N2(n;m)) j2N2(n;m) ! (32) The functions xLP(~j;m) and xHP(~j;m) in (29) are lowpass and highpass `residual'
images determined by the radial frequency extent of the set of windows on the level in question, ie they are de ned in the frequency domain as
^ xLP(~!ij;m) = ( (1 rm(j)) ^x(~!ij) j (m) 0 else (33) and ^ xHP(~!ij;m) = ( 1 rm(L(m)(m)j) ^ x(~!ij) j L (m)(m)(m) 0 else (34)
where the frequency coordinates~!ij are discrete versions of those de ned in (19) ~!ij = (jcosi; jsini) (35)
The validity of (29) is best appreciated by noting that the condition in (28) ensures that the transform coecients are not aliased as a result of the spatial sampling and thus each bandpass subimage can be reconstructed by interpolation via convolution with the 2-d sinc function in (32). The summation properties of the windows then enables the original image to be reconstructed by summationof these bandpass images and the addition of the residualsxLP(~j;m) and xHP(~j;m), the latter replacing the
lowpass and highpass regions of the frequency domain not covered by the window functions (cf (26)).
In terms of minimising redundancy, it is appropriate to choose the spatial sam-pling intervals so as to minimise the size of the folding rectangle in Fig. 5 without violating the inversion condition in (28). However, if based solely on this criterion, the resulting transform would be rather unstructured, in that it would be dicult to relate the range of sampling densities that could be employed. It is likely that such a transform would be awkward to use and in any case would not t into the `local spectra' framework which motivates an MFT approach (cf section 2.3). The latter point might suggest the adoption of the same sampling interval for a given value ofm, thus ensuring a complete representation of the frequency domain at some spatial res-olution on each level. However, this ignores the di ering extent of the windows, both in terms of the size of the spatial regions referred to and the resultant redundancy. As a compromise, the approach adopted here is to choose the sampling densities to be the same along each axis and to be a power of 2, the appropriate power being selected to minimise redundancy and ensure agreement with (28). The result is a transform in which each level has a nite number of coecient subsets, corresponding to dif-ferent sampling densities, which have a well de ned relationship amongst themselves and with other subsets on other levels. In addition, the scheme has the advantage of enabling ecient calculation of the transform coecients via FFT methods. To be more precise, the sampling lattices are de ned such that
N1(n;m) = N2(n;m) = N(n;m) (36)
where
N(n;m) = 2 (n;m)
(n;m) > 0 (37)
and (n;m) is an integer chosen to minimise the number of coordinates in the region
R(k;n;m) without violating (28). As examples presented in the following sections
will demonstrate, this approach enables an appropriate trade-o between minimising redundancy and maintaining a usable structure. Moreover, since the window functions are bandlimited and the number of spatial samples are constrained to be a power of 2, each of the bandpass ltering operations required to calculate the transform coecients can be implemented using a separable radix-2 FFT algorithm in a similar manner to that adopted for the cartesian separable case [11].
3.3 An Example
The purpose of this section is to present an example of an orientation selective MFT based on the class of transforms described above. Suitable values for the various parameters involved are discussed and examples given of the associated 2-d window functions. Examples of the transform are presented for both synthetic and natural images.
It is rst necessary to decide upon appropriate values for the parameters de ning the form of the window functions ^wkm(~!) given by (18). These are determined by two criteria: the required aspect ratio of the spectral envelopes; and the degree of overlap between adjacent envelopes. In the case of the former, the approach adopted here is to de ne the radial and angular bandwidths such that resulting windows are approximately circular in the 2-d plane, ie having an aspect ratio 1. This
is motivated by the `local spectra' methodology underlying the MFT described in section 2.3: to be useful, the spatial response of the windows should coincide as much as possible, implying the use of isotropic (ie circular) envelopes which are independent of orientation1. The factors e ecting the choice of frequency overlap are discussed in
the next section.
To obtain approximate circularity of the windows, the radial and angular band-widths are de ned to be equal at the 3dB point. This corresponds to selecting the resolution parameterss(m) and (m) such that (see Appendix)
Arlns(m) = Aq(m) (38)
where Ar and Aq are chosen according to the required degree of overlap. Plots in
both the spatial and frequency domains for window functions satisfying the above equation are shown in Figs. 6 and 7 forAr=Aq = 1 and (m) and s(m) given by
(m) = 2=K(m) K(m) = 2m+2 (39)
s(m + 1) = s(m)1
2 (40)
Frequency responses of examples of sets of scaled windows are shown in Figs. 8 and 9 for 1m2.
The spatial sampling densities adopted for the above windows are summarised in Table 1. These were selected on the basis of using three densities atm = 1, N=2, N=4, andN=8 respectively for each of the three radial frequencies (the same density being used for all orientations) and an image of sizeN N pixels. Subsequent levels then
have densities corresponding to a reduction by a factor of 2 from the previous level
1This is in contrast with some previous orientation selective schemes, see eg [1], in which the
spectral envelopes of the windows are themselves de ned to be sensitive to a particular orientation, a fact which perhaps further illustrates the di erence in methodology adopted.
-10 0 10 0 -10 0 10 0 -2 0 2 0 -2 0 2 0 (a) (b) (c) (d)
Figure 6: Examples of 2-d window functions wknm(~) satisfying (38) and with Ar =
Aq = 1: (a) Spatial envelope with (m) = =4, s(m) = 2:2, k = 1 and n = 2;
(b) Frequency envelope corresponding to (a); (c) Spatial envelope with(m) = =8, s(m) = 1:48, k = 6 and n = 3; (d) Frequency envelope corresponding to (c).
(c) (d)
(a) (b)
Figure 7: 3-d plots of the windows shown in Fig. 6. 16
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
Figure 8: Set of scaled window functions form = 1 (dotted lines are used alternatively to aid visability) -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3
Figure 9: Set of scaled window functions for m = 2. 17
Level Band 1 (LP) Band 2 Band 3 (HP) m=1 N 8 N 4 N 2 m=2 N 16 N 16 N 8 N 8 N 4 N 4 m=3 N 32 N 32 N 32 N 32 N 16 N 16 N 16 N 16 N 8 N 8 N 8 N 8
Table 1: Sampling densities along each spatial coordinate and within each frequency band for 1 m 3 and image of size N N. The same density is used for all
orientations.
within each of these `base' bands. Appropriate choice of the inner radial frequency (m) on each level such that the window condition in (28) is not violated then ensures an invertible transform.
The transform described above was derived for the `FM Pattern' and `Lena' 256 256 8-bit grey level images shown in Fig. 13. The former is a symmetric
logarithmic FM pattern consisting of frequencies approximately in the range =10 to and its transform coecients for 1 m 3 are shown in Figs. 14-16. These are
displayed as bandpass subimages arranged according to orientation selectivity and radial frequency. This example clearly illustrates both the frequency domain proper-ties of the windows and the sampling scheme detailed in Table 1. The highpass and lowpass residual images,xLP(~j;m) and xHP(~j;m), are shown at the top and
bot-tom respectively of each frequency-orientation display. The coecients corresponding to the `Lena' image are shown in Figs. 17-19. This example illustrates clearly the orientation selectivity of the transform and the change in resolution as a function of m.
3.4 Local Region Analysis and Steerability
As outlined in section 2.3, the use of the MFT is based on the notion of analysing image regions using locally de ned spectra at di erent resolutions. In the case of the cartesian separable form, this is a natural framework in which to use the transform, since each level consists of spectra de ned at some given spatial resolution and thus processing can be based on a level by level approach [13]. However, the situation is less clear for the orientation selective form described above; in this case, each `level' contains of several representations of concentric frequency bands de ned at di erent
m=1 m=2 m=3
composite level p=2
Figure 10: Ideal regions of support for window functions forming a composite level from three adjacent levels.
spatial resolutions, giving a part representation of the frequency domain for di erent sized spatial regions. In order to utilise the latter within the MFT `local spectra' framework, it is necessary to combine coecients from di erent levels in such a way that each of the resulting composite levels provide a set of complete local spectra de ned at the same spatial resolution. These can then be used in much the same way as their cartesian separable counterparts.
Each of the composite levels are formed by combining sets of coecients which have the same spatial sampling density. This is particularly straightforward if the sampling within the transform is based on the scheme described above and results in a representation de ned over the whole of the frequency domain. This can be seen from Fig. 10, which shows the spectral regions of the windows forming a given composite level, and by reference to the sampling details in Table 1. However, although such levels clearly provide the required representational structure in terms of sampling density and frequency coverage, there are two signi cant issues that need to be noted: rst, the log-polar frequency de nition of the windows means that the spatial regions referred to by the coecients within each band on a composite level will be of di ering extent; and secondly, as is readily apparent from Fig. 10, the orientation selectivity of the bands will be disparate and de ned at di erent resolutions.
Somewhat paradoxically, the rst of these issues relates directly to the requirement for orientation selectivity within the transform and the adoption of a log-polar coor-dinate lattice. However, if the frequency separation intervals and window functions are chosen appropriately, it is possible to ensure that there is reasonably good
0 0.5 1 0 1 2 3 4 5 6 0 0.5 1 0 1 2 3 4 5 6 (b) (a)
Figure 11: Interpolated angular function in an orientation 45o derived from sets of 4
and 8 basis functionsqkm() (m = 1 and m = 2).
spondence between the spatial envelopes on each composite level. The di erence in orientation selectivity of the coecients can in contrast be overcome rather e ectively by designing the transform to besteerable, ie to choose the window functions such that it is possible to interpolate new coecients in the required orientations within each of the composite level frequency bands (eg corresponding to the preferred orientations of the outer band). In the present case this can be achieved with negligible error in the interpolation by selecting a suitable angular overlap of the windows in the frequency domain via the parameterAq. To illustrate this, Fig. 11 shows a comparison between
interpolated and required versions the angular functionqkm(i) for intermediate
orien-tations derived from the set of `basis' functions used in the transform forAq = 2. The
interpolation is based on a simple least-squares approximation [25, 26]. As can be seen, there is negligible error between the required functions and their interpolated counterparts. Moreover, given the polar separability of the windows, it is therefore possible to derive appropriate transform coecients in arbitrary orientations with the same degree of error. In other words, one can de ne a composite level with index p given by ^ xc(~i(p);~!j(p);p) =K (m) 1 X k=0 ak(;n;m) Re n ^ x(~i(n;m);k(m);n(m);s(m)) o + 20
0 0.5 1 0 1 2 3 4 5 6 -1 0 1 0 1 2 3 4 5 6 (b) (a)
Figure 12: Even and odd parts of the interpolated and basis functions in Fig. 11b. |K(m) 1 X k=0 bk(;n;m) Im n ^ x(~i(n;m);k(m);n(m);s(m)) o (41) for all n and m such that N(n;m) = N(p), where N(p) N(p) is the number of
spatial samples on the composite level. The interpolation coecientsak(;n;m) and
bk(;n;m) are those corresponding to the interpolated quadrature parts (even and
odd) of a function with frequency response given byqm( ) as shown in Figs. 12a and 12b, ie in the spatial domain the polar separable windows are complex valued with the real and imaginary parts forming a quadrature lter pair [27]. Note in particular that the spectra on the levels de ned in (41) are now de ned in terms of the polar coordinate vectors
~!j(p) = (j(p)cos; j(p)sin) (42)
which are continuous in and where the radial frequencies j are de ned for some
L(p) such that 1(p) < 2(p) < ::: < L(p)(p) and
fj(p)g =f(n;m)g (43)
for (n;m) = n(m)(m), ie the ordered set of radial frequencies corresponding to
the coecients that constitute the composite level. Hence, by choosing Aq = 2 and
a corresponding value ofAr given by (38), a transform can be derived which consists
of local spectra de ned over a range of spatial resolutions and which are steerable in 21
any required orientation. This latter property is clearly related to the recent interest in steerable representations, see eg [23][24], and it is shown in the next section how this property can be utilised in the estimation of linear feature orientation. It is interesting to note moreover, that although not used in the present work, the above composite levels could also be designed to be continuous in radial frequency by an appropriate choice ofAr and using the same interpolation scheme as employed in the
angular direction.
4 Application: Linear Feature Extraction
The transform described in the preceding sections was applied to the problem of extracting linear features, such as line and edge segments, from images. The scheme illustrates the general methodology behind the use of the transform and the results demonstrate its orientation selective properties. Experiments were performed using the orientation and position estimation algorithm described in [13] which is based on a linear phase model of `locally 1-d' regions. For completeness, brief details of the algorithm are given here.
4.1 A Maximum Likelihood Estimator
The basis of the algorithm is the observation that the Fourier transform of a 2-d continuous signal varying only in an orientation
x(~) = x(1cos + 2sin) (44)
is given by
^
x(~!) = (!1sin !2cos)^x(!1cos + !2sin) (45)
where ^x(!) is complex function de ned as ^
x(!) = A(!)expf |! + g (46)
and is the centroid of the 1-d function x() and is a phase constant. In other words, the spectrum of such functions, which correspond to ideal line and edge features for example, is con ned to a line pasing through the origin at an angle perpendicular to that in the spatial domain and having a linear phase variation which is determined by the relative `position' of the feature.
The model underlying the estimation scheme is then based on this linear phase relationship. Ignoring the e ects of any local analysis, it takes the form of assuming that a Fourier representation of some image region containing a nite number of such
features can be written in polar coordinates as ^ x(;) =K 1 X k=0 ^ xk()( k) (47)
where ^xk() is of the form in (46) and de nes the variation along the orientation k
due to thekth feature. Moreover, if these functions are sampled at regular intervals j =j=L, 1 j L, then the samples are assumed to correspond to those from a
normal Markov process such that ^ xk(j) = kexpf| kgx^k(j 1) + kvk(j) j > 1 (48) and ^ xk(1) = kvk(1) (49) where 2 k = (1 2
k) and vk(j) are complex innovations which are normally
dis-tributed with zero mean and unit variance. This clearly incorporates the linear phase relationship into the model since if a feature is present in an orientationk then
k = (j j 1)k (50)
wherek is the spatial centroid of the feature.
The purpose of the estimator is to estimate values for the complex coecients kexpf| kg. Assuming that suitable frequency samples are available, then it is
straightforward to show that a maximum likelihood estimation can be derived from the correlation statistic
Rk = 1L 1 L 1 X j=1 ^ x k(j)^xk(j 1) (51)
where indicates complex conjugate and the resulting estimates are unbiased since
E[Rk j k = 0] = 0 (52)
and
E[Rk j k > 0] = kexpf| kg (53)
The required estimate of k is therefore
~ k = argfRkg (54)
and the magnitude, ~ k =jRkj, represents a certainty measure based on the projection
of the samples ^xk(j) onto the model in (48).
In practice, the estimator is implemented by deriving the samples ^xk(j) in a nite
number of orientations from an appropriate local Fourier representation of a given 23
image region, eg from those within an MFT of the image. The correlation statistics in each of these orientations then provide a description of the region in terms of linear features. This description can then be further analysed using for example some form of principal orientation algorithm to determine the dominant feature in the region. Of course, such an estimation is necessarily an approximation, since implementing any form of local analysis automatically impliesthe use of a discrete Fourier representation and a consequent `smoothing' of the ideal model in (47). Nevertheless, as shown in the next section, if these local representations are derived without excessive bias, then acceptable estimates can be obtained at the appropriate resolutions.
4.2 Implementation and Results
The orientation selective form of the MFT described in this report clearly provides an ideal framework in which to implement the above estimator. Using the composite levels de ned in section 3.4, a set of correlation statistics are derived from each local spectrum in a nite numberof orientations, leading to position and certainty estimates for features in those orientations within the regions referred to by the spectra. Since the latter are de ned on a polar coordinate lattice, the statistics can be easily derived directly from the transform coecients. A subsequent principal orientation analysis then provides estimates of the orientation of the dominant feature(s) within each of the local regions and at the various spatial resolutions.
It is however rst necessary to adapt the de nition of the correlation statistics in (51) to the case in which the radial frequency samples are irregularly spaced, as in the composite level de nition of (43). This takes the following form
R(k;i;p) = 1L(p)L(p) 1 X j=1
jRj(k;i;p)jexpf | j(p)argfRj(k;i;p)gg (55)
where
Rj(k;i;p) = ^x
c(~i(p);!jk;p)^xc(~i(p);!jk;p) (56)
and ^xc(~i(p);!jk;p) are the composite coecients in an orientation k = k=K(p),
0 k < K(p). The irregular frequency sampling is accounted for by the phase
correction term j(p) which is de ned as
j(p) = 2(p) 1(p)
j+1 j(p)
(57) and can be seen as aligning the phases to that ofR1(k;i;p), thus enabling meaningful
averaging via the summation overj in (55). Estimation of the feature centroids with respect to the spatial sampling points ~i(p) can then be derived using (50) and (54)
~
(k;i;p) = argfR(k;i;p)g
2(p) 1(p)
(58) 24
and corresponding certainty measures given by ~
(k;i;p) =jR(k;i;p)j (59)
for each of the spatial resolutions de ned by the indexp. Results of implementing the scheme for the Rose, Lena and Baboon images are shown in Figs. 20-22. These show the feature estimates displayed as graphics generated line segments in the appropriate spatial regions and with orientation, position and weighting de ned by k, ~(k;i;p)
and ~ (k;i;p) respectively. Clearly all the linear features have been detected at the appropriate resolutions for the Rose image and the majority detected for the case of the natural image.
To derive estimates for the dominant feature within a region, a tensor based principle orientation sceme was used in a similar manner to that in [13]. This is based on the representational framework proposed by Knutsson [28] and consists of de ning the following inertia tensor for each of the local regions
T(i;p) = X
k (k;i;p)~u(k)~ (60)
where~u(k) is a unit vector in the orientation k. It is readily shown that the
eigenvec-tors ofT(i;p) de ne the major and minor axis of an ellipse which is oriented according
to the principle orientation of the distribution de ned by the ~ (k;i;p) and whose as-pect ratio represents the anisotropy associated with the distribution. In other words, an eigenvalue analysis ofT(i;p) yields an estimate of the dominant orientation within
a given region and a measure of the degree of dominance. The centroid estimate in this orientation can then be obtained by interpolating the appropraite coecients via (41), thus giving an invariant estimate of the dominant feature.
The above scheme is readily extended to the case of more than one signi cant feature within a region, providing they have suciently di erent orientations with respect to the current angular resolution. This is achieved by applying a simple local maximum detection operation to the distribution ~ (k;i;p) and then using a `local' tensor analysis about each maximum to determine the orientation of each signi cant feature. More precisely, for each maximumatk = k1;k2;:::;kS, a local inertia tensor
is derived according to Ts(i;p) = X k2Ks ~ (k;i;p)~u(k) (61)
whereKs is the set of orientations about ks such that
~ (k + 1;i;p) ~ (k;i;p) ( > 8 k < ks < 8 k ks 8 k 2Ks (62)
for some suitable value of. An eigenvalue analysis of each tensor then yields a nite number of `local' principal orientations s, which can then be used to derive centroid
and position estimates for each signi cant feature in the region. 25
Results of implementing the principal orientation analysis on the linear feature estimates derived for the Rose, Lena and Baboon images are shown in Figs. 23-26. The rst of these shows the result of estimating only the dominant principal orientation within each region for the Rose image. Clearly at the junction points this single feature hypothesis is inappropriate, resulting in very low certainty measures in those regions. In contrast, Fig. 24 shows the estimates derived when two dominant orientations are allowed, ie the two largest local maxima in the certainty distribution are used. The junction points are then well represented by this multiple feature hypothesis. Similar results for the other two images are shown in Figs. 25 and 26.
5 Conclusions
A class of transforms which provide orientation selective frequency representation of local regions over a range of spatial resolutions has been presented. They represent a natural adaptation of the cartesian separable MFT described in [13][11] to the case in which oriented features in images are of particular interest. Moreover, the win-dow functions adopted means that they can be de ned to be invertible and steerable in arbitrary orientations. The results of applying one such transform to linear fea-ture extraction demonstrated its orientation selective properties and these compare favourable with the results achieved in [13].
Appendix
To achieve approximate circularity of the window functions, the radial and angular bandwidths are chosen to be equal at the 3dB point, ie when rm() = qkm() = 1=2.
This implies from (20) and (21) that ln(2=(m)) 2Arlns(m) = ln( 1=(m)) 2Arlns(m) = 4 (63) and ( 2 k(m)) 2Aq(m) = ( 1 k(m)) 2Aq(m) = 4 (64)
where = 2 1 > 0 and = 2 1 > 0 are the 3dB radial and angular
bandwidths respectively. Assuming is small, the circularity condition can then be expressed as
(m) (65)
If is also assumed small, then 2
(m) + =2 (66)
and from (63) (m) + =2 = (m)expfArlns(m)=2g (67) giving = 2(m)(expfArlns(m)=2g 1) (68) From (64) = Aq(m) (69)
and then from (65) and (68)
Aq(m) = 2(expfArlns(m)=2g 1) (70)
which becomes
Arlns(m)=2 = ln(Aq(m)=2 + 1) (71)
and nally for small(m)
Arlns(m) Aq(m) (72)
Acknowledgements
The author wishes to thank Dr. Hans Knutsson and Dr. Roland Wilson, both of whom have inspired this work and made signi cant contribution to its content. Thanks are also due to membersof the Computer Vision Laboratory, Linkoping, and in particular Leif Haglund for his help with the window interpolation. The work was supported in part by the Royal Society and the Royal Swedish Academy of Sciences.
References
[1] J.G.Daugman, \Complete discrete 2-d Gabor transforms by neural networks for image analysis and compression,"IEEE Trans. Acoust., Speech, Signal Process-ing, vol. ASSP-36, pp. 1169-1179, 1988.
[2] A.C. Bovik, M. Clark, and W.S. Geisler, \Multichannel texture analysis using localized spatial lters", IEEE Trans. Patt. Anal. Mach. Intell., vol. 12, pp. 55-73, 1990.
[3] M. Porat and Y.Y.Zeevi, \The generalizedGabor schemeof image representation in biological and machine vision", IEEE Trans. Pattern Anal. Machine Intell., vol. PAMI-10, pp. 452-68, 1988.
[4] P.J. Burt and E.H. Adelson, \The Laplacian pyramid as a compact image code ",IEEE Trans. Comun., vol. COM-31, pp. 532-40, 1983.
[5] E.H. Adelson, E. Simoncelli, and R. Hingorani, \Orthogonal pyramid transforms for image coding", Proc. SPIE Colloq. on Vis. Commun. and IP, Cambridge, MA, 1987.
[6] J-B. Martens, \The Hermite transform", IEEE Trans. Acous. Speech Sig. Proc., vol. ASSP-38, pp. 1595-1618, 1990.
[7] I. Daubechies, \The wavelet transform, time-frequency localization and signal analysis", IEEE Trans. Inform. Th., vol. IT-36, pp. 961-1005, 1990.
[8] S. Mallat, \Multifrequency channel decompositions of images and wavelet mod-els", IEEE Trans. Acoust., Speech, Signal Proc., vol. ASSP-37, pp. 2091-2110, 1989.
[9] R. Wilson and A.D. Calway, \A general multiresolution signal descriptor and its application to image analysis", Proc. EUSIPCO-88, pp. 663-666, Grenoble, 1988.
[10] A.D. Calway and R. Wilson, \A uni ed approach to feature extraction based on an invertible image transform", Proc. 3rd IEE Int. Conf. Image Processing, Warwick, 1989.
[11] A.D. Calway,The Multiresolution Fourier Transform: A General Purpose Tool for Image Analysis, Ph.D. Thesis, Warwick Univ., 1989.
[12] A.D. Calway and R. Wilson, \Curve extraction in images using the multiresolu-tion Fourier transform", Proc. IEEE ICASSP-90, pp. 2129-2132, Albuquerque, 1990.
[13] R.Wilson, A.D.Calway, and E.R.S.Pearson, \A generalized wavelet transform for Fourier analysis: the multiresolution Fourier transform and its application to image and audio signal analysis", in submission.
[14] D.H. Hubel, Eye, Brain, and Vision, New York: Sci. Amer. Lib., 1988.
[15] E. Simoncelli and E.H. Adelson, \Non-separable extensions of quadrature mirror lters to multiple dimensions",Proc. IEEE, vol. 78, pp. 652-664, 1990.
[16] A.B. Watson, \The Cortex transform: Rapid computation of simulated neural images",Comput. Vision, Graph., Image Proc., vol. 39, pp. 311-327, 1987. [17] M.R. Portno , \Time-frequency representation of digital signals and systems
based on short-time Fourier analysis", IEEE Trans. Acous. Speech Sig. Proc., vol. ASSP-28, pp. 55-69, 1980.
[18] G.H. Granlund, \In search of a general picture processing operator", Comput. Graph. Image Proc., vol. 8, pp. 155-173, 1978.
[19] H. Knutsson, Filtering and Reconstruction in Image Processing, Ph.D. Thesis, Linkoping Univ., 1982.
[20] R. Wilson, \Finite prolate spheroidal sequences and their applications, pt I",
IEEE Trans. Patt. Anal. Machine Intell., vol. PAMI-9, pp. 787-795, 1987. [21] R. Wilson, \Uncertainty, eigenvalue problems and lter design", Int. Report
LiTH-ISY-I-0580, Linkoping Univ., 1983.
[22] A.K.Jain,Fundamentals of Digital Image Processing, Englewood Cli s: Prentice-Hall, 1989.
[23] W.T.Freeman and E.H.Adelson, \The design and use of steerable lters for image analysis, enhancement, and wavelet representation", Tech. Report 126a, MIT, Cambridge, Mass., 1990.
[24] M.T.Andersson and H.Knutsson, \Orientation estimation in ambiguous neigh-bourhoods", 7th Scan. Conf. Image Anal, Aalborg, Denmark, 1991.
[25] W.H.Press, B.P.Flannery, S.A.Teukolsky, and W.T.Vetterling, Numerical Recipes in C, CUP, 1988.
[26] L. Haglund, \Hierarchical scale analysis of images using phase description", Link0ping University Tech Report LIU-TEK-LIC-1989:08, Linkoping, 1989. [27] A.Papoulis,Signal Analysis, New York: McGraw-Hill, 1977.
[28] H. Knutsson, \Representing local structure using tensors",Proc. 6th Scan. Conf. on Image Analysis, Oulu, 1989.
Figure 13: Original 2562 pixel 8-bit grey level images: `FM Pattern'; `Rose'; `Lena';
`Baboon'.
Figure 14: Transform coecients for `FM Pattern'm = 1. Horizontal subimages cor-respond to di erent orientations and vertical subimages to di erent radial frequencies. Highpass and lowpass residual images are shown at top and bottom respectively.
Figure 15: Transform coecients for `FM Pattern' m = 2.
Figure 16: Transform coecients for `FM Pattern' m = 3. 31
Figure 17: Transform coecients for `Lena' imagem = 1.
Figure 18: Transform coecients for `Lena' imagem = 2. 32
Figure 19: Transform coecients for `Lena' imagem = 3.
Figure 20: Linear feature estimates for `Rose' image: ulp = 1; ur p = 2; ll p = 3. 33
Figure 21: Linear feature estimates for `Lena' image.
Figure 22: Linear feature estimates for `Baboon' image.
Figure 23: Dominant principal orientation estimates for `Rose' image derived from linear feature estimates in Fig. 20.
Figure 24: Two dominant principal orientation estimates for `Rose' image derived from linear feature estimates in Fig. 20.
Figure 25: Two dominant principal orientation estimates for `Lena' image.
Figure 26: Two dominant principal orientation estimates for `Baboon' image. 36
