在MATLAB中查找两个向量的交点

Question

在MATLAB中查找两个向量的交点

mathmatlab

13

我有一个非常简单的MATLAB问题。找到两个向量的交点最简单的方法是什么？我不熟悉各种MATLAB函数--看起来应该有一个函数可以实现这个功能。

例如，如果我有一个从(0,0)到(6,6)的向量和另一个从(0,6)到(6,0)的向量，我需要确定它们相交于(3,3)。

- Mav3rick

1

你应该在mathoverload.com上询问这个问题。 - Michael Bray

2

@Michael：我想你指的是mathoverflow.net，尽管该网站更偏向于“研究级数学”。像这样更简单的问题可能应该留在SO上（请参见此Meta帖子：http://meta.stackexchange.com/questions/34570/mathoverflow-net-how-can-we-get-members-of-this-site-and-the-so-sf-su-sites-to-s）。 - gnovice

此外，这个问题涉及将解决方案编程到Matlab中，而不是如何数学上解决问题。请参阅我的答案，了解如何在Matlab中通常（针对任何维度的任何数据）解决此问题。 - jdmichal

1

向量没有交点，直线才有！ - Dario

4个回答

8

好的，您有两个不同线上的点，想要找到它们的交点。最简单的方法是找出这两条线的方程，然后计算它们的交点。

一条直线的方程可以表示为y = mx + b，其中m是斜率，b是y轴截距。对于一条线，您有两个点，可以得到两个方程式。因此，您可以解出常量m和b。这将得到以下两个方程：

 0 = 0*m + 1*b  % Using the first point x=y=0 into y=m*x+b
 6 = 6*m + 1*b  % Using the second point x=y=6

或者以矩阵形式表示：

 [ 0 ] = [ 0 1 ]* [ m ]
 [ 6 ]   [ 6 1 ]  [ b ]

