两个三维矩形的交集

Question

两个三维矩形的交集

c#algorithmvisual-studio-2008geometry3d

11

为了在三维空间中获取两个矩形的交线，我将它们转换成平面，并使用它们法向量的叉积来获取交线，然后尝试获取每个矩形线段与交线的交点。

问题是该线段与三条线段平行，并且只与一条线段在NAN、NAN、NAN处相交，这完全是错误的。您能告诉我代码哪里有问题吗？

我使用了此链接中的vector3，并创建了以下平面类：

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace referenceLineAlgorithm
{
struct Line
{

    public Vector3 direction;
    public Vector3 point;

}

struct lineSegment
{

    public Vector3 firstPoint;
    public Vector3 secondPoint;

}

class plane_test
{
    public enum Line3DResult
    {
        Line3DResult_Parallel = 0,
        Line3DResult_SkewNoCross = 1,
        Line3DResult_SkewCross = 2
    };

    #region Fields

    public Vector3 Normal;
    public float D;
    public Vector3[] cornersArray;
    public Vector3 FirstPoint;
    public Vector3 SecondPoint;
    public Vector3 temp;
    public Vector3 normalBeforeNormalization;


    #endregion

    #region constructors

    public plane_test(Vector3 point0, Vector3 point1, Vector3 point2, Vector3 point3)
    {
        Vector3 edge1 = point1 - point0;
        Vector3 edge2 = point2 - point0;
        Normal = edge1.Cross(edge2);
        normalBeforeNormalization = Normal;

        Normal.Normalize();
        D = -Normal.Dot(point0);

        ///// Set the Rectangle corners 
        cornersArray = new Vector3[] { point0, point1, point2, point3 };

    }

    #endregion

    #region Methods
    /// <summary>
    /// This is a pseudodistance. The sign of the return value is
    /// positive if the point is on the positive side of the plane,
    /// negative if the point is on the negative side, and zero if the
    ///  point is on the plane.
    /// The absolute value of the return value is the true distance only
    /// when the plane normal is a unit length vector.
    /// </summary>
    /// <param name="point"></param>
    /// <returns></returns>
    public float GetDistance(Vector3 point)
    {
        return Normal.Dot(point) + D;
    }

    public void Intersection(plane_test SecondOne)
    {
        ///////////////////////////// Get the parallel to the line of interrsection (Direction )
        Vector3 LineDirection = Normal.Cross(SecondOne.Normal);

        float d1 = this.GetDistance(LineDirection);
        float d2 = SecondOne.GetDistance(LineDirection);

        temp = (LineDirection - (this.Normal * d1) - (SecondOne.Normal * d2));

        temp.x = Math.Abs((float)Math.Round((decimal)FirstPoint.x, 2));
        temp.y = Math.Abs((float)Math.Round((decimal)FirstPoint.y, 2));

        Line line;
        line.direction = LineDirection;
        line.point = temp;

        ////////// Line segments 

        lineSegment AB, BC, CD, DA;

        AB.firstPoint = cornersArray[0]; AB.secondPoint = cornersArray[1];
        BC.firstPoint = cornersArray[1]; BC.secondPoint = cornersArray[2];
        CD.firstPoint = cornersArray[2]; CD.secondPoint = cornersArray[3];
        DA.firstPoint = cornersArray[3]; DA.secondPoint = cornersArray[0];

        Vector3 r1 = new Vector3(-1, -1, -1);
        Vector3 r2 = new Vector3(-1, -1, -1);
        Vector3 r3 = new Vector3(-1, -1, -1);
        Vector3 r4 = new Vector3(-1, -1, -1);

        /*
        0,0 |----------------| w,0
            |                |
            |                |
        0,h |________________|  w,h


         */

        IntersectionPointBetweenLines(AB, line, ref r1);
        IntersectionPointBetweenLines(BC, line, ref r2);
        IntersectionPointBetweenLines(CD, line, ref r3);
        IntersectionPointBetweenLines(DA, line, ref r4);

        List<Vector3> points = new List<Vector3>();
        points.Add(r1);
        points.Add(r2);
        points.Add(r3);
        points.Add(r4);
        points.RemoveAll(

           t => ((t.x == -1) && (t.y == -1) && (t.z == -1))


           );

        if (points.Count == 2)
        {
            FirstPoint = points[0];
            SecondPoint = points[1];


        }




    }

    public Line3DResult IntersectionPointBetweenLines(lineSegment first, Line aSecondLine, ref Vector3 result)
    {
        Vector3 p1 = first.firstPoint;
        Vector3 n1 = first.secondPoint - first.firstPoint;


        Vector3 p2 = aSecondLine.point;
        Vector3 n2 = aSecondLine.direction;

        bool parallel = AreLinesParallel(first, aSecondLine);
        if (parallel)
        {

            return Line3DResult.Line3DResult_Parallel;
        }
        else
        {
            float d = 0, dt = 0, dk = 0;
            float t = 0, k = 0;

            if (Math.Abs(n1.x * n2.y - n2.x * n1.y) > float.Epsilon)
            {
                d = n1.x * (-n2.y) - (-n2.x) * n1.y;
                dt = (p2.x - p1.x) * (-n2.y) - (p2.y - p1.y) * (-n2.x);
                dk = n1.x * (p2.x - p1.x) - n1.y * (p2.y - p1.y);
            }
            else if (Math.Abs(n1.z * n2.y - n2.z * n1.y) > float.Epsilon)
            {
                d = n1.z * (-n2.y) - (-n2.z) * n1.y;
                dt = (p2.z - p1.z) * (-n2.y) - (p2.y - p1.y) * (-n2.z);
                dk = n1.z * (p2.z - p1.z) - n1.y * (p2.y - p1.y);
            }
            else if (Math.Abs(n1.x * n2.z - n2.x * n1.z) > float.Epsilon)
            {
                d = n1.x * (-n2.z) - (-n2.x) * n1.z;
                dt = (p2.x - p1.x) * (-n2.z) - (p2.z - p1.z) * (-n2.x);
                dk = n1.x * (p2.x - p1.x) - n1.z * (p2.z - p1.z);
            }

            t = dt / d;
            k = dk / d;

            result = n1 * t + p1;

            // Check if the point on the segmaent or not 
           // if (! isPointOnSegment(first, result))
            //{
               // result = new Vector3(-1,-1,-1);


           // }

            return Line3DResult.Line3DResult_SkewCross;

        }



    }
    private bool AreLinesParallel(lineSegment first, Line aSecondLine)
    {
        Vector3 vector = (first.secondPoint - first.firstPoint);
        vector.Normalize();

        float kl = 0, km = 0, kn = 0;
        if (vector.x != aSecondLine.direction.x)
        {
            if (vector.x != 0 && aSecondLine.direction.x != 0)
            {
                kl = vector.x / aSecondLine.direction.x;
            }
        }
        if (vector.y != aSecondLine.direction.y)
        {
            if (vector.y != 0 && aSecondLine.direction.y != 0)
            {
                km = vector.y / aSecondLine.direction.y;
            }
        }
        if (vector.z != aSecondLine.direction.z)
        {
            if (vector.z != 0 && aSecondLine.direction.z != 0)
            {
                kn = vector.z / aSecondLine.direction.z;
            }
        }

        // both if all are null or all are equal, the lines are parallel
        return (kl == km && km == kn);




    }

    private bool isPointOnSegment(lineSegment segment, Vector3 point)
    {
        //(x - x1) / (x2 - x1) = (y - y1) / (y2 - y1) = (z - z1) / (z2 - z1)
        float component1 = (point.x - segment.firstPoint.x) / (segment.secondPoint.x  - segment.firstPoint.x);
        float component2 = (point.y - segment.firstPoint.y) / (segment.secondPoint.y - segment.firstPoint.y);
        float component3 = (point.z - segment.firstPoint.z) / (segment.secondPoint.z - segment.firstPoint.z); 

        if ((component1 == component2) && (component2 == component3))
        {
            return true;


        }
        else
        {
            return false;

        }

    }

    #endregion
}
}

