【源码】曲线交点计算——基于矢量运算的快速互交点与自交点定位函数
本函数的作用是快速计算曲线的互交点与自交点，互交点是指两条曲线之间的交点，自交点是指该曲线自己与自己的交点，如曲线“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
下载源代码地址：
http://page2.dfpan.com/fs/6lcc3j52627192d9169/
更多精彩文章请关注微信号：