## Abstract

Existing edge detection filters work well on straight edges but make significant errors near sharp corners by producing rounded corners. This is due to the fact that the edge maps produced by these filters are scale variant. We enhance Canny’s optimality criteria to incorporate detection performance near corners as an explicit design objective. The resulting optimal filter, termed ‘Bessel integral filter’, can be derived analytically and exhibits superior performance over recent alternatives, both in terms of numerical accuracy and experimental fidelity. A noise-free localization index is also derived here to account for the detection accuracy of discontinuities forming sharp corners in the absence of noise. We prove here that edges detected by the filters that are not optimal with respect to this noise-free localization index are scale variant. However, the Bessel integral filter proposed here is optimal with respect to the noise-free localization index and therefore it is a scale-invariant filter.

## 1. Introduction

One of the popular approaches in edge detection is to detect discontinuities by convolving a linear filter with an image. It is a common practice in the literature to design a linear filter in one dimension for signals by modelling the edge as a one-dimensional Heaviside function and then extend it for a two-dimensional image by simply adding one dimension to the calculated one-dimensional linear filter. Such an approach in filter design leads to edge detection algorithms, which do not account for two-dimensional features in edges, such as sharp corners and curved discontinuities and therefore may demonstrate a behaviour known as scale-space smoothing causing edges to deform, distort and smooth. Such a behaviour is a consequence of the fact that these filters are scale variant (Lindeberg 1998*a*,*b*). The amount of feature distortion depends on the size (scale) of the filter. This smoothing effect distorts the original shape of objects in images. It is, therefore, important to design filters to detect edges regardless of the filter size (scale). Our contribution in this paper is to introduce a new framework to design scale-invariant filters for images with two-dimensional features. To this end, we exploit some basic principles similar to the ones introduced by Canny (1986). A two-dimensional ideal wedge is defined here to model various types of corner-shaped discontinuities. The problem is then reduced to design a linear filter to detect any point on sharp corners forming discontinuities regardless of the wedge angle and orientation. In his seminal paper, Canny (1986) proposes a framework for computational edge detection by introducing certain criteria such as maximum signal-to-noise ratio (SNR) and localization and least multiplicity of response to derive an optimal filter approximated by a derivative of Gaussian (*DroG*) for edge detection. Deriche proposes a computationally more efficient recursive filter with a constant execution time with infinite support to implement Canny’s *DroG* in Deriche (1987). Tagare & deFigueiredo (1990) exploit the theory of stochastic processes to model the distribution of zero-crossings to calculate an optimal localization index for *DroG*. An effective width is defined for an indefinite length filter in Sarkar & Boyer (1991) to implement a filter that outperforms *DroG*and is computationally more efficient. Canny-like criteria are employed by Petrou & Kittler (1991) and Petrou (1995) to propose an optimal filter for ramp edges by arguing that edges in real world images are closer to ramps. Infinite symmetric exponential filter is also proposed by Shen & Castan (1992) to demonstrate that such a filter has a better localization property than *DroG* filter. A Wiener filter is proposed in Rao & Ben-Arie (1994) by optimizing discriminative signal-to-noise ratio (DSNR) instead of SNR to demonstrate superior performance in terms of Pratt’s figure of merit (Pratt 2007). DSNR is also employed by Wang *et al.* (1996) to propose an optimal filter for bipolar ramp edges demonstrating less distortion and higher figures of merit than the optimal filters for ramp edges proposed in Petrou (1995) and Petrou & Kittler (1991). Demigny (2002) employs Canny’s criteria for the optimal design of filters for edge detection in a discrete setting. Canny-like criteria are employed directly in a two dimensional domain by Jacob & Unser (2004) to design optimal steerable filters initially introduced in Freeman & Adelson (1991) by using Gaussian window functions. The technique of scale multiplication is proposed in Bao & Wu (2005) to employ multiple scales to enhance the localization criterion in the *DroG* filter for edge detection. The rest of the paper consists of the following sections. In §2, we initially define the Canny-like criteria, then derive an optimal filter in a two-dimensional domain and finally discuss its mathematical properties. Section 3 deals with numerical evaluations of Canny-like criteria for the optimal filter proposed here and some popular filters in the literature. Experimental results are presented in §4 and finally the paper concludes in §5.

