## Abstract

Two new two dimensional (2-D) complex operators for estimating the energy and orientation of 2-D oriented patterns are proposed. The starting point for our work is a new 2-D extension of the Teager-Kaiser energy operator incorporating orientation estimation. The first new energy operator is based on partial derivatives and can be considered a local (point-based) estimator. Using a nonlocal (pseudo-differential) operator we derive a second and more general energy operator. A scale invariant nonlocal operator is derived from the recently proposed spiral phase quadrature (or Riesz) transform. The Teager-Kaiser energy operator and the phase congruency local energy are unified in a single equation for both 1-D and 2-D. Robust orientation estimation, important for isotropic demodulation of fringe patterns is demonstrated. Theoretical error analysis of the local operator is greatly simplified by a logarithmic formulation. Experimental results using the operators on noisy images are shown. In the presence of Gaussian additive noise both the local and nonlocal operators give improved performance when compared with a simple gradient based estimator.

©2005 Optical Society of America

## 1. Introduction

The estimation of image feature orientation is important in many areas of image and pattern analysis, especially for orientation adaptive algorithms. Granlund and Knutsson [1] examined the idea of local orientation and its relation to the local simplicity hypothesis and local Fourier properties. Fringe patterns in interferometry are perhaps the archetypal locally simple images (also known as intrinsically one dimensional, i1D, by Krieger and Zetzche [2], as oriented patterns by Kass and Witkin [3] and as ridge systems by Penrose[4]). Although we restrict our coverage to fringe pattern analysis much of the work applies (with careful modification) to more general signals and images. In fringe pattern analysis the local fringe orientation can be utilized in directional filtering to improve fringe quality: spin filtering in optical interferometry [5–7] and fingerprint ridge enhancement by oriented filtering [1, 8, 9] being two pertinent examples. Recently the author proposed an isotropic 2-D demodulation algorithm [10] which relies upon good orientation estimates for its applicability to intricate fringe patterns. In subsequent publications [11, 12] an outline of the orientation estimation problem was given. The major flaw of gradient-based orientation estimators is their failure in the zero gradient regions (namely ridge tops and valley bottoms). The analogous problem is well known in 1-D demodulation, and has been elegantly resolved by the Teager-Kaiser energy operator [13, 14], designated TKEO from hereon. The energy operator has also been generalized to N-dimensions [15–18] as a magnitude without directivity.

Much of the preceding research on orientation has been concerned with either purely gradient based methods, or purely second derivative methods. We propose to show that by redefining 2-D orientation estimation as a quadratic (Volterra) mix of first and second partial derivatives we obtain an operator that gives energy and orientation estimates that are uniform and isotropic, yet requires just one complex filtering operation (in its extreme formulation). By defining the operator in Fourier space a new family of operators with isotropic performance is revealed. The new operators are related to the recently proposed spiral phase quadrature transform [10–12, 19].

The paper is structured as follows:

Section 2 | Reviews the previous approaches and problems of orientation estimation. |

Section 3 | Reassesses the energy operator literature and introduces a new complex 2-D energy operator. |

Section 4 | Establishes the replacement of local (derivative) operators by nonlocal (pseudo-differential) operators in the 2-D energy operator. |

Section 5 | Describes implementation details, both discrete spatial and spectral methods. |

Section 6 | Investigates the estimator performance in the presence of noise. |

Section 7 | Discusses possible applications (such as fingerprint, oriented patterns, interferometry, image analysis), and developments. |

Section 8 | Concludes. |

Section 9 | Appendix. Presents an alternative interpretation of the analytic image formalism presumed in section 2. |

#### Note

We have endeavoured to condense and fairly represent much of the previously unconnected, cross-disciplinary research which underpins our work on the 2-D energy operator. The task has been much harder than originally anticipated and only a fraction of the relevant research is actually referenced in the text. We apologise in advance for any significant papers that may have been overlooked or misinterpreted; regrettably “The truth is rarely pure and never simple” [20].

## 2. Orientation estimation

Estimation methods are categorized as follows in this section:

First order (linear), first degree derivative method (2.2)

Second order (quadratic), first degree derivative methods (2.3),

First order, second degree derivative methods (2.4),

Second order, first and second degree derivative methods (2.5).

#### 2.1 Direction and orientation

In all the following analysis we assume that the principal model for the image *f*(*x, y*) is a single component amplitude modulated – frequency modulated (AM-FM) function as defined by Havlicek [21], corresponding to “locally simple” as defined by Granlund and Knutsson [1], and equivalently “i1D” as defined by Krieger and Zetzche [2]:

Multi-component AM-FM functions are explicitly excluded from consideration. The envelope or amplitude is *b*(*x,y*), and the phase is *ψ*(*x,y*). There is an underlying assumption
here that any low frequency offset is absent, or has been removed by pre-processing. A similar assumption underlies much demodulation research [10], but it is worth noting that bandpass filtering, or scale space filtering is often the preferred preprocessor for images.

It is convenient to consider the function in Eq.(1) as the real part of a hypothetical complex function *f*_{A}
, analogous to the one dimensional analytic signal of Gabor [22]:

In one sense Eq. (2) can be considered as the *a priori* prototype function from which Eq. (1) is derived. Phase-shifting interferometry allows Eq. (2) to be deduced from a sequence of equally phase-shifted images corresponding to Eq. (1). However, we are primarily interested in the situation where we only have a single image corresponding to Eq. (1). We prefer not to digress into such general issues as how Eq. (2) might be derivable from Eq. (1), and defer further discussion to the appendix (Section 9), noting merely that there is still considerable disagreement about how the definition of analyticity should be extended for signals of more the one dimension.

The complex exponential signal of Eq. (2) allows the fringe (or local structure) direction β to be completely defined by, and calculated from, the complex gradient of the phase:

The direction is a cyclic quantity with principal value in the range -*π*
*<*
*β* ≥ *π*. In contrast the orientation is a quantity that is only well defined over half this range (see Jahne [23] for a detailed definition of direction and orientation). The orientation is applicable to patterns which are of the form shown in Eq. (1), wherein the local reversal of the phase leaves the pattern unchanged:

It is well known that the 180° periodicity of orientation can be gracefully incorporated into complex algebra by the consistent use of the double angle formalism [1]. To be more precise all explicit estimators must contain phase functions of (2*β*) and not (*β*) alone. Our analysis follows this formalism.

#### 2.2 Simple Gradient based methods (first order, first degree)

Kass and Witkin [3] proposed a gradient based technique for analysing oriented patterns, but observed difficulties with aggregating orientation information. Shu and Jain [24] continued the gradient approach. Andresen [5] used the gradient approach to design an orientation adaptive filter. Estimation based on the simple gradient is the classic approach proposed by numerous other researchers, and its problems are discussed in appendix A of Larkin [11]. The main problem is that both components of the gradient of *f*(*x,y*) are near zero at the ridge tops and valley bottoms of the pattern. The estimator is

where the subscripts denote partial derivatives with respect to x and y respectively. On the ridges and in the valleys cos(*ψ*) = ±1, sin(*ψ*) = 0, hence the angle estimate loses its dependence on the phase structure and instead follows the envelope gradient :

Consequently the estimator defined by Eq. (5) is not uniform across features as the phase *ψ* varies, and for this reason uniformity is also known as “phase invariance” [1, 25, 26]. In practice Eq.(6) implies high noise sensitivity at extrema because both numerator and denominator are very small. Several authors [27, 28] have proposed noise reduction in gradient maps by smoothing or averaging the intrinsic errors over a small neighborhood. Recently Quiroga et al [29] have proposed a method of regularizing the orientation estimate as part of a direction unwrapping scheme, but the initial estimates still contain high noise values. Ultimately methods based on the gradient alone are destined to fail at extrema; consequently we are encouraged to look beyond first order methods for improvements.

