Subversion Repositories wimsdev

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
10 reyssat 1
 
2
/*** GEOMETRY.C ***/
3
 
4
#include <math.h>
5
#include "vdefs.h"
6
 
7
float deltax, deltay ;
8
int nedges, sqrt_nsites, nvertices ;
9
Freelist efl ;
10
 
11
void
12
geominit(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
 
21
Edge *
22
bisect(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 +
38
dy*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
 
57
Site *
58
intersect(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
 
110
int
111
right_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
 
175
void
176
endpoint(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
 
190
float
191
dist(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
 
200
void
201
makevertex(Site * v)
202
    {
203
    v->sitenbr = nvertices++ ;
204
    out_vertex(v) ;
205
    }
206
 
207
void
208
deref(Site * v)
209
    {
210
    if (--(v->refcnt) == 0 )
211
        {
212
        makefree((Freenode *)v, (Freelist *)&sfl) ;
213
        }
214
    }
215
 
216
void
217
ref(Site * v)
218
    {
219
    ++(v->refcnt) ;
220
    }
221