## 2. Two-dimensional optimal filter and its mathematical properties

Our objective here is to design a filter to detect image discontinuities with numerical cost as low as that of a linear convolution. In an image, such discontinuities may appear as edges along straight lines, edges forming sharp corners and curved edges. In order to account for all kinds of aforementioned edges, let us consider a piecewise constant image (an ideal wedge) shown in figure 1 forming a sharp corner in the centre of coordinate system, point *O*(*φ*≠*π*). In this figure, without loss of generality, we can assume that the *x*-axis is one of the wedge sides to simplify our calculations. Point *O* can also be considered on discontinuity along a straight line for *φ*=*π*. It can also be assumed to be on a curved edge for *φ*=*π* by considering point *O* to be on an infinitesimal line on the curve. Therefore, an optimal filter is sought to detect the discontinuity at point *O* regardless of the wedge angle and orientation.

### (a) Filter derivation

In this paper, by considering an ideal wedge as an input signal shown in figure 1, the first Canny’s (1986) assumption that two-dimensional edges locally have a constant cross section in some direction, is violated (e.g. see point *O* in figure 1). This wedge is defined as a piecewise constant image conveniently described in polar coordinates as
2.1
where *H*,*r* and *θ* are the Heaviside step function, radial and angular coordinates, respectively. A two-dimensional image *g* given in equation (2.1) (or figure 1) is assumed to be the input signal modelling discontinuity with any arbitrary angle and orientation. Let us also assume that the input image is contaminated with a zero-mean Gaussian noise. The detection performance and localization indices (Canny 1986) are then optimized by using the method of Lagrange multipliers to find the optimal filter.

Let us now start by calculating the detection performance and localization indices for point *O* at the centre of coordinate system (0,0) in figure 1. Our first assumption is that the desired filter is linear and space invariant. We further assume the impulse response is isotropic, i.e. in polar coordinates our desired filter is a function of only *r* (radial coordinate). In order to calculate the detection performance, the response of the isotropic filter to the discontinuity at point *O* is calculated by a convolution in Cartesian coordinates evaluated in polar coordinates as
2.2
where *,*U*_{g}, *g* and *h* are the two-dimensional convolution operation, the convolved image, the ideal wedge defined in equation (2.1) and the desired filter, respectively. To detect the edge in the sharp corner *O*, we consider the directional derivative of *U*_{g} along the unit vector ** n** with respect to which axis

*x*(in Cartesian coordinates) has a fixed angle

*θ*

_{n}so that

*θ*

_{n}∈(0,

*φ*). This unit vector

**is represented in polar coordinates as , where and are radial and angular unit vectors in polar coordinates, respectively. The gradient magnitude of the convolved image**

*n**U*

_{g}is examined to find pixels representing discontinuities. Since the relationship between the gradient magnitude and the directional derivative along the unit vector

**at point**

*n**O*is where

*α*is a constant angle between

**and**

*n***∇**

*U*

_{g}|

_{O}, the directional derivative of the convolved image is proportional to its gradient magnitude. The directional derivative of

*U*

_{g}along the unit vector

**at point**

*n**O*in figure 1 is therefore calculated as 2.3 where . By integrating with respect to

*θ*,

*H*

_{g}(0,

*θ*

_{n}) can be written as 2.4 Let us now assume that the ideal wedge shown in figure 1 is contaminated with zero-mean Gaussian noise

*N*(

*r*,

*θ*) with variance . The filter response to noise is therefore written as 2.5 Therefore, the mean-squared response of filter ∂

*h*/∂

*n*to the noise at point

*O*is , where

*E*(.) is the statistical expectation. By integrating with respect to

*θ*, the above equation is rewritten as 2.6 Using equations (2.4) and (2.6), the detection performance is written as 2.7

