segment_delaunay_graph_2_impl.h

很多二维 三维几何计算算法 C++ 类库

// Copyright (c) 2003,2004,2005,2006  INRIA Sophia-Antipolis (France) and
// Notre Dame University (U.S.A.).  All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Segment_Delaunay_graph_2/include/CGAL/Segment_Delaunay_graph_2/Segment_Delaunay_graph_2_impl.h $
// $Id: Segment_Delaunay_graph_2_impl.h 37157 2007-03-16 10:49:14Z afabri $
// 
//
// Author(s)     : Menelaos Karavelas <[e-mail]>


// class implementation continued
//=================================

CGAL_BEGIN_NAMESPACE

//====================================================================
//====================================================================
//                   CONSTRUCTORS
//====================================================================
//====================================================================

// copy constructor
template<class Gt, class ST, class DS, class LTag>
Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::
Segment_Delaunay_graph_2(const Segment_Delaunay_graph_2& other)
  : DG(other.geom_traits())
{
  Segment_Delaunay_graph_2&
    non_const_other = const_cast<Segment_Delaunay_graph_2&>(other);
  copy(non_const_other);
  CGAL_postcondition( is_valid() );
}

// assignment operator
template<class Gt, class ST, class DS, class LTag>
typename Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::Self&
Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::
operator=(const Self& other)
{
  if ( this != &other ) {
    Segment_Delaunay_graph_2&
      non_const_other = const_cast<Segment_Delaunay_graph_2&>(other);
    copy(non_const_other);
  }
  return (*this);
}

//====================================================================
//====================================================================
//                   METHODS FOR INSERTION
//====================================================================
//====================================================================

//--------------------------------------------------------------------
// insertion of first three sites
//--------------------------------------------------------------------

template<class Gt, class ST, class DS, class LTag>
typename Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::Vertex_handle
Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::
insert_first(const Storage_site_2& ss, const Point_2& )
{
  CGAL_precondition( number_of_vertices() == 0 );

  Vertex_handle v = this->_tds.insert_second();
  v->set_site(ss);
  return v;
}

template<class Gt, class ST, class DS, class LTag>
typename Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::Vertex_handle
Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::
insert_second(const Storage_site_2& ss, const Point_2& p)
{
  CGAL_precondition( number_of_vertices() == 1 );
  // p0 is actually a point
  Site_2 p0 = finite_vertices_begin()->site();
  // MK: change the equality test between points by the functor in
  // geometric traits
  Site_2 tp = Site_2::construct_site_2(p);
  if ( same_points(tp,p0) ) {
    // merge info of identical sites
    merge_info(finite_vertices_begin(), ss);
    return finite_vertices_begin();
  }

  return create_vertex_dim_up(ss);
}

template<class Gt, class ST, class DS, class LTag>
typename Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::Vertex_handle
Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::
insert_third(const Storage_site_2& ss, const Point_2& p)
{
  Site_2 t = Site_2::construct_site_2(p);
  return insert_third(t, ss);
}

