前言

这次要介绍边缘检测和霍夫转换一样参考[1]，而边缘检测的Sobel为本次重点，因Sobel运算在往后取得方向和角度是非常方便的一个算式。

---

这次也修正了坐标和矩形相关的程式逻辑如下。

矩形取得的endX和endY为实际索引位置，举例endX公式为，x + width - 1，若width等于0回传-1。矩形高度为实际高度。修正每个画矩形有关的取得公式。DrawLine8bit修正为自动换算画出线。

边缘运算子

简介

[1]边缘运算主要有Sobel、Scharr和Laplace，它们差异在于乘上的"核"不同，这里使用Sobel中的3x3的一阶核运算来当範例。

Sobel垂直方向3x3一阶"核"：。

Sobel水平方向3x3一阶"核"：。

而垂直方向与水平方向计算出来结果映射在二维系座标，可以取得角度、方向和长度等等。假设垂直方向为x，水平方向为y，则依照三角函数性质可得到tan = 对边 / 斜边 = y / x，反三角函数atan则可得到角度，而距离公式使用欧几里得距离公式(L2算法)，sqrt(x^2　+　y^2)，以上为合併后的性质，接着继续看Sobel的运算步骤。

运算步骤

填补图像。取得核。捲积处理。

程式码

Library.h
这里只使用一阶3x3运算，dx和dy为使用开关。

/*Sobel8bit Parameter:src= source of imagepur= purpose of imagewidth= Image widthheight= Image heightdx= x gradient switchdy= y gradient switch*/void Sobel8bit(C_UCHAE* src, int32_t* pur, C_UINT32 width, C_UINT32 height, const bool dx, const bool dy);

Library.cpp

void Library::Sobel8bit(C_UCHAE* src, int32_t* data, C_UINT32 width, C_UINT32 height, const bool dx, const bool dy){// 1. paddingC_UINT32 padWidth = width + 2;C_UINT32 padHeight = height + 2;UCHAE* padSrc = new UCHAE[padWidth * padHeight];ImagePadding8bit(src, padSrc, width, height, 1);// 2. set kernelImage srcImage(padSrc, padWidth, padHeight, MNDT::ImageType::GRAY_8BIT);int32_t kernels[9];SetSobelKernels(kernels, dx, dy);// 3. calculate convolutionfor (UINT32 row = 1; row < padWidth - 1; row++){for (UINT32 col = 1; col < padHeight - 1; col++){int32_t sum = 0;UINT32 kernelsIndex = 0;for (int32_t blockRow = -1; blockRow <= 1; blockRow++){for (int32_t blockCol = -1; blockCol <= 1; blockCol++){UCHAE pix = srcImage.image[row + blockRow][col + blockCol];sum += (static_cast<int32_t>(pix) * kernels[kernelsIndex]);kernelsIndex++;}}*data = sum;data++;}}delete[] padSrc;padSrc = nullptr;}

结果图

24bit垂直x

24bit水平y

边缘检测

简介

边缘检测主要是使用水平和垂直的绝对值(L1算法)合计后，再去做正规化显示即是边缘检测。而边缘检测方式有很多种例如，Laplacian和Canny等等，但其原理大同小异，部分的差别只在于"核"的不同，这边使用上述的Sobel运算子取得边缘。

运算步骤

取得水平和垂直计算的Sobel值。计算x和y绝对值的合，并取得目前最大值。正规化将最大值转为255。

程式码

Library.h

/*SobelEdgeView8bit Parameter:src= source of imagepur= purpose of imagewidth= Image widthheight= Image height*/void SobelEdgeView8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height);

Library.cpp