To calculate the localization index, let us now assume that (∂*H*_{g}/∂*n*)(0,*θ*_{n})=0 where ** n** is the unit vector that has an

*θ*

_{n}angle (

*θ*

_{n}∈(0,

*φ*)) with respect to axis

*x*(in Cartesian coordinates). Generally the condition (∂

*H*

_{g}/∂

*n*)(0,

*θ*

_{n})=0 may not hold for any given isotropic filter. In §2

*b*, it is demonstrated that filters which do not meet the abovementioned condition, produce scale-variant edges. For now, we assume that (∂

*H*

_{g}/∂

*n*)(0,

*θ*

_{n})=0, since we seek to detect point

*O*in figure 1. Therefore, if

*H*

_{n}is the response of the filter ∂

*h*/∂

*n*to the noise, then the filter response at the presence of noise, should have a local maximum in a point with coordinates (

*r*

_{0},

*θ*

_{n}), i.e. 2.8 Let us now evaluate the individual terms in equation (2.8). The McLaurin’s series expansion along the unit vector

**is employed to calculate (∂**

*n**H*

_{g}/∂

*n*)|

_{r0,θn}for a small displacement

*r*

_{0}: 2.9 As indicated before, we assume (∂

*H*

_{g}/∂

*n*)|

_{(0,θn)}=0 therefore equation (2.9) for a small displacement

*r*

_{0}is written as 2.10 The terms with higher orders of

*r*

_{0}are ignored, since the displacement

*r*

_{0}is assumed small. Equation (2.8) is, therefore, written as 2.11 Owing to the presence of zero-mean Gaussian noise, it is reasonable to assume that

*r*

_{0}is a zero-mean random variable whose variance is calculated as the expectation of by using equation (2.11), i.e. 2.12 Let us now evaluate the variance of the Gaussian random variable (∂

*H*

_{n}/∂

*n*)|

_{(r0,θn)}| in polar coordinates. Similar to the calculation method in equation (2.6), the variance of the random variable (∂

*H*

_{n}/∂

*n*)|

_{(r0,θn)}| can be written as 2.13 On the other hand, the denominator in equation (2.12) can be calculated by considering the boundary conditions, , and therefore 2.14 The localization index is the reciprocal of equation (2.12). By using equations (2.12)–(2.14), the localization index can, therefore, be written as equation (2.15) 2.15 The problem of finding a filter

*h*maximizing the product

*ΣΛ*evaluated in equations (2.7) and (2.15) can be reduced to extremizing one of the integral terms and keeping the other terms constant in equations (2.7) and (2.15) by using the method of Lagrange multipliers. We also note that the terms in equations (2.7) and (2.15) do not depend on

*h*itself and consist of only the first and second derivatives of

*h*. To find the optimal filter in the space of admissible functions, the functional to extremize can, therefore, be written with respect to

*P*=d

*h*/d

*r*, i.e. 2.16 where

*Γ*(

*P*,

*P*

_{r},

*r*) is the Lagrangian,

*γ*,

*α*,

*β*,

*ω*are constant coefficients and

*P*

_{r}=

*dP*/d

*r*. We therefore arrive at the Euler–Lagrange equation: 2.17 where

*α*′=

*α*/

*β*,

*γ*′=

*γ*/

*β*and

*ω*′=

*ω*/

*β*. By assuming that

*r*≠0, the above equation is rewritten as 2.18 Equation (2.18) is of the type of modified Bessel differential equation (Polyanin & Zaitsev 2003). Solutions for differential equation (2.18) exist, only if (

*γ*/

*α*)=(

*ω*/

*β*). The necessary condition for functional (2.16) to have a minimizer is that its second variation is positive. The second variation for functional (2.16) can then be calculated as 2.19 where

*Γ*

_{PP}=(∂

^{2}

*Γ*/∂

*P*

^{2}), ,

*Γ*

_{PPr}=(∂

^{2}

*Γ*/∂

*P*∂

*P*

_{r})=0,

*q*(

*r*)=

*δP*(

*r*) and

*q*

_{r}=d

*q*/d

