source: trunk/src2/voronoi/geometry.c @ 666

Last change on this file since 666 was 666, checked in by mmc, 17 years ago

Added voronoi code, which is currently used by the nanovis server
to generate triangular meshes for point clouds.

File size: 4.7 KB
Line 
1
2/*** GEOMETRY.C ***/
3
4#include <math.h>
5#include "vdefs.h"
6
7float deltax, deltay ;
8int nedges, sqrt_nsites, nvertices ;
9Freelist efl ;
10
11void
12geominit(void)
13    {
14    freeinit(&efl, sizeof(Edge)) ;
15    nvertices = nedges = 0 ;
16    sqrt_nsites = sqrt(nsites+4) ;
17    deltay = ymax - ymin ;
18    deltax = xmax - xmin ;
19    }
20
21Edge *
22bisect(Site * s1, Site * s2)
23    {
24    float dx, dy, adx, ady ;
25    Edge * newedge ;
26
27    newedge = (Edge *)getfree(&efl) ;
28    newedge->reg[0] = s1 ;
29    newedge->reg[1] = s2 ;
30    ref(s1) ;
31    ref(s2) ;
32    newedge->ep[0] = newedge->ep[1] = (Site *)NULL ;
33    dx = s2->coord.x - s1->coord.x ;
34    dy = s2->coord.y - s1->coord.y ;
35    adx = dx>0 ? dx : -dx ;
36    ady = dy>0 ? dy : -dy ;
37    newedge->c = s1->coord.x * dx + s1->coord.y * dy + (dx*dx +
38dy*dy) * 0.5 ;
39    if (adx > ady)
40        {
41        newedge->a = 1.0 ;
42        newedge->b = dy/dx ;
43        newedge->c /= dx ;
44        }
45    else
46        {
47        newedge->b = 1.0 ;
48        newedge->a = dx/dy ;
49        newedge->c /= dy ;
50        }
51    newedge->edgenbr = nedges ;
52    out_bisector(newedge) ;
53    nedges++ ;
54    return (newedge) ;
55    }
56
57Site *
58intersect(Halfedge * el1, Halfedge * el2)
59    {
60    Edge * e1, * e2, * e ;
61    Halfedge * el ;
62    float d, xint, yint ;
63    int right_of_site ;
64    Site * v ;
65
66    e1 = el1->ELedge ;
67    e2 = el2->ELedge ;
68    if ((e1 == (Edge*)NULL) || (e2 == (Edge*)NULL))
69        {
70        return ((Site *)NULL) ;
71        }
72    if (e1->reg[1] == e2->reg[1])
73        {
74        return ((Site *)NULL) ;
75        }
76    d = (e1->a * e2->b) - (e1->b * e2->a) ;
77    if ((-1.0e-10 < d) && (d < 1.0e-10))
78        {
79        return ((Site *)NULL) ;
80        }
81    xint = (e1->c * e2->b - e2->c * e1->b) / d ;
82    yint = (e2->c * e1->a - e1->c * e2->a) / d ;
83    if ((e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
84        (e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
85        e1->reg[1]->coord.x < e2->reg[1]->coord.x))
86        {
87        el = el1 ;
88        e = e1 ;
89        }
90    else
91        {
92        el = el2 ;
93        e = e2 ;
94        }
95    right_of_site = (xint >= e->reg[1]->coord.x) ;
96    if ((right_of_site && (el->ELpm == le)) ||
97       (!right_of_site && (el->ELpm == re)))
98        {
99        return ((Site *)NULL) ;
100        }
101    v = (Site *)getfree(&sfl) ;
102    v->refcnt = 0 ;
103    v->coord.x = xint ;
104    v->coord.y = yint ;
105    return (v) ;
106    }
107
108/*** returns 1 if p is to right of halfedge e ***/
109
110int
111right_of(Halfedge * el, Point * p)
112    {
113    Edge * e ;
114    Site * topsite ;
115    int right_of_site, above, fast ;
116    float dxp, dyp, dxs, t1, t2, t3, yl ;
117
118    e = el->ELedge ;
119    topsite = e->reg[1] ;
120    right_of_site = (p->x > topsite->coord.x) ;
121    if (right_of_site && (el->ELpm == le))
122        {
123        return (1) ;
124        }
125    if(!right_of_site && (el->ELpm == re))
126        {
127        return (0) ;
128        }
129    if (e->a == 1.0)
130        {
131        dyp = p->y - topsite->coord.y ;
132        dxp = p->x - topsite->coord.x ;
133        fast = 0 ;
134        if ((!right_of_site & (e->b < 0.0)) ||
135         (right_of_site & (e->b >= 0.0)))
136            {
137            fast = above = (dyp >= e->b*dxp) ;
138            }
139        else
140            {
141            above = ((p->x + p->y * e->b) > (e->c)) ;
142            if (e->b < 0.0)
143                {
144                above = !above ;
145                }
146            if (!above)
147                {
148                fast = 1 ;
149                }
150            }
151        if (!fast)
152            {
153            dxs = topsite->coord.x - (e->reg[0])->coord.x ;
154            above = (e->b * (dxp*dxp - dyp*dyp))
155                    <
156                    (dxs * dyp * (1.0 + 2.0 * dxp /
157                    dxs + e->b * e->b)) ;
158            if (e->b < 0.0)
159                {
160                above = !above ;
161                }
162            }
163        }
164    else  /*** e->b == 1.0 ***/
165        {
166        yl = e->c - e->a * p->x ;
167        t1 = p->y - yl ;
168        t2 = p->x - topsite->coord.x ;
169        t3 = yl - topsite->coord.y ;
170        above = ((t1*t1) > ((t2 * t2) + (t3 * t3))) ;
171        }
172    return (el->ELpm == le ? above : !above) ;
173    }
174
175void
176endpoint(Edge * e, int lr, Site * s)
177    {
178    e->ep[lr] = s ;
179    ref(s) ;
180    if (e->ep[re-lr] == (Site *)NULL)
181        {
182        return ;
183        }
184    out_ep(e) ;
185    deref(e->reg[le]) ;
186    deref(e->reg[re]) ;
187    makefree((Freenode *)e, (Freelist *) &efl) ;
188    }
189
190float
191dist(Site * s, Site * t)
192    {
193    float dx,dy ;
194
195    dx = s->coord.x - t->coord.x ;
196    dy = s->coord.y - t->coord.y ;
197    return (sqrt(dx*dx + dy*dy)) ;
198    }
199
200void
201makevertex(Site * v)
202    {
203    v->sitenbr = nvertices++ ;
204    out_vertex(v) ;
205    }
206
207void
208deref(Site * v)
209    {
210    if (--(v->refcnt) == 0 )
211        {
212        makefree((Freenode *)v, (Freelist *)&sfl) ;
213        }
214    }
215
216void
217ref(Site * v)
218    {
219    ++(v->refcnt) ;
220    }
221
Note: See TracBrowser for help on using the repository browser.