OOFEM 3.0
Loading...
Searching...
No Matches
geotoolbox.C
Go to the documentation of this file.
1/*
2 *
3 * ##### ##### ###### ###### ### ###
4 * ## ## ## ## ## ## ## ### ##
5 * ## ## ## ## #### #### ## # ##
6 * ## ## ## ## ## ## ## ##
7 * ## ## ## ## ## ## ## ##
8 * ##### ##### ## ###### ## ##
9 *
10 *
11 * OOFEM : Object Oriented Finite Element Code
12 *
13 * Copyright (C) 1993 - 2025 Borek Patzak
14 *
15 *
16 *
17 * Czech Technical University, Faculty of Civil Engineering,
18 * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19 *
20 * This library is free software; you can redistribute it and/or
21 * modify it under the terms of the GNU Lesser General Public
22 * License as published by the Free Software Foundation; either
23 * version 2.1 of the License, or (at your option) any later version.
24 *
25 * This program is distributed in the hope that it will be useful,
26 * but WITHOUT ANY WARRANTY; without even the implied warranty of
27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28 * Lesser General Public License for more details.
29 *
30 * You should have received a copy of the GNU Lesser General Public
31 * License along with this library; if not, write to the Free Software
32 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33 */
34
35#include "geotoolbox.h"
36#include "mathfem.h"
37
38#ifdef __OOFEG
39 #include "oofeggraphiccontext.h"
40#endif
41
42namespace oofem {
43//#define GRAPH_DEBUG_PRINT
44#define GRAPH_LENGTH_FRAC 1.e-4
45
46int
47Polygon :: testPoint(double x, double y) const
48{
49 bool oddNODES = false;
50
51 Polygon :: PolygonEdgeIterator it(this);
52 Vertex p1, p2;
53 double x1, x2, y1, y2;
54 while ( it.giveNext(p1, p2) ) {
55 x1 = p1.coords(0);
56 y1 = p1.coords(1);
57 x2 = p2.coords(0);
58 y2 = p2.coords(1);
59
60 if ( ( ( y1 < y ) && ( y2 >= y ) ) ||
61 ( ( y2 < y ) && ( y1 >= y ) ) ) {
62 if ( x1 + ( y - y1 ) / ( y2 - y1 ) * ( x2 - x1 ) < x ) {
63 oddNODES = !oddNODES;
64 }
65 }
66 }
67
68 return oddNODES;
69}
70
71double
72Polygon :: computeVolume() const
73{
74 double area = 0.0;
75 double x1, x2, x3, y1, y2, y3;
76 Vertex p;
77 Polygon :: PolygonVertexIterator it(this);
78
79 it.giveNext(p);
80 x1 = p.coords(0);
81 y1 = p.coords(1);
82 it.giveNext(p);
83 x2 = p.coords(0);
84 y2 = p.coords(1);
85 while ( it.giveNext(p) ) {
86 x3 = p.coords(0);
87 y3 = p.coords(1);
88 area += ( x2 * y3 + x1 * y2 + y1 * x3 - x2 * y1 - x3 * y2 - x1 * y3 );
89 x2 = x3;
90 y2 = y3;
91 }
92
93 return 0.5 * fabs(area);
94}
95
96
97double
98Polygon :: pointDistance(double xp, double yp) const
99{
100 // computes distance of a given point from a closed polygon
101 // distance is signed, for polygon with positive orientation
102 // (anticlockwise vertex ordering) the positive distance is for a point
103 // located outside
104
105 Polygon :: PolygonEdgeIterator it(this);
106 Vertex p1, p2;
107 double x1, x2, y1, y2, d, l, tx, ty, t, nx, ny, dist = 0.0;
108 bool init = true;
109
110 while ( it.giveNext(p1, p2) ) {
111 x1 = p1.coords(0);
112 y1 = p1.coords(1);
113 x2 = p2.coords(0);
114 y2 = p2.coords(1);
115
116 // first check start vertex (end vertex checked by next edge)
117 d = sqrt( ( xp - x1 ) * ( xp - x1 ) + ( yp - y1 ) * ( yp - y1 ) );
118 if ( init ) {
119 dist = d;
120 init = false;
121 } else {
122 dist = min(dist, d);
123 }
124
125 // edge unit tangent
126 l = sqrt( ( x2 - x1 ) * ( x2 - x1 ) + ( y2 - y1 ) * ( y2 - y1 ) );
127 tx = ( x2 - x1 ) / l;
128 ty = ( y2 - y1 ) / l;
129 // projection of position vector (xp-x1,yp-y1) onto edge
130 t = ( xp - x1 ) * tx + ( yp - y1 ) * ty;
131 // test if inside
132 if ( ( t >= 0.0 ) && ( t <= l ) ) {
133 // compute distance from edge
134 nx = ty;
135 ny = -tx;
136 d = fabs( ( xp - x1 ) * nx + ( yp - y1 ) * ny );
137 dist = min(dist, d);
138 }
139 }
140
141 if ( this->testPoint(xp, yp) ) {
142 return -1.0 * dist;
143 } else {
144 return dist;
145 }
146}
147
148
149
150#ifdef __OOFEG
151GraphicObj *
152Polygon :: draw(oofegGraphicContext &gc, bool filled, int layer)
153{
154 GraphicObj *go;
155 LIST ggroup = make_list();
156
157 EASValsSetLayer(layer);
158 //EASValsSetColor(gc.getElementColor());
159
160 if ( filled ) {
161 int count = 0;
162 double xc = 0.0, yc = 0.0;
163 Vertex p;
164 WCRec r [ 3 ];
165 Polygon :: PolygonVertexIterator it(this);
166
167 while ( it.giveNext(p) ) {
168 xc += p.coords(0);
169 yc += p.coords(1);
170 count++;
171 }
172
173 xc /= count;
174 yc /= count;
175
176 EASValsSetFillStyle(FILL_SOLID);
177 r [ 0 ].x = xc;
178 r [ 0 ].y = yc;
179 r [ 0 ].z = 0.0;
180 Polygon :: PolygonEdgeIterator it2(this);
181 Vertex p1, p2;
182 while ( it2.giveNext(p1, p2) ) {
183 r [ 1 ].x = p1.coords(0);
184 r [ 1 ].y = p1.coords(1);
185 r [ 1 ].z = 0.0;
186 r [ 2 ].x = p2.coords(0);
187 r [ 2 ].y = p2.coords(1);
188 r [ 2 ].z = 0.0;
189 go = CreateTriangle3D(r);
190 add_to_tail(ggroup, go);
191 //EMAddGraphicsToModel(ESIModel(), go);
192 EGWithMaskChangeAttributes(COLOR_MASK | FILL_MASK | LAYER_MASK, go);
193 }
194 } else {
195 WCRec p [ 2 ];
196 Polygon :: PolygonEdgeIterator it(this);
197 Vertex p1, p2;
198 EASValsSetFillStyle(FILL_HOLLOW);
199 while ( it.giveNext(p1, p2) ) {
200 p [ 0 ].x = p1.coords(0);
201 p [ 0 ].y = p1.coords(1);
202 p [ 0 ].z = 0.0;
203 p [ 1 ].x = p2.coords(0);
204 p [ 1 ].y = p2.coords(1);
205 p [ 1 ].z = 0.0;
206
207 go = CreateLine3D(p);
208 add_to_tail(ggroup, go);
209 //EMAddGraphicsToModel(ESIModel(), go);
210 EGWithMaskChangeAttributes(COLOR_MASK | FILL_MASK | LAYER_MASK, go);
211 }
212 }
213
214 go = CreateGgroup(ggroup);
215 EMAddGraphicsToModel(ESIModel(), go);
216 //EGWithMaskChangeAttributes(FILL_MASK | LAYER_MASK, go);
217 return go;
218}
219#endif
220
221
222
223
224Graph :: ~Graph()
225{
226 this->clear();
227}
228
229
230void
231Graph :: clear()
232{
233 node *next, *auxs = s, *auxc = c;
234 if ( s ) {
235 do {
236 next = auxs->next;
237 delete auxs;
238 auxs = next;
239 } while ( auxs != s );
240 }
241
242 if ( c ) {
243 do {
244 next = auxc->next;
245 delete auxc;
246 auxc = next;
247 } while ( auxc != c );
248 }
249
250 s = c = NULL;
251}
252
253
254int
255Graph :: testIfCoincident(node *p1, node *p2, node *q1, node *q2, double *alpha1, double *alpha2)
256{
257 /* return codes:
258 * -1 ... lines are not parallel
259 * 0 .... lines do not coincide
260 *
261 * 1 .... lines are the same
262 * 2 .... lines are the same (difffferent orientation)
263 *
264 * 10 ... p1,p2 inside q1,q2 (alpha1 is parametric coord of p1 on q1-q2, alpha2 is param.coord of p2 on q1-q2)
265 * 20 ... q1,q2 inside p1,p2 (alpha1 is parametric coord of q1 on p1-p2, alpha2 is param.coord of q2 on p1-p2)
266 *
267 * 11 ... p1 inside q1,q2 (alpha1 is parametric coord of p1 on q1-q2)
268 * q1 inside p1,p1 (alpha2 is parametric coord of q1 on p1-p2)
269 * 12 ... p1 inside q1,q2 (alpha1 is parametric coord of p1 on q1-q2)
270 * q2 inside p1,p1 (alpha2 is parametric coord of q2 on p1-p2)
271 * 21 ... p2 inside q1,q2 (alpha1 is parametric coord of p2 on q1-q2)
272 * q1 inside p1,p1 (alpha2 is parametric coord of q1 on p1-p2)
273 * 22 ... p2 inside q1,q2 (alpha1 is parametric coord of p2 on q1-q2)
274 * q2 inside p1,p1 (alpha2 is parametric coord of q2 on p1-p2)
275 *
276 *
277 * 1102 ... common part is p1=q1, q2 (alpha1 is param coord of q2 on p1-p2)
278 * 1120 ... common part is p1=q1, p2 (alpha1 is param coord of p2 on q1-q2)
279 * 2102 ... common part is p2=q1, q2 (alpha1 is param coord of q2 on p1-p2)
280 * 2110 ... common part is p2=q1, p1 (alpha1 is param coord of p1 on q1-q2)
281 * 1201 ... common part is p1=q2, q1 (alpha1 is param coord of q1 on p1-p2)
282 * 1220 ... common part is p1=q2, p2 (alpha1 is param coord of p2 on q1-q2)
283 * 2201 ... common part is p2=q2, q1 (alpha1 is param coord of q1 on p1-p2)
284 * 2210 ... common part is p2=q2, p1 (alpha1 is param coord of p1 on q1-q2)
285 *
286 * 110 ... common part is only p1=q1, lines are parallel
287 * 210 ... common part is only p2=q1, lines are parallel
288 * 120 ... common part is only p1=q2, lines are parallel
289 * 220 ... common part is only p2=q2, lines are parallel
290 */
291
292
293 //par = ((p2->x - p1->x)*(q2->y - q1->y) -
294 // (p2->y - p1->y)*(q2->x - q1->x));
295
296 double par = ( ( p1->x - p2->x ) * ( q2->y - q1->y ) -
297 ( p1->y - p2->y ) * ( q2->x - q1->x ) );
298
299#ifdef GRAPH_DEBUG_PRINT
300 fprintf(stderr, "p1 [%e %e] p2 [%e %e] q1 [%e %e] q2 [%e %e]\n", p1->x, p1->y, p2->x, p2->y, q1->x, q1->y, q2->x, q2->y);
301 fprintf(stderr, "par %e\n", par);
302#endif
303 //if (!par) return 0; /* parallel lines */
304 if ( fabs(par) < GT_TEPS ) { /* parallel lines */
305 /* test if identical part exist */
306 double plength = dist(p1->x, p1->y, p2->x, p2->y);
307 double qlength = dist(q1->x, q1->y, q2->x, q2->y);
308 double ptx = ( p2->x - p1->x ) / plength;
309 double pty = ( p2->y - p1->y ) / plength;
310 double pnx = -pty;
311 double pny = ptx;
312 double d1 = ( pnx * ( q1->x - p1->x ) + pny * ( q1->y - p1->y ) );
313 double d2 = ( pnx * ( q2->x - p1->x ) + pny * ( q2->y - p1->y ) );
314#ifdef GRAPH_DEBUG_PRINT
315 fprintf(stderr, "dist (q-line, p-line)=%e\n", d1);
316#endif
317 if ( ( fabs( d1 / min(plength, qlength) ) >= GT_EPS ) && ( fabs( d2 / min(plength, qlength) ) >= GT_EPS ) && ( sgn(d1) * sgn(d2) > 0. ) ) {
318 return 0; // lines are parallel with nonzero distance
319 }
320
321 double qt1 = ( ptx * ( q1->x - p1->x ) + pty * ( q1->y - p1->y ) ) / plength;
322 double qt2 = ( ptx * ( q2->x - p1->x ) + pty * ( q2->y - p1->y ) ) / plength;
323#ifdef GRAPH_DEBUG_PRINT
324 fprintf(stderr, "param of q1 on pline t1=%e\n", qt1);
325 fprintf(stderr, "param of q2 on pline t2=%e\n", qt2);
326#endif
327 if ( ( qt1 < 0. ) && ( qt2 < 0. ) ) {
328 return 0; // lines do not have common part
329 }
330
331 if ( ( qt1 > 1. ) && ( qt2 > 1. ) ) {
332 return 0; // lines do not have common part
333 }
334
335 if ( ( fabs(qt1) < GT_EPS ) && ( fabs(1. - qt2) < GT_EPS ) ) { // lines are the same
336 return 1;
337 }
338
339 if ( ( fabs(qt2) < GT_EPS ) && ( fabs(1. - qt1) < GT_EPS ) ) { // lines are the same but orientation is different
340 return 2;
341 }
342
343 double qtx = ( q2->x - q1->x ) / qlength;
344 double qty = ( q2->y - q1->y ) / qlength;
345 double pt1 = ( qtx * ( p1->x - q1->x ) + qty * ( p1->y - q1->y ) ) / qlength;
346 double pt2 = ( qtx * ( p2->x - q1->x ) + qty * ( p2->y - q1->y ) ) / qlength;
347
348 // remains - partial overlap
349 if ( ( qt1 < -GT_EPS ) && ( qt2 > ( 1. + GT_EPS ) ) ) { // p1,p2 inside q1,q2
350 * alpha1 = pt1;
351 * alpha2 = pt2;
352 return 10;
353 } else if ( ( qt2 < -GT_EPS ) && ( qt1 > ( 1. + GT_EPS ) ) ) { // p1,p2 inside q1,q2
354 * alpha1 = pt1;
355 * alpha2 = pt2;
356 return 10;
357 } else if ( ( qt1 > GT_EPS ) && ( qt1 < ( 1. - GT_EPS ) ) && ( qt2 > GT_EPS ) && ( qt2 < ( 1. - GT_EPS ) ) ) { // q1,q2 inside p1,p2
358 * alpha1 = qt1;
359 * alpha2 = qt2;
360 return 20;
361 } else if ( fabs(qt1) < GT_EPS ) { // q1 concident with p1
362 if ( ( qt2 > GT_EPS ) && ( qt2 < ( 1. - GT_EPS ) ) ) { // q2 inside p1,p2
363 * alpha1 = qt2;
364 return 1102;
365 } else if ( ( pt2 > GT_EPS ) && ( pt2 < ( 1. - GT_EPS ) ) ) { // p2 inside q1,q2
366 * alpha1 = pt2;
367 return 1120;
368 } else { // only q1=p1 but no common part
369 return 110;
370 }
371 } else if ( fabs(1. - qt1) < GT_EPS ) { // q1 coincident with p2
372 if ( ( qt2 > GT_EPS ) && ( qt2 < ( 1. - GT_EPS ) ) ) { // q2 inside p1,p2
373 * alpha1 = qt2;
374 return 2102;
375 } else if ( ( pt1 > GT_EPS ) && ( pt1 < ( 1. - GT_EPS ) ) ) { // p1 inside q1,q2
376 * alpha1 = pt1;
377 return 2110;
378 } else { // only q1=p2 but no common part
379 return 210;
380 }
381 } else if ( fabs(qt2) < GT_EPS ) { // q2 coincident with p1
382 if ( ( qt1 > GT_EPS ) && ( qt1 < ( 1. - GT_EPS ) ) ) { // q1 inside p1,p2
383 * alpha1 = qt1;
384 return 1201;
385 } else if ( ( pt2 > GT_EPS ) && ( pt2 < ( 1. - GT_EPS ) ) ) { // p2 inside q1,q2
386 * alpha1 = pt2;
387 return 1220;
388 } else { // only q2=p1 but no common part
389 return 120;
390 }
391 } else if ( fabs(1. - qt2) < GT_EPS ) { // q2 coincident with p2
392 if ( ( qt1 > GT_EPS ) && ( qt1 < ( 1. - GT_EPS ) ) ) { // q1 inside p1,p2
393 * alpha1 = qt1;
394 return 2201;
395 } else if ( ( pt1 > GT_EPS ) && ( pt1 < ( 1. - GT_EPS ) ) ) { // p1 inside q1,q2
396 * alpha1 = pt1;
397 return 2210;
398 } else { // only q2=p2 but no common part
399 return 220;
400 }
401 } else if ( ( qt1 >= GT_EPS ) && ( qt1 <= ( 1. - GT_EPS ) ) ) { // q1 inside, q2 outside
402 double qlength = dist(q1->x, q1->y, q2->x, q2->y);
403 double qtx = ( q2->x - q1->x ) / qlength;
404 double qty = ( q2->y - q1->y ) / qlength;
405 double pt1 = ( qtx * ( p1->x - q1->x ) + qty * ( p1->y - q1->y ) ) / qlength;
406 double pt2 = ( qtx * ( p2->x - q1->x ) + qty * ( p2->y - q1->y ) ) / qlength;
407 if ( ( pt1 >= GT_EPS ) && ( pt1 <= ( 1. - GT_EPS ) ) ) { // p1 inside
408 * alpha1 = pt1;
409 * alpha2 = qt1;
410 return 11;
411 } else if ( ( pt2 >= GT_EPS ) && ( pt2 <= ( 1. - GT_EPS ) ) ) { // p2 inside
412 * alpha1 = pt2;
413 * alpha2 = qt1;
414 return 21;
415 } else if ( ( pt1 <= GT_EPS || ( 1. - pt1 ) <= GT_EPS ) && ( ( pt2 < 0. ) || ( pt2 > 1.0 ) ) ) { // p1 inside
416 * alpha1 = pt1;
417 * alpha2 = qt1;
418 return 11;
419 } else if ( ( pt2 <= GT_EPS || ( 1. - pt2 ) <= GT_EPS ) && ( ( pt1 < 0. ) || ( pt1 > 1.0 ) ) ) { // p2 inside
420 * alpha1 = pt2;
421 * alpha2 = qt1;
422 return 21;
423 } else {
424 fprintf(stderr, "testIfConcident: unresolved case q1 inside, q2 outside, pt1(t=%e), pt22(t=%e)\n", pt1, pt2);
425 return -1;
426 }
427 } else if ( ( qt2 >= GT_EPS ) && ( qt2 <= ( 1. - GT_EPS ) ) ) { // q2 inside, q1 outside
428 double qlength = dist(q1->x, q1->y, q2->x, q2->y);
429 double qtx = ( q2->x - q1->x ) / qlength;
430 double qty = ( q2->y - q1->y ) / qlength;
431 double pt1 = ( qtx * ( p1->x - q1->x ) + qty * ( p1->y - q1->y ) ) / qlength;
432 double pt2 = ( qtx * ( p2->x - q1->x ) + qty * ( p2->y - q1->y ) ) / qlength;
433 if ( ( pt1 >= GT_EPS ) && ( pt1 <= ( 1. - GT_EPS ) ) ) { // p1 inside
434 * alpha1 = pt1;
435 * alpha2 = qt2;
436 return 12;
437 } else if ( ( pt2 >= GT_EPS ) && ( pt2 <= ( 1. - GT_EPS ) ) ) { // p2 inside
438 * alpha1 = pt2;
439 * alpha2 = qt2;
440 return 22;
441 } else if ( ( pt1 <= GT_EPS || ( 1. - pt1 ) <= GT_EPS ) && ( ( pt2 < 0. ) || ( pt2 > 1.0 ) ) ) { // p1 inside
442 * alpha1 = pt1;
443 * alpha2 = qt2;
444 return 12;
445 } else if ( ( pt2 <= GT_EPS || ( 1. - pt2 ) <= GT_EPS ) && ( ( pt1 < 0. ) || ( pt1 > 1.0 ) ) ) { // p2 inside
446 * alpha1 = pt2;
447 * alpha2 = qt2;
448 return 22;
449 } else {
450 fprintf(stderr, "testIfConcident: unresolved case q2 inside, q1 outside, pt1(t=%e), pt2(t=%e)\n", pt1, pt2);
451 return -1;
452 }
453 } else {
454 fprintf(stderr, "testIfConcident: unresolved case, q1(t=%e), q2(t=%e)\n", qt1, qt2);
455 }
456 }
457
458 return -1;
459}
460
461
462int
463Graph :: testIfIntersect(node *p1, node *p2, node *q1, node *q2,
464 double *alpha_p, double *alpha_q, double *xint, double *yint)
465{
466 /* return codes:
467 * 0 - no intersection detected
468 * 1 - regular intersection found (alpha_q, alpha_q, xint, yint returned)
469 * -1 - unhandled case, internal error
470 *
471 * 11 - p1 q1 coincide, lines intersect there
472 * 22 - p2 q2 coincide, lines intersect there
473 * 12 - p1 q2 coincide, lines intersect there
474 * 21 - p2 q1 coincide, lines intersect there
475 *
476 *
477 * 110 - p1 being intersection on q1,q2 (alpha_q returned)
478 * 120 - p2 being intersection on q1,q2 (alpha_q returned)
479 * 101 - q1 being intersection on p1,p2 (alpha_p returned)
480 * 102 - q2 being intersection on p1,p2 (alpha_p returned)
481 */
482
483
484
485 double x, y, tp, tq, par;
486
487 par = ( ( p1->x - p2->x ) * ( q2->y - q1->y ) -
488 ( p1->y - p2->y ) * ( q2->x - q1->x ) );
489
490#ifdef GRAPH_DEBUG_PRINT
491 fprintf(stderr, "p1 [%e %e] p2 [%e %e] q1 [%e %e] q2 [%e %e]\n", p1->x, p1->y, p2->x, p2->y, q1->x, q1->y, q2->x, q2->y);
492 fprintf(stderr, "par %e\n", par);
493#endif
494 //if (!par) return 0; /* parallel lines */
495 if ( fabs(par) < GT_TEPS ) { /* parallel lines */
496 fprintf(stderr, "testIfIntersect: Parallel lines detected....\n");
497 fprintf(stderr, "p1 [%e %e] p2 [%e %e] q1 [%e %e] q2 [%e %e]\n", p1->x, p1->y, p2->x, p2->y, q1->x, q1->y, q2->x, q2->y);
498 fprintf(stderr, "par %e\n", par);
499 return -1;
500 } else {
501 //tp = ((q1->x - p1->x)*(q2->y - q1->y) - (q1->y - p1->y)*(q2->x - q1->x))/par;
502 //tq = ((p2->y - p1->y)*(q1->x - p1->x) - (p2->x - p1->x)*(q1->y - p1->y))/par;
503
504 tp = ( ( p1->x - q1->x ) * ( q2->y - q1->y ) - ( q2->x - q1->x ) * ( p1->y - q1->y ) ) / par;
505 tq = ( ( p1->x - p2->x ) * ( p1->y - q1->y ) - ( p1->x - q1->x ) * ( p1->y - p2->y ) ) / par;
506
507 //if(tp<0.0 || tp>1.0 || tq<0.0 || tq>1.0) return 0;
508 //if(tp<0.0 || tp>1.0 || tq<0.0 || tq>1.0) return 0;
509 if ( ( tp < -GT_EPS ) || ( tp > ( 1. + GT_EPS ) ) || ( tq < -GT_EPS ) || ( tq > ( 1 + GT_EPS ) ) ) {
510 return 0;
511 } else if ( ( tp >= GT_EPS ) && ( tp <= ( 1. - GT_EPS ) ) && ( tq >= GT_EPS ) && ( tq <= ( 1. - GT_EPS ) ) ) {
512 // regular intersection found
513 x = p1->x + tp * ( p2->x - p1->x );
514 y = p1->y + tp * ( p2->y - p1->y );
515
516 * alpha_p = dist(p1->x, p1->y, x, y) / dist(p1->x, p1->y, p2->x, p2->y);
517 * alpha_q = dist(q1->x, q1->y, x, y) / dist(q1->x, q1->y, q2->x, q2->y);
518 * xint = x;
519 * yint = y;
520
521 return 1;
522 } else if ( ( ( fabs(tp) <= GT_EPS ) || ( fabs(1. - tp) <= GT_EPS ) ) && ( ( tq >= GT_EPS ) && ( tq <= ( 1. - GT_EPS ) ) ) ) {
523 // one of p1 or p2 is intersection inside q1,q2 line
524 if ( fabs(tp) <= GT_EPS ) {
525 // p1 is on q1, q2
526 * alpha_q = dist(q1->x, q1->y, p1->x, p1->y) / dist(q1->x, q1->y, q2->x, q2->y);
527 return 110;
528 } else {
529 // p2 is on q1,q2
530 * alpha_q = dist(q1->x, q1->y, p2->x, p2->y) / dist(q1->x, q1->y, q2->x, q2->y);
531 return 120;
532 }
533 } else if ( ( ( fabs(tq) <= GT_EPS ) || ( fabs(1. - tq) <= GT_EPS ) ) && ( ( tp >= GT_EPS ) && ( tp <= ( 1. - GT_EPS ) ) ) ) {
534 // one of q1 or q2 is intersection inside p1,p2 line
535 if ( fabs(tq) <= GT_EPS ) {
536 // q1 is on q1, q2
537 * alpha_p = dist(p1->x, p1->y, q1->x, q1->y) / dist(p1->x, p1->y, p2->x, p2->y);
538 return 101;
539 } else {
540 // q2 is on q1,q2
541 * alpha_p = dist(p1->x, p1->y, q2->x, q2->y) / dist(p1->x, p1->y, p2->x, p2->y);
542 return 102;
543 }
544 } else if ( ( ( fabs(tq) <= GT_EPS ) || ( fabs(1. - tq) <= GT_EPS ) ) && ( ( fabs(tp) <= GT_EPS ) || ( fabs(1. - tp) <= GT_EPS ) ) ) {
545 // lines intersect at some of its vertices
546 if ( fabs(tq) <= GT_EPS ) {
547 if ( ( fabs(tp) <= GT_EPS ) ) {
548 // q1 p1 coincide, lines intersect there
549 return 11;
550 } else {
551 // q1 p2 coincide, lines intersect there
552 return 21;
553 }
554 } else {
555 if ( ( fabs(tp) <= GT_EPS ) ) {
556 // q2 p1 coincide, lines intersect there
557 return 12;
558 } else {
559 // q2 p2 coincide, lines intersect there
560 return 22;
561 }
562 }
563 } else {
564 fprintf(stderr, "testIfIntersect: unhandled case [tp %e, tq %e, par %e]\n", tp, tq, par);
565 return -1;
566 }
567
568 }
569}
570
571
572
573
574void
575Graph :: clip(Polygon &result, const Polygon &a, const Polygon &b)
576{
577 // create lists s,c from a and b
578 Polygon :: PolygonVertexIterator it( &a);
579 Vertex v;
580 node *nw;
581 int ret;
582 // s,c are input polygons
583 node *auxs = NULL, *auxc = NULL, *is, *ic;
584 node *p1, *p2, *q1, *q2;
585 double xi, yi;
586 double alpha_s, alpha_c, alpha1, alpha2;
587 node *crt;
588
589 result.clear();
590
591 it.init(& a);
592 while ( it.giveNext(v) ) {
593 nw = this->createNode(v.coords(0), v.coords(1), 0, 0, 0, 0, NS_Vertex, 0, 0, 0);
594 if ( c ) {
595 c->prev->next = nw;
596 nw->next = c;
597 nw->prev = c->prev;
598 c->prev = nw;
599 } else {
600 c = nw;
601 c->next = c->prev = c;
602 }
603 }
604
605 it.init(& b);
606 while ( it.giveNext(v) ) {
607 nw = this->createNode(v.coords(0), v.coords(1), 0, 0, 0, 0, NS_Vertex, 0, 0, 0);
608 if ( s ) {
609 s->prev->next = nw;
610 nw->next = s;
611 nw->prev = s->prev;
612 s->prev = nw;
613 } else {
614 s = nw;
615 s->next = s->prev = s;
616 }
617 }
618
619 if ( !s || !c ) {
620 return;
621 }
622
623 //auxs = last_node(s);
624 //createNode(s->x, s->y, 0, auxs, 0, 0, 0, 0, 0, 0.);
625 //auxc = last_node(c);
626 //createNode(c->x, c->y, 0, auxc, 0, 0, 0, 0, 0, 0.);
627
628 // restart:
629
630 auxs = s;
631 do {
632 //for(auxs = s; auxs == s; auxs = auxs->next)
633 if ( auxs->status != NS_Intersection ) {
634 p1 = auxs;
635 p2 = next_node(auxs->next);
636 if ( testCollapsedEdge(p1, p2) ) {
637 continue;
638 }
639
640 auxc = c;
641 do {
642 // for(auxc = c; auxc == c; auxc = auxc->next)
643 if ( auxc->status != NS_Intersection ) {
644 q1 = auxc;
645 q2 = next_node(auxc->next);
646 if ( testCollapsedEdge(q1, q2) ) {
647 continue;
648 }
649
650
651 ret = testIfCoincident(p1, p2, q1, q2, & alpha1, & alpha2);
652 /* return codes:
653 * 0 .... lines are parallel, but do not coincide
654 *
655 * 1 .... lines are the same
656 * 2 .... lines are the same (difffferent orientation)
657 *
658 * 10 ... p1,p2 inside q1,q2 (alpha1 is parametric coord of p1 on q1-q2, alpha2 is param.coord of p2 on q1-q2)
659 * 20 ... q1,q2 inside p1,p2 (alpha1 is parametric coord of q1 on p1-p2, alpha2 is param.coord of q2 on p1-p2)
660 *
661 * 11 ... p1 inside q1,q2 (alpha1 is parametric coord of p1 on q1-q2)
662 * q1 inside p1,p1 (alpha2 is parametric coord of q1 on p1-p2)
663 * 12 ... p1 inside q1,q2 (alpha1 is parametric coord of p1 on q1-q2)
664 * q2 inside p1,p1 (alpha2 is parametric coord of q2 on p1-p2)
665 * 21 ... p2 inside q1,q2 (alpha1 is parametric coord of p2 on q1-q2)
666 * q1 inside p1,p1 (alpha2 is parametric coord of q1 on p1-p2)
667 * 22 ... p2 inside q1,q2 (alpha1 is parametric coord of p2 on q1-q2)
668 * q2 inside p1,p1 (alpha2 is parametric coord of q2 on p1-p2)
669 *
670 * 1102 ... common part is p1=q1, q2 (alpha1 is param coord of q2 on p1-p2)
671 * 1120 ... common part is p1=q1, p2 (alpha1 is param coord of p2 on q1-q2)
672 * 2102 ... common part is p2=q1, q2 (alpha1 is param coord of q2 on p1-p2)
673 * 2110 ... common part is p2=q1, p1 (alpha1 is param coord of p1 on q1-q2)
674 * 1201 ... common part is p1=q2, q1 (alpha1 is param coord of q1 on p1-p2)
675 * 1220 ... common part is p1=q2, p2 (alpha1 is param coord of p2 on q1-q2)
676 * 2201 ... common part is p2=q2, q1 (alpha1 is param coord of q1 on p1-p2)
677 * 2210 ... common part is p2=q2, p1 (alpha1 is param coord of p1 on q1-q2)
678 *
679 * 110 ... common part is only p1=q1, lines are parallel
680 * 210 ... common part is only p2=q1, lines are parallel
681 * 120 ... common part is only p1=q2, lines are parallel
682 * 220 ... common part is only p2=q2, lines are parallel
683 */
684 if ( ret == 0 ) { /*linear are parallel, but do not coincide*/
685 } else if ( ret == 1 ) { /* lines are coincident */
686 /*
687 * q1->status=p1->status=q2->status=p2->status=NS_IntersectionVertex;
688 */
689 merge2vertex(p1, q1);
690 merge2vertex(p2, q2);
691 //vertex2IntersectionVertex(p1,q1,q2);
692 //vertex2IntersectionVertex(p2,q1,q2);
693 //vertex2IntersectionVertex(q1,p1,p2);
694 //vertex2IntersectionVertex(q2,p1,p2);
695
696 //p1->neighbor=q1; q1->neighbor=p1;
697 //p2->neighbor=q2; q2->neighbor=p2;
698 p1->entry = -1;
699 q1->entry = -1;
700 } else if ( ret == 2 ) { /*lines are the same (difffferent orientation)*/
701 //q1->status=p1->status=q2->status=p2->status=NS_IntersectionVertex;
702 merge2vertex(p1, q2);
703 merge2vertex(p2, q1);
704 //vertex2IntersectionVertex(p1,q1,q2);
705 //vertex2IntersectionVertex(p2,q1,q2);
706 //vertex2IntersectionVertex(q1,p1,p2);
707 //vertex2IntersectionVertex(q2,p1,p2);
708
709 //p1->neighbor=q2; q2->neighbor=p1;
710 //p2->neighbor=q1; q1->neighbor=p2;
711 p1->entry = -1;
712 q1->entry = -1;
713 } else if ( ret == 10 ) { /* 10 ... p1,p2 inside q1,q2 (alpha1 is param. of p1 on q1-q2, alpha2 is param of p2 on q1-q2)*/
714 //p1->status = NS_IntersectionVertex;
715 if ( vertex2IntersectionVertex(p1, q1, q2) ) {
716 is = createNode(p1->x, p1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
717 p1->neighbor = is;
718 is->neighbor = p1;
719 insert( is, auxc, next_node(auxc->next) );
720 if ( alpha1 < alpha2 ) {
721 is->entry = -1;
722 }
723
724 testNewIntersectionVertexEdgeCollapse(p1, q1, q2);
725 }
726
727 p1->entry = -1;
728
729 //p2->status = NS_IntersectionVertex;
730 if ( vertex2IntersectionVertex(p2, q1, q2) ) {
731 ic = createNode(p2->x, p2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha2);
732 p2->neighbor = ic;
733 ic->neighbor = p2;
734 insert( ic, auxc, next_node(auxc->next) );
735 if ( alpha1 > alpha2 ) {
736 ic->entry = -1;
737 }
738
739 testNewIntersectionVertexEdgeCollapse(p2, q1, q2);
740 }
741 } else if ( ret == 20 ) { /* 20...q1,q2 inside p1,p2 (alpha1 is parametric of q1 on p1-p2, alpha2 is param of q2 on p1-p2) */
742 //q1->status = NS_IntersectionVertex;
743 if ( vertex2IntersectionVertex(q1, p1, p2) ) {
744 is = createNode(q1->x, q1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
745 q1->neighbor = is;
746 is->neighbor = q1;
747 insert( is, auxs, next_node(auxs->next) );
748 if ( alpha1 < alpha2 ) {
749 is->entry = -1;
750 }
751
752 testNewIntersectionVertexEdgeCollapse(q1, p1, p2);
753 }
754
755 q1->entry = -1;
756
757 //q2->status = NS_IntersectionVertex;
758 if ( vertex2IntersectionVertex(q2, p1, p2) ) {
759 ic = createNode(q2->x, q2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha2);
760 q2->neighbor = ic;
761 ic->neighbor = q2;
762 insert( ic, auxs, next_node(auxs->next) );
763 if ( alpha1 > alpha2 ) {
764 ic->entry = -1;
765 }
766
767 testNewIntersectionVertexEdgeCollapse(q2, p1, p2);
768 }
769 } else if ( ret == 11 ) {
770 /*
771 * 11 ... p1 inside q1,q2 (alpha1 is parametric coord of p1 on q1-q2)
772 * q1 inside p1,p1 (alpha2 is parametric coord of q1 on p1-p2)
773 */
774 //p1->status=NS_IntersectionVertex;
775 if ( vertex2IntersectionVertex(p1, q1, q2) ) {
776 is = createNode(p1->x, p1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
777 p1->neighbor = is;
778 is->neighbor = p1;
779 insert( is, auxc, next_node(auxc->next) );
780 testNewIntersectionVertexEdgeCollapse(p1, q1, q2);
781 }
782
783 p1->entry = -1;
784 //q1->status = NS_IntersectionVertex;
785 if ( vertex2IntersectionVertex(q1, p1, p2) ) {
786 is = createNode(q1->x, q1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha2);
787 q1->neighbor = is;
788 is->neighbor = q1;
789 insert( is, auxs, next_node(auxs->next) );
790 testNewIntersectionVertexEdgeCollapse(q1, p1, p2);
791 }
792
793 q1->entry = -1;
794 } else if ( ret == 12 ) {
795 /*
796 * 12 ... p1 inside q1,q2 (alpha1 is parametric coord of p1 on q1-q2)
797 * q2 inside p1,p1 (alpha2 is parametric coord of q2 on p1-p2)
798 */
799 //p1->status=NS_IntersectionVertex;
800 if ( vertex2IntersectionVertex(p1, q1, q2) ) {
801 is = createNode(p1->x, p1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
802 p1->neighbor = is;
803 is->neighbor = p1;
804 insert( is, auxc, next_node(auxc->next) );
805 testNewIntersectionVertexEdgeCollapse(p1, q1, q2);
806 }
807
808 p1->entry = p1->neighbor->entry = -1;
809
810 //q2->status = NS_IntersectionVertex;
811 if ( vertex2IntersectionVertex(q2, p1, p2) ) {
812 is = createNode(q2->x, q2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha2);
813 q2->neighbor = is;
814 is->neighbor = q2;
815 insert( is, auxs, next_node(auxs->next) );
816 testNewIntersectionVertexEdgeCollapse(q2, p1, p2);
817 }
818 } else if ( ret == 21 ) {
819 /*
820 * 21 ... p2 inside q1,q2 (alpha1 is parametric coord of p2 on q1-q2)
821 * q1 inside p1,p1 (alpha2 is parametric coord of q1 on p1-p2)
822 */
823 //p2->status=NS_IntersectionVertex;
824 if ( vertex2IntersectionVertex(p2, q1, q2) ) {
825 is = createNode(p2->x, p2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
826 p2->neighbor = is;
827 is->neighbor = p2;
828 insert( is, auxc, next_node(auxc->next) );
829 testNewIntersectionVertexEdgeCollapse(p2, q1, q2);
830 }
831
832 //q1->status = NS_IntersectionVertex;
833 if ( vertex2IntersectionVertex(q1, p1, p2) ) {
834 is = createNode(q1->x, q1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha2);
835 q1->neighbor = is;
836 is->neighbor = q1;
837 insert( is, auxs, next_node(auxs->next) );
838 testNewIntersectionVertexEdgeCollapse(q1, p1, p2);
839 }
840
841 q1->entry = q1->neighbor->entry = -1;
842 } else if ( ret == 22 ) {
843 /*
844 * 22 ... p2 inside q1,q2 (alpha1 is parametric coord of p2 on q1-q2)
845 * q2 inside p1,p1 (alpha2 is parametric coord of q2 on p1-p2)
846 */
847 //p2->status=NS_IntersectionVertex;
848 if ( vertex2IntersectionVertex(p2, q1, q2) ) {
849 is = createNode(p2->x, p2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
850 p2->neighbor = is;
851 is->neighbor = p2;
852 insert( is, auxc, next_node(auxc->next) );
853 testNewIntersectionVertexEdgeCollapse(p2, q1, q2);
854 }
855
856 p2->neighbor->entry = -1;
857 //q2->status = NS_IntersectionVertex;
858 if ( vertex2IntersectionVertex(q2, p1, p2) ) {
859 is = createNode(q2->x, q2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha2);
860 q2->neighbor = is;
861 is->neighbor = q2;
862 insert( is, auxs, next_node(auxs->next) );
863 testNewIntersectionVertexEdgeCollapse(q2, p1, p2);
864 }
865
866 q2->neighbor->entry = -1;
867 } else if ( ret == 1102 ) { /*1102 ... common part is p1=q1, q2 (alpha1 is param coord of q2 on p1-p2)*/
868 //p1->status = q1->status = NS_IntersectionVertex;
869 merge2vertex(p1, q1);
870 //vertex2IntersectionVertex(p1,q1,q2);
871 //vertex2IntersectionVertex(q1,p1,p2);
872 //p1->neighbor = q1; q1->neighbor = p1;
873 p1->entry = q1->entry = -1;
874
875 if ( vertex2IntersectionVertex(q2, p1, p2) ) {
876 //q2->status = NS_IntersectionVertex;
877 is = createNode(q2->x, q2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
878 q2->neighbor = is;
879 is->neighbor = q2;
880 insert( is, auxs, next_node(auxs->next) );
881 testNewIntersectionVertexEdgeCollapse(q2, p1, p2);
882 }
883 } else if ( ret == 1120 ) { /*1120 ... common part is p1=q1, p2 (alpha1 is param coord of p2 on q1-q2)*/
884 //p1->status = q1->status = NS_IntersectionVertex;
885 merge2vertex(p1, q1);
886 //vertex2IntersectionVertex(p1,q1,q2);
887 //vertex2IntersectionVertex(q1,p1,p2);
888 //p1->neighbor = q1; q1->neighbor = p1;
889 p1->entry = q1->entry = -1;
890
891 //p2->status = NS_IntersectionVertex;
892 if ( vertex2IntersectionVertex(p2, q1, q2) ) {
893 is = createNode(p2->x, p2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
894 p2->neighbor = is;
895 is->neighbor = p2;
896 insert( is, auxc, next_node(auxc->next) );
897 testNewIntersectionVertexEdgeCollapse(p2, q1, q2);
898 }
899 } else if ( ret == 2102 ) { /*2102 ... common part is p2=q1, q2 (alpha1 is param coord of q2 on p1-p2)*/
900 //p2->status = q1->status = NS_IntersectionVertex;
901 merge2vertex(p2, q1);
902 //vertex2IntersectionVertex(p2,q1,q2);
903 //vertex2IntersectionVertex(q1,p1,p2);
904 //p2->neighbor = q1; q1->neighbor = p2;
905 q1->entry = -1;
906
907 //q2->status = NS_IntersectionVertex;
908 if ( vertex2IntersectionVertex(q2, p1, p2) ) {
909 is = createNode(q2->x, q2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
910 q2->neighbor = is;
911 is->neighbor = q2;
912 insert( is, auxs, next_node(auxs->next) );
913 testNewIntersectionVertexEdgeCollapse(q2, p1, p2);
914 }
915
916 q2->neighbor->entry = -1;
917 } else if ( ret == 2110 ) { /*2110 ... common part is p2=q1, p1 (alpha1 is param coord of p1 on q1-q2*/
918 //p2->status = q1->status = NS_IntersectionVertex;
919 merge2vertex(p2, q1);
920 //vertex2IntersectionVertex(p2,q1,q2);
921 //vertex2IntersectionVertex(q1,p1,p2);
922 //p2->neighbor = q1; q1->neighbor = p2;
923 q1->entry = -1;
924 p1->entry = -1;
925
926 //p1->status = NS_IntersectionVertex;
927 if ( vertex2IntersectionVertex(p1, q1, q2) ) {
928 is = createNode(p1->x, p1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
929 p1->neighbor = is;
930 is->neighbor = p1;
931 insert( is, auxc, next_node(auxc->next) );
932 testNewIntersectionVertexEdgeCollapse(p1, q1, q2);
933 }
934 } else if ( ret == 1201 ) { /*1201 ... common part is p1=q2, q1 (alpha1 is param coord of q1 on p1-p2)*/
935 //p1->status = q2->status = NS_IntersectionVertex;
936 merge2vertex(p1, q2);
937 //vertex2IntersectionVertex(p1,q1,q2);
938 //vertex2IntersectionVertex(q2,p1,p2);
939 //p1->neighbor = q2; q2->neighbor = p1;
940 p1->entry = q1->entry = -1;
941
942 //q1->status = NS_IntersectionVertex;
943 if ( vertex2IntersectionVertex(q1, p1, p2) ) {
944 is = createNode(q1->x, q1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
945 q1->neighbor = is;
946 is->neighbor = q1;
947 insert( is, auxs, next_node(auxs->next) );
948 testNewIntersectionVertexEdgeCollapse(q1, p1, p2);
949 }
950 } else if ( ret == 1220 ) { /*1220 ... common part is p1=q2, p2 (alpha1 is param coord of p2 on q1-q2)*/
951 //p1->status = q2->status = NS_IntersectionVertex;
952 merge2vertex(p1, q2);
953 //vertex2IntersectionVertex(p1,q1,q2);
954 //vertex2IntersectionVertex(q2,p1,p2);
955 //p1->neighbor = q2; q2->neighbor = p1;
956 p1->entry = -1;
957
958 //p2->status = NS_IntersectionVertex;
959 if ( vertex2IntersectionVertex(p2, q1, q2) ) {
960 is = createNode(p2->x, p2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
961 p2->neighbor = is;
962 is->neighbor = p2;
963 insert( is, auxc, next_node(auxc->next) );
964 testNewIntersectionVertexEdgeCollapse(p2, q1, q2);
965 }
966
967 p2->neighbor->entry = -1;
968 } else if ( ret == 2201 ) { /*2201 ... common part is p2=q2, q1 (alpha1 is param coord of q1 on p1-p2)*/
969 //p2->status = q2->status = NS_IntersectionVertex;
970 merge2vertex(p2, q2);
971 //vertex2IntersectionVertex(p2,q1,q2);
972 //vertex2IntersectionVertex(q2,p1,p2);
973 //p2->neighbor = q2; q2->neighbor = p2;
974 q1->entry = -1;
975
976 //q1->status = NS_IntersectionVertex;
977 if ( vertex2IntersectionVertex(q1, p1, p2) ) {
978 is = createNode(q1->x, q1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
979 q1->neighbor = is;
980 is->neighbor = q1;
981 insert( is, auxs, next_node(auxs->next) );
982 testNewIntersectionVertexEdgeCollapse(q1, p1, p2);
983 }
984
985 q1->neighbor->entry = -1;
986 } else if ( ret == 2210 ) { /*2210 ... common part is p2=q2, p1 (alpha1 is param coord of p1 on q1-q2)*/
987 //p2->status = q2->status = NS_IntersectionVertex;
988 merge2vertex(p2, q2);
989 //vertex2IntersectionVertex(p2,q1,q2);
990 //vertex2IntersectionVertex(q2,p1,p2);
991 //p2->neighbor = q2; q2->neighbor = p2;
992 p1->entry = -1;
993
994 //p1->status = NS_IntersectionVertex;
995 if ( vertex2IntersectionVertex(p1, q1, q2) ) {
996 is = createNode(p1->x, p1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha1);
997 p1->neighbor = is;
998 is->neighbor = p1;
999 insert( is, auxc, next_node(auxc->next) );
1000 testNewIntersectionVertexEdgeCollapse(p1, q1, q2);
1001 }
1002
1003 p1->neighbor->entry = -1;
1004 } else if ( ret == 110 ) { /*common part is only p1=q1, lines are parallel*/
1005 //p1->status=q1->status = NS_IntersectionVertex;
1006 merge2vertex(p1, q1);
1007 //vertex2IntersectionVertex(p1,q1,q2);
1008 //vertex2IntersectionVertex(q1,p1,p2);
1009 //p1->neighbor = q1; q1->neighbor = p1;
1010 } else if ( ret == 210 ) { /*common part is only p2=q1, lines are parallel*/
1011 //p2->status=q1->status = NS_IntersectionVertex;
1012 merge2vertex(p2, q1);
1013 //vertex2IntersectionVertex(p2,q1,q2);
1014 //vertex2IntersectionVertex(q1,p1,p2);
1015 //p2->neighbor = q1; q1->neighbor = p2;
1016 } else if ( ret == 120 ) { /*common part is only p1=q2, lines are parallel*/
1017 //p1->status=q2->status = NS_IntersectionVertex;
1018 merge2vertex(p1, q2);
1019 //vertex2IntersectionVertex(p1,q1,q2);
1020 //vertex2IntersectionVertex(q2,p1,p2);
1021 //p1->neighbor = q2; q2->neighbor = p1;
1022 } else if ( ret == 220 ) { /*common part is only p2=q2, lines are parallel*/
1023 //p2->status=q2->status = NS_IntersectionVertex;
1024 merge2vertex(p2, q2);
1025 //vertex2IntersectionVertex(p2,q1,q2);
1026 //vertex2IntersectionVertex(q2,p1,p2);
1027 //p2->neighbor = q2; q2->neighbor = p2;
1028 } else if ( ret == -1 ) {
1029 ret = testIfIntersect(auxs, next_node(auxs->next), auxc, next_node(auxc->next),
1030 & alpha_s, & alpha_c, & xi, & yi);
1031 /* return codes:
1032 * 0 - no intersection detected
1033 * 1 - regular intersection found (alpha_q, alpha_q, xint, yint returned)
1034 * -1 - unhandled case, internal error
1035 *
1036 * 11 - p1 q1 coincide, lines intersect there
1037 * 22 - p2 q2 coincide, lines intersect there
1038 * 12 - p1 q2 coincide, lines intersect there
1039 * 21 - p2 q1 coincide, lines intersect there
1040 *
1041 *
1042 * 110 - p1 being intersection on q1,q2 (alpha_q returned)
1043 * 120 - p2 being intersection on q1,q2 (alpha_q returned)
1044 * 101 - q1 being intersection on p1,p2 (alpha_p returned)
1045 * 102 - q2 being intersection on p1,p2 (alpha_p returned)
1046 */
1047 if ( ret == 0 ) { // no intersection
1048 } else if ( ret == 1 ) {
1049 bool skip = false;
1050 // 1 .... regular intersection
1051 if ( ( p1->status == NS_IntersectionVertex ) || ( p2->status == NS_IntersectionVertex ) ||
1052 ( q1->status == NS_IntersectionVertex ) || ( q2->status == NS_IntersectionVertex ) ) {
1053 bool _p1 = false, _p2 = false, _q1 = false, _q2 = false;
1054 int c1 = 0, c2 = 0;
1055 // check if IV not associated to same line
1056 if ( p1->status == NS_IntersectionVertex ) {
1057 if ( ( _p1 = belongs(p1, q1, q2) ) ) {
1058 c1++;
1059 }
1060 }
1061
1062 if ( p2->status == NS_IntersectionVertex ) {
1063 if ( ( _p2 = belongs(p2, q1, q2) ) ) {
1064 c1++;
1065 }
1066 }
1067
1068 if ( q1->status == NS_IntersectionVertex ) {
1069 if ( ( _q1 = belongs(q1, p1, p2) ) ) {
1070 c2++;
1071 }
1072 }
1073
1074 if ( q2->status == NS_IntersectionVertex ) {
1075 if ( ( _q2 = belongs(q2, p1, p2) ) ) {
1076 c2++;
1077 }
1078 }
1079
1080 if ( ( c1 == 1 ) || ( c2 == 1 ) ) {
1081 skip = true;
1082 } else if ( ( c1 > 1 ) || ( c2 > 1 ) ) {
1083 // q1 q2 is already boundary edge on p1 p2 (or vice versa)
1084 skip = true;
1085 //this->printYourself();
1086 //printf ("c1 %d, c2 %d, p1x %e p1y %e, q1x %e, q1y %e\n", c1,c2,p1->x,p1->y,q1->x, q1->y);
1087 //THROW_GT_EXCEPTIONM ("Graph::clip unhandled case");
1088 }
1089 }
1090
1091 if ( !skip ) {
1092 is = createNode(xi, yi, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha_s);
1093 ic = createNode(xi, yi, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha_c);
1094 is->neighbor = ic;
1095 ic->neighbor = is;
1096 insert( is, auxs, next_node(auxs->next) );
1097 insert( ic, auxc, next_node(auxc->next) );
1098 }
1099 } else if ( ret == 11 ) { /*11 - p1 q1 coincide, lines intersect there*/
1100 //p1->status=q1->status=NS_IntersectionVertex;
1101 merge2vertex(p1, q1);
1102 //vertex2IntersectionVertex(p1,q1,q2);
1103 //vertex2IntersectionVertex(q1,p1,p2);
1104 //p1->neighbor=q1; q1->neighbor=p1;
1105 } else if ( ret == 22 ) { /*22 - p2 q2 coincide, lines intersect there*/
1106 //p2->status=q2->status=NS_IntersectionVertex;
1107 merge2vertex(p2, q2);
1108 //vertex2IntersectionVertex(p2,q1,q2);
1109 //vertex2IntersectionVertex(q2,p1,p2);
1110 //p2->neighbor=q2; q2->neighbor=p2;
1111 } else if ( ret == 12 ) { /*12 - p1 q2 coincide, lines intersect there*/
1112 //p1->status=q2->status=NS_IntersectionVertex;
1113 merge2vertex(p1, q2);
1114 //vertex2IntersectionVertex(p1,q1,q2);
1115 //vertex2IntersectionVertex(q2,p1,p2);
1116 //p1->neighbor=q2; q2->neighbor=p1;
1117 } else if ( ret == 21 ) { /*21 - p2 q1 coincide, lines intersect there*/
1118 //p2->status=q1->status=NS_IntersectionVertex;
1119 merge2vertex(p2, q1);
1120 //vertex2IntersectionVertex(p2,q1,q2);
1121 //vertex2IntersectionVertex(q1,p1,p2);
1122 //p2->neighbor=q1; q1->neighbor=p2;
1123 } else if ( ret == 110 ) { /*110 - p1 being intersection on q1,q2 (alpha_q returned)*/
1124 //p1->status=NS_IntersectionVertex;
1125 if ( vertex2IntersectionVertex(p1, q1, q2) ) {
1126 is = createNode(p1->x, p1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha_c);
1127 p1->neighbor = is;
1128 is->neighbor = p1;
1129 insert( is, auxc, next_node(auxc->next) );
1130 testNewIntersectionVertexEdgeCollapse(p1, q1, q2);
1131 }
1132 } else if ( ret == 120 ) { /*120 - p2 being intersection on q1,q2 (alpha_q returned)*/
1133 //p2->status=NS_IntersectionVertex;
1134 if ( vertex2IntersectionVertex(p2, q1, q2) ) {
1135 is = createNode(p2->x, p2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha_c);
1136 p2->neighbor = is;
1137 is->neighbor = p2;
1138 insert( is, auxc, next_node(auxc->next) );
1139 testNewIntersectionVertexEdgeCollapse(p2, q1, q2);
1140 }
1141 } else if ( ret == 101 ) { /* 101 - q1 being intersection on p1,p2 (alpha_p returned)*/
1142 //q1->status=NS_IntersectionVertex;
1143 if ( vertex2IntersectionVertex(q1, p1, p2) ) {
1144 is = createNode(q1->x, q1->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha_s);
1145 q1->neighbor = is;
1146 is->neighbor = q1;
1147 insert( is, auxs, next_node(auxs->next) );
1148 testNewIntersectionVertexEdgeCollapse(q1, p1, p2);
1149 }
1150 } else if ( ret == 102 ) { /* 102 - q2 being intersection on p1,p2 (alpha_p returned)*/
1151 //q2->status=NS_IntersectionVertex;
1152 if ( vertex2IntersectionVertex(q2, p1, p2) ) {
1153 is = createNode(q2->x, q2->y, 0, 0, 0, 0, NS_Intersection, 0, 0, alpha_s);
1154 q2->neighbor = is;
1155 is->neighbor = q2;
1156 insert( is, auxs, next_node(auxs->next) );
1157 testNewIntersectionVertexEdgeCollapse(q2, p1, p2);
1158 }
1159 } else {
1160 char buff [ 100 ];
1161 sprintf(buff, "Graph::clip: unhandled return code (%d) from testIfIntersect service", ret);
1162 THROW_GT_EXCEPTIONM(buff);
1163 }
1164 } else {
1165 char buff [ 100 ];
1166 sprintf(buff, "Graph::clip: unhandled return code (%d) from testIfCoincident service", ret);
1167 THROW_GT_EXCEPTIONM(buff);
1168 }
1169 }
1170 } while ( ( auxc = auxc->next ) != c );
1171 }
1172 } while ( ( auxs = auxs->next ) != s );
1173
1174#ifdef GRAPH_DEBUG_PRINT
1175 /* print graph */
1176 printYourself();
1177 /* end of debug printing */
1178#endif
1179
1180 // ================END OF WORK ==========================================
1181
1182 // assign status (outer/inner) to graph intersection vertices
1183 int cnt = 0;
1184 bool identical = true;
1185 double xx, yy;
1186 //insidea = e = testPoint(c, s->x, s->y);
1187 auxs = s;
1188 do {
1189 if ( auxs->entry != -1 ) {
1190 identical = false;
1191 }
1192
1193 if ( auxs->status == NS_Intersection ) {
1194 if ( auxs->entry != -1 ) { // boundary
1195 double xn, yn;
1196 if ( auxs->next->status == NS_IntersectionVertex ) {
1197 xn = 0.5 * ( auxs->next->x + auxs->next->neighbor->x );
1198 yn = 0.5 * ( auxs->next->y + auxs->next->neighbor->y );
1199 } else {
1200 xn = auxs->next->x;
1201 yn = auxs->next->y;
1202 }
1203
1204 xx = ( auxs->x + xn ) * 0.5;
1205 yy = ( auxs->y + yn ) * 0.5;
1206 auxs->entry = testPoint(c, xx, yy);
1207 cnt++;
1208 }
1209 } else if ( auxs->status == NS_IntersectionVertex ) {
1210 if ( auxs->entry != -1 ) { //boundary
1211 double xn, yn;
1212 if ( auxs->next->status == NS_IntersectionVertex ) {
1213 xn = 0.5 * ( auxs->next->x + auxs->next->neighbor->x );
1214 yn = 0.5 * ( auxs->next->y + auxs->next->neighbor->y );
1215 } else {
1216 xn = auxs->next->x;
1217 yn = auxs->next->y;
1218 }
1219
1220 xx = ( 0.5 * ( auxs->x + auxs->neighbor->x ) + xn ) * 0.5;
1221 yy = ( 0.5 * ( auxs->y + auxs->neighbor->y ) + yn ) * 0.5;
1222 auxs->entry = testPoint(c, xx, yy);
1223 cnt++;
1224 }
1225 }
1226
1227 auxs = auxs->next;
1228 } while ( auxs != s );
1229
1230
1231 //insideb = e = testPoint(s, c->x, c->y);
1232 do {
1233 if ( auxc->status == NS_Intersection ) {
1234 if ( auxc->entry != -1 ) { // boundary
1235 double xn, yn;
1236 if ( auxc->next->status == NS_IntersectionVertex ) {
1237 xn = 0.5 * ( auxc->next->x + auxc->next->neighbor->x );
1238 yn = 0.5 * ( auxc->next->y + auxc->next->neighbor->y );
1239 } else {
1240 xn = auxc->next->x;
1241 yn = auxc->next->y;
1242 }
1243
1244 xx = ( auxc->x + xn ) * 0.5;
1245 yy = ( auxc->y + yn ) * 0.5;
1246 auxc->entry = testPoint(s, xx, yy);
1247 }
1248 } else if ( auxc->status == NS_IntersectionVertex ) {
1249 if ( auxc->entry != -1 ) { //boundary
1250 double xn, yn;
1251 if ( auxc->next->status == NS_IntersectionVertex ) {
1252 xn = 0.5 * ( auxc->next->x + auxc->next->neighbor->x );
1253 yn = 0.5 * ( auxc->next->y + auxc->next->neighbor->y );
1254 } else {
1255 xn = auxc->next->x;
1256 yn = auxc->next->y;
1257 }
1258
1259 xx = ( 0.5 * ( auxc->x + auxc->neighbor->x ) + xn ) * 0.5;
1260 yy = ( 0.5 * ( auxc->y + auxc->neighbor->y ) + yn ) * 0.5;
1261 auxc->entry = testPoint(s, xx, yy);
1262 }
1263 }
1264
1265 auxc = auxc->next;
1266 } while ( auxc != c );
1267
1268
1269 if ( identical ) {
1270 auxs = s;
1271 do {
1272 v.setCoords(auxs->x, auxs->y);
1273 result.addVertex(v);
1274 auxs = auxs->next;
1275 } while ( auxs != s );
1276 } else if ( cnt == 0 ) {
1277 xx = ( s->x + s->next->x ) * 0.5;
1278 yy = ( s->y + s->next->y ) * 0.5;
1279 if ( testPoint(c, xx, yy) ) {
1280 // s fully inside c
1281 auxs = s;
1282 do {
1283 v.setCoords(auxs->x, auxs->y);
1284 result.addVertex(v);
1285 auxs = auxs->next;
1286 } while ( auxs != s );
1287 } else if ( testPoint( s, 0.5 * ( c->x + c->next->x ), 0.5 * ( c->y + c->next->y ) ) ) {
1288 // c fully in s
1289 auxc = c;
1290 do {
1291 v.setCoords(auxc->x, auxc->y);
1292 result.addVertex(v);
1293 auxc = auxc->next;
1294 } while ( auxc != c );
1295 }
1296 } else {
1297 /* loop over all entry points (non visited intersection) */
1298 //while ((crt = first(s)) != s)
1299 while ( ( crt = first(s) ) ) {
1300 //old = 0;
1301 /* loop over not yet visited intersections */
1302 for ( ; !crt->visited; crt = crt->neighbor ) {
1303 if ( crt->entry == 0 ) {
1304 crt->visited = 1;
1305 crt = crt->neighbor; // if not entry intersection use neighbor
1306 //if (crt->entry ==0 || crt->visited ) break;
1307 if ( crt->entry == 0 ) {
1308 break;
1309 }
1310 } else if ( crt->entry == -1 ) {
1311 crt->visited = 1;
1312 crt->neighbor->visited = 1;
1313 }
1314
1315 crt->visited = 1;
1316 crt->neighbor->visited = 1;
1317
1318
1319 for ( ; ; ) { // loop until next intersection
1320 //nw = create(crt->x, crt->y, old, 0, 0, 0, 0, 0, 0, 0.);
1321 v.setCoords(crt->x, crt->y);
1322 result.addVertex(v);
1323 //old = new;
1324 crt->visited = 1;
1325 crt = crt->next;
1326 //if((crt->status==NS_Intersection) || ((crt->status==NS_IntersectionVertex) && (crt->entry != -1)))
1327 if ( ( crt->status == NS_Intersection ) || ( ( crt->status == NS_IntersectionVertex ) ) ) {
1328 crt->visited = 1;
1329 break;
1330 }
1331 }
1332
1333 //old->nextPoly = root;
1334 //root = old;
1335 }
1336 }
1337 }
1338
1339#ifdef GRAPH_DEBUG_PRINT
1340 /* print graph */
1341 printf("Resulting Graph:\n");
1342 it.init(& result);
1343 while ( it.giveNext(v) ) {
1344 printf("Node [xy %e %e %e]\n", v.coords(0), v.coords(1), 0.0);
1345 }
1346
1347#endif
1348
1349 //printf ("Area=%e\n", result.computeVolume());
1350}
1351
1352
1353
1354Graph :: node *
1355Graph :: prev_node(node *p)
1356{
1357 node *aux = p;
1358 while ( aux && ( aux->status == NS_Intersection ) ) {
1359 aux = aux->prev;
1360 }
1361
1362 return aux;
1363}
1364
1365
1366Graph :: node *
1367Graph :: next_node(node *p)
1368{
1369 node *aux = p;
1370 while ( aux && ( aux->status == NS_Intersection ) ) {
1371 aux = aux->next;
1372 }
1373
1374 return aux;
1375}
1376
1377Graph :: node *
1378Graph :: last_node(node *p)
1379{
1380 node *aux = p;
1381 if ( aux ) {
1382 while ( aux->next ) {
1383 aux = aux->next;
1384 }
1385 }
1386
1387 return aux;
1388}
1389
1390Graph :: node *
1391Graph :: first(node *p)
1392{
1393 node *aux = p;
1394 /*
1395 * if (aux)
1396 * do aux=aux->next;
1397 * while(aux!=p && (aux->status==NS_Vertex || aux->status==NS_Intersection && aux->visited || aux->status==NS_IntersectionVertex && ((aux->entry == -1) || (aux->visited)) ));
1398 * return aux;
1399 */
1400 if ( aux ) {
1401 do {
1402 if ( ( aux->status == NS_Intersection && !aux->visited ) || ( aux->status == NS_IntersectionVertex && ( aux->entry != -1 ) && ( !aux->visited ) ) ) {
1403 return aux;
1404 }
1405
1406 aux = aux->next;
1407 } while ( aux != p );
1408 }
1409
1410 return NULL;
1411}
1412
1413
1414bool
1415Graph :: belongs(node *n, node *v1, node *v2)
1416// tests if node n (assoc to graph1) belongs to line v1,v2 of graph2
1417{
1418 if ( n->status == NS_Vertex ) {
1419 return false;
1420 } else if ( n->neighbor->status == NS_Intersection ) { //Intersection or intersection vertex
1421 if ( ( next_node(n->neighbor) == v2 ) && ( prev_node(n->neighbor) == v1 ) ) {
1422 return true;
1423 }
1424
1425 if ( ( prev_node(n->neighbor) == v2 ) && ( next_node(n->neighbor) == v1 ) ) {
1426 return true;
1427 }
1428 } else { // Vertex to vertex
1429 if ( n->neighbor == v1 || n->neighbor == v2 ) {
1430 return true;
1431 }
1432 }
1433
1434 return false;
1435}
1436
1437double
1438Graph :: dist(double x1, double y1, double x2, double y2)
1439{
1440 return sqrt( ( x1 - x2 ) * ( x1 - x2 ) + ( y1 - y2 ) * ( y1 - y2 ) );
1441}
1442
1443Graph :: node *
1444Graph :: createNode(double x, double y, node *next, node *prev, node *nextPoly,
1445 node *neighbor, nodeStatus st, int entry, int visited, double alpha)
1446{
1447 node *nw = new node;
1448 nw->x = x;
1449 nw->y = y;
1450 nw->next = next;
1451 nw->prev = prev;
1452 if ( prev ) {
1453 nw->prev->next = nw;
1454 }
1455
1456 if ( next ) {
1457 nw->next->prev = nw;
1458 }
1459
1460 nw->nextPoly = nextPoly;
1461 nw->neighbor = neighbor;
1462 nw->status = st;
1463 nw->entry = entry;
1464 nw->visited = visited;
1465 nw->alpha = alpha;
1466 return nw;
1467}
1468
1469void
1470Graph :: insert(node *ins, node *first, node *last)
1471{
1472 node *aux = first;
1473 while ( aux != last && aux->alpha <= ins->alpha ) {
1474 aux = aux->next;
1475 }
1476
1477 ins->next = aux;
1478 ins->prev = aux->prev;
1479 ins->prev->next = ins;
1480 ins->next->prev = ins;
1481}
1482
1483void
1484Graph :: remove(node *rem)
1485{
1486 if ( rem ) {
1487 rem->neighbor->entry = 0;
1488 rem->prev->next = rem->next;
1489 rem->next->prev = rem->prev;
1490 delete rem;
1491 rem = NULL;
1492 }
1493}
1494
1495bool
1496Graph :: intersectionExist(node *p1, node *p2, node *q1, node *q2)
1497{
1498 node *aux = q1->next;
1499 while ( aux != q2 ) {
1500 if ( aux->status == NS_Intersection ) {
1501 if ( ( aux->neighbor->prev == p1 ) && ( aux->neighbor->next == p2 ) ) {
1502 return true;
1503 }
1504 }
1505
1506 aux = aux->next;
1507 }
1508
1509 return false;
1510}
1511
1512
1513bool
1514Graph :: testCollapsedEdge(node *p1, node *p2)
1515{
1516 if ( dist(p1->x, p1->y, p2->x, p2->y) <= GT_EPS ) {
1517 return true;
1518 } else {
1519 return false;
1520 }
1521}
1522
1523void
1524Graph :: removeIntersectionIfExist(node *p1, node *p2, node *q1, node *q2)
1525{
1526 node *anext;
1527 node *aux = q1->next;
1528 while ( aux != q2 ) {
1529 anext = aux->next;
1530 if ( aux->status == NS_Intersection ) {
1531 if ( ( prev_node(aux->neighbor->prev) == p1 ) && ( next_node(aux->neighbor->next) == p2 ) ) {
1532 // remove aux and aux->neighbor
1533 node *neighbor = aux->neighbor;
1534 aux->prev->next = aux->next;
1535 aux->next->prev = aux->prev;
1536 delete aux;
1537
1538 neighbor->prev->next = neighbor->next;
1539 neighbor->next->prev = neighbor->prev;
1540 delete neighbor;
1541 }
1542 }
1543
1544 aux = anext;
1545 }
1546}
1547
1548int
1549Graph :: vertex2IntersectionVertex(node *v, node *l1, node *l2)
1550{
1551 // transforms given vertex as IntersectionVertex on given line
1552 // checks for coorect topology, ie check whether previos intersections of lines from IV to l1,l2 are present
1553 // returns 1 if newly formed
1554 // returns 0 if already exist
1555 node *nv = next_node(v->next);
1556 node *pv = prev_node(v->prev);
1557
1558 if ( v->status == NS_IntersectionVertex ) {
1559 if ( v->neighbor->status == NS_IntersectionVertex ) {
1560 if ( ( v->neighbor == l1 ) || ( v->neighbor == l2 ) ) {
1561 return 0;
1562 } else {
1563 THROW_GT_EXCEPTIONM("Graph::vertex2IntersectionVertex: topology error, neighbor is unrelated intersectionVertex");
1564 }
1565 } else if ( ( prev_node(v->neighbor->prev) == l1 ) && ( next_node(v->neighbor->next) == l2 ) ) {
1566 return 0;
1567 } else if ( ( prev_node(v->neighbor->prev) == l2 ) && ( next_node(v->neighbor->next) == l1 ) ) {
1568 return 0;
1569 } else {
1570 // intersection vertex associated to neiborhing edge -> move it to common vertex
1571 // merge vertex together, remove intersection now associated to IntersectionVertex
1572 if ( ( l1 == prev_node(v->neighbor->prev) ) || ( l1 == next_node(v->neighbor->next) ) ) {
1573 merge2vertex(v, l1);
1574 } else if ( ( l2 == prev_node(v->neighbor->prev) ) || ( l2 == next_node(v->neighbor->next) ) ) {
1575 merge2vertex(v, l2);
1576 } else {
1577 THROW_GT_EXCEPTIONM("Graph::vertex2IntersectionVertex: topology error");
1578 }
1579
1580 return 0;
1581 }
1582 } else {
1583 removeIntersectionIfExist(pv, v, l1, l2);
1584 removeIntersectionIfExist(v, nv, l1, l2);
1585
1586
1587 v->status = NS_IntersectionVertex;
1588 return 1;
1589 }
1590}
1591
1592
1593void
1594Graph :: testNewIntersectionVertexEdgeCollapse(node *v, node *l1, node *l2)
1595{
1596 double v_par, vp_par, vn_par;
1597
1598 if ( l1->status == NS_IntersectionVertex && l1->neighbor->status == NS_Intersection ) {
1599 if ( next_node(l1->neighbor->next) == v ) {
1600 l1->neighbor->entry = -1; // to v
1601 l1->entry = -1; // from l1
1602 } else if ( prev_node(l1->neighbor->prev) == v ) {
1603 l1->entry = -1; // from l1
1604 v->entry = -1; // from v
1605 }
1606 } else if ( l2->status == NS_IntersectionVertex && l2->neighbor->status == NS_Intersection ) {
1607 if ( next_node(l2->neighbor->next) == v ) {
1608 l2->neighbor->entry = -1; // to v
1609 v->neighbor->entry = -1; // to l2
1610 } else if ( prev_node(l2->neighbor->prev) == v ) {
1611 v->entry = -1; // from v
1612 v->neighbor->entry = -1; // to l2
1613 }
1614 } else if ( belongs(prev_node(v->prev), l1, l2) ) {
1615 node *vp = prev_node(v->prev);
1616 vp->entry = -1;
1617
1618 if ( vp->neighbor == l1 ) {
1619 vp_par = 0.;
1620 } else if ( vp->neighbor == l2 ) {
1621 vp_par = 1.;
1622 } else {
1623 vp_par = vp->neighbor->alpha;
1624 }
1625
1626 if ( v->neighbor == l1 ) {
1627 v_par = 0.;
1628 } else if ( v->neighbor == l2 ) {
1629 v_par = 1.;
1630 } else {
1631 v_par = v->neighbor->alpha;
1632 }
1633
1634 if ( v_par < vp_par ) {
1635 v->neighbor->entry = -1;
1636 } else {
1637 vp->neighbor->entry = -1;
1638 }
1639 } else if ( belongs(next_node(v->next), l1, l2) ) {
1640 node *vn = next_node(v->next);
1641 v->entry = -1;
1642
1643 if ( vn->neighbor == l1 ) {
1644 vn_par = 0.;
1645 } else if ( vn->neighbor == l2 ) {
1646 vn_par = 1.;
1647 } else {
1648 vn_par = vn->neighbor->alpha;
1649 }
1650
1651 if ( v->neighbor == l1 ) {
1652 v_par = 0.;
1653 } else if ( v->neighbor == l2 ) {
1654 v_par = 1.;
1655 } else {
1656 v_par = v->neighbor->alpha;
1657 }
1658
1659 if ( v_par < vn_par ) {
1660 v->neighbor->entry = -1;
1661 } else {
1662 vn->neighbor->entry = -1;
1663 }
1664 }
1665}
1666
1667
1668
1669void
1670Graph :: merge2vertex(node *v1, node *v2)
1671{
1672 bool boundary = false;
1673 if ( v1->status == NS_IntersectionVertex && v2->status == NS_IntersectionVertex ) {
1674 if ( v1->neighbor->status == NS_IntersectionVertex && v2->neighbor->status == NS_IntersectionVertex &&
1675 v1->neighbor == v2 && v2->neighbor == v1 ) {
1676 return;
1677 } else if ( v1->neighbor->status == NS_Intersection && v2->neighbor->status == NS_Intersection ) {
1678 remove(v1->neighbor);
1679 remove(v2->neighbor);
1680 v1->neighbor = v2;
1681 v2->neighbor = v1;
1682 v1->status = v2->status = NS_IntersectionVertex;
1683 } else {
1684 THROW_GT_EXCEPTIONM("Graph::merge2vertex: consistency error (one of merged vertices (or both) linked to another vertex already");
1685 }
1686 } else if ( v1->status == NS_IntersectionVertex ) {
1687 if ( v1->neighbor->status == NS_Intersection ) {
1688 boundary = ( v1->neighbor->entry == -1 );
1689 remove(v1->neighbor);
1690 }
1691
1692 v1->neighbor = v2;
1693 v2->neighbor = v1;
1694 v2->status = NS_IntersectionVertex;
1695 if ( boundary ) {
1696 v2->entry = -1;
1697 }
1698 } else if ( v2->status == NS_IntersectionVertex ) {
1699 if ( v2->neighbor->status == NS_Intersection ) {
1700 boundary = ( v2->neighbor->entry == -1 );
1701 remove(v2->neighbor);
1702 }
1703
1704 v2->neighbor = v1;
1705 v1->neighbor = v2;
1706 v1->status = NS_IntersectionVertex;
1707 if ( boundary ) {
1708 v1->entry = -1;
1709 }
1710 } else {
1711 v1->neighbor = v2;
1712 v2->neighbor = v1;
1713 v1->status = v2->status = NS_IntersectionVertex;
1714 }
1715}
1716
1717int
1718Graph :: testPoint(node *s, double x, double y) const
1719{
1720 bool oddNODES = false;
1721 double x1, x2, y1, y2;
1722
1723 if ( s ) {
1724 node *auxs = s;
1725 do {
1726 x1 = auxs->x;
1727 y1 = auxs->y;
1728 x2 = auxs->next->x;
1729 y2 = auxs->next->y;
1730
1731 if ( ( ( y1 < y ) && ( y2 >= y ) ) ||
1732 ( ( y2 < y ) && ( y1 >= y ) ) ) {
1733 if ( x1 + ( y - y1 ) / ( y2 - y1 ) * ( x2 - x1 ) < x ) {
1734 oddNODES = !oddNODES;
1735 }
1736 }
1737
1738 auxs = auxs->next;
1739 } while ( auxs != s );
1740 }
1741
1742 return oddNODES;
1743}
1744
1745
1746void
1747Graph :: printYourself()
1748{
1749 node *auxs, *auxc;
1750 /* print graph */
1751 printf("Graph 1/2:\n");
1752 auxs = s;
1753 do {
1754 if ( auxs->status == NS_Vertex ) {
1755 printf("Vertex [xy %e %e %e] [e=%d]\n", auxs->x, auxs->y, 0.0, auxs->entry);
1756 } else if ( auxs->status == NS_IntersectionVertex ) {
1757 printf("IVertex [xy %e %e %e] [e=%d]\n", auxs->x, auxs->y, 0.0, auxs->entry);
1758 } else {
1759 printf("Intrsc [xy %e %e %e] [e=%d]\n", auxs->x, auxs->y, 0.0, auxs->entry);
1760 }
1761
1762 auxs = auxs->next;
1763 } while ( auxs != s );
1764
1765 /* print graph */
1766 printf("Graph 2/2:\n");
1767 auxc = c;
1768 do {
1769 if ( auxc->status == NS_Vertex ) {
1770 printf("Vertex [xy %e %e %e] [e=%d]\n", auxc->x, auxc->y, 0.0, auxc->entry);
1771 } else if ( auxc->status == NS_IntersectionVertex ) {
1772 printf("IVertex [xy %e %e %e] [e=%d]\n", auxc->x, auxc->y, 0.0, auxc->entry);
1773 } else {
1774 printf("Intrsc [xy %e %e %e] [e=%d]\n", auxc->x, auxc->y, 0.0, auxc->entry);
1775 }
1776
1777 auxc = auxc->next;
1778 } while ( auxc != c );
1779
1780 /* end of debug printing */
1781}
1782
1783
1784void
1785GT_Exception :: print()
1786{
1787 fprintf(stderr, "\nGT_Exception thrown in %s:%d\n", file, line);
1788 if ( msg ) {
1789 fprintf(stderr, "msg: %s\n", msg);
1790 }
1791}
1792} // end namespace oofem
remove
void remove(std::vector< T > &cont, const T &t)
Definition CSG.h:1122
oofem::GT_Exception::msg
const char * msg
Definition geotoolbox.h:248
oofem::GT_Exception::line
int line
Definition geotoolbox.h:249
oofem::GT_Exception::file
const char * file
Definition geotoolbox.h:248
oofem::Graph::belongs
bool belongs(node *n, node *v1, node *v2)
Definition geotoolbox.C:1415
oofem::Graph::testIfIntersect
int testIfIntersect(node *p1, node *p2, node *q1, node *q2, double *alpha_p, double *alpha_q, double *xint, double *yint)
Definition geotoolbox.C:463
oofem::Graph::testPoint
int testPoint(node *poly, double x, double y) const
Definition geotoolbox.C:1718
oofem::Graph::prev_node
node * prev_node(node *p)
Definition geotoolbox.C:1355
oofem::Graph::testIfCoincident
int testIfCoincident(node *p1, node *p2, node *q1, node *q2, double *alpha_1, double *alpha_2)
Definition geotoolbox.C:255
oofem::Graph::printYourself
void printYourself()
Definition geotoolbox.C:1747
oofem::Graph::testCollapsedEdge
bool testCollapsedEdge(node *p1, node *p2)
Definition geotoolbox.C:1514
oofem::Graph::createNode
node * createNode(double x, double y, node *next, node *prev, node *nextPoly, node *neighbor, nodeStatus st, int entry, int visited, double alpha)
Definition geotoolbox.C:1444
oofem::Graph::dist
double dist(double x1, double y1, double x2, double y2)
Definition geotoolbox.C:1438
oofem::Graph::merge2vertex
void merge2vertex(node *v1, node *v2)
Definition geotoolbox.C:1670
oofem::Graph::testNewIntersectionVertexEdgeCollapse
void testNewIntersectionVertexEdgeCollapse(node *v, node *l1, node *l2)
Definition geotoolbox.C:1594
oofem::Graph::insert
void insert(node *ins, node *first, node *last)
Definition geotoolbox.C:1470
oofem::Graph::next_node
node * next_node(node *p)
Definition geotoolbox.C:1367
oofem::Graph::vertex2IntersectionVertex
int vertex2IntersectionVertex(node *v, node *l1, node *l2)
Definition geotoolbox.C:1549
oofem::Graph::NS_IntersectionVertex
@ NS_IntersectionVertex
Definition geotoolbox.h:194
oofem::Graph::removeIntersectionIfExist
void removeIntersectionIfExist(node *p1, node *p2, node *q1, node *q2)
Definition geotoolbox.C:1524
oofem::Graph::first
node * first(node *p)
Definition geotoolbox.C:1391
oofem::Polygon::addVertex
void addVertex(Vertex v)
Definition geotoolbox.h:96
oofem::Polygon::testPoint
int testPoint(double x, double y) const
Definition geotoolbox.C:47
oofem::Vertex::coords
FloatArray coords
Definition geotoolbox.h:63
oofem::Vertex::setCoords
void setCoords(double x, double y)
Definition geotoolbox.h:77
GT_EPS
#define GT_EPS
Definition geotoolbox.h:50
GT_TEPS
#define GT_TEPS
Definition geotoolbox.h:52
THROW_GT_EXCEPTIONM
#define THROW_GT_EXCEPTIONM(m)
Definition geotoolbox.h:260
oofem::min
FloatArrayF< N > min(const FloatArrayF< N > &a, const FloatArrayF< N > &b)
Definition floatarrayf.h:522
oofem::sgn
double sgn(double i)
Returns the signum of given value (if value is < 0 returns -1, otherwise returns 1).
Definition mathfem.h:104
gc
oofem::oofegGraphicContext gc[OOFEG_LAST_LAYER]
oofem::Graph::node::nextPoly
struct node * nextPoly
Definition geotoolbox.h:200
oofem::Graph::node::prev
struct node * prev
Definition geotoolbox.h:199
oofem::Graph::node::neighbor
struct node * neighbor
Definition geotoolbox.h:201
oofem::Graph::node::status
nodeStatus status
Definition geotoolbox.h:203
oofem::Graph::node::next
struct node * next
Definition geotoolbox.h:198
This page is part of the OOFEM-3.0 documentation. Copyright Copyright (C) 1994-2025 Borek Patzak Bořek Patzák
Project e-mail: oofem@fsv.cvut.cz
Generated at for OOFEM by doxygen 1.15.0 written by Dimitri van Heesch, © 1997-2011