Back to index

wims  3.65+svn20090927
geometry.c
Go to the documentation of this file.
00001 
00002 /*** GEOMETRY.C ***/
00003 
00004 #include <math.h>
00005 #include "vdefs.h"
00006 
00007 float deltax, deltay ;
00008 int nedges, sqrt_nsites, nvertices ;
00009 Freelist efl ;
00010 
00011 void
00012 geominit(void)
00013     {
00014     freeinit(&efl, sizeof(Edge)) ;
00015     nvertices = nedges = 0 ;
00016     sqrt_nsites = sqrt(nsites+4) ;
00017     deltay = ymax - ymin ;
00018     deltax = xmax - xmin ;
00019     }
00020 
00021 Edge *
00022 bisect(Site * s1, Site * s2)
00023     {
00024     float dx, dy, adx, ady ;
00025     Edge * newedge ;
00026 
00027     newedge = (Edge *)getfree(&efl) ;
00028     newedge->reg[0] = s1 ;
00029     newedge->reg[1] = s2 ;
00030     ref(s1) ;
00031     ref(s2) ;
00032     newedge->ep[0] = newedge->ep[1] = (Site *)NULL ;
00033     dx = s2->coord.x - s1->coord.x ;
00034     dy = s2->coord.y - s1->coord.y ;
00035     adx = dx>0 ? dx : -dx ;
00036     ady = dy>0 ? dy : -dy ;
00037     newedge->c = s1->coord.x * dx + s1->coord.y * dy + (dx*dx +
00038 dy*dy) * 0.5 ;
00039     if (adx > ady)
00040         {
00041         newedge->a = 1.0 ;
00042         newedge->b = dy/dx ;
00043         newedge->c /= dx ;
00044         }
00045     else
00046         {
00047         newedge->b = 1.0 ;
00048         newedge->a = dx/dy ;
00049         newedge->c /= dy ;
00050         }
00051     newedge->edgenbr = nedges ;
00052     out_bisector(newedge) ;
00053     nedges++ ;
00054     return (newedge) ;
00055     }
00056 
00057 Site *
00058 intersect(Halfedge * el1, Halfedge * el2)
00059     {
00060     Edge * e1, * e2, * e ;
00061     Halfedge * el ;
00062     float d, xint, yint ;
00063     int right_of_site ;
00064     Site * v ;
00065 
00066     e1 = el1->ELedge ;
00067     e2 = el2->ELedge ;
00068     if ((e1 == (Edge*)NULL) || (e2 == (Edge*)NULL))
00069         {
00070         return ((Site *)NULL) ;
00071         }
00072     if (e1->reg[1] == e2->reg[1])
00073         {
00074         return ((Site *)NULL) ;
00075         }
00076     d = (e1->a * e2->b) - (e1->b * e2->a) ;
00077     if ((-1.0e-10 < d) && (d < 1.0e-10))
00078         {
00079         return ((Site *)NULL) ;
00080         }
00081     xint = (e1->c * e2->b - e2->c * e1->b) / d ;
00082     yint = (e2->c * e1->a - e1->c * e2->a) / d ;
00083     if ((e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
00084         (e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
00085         e1->reg[1]->coord.x < e2->reg[1]->coord.x))
00086         {
00087         el = el1 ;
00088         e = e1 ;
00089         }
00090     else
00091         {
00092         el = el2 ;
00093         e = e2 ;
00094         }
00095     right_of_site = (xint >= e->reg[1]->coord.x) ;
00096     if ((right_of_site && (el->ELpm == le)) ||
00097        (!right_of_site && (el->ELpm == re)))
00098         {
00099         return ((Site *)NULL) ;
00100         }
00101     v = (Site *)getfree(&sfl) ;
00102     v->refcnt = 0 ;
00103     v->coord.x = xint ;
00104     v->coord.y = yint ;
00105     return (v) ;
00106     }
00107 
00108 /*** returns 1 if p is to right of halfedge e ***/
00109 
00110 int
00111 right_of(Halfedge * el, Point * p)
00112     {
00113     Edge * e ;
00114     Site * topsite ;
00115     int right_of_site, above, fast ;
00116     float dxp, dyp, dxs, t1, t2, t3, yl ;
00117 
00118     e = el->ELedge ;
00119     topsite = e->reg[1] ;
00120     right_of_site = (p->x > topsite->coord.x) ;
00121     if (right_of_site && (el->ELpm == le))
00122         {
00123         return (1) ;
00124         }
00125     if(!right_of_site && (el->ELpm == re))
00126         {
00127         return (0) ;
00128         }
00129     if (e->a == 1.0)
00130         {
00131         dyp = p->y - topsite->coord.y ;
00132         dxp = p->x - topsite->coord.x ;
00133         fast = 0 ;
00134         if ((!right_of_site & (e->b < 0.0)) ||
00135          (right_of_site & (e->b >= 0.0)))
00136             {
00137             fast = above = (dyp >= e->b*dxp) ;
00138             }
00139         else
00140             {
00141             above = ((p->x + p->y * e->b) > (e->c)) ;
00142             if (e->b < 0.0)
00143                 {
00144                 above = !above ;
00145                 }
00146             if (!above)
00147                 {
00148                 fast = 1 ;
00149                 }
00150             }
00151         if (!fast)
00152             {
00153             dxs = topsite->coord.x - (e->reg[0])->coord.x ;
00154             above = (e->b * (dxp*dxp - dyp*dyp))
00155                     <
00156                     (dxs * dyp * (1.0 + 2.0 * dxp /
00157                     dxs + e->b * e->b)) ;
00158             if (e->b < 0.0)
00159                 {
00160                 above = !above ;
00161                 }
00162             }
00163         }
00164     else  /*** e->b == 1.0 ***/
00165         {
00166         yl = e->c - e->a * p->x ;
00167         t1 = p->y - yl ;
00168         t2 = p->x - topsite->coord.x ;
00169         t3 = yl - topsite->coord.y ;
00170         above = ((t1*t1) > ((t2 * t2) + (t3 * t3))) ;
00171         }
00172     return (el->ELpm == le ? above : !above) ;
00173     }
00174 
00175 void
00176 endpoint(Edge * e, int lr, Site * s)
00177     {
00178     e->ep[lr] = s ;
00179     ref(s) ;
00180     if (e->ep[re-lr] == (Site *)NULL)
00181         {
00182         return ;
00183         }
00184     out_ep(e) ;
00185     deref(e->reg[le]) ;
00186     deref(e->reg[re]) ;
00187     makefree((Freenode *)e, (Freelist *) &efl) ;
00188     }
00189 
00190 float
00191 dist(Site * s, Site * t)
00192     {
00193     float dx,dy ;
00194 
00195     dx = s->coord.x - t->coord.x ;
00196     dy = s->coord.y - t->coord.y ;
00197     return (sqrt(dx*dx + dy*dy)) ;
00198     }
00199 
00200 void
00201 makevertex(Site * v)
00202     {
00203     v->sitenbr = nvertices++ ;
00204     out_vertex(v) ;
00205     }
00206 
00207 void
00208 deref(Site * v)
00209     {
00210     if (--(v->refcnt) == 0 )
00211         {
00212         makefree((Freenode *)v, (Freelist *)&sfl) ;
00213         }
00214     }
00215 
00216 void
00217 ref(Site * v)
00218     {
00219     ++(v->refcnt) ;
00220     }
00221