使用C#实现Hoey Shamos算法

4
好的,我现在正在从我的当前算法中获取正确信息!然而,有700,000个多边形需要检查,速度太慢了!先前的问题已经解决(我的Line2D intersectsWith方法不正确)。
现在问题在于确定我的瓶颈在哪里!这个算法应该是O(nlog-n),所以应该更快。我的intersectsWith方法看起来无法更快,但是我会发布它的代码,以防我错了。
编辑:添加了IComparable接口
我的方法用于读取线段相交。一些代码已被省略以提高可读性。
    public class Line2D : IComparable
    {


    public Line2D(XYPoints p1, XYPoints p2)
    {

    }

    public bool intersectsLine(Line2D comparedLine)
    {

        if ((X2 == comparedLine.X1) && (Y2 == comparedLine.Y1)) return false;
        if ((X1 == comparedLine.X2) && (Y1 == comparedLine.Y2)) return false;

        if (X2 == comparedLine.X1 && Y2 == comparedLine.Y1)
        {
            return false;
        }

        if (X1 == comparedLine.X2 && Y1 == comparedLine.Y2)
        {
            return false;
        }
        double firstLineSlopeX, firstLineSlopeY, secondLineSlopeX, secondLineSlopeY;

        firstLineSlopeX = X2 - X1;
        firstLineSlopeY = Y2 - Y1;

        secondLineSlopeX = comparedLine.getX2() - comparedLine.getX1();
        secondLineSlopeY = comparedLine.getY2() - comparedLine.getY1();

        double s, t;
        s = (-firstLineSlopeY * (X1 - comparedLine.getX1()) + firstLineSlopeX * (getY1() - comparedLine.getY1())) / (-secondLineSlopeX * firstLineSlopeY + firstLineSlopeX * secondLineSlopeY);
        t = (secondLineSlopeX * (getY1() - comparedLine.getY1()) - secondLineSlopeY * (getX1() - comparedLine.getX1())) / (-secondLineSlopeX * firstLineSlopeY + firstLineSlopeX * secondLineSlopeY);

        if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
        {
            return true;
        }

        return false; // No collision 
    }

    int IComparable.CompareTo(object obj)
    {

        //return Y1.GetHashCode();
        Line2D o1 = this;
        Line2D o2 = (Line2D)obj;
        if (o1.getY1() < o2.getY1())
        {
            return -1;
        }
        else if (o1.getY1() > o2.getY2())
        {
            return 1;
        }
        else
        {
            if (o1.getY2() < o2.getY2())
            {
                return -1;
            }
            else if (o1.getY2() > o2.getY2())
            {
                return 1;
            }
            else
            {
                return 0;
            }
        } 
    }


}

我意识到在算法实现中,List 不是最快的算法,但我需要索引!

//Create a new list, sort by Y values.

 List<AlgEvent> SortedList = events.OrderBy(o => o.getY()).ToList();                
 List<Line2D> sweepline = new List<Line2D>();

 for (var g = 0; g < SortedList.Count; g++)
 {
 if (SortedList[g].isStart)
 {
    Line2D nl = SortedList[g].line;
    Line2D above;
    /* Start generating above */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        if (index == -1)
        {
            above = null;
        }
        else
        {
            above = sweepline[index + 1];
        }


    }
    catch (ArgumentOutOfRangeException)
    {
        above = null;
    }
    /* End generating above */
    if (above != null)
    {
        if (above.intersectsLine(nl))
        {
            return true;
        }
    }
    Line2D below;
    /* Start generating below */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        below = sweepline[index - 1];

    }
    catch (ArgumentOutOfRangeException)
    {

        below = null;

    }
    /* End generating below */
    if (below != null)
    {
        if (below.intersectsLine(nl))
        {
            return true;
        }
    }
    sweepline.Add(nl);
    sweepline = sweepline.OrderBy(o => o.getY1()).ToList();

}
else
{
    Line2D nl = SortedList[g].line;
    Line2D above;
    Line2D below;
    /* Start generating above */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        Console.Out.WriteLine("index:" + index);
        //add 1 to get above line
        above = sweepline[index + 1];

    }
    catch (ArgumentOutOfRangeException)
    {

        above = null;

    }
    /* End generating above */
    /* Start generating below */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        below = sweepline[index - 1];

    }
    catch (ArgumentOutOfRangeException)
    {

        below = null;

    }
    /* End generating below */
    sweepline = sweepline.OrderBy(o => o.getY1()).ToList();
    sweepline.Remove(nl);
    if (above != null && below != null)
    {
        if (above.intersectsLine(below))
        {
            return true;
        }
    }
}
Console.WriteLine("");
  }



   } // end numofparts for-loop

   return false;