#### 2.3 Quadratic and tensor based methods (second order, first degree)

Much has been published on second (quadratic) order methods for orientation estimation. Intuitively attractive because of the link between the double-angle formalism and second order forms, the basic idea is presented in the classic work of Granlund and Knutsson [1], and Jahne [23]. A (structure) tensor is constructed from quadratic products of gradient components. The products are actually smoothed (over a local region) before the tensor is constructed. Orientation estimation is accomplished by the direct or iterative solution of the eigenvalue problem for the tensor. We have omitted details of the technique because it is not designed for patterns of the form given in Eq. (1), and is therefore prone to the problem exemplified in Eq. (6), although details are obscured by the smoothing operation. A similar approach, termed “tensor whitening” [30] has been proposed to ameliorate low signal sensitivity of orientation filters.

#### 2.4 First order (linear), second degree derivative methods

A comprehensive review of the application of (purely) second partial derivatives to structure determination in images was given by Danielson *et al*. [31]. It transpires that in 2-D three features have to be disentangled from the second derivatives: energy, orientation, and shape. Although an ingenious method is presented for extricating the orientation the method fails near the zero crossings of Eq.(1) where *all* second derivatives *and* the Hessian are exactly zero (and as a result noise can dominate any estimate). More recently Da Costa *et al*. [32] recognized that gradient-based orientation estimators alone do not operate uniformly over oriented fringe patterns, and that second derivatives are most effective in the regions where first derivatives are least effective. Their proposed solution estimated the likelihood or “coherence” of both estimates and makes a binary decision on which to use. The connection with our work is considered in more detail section 3.

#### 2.5 The “Local Energy”, the 2-D Energy Operator and the Energy Tensor (quadratics of first and second degree derivatives)

In the area of human vision research a number of attempts have been made to develop an energy model related to the local phase coherency of spatial frequency components. In the influential work of Morrone and Burr [33] they note that idea of a “local energy” function was first proposed by Adelson and Bergen [34] in 1985. The Adelson-Bergen local energy function essentially combines directionally filtered quadrature functions in a fashion analogous to that of the human visual system. The connection with our nonlocal energy operator is explored further in section 4. The connection of the “local energy” with the 1-D energy operator of Teager-Kaiser-Maragos-Quatieri-Bovik is reviewed in section 3. There is an even closer parallel with the analytic signal first proposed by Denis Gabor [22]. The related demodulation approach to human vision is well summarized by Daugman [35].

The 2-D energy operator proposed by Yu and Mitra [15] and later by Maragos *et al*. [17, 18] is simply the sum of 2 orthogonal energy operators. The result is an isotropic operator. Unfortunately it is not possible to estimate the orientation of the 2-D frequency because the two (normally positive) orthogonal terms only define orientation within a quadrant.

The idea of combining first and second degree derivatives in a 2-D energy tensor was recently proposed by Felsberg and Granlund [25] and is discussed further in the next section.

## 3. A New Definition of the Two Dimensional Complex Energy Operator

Our initial work has been motivated by a search for a local (or point-wise) orientation estimator that does not require *a priori* knowledge of the local fringe spacing (i.e., it is relatively scale invariant). Clearly any scheme that requires local smoothing is not scale invariant. Another important criterion for an estimator is uniformity: we would like an estimate that is equally reliable at any point in our model, Eq. (1), and not dependent on the sinusoid phase. Reiterating, the term “uniformity” is equivalent to the term “phase invariance” used by other researchers [25].

#### 3.1 Essential Properties of the Teager-Kaiser Energy Operator (TKEO) in One Dimension

It is instructive to briefly review the 1-D TKEO defined by

The main property of the TKEO is that it gives an output proportional to the square of the amplitude and the square of the frequency of an oscillating input signal. If the input signal is compared to a simple harmonic oscillator, the TKEO outputs a measure of the total energy of the oscillator. The first term on the RHS of Eq. (7) is analogous to the kinetic energy whereas the second term is analogous to the potential energy of the oscillator. More details on the properties of the TKEO can be found in a very large literature base (see [13, 14, 17, 18, 36–43] for a small selection).

The main property is best illustrated by the example of the AM-FM signal *g*(*x*) = *b* cos(*ωx*) having modulations varying slowly with respect to the carrier period:

Essentially the TKEO generates two oscillating terms in phase quadrature: the first - related to the gradient - is a sine, the second is a cosine. The operation is equivalent to squaring and adding two quadrature terms, and the oscillation exactly cancels out resulting in a total energy measure *ω*
^{2}
*b*
^{2}. We refer to this property as “uniform estimation” or uniform demodulation. Significantly, operators like the gradient are not uniform demodulators on a sinusoidal signal, but produce a rapidly varying output:

#### 3.2 Extending the Teager-Kaiser Energy Operator to Two Dimensions

We would like to perform the Teager-Kaiser magic in two dimensions; thereby avoiding the non-uniformity problems of both the gradient based methods and the purely second differential methods. The method proposed by Yu and Mitra and, later, Maragos *et al*. [15, 17, 18] simply adds TKEO outputs from orthogonal 1-D operators. From our perspective the problem with this formalism is that it is a positive scalar with no directional information, even though the input is nominally constrained to be highly directional. A first guess for a 2-D TKEO would be to form a vector from the individual x and y TKEO outputs. The main problem is that such a vector can only define orientations within a quadrant and is quite powerless to distinguish between the distinct orientations +45°and -45°. Using the individual x and y TKEO components alone it is not possible to go beyond the quadrant limit. The missing directional information is actually located in the cross derivatives, but the inclusion of cross derivatives would seem to undermine the simplicity and elegance of the TKEO.

Observe that the TKEO in Eq. (7) can be written in simple operator form, if the derivatives are replaced by a general differential operator *D*.

Then the operator form is suggestive of a possible 2-D extension, namely replacing the operator *D* by a two dimensional gradient operator. Our previous work [10] benefiting from the efficiencies inherent in complex image processing has predisposed us toward a complex (rather than the traditional vector) interpretation of the 2-D gradient operator. The complex approach is gaining some popularity in diverse areas of image analysis [44–46] and restoration [47]. We define a complex gradient operator analogous to the complex differential or Cauchy operator of complex analysis (except for a factor of 2 and a sign change):

and apply it, in the form of Eq. (10), to the fringe pattern defined in Eq.(1)

Splitting into real and imaginary parts gives

The combination is clearly a mixture of quadratic forms of first partial derivatives and mixed forms of zeroth and second partial derivatives (including crucial cross terms). But what does it do to a simple linear fringe pattern? In this case a fringe pattern with spatial frequency *q*
_{0}, fringe normal direction *β*, and phase offset *ψ*
_{00}:

$$\mathrm{where}\phantom{\rule{.2em}{0ex}}\left\{{u}_{0}={q}_{0}\phantom{\rule{.2em}{0ex}}\mathrm{cos}\beta ,\phantom{\rule{.2em}{0ex}}{v}_{0}={q}_{0}\phantom{\rule{.2em}{0ex}}\mathrm{sin}\phantom{\rule{.2em}{0ex}}\beta \right\}.$$

The complex differential operations simplify:

In effect the cross terms of the zeroth, first and second derivatives coalesce and eliminate the sinusoidal modulations:

Finally, introducing the oriented fringe parameters of Eq.(13), we obtain a complex estimator that has uniform magnitude across a sinusoidal fringe:

