Matlab交互式生成特定曲線

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);

鼠標拾取即時畫圖：

生成的數據重新畫圖：

自己的工程導入後效果：

生成的數據文件：