void Library::SobelEdgeView8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height){// 1. get sobelint32_t* Gx = new int32_t[width * height];int32_t* Gy = new int32_t[width * height];C_UCHAE* srcEnd = src + width * height;Sobel8bit(src, Gx, width, height, true, false);Sobel8bit(src, Gy, width, height, false, true);// 2. calculate abs and get maxint32_t max = 0;int32_t* data = new int32_t[width * height];int32_t* dataPointer = data;int32_t* GxPointer = Gx;int32_t* GyPointer = Gy;while (src < srcEnd){*dataPointer = abs(*GxPointer) + abs(*GyPointer);max = max < *dataPointer ? *dataPointer : max;dataPointer++;src++;GxPointer++;GyPointer++;}delete[] Gx;Gx = nullptr;delete[] Gy;Gy = nullptr;// 3. normalizationC_UCHAE* purEnd = pur + width * height;dataPointer = data;while (pur < purEnd){*pur = *dataPointer * 255 / max;dataPointer++;pur++;}delete[] data;data = nullptr;}

结果图

霍夫找线

参考[2]的OpenCV代码，主要将最后改为全部重新计算一次在画出结果。

---

简介

霍夫转换是计算x和y在0~180角度的极座标(离散值)，使用三角函数cos　=　邻边　/　斜边，sin　=　对边　/　斜边，而cos乘上斜边得取x，sin乘上斜边得取y，使用图片的位置乘上cos和sin，cos(theta)　*　x取得x在霍夫分布和sin(theta)　*　y取得y分布，两者相加即是霍夫转换的离散值，如下图一，图二则是实际图片输出的分布有兴趣可找函数DrawHoughPolar。

输入资料[x, y]由上而下[7, 7], [5, 5], [3, 3], [1, 1]，可以清楚的看到都相交于弧度2.35(135度)位置，利用此特性统计直线的角度，将可能为直线的角度列出就可以减少计算量。

图一。

图二，整张图片实际画出的线。

初始化

thetaSize：分割数量。PI除theta，theta为每个切割弧度大小。fixRho：三角函数倍率。1除rho，rho为倍率参数。originalR：最大离散值。最大为长宽乘上0.7...(45度)或者高或宽乘上1。maxR：修正的最大离散值。为了算邻近，加上填补和修正索引位置。xAxisOffset：负数的偏移位置。主要用来将离散值分为正和负两部分，所以每一个分割角度的範围是[-originalR,originalR]，。

运算步骤

初始化sin和cos的倍率与弧度分割数量。像素大于0的座标做霍夫转换，并且累加计算出来的离散位置。使用四方邻近法与上下左右的值比对，若该点可能是线则数值会比较大，若不是累加设为0。像素大于0的座标做霍夫转换，若累积不为0则输出白色。(可以改为依照count数量排序后，依照取得前N个再去跑迴圈效率会较高)。

程式码

Library.h

/*HoughLines Parameter:src= source of imagepur= purpose of imagewidth= Image widthheight= Image heightrho= calculation sin and costheta= calculation split parmasthreshold= calculation output parmas*/void HoughLines(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_FLOAT rho, C_FLOAT theta, C_UINT32 threshold);        void HoughLineNeighboursUpdate(UINT32* count, C_UINT32 thetaSize, C_UINT32 maxR, C_UINT32 threshold);

Library.cpp

void Library::HoughLines(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_FLOAT rho, C_FLOAT theta, C_UINT32 threshold){// 1. init paramsC_FLOAT fixRho = 1.0f / rho;C_UINT32 thetaSize = static_cast<UINT32>(MNDT::PI / theta);float* fixSin = new float[thetaSize];float* fixCos = new float[thetaSize];for (UINT32 index = 0; index < thetaSize; index++){fixSin[index] = MNDT::FixValue(sin(theta * index)) * fixRho;fixCos[index] = MNDT::FixValue(cos(theta * index)) * fixRho;}// 2. hough totalImage srcImage(const_cast<UCHAE*>(src), width, height, MNDT::ImageType::GRAY_8BIT);C_DOUBLE originalR = std::max((width + height) * sin(MNDT::PI / 4.0), static_cast<double>(std::max(width, height)));C_UINT32 maxR = static_cast<UINT32>(originalR * fixRho) + 1 + 2;C_UINT32 xAxisOffset = maxR * (thetaSize + 2);UINT32* count = new UINT32[xAxisOffset * 2]{ 0 };for (UINT32 row = 0; row < height; row++){for (UINT32 col = 0; col < width; col++){if (srcImage.image[row][col] > 0){for (UINT32 index = 0; index < thetaSize; index++){int32_t r = static_cast<int32_t>(fixSin[index] * row + fixCos[index] * col);r = r < 0 ? abs(r) + xAxisOffset : r;count[maxR * (index + 1) + r + 1]++;}}}}//// draw hough//DrawHoughPolar(pur//, width, height//, count//, theta, maxR);// 3. update neighboursHoughLineNeighboursUpdate(count, thetaSize, maxR, threshold);// 4. draw lineImage purImage(pur, width, height, MNDT::ImageType::GRAY_8BIT);for (UINT32 row = 0; row < height; row++){for (UINT32 col = 0; col < width; col++){if (srcImage.image[row][col] > 0){for (UINT32 index = 0; index < thetaSize; index++){int32_t r = static_cast<int32_t>(fixSin[index] * row + fixCos[index] * col);r = r < 0 ? abs(r) + xAxisOffset : r;if (count[maxR * (index + 1) + r + 1] > 0){purImage.image[row][col] = 255;}}}}}delete[] fixSin;fixSin = nullptr;delete[] fixCos;fixCos = nullptr;delete[] count;count = nullptr;}void Library::HoughLineNeighboursUpdate(UINT32* count, C_UINT32 thetaSize, C_UINT32 maxR, C_UINT32 threshold){// 4邻比较C_UINT32 xAxisOffset = maxR * (thetaSize + 2);for (UINT32 index = 0; index < thetaSize; index++){C_UINT32 offset = maxR * (index + 1);for (UINT32 countIndex = offset + 1; countIndex < offset + maxR - 1; countIndex++){// x axis -> + and -for (UINT32 axis = 0; axis < 2; axis++){UINT32* nowCountPointer = count + countIndex + xAxisOffset * axis;if (*nowCountPointer < threshold|| *nowCountPointer < *(nowCountPointer - 1)|| *nowCountPointer < *(nowCountPointer + 1)|| *nowCountPointer < *(nowCountPointer - maxR)|| *nowCountPointer < *(nowCountPointer + maxR)){*nowCountPointer = 0;}}}}}

原图

结果

霍夫找圆

参考[3]的OpenCV代码，在用更直观方式实现。

---

简介

圆的方式与直线相同，公式改为，r^2　=　x^2　+　y^2，而在OpenCV当中先统计有可能为圆心的位置，在去判别哪些是真正的圆心。首先找到大于0的像素，依照Sobel运算可取得梯度，而在这里的梯度方向即是切线的法线(与切线垂直)，如此一来只需要往梯度方向和反方向累加通过的可能圆心点即可，如下图一原图，图二为累加的所有路径。

图一

图二

运算步骤

计算垂直x和水平y的Sobel运算子。使用L2正规化梯度偏移量，使用count阵列统计像素大于0从minRadius到maxRadius经过的直线。这里多用一个函数计算点到点之间经过的距离。使用四方邻近法与上下左右的值比对，若该点可能是线则数值会比较大，若不是累加设为0。count大于0的座标从minRadius到maxRadius判断360度的位置是否为大于0的像素，若360度都是即是原心。(可以改为依照count数量排序后，依照取得前N个再去跑迴圈效率会较高)。

程式码

Library.h

/*HoughCircles Parameter:src= source of imagepur= purpose of imagewidth= Image widthheight= Image heightminRadius= min radiusmaxRadius= max radiusthreshold= calculation output parmas*/void HoughCircles(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_FLOAT minRadius, C_FLOAT maxRadius, C_UINT32 threshold);            void HoughCirclesCount(C_UCHAE* src, C_UINT32 width, C_UINT32 height, UINT32* count, C_FLOAT minRadius, C_FLOAT maxRadius);void HoughCirclePointCount(UINT32* count, C_UINT32 width, const Point& p1, const Point& p2);void HoughCircleNeighboursUpdate(UINT32* count, C_UINT32 width, C_UINT32 height, C_UINT32 threshold);

Library.cpp

void Library::HoughCircles(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_FLOAT minRadius, C_FLOAT maxRadius, C_UINT32 threshold){// 1. totalUINT32* count = new UINT32[(width + 2) * (height + 2)]{ 0 };HoughCirclesCount(src, width, height, count, minRadius, maxRadius);// 2. update neighboursHoughCircleNeighboursUpdate(count, width, height, threshold);// 3. drawImage srcImage(const_cast<UCHAE*>(src), width, height, MNDT::ImageType::GRAY_8BIT);Image purImage(pur, width, height, MNDT::ImageType::GRAY_8BIT);for (UINT32 row = 0; row < height; row++){C_UINT32 rowIndex = (row + 1) * width;for (UINT32 col = 0; col < width; col++){C_UINT32 nowIndex = rowIndex + col + 1;if (count[nowIndex] == 0){continue;}for (float radius = minRadius; radius <= maxRadius; radius++){bool nCircle = false;for (float pi = 0; pi <= MNDT::PI * 2; pi += 0.2f){int32_t y = static_cast<int32_t>(static_cast<float>(row) + MNDT::FixValue(sin(pi)) * radius);int32_t x = static_cast<int32_t>(static_cast<float>(col) + MNDT::FixValue(cos(pi)) * radius);if (x < 0 || y < 0|| (unsigned)x >= width || (unsigned)y >= height|| srcImage.image[y][x] == 0){nCircle = true;break;}}if (nCircle){continue;}for (float pi = 0; pi <= MNDT::PI * 2; pi += 0.2f){C_UINT32 y = static_cast<UINT32>(row + MNDT::FixValue(sin(pi)) * radius);C_UINT32 x = static_cast<UINT32>(col + MNDT::FixValue(cos(pi)) * radius);purImage.image[y][x] = 255;}}}}delete[] count;count = nullptr;}void Library::HoughCirclesCount(C_UCHAE* src, C_UINT32 width, C_UINT32 height, UINT32* count, C_FLOAT minRadius, C_FLOAT maxRadius){assert(count != nullptr);// get gradientLibrary lib;int32_t* Gx = new int32_t[width * height];int32_t* Gy = new int32_t[width * height];lib.Sobel8bit(src, Gx, width, height, true, false);lib.Sobel8bit(src, Gy, width, height, false, true);// calculate gradient totalImage srcImage(const_cast<UCHAE*>(src), width, height, MNDT::ImageType::GRAY_8BIT);for (UINT32 row = 0; row < height; row++){C_UINT32 rowIndex = row * width;for (UINT32 col = 0; col < width; col++){C_UINT32 index = rowIndex + col;if (srcImage.image[row][col] > 0){C_FLOAT magnitude = static_cast<float>(sqrt(Gx[index] * Gx[index] + Gy[index] * Gy[index] + 0.0001f));float offsetX = Gx[index] / magnitude;float offsetY = Gy[index] / magnitude;for (UINT32 sign = 0; sign < 2; sign++){float x = col + minRadius * offsetX;float y = row + minRadius * offsetY;for (float radius = minRadius; radius <= maxRadius; radius++){if (x < 0 || y < 0|| x >= width || y >= height){continue;}//MNDT::DrawLine8bit(purImage, Point(col, row), Point(static_cast<UINT32>(x), static_cast<UINT32>(y)));HoughCirclePointCount(count, width, Point(col, row), Point(static_cast<UINT32>(x), static_cast<UINT32>(y)));x += offsetX;y += offsetY;}offsetX = -offsetX;offsetY = -offsetY;}}}}delete[] Gx;Gx = nullptr;delete[] Gy;Gy = nullptr;}void Library::HoughCirclePointCount(UINT32* count, C_UINT32 width, const Point& p1, const Point& p2){C_UINT32 absDiffX = abs(static_cast<int32_t>(p1.X()) - static_cast<int32_t>(p2.X()));C_UINT32 absDiffY = abs(static_cast<int32_t>(p1.Y()) - static_cast<int32_t>(p2.Y()));C_UINT32 base = absDiffX > absDiffY ? absDiffY : absDiffX;C_INT32 diffX = static_cast<int32_t>(p1.X()) - static_cast<int32_t>(p2.X());C_INT32 diffY = static_cast<int32_t>(p1.Y()) - static_cast<int32_t>(p2.Y());C_INT32 baseX = diffX < 0 ? 1 : -1;C_INT32 baseY = diffY < 0 ? 1 : -1;int32_t x = p1.X();int32_t y = p1.Y();for (UINT32 index = 0; index < base; index++){count[static_cast<UINT32>(static_cast<UINT32>(y + 1) * width + x + 1)]++;x += baseX;y += baseY;}}void Library::HoughCircleNeighboursUpdate(UINT32* count, C_UINT32 width, C_UINT32 height, C_UINT32 threshold){// 4邻比较for (UINT32 row = 0; row < height; row++){C_UINT32 rowIndex = (row + 1) * width;for (UINT32 col = 0; col < width; col++){C_UINT32 nowIndex = rowIndex + col + 1;UINT32* nowCountPointer = count + nowIndex;if (*nowCountPointer < threshold|| *nowCountPointer < *(nowCountPointer - 1)|| *nowCountPointer < *(nowCountPointer + 1)|| *nowCountPointer < *(nowCountPointer - width - 2)|| *nowCountPointer < *(nowCountPointer + width + 2)){*nowCountPointer = 0;}}}}

结果图

原图

结果图

结语

这次做的霍夫转换直接把全部的可能输出效率会较差，但对于理解原理没有太大的影响。在[1]有提到分水岭算法它是利用边缘製作出一个标记表，再利用像素差进行合併，由此可知Sobel运算子运用的地方非常广泛，接下来会介绍影像特徵取得，希望在月底能把AdaBoost做完。若有问题或错误欢迎留言指教。

---

C#视窗原始码
C++函数原始码

参考文献

[1]阿洲(2015). OpenCV教学 | 阿洲的程式教学 from: http://monkeycoding.com/?page_id=12 (2018.11.21).
[2]viewcode(2012) 霍夫变换直线检测houghlines及opencv的实现分析 from: https://blog.csdn.net/viewcode/article/details/8090932 (2018.11.21).
[3]zhaocj(2016). Opencv2.4.9源码分析——HoughCircles from: https://blog.csdn.net/zhaocj/article/details/50454847 (2018.11.22)