The magnitude corresponds exactly with previous incarnations of the TKEO, and is proportional to both the amplitude-squared and the frequency-squared of the original cosine. The novel feature is the phase term related to twice the direction angle of the fringe. As many researchers have noted before (in particular Granlund and Knutsson [1] and Jahne [23]) the angle doubling corresponds to the angle *β* being an orientation defined modulo *π*, whilst the complex exponentiation automatically accounts for the phase-wrapping. Equation (16) is really quite remarkable; it represents an isotropic orientation estimator that does not fail at points of zero gradient. Recently Felsberg and Granlund [25] have proposed a 2-D energy tensor that closely resembles our 2-D energy operator, but where the grad operator is replaced by the grad tensor and the second differential operator becomes the Hessian.

It can be shown that our definition of 2-D energy operator also works for complex images. However, the operator must be applied to the real and imaginary parts separately and the results added as follows:

The output is complex. This is structurally consistent with the complex energy operator defined by Hamila [48] for one dimensional signals, which uses symmetric sums of conjugate differential operators to obtain a real output:

$$\mathrm{RHS}\equiv {\left(\frac{d{f}_{R}}{\mathit{dx}}\right)}^{2}-{f}_{R}\frac{{d}^{2}{f}_{R}}{d{x}^{2}}+{\left(\frac{d{f}_{I}}{\mathit{dx}}\right)}^{2}-{f}_{I}\frac{{d}^{2}{f}_{I}}{d{x}^{2}}.$$

An equivalent definition was given by Maragos and Bovik [18] to ensure that their energy operator gives a real output from a complex input. They considered signals of two or more dimensions and obtained a sum of individual components with the result again real. The Maragos approach is suitable whenever the total energy in a multidimensional signal is the quantity of interest. However, in applications where the multi-dimensional structural components are of interest it is desirable to have both directional and magnitude quantities. In the simplest implementation in N dimensions, each of the components can be associated with the components of an N-vector. Such an extension is not proposed by Maragos et al, although it is a possible starting point for a directed energy operator.

An obvious question about our development above is: why not use a vector formalism? The short answer is that the correct combination of partial derivatives cannot be obtained from a simple vector operator, and the elegance of the TKEO formalism is therefore diminished for a vector formalism. It is true that care must be taken to interpret the real and imaginary parts of the complex 2-D energy operator, however there is much to be gained from an operator that fits into current software implementations of complex image processing. Furthermore, it can be shown that the apparent problems with (complex) channel crosstalk [49] do not necessarily require the n-vector formalism for solution.

#### 3.3 Some observation regarding Subharmonicity and the proposed 2-D Energy Operator

In 1-D the TKEO can also be expressed as the Laplacian of a logarithm of a signal, thus going back to Eq. (7) we obtain:

Although insightful, the LHS form of the equation is not particularly useful numerically for values of *g* near zero, owing to the divergence of the logarithm. The above equation was used by Bovik to show that the positivity restriction of the TKEO [40] is equivalent to concavity of the logarithm. More generally the condition defines subharmonic functions in complex variable theory (see Hormander [50] for example):

We can rewrite Eq. (19A) for the 2-D energy operator; again the LHS form is numerically problematic although mathematically concise:

If *g* is real then log *g*
^{2} is also real, giving rise to a well defined operator except at the zero crossings of *g*. The above equation shows that the 2-D energy operator can, in principle, be implemented with just one complex filter; the Hessian-like *D*
^{2} operator, and three other algebraic operations (squaring, multiplication and taking the logarithm). This result is rather surprising when we realize that the quadrature filter approach of Knutsson [1, 51] requires at least three complex linear filters (preferably spaced at 60°intervals). The efficiency is due, in part at least, to the energy operator’s nonlinear structure [42, 43] as a homogeneous quadratic (second order) Volterra filter [52]. Equation (19B) can also be used to represent Felsberg’s energy tensor in terms of a Hessian tensor operating on the logarithm of the signal:

Equation (20) indicates that the positivity of the energy operator is equivalent to the subharmonicity of the logarithm of a signal. In his 1970 book (page 8) on Hardy spaces, Duren [53] observed that if a real-valued function *u*(*z*) of a complex variable *z* is harmonic
in a domain, then |*u*(*z*)|^{p} is subharmonic in the domain if *p* > 0 . Furthermore log^{+}|*f*(*z*)|, the logarithm of the modulus of an analytic function *f*(*z*), is automatically subharmonic. Clarification of these tantalising connections is currently under investigation.

#### 3.4 Almost uniformity of Da Costa scheme

A previous attempt by Da Costa et al [32] to design a uniform orientation estimator that is invariant to phase came close to our energy operator. They recognised that the second derivatives can be nonlinearly combined with pure gradient methods to eliminate null regions.

However instead of utilising the perfect uniformity of the cos^{2} + sin^{2} = 1 squared combination they chose the rather less optimal |cos| + |sin| modulus combination rule which gives non-uniformity between 1 and √2.

## 4. A Nonlocal Complex Energy Operator

The previous section introduced a 2-D energy operator that directly parallels the classical 1-D TKEO. In this section we propose an absolutely scale invariant 2-D energy operator based on pseudo-differential operators. The 1-D parallel will eventually become clear.

#### 4.1 Spiral-Phase-Riesz Transform Formulation

The 2-D operator defined in section 3, (Eq. (10)), neatly combines demodulation and orientation estimation. It can be classified as a nonlinear local (point) operator since it combines the local (or point) derivatives in quadratic products. It is well known that derivative operations tend to enhance noise. Even though the 2-D energy operator works well with perfect fringe patterns, the performance degrades significantly for noisy input.

The complex gradient operator used in Eq. (10) is isotropic and has a Fourier transform multiplier defined as follows, where F{} is the Fourier transform operator:

One interpretation of Eq. (22) is that the complex derivative operator is a spiral phase multiplier in the Fourier domain, with a linearly increasing (conical) radial frequency factor *q*. The increasing gain with frequency *q* is the main cause of the noise enhancement associated with the gradient operator. If we now reconsider Eq. (16) we see that the 2-D energy operator quadratically (*q*
^{2}) enhances the magnitude. The question inevitably arises: if the gradient operator were replaced by another spiral phase operator, but without the radial frequency factor, would it still work? Previously we defined a spiral phase quadrature operator for fringe pattern demodulation [10]. We define a new modified operator but call it *D*_{m}
to show that it is a modified differential operator. In reality it is a type of pseudo-differential operator, better known in the pure-mathematical realm of harmonic analysis as the Riesz transform. The essential Fourier multiplier property is that it has the same Fourier phase as the differential operator, but unit Fourier magnitude as follows:

Note that the operator defined above is slightly different to our “spiral phase quadrature transform” in that it contains an extra imaginary factor. This is purely for the convenience of paralleling the complex differential operator in section 3, and avoiding a superfluous negative sign in the final energy. As we know, this multiplier has an associated spatial domain convolution kernel that is itself a spiral phase, but with a singular inverse square drop-off:

Where the spatial polar coordinates are(*r,θ*). For locally simple fringe patterns (which satisfy quite loose restrictions on the rate of change of frequency and curvature) it can be shown using asymptotic methods [11] that the following approximations are applicable:

Consequently a modified energy operator based on these circular harmonic operators gives the following intriguing result:

The use of higher order spiral phase (Riesz) operators was foreshadowed in the author’s thesis [19]. One interpretation of the above result is as a 2-D dirigible (steerable) analytic signal demodulator, with orientation free of charge. In one dimension the TKEO is considered to parallel the magnitude of the analytic signal, [41] or “local energy” (from the phase congruency perspective [54, 55]). The modified 2-D energy operator in Eq. (26) parallels the differential 2-D energy operator in the same way. The derivative is replaced by the Hilbert transform in 1-D and the complex gradient is replaced by the spiral phase operator (or Riesz transform) in 2-D. It seems rather paradoxical that the phase congruency “*local* energy” is derived from a signal using the Hilbert (1-D) or Riesz (2-D) transform, both of which are archetypal *nonlocal* operators. We prefer to use *operator localization* (rather than *analytic signal localisation*) to characterize the two new energy operators we have presented. Accordingly Eq. (10) defines a local 2-D energy operator, whereas Eq. (26) defines a nonlocal 2-D energy operator.