*r*are from the space of admissible functions. The second variation calculated above is positive if inequalities (2.20) are satisfied for ∀

*r*∈

*R*

^{+}: 2.20 In order for inequalities (2.20) to be true for any

*r*∈

*R*

^{+}, it is required that

*α*>0,

*β*>0. The general solution of equation (2.18) is therefore written as , where

*I*

_{1}(

*r*) and

*K*

_{1}(

*r*) are modified Bessel functions of the first and second kind and degree one (Polyanin & Zaitsev 2003) and

*A*,

*B*and

*C*are constants, respectively. It is noted that

*α*′=

*α*/

*β*is positive since both

*α*and

*β*are positive. This solution is subject to the boundary conditions and

*P*(0)=0. We also note that

*P*(

*r*) is required to be asymmetric along the unit vector

**. Therefore,**

*n**P*(

*r*) along the unit vector

**is written as 2.21**

*n*Figure 2*a* depicts *P*(*r*) for *α*′=0.2 along the unit vector ** n**. It is clear that as ,

*P*(

*r*) approaches to infinity. For implementation purposes, this is not computationally tractable. Therefore, a regularization procedure is proposed here as follows: 2.22 The isotropic filter

*h*is then calculated by integrating the above function with respect to

*r*. 2.23 In this paper, we refer to isotropic filter (2.23) as

*Bessel Integral*filter. Figure 2

*b*,

*c*shows a three-dimensional view of the isotropic filter

*h*and its cross section for

*α*′=0.2 and

*ε*=

**Δ**=1, respectively, where

**Δ**is the sampling distance.

### (b) Mathematical analysis for edge detection filters and algorithms

In the first part of this section, we present some mathematical properties of the Bessel Integral filter. In the second part of §2*b*, we describe general properties of scale-variant and invariant filters to demonstrate that the Bessel integral filter belongs to a family of scale-invariant filters. In the following theorem, we initially prove that the local exterma of the gradient magnitude of the convolved image calculated in equation (2.2) correspond to the discontinuities of the piecewise constant image *g*.

### Theorem 2.1

*Let* *is the image domain) be the image obtained in equation (2.2) by convolving the regularized function h*_{ε} *in equation (2.23) and a piecewise constant image g shown in* figure 1*. In a neighbourhood of point O, if the gradient magnitude of U*_{g} *has a local extermum in (r*_{0}*,θ*_{0}*) so that r*_{0} *is a small displacement from point O and θ*_{0}*≠0,φ, then this local extermum approaches point O, as ε→0.*

A proof for theorem 2.1 is provided in appendix A. In theorem 2.2 presented below, we prove that the gradient magnitude of the image convolved by the input image and the isotropic filter *h* is infinity on the discontinuities of the original image when the regularizing parameter *ε* approaches zero.

### Theorem 2.2

*The gradient magnitude of the convolved image U*_{g}*:Ω→R*^{+} *between the regularized isotropic function h*_{ε} *calculated in equation (2.23) and a piecewise constant image g approaches infinity on the discontinuities of the input image g, if* *.*

A proof for the above theorem is presented in appendix B. Proposition 2.3 presented below readily follows from theorems 2.1 and 2.2. This proposition is the basis for an edge detection algorithm by finding the local maxima of the gradient magnitude of the input image convolved with the Bessel integral filter.

### Proposition 2.3

*The gradient magnitude of the convolved image* *U*_{g}:*Ω*→*R*^{+} *between the regularized isotropic filter* *h*_{ε} *calculated in equation (2.23) and a piecewise constant image g has local maxima on the discontinuities of the input image g, if* .

Given that the gradient magnitude is always positive, the proof of proposition 2.3 is straightforward and is therefore skipped here. According to proposition 2.3, if filter (2.23) is convolved with any input image, the local maxima of the gradient magnitude of the convolved image correspond to discontinuities of the input image. Proposition 2.3 is therefore a fundamental statement for our discontinuity detection algorithm proposed in this paper. The result obtained from theorem 2.1 leads also to proposition 2.4. This proposition rationalizes the assumption made to arrive at equation (2.10) from equation (2.9) for the Bessel integral filter.

