0024914: Micro edge is created during Boolean Operations
[occt.git] / src / Graphic3d / Graphic3d_Strips.cxx
1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
2 //
3 // This file is part of Open CASCADE Technology software library.
4 //
5 // This library is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU Lesser General Public License version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
10 //
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
13
14 #include <Graphic3d_Strips.ixx>
15
16 /*
17                                 TRIANGLE_STRIP
18 */
19
20
21 /*
22       Algorithm to find Triangle-Strips in a triangles graph.
23
24   A triangle graph is a set of vertices (numbered from 1 to nvertices)
25 and of triangles (numbered from 1  to ntriangles), each  triangle is a
26 triplet of vertices.
27
28   A  triangle-strip  is  a  strip of  adjacent  triangles that can  be
29 described as a   list of vertices.  The  strip   v1,v2,v3,v4,v5,v6,...
30 describes the triangles (v1,v2,v3), (v2,v3,v4), (v3,v4,v5), (v4,v5,v6)
31 ...
32
33   Triangle-strips  are an economic  way to pass  around triangles to a
34 hardware shader as they require an average of one vertex per triangle.
35
36   The purpose  of  this algorithm is to  break a triangle graph into a
37 (minimal cardinality) set of triangle strips.
38
39 It is purely topological, triangles are triplets of long integers.
40 There is no limit of size, memory is allocated from free store.
41 The quantity allocated during the algorithm is about
42   66 * t bytes where t is the number of triangles.
43
44 ( the description of the algorithm can be found at the end of the file
45   starting with  "COMMENTS")
46 */
47
48
49 /*******************************/
50 /*  GLOBAL TYPES AND VARIABLES */
51 /*******************************/
52
53 #include <stddef.h>
54 #include <Standard.hxx>
55
56 #define NULLTRIANGLE 0
57 #define DELETED 0
58
59 /* the array of triangles is the basic data structure to folow strips
60    it is allocated from free store */
61 typedef struct {
62   int  v[3];  /* the three vertices of the triangle */
63   int  tn[3]; /* neighbouring triangles */
64   int ivn[3]; /* the index of the neigbouring vertex */
65   int  state; /* the last strip crossing the triangle, 0 means deleted */
66 } triangle;
67
68 triangle *trianglesptr;
69 int TrianglesPtrSize;
70
71 /* remarks about the above structure :
72  Triangles are ranging from 1 two nbtriangles, triangle 0 will always
73  be deleted.
74  A state of 0 means the triangle is deleted from the graph.
75  The vertices are v[0],v[1],v[2]
76  To get the neigbour of the triangle on the other side of the edge
77  v[i],v[j] just pick tn[i+j-1] and ivn[i+j-1].
78  If tn[i+j-1] is 0 there is no neighbour.
79  ivn is the index (0,1,2) of the vertex of the neighbouring triangle
80  tn[] which is not shared with this triangle.
81 */
82
83 /* the number of triangles */
84 int  nbtriangles;
85
86 /* a strip is described by a triangle and the index of the two last
87   NBVERTICES of the triangle in the strip */
88 typedef struct {
89   int  t;     /* the triangle */
90   int iv1,iv2; /* the last NBVERTICES of t in the strip, in 0,1,2 */
91 } stript;
92
93 /* the current strip position is saved in this variable */
94 /* between calls to to GET_VERTEX */
95 static stript current_stript;
96
97 /* this index is used to label the last strip under exploration */
98 static int  last_stript;
99
100 /* tell this dumb compiler that stript_next returns  nothing */
101 void stript_next(stript *st);
102 int  stript_score(stript* pstrip, int *plength);
103
104 /********************************************************************/
105 /*                                                                  */
106 /* STRIPT_INIT : get data and build the triangles array             */
107 /*                                                                  */
108 /********************************************************************/
109
110 void Graphic3d_Strips :: STRIPT_INIT ( const Standard_Integer NBVERTICES,
111                                                                            const TColStd_Array1OfInteger& TABTRIANGLES )  
112 {
113  /* In  order to build  the triangles array  we  will  use a temporary
114 array  : edges. This array  is of  length NBVERTICES.  Each entry is a
115 pointer  to a  list of  structures  of  the  type edge. This structure
116 describes  an edge by : the  second vertex  of the  edge   and the two
117 triangles adjacent to the edge, the starting vertex of the edge is the
118 entry of the array edges. The smallest vertex index of an edge is used
119 to index it in the edges array */
120
121
122   int NBTRIANG = (int) TABTRIANGLES.Length() / 3;
123
124   typedef struct edg {
125     struct edg *next; /* next edge in the list for a vertex */
126     int  v;           /* the second vertex of the edge */
127     int  tn[2];       /* neighbour triangles */
128     int ivn[2];       /* index of third vertex of neighbour triangles */
129   } edge;
130
131   edge **edges;
132   edge *cedge;
133   int  ivert,triang;
134   int  vmin,vmax;
135   int ivthird;
136   int TedgesSize;
137
138   int i,j;
139
140   /* copy the number and initialize a few */
141   nbtriangles = NBTRIANG;
142   last_stript = 1;
143
144   /* allocate the array edges 
145      vertices are ranging from 1 to NBVERTICES */
146   TedgesSize = (NBVERTICES+1) * sizeof(edge*);
147   edges = (edge**) Standard::Allocate(TedgesSize);
148   for (ivert=0;ivert<= NBVERTICES; ivert++) {
149     edges[ivert] = NULL;
150   }
151
152   /* allocate the array triangles from 0 to nbtriangles */
153   TrianglesPtrSize = (nbtriangles+1)*sizeof(triangle);
154   trianglesptr = (triangle*) Standard::Allocate (TrianglesPtrSize);
155   trianglesptr[0].state = DELETED;
156   trianglesptr[0].tn[0] = NULLTRIANGLE;
157   trianglesptr[0].tn[1] = NULLTRIANGLE;
158   trianglesptr[0].tn[2] = NULLTRIANGLE;
159   trianglesptr[0].ivn[0] = 0;
160   trianglesptr[0].ivn[1] = 0;
161   trianglesptr[0].ivn[2] = 0;
162
163
164   /* copy the triangles into the arrays */
165   for (triang=1;triang<=nbtriangles;triang++) {
166
167     /* copy the vertices */
168     trianglesptr[triang].state = 1;
169     for (j=0;j<3;j++) 
170           trianglesptr[triang].v[j] = TABTRIANGLES(3*(triang-1)+j+1);
171       
172     /* insert the edges in the edges array */
173     for (j=0;j<3;j++) {
174       if (trianglesptr[triang].v[j] <= trianglesptr[triang].v[(j+1)%3]) 
175           {
176                 vmin = trianglesptr[triang].v[j];
177                 vmax = trianglesptr[triang].v[(j+1)%3];
178       }
179       else 
180           {
181                 vmax = trianglesptr[triang].v[j];
182                 vmin = trianglesptr[triang].v[(j+1)%3];
183       }
184       ivthird = (j+2)%3;
185
186       /* the edge is inserted in the array at the entry for the 
187          smallest vertex */
188
189       /* first search if there is an entry for this edge */
190       cedge = edges[vmin];
191       while(cedge != NULL) {
192         if (cedge->v == vmax) 
193           break;
194         cedge = cedge->next;
195       }
196       /* if the edge was not found, create it */
197       if (cedge == NULL) {
198         cedge = (edge*) Standard::Allocate (sizeof(edge));
199         cedge->next = edges[vmin];
200         edges[vmin] = cedge;
201         cedge->v = vmax;
202         cedge->tn[0] = triang;
203         cedge->ivn[0] = ivthird;
204         cedge->tn[1] = 0;
205         cedge->ivn[1] = 0;
206       }
207       else {
208         cedge->tn[1] = triang;
209         cedge->ivn[1] = ivthird;
210       }
211     }
212   }
213
214   /* now complete the triangles array (neighbours) using the edges */
215   /* array */
216
217   for (triang=1;triang<=nbtriangles;triang++) {
218     /* on each edge of the triangle : find the neighbour */
219     for (j=0;j<3;j++) {
220       if (trianglesptr[triang].v[j] <= trianglesptr[triang].v[(j+1)%3]) {
221         vmin = trianglesptr[triang].v[j];
222         vmax = trianglesptr[triang].v[(j+1)%3];
223       }
224       else {
225         vmax = trianglesptr[triang].v[j];
226         vmin = trianglesptr[triang].v[(j+1)%3];
227       }
228
229       /* search the entry for the edge */
230       cedge = edges[vmin];
231       while(cedge->v != vmax) {
232         cedge = cedge->next;
233       }
234
235       /* find the neighbouring triangle */
236       i = 0;
237       if (cedge->tn[0] == triang) i = 1;
238       trianglesptr[triang].tn[(2*j)%3] = cedge->tn[i];
239       trianglesptr[triang].ivn[(2*j)%3] = cedge->ivn[i];
240     }
241   }
242
243   /* destroy the edges array which has done it's duty */
244   for (ivert = 1; ivert <= NBVERTICES; ivert++) {
245     while(edges[ivert] != NULL) {
246       cedge = edges[ivert];
247       edges[ivert] = cedge->next;
248       Standard::Free(cedge);
249     }
250   }
251   Standard::Free(edges);
252 }
253
254
255 /********************************************************************/
256 /*                                                                  */
257 /* STRIPT_GET_STRIP : find the next strip                           */
258 /*                                                                  */
259 /********************************************************************/
260
261 void Graphic3d_Strips :: STRIPT_GET_STRIP ( Standard_Integer& NBTRIANGLES,
262                                                                                     Standard_Integer& V1,
263                                                                                         Standard_Integer& V2 )
264 {
265
266   int  btriang;   /* the triangle with the lowest number of neigbours */
267   int  triang;
268   int  tr;
269   int bneib,neib;
270   stript cstrip;   /* the current strip */
271   int   cscore;    /* it's score */
272   int  cleng;      /* it's length */
273   /* the best strip is stored in current_strip */
274   int  blength; /* the best strip length */
275   int  bscore;  /* the best strip score */
276   
277   int i;
278
279   /* first find the triangle with the lowest number of neighbours */
280   btriang = 0;
281   bneib = 4;
282   for (triang=1; triang<=nbtriangles; triang++) {
283     if (trianglesptr[triang].state != 0) {
284       neib = 0;
285       for (i=0;i<3;i++) {
286         tr = trianglesptr[triang].tn[i];
287         if ((tr != 0) && (trianglesptr[tr].state != 0)) {
288           neib++;
289         }
290       }
291       if (neib < bneib) {
292         bneib = neib;
293         btriang = triang;
294         /* a triangle with 0 or one neighbours is fine */
295         if (neib <= 1) break;
296       }
297     }
298   }
299
300   /* if none was found stop the process and free the memory */
301   if (btriang == 0) 
302   {
303     NBTRIANGLES = 0;
304     current_stript.t = 0;
305     Standard::Free(trianglesptr);
306     return;
307   }
308
309   /* now search the best strip from this triangle 
310      the strip with the biggest score.
311      If score are even the biggest length win */
312
313   /* try 0,1,2 */
314   current_stript.t = btriang;
315   current_stript.iv1 = 1;
316   current_stript.iv2 = 2;
317   bscore = stript_score(&current_stript,&blength);
318
319   /* try 1,2,0 */
320   cstrip.t = btriang;
321   cstrip.iv1 = 2;
322   cstrip.iv2 = 0;
323   cscore = stript_score(&cstrip,&cleng);
324   if ((cscore > bscore)  ||
325       ((cscore == bscore) && (cleng > blength))){
326     bscore = cscore;
327     blength = cleng;
328     current_stript.t = cstrip.t;
329     current_stript.iv1 = cstrip.iv1;
330     current_stript.iv2 = cstrip.iv2;
331   }
332
333   /* try 2,0,1 */
334   cstrip.t = btriang;
335   cstrip.iv1 = 0;
336   cstrip.iv2 = 1;
337   cscore = stript_score(&cstrip,&cleng);
338   if ((cscore > bscore) ||
339       ((cscore == bscore) && (cleng > blength))){
340     bscore = cscore;
341     blength = cleng;
342     current_stript.t = cstrip.t;
343     current_stript.iv1 = cstrip.iv1;
344     current_stript.iv2 = cstrip.iv2;
345   }
346
347   /* return the best strip */
348   NBTRIANGLES = blength;
349   triang = current_stript.t;
350   V2 = trianglesptr[triang].v[current_stript.iv1];
351   V1 = trianglesptr[triang].v[3-current_stript.iv1-current_stript.iv2];
352   return;
353 }
354
355
356 /********************************************************************/
357 /*                                                                  */
358 /* STRIPT_GET_VERTEX : get next vertex & triangle in current strip   */
359 /*                                                                  */
360 /********************************************************************/
361
362 void Graphic3d_Strips :: STRIPT_GET_VERTEX ( Standard_Integer& VERTEX,
363                                                                                          Standard_Integer& TRIANGLE )  
364 {
365   int  triang;
366   triang = current_stript.t;
367   /* delete this triangle */
368   trianglesptr[triang].state = 0;
369   TRIANGLE = triang;
370   VERTEX = trianglesptr[triang].v[current_stript.iv2];
371   stript_next(&current_stript);
372   return;
373 }
374
375
376
377 /********************************************************************/
378 /*                                                                  */
379 /* stript_score  : find the start of a strip and it's lenght         */
380 /*                returns the score of the strip                    */
381 /*                                                                  */
382 /********************************************************************/
383
384 int stript_score(stript* pstrip, int *plength)
385 {
386 /* st is set to the beginning of the strip and the length of the strip 
387    is returned. The strip is explored in two directions, if it loops
388    on itself it is detected. */
389 /* the score is a value to optimise. The number of boundary triangles */
390  /* in a strip seems to be a nice choice. */
391   stript cstrip,savstrip;
392   int  length;
393   int  score;
394   int i;
395   int  triang;
396
397   length = 0;
398   score = 0;
399   last_stript++; /* this is used to mark triangles in this strip */
400
401   /* go in the first direction */
402   cstrip.t = pstrip->t;
403   cstrip.iv1 = pstrip->iv1;
404   cstrip.iv2 = pstrip->iv2;
405   while ((cstrip.t != 0) &&                           /* - on a boundary */
406          (trianglesptr[cstrip.t].state != 0) &&          /* - deleted */
407          (trianglesptr[cstrip.t].state != last_stript)) { /* - on the same */
408                                                       /* strip */
409     trianglesptr[cstrip.t].state = last_stript;
410     /* increment the length */
411     length++;
412     /* compute the score */
413     /* increment the score if the triangle has less than three */
414     /* neigbours */
415     for (i=0;i<3;i++)  {
416       triang = trianglesptr[cstrip.t].tn[i];
417       if ((triang == 0) || (trianglesptr[triang].state == 0)) {
418         score++;
419         break;
420       }
421     }
422     /* next in the strip */
423     stript_next(&cstrip);
424   }
425   /* go in the reversed direction */
426   cstrip.t = pstrip->t;
427   cstrip.iv1 = pstrip->iv1;
428   cstrip.iv2 = 3 - pstrip->iv2 - pstrip->iv1;
429   /* save the position of the strip before moving  */
430   savstrip.t = cstrip.t;
431   savstrip.iv1 = cstrip.iv1;
432   savstrip.iv2 = cstrip.iv2;
433   stript_next(&cstrip);
434   while ((cstrip.t != 0) && 
435          (trianglesptr[cstrip.t].state != 0) &&
436          (trianglesptr[cstrip.t].state != last_stript)) {
437     trianglesptr[cstrip.t].state = last_stript;
438     /* save the position of the strip before moving  */
439     savstrip.t = cstrip.t;
440     savstrip.iv1 = cstrip.iv1;
441     savstrip.iv2 = cstrip.iv2;
442     /* increment the length */
443     length++;
444     /* compute the score */
445     /* increment the score if the triangle has less than three */
446     /* neigbours */
447     for (i=0;i<3;i++)  {
448       triang = trianglesptr[cstrip.t].tn[i];
449       if ((triang == 0) || (trianglesptr[triang].state == 0)) {
450         score++;
451         break;
452       }
453     }
454     /* next in the strip */
455     stript_next(&cstrip);
456   }
457
458   /* reverse in the good direction the saved position */
459   pstrip->t = savstrip.t;
460   pstrip->iv1 = savstrip.iv1;
461   pstrip->iv2 = 3 - savstrip.iv1 - savstrip.iv2;
462   *plength = length;
463   return score;
464 }
465
466 /********************************************************************/
467 /*                                                                  */
468 /* stript_next : jump to next triangle in a strip                    */
469 /*                                                                  */
470 /********************************************************************/
471
472 void stript_next(stript *st) 
473 {
474 /* st points toward a triangle and a vertex ordering defining a unique */
475 /* it's content is changed for the next triangle in the strip */
476 /* the triangle may becomes 0 if there was no neighbour */
477
478   int  triang,ntriang;
479   int i,j;
480
481   triang = st->t;
482
483   if (triang == 0) {
484     st->t = 0;
485     st->iv1 = 0;
486     st->iv2 = 0;
487     return;
488   }
489   /* get the neighbouring triangle */
490   i = st->iv1+st->iv2-1;
491   ntriang = trianglesptr[triang].tn[i];
492   /* if there is no neighbour */
493   if (ntriang == 0) {
494     st->t = 0;
495     st->iv1 = 0;
496     st->iv2 = 0;
497     return;
498   }
499   /* compute the new index for the last vertex */
500   j = 0;
501   while(trianglesptr[triang].v[st->iv2] != trianglesptr[ntriang].v[j]) {
502     j++;
503   }
504   st->t = ntriang;
505   st->iv1 = j;
506   st->iv2 = trianglesptr[triang].ivn[i];
507
508   return;
509 }
510
511 /*******************************************************************/
512 /*********                                                **********/
513 /*********              COMMENTS                          **********/
514 /*********                                                **********/
515 /*******************************************************************/
516
517
518 /* 
519   Architecture
520   ************
521
522 The present C implementation was designed to be called from FORTRAN
523 with the following syntaxes :
524
525 1. STRIP_INIT(int * NBVERTICES,int * NBTRIANGLES,int * TABTRIANGLES)
526
527    This is the initiating call where :
528    NBVERTICES is the number of vertices
529    NBTRIANGLES is the number of triangles
530    TABTRIANGLES is the table describing the triangles
531
532    This  function copies  the arguments to  nvertices, ntriangles, and
533 build an inner table of triangles (triangles) were the neighbours of a
534 triangle can be found easily.
535
536
537 IMPORTANT WARNINGS
538   Vertices and Triangles are in the range 1...NBVERTICES
539   1...NBTRIANGLES
540   Double arrays are FORTRAN arrays so the three vertices of triangle I
541   are found in TABTRIANGLES[I+J-1] where J = 0,1,2
542   The FORTRAN array is declared as TABTRIANGLES(3,NBTRIANGLES)
543
544
545   2.  STRIP_GET_STRIP(int * NBTRIANGLES,int * V1,int * V2)
546       STRIP_GET_VERTEX(int * VERTEX,int * TRIANGLE)
547
548       Both functions are used to get the strips. each iteration of the
549      GET_STRIP brings a new strip wre NBTRIANGLES is the number of triangles
550      in the strip (i.e the number  of  vertices is NBTRIANGLES plus two) and  V1, V2
551      are the  two first vertices.  NBTRIANGLES becomes   zero when the last
552      strip has been read.  GET_VERTEX is  used two get the successives
553      vertices of a strip, each  vertex is associated with  a triangle.
554      This start with  the third vertex, the two   first  are given  by
555      GET_STRIP.
556
557      An example of correct call from C to read the strips is
558
559      while(1) {
560        STRIP_GET_STRIP(&nb_triangles,&v[1],&v[2]);
561        if (nb_triangles == 0) break;
562        for (i=3;i <= ( nb_triangles+2 ); i++) {
563          STRIP_GET_VERTEX(&v[i],&t[i]);
564          }
565        }
566
567        Unwise calls to those functions will generally return zero but
568        are not recommanded.
569
570        The data are not checked for coherence, but a zero
571        number of triangle will give a zero number of strips.
572       
573 */
574
575 /*
576 OUTLINE OF THE ALGORITHM
577 ************************
578
579 The algorithm is purely topological.  No coordinates  are given,  it's
580 input are the number  of  vertices, the number  of triangles,  and for
581 each triangle a triplet of vertices indices.
582
583 Let us consider a   triangle T =   (V1,  V2, V3),  this  triangle  has
584 neighbours in the triangle graphs, let us call them T12, T23, T31. T12
585 is the triangle sharing the edge V1,V2 with T, etc...  Of course those
586 three triangles may not exist  in the graph in  this case T is "on the
587 border" or even "in a corner".
588
589 The key remark is that at most three  triangle-strips may cross T they
590 are : T12-T-T23, T23-T-T31,  T31-T-T12. Once three adjacent  triangles
591 are given the entire strip  is uniquely  defined, the orientation of a
592 strip is meaningless as if you reverse it you get the same strip.
593 To describe a current position in a strip you need a triangle and the
594 two last vertices of the triangle in the strip.
595
596 Our algorithm (more precisely heuristic) is the following : 
597 - Find a triangle with the lowest number of neighbours.
598 - List the three (or less) strips crossing this triangle.
599 - Chose the best among them and remove from the graph all the
600 triangles in this strip.
601   The best strip is rated with a "score", the score we used is the
602   number of triangles in the strip which have less than three
603   neighbours  (they are on the "border") in case of score equality the
604   longest strip is selected.
605 - Reiterate the process until there are no triangles left.
606
607 There are no demonstrations of  the optimality  of this algorithm, but
608 it seems to give expected results on regular graphs which are the most
609 commonly fed to it. On a rectangular array of squares, each square cut
610 in  two triangles, it  will generate strips parallels  to  the longest
611 side of the rectangle.
612
613 */
614
615 /*
616 Implementation
617 **************
618
619 First the STRIP_INIT function stores the triangles in a data structure
620 (the triangles array   allocated on free  store),  containing for each
621 triangle it's three neighbours and the third vertex for each neighbour
622 (a zero neighbour is inexistant), a triangle get's also a state 1.  To
623 build this  array a temporary  array "edges" is build giving  for each
624 edge (pair of vertices) the two neghbouring triangles.
625
626 Then at  each  call    of STRIP_GET_STRIP  a  triangle    with minimum
627 neighbours is first chosen.  For   the three possible strips  crossing
628 the triangle the strip_score function is called  which brings back the
629 start of the strip, the length and the score.  The strip-next function
630 is used to jump to  the next triangle  in a strip.   The best strip is
631 chosen, stored in the current_strip  and returned. 
632
633
634 Each   call  to   STRIP_GET_VERTEX  increment  the   the current-strip
635 structure  to the next triangle in   the  strip,  using the strip_next
636 function. The triangle is deleted from the graph and returned.
637
638 The last call  to STRIP_GET_STRIP returns  the  triangles array to the
639 free store.
640
641
642 */
643
644
645 /*
646                                 QUADRANGLE_STRIP
647 */
648
649 /*
650       Algorithm to find Quadragle-Strips in a quadrangles graph.
651
652   A quadrangle graph is a set of vertices (numbered from 1 to nbvertices)
653 and of quadrangles (numbered from 1  to nbquadrangles), each  quadrangle is a
654 quadruplet of vertices.
655
656   A quadrangle-strip is a strip of adjacent quadrangles that can be
657 described as a list of vertices.
658   The  strip v1, v2, v3, v4, v5, v6, v7, v8 ...
659 describes quadrangles (v1, v2, v4, v3), (v4, v3, v5, v6), (v5, v6, v8, v7) ...
660         1-3-5-7
661         | | | |
662         2-4-6-8
663
664   Quadrangle-strips  are an economic  way to pass quadrangles to a
665 hardware renderer as they require an average of two vertex per quadrangle.
666
667   The purpose  of  this algorithm is to  break a quadrangle graph into a
668 (minimal cardinality) set of quadrangle strips.
669
670 It is purely topological, quadrangles are quadruplets of integers.
671 There is no limit of size, memory is allocated from free store.
672 The quantity allocated during the algorithm is about
673   (17*sizeof(int)+align)*q bytes where q is the number of quadrangles
674 and align is system-dependent alignment.
675
676 ( the description of the algorithm can be found at the end of the file
677   starting with  "COMMENTS")
678 */
679
680
681 /*******************************/
682 /*  GLOBAL TYPES AND VARIABLES */
683 /*******************************/
684
685 #define NULLQUADRANGLE 0
686
687 /* the array of quadrangles is the basic data structure to follow strips
688    it is allocated from free store */
689 typedef struct {
690   int  v[4];  /* the four vertices of the quadrangle */
691   int  qn[4]; /* neighbouring quadrangles */
692   int ivn[4][2]; /* the index of two neighbouring vertice [q][v]*/
693   int  state; /* the last strip crossing the quadrangle, 0 means deleted */
694 } quadrangle;
695
696 quadrangle *quadranglesptr;
697 int QuadranglesPtrSize;
698
699 /* the number of quadrangles */
700 int  nbquadrangles;
701
702 /* remarks about the above structure :
703  Quadrangles are ranging from 1 two nbquadrangles, quadrangle 0 will always
704  be deleted.
705  A state of 0 means the quadrangle is deleted from the graph.
706  The NBVERTICES are v[0], v[1], v[2], v[3].
707  To get the neigbour of the quadrangle on the other side of the edge
708  v[i], v[j] just pick qn[i+j-1] and ivn[i+j-1][0], ivn[i+j-1][1].
709  If qn[i+j-1] is 0 there is no neighbour.
710  ivn is the index (0, 1, 2, 3)(0, 1) of two vertice of the neighbouring
711  quadrangle qn[] which are not shared with this quadrangle.
712 */
713
714 /* a strip is described by a quadrangle and the index of the two last
715   NBVERTICES of the quadrangle in the strip */
716 typedef struct {
717   int  q;     /* the quadrangle */
718   int iv2, iv3; /* the last NBVERTICES of q in the strip, in (0, 1, 2, 3) */
719 } stripq;
720
721 /* the current strip position is saved in this variable */
722 /* between calls to to STRIPQ_GET_NEXT */
723  static stripq current_stripq;
724
725 /* this index is used to label the last strip under exploration */
726 static int  last_stripq;
727
728 /* tell this dumb compiler that stripq_next returns  nothing */
729 void stripq_next(stripq *st);
730 int stripq_score(stripq *pstrip, int *plength);
731
732 /********************************************************************/
733 /*                                                                  */
734 /* STRIPQ_INIT : get data and build the quadrangles array           */
735 /*                                                                  */
736 /********************************************************************/
737
738 void Graphic3d_Strips :: STRIPQ_INIT ( const Standard_Integer NBNBVERTICES,
739                                                                            const Standard_Integer NBQUADRANG,
740                                                                            const TColStd_SequenceOfInteger& TABQUADRANGLES )  
741 {
742
743  /* In  order to build  the quadrangles array  we  will  use a temporary
744 array: edges. This array  is of  length NBNBVERTICES.  Each entry is a
745 pointer  to a  list of  structures  of  the  type edge. This structure
746 describes  an edge by: the  second vertex  of the  edge   and the two
747 quadrangles adjacent to the edge, the starting vertex of the edge is the
748 entry of the array edges. The smallest vertex index of an edge is used
749 to index it in the edges array */
750
751   typedef struct edg {
752     struct edg *next; /* next edge in the list for a vertex */
753     int  v;           /* the second vertex of the edge */
754     int  qn[2];       /* neighbour quadrangles */
755     int ivn[2][2];    /* index of two vertice of neighbour quadrangles [q][v]*/
756   } edge;
757
758   edge **edges;
759   edge *cedge;
760   int  ivert, quadrang;
761   int  vmin, vmax;
762   int  iv3, iv4;
763   int QedgesSize;
764
765   int i, j;
766
767   /* copy the number and initialize a few */
768   nbquadrangles = NBQUADRANG;
769   last_stripq = 1;
770
771   /* allocate the array edges 
772      NBVERTICES are ranging from 1 to NBNBVERTICES */
773   QedgesSize = (NBNBVERTICES+1) * sizeof(edge*);
774   edges = (edge**) Standard::Allocate (QedgesSize);
775   for (ivert=0; ivert<= NBNBVERTICES; ivert++) {
776     edges[ivert] = NULL;
777   }
778
779   /* allocate the array quadrangles from 0 to nbquadrangles */
780   QuadranglesPtrSize = (nbquadrangles+1)*sizeof(quadrangle);
781   quadranglesptr = (quadrangle*) Standard::Allocate (QuadranglesPtrSize);
782   quadranglesptr[0].v[0] = 0;
783   quadranglesptr[0].v[1] = 0;
784   quadranglesptr[0].v[2] = 0;
785   quadranglesptr[0].v[3] = 0;
786   quadranglesptr[0].qn[0] = NULLQUADRANGLE;
787   quadranglesptr[0].qn[1] = NULLQUADRANGLE;
788   quadranglesptr[0].qn[2] = NULLQUADRANGLE;
789   quadranglesptr[0].qn[3] = NULLQUADRANGLE;
790   quadranglesptr[0].ivn[0][0] = 0;
791   quadranglesptr[0].ivn[0][1] = 0;
792   quadranglesptr[0].ivn[1][0] = 0;
793   quadranglesptr[0].ivn[1][1] = 0;
794   quadranglesptr[0].ivn[2][0] = 0;
795   quadranglesptr[0].ivn[2][1] = 0;
796   quadranglesptr[0].ivn[3][0] = 0;
797   quadranglesptr[0].ivn[3][1] = 0;
798   quadranglesptr[0].state = DELETED;
799
800   /* copy the quadrangles into the arrays */
801   for (quadrang=1; quadrang<=nbquadrangles; quadrang++)
802   {
803     /* copy the NBVERTICES */
804     quadranglesptr[quadrang].state = 1;
805     for (j=0; j<4; j++)
806       quadranglesptr[quadrang].v[j] = TABQUADRANGLES(4*(quadrang-1)+j+1);
807     /* insert the edges in the edges array */
808     for (j=0; j<4; j++)
809     {
810       if (quadranglesptr[quadrang].v[j] <= quadranglesptr[quadrang].v[(j+1)%4])
811       {
812         vmin = quadranglesptr[quadrang].v[j];
813         vmax = quadranglesptr[quadrang].v[(j+1)%4];
814       }
815       else
816       {
817         vmax = quadranglesptr[quadrang].v[j];
818         vmin = quadranglesptr[quadrang].v[(j+1)%4];
819       }
820       iv3 = (j+2)%4;
821       iv4 = (j+3)%4;
822       /* the edge is inserted in the array at the entry for the 
823          smallest vertex */
824       /* first search if there is an entry for this edge */
825       cedge = edges[vmin];
826       while(cedge != NULL)
827       {
828         if (cedge->v == vmax) 
829           break;
830         cedge = cedge->next;
831       }
832       /* if the edge was not found, create it */
833       if (cedge == NULL) {
834         cedge = (edge*) Standard::Allocate (sizeof(edge));
835         cedge->next = edges[vmin];
836         edges[vmin] = cedge;
837         cedge->v = vmax;
838         cedge->qn[0] = quadrang;
839         cedge->ivn[0][0] = iv3;
840         cedge->ivn[0][1] = iv4;
841         cedge->qn[1] = 0;
842         cedge->ivn[1][0] = 0;
843         cedge->ivn[1][1] = 0;
844       }
845       else
846       {
847         cedge->qn[1] = quadrang;
848         cedge->ivn[1][0] = iv3;
849         cedge->ivn[1][1] = iv4;
850       }
851     }
852   }
853   /* now complete the quadrangles array (neighbours) using the edges array */
854   for (quadrang=1; quadrang<=nbquadrangles; quadrang++)
855   {
856     /* on each edge of the quadrangle: find the neighbour */
857     for (j=0; j<4; j++)
858     {
859       if (quadranglesptr[quadrang].v[j] <= quadranglesptr[quadrang].v[(j+1)%4])
860       {
861         vmin = quadranglesptr[quadrang].v[j];
862         vmax = quadranglesptr[quadrang].v[(j+1)%4];
863       }
864       else
865       {
866         vmax = quadranglesptr[quadrang].v[j];
867         vmin = quadranglesptr[quadrang].v[(j+1)%4];
868       }
869       /* search the entry for the edge */
870       cedge = edges[vmin];
871       while(cedge->v != vmax)
872         cedge = cedge->next;
873       /* find the neighbouring quadrangle */
874       i = 0;
875       if (cedge->qn[0] == quadrang)
876         i = 1;
877       quadranglesptr[quadrang].qn[j] = cedge->qn[i];
878       quadranglesptr[quadrang].ivn[j][0] = cedge->ivn[i][0];
879       quadranglesptr[quadrang].ivn[j][1] = cedge->ivn[i][1];
880     }
881   }
882   /* destroy the edges array which has done it's duty */
883   for (ivert = 1; ivert <= NBNBVERTICES; ivert++)
884   {
885     while(edges[ivert] != NULL)
886     {
887       cedge = edges[ivert];
888       edges[ivert] = cedge->next;
889       Standard::Free(cedge);
890     }
891   }
892   Standard::Free(edges);
893 }
894
895
896 /********************************************************************/
897 /*                                                                  */
898 /* STRIPQ_GET_STRIP : find the next strip                           */
899 /*                                                                  */
900 /********************************************************************/
901
902 void Graphic3d_Strips :: STRIPQ_GET_STRIP ( Standard_Integer& NBQUAD,Standard_Integer& V1,
903                                                                                     Standard_Integer& V2 )
904 {
905   int  bquadrang; /* the quadrangle with the lowest number of neigbours */
906   int  quadrang;
907   int  quad;
908   int  bneib, neib;
909   stripq cstrip;   /* the current strip */
910   int  cscore;    /* it's score */
911   int  cleng;     /* it's length */
912   /* the best strip is stored in current_strip */
913   int  blength;   /* the best strip length */
914   int  bscore;    /* the best strip score */
915   int  i;
916
917   /* first find the quadrangle with the lowest number of neighbours */
918   bquadrang = 0;
919   bneib = 5;
920   for (quadrang=1; quadrang<=nbquadrangles; quadrang++)
921   {
922     if (quadranglesptr[quadrang].state != 0)
923     {
924       neib = 0;
925       for (i=0; i<4; i++)
926       {
927         quad = quadranglesptr[quadrang].qn[i];
928         if ((quad != 0) && (quadranglesptr[quad].state != 0))
929           neib++;
930       }
931       if (neib < bneib)
932       {
933         bneib = neib;
934         bquadrang = quadrang;
935         /* a quadrangle with 0 or one neighbours is fine */
936         if (neib <= 1)
937           break;
938       }
939     }
940   }
941   /* if none was found stop the process and free the memory */
942   if (bquadrang == 0)
943   {
944     NBQUAD = 0;
945     current_stripq.q = 0;
946     Standard::Free(quadranglesptr);
947     return;
948   }
949   /* Now search the best strip from this quadrangle 
950      the strip with the biggest score.
951      If score were even the biggest length win. */
952   /* try 0, 1, 2, 3 */
953   current_stripq.q = bquadrang;
954   current_stripq.iv2 = 2;
955   current_stripq.iv3 = 3;
956   bscore = stripq_score(&current_stripq, &blength);
957   /* try 1, 2, 3, 0 */
958   cstrip.q = bquadrang;
959   cstrip.iv2 = 3;
960   cstrip.iv3 = 0;
961   cscore = stripq_score(&cstrip, &cleng);
962   if ((cscore > bscore)  ||
963       ((cscore == bscore) && (cleng > blength)))
964   {
965     bscore = cscore;
966     blength = cleng;
967     current_stripq.q = cstrip.q;
968     current_stripq.iv2 = cstrip.iv2;
969     current_stripq.iv3 = cstrip.iv3;
970   }
971   /* return the best strip */
972   NBQUAD = blength;
973   quadrang = current_stripq.q;
974
975   V1 = quadranglesptr[quadrang].v[(current_stripq.iv2+2)%4];
976   V2 = quadranglesptr[quadrang].v[(current_stripq.iv3+2)%4];
977   return;
978 }
979
980 /********************************************************************/
981 /*                                                                  */
982 /* STRIPQ_GET_NEXT : get next vertex & quadrangle in current strip  */
983 /*                                                                  */
984 /********************************************************************/
985 void Graphic3d_Strips :: STRIPQ_GET_NEXT ( Standard_Integer& VERTEX1,
986                                                                                    Standard_Integer& VERTEX2,
987                                                                                    Standard_Integer& QUADRANGLE )  
988 {
989   int   quadrang = current_stripq.q;
990   /* delete this quadrangle */
991   quadranglesptr[quadrang].state = 0;
992   QUADRANGLE = quadrang;
993   /* reversed */
994   VERTEX2 = quadranglesptr[quadrang].v[current_stripq.iv2];
995   VERTEX1 = quadranglesptr[quadrang].v[current_stripq.iv3];
996   stripq_next(&current_stripq);
997   return;
998 }
999
1000 /********************************************************************/
1001 /*                                                                  */
1002 /* stripq_score  : find the start of a strip and it's length         */
1003 /*                returns the score of the strip                    */
1004 /*                                                                  */
1005 /********************************************************************/
1006 int stripq_score(stripq *pstrip, int *plength)
1007 {
1008 /* st is set to the beginning of the strip and the length of the strip 
1009    is returned. The strip is explored in two directions, if it loops
1010    on itself it is detected. */
1011 /* The score is a value to optimise. The number of boundary quadrangles */
1012 /* in a strip seems to be a nice choice. */
1013   stripq        cstrip, savstrip;
1014   int   length;
1015   int   score;
1016   int   i;
1017   int   quadrang;
1018
1019   length = 0;
1020   score = 0;
1021   last_stripq++; /* this is used to mark quadrangles in this strip */
1022
1023   /* go forwards till possible... */
1024   cstrip.q = pstrip->q;
1025   cstrip.iv2 = pstrip->iv2;
1026   cstrip.iv3 = pstrip->iv3;
1027   while ((cstrip.q != 0) &&                             /* on a boundary */
1028          (quadranglesptr[cstrip.q].state != 0) &&       /* deleted */
1029          (quadranglesptr[cstrip.q].state != last_stripq))/* on the same strip */
1030   {
1031     quadranglesptr[cstrip.q].state = last_stripq;
1032     /* increment the length */
1033     length++;
1034     /* compute the score */
1035     /* increment the score if the quadrangle has less than four neighbours */
1036     for (i=0; i<4; i++)
1037     {
1038       quadrang = quadranglesptr[cstrip.q].qn[i];
1039       if ((quadrang == 0) || (quadranglesptr[quadrang].state == 0))
1040       {
1041         score++;
1042         break;
1043       }
1044     }
1045     /* next in the strip */
1046     stripq_next(&cstrip);
1047   }
1048   /* turn back... */
1049   cstrip.q = pstrip->q;
1050   cstrip.iv2 = (pstrip->iv2+2)%4;
1051   cstrip.iv3 = (pstrip->iv3+2)%4;
1052   /* ... but save the position of the strip before moving */
1053   savstrip.q = cstrip.q;
1054   savstrip.iv2 = cstrip.iv2;
1055   savstrip.iv3 = cstrip.iv3;
1056   stripq_next(&cstrip);
1057   while ((cstrip.q != 0) && 
1058          (quadranglesptr[cstrip.q].state != 0) &&
1059          (quadranglesptr[cstrip.q].state != last_stripq))
1060   {
1061     quadranglesptr[cstrip.q].state = last_stripq;
1062     /* save the position of the strip each time before moving */
1063     savstrip.q = cstrip.q;
1064     savstrip.iv2 = cstrip.iv2;
1065     savstrip.iv3 = cstrip.iv3;
1066     /* increment the length */
1067     length++;
1068     /* compute the score */
1069     /* increment the score if the quadrangle has less than four neighbours */
1070     for (i=0; i<4; i++)
1071     {
1072       quadrang = quadranglesptr[cstrip.q].qn[i];
1073       if ((quadrang == 0) || (quadranglesptr[quadrang].state == 0))
1074       {
1075         score++;
1076         break;
1077       }
1078     }
1079     /* next in the strip */
1080     stripq_next(&cstrip);
1081   }
1082   /* ... back end reached. Now turn forward again at recent saved position. */
1083   pstrip->q = savstrip.q;
1084   pstrip->iv2 = (savstrip.iv2+2)%4;
1085   pstrip->iv3 = (savstrip.iv3+2)%4;
1086   *plength = length;
1087   return score;
1088 }
1089
1090 /********************************************************************/
1091 /*                                                                  */
1092 /* stripq_next : jump to next quadrangle in a strip                 */
1093 /*                                                                  */
1094 /********************************************************************/
1095
1096 void stripq_next(stripq *st) 
1097 { /* st points toward a quadrangle and a vertex ordering defining a unique. */
1098   /* Its content is changed for the next quadrangle in the strip. */
1099   /* The quadrangle may become 0 if there was no neighbour. */
1100   int   quadrang=st->q; /* current */
1101   int   i=st->iv2;
1102   int   nquadrang=quadranglesptr[quadrang].qn[i];       /* neighbour */
1103
1104   if (!quadrang || !nquadrang)
1105   { /* There is no neighbour on this edge. */
1106     st->q = 0;
1107     st->iv2 = 0;
1108     st->iv3 = 0;
1109   }
1110   else
1111   { /* Compute the new index for the last vertex. */
1112     st->q = nquadrang;
1113     st->iv2 = quadranglesptr[quadrang].ivn[i][0];
1114     st->iv3 = quadranglesptr[quadrang].ivn[i][1];
1115   }
1116 }
1117
1118 /*******************************************************************/
1119 /*********                                                **********/
1120 /*********              COMMENTS                          **********/
1121 /*********                                                **********/
1122 /*******************************************************************/
1123
1124
1125 /* 
1126   Architecture
1127   ************
1128
1129 The present C implementation was designed to be called from FORTRAN
1130 with the following syntaxes:
1131
1132 1. STRIPQ_INIT(int *NBNBVERTICES, int *NBQUADRANGLES, int *TABQUADRANGLES)
1133
1134    This is the initiating call where:
1135    NBNBVERTICES is the number of NBVERTICES
1136    NBQUADRNGLES is the number of quadrangles
1137    TABQUADRANGLES is the table describing quadrangles
1138
1139    This  function copies  its arguments to  nbNBVERTICES, nquadrangles, and
1140 build an inner table of quadrangles, where neighbours of a
1141 quadrangle can be found easily.
1142
1143
1144 IMPORTANT WARNINGS
1145   NBVERTICES and Quadrangles are in the range 1...NBNBVERTICES and 1...NBQUADRANGLES
1146   Double arrays are FORTRAN arrays so the three NBVERTICES of quadrangle I
1147   are found in TABQUADRANGLES[I+J-1] where J = 0, 1, 2, 3.
1148   The FORTRAN array is declared as TABQUADRANGLES(4, NBQUADRANGLES)
1149
1150
1151   2.  STRIPQ_GET_STRIP(int *NBQUAD, int *V1, int *V2)
1152       STRIPQ_GET_NEXT(int *VERTEX1, int *VERTEX2, int *QUADRANGLE)
1153 ???
1154       Both functions are used to get the strips. Each iteration of the
1155      GET_STRIP brings a new strip where NBQUAD is the number of quadrangles
1156      in the strip (NB: number of NBVERTICES would be NBQUADS*2+2) and  V1, V2
1157      are the  two first NBVERTICES.  NBQUADS becomes   zero when the last
1158      strip has been read.  STRIPQ_GET_NEXT is  used to get the successive
1159      NBVERTICES of a strip, each two vertice are associated with next quadrangle.
1160      This start with  the 3d and 4th vertice, the two first are given  by
1161      STRIPQ_GET_STRIP.
1162
1163      An example of correct call from C to read the strips is
1164
1165      while(1)
1166        {
1167          STRIPQ_GET_STRIP(&nbquad, &v[0], &v[1]);
1168        if (nbquad == 0)
1169          break;
1170        for (i=1; i <= nbquad; i++)
1171        {
1172          STRIPQ_GET_NEXT(&v[i*2], &v[i*2+1], &q[i]);
1173        }
1174      }
1175
1176        Unwise calls to those functions will generally return zero but
1177        are not recommanded.
1178
1179        The data are not checked for coherence, but a zero
1180        number of quadrangle will give a zero number of strips.
1181       
1182 */
1183
1184 /*
1185 OUTLINE OF THE ALGORITHM
1186 ************************
1187
1188 The algorithm is purely topological.  No coordinates  are given,  its
1189 input are the number  of  NBVERTICES, the number of quadangles, and for
1190 each quadrangle a quadruplet of NBVERTICES indices.
1191
1192 Let us consider a   quadrangle Q=(V1, V2, V3, V4), this quadrangle has
1193 neighbours in the quadrangle graphs, let us call them Q12, Q23, Q34 and Q41.
1194 Q12 is the quadrangle sharing the edge [V1, V2] with T, etc... Of course those
1195 fouree quadrangles may not exist in the graph, in this case T is "on the
1196 border" or even "in a corner".
1197
1198 The key remark is that at most two quadrangle-strips may cross Q, they
1199 are: Q12-Q-Q34, Q23-Q-Q41. Once four adjacent quadrangles
1200 are given the entire strip is uniquely defined. The orientation of a
1201 strip is meaningless as if you reverse it you would get the same strip.
1202 To describe a current position in a strip you need a quadrangle and the
1203 two last NBVERTICES of the quadrangle in the strip.
1204
1205 Our algorithm (more precisely heuristic) is the following:
1206 - Find a quadrangle with the lowest number of neighbours.
1207 - List the four (or less) strips crossing this quadrangle.
1208 - Chose the best among them and remove from the graph all the
1209 quadrangles in this strip.
1210   The best strip is rated with a "score", the score we used is the
1211   number of quadrangles in the strip which have less than four
1212   neighbours  (they are on the "border") in case of score equality the
1213   longest strip is selected.
1214 - Reiterate the process until there are no quadrangles left.
1215
1216 There are no demonstrations of  the optimality  of this algorithm, but
1217 it seems to give expected results on regular graphs which are the most
1218 commonly fed to it.
1219 */
1220
1221 /*
1222 Implementation
1223 **************
1224
1225 First the STRIPQ_INIT function stores the quadrangles in a data structure
1226 (the quadrangles array allocated on free  store),  containing for each
1227 quadrangle its four neighbours, then third and fourth vertice for each neighbour
1228 (a zero neighbour is inexistant), a quadrangle gets also a state 1.  To
1229 build this  array a temporary  array "edges" is build giving  for each
1230 edge (pair of NBVERTICES) the two neghbouring quadrangles.
1231
1232 Then at  each  call    of STRIPQ_GET_STRIP  a  quadrangle  with minimum
1233 neighbours is first chosen.  For   the four possible strips  crossing
1234 the quadrangle the strip_score function is called  which brings back the
1235 start of the strip, the length and the score.  The strip-next function
1236 is used to jump to  the next quadrangle  in a strip.   The best strip is
1237 chosen, stored in the current_strip  and returned. 
1238
1239
1240 Each   call  to   STRIPQ_GET_NEXT  increments the current-strip
1241 structure  to the next quadrangle in   the  strip,  using the strip_next
1242 function. The quadrangle is deleted from the graph and returned.
1243
1244 The last call  to STRIPQ_GET_STRIP frees  the  quadrangles array to the
1245 free store.
1246
1247
1248 */