/* GTS - Library for the manipulation of triangulated surfaces
* Copyright (C) 1999--2002 Stéphane Popinet
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Library General Public
* License as published by the Free Software Foundation; either
* version 2 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Library General Public License for more details.
*
* You should have received a copy of the GNU Library General Public
* License along with this library; if not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <math.h>
#include "gts.h"
/*#define DEBUG*/
/*#define DEBUG_BOOLEAN*/
/*#define CHECK_ORIENTED*/
#ifdef DEBUG
# include "gts-private.h"
#endif /* DEBUG */
static void surface_inter_destroy (GtsObject * object)
{
GtsSurfaceInter * si = GTS_SURFACE_INTER (object);
gts_object_destroy (GTS_OBJECT (si->s1));
gts_object_destroy (GTS_OBJECT (si->s2));
g_slist_free (si->edges);
(* GTS_OBJECT_CLASS (gts_surface_inter_class ())->parent_class->destroy)
(object);
}
static void surface_inter_class_init (GtsObjectClass * klass)
{
klass->destroy = surface_inter_destroy;
}
static void surface_inter_init (GtsSurfaceInter * si)
{
si->s1 = si->s2 = NULL;
si->edges = NULL;
}
/**
* gts_surface_inter_class:
*
* Returns: the #GtsSurfaceInterClass.
*/
GtsSurfaceInterClass * gts_surface_inter_class (void)
{
static GtsSurfaceInterClass * klass = NULL;
if (klass == NULL) {
GtsObjectClassInfo surface_inter_info = {
"GtsSurfaceInter",
sizeof (GtsSurfaceInter),
sizeof (GtsSurfaceInterClass),
(GtsObjectClassInitFunc) surface_inter_class_init,
(GtsObjectInitFunc) surface_inter_init,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
klass = gts_object_class_new (gts_object_class (), &surface_inter_info);
}
return klass;
}
/* EdgeInter: Header */
typedef struct _EdgeInter EdgeInter;
struct _EdgeInter {
GtsEdge parent;
GtsTriangle * t1, * t2;
};
#define EDGE_INTER(obj) GTS_OBJECT_CAST (obj,\
EdgeInter,\
edge_inter_class ())
#define IS_EDGE_INTER(obj) (gts_object_is_from_class (obj,\
edge_inter_class ()))
static GtsEdgeClass * edge_inter_class (void);
static EdgeInter * edge_inter_new (GtsVertex * v1, GtsVertex * v2,
GtsTriangle * t1, GtsTriangle * t2);
/* EdgeInter: Object */
static GtsEdgeClass * edge_inter_class (void)
{
static GtsEdgeClass * klass = NULL;
if (klass == NULL) {
GtsObjectClassInfo edge_inter_info = {
"EdgeInter",
sizeof (EdgeInter),
sizeof (GtsEdgeClass),
(GtsObjectClassInitFunc) NULL,
(GtsObjectInitFunc) NULL,
(GtsArgSetFunc) NULL,
(GtsArgGetFunc) NULL
};
klass = gts_object_class_new (GTS_OBJECT_CLASS (gts_constraint_class ()),
&edge_inter_info);
}
return klass;
}
static EdgeInter * edge_inter_new (GtsVertex * v1, GtsVertex * v2,
GtsTriangle * t1, GtsTriangle * t2)
{
EdgeInter * object;
object = EDGE_INTER (gts_edge_new (GTS_EDGE_CLASS (edge_inter_class ()),
v1, v2));
object->t1 = t1;
object->t2 = t2;
return object;
}
#ifdef DEBUG
static void write_surface_graph (GtsSurface * s, FILE * fp)
{
GSList * l = NULL;
GtsGraph * g;
static void add_to_list (gpointer data, GSList ** l) {
*l = g_slist_prepend (*l, data);
}
gts_surface_foreach_vertex (s, (GtsFunc) gts_object_reset_reserved, NULL);
gts_surface_foreach_edge (s, (GtsFunc) gts_object_reset_reserved, NULL);
gts_surface_foreach_edge (s, (GtsFunc) add_to_list, &l);
g = gts_segments_graph_new (gts_graph_class (), l);
gts_graph_write_dot (g, fp);
gts_object_destroy (GTS_OBJECT (g));
g_slist_free (l);
}
#endif /* DEBUG */
static GtsPoint * segment_triangle_intersection (GtsSegment * s,
GtsTriangle * t,
GtsPointClass * klass)
{
GtsPoint * A, * B, * C, * D, * E;
gint ABCE, ABCD, ADCE, ABDE, BCDE;
GtsEdge * AB, * BC, * CA;
gdouble a, b, c;
g_return_val_if_fail (s != NULL, NULL);
g_return_val_if_fail (t != NULL, NULL);
g_return_val_if_fail (klass != NULL, NULL);
gts_triangle_vertices_edges (t, NULL,
(GtsVertex **) &A,
(GtsVertex **) &B,
(GtsVertex **) &C,
&AB, &BC, &CA);
D = GTS_POINT (s->v1);
E = GTS_POINT (s->v2);
ABCE = gts_point_orientation_3d_sos (A, B, C, E);
ABCD = gts_point_orientation_3d_sos (A, B, C, D);
if (ABCE < 0 || ABCD > 0) {
GtsPoint * tmpp;
gint tmp;
tmpp = E; E = D; D = tmpp;
tmp = ABCE; ABCE = ABCD; ABCD = tmp;
}
if (ABCE < 0 || ABCD > 0)
return NULL;
ADCE = gts_point_orientation_3d_sos (A, D, C, E);
if (ADCE < 0)
return NULL;
ABDE = gts_point_orientation_3d_sos (A, B, D, E);
if (ABDE < 0)
return NULL;
BCDE = gts_point_orientation_3d_sos (B, C, D, E);
if (BCDE < 0)
return NULL;
a = gts_point_orientation_3d (A, B, C, E);
b = gts_point_orientation_3d (A, B, C, D);
if (a != b) {
c = a/(a - b);
return gts_point_new (klass,
E->x + c*(D->x - E->x),
E->y + c*(D->y - E->y),
E->z + c*(D->z - E->z));
}
/* D and E are contained within ABC */
#ifdef DEBUG
fprintf (stderr,
"segment: %p:%s triangle: %p:%s intersection\n"
"D and E contained in ABC\n",
s, GTS_NEDGE (s)->name, t, GTS_NFACE (t)->name);
#endif /* DEBUG */
g_assert (a == 0.);
return gts_point_new (klass,
(E->x + D->x)/2.,
(E->y + D->y)/2.,
(E->z + D->z)/2.);
}
static gint triangle_triangle_orientation (GtsPoint * p1,
GtsPoint * p2, GtsPoint * p3,
GtsPoint * p4, GtsPoint * p5,
GtsPoint * p6)
{
gint o4 = 0, o5 = 0, o6 = 0;
if (p4 != p1 && p4 != p2 && p4 != p3)
o4 = gts_point_orientation_3d_sos (p1, p2, p3, p4);
if (p5 != p1 && p5 != p2 && p5 != p3)
o5 = gts_point_orientation_3d_sos (p1, p2, p3, p5);
if (o4*o5 < 0)
return 0;
if (p6 != p1 && p6 != p2 && p6 != p3)
o6 = gts_point_orientation_3d_sos (p1, p2, p3, p6);
if (o4*o6 < 0 || o5*o6 < 0)
return 0;
if (o4) return o4;
if (o5) return o5;
g_assert (o6);
return o6;
}
static gint triangle_point_orientation (GtsTriangle * t1,
GtsTriangle * t2,
gint o1,
GtsPoint * p)
{
GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (t1->e1)->v1);
GtsPoint * p2 = GTS_POINT (GTS_SEGMENT (t1->e1)->v2);
GtsPoint * p3 = GTS_POINT (gts_triangle_vertex (t1));
GtsPoint * p4 = GTS_POINT (GTS_SEGMENT (t2->e1)->v1);
GtsPoint * p5 = GTS_POINT (GTS_SEGMENT (t2->e1)->v2);
GtsPoint * p6 = GTS_POINT (gts_triangle_vertex (t2));
gint o = triangle_triangle_orientation (p1, p2, p3, p4, p5, p6);
if (o != 0)
return o;
o = triangle_triangle_orientation (p4, p5, p6, p1, p2, p3);
if (o != 0) {
gint o2 = gts_point_orientation_3d_sos (p4, p5, p6, p);
return - o*o1*o2;
}
return 0;
}
static void add_edge_inter (GtsEdge * e,
GtsTriangle * t,
GtsVertex * v)
{
GtsVertex * ev1 = GTS_SEGMENT (e)->v1, * ev2 = GTS_SEGMENT (e)->v2;
GList * i = GTS_OBJECT (e)->reserved;
GTS_OBJECT (v)->reserved = t;
if (i == NULL) {
GTS_OBJECT (e)->reserved = g_list_prepend (NULL, v);
#ifdef DEBUG
fprintf (stderr, "add_edge_inter: inserting %p (%p,%p)\n", v, e, t);
#endif /* DEBUG */
}
else {
GtsPoint * p1 = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
GtsPoint * p2 = GTS_POINT (GTS_SEGMENT (t->e1)->v2);
GtsPoint * p3 = GTS_POINT (gts_triangle_vertex (t));
gint o1, oref = gts_point_orientation_3d_sos (p1, p2, p3, GTS_POINT (ev1));
o1 = oref;
while (i) {
gint o2 = triangle_point_orientation (t, GTS_OBJECT (i->data)->reserved,
oref, GTS_POINT (ev1));
if (o2 == 0) {
#ifdef DEBUG
g_warning ("add_edge_inter: safe sign evaluation failed\n");
#endif /* DEBUG */
o2 = gts_point_orientation_3d_sos (p1, p2, p3, i->data);
}
if (o1*o2 < 0)
break;
ev1 = i->data;
o1 = o2;
i = i->next;
}
if (i != NULL) {
GList * n = g_list_prepend (NULL, v);
ev2 = i->data;
n->next = i;
n->prev = i->prev;
i->prev = n;
if (n->prev == NULL)
GTS_OBJECT (e)->reserved = n;
else
n->prev->next = n;
}
else {
g_assert (o1*gts_point_orientation_3d_sos (p1, p2, p3, GTS_POINT (ev2))
< 0);
GTS_OBJECT (e)->reserved = g_list_append (GTS_OBJECT (e)->reserved, v);
}
#ifdef DEBUG
fprintf (stderr,
"add_edge_inter: inserting %p (%p,%p) between %p and %p\n",
v, e, t, ev1, ev2);
i = GTS_OBJECT (e)->reserved;
while (i) {
fprintf (stderr, " %p", i->data);
i = i->next;
}
fprintf (stderr, "\n");
#endif /* DEBUG */
}
}
static GtsVertex * intersects (GtsEdge * e,
GtsTriangle * t,
GtsSurface * s)
{
GList * i = GTS_OBJECT (e)->reserved;
GtsVertex * v;
while (i) {
if (GTS_OBJECT (i->data)->reserved == t)
return i->data;
i = i->next;
}
v = GTS_VERTEX (segment_triangle_intersection (GTS_SEGMENT (e), t,
GTS_POINT_CLASS (s->vertex_class)));
if (v != NULL) {
#ifdef DEBUG
if (GTS_IS_NVERTEX (v) && GTS_IS_NEDGE (e) && GTS_IS_NFACE (t) &&
GTS_NVERTEX (v)->name[0] == '\0')
g_snprintf (GTS_NVERTEX (v)->name, GTS_NAME_LENGTH, "%s|%s",
GTS_NEDGE (e)->name, GTS_NFACE (t)->name);
#endif /* DEBUG */
if (s->vertex_class->intersection_attributes)
(*s->vertex_class->intersection_attributes)
(v, GTS_OBJECT (e), GTS_OBJECT (t));
add_edge_inter (e, t, v);
}
return v;
}
/* see figure misc/orientation.fig */
static gint intersection_orientation (GtsTriangle * t1,
GtsEdge * e,
GtsTriangle * t2)
{
GtsVertex * v1, * v2, * v3;
GtsEdge * e2, * e3;
GtsVertex * v4, * v5, * v6;
gts_triangle_vertices_edges (t1, e, &v1, &v2, &v3, &e, &e2, &e3);
gts_triangle_vertices (t2, &v4, &v5, &v6);
return gts_point_orientation_3d_sos (GTS_POINT (v4),
GTS_POINT (v5),
GTS_POINT (v6),
GTS_POINT (v2));
}
#define UPDATE_ORIENTATION if (o > 0) { vi2 = v; e2 = e; } else { vi2 = vi1;\
e2 = e1;\
vi1 = v;\
e1 = e; }
static void intersect_edges (GtsBBox * bb1, GtsBBox * bb2,
GtsSurfaceInter * si)
{
GtsSurface * s1 = GTS_OBJECT (si->s1)->reserved;
GtsTriangle * t1 = GTS_TRIANGLE (bb1->bounded);
GtsTriangle * t2 = GTS_TRIANGLE (bb2->bounded);
GtsVertex * v, * vi1 = NULL, * vi2 = NULL;
GtsEdge * e1 = NULL, * e2 = NULL, * e;
vi1 = intersects (t2->e1, t1, s1);
e1 = t2->e1;
v = intersects (t2->e2, t1, s1);
e = t2->e2;
if (!vi1) {
vi1 = v;
e1 = e;
}
else if (v) {
gint o = intersection_orientation (t2, t2->e2, t1);
UPDATE_ORIENTATION;
}
if (!vi2) {
v = intersects (t2->e3, t1, s1);
e = t2->e3;
if (!vi1) {
vi1 = v;
e1 = e;
}
else if (v) {
gint o = intersection_orientation (t2, t2->e3, t1);
UPDATE_ORIENTATION;
}
}
if (!vi2) {
v = intersects (t1->e1, t2, s1);
e = t1->e1;
if (!vi1) {
vi1 = v;
e1 = e;
}
else if (v) {
gint o = - intersection_orientation (t1, t1->e1, t2);
UPDATE_ORIENTATION;
}
}
if (!vi2) {
v = intersects (t1->e2, t2, s1);
e = t1->e2;
if (!vi1) {
vi1 = v;
e1 = e;
}
else if (v) {
gint o = - intersection_orientation (t1, t1->e2, t2);
UPDATE_ORIENTATION;
}
}
if (!vi2) {
v = intersects (t1->e3, t2, s1);
e = t1->e3;
if (!vi1) {
vi1 = v;
e1 = e;
}
else if (v) {
gint o = - intersection_orientation (t1, t1->e3, t2);
UPDATE_ORIENTATION;
}
}
g_assert ((!vi1 && !vi2) || (vi1 && vi2));
if (vi1) {
GtsEdge * e = GTS_EDGE (edge_inter_new (vi1, vi2, t1, t2));
#ifdef DEBUG
fprintf (stderr, "creating constraint %p: %p->%p: %p/%p\n",
e, vi1, vi2, t1, t2);
#endif /* DEBUG */
gts_surface_add_face (si->s1, GTS_FACE (t1));
gts_surface_add_face (si->s2, GTS_FACE (t2));
si->edges = g_slist_prepend (si->edges, e);
GTS_OBJECT (t1)->reserved = g_slist_prepend (GTS_OBJECT (t1)->reserved, e);
GTS_OBJECT (t2)->reserved = g_slist_prepend (GTS_OBJECT (t2)->reserved, e);
}
}
static GtsSurfaceInter * surface_inter_new (GtsSurfaceInterClass * klass,
GtsSurface * s1,
GtsSurface * s2,
GNode * faces_tree1,
GNode * faces_tree2)
{
GtsSurfaceInter * si;
si = GTS_SURFACE_INTER (gts_object_new (GTS_OBJECT_CLASS (klass)));
si->s1 = gts_surface_new (gts_surface_class (),
s1->face_class,
s1->edge_class,
s1->vertex_class);
GTS_OBJECT (si->s1)->reserved = s1;
si->s2 = gts_surface_new (gts_surface_class (),
s2->face_class,
s2->edge_class,
s2->vertex_class);
GTS_OBJECT (si->s2)->reserved = s2;
gts_bb_tree_traverse_overlapping (faces_tree1, faces_tree2,
(GtsBBTreeTraverseFunc) intersect_edges,
si);
return si;
}
static void free_slist (GtsObject * o)
{
g_slist_free (o->reserved);
o->reserved = NULL;
}
static void free_glist (GtsObject * o)
{
g_list_foreach (o->reserved, (GFunc) gts_object_reset_reserved, NULL);
g_list_free (o->reserved);
o->reserved = NULL;
}
/**
* gts_surface_intersection:
* @s1: a #GtsSurface.
* @s2: a #GtsSurface.
* @faces_tree1: a bounding box tree (see gts_bb_tree_new()) for
* the faces of @s1.
* @faces_tree2: a bounding box tree for the faces of @s2.
*
* Returns: a list of #GtsEdge defining the curve intersection of the
* two surfaces.
*/
GSList * gts_surface_intersection (GtsSurface * s1,
GtsSurface * s2,
GNode * faces_tree1,
GNode * faces_tree2)
{
GtsSurfaceInter * si;
GSList * inter;
g_return_val_if_fail (s1 != NULL, NULL);
g_return_val_if_fail (s2 != NULL, NULL);
g_return_val_if_fail (faces_tree1 != NULL, NULL);
g_return_val_if_fail (faces_tree2 != NULL, NULL);
si = surface_inter_new (gts_surface_inter_class (),
s1, s2, faces_tree1, faces_tree2);
gts_surface_foreach_face (si->s1, (GtsFunc) free_slist, NULL);
gts_surface_foreach_face (si->s2, (GtsFunc) free_slist, NULL);
gts_surface_foreach_edge (si->s1, (GtsFunc) free_glist, NULL);
gts_surface_foreach_edge (si->s2, (GtsFunc) free_glist, NULL);
inter = si->edges;
si->edges = NULL;
gts_object_destroy (GTS_OBJECT (si));
return inter;
}
typedef enum {
INTERIOR = 1 << (GTS_USER_FLAG),
RELEVANT = 1 << (GTS_USER_FLAG + 1)
} CurveFlag;
#define IS_SET(s, f) ((GTS_OBJECT_FLAGS (s) & (f)) != 0)
#define SET(s, f) (GTS_OBJECT_FLAGS (s) |= (f))
#define UNSET(s, f) (GTS_OBJECT_FLAGS (s) &= ~(f))
#define NEXT(s) (GTS_OBJECT (s)->reserved)
#ifdef DEBUG
static void print_segment (GtsSegment * s)
{
fprintf (stderr, "%p: %s->%s ", s,
GTS_NVERTEX (s->v1)->name,
GTS_NVERTEX (s->v2)->name);
if (NEXT (s)) {
GtsSegment * next = NEXT (s);
fprintf (stderr, "next %p: %s->%s\n", next,
GTS_NVERTEX (next->v1)->name,
GTS_NVERTEX (next->v2)->name);
}
else
fprintf (stderr, "next NULL\n");
}
static void write_nodes (GSList * i, GHashTable * hash, guint * nn,
FILE * fp)
{
while (i) {
GtsSegment * s = i->data;
if (!g_hash_table_lookup (hash, s->v1)) {
fprintf (fp, " %u [ label = \"%p\" ];\n", *nn, s->v1);
g_hash_table_insert (hash, s->v1, GUINT_TO_POINTER ((*nn)++));
}
if (!g_hash_table_lookup (hash, s->v2)) {
fprintf (fp, " %u [ label = \"%p\" ];\n", *nn, s->v2);
g_hash_table_insert (hash, s->v2, GUINT_TO_POINTER ((*nn)++));
}
i = i->next;
}
}
static void write_edges (GSList * i, GHashTable * hash,
GtsSurface * surface,
FILE * fp)
{
while (i) {
GtsSegment * s = i->data;
fprintf (fp, " %u -> %u [ label = \"%p:%d\" ];\n",
GPOINTER_TO_UINT (g_hash_table_lookup (hash, s->v1)),
GPOINTER_TO_UINT (g_hash_table_lookup (hash, s->v2)),
s,
gts_edge_face_number (GTS_EDGE (s), surface));
i = i->next;
}
}
static void write_graph (GSList * boundary, GSList * interior,
GtsSurface * surface,
FILE * fp)
{
GHashTable * hash = g_hash_table_new (NULL, NULL);
guint nn = 1;
fprintf (fp, "digraph oriented_curve {\n");
write_nodes (boundary, hash, &nn, fp);
write_nodes (interior, hash, &nn, fp);
write_edges (boundary, hash, surface, fp);
fprintf (fp, " edge [ color = red ];\n");
write_edges (interior, hash, surface, fp);
fprintf (fp, "}\n");
g_hash_table_destroy (hash);
}
static void write_graph1 (GtsSegment * start, GSList * i,
GtsSurface * surface,
FILE * fp)
{
GSList * boundary = NULL, * interior = NULL;
GtsSegment * s = start;
do {
boundary = g_slist_prepend (boundary, s);
s = NEXT (s);
} while (s != start);
while (i) {
if (IS_SET (i->data, INTERIOR))
interior = g_slist_prepend (interior, i->data);
i = i->next;
}
write_graph (boundary, interior, surface, fp);
g_slist_free (boundary);
g_slist_free (interior);
}
static void print_loop (GtsSegment * start, FILE * fp)
{
GtsSegment * s = start;
do {
fprintf (fp, " %p: %p:%s -> %p:%s\n",
s,
s->v1, GTS_NVERTEX (s->v1)->name,
s->v2, GTS_NVERTEX (s->v2)->name);
s = NEXT (s);
} while (s != start && s != NULL);
}
static void draw_vector (GtsPoint * p1, GtsPoint * p2, FILE * fp)
{
gdouble x = p2->x - p1->x;
gdouble y = p2->y - p1->y;
gdouble z = p2->z - p1->z;
fprintf (fp, "VECT 1 3 0 3 0 %g %g %g %g %g %g %g %g %g\n",
p1->x + x - (x - y/2.)/5.,
p1->y + y - (x/2. + y)/5.,
p1->z + z - (x/2. + z)/5.,
p1->x + x,
p1->y + y,
p1->z + z,
p1->x + x - (x + y/2.)/5.,
p1->y + y + (x/2. - y)/5.,
p1->z + z + (x/2. - z)/5.);
fprintf (fp, "VECT 1 2 0 2 0 %g %g %g %g %g %g\n",
p1->x, p1->y, p1->z,
p1->x + x,
p1->y + y,
p1->z + z);
}
static void draw_vector1 (GtsPoint * p1, GtsPoint * p2, GtsPoint * o,
FILE * fp)
{
gdouble x1 = o->x + 0.9*(p1->x - o->x);
gdouble y1 = o->y + 0.9*(p1->y - o->y);
gdouble z1 = o->z + 0.9*(p1->z - o->z);
gdouble x2 = o->x + 0.9*(p2->x - o->x);
gdouble y2 = o->y + 0.9*(p2->y - o->y);
gdouble z2 = o->z + 0.9*(p2->z - o->z);
gdouble x = x2 - x1;
gdouble y = y2 - y1;
gdouble z = z2 - z1;
fprintf (fp, "VECT 1 3 0 3 0 %g %g %g %g %g %g %g %g %g\n",
x1 + x - (x - y/2.)/5.,
y1 + y - (x/2. + y)/5.,
z1 + z - (x/2. + z)/5.,
x1 + x,
y1 + y,
z1 + z,
x1 + x - (x + y/2.)/5.,
y1 + y + (x/2. - y)/5.,
z1 + z + (x/2. - z)/5.);
fprintf (fp, "VECT 1 2 0 2 0 %g %g %g %g %g %g\n",
x1, y1, z1,
x1 + x,
y1 + y,
z1 + z);
}
static void write_segments (GSList * boundary, GSList * interior,
FILE * fp)
{
GSList * i = boundary;
fprintf (fp, "LIST {\n");
while (i) {
GSList * inext = i->next;
GtsSegment * s = i->data;
GtsSegment * next = inext ? inext->data : boundary->data;
GtsVertex * v1, * v2;
if (s->v1 != next->v1 && s->v1 != next->v2) {
v1 = s->v1;
v2 = s->v2;
}
else {
v1 = s->v2;
v2 = s->v1;
}
draw_vector (GTS_POINT (v1), GTS_POINT (v2), fp);
i = inext;
}
i = interior;
while (i) {
GtsSegment * s = i->data;
draw_vector (GTS_POINT (s->v1), GTS_POINT (s->v2), fp);
i = i->next;
}
fprintf (fp, "}\n");
}
static void write_loops (GSList * i, FILE * fp)
{
guint nl = 0;
while (i) {
GtsSegment * start = i->data, * s;
GtsPoint os;
guint n = 0;
fprintf (fp, "(geometry \"loop%d\" = LIST {\n", nl++);
os.x = os.y = os.z = 0.;
s = start;
do {
GtsSegment * next = NEXT (s);
GtsPoint * p;
if (s->v1 != next->v1 && s->v1 != next->v2)
p = GTS_POINT (s->v1);
else
p = GTS_POINT (s->v2);
os.x += p->x; os.y += p->y; os.z += p->z; n++;
s = next;
} while (s != start);
os.x /= n; os.y /= n; os.z /= n;
s = start;
do {
GtsSegment * next = NEXT (s);
if (s->v1 != next->v1 && s->v1 != next->v2)
draw_vector1 (GTS_POINT (s->v1), GTS_POINT (s->v2), &os, fp);
else
draw_vector1 (GTS_POINT (s->v2), GTS_POINT (s->v1), &os, fp);
s = next;
} while (s != start);
fprintf (fp, "})\n");
i = i->next;
}
}
#define NAME(v) (GTS_IS_NVERTEX (v) ? GTS_NVERTEX (v)->name : "")
#endif /* DEBUG */
static GtsSegment * prev_flag (GtsSegment * s, CurveFlag flag)
{
GSList * i = s->v1->segments;
while (i) {
if (i->data != s && IS_SET (i->data, flag))
return i->data;
i = i->next;
}
return NULL;
}
static GtsSegment * next_flag (GtsSegment * s, CurveFlag flag)
{
GSList * i = s->v2->segments;
while (i) {
if (i->data != s && IS_SET (i->data, flag))
return i->data;
i = i->next;
}
return NULL;
}
static GtsSegment * next_interior (GtsVertex * v)
{
GSList * i = v->segments;
while (i) {
GtsSegment * s = i->data;
if (s->v1 == v && IS_SET (s, INTERIOR))
return s;
i = i->next;
}
return NULL;
}
static GtsSegment * prev_interior (GtsVertex * v)
{
GSList * i = v->segments;
while (i) {
GtsSegment * s = i->data;
if (s->v2 == v && IS_SET (s, INTERIOR))
return s;
i = i->next;
}
return NULL;
}
static GtsSegment * reverse (GtsSegment * start,
gboolean interior,
gboolean * isloop)
{
GtsSegment * s = start, * prev = NULL, * rprev = NULL;
GtsSegment * rstart = NULL, * rstart1 = NULL;
do {
GtsSegment * rs;
g_assert (IS_EDGE_INTER (s));
rs = GTS_SEGMENT (edge_inter_new (s->v2, s->v1,
EDGE_INTER (s)->t1, EDGE_INTER (s)->t2));
if (rstart == NULL)
rstart = rs;
else if (rstart1 == NULL)
rstart1 = rs;
if (interior)
SET (rs, INTERIOR);
NEXT (rs) = rprev;
rprev = rs;
prev = s;
s = NEXT (s);
} while (s != NULL && s != start);
if (s == start) {
NEXT (rstart) = rprev;
*isloop = TRUE;
}
else {
NEXT (rstart) = start;
NEXT (prev) = rprev;
*isloop = FALSE;
}
return rstart1;
}
static GSList * interior_loops (GSList * interior)
{
GSList * i = interior;
GSList * loops = NULL;
i = interior;
while (i) {
GtsSegment * s = i->data;
if (IS_SET (s, RELEVANT)) {
GtsSegment * start = s, * end;
do {
GtsSegment * next = next_flag (s, INTERIOR);
UNSET (s, RELEVANT);
end = s;
s = NEXT (s) = next;
} while (s != NULL && s != start);
if (s == start)
loops = g_slist_prepend (loops, start);
else {
GtsSegment * next, * prev;
gboolean isloop;
s = prev_flag (start, INTERIOR);
while (s) {
UNSET (s, RELEVANT);
NEXT (s) = start;
start = s;
s = prev_flag (s, INTERIOR);
}
next = next_flag (end, RELEVANT);
prev = prev_flag (start, RELEVANT);
if (prev != NULL)
SET (start->v1, INTERIOR);
if (next != NULL)
SET (end->v2, INTERIOR);
if (next == NULL && prev == NULL)
loops = g_slist_prepend (loops, start);
else
reverse (start, TRUE, &isloop);
}
}
i = i->next;
}
return loops;
}
#define ORIENTATION(p1,p2,p3,o) (gts_point_orientation_3d (p1, p2, o, p3))
#define ORIENTATION_SOS(p1,p2,p3,o) (gts_point_orientation_3d_sos (p1, p2, o, p3))
#define ORIENTED_VERTICES(s,next,w1,w2) {\
if ((s)->v1 == (next)->v1 || (s)->v1 == (next)->v2) {\
w1 = (s)->v2;\
w2 = (s)->v1;\
}\
else {\
w1 = (s)->v1;\
w2 = (s)->v2;\
}\
}
#if 0
static GtsSegment * segment_intersects (GtsPoint * p1, GtsPoint * p2,
GSList * i,
GtsPoint * o)
{
while (i) {
GtsSegment * s = i->data;
GtsPoint * p3 = GTS_POINT (s->v1);
GtsPoint * p4 = GTS_POINT (s->v2);
if (p3 != p1 && p3 != p2 && p4 != p1 && p4 != p2) {
gdouble o1 = ORIENTATION (p3, p4, p1, o);
gdouble o2 = ORIENTATION (p3, p4, p2, o);
if ((o1 < 0. && o2 > 0.) || (o1 > 0. && o2 < 0.)) {
o1 = ORIENTATION (p1, p2, p3, o);
o2 = ORIENTATION (p1, p2, p4, o);
if ((o1 <= 0. && o2 >= 0.) || (o1 >= 0. && o2 <= 0.))
return s;
}
}
i = i->next;
}
return NULL;
}
#else
static GtsSegment * segment_intersects (GtsPoint * p1, GtsPoint * p2,
GSList * i,
GtsPoint * o)
{
while (i) {
GtsSegment * s = i->data;
GtsPoint * p3 = GTS_POINT (s->v1);
GtsPoint * p4 = GTS_POINT (s->v2);
if (p3 != p1 && p3 != p2 && p4 != p1 && p4 != p2) {
gint o1 = ORIENTATION_SOS (p3, p4, p1, o);
gint o2 = ORIENTATION_SOS (p3, p4, p2, o);
if (o1*o2 < 0) {
o1 = ORIENTATION_SOS (p1, p2, p3, o);
o2 = ORIENTATION_SOS (p1, p2, p4, o);
if (o1*o2 < 0)
return s;
}
}
i = i->next;
}
return NULL;
}
#endif
static gboolean is_inside_wedge (GtsSegment * s1, GtsSegment * s2,
GtsPoint * p, GtsPoint * o)
{
GtsVertex * v1, * v2, * v3;
ORIENTED_VERTICES (s1, s2, v1, v2);
v3 = s2->v1 != v2 ? s2->v1 : s2->v2;
if (ORIENTATION (GTS_POINT (v1), GTS_POINT (v2),
GTS_POINT (v3), o) >= 0.) {
if (ORIENTATION (GTS_POINT (v1), GTS_POINT (v2), p, o) <= 0. ||
ORIENTATION (GTS_POINT (v2), GTS_POINT (v3), p, o) <= 0.)
return FALSE;
}
else if (ORIENTATION (GTS_POINT (v1), GTS_POINT (v2), p, o) <= 0. &&
ORIENTATION (GTS_POINT (v2), GTS_POINT (v3), p, o) <= 0.)
return FALSE;
return TRUE;
}
static GtsSegment * connection (GtsPoint * p,
GSList * interior,
GSList * bloops,
GtsPoint * o)
{
while (bloops) {
GtsSegment * start = bloops->data, * s = start;
do {
GtsSegment * next = NEXT (s);
GtsVertex * v2 = s->v1 == next->v1 || s->v1 == next->v2 ? s->v1 : s->v2;
if (is_inside_wedge (s, next, p, o) &&
!segment_intersects (p, GTS_POINT (v2), interior, o))
return s;
s = next;
} while (s != start);
bloops = bloops->next;
}
return NULL;
}
static gdouble loop_orientation (GtsSegment * start,
GtsPoint * p, GtsPoint * o)
{
GtsSegment * s = start;
gdouble or = 0.;
do {
GtsSegment * next = NEXT (s);
GtsVertex * v1, * v2;
ORIENTED_VERTICES (s, next, v1, v2);
or += ORIENTATION (p, GTS_POINT (v1), GTS_POINT (v2), o);
s = next;
} while (s != start);
#ifdef DEBUG
fprintf (stderr, "loop orientation: %g\n", or);
#endif /* DEBUG */
return or;
}
static void connect_interior_loop (GtsSegment * start,
GSList ** interior,
GSList ** bloops,
GtsSurface * surface,
GtsPoint * o)
{
GtsSegment * s = start, * c = NULL, * next, * s1, * rs1, * rs;
GtsVertex * v, * cv;
gboolean isloop;
do {
if (!(c = connection (GTS_POINT (s->v2), *interior, *bloops, o)))
s = NEXT (s);
} while (s != start && !c);
g_assert (c);
next = NEXT (c);
v = c->v1 == next->v1 || c->v1 == next->v2 ? c->v1 : c->v2;
cv = s->v2;
#ifdef DEBUG
fprintf (stderr, "connecting %p:%s with %p:%s\n",
cv, NAME (cv), v, NAME (v));
fprintf (stderr, " c: %p: %p:%s %p:%s\n", c,
c->v1, NAME (c->v1),
c->v2, NAME (c->v2));
fprintf (stderr, " next: %p: %p:%s %p:%s\n", next,
next->v1, NAME (next->v1),
next->v2, NAME (next->v2));
#endif /* DEBUG */
rs = reverse (s, FALSE, &isloop);
if (isloop) {
if (loop_orientation (rs, GTS_POINT (v), o) < 0.) {
GtsSegment * tmp = s;
s = rs;
rs = tmp;
}
*bloops = g_slist_prepend (*bloops, rs);
}
s1 = GTS_SEGMENT (gts_edge_new (surface->edge_class, v, cv));
rs1 = GTS_SEGMENT (gts_edge_new (surface->edge_class, cv, v));
NEXT (c) = s1;
NEXT (rs1) = next;
*interior = g_slist_prepend (*interior, s1);
NEXT (s1) = NEXT (s);
NEXT (s) = rs1;
}
static GSList * boundary_loops (GSList * boundary)
{
GSList * i = boundary;
GtsSegment * start = i->data;
GSList * loops = NULL;
while (i) {
GtsSegment * s = i->data;
GSList * inext = i->next;
GtsSegment * next = inext ? inext->data : start;
GtsVertex * v = s->v1 == next->v1 || s->v1 == next->v2 ? s->v1 : s->v2;
if (IS_SET (v, INTERIOR)) {
GtsSegment * intprev = prev_interior (v);
NEXT (intprev) = next;
NEXT (s) = next_interior (v);
UNSET (v, INTERIOR);
}
else
NEXT (s) = next;
i = inext;
}
i = boundary;
while (i) {
start = i->data;
if (IS_SET (start, RELEVANT)) {
GtsSegment * s = start;
do {
UNSET (s, RELEVANT);
UNSET (s, INTERIOR);
s = NEXT (s);
} while (s != start);
loops = g_slist_prepend (loops, start);
}
i = i->next;
}
return loops;
}
typedef struct _Ear Ear;
struct _Ear {
GtsVertex * v1, * v2, * v3;
GtsSegment * s1, * s2, * s3;
};
static gboolean point_in_wedge (GtsPoint * p1, GtsPoint * p2, GtsPoint * p3,
GtsPoint * p, gboolean closed, GtsPoint * o)
{
gdouble o1;
if (p == p2 || p == p3)
return FALSE;
o1 = ORIENTATION (p1, p2, p, o);
if ((closed && o1 < 0.) || (!closed && o1 <= 0.)) return FALSE;
o1 = ORIENTATION (p3, p1, p, o);
if ((closed && o1 < 0.) || (!closed && o1 <= 0.)) return FALSE;
return TRUE;
}
#if 0
static gboolean segment_intersects1 (GtsPoint * p1, GtsPoint * p2,
GtsPoint * p3, GtsPoint * p4,
gboolean closed, GtsPoint * o)
{
gdouble o1 = ORIENTATION (p3, p4, p1, o);
gdouble o2 = ORIENTATION (p3, p4, p2, o);
gdouble o3, o4;
if ((closed && ((o1 > 0. && o2 > 0.) || (o1 < 0. && o2 < 0.))) ||
(!closed && ((o1 >= 0. && o2 >= 0.) || (o1 <= 0. && o2 <= 0.))))
return FALSE;
o3 = ORIENTATION (p1, p2, p3, o);
o4 = ORIENTATION (p1, p2, p4, o);
if ((o3 > 0. && o4 > 0.) || (o3 < 0. && o4 < 0.))
return FALSE;
if (closed) return TRUE;
if ((o3 == 0. && o4 > 0.) || (o4 == 0. && o3 > 0.))
return TRUE;
return FALSE;
}
#else
static gboolean segment_intersects1 (GtsPoint * p1, GtsPoint * p2,
GtsPoint * p3, GtsPoint * p4,
gboolean closed, GtsPoint * o)
{
gint o1, o2;
o1 = ORIENTATION_SOS (p3, p4, p1, o);
o2 = ORIENTATION_SOS (p3, p4, p2, o);
if (o1*o2 > 0)
return FALSE;
o1 = ORIENTATION_SOS (p1, p2, p3, o);
o2 = ORIENTATION_SOS (p1, p2, p4, o);
if (o1*o2 > 0)
return FALSE;
return TRUE;
}
#endif
static GtsSegment * triangle_intersects_segments (GtsPoint * p1,
GtsPoint * p2,
GtsPoint * p3,
gboolean closed,
GtsSegment * start,
GtsPoint * o)
{
GtsSegment * s = start;
do {
GtsPoint * p4 = GTS_POINT (s->v1);
GtsPoint * p5 = GTS_POINT (s->v2);
if (p4 == p1) {
if (point_in_wedge (p1, p2, p3, p5, closed, o))
return s;
}
else if (p4 == p2) {
if (point_in_wedge (p2, p3, p1, p5, closed, o))
return s;
}
else if (p4 == p3) {
if (point_in_wedge (p3, p1, p2, p5, closed, o))
return s;
}
else if (p5 == p1) {
if (point_in_wedge (p1, p2, p3, p4, closed, o))
return s;
}
else if (p5 == p2) {
if (point_in_wedge (p2, p3, p1, p4, closed, o))
return s;
}
else if (p5 == p3) {
if (point_in_wedge (p3, p1, p2, p4, closed, o))
return s;
}
else if (segment_intersects1 (p1, p2, p4, p5, closed, o) ||
segment_intersects1 (p2, p3, p4, p5, closed, o) ||
segment_intersects1 (p3, p1, p4, p5, closed, o))
return s;
s = NEXT (s);
} while (s != start);
return NULL;
}
static gboolean new_ear (GtsSegment * s,
Ear * e,
GtsSegment * start,
guint sloppy,
GtsPoint * o)
{
gdouble or;
e->s1 = s;
e->s2 = NEXT (s);
g_return_val_if_fail (e->s2, FALSE);
g_return_val_if_fail (e->s2 != e->s1, FALSE);
ORIENTED_VERTICES (e->s1, e->s2, e->v1, e->v2);
e->v3 = e->s2->v1 != e->v2 ? e->s2->v1 : e->s2->v2;
if (e->v3 == e->v1)
return FALSE;
e->s3 = NEXT (e->s2);
if (gts_segment_connect (e->s3, e->v1, e->v3)) {
if (NEXT (e->s3) != e->s1)
return FALSE;
}
else if (gts_vertices_are_connected (e->v1, e->v3))
return FALSE;
else
e->s3 = NULL;
or = ORIENTATION (GTS_POINT (e->v1), GTS_POINT (e->v2), GTS_POINT (e->v3),o);
switch (sloppy) {
case 0:
if (or <= 0. ||
triangle_intersects_segments (GTS_POINT (e->v1), GTS_POINT (e->v2),
GTS_POINT (e->v3), TRUE, start, o))
return FALSE;
break;
case 1:
if (or < 0. ||
(or > 0. &&
triangle_intersects_segments (GTS_POINT (e->v1), GTS_POINT (e->v2),
GTS_POINT (e->v3), FALSE, start, o)))
return FALSE;
break;
case 2:
if ((or > 0. &&
triangle_intersects_segments (GTS_POINT (e->v1), GTS_POINT (e->v2),
GTS_POINT (e->v3), FALSE, start, o)) ||
(or < 0. &&
triangle_intersects_segments (GTS_POINT (e->v2), GTS_POINT (e->v1),
GTS_POINT (e->v3), FALSE, start, o)))
return FALSE;
break;
case 3:
if (or < 0.)
return FALSE;
break;
}
#ifdef DEBUG
if (or <= 0.)
fprintf (stderr, "or: %g\n", or);
#endif /* DEBUG */
g_assert (or > -1e-6);
return TRUE;
}
static void triangulate_loop (GtsSegment * start,
GtsSurface * surface,
GtsPoint * o)
{
GtsSegment * prev = start, * s;
guint sloppy = 0;
#ifdef DEBUG
guint nt = 0;
#endif /* DEBUG */
s = NEXT (start);
while (NEXT (s) != s) {
GtsSegment * next = NEXT (s);
Ear e;
#ifdef DEBUG
fprintf (stderr, "prev: %p s: %p next: %p\n", prev, s, next);
#endif /* DEBUG */
if (!new_ear (s, &e, start, sloppy, o)) {
if (s == start) {
sloppy++;
#ifdef DEBUG
fprintf (stderr, "sloppy: %u\n", sloppy);
#endif /* DEBUG */
}
prev = s;
s = next;
}
else {
GtsFace * f;
if (!GTS_IS_EDGE (e.s3))
e.s3 = GTS_SEGMENT (gts_edge_new (surface->edge_class, e.v1, e.v3));
f = gts_face_new (surface->face_class,
GTS_EDGE (e.s1), GTS_EDGE (e.s2), GTS_EDGE (e.s3));
gts_surface_add_face (surface, f);
UNSET (e.s1, RELEVANT);
UNSET (e.s1, INTERIOR);
UNSET (e.s2, RELEVANT);
UNSET (e.s2, INTERIOR);
NEXT (prev) = e.s3;
NEXT (e.s3) = NEXT (e.s2);
NEXT (e.s1) = NEXT (e.s2) = NULL;
start = prev;
s = NEXT (prev);
sloppy = 0;
#ifdef DEBUG
{
gchar name[80];
FILE * fp;
fprintf (stderr, " t.%u: (%p:%s,%p:%s,%p:%s)\n",
nt,
e.v1, NAME (e.v1),
e.v2, NAME (e.v2),
e.v3, NAME (e.v3));
sprintf (name, "/tmp/t.%u", nt++);
fp = fopen (name, "wt");
// gts_surface_write (surface, fp);
gts_write_triangle (GTS_TRIANGLE (f), NULL, fp);
// write_graph1 (start, interior, surface, fp);
fclose (fp);
print_loop (start, stderr);
}
#endif /* DEBUG */
}
}
UNSET (s, RELEVANT);
UNSET (s, INTERIOR);
NEXT (s) = NULL;
}
static void check_object (GtsObject * o)
{
g_assert (o->reserved == NULL);
g_assert (o->flags == 0);
}
static void check_boundary (GtsEdge * e, GtsSurface * s)
{
check_object (GTS_OBJECT (e));
check_object (GTS_OBJECT (GTS_SEGMENT (e)->v1));
check_object (GTS_OBJECT (GTS_SEGMENT (e)->v2));
g_assert (gts_edge_face_number (e, s) == 1);
}
static void check_interior (GtsEdge * e, GtsSurface * s)
{
guint n;
check_object (GTS_OBJECT (e));
check_object (GTS_OBJECT (GTS_SEGMENT (e)->v1));
check_object (GTS_OBJECT (GTS_SEGMENT (e)->v2));
n = gts_edge_face_number (e, s);
#ifdef DEBUG
if (n != 2)
gts_surface_print_stats (s, stderr);
#endif /* DEBUG */
g_assert (n == 2);
}
static void check_boundary_interior_triangulation (GSList * boundary,
GSList * interior,
GtsSurface * surface)
{
g_slist_foreach (boundary, (GFunc) check_boundary, surface);
g_slist_foreach (interior, (GFunc) check_interior, surface);
}
static void merge_duplicate (GtsEdge * e)
{
GtsEdge * dup = gts_edge_is_duplicate (e);
g_assert (dup);
gts_edge_replace (dup, e);
gts_object_destroy (GTS_OBJECT (dup));
}
static void triangulate_boundary_interior (GSList * boundary,
GSList * interior,
GtsSurface * s,
GtsPoint * o)
{
GSList * iloops, * bloops, * i;
i = boundary;
while (i) {
SET (i->data, RELEVANT);
i = i->next;
}
i = interior;
while (i) {
SET (i->data, RELEVANT);
SET (i->data, INTERIOR);
i = i->next;
}
iloops = interior_loops (interior);
bloops = boundary_loops (boundary);
i = iloops;
while (i) {
#ifdef DEBUG
fprintf (stderr, "--- interior loop ---\n");
print_loop (i->data, stderr);
#endif /* DEBUG */
connect_interior_loop (i->data, &interior, &bloops, s, o);
i = i->next;
}
#ifdef DEBUG
{
FILE * fp = fopen ("/tmp/bloops", "w");
write_loops (bloops, fp);
fclose (fp);
}
#endif /* DEBUG */
i = bloops;
while (i) {
#ifdef DEBUG
fprintf (stderr, "--- boundary loop ---\n");
print_loop (i->data, stderr);
#endif /* DEBUG */
triangulate_loop (i->data, s, o);
i = i->next;
}
g_slist_foreach (interior, (GFunc) merge_duplicate, NULL);
g_slist_free (iloops);
g_slist_free (bloops);
#ifdef CHECK_ORIENTED
check_boundary_interior_triangulation (boundary, interior, s);
#endif /* CHECK_ORIENTED */
}
static void create_edges (GtsSegment * s, GtsSurface * surface)
{
if (GTS_OBJECT (s)->reserved) {
GList * i = GTS_OBJECT (s)->reserved;
GtsVertex * v1 = i->data;
GTS_OBJECT (s)->reserved = g_list_prepend (i,
gts_edge_new (surface->edge_class, s->v1, v1));
while (i) {
GList * next = i->next;
GtsVertex * v2 = next ? next->data : s->v2;
GTS_OBJECT (i->data)->reserved = NULL;
i->data = gts_edge_new (surface->edge_class, v1, v2);
v1 = v2;
i = next;
}
}
}
static void add_boundary (GtsSegment * s, GtsSegment * next,
GSList ** boundary)
{
if (GTS_OBJECT (s)->reserved == NULL)
*boundary = g_slist_prepend (*boundary, s);
else {
if (s->v2 == next->v2 || s->v2 == next->v1) {
GList * i = g_list_last (GTS_OBJECT (s)->reserved);
while (i) {
*boundary = g_slist_prepend (*boundary, i->data);
i = i->prev;
}
}
else {
GList * i = GTS_OBJECT (s)->reserved;
while (i) {
*boundary = g_slist_prepend (*boundary, i->data);
i = i->next;
}
}
}
}
static void triangulate_face (GtsTriangle * t, GtsSurface * surface)
{
GSList * interior = GTS_OBJECT (t)->reserved;
GSList * boundary = NULL;
GtsSurface * s = gts_surface_new (gts_surface_class (),
surface->face_class,
surface->edge_class,
surface->vertex_class);
gdouble x, y, z;
GtsPoint * p = GTS_POINT (GTS_SEGMENT (t->e1)->v1);
GtsPoint * o;
GTS_OBJECT (t)->reserved = NULL;
gts_triangle_normal (t, &x, &y, &z);
g_assert (x != 0. || y != 0. || z != 0.);
o = gts_point_new (gts_point_class (), p->x + x, p->y + y, p->z + z);
add_boundary (GTS_SEGMENT (t->e3), GTS_SEGMENT (t->e1), &boundary);
add_boundary (GTS_SEGMENT (t->e2), GTS_SEGMENT (t->e3), &boundary);
add_boundary (GTS_SEGMENT (t->e1), GTS_SEGMENT (t->e2), &boundary);
#ifdef DEBUG
{
static guint nt = 0;
char name[80];
FILE * fp;
fprintf (stderr, "%u: triangulating %p\n", nt, t);
if (nt == 28)
fprintf (stderr, "tintin!!!!\n");
sprintf (name, "/tmp/oc.%u", nt++);
fp = fopen (name, "wt");
// write_graph (boundary, interior, s, fp);
write_segments (boundary, interior, fp);
fclose (fp);
}
#endif /* DEBUG */
triangulate_boundary_interior (boundary, interior, s, o);
g_slist_free (interior);
g_slist_free (boundary);
if (GTS_OBJECT (t)->klass->attributes)
gts_surface_foreach_face (s, (GtsFunc) gts_object_attributes, t);
gts_surface_merge (surface, s);
gts_object_destroy (GTS_OBJECT (s));
gts_object_destroy (GTS_OBJECT (o));
}
static void free_edge_list (GtsObject * o)
{
g_list_free (o->reserved);
o->reserved = NULL;
}
/**
* gts_surface_inter_new:
* @klass: a #GtsSurfaceInterClass.
* @s1: a #GtsSurface.
* @s2: a #GtsSurface.
* @faces_tree1: a bounding box tree (see gts_bb_tree_new()) for
* the faces of @s1.
* @faces_tree2: a bounding box tree for the faces of @s2.
* @is_open1: whether @s1 is an "open" surface.
* @is_open2: whether @s2 is an "open" surface.
*
* When triangulating the cut faces, the new faces inherit the
* attributes of these original faces through their attributes()
* method.
*
* Returns: a new #GtsSurfaceInter describing the intersection of @s1
* and @s2.
*/
GtsSurfaceInter * gts_surface_inter_new (GtsSurfaceInterClass * klass,
GtsSurface * s1,
GtsSurface * s2,
GNode * faces_tree1,
GNode * faces_tree2,
gboolean is_open1,
gboolean is_open2)
{
GtsSurfaceInter * si;
GtsSurface * s;
g_return_val_if_fail (klass != NULL, NULL);
g_return_val_if_fail (s1 != NULL, NULL);
g_return_val_if_fail (s2 != NULL, NULL);
g_return_val_if_fail (faces_tree1 != NULL, NULL);
g_return_val_if_fail (faces_tree2 != NULL, NULL);
si = surface_inter_new (klass, s1, s2, faces_tree1, faces_tree2);
gts_surface_foreach_edge (si->s1, (GtsFunc) create_edges, si->s1);
gts_surface_foreach_edge (si->s2, (GtsFunc) create_edges, si->s2);
#ifdef DEBUG
fprintf (stderr, "====== triangulating s1 ======\n");
#endif /* DEBUG */
s = gts_surface_new (gts_surface_class (),
s1->face_class,
s1->edge_class,
s1->vertex_class);
gts_surface_foreach_face (si->s1, (GtsFunc) triangulate_face, s);
gts_surface_foreach_edge (si->s1, (GtsFunc) free_edge_list, NULL);
gts_object_destroy (GTS_OBJECT (si->s1));
si->s1 = s;
GTS_OBJECT (si->s1)->reserved = s1;
#ifdef DEBUG
fprintf (stderr, "====== triangulating s2 ======\n");
#endif /* DEBUG */
s = gts_surface_new (gts_surface_class (),
s2->face_class,
s2->edge_class,
s2->vertex_class);
gts_surface_foreach_face (si->s2, (GtsFunc) triangulate_face, s);
gts_surface_foreach_edge (si->s2, (GtsFunc) free_edge_list, NULL);
gts_object_destroy (GTS_OBJECT (si->s2));
si->s2 = s;
GTS_OBJECT (si->s2)->reserved = s2;
return si;
}
static void check_surface_edge (GtsEdge * e, gpointer * data)
{
gboolean * ok = data[0];
GtsSurface * s = data[1];
GtsSurface * bs = GTS_OBJECT (s)->reserved;
guint nf = gts_edge_face_number (e, s);
if (nf < 1 || nf > 2) {
*ok = FALSE;
g_return_if_fail (nf >= 1 && nf <= 2);
}
if (nf == 1 && gts_edge_face_number (e, bs) == 0) {
*ok = FALSE;
g_return_if_fail (gts_edge_face_number (e, bs) > 0);
}
}
static void mark_edge (GtsObject * o, gpointer data)
{
o->reserved = data;
}
static gint triangle_orientation (GtsTriangle * t, GtsEdge * e)
{
GtsSegment * s = GTS_SEGMENT (t->e1 == e ? t->e2
:
t->e2 == e ? t->e3
:
t->e1);
GtsVertex * v2 = GTS_SEGMENT (e)->v2;
if (s->v1 == v2 || s->v2 == v2)
return 1;
return -1;
}
static gboolean check_orientation (GtsEdge * e, GtsSurface * s)
{
GtsTriangle * t1 = NULL, * t2 = NULL;
GSList * i = e->triangles;
gint o1 = 0, o2 = 0;
while (i) {
if (GTS_IS_FACE (i->data) &&
gts_face_has_parent_surface (i->data, s)) {
if (t1 == NULL) {
t1 = i->data;
o1 = triangle_orientation (t1, e);
}
else if (t2 == NULL) {
t2 = i->data;
o2 = triangle_orientation (t2, e);
g_return_val_if_fail (o1*o2 < 0, FALSE);
}
else
g_assert_not_reached ();
}
i = i->next;
}
g_return_val_if_fail (t1 && t2, FALSE);
return TRUE;
}
static void check_edge (GtsSegment * s, gpointer * data)
{
gboolean * ok = data[0];
GtsSurfaceInter * si = data[1];
gboolean * closed = data[2];
GSList * j;
guint nn = 0;
j = s->v1->segments;
while (j && *ok) {
GtsSegment * s1 = j->data;
if (s1 != s && GTS_OBJECT (s1)->reserved == si) {
if (s1->v2 != s->v1)
*ok = FALSE;
nn++;
}
j = j->next;
}
j = s->v2->segments;
while (j && *ok) {
GtsSegment * s1 = j->data;
if (s1 != s && GTS_OBJECT (s1)->reserved == si) {
if (s1->v1 != s->v2)
*ok = FALSE;
nn++;
}
j = j->next;
}
if (nn != 2)
*closed = FALSE;
if (!check_orientation (GTS_EDGE (s), si->s1))
*ok = FALSE;
if (!check_orientation (GTS_EDGE (s), si->s2))
*ok = FALSE;
}
/**
* gts_surface_inter_check:
* @si: a #GtsSurfaceInter.
* @closed: is set to %TRUE if @si->edges is a closed curve, %FALSE
* otherwise.
*
* Returns: %TRUE if the curve described by @si is an orientable
* manifold, %FALSE otherwise.
*/
gboolean gts_surface_inter_check (GtsSurfaceInter * si,
gboolean * closed)
{
gboolean ok = TRUE;
gpointer data[3];
g_return_val_if_fail (si != NULL, FALSE);
g_return_val_if_fail (closed != NULL, FALSE);
*closed = si->edges ? TRUE : FALSE;
/* mark edges as used by si */
g_slist_foreach (si->edges, (GFunc) mark_edge, si);
data[0] = &ok;
data[1] = si;
data[2] = closed;
g_slist_foreach (si->edges, (GFunc) check_edge, data);
g_slist_foreach (si->edges, (GFunc) gts_object_reset_reserved, NULL);
/* check connectivity of the faces of @si */
if (*closed) {
gpointer data[2];
data[0] = &ok;
data[1] = si->s1;
gts_surface_foreach_edge (si->s1, (GtsFunc) check_surface_edge, data);
data[1] = si->s2;
gts_surface_foreach_edge (si->s2, (GtsFunc) check_surface_edge, data);
}
return ok;
}
/* Given @e and @f returns a #GtsFace compatible with @f and belonging to
@s1 or @s2 */
static GtsFace * next_compatible_face (GtsEdge * e,
GtsFace * f,
GtsSurface * s1,
GtsSurface * s2)
{
GSList * i = e->triangles;
GtsFace * f2 = NULL, * f3 = NULL;
while (i) {
GtsFace * f1 = i->data;
if (f1 != f && GTS_IS_FACE (f1)) {
if (gts_face_has_parent_surface (f1, s1))
return f1;
if (gts_face_has_parent_surface (f1, s2)) {
if (f2 == NULL) f2 = f1;
else if (f3 == NULL) f3 = f1;
else g_assert_not_reached (); /* s2 is a non-manifold surface */
}
}
i = i->next;
}
if (f3 == NULL) {
if (gts_edge_is_boundary (e, s2))
return NULL;
return f2;
}
g_assert (gts_face_has_parent_surface (f, s1));
if (gts_triangles_are_compatible (GTS_TRIANGLE (f), GTS_TRIANGLE (f2), e))
return f2;
return f3;
}
static void walk_faces (GtsEdge * e, GtsFace * f,
GtsSurface * s1,
GtsSurface * s2,
GtsSurface * s)
{
GtsFifo * faces = gts_fifo_new ();
GtsFifo * edges = gts_fifo_new ();
gts_fifo_push (faces, f);
gts_fifo_push (edges, e);
while ((f = gts_fifo_pop (faces)) && (e = gts_fifo_pop (edges))) {
if (!GTS_OBJECT (f)->reserved) {
GtsTriangle * t = GTS_TRIANGLE (f);
GtsFace * f1;
gts_surface_add_face (s, f);
GTS_OBJECT (f)->reserved = s;
if (t->e1 != e && !GTS_OBJECT (t->e1)->reserved &&
(f1 = next_compatible_face (t->e1, f, s1, s2))) {
gts_fifo_push (faces, f1);
gts_fifo_push (edges, t->e1);
}
if (t->e2 != e && !GTS_OBJECT (t->e2)->reserved &&
(f1 = next_compatible_face (t->e2, f, s1, s2))) {
gts_fifo_push (faces, f1);
gts_fifo_push (edges, t->e2);
}
if (t->e3 != e && !GTS_OBJECT (t->e3)->reserved &&
(f1 = next_compatible_face (t->e3, f, s1, s2))) {
gts_fifo_push (faces, f1);
gts_fifo_push (edges, t->e3);
}
}
}
gts_fifo_destroy (faces);
gts_fifo_destroy (edges);
}
/**
* gts_surface_inter_boolean:
* @si: a #GtsSurfaceInter.
* @surface: a #GtsSurface.
* @op: a #GtsBooleanOperation.
*
* Adds to @surface the part of the surface described by @si and @op.
*/
void gts_surface_inter_boolean (GtsSurfaceInter * si,
GtsSurface * surface,
GtsBooleanOperation op)
{
GtsSurface * s = NULL;
gint orient = 1;
GSList * i;
g_return_if_fail (si != NULL);
g_return_if_fail (surface != NULL);
switch (op) {
case GTS_1_OUT_2: s = si->s1; orient = 1; break;
case GTS_1_IN_2: s = si->s1; orient = -1; break;
case GTS_2_OUT_1: s = si->s2; orient = -1; break;
case GTS_2_IN_1: s = si->s2; orient = 1; break;
default: g_assert_not_reached ();
}
/* mark edges as belonging to intersection */
g_slist_foreach (si->edges, (GFunc) mark_edge, si);
i = si->edges;
while (i) {
GtsEdge * e = i->data;
GSList * j = e->triangles;
while (j) {
if (gts_face_has_parent_surface (j->data, s) &&
orient*triangle_orientation (j->data, e) > 0) {
#ifdef DEBUG_BOOLEAN
GtsFace * boundary = gts_edge_is_boundary (e, surface);
g_assert (!boundary || boundary == j->data);
#endif /* DEBUG_BOOLEAN */
walk_faces (e, j->data, s, GTS_OBJECT (s)->reserved, surface);
break;
}
j = j->next;
}
i = i->next;
}
g_slist_foreach (si->edges, (GFunc) gts_object_reset_reserved, NULL);
gts_surface_foreach_face (surface,
(GtsFunc) gts_object_reset_reserved, NULL);
}
static void self_intersecting (GtsBBox * bb1, GtsBBox * bb2,
gpointer * d)
{
GtsTriangle * t1 = bb1->bounded;
GtsTriangle * t2 = bb2->bounded;
if (t1 != t2) {
GtsSegment * s1 = GTS_SEGMENT (t1->e1);
GtsSegment * s2 = GTS_SEGMENT (t1->e2);
GtsSegment * s3 = GTS_SEGMENT (t1->e3);
GtsSegment * s4 = GTS_SEGMENT (t2->e1);
GtsSegment * s5 = GTS_SEGMENT (t2->e2);
GtsSegment * s6 = GTS_SEGMENT (t2->e3);
GtsPoint * pi;
if ((!gts_segments_touch (s4, s1) &&
!gts_segments_touch (s4, s2) &&
!gts_segments_touch (s4, s3) &&
(pi = segment_triangle_intersection (s4, t1, gts_point_class ()))
!= NULL) ||
(!gts_segments_touch (s5, s1) &&
!gts_segments_touch (s5, s2) &&
!gts_segments_touch (s5, s3) &&
(pi = segment_triangle_intersection (s5, t1, gts_point_class ()))
!= NULL) ||
(!gts_segments_touch (s6, s1) &&
!gts_segments_touch (s6, s2) &&
!gts_segments_touch (s6, s3) &&
(pi = segment_triangle_intersection (s6, t1, gts_point_class ()))
!= NULL)) {
GtsBBTreeTraverseFunc func = d[0];
gpointer data = d[1];
gboolean * self_inter = d[2];
gts_object_destroy (GTS_OBJECT (pi));
*self_inter = TRUE;
(* func) (bb1, bb2, data);
}
}
}
/**
* gts_surface_foreach_intersecting_face:
* @s: a #GtsSurface.
* @func: a #GtsBBTreeTraverseFunc.
* @data: user data to pass to @func.
*
* Calls @func for each intersecting pair of faces of @s.
*
* Returns: %TRUE if @func was called at least once, %FALSE otherwise.
*/
gboolean gts_surface_foreach_intersecting_face (GtsSurface * s,
GtsBBTreeTraverseFunc func,
gpointer data)
{
GNode * tree;
gpointer d[3];
gboolean self_inter = FALSE;
g_return_val_if_fail (s != NULL, FALSE);
g_return_val_if_fail (func != NULL, FALSE);
tree = gts_bb_tree_surface (s);
d[0] = func;
d[1] = data;
d[2] = &self_inter;
gts_bb_tree_traverse_overlapping (tree, tree,
(GtsBBTreeTraverseFunc) self_intersecting,
d);
gts_bb_tree_destroy (tree, TRUE);
return self_inter;
}
static void add_intersecting (GtsBBox * bb1, GtsBBox * bb2,
GtsSurface * intersected)
{
gts_surface_add_face (intersected, bb1->bounded);
gts_surface_add_face (intersected, bb2->bounded);
}
/**
* gts_surface_is_self_intersecting:
* @s: a #GtsSurface.
*
* Returns: a new #GtsSurface containing the faces of @s which are
* self-intersecting or %NULL if no faces of @s are self-intersecting.
*/
GtsSurface * gts_surface_is_self_intersecting (GtsSurface * s)
{
GtsSurface * intersected;
g_return_val_if_fail (s != NULL, NULL);
intersected = gts_surface_new (GTS_SURFACE_CLASS (GTS_OBJECT (s)->klass),
s->face_class,
s->edge_class,
s->vertex_class);
if (!gts_surface_foreach_intersecting_face (s,
(GtsBBTreeTraverseFunc) add_intersecting, intersected)) {
gts_object_destroy (GTS_OBJECT (intersected));
intersected = NULL;
}
return intersected;
}
syntax highlighted by Code2HTML, v. 0.9.1