Edge Detection

Edge detection is a problem of fundamental importance in image analysis. In typical images, edges characterize object boundaries and are therefore useful for segmentation, registration, and identification of objects in a scene. In this section, the construction, characteristics, and performance of a number of gradient and zero-crossing edge operators will be presented.

An edge is a jump in intensity. The cross section of an edge has the shape of a ramp. An ideal edge is a discontinuity (i.e., a ramp with an infinite slope). The first derivative assumes a local maximum at an edge. For a continuous image f(x, y), where x and y are the row and column coordinates, respectively, we typically consider the two directional derivatives (∂_x) f(x, y) and (∂_y) f(x, y). Of particular interest in edge detection are two functions that can be expressed in terms of these directional derivatives: the gradient magnitude and orientation. The gradient magnitude is defined as

\(TraditionalForm\`\(\(|\)\(∇\(f(x, 
              y)\)\)\(|\)\)\  = \ \@\(\((∂\_x\ \(f(x, y)\))\)\^2\  + \ \((∂\_y\ \(f(x, y)\))\)\^2\)\),

and the gradient orientation is given by

∠ ∇ f(x, y) = ArcTan( (∂_y) f(x, y) / (∂_x) f(x, 
                  y)).

Local maxima of the gradient magnitude identify edges in f(x, y). When the first derivative achieves a maximum, the second derivative is zero. For this reason, an alternative edge detection strategy is to locate zeros of the second derivatives of f(x, y). The differential operator used in these so-called zero-crossing edge detectors, is the Laplacian

\(TraditionalForm\`Δ\ \(f(x, 
            y)\)\  = \ ∂\_{x, 2}\ \(f(x, 
              y)\)\  + \ ∂\_{y, 2}\ \(f(x, y)\)\).

This loads the package.

<< ImageProcessing`

The following demonstration illustrates the principles of edge detection. We model a one-dimensional edge with a simple smooth function, for example f(t)=Tanh(α t). As the parameter α increases, the slope of the edge gets steeper. Here we use α  =  3.

Plot[1 + Tanh[3 t], {t, -2, 2}, 
        PlotRange -> All];

This evaluates and plots the first and second derivatives of our “edge” function.

\(\(Plot[
        Evaluate[{1 + \ 
              Tanh[3\ t], ∂\_t Tanh[3\ t], ∂\_{t, 2}Tanh[
                3\ t]}], {t, \(-2\), 2}, 
        PlotStyle -> {Hue[0], Hue[1/3], Hue[2/3]}];\)\)

As expected, the result shows that maxima of first derivative (green) and zero-crossings of the second derivative (blue) locate the center of our model edge. It can also be shown that the zero-crossings are independent of the steepness of the transition, while the gradient magnitude is directly related to the edge slope.

In practice, finite difference approximations of first-order directional derivatives are used. These are represented by a pair of masks, say (Subscript[h, x]) and (Subscript[h, y]). Formally these are linear phase FIR filters. A convolution of the image with (Subscript[h, x]) and (Subscript[h, y]) gives two directional derivative images (Subscript[g, x]) and (Subscript[g, y]), respectively. The gradient image is traditionally calculated as  ∇=(Sqrt[((Subscript[g, x])^2)  +  
                ((Subscript[g, y])^2)]), or alternatively using \(TraditionalForm\`∇\ \(\(=\)\(\ \)\(\(|\)\(g\_x\)\(|\)\(\ \)\(+\ \(\(|\)\(g\_y\)\(|\)\)\)\)\)\) [1]. A pixel location is declared an edge location if the value of ∇ (at point x, y) exceeds some threshold. The locations of all edge points constitutes an edge map. The selection of a threshold value is an important design decision that depends on a number of factors, such as image brightness, contrast, level of noise, and even edge direction. Typically, the threshold is selected following an analysis of the gradient image histogram. It is sometimes useful to calculate edge direction information. This is given by ϕ=ArcTan[-((Subscript[g, x])/(Subscript[g, y]))] .

Edge operators.

The method of computing the gradient magnitude is user selectable with option GradientType.

Option of EdgeMagnitude.

A number of commonly used gradient operators are available.

Gradient edge operators.

Here are two common gradient edge operators.

MatrixForm /@ RobertsFilter[]

MatrixForm /@ SobelFilter[]

Note that each of the Sobel edge masks is a combination of a digital differentiator in one of the spatial directions and a smoothing operator in the other.

Here we read one of the example images and consider a subregion of the image.

wing = 
        ImageTake[ImageRead["wing.tif"], -125, 125];

We compute the gradient magnitude and edge direction using the Sobel mask.

{mag, dir} = {EdgeMagnitude[wing, SobelFilter[]], 
          EdgeOrientation[wing, SobelFilter[]]};

Show[GraphicsArray[Graphics /@ {mag, dir}]];

To recover the edges, the gradient image must be segmented using a global or local (i.e., adaptive) threshold operator. The choice of a threshold value determines the resulting segmentation and, therefore, the perceived quality of the edge detector. It is useful to consider the cumulative histogram of the gradient image in selecting an appropriate threshold value. Here we calculate the histogram and cumulative histogram of the gradient magnitude image. In a typical approach, the top 10 percent to 20 percent of the largest gradient values are selected as edge points. This is easily done based on an investigation of the cumulative histogram.

tmp = CumulativeHistogram[ImageHistogram[Round[mag]]];

ListPlot[tmp, PlotJoined -> True, 
        PlotRange -> All];

The threshold can now be estimated visually or even calculated from the cumulative histogram by setting a threshold at the 80% level.

Length[Select[tmp, (#1[[2]] < 0.8 &)]]

We display the edge maps for three selected threshold values.

Show[
        GraphicsArray[(Graphics[Threshold[mag, #1]] &) /@ {% + 6, %, % - 
                6}]];

Based on visual inspection of the edge maps, it seems reasonable to accept the threshold value of  t.

A zero-crossing edge operator was originally proposed by Marr and Hildreth [2]. They suggested that in order to effectively detect intensity changes (edges), the operator needs to have two characteristics. First, it must be a differential operator, taking either a first or second spatial derivative of the image. Second, it should be capable of being tuned to act at any desired scale, so that large filters can be used to detect blurry shadow edges, and small ones can be used to detect sharply focused fine detail in the image. This led to the so-called Laplacian-of-Gaussian edge operator. This is a compound operator, which combines a smoothing operation, using a Gaussian-shaped, linear phase FIR filter, with a differentation operation, using a discrete Laplacian. The edges are identified by the location of zero-crossings (recall that the second derivative changes sign in the vicinity of maxima of the first derivative).

Smoothing and differentiating filters.

LaplacianFilter returns a number of common FIR filter approximations to the Laplacian operator. The default argument value (α = 2) returns the minimum-variance discrete Laplacian subject to the conditions that pixel noise is uncorrelated and has uniform variance [3]. The Laplacian is the lowest order isotropic (i.e., orientation-independent) operator.

MatrixForm  /@  {LaplacianFilter[0. ], LaplacianFilter[0.5], 
        LaplacianFilter[2. ]}

This computes the Laplacian of the example image.

laplacian = EdgeMagnitude[wing, LaplacianFilter[]];

This shows the original and filtered images.

Show[GraphicsArray[Graphics /@ {wing, laplacian}]];

The zero crossings are clearly visible as they occur at the border between bright and dark areas of the Laplacian image. The mid-tones cover areas of small responses to the Laplacian operator, which corresponds to uniform (flat) brightness areas of the original image.

ImageRange[laplacian]

In order to mitigate the increase in pixel noise due to differentiation, the image is filtered with a lowpass filter. This is accomplished by a Gaussian-shaped, linear phase FIR filter. Since convolution is associative and commutative, the two-step sequence can be reduced to one by constructing a compound operator. LoGFilter is a compound operator, whose values are samples of the Laplacian of the bivariate Gaussian function with variance (σ^2).

\(\((g = 
        Exp[\(-\(\(x\^2 + 
                  y\^2\)\/\(2\ σ\^2\)\)\)]; 
Δg\  = \ ∂\_{x, 2}g + ∂\_{y, 2}g\  // Simplify)\)\)

Here we set (σ^2) = 1 and plot both the resulting function as a surface, and the functions profile for y = 0.

Show[
        GraphicsArray[{Plot3D[
              Evaluate[Δg /. {σ -> 1}], {x, -3, 
                3}, {y, -3, 3}, PlotRange -> All, 
              DisplayFunction -> Identity], 
            Plot[Evaluate[Δg /. {σ -> 1, 
                    y -> 0}], {x, -4, 4}, 
              DisplayFunction -> Identity]}]];

Here we filter the example image using LoGFilter.

EdgeMagnitude[wing, LoGFilter[21, 21]];

Shown below are the results of convolving the example image with a Laplacian-of-Gaussian filter followed by zero-crossing detection. The latter is implemented with the function ZeroCrossing.

Show[GraphicsArray[Graphics /@ {%,  ZeroCrossing[%]}]];

[1] Gonzalez, R. C. and Woods, R. E., Digital Image Processing, Addison-Wesley, Reading, MA, 1992.

[2] Marr, D. and Hildreth, E., "Theory of edge detection," Proceedings The Royal Society London, Vol. 207, 1980, 187-217.

[3] Haralick, R. M. and Shapiro, L.G., Computer and Robot Vision, Vol.1, Addison-Wesley, Reading, MA, 1992.