% 已知圆心和半径画圆：参考：https://blog.csdn.net/ZLK961543260/article/details/70216089
 
% 对比三种画圆方法，运算时长对比如下
 
% viscircle:0.22;比较快
 
% rectangle:0.21;比较快
 
% function:0.39:慢
 
[Type Sheet Format]=xlsfinfo('现状OD数据及其他数据.xls');
 
OD=xlsread('现状OD数据及其他数据.xls',Sheet{1});
 
center_area=xlsread('现状OD数据及其他数据.xls',Sheet{2});
 
x_pos=center_area(:,4); %数据第四列为x坐标(米)
 
y_pos=center_area(:,5); %数据第五列为y坐标(米)
 
area=center_area(:,3); %数据第3列为面积(平方米)：注意与x,y的坐标对应
 
r_all=sqrt(area/pi); %计算节点覆盖半径，为后文画圆做准备
 
nanbool=isnan(r_all); %数据1到4行为四个物流园区，通过r为nan，提取出四个物流园区
 
first_node_num=length(find(nanbool==1)); %统计物流园区个数
 
r_index=find(nanbool==0,1); %节点开始的行数index
 
r=r_all(r_index:end); %r存储节点的覆盖半径
 
jam_coef=center_area(r_index:end,6); %数据第六列为交通拥堵系数，只有节点才有,具体在散点图中的标注用text实现，
 
% 具体标注可参考https://zhidao.baidu.com/question/924100748051179779.html
 
% 散点图颜色可参考：https://zhidao.baidu.com/question/521200141403285445.html?qbl=relate_question_1&word=matlabscatter%B1%EA%D7%A2
 
% 具体标注和颜色还未完善
 
% 方法一:viscircle
 
%用viscircle画节点的覆盖范围，用scatter画物流园区和节点
 
x_2_pos=x_pos(r_index:end); %节点x
 
y_2_pos=y_pos(r_index:end); %节点y
 
centers=[x_2_pos,y_2_pos]; %节点圆心位置
 
tic
 
fig1=figure
 
colors = {'b','r','g','y','k'};
 
%定义colors,后文画图方便，也可以只是掉，直接在后文加入'k','g'等，具体可查询matlab画图颜色标记表示
 
viscircles(centers,r,'color',colors{1}); %matlab自带函数：已知圆心和半径画圆
 
hold on
 
scatter(x_2_pos,y_2_pos,'filled');
 
%已知x,y坐标，画散点图命令scatter(x,y),filled表示圆是实心填充
 
hold on %画多个圆时用hold on命令
 
scatter(x_pos(1:4),y_pos(1:4),'filled'); %数字4指的是物流园区的个数，均可以替换成first_node_num
 
%scatter(x_pos(1:first_node_num),y_pos(1:first_node_num),'filled');
 
toc
 
% % 方法二：rectangle
 
% % 用矩形函数rectangle 画圆,rectangle函数一次只能画一个矩形/圆；需要用到循环，与viscircle比较计算时长，tic,toc
 
tic
 
fig2=figure
 
%矩形函数画多个圆只能用for命令进行循环，用矩阵读取会发生错误
 
for i=1:length(r)
 
rectangle('Position',[x_2_pos(i)-r(i),y_2_pos(i)-r(i),2*r(i),2*r(i)],'Curvature',[1,1],'linewidth',1,'EdgeColor','b') %rectangle('position',[x_pos,y_pos,length,width])
 
end %[1,1]表示构造圆/椭圆
 
hold on
 
scatter(x_2_pos,y_2_pos,'filled');
 
hold on
 
scatter(x_pos(1:4),y_pos(1:4),'filled');
 
%scatter(x_pos(1:first_node_num),y_pos(1:first_node_num),'filled');
 
toc
 
% 方法三：构建function
 
tic
 
fig3=figure
 
for j=1:length(r)
 
x=x_2_pos(j);
 
y=y_2_pos(j);
 
r1=r(j);
 
% function [] = plot1( x,y,r1 )
 
theta=0:0.1:2*pi;
 
Circle1=x+r1*cos(theta);
 
Circle2=y+r1*sin(theta);
 
c=[123,14,52]; %color:BGR? RBG？help
 
plot(Circle1,Circle2,'c','linewidth',1);
 
hold on
 
end
 
scatter(x_2_pos,y_2_pos,'filled');
 
hold on
 
scatter(x_pos(1:4),y_pos(1:4),'filled');
 
toc
 
绘图结果：
 
figure1:
 

 
figure2:
 

 
figure3:
 

 
赛题具体地图：
 

 
参考资料：
 
1.https://blog.csdn.net/ZLK961543260/article/details/70216089
 
2.数据来源为2017年研究生数学建模竞赛F题赛题数据
 
来源：https://www.cnblogs.com/feynmania/p/12970209.html

 
MATLAB找圆心之路
 
项目场景：
问题描述：
原因分析：
解决方案：
特征提取：
存在的问题：

 
项目场景：
 
这几天在学习matlab的过程之中在处理图像的过程中遇到了一些问题，如何实现在单一背景之下，对圆心坐标准确的提取。
 
 
 
问题描述：
 
如下图所示如何动态返回圆心坐标：
 
 
 
原因分析：
 
目标颜色与背景色有着较大的差异，我们可以利用这一点，来将圆提取出来，但是其中的一些干扰与误判较难处理。
 
 
解决方案：
 
利用MATLAB解决
 
上位机图像接受
 
首先我们需要在MATLAB中安装摄像头的驱动，这个只需要在附加资源管理器中安装即可
 
 
安装完成之后通过下面代码来测试是否驱动安装成功
 
	videoinput('winvideo')
	win_info=imaqhwinfo('winvideo');
	win_info.DeviceIDs
 
 