对于第一行，可以在MATLAB中计算常数：

 C1 = inv([0 1;6 1]*[1;0]; % m=C1(1) and b=C(2)

现在你已经得到了两条线的方程，可以通过解决以下方程组（这些方程是通过操纵一条线的方程得到的）来解决交点：

 m_1*x-y = -b_1
 m_2*x-y = -b_2

现在我们只需将上述方程组写成矩阵形式并解决：

 [x] = inv [m_1 -1] * [-b_1]
 [y]       [m_2 -1]   [-b_2]

或者用MATLAB语法表示：

 I = inv([m_1 -1; m_2 -1])*[-b_1;-b_2]; % I is the intersection.

注意事项

根据gnovice的评论，如果这些线实际上是线段，则需要检查交点是否在线段的端点之间。
如果两条直线的斜率相等，m_1 = m_2，则要么没有交点，要么有无限多个交点。

- Azim J

3

另外需要说明的一点是：如果这两条线被视为线段，则需要额外检查交点是否位于每条线段的端点之间。 - gnovice

@Amro：你能解释一下为什么应该避免使用_inv_吗？ - Azim J

4

对于方程AX=B，如果A是方阵并且可逆，那么理论上X = inv(A)*B与X = A\B是相同的。但使用反斜杠操作符进行计算更好，因为它们需要更少的计算机时间，更少的内存，并具有更好的错误检测性能。请参考http://www.mathworks.com/access/helpdesk/help/techdoc/math/f4-983672.html和http://www.mathworks.com/access/helpdesk/help/techdoc/math/f4-2224.html获取更多解释。 - Amro

2

对于一般的多维解决方案，实际上你正在解决一系列线性系统。

首先，你需要将方程转化为线性形式：Ax+By=C（根据需要扩展维度）

对于两点：

y - y1 = (y2 - y1) / (x2 - x1) * (x - x1)
y - y1 = (y2 - y1) / (x2 - x1) * x - (y2 - y1) / (x2 - x1) * x1
(y1 - y2) / (x2 - x1) * x + y - y1 = (y1 - y2) / (x2 - x1) * x1
(y1 - y2) / (x2 - x1) * x + y = (y1 - y2) / (x2 - x1) * x1 + y1
(y1 - y2) * x + (x2 - x1) * y = (y1 - y2) * x1 + (x2 - x1) * y1

A = (y1 - y2)
B = (x2 - x1)
C = (y1 - y2) * x1 + (x2 - x1) * y1 = A * x1 + B * y1

举个例子：

x1 = 0, x2 = 6, y1 = 0, y2 = 6
A1 = (0 - 6) = -6
B1 = (6 - 0) = 6
C1 = A1 * 0 + B1 * 0 = 0

x1 = 0, x2 = 6, y1 = 6, y2 = 0
A2 = (6 - 0) = 6
B2 = (6 - 0) = 6
C2 = A2 * 0 + B2 * 6 = 6 * 6 = 36

然后，形成一个矩阵，A、B和C在行中：

[A1 B1 C1]
[A2 B2 C2]

[-6 6 0]
[ 6 6 36]

现在使用Matlab函数rref(matrix)将矩阵化简为阶梯形式：

[ 1 0 3]
[ 0 1 3]

正如你所猜测的那样，最后一列是你的交点。这可以扩展到任意多个维度。如果你的简化阶梯形式除了身份矩阵之外还有其他内容，则你的向量要么没有唯一的交点，要么根据矩阵的形式没有交点。

dim = 2;

% Do other stuff, ending with rref(matrix)

if (matrix(:,1:dim) == eye(dim))
    % Matrix has unique solution.
    solution = (matrix(:,dim+1))'
else
    % No unique solution.
end

在二维中，变化的形式有：

线性解，表示解是形如x + By = C的一条直线：

[ 1 B C]
[ 0 0 0]

无解，表示这些线不相交，其中 C2 <> 0：

[ 1 B C1]
[ 0 0 C2]

- jdmichal

2

其他的结果在我看来比较混乱、啰嗦且不完整。以下是我的两分钱，也可能会让人感到困惑和啰嗦。

如果你确定你的线不是斜交或平行的，那么下面的内容就是你所需要的：

% Let each point be def as a 3x1 array
% Let points defining first line be  : p1, q1
% Let points defining second line be : p2, q2

L = p1-p2;
M = p1-q1;
N = p2-q2;
A = [M N];
T = pinv(A)*L;
h = p1-T(1)*(p1-q1); % h is a 3x1 array representing the actual pt of intersection

是的，Moore-Penrose伪逆是一种强大的工具。该方法的解释是：您想要找到“方向向量”的权重或缩放因子（M和N是方向向量），它们线性组合M和N以给出L。

下面是完整的描述。它提供了一个简单的异常检测方案，并将其处理留给用户。（两条直线算法之间的最小距离来自维基百科；比较方向余弦（DCS）以检查向量态度是常识。）

% Let each point be def as a 3x1 array
% Let points defining first line be : p1, q1
% Let points defining second line be: p2, q2

% There are two conditions that prevent intersection of line segments/lines
% in L3 space. 1. parallel 2. skew-parallel (two lines on parallel planes do not intersect)
% Both conditions need to be identified and handled in a general algorithm.

% First check that lines are not parallel, this is done by comparing DCS of
% the line vectors
% L, M, N ARE DIRECTION VECTORS.
L = p1-p2;
M = p1-q1;
N = p2-q2;

% Calculate a normalized DCS for comparison. If equal, it means lines are parallel.
MVectorMagnitude = sqrt(sum(M.*M,2)); % The rowsum is just a generalization for N-D vectors.
NVectorMagnitude=sqrt(sum(N.*N,2)); % The rowsum is just a generalization for N-D vectors.

if isequal(M/MVectorMagnitude,N/NVectorMagnitude) % Compare the DCS for equality
     fprintf('%s\n', 'lines are parallel. End routine')
end;

% Now check that lines do not exist on parallel planes
% This is done by checking the minimum distance between the two lines. If there's a minimum distance, then the lines are skew.

a1 = dot(M,L); b1 = dot(M,M); c1 = dot(M,N);
a2 = dot(N,L); b2 = dot(N,M); c2 = dot(N,N);

s1 = -(a1*c2 - a2*c1)/(b1*c2-b2*c1);
s2 = -(a1*b2 - a2*b1)/(b1*c2-b2*c1);

Sm = (L + s1*M - s2*N);
s = sqrt(sum(Sm.*Sm,2));

if ~isequal(s,0) % If the minimum distance between two lines is not zero, then the lines do not intersect
    fprintf('%s\n','lines are skew. End routine')
end;

% Here's the actual calculation of the point of intersection of two lines.
A = [M N];
T = pinv(A)*L;
h = p1-T(1)*(p1-q1); % h is a 3x1 array representing the actual pt of intersection.

因此，使用伪逆方法即使当您的M和N向量不平行时也会给出结果（但是不能平行，因为需要inv(A'.A)存在）。您可以使用此方法确定两条平行线或两个平面之间的最小距离-为此，请定义k = p2+T(2)*(p2-q2)，然后所需的距离为h-k。还要注意，如果线是斜的，则h和k是彼此最接近的线上的点。

因此，伪逆和投影空间的使用为我们提供了一个简洁的算法来：

确定两条线（不平行，不斜）的交点
确定两条线（不平行）之间的最小距离
确定两条斜线上彼此最接近的点。

简洁并不意味着时间效率高。很多东西取决于您的精确伪逆函数实现- MATLAB使用svd进行求解，从而达到容差。此外，在更高的维度和更高阶的测量度量（或矢量范数）定义中，某些结果只能近似准确。除了明显的与维度无关的实现外，这可以用于统计回归分析和代数最大化点估计的可能性。

- SkyRat

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- gnovice · Accepted Answer

一种解决方法是使用此教程中推导的方程式，在二维中找到两条直线的交点 (更新：这是互联网档案馆的链接，因为该网站已经不存在)。您可以先创建两个矩阵：一个用于保存线段端点的x坐标，另一个用于保存y坐标。

x = [0 0; 6 6];  %# Starting points in first row, ending points in second row
y = [0 6; 6 0];

然后可以将上述来源的方程式编码如下：

dx = diff(x);  %# Take the differences down each column
dy = diff(y);
den = dx(1)*dy(2)-dy(1)*dx(2);  %# Precompute the denominator
ua = (dx(2)*(y(1)-y(3))-dy(2)*(x(1)-x(3)))/den;
ub = (dx(1)*(y(1)-y(3))-dy(1)*(x(1)-x(3)))/den;

现在你可以计算出两条线的交点：

xi = x(1)+ua*dx(1);
yi = y(1)+ua*dy(1);

对于问题中的示例，上述代码会给出xi = 3和yi = 3，符合预期。如果您想检查交点是否在直线段的端点之间（即它们是有限直线段），只需检查ua和ub值是否都在0和1之间：

isInSegment = all(([ua ub] >= 0) & ([ua ub] <= 1));

我在上面链接的教程中还有几个要点：

如果分母den为0，则这两条线是平行的。
如果ua和ub的方程式的分子和分母都为0，则这两条线是重合的。