src/gts/split.c

   1 /* GTS - Library for the manipulation of triangulated surfaces
   2  * Copyright (C) 1999 Stéphane Popinet
   3  *
   4  * This library is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Library General Public
   6  * License as published by the Free Software Foundation; either
   7  * version 2 of the License, or (at your option) any later version.
   8  *
   9  * This library is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Library General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Library General Public
  15  * License along with this library; if not, write to the
  16  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  17  * Boston, MA 02111-1307, USA.
  18  */
  19
  20 #include <stdlib.h>
  21 #include <string.h>
  22 #include "gts.h"
  23
  24 #define DYNAMIC_SPLIT
  25 #define NEW
  26
  27 /* #define DEBUG
  28    #define DEBUG_HEXPAND
  29    #define DEBUG_EXPAND */
  30
  31 struct _GtsSplitCFace {
  32   GtsFace * f;
  33   GtsTriangle ** a1, ** a2;
  34 };
  35
  36 typedef struct _CFace      CFace;
  37 typedef struct _CFaceClass CFaceClass;
  38
  39 struct _CFace {
  40   GtsObject object;
  41
  42   GtsSplit * parent_split;
  43   GtsTriangle * t;
  44   guint flags;
  45 };
  46 /* the size of the CFace structure must be smaller or equal to the size
  47    of the GtsFace structure as both structures use the same memory location */
  48
  49 struct _CFaceClass {
  50   GtsObjectClass parent_class;
  51 };
  52
  53 #define IS_CFACE(obj) (gts_object_is_from_class (obj, cface_class ()))
  54 #define CFACE(obj)    ((CFace *) obj)
  55 #define CFACE_ORIENTATION(cf) ((cf)->flags & 0x1)
  56 #define CFACE_ORIENTATION_DIRECT(cf) ((cf)->flags |= 0x1)
  57 #define CFACE_VVS(cf)                ((cf)->flags & 0x2)
  58 #define CFACE_VVS_DIRECT(cf)         ((cf)->flags |= 0x2)
  59 #define CFACE_E1                     0x4
  60 #define CFACE_E2                     0x8
  61 #define CFACE_KEEP_VVS               0x10
  62
  63 #define ROTATE_ORIENT(e, e1, e2, e3)  { if (e1 == e) { e1 = e2; e2 = e3; }\
  64                                  else if (e2 == e) { e2 = e1; e1 = e3; }\
  65                                  else g_assert (e3 == e); }
  66 #define SEGMENT_USE_VERTEX(s, v) ((s)->v1 == v || (s)->v2 == v)
  67 #define TRIANGLE_REPLACE_EDGE(t, e, with) { if ((t)->e1 == e)\
  68                                               (t)->e1 = with;\
  69                                             else if ((t)->e2 == e)\
  70                                               (t)->e2 = with;\
  71                                             else {\
  72                                               g_assert ((t)->e3 == e);\
  73                                               (t)->e3 = with;\
  74                                             }\
  75                                           }
  76
  77 #define HEAP_INSERT_OBJECT(h, e) (GTS_OBJECT (e)->reserved =\
  78                                   gts_eheap_insert (h, e))
  79 #define HEAP_REMOVE_OBJECT(h, e) (gts_eheap_remove (h, GTS_OBJECT (e)->reserved),\
  80                                    GTS_OBJECT (e)->reserved = NULL)
  81
  82 static GtsObjectClass * cface_class (void)
  83 {
  84   static GtsObjectClass * klass = NULL;
  85
  86   if (klass == NULL) {
  87     GtsObjectClassInfo cface_info = {
  88       "GtsCFace",
  89       sizeof (CFace),
  90       sizeof (CFaceClass),
  91       (GtsObjectClassInitFunc) NULL,
  92       (GtsObjectInitFunc) NULL,
  93       (GtsArgSetFunc) NULL,
  94       (GtsArgGetFunc) NULL
  95     };
  96     klass = gts_object_class_new (gts_object_class (), &cface_info);
  97     g_assert (sizeof (CFace) <= sizeof (GtsFace));
  98   }
  99
 100   return klass;
 101 }
 102
 103 /* Replace @e with @with for all the triangles using @e but @f.
 104    Destroys @e and removes it from @heap (if not %NULL).
 105    Returns a triangle using e different from f or %NULL. */
 106 static GtsTriangle * replace_edge_collapse (GtsEdge * e,
 107                                             GtsEdge * with,
 108                                             CFace * cf,
 109                                             GtsEHeap * heap
 110 #ifdef DYNAMIC_SPLIT
 111                                             , GtsTriangle *** a1
 112 #endif
 113 #ifdef NEW
 114                                             , guint edge_flag
 115 #endif
 116                                             )
 117 {
 118   GSList * i;
 119   GtsTriangle * rt = NULL;
 120 #ifdef DYNAMIC_SPLIT
 121   guint size;
 122   GtsTriangle ** a;
 123 #endif
 124
 125 #ifdef NEW
 126   i = e->triangles;
 127   e->triangles = NULL;
 128   size = g_slist_length (i)*sizeof (GtsTriangle *);
 129   *a1 = a = g_malloc (size > 0 ? size : sizeof (GtsTriangle *));
 130   while (i) {
 131     GtsTriangle * t = i->data;
 132     GSList * next = i->next;
 133     if (t != ((GtsTriangle *) cf)) {
 134       if (IS_CFACE (t)) {
 135         i->next = e->triangles;
 136         e->triangles = i;
 137         /* set the edge given by edge_flag (CFACE_E1 or CFACE_E2) */
 138         GTS_OBJECT (t)->reserved = GUINT_TO_POINTER (edge_flag);
 139         cf->flags |= CFACE_KEEP_VVS;
 140       }
 141       else {
 142         TRIANGLE_REPLACE_EDGE (t, e, with);
 143         i->next = with->triangles;
 144         with->triangles = i;
 145         rt = t;
 146         *(a++) = t;
 147       }
 148     }
 149     else
 150       g_slist_free_1 (i);
 151     i = next;
 152   }
 153   *a = NULL;
 154   if (!e->triangles) {
 155     if (heap)
 156       HEAP_REMOVE_OBJECT (heap, e);
 157     gts_object_destroy (GTS_OBJECT (e));
 158   }
 159 #else /* not NEW */
 160   i = e->triangles;
 161 #ifdef DYNAMIC_SPLIT
 162   size = g_slist_length (i)*sizeof (GtsTriangle *);
 163   *a1 = a = g_malloc (size > 0 ? size : sizeof (GtsTriangle *));
 164 #endif
 165   while (i) {
 166     GtsTriangle * t = i->data;
 167     GSList * next = i->next;
 168     if (t != ((GtsTriangle *) cf)) {
 169       TRIANGLE_REPLACE_EDGE (t, e, with);
 170       i->next = with->triangles;
 171       with->triangles = i;
 172       rt = t;
 173 #ifdef DYNAMIC_SPLIT
 174       *(a++) = t;
 175 #endif
 176     }
 177     else
 178       g_slist_free_1 (i);
 179     i = next;
 180   }
 181 #ifdef DYNAMIC_SPLIT
 182   *a = NULL;
 183 #endif
 184   if (heap)
 185     HEAP_REMOVE_OBJECT (heap, e);
 186   e->triangles = NULL;
 187   gts_object_destroy (GTS_OBJECT (e));
 188 #endif /* NEW */
 189
 190   return rt;
 191 }
 192
 193 static CFace * cface_new (GtsFace * f,
 194                           GtsEdge * e,
 195                           GtsVertex * v1,
 196                           GtsVertex * v2,
 197                           GtsSplit * vs,
 198                           GtsEHeap * heap,
 199                           GtsEdgeClass * klass
 200 #ifdef DYNAMIC_SPLIT
 201                           , GtsSplitCFace * scf
 202 #endif
 203                           )
 204 {
 205   CFace * cf;
 206   GtsVertex * v;
 207   GtsEdge * e1, * e2, * e3, * vvs;
 208   GSList * i;
 209   GtsTriangle * t, * t1 = NULL, * t2 = NULL;
 210   guint flags;
 211
 212   g_return_val_if_fail (f != NULL, NULL);
 213 #ifndef NEW
 214   g_return_val_if_fail (GTS_IS_FACE (f), NULL);
 215 #endif
 216   g_return_val_if_fail (e != NULL, NULL);
 217   g_return_val_if_fail (vs != NULL, NULL);
 218
 219   t = ((GtsTriangle *) f);
 220   if (heap)
 221     g_return_val_if_fail (!gts_triangle_is_duplicate (t), NULL);
 222
 223 #ifdef NEW
 224   /* get CFACE_E1 and CFACE_E2 info */
 225   flags = GPOINTER_TO_UINT (GTS_OBJECT (f)->reserved);
 226 #endif
 227   GTS_OBJECT_SET_FLAGS (f, GTS_DESTROYED);
 228
 229   i = f->surfaces;
 230   while (i) {
 231     GSList * next = i->next;
 232     gts_surface_remove_face (i->data, f);
 233     i = next;
 234   }
 235   g_slist_free (f->surfaces);
 236
 237   e1 = t->e1; e2 = t->e2; e3 = t->e3;
 238   ROTATE_ORIENT (e, e1, e2, e3);
 239
 240   cf = (CFace *) f;
 241 #ifndef NEW
 242   GTS_OBJECT (cf)->klass = cface_class ();
 243 #else
 244   cf->flags = flags;
 245 #endif
 246   gts_object_init (GTS_OBJECT (cf), cface_class ());
 247   cf->parent_split = vs;
 248
 249   if (SEGMENT_USE_VERTEX (GTS_SEGMENT (e1), v2)) {
 250     CFACE_ORIENTATION_DIRECT (cf); /* v1->v2->v */
 251     e3 = e1; e1 = e2; e2 = e3;
 252   }
 253   v = GTS_SEGMENT (e1)->v1 == v1 ?
 254     GTS_SEGMENT (e1)->v2 : GTS_SEGMENT (e1)->v1;
 255 #ifdef NEW
 256   if ((cf->flags & CFACE_E1) || (cf->flags & CFACE_E2))
 257     g_assert ((vvs = GTS_EDGE (gts_vertices_are_connected (vs->v, v))));
 258   else
 259 #endif
 260   vvs = gts_edge_new (klass, v, vs->v);
 261
 262   t1 = replace_edge_collapse (e1, vvs, cf, heap
 263 #ifdef DYNAMIC_SPLIT
 264                               , &scf->a1
 265 #endif
 266 #ifdef NEW
 267                               , CFACE_E1
 268 #endif
 269                               );
 270   t2 = replace_edge_collapse (e2, vvs, cf, heap
 271 #ifdef DYNAMIC_SPLIT
 272                               , &scf->a2
 273 #endif
 274 #ifdef NEW
 275                               , CFACE_E2
 276 #endif
 277                               );
 278   t = cf->t = t1 ? t1 : t2;
 279   g_assert (t);
 280
 281   /* set up flags necessary to find vvs */
 282   if (t->e1 == vvs) e2 = t->e2;
 283   else if (t->e2 == vvs) e2 = t->e3;
 284   else {
 285     g_assert (t->e3 == vvs);
 286     e2 = t->e1;
 287   }
 288   if (SEGMENT_USE_VERTEX (GTS_SEGMENT (e2), v))
 289     CFACE_VVS_DIRECT (cf);
 290
 291   return cf;
 292 }
 293
 294 static void find_vvs (GtsVertex * vs,
 295                       GtsTriangle * t,
 296                       GtsVertex ** v, GtsEdge ** vvs,
 297                       gboolean orientation)
 298 {
 299   GtsEdge * e1 = t->e1, * e2 = t->e2, * e3 = t->e3, * tmp;
 300
 301   if (SEGMENT_USE_VERTEX (GTS_SEGMENT (e2), vs)) {
 302     tmp = e1; e1 = e2; e2 = e3; e3 = tmp;
 303   }
 304   else if (SEGMENT_USE_VERTEX (GTS_SEGMENT (e3), vs)) {
 305     tmp = e1; e1 = e3; e3 = e2; e2 = tmp;
 306   }
 307   else
 308     g_assert (SEGMENT_USE_VERTEX (GTS_SEGMENT (e1), vs));
 309   if (SEGMENT_USE_VERTEX (GTS_SEGMENT (e2), vs) ||
 310       !gts_segments_touch (GTS_SEGMENT (e1), GTS_SEGMENT (e2))) {
 311     tmp = e1; e1 = e2; e2 = e3; e3 = tmp;
 312     g_assert (gts_segments_touch (GTS_SEGMENT (e1), GTS_SEGMENT (e2)));
 313   }
 314
 315   *vvs = orientation ? e1 : e3;
 316
 317   if (GTS_SEGMENT (*vvs)->v1 != vs) {
 318     g_assert (GTS_SEGMENT (*vvs)->v2 == vs);
 319     *v = GTS_SEGMENT (*vvs)->v1;
 320   }
 321   else
 322     *v = GTS_SEGMENT (*vvs)->v2;
 323 }
 324
 325 static void replace_edge_expand (GtsEdge * e,
 326                                  GtsEdge * with,
 327                                  GtsTriangle ** a,
 328                                  GtsVertex * v)
 329 {
 330   GtsTriangle ** i = a, * t;
 331
 332   while ((t = *(i++))) {
 333 #ifdef DEBUG_EXPAND
 334     g_assert (!IS_CFACE (t));
 335     fprintf (stderr, "replacing %p->%d: e: %p->%d with: %p->%d\n",
 336              t, id (t), e, id (e), with, id (with));
 337 #endif
 338     TRIANGLE_REPLACE_EDGE (t, e, with);
 339     with->triangles = g_slist_prepend (with->triangles, t);
 340     if (GTS_OBJECT (t)->reserved) {
 341       /* apart from the triangles having e as an edge, t is the only
 342          triangle using v */
 343       g_assert (GTS_OBJECT (t)->reserved == v);
 344       GTS_OBJECT (t)->reserved = NULL;
 345     }
 346     else
 347       GTS_OBJECT (t)->reserved = v;
 348   }
 349 }
 350
 351 static void cface_expand (CFace * cf,
 352                           GtsTriangle ** a1,
 353                           GtsTriangle ** a2,
 354                           GtsEdge * e,
 355                           GtsVertex * v1,
 356                           GtsVertex * v2,
 357                           GtsVertex * vs,
 358                           GtsEdgeClass * klass)
 359 {
 360   GtsVertex * v;
 361   GtsEdge * e1, * e2, * vvs;
 362   gboolean orientation;
 363   guint flags;
 364
 365   g_return_if_fail (cf != NULL);
 366   g_return_if_fail (IS_CFACE (cf));
 367   g_return_if_fail (e != NULL);
 368   g_return_if_fail (vs != NULL);
 369
 370   flags = cf->flags;
 371   orientation = CFACE_ORIENTATION (cf);
 372
 373   find_vvs (vs, cf->t, &v, &vvs, CFACE_VVS (cf));
 374
 375 #ifdef NEW
 376   if (flags & CFACE_E1)
 377     e1 = GTS_EDGE (gts_vertices_are_connected (v1, v));
 378   else
 379     e1 = gts_edge_new (klass, v, v1);
 380   if (flags & CFACE_E2)
 381     e2 = GTS_EDGE (gts_vertices_are_connected (v2, v));
 382   else
 383     e2 = gts_edge_new (klass, v, v2);
 384 #else
 385   e1 = gts_edge_new (v, v1);
 386   e2 = gts_edge_new (v, v2);
 387 #endif
 388
 389   replace_edge_expand (vvs, e1, a1, v1);
 390   replace_edge_expand (vvs, e2, a2, v2);
 391
 392 #ifdef NEW
 393   if (!(flags & CFACE_KEEP_VVS)) {
 394     g_slist_free (vvs->triangles);
 395     vvs->triangles = NULL;
 396     gts_object_destroy (GTS_OBJECT (vvs));
 397   }
 398 #else
 399   g_slist_free (vvs->triangles);
 400   vvs->triangles = NULL;
 401   gts_object_destroy (GTS_OBJECT (vvs));
 402 #endif
 403
 404   /* gts_face_new : because I am "creating" a face */
 405   GTS_OBJECT (cf)->klass = GTS_OBJECT_CLASS (gts_face_class ());
 406   gts_object_init (GTS_OBJECT (cf), GTS_OBJECT (cf)->klass);
 407
 408   if (orientation)
 409     gts_triangle_set (GTS_TRIANGLE (cf), e, e2, e1);
 410   else
 411     gts_triangle_set (GTS_TRIANGLE (cf), e, e1, e2);
 412 }
 413
 414 static void split_destroy (GtsObject * object)
 415 {
 416   GtsSplit * vs = GTS_SPLIT (object);
 417   guint i = vs->ncf;
 418   GtsSplitCFace * cf = vs->cfaces;
 419
 420   while (i--) {
 421     if (IS_CFACE (cf->f))
 422       gts_object_destroy (GTS_OBJECT (cf->f));
 423     g_free (cf->a1);
 424     g_free (cf->a2);
 425     cf++;
 426   }
 427   g_free (vs->cfaces);
 428
 429   if (!gts_allow_floating_vertices && vs->v && vs->v->segments == NULL)
 430     gts_object_destroy (GTS_OBJECT (vs->v));
 431
 432   (* GTS_OBJECT_CLASS (gts_split_class ())->parent_class->destroy) (object);
 433 }
 434
 435 static void split_class_init (GtsObjectClass * klass)
 436 {
 437   klass->destroy = split_destroy;
 438 }
 439
 440 static void split_init (GtsSplit * split)
 441 {
 442   split->v1 = split->v2 = NULL;
 443   split->v = NULL;
 444   split->cfaces = NULL;
 445   split->ncf = 0;
 446 }
 447
 448 /**
 449  * gts_split_class:
 450  *
 451  * Returns: the #GtsSplitClass.
 452  */
 453 GtsSplitClass * gts_split_class (void)
 454 {
 455   static GtsSplitClass * klass = NULL;
 456
 457   if (klass == NULL) {
 458     GtsObjectClassInfo split_info = {
 459       "GtsSplit",
 460       sizeof (GtsSplit),
 461       sizeof (GtsSplitClass),
 462       (GtsObjectClassInitFunc) split_class_init,
 463       (GtsObjectInitFunc) split_init,
 464       (GtsArgSetFunc) NULL,
 465       (GtsArgGetFunc) NULL
 466     };
 467     klass = gts_object_class_new (gts_object_class (),
 468                                   &split_info);
 469   }
 470
 471   return klass;
 472 }
 473
 474 #ifdef DEBUG
 475 static gboolean edge_collapse_is_valid (GtsEdge * e)
 476 {
 477   GSList * i;
 478
 479   g_return_val_if_fail (e != NULL, FALSE);
 480
 481   if (gts_segment_is_duplicate (GTS_SEGMENT (e))) {
 482     g_warning ("collapsing duplicate edge");
 483     return FALSE;
 484   }
 485
 486   i = GTS_SEGMENT (e)->v1->segments;
 487   while (i) {
 488     GtsEdge * e1 = i->data;
 489     if (e1 != e && GTS_IS_EDGE (e1)) {
 490       GtsEdge * e2 = NULL;
 491       GSList * j = GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e)->v1 ?
 492         GTS_SEGMENT (e1)->v2->segments : GTS_SEGMENT (e1)->v1->segments;
 493       while (j && !e2) {
 494         GtsEdge * e1 = j->data;
 495         if (GTS_IS_EDGE (e1) &&
 496             (GTS_SEGMENT (e1)->v1 == GTS_SEGMENT (e)->v2 ||
 497              GTS_SEGMENT (e1)->v2 == GTS_SEGMENT (e)->v2))
 498           e2 = e1;
 499         j = j->next;
 500       }
 501       if (e2 && !gts_triangle_use_edges (e, e1, e2)) {
 502         g_warning ("collapsing empty triangle");
 503         return FALSE;
 504       }
 505     }
 506     i = i->next;
 507   }
 508
 509   if (gts_edge_is_boundary (e, NULL)) {
 510     GtsTriangle * t = e->triangles->data;
 511     if (gts_edge_is_boundary (t->e1, NULL) &&
 512         gts_edge_is_boundary (t->e2, NULL) &&
 513         gts_edge_is_boundary (t->e3, NULL)) {
 514       g_warning ("collapsing single triangle");
 515       return FALSE;
 516     }
 517   }
 518   else {
 519     if (gts_vertex_is_boundary (GTS_SEGMENT (e)->v1, NULL) &&
 520         gts_vertex_is_boundary (GTS_SEGMENT (e)->v2, NULL)) {
 521       g_warning ("collapsing two sides of a strip");
 522       return FALSE;
 523     }
 524     if (gts_edge_belongs_to_tetrahedron (e)) {
 525       g_warning ("collapsing tetrahedron");
 526       return FALSE;
 527     }
 528   }
 529
 530   return TRUE;
 531 }
 532 #endif /* DEBUG */
 533
 534 /* Not currently used.  May be useful for some debug code */
 535 #ifdef DEBUG
 536 static void print_split (GtsSplit * vs, FILE * fptr)
 537 {
 538   guint j;
 539   GtsSplitCFace * cf;
 540
 541   g_return_if_fail (vs != NULL);
 542   g_return_if_fail (fptr != NULL);
 543
 544   fprintf (fptr, "%p: v: %p v1: %p v2: %p ncf: %u cfaces: %p\n",
 545            vs, vs->v, vs->v1, vs->v2, vs->ncf, vs->cfaces);
 546   cf = vs->cfaces;
 547   j = vs->ncf;
 548   while (j--) {
 549     fprintf (stderr, "  f: %p a1: %p a2: %p\n",
 550              cf->f, cf->a1, cf->a2);
 551     cf++;
 552   }
 553 }
 554 #endif
 555
 556 /**
 557  * gts_split_collapse:
 558  * @vs: a #GtsSplit.
 559  * @klass: a #GtsEdgeClass.
 560  * @heap: a #GtsEHeap or %NULL.
 561  *
 562  * Collapses the vertex split @vs. Any new edge created during the process will
 563  * be of class @klass. If heap is not %NULL, the new edges will be inserted
 564  * into it and the destroyed edges will be removed from it.
 565  */
 566 void gts_split_collapse (GtsSplit * vs,
 567                          GtsEdgeClass * klass,
 568                          GtsEHeap * heap)
 569 {
 570   GtsEdge * e;
 571   GtsVertex * v, * v1, * v2;
 572   GSList * i, * end;
 573 #ifdef DYNAMIC_SPLIT
 574   GtsSplitCFace * cf;
 575   guint j;
 576 #endif
 577 #ifdef DEBUG
 578   gboolean invalid = FALSE;
 579   static guint ninvalid = 0;
 580 #endif
 581
 582   g_return_if_fail (vs != NULL);
 583   g_return_if_fail (klass != NULL);
 584
 585   v = vs->v;
 586
 587   g_return_if_fail (v->segments == NULL);
 588
 589   /* we don't want to destroy vertices */
 590   gts_allow_floating_vertices = TRUE;
 591
 592   v1 = GTS_SPLIT_V1 (vs);
 593   v2 = GTS_SPLIT_V2 (vs);
 594   g_assert ((e = GTS_EDGE (gts_vertices_are_connected (v1, v2))));
 595
 596 #ifdef DEBUG
 597   fprintf (stderr, "collapsing %p: v1: %p v2: %p v: %p\n", vs, v1, v2, v);
 598   if (!edge_collapse_is_valid (e)) {
 599     char fname[80];
 600     FILE * fptr;
 601     GSList * triangles, * i;
 602
 603     g_warning ("invalid edge collapse");
 604     invalid = TRUE;
 605     sprintf (fname, "invalid.%d", ninvalid);
 606     fptr = fopen (fname, "wt");
 607     gts_write_segment (GTS_SEGMENT (e), GTS_POINT (v), fptr);
 608     triangles = gts_vertex_triangles (v1, NULL);
 609     triangles = gts_vertex_triangles (v2, triangles);
 610     i = triangles;
 611     while (i) {
 612       gts_write_triangle (i->data, GTS_POINT (v), fptr);
 613       i = i->next;
 614     }
 615     g_slist_free (triangles);
 616     fclose (fptr);
 617   }
 618 #endif
 619
 620   i = e->triangles;
 621 #ifdef DYNAMIC_SPLIT
 622   cf = vs->cfaces;
 623   j = vs->ncf;
 624   while (j--) {
 625     g_free (cf->a1);
 626     g_free (cf->a2);
 627     cf++;
 628   }
 629   g_free (vs->cfaces);
 630
 631   vs->ncf = g_slist_length (i);
 632   g_assert (vs->ncf > 0);
 633   cf = vs->cfaces = g_malloc (vs->ncf*sizeof (GtsSplitCFace));
 634 #endif /* DYNAMIC_SPLIT */
 635 #ifdef NEW
 636   while (i) {
 637     cf->f = i->data;
 638     g_assert (GTS_IS_FACE (cf->f));
 639     GTS_OBJECT (cf->f)->klass = GTS_OBJECT_CLASS (cface_class ());
 640     cf++;
 641     i = i->next;
 642   }
 643   i = e->triangles;
 644   cf = vs->cfaces;
 645   while (i) {
 646     cface_new (i->data, e, v1, v2, vs, heap, klass, cf);
 647 #ifdef DEBUG
 648     fprintf (stderr, "cface: %p->%d t: %p->%d a1: ",
 649              cf->f, id (cf->f), CFACE (cf->f)->t, id (CFACE (cf->f)->t));
 650     {
 651       GtsTriangle * t, ** a;
 652       a = cf->a1;
 653       while ((t = *(a++)))
 654         fprintf (stderr, "%p->%d ", t, id (t));
 655       fprintf (stderr, "a2: ");
 656       a = cf->a2;
 657       while ((t = *(a++)))
 658         fprintf (stderr, "%p->%d ", t, id (t));
 659       fprintf (stderr, "\n");
 660     }
 661 #endif
 662     cf++;
 663     i = i->next;
 664   }
 665 #else /* not NEW */
 666   while (i) {
 667     cface_new (i->data, e, v1, v2, vs, heap
 668 #ifdef DYNAMIC_SPLIT
 669                , cf
 670 #endif /* DYNAMIC_SPLIT */
 671                );
 672 #ifdef DYNAMIC_SPLIT
 673     cf->f = i->data;
 674     cf++;
 675 #endif /* DYNAMIC_SPLIT */
 676     i = i->next;
 677   }
 678 #endif /* NEW */
 679   g_slist_free (e->triangles);
 680   e->triangles = NULL;
 681   gts_object_destroy (GTS_OBJECT (e));
 682
 683   gts_allow_floating_vertices = FALSE;
 684
 685   end = NULL;
 686   i = v1->segments;
 687   while (i) {
 688     GtsSegment * s = i->data;
 689     if (s->v1 == v1)
 690       s->v1 = v;
 691     else
 692       s->v2 = v;
 693     end = i;
 694     i = i->next;
 695   }
 696   if (end) {
 697     end->next = v->segments;
 698     v->segments = v1->segments;
 699     v1->segments = NULL;
 700   }
 701
 702   end = NULL;
 703   i = v2->segments;
 704   while (i) {
 705     GtsSegment * s = i->data;
 706     if (s->v1 == v2)
 707       s->v1 = v;
 708     else
 709       s->v2 = v;
 710     end = i;
 711     i = i->next;
 712   }
 713   if (end) {
 714     end->next = v->segments;
 715     v->segments = v2->segments;
 716     v2->segments = NULL;
 717   }
 718
 719 #ifdef DEBUG
 720   if (invalid) {
 721     char fname[80];
 722     FILE * fptr;
 723     GSList * triangles, * i;
 724     GtsSurface * surface = NULL;
 725
 726     sprintf (fname, "invalid_after.%d", ninvalid);
 727     fptr = fopen (fname, "wt");
 728     triangles = gts_vertex_triangles (v, NULL);
 729     i = triangles;
 730     while (i) {
 731       GtsTriangle * t = i->data;
 732       fprintf (stderr, "checking %p->%d\n", t, id (t));
 733       g_assert (GTS_IS_FACE (t));
 734       gts_write_triangle (t, GTS_POINT (v), fptr);
 735       surface = GTS_FACE (t)->surfaces->data;
 736       if (gts_triangle_is_duplicate (t))
 737         fprintf (stderr, "%p->%d is duplicate\n", t, id (t));
 738       if (gts_segment_is_duplicate (GTS_SEGMENT (t->e1)))
 739         fprintf (stderr, "e1 of %p->%d is duplicate\n", t, id (t));
 740       if (gts_segment_is_duplicate (GTS_SEGMENT (t->e2)))
 741         fprintf (stderr, "e2 of %p->%d is duplicate\n", t, id (t));
 742       if (gts_segment_is_duplicate (GTS_SEGMENT (t->e3)))
 743         fprintf (stderr, "e3 of %p->%d is duplicate\n", t, id (t));
 744       i = i->next;
 745     }
 746     fclose (fptr);
 747     g_slist_free (triangles);
 748 #if 0
 749     gts_split_expand (vs, surface);
 750
 751     sprintf (fname, "invalid_after_after.%d", ninvalid);
 752     fptr = fopen (fname, "wt");
 753     triangles = gts_vertex_triangles (v1, NULL);
 754     triangles = gts_vertex_triangles (v2, triangles);
 755     i = triangles;
 756     while (i) {
 757       GtsTriangle * t = i->data;
 758       gts_write_triangle (t, GTS_POINT (v), fptr);
 759       surface = GTS_FACE (t)->surfaces->data;
 760       if (gts_triangle_is_duplicate (t))
 761         fprintf (stderr, "%p->%d is duplicate\n", t, id (t));
 762       if (gts_segment_is_duplicate (GTS_SEGMENT (t->e1)))
 763         fprintf (stderr, "e1 of %p->%d is duplicate\n", t, id (t));
 764       if (gts_segment_is_duplicate (GTS_SEGMENT (t->e2)))
 765         fprintf (stderr, "e2 of %p->%d is duplicate\n", t, id (t));
 766       if (gts_segment_is_duplicate (GTS_SEGMENT (t->e3)))
 767         fprintf (stderr, "e3 of %p->%d is duplicate\n", t, id (t));
 768       i = i->next;
 769     }
 770     fclose (fptr);
 771     g_slist_free (triangles);
 772
 773     exit (1);
 774 #endif
 775     ninvalid++;
 776   }
 777 #endif
 778 }
 779
 780 /**
 781  * gts_split_expand:
 782  * @vs: a #GtsSplit.
 783  * @s: a #GtsSurface.
 784  * @klass: a #GtsEdgeClass.
 785  *
 786  * Expands the vertex split @vs adding the newly created faces to @s. Any
 787  * new edge will be of class @klass.
 788  */
 789 void gts_split_expand (GtsSplit * vs,
 790                        GtsSurface * s,
 791                        GtsEdgeClass * klass)
 792 {
 793   GSList * i;
 794   GtsEdge * e;
 795   GtsVertex * v, * v1, * v2;
 796   gboolean changed = FALSE;
 797   GtsSplitCFace * cf;
 798   guint j;
 799
 800   g_return_if_fail (vs != NULL);
 801   g_return_if_fail (s != NULL);
 802   g_return_if_fail (klass != NULL);
 803
 804   /* we don't want to destroy vertices */
 805   gts_allow_floating_vertices = TRUE;
 806
 807   v1 = GTS_SPLIT_V1 (vs);
 808   v2 = GTS_SPLIT_V2 (vs);
 809   v = vs->v;
 810 #ifdef DEBUG_EXPAND
 811   fprintf (stderr, "expanding %p->%d: v1: %p->%d v2: %p->%d v: %p->%d\n",
 812            vs, id (vs), v1, id (v1), v2, id (v2), v, id (v));
 813 #endif
 814   e = gts_edge_new (klass, v1, v2);
 815   cf = vs->cfaces;
 816   j = vs->ncf;
 817   while (j--) {
 818     cface_expand (CFACE (cf->f), cf->a1, cf->a2, e, v1, v2, v, klass);
 819     gts_surface_add_face (s, cf->f);
 820     cf++;
 821   }
 822
 823   gts_allow_floating_vertices = FALSE;
 824
 825   /* this part is described by figure "expand.fig" */
 826   i = v->segments;
 827   while (i) {
 828     GtsEdge * e1 = i->data;
 829     GtsVertex * with = NULL;
 830     GSList * j = e1->triangles, * next = i->next;
 831     // fprintf (stderr, "e1: %p->%d\n", e1, id (e1));
 832     while (j && !with) {
 833       with = GTS_OBJECT (j->data)->reserved;
 834       j = j->next;
 835     }
 836     if (with) {
 837       j = e1->triangles;
 838       while (j) {
 839         GtsTriangle * t = j->data;
 840         if (GTS_OBJECT (t)->reserved) {
 841           g_assert (GTS_OBJECT (t)->reserved == with);
 842           GTS_OBJECT (t)->reserved = NULL;
 843         }
 844         else
 845           GTS_OBJECT (t)->reserved = with;
 846         j = j->next;
 847       }
 848       if (GTS_SEGMENT (e1)->v1 == v)
 849         GTS_SEGMENT (e1)->v1 = with;
 850       else
 851         GTS_SEGMENT (e1)->v2 = with;
 852
 853       v->segments = g_slist_remove_link (v->segments, i);
 854       i->next = with->segments;
 855       with->segments = i;
 856       changed = TRUE;
 857     }
 858     if (next)
 859       i = next;
 860     else {
 861       /* check for infinite loop (the crossed out case in
 862          figure "expand.fig") */
 863       g_assert (changed);
 864       changed = FALSE;
 865       i = v->segments;
 866     }
 867   }
 868 }
 869
 870 #ifndef DYNAMIC_SPLIT
 871 static void cface_neighbors (GtsSplitCFace * cf,
 872                              GtsEdge * e,
 873                              GtsVertex * v1,
 874                              GtsVertex * v2)
 875 {
 876   GtsTriangle * t = GTS_TRIANGLE (cf->f), ** a;
 877   GtsEdge * e1 = t->e1, * e2 = t->e2, * e3 = t->e3;
 878   GSList * i;
 879   guint size;
 880
 881   ROTATE_ORIENT (e, e1, e2, e3);
 882   if (SEGMENT_USE_VERTEX (GTS_SEGMENT (e1), v2)) {
 883     e3 = e1; e1 = e2; e2 = e3;
 884   }
 885
 886   i = e1->triangles;
 887   size = g_slist_length (i)*sizeof (GtsTriangle *);
 888   a = cf->a1 = g_malloc (size > 0 ? size : sizeof (GtsTriangle *));
 889   while (i) {
 890     if (i->data != t)
 891       *(a++) = i->data;
 892     i = i->next;
 893   }
 894   *a = NULL;
 895
 896   i = e2->triangles;
 897   size = g_slist_length (i)*sizeof (GtsTriangle *);
 898   a = cf->a2 = g_malloc (size > 0 ? size : sizeof (GtsTriangle *));
 899   while (i) {
 900     if (i->data != t)
 901       *(a++) = i->data;
 902     i = i->next;
 903   }
 904   *a = NULL;
 905 }
 906 #endif /*ifndef DYNAMIC_SPLIT */
 907
 908 /**
 909  * gts_split_new:
 910  * @klass: a #GtsSplitClass.
 911  * @v: a #GtsVertex.
 912  * @o1: either a #GtsVertex or a #GtsSplit.
 913  * @o2: either a #GtsVertex or a #GtsSplit.
 914  *
 915  * Creates a new #GtsSplit which would collapse @o1 and @o2 into @v. The
 916  * collapse itself is not performed.
 917  *
 918  * Returns: the new #GtsSplit.
 919  */
 920 GtsSplit * gts_split_new (GtsSplitClass * klass,
 921                           GtsVertex * v,
 922                           GtsObject * o1,
 923                           GtsObject * o2)
 924 {
 925   GtsSplit * vs;
 926   GtsVertex * v1, * v2;
 927 #ifndef DYNAMIC_SPLIT
 928   GtsEdge * e;
 929   GSList * i;
 930   GtsSplitCFace * cf;
 931 #endif
 932
 933   g_return_val_if_fail (klass != NULL, NULL);
 934   g_return_val_if_fail (v != NULL, NULL);
 935   g_return_val_if_fail (GTS_IS_SPLIT (o1) || GTS_IS_VERTEX (o1), NULL);
 936   g_return_val_if_fail (GTS_IS_SPLIT (o2) || GTS_IS_VERTEX (o2), NULL);
 937
 938   vs = GTS_SPLIT (gts_object_new (GTS_OBJECT_CLASS (klass)));
 939   vs->v = v;
 940   vs->v1 = o1;
 941   vs->v2 = o2;
 942   v1 = GTS_SPLIT_V1 (vs);
 943   v2 = GTS_SPLIT_V2 (vs);
 944 #ifdef DYNAMIC_SPLIT
 945   vs->ncf = 0;
 946   vs->cfaces = NULL;
 947 #else
 948   g_assert ((e = GTS_EDGE (gts_vertices_are_connected (v1, v2))));
 949   i = e->triangles;
 950   vs->ncf = g_slist_length (i);
 951   g_assert (vs->ncf > 0);
 952   cf = vs->cfaces = g_malloc (vs->ncf*sizeof (GtsSplitCFace));
 953   while (i) {
 954     cf->f = i->data;
 955     cface_neighbors (cf, e, v1, v2);
 956     i = i->next;
 957     cf++;
 958   }
 959 #endif
 960
 961   return vs;
 962 }
 963
 964 static gboolean
 965 split_traverse_pre_order (GtsSplit *           vs,
 966                           GtsSplitTraverseFunc func,
 967                           gpointer             data)
 968 {
 969   if (func (vs, data))
 970     return TRUE;
 971   if (GTS_IS_SPLIT (vs->v1) &&
 972       split_traverse_pre_order (GTS_SPLIT (vs->v1), func, data))
 973     return TRUE;
 974   if (GTS_IS_SPLIT (vs->v2) &&
 975       split_traverse_pre_order (GTS_SPLIT (vs->v2), func, data))
 976     return TRUE;
 977   return FALSE;
 978 }
 979
 980 static gboolean
 981 split_depth_traverse_pre_order (GtsSplit *             vs,
 982                                 guint                  depth,
 983                                 GtsSplitTraverseFunc   func,
 984                                 gpointer               data)
 985 {
 986   if (func (vs, data))
 987       return TRUE;
 988
 989   depth--;
 990   if (!depth)
 991     return FALSE;
 992
 993   if (GTS_IS_SPLIT (vs->v1) &&
 994       split_depth_traverse_pre_order (GTS_SPLIT (vs->v1), depth, func, data))
 995     return TRUE;
 996   if (GTS_IS_SPLIT (vs->v2) &&
 997       split_depth_traverse_pre_order (GTS_SPLIT (vs->v2), depth, func, data))
 998     return TRUE;
 999   return FALSE;
1000 }
1001
1002 static gboolean
1003 split_traverse_post_order (GtsSplit *           vs,
1004                            GtsSplitTraverseFunc func,
1005                            gpointer             data)
1006 {
1007   if (GTS_IS_SPLIT (vs->v1) &&
1008       split_traverse_post_order (GTS_SPLIT (vs->v1), func, data))
1009     return TRUE;
1010   if (GTS_IS_SPLIT (vs->v2) &&
1011       split_traverse_post_order (GTS_SPLIT (vs->v2), func, data))
1012     return TRUE;
1013   if (func (vs, data))
1014     return TRUE;
1015   return FALSE;
1016 }
1017
1018 static gboolean
1019 split_depth_traverse_post_order (GtsSplit *           vs,
1020                                  guint                depth,
1021                                  GtsSplitTraverseFunc func,
1022                                  gpointer             data)
1023 {
1024   depth--;
1025   if (depth) {
1026     if (GTS_IS_SPLIT (vs->v1) &&
1027         split_depth_traverse_post_order (GTS_SPLIT (vs->v1),
1028                                          depth, func, data))
1029       return TRUE;
1030     if (GTS_IS_SPLIT (vs->v2) &&
1031         split_depth_traverse_post_order (GTS_SPLIT (vs->v2),
1032                                          depth, func, data))
1033       return TRUE;
1034   }
1035   if (func (vs, data))
1036     return TRUE;
1037   return FALSE;
1038 }
1039
1040 /**
1041  * gts_split_traverse:
1042  * @root: the #GtsSplit to start the traversal from.
1043  * @order: the order in which nodes are visited - G_PRE_ORDER or G_POST_ORDER.
1044  * @depth: the maximum depth of the traversal. Nodes below this depth
1045  * will not be visited. If depth is -1 all nodes in the tree are
1046  * visited. If depth is 1, only the root is visited. If depth is 2,
1047  * the root and its children are visited. And so on.
1048  * @func: the function to call for each visited #GtsHSplit.
1049  * @data: user data to pass to the function.
1050  *
1051  * Traverses the #GtsSplit tree having @root as root. Calls @func for each
1052  * #GtsSplit of the tree in the order specified by @order. If order is set
1053  * to G_PRE_ORDER @func is called for the #GtsSplit then its children, if order
1054  * is set to G_POST_ORDER @func is called for the children and then for the
1055  * #GtsSplit.
1056  */
1057 void gts_split_traverse (GtsSplit *           root,
1058                          GTraverseType        order,
1059                          gint                 depth,
1060                          GtsSplitTraverseFunc func,
1061                          gpointer             data)
1062 {
1063   g_return_if_fail (root != NULL);
1064   g_return_if_fail (func != NULL);
1065   g_return_if_fail (order < G_LEVEL_ORDER);
1066   g_return_if_fail (depth == -1 || depth > 0);
1067
1068   switch (order) {
1069   case G_PRE_ORDER:
1070     if (depth < 0)
1071       split_traverse_pre_order (root, func, data);
1072     else
1073       split_depth_traverse_pre_order (root, depth, func, data);
1074     break;
1075   case G_POST_ORDER:
1076     if (depth < 0)
1077       split_traverse_post_order (root, func, data);
1078     else
1079       split_depth_traverse_post_order (root, depth, func, data);
1080     break;
1081   default:
1082     g_assert_not_reached ();
1083   }
1084 }
1085
1086 /**
1087  * gts_split_height:
1088  * @root: a #GtsSplit.
1089  *
1090  * Returns: the maximum height of the vertex split tree having @root as root.
1091  */
1092 guint gts_split_height (GtsSplit * root)
1093 {
1094   guint height = 0, tmp_height;
1095
1096   g_return_val_if_fail (root != NULL, 0);
1097
1098   if (GTS_IS_SPLIT (root->v1)) {
1099     tmp_height = gts_split_height (GTS_SPLIT (root->v1));
1100     if (tmp_height > height)
1101       height = tmp_height;
1102   }
1103   if (GTS_IS_SPLIT (root->v2)) {
1104     tmp_height = gts_split_height (GTS_SPLIT (root->v2));
1105     if (tmp_height > height)
1106       height = tmp_height;
1107   }
1108
1109   return height + 1;
1110 }
1111
1112 #ifndef DYNAMIC_SPLIT
1113 static gboolean list_array_are_identical (GSList * list,
1114                                           gpointer * array,
1115                                           gpointer excluded)
1116 {
1117   while (list) {
1118     gpointer data = list->data;
1119     if (data != excluded) {
1120       gboolean found = FALSE;
1121       gpointer * a = array;
1122
1123       while (!found && *a)
1124         if (*(a++) == data)
1125           found = TRUE;
1126       if (!found)
1127         return FALSE;
1128     }
1129     list = list->next;
1130   }
1131   return TRUE;
1132 }
1133 #endif /* ifndef DYNAMIC_SPLIT */
1134
1135 #ifndef NEW
1136 gboolean gts_split_is_collapsable (GtsSplit * vs)
1137 {
1138   guint i;
1139   GtsSplitCFace * cf;
1140   GtsVertex * v1, * v2;
1141   GtsEdge * e;
1142
1143   g_return_val_if_fail (vs != NULL, FALSE);
1144
1145   v1 = GTS_SPLIT_V1 (vs);
1146   v2 = GTS_SPLIT_V2 (vs);
1147   g_return_val_if_fail ((e = GTS_EDGE (gts_vertices_are_connected (v1, v2))),
1148                         FALSE);
1149
1150 #ifdef DYNAMIC_SPLIT
1151   if (!gts_edge_collapse_is_valid (e))
1152     return FALSE;
1153 #else
1154   i = vs->ncf;
1155   cf = vs->cfaces;
1156   while (i--) {
1157     GtsTriangle * t = GTS_TRIANGLE (cf->f);
1158     GtsEdge * e1 = t->e1, * e2 = t->e2, * e3 = t->e3;
1159
1160     ROTATE_ORIENT (e, e1, e2, e3);
1161     if (SEGMENT_USE_VERTEX (GTS_SEGMENT (e1), v2)) {
1162       e3 = e1; e1 = e2; e2 = e3;
1163     }
1164
1165     if (!list_array_are_identical (e1->triangles, (gpointer *) cf->a1, t))
1166       return FALSE;
1167     if (!list_array_are_identical (e2->triangles, (gpointer *) cf->a2, t))
1168       return FALSE;
1169
1170     cf++;
1171   }
1172 #endif
1173   return TRUE;
1174 }
1175 #endif /* not NEW */
1176
1177 #ifdef DEBUG_HEXPAND
1178 static guint expand_level = 0;
1179
1180 static void expand_indent (FILE * fptr)
1181 {
1182   guint i = expand_level;
1183   while (i--)
1184     fputc (' ', fptr);
1185 }
1186 #endif
1187
1188 /**
1189  * gts_hsplit_force_expand:
1190  * @hs: a #GtsHSplit.
1191  * @hsurface: a #GtsHSurface.
1192  *
1193  * Forces the expansion of @hs by first expanding all its dependencies not
1194  * already expanded.
1195  */
1196 void gts_hsplit_force_expand (GtsHSplit * hs,
1197                               GtsHSurface * hsurface)
1198 {
1199   guint i;
1200   GtsSplitCFace * cf;
1201
1202   g_return_if_fail (hs != NULL);
1203   g_return_if_fail (hsurface != NULL);
1204   g_return_if_fail (hs->nchild == 0);
1205
1206 #ifdef DEBUG_HEXPAND
1207   expand_level += 2;
1208 #endif
1209
1210   if (hs->parent && hs->parent->nchild == 0) {
1211 #ifdef DEBUG_HEXPAND
1212     expand_indent (stderr);
1213     fprintf (stderr, "expand parent %p\n", hs->parent);
1214 #endif
1215     gts_hsplit_force_expand (hs->parent, hsurface);
1216   }
1217
1218   i = GTS_SPLIT (hs)->ncf;
1219   cf = GTS_SPLIT (hs)->cfaces;
1220   while (i--) {
1221     GtsTriangle ** j, * t;
1222
1223     j = cf->a1;
1224     while ((t = *(j++)))
1225       if (IS_CFACE (t)) {
1226 #ifdef DEBUG_HEXPAND
1227         expand_indent (stderr);
1228         fprintf (stderr, "expand a1: cf->f: %p t: %p parent_split: %p\n",
1229                  cf->f,
1230                  t,
1231                  GTS_HSPLIT (CFACE (t)->parent_split));
1232 #endif
1233         gts_hsplit_force_expand (GTS_HSPLIT (CFACE (t)->parent_split),
1234                                  hsurface);
1235 #ifdef DEBUG_HEXPAND
1236         g_assert (!IS_CFACE (t));
1237 #endif
1238       }
1239     j = cf->a2;
1240     while ((t = *(j++)))
1241       if (IS_CFACE (t)) {
1242 #ifdef DEBUG_HEXPAND
1243         expand_indent (stderr);
1244         fprintf (stderr, "expand a2: cf->f: %p t: %p parent_split: %p\n",
1245                  cf->f,
1246                  t,
1247                  GTS_HSPLIT (CFACE (t)->parent_split));
1248 #endif
1249         gts_hsplit_force_expand (GTS_HSPLIT (CFACE (t)->parent_split),
1250                                  hsurface);
1251       }
1252     cf++;
1253   }
1254
1255   gts_hsplit_expand (hs, hsurface);
1256
1257 #ifdef DEBUG_HEXPAND
1258   expand_level -= 2;
1259   expand_indent (stderr);
1260   fprintf (stderr, "%p expanded\n", hs);
1261 #endif
1262 }
1263
1264 static void index_object (GtsObject * o, guint * n)
1265 {
1266   o->reserved = GUINT_TO_POINTER ((*n)++);
1267 }
1268
1269 static void index_face (GtsFace * f, gpointer * data)
1270 {
1271   guint * nf = data[1];
1272
1273   g_hash_table_insert (data[0], f, GUINT_TO_POINTER ((*nf)++));
1274 }
1275
1276 /**
1277  * gts_psurface_write:
1278  * @ps: a #GtsPSurface.
1279  * @fptr: a file pointer.
1280  *
1281  * Writes to @fptr a GTS progressive surface description.
1282  */
1283 void gts_psurface_write (GtsPSurface * ps, FILE * fptr)
1284 {
1285   guint nv = 1;
1286   guint nf = 1;
1287   GHashTable * hash;
1288   gpointer data[2];
1289
1290   g_return_if_fail (ps != NULL);
1291   g_return_if_fail (fptr != NULL);
1292   g_return_if_fail (GTS_PSURFACE_IS_CLOSED (ps));
1293
1294   while (gts_psurface_remove_vertex (ps))
1295     ;
1296
1297   GTS_POINT_CLASS (ps->s->vertex_class)->binary = FALSE;
1298   gts_surface_write (ps->s, fptr);
1299
1300   gts_surface_foreach_vertex (ps->s, (GtsFunc) index_object, &nv);
1301   hash = g_hash_table_new (NULL, NULL);
1302   data[0] = hash;
1303   data[1] = &nf;
1304   gts_surface_foreach_face (ps->s, (GtsFunc) index_face, data);
1305
1306   fprintf (fptr, "%u\n", ps->split->len);
1307   while (ps->pos) {
1308     GtsSplit * vs = g_ptr_array_index (ps->split, --ps->pos);
1309     GtsSplitCFace * scf = vs->cfaces;
1310     GtsVertex * v1, * v2;
1311     guint i = vs->ncf;
1312
1313     fprintf (fptr, "%u %u",
1314              GPOINTER_TO_UINT (GTS_OBJECT (vs->v)->reserved),
1315              vs->ncf);
1316     if (GTS_OBJECT (vs)->klass->write)
1317       (*GTS_OBJECT (vs)->klass->write) (GTS_OBJECT (vs), fptr);
1318     fputc ('\n', fptr);
1319
1320     v1 = GTS_IS_SPLIT (vs->v1) ? GTS_SPLIT (vs->v1)->v : GTS_VERTEX (vs->v1);
1321     GTS_OBJECT (v1)->reserved = GUINT_TO_POINTER (nv++);
1322     v2 = GTS_IS_SPLIT (vs->v2) ? GTS_SPLIT (vs->v2)->v : GTS_VERTEX (vs->v2);
1323     GTS_OBJECT (v2)->reserved = GUINT_TO_POINTER (nv++);
1324
1325     (*GTS_OBJECT (v1)->klass->write) (GTS_OBJECT (v1), fptr);
1326     fputc ('\n', fptr);
1327
1328     (*GTS_OBJECT (v2)->klass->write) (GTS_OBJECT (v2), fptr);
1329     fputc ('\n', fptr);
1330
1331     while (i--) {
1332       CFace * cf = CFACE (scf->f);
1333       GtsTriangle ** a, * t;
1334
1335       fprintf (fptr, "%u %u",
1336                GPOINTER_TO_UINT (g_hash_table_lookup (hash, cf->t)),
1337                cf->flags);
1338       if (GTS_OBJECT_CLASS (ps->s->face_class)->write)
1339         (*GTS_OBJECT_CLASS (ps->s->face_class)->write) (GTS_OBJECT (cf), fptr);
1340       fputc ('\n', fptr);
1341
1342       a = scf->a1;
1343       while ((t = *(a++)))
1344         fprintf (fptr, "%u ",
1345                  GPOINTER_TO_UINT (g_hash_table_lookup (hash, t)));
1346       fprintf (fptr, "\n");
1347
1348       a = scf->a2;
1349       while ((t = *(a++)))
1350         fprintf (fptr, "%u ",
1351                  GPOINTER_TO_UINT (g_hash_table_lookup (hash, t)));
1352       fprintf (fptr, "\n");
1353
1354       g_hash_table_insert (hash, cf, GUINT_TO_POINTER (nf++));
1355
1356       scf++;
1357     }
1358
1359     gts_split_expand (vs, ps->s, ps->s->edge_class);
1360   }
1361
1362   gts_surface_foreach_vertex (ps->s,
1363                               (GtsFunc) gts_object_reset_reserved, NULL);
1364   g_hash_table_destroy (hash);
1365 }
1366
1367 static guint surface_read (GtsSurface * surface,
1368                            GtsFile * f,
1369                            GPtrArray * vertices,
1370                            GPtrArray * faces)
1371 {
1372   GtsEdge ** edges;
1373   guint n, nv, ne, nf;
1374
1375   g_return_val_if_fail (surface != NULL, 1);
1376   g_return_val_if_fail (f != NULL, 1);
1377
1378   if (f->type != GTS_INT) {
1379     gts_file_error (f, "expecting an integer (number of vertices)");
1380     return f->line;
1381   }
1382   nv = atoi (f->token->str);
1383
1384   gts_file_next_token (f);
1385   if (f->type != GTS_INT) {
1386     gts_file_error (f, "expecting an integer (number of edges)");
1387     return f->line;
1388   }
1389   ne = atoi (f->token->str);
1390
1391   gts_file_next_token (f);
1392   if (f->type != GTS_INT) {
1393     gts_file_error (f, "expecting an integer (number of faces)");
1394     return f->line;
1395   }
1396   nf = atoi (f->token->str);
1397
1398   gts_file_next_token (f);
1399   if (f->type == GTS_STRING) {
1400     if (f->type != GTS_STRING) {
1401       gts_file_error (f, "expecting a string (GtsSurfaceClass)");
1402       return f->line;
1403     }
1404     gts_file_next_token (f);
1405     if (f->type != GTS_STRING) {
1406       gts_file_error (f, "expecting a string (GtsFaceClass)");
1407       return f->line;
1408     }
1409     gts_file_next_token (f);
1410     if (f->type != GTS_STRING) {
1411       gts_file_error (f, "expecting a string (GtsEdgeClass)");
1412       return f->line;
1413     }
1414     gts_file_next_token (f);
1415     if (f->type != GTS_STRING) {
1416       gts_file_error (f, "expecting a string (GtsVertexClass)");
1417       return f->line;
1418     }
1419     if (!strcmp (f->token->str, "GtsVertexBinary"))
1420       GTS_POINT_CLASS (surface->vertex_class)->binary = TRUE;
1421     else
1422       gts_file_first_token_after (f, '\n');
1423   }
1424   else
1425     gts_file_first_token_after (f, '\n');
1426
1427   g_ptr_array_set_size (vertices, nv);
1428   g_ptr_array_set_size (faces, nf);
1429   /* allocate nv + 1 just in case nv == 0 */
1430   edges = g_malloc ((ne + 1)*sizeof (GtsEdge *));
1431
1432   n = 0;
1433   while (n < nv && f->type != GTS_ERROR) {
1434     GtsObject * new_vertex =
1435       gts_object_new (GTS_OBJECT_CLASS (surface->vertex_class));
1436
1437     (* GTS_OBJECT_CLASS (surface->vertex_class)->read) (&new_vertex, f);
1438     if (f->type != GTS_ERROR) {
1439       if (!GTS_POINT_CLASS (surface->vertex_class)->binary)
1440         gts_file_first_token_after (f, '\n');
1441       g_ptr_array_index (vertices, n++) = new_vertex;
1442     }
1443     else
1444       gts_object_destroy (new_vertex);
1445   }
1446   if (f->type == GTS_ERROR)
1447     nv = n;
1448   if (GTS_POINT_CLASS (surface->vertex_class)->binary)
1449     gts_file_first_token_after (f, '\n');
1450
1451   n = 0;
1452   while (n < ne && f->type != GTS_ERROR) {
1453     guint p1, p2;
1454
1455     if (f->type != GTS_INT)
1456       gts_file_error (f, "expecting an integer (first vertex index)");
1457     else {
1458       p1 = atoi (f->token->str);
1459       if (p1 == 0 || p1 > nv)
1460         gts_file_error (f, "vertex index `%d' is out of range `[1,%d]'",
1461                         p1, nv);
1462       else {
1463         gts_file_next_token (f);
1464         if (f->type != GTS_INT)
1465           gts_file_error (f, "expecting an integer (second vertex index)");
1466         else {
1467           p2 = atoi (f->token->str);
1468           if (p2 == 0 || p2 > nv)
1469             gts_file_error (f, "vertex index `%d' is out of range `[1,%d]'",
1470                             p2, nv);
1471           else {
1472             GtsEdge * new_edge =
1473               gts_edge_new (surface->edge_class,
1474                             g_ptr_array_index (vertices, p1 - 1),
1475                             g_ptr_array_index (vertices, p2 - 1));
1476
1477             gts_file_next_token (f);
1478             if (f->type != '\n')
1479               if (GTS_OBJECT_CLASS (surface->edge_class)->read)
1480                 (*GTS_OBJECT_CLASS (surface->edge_class)->read)
1481                   ((GtsObject **) &new_edge, f);
1482             gts_file_first_token_after (f, '\n');
1483             edges[n++] = new_edge;
1484           }
1485         }
1486       }
1487     }
1488   }
1489   if (f->type == GTS_ERROR)
1490     ne = n;
1491
1492   n = 0;
1493   while (n < nf && f->type != GTS_ERROR) {
1494     guint s1, s2, s3;
1495
1496     if (f->type != GTS_INT)
1497       gts_file_error (f, "expecting an integer (first edge index)");
1498     else {
1499       s1 = atoi (f->token->str);
1500       if (s1 == 0 || s1 > ne)
1501         gts_file_error (f, "edge index `%d' is out of range `[1,%d]'",
1502                         s1, ne);
1503       else {
1504         gts_file_next_token (f);
1505         if (f->type != GTS_INT)
1506           gts_file_error (f, "expecting an integer (second edge index)");
1507         else {
1508           s2 = atoi (f->token->str);
1509           if (s2 == 0 || s2 > ne)
1510             gts_file_error (f, "edge index `%d' is out of range `[1,%d]'",
1511                             s2, ne);
1512           else {
1513             gts_file_next_token (f);
1514             if (f->type != GTS_INT)
1515               gts_file_error (f, "expecting an integer (third edge index)");
1516             else {
1517               s3 = atoi (f->token->str);
1518               if (s3 == 0 || s3 > ne)
1519                 gts_file_error (f, "edge index `%d' is out of range `[1,%d]'",
1520                                 s3, ne);
1521               else {
1522                 GtsFace * new_face = gts_face_new (surface->face_class,
1523                                                    edges[s1 - 1],
1524                                                    edges[s2 - 1],
1525                                                    edges[s3 - 1]);
1526
1527                 gts_file_next_token (f);
1528                 if (f->type != '\n')
1529                   if (GTS_OBJECT_CLASS (surface->face_class)->read)
1530                     (*GTS_OBJECT_CLASS (surface->face_class)->read)
1531                       ((GtsObject **) &new_face, f);
1532                 gts_file_first_token_after (f, '\n');
1533                 gts_surface_add_face (surface, new_face);
1534                 g_ptr_array_index (faces, n++) = new_face;
1535               }
1536             }
1537           }
1538         }
1539       }
1540     }
1541   }
1542
1543   g_free (edges);
1544
1545   if (f->type == GTS_ERROR) {
1546     gts_allow_floating_vertices = TRUE;
1547     while (nv)
1548       gts_object_destroy (GTS_OBJECT (g_ptr_array_index (vertices, nv-- - 1)));
1549     gts_allow_floating_vertices = FALSE;
1550     return f->line;
1551   }
1552
1553   return 0;
1554 }
1555
1556 /**
1557  * gts_psurface_open:
1558  * @klass: a #GtsPSurfaceClass.
1559  * @s: a #GtsSurface.
1560  * @split_class: a #GtsSplitClass to use for the #GtsSplit.
1561  * @f: a #GtsFile.
1562  *
1563  * Creates a new #GtsPSurface prepared for input from the file @f
1564  * containing a valid GTS representation of a progressive surface. The initial
1565  * shape of the progressive surface is loaded into @s.
1566  *
1567  * Before being usable as such this progressive surface must be closed using
1568  * gts_psurface_close(). While open however, the functions
1569  * gts_psurface_get_vertex_number(), gts_psurface_min_vertex_number() and
1570  * gts_psurface_max_vertex_number() can still be used.
1571  *
1572  * Returns: a new #GtsPSurface or %NULL if there was a format error while
1573  * reading the file, in which case @f contains information about the error.
1574  */
1575 GtsPSurface * gts_psurface_open (GtsPSurfaceClass * klass,
1576                                  GtsSurface * s,
1577                                  GtsSplitClass * split_class,
1578                                  GtsFile * f)
1579 {
1580   GtsPSurface * ps;
1581
1582   g_return_val_if_fail (klass != NULL, NULL);
1583   g_return_val_if_fail (s != NULL, NULL);
1584   g_return_val_if_fail (split_class != NULL, NULL);
1585   g_return_val_if_fail (f != NULL, NULL);
1586
1587   ps = GTS_PSURFACE (gts_object_new (GTS_OBJECT_CLASS (klass)));
1588   ps->s = s;
1589   ps->split_class = split_class;
1590
1591   ps->vertices = g_ptr_array_new ();
1592   ps->faces = g_ptr_array_new ();
1593
1594   if (surface_read (s, f, ps->vertices, ps->faces)) {
1595     ps->s = NULL;
1596     gts_object_destroy (GTS_OBJECT (ps));
1597     return NULL;
1598   }
1599
1600   ps->min = gts_surface_vertex_number (ps->s);
1601   ps->pos = 0;
1602
1603   if (f->type == GTS_INT) {
1604     gint ns = atoi (f->token->str);
1605
1606     if (ns > 0) {
1607       g_ptr_array_set_size (ps->split, ns);
1608       gts_file_first_token_after (f, '\n');
1609     }
1610   }
1611
1612   return ps;
1613 }
1614
1615 /**
1616  * gts_psurface_read_vertex:
1617  * @ps: a #GtsPSurface prealably created with gts_psurface_open().
1618  * @fp: a #GtsFile.
1619  *
1620  * Reads in one vertex split operation from @fp and performs the expansion.
1621  *
1622  * If an error occurs while reading the file, the @error field of @fp is set.
1623  *
1624  * Returns: the newly created #GtsSplit or %NULL if no vertex split could be
1625  * read from @fp.
1626  */
1627 GtsSplit * gts_psurface_read_vertex (GtsPSurface * ps, GtsFile * fp)
1628 {
1629   guint nv, ncf;
1630   GtsSplit * vs, * parent;
1631   GtsSplitCFace * scf;
1632
1633   g_return_val_if_fail (ps != NULL, NULL);
1634   g_return_val_if_fail (fp != NULL, NULL);
1635   g_return_val_if_fail (!GTS_PSURFACE_IS_CLOSED (ps), NULL);
1636
1637   if (ps->pos >= ps->split->len)
1638     return NULL;
1639
1640   if (fp->type == GTS_NONE)
1641     return NULL;
1642   if (fp->type != GTS_INT) {
1643     gts_file_error (fp, "expecting an integer (vertex index)");
1644     return NULL;
1645   }
1646   nv = atoi (fp->token->str);
1647   if (nv == 0 || nv > ps->vertices->len) {
1648     gts_file_error (fp, "vertex index `%d' is out of range `[1,%d]'",
1649                     nv, ps->vertices->len);
1650     return NULL;
1651   }
1652
1653   gts_file_next_token (fp);
1654   if (fp->type != GTS_INT) {
1655     gts_file_error (fp, "expecting an integer (ncf)");
1656     return NULL;
1657   }
1658   ncf = atoi (fp->token->str);
1659
1660   vs = GTS_SPLIT (gts_object_new (GTS_OBJECT_CLASS (ps->split_class)));
1661
1662   vs->v = g_ptr_array_index (ps->vertices, nv - 1);
1663   vs->v1 = vs->v2 = NULL;
1664   vs->cfaces = NULL;
1665   vs->ncf = 0;
1666
1667   gts_file_next_token (fp);
1668   if (fp->type != '\n')
1669     if (GTS_OBJECT (vs)->klass->read)
1670       (* GTS_OBJECT (vs)->klass->read) ((GtsObject **) &vs, fp);
1671   gts_file_first_token_after (fp, '\n');
1672
1673   if (fp->type != GTS_ERROR) {
1674     vs->v1 = gts_object_new (GTS_OBJECT_CLASS (ps->s->vertex_class));
1675     (* GTS_OBJECT_CLASS (ps->s->vertex_class)->read) (&(vs->v1), fp);
1676     if (fp->type != GTS_ERROR) {
1677       vs->v1->reserved = vs;
1678       g_ptr_array_add (ps->vertices, vs->v1);
1679
1680       gts_file_first_token_after (fp, '\n');
1681
1682       vs->v2 = gts_object_new (GTS_OBJECT_CLASS (ps->s->vertex_class));
1683       (*GTS_OBJECT_CLASS (ps->s->vertex_class)->read) (&(vs->v2), fp);
1684       if (fp->type != GTS_ERROR) {
1685         vs->v2->reserved = vs;
1686         g_ptr_array_add (ps->vertices, vs->v2);
1687         gts_file_first_token_after (fp, '\n');
1688       }
1689     }
1690   }
1691
1692   if (fp->type != GTS_ERROR) {
1693     scf = vs->cfaces = g_malloc (sizeof (GtsSplitCFace)*ncf);
1694     while (fp->type != GTS_ERROR && ncf--) {
1695       guint it, flags;
1696       GtsFace * f;
1697       CFace * cf;
1698       GPtrArray * a;
1699
1700       if (fp->type != GTS_INT)
1701         gts_file_error (fp, "expecting an integer (face index)");
1702       else {
1703         it = atoi (fp->token->str);
1704         if (it == 0 || it > ps->faces->len)
1705           gts_file_error (fp, "face index `%d' is out of range `[1,%d]'",
1706                           it, ps->faces->len);
1707         else {
1708           gts_file_next_token (fp);
1709           if (fp->type != GTS_INT)
1710             gts_file_error (fp, "expecting an integer (flags)");
1711           else {
1712             flags = atoi (fp->token->str);
1713             f =
1714               GTS_FACE (gts_object_new (GTS_OBJECT_CLASS (ps->s->face_class)));
1715
1716             gts_file_next_token (fp);
1717             if (fp->type != '\n')
1718               if (GTS_OBJECT (f)->klass->read)
1719                 (*GTS_OBJECT (f)->klass->read) ((GtsObject **) &f, fp);
1720             gts_file_first_token_after (fp, '\n');
1721             if (fp->type != GTS_ERROR) {
1722               scf->f = f;
1723
1724               cf = (CFace *) f;
1725               GTS_OBJECT (cf)->klass = GTS_OBJECT_CLASS (cface_class ());
1726               cf->parent_split = vs;
1727               cf->t = g_ptr_array_index (ps->faces, it - 1);
1728               cf->flags = flags;
1729
1730               a = g_ptr_array_new ();
1731               do {
1732                 if (fp->type != GTS_INT)
1733                   gts_file_error (fp, "expecting an integer (face index)");
1734                 else {
1735                   it = atoi (fp->token->str);
1736                   if (it > ps->faces->len)
1737                     gts_file_error (fp,
1738                                     "face index `%d' is out of range `[1,%d]'",
1739                                     it, ps->faces->len);
1740                   else {
1741                     g_ptr_array_add (a, g_ptr_array_index (ps->faces,
1742                                                            it - 1));
1743                     gts_file_next_token (fp);
1744                   }
1745                 }
1746               } while (fp->type != GTS_ERROR && fp->type != '\n');
1747               gts_file_first_token_after (fp, '\n');
1748               g_ptr_array_add (a, NULL);
1749               scf->a1 = (GtsTriangle **) a->pdata;
1750               g_ptr_array_free (a, FALSE);
1751
1752               if (fp->type != GTS_ERROR) {
1753                 a = g_ptr_array_new ();
1754                 do {
1755                   if (fp->type != GTS_INT)
1756                     gts_file_error (fp, "expecting an integer (face index)");
1757                   else {
1758                     it = atoi (fp->token->str);
1759                     if (it > ps->faces->len)
1760                       gts_file_error (fp,
1761                                    "face index `%d' is out of range `[1,%d]'",
1762                                       it, ps->faces->len);
1763                     else {
1764                       g_ptr_array_add (a, g_ptr_array_index (ps->faces,
1765                                                              it - 1));
1766                       gts_file_next_token (fp);
1767                     }
1768                   }
1769                 } while (fp->type != GTS_ERROR && fp->type != '\n');
1770                 gts_file_first_token_after (fp, '\n');
1771                 g_ptr_array_add (a, NULL);
1772                 scf->a2 = (GtsTriangle **) a->pdata;
1773                 g_ptr_array_free (a, FALSE);
1774
1775                 g_ptr_array_add (ps->faces, f);
1776
1777                 vs->ncf++;
1778                 scf++;
1779               }
1780             }
1781           }
1782         }
1783       }
1784     }
1785   }
1786
1787   if (fp->type != GTS_ERROR) {
1788     if ((parent = GTS_OBJECT (vs->v)->reserved)) {
1789       GTS_OBJECT (vs->v)->reserved = NULL;
1790       if (parent->v1 == GTS_OBJECT (vs->v))
1791         parent->v1 = GTS_OBJECT (vs);
1792       else {
1793         g_assert (parent->v2 == GTS_OBJECT (vs->v));
1794         parent->v2 = GTS_OBJECT (vs);
1795       }
1796     }
1797     g_ptr_array_index (ps->split, ps->pos++) = vs;
1798     gts_split_expand (vs, ps->s, ps->s->edge_class);
1799
1800     return vs;
1801   }
1802
1803   if (vs->v1) gts_object_destroy (vs->v1);
1804   if (vs->v2) gts_object_destroy (vs->v2);
1805   gts_object_destroy (GTS_OBJECT (vs));
1806
1807   return NULL;
1808 }
1809
1810 /**
1811  * gts_psurface_close:
1812  * @ps: a #GtsPSurface prealably created with gts_psurface_open().
1813  *
1814  * Closes a progressive surface.
1815  */
1816 void gts_psurface_close (GtsPSurface * ps)
1817 {
1818   g_return_if_fail (ps != NULL);
1819   g_return_if_fail (!GTS_PSURFACE_IS_CLOSED (ps));
1820
1821   g_ptr_array_free (ps->vertices, TRUE);
1822   g_ptr_array_free (ps->faces, TRUE);
1823   ps->faces = ps->vertices = NULL;
1824
1825   gts_surface_foreach_vertex (ps->s,
1826                               (GtsFunc) gts_object_reset_reserved, NULL);
1827   if (ps->pos > 0)
1828     g_ptr_array_set_size (ps->split, ps->pos);
1829   if (ps->split->len > 1) {
1830     guint i, half = ps->split->len/2, n = ps->split->len - 1;
1831
1832     for (i = 0; i < half; i++) {
1833       gpointer p1 = g_ptr_array_index (ps->split, i);
1834       gpointer p2 = g_ptr_array_index (ps->split, n - i);
1835       g_ptr_array_index (ps->split, n - i) = p1;
1836       g_ptr_array_index (ps->split, i) = p2;
1837     }
1838   }
1839   ps->pos = 0;
1840 }