用C#实现Hoey Shamos算法

好的,我现在从我当前的算法中获取正确的信息! 但是,要检查700,000个多边形,这太慢了! 上一期是固定的(My 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 = 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.getY2()) { return 1; } else { if (o1.getY2()  o2.getY2()) { return 1; } else { return 0; } } } } 

我的算法实现的大部分,我意识到List不是算法的最快,但我需要索引!:

 //Create a new list, sort by Y values. List SortedList = events.OrderBy(o => o.getY()).ToList(); List sweepline = new List(); for (var g = 0; g  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 sweepline = new TreeSet(); for (var g = 0; g  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(""); 

}

    首先,关于线交叉:您不需要实际的交点,只知道它们是否相交。 请参阅http://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/获取可实现此目的的算法。


    关于List实现:

    在使用List s的实现中,您在sweepline上调用indexOf来查找nl 。 这将从开始到结束搜索扫描线。 见List(T).IndexOf 。 如果您使用BinarySearch方法,那么应该大大加快搜索速度。

    List的文档有一段称为性能注意事项。 他们敦促您使用实现IEquatableIComparable的值类型。 所以,你的Line2D应该是一个struct并实现这些接口。

    如果您遵循该建议,从扫描线检索端点应该是O(log n),这足以满足您的需要,并且应该更有效地使用内存。

    插入和删除对于列表是O(n),因为底层数组需要在内存中移动。 SortedSet有更快的插入和删除,但我不太清楚如何在那里找到O(log n)中项目的邻居。 任何人? (另请参阅为什么SortedSet .GetViewBetween不是O(log N)? )


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

    我在用户指南中查找了IndexOf和[i]的性能,它们都被列为O(log n)。 所以这不应该是问题。 它仍然可能稍微快一些,但不超过一个固定因素,调用特定的方法来寻找扫描线上的邻居,即Successor和Predecessor ,它们也是O(log n)。

    所以

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

    我不喜欢他们没有一个不抛出exception的方法,因为抛出exception是非常昂贵的。 你的扫描线通常会非常充实,所以我最好的猜测是找不到它会很少见,并且调用Successor是最有效的方法。 或者,您可以像现在一样继续调用IndexOf ,但在检索[index + 1]之前检查它是否等于Count[index + 1] ,并防止抛出exception:

     [...] int index = sweepline.IndexOf(nl); if( index < sweepline.Count-1 ) { Line2D above = sweepline[index + 1]; if (above.intersectsLine(nl)) { return false; } } [...] 

    本手册的第二章介绍了C5集合的相等性和比较。 在这里,至少你必须实现IEquatableIComparable

    还有一个想法:你报告喂算算法700000行。 你可以从1000,2500,5000,10000行的时间开始,看看算法如何在它们不相交的情况下进行缩放?


    关于如何比较扫描线上的线条:

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

    我认为你的清扫线从下到上运行:

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

    穿过线段的扫描线

    所以我们假设在事件1中将段S1添加到TreeSet中,我们希望将它与S2进行比较,S2现在正在事件2中添加。

    线段可能在某个点交叉,这会改变排序,但算法会在插入之后检查这个,在上面和下面的检查中。 或许更好的称为左右检查,来考虑它。

    无论如何..所以最简单的方法是比较两个线段的底端点。 左边较小,右边较大。 但是,我们需要查看扫描线位置的顺序,从那时起它们可能会改变位置,如图所示。

    因此,我们需要将S2的底部端点与S1上的红点进行比较,看它是否位于该点的左侧或右侧。 它位于左侧,因此S2被认为小于S1。

    通常它比这更简单:如果S1全部位于S2底部端点的左侧,则S1小于S2。 如果S1全部位于S2的底部端点的右侧,则S1大于S2。

    我想你正在寻找类型更好的界面版本:

     public class Line2D : IComparable 

    假设两个属性BottomY是两个Y值中最低的,而BottomX是最低端点的X值,这是一个有点测试的尝试:

     int IComparable.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.Default.Compare( (BottomX - other.X1) * (other.Y2 - other.Y1), (BottomY - other.Y1) * (other.X2 - other.X1) ); } 

    [原始答案]

    我不是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 } 

    对于其余的代码,添加时间测量以找出确切缓慢的事情。 我的猜测是列表管理......不必要的重新分配会大大减慢速度。

    [EDIT1]

    经过对随机线数据的一些研究后,我得出结论:

    1. 如果整个区域的线数太多,则优化没有效果
    2. 如果有更多的小行,那么任何优化都会有更多的加速
    3. 蛮力是T((N*NN)/2) ,对于700K线待处理(不可用),仍然是O(N*N)估计约35小时
    4. 区域细分的优化powershell是T((((N/M)^2)-N)/2) - 优化~O((N/M)^2)其中

      • N是区域行数的最大值
      • M是每个轴的区域划分数,其思想是仅检查穿过某个区域的线(将数据集区域划分为M*M正方形/矩形)。 对于700K线路,最好细分为16x16区域。 测量时间:

        每32K线0.540s每64K线1.950s每128K线7.000s每256K线27.514s

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

    5. 使用多CPU /核心可以实现另一种加速

      算法是完全并行的,4个CPU /核心3.7min/4 -> 56s或者你可以将它移植到GPU ...

    我的优化powershell算法,区域细分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;iap0.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.bx1p0.p[0]) sz.bx1=l->p0.p[0]; if (sz.by0>l->p0.p[1]) sz.by0=l->p0.p[1]; if (sz.by1p0.p[1]) sz.by1=l->p0.p[1]; if (sz.bx0>l->p1.p[0]) sz.bx0=l->p1.p[0]; if (sz.bx1p1.p[0]) sz.bx1=l->p1.p[0]; if (sz.by0>l->p1.p[1]) sz.by0=l->p1.p[1]; if (sz.by1p1.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;ibx0<=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];iabx1bx0) continue; if (a->bx0>b->bx1) continue; if (a->by1by0) 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 pp[0]=x0; pp[1]=y0; pp[2]=0.0; // do not add points out of zone (allmost all duplicit points removal) if (x0bz.bx1) continue; if (y0bz.by1) continue; view.pnt.add(p); } } view.redraw=true; delete lin; } 

    笔记:

    1. List x; 与“无限”大小的T x[]相同
      • x.num; 是Ts x[]字节x[]实际大小x[]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

    [EDIT2]

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

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