---

delaunay.cpp

//$Header: /home/ben/Mapper/c++/RCS/delaunay.cpp,v 6.8 2002/06/14 22:28:47 ben Exp $
// Copyright Benedict Adamson 2002.
// This file is part of Mapper.

// Mapper 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.

// Mapper 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 General Public License for more details.

// You should have received a copy of the GNU General Public License
// along with Mapper; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include <algorithm>
#include <cassert>
#include <cstdlib>

#include "delaunay.h"
#include "matrix2x2.h"



delaunay::triangle::triangle(
   const vector2 &v0,
   const vector2 &v1,
   const vector2 &v2
   )
{
   this->vertex_[0] = v0;
   this->vertex_[1] = v1;
   this->vertex_[2] = v2;
   this->calculate_circumcircle();
   {//Postconditions:
      assert( this->vertex(0) == v0 );
      assert( this->vertex(1) == v1 );
      assert( this->vertex(2) == v2 );
   };
};



delaunay::triangle::triangle(
   const delaunay::triangle::triangle &t
   )
   :
   circumcentre_( t.circumcentre_ ),
   circumradius2( t.circumradius2 )
{
   this->vertex_[0] = t.vertex_[0];
   this->vertex_[1] = t.vertex_[1];
   this->vertex_[2] = t.vertex_[2];
}

const double delaunay::triangle::circumradius(void) const
{
   return sqrt(this->circumradius2);
}



bool delaunay::triangle::circumcircle_contains(
   const vector2 &p
   ) const
{
   const vector2 rp( p - this->circumcentre_ );
   return rp*rp <= this->circumradius2;
}



void delaunay::triangle::calculate_circumcircle(void)
{
   const vector2 &v0( this->vertex(0) );
   const vector2 &v1( this->vertex(1) );
   const vector2 &v2( this->vertex(2) );
   {//Calculate the circumcentre
      const vector2 d01( v1 - v0 );
      const vector2 d02( v2 - v0 );
      const vector2 d12( v2 - v1 );
      const vector2 m01( 0.5*(v0 + v1) );
      const vector2 m02( 0.5*(v0 + v2) );
      const vector2 r01( -d01.y, d01.x );
      const vector2 r02( -d02.y, d02.x );
      const matrix2x2 M(
         r01.x, d02.y,
         d01.x, -r02.y
         );
      const vector2 y( 0.5*d12 );
      const vector2 a( matrix2x2::solve( M, y ) );
      const vector2 c01 = m01 + a[0] * r01;
      const vector2 c02 = m02 + a[1] * r02;
      this->circumcentre_ = 0.5*(c01 + c02);
   };
   {//Calculate the square of the circumradius
      const vector2 r0 = v0 - this->circumcentre_;
      this->circumradius2 = r0*r0;
   };
}



delaunay::triangulation::triangulation(
   const double l,
   const double r,
   const double t,
   const double b
   )
   :
   left( l ),
   right( r ),
   top( t ),
   bottom( b ),
   tl( l, t ),
   tr( r, t ),
   br( r, b ),
   bl( l, b )
{
   this->triangles_.reserve( 2 );
   this->triangles_.push_back( triangle( tl, tr, bl ) );
   this->triangles_.push_back( triangle( bl, tr, br ) );
}



delaunay::triangulation::~triangulation(void)
{
   //Do nothing
}



void delaunay::triangulation::add_point(
   const vector2 &p
   )
{
   {//Preconditions:
      assert( this->left <= p.x && p.x <= this->right );
      assert( this->bottom <= p.y && p.y <= this->top );
      assert( !this->contains( p ) );
   };
   vector< delaunay::triangle > new_triangles;
   vector< delaunay::edge > edges;
   {//heuristic: reserve space
      static const unsigned expected_intersect = 10u;
      new_triangles.reserve( this->triangles_.size() + expected_intersect );
      edges.reserve( expected_intersect * 3u );
   };
   for( vector< delaunay::triangle >::iterator i = this->triangles_.begin();
        i != this->triangles_.end();
        i++ ){//scan existing triangles
      if( (*i).circumcircle_contains( p ) ){
         //triangle will be changed
         edges.push_back( edge( (*i).vertex(0), (*i).vertex(1) ) );
         edges.push_back( edge( (*i).vertex(1), (*i).vertex(2) ) );
         edges.push_back( edge( (*i).vertex(2), (*i).vertex(0) ) );
      }else{
         //triangle not changed
         new_triangles.push_back( *i );
      };
   };
   {//eliminate edges that are duplicated
      vector< delaunay::edge > new_edges;
      new_edges.reserve( edges.size() ); //heuristic
      vector< delaunay::edge >::iterator begin = edges.begin();
      vector< delaunay::edge >::iterator end = edges.end();
      for( vector< delaunay::edge >::iterator i = begin;
           i != end;
           i++ ){
         edge &e( *i );
         bool duplicate = false;
         {//determine whether duplicate
            //Not efficient, but N is usually small
            vector< delaunay::edge >::iterator j = begin;
            while( !duplicate && j != end ){
               if( j != i ){
                  duplicate = (e == *j);
               };//else skip self
               j++;
            };
         };
         if( !duplicate ){
            new_edges.push_back( e );
         };
      };
      edges.swap( new_edges );
   };
   {//create the new triangles
      for( vector< delaunay::edge >::iterator i = edges.begin();
           i != edges.end();
           i++ ){
         edge &e( *i );
         new_triangles.push_back(
            triangle( e.vertex(0), e.vertex(1), p )
            );
      };
   }
   this->triangles_.swap( new_triangles );
   {//Postconditions:
      assert( this->contains( p ) );
   };
}



bool delaunay::triangulation::contains(
   const vector2 &p
   ) const
{
   //assert( this->contains( this->tl ) );
   //assert( this->contains( this->tr ) );
   //assert( this->contains( this->br ) );
   //assert( this->contains( this->bl ) );
   //Indicates whether the given point is a vertex
   //of any of the traingles of the triangulation.
   register bool c = false;
   vector< triangle >::const_iterator i = this->triangles_.begin();
   const vector< triangle >::const_iterator end = this->triangles_.end();
   while( !c && i != end ){
      const triangle &t( *i );
      c = (t.vertex(0) == p) ||
         (t.vertex(1) == p ) ||
         (t.vertex(2) == p);
      i++;
   };
   return c;
}

---