技術標籤：java 根據圓心計算圓弧上點的經緯度matlab比較整體插值和分段插值圓弧裁剪演算法c++已知3個座標點xy畫圓弧

介紹一種圓弧插值方法：

這是實際應用中有人需要這個功能而設計的。實際問題是：已知道路轉彎處三個點，擬合出這個道路的彎。

實際上好多軟體都有這個演算法的。我這裡班門弄斧了。

1、計算圓心座標

根據圓上任意點到圓心的距離相等，列方程；聯立方程組，求解圓心座標(這裡提一下，matlab可以直接解方程)程式碼如下：

%{本函式用於求解平面圓的圓心座標和半徑求解思路：構造圓的方程，解方程組(圓上任一點到圓心的距離相等)輸入：圓上三個點的座標輸出：圓的圓心座標(x0,y0)和半徑r%}function [x0,y0,r] = FindParametersOfACircle(pointsTab)% 定義圓心座標syms x0 y0% 各點到圓心的距離(即半徑r)s1=sqrt((pointsTab(1,1)-x0)^2+(pointsTab(1,2)-y0)^2);s2=sqrt((pointsTab(2,1)-x0)^2+(pointsTab(2,2)-y0)^2); s3=sqrt((pointsTab(3,1)-x0)^2+(pointsTab(3,2)-y0)^2); % 構造方程組eq1=s1-s2;eq2=s2-s3; % 解方程組，得圓心座標[x0,y0]=solve(eq1,eq2);% 輸出待求引數x0 = double(x0);y0 = double(y0);r = sqrt((pointsTab(1,1) - x0)^2+(pointsTab(1,2) - y0)^2);
 

2、計算方位角

插值過程需要計算圓上點的座標。但是問題是：圓是方程，不是函式，必須分段處理，不夠方便。所以將圓引數化，那麼，插值範圍就不再是x或y的最大最小值了，而是方位角的最大最小值。這裡給出已知兩點求方位角的程式碼：

%由兩點計算方位角function [fw] = getFWfrom2points(x0,y0,x1,y1)detx = x1-x0;dety = y1-y0;if detx == 0    if dety >= 0        fw = pi()/2;    else        fw = pi()*3/2;    endelse    if detx > 0 && dety >= 0        fw = atan(dety/detx);    elseif detx < 0 && dety >=0        fw = pi() - atan(-dety/detx);    elseif detx < 0 && dety < 0        fw = pi() + atan(dety/detx);    else        fw = 2*pi() - atan(-dety/detx);    endend
 

3、插值

互動輸入插值點的數量，確定點的圓滑程度。

確定步長，遍歷計算即可。

%{圓弧加密函式輸入：加密區間，即圓弧的起點和終點座標     圓的引數，即圓心座標和半徑     加密點個數num輸出：加密點座標%}function  [pointsTab] = ArcEncryption(x1,y1,x2,y2,x0,y0,r,num)[alpha0] = getFWfrom2points(x0,y0,x1,y1);[alpha1] = getFWfrom2points(x0,y0,x2,y2);detalpha = (alpha1 - alpha0)/num;jsq = 0;for i = alpha0:detalpha:alpha1    jsq = jsq + 1;    asd = cos(i);    asdasd = sin(i);    pointsTab(jsq,1) = r*cos(i) + x0;    pointsTab(jsq,2) = r*sin(i) + y0;end
 

4、主函式如下：

%% main函式clc;clear;close all;% pointsTab裡是圓弧上三個點的座標[x1 y1;x2 y2;x3 y3]% 實際應用中，將這三個點替換為自己需要輸入的點即可pointsTab = [564675.019 3467079.061;564841.178 3467117.161;564960.770 3467217.703];A(:,1) = pointsTab(:,1) - min(pointsTab(:,1));A(:,2) = pointsTab(:,2) - min(pointsTab(:,2));[x0,y0,r] = FindParametersOfACircle(A);% num為插值點的個數% 實際應用中，需要插值多少個點就將數字20改成相應的數值num = 20;[B] = ArcEncryption(A(1,1),A(1,2),A(3,1),A(3,2),x0,y0,r,20);% C矩陣中儲存了插值點的座標% 第一列為x座標，第二列為y座標C(:,1) = B(:,1) + min(pointsTab(:,1));C(:,2) = B(:,2) + min(pointsTab(:,2));plot(C(:,1),C(:,2));hold on;plot(pointsTab(:,1),pointsTab(:,2));axis equal;msgbox '';

結果如下：