Implementation of the Hoi Chamos algorithm using C #

Ok, now I get the correct information from my current algorithm! However, to check 700,000 polygons is just too slow! Previous issue fixed (My Line2D intersects with the wrong way)

Now it is a matter of identifying my bottleneck! This algorithm should be O (nlog-n), so it should be much faster. My intersectsWith method looks as if it couldn’t be faster, however I will send its code if I am wrong

EDIT: Added IComparable interface

My method of reading intersections of line segments. For readability, some code has been omitted.

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; } } } } 

The main part of my implementation of the algorithm, I understand that List is not the fastest for the algorithm, however I need to index !:

 //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; 

==============================================>

UPDATE: September 12th: Implemented TreeSet from C5, implemented IComparable for my classes and slowed it down even more? Am I still indexing it if that matters?

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

Code using 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(""); 

}

+4
source share
2 answers

First, regarding the intersection of lines: you do not need the actual intersection point, only to find out if they intersect. See http://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/ for an algorithm that does just that.


About List implementation:

In your implementation using List s, you call indexOf in sweepline to find nl . It seeks a sweepline from start to finish. See List(T).IndexOf . If you used the BinarySearch method, this should speed up the search significantly.

The documentation list contains a paragraph called performance considerations. They strongly recommend using a value type that implements IEquatable<T> and IComparable<T> . So your Line2D should be struct and implement these interfaces.

If you follow this advice, extracting the endpoint from the sweepline should be O (log n), which is enough for your purpose, and the memory should be used more efficiently.

Insert and delete are O (n) for lists, as a result of which the base array needs to be moved in memory. A SortedSet has faster insertion and deletion, but I don't quite understand how to find the neighbors of elements in O (log n). Anyone? (See also Why SortedSet <T> .GetViewBetween is not O (log N)? )


In any case, C5 TreeSet should solve this problem.

I looked at the performance of IndexOf and [i] in the user manual , and both are listed as O (log n), so this should not be a problem. This is probably somewhat faster, but nothing more than a fixed factor, to name specific methods for finding neighbors on sweepline, i.e. Successor and Predecessor , which are also O (log n).

So,

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

I do not like that they do not have a method that does not throw an exception, since throwing a throw is very expensive. Your scan line will be quite complete, so I think it’s rare to find it, and calling Successor is the most efficient way. Alternatively, you could continue to call indexOf as it is now, but check to see if it is Count minus one before retrieving [index + 1] and prevent the exception from being thrown:

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

Chapter Two of the manual describes equality and comparison for C5 collections. Here, too, at least you should implement IEquatable<T> and IComparable<T> !

Another thought: you are reporting an algorithm feed of 700,000 lines. Could you start by synchronizing, for example, 1000, 2500, 5000, 10000 lines and see how the algorithm scales for cases when they do not overlap?


How to compare lines on sweepline:

You need to find the natural order for Line2Ds in the Sweepline TreeSet, as the CompareTo method will ask you to compare one Line2D with another. One of Line2Ds is already in the Sweepline TreeSet, the other just met and added.

Your sweepline works from bottom to top, I think:

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

Sweepline running through the line segments

So, let the S1 segment be added to the TreeSet at event 1, and we want to compare it with S2, which is added to event 2, right now.

Line segments may intersect at some point, which would change the order, but the algorithm will check this right after inserting them in the checks above and below. Which might be better called left and right checks, think about it.

In any case ... it would be easiest to compare the lower endpoints of both line segments. Less on the left, more on the right. However, we need to look at the ordering in the sweepline position, and since then they may have changed positions, as in the picture.

So, we need to compare the lower endpoint S2 with the red dot on S1 and see if it lies to the left or right of this point. It lies to the left, so S2 is considered less than S1.

This is usually simpler: if all S1 lies to the left of the lower endpoint S2, S1 is less than S2. If all S1 is to the right of the lower endpoint of S2, S1 is greater than S2.

I think you're looking for an interface type version:

 public class Line2D : IComparable<Line2D> 

Assuming the two BottomY properties, the smallest of the two values ​​Y and BottomX , the X value of the smallest endpoint, a somewhat verified attempt:

 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) ); } 
+3
source

[original answer]

I am not C # , but this should speed things up a bit.

  • less heap interception
  • do not calculate the same thing twice
  • avoid all helper calls if you can (remove functions)

the code:

 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 } 

for the rest of your code, add time dimensions to find what exactly slows down. My hunch is list management ... unnecessary redistributions can slow down a lot.

[edit1]

After doing some research on random line data, I did the following:

  • if there are too many rows throughout the territory than no optimizations are performed
  • if there are smaller lines, the more accelerations for any optimizations
  • brute force T((N*NN)/2) , which is still O(N*N) estimate about 35 hours to process 700K lines (unusable)
  • optimized brute force with partitioning of the domain T((((N/M)^2)-N)/2) - optimization ~O((N/M)^2) where

    • N - maximum number of line numbers
    • M - the number of divisions of the area by any axis. The idea is to check only the lines crossing any region (divide the data set area by M*M squares / rectangles). For 700K lines, it is best to subdivide into 16x16 areas. Measured time:

      0.540 s per 32 thousand lines 1.950 s per 64 thousand lines 7.000s per 128K lines 27.514s per 256K lines

    the estimated operating time is 3.7 min per 700 thousand lines (for lines with a maximum length of ~ 10% of the total area) . I think this is better than your 19 minutes.

  • one more speed is possible using several CPU / core

    Algorithm

    fully parallel for 4 CPU / Core 3.7min/4 -> 56s , or you can transfer it to the GPU ...

My optimized brute force algorithm with regional division O (((((N / M) ^ 2) -N) / 2) - optimization

  • get the size of the used area (xmin,xmax,ymin,ymax) O(N)
  • select unit M best I tried for my random 32K-256K datasets was M=16
  • cycle through the entire area of ​​the unit (evenly divided area of ​​the data set)

    create a list of lines intersecting the actual area of ​​the unit, and check the intersection of all lines in this list. If you do not want to duplicate intersections, then discard all intersections outside the current area

my code (I use BDS2006 C ++ and my own lists, so you will need to map it to your code)

 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 pp[0]=x0; pp[1]=y0; pp[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; } 

Notes:

  • List<T> x; matches T x[] with "unlimited" size
    • x.num; - actual size x[] in Ts not Bytes !!! index = <0,x.num-1>
    • x.add(q); adds q to the list at the end
    • x.num=0; clears the list
    • x.allocate(N); allocate space for N items in the list to avoid moving
  • input List<> is view.lin ... contains points p0,p1 each has double p[2] ... x,y
  • output List<> is view.pnt ... contains double p[2] ... x,y

[Edit2]

In addition, I found out that the best performance of the algorithm is higher when M=12+(N>>15)

  • M - the number of areas of division into an axis
  • N - the number of lines to check
+1
source

Source: https://habr.com/ru/post/1499628/


All Articles