static void Main(string[] args)
    {

        //// create the first plane points 
        Vector3 point11 =new Vector3(-255.5f, -160.0f,-1.5f) ;    //0,0
        Vector3 point21 = new Vector3(256.5f, -160.0f, -1.5f);   //0,w
        Vector3 point31 = new Vector3(256.5f, -160.0f, -513.5f); //h,0
        Vector3 point41 = new Vector3(-255.5f, -160.0f, -513.5f); //w,h 

        plane_test plane1 = new plane_test(point11, point21, point41, point31);

        //// create the Second plane points 

        Vector3 point12 = new Vector3(-201.6289f, -349.6289f, -21.5f);
        Vector3 point22 =new Vector3(310.3711f,-349.6289f,-21.5f);
        Vector3 point32 = new Vector3(310.3711f, 162.3711f, -21.5f);
        Vector3 point42 =new Vector3(-201.6289f,162.3711f,-21.5f);
        plane_test plane2 = new plane_test(point12, point22, point42, point32);


        plane2.Intersection(plane1);



    }

这是测试值。最好的问候。

- AMH

2

但我相信交集代码可以简化为其本质，并不使用所有类的完整代码。否则，这些类将非常小，无论如何都可以发布。 - Christian Rau

我不知道该怎么做才能让你帮助我。 - AMH

