svd_voronoi_vertex_sqrt_field_h2.h

来自「CGAL is a collaborative effort of severa」· C头文件 代码 · 共 1,028 行 · 第 1/2 页

// Copyright (c) 2003,2004  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.
//
// $Source: /CVSROOT/CGAL/Packages/Segment_Voronoi_diagram_2/include/CGAL/predicates/Svd_Voronoi_vertex_sqrt_field_H2.h,v $
// $Revision: 1.5 $ $Date: 2004/02/23 22:18:43 $
// $Name:  $
//
// Author(s)     : Menelaos Karavelas <[e-mail]>




#ifndef CGAL_SEGMENT_VORONOI_DIAGRAM_VORONOI_VERTEX_SQRT_FIELD_H2_H
#define CGAL_SEGMENT_VORONOI_DIAGRAM_VORONOI_VERTEX_SQRT_FIELD_H2_H



#include <CGAL/enum.h>
#include <CGAL/determinant.h>
#include <CGAL/predicates/Svd_basic_predicates_H2.h>



CGAL_BEGIN_NAMESPACE




template<class K>
class Svd_voronoi_vertex_sqrt_field_H2
  : public Svd_basic_predicates_H2<K>
{
public:
  typedef Svd_basic_predicates_H2<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;

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

  void
  compute_ppp(const Point_2& p, const Point_2& q, const Point_2& r)
  {
    v_type = PPP;

    RT np = CGAL::square(p.hx()) + CGAL::square(p.hy());
    RT nq = CGAL::square(q.hx()) + CGAL::square(q.hy());
    RT nr = CGAL::square(r.hx()) + CGAL::square(r.hy());

    RT p_xw = p.hx() * p.hw();
    RT q_xw = q.hx() * q.hw();
    RT r_xw = r.hx() * r.hw();

    RT p_yw = p.hy() * p.hw();
    RT q_yw = q.hy() * q.hw();
    RT r_yw = r.hy() * r.hw();

    RT p_hw_sq = CGAL::square( p.hw() );
    RT q_hw_sq = CGAL::square( q.hw() );
    RT r_hw_sq = CGAL::square( r.hw() );

    ux = -det3x3_by_formula(p_yw, np, p_hw_sq,
			    q_yw, nq, q_hw_sq,
			    r_yw, nr, r_hw_sq);

    uy = det3x3_by_formula(p_xw, np, p_hw_sq,
			   q_xw, nq, q_hw_sq,
			   r_xw, nr, r_hw_sq);

    uz = RT(2) * det3x3_by_formula(p_xw, p_yw, p_hw_sq,
				   q_xw, q_yw, q_hw_sq,
				   r_xw, r_yw, r_hw_sq);
  }

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

  void
  compute_pss(const Point_2& p, const Segment_2& q, const Segment_2& r)
  {
    v_type = PSS;

    bool pq = (p == q.source() || p == q.target());
    bool pr = (p == r.source() || p == r.target());

    if ( pq && pr ) {
      ux = p.hx();
      uy = p.hy();
      uz = p.hw();
      return;
    }


    RT a1, b1, c1, a2, b2, c2;
    compute_supporting_line(q, a1, b1, c1);
    compute_supporting_line(r, a2, b2, c2);

    RT c1_ = a1 * p.hx() + b1 * p.hy() + c1 * p.hw();
    RT c2_ = a2 * p.hx() + b2 * p.hy() + c2 * p.hw();

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

    if ( pr ) {
      c2_ = RT(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 ( p == q.target() ) {
	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 ( p == r.source() ) {
	a2 = -a2;  b2 = -b2;  c2_ = -c2_;
      }
    }

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

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

      RT D1D2 = n1 * n2;

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

      ux = J + p.hx() * uz;
      uy = I + p.hy() * uz;

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

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

      RT D1D2 = n1 * n2;

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

      ux = J + p.hx() * uz;
      uy = I + p.hy() * uz;

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

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

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

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

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

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

    RT D1D2 = n1 * n2;

    Sign s_pw = CGAL::sign( p.hw() );

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

    int i_sigma(s_sigma);
    RT sigma(i_sigma);

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

    int i_rho(s_rho);
    RT rho(i_rho);

    RT sqrt_D1D2 = CGAL::sqrt(D1D2);

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

    uz = p.hw() * (-a1a2 - b1b2 + CGAL::sqrt(D1D2));

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


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


  void
  compute_pps(const Point_2& p, const Point_2& q, const Segment_2& r)
  {
    v_type = PPS;

    RT a, b, c;
    compute_supporting_line(r, a, b, c);

    RT c_ = a * p.hx() + b * p.hy() + c * p.hw();
    RT cq_ = a * q.hx() + b * q.hy() + c * q.hw();


    if ( p == r.source() || p == r.target() ) {
      c_ = RT(0);
    }
    if ( q == r.source() || q == r.target() ) {
      cq_ = RT(0);
    }

    Sign s_pw = CGAL::sign( p.hw() );
    Sign s_qw = CGAL::sign( q.hw() );

    Sign s = s_pw * CGAL::sign(c_);

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

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


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

    RT x_ = q.hx() * p.hw() - p.hx() * q.hw();
    RT y_ = q.hy() * p.hw() - p.hy() * q.hw();
    RT n_ = CGAL::square(x_) + CGAL::square(y_);

    // old equality of distcne check
    // ------------------------------
    //    Comparison_result res = CGAL::compare( c_, cq_);

    // checking equality of distance through the quantity X:
    Comparison_result res = CGAL::compare(a * x_, -b * y_);

    if ( res == EQUAL ) {
      RT qw_sq = CGAL::square( q.hw() );
      RT e1 = CGAL::square(c_);
      RT J = nl * (a * n_ + RT(4) * c_ * x_ * q.hw())
	- RT(4) * a * e1 * qw_sq;
      RT I = nl * (b * n_ + RT(4) * c_ * y_ * q.hw())
	- RT(4) * b * e1 * qw_sq;
      RT X = RT(8) * nl * c_ * qw_sq * p.hw();

      ux = J + p.hx() * X;
      uy = I + p.hy() * X;
      uz = p.hw() * X;
      return;
    }

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

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

    ux = J + p.hx() * X - FT(2) * y_ * sqrt_S;
    uy = I + p.hy() * X + FT(2) * x_ * sqrt_S;
    uz = p.hw() * X;
  }


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

  bool
  is_consistent(const Segment_2& p, const Segment_2& q,
		const Segment_2& r, RT a[], RT b[], RT c[]) const
  {
    Line_2 l[3] = {Line_2(a[0], b[0], c[0]), Line_2(a[1], b[1], c[1]), 
		   Line_2(a[2], b[2], c[2])};

    int num_oriented(0);

    if ( is_on_positive_halfspace(l[0], q) ||
	 is_on_positive_halfspace(l[0], r) ) {
      num_oriented++;
    }

    if ( is_on_positive_halfspace(l[1], p) ||
	 is_on_positive_halfspace(l[1], r) ) {
      num_oriented++;
    }

    if ( is_on_positive_halfspace(l[2], p) ||
	 is_on_positive_halfspace(l[2], q) ) {
      num_oriented++;
    }

    return ( num_oriented >= 2 );
  }

  void
  orient_lines(const Segment_2& p, const Segment_2& q,
	       const Segment_2& r, RT a[], RT b[], RT c[]) const 
  {
    Line_2 l[3];
    l[0] = compute_supporting_line(p);
    l[1] = compute_supporting_line(q);
    l[2] = compute_supporting_line(r);
    
    bool is_oriented[3] = {false, false, false};

    if ( is_on_positive_halfspace(l[0], q) ||
	 is_on_positive_halfspace(l[0], r) ) {
      is_oriented[0] = true;
    } else {
      l[0] = opposite_line(l[0]);
      if ( is_on_positive_halfspace(l[0], q) ||
	   is_on_positive_halfspace(l[0], r) ) {
	is_oriented[0] = true;
      } else {
	l[0] = opposite_line(l[0]);
      }
    }

    if ( is_on_positive_halfspace(l[1], p) ||
	 is_on_positive_halfspace(l[1], r) ) {
      is_oriented[1] = true;
    } else {
      l[1] = opposite_line(l[1]);
      if ( is_on_positive_halfspace(l[1], p) ||
	   is_on_positive_halfspace(l[1], r) ) {
	is_oriented[1] = true;
      } else {
	l[1] = opposite_line(l[1]);
      }
    }

    if ( is_on_positive_halfspace(l[2], p) ||
	 is_on_positive_halfspace(l[2], q) ) {
      is_oriented[2] = true;
    } else {
      l[2] = opposite_line(l[2]);
      if ( is_on_positive_halfspace(l[2], p) ||
	   is_on_positive_halfspace(l[2], q) ) {
	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 );

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

    RT 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();
    }


    RT 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();
    }

    RT 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]);
    }

    RT vx = RT(0), vy = RT(0), vw = RT(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);

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

    Sign sgn_dist = Sign(s_vw * CGAL::sign(dist));

    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 Segment_2& p, const Segment_2& q, const Segment_2& r)
  {
    v_type = SSS;