#### 4.2 Unification of Kaiser-Teager energy operator and the “Local Energy”

The Kaiser-Teager energy operator and the phase congruency “local energy” in one dimension can be unified under a generic operator that also unifies the family of two dimensional operators proposed in this paper. The trick is to utilize the fact that applying a 1-D Hilbert transform twice is just the negation operator. This means that the usual equation. for “local energy” [55] which sums the squares of the in-phase signal *g* and quadrature signal *ĝ* can be written in the canonical operator form of Eq. (10) using **S** to represent the Hilbert transform operator:

It is perhaps worth noting that the above LHS form of energy operator, being the sum of squares, is never negative (in contrast to the differential energy operator which is only rarely negative).

#### 4.3 Other Scale Invariant Formulations

The 2-D energy operator of Eq. (26) has flat spectral response, giving a result independent of local frequency. Application of Plancherel/Rayleigh theorem [56] shows immediately that the total energy measure is unaltered by operators with unit-magnitude Fourier multipliers. The operators *D* and *D*
^{2} can be further generalised to homogeneous functions with index *α* :

Accordingly that the gradient based operator can be characterised by *α* = 1, *D*_{gradient}
= *D*
_{1} and the Riesz operator characterised by *α* = 0, *D*_{Riesz}
= *D*
_{0} . Local spectral analysis with the generalized energy operator the gives the following result:

The localisation of the Riesz operator and its alpha generalisation is important in the study of the uncertainty principle. In a fascinating twist it can be shown that in 1-D the Hilbert transform is completely anti-local (an extreme form of nonlocality, page 485 of Havin [57]).

#### 4.4 Radial Spectral Filter Formulations

We may wish to apply our energy operators to images, but the assumption of local simplicity, valid for fringe patterns, is no longer applicable. Felsberg [25] uses a “difference of Gaussians” filter to invoke the Gaussian scale-space guarantee of positive eigenvalues for their energy tensor. Granlund and Knutsson [1] prefer radial lognormal filters to limit the signal spectrum. Havlicek suggests several forms of multi-band filtering [58]. Almost any desired filter can be integrated into a 2-D energy operator by defining and applying the differential operator in the Fourier domain. A radial and polar separable Fourier multiplier could be of the form

where the absolute value |*R*(*q*)| ensures that no sudden sign flipping occurs in the orientation estimate.

## 5. Implementation

#### 5.1 Discretization: Spatial Domain

It is well known that the 1-D energy operator has a very compact and neat discrete formulation:

Only three adjacent samples are needed to implement the operator. In practice it is necessary to remove the DC or background signal, so the above triplet is replaced by a zero mean triplet [59]. An alternative is to replace each sample by a symmetric finite difference [42, 43], necessitating 5 samples in total:

Corresponding attempts to discretize the gradient form of the 2-D energy operator (defined by Eq. (9) and Eq. (11)) cannot achieve isotropic performance with small kernels. The main difficulty is avoiding spectral phase reversals (and spectral magnitude variation) related to the finite difference approximation of differential operators. Although it is possible in principal to implement an effective local estimate of the gradient using just a 3×3 region (and its auto-convolution the 5×5 kernel Hessian operator), we recommend that larger kernels be used to improve the spectral accuracy to the desired response.

Assuming a single frequency input and D operators with the following spectral magnitude (M) and phase (P) response, Eq. (14) and Eq. (15) are modified as follows by the discrete operator response:

The single frequency input results in the following energy operator output

Clearly, if the D operators do not originate from the same source, then magnitude ${\mathbf{M}}_{1}^{2}$ ≠ **M**
_{2}, and phase 2**P**
_{1} ≠ **P**
_{2} mismatch can occur in the RHS of Eq (33), and the uniformity condition will no longer apply. Initial analysis indicates that the simplest 3×3 and 5×5 kernels can produce systematic orientation errors of up to 8°. There is scope for discrete optimisation of the 2-D energy operator but it will not be considered further in this paper. Felsberg et al [26, 60] have recently considered spatial discretization in more detail.

#### 5.2 Discretization: Fourier Domain

Typical computer speeds are now such that a 1024 by 1024 pixel fast Fourier transform (single precision, real data, 2.8Ghz P4 Xeon) takes about 45 milliseconds [61]. Consequently for many applications it is realistic to perform filtering as complex multiplication in the Fourier domain. This allows full control of magnitude and phase over all frequencies. We have found in practice that the following Fourier multipliers give good results:

For much image analysis the image will be split into subbands or scale spaces before application of energy operator as a feature detector [25]. Note that the sub-band filtering and the differential operator filtering can be efficiently combined into one Fourier operation. In the case of the spiral-phase-Riesz implementation of the nonlocal energy operator (where the nonlocal filter kernel is delocalized over the full image area) it is more efficient to implement filtering in the Fourier domain.

## 6. Performance

#### 6.1 Test Pattern Design

NOTE:The following tests are applied to images that can be classified as locally simple, or equivalently i1D (as described in section 2.1). This means that multi-band separation is not required as a pre-processing step for the local and nonlocal 2-D energy operators. Felsberg et al [25, 26, 60, 62] have considered the pre-filtering issues necessary for the application of their closely related energy tensor to more general (i2D) images.

Our objective here is to develop and use a fringe pattern that tests representative properties of orientation estimators. One of the problems with the evaluation of orientation algorithms is how to define a good test image or image sequence to test the properties over a wide range of conditions. As our basis pattern we have chosen a periodic pattern (not unlike that of Quiroga [29]) which evenly fills half the Nyquist frequency range in both spatial directions. We circumvent edge effects by enforcing periodicity in this initial analysis. The pattern contains almost equal components at all frequencies and orientations within the half Nyquist zone, the pattern has circular and hyperbolic regions with undefined orientation (zero phase gradient), where we might expect orientation estimators to perform poorly.

The test pattern is the real part of the following pure phase function defined over a 256×256 pixel image:

In actual tests we quantize this to 8 bit representation typical of many imaging devices

In orientation tests we can compare the estimator with the exact theoretical orientation *β*_{test}
give by

Figure 1 shows the test image and the Fourier transform magnitude of the test image. The maximum horizontal or vertical frequency present in the test image is *π*/2 radians per pixel - precisely half the Nyquist rate.

Felsberg [25] mentions that the energy operator has the elegant property of exactly compensating the aliasing inherent in the gradient squared term alone, and the function times Hessian term alone, thus oversampling is not required in practice. In one dimension it is known that the energy operator can perform remarkably well with significant (3x) undersampling [19, 42, 43], and there is a precise interpretation is in terms of bandpass sampling theory. Nevertheless our test pattern here is 2x oversampled to reduce any printing artifacts that might otherwise appear.

#### 6.2 Systematic Error Theory

In 1-D the error propagation for the TKEO has been analysed in great detail [38, 40]. The estimation error is proportional to the second derivative of the phase (namely the rate of change of the frequency).

Error analysis in 2-D is greatly simplified by the use of Eq. (19B) which uses the Hessian of the logarithm. Amplitude modulation affects are further simplified by the use of the exponential (or attenuation *ρ*) representation, *b* = *e*^{ρ}
, in Eq.(1) so that :

Equation (16) then becomes:

The RHS consists of three terms; the first is the desired 2-D energy operator, the second term is an error proportional to the Hessian of the attenuation, and the third term is an error proportional to the Hessian of the phase. The logarithmic formulation effortlessly removes the complication due to negative amplitude modulation (since the squaring operation eliminates awkward *π* phase flips). The relative errors oscillate at twice the fringe frequency and are inversely proportional to the fringe frequency squared. Significant errors (>10% for example) only occur when the amplitude or frequency change substantially over the period of a single fringe, or when the fringe radius of curvature is less than a few fringe periods.

Systematic error analysis for the nonlocal 2-D energy operators is more involved because there is no simple logarithmic formulation. Nevertheless the asymptotic analysis of the (single) spiral phase quadrature transform [11] indicates that the quadrature phase has an error directly proportional to the curvature of the underlying fringe phase *ψ*. The analysis can be extended to the double spiral phase operator (second order Riesz transform), but this would constitute a research paper on its own, and is not considered further here.

#### 6.3 Energy Operator Performance: High Signal to Noise Conditions

For relatively slowly changing fringe frequencies the 2-D energy operators exhibit very small systematic errors, and their robustness in the presence of noise is of more practical concern. Figure 2 shows the performance of an idealized energy operator as a magnitude (energy) and phase (orientation angle) map. The magnitude is presented with a linear scale (black=zero, white =maximum value). The phase is presented with a linear scale between -*π* (black) and +*π* (white).

By way of comparison we have computed a simple gradient based orientation estimator (*D*_{g}*f*)^{2}, using simple finite differences, and applied it to the 8 bit quantized test pattern of Eq.
(35). The result is shown in Fig. 3. The intrinsic noisiness of the estimate is clear, even when the input is an 8 bit image with SNR of 53dB (as defined by Granlund [1] on page 257):

The characteristic of all gradient based methods is the failure at fringe extrema. In this particular example the increased noise sensitivity is clear around the fringe maxima. The gradient orientation estimator also fails at half the Nyquist frequency (the diamond shaped line in Fig. 3) due to aliasing of the (*D*_{g}*f*_{s}
)^{2} operator.

Figure 4 shows the performance of the local (differential) 2-D energy operator as a magnitude and phase map. Notice how closely it resembles the ideal shown in Fig. 2, the only visible difference being small variations of the zero to 2*π* phase transitions.

Figure 5 shows the performance of the nonlocal (spiral-phase) 2-D energy operator as a magnitude and phase map. Again notice how closely the phase resembles the ideal shown in Fig. 2. The magnitude however is quite different, primarily because it now represents a 2-D signal magnitude and is largely independent of fringe frequency. Observe that the magnitude is quite uniform, except where the original fringe patterns have singular structure. Crucially the orientation phase is close to ideal, with some dimly visible artifacts near the singular points of zero fringe frequency.

#### 6.4 Energy Operator Performance: Low Signal to Noise Conditions

Typical oriented patterns (such as optical interferograms and fingerprints) have relatively low signal to noise ratios. A useful measure of estimator performance is the RMS phase error (in degrees) versus SNR. Such a performance measure was calculated by Granlund ([1], page 258) for the orientation tensor, although their test pattern avoids edge effects by windowing fringe modulation, not by periodic tiling. Their pattern also carefully avoids very low frequency and high curvature effects that occur at the centre and edges of our pattern. This means that our results are not directly comparable to theirs. For large phase errors it is preferable to use period statistical measures which avoid phase wrapping artifacts [63, 64], in the same way that Granlund’s arcsine of the tensor estimate gives properly wrapped phase.

The most visually representative case is for the 10 dB SNR, as defined by Eq. (38). The noise is zero mean Gaussian. We have applied the three estimators to the pattern in Fig. 6 with the following results:

The sequence of results shown in Fig 7 is revealing. The simple gradient phase shown in Fig. 7(b) is still working at 10dB, but is beginning to become overwhelmed by the intrinsic high frequency noise from the spatial derivative. The differential energy operator phase shown in Fig. 7(d) is clearer, but the RMS orientation error is actually only marginally better. The spiral phase (Riesz) energy operator phase shown in Fig. 7(f) is much clearer again as the noise has not been amplified by any spatial derivatives. The results for other SNR ratios are summarized in the following table (Table 1)

The results for the spiral phase operator are directly comparable to those of Granlund, but it should be noted that our test pattern contains features with much higher curvature, which we expect to undermine the initial assumptions of i1D.

#### 6.5 Noise reduction by Low Pass Filtering

The preceding example test patterns have rather slowly varying orientation phase patterns. In 1996 Strobel [65] proposed a method of filtering noisy fringe patterns (especially electronic speckle interferograms) by forming a complex (phasor) image that can be simply low pass filtered. The method elegantly avoids all the phase wrapping problems normally associated with processing of phase images, and has been an inspiration for a systematically complex approach to fringe analysis [10, 66]. Although it is true that many fringe patterns have rather slow variation of orientation over the image, some others (like fingerprints) have rapid variations of orientation (discontinuities, minutiae, folds etc) that would be simply lost if low pass filtered. Our proposed local and nonlocal operators have the potential to track very fast changes in orientation. For this reason phasor low pass filtering is not considered further.

#### 6.6 Energy Operator Performance: Fingerprint Image

As an example of how the new energy operators perform on real images a publicly available digitized fingerprint image from the NIST test data CD-ROM [67] is used, as shown in Fig.8.

The original fingerprint image was first cropped to 768×768 then Fourier down-sampled to 384×384. The only other pre-processing is the removal of spectral components near DC. Figure 9 shows the results of applying the three basic orientation estimators of sections 6.3, and 6.4 The differential energy operator gives considerably better orientation estimation than the simple gradient-squared estimator, which shows substantial high frequency noise. The spiral phase based estimator gives an even more recognizable orientation phase in Fig. 9(f) as one might expect from an operator that does not enhance high frequencies. The magnitude images shown in Figs. 9(a), (c), and (e) all show good discrimination between good signal modulation regions (white) and poor signal modulation (black). The simple gradient-squared magnitude picks out the edges of the saturated ridge structure very clearly especially when compared to the differential operator in Fig. 9(c). The fingerprint orientation estimation using energy operators is surprisingly good considering that the image is almost binary and scarcely resembles the smooth grayscale sinusoid underlying Eq. (1).

## 7. Discussion

#### 7.1 Filtering performance of local and nonlocal operators

Care must be taken interpreting the results of the preceding section. Both the squared-gradient operator and the local 2-D energy operator have differential structures that give them a high pass filtering character. Both are local (point) operators, hence it is not surprising that they appear to enhance noise in noisy fringe maps. Nevertheless it is clear from the tabulated noise data that the local 2-D energy operator has about half the orientation noise of the squared-gradient operator. The noise improvement has been gained by careful use of the second partial derivatives, without loss of signal resolution. In contrast low pass filtering of the squared-gradient operator would reduce noise, but it would also reduce resolution.

Introducing the nonlocal (spiral phase) 2-D energy operator we achieve further diminution of orientation noise at low SNR. The spiral phase (Riesz) structure of the operator has a flat frequency response, neither enhancing nor attenuating the noise spectrum. One could argue that the nonlocal spiral phase kernel in Eq. (24) has a large domain and then it is only to be expected that it tends to average out the effects of noise on the orientation estimate. In fact there is a view that our comparison of the two local estimators (having small spatial kernels) with the nonlocal estimator (having a kernel as big as the image itself) does not make sense because the operator windows are so disparate. The important point, however, is that the nonlocal energy operator has the extraordinary ability to demodulate without explicit low pass filtering and as a consequence the resolution loss associated with low pass filtering is completely avoided simply by using Eq. (26 ). A similar process is well known in one dimension with the nonlocal character of the Hilbert transform.

