Subpixel Localization of Isolated Edges and Streaks in Digital Images

Abstract

Many modern sensing systems rely on the accurate extraction of measurement data from digital images. The localization of edges and streaks in digital images is an important example of this type of measurement, with these techniques appearing in many image processing pipelines. Several approaches attempt to solve this problem at both the pixel level and subpixel level. While the subpixel methods are often necessary for applications requiring best-possible accuracy, they are often susceptible to noise, use iterative methods, or require pre-processing. This work investigates a unified framework for subpixel edge and streak localization using Zernike moments with ramp-based and wedge-based signal models. The method described here is found to outperform the current state-of-the-art for digital images with common signal-to-noise ratios. Performance is demonstrated on both synthetic and real images.

1. Introduction

Digital images frequently contain valuable information about the real-world objects observed by a camera, telescope, or other optical system. This information may be used by a sensing system to understand, interpret, monitor, or analyze the properties of objects contained within the scene. Oftentimes, such systems attempt to distill the dense and complex information content of an image into a sparse set of simple and descriptive primitives—with edges and streaks being especially common examples. Making use of these primitives requires knowledge of their location in the image. Although pixel-level edge/streak localization is often adequate, some applications demand higher accuracy and motivate the need for subpixel localization.

Pixel-level edge localization is a ubiquitous image processing task, with a variety of techniques that can be found in almost every introductory text on image processing. Some popular classical methods are those of Sobel [1], Prewitt [2], Marr–Hildreth [3], and Canny [4], although there are many more. Motivated largely by problems in image segmentation, there has also been recent interest in edge detection and localization using deep learning [5], with notable contemporary examples including DeepEdge [6], DeepContour [7], holistic edge detector (HED) [8], and crisp edge detection (CED) [9].

There are also several methods available for subpixel edge localization. Many of these subpixel methods operate by refining a pixel-level edge guess into a subpixel-level estimate. The approaches for achieving such a subpixel correction vary, but generally belong to one of four different categories: moment-based [10,11], least-squares fitting [12], partial area effect [13], and interpolation [14].

In this work, we treat both edges and streaks as a local (and not global) concept, with each being identifiable by the 2D intensity pattern within a small image patch around a particular image point. Edge points are generally identified by finding pixels possessing a large intensity gradient (attempting to describe image points where there is thought to be an intensity discontinuity). A streak point is identified by finding pixels belonging to a bright (or dark) 1D path against a dark (or bright) background. In either case (edges or streaks), we seek only to localize isolated edge/streak points in this work.

This work presents a localization framework that is equally suitable to finding the subpixel location of edge and streak points within a digital image. Our method belongs to the moment-based category of techniques. We refine recent work on improved edge localization with Zernike moments [15] and then extend this approach to the closely related problem of streak localization. Unlike conventional subpixel edge localization methods using Zernike moments that assume an intensity step function [11], we model the underlying edge intensity function as a linear ramp. We use the same approach to model the underlying streak intensity function as a triangular wedge. This is the first application of Zernike moments (that the authors know of) to the subpixel localization of streaks in a digital image. The framework presented here is computationally efficient, non-iterative, and can be used within most imaging pipelines.

The remainder of this work is organized as follows. Section 2 introduces the coordinate frames and scaling conventions that are used in Section 3 to construct Zernike moments on local image patches. Section 4 describes how to use these Zernike moments for the subpixel localization of both edges and streaks. Performance of this approach is then demonstrated quantitatively on synthetic images (Section 5) and qualitatively on real images (Section 6).

2. Coordinate Frames and Conventions

Suppose that we have a digital image with N rows and M columns, with pixel intensity values stored in a $N\times M$ array (for a monochrome image). Define the u − v coordinate system with the origin in the upper lefthand corner such that pixel centers occur at integer values of u and v. The u-direction is to the right (corresponding to column number) and the v-direction is down (corresponding to the row number). We presume in this work that a different algorithm (e.g., Sobel [1], Canny [4]) has already produced pixel-level estimates for either an edge or streak location. Assuming such an algorithm has detected m such pixel locations, we denote the set of pixel-level guesses as ${\{{\tilde{u}}_i,{\tilde{v}}_i\}}_i^m\subset {{\mathbb{Z}}^{\ast }}^2$ (where ${\mathbb{Z}}^{\ast }$ is the set of non-negative integers).

The algorithms presented in this work use a small image patch (e.g., $5\times 5$ or $7\times 7$) centered about a pixel-level estimate of an edge or streak location to compute a small correction to that feature’s location. The result is subpixel-level localization of a point belonging to an edge or streak. Furthermore, the moment-based methods to be discussed in Section 3 require the signal to be contained within the unit circle. Thus, data within each small image patch must be scaled to lie within the unit circle. For the i-th patch,

$$\begin{array}{c}\overline{u}=\frac{2}{N_p}(u-{\tilde{u}}_i)\overline{v}=\frac{2}{N_p}(v-{\tilde{v}}_i),\end{array}$$

(1)

