Edge, Ridge, and Blob Detection with Symmetric Molecules

TL;DR

This paper introduces a new multiscale, symmetry-based method for detecting and characterizing edges, ridges, and blobs in images, demonstrating robustness to noise and invariance to contrast changes.

Contribution

The paper presents a novel symmetry-based approach using alpha-molecules for feature detection, providing comprehensive geometric characterization and improved robustness over existing methods.

Findings

High accuracy in synthetic image detection tasks

Robustness to noise and contrast variations demonstrated

Effective in real-world applications like blood vessel and cell colony analysis

Abstract

We present a novel approach to the detection and characterization of edges, ridges, and blobs in two-dimensional images which exploits the symmetry properties of directionally sensitive analyzing functions in multiscale systems that are constructed in the framework of alpha-molecules. The proposed feature detectors are inspired by the notion of phase congruency, stable in the presence of noise, and by definition invariant to changes in contrast. We also show how the behavior of coefficients corresponding to differently scaled and oriented analyzing functions can be used to obtain a comprehensive characterization of the geometry of features in terms of local tangent directions, widths, and heights. The accuracy and robustness of the proposed measures are validated and compared to various state-of-the-art algorithms in extensive numerical experiments in which we consider sets of clean and distorted synthetic images that are associated with reliable ground truths. To further demonstrate the applicability, we show how the proposed ridge measure can be used to detect and characterize blood vessels in digital retinal images and how the proposed blob measure can be applied to automatically count the number of cell colonies in a Petri dish.

Figures

Tables

Table 1. Table 1: Edge Detection Performance on Synthetic Images

	No Noise				Medium Noise				Severe Noise
	$\mathrm{FOM}$	${\mathrm{MAE}}_{O_E}$	$\mathrm{TPR}$		$\mathrm{FOM}$	${\mathrm{MAE}}_{O_E}$	$\mathrm{TPR}$		$\mathrm{FOM}$	${\mathrm{MAE}}_{O_E}$	$\mathrm{TPR}$
synthetic image 1
SymFD	0.91	1.79 ∘	100 %		0.89	3.19 ∘	100 %		0.72	4.46 ∘	83 %
Shearlet [74]	0.85	3.78 ∘	100 %		0.82	5.78 ∘	98 %		0.64	8.71 ∘	84 %
PhaseCong [45]	0.91	6.13 ∘	100 %		0.85	8.48 ∘	95 %		0.58	15.59 ∘	72 %
Canny [21]	0.92	n/a	100 %		0.90	n/a	100 %		0.72	n/a	86 %
Sobel [65]	0.93	n/a	100 %		0.73	n/a	83 %		0.11	n/a	17 %
synthetic image 2
SymFD	0.90	4.55 ∘	100 %		0.85	7.43 ∘	96 %		0.64	13.28 ∘	76 %
Shearlet [74]	0.85	6.58 ∘	100 %		0.79	9.67 ∘	96 %		0.55	13.32 ∘	73 %
PhaseCong [45]	0.87	9.46 ∘	97 %		0.72	12.11 ∘	82 %		0.42	17.01 ∘	52 %
Canny [21]	0.91	n/a	100 %		0.87	n/a	99 %		0.65	n/a	79 %
Sobel [65]	0.93	n/a	100 %		0.65	n/a	76 %		0.16	n/a	22 %

Table 2. Table 2: Evaluation of four ridge detection methods on synthetic images.

	No noise					Medium noise					Severe noise
	$\mathrm{FOM}$	${\mathrm{MAE}}_{W_R}$	${\mathrm{MAE}}_{O_R}$	$\mathrm{TPR}$		$\mathrm{FOM}$	${\mathrm{MAE}}_{W_R}$	${\mathrm{MAE}}_{O_R}$	$\mathrm{TPR}$		$\mathrm{FOM}$	${\mathrm{MAE}}_{W_R}$	${\mathrm{MAE}}_{O_R}$	$\mathrm{TPR}$
Synthetic image 3
SymFD	0.90	0.12 px	1.20 ∘	100 %		0.89	0.26 px	2.15 ∘	100 %		0.81	0.84 px	4.80 ∘	94 %
Steger [66]	0.90	0.71 px	n/a	100 %		0.85	1.31 px	n/a	100 %		0.81	2.56 px	n/a	98 %
Frangi [32]	0.80	n/a	1.24 ∘	99 %		0.80	n/a	2.18 ∘	99 %		0.20	n/a	3.26 ∘	23 %
PhaseCong [45]	0.89	n/a	15.67 ∘	99 %		0.81	n/a	31.39 ∘	96 %		0.48	n/a	36.07 ∘	52 %
Synthetic image 4
SymFD	0.88	0.55 px	4.18 ∘	99 %		0.87	0.59 px	4.28 ∘	98 %		0.84	0.77 px	5.12 ∘	97 %
Steger [66]	0.88	1.12 px	n/a	98 %		0.87	1.80 px	n/a	99 %		0.83	1.85 px	n/a	95 %
Frangi [32]	0.85	n/a	4.30 ∘	98 %		0.85	n/a	4.55 ∘	97 %		0.82	n/a	6.94 ∘	97 %
PhaseCong [45]	0.81	n/a	25.71 ∘	97 %		0.82	n/a	29.18 ∘	96 %		0.64	n/a	33.55 ∘	74 %

Table 3. Table 3: Evaluation of two blob detection methods on synthetic images.

	No noise					Medium noise					Severe noise
	$\mathrm{TP}$	$\mathrm{FP}$	${\mathrm{MAE}}_{\mathrm{ctr}}$	${\mathrm{MAE}}_{W_B}$		$\mathrm{TP}$	$\mathrm{FP}$	${\mathrm{MAE}}_{\mathrm{ctr}}$	${\mathrm{MAE}}_{W_B}$		$\mathrm{TP}$	$\mathrm{FP}$	${\mathrm{MAE}}_{\mathrm{ctr}}$	${\mathrm{MAE}}_{W_B}$
Synthetic image 5 (31 blobs)
SymFD	31	0	2.16 px	2.62 px		31	0	2.16 px	2.95 px		31	0	2.09 px	3.13 px
Circular Hough [7, 27]	31	0	0.99 px	0.61 px		30	0	1.82 px	2.33 px		20	0	2.03 px	8.15 px
Synthetic image 6 (200 blobs)
SymFD	200	0	1.10 px	0.87 px		198	0	1.04 px	0.97 px		155	1	1.30 px	1.19 px
Circular Hough [7, 27]	200	0	0.88 px	0.56 px		194	3	0.99 px	1.04 px		47	1	1.10 px	1.66 px

Table 4. Table 4: Evaluation of Orientation Measurements on the REVIEW Database

	VDIS			KPIS			HRIS			CLRIS
	${\mathrm{MAE}}_{O_R}$	$\mathrm{SR}$		${\mathrm{MAE}}_{O_R}$	$\mathrm{SR}$		${\mathrm{MAE}}_{O_R}$	$\mathrm{SR}$		${\mathrm{MAE}}_{O_R}$	$\mathrm{SR}$
Human Observers
O1	2.69 ∘	100 %		2.40 ∘	100 %		2.82 ∘	100 %		1.62 ∘	100 %
O2	2.97 ∘	100 %		2.52 ∘	100 %		2.41 ∘	100 %		1.61 ∘	100 %
O3	3.47 ∘	100 %		3.18 ∘	100 %		2.70 ∘	100 %		1.51 ∘	100 %
Algorithm
SymFD	3.30 ∘	99 %		2.74 ∘	100 %		2.46 ∘	100 %		2.93 ∘	95 %

Table 5. Table 5: Evaluation of Width Measurements on the REVIEW Database

	VDIS			KPIS			HRIS			CLRIS
	${\mathrm{SD}}_{W_R}$	$\mathrm{SR}$		${\mathrm{SD}}_{W_R}$	$\mathrm{SR}$		${\mathrm{SD}}_{W_R}$	$\mathrm{SR}$		${\mathrm{SD}}_{W_R}$	$\mathrm{SR}$
Human Observers
O1	0.54 px	100 %		0.23 px	100 %		0.29 px	100 %		0.57 px	100 %
O2	0.62 px	100 %		0.21 px	100 %		0.26 px	100 %		0.70 px	100 %
O3	0.67 px	100 %		0.23 px	100 %		0.28 px	100 %		0.57 px	100 %
Algorithms
SymFD	1.20 px	99 %		0.41 px	100 %		0.60 px	100 %		1.44 px	95 %
ARIA [9]	0.95 px	99 %		0.29 px	100 %		0.32 px	100 %		0.95 px	100 %
ETOS [10]	0.80 px	100 %		0.36 px	100 %		0.45 px	100 %		0.53 px	100 %
ESP [4]	0.77 px	100 %		0.33 px	100 %		0.42 px	100 %		1.47 px	93 %
Graph [73]	1.43 px	96 %		0.67 px	99 %		0.57 px	100 %		1.78 px	94 %
Gregson [34]	1.49 px	100 %		0.60 px	100 %		2.84 px	100 %		2.84 px	100 %
HHFW [15]	0.88 px	78 %		0.39 px	96 %		0.93 px	88 %		n/a	0 %
1DG [75]	2.11 px	100 %		0.40 px	100 %		4.14 px	100 %		4.14 px	99 %
2DG [52]	1.33 px	77 %		0.34 px	100 %		6.02 px	99 %		6.02 px	27 %

Table 6. Table 6: Accuracy of Cell Colony Counts

	HD			Webcam
	With Mask	Without Mask		With Mask	Without Mask
SymFD	1.95 %	3.57 %		2.23 %	5.95 %
OpenCFU [33]	1.93 %	3.85 %		2.79 %	n/a
ImageJ [16]	6.64 %	80.93 %		11.11 %	n/a
NICE [24]	9.91 %	20.56 %		13.04 %	n/a

Equations

$$A_{s,\alpha }=\left(\begin{array}{cc}s& 0\\ 0& s^{\alpha }\end{array}\right),s>0,\alpha \in [0,1],$$

$$R_{\theta }=\left(\begin{array}{cc}\cos (\theta )& -\sin (\theta )\\ \sin (\theta )& \cos (\theta )\end{array}\right),\theta \in \mathbb{T},$$

$$\mathbb{P}:={\mathbb{R}}_{>0}\times \mathbb{T}\times {\mathbb{R}}^2,$$

$$m_{\lambda }(x)=s_{\lambda }^{(1+\alpha )/2}{g}^{(\lambda )}\left(A_{s_{\lambda },\alpha }{R}_{{\theta }_{\lambda }}(x-y_{\lambda })\right),$$

$$\left|{\partial }^{\rho }\widehat{g^{(\lambda )}}(\xi )\right|\lesssim \min {\left\{1,s_{\lambda }^{-1}+\left|{\xi }_1\right|+s_{\lambda }^{-(1-\alpha )}\left|{\xi }_2\right|\right\}}^M\cdot (1+{\left|\xi \right|}^2{)}^{\frac{-N_1}{2}}\cdot (1+({\xi }_2{)}^2{)}^{\frac{-N_2}{2}}$$

$$(\mathcal{T}f)(g^{(\lambda )},s_{\lambda },{\theta }_{\lambda },y_{\lambda })=⟨f,m_{\lambda }⟩=s_{\lambda }^{(1+\alpha )/2}\underset{{\mathbb{R}}^2}{\int }f(x)\overline{g^{(\lambda )}\left(A_{s_{\lambda },\alpha }{R}_{{\theta }_{\lambda }}(x-y_{\lambda })\right)}\mathrm{d}x.$$

$${\mathrm{PC}}_{{\psi }^{\text{c}}}(f,x)=\frac{\left|\sum _{a\in A}⟨f,{{\psi }^{\text{c}}}_{a,x}{⟩}_{L^2}\right|}{\sum _{a\in A}\left|⟨f,{{\psi }^{\text{c}}}_{a,x}{⟩}_{L^2}\right|+ϵ},$$

$$G_k(x)=(-1{)}^k{H}_k(x)e^{-x^2}=(-1{)}^k{H}_k(x)G_0(x),$$

$$\widehat{G_k}(\xi )=(2\pi i\xi {)}^k\sqrt{\pi }{e}^{-{\pi }^2{\xi }^2}.$$

$$\widehat{\mathcal{H}f}(\xi )=-i\mathrm{sgn}(\xi )\hat{f}(\xi ),$$

$$n^{\ast }=\underset{n\le N}{\mathrm{argmax}}f(x_n).$$

$$f(x)=c_2(x-x_{n^{\ast }}{)}^2+c_1(x-x_{n^{\ast }})+c_0,\text{for x∈{xn∗−1,xn∗,xn∗+1}}.$$

$$x^{\ast }=x_{n^{\ast }}-\frac{c_1}{2c_2}.$$

$$\text{an edge measure}E(f,y)$$

$$\text{a ridge measure}R(f,y)$$

$$\text{and a blob measure}B(f,y)$$

$$O_x(f,y)$$

$$W_x(f,y)$$

$$H_x(f,y)$$

$${\Psi }^{\mathrm{e}}$$

$${\Psi }^{\mathrm{o}}$$

$${\psi }_{j,y}^{\mathrm{o}}(x)=a^j{\psi }^{\mathrm{o}}(a^j(x-y)),{\psi }_{j,y}^{\mathrm{e}}(x)=a^j{\psi }^{\mathrm{e}}(a^j(x-y)),$$

$$K_{{\psi }^{\mathrm{o}}}={\int }_{-\infty }^0{\psi }^{\mathrm{o}}(x)\mathrm{d}x,$$

