cae-tools/poly2d/f2d_tri_holes.cpp

/*
 * f2d_tri_holes.cpp - Triangulate a polygon with holes
 *
 * Written 2013 by Werner Almesberger
 * Copyright 2013 by Werner Almesberger
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 */

/*
 * References:
 * http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Triangulation_2/
 *   Chapter_main.html#Subsection_37.8.2
 * http://www.cgal.org/Manual/latest/examples/Triangulation_2/
 *   polygon_triangulation.cpp
 */

extern "C" {
	#include <assert.h>

	#include "util.h"
	#include "poly2d.h"
}

#if 0
/*
 * @@@ Prevent spurious aborts with
 *
 * terminate called after throwing an instance of 'CGAL::Precondition_exception'
 * what():  CGAL ERROR: precondition violation!
 * Expr: is_simple_2(first, last, traits)
 * File: /usr/include/CGAL/Polygon_2/Polygon_2_algorithms_impl.h
 * Line: 420
 *
 * Note that we also need to check the polygons for simplicity in recurse_area,
 * or this may still lead to assertion failures.
 */

#define CGAL_POLYGON_NO_PRECONDITIONS
#endif

#include "cgal_helper.h"

//#include <vector>
//#include <boost/shared_ptr.hpp>

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
#include <CGAL/Polygon_2.h>


struct FaceInfo2 {
	FaceInfo2(){}
	int nesting_level;

	bool in_domain(void)
	{ 
		return nesting_level % 2 == 1;
	}
};


typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Triangulation_vertex_base_2<K> Vb;
typedef CGAL::Triangulation_face_base_with_info_2<FaceInfo2, K> Fbb;
typedef CGAL::Constrained_triangulation_face_base_2<K, Fbb> Fb;
typedef CGAL::Triangulation_data_structure_2<Vb, Fb> TDS;
typedef CGAL::Exact_predicates_tag Itag;
typedef CGAL::Constrained_Delaunay_triangulation_2<K, TDS, Itag> CDT;
typedef CDT::Point Point;
typedef CGAL::Polygon_2<K> Polygon_2;


#if 0
typedef CGAL::Polygon_with_holes_2<K> Polygon_with_holes;

typedef boost::shared_ptr<Polygon_with_holes> PolygonPtr;

typedef std::vector<PolygonPtr> PolygonPtrVector;

	struct p2d *res = NULL, **last = &res, *np;


static void append_poly(Polygon_2 poly, struct p2d ***last)
{
	**last = P2_to_p2d(poly);
	*last = &(**last)->next;
}
#endif


/* ----- Mark domains ------------------------------------------------------ */


static void mark_domains(CDT &ct, CDT::Face_handle start, int index, 
    std::list<CDT::Edge> &border)
{
	std::list<CDT::Face_handle> queue;

	if (start->info().nesting_level != -1)
		return;

	queue.push_back(start);
	while (!queue.empty()) {
		CDT::Face_handle fh = queue.front();

		queue.pop_front();
		if (fh->info().nesting_level == -1){
			fh->info().nesting_level = index;
			for (int i = 0; i < 3; i++) {
				CDT::Edge e(fh, i);
				CDT::Face_handle n = fh->neighbor(i);
				if (n->info().nesting_level == -1) {
					if (ct.is_constrained(e))
						border.push_back(e);
					else
						queue.push_back(n);
				}
			}
		}
	}
}


static void mark_domains(CDT &cdt)
{
	int index = 0;

	for (CDT::All_faces_iterator it = cdt.all_faces_begin();
	    it != cdt.all_faces_end(); ++it)
		it->info().nesting_level = -1;

	std::list<CDT::Edge> border;
	mark_domains(cdt, cdt.infinite_face(), index++, border);
	while(!border.empty()) {
		CDT::Edge e = border.front();

		border.pop_front();
		CDT::Face_handle n = e.first->neighbor(e.second);
		if (n->info().nesting_level == -1)
			mark_domains(cdt, n, e.first->info().nesting_level+1,
			    border);
	}
}


/* ----- Inser polygon ----------------------------------------------------- */


void insert_polygon(CDT &cdt, const Polygon_2 &polygon)
{
	if (polygon.is_empty())
		return;
	CDT::Vertex_handle v_prev =
	    cdt.insert(*CGAL::cpp0x::prev(polygon.vertices_end()));
	for (Polygon_2::Vertex_iterator vit = polygon.vertices_begin();
	    vit != polygon.vertices_end(); ++vit) {
		CDT::Vertex_handle vh = cdt.insert(*vit);

		cdt.insert_constraint(vh, v_prev);
		v_prev = vh;
	}  
}


static const struct v2d *find_point(const struct p2d *p, double x, double y)
{
	const struct v2d *v = p->v;

	while (v) {
		if (v->x == x && v->y == y)
			break;
		v = v->next;
		if (v == p->v)
			return NULL;
	}
	return v;
}


static struct f2d *make_face(CDT::Finite_faces_iterator fit,
    const struct p2d *p, const struct p2d *holes)
{
	struct f2d *f;
	const struct v2d *v[3];
	const struct p2d *h;
	int i, j;

	for (i = 0; i != 3; i++) {
		Point point = fit->vertex(i)->point();
		const struct v2d *m;

		m = find_point(p, point.x(), point.y());
		if (m) {
			v[i] = m;
			continue;
		}
		for (h = holes; h; h = h->next) {
			m = find_point(h, point.x(), point.y());
			if (!m)
				continue;
			v[i] = m;
			break;
		}
		assert(m);
	}

	f = alloc_type(struct f2d);
	for (i = 0; i != 3; i++) {
		f->x[i] = v[i]->x;
		f->y[i] = v[i]->y;
		j = i+1;
		if (j == 3)
			j = 0;
		f->side[i] = v[i]->next == v[j] || v[j]->next == v[i];
	}
	f->next = NULL;
	
	return f;
}


void f2d_tri_holes_append(const struct p2d *p, const struct p2d *holes,
    struct f2d ***last)
{
	CDT cdt;
	const struct p2d *h;
	struct f2d *f;

	insert_polygon(cdt, p2d_to_P2(p, 0));
	for (h = holes; h; h = h->next) {
		assert(p2d_is_closed(h));
		insert_polygon(cdt, p2d_to_P2(h, 0));
	}
	mark_domains(cdt);

	for (CDT::Finite_faces_iterator fit = cdt.finite_faces_begin();
	    fit != cdt.finite_faces_end(); ++fit)
		if (fit->info().in_domain()) {
			f = make_face(fit, p, holes);
			**last = f;
			*last = &f->next;
		}
}


extern "C" struct f2d *f2d_tri_holes(const struct p2d *p,
    const struct p2d *holes)
{
	struct f2d *res = NULL, **last = &res;

	f2d_tri_holes_append(p, holes, &last);
	return res;
}