动态图像提取
 
   % 查看摄像头支持的图像格式
   d_win_info=imaqhwinfo('winvideo',1);
   d_win_info.SupportedFormats
 
 
在这里我们选用 YUY2_320x180 格式
 
   % 生成对象并同步画面
   obj = videoinput('winvideo',1,'YUY2_320x180');%,
   h = preview(obj);
 
 
至此我们可以把图像提取出来，接下来对图像进行动态处理即可。
 
特征提取：
 
通过遍历图像的每个像素点完成
 
在我们得到图像之后，我们便可以通过一系列操作对图像进行一定的处理，得到自己想要的要的（如下图中运球的位置便是自己想要的）
 
 这里的提取思路也比较简单，首先通过选定特定的阈值，得到二值化图像；对二值化图像进行边缘化提取，所提取到的图像再进行填充、降噪等处理；通过判断填充部分的面积经行圆的判断。
 
 在程序后面可以实时输出圆心的坐标
 
存在的问题：
 
无法提取残缺圆只可以提取边缘闭合的圆

 
clc
 
clear all
 
I=imread('a003.bmp');
 
[m,n,p]=size(I);
 
I=im2bw(I);
 
se=strel('square',3);
 
Ia=imerode(I,se);
 
Iout=I-Ia;
 
figure(1);
 
imshow(Iout);
 
Iout1=Iout;
 
N=1;
 
for i=1:m
 
for j=1:n
 
if Iout(i,j)==1
 
p(N,1)=i;
 
p(N,2)=j;
 
N=N+1;
 
Iout(i,j+5:end)=0;
 
end
 
end
 
end
 
figure(2);
 
imshow(Iout);
 
for i=1:m
 
for j=n:(-1):1
 
if Iout1(i,j)==1
 
p(N,1)=i;
 
p(N,2)=j;
 
N=N+1;
 
Iout1(i,1:j-5)=0;
 
end
 
end
 
end
 
figure(3);
 
imshow(Iout1);
 
x=p(:,1);
 
y=p(:,2);
 
[R,x0,y0]=circ(x,y,N-1);
 
disp(['这个圆的半径是: ' num2str(R),' mm' ]);
 
disp(['这个圆圆心坐标:( ' num2str(x0) ,',',num2str(y0),')']);
 
figure(4)
 
imshow(I);
 
hold on
 
plot(x0,y0,'r+');
 
%%%%%%%%%%%%%%%%%%%%%%%
 
function [R,A,B]=circ(x,y,N) x1 = 0; x2 = 0; x3 = 0; y1 = 0; y2 = 0; y3 = 0; x1y1 = 0; x1y2 = 0; x2y1 = 0; for i = 1 : N x1 = x1 + x(i); x2 = x2 + x(i)*x(i); x3 = x3 + x(i)*x(i)*x(i); y1 = y1 + y(i); y2 = y2 + y(i)*y(i); y3 = y3 + y(i)*y(i)*y(i); x1y1 = x1y1 + x(i)*y(i); x1y2 = x1y2 + x(i)*y(i)*y(i); x2y1 = x2y1 + x(i)*x(i)*y(i); end C = N * x2 - x1 * x1; D = N * x1y1 - x1 * y1; E = N * x3 + N * x1y2 - (x2 + y2) * x1; G = N * y2 - y1 * y1; H = N * x2y1 + N * y3 - (x2 + y2) * y1; a = (H * D - E * G)/(C * G - D * D); b = (H * C - E * D)/(D * D - G * C); c = -(a * x1 + b * y1 + x2 + y2)/N; A = a/(-2); %x 坐标 B = b/(-2); %y 坐标 R = sqrt(a * a + b * b - 4 * c)/2;

 
function result=myfitcicledemo1();
 
y=[223.955000000000,219.390000000000,215.211000000000,211.554000000000,209.053000000000,207.667000000000,207.669000000000,208.964000000000,211.374000000000,215.773000000000,220.426000000000,225.123000000000,230.133000000000,234.219000000000,237.983000000000,241.249000000000,243.601000000000,245.165000000000,245.525000000000,244.674000000000,242.762000000000,240.792000000000,238.179000000000,233.974000000000,229.993000000000,226.155000000000,222.350000000000,218.950000000000,215.010000000000,211.693000000000];
 
x=[142.904000000000,141.692000000000,139.350000000000,135.785000000000,131.562000000000,126.453000000000,121.702000000000,116.802000000000,112.593000000000,108.394000000000,106.064000000000,105.106000000000,105.401000000000,106.680000000000,108.899000000000,112.040000000000,115.709000000000,120.305000000000,124.867000000000,129.748000000000,133.972000000000,136.639000000000,139.097000000000,141.563000000000,142.706000000000,143.078000000000,142.614000000000,141.508000000000,139.195000000000,135.963000000000];
 
[x0,y0,R0,a]=circfit5(x,y);
 
x0
 
y0
 
R0
 
alpha=0:pi/20:2*pi;
 
xx=R0*cos(alpha)+x0;
 
yy=R0*sin(alpha)+y0;
 
plot(y,x,'b+',yy,xx,'r-');
 
p=sqrt((x-x0).^2+(y-y0).^2)-R0
 
function [xc,yc,R,a] = circfit5(x,y)
 
n=length(x); xx=x.*x; yy=y.*y; xy=x.*y;
 
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
 
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
 
a=A\B;
 
xc = -.5*a(1);
 
yc = -.5*a(2);
 
R = sqrt((a(1)^2+a(2)^2)/4-a(3));A