#### 7.2 Signal Bandwidths

A common misapprehension with both the spiral phase quadrature transform and the 2-D energy operator is that the method is limited to narrow bandwidth signals (see for example Felsberg [25] section 2). As can be seen from the (i1D) examples in section 6, the signal covers frequencies from DC right up to half Nyquist. It is more correct to say that the input signal must be *locally*, not globally, narrowband. Nevertheless it is true that in the application to more general images (i2D), the operator must be limited to small bandwidths by pre-filtering.

#### 7.3 Applications: Fingerprint Analysis

An obvious application is in the demodulation of fingerprint patterns. The preliminary results shown in section 6.6 are very encouraging for automation of fingerprint analysis because they are totally scale invariant and do not require the selection of band pass filters. This work represents the second stage in our search for an efficient method to analyse, compress, and accurately resynthesize human fingerprints [68]. To do this we must reliably demodulate a 2-D signal into its amplitude and phase. The first stage is based on the isotropic quadrature transform [10, 11]. The third and final stage will be presented in a forthcoming paper by the author. Daugman originally suggested in 1987 that image demodulation could be used for compact image coding [69]. In his particular case the signal of interest was another biometric: the human iris pattern.

#### 7.4 Applications: Curvature Estimation

Once the orientation *β* has been robustly estimated, the curvature is easily derived by differentiation – the rate of change of orientation (see [1] page 361). Continuing with the complex *D* operator approach, and utilising the compactness of complex algebra [65, 66] we find

The LHS of the above equation defines an effective algorithm for computing the curvature without explicitly having to compensate for wrapping of the orientation angle. The usual approach to curvature estimation typically fails in zero gradient regions (see Koenderink et al [70] or van de Weijer [71] for more detailed discussion), whereas Eq.(39) has no such difficulty as it operates on a uniform orientation estimate.

#### 7.5 Applications: Image Analysis

Felsberg and Granlund [25] have already proposed the use of the 2-D energy tensor for point of interest detection in images and subsequently a new variant they call the Gradient Energy Tensor (GET). Countless image processing algorithms depend on structural estimators; hence whenever robust estimation of orientation or curvature is required the 2-D energy operator may be the most effective method available.

#### 7.6 Alternative Implementations

In the world of fringe analysis the 2-D energy operator can replace a multitude of ad-hoc estimators which predominantly rely on low pass filtering to reduce their inherent noise characteristics. In the world of AM-FM image decomposition the 2-D energy operator now provides orientation information in addition to the energy with little additional computation. The connection between the positivity of an energy operator and the subharmonicity of a function has been noted, but there needs to be more research into the relationship. As in previous publications we have utilised complex algebra to effortlessly combine various operations. It is possible to utilise the more general properties of other algebras such as matrix, tensor, quaternion or Clifford and extend the ideas to higher dimensions (see Rieger *et al*. [72] or Felsberg [73] for recent insights into n-D orientation) . It is even possible to use vector algebra, but the simplicity is then lost. The complex operators are easily implemented in image and signal processing systems which utilise complex algebra; virtually any system using discrete Fourier transforms for example.

## 8. Conclusion

We have presented two new operators for estimating orientation in 2-D fringe patterns. The (local) differential 2-D energy operator generalises the 1-D TKEO, whilst the (nonlocal) spiral phase 2-D energy operator generalises the 1-D phase congruency energy estimator, making it automatically steerable. A further generalization unifies all four operators in one standard form, Eq. (10) and Eq. (27). Preliminary analysis of their performance shows that the operators successfully combine the advantages of first and second order as well as linear and quadratic orientation estimators. Indeed, for pure sinusoidal inputs the operators give perfectly uniform energy and orientation estimates. We are not aware of any other estimators that perform in this way. The operators also work remarkably well for fringe patterns with strong curvature, and this is most likely due to the avoidance of low pass filters in the operator design.

## 9. Appendix: Extension of the Analytic Signal to Images

At present the image processing research community is split on the issue of how to extend the analytic signal from 1-D to 2-D. The situation is similar in the optical community except that many researchers have recently realized the advantages of the Riesz transform approach when applied to optical fringe patterns. We would like to clarify some of the issues that lead us to the fundamental definition of the orientation problem in terms of Eq. (1) and Eq .(2).

Essentially there are two approaches: the orthant definition outlined below, and the Riesz transform (or spiral phase) approach that has been covered in the main body of this paper and others [10, 11, 73, 74]. The crucial point is that the orthant approach enforces the validity of Eq. (1) by multi-band filtering; Eq. (2) then automatically results from the one-sided, or orthant, nature of the filters. Our approach, in contradistinction, must assume the validity of Eq. (1), but never explicitly realize Eq. (2).

#### 9.1 Orthant/Hyperquadrant Analytic Image

The quintessential exposition of the orthant analytic image occurs in the 1996 thesis of Havlicek [21] and a subsequent paper [75]. Havlicek appears to be the first image processing researcher to have considered the Riesz-based definition of the Hilbert transform in higher dimensions; referring to the classic work of Stein [76]. Having discussed possible reasons for using the Riesz transform to define an analytic image Havlicek (page 79 [21]) eschews the following properties:

- Uniqueness,
- Cauchy Riemann relations,
- Isotropy, in exchange for:
- Spectrum only supported on half the frequency space (orthants or hyperquadrants),
- Analytic signal is complex.

Many good reasons are given for the final choice of an orthant analytic image. Havlicek then goes on to carefully define an “adjusted” signum function in such a way to avoid the problems of an orthant defined as the product of several 1-D signum functions. Clearly the *n* -dimensional product contains *n* hyperplanes with zero values. If a significant signal component just happens to occur on one of these hyperplanes, then the essential quadrature property of the signum multiplier fails because the signal is zeroed. An adjustment has to be made on a set of measure zero (the hyperplanes). In contrast the Riesz transform approach only has a discontinuity at a single point (DC), rather than *n* hyperplanes, and thus remains essentially continuous, isotropic, and unique. But the Riesz based analytic image requires three separate components, therefore it cannot be represented by a complex (two component) function.

#### 9.2 Quasi Eigenfunction Approximation

The application of the orthant based analytic image requires multi-band decomposition. The idea of multi-band demodulation is to separate an image into sub-bands each of which must be smaller than a single orthant. The theory of quasi eigen-function approximation (QEA) is developed by Havlicek in Chapter 3 [21] and shows that filtering of the sub-band effects the phase and amplitude modulation of each component of the analytic image in an intuitive and well defined way. Interestingly QEA gives the same results as the application of the asymptotic method of stationary phase [11].

The application of multi-band decomposition to the kind of fringe patterns we consider in this paper appears to be possible. For a closed fringe pattern with almost constant fringe spacing (like a fingerprint, for example) the dominant frequency component forms an annulus. Multi-band filtering would likely require at least three bands. Each band-pass analytic image would provide both phase and amplitude information, thus avoiding the uniformity problem of conventional approach outlined in our main text. The phase gradient of each provides an estimate of the orientation, but the multiple bands need to be combined in some manner to give an overall orientation and energy map corresponding to our energy operators.

## Acknowledgments

Michael Oldfield created the software to implement and test the orientation estimators. Peter Fletcher first pointed out the intrinsic quadrant limitations of traditional 2-D energy operators. Michael Felsberg kindly answered a number of my questions regarding his recent innovation: the energy tensor.

## References and links

**1. **G. H. Granlund and H. Knutsson, *Signal processing for computer vision*, Kluwer, Dordrecht, Netherlands, 1995.

**2. **G. Krieger and C. Zetzche, “Nonlinear image operators for evaluation of local intrinsic dimension,” IEEE Trans. Image Process. **5**, 1026–1042 (1996). [CrossRef] [PubMed]

