conic_x_monotone_arc_2.h

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

                                 alg_t * sl_numer*sl_denom +
                                 alg_s * sl_denom*sl_denom;
    const Algebraic  sl2_denom = sl_denom*sl_denom*sl_denom;

    if (i == 2)
    {

      // Make sure that the denominator is always positive.
      if (CGAL::sign (sl_denom) != NEGATIVE)
      {
        slope_numer = -_two *sl2_numer;
        slope_denom = sl2_denom;
      }
      else
      {
        slope_numer = _two *sl2_numer;
        slope_denom = -sl2_denom;
      }

      return;
    }

    // The third-order derivative is given by:
    //
    //              (2t*alpha - t*beta) * gamma
    //   y''' = -6 ------------------------------
    //                    beta^2 * delta
    //
    const Algebraic  sl3_numer = (_two * alg_r * sl_numer -
                                  alg_t * sl_denom) * sl2_numer;
    const Algebraic  sl3_denom = sl_denom*sl_denom * sl2_denom;

    if (i == 3)
    {
      // Make sure that the denominator is always positive.
      if (CGAL::sign (sl_denom) != NEGATIVE)
      {
        slope_numer = -6 * sl3_numer;
        slope_denom = sl3_denom;
      }
      else
      {
        slope_numer = 6 * sl3_numer;
        slope_denom = -sl2_denom;
      }

      return;
    }

    // \todo Handle higher-order derivatives as well.
    CGAL_assertion (false);
    return;
  }

  /*!
   * Compute the overlap with a given arc, which is supposed to have the same
   * supporting conic curve as this arc.
   * \param arc The given arc.
   * \param overlap Output: The overlapping arc (if any).
   * \return Whether we found an overlap.
   */
  bool _compute_overlap (const Self& arc, Self& overlap) const
  {
    // Check if the two arcs are identical.
    if (equals (arc))
    {
      overlap = arc;
      return (true);
    }

    if (_is_strictly_between_endpoints (arc.left()))
    {
      if (_is_strictly_between_endpoints(arc.right()))
      {
        // Case 1 - *this:     +----------->
        //            arc:       +=====>
        overlap = arc;
        return (true);
      }
      else
      {
        // Case 2 - *this:     +----------->
        //            arc:               +=====>
        overlap = *this;

        if ((overlap._info & IS_DIRECTED_RIGHT) != 0)
          overlap._source = arc.left();
        else
          overlap._target = arc.left();

        return (true);
      }
    }
    else if (_is_strictly_between_endpoints (arc.right()))
    {
      // Case 3 - *this:     +----------->
      //            arc:   +=====>
      overlap = *this;

      if ((overlap._info & IS_DIRECTED_RIGHT) != 0)
        overlap._target = arc.right();
      else
        overlap._source = arc.right();

      return (true);
    }
    else if (arc._is_between_endpoints (this->_source) &&
             arc._is_between_endpoints (this->_target) &&
             (arc._is_strictly_between_endpoints(this->_source) ||
              arc._is_strictly_between_endpoints(this->_target)))
    {
      // Case 4 - *this:     +----------->
      //            arc:   +================>
      overlap = *this;
      return (true);
    }

    // If we reached here, there are no overlaps:
    return (false);
  }

  /*!
   * Intersect the supporing conic curves of this arc and the given arc.
   * \param arc The arc to intersect with.
   * \param inter_list The list of intersection points.
   */
  void _intersect_supporting_conics (const Self& arc,
                                     Intersection_list& inter_list) const
  {
    if (_is_special_segment() && ! arc._is_special_segment())
    {
      // If one of the arcs is a special segment, make sure it is (arc).
      arc._intersect_supporting_conics (*this, inter_list);
      return;
    }

    const int   deg1 = ((this->_info & DEGREE_MASK) == DEGREE_1) ? 1 : 2;
    const int   deg2 = ((arc._info & DEGREE_MASK) == DEGREE_1) ? 1 : 2;
    Nt_traits   nt_traits;
    Algebraic   xs[4];
    int         n_xs = 0;
    Algebraic   ys[4];
    int         n_ys = 0;

    if (arc._is_special_segment())
    {
      // The second arc is a special segment (a*x + b*y + c = 0).
      if (_is_special_segment())
      {
        // Both arc are sepcial segment, so they have at most one intersection
        // point.
        Algebraic   denom = this->_extra_data_P->a * arc._extra_data_P->b -
                            this->_extra_data_P->b * arc._extra_data_P->a;

        if (CGAL::sign (denom) != CGAL::ZERO)
        {
          xs[0] = (this->_extra_data_P->b * arc._extra_data_P->c -
                   this->_extra_data_P->c * arc._extra_data_P->b) / denom;
          n_xs = 1;

          ys[0] = (this->_extra_data_P->c * arc._extra_data_P->a -
                   this->_extra_data_P->a * arc._extra_data_P->c) / denom;
          n_ys = 1;
        }
      }
      else
      {
        // Compute the x-coordinates of the intersection points.
        n_xs = _compute_resultant_roots (nt_traits,
                                         alg_r, alg_s, alg_t,
                                         alg_u, alg_v, alg_w,
                                         deg1,
                                         arc._extra_data_P->a,
                                         arc._extra_data_P->b,
                                         arc._extra_data_P->c,
                                         xs);
        CGAL_assertion (n_xs <= 2);
      
        // Compute the y-coordinates of the intersection points.
        n_ys = _compute_resultant_roots (nt_traits,
                                         alg_s, alg_r, alg_t,
                                         alg_v, alg_u, alg_w,
                                         deg1,
                                         arc._extra_data_P->b,
                                         arc._extra_data_P->a,
                                         arc._extra_data_P->c,
                                         ys);
        CGAL_assertion (n_ys <= 2);
      }
    }
    else
    {
      // Compute the x-coordinates of the intersection points.
      n_xs = _compute_resultant_roots (nt_traits,
                                       this->_r, this->_s, this->_t,
                                       this->_u, this->_v, this->_w,
                                       deg1,
                                       arc._r, arc._s, arc._t,
                                       arc._u, arc._v, arc._w,
                                       deg2,
                                       xs);
      CGAL_assertion (n_xs <= 4);
      
      // Compute the y-coordinates of the intersection points.
      n_ys = _compute_resultant_roots (nt_traits,
                                       this->_s, this->_r, this->_t,
                                       this->_v, this->_u, this->_w,
                                       deg1,
                                       arc._s, arc._r, arc._t,
                                       arc._v, arc._u, arc._w,
                                       deg2,
                                       ys);
      CGAL_assertion (n_ys <= 4);
    }

    // Pair the coordinates of the intersection points. As the vectors of
    // x and y-coordinates are sorted in ascending order, we output the
    // intersection points in lexicographically ascending order.
    unsigned int  mult;
    int           i, j;

    if (arc._is_special_segment())
    {
      if (n_xs == 0 || n_ys == 0)
        return;

      if (n_xs == 1 && n_ys == 1)
      {
        // Single intersection.
        Conic_point_2         ip (xs[0], ys[0]);

        ip.set_generating_conic (_id);
        ip.set_generating_conic (arc._id);

        // In case the other curve is of degree 2, this is a tangency point.
        mult = (deg1 == 1 || _is_special_segment()) ? 1 : 2;
        inter_list.push_back (Intersection_point_2 (ip, mult));
      }
      else if (n_xs == 1 && n_ys == 2)
      {
        Conic_point_2         ip1 (xs[0], ys[0]);

        ip1.set_generating_conic (_id);
        ip1.set_generating_conic (arc._id);

        inter_list.push_back (Intersection_point_2 (ip1, 1));

        Conic_point_2         ip2 (xs[0], ys[1]);

        ip2.set_generating_conic (_id);
        ip2.set_generating_conic (arc._id);

        inter_list.push_back (Intersection_point_2 (ip2, 1));
      }
      else if (n_xs == 2 && n_ys == 1)
      {
        Conic_point_2         ip1 (xs[0], ys[0]);

        ip1.set_generating_conic (_id);
        ip1.set_generating_conic (arc._id);

        inter_list.push_back (Intersection_point_2 (ip1, 1));

        Conic_point_2         ip2 (xs[1], ys[0]);

        ip2.set_generating_conic (_id);
        ip2.set_generating_conic (arc._id);

        inter_list.push_back (Intersection_point_2 (ip2, 1));

      }
      else
      {
        CGAL_assertion (n_xs == 2 && n_ys == 2);

        // The x-coordinates and the y-coordinates are given in ascending
        // order. If the slope of the segment is positive, we pair the
        // coordinates as is - otherwise, we swap the pairs.
        int                   ind_first_y = 0, ind_second_y = 1;

        if (CGAL::sign (arc._extra_data_P->b) == 
            CGAL::sign(arc._extra_data_P->a))
        {
          ind_first_y = 1;
          ind_second_y = 0;
        }

        Conic_point_2         ip1 (xs[0], ys[ind_first_y]);

        ip1.set_generating_conic (_id);
        ip1.set_generating_conic (arc._id);

        inter_list.push_back (Intersection_point_2 (ip1, 1));

        Conic_point_2         ip2 (xs[1], ys[ind_second_y]);

        ip2.set_generating_conic (_id);
        ip2.set_generating_conic (arc._id);

        inter_list.push_back (Intersection_point_2 (ip2, 1));
      }
      
      return;
    }

    for (i = 0; i < n_xs; i++)
    {
      for (j = 0; j < n_ys; j++)
      {
        if (_is_on_supporting_conic (xs[i], ys[j]) &&
            arc._is_on_supporting_conic (xs[i], ys[j]))
        {
          // Create the intersection point and set its generating conics.
          Conic_point_2         ip (xs[i], ys[j]);

          ip.set_generating_conic (_id);
          ip.set_generating_conic (arc._id);

          // Compute the multiplicity of the intersection point.
          if (deg1 == 1 && deg2 == 1)
            mult = 1;
          else
            mult = _multiplicity_of_intersection_point (arc, ip);

          // Insert the intersection point to the output list.
          inter_list.push_back (Intersection_point_2 (ip, mult));
        }
      }
    }

    return;
  }

  /*!
   * Compute the multiplicity of an intersection point.
   * \param arc The arc to intersect with.
   * \param p The intersection point.
   * \return The multiplicity of the intersection point.
   */
  unsigned int _multiplicity_of_intersection_point (const Self& arc,
                                                    const Point_2& p) const
  {
    CGAL_assertion (! _is_special_segment() || ! arc._is_special_segment());

    // Compare the slopes of the two arcs at p, using their first-order
    // partial derivatives.
    Algebraic      slope1_numer, slope1_denom;
    Algebraic      slope2_numer, slope2_denom;

    _derive_by_x_at (p, 1, slope1_numer, slope1_denom);
    arc._derive_by_x_at (p, 1, slope2_numer, slope2_denom);

    if (CGAL::compare (slope1_numer*slope2_denom,
                       slope2_numer*slope1_denom) != EQUAL)
    {
      // Different slopes at p - the mutiplicity of p is 1:
      return (1);
    }

    if (CGAL::sign (slope1_denom) != ZERO &&
        CGAL::sign (slope2_denom) != ZERO)
    {
      // The curves do not have a vertical slope at p.
      // Compare their second-order derivative by x:
      _derive_by_x_at (p, 2, slope1_numer, slope1_denom);
      arc._derive_by_x_at (p, 2, slope2_numer, slope2_denom);
    }
    else
    {
      // Both curves have a vertical slope at p.
      // Compare their second-order derivative by y:
      _derive_by_y_at (p, 2, slope1_numer, slope1_denom);
      arc._derive_by_y_at (p, 2, slope2_numer, slope2_denom);
    }

    if (CGAL::compare (slope1_numer*slope2_denom,
                       slope2_numer*slope1_denom) != EQUAL)
    {
      // Different curvatures at p - the mutiplicity of p is 2:
      return (2);
    }

    // If we reached here, the multiplicity of the intersection point is 3:
    return (3);
  }
  //@}

};

/*!
 * Exporter for x-monotone conic arcs.
 */
template <class Conic_arc_2>
std::ostream& operator<< (std::ostream& os, 
                          const _Conic_x_monotone_arc_2<Conic_arc_2>& arc)
{
  // Output the supporting conic curve.
  os << "{" << CGAL::to_double(arc.r()) << "*x^2 + "
     << CGAL::to_double(arc.s()) << "*y^2 + "
     << CGAL::to_double(arc.t()) << "*xy + " 
     << CGAL::to_double(arc.u()) << "*x + "
     << CGAL::to_double(arc.v()) << "*y + "
     << CGAL::to_double(arc.w()) << "}";

  // Output the endpoints.
  os << " : (" << CGAL::to_double(arc.source().x()) << "," 
     << CGAL::to_double(arc.source().y()) << ") ";

  if (arc.orientation() == CLOCKWISE)
    os << "--cw-->";
  else if (arc.orientation() == COUNTERCLOCKWISE)
    os << "--ccw-->";
  else
    os << "--l-->";

  os << " (" << CGAL::to_double(arc.target().x()) << "," 
     << CGAL::to_double(arc.target().y()) << ")";

  return (os);
}

CGAL_END_NAMESPACE

#endif