Canny边界检测实现

 

Canny's Edge Detector: Implementation

Canny's Edge Detector is well known for its ability to generate single-picel thick continuous edges. It consists of four major steps, which are described below, along with interesting implementation details and outputs. The program has four inputs: Input image I, value of smoothing parameter sigma, high threshold Th and low threshold Tl.

void canny(Image& img, Image& gradientMap, Image& gradientStrength, Image& gradientOrientation, bool useGauss, float sigma, float lowerThreshold, float upperThreshold, Image& result){
  if (useGauss){
   gaussian(img, sigma);    // step 1
   sobelGradientMap(img, gradientMap);
  } else{
   sobelGradientMap(img, gradientMap);
  }
  gradientMapToGradientStrength(gradientMap, gradientStrength);   // step 2
  gradientMapToGradientOrientation(gradientMap, gradientOrientation);
  nonMaximaSurpression(gradientStrength, gradientOrientation);     // step 3
  hysteresisTheshold(gradientStrength, gradientOrientation, lowerThreshold, upperThreshold, result);   //step 4

}

 

Step 1: Generation of Masks:

This module requires the value of sigma as an input and generates x- and y-derivative masks as output. In Canny's method, the masks used are 1st derivative of a Gaussian in x- and y-directions.The masks are also written in two text files for reference.