where $N_p$ is the size of the image patch (e.g., $N_p=5$ is a $5\times 5$ patch). We generally constrain $N_p$ to be an odd integer, such that the pixel-level guess occurs at the center of the patch. This scaling ensures that $\parallel \overline{u}\parallel \le 1$ and $\parallel \overline{v}\parallel \le 1$ for every point within the square patch, and that ${\overline{u}}^2+{\overline{v}}^2\le 1$ within the inscribed circle.

We also find it convenient to define a rotated version of the $\overline{u}$ − $\overline{v}$ coordinate frame with an orientation dictated by the local normal of the edge or streak. Define a frame with coordinate axes ${\overline{u}}^{\prime }$ and ${\overline{v}}^{\prime }$ that are rotated by an angle $\psi$ relative to the unprimed frame (Figure 1) such that the ${\overline{u}}^{\prime }$ direction is parallel with the local edge/streak normal and ${\overline{v}}^{\prime }$ is parallel to the edge/streak tangent. The direction ${\overline{u}}^{\prime }$ is chosen to be positive in the direction from dark to bright for an edge. Alternatively, for streaks, the positive ${\overline{u}}^{\prime }$ direction is chosen to be from the patch center towards the streak’s center. Thus, by construction, the correction from the pixel-level guess to the subpixel streak location is a small positive update along the ${\overline{u}}^{\prime }$ direction. The subpixel update along the ${\overline{u}}^{\prime }$ direction may be either positive or negative for an edge.

3. Computation of Zernike Moments in Digital Images

It is well-established that image moments are a useful tool for compactly describing the shape of the 2D intensity pattern within an image patch using only a small number of parameters. In general, 2D moments are a weighted average of the 2D signal value, with the weights for a particular moment coming from its corresponding basis function. That is, given a basis function $P_{nm}(u,v)$, the corresponding moment of the arbitrary 2D signal $f(u,v)$ is computed as

$$\int \int P_{nm}(\overline{u},\overline{v})f(\overline{u},\overline{v})d\overline{u}d\overline{v}$$

(2)

Since we will be computing moments within small image patches, we have chosen to express all functions in terms of the scaled pixel coordinates $\{\overline{u},\overline{v}\}$ as defined in Equation (1).

The choice of basis functions $P_{nm}$ is somewhat arbitrary, although it is desirable that the chosen set is both complete and orthogonal. In the case of edge or streak localization, we are looking for basis function sets defined within the unit disk. If $P_{nm}$ is chosen to be a polynomial in two variables, there are an infinite number of complete orthogonal sets [16], with the Zernike polynomials being the most commonly used.

3.1. Zernike Polynomials

Zernike polynomials, originally developed to aid in the study of spherical aberrations in optical lenses [17], have since found uses for a broad array of applications [11,15,18,19,20,21,22]. The Zernike polynomials may be written in either Cartesian or polar coordinates, with the polar form being the most commonly used [11],

$$\begin{array}{c}{P}_{nm}(\rho ,\theta )=R_{nm}\left(\rho \right)exp\left(jm\theta \right)\end{array}$$

(3)

where $j=\sqrt{-1}$ and

$${\rho }^2={\overline{u}}^2+{\overline{v}}^2$$

(4)

$$R_{nm}\left(\rho \right)=\sum _{s=0}^{(n-|m\left|\right)/2}\frac{{(-1)}^s(n-s)!{\rho }^{n-2s}}{s!\left(\frac{n+\left|m\right|}{2}-s\right)!\left(\frac{n-\left|m\right|}{2}-s\right)!}$$

(5)

These polynomials form a complete set over a continuous space contained within the unit circle. The 1D radial polynomials, $R_{nm}\left(\rho \right)$, and their corresponding 2D Zernike polynomials, $P_{nm}(u,v)$, may be computed for a few common combinations of n and m,

$$\begin{array}{}\hfill \mathrm{(6)}& R_{00}\left(\rho \right)=1\hfill & \hfill \Rightarrow & P_{00}(\overline{u},\overline{v})=1\hfill \hfill \mathrm{(7)}& R_{11}\left(\rho \right)=\rho \hfill & \hfill \Rightarrow & P_{11}(\overline{u},\overline{v})=\overline{u}+j\overline{v}\hfill \hfill \mathrm{(8)}& R_{20}\left(\rho \right)=2{\rho }^2-1\hfill & \hfill \Rightarrow & P_{20}(\overline{u},\overline{v})=2{\overline{u}}^2+2{\overline{v}}^2-1\hfill \hfill \mathrm{(9)}& R_{22}\left(\rho \right)={\rho }^2\hfill & \hfill \Rightarrow & P_{22}(\overline{u},\overline{v})={\overline{u}}^2+{\overline{v}}^2\hfill \end{array}$$

where the order n and repetition m [15] (or angular dependence [23]) can assume any values that satisfy

$$\begin{array}{c}n\ge \left|m\right|\ge 0,\end{array}$$

(10)

$$\begin{array}{c}n-meven.\end{array}$$

(11)

It is straightforward to show that Zernike polynomials are orthogonal under an $L_2$-inner product,