template<class Gt, class ST, class DS, class LTag>
typename Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::Vertex_handle
Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::
insert_third(const Site_2& t, const Storage_site_2& ss)
{
  CGAL_precondition( number_of_vertices() == 2 );

  // p0 and p1 are actually points
  Vertex_handle v0 = finite_vertices_begin();
  Vertex_handle v1 = ++finite_vertices_begin();
  Site_2 t0 = v0->site();
  Site_2 t1 = v1->site();

  if ( same_points(t, t0) ) {
    // merge info of identical sites
    merge_info(v0, ss);
    return v0;
  }
  if ( same_points(t, t1) ) {
    // merge info of identical sites
    merge_info(v1, ss);
    return v1;
  }

  Vertex_handle v = create_vertex_dim_up(ss);

  Face_handle f(finite_faces_begin());

  Site_2 s1 = f->vertex(0)->site();
  Site_2 s2 = f->vertex(1)->site();
  Site_2 s3 = f->vertex(2)->site();

  Orientation o =
    geom_traits().orientation_2_object()(s1, s2, s3);

  if ( o != COLLINEAR ) {
    if ( o == RIGHT_TURN ) {
      f->reorient();
      for (int i = 0; i < 3; i++) {
	f->neighbor(i)->reorient();
      }
    }
  } else {
    typename Geom_traits::Compare_x_2 compare_x =
      geom_traits().compare_x_2_object();

    Comparison_result xcmp12 = compare_x(s1, s2);
    if ( xcmp12 == SMALLER ) {        // x1 < x2
      Comparison_result xcmp23 = compare_x(s2, s3);
      if ( xcmp23 == SMALLER ) {            // x2 < x3
	flip(f, f->index(v1));
      } else {
	Comparison_result xcmp31 = compare_x(s3, s1);
	if ( xcmp31 == SMALLER ) {          // x3 < x1
	  flip(f, f->index(v0));
	} else {                            // x1 < x3 < x2
	  flip(f, f->index(v)); 
	}
      }
    } else if ( xcmp12 == LARGER ) {  // x1 > x2
      Comparison_result xcmp32 = compare_x(s3, s2);
      if ( xcmp32 == SMALLER ) {            // x3 < x2
	flip(f, f->index(v1));
      } else {
	Comparison_result xcmp13 = compare_x(s1, s3);
	if ( xcmp13 == SMALLER ) {          // x1 < x3
	  flip(f, f->index(v0));
	} else {                            // x2 < x3 < x1
	  flip(f, f->index(v));
	}
      }
    } else {                          // x1 == x2
      typename Geom_traits::Compare_y_2 compare_y =
	geom_traits().compare_y_2_object();

      Comparison_result ycmp12 = compare_y(s1, s2);
      if ( ycmp12 == SMALLER ) {      // y1 < y2
	Comparison_result ycmp23 = compare_y(s2, s3);
	if ( ycmp23 == SMALLER ) {          // y2 < y3
	  flip(f, f->index(v1));
	} else {
	  Comparison_result ycmp31 = compare_y(s3, s1);
	  if ( ycmp31 == SMALLER ) {        // y3 < y1
	    flip(f, f->index(v0));
	  } else {                          // y1 < y3 < y2
	    flip(f, f->index(v));
	  }
	}
      } else if ( ycmp12 == LARGER ) { // y1 > y2
	Comparison_result ycmp32 = compare_y(s3, s2);
	if ( ycmp32 == SMALLER ) {           // y3 < y2
	  flip(f, f->index(v1));
	} else {
	  Comparison_result ycmp13 = compare_y(s1, s3);
	  if ( ycmp13 == SMALLER ) {         // y1 < y3
	    flip(f, f->index(v0));
	  } else {                           // y2 < y3 < y1
	    flip(f, f->index(v));
	  }
	}
      } else {
	// this line should never have been reached
	CGAL_assertion( false );
      }
    }
  }

  return v;
}


template<class Gt, class ST, class DS, class LTag>
typename Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::Vertex_handle
Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::
insert_third(const Storage_site_2& ss, Vertex_handle , Vertex_handle )
{
  CGAL_precondition( number_of_vertices() == 2 );

  //  this can only be the case if the first site is a segment
  CGAL_precondition( dimension() == 1 );

  Vertex_handle v = create_vertex_dim_up(ss);

  Face_circulator fc = incident_faces(v);

  while ( true ) {
    Face_handle f(fc);
    if ( !is_infinite(f) ) {
      flip(f, f->index(v));
      break;
    }
    ++fc;
  }
  
  return v;
}

//--------------------------------------------------------------------
// insertion of a point
//--------------------------------------------------------------------

template<class Gt, class ST, class DS, class LTag>
typename Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::Vertex_handle
Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::
insert_point(const Storage_site_2& ss, const Point_2& p, Vertex_handle vnear)
{
  int n = number_of_vertices();
  if ( n == 0 ) {
    return insert_first(ss, p);
  } else if ( n == 1 ) {
    return insert_second(ss, p);
  } else if ( n == 2 ) {
    return insert_third(ss, p);
  }

  Site_2 t = Site_2::construct_site_2(p);
  return insert_point(ss, t, vnear);
}


