matlab calccircle,Matlab交互式生成特定曲线

工程中，经常要产生测试数据，比如特定包含噪声的曲线。以下是Matlab实现的通过鼠标拾取坐标点，然后生成包含直线和圆弧的x,y坐标的代码。

操作方式：鼠标拾取一个起点，弹出选择菜单，选择直线或者圆弧，选择直线后，鼠标拾取第二个点；选择圆弧之后，鼠标拾取另外两个点，通过三个点拟合圆弧并生成圆弧数据。

主文件：

clc,clear,close all%

figure

axis([0 1000 0 1000])%坐标轴大小

axis square

x_out=[];y_out=[];%保存所有生成的数据

xy0=ginput(1);%起点

temp_xy=xy0;%上一次拾取的点

plot(xy0(1),xy0(2),'*')

axis([0 1000 0 1000]),hold on

axis square

n=4;%曲线的段数

for i=1:n

res=menu('选择','直线','圆弧');

if(res==1)

%disp(['选择的是直线'])

xy1=ginput(1);

[x,y]=drawLine(temp_xy,xy1);

x_out=[x_out x];y_out=[y_out y];

temp_xy=xy1(end,:);

plot(x,y);

else

%disp(['选择的是圆弧'])

%注意圆弧y坐标唯一，即扇形角度不能超过180度

xy1=ginput(2);

[x,y]=getArc(temp_xy,xy1(1,:),xy1(2,:));

x_out=[x_out x];y_out=[y_out y];

temp_xy=xy1(end,:);

end

end

%加噪声

noise=2;%设定噪声幅度

for i=1:length(x_out)

y_out(i)=y_out(i)+random('Normal',-1*noise,noise);

end

xy_file=[x_out' y_out'];%转置为两列

xy_file=sortrows(xy_file,1);%按x坐标(第一列)重排序，放置圆弧产生的顺序颠倒

xy_file=[xy_file(:,1)';xy_file(:,2)']%转置回两行坐标形式

figure

plot(xy_file(1,:),xy_file(2,:))

%写入txt文件，文件名自动递增

count=1;

fileName=['D:\data' num2str(count) '.txt'];

if(exist(fileName,'file'))

count=count+1;

fileName=['D:\data' num2str(count) '.txt'];

end

fid=fopen(fileName,'wt');

fprintf(fid,'%f,%f\n',xy_file);

fclose(fid);

两点生成直线文件：

%两点产生直线数据

function [x,y]=drawLine(A,B)

x=A(1):0.1:B(1);

y=(B(2)-A(2))/(B(1)-A(1))*(x-A(1))+A(2);

end

三点生成圆弧数据文件：

%三点计算圆弧(画图)

%数据A=[1 2];

%B=[5 6];

%C=[3 5];

function [x,y]=getArc(A,B,C)

[c,r]=calcCircle(A,B,C);

a=c(1);

b=c(2);

th=[A;B;C];

th2=[th(:,1)-a th(:,2)-b];

theta=atan2(th2(:,2),th2(:,1));

[theta_max,num_max]=max(theta);

[theta_min,num_min]=min(theta);

t=linspace(theta_min,theta_max,1000);

%t=0:0.1:2*pi;

a=c(1);

b=c(2);

x=r*cos(t)+a;

y=r*sin(t)+b;

plot(x,y,'r-',th(:,1),th(:,2),'o')

%axis equal

end

圆弧拟合文件(注：该文件来源于网络)：

function [centre radius] = calcCircle(pt1, pt2, pt3)

% calcCircle: Fit a circle to a set of 3 points

%

% Inputs:

% pt1, pt2 and pt3 are vectors with 2 elements representing a point

% in 2D Cartesian coordinates.

%

% Returns:

% The centre coordinate (2 elements) and radius of the circle.

% A centre value of [0,0] and radius of -1 if the points are collinear.

%

% Example:

%

% p1 = rand(1,2);

% p2 = rand(1,2);

% p3 = rand(1,2);

%

% [c r] = calcCircle(p1, p2, p3);

%

% figure(1)

% cla

% axis equal

% hold on

% if r ~= -1

% rectangle('Position',[c(1)-r,c(2)-r,2*r,2*r],'Curvature',[1,1],'EdgeColor','g')

% end

% plot(p1(1), p1(2), '*')

% plot(p2(1), p2(2), '*')

% plot(p3(1), p3(2), '*')

%

% for Matlab R13 and up

% version 1.2 (mar 2008)

% Author: Peter Bone (email: peterbone@hotmail.com)

%

% History

% Created: 6th March 2008, version 1.1

% Revisions

% 7th March 2008: Version 1.2 for improved help and usability

% argument checking

if nargin < 3

error('Three input points are required.');

elseif ~isequal(numel(pt1),numel(pt2),numel(pt3),2)

error('The three input points should all have two elements.')

end

pt1 = double(pt1);

pt2 = double(pt2);

pt3 = double(pt3);

epsilon = 0.000000001;

delta_a = pt2 - pt1;

delta_b = pt3 - pt2;

ax_is_0 = abs(delta_a(1)) <= epsilon;

bx_is_0 = abs(delta_b(1)) <= epsilon;

% check whether both lines are vertical - collinear

if ax_is_0 && bx_is_0

centre = [0 0];

radius = -1;

warning([mfilename ':CollinearPoints'],'Points are on a straight line (collinear).');

return

end

% make sure delta gradients are not vertical

% swap points to change deltas

if ax_is_0

tmp = pt2;

pt2 = pt3;

pt3 = tmp;

delta_a = pt2 - pt1;

end

if bx_is_0

tmp = pt1;

pt1 = pt2;

pt2 = tmp;

delta_b = pt3 - pt2;

end

grad_a = delta_a(2) / delta_a(1);

grad_b = delta_b(2) / delta_b(1);

% check whether the given points are collinear

if abs(grad_a-grad_b) <= epsilon

centre = [0 0];

radius = -1;

warning([mfilename ':CollinearPoints'],'Points are on a straight line (collinear).');

return

end

% swap grads and points if grad_a is 0

if abs(grad_a) <= epsilon

tmp = grad_a;

grad_a = grad_b;

grad_b = tmp;

tmp = pt1;

pt1 = pt3;

pt3 = tmp;

end

% calculate centre - where the lines perpendicular to the centre of

% segments a and b intersect.

centre(1) = ( grad_a*grad_b*(pt1(2)-pt3(2)) + grad_b*(pt1(1)+pt2(1)) - grad_a*(pt2(1)+pt3(1)) ) / (2*(grad_b-grad_a));

centre(2) = ((pt1(1)+pt2(1))/2 - centre(1)) / grad_a + (pt1(2)+pt2(2))/2;

% calculate radius

radius = norm(centre - pt1);

部分截图：