$${〈P_{\left(\alpha \right)},P_{\left(\beta \right)}〉}_{L_2}=\underset{u^2+v^2\le 1}{\int \int }{P}_{\left(\alpha \right)}{P}_{\left(\beta \right)}^{\ast }dudv=Q_{\left(\alpha \right)}{\delta }_{\alpha \beta }$$

(12)

where $P_{\left(\alpha \right)}$ and $P_{\left(\beta \right)}$ are two arbitrary polynomials of the set, $P_{\left(\beta \right)}^{\ast }$ is the complex conjugate of $P_{\left(\beta \right)}$, and ${\delta }_{\alpha \beta }$ is the Kronecker delta function. Additionally, $Q_{\left(\alpha \right)}$ is the normalization coefficient and may be computed as [23]

$$\begin{array}{c}{Q}_{nm}=\frac{\pi }{n+1}\end{array}$$

(13)

3.2. Zernike Moments for a Continuous 2D Signal

Zernike moments are formed by using the Zernike polynomials from Equation (3) as the basis functions in the 2D moment equation (Equation (2)). We express such a moment as

$$\begin{array}{c}{Z}_{nm}=\frac{1}{Q_{nm}}\underset{u^2+v^2\le 1}{\int \int }{P}_{nm}(u,v)f(u,v)\mathrm{d}u\mathrm{d}v\end{array}$$

(14)

although we often find that scaling with the normalization coefficient is not required,

$$\begin{array}{c}{A}_{nm}=\underset{u^2+v^2\le 1}{\int \int }{P}_{nm}(u,v)f(u,v)\mathrm{d}u\mathrm{d}v,\end{array}$$

(15)

This, of course, leads to the simple scaling relation

$$\begin{array}{c}{Z}_{nm}=A_{nm}/Q_{nm}\end{array}$$

(16)

3.3. Rotational Properties of Zernike Moments

Zernike moments of repetition $m=0$ are rotationally invariant, as the value of the moment $A_{nm}$ is unaffected by the orientation of the underlying signal relative to the $\overline{u}$ − $\overline{v}$ coordinate system. For other values of m (i.e., for $m>0$), we find that the moment $A_{nm}$ changes as the orientation of the underlying signal changes.

Consider, for example, the moment $A_{nm}$ for a particular image patch as computed in the $\overline{u}$ − $\overline{v}$ frame. Now consider the moment $A_{nm}^{\prime }$ for this same image patch as computed in the ${\overline{u}}^{\prime }$ − ${\overline{v}}^{\prime }$ frame that has been rotated by an angle $\psi$ relative to the unprimed frame (see Figure 1). Noting that ${\theta }^{\prime }=\theta -\psi$, it is clear from Equations (3) and (15) that

$$A_{nm}^{\prime }=A_{nm}exp(-jm\psi )$$

(17)

It is this relation that will ultimately allow us to determine the orientation of an edge or streak from the moment $A_{11}$.

3.4. Zernike Moments for a Digital Image

A digital image, $I(u,v)$, is a quantized representation of the continuous signal $f(u,v)$. The image $I(u,v)$ is presumed to be an array of digital numbers, with integer intensity values (e.g., 0–65,535 for a 16-bit image) occurring at integer values of u and v.

In this case, we approximate the Zernike moment integral from Equation (15) with a double summation. Therefore, assuming a local image patch of size $N_p\times N_p$ centered at a pixel-level edge/streak guess of $\{{\tilde{u}}_i,{\tilde{v}}_i\}$, one may compute the moment as

$$A_{nm}({\tilde{u}}_i,{\tilde{v}}_i)\approx \sum _{k=-p}^p\sum _{s=-p}^{p}I({\tilde{u}}_i+s,{\tilde{v}}_i+k)M_{nm}(p+s,p+k)$$

(18)

where $p=(N_p-1)/2$ is a non-negative integer (since $N_p$ is an odd integer greater than one). The mask $M_{nm}$ is an $N_p\times N_p$ matrix of values found by the integration of $P_{nm}$ over the corresponding pixel and within the patch’s inscribed circle. Values of $M_{11}$ and $M_{20}$ are shown for a $5\times 5$ and $7\times 7$ mask in [15]. It is observed that Equation (18) is simply an image correlation, such that one may compute the moment everywhere in the image according to

$$A_{nm}=M_{nm}\ast I$$

(19)

where ∗ is the 2D correlation operator.

The edge and streak localization methods presented here will ultimately only use the moments $A_{11}$ and $A_{20}$. Of note is that $M_{20}$ is real valued such that

$$A_{20}=M_{20}\ast I$$

(20)

We observe, however, that $M_{11}$ is complex valued,

$$M_{11}=\mathrm{Re}\left[M_{11}\right]+j\mathrm{Im}\left[M_{11}\right]$$

(21)

Fortunately, given the structure of $M_{11}$, one only needs to keep track of the real component in practice since [15]

$$\mathrm{Re}\left[A_{11}\right]=\mathrm{Re}\left[M_{11}\right]\ast I$$

(22)

$$\mathrm{Im}\left[A_{11}\right]=\mathrm{Re}{\left[M_{11}\right]}^T\ast I$$