$${\breve{E}}^{\text{1D}}(f,y)=\frac{\left|\sum _{j\in J}⟨f,{\psi }_{j,y}^{\mathrm{o}}⟩\right|-\sum _{j\in J}\left|⟨f,{\psi }_{j-j^{\mathrm{e}},y}^{\mathrm{e}}⟩\right|-\beta N_J{K}_{{\psi }^{\mathrm{o}}}}{N_J\underset{j\in J}{\max }\left\{\left|⟨f,{\psi }_{j,y}^{\mathrm{o}}⟩\right|\right\}+ϵ},$$

$$E^{\text{1D}}(f,y)=\max \{0,{\breve{E}}^{\text{1D}}(f,y)\}.$$

$${\breve{E}}^{\text{1D}}(g,y^{\ast })=\frac{N_{J}s{K}_{{\psi }^{\mathrm{o}}}-\beta N_J{K}_{{\psi }^{\mathrm{o}}}}{N_{J}s{K}_{{\psi }^{\mathrm{o}}}}=1-\frac{\beta }{s},$$

$${\Psi }_2^{\mathrm{e}}$$

$${\Psi }_2^{\mathrm{o}}$$

$$\widetilde{{\Psi }_2^{\mathrm{e}}}=\left\{c_1{c}_2{\psi }^{\mathrm{e}}(c_1{x}_1){\psi }^{\mathrm{e}}(c_2{x}_2):{\psi }^{\mathrm{e}}\in {\Psi }^{\mathrm{e}},c_1,c_2\in {\mathbb{R}}_{>0}\right\}\subset L^2({\mathbb{R}}^2)\cap L^1({\mathbb{R}}^2)$$

$$SM(g,\alpha ,a,J,\Theta )=\left\{m_{j,\theta ,y}(x)=a^{j(1+\alpha )}g(A_{a^j,\alpha }{R}_{\theta }(x-y)):j\in J,\theta \in \Theta ,y\in {\mathbb{R}}^2\right\}.$$

$$\left|{\partial }^{\rho }\hat{g}(\xi )\right|\lesssim \min {\left\{1,a^{-j}+\left|{\xi }_1\right|+a^{-j(1-\alpha )}\left|{\xi }_2\right|\right\}}^M\cdot (1+{\left|\xi \right|}^2{)}^{\frac{-N_1}{2}}\cdot (1+({\xi }_2{)}^2{)}^{\frac{-N_2}{2}}$$

$$g_k^{(1)}(x)$$

$$g_k^{(3)}(x)$$

$$\left|{\partial }_1^{{\rho }_1}({\xi }_1{)}^k{e}^{-({\xi }_1{)}^2}\right|=\left|P({\xi }_1)e^{-({\xi }_1{)}^2}\right|\lesssim \left({\left|{\xi }_1\right|}^{k+{\rho }_1}+{\left|{\xi }_1\right|}^{\max \{0,k-{\rho }_1\}}\right)e^{-({\xi }_1{)}^2},$$

Full text

1 Introduction

The correct localization of the significant structures in an image as well as the precise characterization of their geometry are two eminent tasks in digital image processing with an overwhelming number of applications. Since the advent of digital image processing, a large body of research has been devoted to the development and analysis of algorithms for the detection and characterization of features such as edges, ridges, and well-defined two-dimensional shapes. Surprisingly, in many practical situations, even highly optimized implementations of popular methods such as the Canny edge detector [21], approaches that are based on directionally sensitive filters [40], or even multiscale systems of anisotropic functions [74] are not always capable of reliably identifying the features in question. This is often the case for images which are heavily distorted by noise, in which different features are strongly overlapping, or where the geometry of the considered features is characterized by a high variation and irregularities such as corner points. It is also worth noting that feature detection in general is a task where, despite recent advances, computers are still often outperformed by humans.

A main difficulty researchers face when developing methods for feature detection is to identify universal and computationally tractable properties that are characteristic of or even unique to points at which a certain feature is localized. In the case of edge detection, the first and most basic observation was that edges are typically associated with changes in contrast and can thus be identified by considering points with large image gradients. This approach led to the development of early edge enhancement filters such as the Roberts filter [63], the Prewitt filter [59], and the Sobel filter [65], which are still widely applied today. However, any procedure merely considering the magnitude of gradient filters as an indicator is highly sensitive to noise and variations in image illumination. Two famous and influential descriptors which already made a big step towards capturing the structural nature of edges were proposed in the 1980’s by Marr and Hildreth [54] and Canny [21]. Marr and Hildreth observed that edges coincide with zero-crossings in the second directional derivative of the image intensities while Canny developed an algorithm for edge detection which identifies points at which the image gradient reaches a local maximum. The Canny edge detector is also based on a gradient filter, and the Marr-Hildreth-operator applies a discrete Laplace filter. However, both methods are not just evaluating whether a magnitude exceeds a certain threshold but consider local structural properties by testing for zero-crossings and local maxima.

A dimensionless approach to identifying features such as edges and ridges may be found in the so-called local energy model, which was proposed in the late 1980’s and postulates that the location of features coincides with points in an image for which the Fourier components are maximally “in phase” [56, 57, 70]. Dimensionless here means that the values yielded by the measure are not describing a physical quantity such as image intensities. To precisely calculate the degree of phase congruency, it is necessary to optimize a certain value over inputs from a continuous interval for every point in an image. This computationally highly inefficient procedure was later significantly simplified by Kovesi, who showed that an equivalent measure can be formulated in terms of convolutions with differently scaled and oriented complex-valued wavelet filters [45, 46]. The resulting phase congruency measure yields a contrast-invariant feature detector that simultaneously detects edges and ridges and can be efficiently implemented to process digital images. It was one of the inspirations for the measures for edge, ridge, and blob detection presented here and will be explicitly defined later.

In recent years, multiscale systems of analyzing functions developed in the field of applied harmonic analysis have become increasingly successful in efficiently representing, approximating, and analyzing two-dimensional digital images. In particular, systems that are based on anisotropic scaling such as curvelets [19, 20], shearlets [51, 49], or contourlets [28] yield provably in-a-sense optimal error decay rates when approximating certain classes of digital images that are governed by curvilinear and anisotropic structures. These constructions have been considered in a wide range of applications in digital image processing including edge and feature detection [74, 50, 31, 43]. Grohs et al. recently introduced the so-called $\alpha$-molecules framework, which allows for a unified treatment of the analysis of most multiscale systems developed in the field of applied harmonic analysis [35]. In the present work, we derive measures for the detection of edges, ridges, and blobs that are inspired by the notion of phase congruency but realized by making full use of the flexibility of modern construction principles for anisotropic analyzing functions provided within the framework of $\alpha$-molecules.

Saying that the Fourier components are in phase at the location of an edge basically means that locally the considered function is purely defined by odd-symmetric sine components while all even-symmetric cosine components vanish. In the first step, we will derive a novel edge measure that utilizes this observation by testing for the local symmetry properties of a function by considering differently scaled and oriented odd- and even-symmetric analyzing functions constructed in the $\alpha$-molecule framework. This measure can then be generalized to detect ridges by interchanging the roles of the odd- and even-symmetric molecules and taking into account the width of the feature. We further show how, with a few modifications regarding the measure and the construction of analyzing functions, similar principles can be applied to the case of blob detection. Each analyzing function in a multiscale system considered in this work is associated with a specific pair of scaling and orientation parameters and thus conveys not only information about the location of features but also about their geometry in terms of feature width and orientation. Finally, we demonstrate how this information can be utilized to obtain measures for local tangent directions, widths, and heights, that yield detailed information about the geometry of an image.

1.1 Contributions

Earlier variants of the edge measure derived in Section 3.2 and the ridge measure derived in Section 3.3 that are based on complex-valued shearlets have already been published previously by the authors of the present work and their co-authors [44, 60, 62]. The main contribution of this work is to provide a comprehensive derivation and description of the proposed measures and bring them into the framework of $\alpha$-molecules which provides enough flexibility to allow for a unified treatment of the cases of edge, ridge, and blob detection. In particular, we demonstrate in Theorem 3.1 that by considering tensor products of derivatives of the one-dimensional Gaussian and their Hilbert transforms, one can obtain a large class of well-behaved $\alpha$-molecules that also provide the necessary symmetry properties. We further show in Section 3.4 that local symmetry properties may be used in a related manner to yield a blob measure. We specifically derive functions that estimate the local tangent orientations of edge contours and ridge centerlines, local diameters of ridges and blobs, and the local contrast of a feature; these are summarized in Section 3.5.

To evaluate the strengths and weaknesses of the proposed measures in comparison with other state-of-the-art methods for feature detection, we performed extensive numerical experiments with respect to synthetic images with reliable ground truths, which may be found in Section 4. This section also contains a figure of merit to measure success of edge detection that is less forgiving of false positives than Pratt’s figure of merit ($\mathrm{PFOM}$) [1]. To further demonstrate the applicability of the proposed framework, we consider the application of the ridge measure in the context of retinal image analysis in Section 5.1 and show how the proposed blob measure can be applied for automatically counting the number of cell colonies in a Petri dish in Section 5.2.

We have also developed a Matlab toolbox that implements all of the proposed measures for two-dimensional digital images. The implementation strategies are briefly explained in Section LABEL:sec:digitalimplementation of the Supplemental Materials, and the toolbox can be downloaded from http://www.math.uni-bremen.de/cda/software.html.

2 Related Work

We begin by giving a short review of the most important tools and concepts that will later be used to derive the respective edge, ridge, and blob measures.

2.1 $\alpha$-Molecules

Active research in the fields of time-frequency analysis and applied harmonic analysis led to the development of a wide variety of function systems and techniques that yield efficient multiscale representations for different classes of functions and signals. Building on results and ideas from classical Fourier analysis and the work of pioneers like Alfred Haar [37], prominent examples include wavelets [25] but also systems based on anisotropic scaling such as ridgelets [17, 18], curvelets [19, 20], shearlets [51, 49], and contourlets [28], which were developed to address weaknesses of the traditional wavelet transform when dealing with anisotropic features such as edges or ridges that can appear in two- or higher-dimensional signals.

In spite of various differences in their construction and their approximation behavior for different classes of signals, many of the aforementioned function systems are constructed by applying translation, scaling, and rotation or shear operators to generator functions, which typically need to satisfy a certain admissibility condition. In [35], so-called $\alpha$-molecules were introduced to provide a unifying framework that includes different constructions of multiscale representations in order to enable the transfer of results from one theory to another, to allow for a categorization of multiscale representations with respect to the approximation behavior and also to facilitate the construction of novel systems.

Systems of $\alpha$-molecules are constructed by translating, scaling, and rotating members of a possibly infinite set of generator functions, where the latter two operations are performed by applying scaling matrices

[TABLE]

and rotation matrices

[TABLE]

to the argument of a generator function, where $\mathbb{T}=\left[-\frac{\pi}{2},\frac{\pi}{2}\right)$. Note that the parameter $\alpha$ in (1) controls the degree of anisotropy and interpolates between isotropic scaling ($\alpha=1$), which is used in the definition of wavelet-based systems, and fully anisotropic scaling ($\alpha=0$), as applied in the definition of ridgelets.

For a fixed parameter $\alpha\in[0,1]$, a single $\alpha$-molecule is thus defined by its corresponding generator and a point $(s,\theta,y)$ in the parameter space

[TABLE]

where $s$ is a scaling, $\theta$ a rotation, $y$ a translation parameter, and $\mathbb{R}_{>0}$ denotes the positive real numbers.

Let $\Lambda$ be an index set and $\alpha\in[0,1]$ be fixed. A family of functions $\left(m_{\lambda}\right)_{\lambda\in\Lambda}$ is called a system of $\alpha$-molecules of order $(L,M,N_{1},N_{2})$ with $L,M,N_{1},N_{2}\in\mathbb{N}_{0}\cup\{\infty\}$, if each molecule can be written as

[TABLE]

where $(s_{\lambda},\theta_{\lambda},y_{\lambda})\in\mathbb{P}$ is a point in the $\alpha$-molecule parameter space and the corresponding generator $g^{(\lambda)}\in L^{2}(\mathbb{R}^{2})$ satisfies

[TABLE]

for all $\left|\rho\right|\leq L$, where $\lesssim$ means that there exists a constant $C>0$ such that the left hand side is smaller than the right hand side times $C$ for all points in the considered domain. If one or more of the parameters $L$, $M$, $N_{1}$, and $N_{2}$ equal $\infty$, then (5) is interpreted to hold for arbitrarily large values of the respective parameter(s).

The authors of [35] aimed to provide a maximum amount of flexibility by allowing each molecule of an $\alpha$-molecule system to be based on a different generator in order to allow shears to be replaced by rotations. In practice, however, systems of $\alpha$-molecules are often constructed from only a small number of generators. For instance, wavelet orthonormal bases of $L^{2}(\mathbb{R}^{2})$ are typically based on four generators, while so-called cone-adapted shearlet systems (see, e.g., [49]) only require three different generators. Note that in the latter case, due to the application of shear matrices instead of rotation matrices, the elements of a cone-adapted shearlet system can not directly be expressed in the form (4). However, it can be shown that certain cone-adapted discrete shearlet systems are equivalent to an $\alpha$-molecule system that is based on an infinite number of generators, all of which satisfy Condition (5) [35].

Condition (5) requires the generators of an $\alpha$-molecule to satisfy certain time-frequency localization properties, where the parameter $L$ describes the degree of spatial localization, $M$ is the number of directional vanishing moments and $N_{1}$ and $N_{2}$ describe the smoothness of a generator.

The transform of a function $f\in L^{2}(\mathbb{R}^{2})$ with respect to an $\alpha$-molecule system $\left(m_{\lambda}\right)_{\lambda\in\Lambda}\subset L^{2}(\mathbb{R}^{2})$ can be defined via the $L^{2}$-inner product:

[TABLE]

In the remainder of this paper, we will often refer to the values $\langle f,m_{\lambda}\rangle$ as the coefficients of an $\alpha$-molecule-based transform of a signal $f$.

