voronoi_vertex_ring_c2.h

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

    Sqrt_1 Zero(RT(0), RT(0), d[0]);
    Sqrt_1 sqrt_D0(RT(0), RT(1), d[0]);

    Sqrt_1 D1 = d[1] + Zero;
    Sqrt_1 D2 = d[2] + Zero;

    Sqrt_3 vz(z[0] * sqrt_D0, z[1] + Zero, z[2] + Zero, Zero, D1, D2);

    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 = Sqrt_3(z[0] * sqrt_D0, z[1] + Zero, z[2] + Zero, Zero, D1, D2);

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

    Sqrt_3 vx, vy, vw;

    vx = Sqrt_3(x[0] * sqrt_D0, x[1] + Zero, x[2] + Zero, Zero, D1, D2);
    vy = Sqrt_3(y[0] * sqrt_D0, y[1] + Zero, y[2] + Zero, Zero, D1, D2);
    vw = Sqrt_3(w[0] * sqrt_D0, w[1] + Zero, w[2] + Zero, Zero, D1, D2);

    Sqrt_1 a1(a[(i_no+1)%3], RT(0), d[0]);
    Sqrt_1 b1(b[(i_no+1)%3], RT(0), d[0]);
    Sqrt_1 c1(c[(i_no+1)%3], RT(0), d[0]);

    Sqrt_3 dist = a1 * vx + b1 * vy + c1 * vw;

    Sign s_vw = CGAL::sign(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 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;

    RT a[3], b[3], c[3];
    RT cx[3], cy[3], cz[3], 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]);
      D[i] = CGAL::square(a[i]) + CGAL::square(b[i]);

      Sqrt_1 Zero(RT(0), RT(0), D[0]);
      Sqrt_1 sqrt_D0(RT(0), RT(1), D[0]);
      Sqrt_1 D1 = Zero + D[1];
      Sqrt_1 D2 = Zero + D[2];

      ux = Sqrt_3(cx[0] * sqrt_D0, cx[1] + Zero, cx[2] + Zero,
		  Zero, D1, D2);
      uy = Sqrt_3(cy[0] * sqrt_D0, cy[1] + Zero, cy[2] + Zero,
		  Zero, D1, D2);
      uz = Sqrt_3(cz[0] * sqrt_D0, cz[1] + Zero, cz[2] + Zero,
		  Zero, D1, D2);      
    }
  }

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

  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);

    } else if ( s1.is_segment() && s2.is_point() && s3.is_point() ) {
      compute_vertex(s2, s3, s1);
      pps_idx = 1;

    } else if ( s1.is_point() && s2.is_segment() && s3.is_point() ) {
      compute_vertex(s3, s1, s2);
      pps_idx = 2;

    } else if ( s1.is_point() && s2.is_point() && s3.is_segment() ) {
      compute_pps(s1, s2, s3);
      pps_idx = 0;

    } else if ( s1.is_point() && s2.is_segment() && s3.is_segment() ) {
      compute_pss(s1, s2, s3);
    } else if ( s1.is_segment() && s2.is_point() && s3.is_segment() ) {
      compute_vertex(s2, s3, s1);
    } else if ( s1.is_segment() && s2.is_segment() && s3.is_point() ) {
      compute_vertex(s3, s1, s2);
    } else {
      compute_sss(s1, s2, s3);
    }

  }

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

  bool is_endpoint_of(const Site_2& p, const Site_2& s) const
  {
    CGAL_precondition( p.is_point() && s.is_segment() );
    return ( same_points(p, s.source_site()) ||
	     same_points(p, s.target_site()) );
  }
  

  //--------------------------------------------------------------------------
  //--------------------------------------------------------------------------
  //                           the orientation test
  //--------------------------------------------------------------------------
  //--------------------------------------------------------------------------

  Orientation
  orientation(const Line_2& l, PPP_Type) const
  {
    Sign s_uz = CGAL::sign(uz_ppp);
    Sign s_l =
      CGAL::sign(l.a() * ux_ppp + l.b() * uy_ppp + l.c() * uz_ppp);

    Sign s = Sign(s_uz * s_l);

    if ( s == ZERO ) { return COLLINEAR; }
    return ( s == POSITIVE ) ? LEFT_TURN : RIGHT_TURN;
  }

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

  Orientation
  orientation(const Line_2& l, PPS_Type) const
  {
    Sign s_uz = CGAL::sign(uz_pps);
    Sign s_l =
      CGAL::sign(l.a() * ux_pps + l.b() * uy_pps + l.c() * uz_pps);

    Sign s = Sign(s_uz * s_l);

    if ( s == ZERO ) { return COLLINEAR; }
    return ( s == POSITIVE ) ? LEFT_TURN : RIGHT_TURN;
  }


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

  // the cases PSS and SSS are identical
  template<class Type>
  Orientation
  orientation(const Line_2& l, Type) const
  {
    Sqrt_1 Zero(RT(0), RT(0), ux.a().c());

    Sqrt_1 a = l.a() + Zero;
    Sqrt_1 b = l.b() + Zero;
    Sqrt_1 c = l.c() + Zero;

    Sign s_uz = CGAL::sign(uz);
    Sign s_l = CGAL::sign(a * ux + b * uy + c * uz);

    Sign s = Sign(s_uz * s_l);

    if ( s == ZERO ) { return COLLINEAR; }
    return ( s == POSITIVE ) ? LEFT_TURN : RIGHT_TURN;
  }
    

  //--------------------------------------------------------------------------
  //--------------------------------------------------------------------------
  //                              the incircle test
  //--------------------------------------------------------------------------
  //--------------------------------------------------------------------------

  //--------------------------------------------------------------------------
  //  the incircle test when the fourth site is a point
  //--------------------------------------------------------------------------

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

  Sign check_easy_degeneracies(const Site_2& t, PPS_Type,
			       bool& use_result) const
  {
    CGAL_precondition( t.is_point() );

    use_result = false;
    if (  ( p_.is_point() && same_points(p_, t) ) ||
	  ( q_.is_point() && same_points(q_, t) ) ||
	  ( r_.is_point() && same_points(r_, t) )  ) {
      use_result = true;
      return ZERO;
    }

    if (  ( p_.is_segment() && is_endpoint_of(t, p_) ) || 
	  ( q_.is_segment() && is_endpoint_of(t, q_) ) ||
	  ( r_.is_segment() && is_endpoint_of(t, r_) )  ) {
      use_result = true;
      return POSITIVE;
    }

    return ZERO;
  }

  inline
  Sign check_easy_degeneracies(const Site_2& t, PSS_Type,
			       bool& use_result) const
  {
    CGAL_precondition( t.is_point() );
    return check_easy_degeneracies(t, PPS_Type(), use_result);
  }

  inline
  Sign check_easy_degeneracies(const Site_2& t, SSS_Type,
			       bool& use_result) const
  {
    CGAL_precondition( t.is_point() );
    use_result = false;
    // ADD THE CASES WHERE t IS AN ENDPOINT OF ONE OF THE SEGMENTS
    return ZERO;
  }

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

  template<class Type>
  inline  
  Sign incircle_p(const Site_2& st, Type type) const
  {
    CGAL_precondition( st.is_point() );

    bool use_result(false);
    Sign s = check_easy_degeneracies(st, type, use_result);
    if ( use_result ) { return s; }

    return incircle_p_no_easy(st, type);
  }

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

  Sign incircle_p(const Site_2& st, PPP_Type) const
  {
    CGAL_precondition( st.is_point() );

    Point_2 t = st.point();

    Oriented_side os =
      side_of_oriented_circle(p_.point(), q_.point(), r_.point(), t);
    if ( os == ON_POSITIVE_SIDE ) { return NEGATIVE; }
    if ( os == ON_NEGATIVE_SIDE ) { return POSITIVE; }
    return ZERO;
  }

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

  Sign incircle_p_no_easy(const Site_2& st, PPS_Type ) const
  {
    CGAL_precondition( st.is_point() );

    Point_2 t = st.point();

    Point_2 pref = p_ref().point();

    Sqrt_1 vx = ux_pps - pref.x() * uz_pps;
    Sqrt_1 vy = uy_pps - pref.y() * uz_pps;

    Sqrt_1 Rs = CGAL::square(vx) + CGAL::square(vy);

    Sqrt_1 Rs1 = CGAL::square(ux_pps - t.x() * uz_pps)
      + CGAL::square(uy_pps - t.y() * uz_pps);

    return CGAL::sign(Rs1 - Rs);
  }

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

  Sign incircle_p_no_easy(const Site_2& st, PSS_Type ) const
  {
    CGAL_precondition( st.is_point() );
    Point_2 t = st.point();

    Sqrt_1 Zero(RT(0), RT(0), ux.a().c());

    Point_2 pref = p_ref().point();

    Sqrt_1 xref = pref.x() + Zero;
    Sqrt_1 yref = pref.y() + Zero;

    Sqrt_3 vx = ux - xref * uz;
    Sqrt_3 vy = uy - yref * uz;

    Sqrt_3 Rs = CGAL::square(vx) + CGAL::square(vy);

    Sqrt_1 tx = t.x() + Zero;
    Sqrt_1 ty = t.y() + Zero;

    Sqrt_3 Rs1 =
      CGAL::square(ux - tx * uz) + CGAL::square(uy - ty * uz);

    Sign s_Q = CGAL::sign(Rs1 - Rs);

    return s_Q;
  }

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

  Sign incircle_p_no_easy(const Site_2& st, SSS_Type ) const
  {
    CGAL_precondition( st.is_point() );

    Point_2 t = st.point();

    Sqrt_1 Zero(RT(0), RT(0), ux.a().c());

    RT a1, b1, c1;
    compute_supporting_line(p_.supporting_site(), a1, b1, c1);

    RT ns = CGAL::square(a1) + CGAL::square(b1);

    Sqrt_1 a = a1 + Zero;
    Sqrt_1 b = b1 + Zero;
    Sqrt_1 c = c1 + Zero;

    Sqrt_1 Ns = ns + Zero;
    Sqrt_3 Ls = CGAL::square(a * ux + b * uy + c * uz);

    Sqrt_1 tx = t.x() + Zero;
    Sqrt_1 ty = t.y() + Zero;

    Sqrt_3 R1s = CGAL::square(ux - tx * uz)
      + CGAL::square(uy - ty * uz);

    return CGAL::sign(R1s * Ns - Ls);
  }



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

  Sign incircle_p(const Site_2& t) const 
  {
    if ( is_degenerate_Voronoi_circle() ) {
      return POSITIVE;
    }

    Sign s(ZERO);
    switch ( v_type ) {
    case PPP:
      s = incircle_p(t, PPP_Type());
      break;
    case PPS:
      s = incircle_p(t, PPS_Type());
      break;
    case PSS:
      s = incircle_p(t, PSS_Type());
      break;
    case SSS:
      s = incircle_p(t, SSS_Type());
      break;
    }

    return s;
  }

  Sign incircle_p_no_easy(const Site_2& t) const 
  {
    Sign s(ZERO);
    switch ( v_type ) {
    case PPP:
      s = incircle_p(t, PPP_Type());
      break;
    case PPS:
      s = incircle_p_no_easy(t, PPS_Type());
      break;
    case PSS:
      s = incircle_p_no_easy(t, PSS_Type());
      break;
    case SSS:
      s = incircle_p_no_easy(t, SSS_Type());
      break;
    }

    return s;
  }

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

  //--------------------------------------------------------------------------
  //  the incircle test when the fourth site is a segment
  //--------------------------------------------------------------------------

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

  Oriented_side
  oriented_side(const Line_2& l, const Point_2& p, PPP_Type) const
  {
    Sign s_uz = CGAL::sign(uz_ppp);

    RT px = uz_ppp * p.x() - ux_ppp;
    RT py = uz_ppp * p.y() - uy_ppp;

    Sign s1 = CGAL::sign(l.b() * px - l.a() * py);

    Sign s = Sign(s_uz * s1);

    if ( s == POSITIVE ) { return ON_POSITIVE_SIDE; }
    if ( s == NEGATIVE ) { return ON_NEGATIVE_SIDE; }
    return ON_ORIENTED_BOUNDARY;
  }

  Oriented_side
  oriented_side(const Line_2& l, const Point_2& p, PPS_Type) const
  {
    Sqrt_1 dx = ux_pps - uz_pps * p.x();
    Sqrt_1 dy = uy_pps - uz_pps * p.y();

    Sign s = Sign(CGAL::sign(uz_pps) *
		  CGAL::sign(dy * l.a() - dx * l.b())
		  );

    if ( s == POSITIVE ) { return ON_POSITIVE_SIDE; }
    if ( s == NEGATIVE ) { return ON_NEGATIVE_SIDE; }
    return ON_ORIENTED_BOUNDARY;
  }

  // the cases PSS and SSS are identical
  template<class Type>
  Oriented_side
  oriented_side(const Line_2& l, const Point_2& p, Type) const
  {
    Sqrt_1 Zero(RT(0), RT(0), ux.a().c());
    Sqrt_1 px = p.x() + Zero;
    Sqrt_1 py = p.y() + Zero;

    Sqrt_3 dx = ux - px * uz;
    Sqrt_3 dy = uy - py * uz;

    Sqrt_1 a = l.a() + Zero;
    Sqrt_1 b = l.b() + Zero;

    Sign s = Sign(CGAL::sign(uz) *
		  CGAL::sign(a * dy - b * dx)
		  );


    if ( s == POSITIVE ) { return ON_POSITIVE_SIDE; }
    if ( s == NEGATIVE ) { return ON_NEGATIVE_SIDE; }
    return ON_ORIENTED_BOUNDARY;
  }

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

  
  Sign incircle(const Line_2& l, PPP_Type) const
  {
    Point_2 pref = p_ref().point();

    RT a1 = CGAL::square(l.a()) + CGAL::square(l.b());
    RT a2 = CGAL::square(ux_ppp - pref.x() * uz_ppp) +
      CGAL::square(uy_ppp - pref.y() * uz_ppp);

    RT a3 =
      CGAL::square(l.a() * ux_ppp + l.b() * uy_ppp + l.c() * uz_ppp);

    Comparison_result cr = CGAL::compare(a3, a1 * a2);

    if ( cr == LARGER ) { return POSITIVE; }