matlab練習程式（多圓交點）

Dsp Tian發表於2014-10-10
Matlab

最近總是對計算幾何方面的程式比較感興趣。

多圓求交點，要先對圓兩兩求交點。

有交點的圓分為相切圓和相交圓。

相切圓求法：

　　1.根據兩圓心求直線

　　2.求公共弦直線方程

　　3.求兩直線交點即兩圓切點。

相交圓求法：

　　1.求公共弦方程直線。

　　2.公共弦直線方程和其中一個圓方程聯立求解即可。

公共弦直線方程就是兩圓方程的差。

結果如下:

matlab程式碼如下：

main.m:

clear all;close all;clc;

n=20;
cic=rand(n,3); %（x,y,r）

hold on;
for i=1:n-1
    for j=i+1:n
        cic1=cic(i,:);
        cic2=cic(j,:);
        p=circleCross(cic1,cic2);
        if ~isempty(p)
            plot(p(:,1),p(:,2),'.');
        end
    end
end

for i=1:n
    theta=0:0.001:2*pi;
    x=cic(i,1)+cic(i,3)*cos(theta);
    y=cic(i,2)+cic(i,3)*sin(theta);
    plot(x,y,'-');
end
axis equal

circleCross.m:

function p=circleCross(cic1,cic2)

    x0=cic1(1);
    y0=cic1(2);
    r0=cic1(3);

    x1=cic2(1);
    y1=cic2(2);
    r1=cic2(3);

    d=sqrt((x0-x1)^2+(y0-y1)^2);    %兩圓心距離

    k1=(y0-y1)/(x0-x1);         %連線兩圓心直線
    b1=y1-k1*x1;

    k2=-1/k1;               %公共弦方程直線
    b2=(r0^2-r1^2-x0^2+x1^2-y0^2+y1^2)/(2*(y1-y0));

    p=[];
    if d==abs(r1-r0) || d==r1+r0        %相切時的交點
        xx=-(b1-b2)/(k1-k2);
        yy=-(-b2*k1+b1*k2)/(k1-k2);
        p=[xx yy];
    elseif abs(r1-r0)<d && d<r1+r0      %相交時的交點
                                        %公共弦方程與其中一個圓的交點
        xx1=(-b2* k2 + x1 + k2 *y1 - sqrt(-b2^2 + r1^2 + k2^2 *r1^2 - 2 *b2* k2* x1 - k2^2* x1^2 + 2*b2*y1 + 2*k2*x1*y1 - y1^2))/(1 + k2^2);
        yy1=k2*xx1+b2;
     
        xx2=(-b2* k2 + x1 + k2 *y1 + sqrt(-b2^2 + r1^2 + k2^2 *r1^2 - 2 *b2* k2* x1 - k2^2* x1^2 + 2*b2*y1 + 2*k2*x1*y1 - y1^2))/(1 + k2^2);  
        yy2=k2*xx2+b2;
        
        p=[xx1 yy1;xx2 yy2];         
    end

end

多圓求交點我只能兩兩比較了，不知道有沒有什麼快速的方法。

相關文章