你能开发这个吗："问题是这条线与三个线段平行，并且只与一个位置（NAN，NAN，NAN）相交，这是完全错误的。如果可以，请告诉我代码中出了什么问题"。 - Ricky Bobby

2个回答

3

该程序计算通过两个矩形的平面的交点线。然后，程序寻找该线与其中一个矩形的边缘之间的交点。如果找到这样的两个点的交点，则返回它们。我不会争论这是否是明智的做法，因为我不知道程序的上下文。

让我们浏览一下代码，看看可能出错的地方。

该程序计算通过两个平面的直线如下：

Vector3 LineDirection = Normal.Cross(SecondOne.Normal);

float d1 = this.GetDistance(LineDirection);
float d2 = SecondOne.GetDistance(LineDirection);

temp = (LineDirection - (this.Normal * d1) - (SecondOne.Normal * d2));

temp.x = Math.Abs((float)Math.Round((decimal)FirstPoint.x, 2));
temp.y = Math.Abs((float)Math.Round((decimal)FirstPoint.y, 2));

Line line;
line.direction = LineDirection;
line.point = temp;

线方向的计算是正确的，但是点的计算是错误的，你可能已经知道了。但我会假装我们有一个有效的点和方向，并继续执行程序的其余部分。

程序调用 AreLinesParallel() 来摆脱与平面线平行的边缘。代码如下：

Vector3 vector = (first.secondPoint - first.firstPoint);
vector.Normalize();

float kl = 0, km = 0, kn = 0;
if (vector.x != aSecondLine.direction.x)
{
    if (vector.x != 0 && aSecondLine.direction.x != 0)
    {
        kl = vector.x / aSecondLine.direction.x;
    }
}
if (vector.y != aSecondLine.direction.y)
{
    if (vector.y != 0 && aSecondLine.direction.y != 0)
    {
        km = vector.y / aSecondLine.direction.y;
    }
}
if (vector.z != aSecondLine.direction.z)
{
    if (vector.z != 0 && aSecondLine.direction.z != 0)
    {
        kn = vector.z / aSecondLine.direction.z;
    }
}

// both if all are null or all are equal, the lines are parallel
return ((kl == km && km == kn));

这段代码主要检查边缘方向的元素除以线路方向的元素是否相等。但是，这是一种不可靠的过程。由于舍入误差，即使AreLinesParallel()声称这些线路实际上并不平行，后续程序仍然可能出现除以零的情况。最好不要使用此过程。

现在来到代码的核心部分，即边缘和线路之间的相交测试：

float d = 0, dt = 0, dk = 0;
float t = 0, k = 0;