### Proposition 2.4

*For the Bessel integral filter and the input image g shown in* figure 1, *equation (2.9) approaches equation (2.10) as* *if* *r*_{0} *is a small displacement.*

The proof of this proposition is straightforward and immediate from the result obtained from theorem 2.1 and is therefore skipped. According to proposition 2.4, localization calculated in equation (2.15) is only true for the Bessel integral filter. However, the localization index for any other isotropic edge detection filter is calculated by using equations (2.8) and (2.9). From these equations, the reciprocal of the square root of variance of *r*_{0} being the localization index is calculated as
2.24
In the absence of noise, a noise-free localization index is derived as
2.25
Equations (2.24) and (2.25), theorem 2.1 and proposition 2.4 lead to lemma 2.5.

### Lemma 2.5

*Let the function* *h*(*r*):*R*^{+}→*R* *continuous in* , *be a two-dimensional isotropic edge detection filter in a polar coordinate system. If the central value of such a filter* *h*(0) *is bounded, the two-dimensional noise-free localization index calculated in equation (2.25) for the wedge angle* *φ*≠*π*, *has a bounded value.*

Lemma 2.5 is proved in appendix C. Lemma 2.5 in fact indicates what kind of filters produces scale-variant edges and it therefore leads to lemma 2.6 describing what property a filter should possess to be able to produce scale-invariant edges. In the proof of lemma 2.5, it is demonstrated that equation (2.25) for an isotropic filter *h*(*r*) (*r* is the radial coordinate in the polar coordinate system) can be written as
2.26
It can also be verified that for a Gaussian filter and in the absence of noise, the localization index calculated in equation (2.26) is proportional to the reciprocal of the standard deviation (scale) of the Gaussian filter, i.e. if the scale of the Gaussian filter increases, the noise-free localization of the detected edge corresponding to the corner point *O* in figure 1 decreases. This phenomenon observed in the case of images with discontinuities containing sharp corners such as the one shown in figure 1 is related to the scale variance property of a Gaussian filter. The conditions under which an isotropic filter detects scale-invariant edges are indicated in lemma 2.6 concluded from the results obtained from theorem 2.1 and lemma 2.5.

### Lemma 2.6

*Let the function* *h*(*r*):*R*^{+}→*R* *continuous in* , *be a two-dimensional isotropic edge detection filter in a polar coordinate system. The localization index in the absence of noise for such a filter is infinity regardless of the value of wedge angle and its orientation as shown in* figure 1, *if* *h*(*r*) *approaches infinity as* .

The proof of this lemma is similar to the proof of theorem 2.1 and can be given by evaluating the localization index calculated in equation (2.26) in lemma 2.5 as and it is therefore skipped in this paper. Lemma 2.6 can be considered as a direct consequence of proposition 2.3 for the Bessel integral filter as a member of a family of filters for which this lemma is applicable.

## 3. Numerical evaluations of detection performance and localization

