在极坐标中，圆的表示方式为：
x=x0+rcosθ
y=y0+rsinθ圆心为 (x0,y0) ， r 为半径，θ 为旋转度数，值范围为 0-359 。
如果给定圆心点和半径，则其它点是否在圆上，就能检测出来了。 在图像中，将每个非0像素点作为圆心点，以一定的半径进行检测， 如果有一个点在圆上，就对这个圆心累加一次。 如果检测到一个圆，那么这个圆心点就累加到最大，成为峰值。 因此，在检测结果中，一个峰值点，就对应一个圆心点。
霍夫圆检测的函数：
skimage.transform.hough_circle(image, radius)radius 是一个数组，表示半径的集合，如 [3，4，5，6]
返回一个3维的数组 (radius index, M, N) ，第一维表示半径的索引，后面两维表示图像的尺寸。
绘制两个圆形，用霍夫圆变换将它们检测出来
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from skimage import draw,transform,feature
In [3]:
img = np.zeros((250, 250,3), dtype=np.uint8)
以半径50画一个圆。
In [4]:
rr, cc = draw.circle_perimeter(60, 60, 50)
以半径60画一个圆。
In [5]:
rr1, cc1 = draw.circle_perimeter(150, 150, 60)
img[cc, rr,:] =255
img[cc1, rr1,:] =255
In [6]:
fig, (ax0,ax1) = plt.subplots(1,2, figsize=(8, 5))
ax0.imshow(img) #显示原图
ax0.set_title('origin image')
hough_radii = np.arange(50, 80, 5) #半径范围
hough_res =transform.hough_circle(img[:,:,0], hough_radii) #圆变换
centers = [] #保存所有圆心点坐标
accums = [] #累积值
radii = [] #半径
for radius, h in zip(hough_radii, hough_res):
#每一个半径值，取出其中两个圆
num_peaks = 2
peaks =feature.peak_local_max(h, num_peaks=num_peaks) #取出峰值
centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius] * num_peaks)
画出最接近的圆。
In [7]:
image =np.copy(img)
for idx in np.argsort(accums)[::-1][:2]:
center_x, center_y = centers[idx]
radius = radii[idx]
cx, cy =draw.circle_perimeter(center_y, center_x, radius)
image[cy, cx] =(255,0,0)
ax1.imshow(image)
ax1.set_title('detected image')
Out[7]:
Text(0.5, 1.0, 'detected image')
结果图如上：原图中的圆用白色绘制，检测出的圆用红色绘制。
In [8]:
import numpy as np
import matplotlib.pyplot as plt
from skimage import data, color,draw,transform,feature,util
裁剪原图片，检测 canny 边缘。
In [9]:
image = util.img_as_ubyte(data.coins()[0:95, 70:370])
edges =feature.canny(image, sigma=3, low_threshold=10, high_threshold=50)
画出最接近的5个圆。
In [10]:
fig, (ax0,ax1) = plt.subplots(1,2, figsize=(8, 5))
ax0.imshow(edges, cmap=plt.cm.gray) #显示canny边缘
ax0.set_title('original iamge')
hough_radii = np.arange(15, 30, 2) #半径范围
hough_res =transform.hough_circle(edges, hough_radii) #圆变换
centers = [] #保存中心点坐标
accums = [] #累积值
radii = [] #半径
for radius, h in zip(hough_radii, hough_res):
#每一个半径值，取出其中两个圆
num_peaks = 2
peaks =feature.peak_local_max(h, num_peaks=num_peaks) #取出峰值
centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius] * num_peaks)
#画出最接近的5个圆
image = color.gray2rgb(image)
for idx in np.argsort(accums)[::-1][:5]:
center_x, center_y = centers[idx]
radius = radii[idx]
cx, cy =draw.circle_perimeter(center_y, center_x, radius)
image[cy, cx] = (255,0,0)
ax1.imshow(image)
ax1.set_title('detected image')
Out[10]:
Text(0.5, 1.0, 'detected image')
椭圆变换是类似的，使用函数为：
skimage.transform.hough_ellipse(img,accuracy, threshold, min_size, max_size)输入参数：
img: 待检测图像。accuracy: 使用在累加器上的短轴二进制尺寸，是一个double型的值，默认为1。thresh: 累加器阈值，默认为4。min_size: 长轴最小长度，默认为4。max_size: 短轴最大长度，默认为None，表示图片最短边的一半。
返回一个 [(accumulator, y0, x0, a, b, orientation)] 数组，
accumulator 表示累加器， (y0,x0) 表示椭圆中心点， (a,b) 分别表示长短轴， orientation 表示椭圆方向。
检测出咖啡图片中的椭圆杯口
In [11]:
import matplotlib.pyplot as plt
from skimage import data,draw,color,transform,feature
加载图片，转换成灰度图并检测边缘，裁剪原图像，不然速度非常慢。
In [12]:
image_rgb = data.coffee()[0:220, 160:420]
image_gray = color.rgb2gray(image_rgb)
edges = feature.canny(image_gray, sigma=2.0, low_threshold=0.55, high_threshold=0.8)
执行椭圆变换。
In [13]:
result =transform.hough_ellipse(edges, accuracy=20, threshold=250,min_size=100, max_size=120)
根据累加器排序。
In [14]:
result.sort(order='accumulator')
估计椭圆参数，排完序后取最后一个。
In [15]:
best = list(result[-1])
yc, xc, a, b = [int(round(x)) for x in best[1:5]]
orientation = best[5]
在原图上画出椭圆。
In [16]:
cy, cx =draw.ellipse_perimeter(yc, xc, a, b, orientation)
在原图中用蓝色表示检测出的椭圆。
In [17]:
image_rgb[cy, cx] = (0, 0, 255)
分别用白色表示 canny 边缘，用红色表示检测出的椭圆，进行对比。
In [18]:
edges = color.gray2rgb(edges)
edges[cy, cx] = (250, 0, 0)
In [19]:
fig2, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(8, 4))
ax1.set_title('Original picture')
ax1.imshow(image_rgb)
Out[19]:
<matplotlib.image.AxesImage at 0x7f6dc49b03e0>
In [20]:
ax2.set_title('Edge (white) and result (red)')
ax2.imshow(edges)
plt.show()
霍夫椭圆变换速度非常慢，应避免图像太大。