if (Math.Abs(n1.x * n2.y - n2.x * n1.y) > float.Epsilon)
{
    d = n1.x * (-n2.y) - (-n2.x) * n1.y;
    dt = (p2.x - p1.x) * (-n2.y) - (p2.y - p1.y) * (-n2.x);
    dk = n1.x * (p2.x - p1.x) - n1.y * (p2.y - p1.y);
}
else if (Math.Abs(n1.z * n2.y - n2.z * n1.y) > float.Epsilon)
{
    d = n1.z * (-n2.y) - (-n2.z) * n1.y;
    dt = (p2.z - p1.z) * (-n2.y) - (p2.y - p1.y) * (-n2.z);
    dk = n1.z * (p2.z - p1.z) - n1.y * (p2.y - p1.y);
}
else if (Math.Abs(n1.x * n2.z - n2.x * n1.z) > float.Epsilon)
{
    d = n1.x * (-n2.z) - (-n2.x) * n1.z;
    dt = (p2.x - p1.x) * (-n2.z) - (p2.z - p1.z) * (-n2.x);
    dk = n1.x * (p2.x - p1.x) - n1.z * (p2.z - p1.z);
}

t = dt / d;
k = dk / d;

result = n1 * t + p1;

这段代码的一个问题是缺少注释来解释算法来源。如果没有文档化的算法可供参考，注释可以包含导致公式的推导。第一个分支处理 (x, y)，第二个处理 (y, z)，第三个处理 (z, x)，所以我认为这些分支在2D中解决了交点，并将这些结果提升到3D中。它通过计算行列式来检查每个2D投影的平行线。我不应该做这种反向工程。

无论如何，这就是产生 NaN 值的代码。没有触发任何三个分支，所以最后d = 0，这会导致除以零。与其依赖于AreLinesParallel()避免除以零，还不如检查实际关键值，即d。

当然，这段代码仍需要更多的工作，因为我们还不知道这些线是否也在3D中相交。而且只有在0 <= t && t <= 1时，该点才位于边缘。随着早期错误的修复，可能会出现更多的漏洞。

- antonakos

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

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

- Ricky Bobby · Accepted Answer

首先，您需要明确一件事：

您所说的3D矩形是指3D平面上的平面矩形（而不是长方体）。

假设您的矩形既不共面也不平行，并且因此存在一个唯一的线D1，它表示每个矩形所描述的平面的交点。

在此假设下，矩形R1和R2相交有4种可能的情况：

enter image description here

（注：有时候D1既不与R1相交也不与R2相交，因此可以稍微旋转R1，R2使D1不总是在平行边上相交，而是在连续边上相交。）

当两个矩形相交时，D1总是在相同的交点上与R1和R2相交（请参见第1和第2张图片）。

您的模型不好，因为您的线不能与同一个矩形的3条线段平行...

如您在这个问题中所询问的：3D lines intersection algorithm，一旦您有了D1（请参见Get endpoints of the line segment defined by the intersection of two rectangles），只需确定它与矩形的每个线段的交点。（需要检查每个矩形的4条线段）

然后检查是否有共同的交点... 如果找到一个交点，那么您的矩形相交。

很抱歉很难直接检查代码，但是我想通过这些信息您应该能够找到错误。

希望对您有所帮助。

编辑：

定义一个矩形需要一个点和2个向量：

R2 {A ,u ,v}
R1 {B, u',v'}

定义由R1和R2描述的平面：P1和P2

P1（和P2）的一个正交向量是n1（和n2）。设n1 = u ^ v和n2 = u' ^ v '，其中：

enter image description here

则

P1: n1.(x-xA,y-yA,z-zA)=0
P2: n2.(x-xB,y-yB,z-zB)=0

那么如果您只是寻找D1，则D1的方程为：

D1: P1^2 + P2 ^2 =0 (x,y,z verify P1 =0  an P2 =0 )

D1 : n1.(x-xA,y-yA,z-zA)^2 + n2.(x-xB,y-yB,z-zB)^2 =0

因此，仅通过矩形的表达式，您就可以获得D1的方程闭式公式。

现在让我们看看交点：

R1中的4个点为：

{ A , A+u , A+v, A+u+v }

按照3D线段交点算法所述进行操作：

D1 inter [A,A+u] = I1
D1 inter [A,A+v] = I2
D1 inter [A+u,A+u+v] = I3
D1 inter [A+v,A+u+v] = I4

(I1、I2、I3、I4 可以为空)

same for D2 you get I1' I2' I3' I4'

如果Ij'=Ik'不为null，则它是交点。

如果你按步骤正确执行，你应该能得到正确的解决方案；除非我没有完全理解问题...