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 == 2p1=[p1,0];p2=[p2,0];p3=[p3,0];elseif num1 ~= 3error('仅支持二维或三维坐标输入');endelseerror('输入坐标的维数不一致');end% 共线检查temp01=p1-p2;temp02=p3-p2;temp03=cross(temp01,temp02);temp=(temp03*temp03')/(temp01*temp01')/(temp02*temp02');if temp < 10^-6error('三点共线, 无法确定圆');endmat1=[p1,1;p2,1;p3,1];% size = 3x4m=+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 = 4x5A=+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:Nsum_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 binaryedge_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:Kp1 = [x_edge(random_point(k,1),1),y_edge(random_point(k,1),1)];% point 1p2 = [x_edge(random_point(k,2),1),y_edge(random_point(k,2),1)];% point 2p3 = [x_edge(random_point(k,3),1),y_edge(random_point(k,3),1)];% point 3try[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];catchcontinue;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 binaryedge_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:Kp1 = [x_edge(random_point(k,1),1),y_edge(random_point(k,1),1)];% point 1p2 = [x_edge(random_point(k,2),1),y_edge(random_point(k,2),1)];% point 2p3 = [x_edge(random_point(k,3),1),y_edge(random_point(k,3),1)];% point 3try[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];catchcontinue;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_pointsx = x_edge(i);y = y_edge(i);if abs((mean(sqrt((x-x0).^2+(y-y0).^2))-r0))> 3*sigma_rx_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所示：

可以看到，综合性的解决方案的效果杠杠的！

matlab圆心提取【你想要的方法这里都有】相关推荐

matlab 提取极值,利用matlab 进行极值统计的一个例子——gev 方法.pdf
利用matlab 进行极值统计的一个例子--gev 方法利用 Matlab 进行极值统计的一个例子--GEV 方法科研菜鸟 /u/sanshiphy 数据和例子均来自于 S. Coles, An ...
matlab提取数据的一部分,matlab如何提取数组中的满足一定范围的一段数据
给定一个数组,如何让matlab生成一个这个数组中的一个随机数? x=[102030];x(randi(length(x)));其中randi(length(x))生成从1~(x的长度)这几个自然数中 ...
MATLAB 颜色提取器 —— APP 版
MATLAB 颜色提取器 -- APP 版日常设计GUI或者APP时,往往需要与颜色打交道,文章链接: MATLAB 如何画出漂亮的图. 在颜色选取上可以通过颜色对照表来选择合适的颜色,有个不方 ...
一步一步教你抓数据——用.net精确提取网站数据的通用方法 [转]
一步一步教你抓数据--用.net精确提取网站数据的通用方法 [转] 2008年02月23日星期六 16:53 具体实现思路: 1 首先用WebClient类下载网页源码 public static ...
access查询出生日期格式转换_从身份证中提取出生日期的3个方法和计算年龄和星座的方法...
在我们日常的工作当中,经常会遇到通过身份证来获取出生年月日的需求,今天就给大家介绍三种可以从身份证中提取出生年月日的方法. 我们都知道身份证不同的区域是有不同的含义的,代表出生年月日的数字是第7位到第 ...
MATLAB求解导弹运动的一些基础方法
导弹运动方程求解数值解离不开MATLAB,本文提出一些基本的方法与技巧. 1.微分方程的求解:ode45.给出例子: function dy = lateral( t,y ) %求解偶然干扰作用下侧向 ...
Matlab中计算程序运行时间的几种方法
平常科研当中,当我们在看文献时,没看到一个优秀的算法时都有想要自己动手编程去实现的愿望,算法好坏可以用代码的运行时间来评估,在MATLAB中大致有以下几种方法来计算程序的运行时间: 1.tic和toc ...
MATLAB#183;提取图像中多个目标
基于matlab工具箱提取图像中的多目标特征(代码如下): 代码前面部分为提取图像的边界信息,调用了后面的遍历函数Pixel_Search,函数实现方法见后~ %%ROI Testing close ...
matlab 输入坐标,matlab中坐标希腊符号的输入方法
希腊字母显示 1.Word中不用公式编辑器输入希腊字母的方法首先你需要先打开一个Office文档,然后在你需要输入希腊字母的时候,关闭中文输入法或将输入法切换为英文状态,然后同时按下Ctrl+Shi ...

matlab圆心提取【你想要的方法这里都有】