circle_segment_2.h

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

   */
  Comparison_result _circs_compare_to_left (const Self& cv, 
                                            const Point_2& p) const
  {
    if (_index() != 0 && _index() == cv._index())
    {
      // Check the case of comparing two circular arcs that originate from the
      // same supporting circle. Their comparison result is not EQUAL only if
      // one is an upper arc and the other is a lower arc.
      if (_is_upper() && ! cv._is_upper())
        return (LARGER);
      else if (! _is_upper() && cv._is_upper())
        return (SMALLER);
      else
        return (EQUAL);
    }

    // We have to compare the slopes of the two supporting circles at p:
    //
    //    p.x() - x0(1)         p.x() - x0(2)
    //   ---------------  and  ---------------
    //    y0(1) - p.y()         y0(2) - p.y()
    //
    // Eventually, we should take the opposite result.
    const CGAL::Sign  sign_numer1 = CGAL::sign (p.x() - x0());
    const CGAL::Sign  sign_denom1 = CGAL::sign (y0() - p.y());
    const CGAL::Sign  sign_numer2 = CGAL::sign (p.x() - cv.x0());
    const CGAL::Sign  sign_denom2 = CGAL::sign (cv.y0() - p.y());

    // Check the case of vertical tangents.
    if (sign_denom1 == ZERO)
    {
      if (sign_denom2 == ZERO)
      {
        if (_is_upper())
        {
          if (cv._is_upper())
          {
            // The two circles have a vertical tangent:
            // The one with a larger radius is above the other.
            return (CGAL::compare (sqr_r(), cv.sqr_r()));
          }
          else
          {
            // The other curve is directed downwards:
            return (LARGER);
          }
        }
        else
        {
          if (cv._is_upper())
          {
            // The other curve is directed upwards:
            return (SMALLER);
          }
          else
          {
            // The two circles have a vertical tangent:
            // The one with a smaller radius is above the other.
            return (CGAL::compare (cv.sqr_r(), sqr_r()));
          }
        }
      }

      // The other arc does not have a vertical tangent.
      return (_is_upper() ? LARGER : SMALLER);
    }
    else if (sign_denom2 == ZERO)
    {
      return (cv._is_upper() ? SMALLER : LARGER);
    }

    // Try to act according to the slope signs.
    CGAL::Sign   sign_slope1;
    CGAL::Sign   sign_slope2;

    if (sign_numer1 == sign_denom1)
      sign_slope1 = POSITIVE;
    else if (sign_numer1 == ZERO)
      sign_slope1 = ZERO;
    else
      sign_slope1 = NEGATIVE;

    if (sign_numer2 == sign_denom2)
      sign_slope2 = POSITIVE;
    else if (sign_numer2 == ZERO)
      sign_slope2 = ZERO;
    else
      sign_slope2 = NEGATIVE;

    if ((sign_slope1 == POSITIVE && sign_slope2 != POSITIVE) ||
        (sign_slope1 == ZERO && sign_slope2 == NEGATIVE))
      return (SMALLER);

    if ((sign_slope2 == POSITIVE && sign_slope1 != POSITIVE) ||
        (sign_slope2 == ZERO && sign_slope1 == NEGATIVE))
      return (LARGER);

    // Compare the slopes of the two tangents to the circles.
    Comparison_result  slope_res;
    
    if (sign_slope1 == ZERO && sign_slope2 == ZERO)
    {
      // Special case were both circles have a horizontal tangent:
      slope_res = EQUAL;
    }
    else
    {
      // Actually compare the slopes.    
      const bool    swap_res = (sign_denom1 != sign_denom2);
      const CoordNT A = (cv.y0() - y0())*p.x() + (y0()*cv.x0() - cv.y0()*x0());
      const CoordNT B = (cv.x0() - x0())*p.y();
     
      slope_res = CGAL::compare (A, B);

      if (slope_res != EQUAL && swap_res)
      {
        // Swap the comparison result, if necessary:
        slope_res = (slope_res == SMALLER) ? LARGER : SMALLER;
      }
    }

    // In case the two circles have different tangent slopes at p, return
    // the opposite of the slope result (since the slope result is the
    // comparison result to the right of the intersection point):
    if (slope_res != EQUAL)
      return ((slope_res == SMALLER) ? LARGER : SMALLER);

    // In this case we have a tangency point at p.
    if (_is_upper())
    {
      if (cv._is_upper())
      {
        // The circle with a larger radius is above the other.
        return (CGAL::compare (sqr_r(), cv.sqr_r()));
      }
      else
      {
        // The other curve is above our curve:
        return (SMALLER);
      }
    }
    else
    {
      if (cv._is_upper())
      {
        // Out curve is above the other curve:
        return (LARGER);
      }
      else
      {
        // The circle with a smaller radius is above the other.
        return (CGAL::compare (cv.sqr_r(), sqr_r()));
      }
    }
  }
  //@}

  /// \name Auxiliary functions for computing intersections.
  //@{

  /*!
   * Compute the intersections between two line segments.
   */
  void _lines_intersect (const Self& cv,
                         Intersection_list& inter_list) const
  {
    // The intersection of the lines:
    //   a1*x + b1*y + c1 = 0   and   a2*x + b2*y + c2 = 0 ,
    // is given by:
    //
    //      b1*c2 - c1*b2     c1*a2 - a1*c2
    //   ( --------------- , --------------- )
    //      a1*b2 - b1*a2     a1*b2 - b1*a2
    //
    unsigned int  mult = 1;
    const NT      denom = a()*cv.b() - b()*cv.a();

    // Make sure the supporting lines are not parallel.
    if (CGAL::sign(denom) == ZERO)
      return;

    const NT      x_numer = b()*cv.c() - c()*cv.b();
    const NT      y_numer = c()*cv.a() - a()*cv.c();
    Point_2       p (x_numer / denom, y_numer / denom);

    inter_list.push_back (Intersection_point_2 (p, mult));
    return;
  }

  /*!
   * Compute the intersections between the supporting circle of (*this) and 
   * the supporting line of the segement cv.
   */
  void _circ_line_intersect (const Self& cv,
                             Intersection_list& inter_list) const
  {
    Point_2       p;
    unsigned int  mult;

    // First check the special cases of vertical and horizontal lines.
    if (cv.is_vertical())
    {
      // The equation of the vertical line is x = -c / a.
      // The y-coordinates of the intersection points are:
      //   y =  y0 +/- sqrt(r^2 - (x - x0)^2)
      //
      const NT   vx = -cv.c() / cv.a();
      const NT   vdisc = sqr_r() - CGAL::square (vx - x0());
      CGAL::Sign sign_vdisc = CGAL::sign (vdisc);
      
      if (sign_vdisc == NEGATIVE)
      {
        // The circle and the vertical line do not intersect.
        return;
      }
      else if (sign_vdisc == ZERO)
      {
        // A single tangency point, given by:
        mult = 2;
        p = Point_2 (vx, y0());
        inter_list.push_back (Intersection_point_2 (p, mult));
        
        return;
      }
      
      // Compute the two intersection points:
      mult = 1;

      p = Point_2 (CoordNT (vx),
                   CoordNT (y0(), -1, vdisc));
      inter_list.push_back (Intersection_point_2 (p, mult));
      
      p = Point_2 (CoordNT (vx),
                   CoordNT (y0(), 1, vdisc));
      inter_list.push_back (Intersection_point_2 (p, mult));

      return;
    }
    else if (CGAL::sign (cv.a()) == ZERO)
    {
      // The equation of the horizontal line is y = -c / b.
      // The y-coordinates of the intersection points are:
      //   x =  x0 +/- sqrt(r^2 - (y - y0)^2)
      //
      const NT   hy = -cv.c() / cv.b();
      const NT   hdisc = sqr_r() - CGAL::square (hy - y0());
      CGAL::Sign sign_hdisc = CGAL::sign (hdisc);
      
      if (sign_hdisc == NEGATIVE)
      {
        // The circle and the vertical line do not intersect.
        return;
      }
      else if (sign_hdisc == ZERO)
      {
        // A single tangency point, given by:
        mult = 2;
        p = Point_2 (x0(), hy);
        inter_list.push_back (Intersection_point_2 (p, mult));

        return;
      }
      
      // Compute the two intersection points:
      mult = 1;

      p = Point_2 (CoordNT (x0(), -1, hdisc),
                   CoordNT (hy));
      inter_list.push_back (Intersection_point_2 (p, mult));

      p = Point_2 (CoordNT (x0(), 1, hdisc),
                   CoordNT (hy));
      inter_list.push_back (Intersection_point_2 (p, mult));

      return;
    }

    // Compute the squared distance between the line and the circle center,
    // inducing the discriminant of the quadratic equations we have to solve.
    const NT   line_factor = CGAL::square(cv.a()) + CGAL::square(cv.b());
    const NT   disc = line_factor*sqr_r() -
                      CGAL::square(cv.a()*x0() + cv.b()*y0() + cv.c());
    CGAL::Sign sign_disc = CGAL::sign (disc);

    if (sign_disc == NEGATIVE)
    {
      // The circle and the line do not intersect:
      return;
    }

    // Compare the square-free part of the solution:
    const NT   aux = cv.b()*x0() - cv.a()*y0();
    const NT   x_base = (aux*cv.b() - cv.a()*cv.c()) / line_factor;
    const NT   y_base = (-aux*cv.a() - cv.b()*cv.c()) / line_factor;

    if (sign_disc == ZERO)
    {
      // A single tangency point, given by:
      mult = 2;
      p = Point_2 (x_base, y_base);
      inter_list.push_back (Intersection_point_2 (p, mult));

      return;
    }

    // We have two intersection points, whose coordinates are one-root numbers.
    bool       minus_root_first = (CGAL::sign(cv.b()) == POSITIVE);
    const NT   x_root_coeff = cv.b() / line_factor;
    const NT   y_root_coeff = cv.a() / line_factor;

    mult = 1;
    if (minus_root_first)
    {
      p = Point_2 (CoordNT (x_base, -x_root_coeff, disc),
                   CoordNT (y_base, y_root_coeff, disc));
      inter_list.push_back (Intersection_point_2 (p, mult));

      p = Point_2 (CoordNT (x_base, x_root_coeff, disc),
                   CoordNT (y_base, -y_root_coeff, disc));
      inter_list.push_back (Intersection_point_2 (p, mult));
    }
    else
    {
      p = Point_2 (CoordNT (x_base, x_root_coeff, disc),
                   CoordNT (y_base, -y_root_coeff, disc));
      inter_list.push_back (Intersection_point_2 (p, mult));

      p = Point_2 (CoordNT (x_base, -x_root_coeff, disc),
                   CoordNT (y_base, y_root_coeff, disc));
      inter_list.push_back (Intersection_point_2 (p, mult));
    }

    return;
  }

  /*!
   * Compute the intersections between two circles.
   */
  void _circs_intersect (const Self& cv,
                         Intersection_list& inter_list) const
  {
    Point_2       p;
    unsigned int  mult;

    // Compute the squared distance between the circle centers, inducing the
    // discriminant of the quadratic equations we have to solve.
    const NT   diff_x = cv.x0() - x0();
    const NT   diff_y = cv.y0() - y0();    
    const NT   sqr_dist = CGAL::square(diff_x) + CGAL::square(diff_y);
    const NT   diff_sqr_rad = sqr_r() - cv.sqr_r();
    const NT   disc = 2*sqr_dist*(sqr_r() + cv.sqr_r()) -
                      (CGAL::square(diff_sqr_rad) + CGAL::square(sqr_dist));
    CGAL::Sign sign_disc = CGAL::sign (disc);

    if (sign_disc == NEGATIVE)
    {
      // The two circles do not intersect.
      return;
    }

    // Compare the square-free part of the solution:
    const NT   x_base = ((x0() + cv.x0()) + diff_x*diff_sqr_rad/sqr_dist) / 2;
    const NT   y_base = ((y0() + cv.y0()) + diff_y*diff_sqr_rad/sqr_dist) / 2;

    if (sign_disc == ZERO)
    {
      // A single tangency point, given by:
      mult = 2;
      p = Point_2 (x_base, y_base);
      inter_list.push_back (Intersection_point_2 (p, mult));

      return;
    }

    // We have two intersection points, whose coordinates are one-root numbers.
    CGAL::Sign sign_diff_y = CGAL::sign (diff_y);
    bool       minus_root_first;

    if (sign_diff_y == ZERO)
      minus_root_first = (CGAL::sign (diff_x) == NEGATIVE);
    else
      minus_root_first = (sign_diff_y == POSITIVE);

    const NT   x_root_coeff = diff_y / (2 * sqr_dist);
    const NT   y_root_coeff = diff_x / (2 * sqr_dist);

    mult = 1;
    if (minus_root_first)
    {
      p = Point_2 (CoordNT (x_base, -x_root_coeff, disc),
                   CoordNT (y_base, y_root_coeff, disc));
      inter_list.push_back (Intersection_point_2 (p, mult));

      p = Point_2 (CoordNT (x_base, x_root_coeff, disc),
                   CoordNT (y_base, -y_root_coeff, disc));
      inter_list.push_back (Intersection_point_2 (p, mult));
    }
    else
    {
      p = Point_2 (CoordNT (x_base, x_root_coeff, disc),
                   CoordNT (y_base, -y_root_coeff, disc));
      inter_list.push_back (Intersection_point_2 (p, mult));

      p = Point_2 (CoordNT (x_base, -x_root_coeff, disc),
                   CoordNT (y_base, y_root_coeff, disc));
      inter_list.push_back (Intersection_point_2 (p, mult));
    }

    return;
  }

  /*!
   * Check if the given point lies on the arc.
   * \pre p lies on the supporting curve.
   */
  bool _is_between_endpoints (const Point_2& p) const
  {
    if (is_linear())
    {
      if (is_vertical())
      {
        // Check if the point is in the y-range of the arc.
        // Note that left() is the lower endpoint and right() is the upper
        // endpoint of the segment in this case. 
        Comparison_result    res = CGAL::compare (p.y(), left().y());

        if (res == SMALLER)
          return (false);
        else if (res == EQUAL)
          return (true);

        return (CGAL::compare (p.y(), right().y()) != LARGER);
      }

      // For non-vertical segments, it is sufficient to check if the point
      // is in the x-range of the arc.
      return (this->is_in_x_range (p));
    }

    // The supporting curve is a circle:
    // Check whether p lies on the upper or on the lower part of the circle.
    Comparison_result   c_res = CGAL::compare (p.y(), y0());

    if ((_is_upper() && c_res == SMALLER) ||
        (! _is_upper() && c_res == LARGER))
    {
      // The point lies on the other half of the circle:
      return (false);
    }

    // Check if the point is in the x-range of the arc.
    return (this->is_in_x_range (p));
  }

  /*!
   * Check if the given point lies in the interior of the arc.
   * \pre p lies on the supporting cu