In this section, the detection performance and localization indices are numerically evaluated with respect to the parameters of the Bessel integral filter derived in the previous section. These indices are evaluated here to compare with some existing filters in the literature. The filter *P*_{ε}(*r*) regularized in equation (2.22) is employed to numerically evaluate indices with respect to the regularizing parameter *ε* and *μ*=1/*α*′. The parameter *μ* is used here instead of *α*^{′} for convenience. The detection performance and localization indices are calculated according to Canny (1986). Figure 3 depicts the detection performance and localization indices and their product with respect to *μ* and the regularizing parameter *ε*. As shown in this figure, by lowering the regularizing parameter *ε*, the detection performance index decreases as depicted in figure 3*a*. The localization index greatly improves as *ε* is lowered (figure 3*b*). The product of these two indices improves as *ε* is lowered (figure 3*c*). We also notice that the ability of the filter for noise removal increases by raising the values for *μ* (and as a result the filter size). The localization, on the other hand, almost remains unchanged, if *μ* is increased (or decreased) as demonstrated in figure 3*b*. We now need to determine the window size according to which the regularized Bessel integral filter is truncated, since this filter in spatial domain is not band limited. Therefore, we define the truncated filter as
3.1
The normalized truncation error *E*(*d*) defined in equation (3.2) as
3.2
We notice that the normalized truncation error *E*(*d*) for *ε*∈[0.01,1] depends only on . For , the error term *E* becomes negligible (*E*=0.0265). Therefore, in this paper, the window size of the filter proposed here is set to . Figure 4 depicts the detection performance, localization indices and their product with respect to the filter size for *DroG*, Demigny and Bessel integral filters with *ε*=0.01, 0.1 and 1. As shown in the figure, the localization of the filter proposed here increases and its detection performance decreases, as lower values for *ε* are selected, which confirms the result shown in figure 3. The localization of the Bessel integral filter for all values of *ε* is superior to the localization of the other two filters (figure 4*b*). Finally, figure 4*c* depicts the product of the detection performance and localization. The Bessel integral filter with *ε*=0.01 has the best overall performance (product of indices) over the whole range of filter sizes.

## 4. Experimental results

Before the experimental results are discussed in this section, there are two implementation issues that require more clarification. Firstly, if the regularizing parameter less than unity (e.g. as low as *ε*=0.01) is required anywhere in this paper, the filter is initially constructed in a higher resolution (e.g. **Δ**=0.01) with the required *ε* and then down-sampled to the normal resolution **Δ**=1. Secondly, similar to other filters, the Bessel integral filter is also convolved with the original image, and then the non-maximum suppression algorithm is applied on the gradient magnitude of the convolved image to detect discontinuities. Let us now compare the algorithm proposed here with *DroG* (Canny 1986), Demigny (2002) and Jacob & Unser (2004) edge detection algorithms. The three Jacob–Unser filters employed in this paper are a wedge filter with wedge angle *π*/2, and two third-order edge detectors with filter smoothing parameters 0.09 and 0.2 (named as edge 1 and edge 2 filters, respectively, in this paper) (see Jacob & Unser (2004) table 1 for more details). The six filters are applied on the synthetic star image shown in figure 6 to compare the accuracy (localization) of the detected edges. To measure the distance between each detected edge pixel from the corresponding real edge pixel on the original image, we use signed distance function (SDF) (Sethian 1996) of the original image of figure 6. To compute an error term representing the average error between the true and the detected edge pixels, let us assume that *z*=*ϕ*(*x*,*y*) is the SDF of the synthetic star image calculated by using the fast marching algorithm (Sethian 1996) and *x*_{i} and *y*_{i} are the *x* and *y* coordinates of the *i*th edge pixel of an edge map produced by one of the algorithms investigated here. The average error term associated with the edge map containing *N*edge pixels is measured by . Same filter size is used for all filters in each measurement. Figure 5 shows the average error term for each algorithm for filter sizes 5, 9, 13, 17, 21, 41 and 101. By changing the filter sizes for each measurement, other filter parameters are also changed accordingly. For example, for *DroG* and Jacob–Unser algorithms, standard deviation=filter size/8, and for the Bessel integral filter *μ*=(filter size)^{2}/64 are calculated. As depicted in figure 5, the average error term increases for all scale-variant filters. However, the average error of edge detection for the Bessel integral filter is unchanged and lowest as implied from theorems 2.1, 2.2, proposition 2.3 and lemma 2.6 proved in §2. The edge maps detected by the six filters for the filter size=41 is depicted in figure 6. As shown in this figure, the sharp corners are detected properly only with the Bessel integral filter regardless of filter size and wedge angle. The other filters with the same size have smoothed the sharp corners. We notice that the central processing unit (CPU) time required to produce edge maps is the same for all filters with the same size. Because all edge detection algorithms investigated here are based on the calculation of the gradient of the convolved image with a linear filter. For example, it takes around 0.1 s for a PC with a 2.6 GHz CPU to run every algorithm discussed in the section with filter size 41×41 in a Matlab environment (v. 7.3) on the original image of figure 6 for edge detection. In order to choose equivalent threshold levels in all algorithms investigated in this paper for a fair comparison, a strip of logarithmically increasing grey scales is added to the bottom of the noisy and real world images used here. This strip of grey scales named as threshold scaling strip in this paper and its cross section are shown in figure 7. The threshold levels in all algorithms are adjusted so that the same number of discontinuities in the strip is detected by all algorithms. Figure 8 shows a real world image of a text book. Noise is introduced to this image in the data acquisition process. The Bessel integral, Jacob–Unser edge 1 and Jacob–Unser wedge filters (size=9×9) are applied to the image of this figure. The threshold levels for all algorithms are adjusted so that the first edge in the threshold scaling strip is detected. A close inspection of this figure indicates that the edge map produced by the Bessel integral filter has the best accuracy (e.g. see ‘t’ in ‘that’, ‘transition’ and ‘gradient’, ‘r’ in ‘gradient’ and ‘propose’, and so on). As the filter size increases, the scale-space distortion of Jacob–Unser filters becomes more severe. However, the Bessel integral filter experiences less scale-space distortion. The Bessel integral and Jacob–Unser edge 1 filters (filter size=13×13) are applied to the noisy image of Cameraman as shown in figure 9*a*. This image is contaminated with the zero-mean Gaussian noise with standard deviation 50. The edge maps produced by the Bessel integral and Jacob–Unser filters are depicted in figure 9*b* and *c*, respectively. As shown in this figure, the Bessel integral filter produces competitive results in comparison with Jacob–Unser filter in the presence of excessive noise. Soft edge maps, as suggested in Martin *et al.* (2004), are produced by applying Bessel integral and Jacob–Unser edge 1 filters (with size 13) on a real world image from Berkeley segmentation data base (Martin *et al*. 2004), as shown in figure 10*a*. There are many details that are detected with higher accuracy by Bessel integral filter. For example, number ‘96’ on the front car and the signs ‘CASTROL’ in the background (figure 10*a*) are completely readable in figure 10*b*. However, it is very difficult (if not impossible) to read these details in the soft edge map produced by Jacob–Unser filter as shown in figure 10*c*. The source code of the edge detection algorithm proposed in this paper is written in Matlab (v. 7.3) and is available in http://users.ecs.soton.ac.uk/sm3/curricul2.htm.

