0 引言
圆已经成为了生活中最常见的几何图形，抬头望去，就能看见非常多由圆构成的东西，如水杯、碗、汤圆、鸡蛋等。离开生活，在各种各样的光学系统中，圆也是最常见的，照相机的光圈，望远镜的物镜，显微镜的照明光等等，都可以由圆这样一个优雅的图形所描述。因此，在图像处理领域，提取圆的参数，例如圆心和半径，已经成为了非常常见的操作。
以matlab为例，其自带的imfindcircles()函数就可以在图像中找到圆并返回参数，但其现有的缺点也制约了圆心提取的效率和可行性。为此，探求更高效率、更具有鲁棒性的圆心提取算法成为了一个热点话题，本文结合我在科研中遇到的复杂圆心提取问题，介绍利用matlab提取圆心的各种方法，并给出一个综合性的解决方案。
1 自带函数
1.1 函数简介
matlab软件有一个内置函数imfindcircles()，可以用于圆心提取，其使用方法如下：
[centers,radii,metric] = imfindcircles(A,radiusRange,'Name','Value')其中，输入参数解释如下：
- A ，输入的图片，即用于检测圆形物体的图像，可以是灰度、真彩色图像或二值图像。
- radiusRange, 要检测圆的半径，或者检测的圆形目标的近似半径，指定为正数。也可以是某个范围,指定为
[rmin rmax]形式的由正整数组成的二元素向量，其中rmin小于rmax。 - Name和value参数，用于设定某些偏好。例如对象极性‘ObjectPolarity’，可以设为'bright'（默认）或者'dark'，表明圆比背景亮还是暗；计算方法'method'，可以设为‘PhaseCode’ （默认）和’TwoStage’；敏感性因子‘Sensitivity’，设定为[0 1]区间的标量值；边缘阈值‘EdgeThreshold’，设定为[0 1]区间的标量值。
输出参数解释如下：
- centers，圆心坐标，返回为 P×2 矩阵，第一列中包含圆心的 x 坐标，第二列中包含 y 坐标。行数 P 是检测到的圆形的个数。centers 根据圆形的强度排序。
- radii，各圆心对应的估计半径，以列向量形式返回。
radii(j)处的半径值对应于以centers(j,:)为中心的圆形。 - metric，圆形的强度，可以认为是正确性或者置信度，以列向量形式返回。
metric(j)处的值对应于以centers(j,:)为中心、半径为radii(j)的圆。
1.2 函数缺点
内置函数imfindcircles()，为简单的图像处理提供了便捷，但是也存在一些显著的问题，使得其难以满足科学研究的需要，主要局限如下：
-
当 radius（或
rmin）的值小于或等于 5 时，imfindcircles的准确度会受到限制。 -
'PhaseCode'和'TwoStage'这两种计算方法检测同心圆的能力有限。同心圆的结果可能因输入图像而异。 -
imfindcircles找不到圆心位于图像区域之外的圆形。 -
imfindcircles会预处理二值（逻辑值）图像以提高结果的准确度。在处理真彩色图像之前，它使用rgb2gray函数将其转换为灰度图像。 -
圆心计算的时间会随着圆半径的范围变大而增大，计算负载较大。
除了使用自带函数之外，回归数学，也还存在很多优质的算法，下面将详细介绍。
2 公式法
所谓公式法，就是给定三个不共线的点，然后求解圆心。具体代码来自此处，如下：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 由圆上三点确定圆心和半径
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT
% p1 : - 第一个点坐标, 行向量 1x3
% p2 : - 第二个点坐标, 行向量 1x3
% p3 : - 第三个点坐标, 行向量 1x3
% 若输入1x2的行向量, 末位自动补0, 变为1x3的行向量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OUTPUT
% pc : - 圆心坐标, 行向量 1x3
% r : - 半径, 标量
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 调用示例1 - 平面上三个点
% [pc1,r1]=points2circle([1,2],[-2,1],[0,-3])
% 调用示例2 - 空间中三个点
% [pc2,r2]=points2circle([1,2,-1],[-2,1,2],[0,-3,-3])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pc,r]=points2circle(p1,p2,p3)
% 输入检查
validateattributes(p1,{'numeric'},{'row'},1);% 行向量
validateattributes(p2,{'numeric'},{'row'},2);
validateattributes(p3,{'numeric'},{'row'},3);
num1=length(p1);num2=length(p2);num3=length(p3);
if (num1 == num2) && (num2 == num3)
if num1 == 2
p1=[p1,0];p2=[p2,0];p3=[p3,0];
elseif num1 ~= 3
error('仅支持二维或三维坐标输入');
end
else
error('输入坐标的维数不一致');
end
% 共线检查
temp01=p1-p2;temp02=p3-p2;
temp03=cross(temp01,temp02);
temp=(temp03*temp03')/(temp01*temp01')/(temp02*temp02');
if temp < 10^-6
error('三点共线, 无法确定圆');
end
mat1=[p1,1;p2,1;p3,1];% size = 3x4
m=+det(mat1(:,2:4));
n=-det([mat1(:,1),mat1(:,3:4)]);
p=+det([mat1(:,1:2),mat1(:,4)]);
q=-det(mat1(:,1:3));
mat2=[[p1*p1';p2*p2';p3*p3'],mat1;2*q,[-m,-n,-p,0]];% size = 4x5
A=+det(mat2(:,2:5));
B=-det([mat2(:,1),mat2(:,3:5)]);
C=+det([mat2(:,1:2),mat2(:,4:5)]);
D=-det([mat2(:,1:3),mat2(:,5)]);
E=+det(mat2(:,1:4));
pc=-[B,C,D]/2/A;
r=sqrt(B^2+C^2+D^2-4*A*E)/2/abs(A);
end
通过代码可以看出，公式法也存在问题，即三点共线时无法确定一个圆，但是对实际图像作边缘检测时，经常会出现很多共线的点。例如，对图1所示的光斑作边缘检测时，会得到图2所示的共线边缘点，此时如何选择会成为一个问题。
3 最小二乘法
在难以选择时，最小二乘法往往可以提供帮助，提供一个全局最优解。详细的数学推导和代码可以参考这篇博客，代码：
function [ p ] = Circle_Fitting( XZ )
% use data XZ to fit a circles whose center and radius are (p(1),p(2)) and
% p(3) respectively
N = size(XZ,1);
x = XZ(:,1);
z = XZ(:,2);
sum_X_Raw = 0;
sum_Z_Raw = 0;
sum_XSquare_Raw = 0;
sum_ZSquare_Raw = 0;
sum_XCube_Raw = 0;
sum_ZCube_Raw = 0;
sum_XZZ_Raw = 0;
sum_XZ_Raw = 0;
sum_XXZ_Raw = 0;
for i=1:N
sum_X_Raw = sum_X_Raw+x(i);
sum_Z_Raw = sum_Z_Raw+z(i);
sum_XSquare_Raw = sum_XSquare_Raw+x(i)*x(i);
sum_ZSquare_Raw = sum_ZSquare_Raw+z(i)*z(i);
sum_XCube_Raw = sum_XCube_Raw+x(i)*x(i)*x(i);
sum_ZCube_Raw = sum_ZCube_Raw+z(i)*z(i)*z(i);
sum_XZ_Raw = sum_XZ_Raw+x(i)*z(i);
sum_XZZ_Raw = sum_XZZ_Raw+x(i)*z(i)*z(i);
sum_XXZ_Raw = sum_XXZ_Raw+x(i)*x(i)*z(i);
end
D = N*sum_XZ_Raw-sum_X_Raw*sum_Z_Raw;
C = N*sum_XSquare_Raw-sum_X_Raw*sum_X_Raw;
E = N*sum_XCube_Raw+N*sum_XZZ_Raw-(sum_XSquare_Raw+sum_ZSquare_Raw)*sum_X_Raw;
G = N*sum_ZSquare_Raw-sum_Z_Raw*sum_Z_Raw;
H = N*sum_ZCube_Raw+N*sum_XXZ_Raw-(sum_XSquare_Raw+sum_ZSquare_Raw)*sum_Z_Raw;
a = (H*D-E*G)/(C*G-D*D);
b = (H*C-E*D)/(D*D-G*C);
c = -((sum_XSquare_Raw+sum_ZSquare_Raw)+a*sum_X_Raw+b*sum_Z_Raw)/N;
p(1) = -0.5*a;
p(2) = -0.5*b;
p(3) = 0.5*sqrt(a*a+b*b-4*c);
end
最小二乘无疑是一个万精油的方法，但也还是有一个问题，就是用于计算的数据如果包含了很多噪点，那么这个方法也会把噪点考虑进去。因此，对原始数据进行去噪是必须的。
4 RANSAC算法
RANSAC算法，即随机采样一致性算法，其实某种程度上来说也算是一种Hough变换的方法，但是其计算圆时可以采用公式法，总体的速度更快。其基本步骤如下：
- 假定总迭代次数为K，初始迭代次数k=0,计数器c=0;
- 从边缘点中随机选取三个点，再利用公式法拟合得到一个圆心坐标和半径。
- 计算剩余点到圆心的距离，并与第1步得到的半径进行比较，如果两者相差较小，例如小于5个像素，既可以认为这个点位于圆上，计数器值c加1。所有点都比较完时，当前迭代次数k+1。
- 重复2-3步，指导k=K时结束迭代。
- 选取计数器最大的那一次迭代，并将对应的圆心坐标和半径作为拟合圆的参数。
代码如下：
clc;
close all;
clear all;
path = 'G:\SegmentationClass\13.tif';
spot = imread(path);
bw_spot = imbinarize(spot); % change the format of sopt image to binary
edge_spot = edge(bw_spot,'canny');% detect the edge of spot
figure; imshow(edge_spot,[]);
[y_edge,x_edge] = find(edge_spot);% coordinates of edge points
N_points = length(y_edge);% number of edge points
K = 60;% iterated times
r_error_th = 15;% the error of estimated r
random_point = randi(N_points,K,3);% generate K rows, 3 columns random ints whose maxmum is N_point
metric = zeros(K,5);
%% step1: estimate initial r,xc,yc
% tic
for k = 1:K
p1 = [x_edge(random_point(k,1),1),y_edge(random_point(k,1),1)];% point 1
p2 = [x_edge(random_point(k,2),1),y_edge(random_point(k,2),1)];% point 2
p3 = [x_edge(random_point(k,3),1),y_edge(random_point(k,3),1)];% point 3
try
[pc,r]=points2circle(p1,p2,p3);
Distance_center_edge = sqrt((x_edge-pc(1)).^2+(y_edge-pc(2)).^2);
r_residual = abs(Distance_center_edge-r);
count = length(find(r_residual<=r_error_th));
accuracy = count/N_points;
var = sqrt(sum(r_residual.^2)/(N_points-1));
metric(k,:) = [r,pc(1),pc(2),accuracy,var];
catch
continue;
end
end
%% step2: denoise
Accuracy = metric(:,4);
[row_accuracy_max,~,] = find(Accuracy==max(Accuracy));
r0 = metric(row_accuracy_max,1);
sigma_r = 5*sqrt(2);
x0 = metric(row_accuracy_max,2);
y0 = metric(row_accuracy_max,3);5 综合性解决方案
下面来看一个实际工程问题，如下提取背景复杂，存在遮挡，且圆心在图像外的圆弧的几何参数。例如，如果图像如图3所示，想要计算白色光斑的参数，应该怎么办？
看到这张图，你会发现，上述的所有办法都失效了！那咋办呢，答案是取各家之长，给出一种综合性的方法。仔细分析你会发现上述方法的优缺点如下：
- 公式法，优点是快，缺点是实际问题中没有那么好的数据。
- RANSAC算法，优点是较快，但容易受到噪声干扰，且选择三个点不合适时，会出现较大的问题。
- 最小二乘法，优点是综合考虑了所有边缘点，能够给出一个全局兼顾的解，但是受到噪声的干扰非常严重。
融合各家之长，给出的解决方案图下：
- 图像预处理。包括二值化，移除背景，边缘检测和边缘点坐标记录。
- 图像去噪。利用RANSAC算法计算圆参数，并去掉偏离很大的点。
- 利用最小二乘法拟合剩余的点，输出圆的参数。
代码如下：
clc;
close all;
clear all;
path = 'G:\pdcFPM\DATA\USAF1951\00\17.tif';
spot = imread(path);
bw_spot = imbinarize(spot); % change the format of sopt image to binary
edge_spot = edge(bw_spot,'canny');% detect the edge of spot
% figure; imshow(edge_spot,[]);
[y_edge,x_edge] = find(edge_spot);% coordinates of edge points
N_points = length(y_edge);% number of edge points
K = 60;% iterated times
r_error_th = 10;% the error of estimated r
random_point = randi(N_points,K,3);% generate K rows, 3 columns random ints whose maxmum is N_point
metric = zeros(K,5);
%% step1: estimate initial r,xc,yc
% tic
for k = 1:K
p1 = [x_edge(random_point(k,1),1),y_edge(random_point(k,1),1)];% point 1
p2 = [x_edge(random_point(k,2),1),y_edge(random_point(k,2),1)];% point 2
p3 = [x_edge(random_point(k,3),1),y_edge(random_point(k,3),1)];% point 3
try
[pc,r]=points2circle(p1,p2,p3);
Distance_center_edge = sqrt((x_edge-pc(1)).^2+(y_edge-pc(2)).^2);
r_residual = abs(Distance_center_edge-r);
count = length(find(r_residual<=r_error_th));
accuracy = count/N_points;
var = sqrt(sum(r_residual.^2)/(N_points-1));
metric(k,:) = [r,pc(1),pc(2),accuracy,var];
catch
continue;
end
end
%% step2: denoise
Accuracy = metric(:,4);
[row_accuracy_max,~,] = find(Accuracy==max(Accuracy));
r0 = metric(row_accuracy_max,1);
sigma_r = 5*sqrt(2);
x0 = metric(row_accuracy_max,2);
y0 = metric(row_accuracy_max,3);
x_new_edge_points = x_edge;
y_new_edge_points = y_edge;
for i = 1:N_points
x = x_edge(i);
y = y_edge(i);
if abs((mean(sqrt((x-x0).^2+(y-y0).^2))-r0))> 3*sigma_r
x_new_edge_points(i)=0;
y_new_edge_points(i)=0;
end
end
mask_signal = find(x_new_edge_points~=0);
x_new_edge_points = x_new_edge_points(mask_signal);
y_new_edge_points = y_new_edge_points(mask_signal);
% scatter(x_new_edge_points,y_new_edge_points);
%% use least square regression algorithm to fit parameters of circle
[para] = Circle_Fitting([x_new_edge_points,y_new_edge_points]);
xc = para(1);
yc = para(2);
radi = para(3);
disp([xc yc radi]);
imshow(bw_spot,[]);
viscircles([xc,yc],radi);运行结果如下图5所示：
可以看到，综合性的解决方案的效果杠杠的！
102
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