============================================

更新:9月12日:

实现了来自C5的TreeSet,为我的类实现了IComparable,并使其变得更慢了?如果这很重要,我仍在对其进行索引。

http://www.itu.dk/research/c5/

使用TreeSet的代码:

TreeSet<Line2D> sweepline = new TreeSet<Line2D>();
for (var g = 0; g < SortedList.Count; g++)
{
if (SortedList[g].isStart)
{
    Line2D nl = SortedList[g].line;
    Line2D above;
    /* Start generating above */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        above = sweepline[index + 1];

    }
    catch (IndexOutOfRangeException)
    {

        above = null;

    }
    /* End generating above */
    if (above != null)
    {
        if (above.intersectsLine(nl))
        {
            return false;
        }
    }
    Line2D below;
    /* Start generating below */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        below = sweepline[index - 1];

    }
    catch (IndexOutOfRangeException)
    {

        below = null;

    }
    /* End generating below */
    if (below != null)
    {
        if (below.intersectsLine(nl))
        {
            return false;
        }
    }
    sweepline.Add(nl);
    //sweepline = sweepline.OrderBy(o => o.getY1()).ToList();

}
else
{
    Line2D nl = SortedList[g].line;
    Line2D above;
    Line2D below;
    /* Start generating above */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //Console.Out.WriteLine("index:" + index);
        //add 1 to get above line
        above = sweepline[index + 1];

    }
    catch (IndexOutOfRangeException)
    {

        above = null;

    }
    /* End generating above */
    /* Start generating below */
    try
    {
        //grab index in sweepline
        int index = sweepline.IndexOf(nl);
        //add 1 to get above line
        below = sweepline[index - 1];

    }
    catch (IndexOutOfRangeException)
    {

        below = null;

    }
    /* End generating below */
    //sweepline = sweepline.OrderBy(o => o.getY1()).ToList();
    sweepline.Remove(nl);
    if (above != null && below != null)
    {
        if (above.intersectsLine(below))
        {
            return false;
        }
    }
}
//Console.WriteLine("");

}


2
你为什么不能自己在C#中实现TreeSet呢? - Zac Howland
那你的意思是我可以扩展OrderedSet类,并添加higher和lower方法? - Evan Parsons
@EvanParsons 我认为Zac的意思是要从头开始实现一个TreeSet。TreeSet和类似的数据结构使用平衡树数据结构(例如红黑树或AVL树),但这些并不容易实现。跳表是一种更简单的数据结构,可以胜任这项工作。 - john
1
@john:是的,实现树形数据结构比试图强制一种不适合算法要求的数据结构与之配合更为简单(方孔遇圆钉)。虽然跳表可以使用(尽管算法代码可能会更难编写),但有序集合则不行。 - Zac Howland
我想我有一个想法。我将使用try-catch块。我将使用列表并使用“indexOf”返回位置。如果它返回-1或不在列表范围内的参数,则将其设置为返回null。 - Evan Parsons
1
一般情况下,您可以使用基于堆的结构来处理事件队列,并使用平衡二叉搜索树来处理Y结构。如果使用列表,会使事情变慢... 如果我没记错的话,所述的nlogn复杂度是使用这些数据结构实现的。 - Kshitij Banerjee
2个回答

3

