Parallel Plane Sweep  0.1
Shared memory multithreaded version of the plane sweep algorithm
parPlaneSweep.cpp
Go to the documentation of this file.
1 /*
2  * The MIT License (MIT)
3  * Copyright (c) <2016> <Mark McKenney>
4  *
5  * Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
6  *
7  * The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
8  *
9  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
10  * */
11 
12 
13 #include <iomanip>
14 #include <omp.h>
15 #include "parPlaneSweep.h"
16 #include "vectorAlEq.h"
17 #include <limits>
18 #include <algorithm>
19 
23 int binarySearchExists( vector< halfsegment > &region, double x, int hi=-1, int lo=0 )
24 {
25  if( hi < 0 )
26  hi = region.size();
27  while( lo < hi ) {
28  int mid = (lo+hi)/2;
29  if( region[mid].dx < x ) lo = mid+1;
30  else if( region[mid].dx > x ) hi = mid;
31  else return mid;
32  }
33  return -1;
34 }
35 
39 int binarySearchSmallestGreater( vector< halfsegment > &region, double x, int hi=-1, int lo=-1 )
40 {
41  if( hi < 0 )
42  hi = region.size();
43  while( lo < hi ) {
44  int mid = (lo+hi)/2;
45  if( region[mid].dx <= x ) lo = mid+1;
46  else if( region[mid].dx > x ) {
47  if( mid > 0 && region[mid-1].dx <= x ) return mid;
48  hi = mid-1;
49  }
50  //else return mid;
51  }
52  return hi;
53 }
54 
58 bool binarySearchHalfsegment(vector< halfsegment > &region,
59  const halfsegment &h, int & mid, int hi=-1, int lo=0 )
60 {
61  if( hi < 0 )
62  hi = region.size();
63  while( lo < hi ) {
64  mid = (lo+hi)/2;
65  if( region[mid] < h ) lo = mid+1;
66  else if( region[mid] == h ) return true;
67  else hi = mid;
68  }
69  mid = lo;
70  return false;
71 }
72 
83 void findIsoBoundaries( vector<halfsegment> &r1, vector<halfsegment> &r2,
84  vector< double> & isoBounds );
85 
86 
95 void createStrips( vector< halfsegment> & region, vector<double> &isoBounds,
96  vector<halfsegment> & rStrips, vector< int > &stripStopIndex );
97 
98 
99 
109 void partialOverlay( const vector<halfsegment> & r1Strips, const vector<halfsegment> & r2Strips,
110  vector<halfsegment> & result,
111  vector< int > &r1StripStopIndex, vector< int > &r2StripStopIndex,
112  const int stripID );
113 
114 
115 
116 
120 bool findIntersectionPoint( const halfsegment & h1, const halfsegment & h2,
121  double & X, double & Y, bool & colinear );
122 
123 
127 bool breakHsegs( const halfsegment &alSeg, halfsegment & origCurr,
128  vector< halfsegment> & brokenSegs, bool & colinear,
129  const bool includeCurrSegInBrokenSegs );
130 
137 void createFinalOverlay( vector<halfsegment> & finalResult,
138  vector< vector<halfsegment> > &resultStrips,
139  const vector<double> &isoBounds );
140 
141 
147 void insertBrokenSegsToActiveListAndDiscoveredQueue( const vector<halfsegment> & brokenSegs,
148  vector<halfsegment> & result,
149  eventQueue& discoveredSegs,
150  activeListVec& activeList,
151  const double eventX,
152  const double eventY );
153 
157 void parallelOverlay( vector<halfsegment> &r1, vector<halfsegment> &r2, vector<halfsegment> &result,
158  int numStrips, int numWorkerThreads )
159 {
160  vector<halfsegment> r1Strips, r2Strips;
161  vector< double > isoBounds;
162  vector< vector< halfsegment> > resultStrips;
163  vector< int > r1StripStopIndex, r2StripStopIndex;
164  result.clear(); // make sure the result vec is clear
165 
166  // set default parallel values
167  if( numStrips < 0 ) {
168  numStrips = omp_get_num_procs();
169  }
170  if( numWorkerThreads > 0 ) {
171  omp_set_num_threads( numWorkerThreads );
172  }
173 
174 
175  int numIsoBounds = numStrips+1;
176 
177  // set up vectors for split points, return strips
178  for( int i = 0; i < numStrips; i++ ) {
179  resultStrips.push_back( vector<halfsegment>() );
180  }
181  for( int i = 0; i < numIsoBounds; i++ ) {
182  isoBounds.push_back( 0 );
183  }
184 
185  // find split points
186  findIsoBoundaries( r1, r2, isoBounds );
187 
188  // split up the regions at the iso boundaries
189  #pragma omp parallel for
190  for( int i = 0; i < 2; i++ ){
191  if( i == 0 ) createStrips( r1, isoBounds, r1Strips, r1StripStopIndex );
192  else createStrips( r2, isoBounds, r2Strips, r2StripStopIndex );
193  }
194  // {
195  // int count, total = 0;
196  // cout << "Number of segmenents in each strip (r1): , ";
197  // for( int i = 0; i < r1StripStopIndex.size( ); i++ ) {
198  // count = r1StripStopIndex[i];
199  // if( i > 0 ) count -= r1StripStopIndex[i-1];
200  // total += count;
201  // cout << count << ",";
202  // }
203  // cout << "\ntotal segments in all strips:, "<< total << endl;
204  // total = 0;
205  // cout << "Number of segmenents in each strip (r2): , ";
206  // for( int i = 0; i < r2StripStopIndex.size( ); i++ ) {
207  // count = r2StripStopIndex[i];
208  // if( i > 0 ) count -= r2StripStopIndex[i-1];
209  // total += count;
210  // cout << count << ",";
211  // }
212  // cout << "\ntotal segments in all strips:, "<< total << endl;
213  // }
214 
215  // do the actual plane sweeps
216 #pragma omp parallel for schedule(dynamic,1)
217  for( int i = 0; i < numStrips; i++ ) {
218  partialOverlay( r1Strips, r2Strips, resultStrips[i], r1StripStopIndex, r2StripStopIndex, i );
219  }
220 
221  // create the final overlay
222  createFinalOverlay( result,
223  resultStrips,
224  isoBounds );
225 
226 }
227 
228 // call the individual overlay functions.
229 // this function just sets up the function calls to overlayPlaneSweep
233 void partialOverlay( const vector<halfsegment> & r1Strips, const vector<halfsegment> & r2Strips,
234  vector<halfsegment> & result,
235  vector< int > &r1StripStopIndex, vector< int > &r2StripStopIndex,
236  const int stripID )
237 {
238 
239  int r1Start, r2Start, r1Stop, r2Stop;
240  if( stripID == 0 ) {
241  r1Start = r2Start = 0;
242  }
243  else {
244  r1Start = r1StripStopIndex[stripID-1];
245  r2Start = r2StripStopIndex[stripID-1];
246  }
247  r1Stop = r1StripStopIndex[stripID];
248  r2Stop = r2StripStopIndex[stripID];
249 
250  overlayPlaneSweep( &(r1Strips[r1Start]), r1Stop-r1Start,
251  &(r2Strips[r2Start]), r2Stop-r2Start,
252  result);
253 }
254 
258 void overlayPlaneSweep( const halfsegment r1[], int r1Size,
259  const halfsegment r2[], int r2Size,
260  vector<halfsegment>& result )
261 {
262  halfsegment currSeg, maxSeg, tmpSeg, *tmpSegPtr;
263  maxSeg.dx = maxSeg.dy = maxSeg.sx = maxSeg.sy = std::numeric_limits<double>::max();
264  double eventX, eventY;
265  activeListVec activeList;
266  eventQueue discoveredSegs;
267  vector< halfsegment > brokenSegs;
268  bool colinearIntersection;
269  int r1Pos = 0;
270  int r2Pos = 0;
271  int segSource;
272  int segIndex;
273  while( discoveredSegs.size() > 0 || r1Pos < r1Size || r2Pos < r2Size ) {
274  // get the next seg
275  // next seg is the least seg from r1, r2, and the discoveredSeg tree (event queue)
276  currSeg = maxSeg;
277  if( r1Pos < r1Size ) {
278  currSeg = r1[r1Pos];
279  segSource = 1;
280  }
281  if( r2Pos < r2Size && r2[r2Pos]< currSeg ) {
282  currSeg = r2[r2Pos];
283  segSource = 2;
284  }
285  if(discoveredSegs.peek( tmpSeg ) ) {
286  if( tmpSeg < currSeg || tmpSeg == currSeg ){
287  currSeg = tmpSeg;
288  segSource = 3;
289  }
290  }
291  // remove the next seg from its source
292  if( segSource == 3 ) {
293  discoveredSegs.pop();
294  }
295  else if( segSource == 2 ){
296  r2Pos++;
297  }
298  else {
299  r1Pos++;
300  }
301 
302  // set current event point.
303  // the activeList compare function uses eventX as its |param| argument
304  eventX = currSeg.dx;
305  activeList.xVal = eventX;
306  eventY = currSeg.dy;
307  // cerr << "========================="<<endl;
308  // cerr << "segSrc: " << segSource<< " currSeg: " << currSeg << endl;
309  // If curr is a left seg, insert it and check for intersections with neighbors
310  // Else it is a right seg, remove it and check its neighbors for intersections
311  if( currSeg.isLeft( ) ) {
312  // initialize the overlap labels
313  currSeg.ola = currSeg.olb = -1;
314  // insert the left seg
315  // use avl_t_insert so we can get its neighbors.
316 
317  bool dup;
318  halfsegment segInAL;
319  activeList.insert( currSeg, dup, segInAL, segIndex );
320  // record the active list seg for replacement later (if there are intersections with it)
321  //activeList.print();
322  // cerr << "al seg: "<<segInAL << "segIndex: " << segIndex << endl;
323  // if a duplicate is in the active list, we get a pointer to the duplicate. So we need to update labels
324  // if insert is successful, we get a pointer to the inserted item
325  if( dup ) {
326  // cerr << "dup"<<endl;
327  // activeList.print();
328  // We found a duplicate in active list. update labels
329  // NOTE: overlapping segs are a special case for labels.
330  // NOTE: overlap labels are altered in the |breakHsegs()| as well,
331  // but that function is not called here (no need to break up segs
332  // if they are equal)
333  segInAL.ola = currSeg.la;
334  segInAL.olb = currSeg.lb;
335  activeList.replace( segInAL, segInAL, segIndex );
336  }
337  else {
338  bool needToRemoveCurr = false;
339  bool hasAbove, hasBelow;
340  // inserted successfully. Need to get neighbors
341  halfsegment belowSegCopy;
342  halfsegment aboveSegCopy;
343  hasBelow =activeList.getBelow( currSeg, belowSegCopy, segIndex );
344  hasAbove = activeList.getAbove( currSeg, aboveSegCopy, segIndex );
345  // if( hasBelow) cerr << "below: "<< belowSegCopy<<endl;
346 
347  // if( hasAbove ) cerr << "above: " << aboveSegCopy<<endl;
348 
349  // do intersections with above and below. Update currSeg along the way
350  brokenSegs.clear();
351 
352  // We have to deal with the below seg first because currSegs labels get changed
353  // based on the below seg. Once we begin dealing with the above seg, the labels
354  // for the currSeg are already computed and will carry over into those calaculations
355  if( hasBelow ) {
356  // update labels:
357  // NOTE: overlapping segs are a special case for labels.
358  // NOTE: overlap labels are altered in the |breakHsegs()| as well,
359  if( currSeg.regionID != belowSegCopy.regionID ) {
360  // if below seg is from opposing region, set overlap labels
361  if( belowSegCopy.dx != belowSegCopy.sx ) { // if seg is not vertical, use la (label above)
362  currSeg.ola = currSeg.olb = belowSegCopy.la;
363  } else { // if below seg is vertical, use lb (the label to the right)
364  currSeg.ola = currSeg.olb = belowSegCopy.lb;
365  }
366  }
367  else if(currSeg.regionID == belowSegCopy.regionID ) { // if below seg is from same region, just extend overlap labels
368  currSeg.ola = currSeg.olb = belowSegCopy.ola;
369  // commented code is for checking against vertical seg below from same regions
370  // this should never happen for well formed regions based on hseg order
371  // if( belowSegCopy.dx != belowSegCopy.sx ) { // if seg is not vertical, use ola (label above)
372  // currSeg.ola = currSeg.olb = belowSegCopy.ola;
373  // } else { // if below seg is vertical, use lb (the label to the right)
374  // currSeg.ola = currSeg.olb = belowSegCopy.olb;
375  // }
376  }
377 
378  // Labels are now computed
379  // Compute the segment intersections:
380  if(breakHsegs( belowSegCopy, currSeg, brokenSegs, colinearIntersection, false ) ){
381  // cerr << "attempt to erase below: " << belowSegCopy << " " <<segIndex << endl;
382  needToRemoveCurr = true;
383  // remove below seg
384  activeList.erase( belowSegCopy, segIndex-1 );
385  segIndex--;
386  }
387 
388  }
389  // compute intersections with above seg:
390  if( hasAbove ) {
391  if( breakHsegs( aboveSegCopy, currSeg, brokenSegs, colinearIntersection, false ) ) {
392  //cerr << "attempt to erase above: " << aboveSegCopy << " " <<segIndex << endl;
393 
394  needToRemoveCurr = true;
395  // remove above seg
396  activeList.erase( aboveSegCopy, segIndex +1 );
397  }
398  }
399 
400  // Update the seg inserted into the active list this round
401  // currSeg is the result of intersecting that seg with its neighbors
402  // cerr << "about to rep: " << segInAL << endl <<currSeg << endl;
403  activeList.replace( segInAL, currSeg, segIndex );
404  // Insert all the broken up segs into thier various data structures
406  discoveredSegs, activeList,
407  eventX, eventY );
408 
409  }
410  }
411  else {
412  // This is a right halfsegment.
413  // find its brother (left halfsegment) in the active list,
414  // remove it, and check its neighbors for intersections.
415  currSeg = currSeg.getBrother();
416  // cerr <<"brother: " << currSeg << endl;
417  halfsegment theALseg;
418  int segIndex;
419  if( activeList.exists( currSeg, theALseg, segIndex ) ){
420  // cerr << "found brother"<<endl;
421  // We found the halfsegment in the active list.
422  // Its possible we don't find one, since a right halfsegment may be in r1 or r2
423  // whose brother (left halfsegment) was broken due to an intersection
424  result.push_back( theALseg );
425  result.push_back( theALseg.getBrother( ) );
426  // The copy of the seg in the active list has the appropriate labels, so copy it over
427  currSeg = theALseg;
428  // find the neighbors
429  halfsegment belowCopy, aboveCopy;
430  bool isAbove, isBelow;
431  isBelow = activeList.getBelow( currSeg, belowCopy, segIndex );
432  isAbove = activeList.getAbove( currSeg, aboveCopy, segIndex );
433  halfsegment origAbove(aboveCopy);
434  // delete the segment
435  if( isAbove && isBelow ) {
436  // cerr << "check above/below for inters"<<endl;
437  brokenSegs.clear();
438  if(breakHsegs( belowCopy, aboveCopy, brokenSegs, colinearIntersection, true ) ){
439  // if we got here, we are done updating labels, so we
440  // can kill iterators in the active list by deleting
441  //activeList.replace( tmpSeg, aboveCopy);
442  //activeList.erase( currSeg );
443  //alHSegPtr = NULL;
444  // remove below seg and above seg
445  activeList.erase( belowCopy, segIndex-1 );
446  segIndex--;
447  activeList.erase( origAbove, segIndex+1 );
448  activeList.erase( currSeg, segIndex );
449  // add broken segs to output and discovered list
451  discoveredSegs, activeList,
452  eventX, eventY );
453  } else {
454  // delete the currseg
455  activeList.erase( currSeg, segIndex );
456  }
457  } else {
458  // jsut delete the seg from the active list
459  // cerr << "just erase: " << currSeg << endl;
460  activeList.erase( currSeg, segIndex );
461  }
462  // activeList.print();
463 
464 
465 
466  }
467  }
468 
469  }
470  // sort the result
471 
472  std::sort( result.begin(), result.end() );
473  // cerr << "ps result: " <<endl;
474  //for( int i = 0; i <result.size(); i++ )
475  // if( result[i].isLeft() )
476  // cerr << result[i]<<endl;
477 
478 }
479 
480 void insertBrokenSegsToActiveListAndDiscoveredQueue( const vector<halfsegment> & brokenSegs,
481  vector<halfsegment> & result,
482  eventQueue& discoveredSegs,
483  activeListVec& activeList,
484  const double eventX,
485  const double eventY )
486 {
487  for( int i = 0; i < brokenSegs.size(); i++ ){
488  if( (brokenSegs[i].dx != brokenSegs[i].sx && ( brokenSegs[i].dx <= eventX && brokenSegs[i].sx <= eventX))
489  || (brokenSegs[i].dx == brokenSegs[i].sx && ( brokenSegs[i].dy <= eventY && brokenSegs[i].sy <= eventY)) ){
490  // If the seg is behind sweep line, just put it in the output.
491  // The seg is behind teh sweep line if it is not vertical and dx,sx <= eventX
492  // OR the seg is vertical and dy,sy <= eventY
493  result.push_back( brokenSegs[i] );
494  }
495  else if ( !brokenSegs[i].isLeft()
496  || brokenSegs[i].dx > eventX
497  || ( brokenSegs[i].dx == eventX && brokenSegs[i].dy > eventY ) ) {
498  // If the seg is ahead of the sweep line, or a right halfegment,
499  // it goes in the event queue (discovered list)
500  // The seg is ahead of the sweep line if dx > eventX or dx = eventx and dy > eventY
501  discoveredSegs.insert( brokenSegs[i] );
502  }
503  else {
504  // If we get here, the seg spans the sweep line and is a left halfsegment
505  // The seg goes back into the active list.
506  bool tmp;
507  halfsegment tmph;
508  int tmpi;
509  activeList.insert( brokenSegs[i], tmp,tmph,tmpi );
510  }
511  }
512 }
513 
514 bool breakHsegs( const halfsegment &alSeg, halfsegment & origCurr,
515  vector< halfsegment> & brokenSegs, bool & colinear,
516  const bool includeCurrSegInBrokenSegs)
517 {
518  halfsegment tmpSeg;
519  halfsegment curr = origCurr; // copy
520  const halfsegment h2 = alSeg; // rename
521  bool foundIntersection;
522  // get the intersecion point
523  double X,Y;
524  if( foundIntersection = findIntersectionPoint( h2, curr, X, Y, colinear ) ) {
525  if( colinear ) {
526  // If the segs are colinear, their intersection can have at most 3 components
527  // 1) one seg begins to the left (or below) the other. (non overlapping part)
528  // 2) the portion of the segs that overlap
529  // 3) one seg ends to the right (or above) the other. (non overlapping part)
530  // at most, all three portions exist. At the least, portion (2) will exist
531  if( h2.dx < curr.dx || h2.dy < curr.dy ) { // build the first part (1)
532  tmpSeg = h2;
533  tmpSeg.sx = curr.dx;
534  tmpSeg.sy = curr.dy;
535  brokenSegs.push_back( tmpSeg ); // THIS IS NEW Neighbor SEG from AVL tree
536  brokenSegs.push_back( tmpSeg.getBrother() );
537  }
538  // build the middle part (2)
539  tmpSeg = curr;
540  if( curr.sx > h2.sx || (curr.sx == h2.sx && curr.sy > h2.sy )) {
541  // curr extends beyond h2
542  tmpSeg.sx = h2.sx;
543  tmpSeg.sy = h2.sy;
544  }
545  tmpSeg.ola = h2.la; // Set the overlapping labels.
546  tmpSeg.olb = h2.lb;
547  brokenSegs.push_back( tmpSeg.getBrother() );
548  if( includeCurrSegInBrokenSegs ) { // don't need to return currSeg for inserting left hseg
549  brokenSegs.push_back( tmpSeg );
550  }
551  origCurr = tmpSeg; // UPDATE |origCurr| so that the overlap labels get updated
552  // it is pass by reference
553  // build last part (3)
554  if( curr.sx != h2.sx || curr.sy != h2.sy ) {
555  // h2 extends past curr
556  tmpSeg = h2;
557  tmpSeg.dx = curr.sx;
558  tmpSeg.dy = curr.sy;
559  if( curr.sx > h2.sx || (curr.sx == h2.sx && curr.sy > h2.sy )) {
560  // curr extends past h2
561  tmpSeg = curr;
562  tmpSeg.dx = h2.sx;
563  tmpSeg.dy = h2.sy;
564  }
565  brokenSegs.push_back( tmpSeg ); // some future seg
566  brokenSegs.push_back( tmpSeg.getBrother() );
567  }
568  }
569  else {
570  // regular intersection
571  // split up curr
572  if( (X == curr.dx && Y == curr.dy ) ||
573  (X == curr.sx && Y == curr.sy ) ) {
574  // If this is an end point intersection, do not break the seg
575  // we do not remove the curr seg from the active list, so no need
576  // to add to broken segs to get reinserted. Active list curr seg
577  // gets updated in the plane sweep function after all intersections
578  // with currSegs neighbors have been computed
579  }
580  else {
581  // If the intersection is on the interior of this seg, break it up
582  tmpSeg = curr;
583  tmpSeg.sx = X;
584  tmpSeg.sy = Y;
585  brokenSegs.push_back( tmpSeg.getBrother() );
586  if( includeCurrSegInBrokenSegs ) { // don't need to return currSeg for inserting left hseg
587  brokenSegs.push_back( tmpSeg );
588  }
589  origCurr = tmpSeg; // UPDATE |origCurr|. it is passed by reference
590  tmpSeg = curr;
591  tmpSeg.dx = X;
592  tmpSeg.dy = Y;
593  brokenSegs.push_back( tmpSeg );
594  brokenSegs.push_back( tmpSeg.getBrother() );
595  }
596  // split up h2. Identical code to above. see those comments
597  if( (X == h2.dx && Y == h2.dy ) ||
598  (X == h2.sx && Y == h2.sy ) ) {
599  // If this is an end point intersection, do not break the seg
600  // we will reinsert this seg, but its brother remains the same,
601  // so we don't need to put the brother in
602  // reinsertion is just for convienience
603  // we remove the seg in the plane sweep function so
604  // we don't have to keep track of the aboves belows.
605  // In the case of colinears, seg must be removed anyway, so we just
606  // always remove the above/below segs and reinsert them if needed.
607  brokenSegs.push_back( h2 );
608  }
609  else {
610  // interior intersection, split
611  tmpSeg = h2;
612  tmpSeg.sx = X;
613  tmpSeg.sy = Y;
614  brokenSegs.push_back( tmpSeg ); // THIS IS THE UPDATED H2
615  brokenSegs.push_back( tmpSeg.getBrother() );
616  tmpSeg = h2;
617  tmpSeg.dx = X;
618  tmpSeg.dy = Y;
619  brokenSegs.push_back( tmpSeg );
620  brokenSegs.push_back( tmpSeg.getBrother() );
621  }
622  }
623  }
624  // cerr << "------- broken segs -----" << endl;
625  // for( int i =0; i < brokenSegs.size(); i++ )
626  // cerr << brokenSegs[i]<<endl;
627  // cerr << "------- -----" << endl;
628 
629  return foundIntersection;
630 }
631 
632 bool findIntersectionPoint( const halfsegment & h1, const halfsegment & h2, double & X, double & Y, bool & colinear )
633 {
634  // if colinear, intersection point is h2 dominating (since h2 will be curr seg from overlay)
635  colinear = false;
636  if( h1.colinear( h2 ) ) {
637  X = h2.dx;
638  Y = h2.dy;
639  colinear = true;
640  return true;
641  }
642  // if they share an end point, there is nothing to do
643  if( (h1.dx == h2.dx && h1.dy == h2.dy ) ||
644  (h1.sx == h2.sx && h1.sy == h2.sy ) ) {
645  X = std::numeric_limits<double>::max();
646  Y = std::numeric_limits<double>::max();
647  return false;
648  }
649 
650  // find intersection point
651  double x1 = h1.dx;
652  double y1 = h1.dy;
653  double x2 = h1.sx;
654  double y2 = h1.sy;
655  double x3 = h2.dx;
656  double y3 = h2.dy;
657  double x4 = h2.sx;
658  double y4 = h2.sy;
659 
660  double denom = ((y4-y3)*(x2-x1)) - ((x4-x3)*(y2-y1));
661  double ua = ((x4-x3)*(y1-y3)) - ((y4-y3)*(x1-x3));
662  double ub = ((x2-x1)*(y1-y3)) - ((y2-y1)*(x1-x3));
663 
664  ua = ua/denom;
665  ub = ub/denom;
666  // if ua and ub are between 0 and 1 inclusive, we have an intersection
667  // in at least 1 interior.
668  // end point intersections are handled above
669  if( 0.0 <= ua && ua <= 1.0 && 0.0 <= ub && ub <= 1.0 ) {
670  X = x1 + ( ua*(x2-x1) );
671  Y = y1 + ( ua*(y2-y1) );
672  return true;
673  }
674  else {
675  X = std::numeric_limits<double>::max();
676  Y = std::numeric_limits<double>::max();
677  }
678  return false;
679 }
680 
681 bool hsegIDSort( const halfsegment & h1, const halfsegment & h2 )
682 {
683  return h1.stripID < h2.stripID || (h1.stripID == h2.stripID && h1 < h2 );
684 }
685 
686 void createStrips( vector< halfsegment> & region, vector<double> &isoBounds,
687  vector<halfsegment> & rStrips, vector< int > &stripStopIndex )
688 {
689  halfsegment workSeg;
690  int startBound = 0;
691  // grab a seg, break it on each strip that it crosses, put it in the strips
692  for( int i = 0; i< region.size(); i++ ) {
693  // only need to worry about lefties
694  if( region[i].isLeft( ) ) {
695  workSeg = region[i];
696  for( int j = startBound; j < isoBounds.size()-1; j++ ){
697  if( workSeg.dx > isoBounds[j+1] ) {
698  // we are done with this strip (hseg ordering)
699  startBound++;
700  continue;
701  }
702  // check if we cross the boundary.
703  // remember, we won't have any seg end on a boundary unless we split it
704  else if( workSeg.dx >= isoBounds[j] && workSeg.sx < isoBounds[j+1] ) {
705  workSeg.stripID = j;
706  rStrips.push_back( workSeg );
707  rStrips.push_back( workSeg.getBrother( ) );
708  break; // done with this seg
709  }
710  // otherwise, we cross a boundary
711  else {
712  // get the y value at iso bound
713  halfsegment lhs = workSeg;
714  // make the segs split at isoBounds[j]
715  lhs.sy = workSeg.dy = workSeg.getYvalAtX( isoBounds[j+1] );
716  lhs.sx = workSeg.dx = isoBounds[j+1];
717  lhs.stripID = j;
718  rStrips.push_back( lhs );
719  rStrips.push_back( lhs.getBrother() );
720  }
721  }
722  }
723  }
724  // deallocate input segments (to save memory!)
725  // use the swap to local var trick!
726  //{
727  // vector< halfsegment > localOne;
728  // region.swap(localOne);
729  //}
730 
731  // sort the strips by stripID, then by hseg order
732 
733  std::sort( rStrips.begin(), rStrips.end(), hsegIDSort);
734  // find the startIndex for each strip
735  for( int i = 0; i < isoBounds.size()-1; i++ )
736  stripStopIndex.push_back( std::numeric_limits<int>::min() );
737  for( int i = 0; i < rStrips.size(); i++ ){
738  if( i >= stripStopIndex[ rStrips[i].stripID ] )
739  stripStopIndex[ rStrips[i].stripID ] = i+1;
740  }
741  // remove any remaining min vals
742  int prevVal = 0;
743  for( int i = 0; i < stripStopIndex.size(); i++ ) {
744  if( stripStopIndex[i] == std::numeric_limits<int>::min() ) {
745  stripStopIndex[i] = prevVal;
746  }
747  prevVal = stripStopIndex[i];
748  }
749 
750 #ifdef DEBUG_PRINT
751 #pragma omp critical
752  {
753  /* cerr << "strips: "<<endl;
754  for( int i = 0; i < rStrips.size(); i++ ) {
755  cerr <<"[(" << rStrips[i].dx << "," << rStrips[i].dy << ")(" << rStrips[i].sx << "," << rStrips[i].sy << ") " << rStrips[i].stripID << ", " << rStrips[i].la << ", " << rStrips[i].lb << "]" << endl;
756  }*/
757  cerr<< "stopIndex: " ;
758  for( int i = 0; i < stripStopIndex.size(); i++ )
759  cerr << stripStopIndex[i] << " ";
760  cerr << endl;
761  }
762 #endif
763 }
764 
765 void findIsoBoundaries( vector<halfsegment> &r1, vector<halfsegment> &r2,
766  vector< double> & isoBounds )
767 {
768 
769  // set extrema for isobounds
770  isoBounds[0] = std::numeric_limits<double>::max() *-1;
771  isoBounds[ isoBounds.size()-1 ] = std::numeric_limits<double>::max();
772 
773  // find min/max X
774  double minX = std::numeric_limits<double>::max();
775  double maxX = std::numeric_limits<double>::max()*-1;
776  if( isoBounds.size() == 2){
777  return;
778  }
779  for( int i = 0; i < r1.size(); i++ ){
780  if( r1[i].dx < minX ) minX = r1[i].dx;
781  if( r1[i].dx > maxX ) maxX = r1[i].dx;
782  if( r1[i].sx < minX ) minX = r1[i].sx;
783  if( r1[i].sx > maxX ) maxX = r1[i].sx;
784  }
785  for( int i = 0; i < r2.size(); i++ ){
786  if( r2[i].dx < minX ) minX = r2[i].dx;
787  if( r2[i].dx > maxX ) maxX = r2[i].dx;
788  if( r2[i].sx < minX ) minX = r2[i].sx;
789  if( r2[i].sx > maxX ) maxX = r2[i].sx;
790  }
791  // calc the middle iso bounds (not the extrema).
792  // need numIsoBounds-2 values spaced evenly between minX and maxX
793  double prevIsoVal = minX;
794  double stripWidth = (maxX-minX) / (isoBounds.size()-1);
795  for( int i = 1; i < isoBounds.size()-1; i++ ) {
796  prevIsoVal += stripWidth;
797  isoBounds[i] = prevIsoVal;
798  }
799 
800 #ifdef DEBUG_PRINT
801  cerr<< "iso Bounds: " << endl;
802  cerr << minX <<","<< maxX<<","<< r1.size() << ","<<r2.size() <<endl;
803  for( int i= 0; i < isoBounds.size(); i++ ){
804  cerr << isoBounds[i] << ", ";
805  }
806  cerr << endl;
807 #endif
808 
809  // make sure we don't have an iso boundary on an endpoint
810  for( int i = 1; i < isoBounds.size()-1; i++ ) {
811  int r1Index, r2Index;
812  r1Index = binarySearchExists( r1, isoBounds[i] );
813  r2Index = binarySearchExists( r2, isoBounds[i] );
814  // move forward to find highest segwith same domX
815  if( r1Index >= 0 || r2Index >= 0 ){
816  // get smallest greater
817  r1Index = binarySearchSmallestGreater( r1, isoBounds[i] );
818  r2Index = binarySearchSmallestGreater( r2, isoBounds[i] );
819  // cerr << "riInd, r2Ind " << r1Index << ", " << r2Index << endl;
820  double xVal = isoBounds[i+1];
821  if( r1Index < r1.size() && r1[r1Index].dx < xVal ) {
822  xVal = r1[r1Index].dx;
823  }
824  if( r2Index < r2.size() && r2[r2Index].dx < xVal ) {
825  xVal = r2[r2Index].dx;
826  }
827  if( xVal == std::numeric_limits<double>::max() ){
828  cerr << "Did not find xVal in splits" << endl;
829  exit(-1 );
830  }
831  // find point between isoX and xVal
832  isoBounds[i] = (isoBounds[i]+xVal)/2.0;
833  }
834  }
835 
836 #ifdef DEBUG_PRINT
837  cerr<< "iso Bounds: " << endl;
838  for( int i= 0; i < isoBounds.size(); i++ ){
839  cerr << isoBounds[i] << ", ";
840  }
841  cerr << endl;
842 #endif
843 
844 
845 }
846 
847 
848 void createFinalOverlay( vector<halfsegment> & finalResult,
849  vector< vector<halfsegment> > &resultStrips,
850  const vector<double> &isoBounds )
851 {
852  halfsegment curr;
853  int broIndex;
854  int currBound;
855  int currStrip;
856  int currIndex;
857  // each strip is sorted. Join segs with artifical breaks
858  // easiest way: for each isobound
859  // 1) grab a non-invalidated seg
860  // 2) if it ends on an iso bound, check if it is the only one at that point
861  // 3) if so, find the seg and its brother from the next strip, join, invalidate all others
862  // keep the earliest left seg as the valid one
863  // 4) if the new seg ends on the next iso bound, repeat, building the seg across all bounds.
864  // 5) once seg is finished, put it and its brother in results
865  for( int i = 0; i <resultStrips.size(); i++ ) {
866  for( int j = 0; j < resultStrips[i].size(); j++ ) {
867  curr = resultStrips[i][j];
868  // only process left hsegs that are valid;
869  if( curr.isLeft() && curr.la != curr.lb ) {
870  // if the left end point ends on an iso boundary, merge it
871  currIndex = j;
872  currBound = i+1;
873  currStrip = i;
874  if( curr.sx == isoBounds[currBound] ) {
875  // find its brother
876  if( ! binarySearchHalfsegment(resultStrips[ currStrip ], curr.getBrother(), broIndex ) ) {
877  cerr << "did not find brother 1" <<endl; exit( -1 );
878  }
879  }
880  bool invalidateLast=false;
881  while( curr.sx == isoBounds[currBound] ) {
882  // check the brother's neighbors
883  if( (broIndex+1 < resultStrips[currStrip].size()
884  && resultStrips[currStrip][broIndex+1].dx == curr.sx
885  && resultStrips[currStrip][broIndex+1].dy == curr.sy )
886  || (broIndex > 0
887  && resultStrips[currStrip][broIndex-1].dx == curr.sx
888  && resultStrips[currStrip][broIndex-1].dy == curr.sy ) ) {
889  // multiple segs cross here. We are done with this seg
890  break;
891  } else {
892  invalidateLast = true; // record that we need to invalidate the last part
893  // we have to join this seg with its counterpart in the next strip
894  // invalidate the brother
895  resultStrips[currStrip][broIndex].la = resultStrips[currStrip][broIndex].lb = -1;
896  // invalidate the curr
897  resultStrips[currStrip][currIndex].la = resultStrips[currStrip][currIndex].lb = -1;
898  // find the seg in the next strip
899  binarySearchHalfsegment(resultStrips[ currStrip+1 ], resultStrips[currStrip][broIndex], currIndex );
900  // find its brother
901  binarySearchHalfsegment(resultStrips[ currStrip+1 ], resultStrips[currStrip+1][currIndex].getBrother(), broIndex );
902  // update curr with that segs sub point
903  curr.sx = resultStrips[ currStrip+1 ][currIndex].sx;
904  curr.sy = resultStrips[ currStrip+1 ][currIndex].sy;
905  currStrip++; // we are now in the next strip
906  currBound++; // we are now looking at the next bound
907  }
908  }
909  if( invalidateLast ) {
910  // invalidate the brother
911  resultStrips[currStrip][broIndex].la = resultStrips[currStrip][broIndex].lb = -1;
912  // invalidate the curr
913  resultStrips[currStrip][currIndex].la = resultStrips[currStrip][currIndex].lb = -1;
914 
915  }
916  // add the seg and its brother to the output
917  finalResult.push_back( curr );
918  finalResult.push_back( curr.getBrother() );
919  }
920  }
921  }
922 
923 }
924 
925 
int ola
overlap labels.
Definition: parPlaneSweep.h:87
bool colinear(const halfsegment &rhs) const
A vector-based implementation of a plane sweep event queue.
Definition: vectorAlEq.h:42
holds a single halfsegment, along with labeling information and information indicating if the opposin...
Definition: parPlaneSweep.h:76
bool hsegIDSort(const halfsegment &h1, const halfsegment &h2)
int binarySearchExists(vector< halfsegment > &region, double x, int hi=-1, int lo=0)
int regionID
The region the seg belongs to.
Definition: parPlaneSweep.h:85
void createStrips(vector< halfsegment > &region, vector< double > &isoBounds, vector< halfsegment > &rStrips, vector< int > &stripStopIndex)
bool findIntersectionPoint(const halfsegment &h1, const halfsegment &h2, double &X, double &Y, bool &colinear)
void partialOverlay(const vector< halfsegment > &r1Strips, const vector< halfsegment > &r2Strips, vector< halfsegment > &result, vector< int > &r1StripStopIndex, vector< int > &r2StripStopIndex, const int stripID)
void parallelOverlay(vector< halfsegment > &r1, vector< halfsegment > &r2, vector< halfsegment > &result, int numStrips, int numWorkerThreads)
double getYvalAtX(const double x) const
void insertBrokenSegsToActiveListAndDiscoveredQueue(const vector< halfsegment > &brokenSegs, vector< halfsegment > &result, avl_table *discoveredSegs, avl_table *activeList, const double eventX, const double eventY)
bool pop()
Definition: vectorAlEq.h:89
bool getBelow(const halfsegment &h1, halfsegment &theBelow, const int index)
Definition: vectorAlEq.h:391
void createFinalOverlay(vector< halfsegment > &finalResult, vector< vector< halfsegment > > &resultStrips, const vector< double > &isoBounds)
bool exists(const halfsegment &h1, halfsegment &theCopy, int &index)
Definition: vectorAlEq.h:280
int binarySearchSmallestGreater(vector< halfsegment > &region, double x, int hi=-1, int lo=-1)
int stripID
strip ID
Definition: parPlaneSweep.h:83
bool isLeft() const
double dx
dominating and submissive points
Definition: parPlaneSweep.h:79
void overlayPlaneSweep(const halfsegment r1[], int r1Size, const halfsegment r2[], int r2Size, vector< halfsegment > &result)
bool binarySearchHalfsegment(vector< halfsegment > &region, const halfsegment &h, int &mid, int hi=-1, int lo=0)
int size()
Definition: vectorAlEq.h:99
A vector-based implementation of a plane sweep active list.
Definition: vectorAlEq.h:130
int la
label above, label below
Definition: parPlaneSweep.h:81
void replace(const halfsegment &h1, const halfsegment &newH1)
Definition: vectorAlEq.h:309
void insert(const halfsegment &h1, bool &duplicate, halfsegment &theDup, int &segIndex)
Definition: vectorAlEq.h:244
void findIsoBoundaries(vector< halfsegment > &r1, vector< halfsegment > &r2, vector< double > &isoBounds)
bool breakHsegs(const halfsegment &alSeg, halfsegment &origCurr, vector< halfsegment > &brokenSegs, bool &colinear, const bool includeCurrSegInBrokenSegs)
bool getAbove(const halfsegment &h1, halfsegment &theAbove, const int index)
Definition: vectorAlEq.h:350
void insert(const halfsegment &h1)
Definition: vectorAlEq.h:55
double xVal
The current position of the sweep line.
Definition: vectorAlEq.h:137
halfsegment getBrother() const
void erase(const halfsegment &h1, const int index)
Definition: vectorAlEq.h:423
bool peek(halfsegment &h1)
Definition: vectorAlEq.h:77