Crust

How to generate a crust ?

Input: a set P of points in the plane

Output: a set S of segments with endpoints in P

1. Construct the Voronoi diagram VD of P.
2. Construct the Delaunay triangulation DT of P+V, where V are the vertices of VD.
3. S:= all edges of DT with endpoints in P.

// examples/Tutorial_SCG99/crust.C
// -------------------------------
#include "tutorial.h"
#include <CGAL/IO/leda_window.h>
#include <CGAL/IO/Ostream_iterator.h>
#include <vector>

struct Tmp_point : public Point {
    static bool init;
    bool delaunay;
    Tmp_point() : delaunay( init) {}
    Tmp_point( const Point& p) : Point(p), delaunay( init) {}
};
bool Tmp_point::init = true;

struct Traits : public EucliTraits {
    typedef Tmp_point Point;
};
typedef CGAL::Triangulation_vertex_base_2<Traits>                  Vb;
typedef CGAL::Triangulation_face_base_2<Traits>                    Fb;
typedef CGAL::Triangulation_default_data_structure_2<Traits,Vb,Fb> Tds;
typedef CGAL::Delaunay_triangulation_2<Traits,Tds>                 Delaunay;


template <class InputIterator, class OutputIterator>
void crust( InputIterator first, InputIterator last, OutputIterator out) {

//generate DT of point set
    Delaunay dt;
    Tmp_point::init = true;
    dt.insert( first, last);


//get voronoi vertices from the DT
    vector<Point> voronoi_vertices;
    voronoi_vertices.reserve( dt.number_of_faces());
    Delaunay::Face_iterator f = dt.faces_begin();
    while ( f != dt.faces_end()) {
	voronoi_vertices.push_back( dt.dual( f));
	++f;
    }


//add the VD vertices into the DT

    Tmp_point::init = false;
    dt.insert( voronoi_vertices.begin(), voronoi_vertices.end());


//find the edges with both endpoints in P

    Delaunay::Edge_iterator e = dt.edges_begin();
    while ( e != dt.edges_end()) {
	Delaunay::Face_handle face = (*e).first;
	int index = (*e).second;
	if ( face->vertex( dt.cw(  index))->point().delaunay &&
	     face->vertex( dt.ccw( index))->point().delaunay)
	    *out++ = dt.segment(e);
        ++e;
    }
}

int main () {
    Delaunay_triangulation  dt;
    leda_window* window = CGAL::create_and_display_demo_window();
    vector<Point> points;
    while (true) {
	Point p;
	*window >> p;
	if ( ! (*window))
	    break;
	points.push_back( p);
	window->clear();
	*window << CGAL::BLACK;
	std::copy( points.begin(), points.end(), 
		   CGAL::Ostream_iterator<Point, leda_window>(*window));
	*window << CGAL::GREEN;
	crust( points.begin(), points.end(), 
	       CGAL::Ostream_iterator<Segment, leda_window>(*window));
    }
    delete window;
    return 0;
}