(23)

Thus, we may compute all the necessary moments through three simple image correlations (which, in practice, only need to be computed at the pixel-level edge or streak locations and not at every point in the image).

4. Moment-Based Edge and Streak Localization

The same procedure may be used for both edge and streak localization. In both cases, the image data in a small $N_p\times N_p$ image patch around a pixel-level edge/streak guess is scaled according to Equation (1) and the Zernike moments $A_{11}$ and $A_{20}$ are computed (Equations (20), (22) and (23)). These moments are used to compute the edge/streak orientation ($\psi$) and the distance along this direction by which the pixel-level edge/streak guess should be adjusted (ℓ). Consequently, both the edge and the streak are corrected to subpixel accuracy by

$$\begin{array}{c}\left[\begin{array}{c}\overline{u}\\ \overline{v}\end{array}\right]=\ell \left[\begin{array}{c}cos\psi \\ sin\psi \end{array}\right]\end{array}$$

(24)

which, after rearranging Equation (1), yields the correction we seek in practice

$$\begin{array}{c}\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}\tilde{u}\\ \tilde{v}\end{array}\right]+\frac{N_p\ell }{2}\left[\begin{array}{c}cos\psi \\ sin\psi \end{array}\right]\end{array}$$

(25)

The orientation $\psi$ of both edges and streaks is found in the same way and using the same equation. The difference between the edge and streak correction is simply how the Zernike moments are used to compute ℓ.

4.1. Computing Edge or Streak Orientation

Determining the normal direction to an edge or streak is achieved in the exact same manner, with the final equation being equivalent for both. By construction, and as can be seen from Figure 1, the intensity value is only a function of ${\overline{u}}^{\prime }$ (i.e., not a function of ${\overline{v}}^{\prime }$) for both the edge and the streak. We see immediately from the form of $P_{11}$ in Equation (6) that

$$\mathrm{Im}\left[A_{11}^{\prime }\right]=0$$

(26)

Thus, recalling that $exp(-jm\psi )=cos\left(m\psi \right)-jsin\left(m\psi \right)$, we may rewrite Equation (17) as (for $m=n=1$)

$$A_{11}^{\prime }=A_{11}\left[cos\left(\psi \right)-jsin\left(\psi \right)\right]$$

(27)

such that

$$\mathrm{Re}\left[A_{11}^{\prime }\right]=\mathrm{Re}\left[A_{11}\right]cos\left(\psi \right)+\mathrm{Im}\left[A_{11}\right]sin\left(\psi \right)$$

(28)

$$\mathrm{Im}\left[A_{11}^{\prime }\right]=\mathrm{Im}\left[A_{11}\right]cos\left(\psi \right)-\mathrm{Re}\left[A_{11}\right]sin\left(\psi \right)=0$$

(29)

The streak orientation may be found using the equation for the imaginary component of $A_{11}^{\prime }$. Observing that

$$\mathrm{Im}\left[A_{11}\right]cos\left(\psi \right)=\mathrm{Re}\left[A_{11}\right]sin\left(\psi \right)$$

(30)

we find that the orientation of the streak is computed in terms of the moment $A_{11}$ (computed using Equations (22) and (23)) as

$$\psi =\mathrm{atan}2\left(\mathrm{Im}\left[A_{11}\right],\mathrm{Re}\left[A_{11}\right]\right)$$

(31)

This relation has been known for some time for edges [11,15]. Although obvious within the present framework, this represents the first extension of Equation (31) to streaks (of which the authors are aware).

4.2. Computing ℓ for Edges

An edge is generally understood to describe a discontinuity in image intensity in one direction, with little intensity change in the direction orthogonal to this discontinuity. Real grayscale images, however, rarely exhibit a true intensity discontinuities. Instead, image blur and pixel quantization cause the intensity change rapidly over a small distance (a few pixels). Thus, we seek areas of high intensity gradient rather than true discontinuities. It has long been known [10] that using a step function for the edge model within a moment-based subpixel edge localization algorithm produces a biased edge update if the image is blurred. This was one of the motivations for introducing a ramp edge model in [15].

In many practical image processing problems, the point spread function (PSF) due to camera defocus and other optical effects is well modeled as a 2D Gaussian [24]. Consequently, the streak associated with a crisp edge (a true discontinuity) may be blurred according to

$$I_{blur}\approx K_G\ast I$$

(32)

where I is the perfectly crisp image, $K_G$ is the Gaussian kernel, and $I_{blur}$ is the blurred image. The one-dimensional intensity profile taken perpendicular to the edge is sometimes referred to as the edge spread function (ESF), which will generally take the shape of a sigmoid function. To avoid the mathematical complexities of the sigmoid function within the Zernike moment integrals, it was observed in [15] that a linear ramp provides an adequate engineering approximation for most practical cases. The objective, therefore, is to relate the width of the linear ramp ($2w$ full-width, see Figure 2) with the width of the Gaussian kernel approximating the camera PSF ($\sigma$). We do this using the linear relationship

$$w\approx k_{edge}\sigma$$

(33)

