1 /* GTS - Library for the manipulation of triangulated surfaces
2 * Copyright (C) 1999 Stéphane Popinet
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.
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.
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.
24 * gts_vertex_encroaches_edge:
28 * Returns: %TRUE if @v is strictly contained in the diametral circle of @e,
31 gboolean
gts_vertex_encroaches_edge (GtsVertex
* v
, GtsEdge
* e
)
33 GtsPoint
* p
, * p1
, * p2
;
35 g_return_val_if_fail (v
!= NULL
, FALSE
);
36 g_return_val_if_fail (e
!= NULL
, FALSE
);
39 p1
= GTS_POINT (GTS_SEGMENT (e
)->v1
);
40 p2
= GTS_POINT (GTS_SEGMENT (e
)->v2
);
42 if ((p1
->x
- p
->x
)*(p2
->x
- p
->x
) + (p1
->y
- p
->y
)*(p2
->y
- p
->y
) < 0.0)
48 * gts_edge_is_encroached:
50 * @s: a #GtsSurface describing a (constrained) Delaunay triangulation.
51 * @encroaches: a #GtsEncroachFunc.
52 * @data: user data to be passed to @encroaches.
54 * Returns: a #GtsVertex belonging to @s and encroaching upon @e
55 * (as defined by @encroaches) or %NULL if there is none.
57 GtsVertex
* gts_edge_is_encroached (GtsEdge
* e
,
59 GtsEncroachFunc encroaches
,
64 g_return_val_if_fail (e
!= NULL
, NULL
);
65 g_return_val_if_fail (s
!= NULL
, NULL
);
66 g_return_val_if_fail (encroaches
!= NULL
, NULL
);
70 GtsFace
* f
= i
->data
;
71 if (GTS_IS_FACE (f
) && gts_face_has_parent_surface (f
, s
)) {
72 GtsVertex
* v
= gts_triangle_vertex_opposite (GTS_TRIANGLE (f
), e
);
73 if ((* encroaches
) (v
, e
, s
, data
))
82 #define ALREADY_ENCROACHED(c) (GTS_OBJECT (c)->reserved)
84 static void vertex_encroaches (GtsVertex
* v
,
87 GtsEncroachFunc encroaches
,
90 GSList
* triangles
, * i
;
92 g_return_if_fail (v
!= NULL
);
93 g_return_if_fail (surface
!= NULL
);
94 g_return_if_fail (encroached
!= NULL
);
95 g_return_if_fail (encroaches
!= NULL
);
97 i
= triangles
= gts_vertex_triangles (v
, NULL
);
99 GtsFace
* f
= i
->data
;
100 if (GTS_IS_FACE (f
) && gts_face_has_parent_surface (f
, surface
)) {
101 GtsEdge
* e
= gts_triangle_edge_opposite (i
->data
, v
);
102 if (!ALREADY_ENCROACHED (e
) &&
103 GTS_IS_CONSTRAINT (e
) &&
104 (* encroaches
) (v
, e
, surface
, data
)) {
105 gts_fifo_push (encroached
, e
);
106 ALREADY_ENCROACHED (e
) = encroached
;
111 g_slist_free (triangles
);
114 static void make_encroached_fifo (GtsEdge
* e
, gpointer
* datas
)
116 GtsFifo
* fifo
= datas
[0];
117 GtsSurface
* s
= datas
[1];
118 GtsEncroachFunc encroaches
= (GtsEncroachFunc
) datas
[2];
119 gpointer data
= datas
[3];
121 if (GTS_IS_CONSTRAINT (e
) &&
122 gts_edge_is_encroached (e
, s
, encroaches
, data
)) {
123 gts_fifo_push (fifo
, e
);
124 ALREADY_ENCROACHED (e
) = fifo
;
128 #define SQUARE_ROOT_TWO 1.41421356237309504880168872420969807856967187
129 #define DISTANCE_2D(v1, v2) (sqrt ((GTS_POINT (v2)->x - GTS_POINT (v1)->x)*\
130 (GTS_POINT (v2)->x - GTS_POINT (v1)->x) +\
131 (GTS_POINT (v2)->y - GTS_POINT (v1)->y)*\
132 (GTS_POINT (v2)->y - GTS_POINT (v1)->y)))
134 /* finds where to split the given edge to avoid infinite cycles. (see
135 Shewchuk's thesis for details */
136 static GtsVertex
* split_edge (GtsEdge
* e
,
137 GtsSurface
* surface
)
139 GSList
* i
= e
->triangles
;
142 /* look for constraints touching e */
144 GtsTriangle
* t
= i
->data
;
145 if (GTS_IS_FACE (t
) &&
146 gts_face_has_parent_surface (GTS_FACE (t
), surface
)) {
148 if (t
->e1
== e
) { e1
= t
->e2
; e2
= t
->e3
; }
149 else if (t
->e2
== e
) { e1
= t
->e1
; e2
= t
->e3
; }
150 else { e1
= t
->e1
; e2
= t
->e2
; }
151 if (GTS_IS_CONSTRAINT (e1
) && !GTS_IS_CONSTRAINT (e2
))
153 else if (GTS_IS_CONSTRAINT (e2
) && !GTS_IS_CONSTRAINT (e1
))
159 /* use power of two concentric shells */
160 GtsVertex
* v1
= GTS_SEGMENT (e
)->v1
;
161 GtsVertex
* v2
= GTS_SEGMENT (e
)->v2
;
162 gdouble l
= DISTANCE_2D (v1
, v2
);
163 gdouble nearestpower
= 1., split
;
165 while (l
> SQUARE_ROOT_TWO
*nearestpower
)
167 while (l
< SQUARE_ROOT_TWO
*nearestpower
/2.)
169 split
= nearestpower
/l
/2.;
171 if (GTS_SEGMENT (c
)->v1
== v2
|| GTS_SEGMENT (c
)->v2
== v2
)
173 return gts_vertex_new (surface
->vertex_class
,
174 (1. - split
)*GTS_POINT (v1
)->x
+
175 split
*GTS_POINT (v2
)->x
,
176 (1. - split
)*GTS_POINT (v1
)->y
+
177 split
*GTS_POINT (v2
)->y
,
178 (1. - split
)*GTS_POINT (v1
)->z
+
179 split
*GTS_POINT (v2
)->z
);
182 return gts_segment_midvertex (GTS_SEGMENT (e
), surface
->vertex_class
);
185 static gint
split_encroached (GtsSurface
* surface
,
186 GtsFifo
* encroached
,
188 GtsEncroachFunc encroaches
,
193 while (steiner_max
-- != 0 && (s
= gts_fifo_pop (encroached
))) {
194 GtsVertex
*add_vertex_returned
;
195 GtsVertex
* v
= split_edge (GTS_EDGE (s
), surface
);
196 GtsFace
* boundary
= gts_edge_is_boundary (GTS_EDGE (s
), surface
);
197 GtsFace
* f
= boundary
;
199 GtsEdge
* e1
= GTS_EDGE (gts_object_clone (GTS_OBJECT (s
)));
200 GtsEdge
* e2
= GTS_EDGE (gts_object_clone (GTS_OBJECT (s
)));
202 GTS_SEGMENT (e1
)->v1
= s
->v1
;
203 s
->v1
->segments
= g_slist_prepend (s
->v1
->segments
, e1
);
204 GTS_SEGMENT (e1
)->v2
= v
;
205 v
->segments
= g_slist_prepend (v
->segments
, e1
);
207 GTS_SEGMENT (e2
)->v1
= v
;
208 v
->segments
= g_slist_prepend (v
->segments
, e2
);
209 GTS_SEGMENT (e2
)->v2
= s
->v2
;
210 s
->v2
->segments
= g_slist_prepend (s
->v2
->segments
, e2
);
212 GtsEdge
* e1
= gts_edge_new (GTS_EDGE_CLASS (GTS_OBJECT (s
)->klass
),
214 GtsEdge
* e2
= gts_edge_new (GTS_EDGE_CLASS (GTS_OBJECT (s
)->klass
),
218 GTS_OBJECT (s
)->klass
= GTS_OBJECT_CLASS (surface
->edge_class
);
221 f
= gts_edge_has_parent_surface (GTS_EDGE (s
), surface
);
222 g_assert (f
!= NULL
);
223 add_vertex_returned
= gts_delaunay_add_vertex_to_face (surface
, v
, f
);
224 g_assert (add_vertex_returned
== NULL
);
227 gts_object_destroy (GTS_OBJECT (s
));
229 vertex_encroaches (v
, surface
, encroached
, encroaches
, data
);
231 if (gts_edge_is_encroached (e1
, surface
, encroaches
, data
)) {
232 gts_fifo_push (encroached
, e1
);
233 ALREADY_ENCROACHED (e1
) = encroached
;
235 if (gts_edge_is_encroached (e2
, surface
, encroaches
, data
)) {
236 gts_fifo_push (encroached
, e2
);
237 ALREADY_ENCROACHED (e2
) = encroached
;
245 * gts_delaunay_conform:
246 * @surface: a #GtsSurface describing a constrained Delaunay triangulation.
247 * @steiner_max: maximum number of Steiner points.
248 * @encroaches: a #GtsEncroachFunc.
249 * @data: user-data to pass to @encroaches.
251 * Recursively split constraints of @surface which are encroached by
252 * vertices of @surface (see Shewchuk 96 for details). The split
253 * constraints are destroyed and replaced by a set of new constraints
254 * of the same class. If gts_vertex_encroaches_edge() is used for
255 * @encroaches, the resulting surface will be Delaunay conforming.
257 * If @steiner_max is positive or nul, the recursive splitting
258 * procedure will stop when this maximum number of Steiner points is
259 * reached. In that case the resulting surface will not necessarily be
260 * Delaunay conforming.
262 * Returns: the number of remaining encroached edges. If @steiner_max
263 * is set to a negative value and gts_vertex_encroaches_edge() is used
264 * for @encroaches this should always be zero.
266 guint
gts_delaunay_conform (GtsSurface
* surface
,
268 GtsEncroachFunc encroaches
,
271 GtsFifo
* encroached
;
273 guint encroached_number
;
275 g_return_val_if_fail (surface
!= NULL
, 0);
276 g_return_val_if_fail (surface
!= NULL
, 0);
277 g_return_val_if_fail (encroaches
!= NULL
, 0);
279 datas
[0] = encroached
= gts_fifo_new ();
281 datas
[2] = encroaches
;
283 gts_surface_foreach_edge (surface
, (GtsFunc
) make_encroached_fifo
, datas
);
285 split_encroached (surface
,
289 gts_fifo_foreach (encroached
, (GtsFunc
) gts_object_reset_reserved
, NULL
);
290 encroached_number
= gts_fifo_size (encroached
);
291 gts_fifo_destroy (encroached
);
292 return encroached_number
;
295 #define EHEAP_PAIR(f) (GTS_OBJECT (f)->reserved)
297 static void heap_surface_add_face (GtsSurface
* s
, GtsFace
* f
)
299 GtsEHeap
* heap
= GTS_OBJECT (s
)->reserved
;
300 gdouble key
= gts_eheap_key (heap
, f
);
303 EHEAP_PAIR (f
) = gts_eheap_insert_with_key (heap
, f
, key
);
305 if (GTS_SURFACE_CLASS (GTS_OBJECT (s
)->klass
->parent_class
)->add_face
)
306 (* GTS_SURFACE_CLASS (GTS_OBJECT (s
)->klass
->parent_class
)->add_face
)
310 static void heap_surface_remove_face (GtsSurface
* s
, GtsFace
* f
)
312 GtsEHeap
* heap
= GTS_OBJECT (s
)->reserved
;
315 gts_eheap_remove (heap
, EHEAP_PAIR (f
));
317 if (GTS_SURFACE_CLASS (GTS_OBJECT (s
)->klass
->parent_class
)->remove_face
)
318 (* GTS_SURFACE_CLASS (GTS_OBJECT (s
)->klass
->parent_class
)->remove_face
)
322 static void heap_surface_class_init (GtsSurfaceClass
* klass
)
324 klass
->add_face
= heap_surface_add_face
;
325 klass
->remove_face
= heap_surface_remove_face
;
328 static GtsObjectClass
* heap_surface_class_new (GtsObjectClass
* parent_class
)
330 GtsObjectClassInfo heap_surface_info
;
332 heap_surface_info
= parent_class
->info
;
333 heap_surface_info
.class_init_func
= (GtsObjectClassInitFunc
)
334 heap_surface_class_init
;
335 return gts_object_class_new (parent_class
,
339 static void make_face_heap (GtsFace
* f
, GtsEHeap
* heap
)
341 gdouble key
= gts_eheap_key (heap
, f
);
344 EHEAP_PAIR (f
) = gts_eheap_insert_with_key (heap
, f
, key
);
348 * gts_delaunay_refine:
349 * @surface: a #GtsSurface describing a conforming Delaunay triangulation.
350 * @steiner_max: maximum number of Steiner points.
351 * @encroaches: a #GtsEncroachFunc.
352 * @encroach_data: user-data to pass to @encroaches.
353 * @cost: a #GtsKeyFunc used to sort the faces during refinement.
354 * @cost_data: user-data to pass to @cost.
356 * An implementation of the refinement algorithm described in Ruppert
357 * (1995) and Shewchuk (1996).
359 * Returns: the number of unrefined faces of @surface left. Should be zero
360 * if @steiner_max is set to a negative value.
362 guint
gts_delaunay_refine (GtsSurface
* surface
,
364 GtsEncroachFunc encroaches
,
365 gpointer encroach_data
,
369 GtsObjectClass
* heap_surface_class
;
370 GtsObjectClass
* original_class
;
372 GtsFifo
* encroached
;
374 guint unrefined_number
;
376 g_return_val_if_fail (surface
!= NULL
, 0);
377 g_return_val_if_fail (encroaches
!= NULL
, 0);
378 g_return_val_if_fail (cost
!= NULL
, 0);
380 original_class
= GTS_OBJECT (surface
)->klass
;
381 heap_surface_class
= heap_surface_class_new (original_class
);
382 GTS_OBJECT (surface
)->klass
= heap_surface_class
;
384 heap
= gts_eheap_new (cost
, cost_data
);
385 gts_surface_foreach_face (surface
, (GtsFunc
) make_face_heap
, heap
);
386 encroached
= gts_fifo_new ();
388 GTS_OBJECT (surface
)->reserved
= heap
;
390 while (steiner_max
-- != 0 && (f
= gts_eheap_remove_top (heap
, NULL
))) {
391 GtsVertex
*add_vertex_returned
;
393 GTS_VERTEX (gts_triangle_circumcircle_center (GTS_TRIANGLE (f
),
394 GTS_POINT_CLASS (surface
->vertex_class
)));
395 EHEAP_PAIR (f
) = NULL
;
396 g_assert (c
!= NULL
);
397 add_vertex_returned
= gts_delaunay_add_vertex (surface
, c
, f
);
398 g_assert (add_vertex_returned
== NULL
);
400 vertex_encroaches (c
, surface
, encroached
, encroaches
, encroach_data
);
401 if (!gts_fifo_is_empty (encroached
)) {
402 gts_delaunay_remove_vertex (surface
, c
);
403 steiner_max
= split_encroached (surface
,
411 unrefined_number
= gts_eheap_size (heap
);
412 gts_eheap_foreach (heap
, (GFunc
) gts_object_reset_reserved
, NULL
);
413 gts_eheap_destroy (heap
);
415 gts_fifo_foreach (encroached
, (GtsFunc
) gts_object_reset_reserved
, NULL
);
416 gts_fifo_destroy (encroached
);
418 GTS_OBJECT (surface
)->klass
= original_class
;
419 GTS_OBJECT (surface
)->reserved
= NULL
;
420 g_free (heap_surface_class
);
422 return unrefined_number
;