本函式的作用是快速計算曲線的互交點與自交點，互交點是指兩條曲線之間的交點，自交點是指該曲線自己與自己的交點，如曲線“8”的中間點即為自交點。該函式編寫全部採用向量運算，沒有使用任何for迴圈，提高了程式碼的執行效率。

P = INTERX(L1,L2) 返回曲線L1與L2的交點，曲線的描述方式為2行形式的矩陣，第一行為x座標，第二行為y座標。

如果兩條曲線之間沒有交點，則P返回空矩陣。

P = INTERX(L1) 返回曲線L1的自交點，但不包括切線的交點。

示例如下：

t = linspace(0,2*pi);

r1 = sin(4*t)+2; x1 = r1.*cos(t); y1 = r1.*sin(t);

r2 = sin(8*t)+2; x2 = r2.*cos(t); y2 = r2.*sin(t);

P = InterX([x1;y1],[x2;y2]);

plot(x1,y1,x2,y2,P(1,:),P(2,:),‘ro’)

完整MATLAB程式碼如下：

function P = InterX(L1,varargin)

%INTERX Intersection of curves

% P = INTERX(L1,L2) returns the intersection points of two curves L1

% and L2. The curves L1,L2 can be either closed or open and are described

% by two-row-matrices, where each row contains its x- and y- coordinates.

% The intersection of groups of curves (e.g. contour lines, multiply

% connected regions etc) can also be computed by separating them with a

% column of NaNs as for example

%

% L = [x11 x12 x13 … NaN x21 x22 x23 …;

% y11 y12 y13 … NaN y21 y22 y23 …]

%

% P has the same structure as L1 and L2, and its rows correspond to the

% x- and y- coordinates of the intersection points of L1 and L2. If no

% intersections are found, the returned P is empty.

%

% P = INTERX(L1) returns the self-intersection points of L1. To keep

% the code simple, the points at which the curve is tangent to itself are

% not included. P = INTERX(L1,L1) returns all the points of the curve

% together with any self-intersection points.

%

% Example:

% t = linspace(0,2*pi);

% r1 = sin(4*t)+2; x1 = r1.*cos(t); y1 = r1.*sin(t);

% r2 = sin(8*t)+2; x2 = r2.*cos(t); y2 = r2.*sin(t);

% P = InterX([x1;y1],[x2;y2]);

% plot(x1,y1,x2,y2,P(1,:),P(2,:),‘ro’)

% Author : NS

% Version: 3.0, 21 Sept. 2010

% Two words about the algorithm: Most of the code is self-explanatory.

% The only trick lies in the calculation of C1 and C2. To be brief, this

% is essentially the two-dimensional analog of the condition that needs

% to be satisfied by a function F(x) that has a zero in the interval

% [a,b], namely

% F(a)*F(b) <= 0

% C1 and C2 exactly do this for each segment of curves 1 and 2

% respectively. If this condition is satisfied simultaneously for two

% segments then we know that they will cross at some point.

% Each factor of the ‘C’ arrays is essentially a matrix containing

% the numerators of the signed distances between points of one curve

% and line segments of the other.

%...Argument checks and assignment of L2

error(nargchk(1,2,nargin));

if nargin == 1,

    L2 = L1;    hF = @lt;   %...Avoid the inclusion of common points

else

    L2 = varargin{1}; hF = @le;

end

   

%...Preliminary stuff

x1  = L1(1,:)';  x2 = L2(1,:);

y1  = L1(2,:)';  y2 = L2(2,:);

dx1 = diff(x1); dy1 = diff(y1);

dx2 = diff(x2); dy2 = diff(y2);



%...Determine 'signed distances'   

S1 = dx1.*y1(1:end-1) - dy1.*x1(1:end-1);

S2 = dx2.*y2(1:end-1) - dy2.*x2(1:end-1);



C1 = feval(hF,D(bsxfun(@times,dx1,y2)-bsxfun(@times,dy1,x2),S1),0);

C2 = feval(hF,D((bsxfun(@times,y1,dx2)-bsxfun(@times,x1,dy2))',S2'),0)';



%...Obtain the segments where an intersection is expected

[i,j] = find(C1 & C2); 

if isempty(i),P = zeros(2,0);return; end;



%...Transpose and prepare for output

i=i'; dx2=dx2'; dy2=dy2'; S2 = S2';

L = dy2(j).*dx1(i) - dy1(i).*dx2(j);

i = i(L~=0); j=j(L~=0); L=L(L~=0);  %...Avoid divisions by 0



%...Solve system of eqs to get the common points

P = unique([dx2(j).*S1(i) - dx1(i).*S2(j), ...

            dy2(j).*S1(i) - dy1(i).*S2(j)]./[L L],'rows')';

          

function u = D(x,y)

    u = bsxfun(@minus,x(:,1:end-1),y).*bsxfun(@minus,x(:,2:end),y);

end

end

下載原始碼地址：

更多精彩文章請關注微訊號：