圆检测的最小二乘法
Least squares for circle detection
我正在尝试对数据点的子集使用最小二乘优化从激光扫描中执行圆检测。由于仅对圆的一部分进行测量,因此最小二乘法 returns 是错误的结果,报告的圆比实际距离激光更近。
算法的结果如图所示。散点表示激光测量值,圆圈以算法返回的点为中心。灰色半透明形状表示机器人进行扫描(激光在该形状的左侧和右侧)。
我只对已知半径 RR 的圆的局部坐标感兴趣。
PS。我假设扫描被分成簇(self.clusters[i] 是一个簇),它们是 [x,y] 激光点的列表
def circle(x, scan):
xc, yc = x
f = sqrt((scan[:,0] - xc)**2 + (scan[:,1] - yc)**2) - RR
return f
def optimize_detect_circles(self):
centre = [1,1]
for i in range(0, self.number_of_clusters):
range_points = np.array(self.clusters[i])
sol = optimize.root(circle, centre, args=(range_points), method='lm')
self.circle_candidates.append(sol.x)
print sol.x
图片如下:
1,1
离正确值太远;你很可能陷入局部最优。
尝试从离真实中心更近的点开始。您可以通过首先将一条直线拟合到您的集群来找到它;然后将点分成两半,根据它们投射到线的哪一半;接下来适合 two 行,每行到新的两个子集群之一;求两条垂线在中点的交点。
这是基于您的集群是跨度不大于 180 度的弧形,它们看起来确实如此。如果不是,只需重复细分,得到四个和弦而不是两个。
您可以使用 Circular Hough Transform 找到圆 - 如果您事先知道圆的半径,这将特别容易。
从 skimage 上有关此过程的文档中大量借用,我整理了以下可运行的代码,但可能需要进行一些调整才能找到圆圈:
import numpy as np
import matplotlib.pyplot as plt
import skimage
from skimage import data, filter, io
from skimage.transform import hough_circle
from skimage.feature import peak_local_max
from skimage import data, color
from skimage.draw import circle_perimeter
theImage = np.sum(io.imread("w1s31.png"),2)/4 # Image to greyscale
hough_radii = np.arange(61, 69, 2) # These are the radii to search for
hough_res = hough_circle(theImage, hough_radii)
centers = []
accums = []
radii = []
for radius, h in zip(hough_radii, hough_res):
# For each radius, extract two circles
num_peaks = 2
peaks = peak_local_max(h, num_peaks=num_peaks)
centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius] * num_peaks)
# Draw the most prominent 5 circles
image = color.gray2rgb(theImage)
for idx in np.argsort(accums)[::-1][:5]:
center_x, center_y = centers[idx]
radius = radii[idx]
cx, cy = circle_perimeter(center_y, center_x, radius)
theImage[cy, cx] = (1)
plt.imshow(theImage, cmap=plt.cm.gray)
plt.show()
这是 this 论文中实现的圆检测的一小段。
该算法始终有效,并且不会以线性复杂度迭代。
var add = function(a,b){ return a+b; };
var add_uu = function(a,b){ return a+b[0]*b[0]; };
var add_uv = function(a,b){ return a+b[0]*b[1]; };
var add_vv = function(a,b){ return a+b[1]*b[1]; };
var add_uuu = function(a,b){ return a+b[0]*b[0]*b[0]; };
var add_vvv = function(a,b){ return a+b[1]*b[1]*b[1]; };
var add_uvv = function(a,b){ return a+b[0]*b[1]*b[1]; };
var add_uuv = function(a,b){ return a+b[0]*b[0]*b[1]; };
var getx = function(e){ return e[0]; };
var gety = function(e){ return e[1]; };
$(document).ready(function () {
var paper = Raphael("canvas");
var points=[];
$("#canvas").mousedown(function (e) {
var x = e.offsetX;
var y = e.offsetY;
points.push([x,y]);
paper.circle(x, y, 1);
});
$("#clear").click(function(){
paper.clear();
points = [];
});
$("#go").click(function(){
var N = points.length;
var xb = points.map(getx).reduce(add,0) / N;
var yb = points.map(gety).reduce(add,0) / N;
var u = points.map(function(e){return [e[0]-xb,e[1]-yb];});
var a1 = u.reduce(add_uu,0);
var b1 = u.reduce(add_uv,0);
var c1 = 0.5*(u.reduce(add_uuu,0) + u.reduce(add_uvv,0));
var a2 = u.reduce(add_uv,0);
var b2 = u.reduce(add_vv,0);
var c2 = 0.5*(u.reduce(add_vvv,0) + u.reduce(add_uuv,0));
var q = a2/a1;
var vc = (c2-q*c1)/(b2-q*b1);
var uc = (c1-b1*vc)/a1;
var r = Math.sqrt(uc*uc+vc*vc+(a1+b2)/N);
var x = uc+xb;
var y = vc+yb;
paper.circle(x, y, r).attr({"stroke":"#f00","stroke-width":2});
var e = points.reduce(function(p,c) {
var t = r*r - (c[0]-x)*(c[0]-x) - (c[1]-y)*(c[1]-y);
return p+t*t;
},0);
console.log("Residue = " + e);
});
});
#canvas {
width: 600px;
height: 400px;
border: 2px dotted #ccc;
cursor: crosshair;
}
<script src="http://cdnjs.cloudflare.com/ajax/libs/raphael/2.1.0/raphael-min.js"></script>
<script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script>
<button id="clear">Clear</button>
<button id="go">Compute circle</button>
<div id="canvas"></div>
我正在尝试对数据点的子集使用最小二乘优化从激光扫描中执行圆检测。由于仅对圆的一部分进行测量,因此最小二乘法 returns 是错误的结果,报告的圆比实际距离激光更近。
算法的结果如图所示。散点表示激光测量值,圆圈以算法返回的点为中心。灰色半透明形状表示机器人进行扫描(激光在该形状的左侧和右侧)。
我只对已知半径 RR 的圆的局部坐标感兴趣。
PS。我假设扫描被分成簇(self.clusters[i] 是一个簇),它们是 [x,y] 激光点的列表
def circle(x, scan):
xc, yc = x
f = sqrt((scan[:,0] - xc)**2 + (scan[:,1] - yc)**2) - RR
return f
def optimize_detect_circles(self):
centre = [1,1]
for i in range(0, self.number_of_clusters):
range_points = np.array(self.clusters[i])
sol = optimize.root(circle, centre, args=(range_points), method='lm')
self.circle_candidates.append(sol.x)
print sol.x
图片如下:
1,1
离正确值太远;你很可能陷入局部最优。
尝试从离真实中心更近的点开始。您可以通过首先将一条直线拟合到您的集群来找到它;然后将点分成两半,根据它们投射到线的哪一半;接下来适合 two 行,每行到新的两个子集群之一;求两条垂线在中点的交点。
这是基于您的集群是跨度不大于 180 度的弧形,它们看起来确实如此。如果不是,只需重复细分,得到四个和弦而不是两个。
您可以使用 Circular Hough Transform 找到圆 - 如果您事先知道圆的半径,这将特别容易。
从 skimage 上有关此过程的文档中大量借用,我整理了以下可运行的代码,但可能需要进行一些调整才能找到圆圈:
import numpy as np
import matplotlib.pyplot as plt
import skimage
from skimage import data, filter, io
from skimage.transform import hough_circle
from skimage.feature import peak_local_max
from skimage import data, color
from skimage.draw import circle_perimeter
theImage = np.sum(io.imread("w1s31.png"),2)/4 # Image to greyscale
hough_radii = np.arange(61, 69, 2) # These are the radii to search for
hough_res = hough_circle(theImage, hough_radii)
centers = []
accums = []
radii = []
for radius, h in zip(hough_radii, hough_res):
# For each radius, extract two circles
num_peaks = 2
peaks = peak_local_max(h, num_peaks=num_peaks)
centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
radii.extend([radius] * num_peaks)
# Draw the most prominent 5 circles
image = color.gray2rgb(theImage)
for idx in np.argsort(accums)[::-1][:5]:
center_x, center_y = centers[idx]
radius = radii[idx]
cx, cy = circle_perimeter(center_y, center_x, radius)
theImage[cy, cx] = (1)
plt.imshow(theImage, cmap=plt.cm.gray)
plt.show()
这是 this 论文中实现的圆检测的一小段。
该算法始终有效,并且不会以线性复杂度迭代。
var add = function(a,b){ return a+b; };
var add_uu = function(a,b){ return a+b[0]*b[0]; };
var add_uv = function(a,b){ return a+b[0]*b[1]; };
var add_vv = function(a,b){ return a+b[1]*b[1]; };
var add_uuu = function(a,b){ return a+b[0]*b[0]*b[0]; };
var add_vvv = function(a,b){ return a+b[1]*b[1]*b[1]; };
var add_uvv = function(a,b){ return a+b[0]*b[1]*b[1]; };
var add_uuv = function(a,b){ return a+b[0]*b[0]*b[1]; };
var getx = function(e){ return e[0]; };
var gety = function(e){ return e[1]; };
$(document).ready(function () {
var paper = Raphael("canvas");
var points=[];
$("#canvas").mousedown(function (e) {
var x = e.offsetX;
var y = e.offsetY;
points.push([x,y]);
paper.circle(x, y, 1);
});
$("#clear").click(function(){
paper.clear();
points = [];
});
$("#go").click(function(){
var N = points.length;
var xb = points.map(getx).reduce(add,0) / N;
var yb = points.map(gety).reduce(add,0) / N;
var u = points.map(function(e){return [e[0]-xb,e[1]-yb];});
var a1 = u.reduce(add_uu,0);
var b1 = u.reduce(add_uv,0);
var c1 = 0.5*(u.reduce(add_uuu,0) + u.reduce(add_uvv,0));
var a2 = u.reduce(add_uv,0);
var b2 = u.reduce(add_vv,0);
var c2 = 0.5*(u.reduce(add_vvv,0) + u.reduce(add_uuv,0));
var q = a2/a1;
var vc = (c2-q*c1)/(b2-q*b1);
var uc = (c1-b1*vc)/a1;
var r = Math.sqrt(uc*uc+vc*vc+(a1+b2)/N);
var x = uc+xb;
var y = vc+yb;
paper.circle(x, y, r).attr({"stroke":"#f00","stroke-width":2});
var e = points.reduce(function(p,c) {
var t = r*r - (c[0]-x)*(c[0]-x) - (c[1]-y)*(c[1]-y);
return p+t*t;
},0);
console.log("Residue = " + e);
});
});
#canvas {
width: 600px;
height: 400px;
border: 2px dotted #ccc;
cursor: crosshair;
}
<script src="http://cdnjs.cloudflare.com/ajax/libs/raphael/2.1.0/raphael-min.js"></script>
<script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script>
<button id="clear">Clear</button>
<button id="go">Compute circle</button>
<div id="canvas"></div>