Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members  

gugr_basics.cpp

Go to the documentation of this file.
00001 /* LIBGUL - Geometry Utility Library
00002  * Copyright (C) 1998-1999 Norbert Irmer
00003  *
00004  * This library is free software; you can redistribute it and/or
00005  * modify it under the terms of the GNU Library General Public
00006  * License as published by the Free Software Foundation; either
00007  * version 2 of the License, or (at your option) any later version.
00008  *
00009  * This library is distributed in the hope that it will be useful,
00010  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * Library General Public License for more details.
00013  *
00014  * You should have received a copy of the GNU Library General Public
00015  * License along with this library; if not, write to the
00016  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
00017  * Boston, MA 02111-1307, USA.
00018  */
00019 
00020 #include "stdafx.h"
00021 
00022 #include <iostream>
00023 
00024 #include "gul_types.h"
00025 #include "gul_error.h"
00026 #include "guar_exact.h"
00027 #include "gugr_basics.h"
00028 
00029 namespace gugr {
00030 
00031 using gul::AllocError;
00032 
00033 
00034 GULAPI line::line( const point2<rational>& a, const point2<rational>& b )
00035 {
00036   rational dx,dy,adx,ady; // the constructor initializes these variables to 0
00037                           // (so adx becomes 0, when dx=0 and dy!=0 below)
00038   dx = b.x - a.x;
00039   dy = b.y - a.y;
00040   if( !dx.is_zero() ) ady = dy/dx; else ady = rational(ULong(1));
00041   if( !dy.is_zero() ) { 
00042     if( !dx.is_zero() ) 
00043       adx = rational(ULong(1))/ady;
00044   }
00045   else 
00046     adx = rational(ULong(1));
00047   m_rep = new line_rep( a, adx, ady );
00048 }
00049 
00050 bool intersect_nonvert( const line& a, const line& b, point2<rational>& I )
00051 {
00052   const rational& x0 = a.v().x;
00053   const rational& y0 = a.v().y;
00054   const rational& dy = a.dy();
00055   const rational& X0 = b.v().x;
00056   const rational& Y0 = b.v().y;
00057   const rational& DY = b.dy();
00058   rational d;
00059   int sign;
00060 
00061   if( a.IsVertical() || b.IsVertical() ) return false;
00062 
00063   d = DY - dy;
00064   sign = d.test();
00065   if( sign == 0 ) return false;
00066   
00067   I.x = (y0 - Y0 + X0*DY - x0*dy) / d;
00068   I.y = y0 + (I.x-x0)*dy;
00069 
00070   return true;
00071 }
00072 
00073 bool intersect_vert( const rational& x0, const line& b, point2<rational>& I )
00074 {
00075   const rational& X0 = b.v().x;
00076   const rational& Y0 = b.v().y;
00077   const rational& DY = b.dy();
00078 
00079   if( b.IsVertical() ) return false;
00080   if( b.IsHorizontal() )
00081     I.y = Y0;
00082   else
00083     I.y = Y0 + (x0 - X0)*DY;
00084   I.x = x0;
00085 
00086   return true;
00087 }
00088 
00089 bool intersect_horiz( const rational& y0, const line& b, point2<rational>& I )
00090 {
00091   const rational& X0 = b.v().x;
00092   const rational& Y0 = b.v().y;
00093   const rational& DX = b.dx();
00094 
00095   if( b.IsHorizontal() ) return false;
00096   if( b.IsVertical() )
00097     I.x = X0;
00098   else
00099     I.x = X0 + (y0 - Y0)*DX;
00100   I.y = y0;
00101 
00102   return true;
00103 }
00104 
00105 bool intersect( const line& a, const line& b, point2<rational>& I )
00106 {
00107   if( a.IsVertical() )
00108   {
00109     if( b.IsVertical() ) return false;
00110     return intersect_vert( a.v().x, b, I ); 
00111   }
00112   else if( a.IsHorizontal() )
00113   {
00114     if( b.IsHorizontal() ) return false;
00115     return intersect_horiz( a.v().y, b, I ); 
00116   }
00117 
00118   if( b.IsVertical() )
00119     return intersect_vert( b.v().x, a, I ); 
00120   else if( b.IsHorizontal() )
00121     return intersect_horiz( b.v().y, a, I ); 
00122   
00123   return intersect_nonvert( a, b, I );
00124 }
00125 
00126 
00127 /* ------------------------------------------------------------------------
00128   Appends a edge to the left chain of a monotone polygon. If this closes
00129   the polygon the function returns true
00130 ------------------------------------------------------------------------ */
00131 bool monpoly_info::AppendLeft( const graph_vertex& v1, const graph_vertex& v2 )
00132 {
00133   graph_vertex *vl1,*vl2;
00134   graph_edge *e;
00135   
00136   if( Vl.IsEmpty() )
00137   {
00138     vl2 = new graph_vertex(v2.v,NULL,v2.reserved[1]);
00139     Vl.AppendLeft(vl2);
00140   }
00141   else
00142     vl2 = Vl.head;
00143     
00144   e = new graph_edge();
00145   E.Append(e);
00146     
00147   vl1 = new graph_vertex(v1.v,e,v1.reserved[1]);
00148   Vl.AppendLeft(vl1);
00149  
00150   e->v[0] = vl1;
00151   e->v[1] = vl2;    
00152                                         /* if v is endpoint of right chain: */
00153   if((!Vr.IsEmpty())&&(Vr.head->v == v1.v)) 
00154     return(true);
00155 
00156   return(false);
00157 }
00158 
00159 /* ------------------------------------------------------------------------
00160   Appends a vertex to the right chain of a monotone polygon. If this closes
00161   the polygon the function returns true
00162 -------------------------------------------------------------------------- */
00163 bool monpoly_info::AppendRight( const graph_vertex& v1, const graph_vertex& v2 )
00164 {
00165   graph_vertex *vr1,*vr2;
00166   graph_edge *e;
00167   
00168   if( Vr.IsEmpty() )
00169   {
00170     vr2 = new graph_vertex(v2.v,NULL,v2.reserved[1]);
00171     Vr.AppendLeft(vr2);
00172   }
00173   else
00174     vr2 = Vr.head;
00175     
00176   e = new graph_edge();
00177   E.Append(e);
00178     
00179   vr1 = new graph_vertex(v1.v,e,v1.reserved[1]);
00180   Vr.AppendLeft(vr1);
00181  
00182   e->v[0] = vr1;
00183   e->v[1] = vr2;    
00184   
00185   if((!Vl.IsEmpty())&&(Vl.head->v == v1.v)) /*if v is endpoint of left chain*/
00186     return(true);
00187 
00188   return(false);
00189 }
00190 
00191 /*-------------------------------------------------------------------------
00192   Converts simple polygon to a DCEL (doubly connected edge list, winged edge
00193   list)
00194 --------------------------------------------------------------------------*/
00195 template< class T >
00196 void PolygonToGraph( 
00197       const int n, gul::Ptr< gul::point2<T> > P,
00198       const T orgx, const T orgy,
00199       const T scalex, const T scaley,
00200       int fleft, int fright,
00201       graph_edge_list *E, graph_vertex_list *V )
00202 {
00203   graph_vertex  *v, *vend, vdummy;
00204   graph_edge    *e = 0, *eend, edummy;
00205   int i;
00206   T x,y;
00207   graph_vertex_list hv;
00208   graph_edge_list he;
00209   unsigned long *cbuf = (unsigned long *)alloca(gul::rtr<T>::mantissa_length());
00210   
00211   he.head = eend = &edummy;
00212   hv.head = vend = &vdummy;
00213    
00214   for( i=0; i<n; i++ )
00215   {
00216     x = (P[i].x - orgx)/scalex;
00217     y = (P[i].y - orgy)/scaley;
00218 
00219     e = new graph_edge();
00220     he.Append(e);
00221 
00222     v = new graph_vertex( point2<rational>(
00223          rational(coord2int(x,cbuf),cbuf),
00224          rational(coord2int(y,cbuf),cbuf) ), e );
00225     hv.Append(v);
00226     
00227     e->f[0].set_i(fleft);
00228     e->f[1].set_i(fright);
00229 
00230     e->v[0] = v;
00231     e->e[0] = e->next;
00232 
00233     e->next->v[1] = v;
00234     e->next->e[1] = e;
00235   }
00236   e->v[1] = vend->prev;
00237   e->e[1] = eend->prev;
00238   eend->prev->e[0] = e;
00239 
00240   vend->prev->next = V->head;
00241   eend->prev->next = E->head;
00242   if( V->head != NULL ) V->head->prev = vend->prev;
00243   if( E->head != NULL ) E->head->prev = eend->prev;
00244 
00245   E->head = he.head; E->nElems += n;
00246   V->head = hv.head; V->nElems += n;
00247   he.ReleaseElems();
00248   hv.ReleaseElems();
00249 }
00250 // template instantiation
00251 template void PolygonToGraph( 
00252       const int n, gul::Ptr< gul::point2<float> > P,
00253       const float orgx, const float orgy,
00254       const float scalex, const float scaley,
00255       int fleft, int fright,
00256       graph_edge_list *E, graph_vertex_list *V );
00257 template void PolygonToGraph( 
00258       const int n, gul::Ptr< gul::point2<double> > P,
00259       const double orgx, const double orgy,
00260       const double scalex, const double scaley,
00261       int fleft, int fright,
00262       graph_edge_list *E, graph_vertex_list *V );
00263 
00264 
00265 
00266 template< class T >
00267 void PolygonToGraphNV( 
00268       const int n, gul::Ptr< gul::point2<T> > P,
00269       const T orgx, const T orgy,
00270       const T scalex, const T scaley,
00271       int fleft, int fright,
00272       graph_edge_list *E, graph_vertex_list *V )
00273 {
00274   graph_vertex  *v, *vend, vdummy;
00275   graph_edge    *e = 0, *eend, edummy;
00276   int i;
00277   T x,y;
00278   graph_vertex_list hv;
00279   graph_edge_list he;
00280   unsigned long *cbuf = (unsigned long *)alloca(gul::rtr<T>::mantissa_length());
00281   
00282   he.head = eend = &edummy;
00283   hv.head = vend = &vdummy;
00284    
00285   for( i=0; i<n; i++ )
00286   {
00287     x = (P[i].x - orgx)/scalex;
00288     y = (P[i].y - orgy)/scaley;
00289 
00290     e = new graph_edge();
00291     he.Append(e);
00292 
00293     v = new graph_vertex( point2<rational>(
00294          rational(coord2int(x,cbuf),cbuf),
00295          rational(coord2int(y,cbuf),cbuf) ), e );
00296     hv.Append(v);
00297     v->v.m_rep->reserved.i = i;  // enumber vertices
00298 
00299     e->f[0].set_i(fleft);
00300     e->f[1].set_i(fright);
00301 
00302     e->v[0] = v;
00303     e->e[0] = e->next;
00304 
00305     e->next->v[1] = v;
00306     e->next->e[1] = e;
00307   }
00308   e->v[1] = vend->prev;
00309   e->e[1] = eend->prev;
00310   eend->prev->e[0] = e;
00311 
00312   vend->prev->next = V->head;
00313   eend->prev->next = E->head;
00314   if( V->head != NULL ) V->head->prev = vend->prev;
00315   if( E->head != NULL ) E->head->prev = eend->prev;
00316 
00317   E->head = he.head; E->nElems += n;
00318   V->head = hv.head; V->nElems += n;
00319   he.ReleaseElems();
00320   hv.ReleaseElems();
00321 }
00322 // template instantiation
00323 template void PolygonToGraphNV( 
00324       const int n, gul::Ptr< gul::point2<float> > P,
00325       const float orgx, const float orgy,
00326       const float scalex, const float scaley,
00327       int fleft, int fright,
00328       graph_edge_list *E, graph_vertex_list *V );
00329 template void PolygonToGraphNV( 
00330       const int n, gul::Ptr< gul::point2<double> > P,
00331       const double orgx, const double orgy,
00332       const double scalex, const double scaley,
00333       int fleft, int fright,
00334       graph_edge_list *E, graph_vertex_list *V );
00335 
00336 
00337 /* ----------------------------------------------------------------
00338    gets the edges incident to vertex 'v'
00339 ------------------------------------------------------------------ */
00340 
00341 int IncidentEdges( const graph_vertex *v, graph_edge **A )
00342 {
00343   graph_edge *a0,*a;
00344   int i;
00345   
00346   a = A[0] = a0 = v->e;
00347   if( !a ) return 0;
00348 
00349   if( a->v[0] == v )
00350     a = a->e[0];
00351   else
00352     a = a->e[1];
00353 
00354   i = 1;
00355   while( a != a0 )
00356   {
00357     A[i] = a;
00358     i++;
00359 
00360     if( a->v[0] == v )
00361       a = a->e[0];
00362     else
00363       a = a->e[1];
00364   }
00365 
00366   return(i);
00367 }
00368 
00369 /* ----------------------------------------------------------------
00370    gets the edges incident to edge 'e' at vertex 'v'
00371 ------------------------------------------------------------------ */
00372 
00373 int EdgeCycle( graph_edge *e, graph_vertex *v, graph_edge **A )
00374 {
00375   graph_edge *a0,*a;
00376   int i;
00377   
00378   A[0] = a0 = a = e;
00379 
00380   if( a->v[0] == v )
00381     a = a->e[0];
00382   else
00383     a = a->e[1];
00384 
00385   i = 1;
00386   while( a != a0 )
00387   {
00388     A[i] = a;
00389     i++;
00390 
00391     if( a->v[0] == v )
00392       a = a->e[0];
00393     else
00394       a = a->e[1];
00395   }
00396 
00397   return(i);
00398 }
00399 
00400 /* --------------------------------------------------------------------
00401   Orient edge, so that y2 > y1 or x2 < x1 if y2 = y1
00402 --------------------------------------------------------------------- */
00403 
00404 void OrientEdge( graph_edge *e )
00405 {
00406   graph_edge *tmpe;
00407   graph_vertex *tmpv;
00408   gul::ptr_int_union tmpf;
00409   int sign;
00410   
00411   if( e != NULL )
00412   {
00413     sign = compare( e->v[0]->v.v().y, e->v[1]->v.v().y );
00414     if( !sign ) sign = -compare( e->v[0]->v.v().x, e->v[1]->v.v().x );
00415 
00416     if( sign > 0 )  /* v1 > v2 */
00417     {
00418       tmpv = e->v[0];
00419       e->v[0] = e->v[1];
00420       e->v[1] = tmpv;
00421 
00422       tmpe = e->e[0];
00423       e->e[0] = e->e[1];
00424       e->e[1] = tmpe;
00425 
00426       tmpf = e->f[0];
00427       e->f[0] = e->f[1];
00428       e->f[1] = tmpf;
00429     }
00430   }
00431 }
00432 
00433 /* --------------------------------------------------------------------
00434   Orient edges, so that y2 > y1 or x2 < x1 if y2 = y1
00435 --------------------------------------------------------------------- */
00436 
00437 void OrientEdges( graph_edge *E )
00438 {
00439   graph_edge *e, *tmpe;
00440   graph_vertex *tmpv;
00441   gul::ptr_int_union tmpf;
00442   int sign;
00443   
00444   e = E;
00445   while( e != NULL )
00446   {
00447     sign = compare( e->v[0]->v.v().y, e->v[1]->v.v().y );
00448     if( !sign ) sign = -compare( e->v[0]->v.v().x, e->v[1]->v.v().x );
00449 
00450     if( sign > 0 )  /* v1 > v2 */
00451     {
00452       tmpv = e->v[0];
00453       e->v[0] = e->v[1];
00454       e->v[1] = tmpv;
00455 
00456       tmpe = e->e[0];
00457       e->e[0] = e->e[1];
00458       e->e[1] = tmpe;
00459 
00460       tmpf = e->f[0];
00461       e->f[0] = e->f[1];
00462       e->f[1] = tmpf;
00463     }
00464 
00465     e = e->next;
00466   }
00467 }
00468 
00469 
00470 
00471 
00472 /* ------------------------------------------------------------------------
00473   Determine on which side of AB the point C is
00474 -------------------------------------------------------------------------- */
00475 
00476 int DetermineOrientation( const point2<rational>& A, const point2<rational>& B, 
00477                           const point2<rational>& C )
00478 {
00479   const Interval& ix1 = A.x.get_bounds();
00480   const Interval& iy1 = A.y.get_bounds();
00481   const Interval& ix2 = B.x.get_bounds();
00482   const Interval& iy2 = B.y.get_bounds();
00483   const Interval& ix = C.x.get_bounds();
00484   const Interval& iy = C.y.get_bounds();
00485   Interval i;
00486 
00487   i = (ix2-ix1)*(iy-iy1) - (ix-ix1)*(iy2-iy1);
00488   if( i.m_low > 0.0 ) return(+1);
00489   else if( i.m_high < 0.0 ) return(-1);
00490 
00491   const rational& x1 = A.x;
00492   const rational& y1 = A.y;
00493   const rational& x2 = B.x;
00494   const rational& y2 = B.y;
00495   const rational& x = C.x;
00496   const rational& y = C.y;
00497   rational det;
00498 
00499   det = (x2-x1)*(y-y1) - (x-x1)*(y2-y1);
00500   return( det.test() );
00501 }
00502 
00503 int DetermineOrientationPV( const point2<rational>& A, const point2<rational>& AB, 
00504                             const point2<rational>& C )
00505 {
00506   const Interval& ix1 = A.x.get_bounds();
00507   const Interval& iy1 = A.y.get_bounds();
00508   const Interval& ix1x2 = AB.x.get_bounds();
00509   const Interval& iy1y2 = AB.y.get_bounds();
00510   const Interval& ix = C.x.get_bounds();
00511   const Interval& iy = C.y.get_bounds();
00512   Interval i;
00513 
00514   i = ix1x2*(iy-iy1) - (ix-ix1)*iy1y2;
00515   if( i.m_low > 0.0 ) return(+1);
00516   else if( i.m_high < 0.0 ) return(-1);
00517 
00518   const rational& x1 = A.x;
00519   const rational& y1 = A.y;
00520   const rational& x1x2 = AB.x;
00521   const rational& y1y2 = AB.y;
00522   const rational& x = C.x;
00523   const rational& y = C.y;
00524   rational det;
00525 
00526   det = x1x2*(y-y1) - (x-x1)*y1y2;
00527   return( det.test() );
00528 }
00529 
00530 /* --------------------------------------------------------------------
00531   Compares angle(AB,AC) to angle(AB,AD)
00532 
00533 
00534   Remarks:
00535 
00536   code_x = 3 if X is to the left of AB
00537            2 if X is on AB
00538            1 if X is to the right of AB
00539            0 if it was not necessary to calculate the cross-product
00540   match = 1 if C lies on AB
00541         = 2 if D lies on AB
00542         = 0 if both C and D are not on AB             
00543 ---------------------------------------------------------------------- */
00544 
00545 int CompareAngles(
00546   const point2<rational>& A, const point2<rational>& B,
00547   const point2<rational>& C, const point2<rational>& D,
00548   int *match, int *code_c, int *code_d )
00549 {
00550   const Interval& iAx = A.x.get_bounds();
00551   const Interval& iAy = A.y.get_bounds();
00552   const Interval& iBx = B.x.get_bounds();
00553   const Interval& iBy = B.y.get_bounds();
00554   const Interval& iCx = C.x.get_bounds();
00555   const Interval& iCy = C.y.get_bounds();
00556   const Interval& iDx = D.x.get_bounds();
00557   const Interval& iDy = D.y.get_bounds();
00558   Interval iABx,iABy, iACx,iACy, iADx,iADy;
00559   Interval iABcAC,iABdAC; 
00560   Interval iABcAD,iABdAD;
00561   bool flag1 = false, flag2 = false;
00562   int sideC,sideD,sign;
00563   const rational& rAx = A.x;
00564   const rational& rAy = A.y;
00565   const rational& rBx = B.x;
00566   const rational& rBy = B.y;
00567   const rational& rCx = C.x;
00568   const rational& rCy = C.y;
00569   const rational& rDx = D.x;
00570   const rational& rDy = D.y;
00571   rational rABx,rABy, rACx,rACy, rADx,rADy;
00572   rational rABcAC, rABdAC;   
00573   rational rABcAD, rABdAD;
00574   Interval iLambda;
00575   rational rLambda;   
00576 
00577   *match = 0;
00578   *code_c = 0;
00579   *code_d = 0;
00580   
00581   iABx = iBx - iAx;
00582   iABy = iBy - iAy;
00583   iACx = iCx - iAx;
00584   iACy = iCy - iAy;
00585   
00586   iABcAC = iABx*iACy - iABy*iACx;
00587   
00588   if( iABcAC.m_low > 0.0 ) sideC = +1;
00589   else if( iABcAC.m_high < 0.0 ) sideC = -1;
00590   else
00591   {
00592     flag1 = true;
00593     rABx = rBx - rAx;
00594     rABy = rBy - rAy;
00595     rACx = rCx - rAx;
00596     rACy = rCy - rAy;
00597 
00598     rABcAC = rABx*rACy - rABy*rACx;
00599     sideC = rABcAC.test();
00600   }
00601 
00602   *code_c = sideC + 2;
00603   
00604   if( sideC == 0 )    /* sideC == 0 => check orientation with dot product */
00605   {
00606     iABdAC = iABx*iACx + iABy*iACy;
00607     if( iABdAC.m_low > 0.0 )        /* 0=angle(AB,AC) < angle(AB,AD) */
00608     { *match = 1; return(-1); }   
00609     else if( iABdAC.m_high >= 0.0 )
00610     {    
00611       rABdAC = rABx*rACx + rABy*rACy;
00612       if( rABdAC.test() > 0.0 )     /* 0=angle(AB,AC) < angle(AB,AD) */
00613       { *match = 1; return(-1); } 
00614     }
00615     /* now angle(AB,AC) == 180 if sideC == 0 */
00616   }
00617  
00618   iADx = iDx - iAx;
00619   iADy = iDy - iAy;
00620 
00621   iABcAD = iABx*iADy - iABy*iADx;
00622   
00623   if( iABcAD.m_low > 0.0 ) sideD = +1;
00624   else if( iABcAD.m_high < 0.0 ) sideD = -1;
00625   else
00626   {
00627     flag2 = true;
00628     rADx = rDx - rAx;
00629     rADy = rDy - rAy;
00630     if( !flag1 )
00631     {
00632       rABx = rBx - rAx;
00633       rABy = rBy - rAy;
00634     }
00635     rABcAD = rABx*rADy - rABy*rADx;
00636     sideD = rABcAD.test();
00637   }
00638 
00639   *code_d = sideD + 2;
00640   
00641   if( sideD == 0 ) /* sideD == 0 => check orientation with dot product */
00642   {
00643     if( sideC == 0 ) { *match = 2; return(+1); }
00644 
00645     iABdAD = iABx*iADx + iABy*iADy;
00646     if( iABdAD.m_low > 0.0 )      /* 0=angle(AB,AD) < angle(AB,AC) */
00647     { *match = 2; return(+1); }   
00648     else if( iABdAD.m_high >= 0.0 ) /* exact calculation if necessary */
00649     {    
00650       rABdAD = rABx*rADx + rABy*rADy;
00651       if( rABdAD.test() > 0.0 )     /* 0=angle(AB,AD) < angle(AB,AC) */
00652       { *match = 2; return(+1); } 
00653     }
00654     /* now angle(AB,AD) == 180 if sideD == 0 */
00655 
00656     if( sideC > 0 ) return(-1);
00657     return(+1);
00658   }
00659 
00660   if( sideC == 0 )
00661   {
00662     if( sideD > 0 ) return(+1);
00663     return(-1);
00664   }
00665   
00666   if( sideC > sideD ) return(-1);
00667   else if( sideC < sideD ) return(+1);
00668 
00669   /* intersect Line(C,AB) and Line(A,AD) */
00670   
00671   iLambda = (iADx*iACy - iACx*iADy)/(iABx*iADy - iABy*iADx);
00672   if( iLambda.m_low > 0.0 )
00673   {
00674      if( sideC > 0 ) return(+1);  /* angle(AB,AC) > angle(AB,AD) */
00675      return(-1);
00676   }
00677   if( iLambda.m_high < 0.0 )
00678   {
00679     if( sideC > 0 ) return(-1);  /* angle(AB,AC) < angle(AB,AD) */
00680     return(+1);
00681   }
00682   
00683   if( !flag1 )
00684   {
00685     rACx = rCx - rAx;
00686     rACy = rCy - rAy;
00687   }
00688   if( !flag2 )
00689   {
00690     if( !flag1 )
00691     {
00692       rABx = rBx - rAx;
00693       rABy = rBy - rAy;
00694     }
00695     rADx = rDx - rAx;
00696     rADy = rDy - rAy;
00697   }
00698   
00699   rLambda = (rADx*rACy - rACx*rADy)/(rABx*rADy - rABy*rADx);
00700   sign = rLambda.test();
00701 
00702   if( sign > 0 )
00703   {
00704      if( sideC > 0 ) return(+1);  /* angle(AB,AC) > angle(AB,AD) */
00705      return(-1);
00706   }
00707   if( sideC > 0 ) return(-1);  /* angle(AB,AC) < angle(AB,AD) */
00708   return(+1);
00709 }
00710 
00711 /* --------------------------------------------------------------------
00712   Find the insertion positon for a oriented line segment in an edge 
00713   cycle of a DCEL 
00714 
00715   Remarks:
00716 
00717   code_edge = 1 if angle(ab,edge) > 180
00718             = 2 if angle(ab,edge) = 0 or 180
00719             = 3 if angle(ab,edge) < 180
00720   match = 1 if ab and left are colinear 
00721         = 2 if ab and right are colinear
00722         = 0 if none of both edges lies on ab         
00723 ---------------------------------------------------------------------*/
00724 
00725 int EdgeInsertPosition(
00726          const vertex& a, const vertex& b,
00727          int nE, graph_edge **E,
00728          int *eleft, int *eright,
00729          int *code_left, int *code_right )
00730 {
00731   vertex v1,v2,v3; 
00732   int left,right,mid,sign,s1,s2,s3,tmp;
00733   int match = 0;
00734   
00735   left = 0;
00736   right = nE-1;
00737 
00738   if( right > 0 )
00739   {
00740     if( E[left]->v[0]->v == a ) v1 =  E[left]->v[1]->v;
00741     else v1 = E[left]->v[0]->v;
00742     if( E[right]->v[0]->v == a ) v2 =  E[right]->v[1]->v;
00743     else v2 = E[right]->v[0]->v;
00744 
00745     sign = CompareAngles( a.v(), b.v(), v1.v(), v2.v(),
00746                           &match, &s1, &s2 );
00747     if( sign > 0 )
00748     {
00749       if( match )
00750       {
00751         if( right-1 != left ) s1 = 0;
00752         left = right-1;
00753       }
00754       else
00755       {
00756         mid = (left+right)/2;
00757 
00758         while( mid != left ) 
00759         {
00760           if( E[mid]->v[0]->v == a ) v3 =  E[mid]->v[1]->v;
00761           else v3 = E[mid]->v[0]->v;
00762 
00763           sign = CompareAngles( a.v(), b.v(), v1.v(), v3.v(), 
00764                                 &match, &s1, &s3 );
00765           if( match )
00766           {
00767             right = mid;
00768             s2 = s3;
00769             if( mid-1 != left ) s1 = 0;
00770             left = mid-1;
00771             break;
00772           } 
00773 
00774           if( sign > 0 )
00775           {
00776             right = mid;
00777             s2 = s3;
00778           }
00779           else if( sign < 0 )
00780           {
00781             left = mid;
00782             v1 = v3;
00783             s1 = s3;
00784           }
00785           else 
00786             throw InternalError();   /* this cannot be */
00787 
00788           mid = (left+right)/2;
00789         }
00790       }
00791     }
00792     else
00793     {
00794       gul::Assert< gul::InternalError >( ndebug ||
00795                         ((sign != 0) && (left == 0) && (right == nE-1)) );
00796  
00797       tmp = left; left = right; right = tmp;
00798       tmp = s1; s1 = s2; s2 = tmp;
00799     }
00800   }
00801   else
00802   {
00803     right = left;
00804     s1 = s2 = 0;
00805   }
00806 
00807   if( s1 == 0 )
00808   {
00809     if( E[left]->v[0]->v == a ) v1 =  E[left]->v[1]->v;
00810     else v1 = E[left]->v[0]->v;
00811 
00812     s1 = DetermineOrientation( a.v(), b.v(), v1.v() ) + 2;
00813   }
00814   if( s2 == 0 )
00815   {
00816     if( right != left )
00817     {
00818       if( E[right]->v[0]->v == a ) v2 =  E[right]->v[1]->v;
00819       else v2 = E[right]->v[0]->v;
00820 
00821       s2 = DetermineOrientation( a.v(), b.v(), v2.v() ) + 2;
00822     }
00823     else
00824       s2 = s1;
00825   }
00826 
00827   *eleft = left;
00828   *eright = right; 
00829   *code_left = s1;
00830   *code_right = s2;
00831   return(match);
00832 }
00833 
00834 /* ----------------------------------------------------------------------
00835   Connects a edge at its ends into a graph
00836 -----------------------------------------------------------------------*/
00837 
00838 void ConnectEdge( graph_edge *e, int maxE )
00839 {
00840   graph_edge *left, *right, **E;
00841   graph_vertex *v0 = e->v[0], *v1 = e->v[1];
00842   int ileft, iright, codel, coder, nE;
00843  
00844   E = (graph_edge **)alloca( sizeof(graph_edge *) * maxE );
00845   if( E == NULL ) throw AllocError();
00846  
00847   nE = IncidentEdges( v0, E );
00848   if( nE > 0 )
00849   {
00850     EdgeInsertPosition( v0->v, v1->v, nE, E, &ileft, &iright,
00851                         &codel, &coder );
00852     left = E[ileft];  right = E[iright];
00853 
00854     if( left->v[0] == v0 )
00855       left->e[0] = e;
00856     else
00857       left->e[1] = e;
00858 
00859     e->e[0] = right; 
00860   }
00861   else
00862   {
00863     v0->e = e;    
00864     e->e[0] = e;
00865   }
00866 
00867   nE = IncidentEdges( v1, E );
00868   if( nE > 0 )
00869   {
00870     EdgeInsertPosition( v1->v, v0->v, nE, E, &ileft, &iright,
00871                         &coder, &codel );
00872     left = E[ileft];  right = E[iright];
00873 
00874     if( left->v[0] == v1 )
00875       left->e[0] = e;
00876     else
00877       left->e[1] = e;
00878 
00879     e->e[1] = right; 
00880   }
00881   else
00882   {
00883     v1->e = e;    
00884     e->e[1] = e;
00885   }
00886 }
00887 
00888 /* ----------------------------------------------------------------------
00889   Disconnect the ends of an edge from a graph
00890 -----------------------------------------------------------------------*/
00891 
00892 int DisconnectEdge( graph_edge *e, int maxE )
00893 {
00894   graph_edge *left, *right, **E;
00895   graph_vertex *v;
00896   int n,isolated;
00897  
00898   E = (graph_edge **)alloca( sizeof(graph_edge *) * maxE );
00899   if( E == NULL ) throw AllocError();
00900 
00901   isolated = 0;
00902   v = e->v[0];
00903  
00904   n = EdgeCycle( e, v, E );
00905   if( n > 1 )
00906   {
00907     right = E[1]; 
00908     left = E[n-1];
00909 
00910     if( left->v[0] == v )
00911       left->e[0] = right;
00912     else 
00913       left->e[1] = right;
00914 
00915     if( v->e == e ) v->e = left;
00916   }
00917   else
00918     isolated |= 1;
00919 
00920   v = e->v[1];
00921  
00922   n = EdgeCycle( e, v, E );
00923   if( n > 1 )
00924   {
00925     right = E[1]; 
00926     left = E[n-1];
00927 
00928     if( left->v[0] == v ) 
00929       left->e[0] = right;
00930     else 
00931       left->e[1] = right;
00932 
00933     if( v->e == e ) v->e = left;
00934   }
00935   else
00936     isolated |= 2;
00937 
00938   return isolated;
00939 }
00940 
00941 }

Generated on Mon Jan 21 04:17:28 2002 for GUL 0.6 - Geometry Utility Library by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001