2.2 Edge and Ridge Detection Via Phase Congruency

Two classes of feature detection algorithms served as the inspiration for the methods in this paper: phase congruency and shearlet-based methods. For a more detailed exposition about the development of feature detection based on phase congruency, see [61, Section 4.2.2]. The key idea underlying using phase congruency to detect features is that the inner products of even- or odd-symmetric functions centered at even- or odd-symmetric features have certain traits, like having similar values, or having similarly large values. The original phase congruency measures used complex exponentials as the analyzing functions  [3, 57, 56, 70, 45]. However, such an approach attempts to define a local property based on the coefficients of the discrete Fourier transform, which – when viewed independently – only provide information about which frequency components are present in a signal but not where they occur. A possible solution was given by Kovesi, who proposed extracting the local phase congruency by using sets of differently scaled complex-valued wavelets [45]. Similar to the use of sines and cosines in the local energy approach, the real and imaginary parts of these wavelets formed Hilbert transform pairs, with the real part being even- and the imaginary part being odd-symmetric. By considering multiple scalings, Kovesi was able to combine phase information from different parts of the frequency domain and to eventually derive a dimensionless and easy-to-compute complex wavelet-based phase congruency measure for a one-dimensional function $f\in L^{2}(\mathbb{R})$ and a location $x$, namely

[TABLE]

where ${\psi^{\text{c}}}\in L^{2}(\mathbb{R})$ is a complex-valued wavelet with certain symmetry properties, $A\subset\mathbb{R}_{>0}$ a set of scaling parameters, ${\psi^{\text{c}}}_{a,x}=a^{-\frac{1}{2}}\psi\left(\frac{\cdot-x}{a}\right)$ (i.e. ${\psi^{\text{c}}}$ centered at $x$ and scaled by $a$) and $\epsilon>0$ prevents division by zero. The quantity $\operatorname{PC}_{\psi^{\text{c}}}(f,x)$ approaches $1$ if the phases of the local frequency components are in congruency and approaches [math] if they are widely spread. More details about the generalization of this measure to two dimensions and its implementation may be found here [45, 46, 47].

The measure (7) can effortlessly be computed and is by construction contrast-invariant. Features such as edges and ridges are picked up solely by the traces their structure leaves in the phases of the frequency representation of a signal rather than the local magnitude of contrast. Furthermore, as Kovesi points out in his original work, it is an interesting property of the phase congruency measure that it does not distinguish edges and ridges.

While contrast-invariance is indeed a highly desirable feature in many edge and ridge detection tasks, localizing a property of the frequency components in the time representation of a signal comes at a price. Indeed, it was already observed in [67] and [60] that replacing isotropic wavelets with anisotropically scaled analyzing functions in (7) deteriorates the detection of features in images rather than improving it. This should come as a surprise, as one would assume that anisotropically scaled elements would be better suited for detecting typically anisotropic features such as edges or ridges. Finding a way to bring together the intuitions behind the phase congruency measure $\operatorname{PC}_{\psi^{\text{c}}}$ and modern constructions of anisotropic analyzing functions was in fact one of the main factors motivating the present work.

2.3 Derivative of Gaussian Wavelets

In order to detect features by testing for local symmetry properties, we require even- and odd-symmetric real-valued functions with good localization in both the time and the frequency domain. Throughout this paper, we call a function $f$ even-symmetric, if it satisfies $f(x)=f(-x)$ for all $x$, odd-symmetric, if it holds for all $x$ that $f(x)=-f(-x)$, and symmetric if one of the above is true.

In the remainder of this work, we will consider two simple tools to construct one-dimensional even- and odd-symmetric generators with these properties, namely taking derivatives of the Gaussian as well as applying the so-called Hilbert transform. We note that such functions have already appeared in the literature in other multi-scale edge and ridge detection algorithms [53, 68, 43]. A comprehensive discussion of the harmonic analysis used in this section may be found in [12]. Note that throughout this paper, we denote the Fourier transform of a function $f\in L^{1}(\mathbb{R}^{d})$ with $\hat{f}({\xi})=\int_{\mathbb{R}}f(x)e^{-2\pi{i\mkern 1.0mu}\langle{x},{\xi}\rangle_{\ell^{2}}}\mathrm{d}{x}$ and consider the usual extension to $L^{2}(\mathbb{R}^{d})$. For $k\in\mathbb{N}_{0}$, let us denote the $k$-th derivative of the unnormalized Gaussian by

[TABLE]

where $H_{k}$ is the $k$-th Hermite polynomial (see, e.g., [69]). Given standard results from harmonic analysis, it follows that

[TABLE]

Due to $G_{0}$ being even-symmetric, the $k$-th derivative $G_{k}$ is odd-symmetric if $k$ is odd and even-symmetric if $k$ is even. Consequently, $\widehat{G_{k}}$ is purely imaginary for odd $k$ and purely real if $k$ is even. Furthermore, the Gaussian lies in the Schwartz space of rapidly decreasing functions. That is, the Gaussian is infinitely differentiable and all of its derivatives decay faster than the inverse of any polynomial. With the zero-mean property $\widehat{G_{k}}(0)=\int_{\mathbb{R}}G_{k}(x)\mathrm{d}x=0$, this implies that $G_{k}$ satisfies the wavelet admissibility condition (see, e.g., [25]) for any integer $k>0$. Note that, due to the zero-mean property, the functions $G_{k}$ are orthogonal to constant functions and the associated wavelet coefficients are hence invariant under the addition of constants.

The Hilbert transform is a bounded linear operator on $L^{2}(\mathbb{R})$ and defined for functions $f\in L^{2}(\mathbb{R})$ by

[TABLE]

where $\operatorname{sgn}$ denotes the sign function. It is easy to see from (10) that if $\hat{f}$ is purely real, $\widehat{\mathcal{H}f}$ is purely imaginary and vice versa. This implies that for any even-symmetric function $f$, its Hilbert transform $\mathcal{H}f$ is odd-symmetric, while any odd-symmetric function will become even-symmetric under the Hilbert transform. Further, the Hilbert transform leaves the magnitude of the Fourier transform invariant. In particular, this implies that the zero-mean property as well as the wavelet admissibility condition are preserved by the Hilbert transform and that for any function $f\in L^{2}(\mathbb{R})$, it holds that $\|f\|_{2}=\|\mathcal{H}f\|_{2}$. Note that the latter does not hold for the $L^{1}$-norm, where for a Schwartz function $f$, $\mathcal{H}f\in L^{1}$ if and only if it has the zero-mean property. Two examples of even- and odd-symmetric wavelets based on derivatives of the Gaussian are plotted along with their respective Hilbert transforms in the Supplementary Materials (see Figures LABEL:fig:dog1, and LABEL:fig:dog2).

2.4 Maximum Point Estimation From Discrete Samples

Each two-dimensional analyzing function within a system of $\alpha$-molecules is associated with a certain scaling and rotation parameter. These parameters can be used to obtain first estimates of the local tangent direction and width of a feature by considering the most significant analyzing function (i.e. the member of a system of functions yielding the largest coefficient). To obtain more precise measurements, we consider a refinement procedure which allows us to estimate the maximum point of a function defined in the continuum from only a few discrete samples.

Let $f\in C^{2}(\mathbb{R})$, $\{x_{n}\}_{n=1}^{N}\subset\mathbb{R}$ be a strictly increasing sequence of $N\in\mathbb{N}$ sampling points (i.e. $x_{n}<x_{n+1}$ for all $1\leq n\leq N-1$) and denote

[TABLE]

If $1<n^{*}<N$, the function $f$ has at least one local maximum point in the interval $[x_{n^{*}-1},x_{n^{*}+1}]$, which is in real-world applications often taken as a best guess for the global maximum point of $f$. One method of refining the estimate $x_{n^{*}}$ for a local maximum point is to assume that $f$ can be approximated on $[x_{n^{*}-1},x_{n^{*}+1}]$ by the unique parabola fit through the points $(x_{n^{*}-1},f(x_{n^{*}-1}))$,$(x_{n^{*}},f(x_{n^{*}}))$, and $(x_{n^{*}+1},f(x_{n^{*}+1}))$, that is

[TABLE]

The maximum point of the parabola (12), denoted as $x^{*}$, can now be computed by

[TABLE]

This approach appears in [74], where it is used to refine orientation estimation in a certain shearlet-based algorithm. A short review of this and other methods for estimating maximal points from discrete samples in the context of edge and line detection may be found in [8].

Figure 1 illustrates how (13) can be used to yield refined estimates of local maxima from a finite number of samples. In this paper, we will use (13) to improve on estimates of local tangent orientations and widths of features from finite samples of coefficients associated with differently scaled and oriented $\alpha$-molecules.

3 Symmetric Molecule-Based Feature Detection

In this section, we will estimate the likelihood of a certain feature being centered at a point $y\in\mathbb{R}^{2}$ in the image domain for a given image $f\in L^{2}(\mathbb{R}^{2})$ by considering the transform of $f$ with respect to systems of even- and odd-symmetric $\alpha$-molecules. We will furthermore demonstrate how additional information about the detected features, namely the local tangent orientations of edges and ridges as well as local widths (diameters) of ridges and blobs and heights (contrasts) can be extracted. Formally, we will define three mappings that serve as feature detectors:

[TABLE]

Furthermore, we will define functions that extract local properties of a feature centered at a given point in the image domain, namely

[TABLE]

The measures $\operatorname{E}(f,y)$, and $\operatorname{H}_{\mathrm{E}}(f,y)$ will be defined in Section 3.2; $\operatorname{R}(f,y)$, $\operatorname{W}_{\mathrm{R}}(f,y)$, and $\operatorname{H}_{\mathrm{R}}(f,y)$ in Section 3.3; $\operatorname{B}(f,y)$, $\operatorname{W}_{\mathrm{B}}(f,y)$, and $\operatorname{H}_{\mathrm{B}}(f,y)$ in Section 3.4; and $\operatorname{O}_{\mathrm{E}}(f,y)$, and $\operatorname{O}_{\mathrm{R}}(f,y)$ in Section 3.5.

Before turning to the derivation of the measures in (14) to (19), we would like to explain the basic principle of symmetric molecule-based feature detection in the one-dimensional setting, which we emphasize is distinct from phase congruency, although it is inspired by it. Figure 2a depicts what we consider an ideal one-dimensional edge, namely a step discontinuity preceded and followed by two distinct constant functions. An ideal ridge, which could be described as a short constant function interrupting a constant baseline, is shown in Figure 2b. The features plotted in Figures 2a and 2b are both centered at the point $1/2$ and the main question we aim to answer is how this point could be discriminated from its neighborhood and, in general, from other points not lying in the center of a ridge, an edge, or another kind of feature. Note that the notions of ridges and blobs coincide in the one-dimensional setting.

One way of looking at this problem is to observe that the function plotted in Figure 2a is odd-symmetric around $1/2$, while the function shown in Figure 2b is even-symmetric around $1/2$. For edges, we can additionally observe a certain kind of self-similarity. When restricting the function plotted in Figure 2a to any neighborhood around the point $1/2$ and then scaling it back to the full interval $[0,1]$, the function will remain unchanged. Ridges on the other hand are not defined by a single step-discontinuity but by two step-discontinuities which are in close proximity of each other. Any ridge can thus be associated with a width (i.e. the distance between its step-discontinuities) and is not invariant under the restricting and rescaling operation described above. In the one-dimensional setting, we will only deal with the case of edge detection but this issue will be taken into consideration when defining the two-dimensional ridge measure in Section 3.3.

To summarize: In order to detect edges, we aim to test each point in the domain of a function for the scale-invariant symmetry and self-similarity properties described in the preceding paragraph. If a point is at the location of an edge, restrictions of the analyzed function to reasonably small neighborhoods around this point should look roughly the same and, in particular, be odd-symmetric when centered at the analyzed point. A simple and efficient way of testing for these properties is to consider the coefficients associated with differently scaled even- and odd-symmetric wavelets that are centered around the point of interest.

Figures 2c and 2d show three differently scaled versions of the odd-symmetric first and the even-symmetric second derivative of the Gaussian (cf. (8) and (9)), all of which are $L^{1}$-normalized to one. Figures 2e and 2f depict the corresponding inner products with the idealized one-dimensional edge as a function of the translation parameter $t$. For the remainder of this paper, we shall call inner products of a function $f$ with an even-symmetric wavelet or $\alpha$-molecule the even-symmetric coefficients of $f$ and the inner products with an odd-symmetric wavelet or $\alpha$-molecule the odd-symmetric coefficients. Independently of scale, the even-symmetric coefficients shown in Figure 2f go to zero at the location of the step-discontinuity, as the even-symmetric wavelets are orthogonal to the locally odd-symmetric function centered at the point $1/2$. Furthermore, the odd-symmetric coefficients all peak in the point $1/2$ with a maximum value of approximately $1/2$. This is due to the fact that the jump size of the singularity is one and because the first derivative of the Gaussian is positive on $\mathbb{R}_{<0}$. Consequently, the integral of the $L^{1}$-normalized first derivative of the Gaussian over $\mathbb{R}_{<0}$ is $\frac{1}{2}$.

In summary, Figure 2 illustrates that at the location of an ideal edge, $L^{1}$-normalized symmetric wavelets show scale invariant behavior in the sense that odd-symmetric coefficients remain constant and locally maximal (with respect to translation) under scaling, while even-symmetric coefficients remain fixed at zero. In particular, this behavior is invariant of the jump size of the singularity and thus unaffected by changes in contrast.

To translate the above observations into an edge measure for functions in $L^{2}(\mathbb{R})$, we first use the derivatives of the Gaussian and the Hilbert transform defined in Section 2.3 to define sets of odd- and even-symmetric $L^{1}$-normalized wavelets:

[TABLE]

