approx_offset_base_2.h

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

        n_segments = 1;
      }
      else
      {
        abs_delta_x = (sign_delta_x == POSITIVE) ? delta_x : -delta_x;
        abs_delta_y = (sign_delta_y == POSITIVE) ? delta_y : -delta_y;

        // In this general case, the length d of the current edge is usually
        // an irrational number.
        // Compute the upper bound for the approximation error for d.
        // This bound is given by:
        //
        //                           |  (d - delta_y)  |
        //     bound = 2 * d * eps * | --------------- |
        //                           |     delta_x     |
        //
        const double    dd = CGAL::sqrt (CGAL::to_double (sqr_d));
        const double    derr_bound =
          2 * dd * _eps * ::fabs ((dd - CGAL::to_double (delta_y)) /
                                  CGAL::to_double (delta_x));

        err_bound = NT (derr_bound);
       
        // Compute an approximation for d (the squared root of sqr_d).
        int             denom = _inv_sqrt_eps;
        const int       max_int = (1 << (8*sizeof(int) - 2));

        while (static_cast<double>(max_int) / denom < dd)
        {
          denom /= 2;
        }

        app_d = NT (static_cast<int> (dd * denom + 0.5)) /
                NT (denom);
        app_err = sqr_d - CGAL::square (app_d);

        while (CGAL::compare (CGAL::abs (app_err), 
                              err_bound) == CGAL::LARGER ||
               CGAL::compare (app_d, abs_delta_x) != LARGER ||
               CGAL::compare (app_d, abs_delta_y) != LARGER)
        {
          app_d = (app_d + sqr_d/app_d) / 2;
          app_err = sqr_d - CGAL::square (app_d);
        }

        sign_app_err = CGAL::sign (app_err);

        if (sign_app_err == CGAL::ZERO)
        {
          // In this case d is a rational number, and we should shift the
          // both edge endpoints by (r * delta_y / d, -r * delta_x / d) to
          // obtain the offset points op1 and op2.
          const NT   trans_x = r * delta_y / app_d;
          const NT   trans_y = r * (-delta_x) / app_d;

          op1 = Point_2 (x1 + trans_x, y1 + trans_y);
          op2 = Point_2 (x2 + trans_x, y2 + trans_y);

          seg1 = X_monotone_curve_2 (op1, op2);
          dir_right1 = (sign_delta_x == CGAL::POSITIVE);

          n_segments = 1;
        }
        else
        {
          // Act according to the sign of delta_x.
          if (sign_delta_x == CGAL::NEGATIVE)
          {
            // x1 > x2, so we take a lower approximation for the squared root.
            if (sign_app_err == CGAL::NEGATIVE)
              app_d = sqr_d / app_d;
          }
          else
          {
            // x1 < x2, so we take an upper approximation for the squared root.
            if (sign_app_err == CGAL::POSITIVE)
              app_d = sqr_d / app_d;
          }

          // If theta is the angle that the vector (delta_x, delta_y) forms
          // with the x-axis, the perpendicular vector forms an angle of
          // phi = theta - PI/2, and we can approximate tan(phi/2) from below
          // and from above using:
          lower_tan_half_phi = (app_d - delta_y) / (-delta_x);
          upper_tan_half_phi = (-delta_x) / (app_d + delta_y);

          // Translate (x1, y1) by (r*cos(phi-), r*sin(phi-)) and create the
          // first offset point.
          // If tan(phi/2) = t is rational, then sin(phi) = 2t/(1 + t^2)
          // and cos(phi) = (1 - t^2)/(1 + t^2) are also rational.
          sqr_tan_half_phi = CGAL::square (lower_tan_half_phi);
          sin_phi = 2 * lower_tan_half_phi / (1 + sqr_tan_half_phi);
          cos_phi = (1 - sqr_tan_half_phi) / (1 + sqr_tan_half_phi);
          
          op1 = Point_2 (x1 + r*cos_phi, y1 + r*sin_phi);
          
          // Translate (x2, y2) by (r*cos(phi+), r*sin(phi+)) and create the
          // second offset point.
          sqr_tan_half_phi = CGAL::square (upper_tan_half_phi);
          sin_phi = 2 * upper_tan_half_phi / (1 + sqr_tan_half_phi);
          cos_phi = (1 - sqr_tan_half_phi) / (1 + sqr_tan_half_phi);
          
          op2 = Point_2 (x2 + r*cos_phi, y2 + r*sin_phi);
          
          // Compute the line l1 tangent to the circle centered at (x1, y1)
          // with radius r at the approximated point op1.
          l1 = f_perp_line (f_line (*curr, op1), op1);
        
          // Compute the line l2 tangent to the circle centered at (x2, y2)
          // with radius r at the approximated point op2.
          l2 = f_perp_line (f_line (*next, op2), op2);

          // Intersect the two lines. The intersection point serves as a common
          // end point for the two line segments we are about to introduce.
          obj = f_intersect (l1, l2);
        
          assign_success = CGAL::assign (mid_p, obj);
          CGAL_assertion (assign_success);

          // Create the two segments [op1 -> p_mid] and [p_min -> op2].
          seg1 = X_monotone_curve_2 (op1, mid_p);
          dir_right1 = (f_comp_xy (op1, mid_p) == CGAL::SMALLER);

          seg2 = X_monotone_curve_2 (mid_p, op2);
          dir_right2 = (f_comp_xy (mid_p, op2) == CGAL::SMALLER);

          n_segments = 2;
        }
      }

      if (curr == first)
      {
	// This is the first edge we visit -- store op1 for future use.
	first_op = op1;
      }
      else
      {
        // Connect prev_op and op1 with a circular arc, whose supporting circle
        // is (x1, x2) with radius r.
        arc = Curve_2 (*curr, r, CGAL::COUNTERCLOCKWISE,
                       Tr_point_2 (prev_op.x(), prev_op.y()),
                       Tr_point_2 (op1.x(), op1.y()));

        // Subdivide the arc into x-monotone subarcs and insert them to the
        // convolution cycle.
        xobjs.clear();
        f_make_x_monotone (arc, std::back_inserter(xobjs));

        for (xobj_it = xobjs.begin(); xobj_it != xobjs.end(); ++xobj_it)
        {
          assign_success = CGAL::assign (xarc, *xobj_it);
          CGAL_assertion (assign_success);

          *oi = Labeled_curve_2 (xarc,
                                 X_curve_label (xarc.is_directed_right(),
                                                cycle_id,
                                                curve_index));
          ++oi;
          curve_index++;
        }
      }

      // Append the offset segment(s) to the convolution cycle.
      CGAL_assertion (n_segments == 1 || n_segments == 2);

      *oi = Labeled_curve_2 (seg1,
                             X_curve_label (dir_right1,
                                            cycle_id,
                                            curve_index));
      oi++;
      curve_index++;

      if (n_segments == 2)
      {
        *oi = Labeled_curve_2 (seg2,
                               X_curve_label (dir_right2,
                                              cycle_id,
                                              curve_index));
        ++oi;
        curve_index++;
      }
    
      // Proceed to the next polygon vertex.
      prev_op = op2;
      curr = next;
    
    } while (curr != first);

    // Close the convolution cycle by creating the final circular arc,
    // centered at the first vertex.
    arc = Curve_2 (*first, r, CGAL::COUNTERCLOCKWISE,
                   Tr_point_2 (op2.x(), op2.y()),
                   Tr_point_2 (first_op.x(), first_op.y()));

    // Subdivide the arc into x-monotone subarcs and insert them to the
    // convolution cycle.
    bool           is_last;

    xobjs.clear();
    f_make_x_monotone (arc, std::back_inserter(xobjs));

    xobj_it = xobjs.begin();
    while (xobj_it != xobjs.end())
    {
      assign_success = CGAL::assign (xarc, *xobj_it);
      CGAL_assertion (assign_success);
      
      ++xobj_it;
      is_last = (xobj_it == xobjs.end());

      *oi = Labeled_curve_2 (xarc,
                             X_curve_label (xarc.is_directed_right(),
                                            cycle_id,
                                            curve_index,
                                            is_last));
      curve_index++;
    }

    return (oi);
  }

};


CGAL_END_NAMESPACE

#endif