where $k_{edge}$ is the scaling we seek. In [15], it was suggested to select $k_{edge}=1.66$. We performed a more comprehensive study and found that choosing $k_{edge}=1.80$ produced superior performance, especially as the SNR became very large. In general, we found reduced sensitivity to the choice of $k_{edge}$ as the images became noisier (lower SNR).

Therefore, we choose to model an edge as a ramp, whose intensity changes linearly between a background intensity (h) and a foreground intensity ($h+k$). The midpoint of this transition is defined to occur at a distance ℓ from the image patch center and has a width of $2w$. Since we are using Zernike moments, we define all these quantities within the unit disk (Figure 2). By choosing to define the edge in the ${\overline{u}}^{\prime }$ − ${\overline{v}}^{\prime }$ frame, it is straightforward to write the intensity as a function of ${\overline{u}}^{\prime }$ only,

$$I_{edge}({\overline{u}}^{\prime },{\overline{v}}^{\prime })=\left\{\begin{array}{cc}h\hfill & {\overline{u}}^{\prime }\le \ell -w\hfill \\ h+k\frac{{\overline{u}}^{\prime }-(\ell -w)}{2w}\hfill & \ell -w<{\overline{u}}^{\prime }\le \ell +w\hfill \\ h+k\hfill & \ell +w<{\overline{u}}^{\prime }\hfill \end{array}\right.$$

(34)

Using this ramp edge model, it is possible to analytically solve the double integral in the moment equation from Equation (15) in the edge-aligned (i.e., primed) frame. We do this for the moments $A_{11}^{\prime }$ and $A_{20}^{\prime }$, leading to

$$\begin{array}{c}\begin{array}{c}\hfill A_{11}^{\prime }=\frac{k}{24w}[3arcsin{\ell }_2+(5-2{\ell }_2^2){\ell }_2{B}_2\\ \hfill -3arcsin{\ell }_1-(5-2{\ell }_1^2){\ell }_1{B}_1]\end{array}\end{array}$$

(35)

$$\begin{array}{c}{A}_{20}^{\prime }=A_{20}=\frac{k}{15w}\left[B_1^5-B_2^5\right]\end{array}$$

(36)

where

$$\begin{array}{c}{\ell }_1=\ell -w\mathrm{and}{\ell }_2=\ell +w\end{array}$$

(37)

and

$$\begin{array}{c}{B}_1=\sqrt{1-{\ell }_1^2}\mathrm{and}{B}_2=\sqrt{1-{\ell }_2^2}\end{array}$$

(38)

Looking at the expressions for $A_{11}^{\prime }$ and $A_{20}^{\prime }$, it is immediately evident that the the intensity-dependent variable k (which describes the magnitude of the intensity change across the edge) cancels out if one considers the ratio, $Q_E$

$$\begin{array}{c}{Q}_E=\frac{A_{20}^{\prime }}{A_{11}^{\prime }}=\frac{A_{20}}{A_{11}^{\prime }}\end{array}$$

(39)

In many cases, the edge width w is known (e.g., from the imaging system point spread function), such that $Q_E$ is a function of only ℓ. Although the analytic expression for $Q_E$ is rather cumbersome, it was found in [15] that

$$\begin{array}{c}{Q}_E\approx \ell \left[1-(1+\ell /2)w^2\right]\end{array}$$

(40)

which may be rearranged to solve for the unknown ℓ

$$\begin{array}{c}{\widehat{\ell }}_E\approx \frac{1-w^2-\sqrt{{(w^2-1)}^2-2w^2{Q}_E}}{w^2}\end{array}$$

(41)

Note that the ratio $Q_E$ is easy to compute in practice from the raw image moments found in a digital image,

$$\begin{array}{c}{Q}_E=\frac{A_{20}^{\prime }}{A_{11}^{\prime }}=\frac{A_{20}}{A_{11}^{\prime }}=\frac{A_{20}}{\mathrm{Re}\left[A_{11}\right]cos\psi +\mathrm{Im}\left[A_{11}\right]sin\psi }\end{array}$$

(42)

where $A_{20}$ from Equation (20), $\mathrm{Re}\left[A_{11}\right]$ from Equation (22), $\mathrm{Im}\left[A_{11}\right]$ from Equation (23), and $\psi$ from Equation (31). Thus, with w and $Q_E$ known, Equation (41) may be used to solve for ${\widehat{\ell }}_E$ for a given image patch.

4.3. Computing ℓ for Streaks

As a natural extension to the ideal step-function ESF, we model the ideal line spread function (LSF) as an impulse (where the LSF is defined as the 1D intensity profile perpendicular to the streak). As before, the perfectly crisp image is blurred with a Gaussian kernel, thus spreading out the line intensity, with the resulting LSF being Gaussian PDF. Rather than deal with the mathematical complexities of the Gaussian PDF, we choose to model the streak PSF as a wedge. To make practical use of the wedge model, it is necessary to determine the relationship between the wedge width (w, see Figure 3) and the Gaussian kernel width ($\sigma$),

$$w\approx k_{streak}\sigma$$

(43)

where $k_{streak}$ is the parameter we seek. We found that choosing $k_{streak}=0.90$ provided the best results, with low SNR images exhibiting less sensitivity to the exact choice of this parameter.

The small image patch centered about the pixel-level guess is assumed to have a constant background intensity of h and contain a streak of intensity $h+k$. The wedge has a full width of $2w$ with a peak intensity occurring at a distance ℓ from the image patch (or disk) center. The sides of the wedge are linear ramps transitioning between the background and the streak’s ridgeline. This is shown pictorially on the unit disk in Figure 3.

As with the edge, we choose to define the streak model in the ${\overline{u}}^{\prime }$ − ${\overline{v}}^{\prime }$ frame such that the intensity is a function of ${\overline{u}}^{\prime }$ only (and not a function of ${\overline{v}}^{\prime }$),

$$I_{streak}({\overline{u}}^{\prime },{\overline{v}}^{\prime })=\left\{\begin{array}{cc}h\hfill & {\overline{u}}^{\prime }\le \ell -w\hfill \\ h+k\frac{{\overline{u}}^{\prime }-(\ell -w)}{w}\hfill & \ell -w<{\overline{u}}^{\prime }\le \ell \hfill \\ h+k-k\frac{{\overline{u}}^{\prime }-\ell }{w}\hfill & \ell <{\overline{u}}^{\prime }\le \ell +w\hfill \\ h\hfill & \ell +w<{\overline{u}}^{\prime }\hfill \end{array}\right.$$

(44)

The analytical value of $A_{11}^{\prime }$ and $A_{20}^{\prime }$ may be found by evaluating the double integral from Equation (15)

$$\begin{array}{c}\begin{array}{c}\hfill A_{11}^{\prime }=\frac{k}{12w}[6arcsin\ell -3arcsin{\ell }_1-3arcsin{\ell }_2\\ \hfill +2C\ell (5-2{\ell }^2)\\ \hfill -B_1{\ell }_1(5-2{\ell }_1^2)\\ \hfill -B_2{\ell }_2(5-2{\ell }_2^2)]\end{array}\end{array}$$

(45)

and

$$\begin{array}{c}\begin{array}{c}\hfill A_{20}^{{}^{\prime }}=\frac{2k}{15w}\left[B_1^5+B_2^5-2C^5\right],\end{array}\end{array}$$

(46)

where $B_1$ and $B_2$ are from Equation (38) and

$$\begin{array}{c}C=\sqrt{1-{\ell }^2}\end{array}$$

(47)

The ratio of $A_{20}^{\prime }$ to $A_{11}^{{}^{\prime }}$ eliminates k, thus providing a function of only ℓ and w,

$$\begin{array}{c}{Q}_S(\ell ,w)=\frac{A_{20}^{\prime }}{A_{11}^{\prime }}=\frac{A_{20}}{A_{11}^{\prime }}\end{array}$$

(48)

Assuming the streak width w is known, we seek to rearrange $Q_S$ to solve for the unknown ℓ. The complicated form of $Q_S$ after substitution of Equation (45) and Equation (46) makes finding an analytic solution difficult for arbitrary values of w and ℓ. Fortunately, it is straightforward to find an approximation that is good enough for most practical image processing applications.

We know that streaks are thin, so it is instructive to explore what happens to $Q_S$ as $w\to 0$. We find that the limit does permit a simple analytic solution,

$$\begin{array}{c}{Q}_{S_0}=\underset{w\to 0}{lim}{Q}_S=\frac{4{\ell }_{S_0}^2-1}{3{\ell }_{S_0}}\end{array}$$

(49)

which may be solved for for ${\ell }_{S_0}$

$$\begin{array}{c}{\widehat{\ell }}_{S_0}=\frac{3}{8}{Q}_{S_0}\pm \sqrt{{\left(\frac{3}{8}{Q}_{S_0}\right)}^2+\frac{1}{4}}\end{array}$$

(50)

To choose the correct root, observe that

$$\begin{array}{c}\underset{{\ell }_{S_0}\to 0^{+}}{lim}{Q}_{S_0}=\underset{{\ell }_{S_0}\to 0^{+}}{lim}\frac{4{\ell }_{S_0}^2-1}{3{\ell }_{S_0}}=-\infty \end{array}$$

(51)

where we know to choose the right limit since ${\ell }_S\ge 0$ by construction. Thus we seek the root that is approximately zero when $Q_{S_0}$ is a large negative number, which only happens when the plus sign is chosen in Equation (50). Therefore,

$$\begin{array}{c}{\widehat{\ell }}_{S_0}=\frac{3}{8}{Q}_{S_0}+\sqrt{{\left(\frac{3}{8}{Q}_{S_0}\right)}^2+\frac{1}{4}}\end{array}$$

(52)

This analytic result may be generalized to the situation where $w>0$, which does not appear to permit an exact analytic solution. Therefore, we write a parameterized expression for ${\widehat{\ell }}_S$ that simplifies exactly to the form of Equation (52) when $w=0$ and fit the parameters in a least squares sense. Using this approach, consider a model of the form

$$\begin{array}{c}\begin{array}{c}\hfill {\widehat{\ell }}_S\approx \sqrt{a_1^2{Q}_S^2+a_2{Q}_{S}w+a_3{w}^2+a_{4}w+a_5}+a_6{Q}_S+a_{7}w+a_8\end{array}\end{array}$$

(53)

We found the terms associated with $a_1$, $a_3$, $a_5$, and $a_6$ to dominate the estimate of ${\ell }_S$, with the remaining terms contributing relatively little. Furthermore, it was found that $a_1\approx a_6$ regardless of the test set-up. Therefore, discarding the unimportant terms and letting $a_1=a_6$, we performed a three parameter fit for the streak correction of the form

$$\begin{array}{c}\begin{array}{c}\hfill {\widehat{\ell }}_S\approx \sqrt{a_1^2{Q}_S^2+a_3{w}^2+a_5}+a_1{Q}_S.\end{array}\end{array}$$

(54)

The result of the least squares fit found the values of $a_1$ and $a_5$ to exactly match the analytically derived coefficients for ${\ell }_{S_0}$ in Equation (52)

$$\begin{array}{c}{a}_1=3/8,a_5=1/4\end{array}$$

(55)

and empirically found that

$$\begin{array}{c}{a}_3\approx -1/10\end{array}$$

(56)

Therefore, we may write the empirically derived expression for the streak update for arbitrary w as

$$\begin{array}{c}{\widehat{\ell }}_S=\frac{3}{8}{Q}_S+\sqrt{{\left(\frac{3}{8}{Q}_S\right)}^2-\frac{1}{10}{w}^2+\frac{1}{4}}\end{array}$$

(57)

As with the edge update, note that the streak ratio $Q_S$ is easy to compute in practice from the raw image moments found in a digital image,

$$\begin{array}{c}{Q}_S=\frac{A_{20}^{\prime }}{A_{11}^{\prime }}=\frac{A_{20}}{A_{11}^{\prime }}=\frac{A_{20}}{\mathrm{Re}\left[A_{11}\right]cos\psi +\mathrm{Im}\left[A_{11}\right]sin\psi }\end{array}$$

(58)

where $A_{20}$ from Equation (20), $\mathrm{Re}\left[A_{11}\right]$ from Equation (22), $\mathrm{Im}\left[A_{11}\right]$ from Equation (23), and $\psi$ from Equation (31). Thus, with w and $Q_S$ known, Equation (57) may be used to solve for ${\widehat{\ell }}_S$ for a given image patch.

Observe that $Q_E$ and $Q_S$ are the same moment ratio, $A_{20}^{\prime }/A_{11}^{\prime }$; hence, the equations to compute these ratios from the raw image moments are the same (compare Equation (42) and Equation (58)). What differs is the assumption of the underlying signal (a ramp or an edge), leading to a different relationship (Equation (41) or Equation (57)) between the moment ratio and the subpixel location of the edge or streak.

5. Numerical Validation on Synthetic Images

The performance of the edge and streak localization methods presented in this work were quantitatively evaluated using synthetic images. We find synthetic images to be especially useful in this context since the true continuous location of every image feature is known. The perfectly known continuous underlying signal may be blurred to simulate camera defocus and quantized (both spatially and in intensity) to simulate differing image resolutions. Further, noise may be added with a prescribed intensity, allowing the unambiguous evaluation of performance as a function of signal-to-noise ratio (SNR). This is important, as the localization of edges and streaks is known to become more challenging as SNR decreases [25,26]. Of particular note is that our new streak localization method works for 1D paths of arbitrary shape, whereas most existing streak detection algorithms—especially for faint (low SNR) streaks—presume the streaks are straight lines.

For the examples presented here, perfect images were blurred by using a Gaussian point spread function (PSF). After blurring, zero-mean Gaussian noise was added to achieve the specified SNR.

5.1. Synthetic Images with Edges

5.1.1. Ideal Edge Localization Performance

It is important to quantify the error associated with the approximations used to arrive at the analytic edge update given in Equation (41). Therefore, as a bounding case, suppose that we perfectly compute the Zernike moments for a noise-free continuous signal. In this situation, the error in ${\widehat{\ell }}_E$ is given by the contours in Figure 4 for different situations. These contours visually demonstrate the performance improvement afforded by switching from the step-function edge model (red contours) to the ramp edge model (black contours). The results shown here are identical to the observations of Christian in [15].

5.1.2. Digital Image Edge Localization Performance

Our method performed favorably to other existing techniques when processing synthetic digital imagery. This was assessed through a Monte Carlo analysis where we evaluated performance of different algorithms for images having varying amounts of blur and noise. Figure 5 shows edge localization error with our technique (black contours) compared against the moment-based solution with a step-function edge model [11] and the partial area effect (PAE) [13]. Results for both of the two moment-based methods shown here assume a $5\times 5$ pixel mask.

Note that the PAE method from [13] was chosen as one of the two comparison methods in Figure 5 since this represents the current state-of-the-art. Indeed, this method has recently been used for the subpixel localization of edges in a wide variety of applications [27,28,29].

We observe that the PAE algorithm produced nearly perfect edge localization in cases with no noise (infinite SNR; off the right-hand side of Figure 5). The Zernike moment methods tended to perform better than the PAE method as noise increased (as SNR decreased; towards the left-hand side of Figure 5). The method presented in this work outperforms the PAE method for most real-life SNR values.

Example performance of our subpixel edge localization algorithm in different noise/blur cases is shown in Figure 6. This example shows localization of the edge of a circle. Clear improvement is evident in all cases, as the algorithm moves the pixel-level edge guess (red ×) towards the true edge location (black line). We know the true edge location since these are synthetic images.

5.2. Synthetic Images with Streaks

5.2.1. Ideal Streak Localization Performance

As with the case of edges, we begin the numerical assessment of our subpixel streak localization method by considering the case of a continuous signal. This allows us to directly quantify the error associated with the approximations used to arrive at the analytic expression in Equation (57). We considered all reasonably plausible combinations of streak location (ℓ) and streak width (w) and produced contours of errors in the estimate ${\widehat{\ell }}_S$, as shown in Figure 7. These errors are low enough to be negligible when applied to a pixelated image.

5.2.2. Digital Image Streak Localization Performance

Our Zernike moment method also performed well in the subpixel localization of streaks. We performed a Monte Carlo analysis where streak localization error was recorded for varying amounts of image blur and noise. The results are shown as contours in Figure 8. As expected, localization performance decreases with increased noise and blur.

Example performance of our subpixel streak localization algorithm in different noise/blur cases is shown in Figure 9. This example shows localization of a circular streak. Clear improvement is evident in all cases, as the algorithm moves the pixel-level streak guess (red ×) towards the true streak location (black line). We know the true streak location since these are synthetic images.

6. Validation on Real Data

After confirming that estimated edge and streak locations agree with the truth in simulated images, we apply our method to real digital images. As these real-world images do not provide perfect subpixel knowledge of the edge or streak location, verifying results is from visual inspection and is largely qualitative.

It is important to remember that the algorithm presented here only performs the subpixel localization (i.e., correction) on pixel-level location guesses (e.g., using Sobel [1], Canny [4], or other method); any edges or streaks that the higher-level algorithm fails to identify will not contribute to the final result. Note that these pixel-level guesses may be found automatically or manually. Regardless of how they are found, the subpixel correction discussed in this manuscript is automatic.

This section includes a number of example images with the accompanying results from the methods proposed in this paper. These examples show the raw image on the left-most frame, followed by two sections of the image in grayscale containing streaks or edges of interest (center and right frame). We highlight performance by progressively zooming in on a specific portion of the image (moving left to right), with blue boxes indicating the region-of-interest for the subsequent frame. The middle frame of each example only shows the subpixel estimate overlay (green dots). The right frame of each example shows both the pixel-level guess overlay (red ×) and the subpixel estimate overlay (green dots). The right frame also shows the edge or streak estimates connected by a line to help illustrate the improvement in smoothness naturally produced by the subpixel correction.

Figure 10 shows an application to natural disaster management that illustrates the difference in the shores of the Mississippi River in the aftermath of a flood (bottom) and its normal banks (top). Figure 11 and Figure 12 show the application of the proposed technique for the subpixel localization of common road surface markings (e.g., pedestrian crosswalk markings, lane markings). Figure 13, Figure 14 and Figure 15 show various applications to space exploration. Finally, Figure 16 and Figure 17 highlight the potential use of this method in medical imaging (e.g., tracing the routes of blood vessels in a retinal scan, microscope imaging of tumors). The diversity of example images is intended to emphasize that the techniques presented in this manuscript are application agnostic and can be applied to a wide variety of image processing tasks.

7. Conclusions

Many modern sensing systems rely on the accurate extraction of measurement data from digital images. The localization of edges and streaks in digital images is an important example of this type of measurement, with these techniques appearing in many image processing pipelines.

Zernike moments are powerful tools in image processing and have been used for subpixel edge localization for over 25 years. In this manuscript, we describe a new way to exploit Zernike moment data to produce subpixel edge estimates, resulting in improved localization performance relative to earlier techniques using Zernike moments to achieve this same task. We also show how this same framework can be extended to the task of subpixel localization of streak. As far as the authors know, this represents the first application of Zernike moments to subpixel streak localization.

Correcting a pixel-level guess of either an edge or a streak requires use of only two Zernike moments ($A_{11}$ and $A_{20}$), with both of these moments being computed over a small image patch centered about the pixel-level guess. One of the principal innovations of this work is the use of a linear ramp (for an edge) or triangular wedge (for a streak) signal model. These simplified models make it possible to refine the pixel-level guess to subpixel accuracy using an analytic function of these two moments and knowledge of the edge/streak width. Furthermore, we show this new method to be tolerant to noise and to outperform many existing methods. Performance was quantitatively evaluated on synthetic images (localization error less than 0.1 pixel for both edges and streaks) and qualitatively evaluated on real images. Applications were shown for remote sensing, localization of road markings, space exploration, and medical imaging.