## 5. Conclusion

A scale-invariant filter termed as Bessel integral filter is derived here by defining the Canny-like criteria in a two-dimensional image domain. The Bessel integral filter proposed here is singular at the centre. A regularization method for the filter is therefore proposed here to improve the robustness of this edge detection method in the presence of noise. It is demonstrated that the product of the detection performance and localization indices of the filter proposed here is higher than some recent filters in the literature. It is also analytically proved in this paper that the gradient amplitude of the filtered image has local maxima on the discontinuities of the input image, when the regularizing parameter of the Bessel integral filter approaches zero. This property of the Bessel integral filter is the main reason for its scale-invariance property in edge detection. The condition under which a filter may be scale variant or invariant is presented here. Two issues, not discussed in this paper, are open for further research: (i) the possibility of employing the filter proposed here in a class of steerable filters as in Jacob & Unser (2004); (ii) the filter proposed here is designed in a continuous domain. A discrete treatment as in Demigny (2002) for the filtering design introduced here would lead to a discrete filter equivalent to the Bessel integral filter derived here in a continuous domain.

## Acknowledgements

This work was supported in part by the IST program of the European Community, under the PASCAL2 Network of Excellence, the IST-2007-216886 and PinView project with grant number 216529.

- Received September 10, 2010.
- Accepted December 15, 2010.

- This journal is © 2011 The Royal Society

## Appendix A

### Proof for theorem 2.1.

Let us consider the piecewise constant image *g* as shown in figure 1. Let us further assume that the local extermum of *H*_{g}=∂*U*_{g}/∂*n* corresponding to point *O* in the original image is at coordinate (*r*_{0},*θ*_{0}), where *r*_{0} is a small displacement, *θ*_{0}≠0,*φ* and . By using the McLaurin’s series, we can write:
A1
where is perpendicular to ** n**. The higher order terms are ignored owing to the fact that