Note that $\|\mathcal{H}G_{k}\|_{1}<\infty$ for all $k\in\mathbb{N}$ is a consequence of $G_{k}$ being a Schwartz function and the zero-mean property $\widehat{G_{k}}(0)=0$. To see this, one can split the integral of the Fourier representation $\mathcal{H}G_{k}(x)=\int_{\mathbb{R}}\widehat{\mathcal{H}G_{k}}(\xi)e^{2\pi{i\mkern 1.0mu}x\xi}\mathrm{d}\xi$ at the origin and apply integration by parts to show that $\mathcal{H}G_{k}$ is dominated by $Cx^{-2}$ for some $C>0$ and hence in $L^{1}(\mathbb{R})$.

Let $\psi^{\mathrm{o}}\in\Psi^{\mathrm{o}}$, $\psi^{\mathrm{e}}\in\Psi^{\mathrm{e}}$, $a>1$ be fixed, $J=\{j_{n}\}_{n=1}^{N_{J}}\subset\mathbb{Z}$ be an increasing sequence of $N_{J}\in\mathbb{N}$ scaling parameters and denote $L^{1}$-normalized dilates and shifts of the generating wavelets $\psi^{\mathrm{o}}$ and $\psi^{\mathrm{e}}$ with

[TABLE]

for scaling parameters $j\in J$ and translation parameters $y\in\mathbb{R}$. Note that $f\mapsto a^{j}f(a^{j}\cdot)$ defines an isometric mapping of $L^{1}$. As already discussed earlier in the case of the first derivative of the Gaussian, the odd-symmetric coefficient at the location of an ideal edge is, independently of scaling, fully determined by the jump size of the discontinuity and the integral of the generating wavelet $\psi^{\mathrm{o}}$ over $\mathbb{R}_{\leq 0}$. To be precise, the value

[TABLE]

is equal to the odd-symmetric coefficients at the location of an ideal edge with jump size one. With a soft-thresholding parameter $\beta>0$ that corresponds to the minimal jump size that is required for the detection of an edge and a fixed scaling offset for the even-symmetric wavelets $j^{\mathrm{e}}\in\mathbb{R}_{>0}$, a one-dimensional edge measure for a function $f\in L^{2}(\mathbb{R})$ and $y\in\mathbb{R}$ is given by

[TABLE]

where $\epsilon>0$ prevents division by zero. The measure $\breve{\operatorname{E}}^{\text{1D}}$ is bounded above by one but to ensure a mapping to $[0,1]$, we finally set

[TABLE]

For example, if a function $g$ has an ideal edge at a point $y^{*}$ with jump size $s$, then all odd-symmetric coefficients will virtually111Note that the derivatives of the Gaussian are not compactly supported in the time domain and that an ideal edge that separates two different constant functions covering the whole real line is not in $L^{2}(\mathbb{R})$ (cf. Figure 2a). be equal to $K_{\psi^{\mathrm{o}}}$ times the jump size $s$, and all even-symmetric coefficients will be zero. That is, when disregarding $\epsilon$, the measure $\breve{\operatorname{E}}^{\text{1D}}(g,y^{*})$ reduces to

[TABLE]

illustrating that the jump size of edges that can be detected by the measure $\operatorname{E}^{\text{1D}}$ is indeed bounded below by the parameter $\beta$.

The term $\beta K_{\psi^{\mathrm{o}}}$ in (22) can also be seen as a soft-thresholding parameter that is applied to all odd-symmetric coefficients and thereby implicitly denoises the analyzed signal $f$. In practice, the parameter $\beta$ is typically chosen as a function of the expected level of noise present in a given signal or image.

Apart from the lower bound on $s$ introduced by $\beta$, the measure $\operatorname{E}^{\text{1D}}$ is by construction contrast-invariant. Due to normalization with respect to the odd-symmetric coefficient with the greatest magnitude, $\operatorname{E}^{\text{1D}}$ is independent of the jump size of the discontinuity associated with an edge. In particular, regardless of the jump size, $\operatorname{E}^{\text{1D}}$ will be close to one at the location of an ideal edge and only diminish at a point $y$ if the analyzed function has even-symmetric components around $y$, or if the odd-symmetric coefficients at $y$ have a strong variance with respect to scaling. $\operatorname{E}^{\text{1D}}$ thus provides a continuous measure of the structural similarity of the neighborhood of a given point with an ideal edge, which sets it categorically apart from most gradient-based approaches to edge detection that are measuring jump sizes.

Broadly speaking, the odd-symmetric coefficients in (22) serve as evidence in favor of – while the even-symmetric coefficients serve as evidence against – the presence of an edge at a given location. To suppress the emergence of false positives in the extended neighborhood of an actual edge, it has proven useful to consider even-symmetric wavelets that have lower frequencies than their odd-symmetric counterparts, which can be achieved by setting the scale offset parameter $j^{\mathrm{e}}$ to any real value greater than [math]. By doing so, the even-symmetric coefficients have a slower decay on the left and right side of an edge and can thereby better suppress the faster decaying odd-symmetric coefficients in the neighborhood of an edge. The effect of the scale offset parameter is illustrated in Figure 3, which visualizes the one-dimensional edge measure $\operatorname{E}^{\text{1D}}$ computed on a piecewise polynomial signal with additive Gaussian noise. The odd- and even-symmetric wavelets were chosen to be the first derivative of the Gaussian and its Hilbert transform, that is $\psi^{\mathrm{o}}=\frac{G_{1}}{\|G_{1}\|_{1}}$ and $\psi^{\mathrm{e}}=\frac{\mathcal{H}G_{1}}{\|\mathcal{H}G_{1}\|_{1}}$. All numerical results in this manuscript were in fact obtained by considering even- and odd-symmetric generators that form Hilbert transform pairs, that is, up to normalization $\psi^{\mathrm{o}}=\mathcal{H}\psi^{\mathrm{e}}$. Such pairs of wavelets also play an important role in the dual-tree complex wavelet transform [64]. One advantage of considering Hilbert transform pairs is that the respective functions only differ with respect to their symmetry properties in the sense that they differ by a $90^{\circ}$ phase shift but have the same magnitude response in the Fourier domain. The other parameters were selected as follows: $a=2$, $N_{J}=4$, $\beta=0.03$, $j^{\mathrm{e}}=0$ for the plot in Figure 3b and $j^{\mathrm{e}}=2$ for the plot in Figure 3c.

The main difference between edges and ridges is that edges are locally odd-symmetric features while ridges are locally even-symmetric features. Thus, to obtain a one-dimensional ridge measure, it almost suffices to simply switch the roles of the even- and odd-symmetric wavelets in (22). However, due to the fact that ridges, in contrast to edges, are associated with widths, the even-symmetric coefficients are not scaling invariant at the location of an ideal ridge as it is the case for odd-symmetric coefficients at the location of an ideal edge. As the one-dimensional edge measure was mainly derived to illustrate the general concept behind this work, a ridge measure will only be rigorously defined for the two-dimensional case in Section 3.3.

3.1 Symmetric $\alpha$-Molecules for Edge, Ridge, and Blob Detection in 2D

We now generalize the one-dimensional edge measure $\operatorname{E}^{\text{1D}}$ to edges, ridges, and blobs in the two-dimensional setting (cf. Figure 4). With small adjustments, the symmetry and self-similarity properties that characterize edges and ridges in the one-dimensional setting can also be utilized to extract features from two-dimensional signals. A function describing a two-dimensional image is locally odd-symmetric at a point lying on an edge with respect to the tangent of the edge contour. On the other hand, the tangents of the contour of a ridge define axes of symmetry associated with even-symmetry while a function is locally point symmetric at the center of a blob. In order to detect two-dimensional edges and ridges, we will consider generating functions in $L^{2}(\mathbb{R}^{2})\cap L^{1}(\mathbb{R}^{2})$ that are defined as tensor products of the Gaussian $G_{0}$ with one-dimensional odd- or even-symmetric wavelets from the sets $\Psi^{\mathrm{e}}$ and $\Psi^{\mathrm{o}}$. This way, we will be able to test for, in a sense, one-dimensional symmetry properties that are present in two-dimensional functions along specific directions. In the case of blob detection, we will consider tensor products of two even-symmetric wavelets. To this end, we define the sets

[TABLE]

that contain two-dimensional wavelet generators that are based on even- and odd-symmetric one-dimensional wavelets and the Gaussian $G_{0}$, while the generators in

[TABLE]

are even-symmetric and have vanishing moments in both dimensions. The parameters $c_{1}$ and $c_{2}$ can be used to adjust the scale of the even- and odd-symmetric generators as well as their aspect ratio. Examples of two-dimensional symmetric generators contained in eqs. 25, 26 and 27 are depicted in Figures LABEL:fig:genodd2d to LABEL:fig:geneventilde2d.

Analogous to the one-dimensional case, we require sets of differently scaled analyzing functions to test for scale-invariant local symmetry properties. However, the local symmetry at the location of a two-dimensional feature is typically only given with respect to a specific symmetry axis. Thus, the analysis of two-dimensional signals cannot solely be based on differently scaled wavelet-like functions but also requires a means of changing their preferred orientation. A natural way of obtaining systems of dilations, rotations, and shifts of a single two-dimensional generator is to consider the $\alpha$-molecule framework [35] briefly described in Section 2.1. Let $\alpha\in[0,1]$ be fixed, $g\in\Psi^{\mathrm{e}}_{2}\cup\Psi^{\mathrm{o}}_{2}\cup\widetilde{\Psi^{\mathrm{e}}_{2}}$ be a symmetric generator, $a>1$, $J\subset\mathbb{Z}$ be a set of scaling parameters, and $\Theta\subset\left[-\frac{\pi}{2},\frac{\pi}{2}\right)$ be a set of orientation parameters, then the corresponding system of $L^{1}$-normalized symmetric molecules (noting the $L^{1}$-isometric dilation) is given by

[TABLE]

If further for all $\left|\rho\right|\leq L$ with $L,M,N_{1},N_{2}\in\mathbb{N}_{0}\cup\{\infty\}$

[TABLE]

holds with the implicit constants independent of the choice of parameters $j\in J$, then we call the associated set $SM(g,\alpha,a,J,\Theta)$ a system of $L^{1}$-normalized symmetric $\alpha$-molecules of order $(L,M,N_{1},N_{2})$.

The parameter $\alpha$ describes the degree of anisotropy introduced to a system by the scaling matrix $A_{a^{j},\alpha}$ and should be chosen as a function of the expected level of noise and the regularity of the structures that are to be detected. If a signal was perturbed with severe noise but the contours of its edges or ridges are smooth, introducing a high degree of anisotropic scaling by choosing $\alpha$ close to zero can significantly improve the detection performance. However, setting $\alpha$ close to zero can lead to problems at locations where the contour of a feature is irregular and can thus not be well approximated by a linear function. We emphasize that the elements in (28) are $L^{1}$-normalized, while classical systems of $\alpha$-molecules are normalized with respect to the $L^{2}$-norm.

The following proposition specifies the order of symmetric molecule systems obtained from generators defined by eqs. 26, 25 and 27 in terms of a localization parameter $L$, a vanishing moments parameter $M$ and smoothness parameters $N_{1}$ and $N_{2}$ (cf. (5)).

Theorem 3.1**.**

Let $k\in\mathbb{N}$ and denote four types of generators by

[TABLE]

*Let $\alpha\in[0,1]$, $a>1$, $J\subset\mathbb{Z}$ be a set of scaling parameters, and $\Theta\subset[-\frac{\pi}{2},\frac{\pi}{2})$ a set of rotation parameters then it holds that $SM(\|g_{k}^{(i)}\|_{1}^{-1}g_{k}^{(i)},\alpha,a,J,\Theta)$ is of order $(L,\max\{0,k-L\},\infty,\infty)$, where $L\in\mathbb{N}_{0}\cup\{\infty\}$ for $i\in\{1,3\}$, and $L\in\{0,\ldots,k-1\}$ for $i\in\{2,4\}$. *

Proof.

We need to verify that (29) holds. Let $\rho\in\mathbb{N}_{0}^{2}$ and observe that

[TABLE]

where $P(\xi_{1})$ is a polynomial of degree $k+\rho_{1}$ whose growth is asymptotically bounded by $\left|\xi_{1}\right|^{k+\rho_{1}}$ and whose decay towards zero is asymptotically bounded by $\left|\xi_{1}\right|^{\max\{0,k-\rho_{1}\}}$. Let $N_{1},N_{2}\in\mathbb{N}_{0}$. We use that $e^{-x^{2}}\lesssim(1+\left|x\right|)^{-N}$ for all $N\in\mathbb{N}_{0}$ and write

[TABLE]

independently of the choice of parameters $a>1$, $j\in\mathbb{Z}$ and $\alpha\in[0,1]$. This shows (29) and a similar argument leads to the same properties for the generator $g_{k}^{(3)}$.

For $g_{k}^{(2)}$, we consider $\left|\rho\right|<k$, $N_{1},N_{2}\in\mathbb{N}_{0}$ and write

[TABLE]

where $\partial_{1}^{\rho_{1}}\left|\xi_{1}\right|(\xi_{1})^{k-1}e^{-\pi^{2}(\xi_{1})^{2}}$ only exists for $\rho_{1}<k$. By repeating the same steps as in the case of $\widehat{g_{k}^{(1)}}$, we can conclude that

[TABLE]

independently of the choice of parameters $a>1$, $j\in\mathbb{Z}$ and $\alpha\in[0,1]$. The same argument can be used to show that $SM(\|g_{k}^{(4)}\|_{1}^{-1}g_{k}^{(4)},\alpha,a,J,\Theta)$ is of order $(L,\max\{0,k-L\},\infty,\infty)$ for $L\in\{0,\ldots,k-1\}$. ∎

