gts/psurface.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 <math.h>
  22 #include "gts.h"
  23
  24 #define HEAP_INSERT_OBJECT(h, e) (GTS_OBJECT (e)->reserved =\
  25                                    gts_eheap_insert (h, e))
  26 #define HEAP_REMOVE_OBJECT(h, e) (gts_eheap_remove (h, GTS_OBJECT (e)->reserved),\
  27                                   GTS_OBJECT (e)->reserved = NULL)
  28
  29 static void psurface_destroy (GtsObject * object)
  30 {
  31   GtsPSurface * ps = GTS_PSURFACE (object);
  32   guint i;
  33
  34   if (!GTS_PSURFACE_IS_CLOSED (ps))
  35     gts_psurface_close (ps);
  36
  37   for (i = 0; i < ps->split->len; i++)
  38     if (g_ptr_array_index (ps->split, i))
  39       gts_object_destroy (GTS_OBJECT (g_ptr_array_index (ps->split, i)));
  40   g_ptr_array_free (ps->split, TRUE);
  41
  42   (* GTS_OBJECT_CLASS (gts_psurface_class ())->parent_class->destroy) (object);
  43 }
  44
  45 static void psurface_class_init (GtsObjectClass * klass)
  46 {
  47   klass->destroy = psurface_destroy;
  48 }
  49
  50 static void psurface_init (GtsPSurface * psurface)
  51 {
  52   psurface->s = NULL;
  53   psurface->split = g_ptr_array_new ();
  54   psurface->split_class = gts_split_class ();
  55   psurface->pos = psurface->min = 0;
  56   psurface->vertices = psurface->faces = NULL;
  57 }
  58
  59 /**
  60  * gts_psurface_class:
  61  *
  62  * Returns: the #GtsPSurfaceClass.
  63  */
  64 GtsPSurfaceClass * gts_psurface_class (void)
  65 {
  66   static GtsPSurfaceClass * klass = NULL;
  67
  68   if (klass == NULL) {
  69     GtsObjectClassInfo psurface_info = {
  70       "GtsPSurface",
  71       sizeof (GtsPSurface),
  72       sizeof (GtsPSurfaceClass),
  73       (GtsObjectClassInitFunc) psurface_class_init,
  74       (GtsObjectInitFunc) psurface_init,
  75       (GtsArgSetFunc) NULL,
  76       (GtsArgGetFunc) NULL
  77     };
  78     klass = gts_object_class_new (gts_object_class (),
  79                                   &psurface_info);
  80   }
  81
  82   return klass;
  83 }
  84
  85 static GtsVertex * edge_collapse (GtsPSurface * ps,
  86                                   GtsEdge * e,
  87                                   GtsEHeap * heap,
  88                                   GtsCoarsenFunc coarsen_func,
  89                                   gpointer coarsen_data,
  90                                   gdouble maxcosine2)
  91 {
  92   GtsVertex  * v1 = GTS_SEGMENT (e)->v1, * v2 = GTS_SEGMENT (e)->v2, * mid;
  93   GtsSplit * vs;
  94   GtsObject * o1, * o2;
  95
  96   /* if the edge is degenerate (i.e. v1 == v2), destroy and return */
  97   if (v1 == v2) {
  98     gts_object_destroy (GTS_OBJECT (e));
  99     return NULL;
 100   }
 101
 102   if (!gts_edge_collapse_is_valid (e) ||
 103       /* check that a non-manifold edge is not a contact edge */
 104       (g_slist_length (e->triangles) > 2 && gts_edge_is_contact (e) > 1)) {
 105     GTS_OBJECT (e)->reserved =
 106       gts_eheap_insert_with_key (heap, e, G_MAXDOUBLE);
 107     return NULL;
 108   }
 109
 110   mid = (*coarsen_func) (e, ps->s->vertex_class, coarsen_data);
 111
 112   if (gts_edge_collapse_creates_fold (e, mid, maxcosine2)) {
 113     GTS_OBJECT (e)->reserved =
 114       gts_eheap_insert_with_key (heap, e, G_MAXDOUBLE);
 115     gts_object_destroy (GTS_OBJECT (mid));
 116     return NULL;
 117   }
 118
 119   if (GTS_OBJECT (v1)->reserved)
 120     o1 = GTS_OBJECT (v1)->reserved;
 121   else
 122     o1 = GTS_OBJECT (v1);
 123   if (GTS_OBJECT (v2)->reserved)
 124     o2 = GTS_OBJECT (v2)->reserved;
 125   else
 126     o2 = GTS_OBJECT (v2);
 127   vs = gts_split_new (ps->split_class, mid, o1, o2);
 128   gts_split_collapse (vs, ps->s->edge_class, heap);
 129   GTS_OBJECT (vs->v)->reserved = vs;
 130   g_ptr_array_add (ps->split, vs);
 131
 132   return mid;
 133 }
 134
 135 static void update_2nd_closest_neighbors (GtsVertex * v, GtsEHeap * heap)
 136 {
 137   GSList * i = v->segments;
 138   GSList * list = NULL;
 139
 140   while (i) {
 141     GtsSegment * s = i->data;
 142     if (GTS_IS_EDGE (s)) {
 143       GtsVertex * v1 = s->v1 == v ? s->v2 : s->v1;
 144       GSList * j = v1->segments;
 145       while (j) {
 146         GtsSegment * s1 = j->data;
 147         if (GTS_IS_EDGE (s1) && !g_slist_find (list, s1))
 148           list = g_slist_prepend (list, s1);
 149         j = j->next;
 150       }
 151     }
 152     i = i->next;
 153   }
 154
 155   i = list;
 156   while (i) {
 157     GtsEdge * e = i->data;
 158     if (GTS_OBJECT (e)->reserved)
 159       HEAP_REMOVE_OBJECT (heap, e);
 160     HEAP_INSERT_OBJECT (heap, e);
 161     i = i->next;
 162   }
 163
 164   g_slist_free (list);
 165 }
 166
 167 static gdouble edge_length2 (GtsEdge * e)
 168 {
 169   return gts_point_distance2 (GTS_POINT (GTS_SEGMENT (e)->v1),
 170                               GTS_POINT (GTS_SEGMENT (e)->v2));
 171 }
 172
 173 static void create_heap_coarsen (GtsEdge * e, GtsEHeap * heap)
 174 {
 175   HEAP_INSERT_OBJECT (heap, e);
 176 }
 177
 178 /* #define DEBUG_FOLD */
 179 /* #define DEBUG_CONTACT_VERTEX */
 180
 181 #ifdef DEBUG_FOLD
 182 static void check_fold (GtsTriangle * t, gdouble * maxcosine2)
 183 {
 184   GtsEdge * e1 = t->e1, * e2 = t->e2, * e3 = t->e3;
 185
 186
 187   if (gts_triangles_are_folded (e1->triangles,
 188                                 GTS_SEGMENT (e1)->v1,
 189                                 GTS_SEGMENT (e1)->v2,
 190                                 *maxcosine2) ||
 191       gts_triangles_are_folded (e2->triangles,
 192                                 GTS_SEGMENT (e2)->v1,
 193                                 GTS_SEGMENT (e2)->v2,
 194                                 *maxcosine2) ||
 195       gts_triangles_are_folded (e3->triangles,
 196                                 GTS_SEGMENT (e3)->v1,
 197                                 GTS_SEGMENT (e3)->v2,
 198                                 *maxcosine2)) {
 199     fprintf (stderr, "triangle %p:(%p,%p,%p) is folded\n", t, e1, e2, e3);
 200     g_assert_not_reached ();
 201   }
 202 }
 203 #endif
 204
 205 /**
 206  * gts_psurface_new:
 207  * @klass: a #GtsPSurfaceClass.
 208  * @surface: a #GtsSurface.
 209  * @split_class: a #GtsSplitClass to use for the new progressive surface.
 210  * @cost_func: cost function for the edge collapse algorithm.
 211  * @cost_data: data to pass to @cost_func.
 212  * @coarsen_func: the function returning the vertex replacement for the edge
 213  * collapse.
 214  * @coarsen_data: data to pass to @coarsen_func.
 215  * @stop_func: the function to call to decide whether to stop the coarsening
 216  * process.
 217  * @stop_data: data to pass to @stop_func.
 218  * @minangle: the minimum angle allowable between two neighboring triangles.
 219  * This is used to avoid introducing folds in the mesh during simplification.
 220  *
 221  * This function works in exactly the same way as the
 222  * gts_surface_coarsen() function, except that the history of edge
 223  * collapse is saved in an array of #GtsSplit objects. This allows for
 224  * dynamic continuous multiresolution control of the input @surface.
 225  *
 226  * Returns: a new progressive surface.
 227  */
 228 GtsPSurface * gts_psurface_new (GtsPSurfaceClass * klass,
 229                                 GtsSurface * surface,
 230                                 GtsSplitClass * split_class,
 231                                 GtsKeyFunc cost_func,
 232                                 gpointer cost_data,
 233                                 GtsCoarsenFunc coarsen_func,
 234                                 gpointer coarsen_data,
 235                                 GtsStopFunc stop_func,
 236                                 gpointer stop_data,
 237                                 gdouble minangle)
 238 {
 239   GtsPSurface * psurface;
 240   GtsEHeap * heap;
 241   GtsEdge * e;
 242   gdouble top_cost, maxcosine2;
 243   guint i;
 244
 245   g_return_val_if_fail (klass != NULL, NULL);
 246   g_return_val_if_fail (surface != NULL, NULL);
 247   g_return_val_if_fail (split_class != NULL, NULL);
 248   g_return_val_if_fail (stop_func != NULL, NULL);
 249
 250   psurface = GTS_PSURFACE (gts_object_new (GTS_OBJECT_CLASS (klass)));
 251   psurface->s = surface;
 252   psurface->split_class = split_class;
 253
 254   if (cost_func == NULL)
 255     cost_func = (GtsKeyFunc) edge_length2;
 256   if (coarsen_func == NULL)
 257     coarsen_func = (GtsCoarsenFunc) gts_segment_midvertex;
 258
 259   heap = gts_eheap_new (cost_func, cost_data);
 260   maxcosine2 = cos (minangle); maxcosine2 *= maxcosine2;
 261
 262   gts_eheap_freeze (heap);
 263   gts_surface_foreach_edge (surface, (GtsFunc) create_heap_coarsen, heap);
 264   gts_eheap_thaw (heap);
 265   /* we want to control edge destruction manually */
 266   gts_allow_floating_edges = TRUE;
 267   while ((e = gts_eheap_remove_top (heap, &top_cost)) &&
 268          (top_cost < G_MAXDOUBLE) &&
 269          !(*stop_func) (top_cost, gts_eheap_size (heap) -
 270                         gts_edge_face_number (e, surface), stop_data)) {
 271     GtsVertex * v = edge_collapse (psurface, e, heap,
 272                                    coarsen_func, coarsen_data, maxcosine2);
 273     if (v != NULL) {
 274       update_2nd_closest_neighbors (v, heap);
 275 #ifdef DEBUG_FOLD
 276       {
 277         GSList * triangles = gts_vertex_triangles (v, NULL), * i;
 278         fprintf (stderr, "\n---- Check for folds ----\n%p: ", v);
 279         i = triangles;
 280         while (i) {
 281           GtsTriangle * t = i->data;
 282           fprintf (stderr, "%p:(%p,%p,%p) ", t, t->e1, t->e2, t->e3);
 283           i = i->next;
 284         }
 285         fprintf (stderr, "\n");
 286         g_slist_free (triangles);
 287         gts_surface_foreach_face (surface, (GtsFunc) check_fold, &maxcosine2);
 288       }
 289 #endif
 290 #ifdef DEBUG_CONTACT_VERTEX
 291       if (gts_vertex_is_contact (v, FALSE) != 1) {
 292         FILE * fptr = fopen ("after", "wt");
 293         GSList * triangles = gts_vertex_triangles (v, NULL), * i;
 294
 295         fprintf (stderr, "collapse of %p created a contact vertex\n", e);
 296
 297         fprintf (fptr,
 298                  "(geometry \"sphere\" { = SPHERE 0.1 0. 0. 0. })\n"
 299                  "(normalization \"sphere\" none)\n");
 300         i = triangles;
 301         while (i) {
 302           gts_write_triangle (i->data, GTS_POINT (v), fptr);
 303           i = i->next;
 304         }
 305         g_assert_not_reached ();
 306       }
 307 #endif
 308     }
 309   }
 310   gts_allow_floating_edges = FALSE;
 311
 312   /* set reserved field of remaining edges back to NULL */
 313   if (e) GTS_OBJECT (e)->reserved = NULL;
 314   gts_eheap_foreach (heap, (GFunc) gts_object_reset_reserved, NULL);
 315
 316   gts_eheap_destroy (heap);
 317
 318   psurface->pos = psurface->split->len;
 319   psurface->min = gts_surface_vertex_number (psurface->s);
 320
 321   /* set reserved field of vertices (used to build the hierarchy)
 322      back to NULL */
 323   for (i = 0; i < psurface->split->len; i++) {
 324     GtsSplit * vs = g_ptr_array_index (psurface->split, i);
 325     gts_object_reset_reserved (GTS_OBJECT (vs->v));
 326   }
 327
 328   return psurface;
 329 }
 330
 331 /**
 332  * gts_psurface_add_vertex:
 333  * @ps: a #GtsPSurface.
 334  *
 335  * Adds a vertex to the progressive surface @ps by expanding the next
 336  * available #GtsSplit.
 337  *
 338  * Returns: the expanded #GtsSplit or %NULL if all the #GtsSplit have already
 339  * been expanded.
 340  */
 341 GtsSplit * gts_psurface_add_vertex (GtsPSurface * ps)
 342 {
 343   GtsSplit * vs;
 344
 345   g_return_val_if_fail (ps != NULL, NULL);
 346   g_return_val_if_fail (GTS_PSURFACE_IS_CLOSED (ps), NULL);
 347
 348   if (ps->pos == 0)
 349     return NULL;
 350
 351   vs = g_ptr_array_index (ps->split, --ps->pos);
 352   gts_split_expand (vs, ps->s, ps->s->edge_class);
 353
 354   return vs;
 355 }
 356
 357 /**
 358  * gts_psurface_remove_vertex:
 359  * @ps: a #GtsPSurface.
 360  *
 361  * Removes one vertex from the progressive surface @ps by collapsing the first
 362  * available #GtsSplit.
 363  *
 364  * Returns: the collapsed #GtsSplit or %NULL if all the #GtsSplit have already
 365  * been collapsed.
 366  */
 367 GtsSplit * gts_psurface_remove_vertex (GtsPSurface * ps)
 368 {
 369   GtsSplit * vs;
 370
 371   g_return_val_if_fail (ps != NULL, NULL);
 372   g_return_val_if_fail (GTS_PSURFACE_IS_CLOSED (ps), NULL);
 373
 374   if (ps->pos == ps->split->len)
 375     return NULL;
 376
 377   vs = g_ptr_array_index (ps->split, ps->pos++);
 378   gts_split_collapse (vs, ps->s->edge_class, NULL);
 379
 380   return vs;
 381 }
 382
 383 /**
 384  * gts_psurface_max_vertex_number:
 385  * @ps: a #GtsPSurface.
 386  *
 387  * Returns: the maximum number of vertices of @ps i.e. the number of vertices
 388  * if all the #GtsSplit were expanded.
 389  */
 390 guint gts_psurface_max_vertex_number (GtsPSurface * ps)
 391 {
 392   g_return_val_if_fail (ps != NULL, 0);
 393
 394   return ps->min + ps->split->len;
 395 }
 396
 397 /**
 398  * gts_psurface_min_vertex_number:
 399  * @ps: a #GtsPSurface.
 400  *
 401  * Returns: the minimum number of vertices of @ps i.e. the number of vertices
 402  * if all the #GtsSplit were collapsed.
 403  */
 404 guint gts_psurface_min_vertex_number (GtsPSurface * ps)
 405 {
 406   g_return_val_if_fail (ps != NULL, 0);
 407
 408   return ps->min;
 409 }
 410
 411 /**
 412  * gts_psurface_set_vertex_number:
 413  * @ps: a #GtsPSurface.
 414  * @n: a number of vertices.
 415  *
 416  * Performs the required number of collapses or expansions to set the number
 417  * of vertices of @ps to @n.
 418  */
 419 void gts_psurface_set_vertex_number (GtsPSurface * ps, guint n)
 420 {
 421   g_return_if_fail (ps != NULL);
 422   g_return_if_fail (GTS_PSURFACE_IS_CLOSED (ps));
 423
 424   n = ps->min + ps->split->len - n;
 425   while (ps->pos > n && gts_psurface_add_vertex (ps))
 426     ;
 427   while (ps->pos < n && gts_psurface_remove_vertex (ps))
 428     ;
 429 }
 430
 431 /**
 432  * gts_psurface_get_vertex_number:
 433  * @ps: a #GtsPSurface.
 434  *
 435  * Returns: the current number of vertices of @ps.
 436  */
 437 guint gts_psurface_get_vertex_number (GtsPSurface * ps)
 438 {
 439   g_return_val_if_fail (ps != NULL, 0);
 440
 441   if (!GTS_PSURFACE_IS_CLOSED (ps))
 442     return ps->min + ps->pos;
 443   else
 444     return ps->min + ps->split->len - ps->pos;
 445 }
 446
 447 /**
 448  * gts_psurface_foreach_vertex:
 449  * @ps: a #GtsPSurface.
 450  * @func: a function to call for each vertex of @ps.
 451  * @data: data to be passed to @func.
 452  *
 453  * Calls @func for each (potential) vertex of @ps, whether actually used
 454  * or not. The vertices are called in the order they were created during the
 455  * edge collapse operation.
 456  */
 457 void gts_psurface_foreach_vertex (GtsPSurface * ps,
 458                                   GtsFunc func,
 459                                   gpointer data)
 460 {
 461   guint i;
 462
 463   g_return_if_fail (ps != NULL);
 464   g_return_if_fail (func != NULL);
 465   g_return_if_fail (GTS_PSURFACE_IS_CLOSED (ps));
 466
 467   for (i = 0; i < ps->split->len; i++) {
 468     GtsSplit * vs = g_ptr_array_index (ps->split, i);
 469     (*func) (vs->v, data);
 470   }
 471 }