*r*

_{0}is a small displacement. Since

*H*

_{g}(

*r*

_{0},

*θ*

_{0}) is a local extermum, (∂

*H*

_{g}/∂

*n*)|

_{(r0,θ0)}=0. Equation (A1) can therefore be used to calculate

*r*

_{0}the distance between the local extermum and point

*O*, i.e. A2 From equation (2.2), straightforward calculations show that A3 where

*g*is defined in equation (2.1). By integrating the above equation with respect to

*θ*, and then integrating by parts with respect to

*r*for the first term in the above equation and finally considering the asymptotic behaviour of the Bessel function, we can write: A4 It is also straightforward to show that A5

By substituting equations (A4), (A5) and (2.14) into equation (A2), we can write:
A6
where . Equation (A6) is indeterminate as *ε*→0. By using l’Hôpital rule for intermediate terms, we can write equation (A6) as
A7
Therefore, the distance between point *O* and the local extermum *H*_{g}(*r*_{0},*θ*_{0}) approaches zero, as , i.e. *H*_{g}(0,*θ*_{0}) is a local extermum as . Therefore |**∇***U*_{g}| has a local extermum in point *O*. The above argument is also correct for edges along straight lines for *φ*=*π*. Any edge point along curved boundaries can be considered a point on an infinitesimal line for which the above argument also applies. ■

## Appendix B

### Proof for theorem 2.2.

Let us consider the piecewise constant image shown in figure 1 described by equation (2.1) in polar coordinates as an input image *g* in equation (2.2). In figure 1, point *O* is placed on the image discontinuity forming a sharp corner. We aim to show that |**∇***U*_{g}|_{O} approaches infinity as *ε*→0. By using equation (2.2), let us evaluate *H*_{g}=(∂*U*_{g}/∂*n*)|_{(0,θn)}, where , and *θ*_{n}, and are any arbitrary angle between 0 and *φ*, and the unit vectors (radial and angular unit vectors) in polar coordinates, respectively. The Cartesian convolution evaluated in polar coordinates is written as
B1
On the other hand, ∂*h*/∂*n* can be calculated as
B2
where *P*_{ε}(*r*) is defined in equation (2.22). By substituting equation (B2) into equation (B1), then integrating by parts with respect to *r* and finally integrating with respect to *θ*, we can write:
B3
where . Using the asymptotic behaviour of *K*_{1}(*r*) as and (e.g. Polyanin & Zaitsev 2003; Arfken & Weber 2005) and l’Hôpital’s rule for indeterminate terms, equation (B3) can be rewritten as
B4
As , the integration term in the above equation approaches infinity. Since (∂*U*_{g}/∂*n*)(0,0)=**∇***U*_{g}|_{O}⋅** n**, then . We notice that for

*φ*=

*π*, the wedge shown in figure 1 becomes an edge along a straight line on which point

*O*can be regarded as an arbitrary point. It is also noted that if point

*O*is placed on some part of the edge (contour) that can be regarded as a regular curve, then the contour at point

*O*can be approximated as an infinitesimal line passing through point

*O*for both of which the above argument applies. ■

## Appendix C

### Proof for lemma 2.5.

In the absence of noise, the noise-free localization index is calculated according to equation (2.25). We would like to evaluate the noise-free localization by calculating the terms (∂*H*_{g}/∂*n*)|_{(0,θn)} and (∂^{2}*H*_{g}/∂*n*^{2})|_{(0,θn)}:
C1
and
C2
where *h*(*r*) is only a function of radial coordinate in a polar coordinate system. By replacing (C1) and (C2) into equation (2.25), the noise-free localization index is calculated as
C3
Since *h*(0) is bounded, (d*h*/d*r*)|_{(r=0)}=0, and the wedge angle *φ*≠*π*, then the noise-free localization index evaluated in (C3) is bounded. ■