Theorem 3.1 illustrates that symmetric molecule systems based on derivatives of the Gaussian are well-localized and smooth in the time domain while the number of vanishing moments increases with the number of derivatives. However, Theorem 3.1 also shows that changing the symmetry properties of a generator by applying the Hilbert transform reduces smoothness in the Fourier domain and thus localization in the time domain.

3.2 Edge Detection

In the neighborhood of a point $y$ that lies on an ideal edge (cf. Figure 4a), the image is locally odd-symmetric with respect to a symmetry axis defined by the tangent of the edge contour. The measure $\operatorname{E}^{\text{1D}}$ can thus be generalized to yield a two-dimensional edge measure by restricting the analysis to symmetric molecules whose direction of vanishing moments is orthogonal to the tangent of the edge contour going through $y$. Let $g^{\mathrm{o}}\in\Psi^{\mathrm{o}}_{2}$, $g^{\mathrm{e}}\in\Psi^{\mathrm{e}}_{2}$ be a pair of odd- and even-symmetric generators, $\alpha\in[0,1]$, $a>1$, $J=\{j_{n}\}_{n=1}^{N_{J}}\subset\mathbb{Z}$ be an increasing sequence of $N_{J}\in\mathbb{N}$ scaling parameters, and $\Theta\subset\mathbb{T}$ be a set of orientation parameters. For the remainder of this section, we will denote elements from the symmetric molecule system $SM(g^{\mathrm{o}},\alpha,a,J,\Theta)$ with $m^{\mathrm{o}}_{j,\theta,y}$ and elements from the system $SM(g^{\mathrm{e}},\alpha,a,J,\Theta)$ with $m^{\mathrm{e}}_{j,\theta,y}$ (cf. (28)). For an image defined by a two-dimensional function $f\in L^{2}(\mathbb{R}^{2})$, we first denote the most significant scaling and orientation parameters at a point $y\in\mathbb{R}^{2}$ by

[TABLE]

Then, with a soft-thresholding parameter $\beta>0$ and a possibly non-integer scaling offset parameter $j^{\mathrm{e}}\in\mathbb{R}_{>0}$, a two-dimensional edge measure is given by

[TABLE]

where $\epsilon>0$ prevents division by zero and $K_{\psi^{\mathrm{o}}}$ (as defined in (21) for the one-dimensional wavelet $\psi^{\mathrm{o}}$ appearing in the definition of $g^{\mathrm{o}}$ (25)) corresponds to the odd-symmetric coefficient at the location of an ideal edge with jump-size one. As in the one-dimensional case, we finally set

[TABLE]

to ensure that the measure maps to $[0,1]$.

As we already pointed out during the derivation of the one-dimensional edge measure $\operatorname{E}^{\text{1D}}(f,y)$ earlier in this section (cf. (22)), the measure (44) is partly based on the observation that the coefficients obtained from an odd-symmetric and $L^{1}$-normalized analyzing function are invariant to scaling if it is spatially and directionally aligned with an ideal edge. In this regard, we would like to emphasize an interesting connection to other shearlet-based approaches to edge detection, which utilize the special decay behavior of $L^{2}$-normalized shearlets. It was shown in [36] that the coefficents of band-limited shearlets at a point that lies on a smooth boundary curve, and whose local tangent direction is aligned with the considered shearlet, decay with $O(s^{-3/4})$ for $s\to\infty$, where $s$ denotes the scaling parameter (cf. (4)). Similar results have also been shown in the case of compactly supported shearlet generators [50]. This is related to the measures proposed here in the sense that, when considering a two-dimensional generator $\psi\in L^{1}(\mathbb{R}^{2})\cap L^{2}(\mathbb{R}^{2})$ and parabolic scaling (i.e., $\alpha=1/2$), it holds that

[TABLE]

In other words, considering $L^{1}$- instead of $L^{2}$-normalization precisely neutralizes the decay of shearlet coefficients at points that lie on smooth boundary curves.

Further note that, since $\|m^{\mathrm{o}}\|_{L^{1}}=1$ for all $m^{\mathrm{o}}\in SM(g^{\mathrm{o}},\alpha,a,J,\Theta)$, the coefficient associated with the most significant scale parameter $j^{*}(y)$ at a point $y$ can directly be related to the contrast (i.e, the height) of an edge going through $y$. By considering $K_{\psi^{\mathrm{o}}}$ (cf. (21)), which denotes the value of an odd-symmetric coefficient at the location of an ideal edge with jump size one, we define the local height measure for edges as

[TABLE]

3.3 Ridge Detection

For ridge detection, we can exploit the fact that in the neighborhood of a point that lies on the centerline of an ideal ridge (cf. Figure 4b), the image is locally even-symmetric with respect to a symmetry axis defined by the tangent of the centerline. A two-dimensional ridge measure can in principle be obtained by simply interchanging the roles of the odd- and even-symmetric molecules $m^{\mathrm{o}}$ and $m^{\mathrm{e}}$ in the definition of the edge measure $\operatorname{E}(f,y)$. However, to fully retain contrast invariance, it is necessary to also take into account the width of the ridge, that is, the distance between the two jump singularities that define a ridge in an idealized setting. Contrary to the case of edge detection, even-symmetric coefficients of $L^{1}$-normalized molecules that are centered around the location of a ridge are not invariant to scaling but depend on the scaling parameter as well as the width of the respective ridge. This behavior is illustrated with respect to ridges of three different widths in Figure 5.

For an even-symmetric generator $g^{\mathrm{e}}\in\Psi^{\mathrm{e}}_{2}$, let us denote the coefficient at the location of an ideal ridge whose tangent direction in $y$ agrees with the orientation of $g^{\mathrm{e}}$ and that has width $2r$ and height $1$ by

[TABLE]

where $\psi^{\mathrm{e}}$ is the one-dimensional even-symmetric wavelet associated with the two-dimensional generator $g^{\mathrm{e}}$ (cf. (26)). We furthermore define the radius of $g^{\mathrm{e}}$ as half of the width of the ideal ridge which yields the most significant coefficient, that is,

[TABLE]

and extend this notion to arbitrarily shifted and rotated molecules by setting $r_{m^{\mathrm{e}}_{0,\theta,y}}=r_{m^{\mathrm{e}}_{0,0,0}}=r_{g^{\mathrm{e}}}$. The radius of a dilated molecule can directly be computed as a function of the scale parameter and the underlying generator, that is, $r_{m^{\mathrm{e}}_{j,\theta,y}}=a^{-j}r_{m^{\mathrm{e}}_{0,0,0}}=a^{-j}r_{g^{\mathrm{e}}}$. Note that in most practical cases, the radius $r_{g^{\mathrm{e}}}$ is equivalent to the distance between the origin and the first zero crossing of the associated one-dimensional wavelet $\psi^{\mathrm{e}}$ (cf. Figure 5a).

Let $g^{\mathrm{e}}\in\Psi^{\mathrm{e}}_{2}$, $g^{\mathrm{o}}\in\Psi^{\mathrm{o}}_{2}$ be a pair of odd-symmetric and even-symmetric generators, $\alpha\in[0,1]$, $a>1$, $J=\{j_{n}\}_{n=1}^{N_{J}}\subset\mathbb{Z}$ be an increasing sequence of $N_{J}\in\mathbb{N}$ scaling parameters, and $\Theta\subset\mathbb{T}$ be a set of orientation parameters. Furthermore let $f\in L^{2}(\mathbb{R}^{2})$ be a two-dimensional image. Analogous to the edge detection case, we denote the scaling and rotation parameters associated with the most significant even-symmetric molecule at a point $y\in\mathbb{R}^{2}$ in the image plane by

[TABLE]

For each point $y$, we first estimate the width of an assumed ideal ridge whose centerline passes through $y$ and whose tangent direction at $y$ agrees with the most significant orientation parameter $\theta^{*}(y)$. We use $j^{-}(y)$ and $j^{+}(y)$ to denote the scaling parameters in the strictly increasing sequence $J$ that precede, respectively succeed, the most significant scale $j^{*}(y)$. The width measure is obtained by computing a refinement of the most significant scaling parameter as the maximum point of a parabola fit through the points

[TABLE]

Following (13), we define the local width measure as

[TABLE]

where the values $c_{1}$ and $c_{2}$ are chosen according to (12). If $j^{*}(y)$ is the first or last index in $J$, we set $\operatorname{W}_{\mathrm{R}}(f,y)=2\left(a^{-j^{*}(y)}r_{g^{\mathrm{e}}}\right)$. Under the assumption that an ideal ridge with width $\operatorname{W}_{\mathrm{R}}(f,y)$, height $1$, and a local tangent direction that agrees with $\theta^{*}(y)$ passes through $y$, we can directly compute the associated even-symmetric coefficients of differently dilated molecules with rotation parameter $\theta^{*}(y)$ and denote

[TABLE]

Based on the estimated width $\operatorname{W}_{\mathrm{R}}(f,y)$, we can further define a measure that computes the height of an assumed ideal ridge whose centerline passes through $y$ as the ratio between the most significant coefficient and the value $K_{g^{\mathrm{e}}}(f,j^{*}(y),y)$, that is,

[TABLE]

A two-dimensional ridge measure can now be defined by carefully modifying the edge measure $\operatorname{E}(f,y)$. To detect locally even-symmetric features instead of odd-symmetric ones, the roles of the even-symmetric molecules $m^{\mathrm{e}}$ and odd-symmetric molecules $m^{\mathrm{o}}$ need to be interchanged. While the odd-symmetric coefficients at the location of an ideal edge are invariant to scaling, this is not the case for the even-symmetric coefficients at the location of an ideal ridge. This can be taken into account by adjusting the normalizing factor for every scaling parameter $j\in J$ using the values $\operatorname{H}_{\mathrm{R}}(f,y)K_{g^{\mathrm{e}}}(f,j,y)$, which correspond to the even-symmetric coefficients at the location of an ideal ridge with width $\operatorname{W}_{\mathrm{R}}(f,y)$ and height $\operatorname{H}_{\mathrm{R}}(f,y)$. Finally, when using strongly oscillating generators, even the sign of the even-symmetric coefficients at the location of an ideal ridge is not invariant to scaling. This can be addressed by considering the expected sign $\operatorname{sgn}(\operatorname{H}_{\mathrm{R}}(f,y)K_{g^{\mathrm{e}}}(f,j,y))$ for each scale parameter when taking the sum of the even-symmetric coefficients. To simplify our notation, we write $h(f,j,y)=\operatorname{H}_{\mathrm{R}}(f,y)K_{g^{\mathrm{e}}}(f,j,y)$ and a two-dimensional ridge measure is now given by

[TABLE]

where $\epsilon>0$ prevents division by zero, $\beta>0$ is a soft-thresholding parameter, and $j^{\mathrm{o}}\in\mathbb{R}_{>0}$ denotes the scaling offset between even- and odd-symmetric molecules. As in the case of edge detection, we finally set

[TABLE]

3.4 Blob Detection

The measures $\operatorname{E}(f,y)$ and $\operatorname{R}(f,y)$ detect features that exhibit symmetry properties that are locally defined by a singular symmetry axis, namely the normal of the tangent of an edge contour or the centerline of a ridge. A similar approach can be used to detect features in two-dimensional images that have less anisotropic symmetry properties in the sense that locally, the center of such features is a point of symmetry for more than one or even all possible directions. Let us consider an idealized blob as a filled square, which has four axes of even symmetry, or a filled circle, for which each line that passes through its center is an axis of even symmetry (cf. Figure 4c).

To detect such blobs, we will use a generator $\widetilde{g^{\mathrm{e}}}\in\widetilde{\Psi^{\mathrm{e}}_{2}}$ that is defined by the tensor product of an even-symmetric wavelet with itself (cf. (27)). Like a square, such a generator has by construction four axes of symmetry and is invariant to rotations by an angle of $\frac{\pi}{2}$. Let us denote the coefficient at the location of an ideal square-shaped blob with side length (width) $2r$ and height $1$ by

[TABLE]

and retain the definition of the radius of $\widetilde{g^{\mathrm{e}}}$ as half of the side length of the square-shaped blob which maximizes its response, that is,

[TABLE]

Let $g^{\mathrm{o}}\in\Psi^{\mathrm{o}}_{2}$ be an odd-symmetric generator, $a>1$, $J=\{j_{n}\}_{n=1}^{N_{J}}\subset\mathbb{Z}$ be an increasing sequence of $N_{J}\in\mathbb{N}$ scaling parameters, and $\Theta\subset\mathbb{T}$ be a set of orientation parameters. Due to the isotropic nature of the features we aim to detect, and contrary to the cases of edge and ridge detection, we set $\alpha=1$ and denote elements from the molecule system $SM(\widetilde{g^{\mathrm{e}}},1,a,J,\Theta)$ with $\widetilde{m^{\mathrm{e}}}_{j,\theta,y}$ and elements from the system $SM(g^{\mathrm{o}},1,a,J,\Theta)$ with $m^{\mathrm{o}}_{j,\theta,y}$. For a two-dimensional image $f\in L^{2}(\mathbb{R}^{2})$, we again denote the most significant pair of scaling and rotation parameters with respect to the even-symmetric generator $\widetilde{g^{\mathrm{e}}}$ at a point $y\in\mathbb{R}^{2}$ with

[TABLE]

As in the case of ridge detection, the even-symmetric coefficients at the center of a blob are not invariant to scaling but depend on both the scaling parameter $j$ and the width of the respective blob. For each point $y$, we again estimate the width of an assumed ideal blob by fitting a parabola through the points

[TABLE]

where $j^{-}(y)$ denotes the predecessor and $j^{+}(y)$ the successor of $j^{*}(y)$ in the sequence $J$ and compute

[TABLE]

