1 /*
2 ** License Applicability. Except to the extent portions of this file are
3 ** made subject to an alternative license as permitted in the SGI Free
4 ** Software License B, Version 1.1 (the "License"), the contents of this
5 ** file are subject only to the provisions of the License. You may not use
6 ** this file except in compliance with the License. You may obtain a copy
7 ** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
8 ** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
9 **
10 ** http://oss.sgi.com/projects/FreeB
11 **
12 ** Note that, as provided in the License, the Software is distributed on an
13 ** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
14 ** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
15 ** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
16 ** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
17 **
18 ** Original Code. The Original Code is: OpenGL Sample Implementation,
19 ** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
20 ** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
21 ** Copyright in any portions created by third parties is as indicated
22 ** elsewhere herein. All Rights Reserved.
23 **
24 ** Additional Notice Provisions: The application programming interfaces
25 ** established by SGI in conjunction with the Original Code are The
26 ** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
27 ** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
28 ** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
29 ** Window System(R) (Version 1.3), released October 19, 1998. This software
30 ** was created using the OpenGL(R) version 1.2.1 Sample Implementation
31 ** published by SGI, but has not been independently verified as being
32 ** compliant with the OpenGL(R) version 1.2.1 Specification.
33 **
34 */
35 /*
36 */
37 
38 #include "gluos.h"
39 //#include <stdlib.h>
40 //#include <stdio.h>
41 //#include <GL/gl.h>
42 
43 //#include "glimports.h"
44 #include "zlassert.h"
45 
46 #include "monoChain.h"
47 #include "quicksort.h"
48 #include "searchTree.h"
49 #include "polyUtil.h"
50 
51 #ifndef max
52 #define max(a,b) ((a>b)? a:b)
53 #endif
54 #ifndef min
55 #define min(a,b) ((a>b)? b:a)
56 #endif
57 
58 extern Int isCusp(directedLine *v);
59 extern Int deleteRepeatDiagonals(Int num_diagonals, directedLine** diagonal_vertices, directedLine** new_vertices);
60 
61 //for debug purpose only
62 #if 0 // UNUSED
63 static void drawDiagonals(Int num_diagonals, directedLine** diagonal_vertices)
64 {
65   Int i;
66   for(i=0; i<num_diagonals; i++)
67     {
68       glBegin(GL_LINE);
69       glVertex2fv(diagonal_vertices[2*i]->head());
70       glVertex2fv(diagonal_vertices[2*i+1]->head());
71       glEnd();
72     }
73 }
74 #endif
75 
76 /*given (x_1, y_1) and (x_2, y_2), and y
77  *return x such that (x,y) is on the line
78  */
79 inline Real intersectHoriz(Real x1, Real y1, Real x2, Real y2, Real y)
80 {
81   return ((y2==y1)? (x1+x2)*0.5 : x1 + ((y-y1)/(y2-y1)) * (x2-x1));
82 }
83 
84 //compare the heads of the two chains
85 static int compChainHeadInY(monoChain* mc1, monoChain* mc2)
86 {
87   return compV2InY(mc1->getHead()->head(), mc2->getHead()->head());
88 }
89 
90 monoChain::monoChain(directedLine* cHead, directedLine* cTail)
91 {
92   chainHead = cHead;
93   chainTail = cTail;
94   next = this;
95   prev = this;
96 
97   nextPolygon = NULL;
98 
99   //compute bounding box
100   directedLine* temp;
101   minX = maxX = chainTail->head()[0];
102   minY = maxY = chainTail->head()[1];
103 
104   for(temp=chainHead; temp!=cTail; temp = temp->getNext())
105     {
106       if(temp->head()[0] < minX)
107 	minX = temp->head()[0];
108       if(temp->head()[0] > maxX)
109 	maxX = temp->head()[0];
110 
111       if(temp->head()[1] < minY)
112 	minY = temp->head()[1];
113       if(temp->head()[1] > maxY)
114 	maxY = temp->head()[1];
115     }
116 
117   //check whether the chain is increasing or decreasing
118   if(chainHead->compInY(chainTail) <0)
119     isIncrease = 1;
120   else
121     isIncrease = 0;
122 
123   //initilize currrent, this is used for accelerating search
124   if(isIncrease)
125     current = chainHead;
126   else
127     current = chainTail;
128 
129   isKey = 0;
130   keyY = 0;
131 }
132 
133 //insert a new line between prev and this
134 void monoChain::insert(monoChain* nc)
135 {
136   nc->next = this;
137   nc->prev = prev;
138   prev->next = nc;
139   prev = nc;
140 }
141 
142 void monoChain::deleteLoop()
143 {
144   monoChain *temp, *tempNext;
145   prev->next = NULL;
146   for(temp=this; temp != NULL; temp = tempNext)
147     {
148       tempNext = temp->next;
149       delete temp;
150     }
151 }
152 
153 void monoChain::deleteLoopList()
154 {
155   monoChain *temp, *tempNext;
156   for(temp=this; temp != NULL; temp = tempNext)
157     {
158       tempNext = temp->nextPolygon;
159       temp->deleteLoop();
160     }
161 }
162 
163 Int monoChain::toArraySingleLoop(monoChain** array, Int index)
164 {
165   monoChain *temp;
166   array[index++] = this;
167   for(temp = next; temp != this; temp = temp->next)
168     {
169       array[index++] = temp;
170     }
171   return index;
172 }
173 
174 monoChain** monoChain::toArrayAllLoops(Int& num_chains)
175 {
176   num_chains = numChainsAllLoops();
177   monoChain **ret =  (monoChain**) malloc(sizeof(monoChain*) * num_chains);
178   assert(ret);
179   monoChain *temp;
180   Int index = 0;
181   for(temp = this; temp != NULL; temp=temp->nextPolygon){
182     index = temp->toArraySingleLoop(ret, index);
183   }
184   return ret;
185 }
186 
187 Int monoChain::numChainsSingleLoop()
188 {
189   Int ret=0;
190   monoChain* temp;
191   if(next == this) return 1;
192   ret = 1;
193   for(temp=next; temp != this; temp = temp->next)
194     ret++;
195   return ret;
196 }
197 
198 Int monoChain::numChainsAllLoops()
199 {
200   Int ret=0;
201   monoChain *temp;
202   for(temp =this; temp != NULL; temp = temp->nextPolygon)
203     ret += temp->numChainsSingleLoop();
204   return ret;
205 }
206 
207 //update 'current'
208 Real monoChain::chainIntersectHoriz(Real y)
209 {
210   directedLine* temp;
211   if(isIncrease)
212     {
213       for(temp= current; temp != chainTail; temp = temp->getNext())
214 	{
215 	  if(temp->head()[1] > y)
216 	    break;
217 	}
218       current = temp->getPrev();
219     }
220   else
221     {
222       for(temp = current; temp != chainHead; temp = temp->getPrev())
223 	{
224 	  if(temp->head()[1] > y)
225 	    break;
226 	}
227       current = temp->getNext();
228     }
229   return intersectHoriz(current->head()[0], current->head()[1], current->tail()[0], current->tail()[1], y);
230 }
231 
232 monoChain* directedLineLoopToMonoChainLoop(directedLine* loop)
233 {
234   directedLine *temp;
235   monoChain *ret=NULL;
236 
237   //find the first cusp
238   directedLine *prevCusp=NULL;
239   directedLine *firstCusp;
240 
241   if(isCusp(loop))
242     prevCusp = loop;
243   else
244     {
245       for(temp = loop->getNext(); temp != loop; temp = temp->getNext())
246 	if(isCusp(temp))
247 	  break;
248       prevCusp = temp;
249     }
250   firstCusp = prevCusp;
251 //printf("first cusp is (%f,%f), (%f,%f), (%f,%f)\n", prevCusp->getPrev()->head()[0], prevCusp->getPrev()->head()[1], prevCusp->head()[0], prevCusp->head()[1], prevCusp->tail()[0], prevCusp->tail()[1]);
252 
253   for(temp = prevCusp->getNext(); temp != loop; temp = temp->getNext())
254     {
255       if(isCusp(temp))
256 	{
257 //printf("the cusp is (%f,%f), (%f,%f), (%f,%f)\n", temp->getPrev()->head()[0], temp->getPrev()->head()[1], temp->head()[0], temp->head()[1], temp->tail()[0], temp->tail()[1]);
258 	  if(ret == NULL)
259 	    {
260 	      ret = new monoChain(prevCusp, temp);
261 	    }
262 	  else
263 	    ret->insert(new monoChain(prevCusp, temp));
264 	  prevCusp = temp;
265 	}
266     }
267   assert(ret);
268   ret->insert(new monoChain(prevCusp, firstCusp));
269 
270   return ret;
271 }
272 
273 monoChain* directedLineLoopListToMonoChainLoopList(directedLine* list)
274 {
275   directedLine* temp;
276   monoChain* mc;
277   monoChain* mcEnd;
278   mc = directedLineLoopToMonoChainLoop(list);
279   mcEnd = mc;
280   for(temp = list->getNextPolygon(); temp != NULL; temp = temp->getNextPolygon())
281     {
282       monoChain *newLoop = directedLineLoopToMonoChainLoop(temp);
283       mcEnd->setNextPolygon(newLoop);
284       mcEnd = newLoop;
285     }
286   return mc;
287 }
288 
289 /*compare two edges of a polygon.
290  *edge A < edge B if there is a horizontal line so that the intersection
291  *with A is to the left of the intersection with B.
292  *This function is used in sweepY for the dynamic search tree insertion to
293  *order the edges.
294  * Implementation: (x_1,y_1) and (x_2, y_2)
295  */
296 static Int compEdges(directedLine *e1, directedLine *e2)
297 {
298   Real* head1 = e1->head();
299   Real* tail1 = e1->tail();
300   Real* head2 = e2->head();
301   Real* tail2 = e2->tail();
302 /*
303   Real h10 = head1[0];
304   Real h11 = head1[1];
305   Real t10 = tail1[0];
306   Real t11 = tail1[1];
307   Real h20 = head2[0];
308   Real h21 = head2[1];
309   Real t20 = tail2[0];
310   Real t21 = tail2[1];
311 */
312   Real e1_Ymax, e1_Ymin, e2_Ymax, e2_Ymin;
313 /*
314   if(h11>t11) {
315     e1_Ymax= h11;
316     e1_Ymin= t11;
317   }
318   else{
319     e1_Ymax = t11;
320     e1_Ymin = h11;
321   }
322 
323   if(h21>t21) {
324     e2_Ymax= h21;
325     e2_Ymin= t21;
326   }
327   else{
328     e2_Ymax = t21;
329     e2_Ymin = h21;
330   }
331 */
332 
333   if(head1[1]>tail1[1]) {
334     e1_Ymax= head1[1];
335     e1_Ymin= tail1[1];
336   }
337   else{
338     e1_Ymax = tail1[1];
339     e1_Ymin = head1[1];
340   }
341 
342   if(head2[1]>tail2[1]) {
343     e2_Ymax= head2[1];
344     e2_Ymin= tail2[1];
345   }
346   else{
347     e2_Ymax = tail2[1];
348     e2_Ymin = head2[1];
349   }
350 
351 
352   /*Real e1_Ymax = max(head1[1], tail1[1]);*/ /*max(e1->head()[1], e1->tail()[1]);*/
353   /*Real e1_Ymin = min(head1[1], tail1[1]);*/ /*min(e1->head()[1], e1->tail()[1]);*/
354   /*Real e2_Ymax = max(head2[1], tail2[1]);*/ /*max(e2->head()[1], e2->tail()[1]);*/
355   /*Real e2_Ymin = min(head2[1], tail2[1]);*/ /*min(e2->head()[1], e2->tail()[1]);*/
356 
357   Real Ymax = min(e1_Ymax, e2_Ymax);
358   Real Ymin = max(e1_Ymin, e2_Ymin);
359 
360   Real y = 0.5*(Ymax + Ymin);
361 
362 /*  Real x1 = intersectHoriz(e1->head()[0], e1->head()[1], e1->tail()[0], e1->tail()[1], y);
363   Real x2 = intersectHoriz(e2->head()[0], e2->head()[1], e2->tail()[0], e2->tail()[1], y);
364 */
365 /*
366   Real x1 = intersectHoriz(h10, h11, t10, t11, y);
367   Real x2 = intersectHoriz(h20, h21, t20, t21, y);
368 */
369   Real x1 = intersectHoriz(head1[0], head1[1], tail1[0], tail1[1], y);
370   Real x2 = intersectHoriz(head2[0], head2[1], tail2[0], tail2[1], y);
371 
372   if(x1<= x2) return -1;
373   else return 1;
374 }
375 
376 Int  compChains(monoChain* mc1, monoChain* mc2)
377 {
378   Real y;
379   assert(mc1->isKey || mc2->isKey);
380   if(mc1->isKey)
381     y = mc1->keyY;
382   else
383     y = mc2->keyY;
384   directedLine *d1 = mc1->find(y);
385   directedLine *d2 = mc2->find(y);
386   mc2->find(y);
387 //  Real x1 = mc1->chainIntersectHoriz(y);
388 //  Real x2 = mc2->chainIntersectHoriz(y);
389   return   compEdges(d1, d2);
390 }
391 
392 //this function modifies current for efficiency
393 directedLine* monoChain::find(Real y)
394 {
395   directedLine *ret;
396   directedLine *temp;
397   assert(current->head()[1] <= y);
398   if(isIncrease)
399     {
400       assert(chainTail->head()[1] >=y);
401       for(temp=current; temp!=chainTail; temp = temp->getNext())
402 	{
403 	  if(temp->head()[1] > y)
404 	    break;
405 	}
406       current = temp->getPrev();
407       ret = current;
408     }
409   else
410     {
411       for(temp=current; temp != chainHead; temp = temp->getPrev())
412 	{
413 	  if(temp->head()[1] > y)
414 	    break;
415 	}
416       current = temp->getNext();
417       ret = temp;
418     }
419   return ret;
420 }
421 
422 void monoChain::printOneChain()
423 {
424   directedLine* temp;
425   for(temp = chainHead; temp != chainTail; temp = temp->getNext())
426     {
427       printf("(%f,%f) ", temp->head()[0], temp->head()[1]);
428     }
429   printf("(%f,%f) \n", chainTail->head()[0], chainTail->head()[1]);
430 }
431 
432 void monoChain::printChainLoop()
433 {
434   monoChain* temp;
435   this->printOneChain();
436   for(temp = next; temp != this; temp = temp->next)
437     {
438       temp->printOneChain();
439     }
440   printf("\n");
441 }
442 
443 void monoChain::printAllLoops()
444 {
445   monoChain* temp;
446   for(temp=this; temp != NULL; temp = temp->nextPolygon)
447     temp->printChainLoop();
448 }
449 
450 //return 1 if error occures
451 Int MC_sweepY(Int nVertices, monoChain** sortedVertices, sweepRange** ret_ranges)
452 {
453   Int i;
454   Real keyY;
455   Int errOccur=0;
456 //printf("enter MC_sweepY\n");
457 //printf("nVertices=%i\n", nVertices);
458   /*for each vertex in the sorted list, update the binary search tree.
459    *and store the range information for each vertex.
460    */
461   treeNode* searchTree = NULL;
462 //printf("nVertices=%i\n", nVertices);
463   for(i=0; i<nVertices; i++)
464     {
465       monoChain* vert = sortedVertices[i];
466       keyY = vert->getHead()->head()[1]; //the sweep line
467       directedLine *dline = vert->getHead();
468       directedLine *dlinePrev = dline->getPrev();
469       if(isBelow(dline, dline) && isBelow(dline, dlinePrev))
470 	{
471 //printf("case 1\n");
472 	  //this<v and prev < v
473 	  //delete both edges
474 	  vert->isKey = 1;
475 	  vert->keyY = keyY;
476 	  treeNode* thisNode = TreeNodeFind(searchTree, vert, (Int (*) (void *, void *))compChains);
477 	  vert->isKey = 0;
478 
479 	  vert->getPrev()->isKey = 1;
480 	  vert->getPrev()->keyY = keyY;
481 	  treeNode* prevNode = TreeNodeFind(searchTree, vert->getPrev(), (Int (*) (void *, void *))compChains);
482 	  vert->getPrev()->isKey = 0;
483 
484 	  if(cuspType(dline) == 1)//interior cusp
485 	    {
486 
487 	      treeNode* leftEdge = TreeNodePredecessor(prevNode);
488 	      treeNode* rightEdge = TreeNodeSuccessor(thisNode);
489 	      if(leftEdge == NULL ||  rightEdge == NULL)
490 		{
491 		  errOccur = 1;
492 		  goto JUMP_HERE;
493 		}
494 
495 	      directedLine* leftEdgeDline = ((monoChain* ) leftEdge->key)->find(keyY);
496 
497 
498 
499 	      directedLine* rightEdgeDline = ((monoChain* ) rightEdge->key)->find(keyY);
500 
501 	      ret_ranges[i] = sweepRangeMake(leftEdgeDline, 1, rightEdgeDline, 1);
502 	    }
503 	  else /*exterior cusp*/
504 	    {
505 	      ret_ranges[i] = sweepRangeMake( dline, 1, dlinePrev, 1);
506 	    }
507 
508 	  searchTree = TreeNodeDeleteSingleNode(searchTree, thisNode);
509 	  searchTree = TreeNodeDeleteSingleNode(searchTree, prevNode);
510 
511 	}
512       else if(isAbove(dline, dline) && isAbove(dline, dlinePrev))
513 	{
514 //printf("case 2\n");
515 	  //insert both edges
516 	  treeNode* thisNode = TreeNodeMake(vert);
517 	  treeNode* prevNode = TreeNodeMake(vert->getPrev());
518 
519 	  vert->isKey = 1;
520           vert->keyY = keyY;
521 	  searchTree = TreeNodeInsert(searchTree, thisNode, (Int (*) (void *, void *))compChains);
522           vert->isKey = 0;
523 
524           vert->getPrev()->isKey = 1;
525           vert->getPrev()->keyY = keyY;
526 	  searchTree = TreeNodeInsert(searchTree, prevNode, (Int (*) (void *, void *))compChains);
527           vert->getPrev()->isKey = 0;
528 
529 	  if(cuspType(dline) == 1) //interior cusp
530 	    {
531 //printf("cuspType is 1\n");
532 	      treeNode* leftEdge = TreeNodePredecessor(thisNode);
533 	      treeNode* rightEdge = TreeNodeSuccessor(prevNode);
534               if(leftEdge == NULL || rightEdge == NULL)
535 		{
536 		  errOccur = 1;
537 		  goto JUMP_HERE;
538 		}
539 //printf("leftEdge is %i, rightEdge is %i\n", leftEdge, rightEdge);
540 	      directedLine* leftEdgeDline = ((monoChain*) leftEdge->key)->find(keyY);
541 	      directedLine* rightEdgeDline = ((monoChain*) rightEdge->key)->find(keyY);
542 	      ret_ranges[i] = sweepRangeMake( leftEdgeDline, 1, rightEdgeDline, 1);
543 	    }
544 	  else //exterior cusp
545 	    {
546 //printf("cuspType is not 1\n");
547 	      ret_ranges[i] = sweepRangeMake(dlinePrev, 1, dline, 1);
548 	    }
549 	}
550       else
551 	{
552 //printf("%i,%i\n", isAbove(dline, dline), isAbove(dline, dlinePrev));
553 	  errOccur = 1;
554 	  goto JUMP_HERE;
555 
556 	  fprintf(stderr, "error in MC_sweepY\n");
557 	  exit(1);
558 	}
559     }
560 
561  JUMP_HERE:
562   //finally clean up space: delete  the search tree
563   TreeNodeDeleteWholeTree(searchTree);
564   return errOccur;
565 }
566 
567 void MC_findDiagonals(Int total_num_edges, monoChain** sortedVertices,
568 		   sweepRange** ranges, Int& num_diagonals,
569 		   directedLine** diagonal_vertices)
570 {
571   Int i,j,k;
572   k=0;
573   //reset 'current' of all the monoChains
574   for(i=0; i<total_num_edges; i++)
575     sortedVertices[i]->resetCurrent();
576 
577   for(i=0; i<total_num_edges; i++)
578     {
579       directedLine* vert = sortedVertices[i]->getHead();
580       directedLine* thisEdge = vert;
581       directedLine* prevEdge = vert->getPrev();
582       if(isBelow(vert, thisEdge) && isBelow(vert, prevEdge) && compEdges(prevEdge, thisEdge)<0)
583 	{
584 	  //this is an upward interior cusp
585 	  diagonal_vertices[k++] = vert;
586 
587 	  directedLine* leftEdge = ranges[i]->left;
588 	  directedLine* rightEdge = ranges[i]->right;
589 
590 	  directedLine* leftVert = leftEdge;
591 	  directedLine* rightVert = rightEdge->getNext();
592 	  assert(leftVert->head()[1] >= vert->head()[1]);
593 	  assert(rightVert->head()[1] >= vert->head()[1]);
594 	  directedLine* minVert = (leftVert->head()[1] <= rightVert->head()[1])?leftVert:rightVert;
595 	  Int found = 0;
596 	  for(j=i+1; j<total_num_edges; j++)
597 	    {
598 	      if(sortedVertices[j]->getHead()->head()[1] > minVert->head()[1])
599 		break;
600 
601 	      if(sweepRangeEqual(ranges[i], ranges[j]))
602 		{
603 		  found = 1;
604 		  break;
605 		}
606 	    }
607 
608 	  if(found)
609 	    diagonal_vertices[k++] = sortedVertices[j]->getHead();
610 	  else
611 	    diagonal_vertices[k++] = minVert;
612 	}
613       else if(isAbove(vert, thisEdge) && isAbove(vert, prevEdge) && compEdges(prevEdge, thisEdge)>0)
614 	{
615 	  //downward interior cusp
616 	  diagonal_vertices[k++] = vert;
617 	  directedLine* leftEdge = ranges[i]->left;
618 	  directedLine* rightEdge = ranges[i]->right;
619 	  directedLine* leftVert = leftEdge->getNext();
620 	  directedLine* rightVert = rightEdge;
621 	  assert(leftVert->head()[1] <= vert->head()[1]);
622 	  assert(rightVert->head()[1] <= vert->head()[1]);
623 	  directedLine* maxVert = (leftVert->head()[1] > rightVert->head()[1])? leftVert:rightVert;
624 	  Int found=0;
625 	  for(j=i-1; j>=0; j--)
626 	    {
627 	      if(sortedVertices[j]->getHead()->head()[1] < maxVert->head()[1])
628 		break;
629 	      if(sweepRangeEqual(ranges[i], ranges[j]))
630 		{
631 		  found = 1;
632 		  break;
633 		}
634 	    }
635 	  if(found)
636 	    diagonal_vertices[k++] = sortedVertices[j]->getHead();
637 	  else
638 	    diagonal_vertices[k++] = maxVert;
639 	}
640     }
641   num_diagonals = k/2;
642 }
643 
644 
645 
646 
647 directedLine* MC_partitionY(directedLine *polygons, sampledLine **retSampledLines)
648 {
649 //printf("enter mc_partitionY\n");
650   Int total_num_chains = 0;
651   monoChain* loopList = directedLineLoopListToMonoChainLoopList(polygons);
652   monoChain** array = loopList->toArrayAllLoops(total_num_chains);
653 
654   if(total_num_chains<=2) //there is just one single monotone polygon
655     {
656       loopList->deleteLoopList();
657       free(array);
658       *retSampledLines = NULL;
659       return polygons;
660     }
661 
662 //loopList->printAllLoops();
663 //printf("total_num_chains=%i\n", total_num_chains);
664   quicksort( (void**)array, 0, total_num_chains-1, (Int (*)(void*, void*))compChainHeadInY);
665 //printf("after quicksort\n");
666 
667   sweepRange** ranges = (sweepRange**)malloc(sizeof(sweepRange*) * (total_num_chains));
668   assert(ranges);
669 
670   if(MC_sweepY(total_num_chains, array, ranges))
671     {
672       loopList->deleteLoopList();
673       free(array);
674       *retSampledLines = NULL;
675       return NULL;
676     }
677 //printf("after MC_sweepY\n");
678 
679 
680   Int num_diagonals;
681   /*number diagonals is < total_num_edges*total_num_edges*/
682   directedLine** diagonal_vertices = (directedLine**) malloc(sizeof(directedLine*) * total_num_chains*2/*total_num_edges*/);
683   assert(diagonal_vertices);
684 
685 //printf("before call MC_findDiagonales\n");
686 
687   MC_findDiagonals(total_num_chains, array, ranges, num_diagonals, diagonal_vertices);
688 //printf("after call MC_findDia, num_diagnla=%i\n", num_diagonals);
689 
690   directedLine* ret_polygons = polygons;
691   sampledLine* newSampledLines = NULL;
692   Int i,k;
693 
694   num_diagonals=deleteRepeatDiagonals(num_diagonals, diagonal_vertices, diagonal_vertices);
695 
696 
697 
698 //drawDiagonals(num_diagonals, diagonal_vertices);
699 //printf("diagoanls are \n");
700 //for(i=0; i<num_diagonals; i++)
701 //  {
702 //    printf("(%f,%f)\n", diagonal_vertices[2*i]->head()[0], diagonal_vertices[2*i]->head()[1]);
703 //    printf("**(%f,%f)\n", diagonal_vertices[2*i+1]->head()[0], diagonal_vertices[2*i+1]->head()[1]);
704 //  }
705 
706   Int *removedDiagonals=(Int*)malloc(sizeof(Int) * num_diagonals);
707   for(i=0; i<num_diagonals; i++)
708     removedDiagonals[i] = 0;
709 //  printf("first pass\n");
710 
711 
712  for(i=0,k=0; i<num_diagonals; i++,k+=2)
713     {
714 
715 
716       directedLine* v1=diagonal_vertices[k];
717       directedLine* v2=diagonal_vertices[k+1];
718       directedLine* ret_p1;
719       directedLine* ret_p2;
720 
721       /*we ahve to determine whether v1 and v2 belong to the same polygon before
722        *their structure are modified by connectDiagonal().
723        */
724 /*
725       directedLine *root1 = v1->findRoot();
726       directedLine *root2 = v2->findRoot();
727       assert(root1);
728       assert(root2);
729 */
730 
731 directedLine*  root1 = v1->rootLinkFindRoot();
732 directedLine*  root2 = v2->rootLinkFindRoot();
733 
734       if(root1 != root2)
735 	{
736 
737 	  removedDiagonals[i] = 1;
738 	  sampledLine* generatedLine;
739 
740 
741 
742 	  v1->connectDiagonal(v1,v2, &ret_p1, &ret_p2, &generatedLine, ret_polygons);
743 
744 
745 
746 	  newSampledLines = generatedLine->insert(newSampledLines);
747 /*
748 	  ret_polygons = ret_polygons->cutoffPolygon(root1);
749 
750 	  ret_polygons = ret_polygons->cutoffPolygon(root2);
751 	  ret_polygons = ret_p1->insertPolygon(ret_polygons);
752 root1->rootLinkSet(ret_p1);
753 root2->rootLinkSet(ret_p1);
754 ret_p1->rootLinkSet(NULL);
755 ret_p2->rootLinkSet(ret_p1);
756 */
757 	  ret_polygons = ret_polygons->cutoffPolygon(root2);
758 
759 
760 
761 root2->rootLinkSet(root1);
762 ret_p1->rootLinkSet(root1);
763 ret_p2->rootLinkSet(root1);
764 
765        /*now that we have connected the diagonal v1 and v2,
766         *we have to check those unprocessed diagonals which
767         *have v1 or v2 as an end point. Notice that the head of v1
768         *has the same coodinates as the head of v2->prev, and the head of
769         *v2 has the same coordinate as the head of v1->prev.
770         *Suppose these is a diagonal (v1, x). If (v1,x) is still a valid
771         *diagonal, then x should be on the left hand side of the directed line:        *v1->prev->head -- v1->head -- v1->tail. Otherwise, (v1,x) should be
772         *replaced by (v2->prev, x), that is, x is on the left of
773         * v2->prev->prev->head, v2->prev->head, v2->prev->tail.
774         */
775         Int ii, kk;
776         for(ii=0, kk=0; ii<num_diagonals; ii++, kk+=2)
777 	  if( removedDiagonals[ii]==0)
778 	    {
779 	      directedLine* d1=diagonal_vertices[kk];
780 	      directedLine* d2=diagonal_vertices[kk+1];
781 	      /*check d1, and replace diagonal_vertices[kk] if necessary*/
782 	      if(d1 == v1) {
783 		/*check if d2 is to left of v1->prev->head:v1->head:v1->tail*/
784 		if(! pointLeft2Lines(v1->getPrev()->head(),
785 				     v1->head(), v1->tail(), d2->head()))
786 		  {
787 /*
788 		    assert(pointLeft2Lines(v2->getPrev()->getPrev()->head(),
789 					   v2->getPrev()->head(),
790 					   v2->getPrev()->tail(), d2->head()));
791 */
792 		    diagonal_vertices[kk] = v2->getPrev();
793 		  }
794 	      }
795 	      if(d1 == v2) {
796 		/*check if d2 is to left of v2->prev->head:v2->head:v2->tail*/
797 		if(! pointLeft2Lines(v2->getPrev()->head(),
798 				     v2->head(), v2->tail(), d2->head()))
799 		  {
800 /*
801 		    assert(pointLeft2Lines(v1->getPrev()->getPrev()->head(),
802 					   v1->getPrev()->head(),
803 					   v1->getPrev()->tail(), d2->head()));
804 */
805 		    diagonal_vertices[kk] = v1->getPrev();
806 		  }
807 	      }
808 	      /*check d2 and replace diagonal_vertices[k+1] if necessary*/
809 	      if(d2 == v1) {
810 		/*check if d1 is to left of v1->prev->head:v1->head:v1->tail*/
811 		if(! pointLeft2Lines(v1->getPrev()->head(),
812 				     v1->head(), v1->tail(), d1->head()))
813 		  {
814 /*		    assert(pointLeft2Lines(v2->getPrev()->getPrev()->head(),
815 					   v2->getPrev()->head(),
816 					   v2->getPrev()->tail(), d1->head()));
817 */
818 		    diagonal_vertices[kk+1] = v2->getPrev();
819 		  }
820 	      }
821 	      if(d2 == v2) {
822 		/*check if d1 is to left of v2->prev->head:v2->head:v2->tail*/
823 		if(! pointLeft2Lines(v2->getPrev()->head(),
824 				     v2->head(), v2->tail(), d1->head()))
825 		  {
826 /*		    assert(pointLeft2Lines(v1->getPrev()->getPrev()->head(),
827 					   v1->getPrev()->head(),
828 					   v1->getPrev()->tail(), d1->head()));
829 */
830 		    diagonal_vertices[kk+1] = v1->getPrev();
831 		  }
832 	      }
833 	    }
834 }/*end if (root1 not equal to root 2)*/
835 }
836 
837   /*second pass,  now all diagoals should belong to the same polygon*/
838 //printf("second pass: \n");
839 
840 //  for(i=0; i<num_diagonals; i++)
841 //    printf("%i ", removedDiagonals[i]);
842 
843 
844   for(i=0,k=0; i<num_diagonals; i++, k += 2)
845     if(removedDiagonals[i] == 0)
846       {
847 
848 
849 	directedLine* v1=diagonal_vertices[k];
850 	directedLine* v2=diagonal_vertices[k+1];
851 
852 
853 
854 	directedLine* ret_p1;
855 	directedLine* ret_p2;
856 
857 	/*we ahve to determine whether v1 and v2 belong to the same polygon before
858 	 *their structure are modified by connectDiagonal().
859 	 */
860 	directedLine *root1 = v1->findRoot();
861 /*
862 	directedLine *root2 = v2->findRoot();
863 
864 
865 
866 	assert(root1);
867 	assert(root2);
868 	assert(root1 == root2);
869   */
870 	sampledLine* generatedLine;
871 
872 
873 
874 	v1->connectDiagonal(v1,v2, &ret_p1, &ret_p2, &generatedLine, ret_polygons);
875 	newSampledLines = generatedLine->insert(newSampledLines);
876 
877 	ret_polygons = ret_polygons->cutoffPolygon(root1);
878 
879 	ret_polygons = ret_p1->insertPolygon(ret_polygons);
880 
881 	ret_polygons = ret_p2->insertPolygon(ret_polygons);
882 
883 
884 
885 	for(Int j=i+1; j<num_diagonals; j++)
886 	  {
887 	    if(removedDiagonals[j] ==0)
888 	      {
889 
890 		directedLine* temp1=diagonal_vertices[2*j];
891 		directedLine* temp2=diagonal_vertices[2*j+1];
892                if(temp1==v1 || temp1==v2 || temp2==v1 || temp2==v2)
893 		if(! temp1->samePolygon(temp1, temp2))
894 		  {
895 		    /*if temp1 and temp2 are in different polygons,
896 		     *then one of them must be v1 or v2.
897 		     */
898 
899 
900 
901 		    assert(temp1==v1 || temp1 == v2 || temp2==v1 || temp2 ==v2);
902 		    if(temp1==v1)
903 		      {
904 			diagonal_vertices[2*j] = v2->getPrev();
905 		      }
906 		    if(temp2==v1)
907 		      {
908 			diagonal_vertices[2*j+1] = v2->getPrev();
909 		      }
910 		    if(temp1==v2)
911 		      {
912 			diagonal_vertices[2*j] = v1->getPrev();
913 		      }
914 		    if(temp2==v2)
915 		      {
916 			diagonal_vertices[2*j+1] = v1->getPrev();
917 		      }
918 		  }
919 	      }
920 	  }
921 
922       }
923 
924 
925   //clean up
926   loopList->deleteLoopList();
927   free(array);
928   free(ranges);
929   free(diagonal_vertices);
930   free(removedDiagonals);
931 
932   *retSampledLines = newSampledLines;
933   return ret_polygons;
934 }
935 
936 
937