【源码】快速而稳健的两条曲线交点计算

此函数计算两条曲线相交的（x，y）位置。

This function computes the (x,y) locations where two curves intersect.

曲线可以用NaN断开或具有垂直分段。

The curves can be broken with NaNs or have vertical segments.

该函数的运算也非常快（至少我认为在典型应用程序的数据处理上是很快的）。

It is also very fast (at least on data that represents what I think is a typical application).

部分代码如下：

function [x0,y0,iout,jout] = intersections(x1,y1,x2,y2,robust)

%INTERSECTIONS Intersections of curves.

% Computes the (x,y) locations where two curves intersect. The curves

% can be broken with NaNs or have vertical segments.

%

% Example:

% [X0,Y0] = intersections(X1,Y1,X2,Y2,ROBUST);

%

% where X1 and Y1 are equal-length vectors of at least two points and

% represent curve 1. Similarly, X2 and Y2 represent curve 2.

% X0 and Y0 are column vectors containing the points at which the two

% curves intersect.

%

% ROBUST (optional) set to 1 or true means to use a slight variation of the

% algorithm that might return duplicates of some intersection points, and

% then remove those duplicates. The default is true, but since the

% algorithm is slightly slower you can set it to false if you know that

% your curves don’t intersect at any segment boundaries. Also, the robust

% version properly handles parallel and overlapping segments.

%

% The algorithm can return two additional vectors that indicate which

% segment pairs contain intersections and where they are:

%

% [X0,Y0,I,J] = intersections(X1,Y1,X2,Y2,ROBUST);

%

% For each element of the vector I, I(k) = (segment number of (X1,Y1)) +

% (how far along this segment the intersection is). For example, if I(k) =

% 45.25 then the intersection lies a quarter of the way between the line

% segment connecting (X1(45),Y1(45)) and (X1(46),Y1(46)). Similarly for

% the vector J and the segments in (X2,Y2).

%

% You can also get intersections of a curve with itself. Simply pass in

% only one curve, i.e.,

%

% [X0,Y0] = intersections(X1,Y1,ROBUST);

%

% where, as before, ROBUST is optional.

% Version: 2.0, 25 May 2017

% Author: Douglas M. Schwarz

% Email: dmschwarz=ieeeorg, dmschwarz=urgradrochester*edu

% Real_email = regexprep(Email,{’=’,’*’},{’@’,’.’})

% Theory of operation:

%

% Given two line segments, L1 and L2,

%

% L1 endpoints: (x1(1),y1(1)) and (x1(2),y1(2))

% L2 endpoints: (x2(1),y2(1)) and (x2(2),y2(2))

%

% we can write four equations with four unknowns and then solve them. The

% four unknowns are t1, t2, x0 and y0, where (x0,y0) is the intersection of

% L1 and L2, t1 is the distance from the starting point of L1 to the

% intersection relative to the length of L1 and t2 is the distance from the

% starting point of L2 to the intersection relative to the length of L2.

%

% So, the four equations are

%

% (x1(2) - x1(1))*t1 = x0 - x1(1)

% (x2(2) - x2(1))*t2 = x0 - x2(1)

% (y1(2) - y1(1))*t1 = y0 - y1(1)

% (y2(2) - y2(1))*t2 = y0 - y2(1)

%

% Rearranging and writing in matrix form,

%

% [x1(2)-x1(1) 0 -1 0; [t1; [-x1(1);

% 0 x2(2)-x2(1) -1 0; * t2; = -x2(1);

% y1(2)-y1(1) 0 0 -1; x0; -y1(1);

% 0 y2(2)-y2(1) 0 -1] y0] -y2(1)]

%

% Let’s call that A*T = B. We can solve for T with T = A\B.

%

% Once we have our solution we just have to look at t1 and t2 to determine

% whether L1 and L2 intersect. If 0 <= t1 < 1 and 0 <= t2 < 1 then the two

% line segments cross and we can include (x0,y0) in the output.

%

% In principle, we have to perform this computation on every pair of line

% segments in the input data. This can be quite a large number of pairs so

% we will reduce it by doing a simple preliminary check to eliminate line

% segment pairs that could not possibly cross. The check is to look at the

% smallest enclosing rectangles (with sides parallel to the axes) for each

% line segment pair and see if they overlap. If they do then we have to

% compute t1 and t2 (via the A\B computation) to see if the line segments

% cross, but if they don’t then the line segments cannot cross. In a

% typical application, this technique will eliminate most of the potential

% line segment pairs.

% Input checks.

if verLessThan(‘matlab’,‘7.13’)

error(nargchk(2,5,nargin)) %#ok<NCHKN>

else

narginchk(2,5)

end

…

源码下载地址：

http://page2.dfpan.com/fs/6lcaj2021729416d2b4/

更多精彩文章请关注微信号：