**3. **M. Kass and A. Witkin, “Analyzing oriented patterns,” CVGIP **37**, 362–385 (1987).

**4. **R. Penrose, “The topology of ridge systems,” Ann. Hum. Genet.,Lond. **42**, 435–444 (1979). [CrossRef]

**5. **K. Andresen and Q. Yu, “Robust Phase Unwrapping By Spin Filtering Combined With a Phase Direction Map,” Optik **94**, pp.145–149 (1993).

**6. **Q. Yu and K. Andresen, “Fringe-orientation maps and fringe skeleton extraction by the two-dimensional derivative-sign binary-fringe method,” Appl. Opt. **33**, 6873–6878 (1994). [CrossRef] [PubMed]

**7. **J. L. Marroquin, R. Rodriguez-Vera, and M. Servin, “Local phase from local orientation by solution of a sequence of linear systems,” J. Opt. Soc. Am. A **15**, 1536–1544 (1998). [CrossRef]

**8. **H. Knutsson, R. Wilson, and G. H. Granlund, “Anisotropic Non-Stationary Image Estimation and its Applications ┅ Part I: Restoration of Noisy Images,” IEEE Trans. Commun. **31**, 388╌397 (1983). [CrossRef]

**9. **L. Hong, Y. F. Wan, and A. Jain, “Fingerprint Image Enhancement - Algorithm and Performance Evaluation,” IEEE Trans. Pattern Anal. Mach. Intell. **20**, 777–789 (1998). [CrossRef]

**10. **K. G. Larkin, D. Bone, and M. A. Oldfield, “Natural demodulation of two-dimensional fringe patterns: I. General background to the spiral phase quadrature transform,” J. Opt. Soc. Am. A **18**, pp.1862–1870 (2001). [CrossRef]

**11. **K. G. Larkin, “Natural demodulation of two-dimensional fringe patterns: II. Stationary phase analysis of the spiral phase quadrature transform.,” J. Opt. Soc. Am. A **18**, pp.1871–1881 (2001). [CrossRef]

**12. **K. G. Larkin, “Natural demodulation of 2D fringe patterns,” Fringe′01 - The Fourth International Workshop on Automatic Processing of Fringe Patterns, Bremen, Germany, (2001), Elsevier, The Data Science Library, eds W. Juptner and W. Osten, ISBN : 2-84299-318-7

**13. **H. M. Teager and S. M. Teager, “Evidence for nonlinear sound production mechanisms in the vocal tract,” *Speech production and speech modelling*, ed. Hardcastle, W. J. and A. Marchal, (France: NATO Advanced Study Institute, Series D, 1989) pp. 55.

**14. **J. F. Kaiser, “On a simple algorithm to calculate the 'energy' of a signal,” Proc IEEE Int. Conf. Acoust. Speech, Signal Processing, Albuquerque, NM , (1990), pp. 381–384.

**15. **T.-H. Yu and S. K. Mitra, “A novel nonlinear filter for image enhancement,” Image Processing algorithms and Techniques II, Proc. SPIE **1452** ,(1991), pp. 303–309.

**16. **S. K. Mitra, H. Li, I-S. Lin, and T-H. Yu, “A new class of nonlinear filters for image enhancement,” Int. Conf. Acoustics, Speech, and Signal Processing, Toronto, Canada , (1991), pp. 2525–2528.

**17. **P. Maragos, A. C. Bovik, and T. F. Quatieri, “A multidimensional energy operator for image processing,” SPIE Conference on Visual Communications and Image Processing, Boston, MA , (1992), pp. 177–186.

**18. **P. Maragos and A. C. Bovik, “Image Demodulation Using Multidimensional Energy Separation” J. Opt. Soc. Am. A **12**, pp.1867–1876 (1995). [CrossRef]

**19. **K. G. Larkin, “Topics in Multi-dimensional Signal Demodulation,” PhD thesis. Dept. of Physical Optics, University of Sydney , 2001. http://setis.library.usyd.edu.au/~thesis/adt-NU/public/

**20. **O. Wilde, “Act 1,” *The Importance of Being Earnest*, 1899.

**21. **J. P. Havlicek, “AM-FM Image models,” PhD thesis. University of Texas , 1996.
http://hotnsour.ou.edu/joebob/PdfPubs/JPHavlicekDiss.pdf

**22. **D. Gabor, “Theory of communications,” Journal of the IEE , **93**, 429–457 (1947).

**23. **B. Jahne, *Practical handbook on Image processing for Scientific applications*, CRC Press, Boca Raton, Florida , 1997.

**24. **C. F. Shu and R. C. Jain, “Direct Estimation and Error Analysis For Oriented Patterns,” CVGIP-Image Understanding **58**, 383–398 (1993). [CrossRef]

**25. **M. Felsberg and G. H. Granlund, “Detection using Channel Clustering and the 2D Energy Tensor,” Pattern Recognition: 26th DAGM Symposium, Tübingen, Germany , (2004), pp. 103–110.

**26. **M. Felsberg and E. Jonsson, “Energy tensors: Quadratic phase invariant image operators,” DAGM Symposium, Mustererkennung, Wien , (2005),

**27. **H. Canabal, J. A. Quiroga, and E. Bernabeu, “Automatic processing in moire deflectometry by local fringe direction calculation,” Appl. Opt. **37**, 5894–5901 (1998). [CrossRef]

**28. **X. Zhou, J. P. Baird, and J. F. Arnold, “Fringe-Orientation Estimation by use of a Gaussian Gradient Filter and Neighboring-Direction Averaging,” Appl. Opt. **38**, 795–804 (1999). [CrossRef]

**29. **J. A. Quiroga, M. Servin, and F. Cuevas, “Modulo two pi fringe orientation angle estimation by phase unwrapping with a regularized phase tracking algorithm,” J. Opt. Soc. Am. A **19**, 1524–1531 (2002). [CrossRef]

**30. **H. Knutsson and M. Andersson, “Robust N-Dimensional Orientation Estimation using Quadrature Filters and Tensor Whitening,” ICASSP ′94, Adelaide, Australia , (1994). http://www.cvl.isy.liu.se/ScOut/Publications/PaperInfo/ka94.html

**31. **P.-E. Danielsson, Q. Lin, and Q.-Z. Ye, “Efficient detection of second degree variations in 2D and 3D images,” Journal of Vis. Commun. Image Represent. **12**, 255–305 (2001). [CrossRef]

**32. **J. P. Da Costa, F. Le Pouliquen, C. Germain, and P. Baylou, “New operators for optimized orientation estimation,” ICIP 2001, Thessaloniki, Greece , (2001).

**33. **M. C. Morrone and D. C. Burr, “Feature detection in human vision: a phase-dependent energy model,” Proceedings of the Royal Society of London, B **235**, 221–245 (1988). [CrossRef]

**34. **E. H. Adelson and J. R. Bergen, “Spatiotemporal energy models for the perception of motion,” J. Opt. Soc. Am. A **2**, 284–299 (1985). [CrossRef] [PubMed]

**35. **J. G. Daugman and C. J. Downing, “Demodulation, predictive coding, and spatial vision,” J. Opt. Soc. Am. A **12**, 641–660 (1995). [CrossRef]

**36. **P. Maragos, T. F. Quatieri, and J. F. Kaiser, “Speech nonlinearities, modulations, and energy operators,” Proc IEEE Int. Conf. ASSP, Toronto, Canada , (1991), 421–424.

**37. **P. Maragos, T. F. Quatieri, and J. F. Kaiser, “On separating amplitude from frequency modulations using energy operators,” Proc IEEE Int. Conf. ASSP, San Francisco, CA , (1992), 1–4.