template<class Gt, class ST, class DS, class LTag>
typename Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::Vertex_handle
Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::
insert_point(const Storage_site_2& ss, const Site_2& t,
	     Vertex_handle vnear)
{
  CGAL_precondition( t.is_point() );
  CGAL_assertion( number_of_vertices() > 2 );

  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  // MK::ERROR: I need to write a insert_point_no_search method that
  // does not search for the nearest neighbor; this should be used by
  // insert_point. Below the first version of the code is correct. The
  // second is what the insert_point method should do before calling
  // insert_point_no_search.

  // first find the nearest neighbor
#if 1
  Vertex_handle  vnearest = nearest_neighbor( t, vnear );
#else
  Vertex_handle vnearest;
  if ( vnear == Vertex_handle() ) {
    vnearest = nearest_neighbor( t, vnear );
  } else {
    vnearest = vnear;
  }
#endif

  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************
  //*********************************************************************


  // check is the point has already been inserted or lies on
  // a segment already inserted
  Arrangement_type at_res = arrangement_type(t, vnearest);
  if ( vnearest->is_point() ) {
    if ( at_res == AT2::IDENTICAL ) {
      // merge info of identical sites
      merge_info(vnearest, ss);
      return vnearest;
    }
  } else {
    CGAL_assertion( vnearest->is_segment() );
    CGAL_assertion( at_res != AT2::TOUCH_1 );
    CGAL_assertion( at_res != AT2::TOUCH_2 );
    CGAL_assertion( at_res == AT2::DISJOINT || at_res == AT2::INTERIOR );
    if ( at_res == AT2::INTERIOR ) {
      CGAL_assertion( t.is_input() );

      Vertex_triple vt = insert_exact_point_on_segment(ss, t, vnearest);
      return vt.first;
    } else {
      // the point to be inserted does not belong to the interior of a
      // segment
      CGAL_assertion( at_res == AT2::DISJOINT );
    }
  }

  return insert_point2(ss, t, vnearest);
}

template<class Gt, class ST, class DS, class LTag>
typename Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::Vertex_handle
Segment_Delaunay_graph_2<Gt,ST,DS,LTag>::
insert_point2(const Storage_site_2& ss, const Site_2& t,
	      Vertex_handle vnearest)
{
  CGAL_precondition( t.is_point() );
  CGAL_assertion( number_of_vertices() > 2 );

  CGAL_expensive_precondition
    ( nearest_neighbor(t, vnearest) == vnearest );

  // find the first conflict

#if !defined(CGAL_NO_ASSERTIONS) && !defined(NDEBUG)
  // verify that there are no intersections...
  Vertex_circulator vc = incident_vertices(vnearest);
  Vertex_circulator vc_start = vc;
  do {
    Vertex_handle vv(vc);
    Arrangement_type at_res = arrangement_type(t, vv);

    CGAL_assertion( at_res == AT2::DISJOINT );
    ++vc;
  } while ( vc != vc_start );
#endif

  // first look for conflict with vertex
  Face_circulator fc_start = incident_faces(vnearest);
  Face_circulator fc = fc_start;
  Face_handle start_f;
  Sign s;

  std::map<Face_handle,Sign> sign_map;

  do {
    Face_handle f(fc);

    s = incircle(f, t);

    sign_map[f] = s;

    if ( s == NEGATIVE ) {
      start_f = f;
      break;
    }
    ++fc;
  } while ( fc != fc_start );

  // we are not in conflict with a Voronoi vertex, so we have to
  // be in conflict with the interior of a Voronoi edge
  if ( s != NEGATIVE ) {
    Edge_circulator ec_start = incident_edges(vnearest);
    Edge_circulator ec = ec_start;

    bool interior_in_conflict(false);
    Edge e;
    do {
      e = *ec;

      Sign s1 = sign_map[e.first];
      Sign s2 = sign_map[e.first->neighbor(e.second)];

      if ( s1 == s2 ) {
	interior_in_conflict = edge_interior(e, t, s1);
      } else {
	// It seems that there was a problem here when one of the
	// signs was positive and the other zero. In this case we
	// still check pretending that both signs where positive
	interior_in_conflict = edge_interior(e, t, POSITIVE);
      }

      if ( interior_in_conflict ) { break; }
      ++ec;
    } while ( ec != ec_start );

    sign_map.clear();

    CGAL_assertion( interior_in_conflict );

    return insert_degree_2(e, ss);
  }


  // we are in conflict with a Voronoi vertex; start from that and 
  // find the entire conflict region and then repair the diagram
  List l;
  Face_map fm;

  Triple<bool, Vertex_handle, Arrangement_type>
    vcross(false, Vertex_handle(), AT2::DISJOINT);

  // MK:: NEED TO WRITE A FUNCTION CALLED find_conflict_region WHICH
  // IS GIVEN A STARTING FACE, A LIST, A FACE MAP, A VERTEX MAP AND A
  // LIST OF FLIPPED EDGES AND WHAT IS DOES IS INITIALIZE THE CONFLICT 
  // REGION AND EXPANDS THE CONFLICT REGION.
  initialize_conflict_region(start_f, l);