To generate the masks, the first step is a reasonable computation of the mask size. The mask size should not be too large compared to the lobes of the mask, otherwise it will result in unnecessary computational overhead during convolution. At the same time, the mask size should not be so small to loose the characteristics of primary lobes. We chose mask size based on analysing the Gaussian and applying a threshold T. This idea is explained in the figure. T is a real-number between 0 and 1. Width of the mask is derived based on T as follows:
First, size of half mask, sHalf is computed by finding the point on the curve where the Gaussian value drops below T, i.e. solving exp(-x^2/(2*sigma^2) = T. This gives sHalf = round(sqrt(-log(T) * 2 * sigma^2)). The mask size is then 2*sHalf + 1 to incorporate both positive and negative sides of the mask. We also put a lower limit on sigma of 0.5, because below that the mask size came out to be less than 3, which was not reasonable for finding fx and fy. The values of mask size obtained by this method for various sigma values are shown below:

Sigma	Size of Mask
0.5	3x3
1	5x5
2	9x9
3	13x13
4	19x19

In our notation, x-direction is taken to be vertically downward, and y-direction horizontally rightward. This means that the lobes of x-mask are vertically aligned, while those of y-mask are horizontally aligned. Once the masks are generated, they are scaled and rounded off so that convolution is done with integer values rather than floats. The scale factor is saved, because gradient magnitude is later scaled down (after convolution) by the same factor. Scaled masks for sigma = 1 are shown below.

x-mask 15 69 114 69 15 35 155 255 155 35 0 0 0 0 0 -35 -155 -255 -155 -35 -15 -69 -114 -69 -15

y-mask 15 35 0 -35 -15 69 155 0 155 69 114 255 0 -255 -114 69 155 0 155 69 15 35 0 35 15

Step 2: Applying Masks to Images

The masks are applied to the images using convolution. The result is then scaled down by the same factor which was used to scale up the masks. To write output to image files, the min and max values are scaled to 0 and 255 respectively. Results on a sample image are shown below.

sigma = 0.5, fx and fy

sigma = 2, fx and fy

The effect of increasing sigma is obvious from the convolved images: the gradients are much smoother. It should be noted that horizontally running edges are identified in fx, because the gradient change occurs when moving along the x (downward) direction. Simlarly vertical edges are seen in fy because our definition of y-axis is along the column direction.

Gradient magnitude, M, is computed from fx and fy images, using the magnitude formula:

M for sigma = 0.5

M for sigma = 2

Again, the effect of increasing sigma is dramatically highlighted in these two examples.

Gradient Direction, Phi, is computed using atan2 function. The difference between atan and atan2 functions is that atan returns output range from -pi/2 to pi/2, whereas atan2 returns in the range of -pi to pi. Thus atan2 is preferred because that is the real range of the possible directions of gradient in an image. To keep things simple in our code, we converted the angle returned by atan2 function to degrees and added 180 to get an output range of 0-360 degrees.

Visualizing gradient direction in a meaningful fashion is difficult. Raw gradient direction, if mapped to the range of 0-255 looks as follows:

This image is difficult to interpret, even though it can be seen that linear structures like pillars, have similar gradient direction values. Also, the ropes from which the flags are hanging show similar values of gradient direction. A better visualization technique is to view only those values of direction for which gradient magnitude is high enough. In this method, we made all values in the gradient direction image -1 where gradient magnitude was below a certain threshold (t = 10). Thus, the output is visible only at significant edges. The image was then false colored to get the following:

This output is much more interesting than the previous one. Edges with similar direction are colored similarly. Each directions has two colors associated with it (because 180 apart edges run in the same direction). Thus the string of the flags, running at approximately -45 degrees has yellow at one side and blue at the other.

The effect of this method of visualizing gradient direction can be used on a circle image to show edges in different directions have different values of Phi. For a white circle on a black background, the gradient direction looks like this:
 
This clearly indicates that Phi is different for edges running in different directions.

 

Step 3: Non Maxima Suppression:

Non maxima supression step makes all edges in M one pixel thick. This is an important step in Canny's algorithm, which distinguishes it from other algorithms. The first step is to quantize gradient direction into just four directions. In our implementation, the following values were assigned during quantization

Angle	Value assigned
0 to 22.5 157.5 to 202.5 337.5 to 360	0
22.5 to 67.5 202.5 to 247.5	1
67.5 to 112.5 247.5 to 292.5	2
112.5 to 157.5 292.5 to 337.5	3

This assignment follows simply from observing that out of the 8 discretized directions, 4 opposing pairs are unique in terms of picking neighborhood pairs in a straight line around each edge point.

Quantized direction for the mecca06 image, visualized using similar color coding as above is shown below:

Here, dark blue shows no edge, and four other colors show four quantized directions. It is clear from this image that quantization has worked properly.

The next step is to pick two neighbors of each edge point along the gradient direction. This is because gradient direction is perpendicular to the edge, and therefore, this is the direction in which we are searching for edge points. For example, for quantized direction 0 (in the table above), the gradient direction could be less than 0 degrees, meaning the edge is a horizontal edge. Therefore the two neighbors that need to be picked for comparison are the north and south neighbors, denoted in code by (r-1, c) and (r+1,c) respectively. If the edge point (r,c) is greater than both these neighbors, then it is maintained in M otherwise it is made zero.

Output before and after nms for sigma = 2 is shown below:

Notice how each edge is 1 pixel thick now, but has different values of gradient magnitude (different gray level in the image) based on its strength.

 

for (unsigned y = 0; y < height; ++y){
    for (unsigned x = 0; x < width; ++x){

   // 方向在<0, PI>，区域0: <0, PI/8>，区域1:<PI/8, 3*PI/8>，区域2:<3*PI/8, 5*PI/8>，区域3:<5*PI/8, 7*PI/8>，

   //区域4:<7*PI/8，PI>
     short direction = (*orientation < 0) ? (short)((* orientation + PI) / PI_4 + 0.5) : (short)(*orientation / PI_4 + 0.5);
     switch(direction){
     case 0:
     case 4:
      strengthBefore = (x != 0) ? *(strength - 1) : 0;
      strengthAfter = (x != lastW) ? *(strength + 1) : 0;
      break;
     case 1:
      strengthBefore = (x != 0 && y != 0) ? *(strength - width - 1) : 0;
      strengthAfter = (x != lastW && y != lastH) ? *(strength + width + 1) : 0;
      break;
     case 2:
      strengthBefore = (y != 0) ? *(strength - width) : 0;
      strengthAfter = (y != lastH) ? *(strength + width) : 0;
      break;
     case 3:
      strengthBefore = (x != 0 && y != lastH) ? *(strength + width - 1) : 0;
      strengthAfter = (x != lastW && y != 0) ? *(strength - width + 1) : 0;
      break;
     }
     *resultPtr = (*strength < strengthBefore || *strength <= strengthAfter) ? 0 : *strength;
     ++resultPtr;
     ++strength;
     ++orientation;
    }
   }

 

Step 4: Hysteresis Thresholding

The final step in Canny's edge detection algorithm is to apply two thresholds to follow edges. Since edges are to be followed recursively by looking at neighbors, we first made the border pixels zero, so that finding neighbors does not go out of bounds of the image. Next, the image is scanned from left to right, top to bottom. The first pixel in non-maxima suppressed magnitude image which is above a certain threshold, Th, is declared an edge. Then all its neighbors are recursively followed, and those above threshold, Tl, are marked as an edge. A visited map is also maintained so that recursion does not loop infinitely. Thus there are really two stopping conditions: if a neighbor is below Tl, we wont recurse on it; also, if a neighbor has already been visited, then we won't recurse on it.

The output can be animated to highlight the effect of edge-following algorithm. Click here to see an mpeg movie. In this video, disconnected strings of edges are shown in the order they were encountered during the left-to-right-top-to-bottom scan. This movie was generated using Th = 50, Tl = 10. If Th is increased but Tl remains the same, lesser edge components will be detected, but their lengths will be the same. If Tl is increased but Th is the same, then the same number of edge components will be detected but their lengths will be shorter. This is supported by the following examples:

Th = 50, Tl = 10	Th = 100, Tl = 10 (incr. in Th only)	Th = 50, Tl = 40 (incr. in Tl only)

Results on other images:

Original Image	Parameters	Edge Map	Non Maxima Suppressed output
	Th = 80  Tl = 30  sigma = 1
	Th = 80 Tl = 30 Sigma = 1

In the 2nd image, note that it is very difficult to obtain the complete closed edge of the egg, because the edge in the middle portion is very weak (shown by the NMS output).