voronoi_vertex_sqrt_field_c2.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/Voronoi_vertex_sqrt_field_C2.h $
// $Id: Voronoi_vertex_sqrt_field_C2.h 37237 2007-03-19 07:47:35Z afabri $
// 
//
// Author(s)     : Menelaos Karavelas <[e-mail]>




#ifndef CGAL_SEGMENT_DELAUNAY_GRAPH_2_VORONOI_VERTEX_SQRT_FIELD_C2_H
#define CGAL_SEGMENT_DELAUNAY_GRAPH_2_VORONOI_VERTEX_SQRT_FIELD_C2_H



#include <CGAL/Segment_Delaunay_graph_2/Basic_predicates_C2.h>
#include <CGAL/Segment_Delaunay_graph_2/Are_same_points_C2.h>
#include <CGAL/Segment_Delaunay_graph_2/Are_same_segments_C2.h>


CGAL_BEGIN_NAMESPACE

CGAL_SEGMENT_DELAUNAY_GRAPH_2_BEGIN_NAMESPACE


template<class K>
class Voronoi_vertex_sqrt_field_C2
  : public Basic_predicates_C2<K>
{
public:
  typedef Basic_predicates_C2<K> Base;

  typedef enum {PPP = 0, PPS, PSS, SSS} vertex_t;
  struct PPP_Type {};
  struct PPS_Type {};
  struct PSS_Type {};
  struct SSS_Type {};

  typedef typename Base::Point_2             Point_2;
  typedef typename Base::Segment_2           Segment_2;
  typedef typename Base::Line_2              Line_2;
  typedef typename Base::Site_2              Site_2;
  typedef typename Base::FT                  FT;
  typedef typename Base::RT                  RT;

  typedef typename Base::Homogeneous_point_2 Homogeneous_point_2;

  typedef typename Base::Orientation         Orientation;
  typedef typename Base::Comparison_result   Comparison_result;
  typedef typename Base::Oriented_side       Oriented_side;
  typedef typename Base::Sign                Sign;

private:
  typedef Are_same_points_C2<K>    Are_same_points_2;
  typedef Are_same_segments_C2<K>  Are_same_segments_2;

  typedef typename K::Intersections_tag ITag;

  Are_same_points_2     same_points;
  Are_same_segments_2   same_segments;


private:
  //--------------------------------------------------------------------------

  void
  compute_ppp(const Site_2& sp, const Site_2& sq, const Site_2& sr)
  {
    CGAL_precondition( sp.is_point() && sq.is_point() &&
		       sr.is_point() );

    Point_2 p = sp.point(), q = sq.point(), r = sr.point();

    v_type = PPP;

    FT np = CGAL::square(p.x()) + CGAL::square(p.y());
    FT nq = CGAL::square(q.x()) + CGAL::square(q.y());
    FT nr = CGAL::square(r.x()) + CGAL::square(r.y());

    ux = np * (q.y() - r.y()) + nq * (r.y() - p.y()) + nr * (p.y() - q.y());
    uy = -(np * (q.x() - r.x()) + nq * (r.x() - p.x()) + nr * (p.x() - q.x()));
    uz = FT(2) * ( (q.x() * r.y() - r.x() * q.y()) +
		   (r.x() * p.y() - p.x() * r.y()) +
		   (p.x() * q.y() - q.x() * p.y()) );
  }

  //--------------------------------------------------------------------------

  void
  compute_pss(const Site_2& p, const Site_2& q, const Site_2& r)
  {
    CGAL_precondition( p.is_point() && q.is_segment() &&
		       r.is_segment() );

    v_type = PSS;

    bool pq =
      same_points(p, q.source_site()) || same_points(p, q.target_site());
    bool pr =
      same_points(p, r.source_site()) || same_points(p, r.target_site());

    Point_2 pp = p.point();

    if ( pq && pr ) {
      ux = pp.x();
      uy = pp.y();
      uz = FT(1);
      return;
    }


    FT a1, b1, c1, a2, b2, c2;
    compute_supporting_line(q.supporting_site(), a1, b1, c1);
    compute_supporting_line(r.supporting_site(), a2, b2, c2);

    FT c1_ = a1 * pp.x() + b1 * pp.y() + c1;
    FT c2_ = a2 * pp.x() + b2 * pp.y() + c2;

    if ( pq ) {
      c1_ = FT(0);
    }

    if ( pr ) {
      c2_ = FT(0);
    }

    Sign sgn_c1_ = CGAL::sign(c1_);
    Sign sgn_c2_ = CGAL::sign(c2_);

    if ( sgn_c1_ == NEGATIVE ) {
      a1 = -a1;  b1 = -b1;  c1_ = -c1_;
    } else if ( sgn_c1_ == ZERO ) {

      CGAL_assertion( pq );
      if ( same_points(p, q.target_site()) ) {
	a1 = -a1;  b1 = -b1;  c1_ = -c1_;
      }
    }

    if ( sgn_c2_ == NEGATIVE ) {
      a2 = -a2;  b2 = -b2;  c2_ = -c2_;
    } else if ( sgn_c2_ == ZERO ) {

      CGAL_assertion( pr );
      if ( same_points(p, r.source_site()) ) {
	a2 = -a2;  b2 = -b2;  c2_ = -c2_;
      }
    }

    if ( pq ) {
      FT J = a1 * c2_;
      FT I = b1 * c2_;

      FT n1 = CGAL::square(a1) + CGAL::square(b1);
      FT n2 = CGAL::square(a2) + CGAL::square(b2);

      FT D1D2 = n1 * n2;

      uz = -a1 * a2 - b1 * b2 + CGAL::sqrt(D1D2);

      ux = J + pp.x() * uz;
      uy = I + pp.y() * uz;

    } else if ( pr ) {
      FT J = a2 * c1_;
      FT I = b2 * c1_;

      FT n1 = CGAL::square(a1) + CGAL::square(b1);
      FT n2 = CGAL::square(a2) + CGAL::square(b2);

      FT D1D2 = n1 * n2;

      uz = -a1 * a2 - b1 * b2 + CGAL::sqrt(D1D2);

      ux = J + pp.x() * uz;
      uy = I + pp.y() * uz;

    } else {
      Line_2 lq(a1, b1, c1_);
      Line_2 lr(a2, b2, c2_);
      compute_pll(pp, lq, lr);
    }
  }


  void
  compute_pll(const Point_2& p, const Line_2& lq, const Line_2& lr)
  {
    FT a1 = lq.a(), b1 = lq.b(), c1_ = lq.c();
    FT a2 = lr.a(), b2 = lr.b(), c2_ = lr.c();

    CGAL_precondition( c1_ >= FT(0) );
    CGAL_precondition( c2_ >= FT(0) );

    FT n1 = CGAL::square(a1) + CGAL::square(b1);
    FT n2 = CGAL::square(a2) + CGAL::square(b2);

    FT I = b1 * c2_ + b2 * c1_;
    FT J = a1 * c2_ + a2 * c1_;

    FT c1c2 = FT(2) * c1_ * c2_;
    FT a1a2 = a1 * a2;
    FT b1b2 = b1 * b2;

    FT D1D2 = n1 * n2;

    // compute sigma
    FT sigma_expr = b1 * CGAL::sqrt(n2) - b2 * CGAL::sqrt(n1);
    Sign s_sigma = CGAL::sign(sigma_expr);

    int i_sigma(s_sigma);
    FT sigma(i_sigma);

    // compute rho
    FT rho_expr = a1 * CGAL::sqrt(n2) - a2 * CGAL::sqrt(n1);
    Sign s_rho = CGAL::sign(rho_expr);

    int i_rho(s_rho);
    FT rho(i_rho);

    FT sqrt_D1D2 = CGAL::sqrt(D1D2);

    FT A = a1a2 - b1b2;
    FT u1 = c1c2 * (sqrt_D1D2 + A);
    FT u2 = c1c2 * (sqrt_D1D2 - A);

    uz = -a1a2 - b1b2 + CGAL::sqrt(D1D2);

    ux = J + p.x() * uz + sigma * CGAL::sqrt(u1);
    uy = I + p.y() * uz - rho * CGAL::sqrt(u2);
  }


  //--------------------------------------------------------------------------


  void
  compute_pps(const Site_2& p, const Site_2& q, const Site_2& r)
  {
    CGAL_precondition( p.is_point() && q.is_point() &&
		       r.is_segment() );
    v_type = PPS;

    FT a, b, c;
    compute_supporting_line(r.supporting_site(), a, b, c);

    Point_2 pp = p.point(), qq = q.point();

    FT c_ = a * pp.x() + b * pp.y() + c;
    FT cq_ = a * qq.x() + b * qq.y() + c;

    if ( same_points(p, r.source_site()) ||
	 same_points(p, r.target_site()) ) {
      c_ = FT(0);
    }
    if ( same_points(q, r.source_site()) ||
	 same_points(q, r.target_site()) ) {
      cq_ = FT(0);
    }

    Sign s = CGAL::sign(c_);

    if ( s == NEGATIVE ) {
      a = -a;  b = -b;  c = -c;  c_ = -c_;  cq_ = -cq_;
    } else if ( s == ZERO ) {
      Sign s1 = CGAL::sign(cq_);

      CGAL_assertion( s1 != ZERO );
      if ( s1 == NEGATIVE ) {
	a = -a;  b = -b;  c = -c;  c_ = -c_;  cq_ = -cq_;
      }
    }

    FT nl = CGAL::square(a) + CGAL::square(b);

    FT x_ = qq.x() - pp.x();
    FT y_ = qq.y() - pp.y();
    FT n_ = CGAL::square(x_) + CGAL::square(y_);


    Comparison_result res = CGAL::compare( c_, cq_ );

    if ( res == EQUAL ) {
      FT e1 = CGAL::square(c_);
      FT J = nl * (a * n_ + FT(4) * c_ * x_) - FT(4) * a * e1;
      FT I = nl * (b * n_ + FT(4) * c_ * y_) - FT(4) * b * e1;
      FT X = FT(8) * nl * c_;

      ux = J + pp.x() * X;
      uy = I + pp.y() * X;
      uz = X;
      return;
    }

    FT e1 = a * x_ + b * y_;
    FT e2 = b * x_ - a * y_;
    FT e3 = n_ * e1;
    FT e4 = FT(2) * c_ * e2;

    FT X = FT(2) * CGAL::square(e1);
    FT I = b * e3 + x_ * e4;
    FT J = a * e3 - y_ * e4;
    FT sqrt_S = CGAL::sqrt(n_ * nl * c_ * cq_);

    ux = J + pp.x() * X - FT(2) * y_ * sqrt_S;
    uy = I + pp.y() * X + FT(2) * x_ * sqrt_S;
    uz = X;
  }


  //--------------------------------------------------------------------------
  bool check_if_exact(const Site_2& , unsigned int ,
		      const Tag_false&) const
  {
    return true;
  } 

  bool check_if_exact(const Site_2& s, unsigned int i,
		      const Tag_true&) const
  {
    return s.is_input(i);
  }

  // determines of the segment s is on the positive halfspace as
  // defined by the supporting line of the segment supp; the line l
  // is supposed to be the supporting line of the segment supp and we
  // pass it so that we do not have to recompute it
  bool
  is_on_positive_halfspace(const Site_2& supp,
			   const Site_2& s, const Line_2& l) const
  {
    CGAL_precondition( supp.is_segment() && s.is_segment() );

    if ( same_segments(supp.supporting_site(),
		       s.supporting_site()) ) {
      return false;
    }

    if ( same_points(supp.source_site(), s.source_site()) ||
	 same_points(supp.target_site(), s.source_site()) ) {
      return oriented_side_of_line(l, s.target()) == ON_POSITIVE_SIDE;
    }

    if ( same_points(supp.source_site(), s.target_site()) ||
	 same_points(supp.target_site(), s.target_site()) ) {
      return oriented_side_of_line(l, s.source()) == ON_POSITIVE_SIDE;
    }

    ITag itag;

    if ( !check_if_exact(s, 0, itag) &&
	 same_segments(supp.supporting_site(),
		       s.crossing_site(0)) ) {
      return oriented_side_of_line(l, s.target()) == ON_POSITIVE_SIDE;
    }

    if ( !check_if_exact(s, 1, itag) &&
	 same_segments(supp.supporting_site(),
		       s.crossing_site(1)) ) {
      return oriented_side_of_line(l, s.source()) == ON_POSITIVE_SIDE;
    }

    return Base::is_on_positive_halfspace(l, s.segment());
  }

  //--------------------------------------------------------------------------

  void
  orient_lines(const Site_2& sp, const Site_2& sq,
	       const Site_2& sr, FT a[], FT b[], FT c[]) const 
  {
    CGAL_precondition( sp.is_segment() && sq.is_segment() &&
		       sr.is_segment() );

    Line_2 l[3];
    l[0] = compute_supporting_line(sp.supporting_site());
    l[1] = compute_supporting_line(sq.supporting_site());
    l[2] = compute_supporting_line(sr.supporting_site());
    
    bool is_oriented[3] = {false, false, false};

    if ( is_on_positive_halfspace(sp, sq, l[0]) ||
    	 is_on_positive_halfspace(sp, sr, l[0]) ) {
      is_oriented[0] = true;
    } else {
      
      l[0] = opposite_line(l[0]);
      if ( is_on_positive_halfspace(sp, sq, l[0]) ||
      	   is_on_positive_halfspace(sp, sr, l[0]) ) {
	is_oriented[0] = true;
      } else {
	l[0] = opposite_line(l[0]);
      }
    }

    if ( is_on_positive_halfspace(sq, sp, l[1]) ||
	 is_on_positive_halfspace(sq, sr, l[1]) ) {
      is_oriented[1] = true;
    } else {
      l[1] = opposite_line(l[1]);
      if ( is_on_positive_halfspace(sq, sp, l[1]) ||
	   is_on_positive_halfspace(sq, sr, l[1]) ) {
	is_oriented[1] = true;
      } else {
	l[1] = opposite_line(l[1]);
      }
    }

    if ( is_on_positive_halfspace(sr, sp, l[2]) ||
	 is_on_positive_halfspace(sr, sq, l[2]) ) {
      is_oriented[2] = true;
    } else {
      l[2] = opposite_line(l[2]);
      if ( is_on_positive_halfspace(sr, sp, l[2]) ||
	   is_on_positive_halfspace(sr, sq, l[2]) ) {
	is_oriented[2] = true;
      } else {
	l[2] = opposite_line(l[2]);
      }
    }

    if ( is_oriented[0] && is_oriented[1] && is_oriented[2] ) {
      for (int i = 0; i < 3; i++) {
	a[i] = l[i].a();
	b[i] = l[i].b();
	c[i] = l[i].c();
      }
      return;
    }

    int i_no(-1);
    for (int i = 0; i < 3; i++) {
      if ( !is_oriented[i] ) {
	i_no = i;
	CGAL_assertion( is_oriented[(i+1)%3] && is_oriented[(i+2)%3] );
	break;
      }
    }

    CGAL_assertion( i_no != -1 );

    FT sqrt_d[3];
    for (int i = 0; i < 3; i++) {
      FT d1 = CGAL::square(l[i].a()) + CGAL::square(l[i].b());
      sqrt_d[i] = CGAL::sqrt(d1);
    }

    FT z[3];
    for (int i = 0; i < 3; i++) {
      z[i] = l[(i+1)%3].a() * l[(i+2)%3].b()
	- l[(i+2)%3].a() * l[(i+1)%3].b();
    }


    FT vz = z[0] * sqrt_d[0] + z[1] * sqrt_d[1] + z[2] * sqrt_d[2];

    Sign s_minus_vz = CGAL::sign(vz);

    CGAL_assertion( s_minus_vz != ZERO );

    if ( s_minus_vz == NEGATIVE ) {
      l[i_no] = opposite_line(l[i_no]);

      for (int i = 0; i < 3; i++) {
	a[i] = l[i].a();
	b[i] = l[i].b();
	c[i] = l[i].c();
      }
      return;
    }

    // now we have to check if the other orientation of l[i_no]
    // corresponds to a CCW triangle as well.
    z[(i_no+1)%3] = -z[(i_no+1)%3];
    z[(i_no+2)%3] = -z[(i_no+2)%3];

    vz = z[0] * sqrt_d[0] + z[1] * sqrt_d[1] + z[2] * sqrt_d[2];

    Sign s_minus_vz_2 = CGAL::sign(vz);

    CGAL_assertion( s_minus_vz_2 != ZERO );

    if ( s_minus_vz_2 == NEGATIVE ) {
      // the other orientation does not correspond to a CCW triangle.
      for (int i = 0; i < 3; i++) {
	a[i] = l[i].a();
	b[i] = l[i].b();
	c[i] = l[i].c();
      }
      return;
    }

    // now compute the Voronoi vertex;
    for (int i = 0; i < 3; i++) {
      a[i] = l[i].a();
      b[i] = l[i].b();
      c[i] = l[i].c();
    }

    FT x[3], y[3], w[3];
    for (int i = 0; i < 3; i++) {
      x[i] = c[(i+1)%3] * b[(i+2)%3] - c[(i+2)%3] * b[(i+1)%3];
      y[i] = -(c[(i+1)%3] * a[(i+2)%3] - c[(i+2)%3] * a[(i+1)%3]);
      w[i] = -(a[(i+1)%3] * b[(i+2)%3] - a[(i+2)%3] * b[(i+1)%3]);
    }

    FT vx = FT(0), vy = FT(0), vw = FT(0);

    for (int i = 0; i < 3; i++) {
      vx += x[i] * sqrt_d[i];
      vy += y[i] * sqrt_d[i];
      vw += w[i] * sqrt_d[i];
    }

    Sign s_vw = CGAL::sign(vw);

    FT dist =
      a[(i_no+1)%3] * vx + b[(i_no+1)%3] * vy + c[(i_no+1)%3] * vw;
    

#ifdef CGAL_CFG_NO_OPERATOR_TIMES_FOR_SIGN
    Sign sgn_dist = CGAL::Sign(s_vw * CGAL::sign(dist));
#else
    Sign sgn_dist = s_vw * CGAL::sign(dist);
#endif

    CGAL_assertion( sgn_dist != ZERO );

    if ( sgn_dist == NEGATIVE ) {
      a[i_no] = -a[i_no];
      b[i_no] = -b[i_no];
      c[i_no] = -c[i_no];
    }
  }



  void
  compute_sss(const Site_2& p, const Site_2& q, const Site_2& r)
  {
    CGAL_precondition( p.is_segment() && q.is_segment() &&
		       r.is_segment() );

    v_type = SSS;

    FT a[3], b[3], c[3];
    FT cx[3], cy[3], cz[3], sqrt_D[3];

    orient_lines(p, q, r, a, b, c);

    for (int i = 0; i < 3; i++) {
      cx[i] = c[(i+1)%3] * b[(i+2)%3] - c[(i+2)%3] * b[(i+1)%3];
      cy[i] = -(c[(i+1)%3] * a[(i+2)%3] - c[(i+2)%3] * a[(i+1)%3]);
      cz[i] = -(a[(i+1)%3] * b[(i+2)%3] - a[(i+2)%3] * b[(i+1)%3]);

      FT d = CGAL::square(a[i]) + CGAL::square(b[i]);
      sqrt_D[i] = CGAL::sqrt(d);
    }

    ux = cx[0] * sqrt_D[0] + cx[1] * sqrt_D[1] + cx[2] * sqrt_D[2];
    uy = cy[0] * sqrt_D[0] + cy[1] * sqrt_D[1] + cy[2] * sqrt_D[2];
    uz = cz[0] * sqrt_D[0] + cz[1] * sqrt_D[1] + cz[2] * sqrt_D[2];
  }

  //--------------------------------------------------------------------------

  void
  compute_vertex(const Site_2& s1, const Site_2& s2, const Site_2& s3)
  {
    if ( s1.is_point() && s2.is_point() && s3.is_point() ) {
      compute_ppp(s1, s2, s3);