OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
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 - 2013 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 
42 namespace oofem {
43 //#define GRAPH_DEBUG_PRINT
44 #define GRAPH_LENGTH_FRAC 1.e-4
45 
46 int
47 Polygon :: 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 
71 double
72 Polygon :: 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 
97 double
98 Polygon :: 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
151 GraphicObj *
152 Polygon :: 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 
224 Graph :: ~Graph()
225 {
226  this->clear();
227 }
228 
229 
230 void
231 Graph :: 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 
254 int
255 Graph :: 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 
462 int
463 Graph :: 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 
574 void
575 Graph :: 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 
1354 Graph :: node *
1355 Graph :: 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 
1366 Graph :: node *
1367 Graph :: 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 
1377 Graph :: node *
1378 Graph :: 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 
1390 Graph :: node *
1391 Graph :: 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 
1414 bool
1415 Graph :: 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 
1437 double
1438 Graph :: dist(double x1, double y1, double x2, double y2)
1439 {
1440  return sqrt( ( x1 - x2 ) * ( x1 - x2 ) + ( y1 - y2 ) * ( y1 - y2 ) );
1441 }
1442 
1443 Graph :: node *
1444 Graph :: 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 
1469 void
1470 Graph :: 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 
1483 void
1484 Graph :: 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 
1495 bool
1496 Graph :: 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 
1513 bool
1514 Graph :: 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 
1523 void
1524 Graph :: 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 
1548 int
1549 Graph :: 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 
1593 void
1594 Graph :: 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 
1669 void
1670 Graph :: 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 
1717 int
1718 Graph :: 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 
1746 void
1747 Graph :: 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 
1784 void
1785 GT_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
GT_TEPS
#define GT_TEPS
Definition: geotoolbox.h:52
THROW_GT_EXCEPTIONM
#define THROW_GT_EXCEPTIONM(m)
Definition: geotoolbox.h:260
oofem::Polygon::computeVolume
double computeVolume() const
Definition: geotoolbox.C:72
GT_EPS
#define GT_EPS
Definition: geotoolbox.h:50
oofem::Polygon::pointDistance
double pointDistance(double x, double y) const
Definition: geotoolbox.C:98
oofem::sgn
double sgn(double i)
Returns the signum of given value (if value is < 0 returns -1, otherwise returns 1) ...
Definition: mathfem.h:91
oofem::Graph::node::neighbor
struct node * neighbor
Definition: geotoolbox.h:201
oofem::Graph::merge2vertex
void merge2vertex(node *v1, node *v2)
Definition: geotoolbox.C:1670
gc
oofem::oofegGraphicContext gc[OOFEG_LAST_LAYER]
oofem::Graph::printYourself
void printYourself()
Definition: geotoolbox.C:1747
oofem::Graph::prev_node
node * prev_node(node *p)
Definition: geotoolbox.C:1355
oofem::Graph::node::status
nodeStatus status
Definition: geotoolbox.h:203
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::Vertex::coords
FloatArray coords
Definition: geotoolbox.h:63
oofem::Graph::node::next
struct node * next
Definition: geotoolbox.h:198
oofem::Graph::dist
double dist(double x1, double y1, double x2, double y2)
Definition: geotoolbox.C:1438
oofem::Graph::vertex2IntersectionVertex
int vertex2IntersectionVertex(node *v, node *l1, node *l2)
Definition: geotoolbox.C:1549
oofem::GT_Exception::print
void print()
Definition: geotoolbox.C:1785
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::Graph::clip
void clip(Polygon &result, const Polygon &a, const Polygon &b)
Definition: geotoolbox.C:575
oofem::Graph::last_node
node * last_node(node *p)
Definition: geotoolbox.C:1378
oofem::Graph::testPoint
int testPoint(node *poly, double x, double y) const
Definition: geotoolbox.C:1718
oofem::Graph::insert
void insert(node *ins, node *first, node *last)
Definition: geotoolbox.C:1470
oofem::Polygon
Class representing 2D polygon.
Definition: geotoolbox.h:91
oofem::Vertex::setCoords
void setCoords(double x, double y)
Definition: geotoolbox.h:77
oofem::Graph::intersectionExist
bool intersectionExist(node *p1, node *p2, node *q1, node *q2)
Definition: geotoolbox.C:1496
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::node::nextPoly
struct node * nextPoly
Definition: geotoolbox.h:200
oofem::Vertex
Class representing vertex.
Definition: geotoolbox.h:60
oofem::Graph::remove
void remove(node *n)
Definition: geotoolbox.C:1484
oofem::min
int min(int i, int j)
Returns smaller value from two given decimals.
Definition: mathfem.h:59
oofem::Graph::testNewIntersectionVertexEdgeCollapse
void testNewIntersectionVertexEdgeCollapse(node *v, node *l1, node *l2)
Definition: geotoolbox.C:1594
oofem::Polygon::testPoint
int testPoint(double x, double y) const
Definition: geotoolbox.C:47
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)
Create new node struct.
Definition: geotoolbox.C:1444
oofem::Polygon::draw
GraphicObj * draw(oofegGraphicContext &, bool filled, int layer=OOFEG_DEBUG_LAYER)
Definition: geotoolbox.C:152
oofem::Graph::testCollapsedEdge
bool testCollapsedEdge(node *p1, node *p2)
Definition: geotoolbox.C:1514
oofem
the oofem namespace is to define a context or scope in which all oofem names are defined.
oofem::Polygon::addVertex
void addVertex(Vertex v)
Definition: geotoolbox.h:96
oofem::Graph::belongs
bool belongs(node *n, node *v1, node *v2)
Definition: geotoolbox.C:1415
oofem::Graph::clear
void clear()
Definition: geotoolbox.C:231
oofem::Graph::node::prev
struct node * prev
Definition: geotoolbox.h:199
oofem::Polygon::PolygonEdgeIterator::giveNext
int giveNext(Vertex &p1, Vertex &p2)
Definition: geotoolbox.h:122
oofem::Graph::next_node
node * next_node(node *p)
Definition: geotoolbox.C:1367
This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:28 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011