Canny边缘检测算子是澳洲计算机科学家约翰·坎尼（John F. Canny）于1986年开发出来的一个多级边缘检测算法。
算法步骤
1. 降噪
任何边缘检测算法都不可能在未经处理的原始数据上很好地处理，所以第一步是对原始数据与高斯平滑模板作卷积，得到的图像与原始图像相比有些轻微的模糊（blurred）。这样，单独的一个像素噪声在经过高斯平滑的图像上变得几乎没有影响。
2. 寻找图像中的亮度梯度
图像中的边缘可能会指向不同的方向，所以Canny算法使用4个mask检测水平、垂直以及对角线方向的边缘。原始图像与每个mask所作的卷积都存储起来。对于每个点我们都标识在这个点上的最大值以及生成的边缘的方向。这样我们就从原始图像生成了图像中每个点亮度梯度图以及亮度梯度的方向。
3. 在图像中跟踪边缘
较高的亮度梯度比较有可能是边缘，但是没有一个确切的值来限定多大的亮度梯度是边缘多大又不是，所以Canny使用了滞后阈值。
滞后阈值需要两个阈值——高阈值与低阈值。假设图像中的重要边缘都是连续的曲线，这样我们就可以跟踪给定曲线中模糊的部分，并且避免将没有组成曲线的噪声像素当成边缘。所以我们从一个较大的阈值开始，这将标识出我们比较确信的真实边缘，使用前面导出的方向信息，我们从这些真正的边缘开始在图像中跟踪整个的边缘。在跟踪的时候，我们使用一个较小的阈值，这样就可以跟踪曲线的模糊部分直到我们回到起点。
一旦这个过程完成，我们就得到了一个二值图像，每点表示是否是一个边缘点。
一个获得亚像素精度边缘的改进实现是在梯度方向检测二阶方向导数的过零点
它在梯度方向的三阶方向导数满足符号条件
其中Lx，Ly，··· ，Lyyy表示用高斯核平滑原始图像得到的尺度空间表示{\displaystyle L}L计算得到的偏导数。用这种方法得到的边缘片断是连续曲线，这样就不需要另外的边缘跟踪改进。滞后阈值也可以用于亚像素边缘检测。
代码实现
C
该程序读取每像素8bit的灰度BMP文件，并将结果保存到`out.bmp’。用-lm编译。
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <string.h>
#include <stdbool.h>
#include <assert.h>#define MAX_BRIGHTNESS 255// C99 doesn't define M_PI (GNU-C99 does)
#define M_PI 3.14159265358979323846264338327/** Loading part taken from* http://www.vbforums.com/showthread.php?t=261522* BMP info:* http://en.wikipedia.org/wiki/BMP_file_format** Note: the magic number has been removed from the bmpfile_header_t* structure since it causes alignment problems* bmpfile_magic_t should be written/read first* followed by the* bmpfile_header_t* [this avoids compiler-specific alignment pragmas etc.]*/typedef struct {
uint8_t magic[2];
} bmpfile_magic_t;typedef struct {
uint32_t filesz;uint16_t creator1;uint16_t creator2;uint32_t bmp_offset;
} bmpfile_header_t;typedef struct {
uint32_t header_sz;int32_t width;int32_t height;uint16_t nplanes;uint16_t bitspp;uint32_t compress_type;uint32_t bmp_bytesz;int32_t hres;int32_t vres;uint32_t ncolors;uint32_t nimpcolors;
} bitmap_info_header_t;typedef struct {
uint8_t r;uint8_t g;uint8_t b;uint8_t nothing;
} rgb_t;// Use short int instead `unsigned char' so that we can
// store negative values.
typedef short int pixel_t;pixel_t *load_bmp(const char *filename,bitmap_info_header_t *bitmapInfoHeader)
{
FILE *filePtr = fopen(filename, "rb");if (filePtr == NULL) {
perror("fopen()");return NULL;}bmpfile_magic_t mag;if (fread(&mag, sizeof(bmpfile_magic_t), 1, filePtr) != 1) {
fclose(filePtr);return NULL;}// verify that this is a bmp file by check bitmap id// warning: dereferencing type-punned pointer will break// strict-aliasing rules [-Wstrict-aliasing]if (*((uint16_t*)mag.magic) != 0x4D42) {
fprintf(stderr, "Not a BMP file: magic=%c%c\n",mag.magic[0], mag.magic[1]);fclose(filePtr);return NULL;}bmpfile_header_t bitmapFileHeader; // our bitmap file header// read the bitmap file headerif (fread(&bitmapFileHeader, sizeof(bmpfile_header_t),1, filePtr) != 1) {
fclose(filePtr);return NULL;}// read the bitmap info headerif (fread(bitmapInfoHeader, sizeof(bitmap_info_header_t),1, filePtr) != 1) {
fclose(filePtr);return NULL;}if (bitmapInfoHeader->compress_type != 0)fprintf(stderr, "Warning, compression is not supported.\n");// move file point to the beginning of bitmap dataif (fseek(filePtr, bitmapFileHeader.bmp_offset, SEEK_SET)) {
fclose(filePtr);return NULL;}// allocate enough memory for the bitmap image datapixel_t *bitmapImage = malloc(bitmapInfoHeader->bmp_bytesz *sizeof(pixel_t));// verify memory allocationif (bitmapImage == NULL) {
fclose(filePtr);return NULL;}// read in the bitmap image datasize_t pad, count=0;unsigned char c;pad = 4*ceil(bitmapInfoHeader->bitspp*bitmapInfoHeader->width/32.) - bitmapInfoHeader->width;for(size_t i=0; i<bitmapInfoHeader->height; i++){
for(size_t j=0; j<bitmapInfoHeader->width; j++){
if (fread(&c, sizeof(unsigned char), 1, filePtr) != 1) {
fclose(filePtr);return NULL;}bitmapImage[count++] = (pixel_t) c;}fseek(filePtr, pad, SEEK_CUR);}// If we were using unsigned char as pixel_t, then:// fread(bitmapImage, 1, bitmapInfoHeader->bmp_bytesz, filePtr);// close file and return bitmap image datafclose(filePtr);return bitmapImage;
}// Return: true on error.
bool save_bmp(const char *filename, const bitmap_info_header_t *bmp_ih,const pixel_t *data)
{
FILE* filePtr = fopen(filename, "wb");if (filePtr == NULL)return true;bmpfile_magic_t mag = {
{
0x42, 0x4d}};if (fwrite(&mag, sizeof(bmpfile_magic_t), 1, filePtr) != 1) {
fclose(filePtr);return true;}const uint32_t offset = sizeof(bmpfile_magic_t) +sizeof(bmpfile_header_t) +sizeof(bitmap_info_header_t) +((1U << bmp_ih->bitspp) * 4);const bmpfile_header_t bmp_fh = {
.filesz = offset + bmp_ih->bmp_bytesz,.creator1 = 0,.creator2 = 0,.bmp_offset = offset};if (fwrite(&bmp_fh, sizeof(bmpfile_header_t), 1, filePtr) != 1) {
fclose(filePtr);return true;}if (fwrite(bmp_ih, sizeof(bitmap_info_header_t), 1, filePtr) != 1) {
fclose(filePtr);return true;}// Palettefor (size_t i = 0; i < (1U << bmp_ih->bitspp); i++) {
const rgb_t color = {
(uint8_t)i, (uint8_t)i, (uint8_t)i};if (fwrite(&color, sizeof(rgb_t), 1, filePtr) != 1) {
fclose(filePtr);return true;}}// We use int instead of uchar, so we can't write img// in 1 call any more.// fwrite(data, 1, bmp_ih->bmp_bytesz, filePtr);// Padding: http://en.wikipedia.org/wiki/BMP_file_format#Pixel_storagesize_t pad = 4*ceil(bmp_ih->bitspp*bmp_ih->width/32.) - bmp_ih->width;unsigned char c;for(size_t i=0; i < bmp_ih->height; i++) {
for(size_t j=0; j < bmp_ih->width; j++) {
c = (unsigned char) data[j + bmp_ih->width*i];if (fwrite(&c, sizeof(char), 1, filePtr) != 1) {
fclose(filePtr);return true;}}c = 0;for(size_t j=0; j<pad; j++)if (fwrite(&c, sizeof(char), 1, filePtr) != 1) {
fclose(filePtr);return true;}}fclose(filePtr);return false;
}// if normalize is true, map pixels to range 0..MAX_BRIGHTNESS
void convolution(const pixel_t *in, pixel_t *out, const float *kernel,const int nx, const int ny, const int kn,const bool normalize)
{
assert(kn % 2 == 1);assert(nx > kn && ny > kn);const int khalf = kn / 2;float min = FLT_MAX, max = -FLT_MAX;if (normalize)for (int m = khalf; m < nx - khalf; m++)for (int n = khalf; n < ny - khalf; n++) {
float pixel = 0.0;size_t c = 0;for (int j = -khalf; j <= khalf; j++)for (int i = -khalf; i <= khalf; i++) {
pixel += in[(n - j) * nx + m - i] * kernel[c];c++;}if (pixel < min)min = pixel;if (pixel > max)max = pixel;}for (int m = khalf; m < nx - khalf; m++)for (int n = khalf; n < ny - khalf; n++) {
float pixel = 0.0;size_t c = 0;for (int j = -khalf; j <= khalf; j++)for (int i = -khalf; i <= khalf; i++) {
pixel += in[(n - j) * nx + m - i] * kernel[c];c++;}if (normalize)pixel = MAX_BRIGHTNESS * (pixel - min) / (max - min);out[n * nx + m] = (pixel_t)pixel;}
}/** gaussianFilter:* http://www.songho.ca/dsp/cannyedge/cannyedge.html* determine size of kernel (odd #)* 0.0 <= sigma < 0.5 : 3* 0.5 <= sigma < 1.0 : 5* 1.0 <= sigma < 1.5 : 7* 1.5 <= sigma < 2.0 : 9* 2.0 <= sigma < 2.5 : 11* 2.5 <= sigma < 3.0 : 13 ...* kernelSize = 2 * int(2*sigma) + 3;*/
void gaussian_filter(const pixel_t *in, pixel_t *out,const int nx, const int ny, const float sigma)
{
const int n = 2 * (int)(2 * sigma) + 3;const float mean = (float)floor(n / 2.0);float kernel[n * n]; // variable length arrayfprintf(stderr, "gaussian_filter: kernel size %d, sigma=%g\n",n, sigma);size_t c = 0;for (int i = 0; i < n; i++)for (int j = 0; j < n; j++) {
kernel[c] = exp(-0.5 * (pow((i - mean) / sigma, 2.0) +pow((j - mean) / sigma, 2.0)))/ (2 * M_PI * sigma * sigma);c++;}convolution(in, out, kernel, nx, ny, n, true);
}/** Links:* http://en.wikipedia.org/wiki/Canny_edge_detector* http://www.tomgibara.com/computer-vision/CannyEdgeDetector.java* http://fourier.eng.hmc.edu/e161/lectures/canny/node1.html* http://www.songho.ca/dsp/cannyedge/cannyedge.html** Note: T1 and T2 are lower and upper thresholds.*/
pixel_t *canny_edge_detection(const pixel_t *in,const bitmap_info_header_t *bmp_ih,const int tmin, const int tmax,const float sigma)
{
const int nx = bmp_ih->width;const int ny = bmp_ih->height;pixel_t *G = calloc(nx * ny * sizeof(pixel_t), 1);pixel_t *after_Gx = calloc(nx * ny * sizeof(pixel_t), 1);pixel_t *after_Gy = calloc(nx * ny * sizeof(pixel_t), 1);pixel_t *nms = calloc(nx * ny * sizeof(pixel_t), 1);pixel_t *out = malloc(bmp_ih->bmp_bytesz * sizeof(pixel_t));if (G == NULL || after_Gx == NULL || after_Gy == NULL ||nms == NULL || out == NULL) {
fprintf(stderr, "canny_edge_detection:"" Failed memory allocation(s).\n");exit(1);}gaussian_filter(in, out, nx, ny, sigma);const float Gx[] = {
-1, 0, 1,-2, 0, 2,-1, 0, 1};convolution(out, after_Gx, Gx, nx, ny, 3, false);const float Gy[] = {
1, 2, 1,0, 0, 0,-1,-2,-1};convolution(out, after_Gy, Gy, nx, ny, 3, false);for (int i = 1; i < nx - 1; i++)for (int j = 1; j < ny - 1; j++) {
const int c = i + nx * j;// G[c] = abs(after_Gx[c]) + abs(after_Gy[c]);G[c] = (pixel_t)hypot(after_Gx[c], after_Gy[c]);}// Non-maximum suppression, straightforward implementation.for (int i = 1; i < nx - 1; i++)for (int j = 1; j < ny - 1; j++) {
const int c = i + nx * j;const int nn = c - nx;const int ss = c + nx;const int ww = c + 1;const int ee = c - 1;const int nw = nn + 1;const int ne = nn - 1;const int sw = ss + 1;const int se = ss - 1;const float dir = (float)(fmod(atan2(after_Gy[c],after_Gx[c]) + M_PI,M_PI) / M_PI) * 8;if (((dir <= 1 || dir > 7) && G[c] > G[ee] &&G[c] > G[ww]) || // 0 deg((dir > 1 && dir <= 3) && G[c] > G[nw] &&G[c] > G[se]) || // 45 deg((dir > 3 && dir <= 5) && G[c] > G[nn] &&G[c] > G[ss]) || // 90 deg((dir > 5 && dir <= 7) && G[c] > G[ne] &&G[c] > G[sw])) // 135 degnms[c] = G[c];elsenms[c] = 0;}// Reuse array// used as a stack. nx*ny/2 elements should be enough.int *edges = (int*) after_Gy;memset(out, 0, sizeof(pixel_t) * nx * ny);memset(edges, 0, sizeof(pixel_t) * nx * ny);// Tracing edges with hysteresis . Non-recursive implementation.size_t c = 1;for (int j = 1; j < ny - 1; j++)for (int i = 1; i < nx - 1; i++) {
if (nms[c] >= tmax && out[c] == 0) {
// trace edgesout[c] = MAX_BRIGHTNESS;int nedges = 1;edges[0] = c;do {
nedges--;const int t = edges[nedges];int nbs[8]; // neighboursnbs[0] = t - nx; // nnnbs[1] = t + nx; // ssnbs[2] = t + 1; // wwnbs[3] = t - 1; // eenbs[4] = nbs[0] + 1; // nwnbs[5] = nbs[0] - 1; // nenbs[6] = nbs[1] + 1; // swnbs[7] = nbs[1] - 1; // sefor (int k = 0; k < 8; k++)if (nms[nbs[k]] >= tmin && out[nbs[k]] == 0) {
out[nbs[k]] = MAX_BRIGHTNESS;edges[nedges] = nbs[k];nedges++;}} while (nedges > 0);}c++;}free(after_Gx);free(after_Gy);free(G);free(nms);return out;
}int main(const int argc, const char ** const argv)
{
if (argc < 2) {
printf("Usage: %s image.bmp\n", argv[0]);return 1;}static bitmap_info_header_t ih;const pixel_t *in_bitmap_data = load_bmp(argv[1], &ih);if (in_bitmap_data == NULL) {
fprintf(stderr, "main: BMP image not loaded.\n");return 1;}printf("Info: %d x %d x %d\n", ih.width, ih.height, ih.bitspp);const pixel_t *out_bitmap_data =canny_edge_detection(in_bitmap_data, &ih, 45, 50, 1.0f);if (out_bitmap_data == NULL) {
fprintf(stderr, "main: failed canny_edge_detection.\n");return 1;}if (save_bmp("out.bmp", &ih, out_bitmap_data)) {
fprintf(stderr, "main: BMP image not saved.\n");return 1;}free((pixel_t*)in_bitmap_data);free((pixel_t*)out_bitmap_data);return 0;
}
Go
Library:Imger
Imger是一个包含了Go编写的图像处理算法集合的库。
注意，在Linux上，示例图像文件名的扩展名需要从.PNG更改为.png，以便库能够识别它。
package mainimport (ed "github.com/Ernyoke/Imger/edgedetection""github.com/Ernyoke/Imger/imgio""log"
)func main() {
img, err := imgio.ImreadRGBA("Valve_original_(1).png")if err != nil {
log.Fatal("Could not read image", err)}cny, err := ed.CannyRGBA(img, 15, 45, 5)if err != nil {
log.Fatal("Could not perform Canny Edge detection")}err = imgio.Imwrite(cny, "Valve_canny_(1).png")if err != nil {
log.Fatal("Could not write Canny image to disk")}
}
PHP
// input: r,g,b in range 0..255
function RGBtoHSV($r, $g, $b) {
$r = $r/255.; // convert to range 0..1$g = $g/255.;$b = $b/255.;$cols = array("r" => $r, "g" => $g, "b" => $b);asort($cols, SORT_NUMERIC);$min = key(array_slice($cols, 1)); // "r", "g" or "b"$max = key(array_slice($cols, -1)); // "r", "g" or "b"// hueif($cols[$min] == $cols[$max]) {
$h = 0;} else {
if($max == "r") {
$h = 60. * ( 0 + ( ($cols["g"]-$cols["b"]) / ($cols[$max]-$cols[$min]) ) );} elseif ($max == "g") {
$h = 60. * ( 2 + ( ($cols["b"]-$cols["r"]) / ($cols[$max]-$cols[$min]) ) );} elseif ($max == "b") {
$h = 60. * ( 4 + ( ($cols["r"]-$cols["g"]) / ($cols[$max]-$cols[$min]) ) );}if($h < 0) {
$h += 360;}}// saturationif($cols[$max] == 0) {
$s = 0;} else {
$s = ( ($cols[$max]-$cols[$min])/$cols[$max] );$s = $s * 255;}// lightness$v = $cols[$max];$v = $v * 255;return(array($h, $s, $v));
}$filename = "image.png";
$dimensions = getimagesize($filename);
$w = $dimensions[0]; // width
$h = $dimensions[1]; // height$im = imagecreatefrompng($filename);for($hi=0; $hi < $h; $hi++) {
for($wi=0; $wi < $w; $wi++) {
$rgb = imagecolorat($im, $wi, $hi);$r = ($rgb >> 16) & 0xFF;$g = ($rgb >> 8) & 0xFF;$b = $rgb & 0xFF;$hsv = RGBtoHSV($r, $g, $b);// compare pixel below with current pixel$brgb = imagecolorat($im, $wi, $hi+1);$br = ($brgb >> 16) & 0xFF;$bg = ($brgb >> 8) & 0xFF;$bb = $brgb & 0xFF;$bhsv = RGBtoHSV($br, $bg, $bb);// if difference in hue > 20, edge is detectedif($hsv[2]-$bhsv[2] > 20) {
imagesetpixel($im, $wi, $hi, imagecolorallocate($im, 255, 0, 0));} else {
imagesetpixel($im, $wi, $hi, imagecolorallocate($im, 0, 0, 0));}}}header('Content-Type: image/jpeg');
imagepng($im);
imagedestroy($im);
Python
在Python中，Canny边缘检测通常将使用scikit-image或OpenCV-Python进行。这是使用numpy / scipy的方法：
#!/bin/python
import numpy as np
from scipy.ndimage.filters import convolve, gaussian_filter
from scipy.misc import imread, imshowdef CannyEdgeDetector(im, blur = 1, highThreshold = 91, lowThreshold = 31):im = np.array(im, dtype=float) #Convert to float to prevent clipping values#Gaussian blur to reduce noiseim2 = gaussian_filter(im, blur)#Use sobel filters to get horizontal and vertical gradientsim3h = convolve(im2,[[-1,0,1],[-2,0,2],[-1,0,1]]) im3v = convolve(im2,[[1,2,1],[0,0,0],[-1,-2,-1]])#Get gradient and directiongrad = np.power(np.power(im3h, 2.0) + np.power(im3v, 2.0), 0.5)theta = np.arctan2(im3v, im3h)thetaQ = (np.round(theta * (5.0 / np.pi)) + 5) % 5 #Quantize direction#Non-maximum suppressiongradSup = grad.copy()for r in range(im.shape[0]):for c in range(im.shape[1]):#Suppress pixels at the image edgeif r == 0 or r == im.shape[0]-1 or c == 0 or c == im.shape[1] - 1:gradSup[r, c] = 0continuetq = thetaQ[r, c] % 4if tq == 0: #0 is E-W (horizontal)if grad[r, c] <= grad[r, c-1] or grad[r, c] <= grad[r, c+1]:gradSup[r, c] = 0if tq == 1: #1 is NE-SWif grad[r, c] <= grad[r-1, c+1] or grad[r, c] <= grad[r+1, c-1]:gradSup[r, c] = 0if tq == 2: #2 is N-S (vertical)if grad[r, c] <= grad[r-1, c] or grad[r, c] <= grad[r+1, c]:gradSup[r, c] = 0if tq == 3: #3 is NW-SEif grad[r, c] <= grad[r-1, c-1] or grad[r, c] <= grad[r+1, c+1]:gradSup[r, c] = 0#Double thresholdstrongEdges = (gradSup > highThreshold)#Strong has value 2, weak has value 1thresholdedEdges = np.array(strongEdges, dtype=np.uint8) + (gradSup > lowThreshold)#Tracing edges with hysteresis #Find weak edge pixels near strong edge pixelsfinalEdges = strongEdges.copy()currentPixels = []for r in range(1, im.shape[0]-1):for c in range(1, im.shape[1]-1): if thresholdedEdges[r, c] != 1:continue #Not a weak pixel#Get 3x3 patch localPatch = thresholdedEdges[r-1:r+2,c-1:c+2]patchMax = localPatch.max()if patchMax == 2:currentPixels.append((r, c))finalEdges[r, c] = 1#Extend strong edges based on current pixelswhile len(currentPixels) > 0:newPix = []for r, c in currentPixels:for dr in range(-1, 2):for dc in range(-1, 2):if dr == 0 and dc == 0: continuer2 = r+drc2 = c+dcif thresholdedEdges[r2, c2] == 1 and finalEdges[r2, c2] == 0:#Copy this weak pixel to final resultnewPix.append((r2, c2))finalEdges[r2, c2] = 1currentPixels = newPixreturn finalEdgesif __name__=="__main__":im = imread("test.jpg", mode="L") #Open image, convert to greyscalefinalEdges = CannyEdgeDetector(im)imshow(finalEdges)