四边形内一点的相对位置
Relative position of a point within a quadrilateral
我正在尝试找到最简单的方法来确定一个点在四边形中的相对位置。已知的是(见图)点1、2、3、4、5在xy坐标系中的位置:x1,y1,x2,y2,x3,y3,x4,y4,x5,y5.
也已知点1、2、3、4在ξ-η坐标系中的位置(见图)
根据这些数据,我想确定点 5 的 ξ 和 η 是多少。
结果
感谢所有回答的人!我发现@dbc 和@agentp 的解决方案相似。我还发现这个解决方案比@MBo 的透视变换解决方案更好,因为我不必计算矩阵的逆 (Ax=B --> x=inv(A)*B).
我得到以下结果:
u = 0.5 * (ξ + 1)
v = 0.5 * (η + 1)
在我的例子中,所有点都在矩形内,因此 u>0 和 v>0。
您需要计算一个透视变换矩阵,将源四边形的 4 个点映射到目标四边形的 4 个点 (example) (more mathemathics),然后将此变换应用于第 5 个点的坐标(将矩阵乘以坐标向量)
这看起来像是标准的有限元参数化
(这个问题没有具体cify 特定映射,但我想有人可能对这个 specific 案例感兴趣)
{x, y} == (
(1 - eta) (1 - ci) {p1x, p1y} +
(1 - eta) (1 + ci) {p2x, p2y} +
(1 + eta) (1 + ci) {p3x, p3y} +
(1 + eta) (1 - ci) {p4x, p4y} )/4
这可以用 {eta,ci} 的封闭形式求解,但是表达式对于 post.
来说相当笨拙
在实践中,计算这些常量:
ax = p1x + p2x + p3x + p4x
bx = p1x - p2x - p3x + p4x
cx = p1x + p2x - p3x - p4x
dx = p1x - p2x + p3x - p4x
ay = p1y + p2y + p3y + p4y
by = p1y - p2y - p3y + p4y
cy = p1y + p2y - p3y - p4y;
dy = p1y - p2y + p3y - p4y;
求解此二次方程 eta
:
(ax by - bx ay) - 4 (by x - bx y) +
eta (dx ay - cx by + bx cy - ax dy + 4 (x dy - dx y)) +
eta^2 (cx dy - dx cy) == 0
然后得到 ci
为:
ci = ((-ax + eta cx + 4 x)/(-bx + eta dx))
如果多边形不太扭曲,只有一种解决方案将满足 -1<eta<1
和 -1<ci<1
你这里有一个 2d bilinear blended surface。为简单起见,让我们将其坐标更改为从零到一:
u = 0.5 * (ξ + 1)
v = 0.5 * (η + 1)
在那种情况下,表面评估器可以表示为
F(u, v) = P1 + u * (P2 - P1) + v * ((P4 + u * (P3 - P4)) - (P1 + u * (P2 - P1)))
即对于给定的u
,构造一条通过以下两点的直线:
Pv0 = P1 + u * (P2 - P1);
Pv1 = P4 + u * (P3 - P4);
然后在给定的 v
之间进行插值
F(u, v) = Pv0 + v * (Pv1 - Pv0)
您要寻找的是满足 F(u, v) = P5
的值 (u,v)
。当从 Pv0
到 Pv1
的线穿过 P5
时,给定 u
会发生这种情况,这将在 P5 - Pv0
平行于 Pv1 - Pv0
时发生-- 即当它们的 2d cross 为零时:
cross2d(P5 - Pv0, Pv1 - Pv0) = 0
⇒
cross2d(P5 - (P1 + u * (P2 - P1)),
P4 + u * (P3 - P4) - (P1 + u * (P2 - P1))) = 0
现在,两个二维向量 A ⨯ B
的 2d cross 由 Ax*By - Ay*Bx
给出,因此等式变为
(x5 - (x1 + u * (x2 - x1))) * (y4 + u * (y3 - y4) - (y1 + u * (y2 - y1))) - (y5 - (y1 + u * (y2 - y1))) * (x4 + u * (x3 - x4) - (x1 + u * (x2 - x1))) = 0
Expanding this expression out and collecting collecting together terms in u
,我们得到
u^2 * (x1*y3 - x1*y4 - x2*y3 + x2*y4 + (-x3)*y1 + x3*y2 + x4*y1 - x4*y2)
+ u * (-x1*y3 + 2*x1*y4 - x1*y5 - x2*y4 + x2*y5 + x3*y1 - x3*y5 - 2*x4*y1 + x4*y2 + x4*y5 + x5*y1 - x5*y2 + x5*y3 - x5*y4)
+ (-x1*y4 + x1*y5 + x4*y1 - x4*y5 - x5*y1 + x5*y4)
= 0
现在这是 u
上的二次方程,并且可以是 solved。请注意,在四边形的顶部和底部边缘平行的情况下,二次方程将转化为线性方程;你的二次方程求解器必须处理这个问题。
double a = (x1 * y3 - x1 * y4 - x2 * y3 + x2 * y4 + (-x3) * y1 + x3 * y2 + x4 * y1 - x4 * y2);
double b = (-x1 * y3 + 2 * x1 * y4 - x1 * y5 - x2 * y4 + x2 * y5 + x3 * y1 - x3 * y5 - 2 * x4 * y1 + x4 * y2 + x4 * y5 + x5 * y1 - x5 * y2 + x5 * y3 - x5 * y4);
double c = (-x1 * y4 + x1 * y5 + x4 * y1 - x4 * y5 - x5 * y1 + x5 * y4);
double[] solutions = Quadratic.Solve(a, b, c);
解决方案可能不止一种。对于退化四边形也可能没有解决方案。
求解 u
的值后,找到等效的 v
很简单。给定积分
Pv0 = P1 + u * (P2 - P1);
Pv1 = P4 + u * (P3 - P4);
你寻求 v
这样
v * (Pv1 - Pv0) = P5 - Pv0;
选择坐标索引 0 或 1,使 |(Pv1 - Pv0)[index]|
最大化 。 (如果两个坐标几乎都为零,那么放弃——这个具体u
无解。然后设置
v = (P5 - Pv0)[index] / (Pv1 - Pv0)[index];
最后,如果您有多个解决方案,请选择混合 [u, v]
边界内的解决方案。然后最后设置
ξ = 2 * u - 1;
η = 2 * v - 1;
参考@blaz 的自我回答(请为@blaze、@dbc 和@agentp 的回答投票)
对于不愿意手工复制公式的人,这里是 C# 代码的公式:
double v_sqrt = Math.Sqrt(
4 * (
(x3 - x4) * (y1 - y2) - (x1 - x2) * (y3 - y4)) * (x4 * (-1 * y + y1) + x1 * (y - y4) + x * (-1 * y1 + y4)) +
Math.Pow(
(x3 * y - x4 * y - x3 * y1 + 2 * x4 * y1 - x4 * y2 + x1 * (y + y3 - 2 * y4) + x2 * (-1 * y + y4) + x * (-1 * y1 + y2 - y3 + y4))
, 2)
);
double u_sqrt = Math.Sqrt(
4 * ((x3 - x4) * (y1 - y2) - (x1 - x2) * (y3 - y4))
* (
x4 * (-1 * y + y1) + x1 * (y - y4) + x * (-1 * y1 + y4)
) +
Math.Pow(
(x3 * y - x4 * y - x3 * y1 + 2 * x4 * y1 - x4 * y2 + x1 * (y + y3 - 2 * y4) + x2 * (-1 * y + y4) + x * (-1 * y1 + y2 - y3 + y4))
, 2)
);
double k = 1 / (2 * ((x3 - x4) * (y1 - y2) - (x1 - x2) * (y3 - y4)));
double l = 1 / (2 * ((x1 - x4) * (y2 - y3) - (x2 - x3) * (y1 - y4)));
///////////////////////////////////////////////////////////////////////////////////////////////
double v1 = l *
(x2 * y - x3 * y + x4 * y + x * y1 - 2 * x2 * y1 + x3 * y1 - x * y2 - x4 * y2 + x * y3 - x1 * (y - 2 * y2 + y3) - x * y4 + x2 * y4 +
v_sqrt);
///////////////////////////////////////////////////////////////////////////////////////////////
double u1 = -1 * k *
(-x2 * y + x3 * y - x * y1 - x3 * y1 + 2 * x4 * y1 + x * y2 - x4 * y2 - x * y3 + x1 * (y + y3 - 2 * y4) + x * y4 + x2 * y4 +
u_sqrt);
double v2 = -1 * l *
(x1 * y + x3 * y - x4 * y - x * y1 - 2 * x3 * y1 + x * y2 - -2 * x1 * y2 + x4 * y2 - x * y3 + x1 * y3 + x * y4 - x2 * (y - 2 * y1 + y4) +
v_sqrt);
/////////////////////////////////////////////////////////////////////////////////////////////////
double u2 = k *
(x2 * y - x3 * y + x4 * y + x * y1 + x3 * y1 - 2 * x4 * y1 - x * y2 + x4 * y2 + x * y3 - x1 * (y + y3 - 2 * y4) - x * y4 - x2 * y4 +
u_sqrt);
在大多数情况下它是 u1 和 v1 所以应该不需要计算其他的。
我用它来校准 Pegasus Air-Pen 设备(超声波笔)在 sheet 纸上的坐标。如果您的点 1 到点 5 的坐标也 >= 0,它的效果最好。
很抱歉 post 将此作为答案,但评论太长了,我认为这对 post 是一个有价值的帮助我.
我正在尝试找到最简单的方法来确定一个点在四边形中的相对位置。已知的是(见图)点1、2、3、4、5在xy坐标系中的位置:x1,y1,x2,y2,x3,y3,x4,y4,x5,y5.
也已知点1、2、3、4在ξ-η坐标系中的位置(见图)
根据这些数据,我想确定点 5 的 ξ 和 η 是多少。
结果
感谢所有回答的人!我发现@dbc 和@agentp 的解决方案相似。我还发现这个解决方案比@MBo 的透视变换解决方案更好,因为我不必计算矩阵的逆 (Ax=B --> x=inv(A)*B).
我得到以下结果:
u = 0.5 * (ξ + 1)
v = 0.5 * (η + 1)
在我的例子中,所有点都在矩形内,因此 u>0 和 v>0。
您需要计算一个透视变换矩阵,将源四边形的 4 个点映射到目标四边形的 4 个点 (example) (more mathemathics),然后将此变换应用于第 5 个点的坐标(将矩阵乘以坐标向量)
这看起来像是标准的有限元参数化 (这个问题没有具体cify 特定映射,但我想有人可能对这个 specific 案例感兴趣)
{x, y} == (
(1 - eta) (1 - ci) {p1x, p1y} +
(1 - eta) (1 + ci) {p2x, p2y} +
(1 + eta) (1 + ci) {p3x, p3y} +
(1 + eta) (1 - ci) {p4x, p4y} )/4
这可以用 {eta,ci} 的封闭形式求解,但是表达式对于 post.
来说相当笨拙在实践中,计算这些常量:
ax = p1x + p2x + p3x + p4x
bx = p1x - p2x - p3x + p4x
cx = p1x + p2x - p3x - p4x
dx = p1x - p2x + p3x - p4x
ay = p1y + p2y + p3y + p4y
by = p1y - p2y - p3y + p4y
cy = p1y + p2y - p3y - p4y;
dy = p1y - p2y + p3y - p4y;
求解此二次方程 eta
:
(ax by - bx ay) - 4 (by x - bx y) +
eta (dx ay - cx by + bx cy - ax dy + 4 (x dy - dx y)) +
eta^2 (cx dy - dx cy) == 0
然后得到 ci
为:
ci = ((-ax + eta cx + 4 x)/(-bx + eta dx))
如果多边形不太扭曲,只有一种解决方案将满足 -1<eta<1
和 -1<ci<1
你这里有一个 2d bilinear blended surface。为简单起见,让我们将其坐标更改为从零到一:
u = 0.5 * (ξ + 1)
v = 0.5 * (η + 1)
在那种情况下,表面评估器可以表示为
F(u, v) = P1 + u * (P2 - P1) + v * ((P4 + u * (P3 - P4)) - (P1 + u * (P2 - P1)))
即对于给定的u
,构造一条通过以下两点的直线:
Pv0 = P1 + u * (P2 - P1);
Pv1 = P4 + u * (P3 - P4);
然后在给定的 v
F(u, v) = Pv0 + v * (Pv1 - Pv0)
您要寻找的是满足 F(u, v) = P5
的值 (u,v)
。当从 Pv0
到 Pv1
的线穿过 P5
时,给定 u
会发生这种情况,这将在 P5 - Pv0
平行于 Pv1 - Pv0
时发生-- 即当它们的 2d cross 为零时:
cross2d(P5 - Pv0, Pv1 - Pv0) = 0
⇒
cross2d(P5 - (P1 + u * (P2 - P1)),
P4 + u * (P3 - P4) - (P1 + u * (P2 - P1))) = 0
现在,两个二维向量 A ⨯ B
的 2d cross 由 Ax*By - Ay*Bx
给出,因此等式变为
(x5 - (x1 + u * (x2 - x1))) * (y4 + u * (y3 - y4) - (y1 + u * (y2 - y1))) - (y5 - (y1 + u * (y2 - y1))) * (x4 + u * (x3 - x4) - (x1 + u * (x2 - x1))) = 0
Expanding this expression out and collecting collecting together terms in u
,我们得到
u^2 * (x1*y3 - x1*y4 - x2*y3 + x2*y4 + (-x3)*y1 + x3*y2 + x4*y1 - x4*y2)
+ u * (-x1*y3 + 2*x1*y4 - x1*y5 - x2*y4 + x2*y5 + x3*y1 - x3*y5 - 2*x4*y1 + x4*y2 + x4*y5 + x5*y1 - x5*y2 + x5*y3 - x5*y4)
+ (-x1*y4 + x1*y5 + x4*y1 - x4*y5 - x5*y1 + x5*y4)
= 0
现在这是 u
上的二次方程,并且可以是 solved。请注意,在四边形的顶部和底部边缘平行的情况下,二次方程将转化为线性方程;你的二次方程求解器必须处理这个问题。
double a = (x1 * y3 - x1 * y4 - x2 * y3 + x2 * y4 + (-x3) * y1 + x3 * y2 + x4 * y1 - x4 * y2);
double b = (-x1 * y3 + 2 * x1 * y4 - x1 * y5 - x2 * y4 + x2 * y5 + x3 * y1 - x3 * y5 - 2 * x4 * y1 + x4 * y2 + x4 * y5 + x5 * y1 - x5 * y2 + x5 * y3 - x5 * y4);
double c = (-x1 * y4 + x1 * y5 + x4 * y1 - x4 * y5 - x5 * y1 + x5 * y4);
double[] solutions = Quadratic.Solve(a, b, c);
解决方案可能不止一种。对于退化四边形也可能没有解决方案。
求解 u
的值后,找到等效的 v
很简单。给定积分
Pv0 = P1 + u * (P2 - P1);
Pv1 = P4 + u * (P3 - P4);
你寻求 v
这样
v * (Pv1 - Pv0) = P5 - Pv0;
选择坐标索引 0 或 1,使 |(Pv1 - Pv0)[index]|
最大化 。 (如果两个坐标几乎都为零,那么放弃——这个具体u
无解。然后设置
v = (P5 - Pv0)[index] / (Pv1 - Pv0)[index];
最后,如果您有多个解决方案,请选择混合 [u, v]
边界内的解决方案。然后最后设置
ξ = 2 * u - 1;
η = 2 * v - 1;
参考@blaz 的自我回答(请为@blaze、@dbc 和@agentp 的回答投票)
对于不愿意手工复制公式的人,这里是 C# 代码的公式:
double v_sqrt = Math.Sqrt(
4 * (
(x3 - x4) * (y1 - y2) - (x1 - x2) * (y3 - y4)) * (x4 * (-1 * y + y1) + x1 * (y - y4) + x * (-1 * y1 + y4)) +
Math.Pow(
(x3 * y - x4 * y - x3 * y1 + 2 * x4 * y1 - x4 * y2 + x1 * (y + y3 - 2 * y4) + x2 * (-1 * y + y4) + x * (-1 * y1 + y2 - y3 + y4))
, 2)
);
double u_sqrt = Math.Sqrt(
4 * ((x3 - x4) * (y1 - y2) - (x1 - x2) * (y3 - y4))
* (
x4 * (-1 * y + y1) + x1 * (y - y4) + x * (-1 * y1 + y4)
) +
Math.Pow(
(x3 * y - x4 * y - x3 * y1 + 2 * x4 * y1 - x4 * y2 + x1 * (y + y3 - 2 * y4) + x2 * (-1 * y + y4) + x * (-1 * y1 + y2 - y3 + y4))
, 2)
);
double k = 1 / (2 * ((x3 - x4) * (y1 - y2) - (x1 - x2) * (y3 - y4)));
double l = 1 / (2 * ((x1 - x4) * (y2 - y3) - (x2 - x3) * (y1 - y4)));
///////////////////////////////////////////////////////////////////////////////////////////////
double v1 = l *
(x2 * y - x3 * y + x4 * y + x * y1 - 2 * x2 * y1 + x3 * y1 - x * y2 - x4 * y2 + x * y3 - x1 * (y - 2 * y2 + y3) - x * y4 + x2 * y4 +
v_sqrt);
///////////////////////////////////////////////////////////////////////////////////////////////
double u1 = -1 * k *
(-x2 * y + x3 * y - x * y1 - x3 * y1 + 2 * x4 * y1 + x * y2 - x4 * y2 - x * y3 + x1 * (y + y3 - 2 * y4) + x * y4 + x2 * y4 +
u_sqrt);
double v2 = -1 * l *
(x1 * y + x3 * y - x4 * y - x * y1 - 2 * x3 * y1 + x * y2 - -2 * x1 * y2 + x4 * y2 - x * y3 + x1 * y3 + x * y4 - x2 * (y - 2 * y1 + y4) +
v_sqrt);
/////////////////////////////////////////////////////////////////////////////////////////////////
double u2 = k *
(x2 * y - x3 * y + x4 * y + x * y1 + x3 * y1 - 2 * x4 * y1 - x * y2 + x4 * y2 + x * y3 - x1 * (y + y3 - 2 * y4) - x * y4 - x2 * y4 +
u_sqrt);
在大多数情况下它是 u1 和 v1 所以应该不需要计算其他的。
我用它来校准 Pegasus Air-Pen 设备(超声波笔)在 sheet 纸上的坐标。如果您的点 1 到点 5 的坐标也 >= 0,它的效果最好。
很抱歉 post 将此作为答案,但评论太长了,我认为这对 post 是一个有价值的帮助我.