29 if( region[mid].dx < x ) lo = mid+1;
30 else if( region[mid].dx > x ) hi = mid;
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;
59 const halfsegment &h,
int & mid,
int hi=-1,
int lo=0 )
65 if( region[mid] < h ) lo = mid+1;
66 else if( region[mid] == h )
return true;
84 vector< double> & isoBounds );
95 void createStrips( vector< halfsegment> & region, vector<double> &isoBounds,
96 vector<halfsegment> & rStrips, vector< int > &stripStopIndex );
109 void partialOverlay(
const vector<halfsegment> & r1Strips,
const vector<halfsegment> & r2Strips,
110 vector<halfsegment> & result,
111 vector< int > &r1StripStopIndex, vector< int > &r2StripStopIndex,
121 double & X,
double & Y,
bool & colinear );
128 vector< halfsegment> & brokenSegs,
bool & colinear,
129 const bool includeCurrSegInBrokenSegs );
138 vector< vector<halfsegment> > &resultStrips,
139 const vector<double> &isoBounds );
148 vector<halfsegment> & result,
152 const double eventY );
157 void parallelOverlay( vector<halfsegment> &r1, vector<halfsegment> &r2, vector<halfsegment> &result,
158 int numStrips,
int numWorkerThreads )
160 vector<halfsegment> r1Strips, r2Strips;
161 vector< double > isoBounds;
162 vector< vector< halfsegment> > resultStrips;
163 vector< int > r1StripStopIndex, r2StripStopIndex;
167 if( numStrips < 0 ) {
168 numStrips = omp_get_num_procs();
170 if( numWorkerThreads > 0 ) {
171 omp_set_num_threads( numWorkerThreads );
175 int numIsoBounds = numStrips+1;
178 for(
int i = 0; i < numStrips; i++ ) {
179 resultStrips.push_back( vector<halfsegment>() );
181 for(
int i = 0; i < numIsoBounds; i++ ) {
182 isoBounds.push_back( 0 );
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 );
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 );
233 void partialOverlay(
const vector<halfsegment> & r1Strips,
const vector<halfsegment> & r2Strips,
234 vector<halfsegment> & result,
235 vector< int > &r1StripStopIndex, vector< int > &r2StripStopIndex,
239 int r1Start, r2Start, r1Stop, r2Stop;
241 r1Start = r2Start = 0;
244 r1Start = r1StripStopIndex[stripID-1];
245 r2Start = r2StripStopIndex[stripID-1];
247 r1Stop = r1StripStopIndex[stripID];
248 r2Stop = r2StripStopIndex[stripID];
251 &(r2Strips[r2Start]), r2Stop-r2Start,
260 vector<halfsegment>& result )
263 maxSeg.
dx = maxSeg.
dy = maxSeg.
sx = maxSeg.
sy = std::numeric_limits<double>::max();
264 double eventX, eventY;
267 vector< halfsegment > brokenSegs;
268 bool colinearIntersection;
273 while( discoveredSegs.
size() > 0 || r1Pos < r1Size || r2Pos < r2Size ) {
277 if( r1Pos < r1Size ) {
281 if( r2Pos < r2Size && r2[r2Pos]< currSeg ) {
285 if(discoveredSegs.
peek( tmpSeg ) ) {
286 if( tmpSeg < currSeg || tmpSeg == currSeg ){
292 if( segSource == 3 ) {
293 discoveredSegs.
pop();
295 else if( segSource == 2 ){
305 activeList.
xVal = eventX;
313 currSeg.
ola = currSeg.
olb = -1;
319 activeList.
insert( currSeg, dup, segInAL, segIndex );
333 segInAL.
ola = currSeg.
la;
334 segInAL.
olb = currSeg.
lb;
335 activeList.
replace( segInAL, segInAL, segIndex );
338 bool needToRemoveCurr =
false;
339 bool hasAbove, hasBelow;
343 hasBelow =activeList.
getBelow( currSeg, belowSegCopy, segIndex );
344 hasAbove = activeList.
getAbove( currSeg, aboveSegCopy, segIndex );
361 if( belowSegCopy.
dx != belowSegCopy.
sx ) {
362 currSeg.
ola = currSeg.
olb = belowSegCopy.
la;
364 currSeg.
ola = currSeg.
olb = belowSegCopy.
lb;
368 currSeg.
ola = currSeg.
olb = belowSegCopy.
ola;
380 if(
breakHsegs( belowSegCopy, currSeg, brokenSegs, colinearIntersection,
false ) ){
382 needToRemoveCurr =
true;
384 activeList.
erase( belowSegCopy, segIndex-1 );
391 if(
breakHsegs( aboveSegCopy, currSeg, brokenSegs, colinearIntersection,
false ) ) {
394 needToRemoveCurr =
true;
396 activeList.
erase( aboveSegCopy, segIndex +1 );
403 activeList.
replace( segInAL, currSeg, segIndex );
406 discoveredSegs, activeList,
419 if( activeList.
exists( currSeg, theALseg, segIndex ) ){
424 result.push_back( theALseg );
430 bool isAbove, isBelow;
431 isBelow = activeList.
getBelow( currSeg, belowCopy, segIndex );
432 isAbove = activeList.
getAbove( currSeg, aboveCopy, segIndex );
435 if( isAbove && isBelow ) {
438 if(
breakHsegs( belowCopy, aboveCopy, brokenSegs, colinearIntersection,
true ) ){
445 activeList.
erase( belowCopy, segIndex-1 );
447 activeList.
erase( origAbove, segIndex+1 );
448 activeList.
erase( currSeg, segIndex );
451 discoveredSegs, activeList,
455 activeList.
erase( currSeg, segIndex );
460 activeList.
erase( currSeg, segIndex );
472 std::sort( result.begin(), result.end() );
481 vector<halfsegment> & result,
485 const double eventY )
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)) ){
493 result.push_back( brokenSegs[i] );
495 else if ( !brokenSegs[i].isLeft()
496 || brokenSegs[i].dx > eventX
497 || ( brokenSegs[i].dx == eventX && brokenSegs[i].dy > eventY ) ) {
501 discoveredSegs.
insert( brokenSegs[i] );
509 activeList.
insert( brokenSegs[i], tmp,tmph,tmpi );
515 vector< halfsegment> & brokenSegs,
bool & colinear,
516 const bool includeCurrSegInBrokenSegs)
521 bool foundIntersection;
531 if( h2.
dx < curr.
dx || h2.
dy < curr.
dy ) {
535 brokenSegs.push_back( tmpSeg );
540 if( curr.
sx > h2.
sx || (curr.
sx == h2.
sx && curr.
sy > h2.
sy )) {
548 if( includeCurrSegInBrokenSegs ) {
549 brokenSegs.push_back( tmpSeg );
554 if( curr.
sx != h2.
sx || curr.
sy != h2.
sy ) {
559 if( curr.
sx > h2.
sx || (curr.
sx == h2.
sx && curr.
sy > h2.
sy )) {
565 brokenSegs.push_back( tmpSeg );
572 if( (X == curr.
dx && Y == curr.
dy ) ||
573 (X == curr.
sx && Y == curr.
sy ) ) {
586 if( includeCurrSegInBrokenSegs ) {
587 brokenSegs.push_back( tmpSeg );
593 brokenSegs.push_back( tmpSeg );
597 if( (X == h2.
dx && Y == h2.
dy ) ||
598 (X == h2.
sx && Y == h2.
sy ) ) {
607 brokenSegs.push_back( h2 );
614 brokenSegs.push_back( tmpSeg );
619 brokenSegs.push_back( tmpSeg );
629 return foundIntersection;
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();
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));
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) );
675 X = std::numeric_limits<double>::max();
676 Y = std::numeric_limits<double>::max();
686 void createStrips( vector< halfsegment> & region, vector<double> &isoBounds,
687 vector<halfsegment> & rStrips, vector< int > &stripStopIndex )
692 for(
int i = 0; i< region.size(); i++ ) {
694 if( region[i].isLeft( ) ) {
696 for(
int j = startBound; j < isoBounds.size()-1; j++ ){
697 if( workSeg.
dx > isoBounds[j+1] ) {
704 else if( workSeg.
dx >= isoBounds[j] && workSeg.
sx < isoBounds[j+1] ) {
706 rStrips.push_back( workSeg );
716 lhs.
sx = workSeg.
dx = isoBounds[j+1];
718 rStrips.push_back( lhs );
733 std::sort( rStrips.begin(), rStrips.end(),
hsegIDSort);
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;
743 for(
int i = 0; i < stripStopIndex.size(); i++ ) {
744 if( stripStopIndex[i] == std::numeric_limits<int>::min() ) {
745 stripStopIndex[i] = prevVal;
747 prevVal = stripStopIndex[i];
757 cerr<<
"stopIndex: " ;
758 for(
int i = 0; i < stripStopIndex.size(); i++ )
759 cerr << stripStopIndex[i] <<
" ";
766 vector< double> & isoBounds )
770 isoBounds[0] = std::numeric_limits<double>::max() *-1;
771 isoBounds[ isoBounds.size()-1 ] = std::numeric_limits<double>::max();
774 double minX = std::numeric_limits<double>::max();
775 double maxX = std::numeric_limits<double>::max()*-1;
776 if( isoBounds.size() == 2){
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;
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;
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;
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] <<
", ";
810 for(
int i = 1; i < isoBounds.size()-1; i++ ) {
811 int r1Index, r2Index;
815 if( r1Index >= 0 || r2Index >= 0 ){
820 double xVal = isoBounds[i+1];
821 if( r1Index < r1.size() && r1[r1Index].dx < xVal ) {
822 xVal = r1[r1Index].dx;
824 if( r2Index < r2.size() && r2[r2Index].dx < xVal ) {
825 xVal = r2[r2Index].dx;
827 if( xVal == std::numeric_limits<double>::max() ){
828 cerr <<
"Did not find xVal in splits" << endl;
832 isoBounds[i] = (isoBounds[i]+xVal)/2.0;
837 cerr<<
"iso Bounds: " << endl;
838 for(
int i= 0; i < isoBounds.size(); i++ ){
839 cerr << isoBounds[i] <<
", ";
849 vector< vector<halfsegment> > &resultStrips,
850 const vector<double> &isoBounds )
865 for(
int i = 0; i <resultStrips.size(); i++ ) {
866 for(
int j = 0; j < resultStrips[i].size(); j++ ) {
867 curr = resultStrips[i][j];
874 if( curr.
sx == isoBounds[currBound] ) {
877 cerr <<
"did not find brother 1" <<endl; exit( -1 );
880 bool invalidateLast=
false;
881 while( curr.
sx == isoBounds[currBound] ) {
883 if( (broIndex+1 < resultStrips[currStrip].size()
884 && resultStrips[currStrip][broIndex+1].dx == curr.
sx
885 && resultStrips[currStrip][broIndex+1].dy == curr.
sy )
887 && resultStrips[currStrip][broIndex-1].dx == curr.
sx
888 && resultStrips[currStrip][broIndex-1].dy == curr.
sy ) ) {
892 invalidateLast =
true;
895 resultStrips[currStrip][broIndex].la = resultStrips[currStrip][broIndex].lb = -1;
897 resultStrips[currStrip][currIndex].la = resultStrips[currStrip][currIndex].lb = -1;
901 binarySearchHalfsegment(resultStrips[ currStrip+1 ], resultStrips[currStrip+1][currIndex].getBrother(), broIndex );
903 curr.
sx = resultStrips[ currStrip+1 ][currIndex].sx;
904 curr.
sy = resultStrips[ currStrip+1 ][currIndex].sy;
909 if( invalidateLast ) {
911 resultStrips[currStrip][broIndex].la = resultStrips[currStrip][broIndex].lb = -1;
913 resultStrips[currStrip][currIndex].la = resultStrips[currStrip][currIndex].lb = -1;
917 finalResult.push_back( curr );
bool colinear(const halfsegment &rhs) const
A vector-based implementation of a plane sweep event queue.
holds a single halfsegment, along with labeling information and information indicating if the opposin...
bool hsegIDSort(const halfsegment &h1, const halfsegment &h2)
int binarySearchExists(vector< halfsegment > ®ion, double x, int hi=-1, int lo=0)
int regionID
The region the seg belongs to.
void createStrips(vector< halfsegment > ®ion, 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 getBelow(const halfsegment &h1, halfsegment &theBelow, const int index)
void createFinalOverlay(vector< halfsegment > &finalResult, vector< vector< halfsegment > > &resultStrips, const vector< double > &isoBounds)
bool exists(const halfsegment &h1, halfsegment &theCopy, int &index)
int binarySearchSmallestGreater(vector< halfsegment > ®ion, double x, int hi=-1, int lo=-1)
double dx
dominating and submissive points
void overlayPlaneSweep(const halfsegment r1[], int r1Size, const halfsegment r2[], int r2Size, vector< halfsegment > &result)
bool binarySearchHalfsegment(vector< halfsegment > ®ion, const halfsegment &h, int &mid, int hi=-1, int lo=0)
A vector-based implementation of a plane sweep active list.
int la
label above, label below
void replace(const halfsegment &h1, const halfsegment &newH1)
void insert(const halfsegment &h1, bool &duplicate, halfsegment &theDup, int &segIndex)
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)
void insert(const halfsegment &h1)
double xVal
The current position of the sweep line.
halfsegment getBrother() const
void erase(const halfsegment &h1, const int index)
bool peek(halfsegment &h1)