where the values $c_{1}$ and $c_{2}$ are chosen according to (12). Under the assumption that an ideal square-shaped blob with width $\operatorname{W}_{\mathrm{B}}(f,y)$ and height $1$ is centered on $y$, we can also directly compute the associated even-symmetric coefficients of differently dilated molecules and denote

[TABLE]

The height of the blob can then be computed by

[TABLE]

We are now ready to define a measure for blob detection that can be understood as a generalization of the ridge measure $\operatorname{R}(f,y)$ in the sense that we aim to detect features that are associated with more than one axis of even symmetry. In particular, we will take into account that axes of even symmetry induce rotation invariance. A perfect square, for instance, has four symmetry axes and is invariant to rotations by $\frac{\pi}{2}$, while a pentagon has five axes of even symmetry and is invariant to rotations by $\frac{2\pi}{5}$. Let for each point $y\in\mathbb{R}^{2}$ in the image plane $\Theta^{*}(y)\subset\Theta$ denote the subset of rotation parameters that correspond to the expected rotation invariance properties of the blob centered on $y$. The subset $\Theta^{*}(y)$ thus depends on the most significant rotation parameter $\theta^{*}(y)$ and on the symmetry and rotation invariance properties of the specific shape we aim to detect. Note that in the case of circles, $\Theta^{*}(y)=\Theta$ independently of the most significant scale paramter $j^{*}(y)$ for all points $y$. With a soft-thresholding parameter $\beta>0$ and a possibly non-integer scaling offset $j^{\mathrm{o}}\in\mathbb{R}_{>0}$, the blob measure is given by

[TABLE]

where $\epsilon>0$ prevents division by zero. We define

[TABLE]

3.5 Estimation of Local Tangent Directions for Edges and Ridges

The width estimates $\operatorname{W}_{\mathrm{R}}(f,y)$ and $\operatorname{W}_{\mathrm{B}}(f,y)$ exploit the fact that even-symmetric generators from the sets $\Psi^{\mathrm{e}}_{2}$ and $\widetilde{\Psi^{\mathrm{e}}_{2}}$ are themselves associated with a width which, in combination with the scale parameter $j$, also determines the width of the elements in the associated system of symmetric molecules. By considering the orientation parameter $\theta$, a similar approach can be used to estimate local tangent directions of anisotropic features such as edges or ridges. Let us first consider the case of edges and a system $SM(g^{\mathrm{o}},\alpha,a,J,\Theta)$ of odd-symmetric molecules with $g^{\mathrm{o}}\in\Psi^{\mathrm{o}}_{2}$, $\alpha\in[0,1]$, $a>1$, an increasing sequence $J=\{j_{n}\}_{n=1}^{N_{J}}\subset\mathbb{Z}$ of $N_{J}\in\mathbb{N}$ scaling parameters, and a clockwise (or counter-clockwise) ordered sequence $\Theta=\{\theta_{n}\}_{n=1}^{N_{\Theta}}\subset\mathbb{T}$ of $N_{\Theta}\in\mathbb{N}$ orientation parameters on the torus. Each odd-symmetric molecule $m^{\mathrm{o}}_{j,\theta,y}\in SM(g^{\mathrm{o}},\alpha,a,J,\Theta)$ (cf. (28)) has a preferred orientation defined by the rotation parameter $\theta$ and the preferred orientation of the generator $g^{\mathrm{o}}$. At the location of an ideal edge, the odd-symmetric coefficient is maximized with respect to rotation if the direction of the one-dimensional wavelet $\psi^{\mathrm{o}}$ used in the construction of the generator $g^{\mathrm{o}}$ is orthogonal to the local tangent of the edge. For an image $f\in L^{2}(\mathbb{R}^{2})$ and a point $y\in\mathbb{R}^{2}$, let us again denote the most significant scaling and orientation parameters by

[TABLE]

Analogous to the width measures $\operatorname{W}_{\mathrm{R}}(f,y)$ and $\operatorname{W}_{\mathrm{B}}(f,y)$, we can use $\theta^{*}(y)$ as a first approximation of the tangent direction of a potential edge going through the point $y$. Again, this approximation can be refined by computing the maximum point of a parabola fit through the points $\left(\theta^{-}(y),\langle f,m^{\mathrm{e}}_{j^{*}(y),\theta^{-}(y),y}\rangle\right)$, $\left(\theta^{*}(y),\langle f,m^{\mathrm{e}}_{j^{*}(y),\theta^{*}(y),y}\rangle\right)$, and $\left(\theta^{+}(y),\langle f,m^{\mathrm{e}}_{j^{*}(y),\theta^{+}(y),y}\rangle\right)$, where $\theta^{-}(y)$ denotes the preceding and $\theta^{+}(y)$ the succeeding orientation parameter of $\theta^{*}(y)$ in the sequence $\Theta$. The local tangent orientation can then be estimated by applying (13), that is,

[TABLE]

where the values $c_{1}$ and $c_{2}$ are chosen as in (12). To obtain a function $\operatorname{O}_{\mathrm{R}}(f,y)$ that estimates local tangent orientation of ridges, it suffices to replace the odd-symmetric generator $g^{\mathrm{o}}$ with an even-symmetric generator $g^{\mathrm{e}}\in\Psi^{\mathrm{e}}_{2}$.

4 Performance Evaluation on Synthetic Images

To evaluate the performance of the proposed measures, we first consider sets of synthetic digital grayscale images for which the locations and properties of features are unambiguously defined. For each type of feature, we created two test images of varying complexity that contain exemplary instances of the respective feature. To test the stability of the proposed measures in adverse conditions, we furthermore applied a combination of Gaussian and Poisson (shot) noise to each test image with three different degrees of severity (no noise, medium noise, severe noise). The Supplemental Materials contain visualizations of all six synthetic images, their distorted versions as well as the corresponding ground truths in Section LABEL:sec:evalsyn_supp and an explanation of the Matlab variables used in Section LABEL:sec:digitalimplementation.

In all cases, both the synthetic image and the ground truth are based on a geometry that was defined semi-randomly on the unit square and then discretized on a grid of $768\times 768$ pixels. Hence, the ground truth is not directly tied to the discrete digital image and therefore less likely to be biased towards any of the different algorithms we consider in our numerical experiments. For all three features, the feature location ground truth is stored as a binary image where white pixels denote the contours of edges, centerlines of ridges, and centers of blobs, respectively. In the case of edges and ridges, the ground truth also contains the orientation of the tangent for each pixel that lies on the contour of an edge or the centerline of a ridge. The same holds for local feature widths in the case of ridges and blobs. The height measures $\operatorname{H}_{\mathrm{R}}(f,y)$, and $\operatorname{H}_{\mathrm{B}}(f,y)$ are essentially defined as functions of the width measures $\operatorname{W}_{\mathrm{R}}(f,y)$, and $\operatorname{W}_{\mathrm{B}}(f,y)$ (cf. (54), respectively (63)). The validity of a height measure therefore almost exclusively depends on the precision of the associated width measure. For brevity, the measures $\operatorname{H}_{\mathrm{R}}(f,y)$, and $\operatorname{H}_{\mathrm{B}}(f,y)$ were thus omitted from the numerical evaluation.

We consider different metrics to evaluate the results obtained by a feature detection algorithm with respect to the corresponding ground truth. For a given image, let $N_{\mathrm{gt}}$ denote the number of pixels that are marked as feature points in the ground truth and $N_{\mathrm{det}}$ the number of points that were detected as feature points by the algorithm under evaluation. We denote the set of feature points in the ground truth with $P_{\mathrm{gt}}=\{p_{n}\}_{n=1}^{N_{\mathrm{gt}}}\subset\mathbb{N}^{2}$ and the set of detected feature points with $P_{\mathrm{det}}=\{q_{n}\}_{n=1}^{N_{\mathrm{det}}}\subset\mathbb{N}^{2}$. Furthermore, let $D_{\mathrm{det}}=\{\min_{p\in P_{\mathrm{gt}}}\|q-p\|\colon q\in P_{\mathrm{det}}\}$ denote the set of all Euclidean distances from pixel coordinates in $P_{\mathrm{det}}$ to the set $P_{\mathrm{gt}}$, and $D_{\mathrm{gt}}=\{\min_{q\in P_{\mathrm{det}}}\|p-q\|\colon p\in P_{\mathrm{gt}}\}$ denote the set of all distances from points in $P_{\mathrm{gt}}$ to the set $P_{\mathrm{det}}$. A figure of merit that evaluates how successful an algorithm was at correctly detecting the location of features can then be defined by

[TABLE]

where $d_{n}\in D=D_{\mathrm{gt}}\cup D_{\mathrm{det}}$ denotes the $n$-th largest element in the union of both sets of distances, and $\gamma>0$. All $\mathrm{FOM}$-values reported in this section were computed with $\gamma=\frac{1}{9}$. Note that for $D=D_{\mathrm{det}}$ and accordingly adjusted summation bounds (i.e., replacing $N_{\mathrm{gt}}$ with $N_{\mathrm{det}}$), (68) yields exactly the definition of Pratt’s figure of merit ($\mathrm{PFOM}$) [1]. The main reason for using a slightly adjusted version here is that the original definition of $\mathrm{PFOM}$ tends to be too forgiving of false positives in the case that $N_{\mathrm{det}}<N_{\mathrm{gt}}$.

For a fixed radius $r>0$, we define the set of successfully detected points in the ground truth (true positives) as $P_{\mathrm{tp}}=\{p\in P_{\mathrm{gt}}\colon\min_{q\in P_{\mathrm{det}}}\|p-q\|\leq r\}$ and let $N_{\mathrm{tp}}\in\mathbb{N}$ denote the number of points in $P_{\mathrm{tp}}$. To evaluate width and tangent orientation measures, we compute the mean of the set of absolute errors obtained by comparing the corresponding ground truth values at points in $P_{\mathrm{tp}}$ with the values yielded by an algorithm at the nearest neighbors of the ground truth points in the set $P_{\mathrm{det}}$. Formally, we define the mean absolute error ($\mathrm{MAE}$) of an arbitrary measure $\operatorname{M}(f,{y})$ computed on a digital image $f$ as

[TABLE]

where $p_{n}\in P_{\mathrm{tp}}$ and $\operatorname{M}_{\mathrm{gt}}(f,p_{n})$ denotes the ground truth of the estimated property at the pixel $p_{n}$. Note that in the case of tangent orientation measures, the error is computed on the torus. Finally, we define the number of true positives ($\mathrm{TP}$) and the ratio of successfully detected feature points ($\mathrm{TPR}$) as

[TABLE]

respectively.

Throughout this section, the results of SymFD (the proposed algorithm) are not only compared to the ground truth but also to the results of other well-known feature detection algorithms. An overview of all considered algorithms is given in Table LABEL:tab:synthalgorithms in the Supplementary Materials, which lists for each method the types of features it can detect, the number of parameters of the algorithm and the average time it took to process one synthetic digital image. The average processing time was measured on an Intel Core i7-4790 CPU clocked at 3.60 GHz.

The parameters of each method were separately optimized for each test image and each noise configuration. In the cases of edge and ridge detection, the optimization was carried out in terms of the $\mathrm{FOM}$-value. For blob detection, the goal was to maximize the number of correctly detected blobs while keeping the number of false positives close to zero.

Note that SymFD also implements two post-processing steps, namely hard thresholding and morphological thinning to obtain thinned binary feature maps. These steps, along with the overall implementation details of SymFD, are explained in detail in Section LABEL:sec:digitalimplementation in the Supplementary Materials.

All Matlab scripts and configurations that were used to obtain the numerical results presented in this section can be downloaded from http://www.math.uni-bremen.de/cda/software.html.

4.1 Edge Detection

The synthetic images for the evaluation the edge measure $\operatorname{E}(f,y)$ and the tangent orientation measure $\operatorname{O}_{\mathrm{E}}(f,y)$ are based on shapes that are defined by closed smooth spline curves and polygons. The main advantage of considering splines and polygons is that in both cases, with the exception of corner points in a polygon, the tangent can easily be computed for every point on the respective curve. The binary images that represent the feature location ground truth were obtained by marking each pixel which intersects with any of the defined curves. Each polygon and each spline curve is constructed from a finite sequence of points. These point sequences were predefined in the case of synthetic image 1 and semi-randomly selected in synthetic image 2. The shapes in both images were filled using a range of constant values and different types of gradients to ensure a certain amount of contrast variation. Furthermore, both images were overlaid with a smooth surface to test the capability of an algorithm to properly distinguish between sharp and smooth transitions. The synthetic images, their distorted versions, and the corresponding ground truths are compiled in Figure LABEL:fig:synthimages_edges in the Supplementary Materials.

We compare the performance of SymFD with results obtained from the Canny edge detector [21], the Sobel edge filter [65], the phase congruency measure [45], and a shearlet-based edge detector [74], which can be seen as a shearlet-based generalization of Canny. In the case of the Canny and Sobel edge detectors, we applied the implementation provided by Matlab’s edge function. The phase congruency measure was computed with the phasecong3 function [47] and an implementation of the shearlet-based edge detector was kindly provided by the authors of [74]. A visual comparison of the results yielded by the different methods when processing synthetic image 2 with severe distortions may be found in Figure 6. The SymFD edge measure $\operatorname{E}(f,y)$ and the tangent orientation measure $\operatorname{O}_{\mathrm{E}}(f,y)$ after post-processing are plotted in Figure LABEL:fig:synthimages_edges_SFD in the Supplementary Materials in the case of synthetic image 1 distorted with medium noise.

