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
#include
#include
#include
#include
#include
#include
#include
#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_twidth;
int32_theight;
uint16_t nplanes;
uint16_t bitspp;
uint32_t compress_type;
uint32_t bmp_bytesz;
int32_thres;
int32_tvres;
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 header
if (fread(&bitmapFileHeader, sizeof(bmpfile_header_t),
1, filePtr) != 1) {
fclose(filePtr);
return NULL;
}
// read the bitmap info header
if (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 data
if (fseek(filePtr, bitmapFileHeader.bmp_offset, SEEK_SET)) {
fclose(filePtr);
return NULL;
}
// allocate enough memory for the bitmap image data
pixel_t *bitmapImage = malloc(bitmapInfoHeader->bmp_bytesz *
sizeof(pixel_t));
// verify memory allocation
if (bitmapImage == NULL) {
fclose(filePtr);
return NULL;
}
// read in the bitmap image data
size_t pad, count=0;
unsigned char c;
pad = 4*ceil(bitmapInfoHeader->bitspp*bitmapInfoHeader->width/32.) - bitmapInfoHeader->width;
for(size_t i=0;
iheight;
i++){
for(size_t j=0;
jwidth;
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 data
fclose(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;
}
// Palette
for (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_storage
size_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 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 array
fprintf(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 deg
nms[c] = G[c];
else
nms[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 edges
out[c] = MAX_BRIGHTNESS;
int nedges = 1;
edges[0] = c;
do {
nedges--;
const int t = edges[nedges];
int nbs[8];
// neighbours
nbs[0] = t - nx;
// nn
nbs[1] = t + nx;
// ss
nbs[2] = t + 1;
// ww
nbs[3] = t - 1;
// ee
nbs[4] = nbs[0] + 1;
// nw
nbs[5] = nbs[0] - 1;
// ne
nbs[6] = nbs[1] + 1;
// sw
nbs[7] = nbs[1] - 1;
// se
for (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 = https://www.it610.com/article/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 =
https://www.it610.com/article/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 main
import (
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"
// hue
if($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;
}
}
// saturation
if($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 detected
if($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
【算法|Canny边缘检测算法(Canny edge detector)】在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, imshow
def 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 noise
im2 = gaussian_filter(im, blur)
#Use sobel filters to get horizontal and vertical gradients
im3h = 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 direction
grad = 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 suppression
gradSup = grad.copy()
for r in range(im.shape[0]):
for c in range(im.shape[1]):
#Suppress pixels at the image edge
if r == 0 or r == im.shape[0]-1 or c == 0 or c == im.shape[1] - 1:
gradSup[r, c] = 0
continue
tq = thetaQ[r, c] % 4
if 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] = 0
if tq == 1: #1 is NE-SW
if grad[r, c] <= grad[r-1, c+1] or grad[r, c] <= grad[r+1, c-1]:
gradSup[r, c] = 0
if 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] = 0
if tq == 3: #3 is NW-SE
if grad[r, c] <= grad[r-1, c-1] or grad[r, c] <= grad[r+1, c+1]:
gradSup[r, c] = 0
#Double threshold
strongEdges = (gradSup > highThreshold)
#Strong has value 2, weak has value 1
thresholdedEdges = np.array(strongEdges, dtype=np.uint8) + (gradSup > lowThreshold)
#Tracing edges with hysteresis
#Find weak edge pixels near strong edge pixels
finalEdges = 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 pixels
while 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: continue
r2 = r+dr
c2 = c+dc
if thresholdedEdges[r2, c2] == 1 and finalEdges[r2, c2] == 0:
#Copy this weak pixel to final result
newPix.append((r2, c2))
finalEdges[r2, c2] = 1
currentPixels = newPix
return finalEdges
if __name__=="__main__":
im = imread("test.jpg", mode="L") #Open image, convert to greyscale
finalEdges = CannyEdgeDetector(im)
imshow(finalEdges)
推荐阅读
- 人工智能|干货！人体姿态估计与运动预测
- 分析COMP122 The Caesar Cipher
- 技术|为参加2021年蓝桥杯Java软件开发大学B组细心整理常见基础知识、搜索和常用算法解析例题（持续更新...)
- C语言学习(bit)|16.C语言进阶——深度剖析数据在内存中的存储
- Python机器学习基础与进阶|Python机器学习--集成学习算法--XGBoost算法
- 数据结构与算法|【算法】力扣第 266场周赛
- 数据结构和算法|LeetCode 的正确使用方式
- leetcode|今天开始记录自己的力扣之路
- 人工智能|【机器学习】深度盘点（详细介绍 Python 中的 7 种交叉验证方法！）
- 网络|简单聊聊压缩网络