首先,关于线段相交:你不需要知道实际的交点,只需要知道它们是否相交即可。参见http://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/,这里有一个可以实现此功能的算法。


关于List实现:
在使用List时,您调用indexOf在扫描线上查找nl。这会从头到尾搜索扫描线。请参见List(T).IndexOf。如果您使用BinarySearch方法,应该可以大大加快搜索速度。 List的文档有一个名为性能考虑的段落。他们建议您使用实现IEquatable<T>IComparable<T>的值类型。因此,您的Line2D可能应该是一个struct并实现这些接口。
如果您遵循这个建议,从扫描线中检索端点的时间复杂度应该是O(log n),这对您的目的已经足够,并且内存应该更有效地使用。

对于列表来说,插入和删除的时间复杂度为O(n),因为底层数组需要在内存中移动。一个SortedSet具有更快的插入和删除速度,但我不太清楚如何在其中以O(log n)的时间复杂度找到一个元素的相邻元素。有人知道吗?(另请参见为什么SortedSet<T>.GetViewBetween不是O(log N)?


无论如何,C5的TreeSet应该可以解决这个问题。

我查阅了用户指南中IndexOf和[i]的性能,它们都被列为O(log n)。所以这不应该是问题所在。仍然可能更快一些,但只是一个固定因素,调用sweepline上找到邻居的特定方法,即SuccessorPredecessor,它们也是O(log n)。

所以

[...]
try 
{
    Line2D above = sweepline.Successor(nl);
    if (above.intersectsLine(nl))
    {
        return false;
    }
}
catch (NoSuchItemException ignore) { }
[...]

我不喜欢他们没有不抛出异常的方法,因为抛出异常非常昂贵。通常情况下,您的扫描线将会相当满,所以我最好的猜测是找不到一个将是罕见的情况,并调用Successor是最有效的方法。或者,您可以像现在一样不断调用IndexOf,但在检索[index + 1]之前先检查它是否等于Count减一,并完全防止抛出异常:
[...]
int index = sweepline.IndexOf(nl);
if( index < sweepline.Count-1 )
{
    Line2D above = sweepline[index + 1];
    if (above.intersectsLine(nl))
    {
        return false;
    }
}
[...]

手册的第二章介绍了C5集合的相等性和比较。在这里,至少您必须实现IEquatable<T>IComparable<T>

还有一个想法:您报告称算法处理了700000行。能否从例如1000、2500、5000、10000行开始计时,并观察算法在它们不相交的情况下如何扩展?


如何比较扫描线上的线段:

你需要为Sweepline TreeSet中的Line2Ds找到某种自然排序,因为CompareTo方法要求你将一个Line2D与另一个进行比较。其中一个Line2D已经存在于Sweepline TreeSet中,另一个刚刚被遇到并正在添加。

我认为你的扫描线是从下往上运行的:

List<AlgEvent> SortedList = events.OrderBy(o => o.getY()).ToList();

Sweepline running through the line segments

假设在事件1中向TreeSet添加了线段S1,我们希望将其与正在事件2中添加的S2进行比较。
这些线段可能会在某一点相交,这会改变它们的顺序,但是算法会在插入它们后立即检查这一点,在上面和下面的检查中。考虑到这一点,也许最好称之为左侧和右侧的检查。
总之...所以最简单的方法是比较两条线段的底部端点。左边是较小的,右边是较大的。然而,我们需要查看扫描线位置处的排序情况,因为它们可能已经改变了位置,就像图片中一样。
因此,我们需要将S2的底部端点与S1上的红点进行比较,并查看它是否位于该点的左侧或右侧。它位于左侧,因此认为S2比S1小。
通常情况下,这比较简单:如果S1的所有部分都位于S2的底部端点的左侧,则S1比S2小。如果S1的所有部分都位于S2的底部端点的右侧,则S1比S2大。
我认为你正在寻找接口的类型安全版本:
public class Line2D : IComparable<Line2D>

假设有两个属性BottomY,它们是两个Y值中的最小值,以及BottomX,它是最低端点的X值,以下是一个略有测试的尝试:
int IComparable<Line2D>.CompareTo(Line2D other)
{
    if( BottomY < other.BottomY )
    {
        return -other.CompareTo(this);
    }

    // we're the segment being added to the sweepline
    if( BottomX >= other.X1 && BottomX >= other.X2 )
    {
        return 1;
    }
    if( BottomX <= other.X1 && BottomX <= other.X2 )
    {
        return -1;
    }

    if( other.Y2 == other.Y1 )
    {
        // Scary edge case: horizontal line that we intersect with. Return 0?
        return 0;
    }

    // calculate the X coordinate of the intersection of the other segment
    // with the sweepline
    // double redX = other.X1 + 
    //    (BottomY - other.Y1) * (other.X2 - other.X1) / (other.Y2 - other.Y1);
    //
    // return BottomX.CompareTo(redX);

    // But more efficient, and more along the lines of the orientation comparison:
    return Comparer<Double>.Default.Compare(
        (BottomX - other.X1) * (other.Y2 - other.Y1),
        (BottomY - other.Y1) * (other.X2 - other.X1) );

}

我可以把数据限制在每次只读取 100,000 条。这比每次读取 700,000 条并记录时间更容易和更快。 - Evan Parsons
我正在实现IComparable时遇到了困难。目前它返回了600,000个“无效”多边形,这绝对是错误的。我将用修订后的Line2D类更新原始帖子。 - Evan Parsons
我已经添加了一些关于如何比较Line2Ds的想法。 - flup
我仍然感到很困惑,需要一些时间来理解并报告回来,再次感谢! - Evan Parsons
我以为我的线是从左到右扫过去的。因为在从左到右扫描时,它们会遇到Y,所以我按Y来组织它们。这是我第一次看到底部的Y,在这个算法中使用了另一个底部的Y。我还没有实现IEquatable,但我已经重写了getHashCode和Equals方法。 - Evan Parsons
显示剩余4条评论

1

[原回答]

我不是C#用户,但这应该会加快一些速度。

  • 减少堆垃圾
  • 不要重复计算相同的东西
  • 尽可能避免所有子调用(删除获取函数)

代码:

public bool intersectsLine(const Line2D &comparedLine)
    {
    if ((X2==comparedLine.X1)&&(Y2==comparedLine.Y1)) return false;
    if ((X1==comparedLine.X2)&&(Y1==comparedLine.Y2)) return false;
    double dx1,dy1,dx2,dy2;
    dx1 = X2 - X1;
    dy1 = Y2 - Y1;
    dx2 = comparedLine.X2 - comparedLine.X1;
    dy2 = comparedLine.Y2 - comparedLine.Y1;
    double s,t,ax,ay,b;
    ax=X1-comparedLine.X1;
    ay=Y1-comparedLine.Y1;
    b=1.0/(-(dx2*dy1)+(dx1*dy2));
    s = (-(dy1*ax)+(dx1*ay))*b;
    t = ( (dx2*ay)-(dy2*ax))*b;
    if ((s>=0)&&(s<=1)&&(t>=0)&&(t<=1)) return true;
    return false; // No collision
    }

对于你的代码的其余部分,添加时间测量以找出究竟是什么使得速度变慢。我猜测是列表管理...不必要的重新分配可能会严重拖慢速度。

[编辑1]

在对随机行数据进行一些研究后,我得出了以下结论:

  1. 如果整个区域有太多的线,那么就不存在任何优化。
  2. 如果小线条数量更多,则会有更多的速度提升。
  3. 暴力破解是 T((N*N-N)/2),仍然是 O(N*N)。估计处理700K行需要约35小时(无法使用)。
  4. 带有区域划分优化的暴力破解是 T((((N/M)^2)-N)/2) - 优化 ~O((N/M)^2),其中:

    • N 是区域线数的最大值;
    • M 是每个轴上区域分割的数目。

      思路是只检查跨越某个区域的线(将数据集区域划分为 M*M 的正方形/矩形)。对于700K行,最佳划分是 16x16 区域。测量时间:

      32K 行每0.540秒 64K 行每1.950秒 128K 行每7.000秒 256K 行每27.514秒

    估计运行时间为每 700K 行 3.7 分钟(对于每行长度不超过整个区域的 10%)。我认为这比你的 19 分钟要好。

  5. 还可以使用多个 CPU/核心进行加速

    算法完全可并行化,对于4个 CPU/核心 3.7min/4 -> 56s 或者您可以将其移植到 GPU 上...

我的优化暴力算法与面积细分 O((((N/M)^2)-N)/2) - 优化

  1. 获取已使用区域的大小(xmin,xmax,ymin,ymax)O(N)
  2. 选择细分M。我尝试了最佳随机数据集 32K-256K 的行数为 M=16
  3. 循环遍历所有的细分区域(均匀细分数据集区域)

    创建跨越实际细分区域的线列表,并对该列表中的所有线段进行交点检查。如果不想有重复交点,则丢弃当前区域之外的所有交点

我的代码(我正在使用 BDS2006 C++ 和自己的列表,因此您需要对其进行移植以使其与您的代码兼容)

void Twin_GLView2D::main_intersect(int M=16)
{
int ia,ib,i,j,N;
double zero=1e-6;
glview2D::_lin *l;
glview2D::_pnt p;

struct _line
    {
    double bx0,by0,bx1,by1;     // bounding rectangle
    double x0,y0,dx,dy;         // precomputed params
    } *lin,*a,*b;

struct _siz
    {
    double bx0,bx1,by0,by1;     // zone bounding rectangle
    } sz,bz;
List<_line*> zone;

// load and precompute lines
N=view.lin.num;
lin=new _line[N];
if (lin==NULL) return;
for (a=lin,l=view.lin.dat,ia=0;ia<N;ia++,a++,l++)
    {
    // line ...
    if (l->p0.p[0]<=l->p1.p[0]) { a->bx0=l->p0.p[0]; a->bx1=l->p1.p[0]; }
    else                        { a->bx0=l->p1.p[0]; a->bx1=l->p0.p[0]; }
    if (l->p0.p[1]<=l->p1.p[1]) { a->by0=l->p0.p[1]; a->by1=l->p1.p[1]; }
    else                        { a->by0=l->p1.p[1]; a->by1=l->p0.p[1]; }
    a->x0=l->p0.p[0]; a->dx=l->p1.p[0]-l->p0.p[0];
    a->y0=l->p0.p[1]; a->dy=l->p1.p[1]-l->p0.p[1];
    // global image size for zone subdivision
    if (!ia)
        {
        sz.bx0=l->p0.p[0];
        sz.by0=l->p0.p[1];
        sz.bx1=sz.bx0;
        sz.by1=sz.by0;
        }
    if (sz.bx0>l->p0.p[0]) sz.bx0=l->p0.p[0];
    if (sz.bx1<l->p0.p[0]) sz.bx1=l->p0.p[0];
    if (sz.by0>l->p0.p[1]) sz.by0=l->p0.p[1];
    if (sz.by1<l->p0.p[1]) sz.by1=l->p0.p[1];
    if (sz.bx0>l->p1.p[0]) sz.bx0=l->p1.p[0];
    if (sz.bx1<l->p1.p[0]) sz.bx1=l->p1.p[0];
    if (sz.by0>l->p1.p[1]) sz.by0=l->p1.p[1];
    if (sz.by1<l->p1.p[1]) sz.by1=l->p1.p[1];
    }
// process lines by zonal subdivision
zone.allocate(N);
view.pnt.num=0; view.pnt.allocate(view.lin.num);
sz.bx1-=sz.bx0; sz.bx1/=double(M);
sz.by1-=sz.by0; sz.by1/=double(M);
for (bz.by0=sz.by0,bz.by1=sz.by0+sz.by1,i=0;i<M;i++,bz.by0+=sz.by1,bz.by1+=sz.by1)
for (bz.bx0=sz.bx0,bz.bx1=sz.bx0+sz.bx1,j=0;j<M;j++,bz.bx0+=sz.bx1,bz.bx1+=sz.bx1)
    {
    // create list of lines for actual zone only
    zone.num=0;         // clear zone list
    for (a=lin,ia=   0;ia<N;ia++,a++)
     if ((a->bx0<=bz.bx1)&&(a->bx1>=bz.bx0))
      if ((a->by0<=bz.by1)&&(a->by1>=bz.by0))
       zone.add(a); // add line to zone list
    // check for intersection within zone only
    // O((((N/M)^2)-N)/2) - optimizations
    for (ia=   0,a=zone.dat[ia];ia<zone.num;ia++,a=zone.dat[ia])
    for (ib=ia+1,b=zone.dat[ib];ib<zone.num;ib++,b=zone.dat[ib])
        {
        // discart lines with non intersecting bound rectangles
        if (a->bx1<b->bx0) continue;
        if (a->bx0>b->bx1) continue;
        if (a->by1<b->by0) continue;
        if (a->by0>b->by1) continue;
        // 2D lines a,b intersect ?
        double x0,y0,x1,y1,t0,t1;
        // compute intersection
        t1=divide(a->dx*(a->y0-b->y0)+a->dy*(b->x0-a->x0),(a->dx*b->dy)-(b->dx*a->dy));
        x1=b->x0+(b->dx*t1);
        y1=b->y0+(b->dy*t1);
        if (fabs(a->dx)>=fabs(a->dy)) t0=divide(b->x0-a->x0+(b->dx*t1),a->dx);
        else                          t0=divide(b->y0-a->y0+(b->dy*t1),a->dy);
        x0=a->x0+(a->dx*t0);
        y0=a->y0+(a->dy*t0);
        // check if intersection exists
        if (fabs(x1-x0)>zero) continue;
        if (fabs(y1-y0)>zero) continue;
        if ((t0<0.0)||(t0>1.0)) continue;
        if ((t1<0.0)||(t1>1.0)) continue;
        // if yes add point
        p.p[0]=x0;
        p.p[1]=y0;
        p.p[2]=0.0;
        // do not add points out of zone (allmost all duplicit points removal)
        if (x0<bz.bx0) continue;
        if (x0>bz.bx1) continue;
        if (y0<bz.by0) continue;
        if (y0>bz.by1) continue;
        view.pnt.add(p);
        }
    }
view.redraw=true;
delete lin;
}

注意:

  1. List<T> x;T x[] 是相同的,只是大小为“无限制”
    • x.num;x[] 中 T 类型元素的实际大小,而不是字节数!index = <0,x.num-1>
    • x.add(q);q 添加到列表末尾
    • x.num=0; 清空列表
    • x.allocate(N); 为列表中的 N 个项目分配空间,以避免重新定位
  2. 输入的 List<>view.lin,包含点 p0,p1,每个点都有 double p[2] ... x,y
  3. 输出的 List<>view.pnt,包含 double p[2] ... x,y

[编辑2]

此外,我发现上述算法的最佳性能是当 M=12+(N>>15)

  • M 是每个轴的细分区域数
  • N 是要检查的线条数量

你可能已经找到了加速的方法。你提供的代码并没有比我的更快,但是当我强制我的“intersectsLine”方法返回true或false时,它只需要19分钟就能完成,而不是超时。有没有更快的方法来检测是否发生了交叉?(不需要精确的点) - Evan Parsons
很难在不了解数据集属性的情况下说出来,通常聪明的重新排序、排序和/或删除不重要的数据会导致速度大幅提高。 - Spektre
今晚我将尝试您的暴力方法。我应该提到,它有700,000个多边形,每个多边形之间有5到4,500条线。 - Evan Parsons
我从一个二进制文件中读取了七十万的数据,没有内存问题可言。 - Evan Parsons
很抱歉,我不能分享数据,但是我已经更新了我的代码以回答这个问题:https://dev59.com/2HXZa4cB1Zd3GeqPAv5o。 - Evan Parsons
显示剩余4条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接