The detection performance of the considered algorithms for both synthetic images and three different levels of noise are summarized in Table 1. The edge detection performance was measured in terms of the $\mathrm{FOM}$. Estimates of the local tangent directions are compared to the ground truth via the $\mathrm{MAE}_{\operatorname{O}_{\mathrm{E}}}$. Note that the $\mathrm{TPR}$ is equal to the percentage of ground truth feature points that were used in the computation of the $\mathrm{MAE}_{\operatorname{O}_{\mathrm{E}}}$ (cf. (69) and (70)). The radius defining the set of true positives $P_{\mathrm{tp}}$ was set to three pixels. Note that out of the algorithms included in the evaluation, only SymFD, the phase congruency measure, and the shearlet-based edge detector provide estimates of tangent orientations.

4.2 Ridge Detection

Similar to the case of edge detection, we generated synthetic digital images that contain exemplary ridges whose centerlines are based on smooth spline curves. The sequences of points defining the spline curves were again predefined in synthetic image 3 and chosen in a semi-random fashion in synthetic image 4. To obtain ridges, the spline curves were thickened by a value ranging from three to ten pixels. In the case of synthetic image 3, the ridge-shapes were filled using different constant values and different types of gradients, while in synthetic image 4, all ridge-shapes are filled with the same value. To test the capability of an algorithm to correctly differentiate between edges and ridges, both synthetic images also contain filled shapes that are based on closed curves. Furthermore, both images were overlaid with a smooth surface. Visualizations of the noise-free and distorted versions of synthetic image 3 and synthetic image 4 as well as the corresponding ground truths may be found in Figure LABEL:fig:synthimages_ridges in the Supplementary Materials.

Detection results for a noisy version of synthetic image 3 yielded by the ridge measure $\operatorname{R}(f,y)$, the tangent orientation measure $\operatorname{O}_{\mathrm{R}}(f,y)$, and the width measure $\operatorname{W}_{\mathrm{R}}(f,y)$ are depicted in Figure 7. The ridge detection performance of SymFD is furthermore compared with results obtained from the phase congruency measure [45], the vessel enhancement filter proposed by Frangi et al. [32], and a widely used scale space-based ridge detection method proposed by Steger [66]. For the phase congruency measure, we applied the same implementation already used during the evaluation of different edge detection methods. An implementation of the Frangi vesselness measure may be found on the Matlab Central File Exchange [48]. To obtain binary detection results, the feature maps computed by the Frangi vesselness measure were further processed with the same thresholding and thinning methods that are implemented in SymFD. For the Steger ridge detection method, we used an implementation provided by a freely available ImageJ plugin [71].

The ridge detection performances of SymFD, the phase congruency measure, the Frangi vessel enhancement filter, and the method proposed by Steger are summarized in Table 2. The similarity between the ridge location ground truth and the map of feature points detected by a single algorithm is again measured in terms of the $\mathrm{FOM}$. Estimates of the local tangent directions and ridge widths are compared to the ground truth via the $\mathrm{MAE}_{\operatorname{O}_{\mathrm{R}}}$ and the $\mathrm{MAE}_{\operatorname{W}_{\mathrm{R}}}$, respectively. The radius defining the set of true positives $P_{\mathrm{tp}}$ was set to three pixels. Similar to the evaluation of edge detection methods, errors in width and tangent direction measurements are only computed for feature points in the set $P_{\mathrm{tp}}$. Out of the algorithms included in the evaluation, SymFD, the phase congruency measure and the implementation applied to obtain the Frangi vesselness measure also provide estimates of local tangent directions. Furthermore, only SymFD and the Steger ridge detector yield estimates of local ridge widths. A visual comparison of the results yielded by the different methods when processing synthetic image 4 with severe distortions may be found in Figure LABEL:fig:synthimages_ridges_comparison in the Supplementary Materials.

4.3 Blobs

To evaluate the proposed blob measure $\operatorname{B}(f,y)$ and width estimates $\operatorname{W}_{\mathrm{B}}(f,y)$, we consider two synthetic images that show aggregations of filled circles with varying diameters and contrasts. Synthetic image 5 consists of 31 large circles with diameters ranging from 30 to 50 pixels while synthetic image 6 contains a total of 200 small blobs whose diameters are between 7 and 13 pixels. In both cases, the center points of the circles were semi-randomly selected in the sense that two center points have at least a distance of 20 pixels in synthetic image 6 and 100 pixels in synthetic image 5. Noisy versions of synthetic image 5 and synthetic image 6 and the corresponding ground truths are visualized in Figure LABEL:fig:synthimages_edges in the Supplementary Materials.

The circle detection performance of SymFD is compared to the circular Hough transform [7, 27], which can be computed in Matlab via the function imfindcircles. To compare detection results of an algorithm with the ground truth, we consider the number $\mathrm{TP}$ of true positives with respect to a radius of six pixels, as defined in (70), and the number $\mathrm{FP}$ of false positives. Here, $\mathrm{FP}$ is equal to the number of detected blob center points which are not in the vicinity of an actual center of a blob. That is, $\mathrm{FP}$ is equal to the number of elements in the set $P_{\mathrm{fp}}=\{p_{1}\in P_{\mathrm{det}}\colon\min_{p_{2}\in P_{\mathrm{gt}}}\|p_{1}-p_{2}\|>6\}$. We further consider the mean absolute errors in terms of center point location and blob width denoted by $\mathrm{MAE}_{\mathrm{ctr}}$ and $\mathrm{MAE}_{\operatorname{W}_{\mathrm{B}}}$, respectively. The results of the numerical evaluation of the blob detection performance of SymFD and the circular Hough transform are summarized in Table 3.

A visual comparison of the results yielded by SymFD and the circular Hough transform in the cases of synthetic image 5 and synthetic image 6 distorted by severe Gaussian and Poisson noise is presented in Figure 8, and SymFD blob detection results with width error map for synthetic image 6 distorted by medium noise are shown in Figure LABEL:fig:synthimages_blobs_SFD in the Supplementary Materials.

5 Applications

Previous versions of $\operatorname{E}(f,y)$ and $\operatorname{R}(f,y)$ that are based on complex-valued shearlets [60, 44, 62] were already successfully applied in different feature extraction tasks such as the detection and characterization of flame fronts [62, 44], the detection of borders of tidal flats in the Wadden Sea [44], the extraction of fracture-traces in rock masses [13], the detection of channel boundaries in seismic data [41, 42], and the automated detection of the boundaries of touching cells in scanning electron (SEM) images [55].

To further evaluate the applicability of the proposed measures, we consider two tasks which require the correct localization and characterization of certain morphological components in biomedical images. In Section 5.1, we demonstrate how the width measure $\operatorname{W}_{\mathrm{R}}(f,y)$ and estimates of local tangent directions yielded by $\operatorname{O}_{\mathrm{R}}(f,y)$ can be used in combination with the ridge measure $\operatorname{R}(f,y)$ to obtain characterizations of the geometry of retinal blood vessels. In particular, the blood vessel profiles yielded by SymFD are validated and compared to other state-of-the-art algorithms on the publicly available REVIEW database [5]. In Section 5.2, we investigate the applicability of the blob measure $\operatorname{B}(f,y)$ for automatically counting the number of grown cell colonies in a Petri dish.

5.1 Determining Widths and Orientations of Retinal Blood Vessels

Retinal imaging is one of only a few existing techniques which provides non-invasive observations of the human circulatory system. Variations in the geometry of the retinal vasculature can furthermore be linked to a wide range of ocular diseases and systemic diseases with a vascular component such as diabetes [23, 11], Alzheimer’s disease [72], or high blood pressure [39]. In recent years, many methods have been proposed for automatically extracting and analyzing the geometry of blood vessels in digital retinal images [2, 58, 29].

To give an example of how SymFD can be applied in the context of retinal image analysis, we consider the REVIEW retinal vessel reference dataset [5]. For each image in the REVIEW database, three independent experts were asked to manually define vessel edge profiles on a number of preselected segments. Each profile consists of a pair of opposite edge points that indicate the boundaries of a blood vessel. A single profile thus not only contains information about the location of a retinal blood vessel but provides a complete characterization in terms of the local width and the local normal direction of the centerline. Two examples of images in the REVIEW database are shown along with manually marked vessel profiles in Figure LABEL:fig:review_gt in the Supplemental Materials.

The REVIEW database consists of four subsets of digital images that aim to provide meaningful benchmarks for different types of challenges in retinal image processing. The high resolution image set (HRIS) contains four images that show severe cases of diabetic retinopathy. Each image in the HRIS data set has a resolution of $896\times 609\times 3$ pixels but has previously been sub-sampled by a factor of four such that, disregarding human error, the vessel widths are known to a precision of a quarter of a pixel. The vascular disease image set (VDIS) contains eight retinal images of size $1360\times 1024\times 3$ that were randomly selected from a database of images of patients attending a retinopathy clinic. The central light reflex image set (CLRIS) consists of only two images of a resolution of $2160\times 1440\times 3$ pixels that show early atherosclerotic changes in retinal vessels. Finally, the kick point image set (KPIS) contains two sub-sampled images where the normal direction of the centerline of a given location was computationally determined from the neighboring pixels on the centerline and the human observers where only asked to identify the kick-points determining the width. In total, the REVIEW database contains 5066 manually marked vessel profiles.

The diameters of blood vessels that are visible in the images of the database roughly range from $2$ to $20$ pixels. To obtain a precise detection and characterization of both thin and wide blood vessels, we use not a single pair but two pairs of systems of even- and odd-symmetric molecules that cover different but overlapping regions in the frequency domain. Both pairs of systems are based on the even-symmetric one-dimensional wavelet $\psi^{\mathrm{e}}=\frac{\mathcal{H}G_{1}}{\|\mathcal{H}G_{1}\|_{L^{1}}}$ and the odd-symmetric wavelet $\psi^{\mathrm{o}}=\frac{G_{1}}{\|G_{1}\|_{L^{1}}}$ and cover two octaves with $\text{{scalesPerOctave}}=4$, $\alpha=1$, and $\text{{nOrientations}}=16$. The first pair of systems aims to detect of thin vessel and is defined by the parameters $\text{{maxFeatureWidth}}_{1}=8$, $\text{{maxFeatureLength}}_{1}=24$, and $\text{{minFeatureWidth}}_{1}=2$. The second pair is capable of detecting vessels with a diameter of up to 24 pixels and is defined by the parameters $\text{{maxFeatureWidth}}_{2}=24$, $\text{{maxFeatureLength}}_{2}=72$, and $\text{{minFeatureWidth}}_{2}=6$. Each image in the database is processed by computing the ridge measure $\operatorname{R}(f,y)$, the width measure $\operatorname{W}_{\mathrm{R}}(f,y)$ and the tangent direction measure $\operatorname{O}_{\mathrm{R}}(f,y)$ for both pairs of systems of even- and odd-symmetric molecules with parameters $j^{\mathrm{o}}=0$, and $\beta=2$. In each case, the computations are carried out solely on the green channel of the processed RGB image, in which the blood vessels typically have the most significant contrast. For every pixel location $y$, we keep the values yielded by the pair of systems associated with the larger ridge measure $\operatorname{R}(f,y)$. We further restrict the analysis to ridges with negative contrast, that is, to locations where $\operatorname{H}_{\mathrm{R}}(f,y)\leq 0$. As the width measure $\operatorname{W}_{\mathrm{R}}(f,y)$ is not only based on the even-symmetric coefficient on the most significant scale but also on the coefficients with respect to the preceding and succeeding scaling parameters, measurements where the first or last scales are most significant are not considered. The obtained feature maps in the case of image 7 from the VDIS dataset are plotted in Figures 9b to 9d.

To validate the obtained measurements, the width and tangent direction maps yielded by SymFD are compared to the widths and orientations of the manually marked profiles. For each vessel edge profile in the ground truth, we consider the difference with respect to the width measure $\operatorname{W}_{\mathrm{R}}(f,y)$ and the orientation measure $\operatorname{O}_{\mathrm{R}}(f,y)$ at the nearest location $y$ in a two-pixel radius for which the ridge measure $\operatorname{R}(f,y)$ is greater than zero. If no such pixel exists, the corresponding vessel edge profile counts as undetected. Visualizations of the tangent direction error and the width error for measurements obtained from SymFD are visualized in the case of image 7 from the VDIS dataset in Figures 9e and 9f. For the tangent orientation estimates $\operatorname{O}_{\mathrm{R}}(f,y)$, we compute the mean absolute error $\mathrm{MAE}_{\operatorname{O}_{\mathrm{R}}}$ on the torus with respect to all manually marked profiles for each of the four subsets of images. The respective values along with the success rate $\mathrm{SR}$, which denotes the percentage of successfully matched vessel profiles, are reported in Table 4.

For the validation of the width measurements $\operatorname{W}_{\mathrm{R}}(f,y)$, we follow the procedure proposed in [4, 9, 73] and consider the standard deviation $\mathrm{SD}_{\operatorname{W}_{\mathrm{R}}}$ of differences between the measure $\operatorname{W}_{\mathrm{R}}(f,y)$ and the ground truth vessel diameters instead of the mean absolute error. This way, the comparison is independent of possibly different implicit definitions of vessel widths. We further compare SymFD to a number of classical and state-of-the-art methods for the extraction of vessel diameters from digital retinal images. In [75], vessel widths are obtained by fitting a one-dimensional Gaussian (1DG) model to the vessel intensity profiles. A similar approach using a two-dimensional Gaussian (2DG) model was proposed in [52]. The half height full width (HHFW) model defines the width of a vessel as the distance between the points at which the vessel intensity profile reaches half of the maximum intensity left and right of an assumed center point [15]. Gregson et al. proposed to obtain vessel widths by fitting a rectangle to the intensity profile such that the area under the rectangle is equal to the area under the intensity profile [34]. A graph-based edge segmentation technique for measuring vessel widths was proposed in [73], while the extraction of segment profiles (ESP) algorithm proposed in [4] is based on an active contour model. The automated retinal image analyzer (ARIA) employs a wavelet-based approach for vessel segmentation which is followed by a refinement of edge locations [9]. Finally, edge tracking based on orientation scores (ETOS) is an algorithm that yields state-of-the-art vessel width measurements by simultaneously tracking both edges of a vessel in the coupled space of positions and orientations [10]. The accuracy of vessel width measurements on images from the REVIEW database obtained by SymFD and all aforementioned methods in terms of the standard deviation $\mathrm{SD}_{\operatorname{W}_{\mathrm{R}}}$ and success rate $\mathrm{SR}$ are compiled in Table 5.

