in: M. van Ginkel, P.W. Verbeek, and L.J. van Vliet, Multi-orientation estimation: Selectivity and localization, in: H.E. Bal, H. Corporaal, P.P. Jonker, J.F.M. Tonino (eds.), ASCI’97, Proc. 3rd Annual Conference of the Advanced School for Computing and Imaging
(Heijen, NL, June 2-4), ASCI, Delft, 1997, 99-105.
Multi-Orientation Estimation:
Selectivity and Localization
M. van Ginkel, P.W. Verbeek, L.J. van Vliet
Pattern Recognition Group
Faculty of Applied Physics,
Delft University of Technology,
Lorentzweg 1, 2628 CJ Delft, The Netherlands,
{michael,piet,lucas}@ph.tn.tudelft.nl
Keywords: lines, edges, orientation, texture
Abstract
Filtering of an image with rotated versions of an orientation selective filter yields a set of im-ages which can be stacked to form an orienta-tion space. Orientaorienta-tion space provides a means of analyzing overlapping and touching anisotropic textures. A set of rotated kth order directional derivatives yields a discrete orientation space, which allows interpolation. Next we apply a de-convolution scheme that results in improved ori-entation selectivity. This scheme allows decom-position of noisy multi-orientation patterns.
1
Introduction
Images are composed of various types of struc-tures, such as lines, edges and textures. These in turn can be characterized by various proper-ties, the most important being intensity, scale and orientation. In this paper we focus on improved orientation selectivity for multi-orientation esti-mation.
An image can often locally be modeled as a weighted sum of translation invariant patterns or (paintbrush) strokes. Each stroke has a one-dimensional intensity profile and a typical orien-tation: the profile orientation across the stroke, see figure 1.
We define directions as angles in the interval (0, 2π), thus making a distinction between an ar-row pointing to the left or right. When we refer to orientation (0, π) we make no such distinction. For many applications a reliable estimator for a single locally dominant orientation suffices, Kass & Witkin [4], Haglund [3] and Van Vliet & Ver-beek [9]. There are a number of applications for
x φ
Figure 1: An oriented pattern
which an estimate of the local orientation will not be sufficient. Textures can often be characterized by a number of overlapping patterns with a dif-ferent orientation. At boundaries between two single orientation regions, an estimator for a lo-cally dominant orientation will give the wrong an-swer. In such cases a more advanced analysis is required.
For another important property, scale, an ex-tensive framework has been created over the last decade, initiated by Witkin [12] and Koen-derink [6]. In the scale space paradigm, scale is explicitly dealt with by embedding an image in a new image with one extra dimension, the scale dimension. This image is generated from the original by iteratively applying an isotropic (e.g. Gaussian) filter.
In order to deal with anisotropy an image or scale space can be replaced by a stack of direction-ally filtered versions of the original. The resulting orientation space or orientation+scale space is pe-riodic along the orientation axis with a period of π. Orientation space and scale space differ in two aspects; the original image is part of scale space, but not of orientation space; scale space has a hierarchical structure, whereas orientation space does not. As shown in Van Vliet and Verbeek [10], such an orientation space can also be used to
per-form a segmentation of overlapping objects. The main objective of this paper is to establish how to generate an orientation space representation of an image with a high orientation selectivity.
Interesting related work on orientation estima-tion and juncestima-tion classificaestima-tion has been done by Andersson [1] and Michaelis [7].
2
Constructing an orientation
space
Orientation space is created by applying ro-tated versions of some orientation selective filter to the image. Initially the only constraint on the choice of filter is that it allows proper sampling of orientation space. The filters we use are all two dimensional and operate in the (x, y) plane. Yet, most of the ensuing (Fourier) analysis takes place in the orientation dimension and will be one-dimensional.
2.1
Directional derivatives
We will consider the possibility of using a Gaussian directional derivative of order k as ori-entation selective filter. In section 2.3 we will use the results of this section to derive a filter with a better orientation selectivity.
The Gaussian is only used for regularization of the derivative operators. Therefore the σ, which specifies the width of the Gaussian, is not of in-terest in what follows. The Gaussian derivative filter of order k in the φ direction is given by:
G(k)(x, y, φ) = (cos φ ∂ ∂x+ sin φ ∂ ∂y) kexp( −x 2+ y2 2σ2 ) (1)
Next consider an image I(x, y) that consists of a sum of strokes fi, with corresponding orientations θi. Thus I(x, y) can locally be described by:
I(x, y) =X i
fi(x cos θi+ y sin θi) (2)
The response of the kth order directional Gaus-sian derivative to the stroke fi when φ = θi, is denoted by fik,max. The response of the direc-tional derivative to the image I(x, y) can then be written as:
I(k)(x, y, φ) = G(k)(x, y, φ)∗(x,y)I(x, y)
=X
i
cosk(φ− θi)fik,max(x, y)(3)
Where ∗(x,y) denotes convolution in the (x, y) plane. This can be rewritten as:
I(k)(x, y, φ) =X i
coskφ∗(φ)fik,max(x, y)δ(φ− θi) (4)
Convolution in the φ-dimension is denoted by ∗(φ). The last equation has a simple interpreta-tion. In orientation space each stroke gives an impulse in the φ-dimension convolved with the following kernel:
a(k)(φ) = coskφ (5) A property of this kernel is that the orientation selectivity increases with increasing order k. This means that a specific selectivity can be chosen for a particular application.
A less desirable property is that not only the orientation selectivity changes with the or-der k, but the radial frequency sensitivity as well. Therefore we will follow Knutsson [5] and decom-pose the filter into an angular and a radial part. The Fourier transform of a directional Gaussian derivative using polar coordinates ω and θ is given by: e G(k)(ω, θ, φ) = (jω)kexp(−1 2σ 2ω2) cosk(φ − θ) (6) Although separable, both the angular response, cosk(φ− θ), and the radial response depend on k. Since we wish to vary only the orientation selec-tivity of the filter, we will keep the radial response fixed at k = 2 by multiplying eG(k)by ω−k+2. This is a rather arbitrary choice, but will suffice for the purpose of this paper. The actual choice of radial function depends on the application. The Fourier transform of the resulting filter A is as follows:
e A(k)(ω, θ, φ) = ω−k+2Ge(k)(ω, θ, φ) = jkω2exp(−1 2σ 2ω2) cosk(φ − θ) (7)
2.2
Sampling the orientation space
In practice only a finite number of filters can be used. Therefore the orientation space must be sampled. In this section we derive how many filters are needed for a given order k to allow re-construction.
Given a certain filter, Freeman & Adelson [2] and Perona [8] address the problem of creating a set of filters that allows interpolation and mini-mizes the error between this set and the filter re-sponse obtained by rotation of the original filter. They also present the filter constraints required to allow interpolation.
The following analysis is essentially identical to the derivations in the references given above and shows that our filter allows interpolation. The Fourier coeffients of cos(φ) are zero, except for ωφ =±1. Repeated convolution (k − 1 times) of
the Fourier series representation of cos φ by itself, yields the Fourier coefficientsea(k)(ω
φ) of a(k)(φ), which are then given by Pascal’s triangle:
e a(k)(ωφ) = 1 2 k k 1 2(k + ωφ) k odd, ωφ odd≤ k k even, ωφ even≤ k 0 elsewhere (8)
The coefficients are depicted in figure 2a for order 10. The filter is bandlimited and the co-efficients that correspond to frequencies higher than kωφ are zero. The Nyquist sampling the-orem states that we need 2k + 1 samples on the interval (0, 2π) to allow reconstruction, or equiva-lently 2k + 1 filters. The number of filters needed can be reduced to k + 1, by noting the following symmetry:
a(k)(φ + π) = (−1)ka(k)(φ) (9)
2.3
Improving the angular
resolu-tion
The kernel a(k)(φ) can be well approximated by a Gaussian shape with a variance that de-creases as 1k. Therefore the peak width decreases as √1
k. Since the kernels are periodic the notion of peak width only makes sense on the interval (−π 2, π 2). The 1/ √ k behaviour is disappointing, because one would wish that doubling the num-ber of filters would halve the peak width. We will now show how this can be achieved.
Equation 4 shows that each stroke gives an im-pulse convolved with the kernel a(k)(φ). To im-prove the orientation selectivity we can simply de-convolve the resulting signal, because the blurring kernel is exactly known. As equation 8 shows, some Fourier series coefficients are zero, mean-ing that only a partial deconvolution can be per-formed. For order k the amplitude of the Fourier coefficients can be flattened (deconvolved) to:
eb(k)(ωφ) = 1 k+1 k odd, ωφ odd≤ k k even, ωφ even≤ k 0 elsewhere (10)
In figure 2b the coefficients eb(k)(ω
φ) are depicted for k = 10. The resulting response in orientation space is given by:
b(k)(φ) = sin((k + 1)φ)
sin φ (11)
In figure 2c the kernels a(k)(φ) and b(k)(φ) are depicted for k = 10. It is clear that the peak width is indeed reduced, but at the cost of adding
some side lobes. These lobes will cause problems if the amplitude of the different strokes differs too much. The first zerocrossing of b(k)(φ) lies at π/(k + 1). This shows that the peak width indeed decreases approximately as 1/k.
The response along the φ axis in orientation space to a stroke is equivalent to the angular re-sponse of the filter that is used. So instead of performing the deconvolution, it is also possible to use a filter that uses equation 11 as its angular response. The Fourier transform of the filter cre-ated by replacing the angular part of equation 7 by b(k)(φ) is given by: e B(k)(ω, θ, φ) = jkω2exp(−1 2σ 2ω2)sin((k + 1)(φ− θ)) sin(φ− θ) (12) In figure 3 the A and B filters are shown in the frequency domain as well as in the spatial domain.
2.4
Measuring the rms orientation
intensity
The filters A and B do not measure the rms ori-entation intensity, but respond to either locally symmetric or antisymmetric structures. For tex-ture analysis we do not wish to discriminate be-tween the two. The rms intensity is measured by the following procedure; First the filter result is squared, yielding an estimate of the local orien-tation energy, but with an image containing high frequency signals superimposed on it. By adding the square of its quadrature counterpart [5] the high frequency signals will cancel. Creating a quadrature counterpart for filter B results in a filter that is discontinuous for even k and non-differentiable for odd k. Therefore it is not pos-sible to use a quadrature filter approach. Instead the image will be smoothed by a Gaussian to re-move the high frequency signals. The square root is taken of the smoothed image to yield the final estimate of the rms orientation intensity. Any in-terpolation must take place before the squaring operation, otherwise the necessary conditions to allow interpolation will no longer hold. All these operations take place in the (x, y) planes of ori-entation space.
3
Single Peak detection
In order to perform an analysis in orientation space, it is necessary to detect peaks in the φ di-mension at every (x, y) position. We assume that the peaks are non-overlapping and in the exper-iments we have ensured this by choosing a suffi-ciently high order k. The shape of the peaks is
-1.5 -1 -0.5 0.5 1 1.5 -0.2 0.2 0.4 0.6 0.8 1 φ ω ωφ −π− − 2 π2 0 0 0 0 252 210 120 45 10 1 (a) (b) (c)
Figure 2: a.) Fourier coefficients ea
(10)
(ωφ). The numbers show relative amplitude. b.) Fourier coefficients
eb(10)(ωφ). c.) The kernels a(10)(φ), dashed, and b(10)(φ), solid line
(a) A(9)(x, y) (b) eA(9)(ω
x, ωy) (c) B(9)(x, y) (d) eB(9)(ωx, ωy)
Figure 3: Spatial and frequency domain versions of the filters A and B
known and given by |bk(φ)|. Since we are not interested in the side lobes of this function, a simpler Gaussian shaped function will be used in the fit procedure. Because a fit has to be per-formed for every (x, y) position, traditional iter-ative methods such as Levenberg-Marquardt re-quire too much computation time. In this section we describe a fast alternative for the fit procedure. The σ of the Gaussian can be determined in ad-vance since the width of the peaks is known. The Gaussian was fitted to a peak for order 9 yield-ing σ = 0.160286 radians. The fitted Gaussian coincides almost perfectly with the main lobe. Since the peak width decreases approximately as 1/(k + 1), the following relation gives the σ for any k:
σ = 1.60286
k + 1 (radians) (13) To estimate the location and amplitude of the Gaussian the following error measure will be min-imized:
|f − Ag(ψ)|2 (14)
Where f is the orientation response (φ dimen-sion), g a Gaussian positioned at ψ with ampli-tude factor A and the σ as given by equation 13. The inner product between functions m and n is defined as:
(b) (a)
Figure 4: a.) Signal and inner product between shifted copies of the template and the signal (Gaus-sian filtering). b.) Signal and the detected peak.
m(φ)• n(φ) = Z
φ
m(φ)n(φ)dφ (15)
The integral is performed over one period of φ. Equation 14 can be written as:
|f|2+ |Ag|2
− 2Af • g (16)
Only the last term is a function of ψ. Minimiz-ing equation 14 with respect to ψ is equivalent to finding maxψ(f• g). Due to the symmetry of the Gaussian, convolving f with the Gaussian kernel yields f• g for every value of ψ. The position of the maximum of the resulting signal is the ψ that maximizes f•g. A typical f is shown in figure 4a. In the same figure f• g is shown for every ψ.
The next step is to fit the A parameter. Setting the deritivative of equation 16 with respect to A to zero and solving the resulting equation yields:
φ y φ x x y C A D B (b) (c) (d) (e) (f) (g) (h) (i) slice CD slice AB orientation (a) φ x y A D C B image original space k=19
Figure 5: a) Schematic representation of orientation space. An x, φ slice as indicated is used in (d-f). b) A noise free image containing two superimposed patterns. The intersection of the slice with the (x, y) plane is indicated by line AB. The intersection with a second slice is indicated by line CD. c) Noisy version of (b). d) The slice indicated in (a) for the rms orientation intensity measured by filter A in image (b). e) Same as (d), but using filter B. f) Same as (e), but applied to image (c). (g-i) Same as (d-f) but using the other slice as indicated by line CD. Note the structures next to (f) indicating where in the slice each pattern gives a response.
A = f• g
|g|2 (17)
In figure 4b the peak as detected by the algo-rithm, Ag, is shown plotted over the original sig-nal f . To improve the fit it is necessary to remove the background offset before using the algorithm as discussed above. Since each parameter is op-timized individually, there is no guarantee that a global optimum will be found. Maximizing f• g will first yield the peaks with the largest volume. Such a peak is either a true peak or may consist of several overlapping peaks (This does not hap-pen in the test images). Therefore a quality mea-sure is needed to assess whether our assumption of isolated peaks holds. We will use the following quality measure:
Q = f W • g
|fW| (18)
Where fW is a windowed version of f with the window’s center at the position of the detected peak. This is to prevent other peaks from influ-encing the quality measure. The window should be about as wide as the Gaussian, for example 4σ. Q lies in the interval [−1, 1], with 1 indicating a
high quality. By subtracting the peak that was found from f , a second peak can now be detected by repeating the procedure.
Since Gaussian filtering can be implemented very efficiently [13], the entire procedure is very fast.
4
Experiments
In figure 5 we show the results of applying the fil-ters proposed in this paper to an image consisting of two overlapping patterns. Each of the patterns has been created by calculating the Euclidean dis-tance to some point in the image and computing the sine of the result. The amplitude of each indi-vidual pattern is one. The second test image was created by adding Gaussian noise with a variance of one to the first test image (SNR=6dB). We have computed the orientation space for both im-ages using both filter A and filter B for k = 19. We show two slices of the orientation space for the following combinations: noise free image and filter A, noise free image and filter B, noisy image and filter B. The slices show the rms orientation intensity, using σ = 5 for the smoothing.
The resulting slices show that the orientation selectivity is indeed improved. It is also clear that the scheme is quite robust with respect to noise. Furthermore, close examination of the
(a) (b)
(c) (d)
Figure 6: a) Three touching patterns: left 51 degrees with respect to x axis, middle 36, right 21. Image size 224 by 224 pixels, diameter of circle 128 pixels. Gaussian noise with variance 0.2 was added to the image. b) x, φ slice of orientation space halfway along the y axis. Parameters used: k = 27, σfilter= 1.5, σsmoothing= 5. c)
Magnitude of second detected peak divided by magnitude of first detected peak for each x, y position. d) Detected boundaries superimposed on (a).
mation present along the orientation axis, and it is therefore possible to interpolate in (the origi-nal) orientation space to get improved accuracy.
The next experiment concerns touching orien-tation fields. A test image was created contain-ing a circular region in the middle and two regions around it separated by a second vertical boundary in the middle of the image. The circular region has an orientation of 36 degrees with respect to the x axis, the left region 51 degrees, and the right region 21 degrees. The amplitude of the patterns is again one, and Gaussian noise with a variance of 0.2 (SNR=13dB) was added.
Using the method described in section 3, the dominant peak in orientation space was found for each x, y position. In figure 6c the magnitude of the second peak divided by the magnitude of the first is shown. When two orientation fields are present, as at the boundaries of the regions, the two peaks found will have approximately the same magnitude and hence the quotient will ap-proach one. Indeed, the quotient is maximal at
the boundaries as the results in figure 6c show. Using a watershed [11] based technique, the ridges in figure 6c are detected and shown superimposed on the original image in figure 6d. They coincide almost perfectly with the actual boundaries.
Since the parameters were chosen to ensure that there are no overlapping peaks in orienta-tion space, it was not necessary to use the quality measure as defined in section 3. The location of the maxima is not determined with subpixel accu-racy in the current implementation. Subtracting a slightly misaligned peak results in larger resid-ual flanks. The orientation of the center region in figure 6a corresponds exactly with a pixel po-sition, whereas the other patterns do not. As can be seen in figure 6c, the residual peaks are indeed larger for the other two regions.
The location of the boundaries found in the second example is very accurate. However, as the difference in angle between the different regions decreases, a higher order filter will be required to distinguish between the regions. But, as the
wide narrow
Figure 7: Localization depends on orientation with respect to the boundary
filter becomes narrower in the frequency domain, the filter will become wider in the spatial domain. This will have a corresponding adverse effect on localization. Due to the shape of the filter as shown in figure 3c, localization depends on the orientation of the boundary with respect to the oriented pattern, see figure 7. Fortunately, the situation where the pattern(s) and the boundary have a similar orientation (allowing good localiza-tion) occurs most frequently in real world images.
5
Conclusions
We have shown that the peak width of the re-sponse in orientation space can be made to de-pend linearly on the number of filters k. This occurs at the cost of introducing some side lobes. Simple experiments with this scheme show that it works well for noisy images.
The method used to detect the boundaries is still rather crude and more advanced techniques to analyze orientation space must be developed. Real orientation fields may not only differ in ori-entation, but also in amplitude and frequency content. Furthermore the interplay between filter order, frequency content, attainable orientation selectivity and localization and its consequences must be investigated.
Acknowledgements
This work was partially supported by the Rolling Grants program 94RG12 of the Netherlands Or-ganization for Fundamental Research of Matter (FOM) and by the Royal Dutch Academy of Sci-ences (KNAW).
References
[1] M. Andersson, Controllable Multidimen-sional Filters and Models in Low Level Com-puter Vision, PhD thesis, Link¨oping Univer-sity, Sweden, 1992.
[2] W.T. Freeman and E.H. Adelson, The De-sign and Use of Steerable filters, IEEE trans-actions on Pattern Analysis and Machine
In-telligence, vol. 13, no. 9, September 1991, pp. 891-906.
[3] L. Haglund, Adaptive Multidimensional Fil-tering, PhD thesis, Link¨oping University, Sweden, 1992.
[4] M. Kass and A. Witkin, Analyzing Oriented Patterns, Computer Vision, Graphics and Image Processing, vol. 37, 1987, pp. 362-385 [5] H. Knutsson, Filtering and Reconstruction in Image Processing, PhD thesis, Link¨oping University, Sweden, 1982.
[6] J.J. Koenderink, The Structure of Images, Biological Cybernetics, vol. 50, 1984, pp. 363-370.
[7] M. Michaelis and G. Sommer, Junction Clas-sification by multiple orientation detection, in: Jan-Olof Eklundh (ed.) ECCV ’94, Third European Conference on Computer Vision, Stockholm, Sweden, May 1994, pp. 101-108. [8] P. Perona, Deformable Kernels for Early Vi-sion, IEEE transactions on Pattern Analysis and Machine Intelligence, vol. 17, no. 5, May 1995, pp. 488-499.
[9] L.J. van Vliet and P.W. Verbeek, Estima-tors for Orientation and Anisotropy in Dig-itized Images, in: J. van Katwijk, J.J. Gerbrands, M.R. van Steen, J.F.M. Tonino (eds.), ASCI’95, Proc. First Annual Confer-ence of the Advanced School for Computing and Imaging (Heijen, NL, May 16-18), ASCI, Delft, 1995, pp. 442-450.
[10] L.J. van Vliet and P.W. Verbeek, Segmen-tation of overlapping objects, in: L.J. van Vliet, I.T. Young (eds.), Abstracts of the ASCI Imaging Workshop 1995, Venray, The Netherlands, October 1995, pp. 5-6.
[11] B.J.H. Verwer, L.J. van Vliet and P.W. Ver-beek, Binary and Grey-value skeletons: Met-rics and Algorithms, International Journal of Pattern Recognition and Artificial Intelli-gence, vol. 7, no. 5, 1993, pp. 1287-1308. [12] A. Witkin, Scale space filtering, Proc. Int.
Joint Conf. on Artif. Intell., Karlsruhe, Ger-many, 1983, pp. 1019-1023.
[13] I.T. Young, L.J. van Vliet, Recursive imple-mentation of the Gaussian filter, Signal Pro-cessing, vol. 44, no. 2, 1995, 139-151.
