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

gugr_split.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 "guma_sorting.h"
00027 #include "guar_exact.h"
00028 #include "gugr_basics.h"
00029 #include "gugr_split.h"
00030 #include "gul_io.h"
00031 #include "gul_std.h"
00032 #include "gugr_planesweep.h"
00033 #include "gugr_io.h"
00034 
00035 namespace gugr {
00036 
00037 // using namespace guar;
00038 using gul::AllocError;
00039 using gul::Ptr;
00040 
00041 /* ------------------------------------------------------------------------
00042   Sort "Cut records" according to their value
00043 ------------------------------------------------------------------------- */
00044 
00045 int compare_cut_records( void *p1, void *p2, void *data )
00046 {
00047   rational &r1 = ((cut_record *)p1)->val, &r2 = ((cut_record *)p2)->val;
00048   int sign;
00049   
00050   sign = compare( r1, r2 );
00051   return(sign);
00052 }
00053 
00054 void sort_cut_records( int ns, Ptr< cut_info > S )
00055 {
00056   cut_record *rec, **tmpR;
00057   int i,j,max_nR;
00058   
00059   max_nR = 0;
00060   for( i = 1; i <= ns; i++ )
00061     if( S[i].R.nElems > max_nR ) max_nR = S[i].R.nElems;
00062 
00063   if( max_nR > 0 )
00064   {
00065     tmpR = (cut_record **)alloca( sizeof(cut_record *)*max_nR );
00066     if( !tmpR ) throw AllocError();
00067   }
00068   else tmpR = 0;  // make compiler happy :)
00069     
00070   for( i = 1; i <= ns; i++ )
00071   {
00072     if( S[i].R.nElems > 1 )
00073     {
00074       rec = S[i].R.head;
00075       for( j = 0; j < S[i].R.nElems; j++ )
00076       {
00077         tmpR[j] = rec;
00078         rec = rec->next;
00079       } 
00080       guma::PtrHeapSort( S[i].R.nElems,(void **)tmpR, compare_cut_records,NULL);
00081 
00082       S[i].R.head = tmpR[0];
00083       tmpR[0]->prev = NULL;
00084       for( j = 1; j < S[i].R.nElems; j++ )
00085       {
00086         tmpR[j-1]->next = tmpR[j];
00087         tmpR[j]->prev = tmpR[j-1];
00088       }
00089       tmpR[S[i].R.nElems-1]->next = NULL;
00090     }
00091   }
00092 }
00093 
00094 
00095 void IntersectWithAscendant( 
00096                      const point2<rational>& A, const point2<rational>& B,
00097                      const graph_edge_list& E, const graph_vertex_list& V,
00098                      line& L, cut_record_list& R )
00099 {
00100   rational y;
00101   graph_edge *e;
00102   graph_vertex *v;
00103   int i,k,sign;
00104   cut_record *rec, **tmpR;
00105   point2<rational> AB,X;
00106 
00107   AB.x = B.x - A.x;
00108   AB.y = B.y - A.y;
00109   L = line(A,B);
00110 
00111   v = V.head;
00112   while( v != NULL )
00113   {
00114     sign = DetermineOrientationPV( A, AB, v->v.v() );
00115     v->reserved[1].set_i(sign);
00116     if( sign == 0 )      // new record, when on diagonal
00117     {
00118       rec = new cut_record();
00119       R.Append( rec );      
00120       rec->v = v;
00121       rec->e = NULL;
00122       rec->val = v->v.v().y;
00123     }
00124     v = v->next;
00125   }
00126 
00127   e = E.head;
00128   while( e != NULL )  
00129   {
00130     i = e->v[0]->reserved[1].i;
00131     k = e->v[1]->reserved[1].i;
00132 
00133     if( k*i < 0 ) // edge intersects diagonal 
00134     {
00135       if( e->l.IsNULL() )
00136         e->CalcLine();
00137 
00138       intersect( e->l, L, X );
00139 
00140       rec = new cut_record();
00141       R.Append( rec );
00142       rec->v = new graph_vertex( X, NULL );
00143       rec->v->reserved[1].set_i(0);  // on diagonal
00144       rec->e = e;          
00145       rec->val = X.y;
00146     }
00147     e = e->next;
00148   }
00149                                          /* --- sort cut_records (on y) ---- */
00150   tmpR = (cut_record **)alloca( sizeof(cut_record *)*R.nElems );
00151   if( !tmpR ) throw AllocError();
00152 
00153   rec = R.head;
00154   for( i = 0; i < R.nElems; i++ )
00155   {
00156     tmpR[i] = rec;
00157     rec = rec->next;
00158   } 
00159   guma::PtrHeapSort( R.nElems,(void **)tmpR, compare_cut_records,NULL);
00160 
00161   R.head = tmpR[0];
00162   tmpR[0]->prev = NULL;
00163   for( i = 1; i < R.nElems; i++ )
00164   {
00165     tmpR[i-1]->next = tmpR[i];
00166     tmpR[i]->prev = tmpR[i-1];
00167   }
00168   tmpR[R.nElems-1]->next = NULL;
00169 }
00170 
00171 
00172 /*--------------------------------------------------------------------------
00173   inserts a diagonal into the graph of a quad
00174   uses the reserved[1] field of vertices 
00175 ---------------------------------------------------------------------------*/
00176 void InsertDiagonal( 
00177               graph_edge_list *E, graph_vertex_list *V,
00178               graph_vertex *P11, graph_vertex *P22 )
00179 
00180 {
00181   cut_record_list Recs;
00182   cut_record *rec;
00183   int nhE,L,R,match,cR,cL,ivh,ivl,maxE;
00184   graph_edge *e,*eh,*el,**hE;
00185   graph_vertex *lv;
00186   vertex *lowerleft = &P11->v;
00187   vertex *upperright = &P22->v;
00188   line lin;
00189 
00190   IntersectWithAscendant( P11->v.v(), P22->v.v(), *E, *V, lin, Recs );
00191 
00192   maxE = 3*E->nElems + V->nElems;
00193   hE = (graph_edge **)alloca( sizeof(graph_edge *)*maxE );
00194   if( hE == NULL ) throw AllocError();
00195 
00196   rec = Recs.head;
00197   lv = rec->v;
00198   gul::Assert<InternalError>( ndebug || (lv == P11) );
00199   rec = rec->next;
00200 
00201   while( rec != NULL )
00202   {
00203     if( rec->e == NULL )    /* vertex on diagonal */
00204     {
00205       nhE = IncidentEdges( rec->v, hE );
00206       match = EdgeInsertPosition( rec->v->v, *lowerleft, nhE, hE, 
00207                                   &L, &R, &cL, &cR );
00208       if( cR != 2 ) // have to create a new edge
00209       {
00210         e = new graph_edge();
00211         E->Append(e);
00212 
00213         e->v[0] = lv;
00214         e->v[1] = rec->v;
00215 
00216         if( hE[L]->v[0] == rec->v ) // connect upper end of new edge
00217         {
00218           e->e[1] = hE[L]->e[0];
00219           hE[L]->e[0] = e;
00220           e->f[0] = e->f[1] = hE[L]->f[0];
00221         }
00222         else
00223         {
00224           e->e[1] = hE[L]->e[1];
00225           hE[L]->e[1] = e;
00226           e->f[0] = e->f[1] = hE[L]->f[1];
00227         }
00228 
00229         nhE = IncidentEdges( lv, hE );
00230         match = EdgeInsertPosition( lv->v, *upperright, nhE, hE, 
00231                                     &L, &R, &cL, &cR ); 
00232         if( hE[L]->v[0] == lv ) // connect lower end of new edge
00233         {
00234           e->e[0] = hE[L]->e[0];
00235           hE[L]->e[0] = e;
00236         }
00237         else
00238         {
00239           e->e[0] = hE[L]->e[1];
00240           hE[L]->e[1] = e;
00241         }
00242       }
00243     }
00244     else           /* diagonal intersects an edge */
00245     {
00246       el = rec->e;
00247 
00248       eh = new graph_edge();  // new edge for half of edge above diagonal  
00249       E->Append(eh);
00250       rec->v->e = eh;         // give new vertex an incident edge
00251       V->Append(rec->v);      // append it to vertex list        
00252 
00253       if( el->v[0]->reserved[1].i > 0 ) // part with v0 above diagonal 
00254       {
00255         ivh = 0;
00256         ivl = 1;
00257       }
00258       else                                 // part with v1 above diagonal
00259       {
00260         ivh = 1;
00261         ivl = 0;
00262       }
00263       eh->v[ivl] = rec->v;
00264       eh->v[ivh] = el->v[ivh];
00265       eh->f[0] = el->f[0];
00266       eh->f[1] = el->f[1];
00267         
00268       nhE = EdgeCycle( el, el->v[ivh], hE );
00269       if( nhE > 1 )
00270       {
00271         eh->e[ivh] = el->e[ivh];
00272         if( hE[nhE-1]->e[0] == el ) 
00273           hE[nhE-1]->e[0] = eh; 
00274         else 
00275           hE[nhE-1]->e[1] = eh;
00276       }
00277       else
00278         eh->e[ivh] = eh;
00279         
00280       el->v[ivh]->e = eh;      // eh replaces original edge in upper edge cycle
00281       el->v[ivh] = rec->v;     // shorten intersecting edge to lower half
00282 
00283                                    /* ---- create new edge for diagonal ---- */
00284       e = new graph_edge();
00285       E->Append(e);
00286       e->v[0] = lv;
00287       e->v[1] = rec->v;
00288       e->f[0] = e->f[1] = el->f[ivl];
00289       
00290       nhE = IncidentEdges( lv, hE );
00291       match = EdgeInsertPosition( lv->v, *upperright, nhE, hE, 
00292                                     &L, &R, &cL, &cR ); 
00293       if( hE[L]->v[0] == lv ) // connect lower end of new edge
00294       {
00295         e->e[0] = hE[L]->e[0];
00296         hE[L]->e[0] = e;
00297       }
00298       else
00299       {
00300         e->e[0] = hE[L]->e[1];
00301         hE[L]->e[1] = e;
00302       }
00303       eh->e[ivl] = e;
00304       e->e[1] = el;
00305       el->e[ivh] = eh;
00306     }
00307     lv = rec->v;
00308     rec = rec->next;
00309   }
00310   gul::Assert<InternalError>( ndebug || (lv == P22) );
00311 
00312   Recs.DeleteElems();
00313 }
00314 
00315 
00316 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00317 
00318   Functions for splitting a graph into horizontal strips
00319 
00320 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00321 
00322 /* -------------------------------------------------------------------------
00323   Calculates the intersection of a graph with horizontal cutting lines 
00324 -------------------------------------------------------------------------- */
00325 
00326 
00327 /*------------------------------------------------------------
00328   found = true/false
00329   result = interval
00330 -------------------------------------------------------------*/
00331 
00332 int FindCutInfo( rational& x, int nA, const Ptr< cut_info >& A,
00333                  bool& found )
00334 {
00335   int mid,right,left,sign;
00336 
00337   if( nA == 0 ) {found = false; return -1;}
00338 
00339   sign = compare( A[0].val, x );
00340   if( sign >= 0 )
00341   {
00342     if( sign > 0 ) {found = false; return -1;}
00343     found = true; return 0;
00344   }
00345 
00346   sign = compare( A[nA-1].val, x );
00347   if( sign <= 0 )
00348   {
00349     if( sign < 0 ) {found = false; return nA-1;}
00350     found = true; return nA-1;
00351   }
00352 
00353   left = 0;
00354   right = nA;
00355   mid = (left+right)/2;
00356 
00357   while( (mid!=left) && (mid!=right) )
00358   {  
00359     sign = compare( A[mid].val, x );
00360 
00361     if( sign == 0  )
00362       break;
00363     else if( sign > 0 )
00364     {
00365       right = mid;
00366       mid = (left+right)>>1;
00367     }
00368     else
00369     {
00370       left = mid;
00371       mid = (left+right)>>1;
00372     }
00373   }
00374 
00375   found = (sign == 0);
00376   return mid;
00377 }
00378 
00379 
00380 
00381 void IntersectWithHorizontals( 
00382           graph_edge_list *E, graph_vertex_list *V,
00383           const int nS, Ptr< cut_info > S )
00384 {
00385   int ns = nS-1, i,j,k,nA;
00386   graph_vertex *v,*v_next;
00387   graph_edge *e, *e_next;
00388   rational y,x0,y0,x,dx;
00389   bool rz;
00390   cut_record *rec;
00391   Ptr< cut_info > A;
00392 
00393   v = V->head;
00394 
00395   nA = nS-1;            
00396   A.use_pointer( &S[1], nS-1 );  // used for binary search
00397 
00398   while( v != NULL )
00399   {
00400     v_next = v->next;
00401 
00402     y = v->v.v().y;
00403 
00404     i = FindCutInfo( y, nA, A, rz ) + 1;
00405  
00406     v->reserved[0].set_i(i);
00407     v->reserved[1].set_i(rz);
00408     
00409     V->Remove( v );
00410     S[i].V.Append( v );
00411 
00412     if( rz )
00413     { 
00414       rec = new cut_record();
00415       S[i].R.Append( rec );
00416            
00417       rec->v = v;
00418       rec->e = NULL;
00419       rec->val = v->v.v().x;
00420     }
00421 
00422     v = v_next;
00423   }
00424 
00425   e = E->head;
00426 
00427   while( e != NULL )  
00428   {
00429     e_next = e->next;
00430 
00431     i = e->v[0]->reserved[0].i;
00432     k = e->v[1]->reserved[0].i;
00433 
00434     E->Remove( e );
00435     S[k].E.Append( e );
00436 
00437     if( i !=  k )
00438     {
00439       if( e->l.IsNULL() )
00440       {
00441         e->CalcLine();
00442       }
00443 
00444       x0 = e->l.v().x;
00445       y0 = e->l.v().y;
00446     
00447       dx = e->l.dx();
00448 
00449       if( e->v[1]->reserved[1].i ) k--;
00450 
00451       for( j = i+1; j <= k; j++ )
00452       {
00453         if( !dx.is_zero() )
00454           x = x0 + (S[j].val-y0)*dx;
00455         else
00456           x = x0;
00457           
00458         rec = new cut_record();
00459         S[j].R.Append( rec );
00460               
00461         rec->v = NULL;
00462         rec->e = e;          
00463         rec->val = x;
00464       }
00465     }
00466 
00467     e = e_next;
00468   }
00469 
00470   sort_cut_records( ns, S );
00471 }
00472 
00473 /* -------------------------------------------------------------------------
00474   Divide graph into an independant upper and lower part
00475 -------------------------------------------------------------------------- */
00476  
00477 void DivideHorizontalStrips( 
00478               cut_info *Sa, cut_info *Sb, int maxE,
00479               const rational& minx, const rational maxx,
00480               const rational& far_left, const rational& far_right,
00481               graph_vertex **minSaV, graph_vertex **minSbV,
00482               graph_vertex **maxSaV, graph_vertex **maxSbV )
00483 {
00484   graph_edge **E;
00485   cut_record *rec;
00486   int nE,L1,L2,R1,R2,match1,match2,cL1,cR1,cL2,cR2;
00487   vertex farleft,farright;
00488   graph_edge *ah,*el,*bl,*bh,*er,*al, *e, *eu, *ew;
00489   graph_vertex *v, *lva, *lvb, *vu, *vl, *delv;
00490   int fa,fb,nfa,nfb;
00491   line hori;
00492 
00493   farleft = vertex( point2<rational>(far_left,Sa->val) );
00494   farright = vertex( point2<rational>(far_right,Sa->val) );
00495   hori = line( farleft.v(), rational(ULong(1)), rational() );
00496  
00497   E = (graph_edge **)alloca( sizeof(graph_edge *)*maxE );
00498   if( E == NULL ) throw AllocError();
00499   
00500   rec = Sa->R.head;
00501   fa=fb=nfa=nfb=0;
00502                                    /* new vertex for strip above cutting line*/
00503   v = new graph_vertex( point2<rational>(minx,Sa->val), NULL );
00504   Sa->V.Append(v);
00505   lva = *minSaV = v;
00506                                   /* new vertex for strip below cutting line */
00507   v = new graph_vertex(lva->v, NULL);     
00508   Sb->V.Append(v);
00509   lvb = *minSbV = v;
00510   
00511   delv = 0;
00512   while( rec != NULL )
00513   {
00514     if( rec->v != NULL )    /* vertex on cutting line */
00515     {
00516       nE = IncidentEdges( rec->v, E );
00517 
00518       match1 = EdgeInsertPosition( rec->v->v, farleft, nE, E, 
00519                                    &L1, &R1, &cL1, &cR1 );
00520       cL1 = 2-(cL1-2);  cR1 = 2-(cR1-2);
00521       
00522       if( cL1 == 1 )      /* special case: E[L1] below cutting line */
00523       {
00524         bh= E[L1];
00525         er = al = ah = NULL;
00526         
00527         if( cR1 == 2 )  
00528         {
00529           el = E[R1];
00530           bl = E[R1]->e[0];
00531         }
00532         else
00533         {
00534           el = NULL;
00535           bl = E[R1];
00536         }
00537       }
00538       else if( cL1 == 2 ) /* special case: E[L1] on cutting line */
00539       {
00540         er = E[L1];
00541         al = ah = NULL;
00542 
00543         if( cR1 == 2 )  
00544         {
00545           el = E[R1];
00546           if( el->e[0] == er )
00547             bl = NULL;
00548           else
00549             bl = el->e[0];
00550         }
00551         else
00552         {
00553           bl = E[R1];
00554           el = NULL;
00555         }
00556         bh = bl;
00557         if( bh != NULL )
00558           while( bh->e[1] != er ) bh = bh->e[1];
00559       }
00560       else                   /* E[L1] above cutting line */
00561       {
00562         ah = E[L1];
00563         
00564         if( cR1 == 3 )  /* E[R1] above cutting line */
00565         {
00566           el = bl = bh = er = NULL;
00567           al = E[R1];
00568         }
00569         else if( (cR1 == 2) && (!match1) ) /* E[R1] on cutting line     */
00570         {                                  /* but different orientation */
00571           el = bl = bh = NULL;
00572           er = E[R1];
00573           al = E[R1]->e[1];     /* using edge orientatation and
00574                                    that cutting line is horizontal ! */
00575         }
00576         else
00577         {
00578           match2 = EdgeInsertPosition( rec->v->v, farright, nE, E, 
00579                                        &L2, &R2, &cL2, &cR2 );
00580           if( match1 )   /* E[R1] on cutting line, same orientation */
00581           {
00582             el = E[R1];
00583             if( cL2 == 1 ) 
00584             {
00585               bl = E[R1]->e[0];
00586               bh = E[L2];
00587             }
00588             else
00589               bl = bh = NULL;
00590           }
00591           else
00592           {
00593             el = NULL;            
00594             bl = E[R1];
00595             bh = E[L2];
00596           } 
00597          
00598           if( cR2 == 2 )  /* R2 on cutting line */
00599           {
00600             er = E[R2];
00601             al = E[R2]->e[1];
00602           }
00603           else
00604           {
00605             er = NULL;
00606             al = E[R2];
00607           }
00608         }
00609       }
00610      
00611       if( el != NULL )
00612       {
00613         fb = el->f[0].i;
00614         fa = el->f[1].i;
00615       }
00616 
00617       if( al != NULL )
00618       {
00619                                       /* create new vertex: */
00620         v = new graph_vertex( ah->v[0]->v, ah );
00621         Sa->V.Append(v);
00622                 
00623         ah->e[0] = al;   /* connect upper edges to new vertex */
00624         e = al;           
00625         do
00626         {
00627           e->v[0] = v;
00628           e = e->e[0];
00629         }
00630         while( e != al );
00631 
00632         if( el == NULL )
00633           fa = ah->f[0].i;
00634                           /* create new edge */
00635         e = new graph_edge();
00636         Sa->E.Append(e);
00637 
00638         e->l = hori;
00639           
00640         e->f[0].set_i(fb);
00641         e->f[1].set_i(fa);
00642 
00643         e->e[0] = al;
00644         ah->e[0] = e;
00645 
00646         if( lva->e != NULL )
00647         {
00648           e->e[1] = lva->e->e[0];
00649           lva->e->e[0] = e;
00650         }
00651         else
00652         {
00653           e->e[1] = e;
00654           lva->e = e;
00655         }
00656         v->e = e;
00657       
00658         e->v[0] = v;
00659         e->v[1] = lva;
00660           
00661         lva = v;
00662         nfa = nfb = al->f[1].i;
00663       }
00664 
00665       if( bl != NULL )
00666       {
00667                                  /* create new vertex: */
00668         v = new graph_vertex( bl->v[1]->v, bh );
00669         Sb->V.Append(v);
00670 
00671         bh->e[1] = bl;   /* connect edges below to new vertex */
00672         e = bl;           
00673         do
00674         {
00675           Sa->E.Remove(e);
00676           Sb->E.Append(e);
00677           e->v[1] = v;
00678           e = e->e[1];
00679         }
00680         while( e != bl );
00681 
00682         if( el == NULL )
00683           fb = bl->f[0].i;
00684                                  /* create new edge */
00685         e = new graph_edge();
00686         Sb->E.Append(e);
00687 
00688         e->l = hori;
00689         
00690         e->f[0].set_i(fb);
00691         e->f[1].set_i(fa);
00692 
00693         e->e[0] = bl;
00694         bh->e[1] = e;
00695 
00696         if( lvb->e != NULL )
00697         {
00698           e->e[1] = lvb->e->e[1];
00699           lvb->e->e[1] = e;
00700         }
00701         else
00702         {
00703           e->e[1] = e;
00704           lvb->e = e;
00705         } 
00706       
00707         e->v[0] = v;
00708         e->v[1] = lvb;
00709         
00710         lvb = v;
00711         nfa = nfb = bh->f[1].i;
00712       }
00713         
00714       if( er != NULL )
00715       {
00716         fa = er->f[1].i;
00717         fb = er->f[0].i;
00718       }
00719       else
00720       {
00721         fa = nfa;
00722         fb = nfb;
00723       }
00724 
00725       Sa->V.Remove( rec->v );            /* remove original vertex */
00726       delete delv;               /* delayed deletion of vertex from original graph */
00727       delv = rec->v;
00728       
00729       if( el != NULL ) {          /* remove original horizontal edge */
00730         Sa->E.Remove(el);
00731         delete el;
00732       }
00733     }
00734     else           /* cutting line intersects edge */
00735     {
00736       ew = rec->e;
00737                                   /* new vertex for strip above cutting line */
00738       vu = new graph_vertex( point2<rational>(rec->val,Sa->val), NULL );
00739       Sa->V.Append(vu);
00740                                 /* new vertex for strip below cutting line */
00741       vl = new graph_vertex( vu->v, NULL );
00742       Sb->V.Append(vl);
00743                             /* new edge for upper halve of intersecting edge */
00744       eu = new graph_edge();
00745       Sa->E.Append(eu);
00746 
00747       eu->l = ew->l;
00748 
00749       eu->f[0] = ew->f[0];
00750       eu->f[1] = ew->f[1];
00751       eu->v[0] = vu;
00752       eu->v[1] = ew->v[1];
00753 
00754       eu->e[0] = NULL;
00755       if( ew->e[1] != ew )
00756       {
00757         eu->e[1] = ew->e[1];
00758         nE = EdgeCycle( ew, ew->v[1], E ); e = E[nE-1];
00759         if( e->e[0] == ew ) e->e[0] = eu; else e->e[1] = eu;
00760       }
00761       else eu->e[1] = eu;
00762       ew->v[1]->e = eu;
00763         
00764                      /* shorten intersecting edge and add it to lower strip */
00765       Sa->E.Remove(ew);
00766       Sb->E.Append(ew);
00767 
00768       ew->v[1] = vl;
00769       ew->e[1] = NULL;
00770       
00771       vu->e = eu;         /* set vertex pointers to an incident edge */
00772       vl->e = ew;
00773                               /* create new horizontal edge for upper strip */
00774       e = new graph_edge();
00775       Sa->E.Append(e);
00776 
00777       e->l = hori;
00778         
00779       e->f[0] = e->f[1] = ew->f[0]; /* same face below and above */
00780 
00781       e->e[0] = eu;
00782       eu->e[0] = e;
00783 
00784       if( lva->e != NULL )
00785       {
00786         e->e[1] = lva->e->e[0];
00787         lva->e->e[0] = e;
00788       }
00789       else
00790       {
00791         e->e[1] = e;
00792         lva->e = e;
00793       }
00794       vu->e = e;
00795       
00796       e->v[0] = vu;
00797       e->v[1] = lva;
00798                               /* create new horizontal edge for lower strip */
00799       e = new graph_edge();
00800       Sb->E.Append(e);
00801       
00802       e->l = hori;
00803       
00804       e->f[0] = e->f[1] = ew->f[0]; /* same face below and above */
00805 
00806       e->e[0] = ew;
00807       ew->e[1] = e;
00808 
00809       if( lvb->e != NULL )
00810       {
00811         e->e[1] = lvb->e->e[1];
00812         lvb->e->e[1] = e;
00813       }
00814       else
00815       {
00816         e->e[1] = e;
00817         lvb->e = e;
00818       }
00819       vl->e = ew;
00820       
00821       e->v[0] = vl;
00822       e->v[1] = lvb;
00823 
00824       lva = vu;
00825       lvb = vl;
00826       fa = fb = ew->f[1].i;
00827     }
00828 
00829     rec = rec->next;
00830   }
00831   delete delv;
00832                                   /* new vertex for strip above cutting line */
00833   v = new graph_vertex( point2<rational>(maxx,Sa->val), NULL );
00834   *maxSaV = v;
00835   Sa->V.Append(v);
00836                                /* create new horizontal edge for upper strip */
00837   e = new graph_edge();
00838   Sa->E.Append(e);
00839   
00840   e->l = hori;
00841     
00842   e->f[0].set_i(fb);
00843   e->f[1].set_i(fa);
00844 
00845   e->e[0] = e;
00846 
00847   if( lva->e != NULL )
00848   {
00849     e->e[1] = lva->e->e[0];
00850     lva->e->e[0] = e;
00851   }
00852   else
00853   {
00854     e->e[1] = e;
00855     lva->e = e;
00856   }
00857   v->e = e;
00858       
00859   e->v[0] = v;
00860   e->v[1] = lva;
00861                                   /* new vertex for strip below cutting line */
00862   v = new graph_vertex( v->v, NULL );
00863   *maxSbV = v;
00864   Sb->V.Append(v);
00865                                 /* create new horizontal edge for lower strip */
00866   e = new graph_edge();
00867   Sb->E.Append(e);
00868 
00869   e->l = hori;
00870   
00871   e->f[0].set_i(fb);
00872   e->f[1].set_i(fa);
00873 
00874   e->e[0] = e;
00875 
00876   if( lvb->e != NULL )
00877   {
00878     e->e[1] = lvb->e->e[1];
00879     lvb->e->e[1] = e;
00880   }
00881   else
00882   {
00883     e->e[1] = e;
00884     lvb->e = e;
00885   } 
00886   v->e = e;
00887       
00888   e->v[0] = v;
00889   e->v[1] = lvb; 
00890 }
00891 
00892 
00893 
00894 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00895 
00896   Functions for splitting a graph into vertical strips
00897 
00898 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00899 
00900 /* -------------------------------------------------------------------------
00901   Calculates the intersection of a graph with vertical cutting lines 
00902 -------------------------------------------------------------------------- */
00903 
00904 void IntersectWithVerticals( 
00905           graph_edge_list *E, graph_vertex_list *V,
00906           const int nS, Ptr< cut_info > S )
00907 {
00908   int ns = nS-1, i,j,k,ii,ik,sign,nA;
00909   graph_vertex *v,*v_next;
00910   graph_edge *e, *e_next;
00911   rational y,x0,y0,x,dx,dy;
00912   bool rz;
00913   cut_record *rec;
00914   Ptr< cut_info > A;
00915 
00916   v = V->head;
00917 
00918   nA = nS-1;            
00919   A.use_pointer( &S[1], nS-1 );  // used for binary search
00920 
00921   while( v != NULL )
00922   {
00923     v_next = v->next;
00924 
00925     x = v->v.v().x;
00926   
00927     i = FindCutInfo( x, nA, A, rz ) + 1;
00928 
00929     v->reserved[0].set_i(i);
00930     v->reserved[1].set_i(rz);
00931     
00932     V->Remove( v );
00933     S[i].V.Append( v );
00934 
00935     if( rz )
00936     { 
00937       rec = new cut_record();
00938       S[i].R.Append( rec );
00939            
00940       rec->v = v;
00941       rec->e = NULL;
00942       rec->val = v->v.v().y;
00943     }
00944 
00945     v = v_next;
00946   }
00947 
00948   e = E->head;
00949 
00950   while( e != NULL )  
00951   {
00952     e_next = e->next;
00953 
00954     sign = compare(e->v[1]->v.v().x, e->v[0]->v.v().x);
00955     if( sign >= 0 ) {
00956       ii = 0;
00957       ik = 1;
00958     } else {
00959       ii = 1;
00960       ik = 0;
00961     }    
00962 
00963     i = e->v[ii]->reserved[0].i;
00964     k = e->v[ik]->reserved[0].i;
00965 
00966     E->Remove( e );
00967     S[k].E.Append( e );
00968 
00969     if( i !=  k )
00970     {
00971       if( e->l.IsNULL() )
00972       {
00973         e->CalcLine();
00974       }
00975 
00976       x0 = e->l.v().x;
00977       y0 = e->l.v().y;
00978     
00979       dy = e->l.dy();
00980 
00981       if( e->v[ik]->reserved[1].i ) k--;
00982 
00983       for( j = i+1; j <= k; j++ )
00984       {
00985         if( !dy.is_zero() )
00986           y = y0 + (S[j].val-x0)*dy;
00987         else
00988           y = y0;
00989           
00990         rec = new cut_record();
00991         S[j].R.Append( rec );
00992               
00993         rec->v = NULL;
00994         rec->e = e;          
00995         rec->val = y;
00996       }
00997     }
00998 
00999     e = e_next;
01000   }
01001 
01002   sort_cut_records( ns, S );
01003 }
01004 
01005 /* -------------------------------------------------------------------------
01006   Divide graph into an independant left and right part
01007 -------------------------------------------------------------------------- */
01008  
01009 void DivideVerticalStrips( 
01010               cut_info *Sb, cut_info *Sa, int maxE,
01011               const rational& miny, const rational maxy,
01012               const rational& far_bottom, const rational& far_top,
01013               graph_vertex **minSbV, graph_vertex **minSaV,
01014               graph_vertex **maxSbV, graph_vertex **maxSaV )
01015 {
01016   graph_edge **E;
01017   cut_record *rec;
01018   int nE,L1,L2,R1,R2,match1,match2,cL1,cR1,cL2,cR2,iv,sign,ew_ori;
01019   vertex farbottom,fartop;
01020   graph_edge *ah,*el,*bl,*bh,*er,*al, *e, *eu, *ew;
01021   graph_vertex *v, *lva, *lvb, *vu, *vl, *delv;
01022   int fa,fb,nfa,nfb;
01023   line verti;
01024 
01025   farbottom = vertex( point2<rational>(Sb->val,far_bottom) );
01026   fartop = vertex( point2<rational>(Sb->val,far_top) );
01027   verti = line( farbottom.v(), rational(), rational(ULong(1)) );
01028  
01029   E = (graph_edge **)alloca( sizeof(graph_edge *)*maxE );
01030   if( E == NULL ) throw AllocError();
01031   
01032   rec = Sb->R.head;
01033   fa=fb=nfa=nfb=0;
01034                                 /* new vertex for strip left from cutting line*/
01035   v = *minSaV = new graph_vertex( point2<rational>(Sb->val,miny), NULL );
01036   Sa->V.Append(v);
01037   lva = v;
01038                               /* new vertex for strip right from cutting line */
01039   v = *minSbV = new graph_vertex(lva->v, NULL);     
01040   Sb->V.Append(v);
01041   lvb = v;
01042 
01043   delv = 0;  
01044   while( rec != NULL )
01045   {
01046     if( rec->v != NULL )    /* vertex on cutting line */
01047     {
01048       nE = IncidentEdges( rec->v, E );
01049 
01050       match1 = EdgeInsertPosition( rec->v->v, farbottom, nE, E, 
01051                                    &L1, &R1, &cL1, &cR1 );
01052       cL1 = 2-(cL1-2);  cR1 = 2-(cR1-2);
01053      
01054       if( cL1 == 1 )      /* special case: E[L1] left from cutting line */
01055       {
01056         bh = E[L1];
01057         er = al = ah = NULL;
01058         
01059         if( cR1 == 2 )  
01060         {
01061           el = E[R1];
01062           bl = E[R1]->e[1];
01063         }
01064         else
01065         {
01066           el = NULL;
01067           bl = E[R1];
01068         }
01069       }
01070       else if( cL1 == 2 ) /* special case: E[L1] on cutting line */
01071       {
01072         er = E[L1];
01073         al = ah = NULL;
01074         
01075         if( cR1 == 2 )  
01076         {
01077           el = E[R1];
01078           if( el->e[1] == er )
01079             bl = NULL;
01080           else
01081             bl = el->e[1];
01082         }
01083         else
01084         {
01085           bl = E[R1];
01086           el = NULL;
01087         }        
01088         bh = bl;
01089         if( bh != NULL )
01090         {
01091           e = (bh->v[0] == rec->v ? bh->e[0] : bh->e[1]);
01092           while( e != er ) 
01093           {
01094             bh = e;
01095             e = (bh->v[0] == rec->v ? bh->e[0] : bh->e[1]);
01096           }
01097         }
01098       }
01099       else    /* E[L1] right from cutting line */
01100       {
01101         ah = E[L1];
01102         
01103         if( cR1 == 3 )  /* E[R1] right from cutting line */
01104         {
01105           el = bl = bh = er = NULL;
01106           al = E[R1];
01107         }
01108         else if( (cR1 == 2) && (!match1) ) /* E[R1] on cutting line     */
01109         {                                  /* but different orientation */
01110           el = bl = bh = NULL;
01111           er = E[R1];
01112           al = E[R1]->e[0];     /* using edge orientatation and
01113                                    that cutting line is vertical ! */
01114         }
01115         else
01116         {
01117           match2 = EdgeInsertPosition( rec->v->v, fartop, nE, E, 
01118                                        &L2, &R2, &cL2, &cR2 );
01119           if( match1 )   /* E[R1] on cutting line, same orientation */
01120           {
01121             el = E[R1];
01122             if( cL2 == 1 ) 
01123             {
01124               bl = E[R1]->e[1];
01125               bh = E[L2];
01126             }
01127             else
01128               bl = bh = NULL;
01129           }
01130           else
01131           {
01132             el = NULL;            
01133             bl = E[R1];
01134             bh = E[L2];
01135           } 
01136          
01137           if( cR2 == 2 )  /* R2 on cutting line */
01138           {
01139             er = E[R2];
01140             al = E[R2]->e[0];
01141           }
01142           else
01143           {
01144             er = NULL;
01145             al = E[R2];
01146           }
01147         }
01148       }
01149      
01150       if( el != NULL )
01151       {
01152         fb = el->f[1].i;
01153         fa = el->f[0].i;
01154       }
01155 
01156       if( al != NULL )
01157       {
01158                                       /* create new vertex: */
01159         iv = (ah->v[0] == rec->v ? 0 : 1); 
01160         v = new graph_vertex( rec->v->v, ah );
01161         Sa->V.Append(v);
01162                 
01163         ah->e[iv] = al;   /* connect edges to the left to new vertex */
01164         e = al;           
01165         do
01166         {
01167           Sb->E.Remove(e);
01168           Sa->E.Append(e);
01169           if( e->v[0] == rec->v )
01170           { e->v[0] = v;
01171             e = e->e[0];
01172           } else
01173           { e->v[1] = v;
01174             e = e->e[1];
01175           }
01176         }
01177         while( e != al );
01178 
01179         if( el == NULL )
01180           fa = ah->f[iv].i;
01181                           /* create new edge */
01182         e = new graph_edge();
01183         Sa->E.Append(e);
01184 
01185         e->l = verti;
01186           
01187         e->f[0].set_i(fa);
01188         e->f[1].set_i(fb);
01189 
01190         e->e[1] = al;
01191         ah->e[iv] = e;
01192 
01193         if( lva->e != NULL )
01194         {
01195           if( lva->e->v[0] == lva ) {
01196             e->e[0] = lva->e->e[0];
01197             lva->e->e[0] = e;
01198           } else {
01199             e->e[0] = lva->e->e[1];
01200             lva->e->e[1] = e;
01201           }
01202         }
01203         else
01204         {
01205           e->e[0] = e;
01206           lva->e = e;
01207         }
01208         v->e = e;
01209       
01210         e->v[1] = v;
01211         e->v[0] = lva;
01212           
01213         lva = v;
01214         nfa = nfb = (al->v[0] == v ? al->f[1].i : al->f[0].i);
01215       }
01216 
01217       if( bl != NULL )
01218       {
01219                                  /* create new vertex: */
01220         v = new graph_vertex( rec->v->v, bh );
01221         Sb->V.Append(v);
01222                             /* connect edges on the right to new vertex */
01223         if( bh->v[0] == rec->v )
01224           bh->e[0] = bl;   
01225         else        
01226           bh->e[1] = bl; 
01227         e = bl;           
01228         do
01229         {
01230           if( e->v[0] == rec->v )
01231           { e->v[0] = v;
01232             e = e->e[0];       
01233           } else {
01234             e->v[1] = v;
01235             e = e->e[1];
01236           }
01237         }
01238         while( e != bl );
01239 
01240         if( el == NULL )
01241           if( bl->v[1] == v )
01242             fb = bl->f[0].i;
01243           else
01244             fb = bl->f[1].i;
01245                                  /* create new edge */
01246         e = new graph_edge();
01247         Sb->E.Append(e);
01248 
01249         e->l = verti;
01250         
01251         e->f[0].set_i(fa);
01252         e->f[1].set_i(fb);
01253 
01254         e->e[1] = bl;
01255         if( bh->v[0] == v ) bh->e[0] = e; else bh->e[1] = e;
01256 
01257         if( lvb->e != NULL )
01258         {
01259           if( lvb->e->v[0] == lvb ) {
01260             e->e[0] = lvb->e->e[0];
01261             lvb->e->e[0] = e;
01262           } else {
01263             e->e[0] = lvb->e->e[1];
01264             lvb->e->e[1] = e;
01265           }
01266         }
01267         else
01268         {
01269           e->e[0] = e;
01270           lvb->e = e;
01271         } 
01272       
01273         e->v[1] = v;
01274         e->v[0] = lvb;
01275         
01276         lvb = v;
01277         nfa = nfb = (bh->v[1] == v ? bh->f[1].i : bh->f[0].i);
01278       }
01279 
01280       Sb->V.Remove( rec->v );            /* remove original vertex */
01281       delete delv;                /* delayed deletion of vertex from original graph */
01282       delv = rec->v;
01283         
01284       if( er != NULL )
01285       {
01286         fa = er->f[0].i;
01287         fb = er->f[1].i;
01288       }
01289       else
01290       {
01291         fa = nfa;
01292         fb = nfb;
01293       } 
01294 
01295       if( el != NULL ) {                 /* remove original vertical edge */
01296         Sb->E.Remove(el); 
01297         delete el;
01298       }
01299     }
01300     else           /* cutting line intersects edge */
01301     {
01302       ew = rec->e;
01303                               /* new vertex for strip right from cutting line */
01304       vu = new graph_vertex( point2<rational>(Sb->val,rec->val), NULL );
01305       Sb->V.Append(vu);
01306                              /* new vertex for strip left from cutting line */
01307       vl = new graph_vertex( vu->v, NULL );
01308       Sa->V.Append(vl);
01309                             /* new edge for right halve of intersecting edge */
01310       eu = new graph_edge();
01311       Sb->E.Append(eu);
01312 
01313       eu->l = ew->l;
01314 
01315       eu->f[0] = ew->f[0];
01316       eu->f[1] = ew->f[1];
01317 
01318       sign = compare(ew->v[1]->v.v().x, ew->v[0]->v.v().x);
01319       if( sign > 0 ) ew_ori = 1; else ew_ori = 0; 
01320 
01321       if( ew_ori )                    /*** ew oriented from left to right ***/
01322       {
01323         eu->v[0] = vu;
01324         eu->v[1] = ew->v[1];
01325 
01326         eu->e[0] = NULL;
01327         if( ew->e[1] != ew )
01328         {
01329           eu->e[1] = ew->e[1];
01330           nE = EdgeCycle( ew, ew->v[1], E ); e = E[nE-1];
01331           if( e->e[0] == ew ) e->e[0] = eu; else e->e[1] = eu;
01332         }
01333         else eu->e[1] = eu;
01334         ew->v[1]->e = eu;
01335 
01336                      /* shorten intersecting edge and add it to left strip */
01337         Sb->E.Remove(ew);
01338         Sa->E.Append(ew);
01339 
01340         ew->v[1] = vl;
01341         ew->e[1] = NULL;
01342       
01343         vu->e = eu;         /* set vertex pointers to an incident edge */
01344         vl->e = ew;
01345                               /* create new vertical edge for right strip */
01346         e = new graph_edge();
01347         Sb->E.Append(e);
01348 
01349         e->l = verti;
01350         
01351         e->f[0] = e->f[1] = ew->f[1]; /* same face on left and right side */
01352                                       /* of edge */
01353         eu->e[0] = e;
01354       }
01355       else                            /*** ew oriented from right to left ***/
01356       {
01357         eu->v[1] = vu;
01358         eu->v[0] = ew->v[0];
01359 
01360         eu->e[1] = NULL;
01361         if( ew->e[0] != ew )
01362         {
01363           eu->e[0] = ew->e[0];
01364           nE = EdgeCycle( ew, ew->v[0], E ); e = E[nE-1];
01365           if( e->e[0] == ew ) e->e[0] = eu; else e->e[1] = eu;
01366         }
01367         else eu->e[0] = eu;
01368         ew->v[0]->e = eu;
01369 
01370                      /* shorten intersecting edge and add it to lower strip */
01371         Sb->E.Remove(ew);
01372         Sa->E.Append(ew);
01373 
01374         ew->v[0] = vl;
01375         ew->e[0] = NULL;
01376       
01377         vu->e = eu;         /* set vertex pointers to an incident edge */
01378         vl->e = ew;
01379                               /* create new vertical edge for right strip */
01380         e = new graph_edge();
01381         Sb->E.Append(e);
01382 
01383         e->l = verti;
01384         
01385         e->f[0] = e->f[1] = ew->f[0]; /* same face on left and right side */
01386                                       /* of edge */
01387         eu->e[1] = e;
01388       } 
01389 
01390       e->e[1] = eu;
01391 
01392       if( lvb->e != NULL )
01393       {
01394         if( lvb->e->v[0] == lvb ) {
01395           e->e[0] = lvb->e->e[0];
01396           lvb->e->e[0] = e;
01397         } else {
01398           e->e[0] = lvb->e->e[1];
01399           lvb->e->e[1] = e;
01400         }
01401       }
01402       else
01403       {
01404         e->e[0] = e;
01405         lvb->e = e;
01406       }
01407       vu->e = eu;
01408 
01409       e->v[1] = vu;
01410       e->v[0] = lvb;
01411                               /* create new vertical edge for left strip */
01412       e = new graph_edge();
01413       Sa->E.Append(e);
01414       
01415       e->l = verti;
01416       
01417       if( ew_ori )                  /*** ew oriented from left to right ***/
01418       {
01419         e->f[0] = e->f[1] = ew->f[1]; /* same face on left and right side */
01420 
01421         e->e[1] = ew;
01422         ew->e[1] = e;
01423 
01424         fa = fb = ew->f[0].i;      /* for next intersection */
01425       }
01426       else                          /*** ew oriented from left to right ***/
01427       {
01428         e->f[0] = e->f[1] = ew->f[0]; /* same face on left and right side */
01429 
01430         e->e[1] = ew;
01431         ew->e[0] = e;
01432 
01433         fa = fb = ew->f[1].i;       /* for next intersection */
01434       }
01435       
01436       if( lva->e != NULL )
01437       {
01438         if( lva->e->v[0] == lva ) {
01439           e->e[0] = lva->e->e[0];
01440           lva->e->e[0] = e;
01441         } else {
01442           e->e[0] = lva->e->e[1];
01443           lva->e->e[1] = e;
01444         }
01445       }
01446       else
01447       {
01448         e->e[0] = e;
01449         lva->e = e;
01450       }
01451       vl->e = e;
01452       
01453       e->v[1] = vl;
01454       e->v[0] = lva;
01455 
01456       lvb = vu;
01457       lva = vl;
01458     }
01459 
01460     rec = rec->next;
01461   }
01462   delete delv;
01463                                /* new vertex for strip left from cutting line */
01464   v = new graph_vertex( point2<rational>(Sb->val,maxy), NULL );
01465   *maxSaV = v;
01466   Sa->V.Append(v);
01467                                /* create new vertical edge for left strip */
01468   e = new graph_edge();
01469   Sa->E.Append(e);
01470   
01471   e->l = verti;
01472     
01473   e->f[0].set_i(fa);
01474   e->f[1].set_i(fb);
01475 
01476   e->e[1] = e;
01477 
01478   if( lva->e != NULL )
01479   {
01480     if( lva->e->v[0] == lva ) {
01481       e->e[0] = lva->e->e[0];
01482       lva->e->e[0] = e;
01483     } else {
01484       e->e[0] = lva->e->e[1];
01485       lva->e->e[1] = e;
01486     }
01487   }
01488   else
01489   {
01490     e->e[0] = e;
01491     lva->e = e;
01492   }
01493   v->e = e;
01494       
01495   e->v[1] = v;
01496   e->v[0] = lva;
01497                              /* new vertex for strip right from cutting line */
01498   v = new graph_vertex( v->v, NULL );
01499   *maxSbV = v;
01500   Sb->V.Append(v);
01501                                 /* create new vertical edge for right strip */
01502   e = new graph_edge();
01503   Sb->E.Append(e);
01504 
01505   e->l = verti;
01506   
01507   e->f[0].set_i(fa);
01508   e->f[1].set_i(fb);
01509 
01510   e->e[1] = e;
01511 
01512   if( lvb->e != NULL )
01513   {
01514     if( lvb->e->v[0] == lvb ) {
01515       e->e[0] = lvb->e->e[0];
01516       lvb->e->e[0] = e;
01517     } else {
01518       e->e[0] = lvb->e->e[1];
01519       lvb->e->e[1] = e;
01520     }
01521   }
01522   else
01523   {
01524     e->e[0] = e;
01525     lvb->e = e;
01526   } 
01527   v->e = e;
01528   
01529   e->v[1] = v;
01530   e->v[0] = lvb; 
01531 }
01532 
01533 /* ---------------------------------------------------------------------------
01534   splits a graph into four independent parts
01535 ---------------------------------------------------------------------------- */
01536 template< class T >
01537 void SplitGraph( GraphInfo *G, T xmed, T ymed,
01538                  GraphConvInfo<T> *Gconv, GraphInfo4& sub )
01539 {
01540   graph_vertex *Vy1Sl, *Vy1Sr, *Vy2Sl, *Vy2Sr, *Vx1Sa, *Vx1Sb, *Vx2Sa, *Vx2Sb;
01541   graph_vertex *hvLeftT, *hvRightT;
01542   graph_edge *e, *el, **hE;
01543   graph_vertex *v, *v0;
01544   int nhE, nMaxE, i, j;
01545   rational X1,X2,X3;
01546   rational Y1,Y2,Y3;
01547   unsigned long *cbuf = (unsigned long *)alloca(gul::rtr<double>::mantissa_length());
01548   T y2,x2;
01549   Ptr< gugr::cut_info > CutsX, CutsY;
01550 
01551   /*
01552   cout << "SPLITGRAPH: input graph\n";
01553   cout << "**********************\n";
01554   gugr::Dump<T>::dump_vertices( G->V.head );
01555   gugr::Dump<T>::dump_edges( G->E.head );
01556   */
01557 
01558   rational& MinX = Gconv->MinX;    // short aliases
01559   rational& MaxX = Gconv->MaxX;
01560   rational& FarMinX = Gconv->FarMinX;
01561   rational& FarMaxX = Gconv->FarMaxX;
01562   rational& MinY = Gconv->MinY;
01563   rational& MaxY = Gconv->MaxY;
01564   rational& FarMinY = Gconv->FarMinY;
01565   rational& FarMaxY = Gconv->FarMaxY;
01566  
01567   CutsX.reserve_place( reserve_stack(cut_info,4), 4 );
01568   CutsY.reserve_place( reserve_stack(cut_info,4), 4 );
01569 
01570   x2 = (xmed - Gconv->minx)/Gconv->scale;
01571   X2 = rational(coord2int(x2,cbuf),cbuf);
01572   y2 = (ymed - Gconv->miny)/Gconv->scale;
01573   Y2 = rational(coord2int(y2,cbuf),cbuf);
01574 
01575   CutsX[1].val = X1 = G->P[0][0]->v.v().x;
01576   CutsX[2].val = X2;
01577   CutsX[3].val = X3 = G->P[0][1]->v.v().x;
01578 
01579   CutsY[1].val = Y1 = G->P[0][0]->v.v().y;
01580   CutsY[2].val = Y2;
01581   CutsY[3].val = Y3 = G->P[1][0]->v.v().y;
01582 
01583   nMaxE = 10*G->E.nElems;      // this is a hack, but its enough even in the
01584                                // in the worst case
01585   hE = (graph_edge **)alloca( sizeof(graph_edge *) * nMaxE );
01586   
01587                                            // split graph into vertical strips
01588   IntersectWithVerticals( &G->E, &G->V, 4, CutsX );
01589 
01590                                 // isolate first vertical strip and discard it
01591   DivideVerticalStrips( &CutsX[3], &CutsX[2], nMaxE,
01592                         MinY, MaxY, FarMinY, FarMaxY,
01593                         &Vy1Sr, &Vy1Sl, &Vy2Sr, &Vy2Sl );  
01594 
01595   /*
01596   cout << "after disconnecting X-Strip " << 3 << " and " << 2 << "\n";
01597   cout << "*******************************************************\n";
01598   for( int k = 3; k >= 0; k-- )
01599   {
01600     cout << "X-Strip " << k << "\n";
01601     gugr::Dump<T>::dump_vertices( CutsX[k].V.head );
01602     gugr::Dump<T>::dump_edges( CutsX[k].E.head );
01603   }
01604   */
01605 
01606   CutsX[3].R.DeleteElems();
01607   CutsX[3].E.DeleteElems();
01608   CutsX[3].V.DeleteElems();
01609 
01610   for( i = 2; i > 0; i-- )
01611   {
01612     DivideVerticalStrips( &CutsX[i], &CutsX[i-1], nMaxE,
01613                           MinY, MaxY, FarMinY, FarMaxY,
01614                           &Vy1Sr, &Vy1Sl, &Vy2Sr, &Vy2Sl );
01615 
01616     /*
01617     cout << "after disconnecting X-Strip " << i << " and " << i-1 << "\n";
01618     cout << "*********************************************************\n";
01619     for( int k = i; k >= 0; k-- )
01620     {
01621       cout << "X-Strip " << k << "\n";
01622       gugr::Dump<T>::dump_vertices( CutsX[k].V.head );
01623       gugr::Dump<T>::dump_edges( CutsX[k].E.head );
01624     }
01625     */
01626                                // split vertical strips into horizontal strips 
01627     IntersectWithHorizontals( &CutsX[i].E, &CutsX[i].V, 4, CutsY );
01628 
01629                                          // isolate first strip and discard it
01630     DivideHorizontalStrips( &CutsY[3], &CutsY[2], nMaxE,
01631                             MinX, MaxX, FarMinX, FarMaxX,
01632                             &Vx1Sa, &Vx1Sb, &Vx2Sa, &Vx2Sb );   
01633 
01634     /*
01635     cout << "after disconnecting Y-Strip " << 3 << " and " << 2 << "\n";
01636     cout << "*******************************************************\n";
01637     for( int k = 3; k >= 0; k-- )
01638     {
01639       cout << "Y-Strip " << k << "\n";
01640       gugr::Dump<T>::dump_vertices( CutsY[k].V.head );
01641       gugr::Dump<T>::dump_edges( CutsY[k].E.head );
01642     }
01643     */
01644   
01645     CutsY[3].R.DeleteElems();
01646     CutsY[3].E.DeleteElems();
01647     CutsY[3].V.DeleteElems();
01648 
01649     hvLeftT = Vx1Sb; hvRightT = Vx2Sb;  // for next strip
01650     
01651                                // divide vertical strips into horizontal strips    
01652     for( j = 2; j > 0; j-- )
01653     {
01654       DivideHorizontalStrips( &CutsY[j], &CutsY[j-1], nMaxE,
01655                               MinX, MaxX, FarMinX, FarMaxX,
01656                               &Vx1Sa, &Vx1Sb, &Vx2Sa, &Vx2Sb );
01657 
01658       /* 
01659       cout << "after disconnecting Y-Strip " << j << " and " << j-1 << "\n";
01660       cout << "*********************************************************\n";
01661       for( int k = j; k >= 0; k-- )
01662       {
01663         cout << "Y-Strip " << k << "\n";
01664         gugr::Dump<T>::dump_vertices( CutsY[k].V.head );
01665         gugr::Dump<T>::dump_edges( CutsY[k].E.head );
01666       }
01667       */
01668       
01669                               // remove help edges from the corner edge cycles
01670       v0 = hvLeftT;
01671       e = v0->e;
01672       sub[2*(j-1)+i-1]->face = e->e[0]->f[1].i; // face index of quad, if it is
01673       v = e->v[0];                              // not divided
01674       sub[2*(j-1)+i-1]->P[1][0] = v;
01675       nhE = EdgeCycle( e, v, hE );
01676       el = hE[nhE-1];
01677       if( el->v[0] == v ) el->e[0] = e->e[0]; else el->e[1] = e->e[0];
01678       v->e = el;
01679       CutsY[j].V.Remove(v0); CutsY[j].E.Remove(e);
01680       delete v0; delete e;
01681       
01682       v0 = hvRightT;
01683       e = v0->e;
01684       v = e->v[1];
01685       sub[2*(j-1)+i-1]->P[1][1] = v;
01686       nhE = EdgeCycle( e, v, hE );
01687       el = hE[nhE-1];
01688       if( el->v[0] == v ) el->e[0] = e->e[1]; else el->e[1] = e->e[1];
01689       v->e = el;
01690       CutsY[j].V.Remove(v0); CutsY[j].E.Remove(e);
01691       delete v0; delete e;
01692  
01693       v0 = Vx2Sa;
01694       e = v0->e;
01695       v = e->v[1];
01696       sub[2*(j-1)+i-1]->P[0][1] = v;
01697       nhE = EdgeCycle( e, v, hE );
01698       el = hE[nhE-1];
01699       if( el->v[0] == v ) el->e[0] = e->e[1]; else el->e[1] = e->e[1];
01700       v->e = el;
01701       CutsY[j].V.Remove(v0); CutsY[j].E.Remove(e);
01702       delete v0; delete e;
01703  
01704       v0 = Vx1Sa;
01705       e = v0->e;
01706       v = e->v[0];
01707       sub[2*(j-1)+i-1]->P[0][0] = v;
01708       nhE = EdgeCycle( e, v, hE );
01709       el = hE[nhE-1];
01710       if( el->v[0] == v ) el->e[0] = e->e[0]; else el->e[1] = e->e[0];
01711       v->e = el;
01712       CutsY[j].V.Remove(v0); CutsY[j].E.Remove(e);
01713       delete v0; delete e;
01714     
01715       hvLeftT = Vx1Sb;  hvRightT = Vx2Sb;   // for next loop
01716 
01717       CutsY[j].R.DeleteElems();             // free cut records only
01718       sub[2*(j-1)+i-1]->E += CutsY[j].E;
01719       sub[2*(j-1)+i-1]->V += CutsY[j].V;
01720       // CutsY[j].E.ReleaseElems();            // clear for next strip
01721       // CutsY[j].V.ReleaseElems();
01722     }
01723 
01724     CutsY[0].R.DeleteElems();
01725     CutsY[0].E.DeleteElems();
01726     CutsY[0].V.DeleteElems();
01727 
01728     CutsX[i].R.DeleteElems();
01729     CutsX[i].E.DeleteElems();
01730     CutsX[i].V.DeleteElems();
01731   }
01732   
01733   CutsX[0].R.DeleteElems();
01734   CutsX[0].E.DeleteElems();
01735   CutsX[0].V.DeleteElems();
01736 }
01737 template 
01738 void SplitGraph( GraphInfo *G, float xmed, float ymed,
01739                  GraphConvInfo<float> *Gconv, GraphInfo4& sub );
01740 template 
01741 void SplitGraph( GraphInfo *G, double xmed, double ymed,
01742                  GraphConvInfo<double> *Gconv, GraphInfo4& sub );
01743 
01744 /* ---------------------------------------------------------------------------
01745   inits a graph info
01746 ---------------------------------------------------------------------------- */
01747 template< class T >
01748 void InitGraph( 
01749        const rational& X1, const rational& X2,
01750        const rational& Y1, const rational& Y2,
01751        GraphConvInfo<T> *Gconv, graph_edge_list *E, graph_vertex_list *V,
01752        GraphInfo *G )
01753 {
01754 
01755   graph_vertex *Vy1Sl, *Vy1Sr, *Vy2Sl, *Vy2Sr, *Vx1Sa, *Vx1Sb, *Vx2Sa, *Vx2Sb;
01756   graph_vertex *hvLeftT, *hvRightT;
01757   graph_edge *e, *el, **hE;
01758   graph_vertex *v, *v0;
01759   int nhE, nMaxE;
01760   Ptr< cut_info > CutsX, CutsY;
01761 
01762   rational& MinX = Gconv->MinX;    // short aliases
01763   rational& MaxX = Gconv->MaxX;
01764   rational& FarMinX = Gconv->FarMinX;
01765   rational& FarMaxX = Gconv->FarMaxX;
01766   rational& MinY = Gconv->MinY;
01767   rational& MaxY = Gconv->MaxY;
01768   rational& FarMinY = Gconv->FarMinY;
01769   rational& FarMaxY = Gconv->FarMaxY;
01770  
01771   CutsX.reserve_place( reserve_stack(cut_info,3), 3 );
01772   CutsY.reserve_place( reserve_stack(cut_info,3), 3 );
01773 
01774   CutsX[1].val = X1;
01775   CutsX[2].val = X2;
01776 
01777   CutsY[1].val = Y1;
01778   CutsY[2].val = Y2;
01779 
01780   nMaxE = 10 * E->nElems;      // this is a hack, but its enough even in the
01781                                // worst case
01782   hE = (graph_edge **)alloca( sizeof(graph_edge *) * nMaxE );
01783   
01784   /*
01785   cout << "input graph\n";
01786   cout << "***********\n";
01787   gugr::Dump<T>::dump_vertices( V->head );
01788   gugr::Dump<T>::dump_edges( E->head );
01789   */
01790                                            // split graph into vertical strips
01791   IntersectWithVerticals( E, V, 3, CutsX );
01792 
01793                                 // isolate first vertical strip and discard it
01794   DivideVerticalStrips( &CutsX[2], &CutsX[1], nMaxE,
01795                         MinY, MaxY, FarMinY, FarMaxY,
01796                         &Vy1Sr, &Vy1Sl, &Vy2Sr, &Vy2Sl );  
01797   /*
01798   cout << "after disconnecting X-Strip " << 2 << " and " << 1 << "\n";
01799   cout << "*******************************************************\n";
01800   for( int k = 2; k >= 0; k-- )
01801   {
01802     cout << "X-Strip " << k << "\n";
01803     gugr::Dump<T>::dump_vertices( CutsX[k].V.head );
01804     gugr::Dump<T>::dump_edges( CutsX[k].E.head );
01805   }
01806   */
01807 
01808   CutsX[2].R.DeleteElems();
01809   CutsX[2].E.DeleteElems();
01810   CutsX[2].V.DeleteElems();
01811 
01812   DivideVerticalStrips( &CutsX[1], &CutsX[0], nMaxE,
01813                         MinY, MaxY, FarMinY, FarMaxY,
01814                         &Vy1Sr, &Vy1Sl, &Vy2Sr, &Vy2Sl );
01815   /*
01816   cout << "after disconnecting X-Strip " << 1 << " and " << 0 << "\n";
01817   cout << "*******************************************************\n";
01818   for( int k = 1; k >= 0; k-- )
01819   {
01820     cout << "X-Strip " << k << "\n";
01821     gugr::Dump<T>::dump_vertices( CutsX[k].V.head );
01822     gugr::Dump<T>::dump_edges( CutsX[k].E.head );
01823   }
01824   */
01825                                // split vertical strips into horizontal strips 
01826   IntersectWithHorizontals( &CutsX[1].E, &CutsX[1].V, 3, CutsY );
01827 
01828                                          // isolate first strip and discard it
01829   DivideHorizontalStrips( &CutsY[2], &CutsY[1], nMaxE,
01830                           MinX, MaxX, FarMinX, FarMaxX,
01831                           &Vx1Sa, &Vx1Sb, &Vx2Sa, &Vx2Sb );   
01832   /*
01833   cout << "after disconnecting Y-Strip " << 2 << " and " << 1 << "\n";
01834   cout << "*******************************************************\n";
01835   for( int k = 2; k >= 0; k-- )
01836   {
01837     cout << "Y-Strip " << k << "\n";
01838     gugr::Dump<T>::dump_vertices( CutsY[k].V.head );
01839     gugr::Dump<T>::dump_edges( CutsY[k].E.head );
01840   }
01841   */
01842   CutsY[2].R.DeleteElems();
01843   CutsY[2].E.DeleteElems();
01844   CutsY[2].V.DeleteElems();
01845 
01846   hvLeftT = Vx1Sb; hvRightT = Vx2Sb;  // for next strip
01847     
01848                                // divide vertical strips into horizontal strips
01849   DivideHorizontalStrips( &CutsY[1], &CutsY[0], nMaxE,
01850                           MinX, MaxX, FarMinX, FarMaxX,
01851                           &Vx1Sa, &Vx1Sb, &Vx2Sa, &Vx2Sb );
01852   /*
01853   cout << "after disconnecting Y-Strip " << 1 << " and " << 0 << "\n";
01854   cout << "*******************************************************\n";
01855   for( int k = 1; k >= 0; k-- )
01856   {
01857     cout << "Y-Strip " << k << "\n";
01858     gugr::Dump<T>::dump_vertices( CutsY[k].V.head );
01859     gugr::Dump<T>::dump_edges( CutsY[k].E.head );
01860   }
01861   */
01862 
01863   v0 = hvLeftT;           // remove help edges from the corner edge cycles
01864   e = v0->e;
01865   G->face = e->e[0]->f[1].i;        // face index of quad, if it is
01866   v = e->v[0];                    // not divided
01867   G->P[1][0] = v;
01868   nhE = EdgeCycle( e, v, hE );
01869   el = hE[nhE-1];
01870   if( el->v[0] == v ) el->e[0] = e->e[0]; else el->e[1] = e->e[0];
01871   v->e = el;
01872   CutsY[1].V.Remove(v0); CutsY[1].E.Remove(e);
01873   delete v0; delete e;
01874       
01875   v0 = hvRightT;
01876   e = v0->e;
01877   v = e->v[1];
01878   G->P[1][1] = v;
01879   nhE = EdgeCycle( e, v, hE );
01880   el = hE[nhE-1];
01881   if( el->v[0] == v ) el->e[0] = e->e[1]; else el->e[1] = e->e[1];
01882   v->e = el;
01883   CutsY[1].V.Remove(v0); CutsY[1].E.Remove(e);
01884   delete v0; delete e;
01885  
01886   v0 = Vx2Sa;
01887   e = v0->e;
01888   v = e->v[1];
01889   G->P[0][1] = v;
01890   nhE = EdgeCycle( e, v, hE );
01891   el = hE[nhE-1];
01892   if( el->v[0] == v ) el->e[0] = e->e[1]; else el->e[1] = e->e[1];
01893   v->e = el;
01894   CutsY[1].V.Remove(v0); CutsY[1].E.Remove(e);
01895   delete v0; delete e;
01896  
01897   v0 = Vx1Sa;
01898   e = v0->e;
01899   v = e->v[0];
01900   G->P[0][0] = v;
01901   nhE = EdgeCycle( e, v, hE );
01902   el = hE[nhE-1];
01903   if( el->v[0] == v ) el->e[0] = e->e[0]; else el->e[1] = e->e[0];
01904   v->e = el;
01905   CutsY[1].V.Remove(v0); CutsY[1].E.Remove(e);
01906   delete v0; delete e;
01907 
01908   CutsY[1].R.DeleteElems();             // free cut records only
01909 
01910   G->E += CutsY[1].E;
01911   G->V += CutsY[1].V;
01912 
01913   CutsY[0].R.DeleteElems();
01914   CutsY[0].E.DeleteElems();
01915   CutsY[0].V.DeleteElems();
01916 
01917   CutsX[1].R.DeleteElems();
01918   CutsX[1].E.DeleteElems();
01919   CutsX[1].V.DeleteElems();
01920   
01921   CutsX[0].R.DeleteElems();
01922   CutsX[0].E.DeleteElems();
01923   CutsX[0].V.DeleteElems();
01924 }
01925 
01926 /* ---------------------------------------------------------------------------
01927   inits a graph info
01928 ---------------------------------------------------------------------------- */
01929 template
01930 void InitGraph( 
01931        const rational& X1, const rational& X2,
01932        const rational& Y1, const rational& Y2,
01933        GraphConvInfo<float> *Gconv, graph_edge_list *E, graph_vertex_list *V,
01934        GraphInfo *G );
01935 template
01936 void InitGraph( 
01937        const rational& X1, const rational& X2,
01938        const rational& Y1, const rational& Y2,
01939        GraphConvInfo<double> *Gconv, graph_edge_list *E, graph_vertex_list *V,
01940        GraphInfo *G );
01941 
01942 }

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