5.2 Automated Cell Colony Counting

Counting the number of cell colonies growing on an agar disk is a common but time-consuming task in many areas of biomedical research. To evaluate the applicability of the blob measure $\operatorname{B}(f,y)$ for automatically obtaining the number of colonies from digital images, we consider a set of pictures which was originally developed for the evaluation of the OpenCFU software package [33]. The set consists of digital images of 19 different plates containing between 10 and 1000 Staphylococcus aureus colonies and can be downloaded from http://opencfu.sourceforge.net/. Each plate was photographed once with a high-definition camera and once with a low-cost webcam. To provide a reference for the evaluation of automated counting methods, the number of colonies on each plate was independently counted by seven trained humans. For a single plate, the ground truth was then defined as the median of the numbers reported by the human experts. Examples of the considered images are displayed in Figure LABEL:fig:cellcol_input in the Supplemental Materials.

To automatically count the number of cell colonies in the high-definition images, all 19 images were processed with SymFD by computing the blob measure $\operatorname{B}(f,y)$ on the blue color channels using the parameters $\psi^{\mathrm{e}}=\frac{\mathcal{H}G_{1}}{\|\mathcal{H}G_{1}\|_{L^{1}}}$, $j^{\mathrm{o}}=-1$, $\text{{maxFeatureWidth}}_{\text{hd}}=20$, $\text{{maxFeatureLength}}_{\text{hd}}=20$, $\text{{minFeatureWidth}}_{\text{hd}}=10$, $\alpha=1$, $\text{{scalesPerOctave}}=3$, $\text{{nOrientations}}=16$, and $\beta_{\text{hd}}=15$. The obtained feature maps were then transformed into binary images by applying the threshold $0.03$. Centers of cell colonies were eventually obtained by computing the centroids of all connected components in the binary image via Matlab’s regionprops function. For processing the webcam images, the parameters specifying the effective support of the applied symmetric molecules had to be adapted to accommodate for the lower resolution and were set to $\text{{maxFeatureWidth}}_{\text{wc}}=8$, $\text{{maxFeatureLength}}_{\text{wc}}=8$, and $\text{{minFeatureWidth}}_{\text{wc}}=4$. Furthermore, mostly due to horizontal stripe artifacts the contrasts of the cell colonies are slightly lower in the images captured by a webcam than in their high-definition counterparts. Thus, for processing the webcam images, the parameter specifying the minimal contrast was chosen as $\beta_{\text{wc}}=6$. Examples of cell colony detection results obtained by SymFD for one high-definition image and one webcam image are visualized in Figure 10.

To validate the accuracy of the cell colony counts obtained by SymFD, we follow the procedure proposed in [33] and consider the median relative deviation from the counts provided by the human experts across all 19 images. We further compare the accuracy of SymFD with three other methods for the automated detection and counting of cell colonies. The NIST’s Integrated Colony Enumerator (NICE) [24] is a freely available software package based on a combination of extended-minima transforms and thresholding operations. We also consider an ImageJ plugin developed by Cai et al. [16] that is based on adaptive thresholding, the application of a watershed algorithm, and particle filtering. The results obtained by SymFD are compared to OpenCFU which is a highly versatile and efficiently implemented open source library for the detection of cell colonies [33]. The median relative deviations from human counts for SymFD and the aforementioned methods with respect to high-definition and webcam images of all 19 plates are compiled in Table 6. The cell colony counts for OpenCFU, NICE, and the ImageJ plugin proposed by Cai et al. were kindly provided by the author of [33]. In particular NICE and the ImageJ macro are sensitive to artifacts in the background and at the border of the Petri dish. To facilitate a fair comparison, Table 6 also reports the median relative errors when using an additional foreground mask. To illustrate the robustness of the considered methods with respect to images of agar plates with both small and large numbers of cell colonies, semilog graphs plotting for each algorithm the relative deviation for a single plate against the number of cell colonies in the corresponding plate are shown Figure 11.

6 Discussion of Numerical Results

The results of the numerical evaluation on synthetic images summarized in Tables 1, 2, and 3 indicate that the detection accuracy of the proposed measures $\operatorname{E}(f,y)$, $\operatorname{R}(f,y)$, and $\operatorname{B}(f,y)$ at least matches the performance of some of the most popular state-of-the-art algorithms for edge, ridge, and blob detection. In particular, all three measures show a high robustness in the presence of noise. Maybe the most significant instance of this may be found in the case of blob detection, where the circular Hough transform clearly fails to reliably detect filled circles in severely distorted images, while SymFD still yields a perfect detection result in the case of synthetic image 5 and correctly identifies three quarters of the blobs in synthetic image 6 (see Table 3 and Figure 8). However, the experiments also show that regarding the mere detection of features, SymFD does not provide a significant advantage over some of the already established methods for the considered test inputs. In the case of edge detection, the well-known Canny edge detector [21] matches or even slightly outperforms the detection rate of SymFD for all of the six test images (see Table 1) while the same is true in the case of ridge detection for the Steger algorithm [66] (see Table 2).

One of the most significant benefits of the proposed framework is that it not only allows for the detection of different types of features but also yields a comprehensive characterization of their geometry in terms of local tangent directions and diameters. The respective mean absolute errors reported in Tables 1 and 2 indicate a comparatively high accuracy of the tangent orientation measures $\operatorname{O}_{\mathrm{E}}(f,y)$ and $\operatorname{O}_{\mathrm{R}}(f,y)$, and the width measure $\operatorname{W}_{\mathrm{R}}(f,y)$. However, the experiments also show that for all considered test images, the width and orientation measures are sensitive to noise, where the most significant decrease in accuracy may be found in the case of $\operatorname{O}_{\mathrm{E}}(f,y)$ and synthetic image 2 (see Table 1). The visualization of the local error of the tangent direction estimates shown in Figure LABEL:fig:synthimages_edges_SFD_orierror in the Supplemental Materials furthermore reveals that the accuracy of the measure $\operatorname{O}_{\mathrm{E}}(f,y)$ deteriorates in the proximity of points on which the edge curve is not smooth (e.g. corner points). In the case of blob detection, the errors of the width measure $\operatorname{W}_{\mathrm{B}}(f,y)$ summarized in Table 3 and visualized in Figure LABEL:fig:synthimages_blobs_SFD_error_width in the Supplemental Materials are relatively large, even in comparison with the results yielded by the circular Hough transform. This can in part be explained by the fact that the two-dimensional generators from the set $\widetilde{\Psi^{\mathrm{e}}_{2}}$ are based on tensor products of one-dimensional even-symmetric wavelets and thus designed to detect square-shaped rather than circular blobs, which is also reflected in the definition of the function $\widetilde{K}_{\widetilde{g^{\mathrm{e}}}}(r)$ in Section 3.4 (see (57)). We suspect that the accuracy of the measure $\operatorname{W}_{\mathrm{B}}(f,y)$ with respect to filled circles would significantly be improved by using circular two-dimensional generators and adjusting the definition of $\widetilde{K}_{\widetilde{g^{\mathrm{e}}}}(r)$ accordingly. A first choice for such a generator could be the two-dimensional Mexican hat wavelet. We further suspect that using a generator that actually matches the symmetry properties and shape of the objects that are to be detected would also have a positive effect on the overall accuracy of the measure $\operatorname{B}(f,y)$.

The numerical experiments conducted in Section 5.1 show that the measures $\operatorname{R}(f,y)$, $\operatorname{O}_{\mathrm{R}}(f,y)$, and $\operatorname{W}_{\mathrm{R}}(f,y)$ can be used to reliably describe the geometry of blood vessels in different types of retinal images. In particular, the tangent direction estimates yielded by $\operatorname{O}_{\mathrm{R}}(f,y)$ are highly accurate in the sense that their mean deviation from the ground truth is only slightly above the variation within the group of human experts (see Table 4). While the width measure $\operatorname{W}_{\mathrm{R}}(f,y)$ yields more accurate estimates of local vessel diameters than many of the classical approaches it is also clearly outperformed by some of the more recent algorithms that were specifically developed for processing retinal images, such as ARIA [9] and ETOS [10] (see Table 5).

The applicability of the blob measure $\operatorname{B}(f,y)$ for automatically counting the number of cell colonies on an agar plate was investigated in Section 5.2. It has been suggested that in practice, an average relative deviation of less than $20\,\%$ will often be negligible in comparison to noise introduced by other experimental factors [33]. When using an additional foreground mask, SymFD achieves a median relative deviation from the ground truth of about $2\,\%$ on a set of 19 high-definition images but also on a set of 19 low-quality webcam images (cf. Table 6). This accuracy is comparable to the accuracy of the widely used OpenCFU software package [33] and suggests a high reliability of the results obtained by SymFD when applied for cell colony detection. The values reported in Table 6 also show only a slight decrease in accuracy for SymFD when omitting the foreground mask, which indicates that the measure $\operatorname{B}(f,y)$ is not sensitive to structures in the image background or at the border of a Petri dish. Similar to other methods, SymFD tends to slightly overestimate the number of cell colonies on sparse plates while underestimating the cell colony count for high-density plates (cf. Figure 11). This is mostly due to the high significance of false positives when only a few colonies are present and the difficulty of correctly identifying single colonies in dense clusters which often occur on high-density plates.

The numerical experiments described in Sections 4 and 5 suggest that the proposed measures indeed define a highly versatile and powerful framework that can be reliably applied in a wide range of different feature detection and characterization tasks. There is, however, a trade-off for this high degree of flexibility, namely the relatively large number of parameters that are required to be selected when configuring SymFD for a specific task. Furthermore, evaluating the proposed edge, ridge, and blob measures on a digital image often requires the computation of more than 100 convolutions with digital symmetric molecule filters. While an average execution time of slightly less than 20 seconds is still acceptable in most tasks, SymFD is currently far from being applicable in real-time applications.

7 Conclusion

In the present work, we propose a novel framework for the detection and characterization of features such as edges, ridges, and blobs in two-dimensional images. Drawing inspiration from investigations of the phase congruency property [57, 56, 45], the developed measures exploit the special symmetry properties of directionally sensitive analyzing functions that are based on tensor products of $L^{1}$-normalized derivatives of the one-dimensional Gaussian and their Hilbert transforms. Differently scaled and oriented two-dimensional analyzing functions are then constructed within the recently introduced $\alpha$-molecule framework [35], which provides a maximum degree of flexibility. Combining both old and new ideas from applied harmonic analysis and computational vision, the developed measures are by construction contrast-invariant and furthermore capable of fully characterizing the geometry of features in terms of local tangent directions and local widths.

By considering pairs of even- and odd-symmetric analyzing functions that are sensitive to scaling and orientation, the measures $\operatorname{E}(f,y)$, $\operatorname{R}(f,y)$, and $\operatorname{B}(f,y)$ are strongly inspired by functional properties that are known to be exhibited by neurons in the early human visual system [26]. Another property that has widely been observed in neural populations is that the response of single neurons is often normalized in the sense that it is divided by a weighted sum of the responses of neighboring neurons [38, 22, 14]. This principle is often termed divisive normalization and closely related to the normalization procedure which ensures contrast invariance in the proposed measures (cf. (22), (44), (55), and (64)).

The implementations of the developed edge, ridge, and blob measures for digital images could also be modeled as artificial neural networks, where the applied even- and odd-symmetric digital filters define a convolutional layer, the absolute value is used as a non-linearity and max-pooling as well as divisive normalization are applied to obtain the final normalized values. It would be interesting to see whether a neural network architecture that is based on similar design principles would yield computations that are related to the formulas proposed in this work when trained for the detection of edges, ridges, or blobs. On one hand, this could lead to new ideas to improve the proposed measures in order to to better handle difficulties such as corner points, intersections of ridges, or densely grouped blobs. On the other hand, neural network architectures such as [6], which have already been trained with great success for similar tasks, could potentially also be further improved by considering the computational principles behind the proposed measures. Ideally, such work would lead to hybrid approaches that combine the best of both worlds in the sense that they yield state-of-the-art performance results for specific applications that are nowadays often only achievable with learning-based methods while also being robust in applications where pure learning is not possible due to a lack of training data.

We have demonstrated that the proposed measures can yield state-of-the-art detection performances when considering sets of both clean and distorted artificial images that are associated with reliable ground truths. Furthermore, we have given examples of how the developed feature detectors can be applied in tasks like the characterization and detection of blood vessels in retinal images, or the automated counting of cell colonies in a Petri dish. The obtained experimental results are promising and suggest a high applicability of the proposed measures in a wide range of diverse applications. However, they also reveal a number of open issues that leave ample room for future improvements.

The present work only considers two-dimensional images. However, similar measures and construction principles of even- and odd-symmetric generators could also be applied in a three-dimensional setting to detect and characterize surfaces, planes, filaments, or filled three-dimensional objects. The definition as well as the mathematical and experimental analysis of such measures could also be a fruitful topic of future research.

8 Acknowledgments

We thank J.-L. Starck for suggesting that we extend the framework to the case of blob detection and M. Schäfer for insightful discussions on the subject of $\alpha$-molecules. We are grateful to Q. Geissmann for kindly sharing a dataset of images of cell colonies and the corresponding ground truth counts. We furthermore thank D. Labate for sharing the Matlab implementation of his shearlet-based edge detector.

R. Reisenhofer acknowledges financial support from the Austrian Science Fund (FWF) under grant number M 2528.