**38. **P. Maragos, J. F. Kaiser, and T. F. Quatieri, “On amplitude and frequency demodulation using energy operators,” IEEE Trans. Sig. Process . **41**, 1532–1550 (1993). [CrossRef]

**39. **A. C. Bovik, P. Maragos, and T. F. Quatieri, “AM-FM energy detection and separation in noise using multiband energy operators,” IEEE Trans. Sig. Process. **41**, 3245–3265 (1993). [CrossRef]

**40. **A. C. Bovik and P. Maragos, “Conditions for positivity of an energy operator,” IEEE Trans. Sig. Process. **42**, 469–471 (1994). [CrossRef]

**41. **A. Potamianos and P. Maragos, “A Comparison of the Energy Operator and the Hilbert Transform Approach to Signal and Speech Demodulation,” Signal Process. **37**, 95–120 (1994). [CrossRef]

**42. **K. G. Larkin, “Efficient nonlinear algorithm for envelope detection in white light interferometry,” J. Opt. Soc. Am. A **13**, 832–843 (1996). [CrossRef]

**43. **K. G. Larkin, “Efficient Demodulator for Bandpass Sampled AM Signals,” Electron. Lett. **32**, 101–102 (1996). [CrossRef]

**44. **P. A. Fletcher and K. G. Larkin, “Direct Embedding and Detection of RST Invariant Watermarks,” IH2002, Fifth International Workshop on Information Hiding, Noordwijkerhout, The Netherlands, (2002), 129–144.

**45. **J. Bigun, T. Bigun, and K. Nilsson., “Recognition by symmetry derivatives and the generalized structure tensor,” IEEE Trans. Pattern Anal. Mach. Intell. **26**, (2004). [CrossRef] [PubMed]

**46. **B. Forster, T. Blu, and M. Unser., “A New Family of Complex Rotation-Covariant Multiresolution Bases in 2D,” Proceedings of the SPIE Conference on Mathematical Imaging: Wavelet Applications in Signal and Image Processing X, San Diego CA, USA , (2003), 475–479.

**47. **M. R. Arnison, K. G. Larkin, C. J. R. Sheppard, and N. I. Smith, et al., “Linear phase imaging using differential interference contrast microscopy,” J. Microsc. **214**, 7–12 (2004). [CrossRef] [PubMed]

**48. **R. Hamila, J. Astola, M. A. Cheikh, and M. Gabbouj, et al., “Teager energy and the ambiguity function,” IEEE Trans. Signal Process. **47**, 260–262 (1999). [CrossRef]

**49. **M. Servin, J. A. Quiroga, and J. L. Marroquin, “General n-dimensional quadrature transform and its application to interferogram demodulation,” J. Opt. Soc. Am. A **20**, 925–934 (2003). [CrossRef]

**50. **L. Hormander, *Introduction to complex analysis in several variables*, North Holland, Amsterdam, 1973.

**51. **H. Knutsson, “Filtering and Reconstruction in Image Processing,” PhD. Linkoping University , 1982.

**52. **S. Thurnhofer, “Two-dimensional Teager filters,” *Nonlinear Image Processing*, ed. S. K. Mitra and G. L. Sicuranza, (San Diego: Academic Press, 2001) 167–202. [CrossRef]

**53. **P. L. Duren, *Theory of Hp Spaces*, Dover Publication, Mineola, NY, 2000.

**54. **S. Venkatesh and R. Owens, “On the classification of image features,” Pattern Recogn. Lett. **11**, 339–349 (1990). [CrossRef]

**55. **P. Kovesi, “Image Features From Phase Congruency,” Videre: A Journal of Computer Vision Research, MIT Press **1**, (1999). http://mitpress.mit.edu/e-journals/Videre/001/v13.html

**56. **R. N. Bracewell, *The Fourier transform and its applications*, McGraw Hill, New York, 1978.

**57. **V. Havin and B. Joricke, *The uncertainty principle in harmonic analysis*, Springer-Verlag, Berlin, 1994.

**58. **J. P. Havlicek, D. S. Harding, and A. C. Bovik, “Multicomponent Multidimensional Signals,” Multidimens. Syst. Signal Process . **9**, 391–398 (1998). [CrossRef]

**59. **A. C. Bovik, J. Havlicek, M. Desai, and D. Harding, “Limits on discrete modulated signals,” IEEE Trans. Signal Process . **45**, 867–879 (1997). [CrossRef]

**60. **M. Felsberg and U. Köthe, “Get: The connection between monogenic scale-space and gaussian derivatives,” Scale Space and PDE Methods in Computer Vision, ed. R. Kimmel, N. Sochen, and J. Weickert, LNCS, Springer, 2005) 3459: 192–203. [CrossRef]

**61. **
FFTW: Fastest Fourier Transform in the West., “FFT Benchmarks,” (2005). http://www.fftw.org/speed/

**62. **M. Felsberg, “The GET Operator,” Dept. EE, Linkoping University, 2004.http://www.cvl.isy.liu.se/ScOut/Publications/

**63. **K. G. Larkin, “A self-calibrating phase-shifting algorithm based on the natural demodulation of two-dimensional fringe patterns,” Opt. Express **9**, 236–253 (2001). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-5-236 [CrossRef] [PubMed]

**64. **M. Alonso and G. W. Forbes, “Measures of spread for periodic distributions and the associated uncertainty relations,” Am. J. Phys. **69**, 340–347 (2000).

**65. **B. Strobel, “Processing of Interferometric Phase Maps As Complex-Valued Phasor Images,” Appl. Opt. **35**, 2192–2198 (1996). [CrossRef] [PubMed]

**66. **J. Burke and H. Helmers, “Complex Division as a Common Basis for Calculating Phase Differences in Electronic Speckle Pattern Interferometry in One Step,” Appl. Opt. **37**, 2589–2590 (1998). [CrossRef]

**67. **
NIST Image Group's Fingerprint Research, Fingerprint Test Data on CD-ROM. http://www.itl.nist.gov/iad/894.03/fing/fing.html

**68. **J. M. Huntley, “Fringe analysis today and tomorrow,” Speckle Metrology 2003, Trondheim, Norway, (2003), 167–174.

**69. **J. G. Daugman, “Image analysis and compact coding by oriented Gabor primitives,” Image Understanding and Man Machine interface , 19–30, (1987).

**70. **J. J. Koenderink and W. Richards, “Two-dimensional curvature operators,” J. Opt. Soc. Am. A **5**, 1136–1141 (1988). [CrossRef]

**71. **J. V. d. Weijer, L. J. v. Vliet, P. W. Verbeek, and M. v. Ginkel, “Curvature estimation in oriented patterns using curvilinear models applied to gradient vector fields,” IEEE Trans. Pattern Anal. Mach. Intell. **23**, 1035–1042 (2001). [CrossRef]

**72. **B. Rieger and J. van Vliet Lucas, 22(6), 2004, “A systematic approach to nD orientation representation,” Image Vision Comput. 22, 453–459 (2004). [CrossRef]

**73. **M. Felsberg, “Low-Level Image Processing with the Structure Multivector,” PhD thesis. Christian-Albrechts-University of Kiel , 2002. http://www.isy.liu.se/~mfe/Diss.ps.gz

**74. **M. Felsberg and G. Sommer, “The monogenic signal,” Institut fur informatik, Chritian-Albrechts-Universitat, 2001.

**75. **J. P. Havlicek, J. W. Havlicek, and A. C. Bovik, “The Analytic Image,” IEEE International Conference on Image Processing, Santa Barbara,California, (1997), 446–449.

**76. **E. M. Stein, *Singular integrals and differentiability properties of functions*, Princeton University Press, Princeton, N.J., 1970.