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: [email protected]) % % 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);
鼠標拾取即時畫圖:
生成的數據重新畫圖:
自己的工程導入後效果:
生成的數據文件: