gts/vopt.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 "gts.h"
  21
  22 /* #define DEBUG_VOPT */
  23
  24 /* compute the normal (nx, ny, nz) as the cross-product of the first two
  25    oriented edges and the norm nt = |t| as (v1xv2).v3 */
  26 static void triangle_normal (GtsTriangle * t,
  27                              gdouble * nx,
  28                              gdouble * ny,
  29                              gdouble * nz,
  30                              gdouble * nt)
  31 {
  32   GtsPoint * p1, * p2 = NULL, * p3 = NULL;
  33   gdouble x1, y1, z1, x2, y2, z2;
  34
  35   g_return_if_fail (t != NULL);
  36
  37   p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
  38   if (GTS_SEGMENT (t->e1)->v1 == GTS_SEGMENT (t->e2)->v1) {
  39     p2 = GTS_POINT (GTS_SEGMENT (t->e2)->v2);
  40     p3 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
  41   }
  42   else if (GTS_SEGMENT (t->e1)->v2 == GTS_SEGMENT (t->e2)->v2) {
  43     p2 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
  44     p3 = GTS_POINT (GTS_SEGMENT (t->e2)->v1);
  45   }
  46   else if (GTS_SEGMENT (t->e1)->v1 == GTS_SEGMENT (t->e2)->v2) {
  47     p2 = GTS_POINT (GTS_SEGMENT (t->e2)->v1);
  48     p3 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
  49   }
  50   else if (GTS_SEGMENT (t->e1)->v2 == GTS_SEGMENT (t->e2)->v1) {
  51     p2 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
  52     p3 = GTS_POINT (GTS_SEGMENT (t->e2)->v2);
  53   }
  54   else
  55     g_assert_not_reached ();
  56
  57   x1 = p2->x - p1->x;
  58   y1 = p2->y - p1->y;
  59   z1 = p2->z - p1->z;
  60
  61   x2 = p3->x - p1->x;
  62   y2 = p3->y - p1->y;
  63   z2 = p3->z - p1->z;
  64
  65   *nt = ((p1->y*p2->z - p1->z*p2->y)*p3->x +
  66          (p1->z*p2->x - p1->x*p2->z)*p3->y +
  67          (p1->x*p2->y - p1->y*p2->x)*p3->z);
  68   *nx = y1*z2 - z1*y2;
  69   *ny = z1*x2 - x1*z2;
  70   *nz = x1*y2 - y1*x2;
  71 }
  72
  73 static void boundary_preservation (GtsEdge * edge,
  74                                    GtsFace * f,
  75                                    GtsVector e1, GtsVector e2,
  76                                    GtsMatrix * H, GtsVector c)
  77 {
  78   GtsTriangle * t = GTS_TRIANGLE (f);
  79   GtsEdge * edge2;
  80   GtsVertex * v1 = GTS_SEGMENT (edge)->v1, * v2 = GTS_SEGMENT (edge)->v2;
  81   GtsPoint * p1, * p2;
  82   GtsVector e, e3;
  83
  84   /* find orientation of segment */
  85   edge2 = edge == t->e1 ? t->e2 : edge == t->e2 ? t->e3 : t->e1;
  86   if (v2 != GTS_SEGMENT (edge2)->v1 && v2 != GTS_SEGMENT (edge2)->v2) {
  87     v2 = v1; v1 = GTS_SEGMENT (edge)->v2;
  88   }
  89   p1 = GTS_POINT (v1);
  90   p2 = GTS_POINT (v2);
  91
  92   e[0] = p2->x - p1->x;
  93   e[1] = p2->y - p1->y;
  94   e[2] = p2->z - p1->z;
  95
  96   e1[0] += e[0];
  97   e1[1] += e[1];
  98   e1[2] += e[2];
  99
 100   e3[0] = p2->y*p1->z - p2->z*p1->y;
 101   e3[1] = p2->z*p1->x - p2->x*p1->z;
 102   e3[2] = p2->x*p1->y - p2->y*p1->x;
 103
 104   e2[0] += e3[0];
 105   e2[1] += e3[1];
 106   e2[2] += e3[2];
 107
 108   H[0][0] += e[1]*e[1] + e[2]*e[2];
 109   H[0][1] -= e[0]*e[1];
 110   H[0][2] -= e[0]*e[2];
 111   H[1][0] = H[0][1];
 112   H[1][1] += e[0]*e[0] + e[2]*e[2];
 113   H[1][2] -= e[1]*e[2];
 114   H[2][0] = H[0][2];
 115   H[2][1] = H[1][2];
 116   H[2][2] += e[0]*e[0] + e[1]*e[1];
 117
 118   c[0] += e[1]*e3[2] - e[2]*e3[1];
 119   c[1] += e[2]*e3[0] - e[0]*e3[2];
 120   c[2] += e[0]*e3[1] - e[1]*e3[0];
 121 }
 122
 123 static gdouble boundary_cost (GtsEdge * edge,
 124                               GtsFace * f,
 125                               GtsVertex * v)
 126 {
 127   GtsTriangle * t = GTS_TRIANGLE (f);
 128   GtsEdge * edge2;
 129   GtsVertex * v1 = GTS_SEGMENT (edge)->v1, * v2 = GTS_SEGMENT (edge)->v2;
 130   GtsPoint * p1, * p2;
 131   GtsVector e;
 132   GtsPoint * p = GTS_POINT (v);
 133
 134   /* find orientation of segment */
 135   edge2 = edge == t->e1 ? t->e2 : edge == t->e2 ? t->e3 : t->e1;
 136   if (v2 != GTS_SEGMENT (edge2)->v1 && v2 != GTS_SEGMENT (edge2)->v2) {
 137     v2 = v1; v1 = GTS_SEGMENT (edge)->v2;
 138   }
 139   p1 = GTS_POINT (v1);
 140   p2 = GTS_POINT (v2);
 141
 142   e[0] = (p2->y - p1->y)*(p->z - p2->z) - (p2->z - p1->z)*(p->y - p2->y);
 143   e[1] = (p2->z - p1->z)*(p->x - p2->x) - (p2->x - p1->x)*(p->z - p2->z);
 144   e[2] = (p2->x - p1->x)*(p->y - p2->y) - (p2->y - p1->y)*(p->x - p2->x);
 145
 146   return e[0]*e[0] + e[1]*e[1] + e[2]*e[2];
 147 }
 148
 149 static gdouble edge_boundary_cost (GtsEdge * e, GtsVertex * v)
 150 {
 151   gdouble cost = 0.;
 152   GSList * i;
 153
 154   i = GTS_SEGMENT (e)->v1->segments;
 155   while (i) {
 156     GtsFace * f;
 157     if (GTS_IS_EDGE (i->data) &&
 158         (f = gts_edge_is_boundary (i->data, NULL)))
 159       cost += boundary_cost (i->data, f, v);
 160     i = i->next;
 161   }
 162   i = GTS_SEGMENT (e)->v2->segments;
 163   while (i) {
 164     GtsFace * f;
 165     if (i->data != e &&
 166         GTS_IS_EDGE (i->data) &&
 167         (f = gts_edge_is_boundary (i->data, NULL)))
 168       cost += boundary_cost (i->data, f, v);
 169     i = i->next;
 170   }
 171
 172   return cost/4.;
 173 }
 174
 175 static gdouble edge_volume_cost (GtsEdge * e, GtsVertex * v)
 176 {
 177   GSList * i, * triangles;
 178   gdouble n1, n2, n3, nt;
 179   gdouble cost = 0.0, a;
 180
 181   triangles = gts_vertex_triangles (GTS_SEGMENT (e)->v1, NULL);
 182   triangles = gts_vertex_triangles (GTS_SEGMENT (e)->v2, triangles);
 183
 184   i = triangles;
 185   while (i) {
 186     if (GTS_IS_FACE (i->data)) {
 187       triangle_normal (i->data, &n1, &n2, &n3, &nt);
 188       a = GTS_POINT (v)->x*n1 +
 189         GTS_POINT (v)->y*n2 +
 190         GTS_POINT (v)->z*n3 - nt;
 191       cost += a*a;
 192     }
 193     i = i->next;
 194   }
 195   g_slist_free (triangles);
 196
 197   return cost/36.;
 198 }
 199
 200 static gdouble edge_shape_cost (GtsEdge * e, GtsVertex * v)
 201 {
 202   GSList * list, * i;
 203   GtsVertex
 204     * v1 = GTS_SEGMENT (e)->v1,
 205     * v2 = GTS_SEGMENT (e)->v2;
 206   gdouble cost = 0.;
 207
 208   list = gts_vertex_neighbors (v1, NULL, NULL);
 209   list = gts_vertex_neighbors (v2, list, NULL);
 210   i = list;
 211   while (i) {
 212     GtsPoint * p = i->data;
 213     if (p != GTS_POINT (v1) && p != GTS_POINT (v2))
 214       cost += gts_point_distance2 (p, GTS_POINT (v));
 215     i = i->next;
 216   }
 217   g_slist_free (list);
 218
 219   return cost;
 220 }
 221
 222 /**
 223  * gts_volume_optimized_vertex:
 224  * @edge: a #GtsEdge.
 225  * @klass: a #GtsVertexClass to be used for the new vertex.
 226  * @params: a #GtsVolumeOptimizedParms.
 227  *
 228  * Returns: a #GtsVertex which can be used to replace @edge for an
 229  * edge collapse operation. The position of the vertex is optimized in
 230  * order to minimize the changes in area and volume for the surface
 231  * using @edge. The volume enclosed by the surface is locally
 232  * preserved. For more details see "Fast and memory efficient
 233  * polygonal simplification" (1998) and "Evaluation of memoryless
 234  * simplification" (1999) by Lindstrom and Turk.
 235  */
 236 GtsVertex * gts_volume_optimized_vertex (GtsEdge * edge,
 237                                          GtsVertexClass * klass,
 238                                          GtsVolumeOptimizedParams * params)
 239 {
 240   GSList * triangles, * i;
 241   gdouble sn1 = 0., sn2 = 0., sn3 = 0.;
 242   gdouble sn11 = 0., sn22 = 0., sn33 = 0.;
 243   gdouble sn12 = 0., sn13 = 0., sn23 = 0.;
 244   gdouble st = 0., stn1 = 0., stn2 = 0., stn3 = 0.;
 245   gdouble n1, n2, n3, nt;
 246   GtsMatrix * A, * Ai;
 247   GtsVector A1, b;
 248   GtsVector e1 = {0., 0., 0.}, e2 = {0., 0., 0.};
 249   GtsMatrix * Hb;
 250   GtsVector cb = {0., 0., 0.};
 251   GtsVertex * v;
 252   GtsVertex * v1, * v2;
 253   guint n = 0, nb = 0;
 254 #ifdef DEBUG_VOPT
 255   guint nold = 0;
 256 #endif
 257
 258   g_return_val_if_fail (edge != NULL, NULL);
 259   g_return_val_if_fail (klass != NULL, NULL);
 260   g_return_val_if_fail (params != NULL, NULL);
 261
 262   A = gts_matrix_zero (NULL);
 263   Hb = gts_matrix_zero (NULL);
 264   v1 = GTS_SEGMENT (edge)->v1;
 265   v2 = GTS_SEGMENT (edge)->v2;
 266
 267   /* boundary preservation */
 268   i = v1->segments;
 269   while (i) {
 270     GtsEdge * edge1 = i->data;
 271     GtsFace * f;
 272     if (GTS_IS_EDGE (edge1) &&
 273         (f = gts_edge_is_boundary (edge1, NULL))) {
 274       boundary_preservation (edge1, f, e1, e2, Hb, cb);
 275       nb++;
 276     }
 277     i = i->next;
 278   }
 279   i = v2->segments;
 280   while (i) {
 281     GtsEdge * edge1 = i->data;
 282     GtsFace * f;
 283     if (edge1 != edge &&
 284         GTS_IS_EDGE (edge1) &&
 285         (f = gts_edge_is_boundary (edge1, NULL))) {
 286       boundary_preservation (edge1, f, e1, e2, Hb, cb);
 287       nb++;
 288     }
 289     i = i->next;
 290   }
 291   if (nb > 0) {
 292     GtsMatrix * H = gts_matrix_new (
 293                e1[2]*e1[2] + e1[1]*e1[1], - e1[0]*e1[1], - e1[0]*e1[2], 0.,
 294                - e1[0]*e1[1], e1[2]*e1[2] + e1[0]*e1[0], - e1[1]*e1[2], 0.,
 295                - e1[0]*e1[2], - e1[1]*e1[2], e1[1]*e1[1] + e1[0]*e1[0], 0.,
 296                0., 0., 0., 0.);
 297     GtsVector c;
 298
 299     c[0] = e1[1]*e2[2] - e1[2]*e2[1];
 300     c[1] = e1[2]*e2[0] - e1[0]*e2[2];
 301     c[2] = e1[0]*e2[1] - e1[1]*e2[0];
 302     n = gts_matrix_quadratic_optimization (A, b, n, H, c);
 303     gts_matrix_destroy (H);
 304   }
 305
 306   g_assert (n <= 2);
 307
 308 #ifdef DEBUG_VOPT
 309   if (n != nold) {
 310     fprintf (stderr, "--- boundary preservation ---\n");
 311     gts_matrix_print (A, stderr);
 312     gts_vector_print (b, stderr);
 313     nold = n;
 314   }
 315 #endif
 316
 317   /* volume preservation */
 318   triangles = gts_vertex_triangles (v1, NULL);
 319   triangles = gts_vertex_triangles (v2, triangles);
 320
 321   i = triangles;
 322   while (i) {
 323     if (GTS_IS_FACE (i->data)) {
 324       triangle_normal (i->data, &n1, &n2, &n3, &nt);
 325       sn1 += n1; sn2 += n2; sn3 += n3;
 326       sn11 += n1*n1; sn22 += n2*n2; sn33 += n3*n3;
 327       sn12 += n1*n2; sn13 += n1*n3; sn23 += n2*n3;
 328       st += nt; stn1 += nt*n1; stn2 += nt*n2; stn3 += nt*n3;
 329     }
 330     i = i->next;
 331   }
 332   g_slist_free (triangles);
 333
 334   A1[0] = sn1; A1[1] = sn2; A1[2] = sn3;
 335   n = gts_matrix_compatible_row (A, b, n, A1, st);
 336
 337 #ifdef DEBUG_VOPT
 338   if (n != nold) {
 339     fprintf (stderr, "--- volume preservation ---\n");
 340     gts_matrix_print (A, stderr);
 341     gts_vector_print (b, stderr);
 342     nold = n;
 343   }
 344 #endif
 345
 346 #if 1 /* Weighted average of volume and boundary optimization */
 347   if (n < 3) {
 348     /* volume optimization and boundary optimization */
 349     GtsMatrix * H = gts_matrix_new (sn11, sn12, sn13, 0.,
 350                                     sn12, sn22, sn23, 0.,
 351                                     sn13, sn23, sn33, 0.,
 352                                     0., 0., 0., 0.);
 353     GtsVector c;
 354     gdouble le = 9.*params->boundary_weight*
 355       gts_point_distance2 (GTS_POINT (v1),
 356                            GTS_POINT (v2));
 357     guint i, j;
 358
 359     c[0] = - stn1; c[1] = - stn2; c[2] = - stn3;
 360     if (nb > 0)
 361       for (i = 0; i < 3; i++) {
 362         for (j = 0; j < 3; j++)
 363           H[i][j] = params->volume_weight*H[i][j] + le*Hb[i][j];
 364         c[i] = params->volume_weight*c[i] + le*cb[i];
 365       }
 366     n = gts_matrix_quadratic_optimization (A, b, n, H, c);
 367     gts_matrix_destroy (H);
 368   }
 369
 370 #ifdef DEBUG_VOPT
 371   if (n != nold) {
 372     fprintf (stderr, "--- volume and boundary optimization ---\n");
 373     gts_matrix_print (A, stderr);
 374     gts_vector_print (b, stderr);
 375     nold = n;
 376   }
 377 #endif
 378
 379   if (n < 3) {
 380     /* triangle shape optimization */
 381     gdouble nv = 0.0;
 382     GtsMatrix * H;
 383     GtsVector c = {0., 0., 0.};
 384     GSList * list, * i;
 385
 386     list = gts_vertex_neighbors (v1, NULL, NULL);
 387     list = gts_vertex_neighbors (v2, list, NULL);
 388
 389     i = list;
 390     while (i) {
 391       GtsPoint * p1 = i->data;
 392       if (p1 != GTS_POINT (v1) && p1 != GTS_POINT (v2)) {
 393         nv += 1.0;
 394         c[0] -= p1->x;
 395         c[1] -= p1->y;
 396         c[2] -= p1->z;
 397       }
 398       i = i->next;
 399     }
 400     g_slist_free (list);
 401
 402     H = gts_matrix_new (nv, 0., 0., 0.,
 403                         0., nv, 0., 0.,
 404                         0., 0., nv, 0.,
 405                         0., 0., 0., 0.);
 406     n = gts_matrix_quadratic_optimization (A, b, n, H, c);
 407     gts_matrix_destroy (H);
 408   }
 409
 410 #ifdef DEBUG_VOPT
 411   if (n != nold) {
 412     fprintf (stderr, "--- triangle shape optimization ---\n");
 413     gts_matrix_print (A, stderr);
 414     gts_vector_print (b, stderr);
 415     nold = n;
 416   }
 417 #endif
 418 #else /* Weighted average of volume, boundary and shape optimization */
 419   if (n < 3) {
 420     /* volume optimization, boundary and shape optimization */
 421     GtsMatrix * H;
 422     GtsVector c;
 423     gdouble l2 = gts_point_distance2 (GTS_POINT (v1),
 424                                       GTS_POINT (v2));
 425     gdouble wv = params->volume_weight/32.;
 426     gdouble wb = params->boundary_weight/4.*l2;
 427     gdouble ws = params->shape_weight*l2*l2;
 428
 429     gdouble nv = 0.0;
 430     GtsVector cs = {0., 0., 0.};
 431     GSList * list, * i;
 432
 433     list = gts_vertex_neighbors (v1, NULL, NULL);
 434     list = gts_vertex_neighbors (v2, list, NULL);
 435
 436     i = list;
 437     while (i) {
 438       GtsPoint * p1 = i->data;
 439       if (p1 != GTS_POINT (v1) && p1 != GTS_POINT (v2)) {
 440         nv += 1.0;
 441         cs[0] -= p1->x;
 442         cs[1] -= p1->y;
 443         cs[2] -= p1->z;
 444       }
 445       i = i->next;
 446     }
 447     g_slist_free (list);
 448
 449     H = gts_matrix_new (wv*sn11 + wb*Hb[0][0] + ws*nv,
 450                         wv*sn12 + wb*Hb[0][1],
 451                         wv*sn13 + wb*Hb[0][2],
 452                         wv*sn12 + wb*Hb[1][0],
 453                         wv*sn22 + wb*Hb[1][1] + ws*nv,
 454                         wv*sn23 + wb*Hb[1][2],
 455                         wv*sn13 + wb*Hb[2][0],
 456                         wv*sn23 + wb*Hb[2][1],
 457                         wv*sn33 + wb*Hb[2][2] + ws*nv);
 458
 459     c[0] = - wv*stn1 + wb*cb[0] + ws*cs[0];
 460     c[1] = - wv*stn2 + wb*cb[1] + ws*cs[1];
 461     c[2] = - wv*stn3 + wb*cb[2] + ws*cs[2];
 462
 463     n = gts_matrix_quadratic_optimization (A, b, n, H, c);
 464     gts_matrix_destroy (H);
 465   }
 466
 467 #ifdef DEBUG_VOPT
 468   if (n != nold) {
 469     fprintf (stderr, "--- volume, boundary and shape optimization ---\n");
 470     gts_matrix_print (A, stderr);
 471     gts_vector_print (b, stderr);
 472     nold = n;
 473   }
 474 #endif
 475 #endif /* Weighted average of volume, boundary and shape optimization */
 476
 477   g_assert (n == 3);
 478   Ai = gts_matrix3_inverse (A);
 479   g_assert (Ai != NULL);
 480
 481   v = gts_vertex_new (klass,
 482                       Ai[0][0]*b[0] + Ai[0][1]*b[1] + Ai[0][2]*b[2],
 483                       Ai[1][0]*b[0] + Ai[1][1]*b[1] + Ai[1][2]*b[2],
 484                       Ai[2][0]*b[0] + Ai[2][1]*b[1] + Ai[2][2]*b[2]);
 485
 486   gts_matrix_destroy (A);
 487   gts_matrix_destroy (Ai);
 488   gts_matrix_destroy (Hb);
 489
 490   return v;
 491 }
 492
 493 /**
 494  * gts_volume_optimized_cost:
 495  * @e: a #GtsEdge.
 496  * @params: a #GtsVolumeOptimizedParams.
 497  *
 498  * Returns: the cost for the collapse of @e as minimized by the function
 499  * gts_volume_optimized_vertex().
 500  */
 501 gdouble gts_volume_optimized_cost (GtsEdge * e,
 502                                    GtsVolumeOptimizedParams * params)
 503 {
 504   GtsVertex * v;
 505   gdouble cost;
 506   gdouble length2;
 507
 508   g_return_val_if_fail (e != NULL, G_MAXDOUBLE);
 509   g_return_val_if_fail (params != NULL, G_MAXDOUBLE);
 510
 511   v = gts_volume_optimized_vertex (e, gts_vertex_class (), params);
 512
 513   length2 = gts_point_distance2 (GTS_POINT (GTS_SEGMENT (e)->v1),
 514                                  GTS_POINT (GTS_SEGMENT (e)->v2));
 515   cost =
 516     params->volume_weight*edge_volume_cost (e, v) +
 517     params->boundary_weight*length2*edge_boundary_cost (e, v) +
 518     params->shape_weight*length2*length2*edge_shape_cost (e, v);
 519   gts_object_destroy (GTS_OBJECT (v));
 520
 521   return cost;
 522 }