In 5.1, we discussed the topic of edge detection from the point-of-view of derivative filters and implemented edge detectors by treating them as modified gradient operators. In this section we take an alternate view of the same problem, and derive a multiscale edge detection algorithm that draws upon wavelet theory. Mallat and Zhong published some of the pioneering work in this field in a 1992 paper that presented a numerical approach for the characterization of signals in terms of their multiscale edges9. In this section we will show, with some amount of hand-waving, that classical edge detectors can be interpreted as discretized wavelet transforms, where the convolution of the image with an edge detection kernel (e.g., our old friend the Sobel kernel) in fact produces the coefficients of the wavelet transform. Readers interested in a more rigorous treatment of the subject should consult the references. We will improve the performance of classical edge detectors by incorporating scale using a wavelet-based model. MATLAB functions for multiscale edge detection are presented, and a fixed-point edge detection framework suitable for embedded DSP deployment and implemented on the C6416 is also given.
As shown in 5.1, detecting edges boils down to identifying the local irregularities in an image. Perceived edges correspond to those locations where the image undergoes a sharp variation in brightness, which correlates to the maxima of the first derivative of the image. From the multi-level pyramid wavelet decomposition tree, it can be seen that performing the wavelet transform on an image is equivalent to successively smoothing the image and then differentiating it. At each level, the previous level’s approximation coefficients are passed through complementary high-pass filters. The output of these high-pass filters are the LL, HL, and HH subbands, and a quick glance at Figures 6-5 and 6-8 confirms that the detail coefficient matrices can be thought of as edge enhanced images pertaining to the edges in the horizontal, vertical, and diagonal directions, respectively. This is evident from the definition of the QMF bank – the high-pass filters produce the detail coefficients and we showed in 4.1.2 and then again in 5.1 that convolution kernels that accentuate edges are high-pass filters. Smoothing and then differentiating the edges results in peaks, corresponding to bright regions in the detail coefficient images in 6-5b and 6-8b, and these peaks can be subsequently segmented by picking out the local maxima of the absolute value of the derivative images.
The Canny edge detector, also introduced in 5.1, uses the first derivative of a Gaussian as the edge-enhancement filter and Mallat generalized Canny’s method, linking it to the wavelet transform in [11]. The Canny edge detector uses the modulus of the gradient vector V’f(x,yj as a low-pass filtered version of the input image. Points in the image are defined as an edge if the modulus of ¥f(x,y) are locally maximum in a direction parallel to Vf(x,y). Mallat showed that the wavelet transform is proportional to the gradient of a smoothed version of f(x,y), as the detail coefficients are derivatives of the approximation coefficients, which in turn is a low-pass filtered or smoothed version of f(x,y). We can therefore implement edge detection by looking at the modulus (absolute value) of the wavelet detail coefficients, and finding the local maxima of the modulus, or the wavelet modulus maxima. Multiscale edge detection is based on analysis and segmentation of the wavelet modulus maxima.
Mallat’s algorithm makes use of the modulus maxima of the wavelet transform and the angle of the DWT8. In addition, only the horizontal and vertical wavelet detail coefficients are used in the computation, and as a result the entire process can be made more efficient. By not computing the diagonal detail coefficients, we lessen the burden on the algorithm at the cost of perfect reconstruction. If we define WH (x, y, 2J) to be the horizontal detail coefficients at scale j and Wv (x, y, 21) the vertical detail coefficients at scale j, then the wavelet transform modulus M(x, y, 2′) at scale j is
The angle Af(x, y, 2s) at scale j is then
Multiscale edges are points (x’, y’) where M(x’, y’, 21) is larger than its two neighbors in the direction given by A f{x’, y’,2j). In the digital realm, just like the Canny edge detector, A fix’, y’,21) is typically quantized to the set of eight possible directions corresponding to the eight neighbors of the pixel flx’,y’)-
The algorithm we implement on the TI C6416 DSP is based on the Mallat’s algorithm, but uses a simpler segmentation scheme more suitable for an embedded platform. Algorithm 6-1 is a very high-level description of a multiscale edge detection scheme that uses two thresholds to segment a multiscale edge-enhanced image.
Algorithm 6-1: Multiscale Edge Detector
INPUT: Image I, global threshold T, number decomposition levels N OUTPUT: binary edge image
for each wavelet decomposition level j = 1 … N
compute DWT coefficients at level j end
let Dn be the detail coefficients at final level N compute local modulus maxima of DN compute global maxima threshold based on scale N if (local modulus maxima > T) then mark as edge
Of course, with an algorithm description as high-level as Algorithm 6-1, the devil is always in the details. In subsequent sections we elaborate on the specifics and implement an image processing scheme based on Algorithm 61 to edge enhance and subsequently segment images on the C6416 DSP. However, instead of the dyadic 2D DWT derived in 6.1.1, we replace it with a version of the wavelet transform that goes by the moniker "algorithme a trous".
The Undecimated Wavelet Transform
Mallat’s fast multiscale edge detector is implemented using an undecimated form of the wavelet transform known as the "algorithme a trous"12. This fast dyadic wavelet transform is similar to that depicted in Figure 6-2, except no downsampling is performed. As a consequence, the transform is said to be redundant, or overcomplete, in the sense that superfluous coefficients are retained and successive coefficient matrices are the same size as the input. Recall that with the traditional dyadic DWT, the output from each filter is one-half the length of the input, due to the decimation step. With each decomposition level in the algorithme a trous, the low-pass wavelet filter hLP and high-pass wavelet filter hHP are "stretched" in a sense, and the convolution is performed without any subsampling. This rescaling of the wavelet filters is achieved by inserting 2J-1 zeros between every filter coefficient (where j is the current decomposition level), hence the name algorithme a trous which means the "holes algorithm" in French. The analysis filter bank for the a trous algorithm is shown in Figure 6-9.
The a trous algorithm outputs N+l images of the same size, where N is the number of decomposition levels – a single approximation image plus N detail images. This algorithm is very attractive for a number of reasons: it is translation-invariant (a shift in the input simply shifts the coefficients), the discrete transform values are known exactly at every pixel location without the need of any interpolation, and correlation across scales is easily exploited due to its inherent structure. Moreover, an even simpler form of the algorithm exists: through careful choice of the filter hw we calculate the detail coefficients as differences between consecutive approximation coefficients13. This scheme is illustrated graphically in Figure 6-10, where Aj and Dj denote the approximation and detail coefficient matrices, respectively. With this formulation, synthesis is merely the sum of all detail coefficient matrices and the final approximation coefficient matrix. For an N-level decomposition, the reconstruction formula is then
Figure 6-11. Analysis filter bank for the a trous algorithm. With each level, the filters expand by inserting zeros between every coefficient. In contrast to the left-hand side of Figure 6-2, note the lack of decimators.
Figure 6-12. Simplified a trous algorithm. Approximation coefficient matrices at level j are denoted by Aj and detail coefficient matrices at level j are denoted by Dj.
Edge Detection with the Undecimated Wavelet Transform
Thus far we have introduced the basics behind multiscale edge detection and have explained the mechanics behind the undecimated wavelet transform. In this section we put the two together and derive a specific edge detection algorithm modeled on the basis of Algorithm 6-1 that can be parameterized by the number of wavelet decompositions.
Recall that there are many choices for wavelet basis functions, and the selection of which basis to use for edge detection is quite important. In particular, the mother wavelet should be symmetric. Symmetry is important in this context because we are essentially differentiating a smoothed image, and therefore a lack of symmetry implies that the location of edges will change as we successively smooth and differentiate the image. With the appropriate choice of wavelet, the locations of edges correspond to the modulus maxima of the wavelet transform at a given scale.
The cubic B-spline (B3) wavelet fits the bill perfectly – it is symmetric and easily derived. A spline of degree 3 is the convolution of four box functions, and in general the spline @N-i(t) of degree N-l is the convolution of N box functions1. Hence to generate the B3 lowpass filter coefficients we convolve (V2, V2) with itself four times:
To generate a low-pass filter kernel hLF suitable for use with the simplified a trous algorithm we extend the above to two dimensions as so:
This kernel is attractive because we can replace the expensive division operations by left bit-shifts since the denominators are all powers of two. We now have enough to prototype the undecimated wavelet transform in MATLAB and see first-hand how it can be used to implement an edge detector. Listing 6-4 is a MATLAB implementation of the simplified a trous algorithm, a trous dwt, which does not utilize any functions from the Wavelet Toolbox. The function requires the caller to provide the image and number of decomposition levels, and returns the final approximation image and a three-dimensional array where each "plane" consists of the detail image at that particular scale. We have not encountered higher-dimensional arrays in MATLAB until now, but they work just as one might expect. For example, if one were to perform a three-level decomposition of an image I using a_trous_dwt, accessing the second detail coefficient matrix is accomplished as:
Listing 6-4: MATLAB function for the simplified a trous algorithm.
function [A, D] = a_trous_dwt(I, N)
The simplicity of the reconstruction formula coupled with MATLAB’s inherent support of higher dimensional matrices leads to a very simple implementation of the inverse a trous wavelet transform, given in Listing 6-5.
Listing 6-5: a_trous_idwt .m, a MATLAB function for the inverse simplified a trous algorithm.
Now that we have readied an undecimated DWT, we can proceed to graft an edge detector on top of it. As in the case of the classical edge detectors, a thresholding strategy is needed, for even though the wavelet modulus maxima is a decent approximation to the edges of an image, there will exist numerous false positives even in the absence of noise for any reasonably complex image. There are numerous thresholding schemes, ranging from a simple fixed global threshold to more sophisticated multiscale correlation thresholds, as described in [14]. The algorithm we port onto the DSP segments the final detail coefficient image using the two thresholds mentioned in Algorithm 6-1. An estimate of the wavelet modulus maxima that does not factor into account the angle of the DWT is first computed by comparing the wavelet modulus to the local average. This approximation to the wavelet modulus maxima is then compared to a global threshold dynamically calculated from the coefficients of the estimated modulus of the detail coefficients. Listing 6-6 is the MATLAB function that implements this scheme, and forms the basis of the fixed-point C code that we will port to the DSP platform in the next section.
Listing 6-6: edgedet. m, a multiscale edge detector.
The variable D is the absolute value of the highest level detail coefficient matrix and the subject of the segmentation process that produces a binary image. This makes intuitive sense, as the most negative and most positive coefficients correspond to the stronger edges, and coefficients in this high-pass filtered image close to zero are locations where the approximation image is of constant brightness. The majority of the edge detection work is performed in the third statement of edgedet.m, and warrants a detailed explanation. The call to filter2 (ones (3) /9, D) returns an image of the same size as D where pixel values are the local average of the surrounding 3×3 neighborhood. The left-hand side of the statement, (D>filter2 (ones (3) /9, D), returns a logical (boolean) coefficient matrix that will have a value of 1 wherever the modulus of the detail coefficient matrix is larger than the average of its neighbors, and 0 otherwise. Hence, we approximate the wavelet modulus maxima without the need for computing an expensive arctangent as Mallat’s algorithm does. The right-hand side of the final statement, (D>mean2 (D) ), returns another logical matrix with Is wherever the modulus of the detail coefficient matrix is larger than the mean of the entire detail image. This is the global threshold mentioned in Algorithm 6-1 and here we dynamically calculate it. The result of multiplying these two matrices together in an element-by-element fashion (not matrix multiplication) is in effect a logical "and" operation that yields the final binary segmented image.
What makes this algorithm so powerful is that we now have the ability to control the behavior of the edge detection scheme by incorporating scale into the overall framework. With classical edge detectors there is no concept of scale, but with a wavelet model we are free to adjust the scale as desired by controlling the number of wavelet decompositions. Scale controls the significance of the detected edges, as illustrated in Figure 6-13. A careful balance must be struck between the image resolution and the scale of edges: high resolution and small scale (reduced number of wavelet decompositions) results in comparatively noisier and more discontinuous edges, as can be seen from Figure 6-13 when the number of wavelet decompositions is 1. The "peppers" image, in particular, suffers from this problem. However, low resolution coupled with large scale may result in undetected edges.
Sharper edges are more likely to be retained by the wavelet transform across subsequent scales, whereas less significant edges are attenuated as the scale increases. Again referring to Figure 6-13, as the number of wavelet decompositions increases, the boundaries between the objects in the peppers image and the structure present in the butterfly’s wings become more pronounced, due to the coarser resolution which has the effect of smoothing out the irregularities in the image. This comes at a price however, because even though the noise and details within the image decrease, the exact location of the edges now becomes more difficult to discern.
State-of-the-art multiscale segmentation schemes are more sophisticated than that employed in edgedet .m. They incorporate scale into the criteria for discriminating between real structure and noise. In general, wavelet transform coefficients are likely to increase across scale wherever there are significant edges and likewise, the coefficients are equally likely to decrease across scale at location of noise and details. In fact, we will leverage this behavior during the development of a wavelet-based denoising algorithm in the 6.3.
Figure 6-13. Wavelet-based edge detection using edgedet .m. (a) "Peppers" image, with segmented binary edge images using 1, 3, and 6 decomposition levels, (b) "Monarch" image, again with 1, 3, and 6 decomposition levels.
