ボロノイ図とその応用

ちょっとした歴史

//   // kInfPoint -  , kNegInfPoint -   Point2D Line2D::GetIntersection(const Line2D& second_line) const { double cross_prod_norms = Vector2D(this->A, this->B).OrientedCCW(Vector2D(second_line.A, second_line.B)); Point2D intersect_point; if (fabs(cross_prod_norms) <= EPS) /* A1 / A2 == B1 / B2 */ { if (fabs(this->B * second_line.C - second_line.B * this->C) <= EPS) /* .. == C1 / C2 */ { intersect_point = kNegInfPoint2D; } else { intersect_point = kInfPoint2D; } } else { double res_x = (second_line.C * this->B - this->C * second_line.B) / cross_prod_norms; double res_y = (second_line.A * this->C - this->A * second_line.C) / cross_prod_norms; intersect_point = Point2D(res_x, res_y); } return intersect_point; } //   Point2D Segment2D::GetIntersection(const Segment2D& second_seg) const { Line2D first_line(*this); Line2D second_line(second_seg); Point2D intersect_point = first_line.GetIntersection(second_line); if (intersect_point == kNegInfPoint2D) { if (this->Contains(second_seg.b)) { intersect_point = second_seg.b; } else if (this->Contains(second_seg.a)) { intersect_point = second_seg.a; } else if (second_seg.Contains(this->b)) { intersect_point = this->b; } else if (second_seg.Contains(this->a)) { intersect_point = this->a; } else { intersect_point = kInfPoint2D; } } else if (!(this->Contains(intersect_point) && second_seg.Contains(intersect_point))) { intersect_point = kInfPoint2D; } return intersect_point; } 

1. 最大反復回数が経過するまで（2 *（n + m）以下であることが証明されます）、次の手順を実行します。
   a）。 1番目と2番目のポリゴンの現在のエッジを取得します。
   b） それらが交差する場合-交差点からの点を考慮します：交差点ポリゴンに追加しただけで、追加を無視する場合があります。そうでない場合-最初の交差点でない場合（つまり、すでに円を作っている）、交差点、それ以外の場合はアルゴリズムを終了します-交差点に到達しました。
   c）。 次に（現在のエッジが交差するかどうかに関係なく）、 モーション関数を呼び出します。この関数は、最初の（または2番目の）ポリゴンのウィンドウを1エッジ前方に移動し、交差ポリゴンに頂点を追加できるようにします。 主な行動は、 運動で正確に行われます。
2. 交差点がなかった場合、あるポリゴンが別のポリゴンの内側にあるかどうかを判断します（ポイントがOの後ろの凸ポリゴンに属しているかどうかを確認します（log（ポリゴンの頂点の数））- [1] 、59-60ページ）。 そうであれば、それを返し、そうでなければ空の交差点を返します。

 //   ,       Convex2D Rectangle::GetIntersectionalConvex2D(const Point2D& cur_point, const Line2D& halfplane) const { vector<Point2D> convex_points; Segment2D cur_side; Point2D intersection_point; for (int i = 0, sz = vertices_.size(); i < sz; ++i) { int j = (i + 1) % sz; cur_side = Segment2D(vertices_[i], vertices_[j]); intersection_point = halfplane.GetIntersection(cur_side); if (intersection_point != kInfPoint2D) convex_points.push_back(intersection_point); if (halfplane.Sign(cur_point) == halfplane.Sign(vertices_[i])) convex_points.push_back(vertices_[i]); } Convex2D result_polygon(MakeConvexHullJarvis(convex_points)); return result_polygon; } //       NumOfCase EdgesCaseNum(const Segment2D& first_edge, const Segment2D& second_edge) { bool first_looks_at_second = first_edge.LooksAt(second_edge); bool second_looks_at_first = second_edge.LooksAt(first_edge); if (first_looks_at_second && second_looks_at_first) { return NumOfCase::kBothLooks; } else if (first_looks_at_second) { return NumOfCase::kFirstLooksAtSecond; } else if (second_looks_at_first) { return NumOfCase::kSecondLooksAtFirst; } else { return NumOfCase::kBothNotLooks; } } // ,    "" WhichEdge WhichEdgeIsInside(const Segment2D& first_edge, const Segment2D& second_edge) { double first_second_side = Vector2D(second_edge).OrientedCCW(Vector2D(second_edge.a, first_edge.b)); double second_first_side = Vector2D(first_edge).OrientedCCW(Vector2D(first_edge.a, second_edge.b)); if (first_second_side < 0) { return WhichEdge::kSecondEdge; } else if (second_first_side < 0) { return WhichEdge::kFirstEdge; } else { return WhichEdge::Unknown; } } //   WhichEdge MoveOneOfEdges(const Segment2D& first_edge, const Segment2D& second_edge, Convex2D& result_polygon) { WhichEdge now_inside = WhichEdgeIsInside(first_edge, second_edge); NumOfCase case_num = EdgesCaseNum(first_edge, second_edge); WhichEdge which_edge_is_moving; switch (case_num) { case NumOfCase::kBothLooks: { if (now_inside == WhichEdge::kFirstEdge) { which_edge_is_moving = WhichEdge::kSecondEdge; } else { which_edge_is_moving = WhichEdge::kFirstEdge; } break; } case NumOfCase::kFirstLooksAtSecond: { which_edge_is_moving = WhichEdge::kFirstEdge; break; } case NumOfCase::kSecondLooksAtFirst: { which_edge_is_moving = WhichEdge::kSecondEdge; break; } case NumOfCase::kBothNotLooks: { if (now_inside == WhichEdge::kFirstEdge) { which_edge_is_moving = WhichEdge::kSecondEdge; } else { which_edge_is_moving = WhichEdge::kFirstEdge; } break; } } if (result_polygon.Size() != 0 && (case_num == NumOfCase::kFirstLooksAtSecond || case_num == NumOfCase::kSecondLooksAtFirst)) { Point2D vertex_to_add; if (now_inside == WhichEdge::kFirstEdge) { vertex_to_add = first_edge.b; } else if (now_inside == WhichEdge::kSecondEdge) { vertex_to_add = second_edge.b; } else { if (case_num == NumOfCase::kFirstLooksAtSecond) vertex_to_add = first_edge.b; // ?! else vertex_to_add = second_edge.b; } if (vertex_to_add != result_polygon.GetCurVertex()) result_polygon.AddVertex(vertex_to_add); } return which_edge_is_moving; } //    (    ,      ) Convex2D Convex2D::GetIntersectionalConvex(Convex2D& second_polygon) { Convex2D result_polygon; size_t max_iter = 2 * (this->Size() + second_polygon.Size()); Segment2D cur_fp_edge; // current first polygon edge Segment2D cur_sp_edge; // current second polygon edge Point2D intersection_point; bool no_intersection = true; WhichEdge moving_edge = WhichEdge::Unknown; for (size_t i = 0; i < max_iter; ++i) { cur_fp_edge = this->GetCurEdge(); cur_sp_edge = second_polygon.GetCurEdge(); intersection_point = cur_fp_edge.GetIntersection(cur_sp_edge); if (intersection_point != kInfPoint2D) { if (result_polygon.Size() == 0) { no_intersection = false; result_polygon.AddVertex(intersection_point); } else if (intersection_point != result_polygon.GetCurVertex()) { if (intersection_point == result_polygon.vertices_[0]) { break; // we already found the intersection polygon } else { result_polygon.AddVertex(intersection_point); } } } moving_edge = MoveOneOfEdges(cur_fp_edge, cur_sp_edge, result_polygon); if (moving_edge == WhichEdge::kFirstEdge) { this->MoveCurVertex(); } else { second_polygon.MoveCurVertex(); } } if (no_intersection == true) { if (second_polygon.Contains(this->GetCurVertex())) { result_polygon = *this; } else if (this->Contains(second_polygon.GetCurVertex())) { result_polygon = second_polygon; } } return result_polygon; } 

 Convex2D GetHalfPlanesIntersection(const Point2D& cur_point, const vector<Line2D>& halfplanes, const Rectangle& border_box) { if (halfplanes.size() == 1) { Convex2D cur_convex(border_box.GetIntersectionalConvex2D(cur_point, halfplanes[0])); return cur_convex; } else { int middle = halfplanes.size() >> 1; vector<Line2D> first_half(halfplanes.begin(), halfplanes.begin() + middle); vector<Line2D> second_half(halfplanes.begin() + middle, halfplanes.end()); Convex2D first_convex(GetHalfPlanesIntersection(cur_point, first_half, border_box)); Convex2D second_convex(GetHalfPlanesIntersection(cur_point, second_half, border_box)); return first_convex.GetIntersectionalConvex(second_convex); } } 

 //      Voronoi2DLocus VoronoiDiagram2D::MakeVoronoi2DLocus(const Point2D& site, const vector<Point2D>& points, const Rectangle& border_box) { Voronoi2DLocus cur_locus; vector<Line2D> halfplanes; for (auto cur_point : points) { if (cur_point != site) { Segment2D cur_seg(site, cur_point); Line2D cur_halfplane(cur_seg.GetCenter(), cur_seg.NormalVec()); halfplanes.push_back(cur_halfplane); } } *cur_locus.region_ = GetHalfPlanesIntersection(site, halfplanes, border_box); cur_locus.site_ = site; return cur_locus; } //    VoronoiDiagram2D VoronoiDiagram2D::MakeVoronoiDiagram2DHalfPlanes(const vector<Point2D>& points, const Rectangle& border_box) { Voronoi2DLocus cur_locus; for (auto cur_point : points) { cur_locus = MakeVoronoi2DLocus(cur_point, points, border_box); this->diagram_.push_back(cur_locus); } return *this; } 

 VoronoiDiagram2D VoronoiDiagram2D::MakeVoronoiDiagram2DFortune(const vector<Point2D>& points, const Rectangle& border_box) { priority_queue<Event> events_queue(points.begin(), points.end()); shared_ptr<Event> cur_event; BeachSearchTree beach_line; DCEL edges; while (!events_queue.empty()) { cur_event = make_shared<Event>(events_queue.top()); shared_ptr<const PointEvent> is_point_event(dynamic_cast<const PointEvent *>(cur_event.get())); if (is_point_event) { events_queue.pop(); beach_line.HandlePointEvent(*is_point_event, border_box, events_queue, edges); } else { shared_ptr<const CircleEvent> is_circle_event(dynamic_cast<const CircleEvent *>(cur_event.get())); events_queue.pop(); beach_line.HandleCircleEvent(*is_circle_event, border_box, events_queue, edges); } } edges.Finish(border_box); this->dcel_ = edges; return *this; } 

 pair<Point2D, Point2D> Arch::GetIntersection(const Arch& second_arch, double line_pos) const { pair<Point2D, Point2D> intersect_points; double p1 = 2 * (line_pos - this->focus_->y); double p2 = 2 * (line_pos - second_arch.focus_->y); // is not 0.0, because line moved down if (fabs(p1) <= EPS) { intersect_points.first = this->GetIntersection(Ray2D(*this->focus_, Point2D(this->focus_->x, this->focus_->y + 1))); intersect_points.second = intersect_points.first; } else { // solving the equation double a1 = 1 / p1; double a2 = 1 / p2; double a = a2 - a1; double b1 = -this->focus_->x / p1; double b2 = -second_arch.focus_->x / p2; double b = b2 - b1; double c1 = pow(this->focus_->x, 2) + pow(this->focus_->y, 2) - pow(line_pos, 2) / p1; double c2 = pow(second_arch.focus_->x, 2) + pow(second_arch.focus_->y, 2) - pow(line_pos, 2) / p1; double c = c2 - c1; double D = pow(b, 2) - 4 * a * c; if (D < 0) { intersect_points = make_pair(kInfPoint2D, kInfPoint2D); } else if (fabs(D) <= EPS) { double x = -b / (2 * a); double y = a1 * pow(x, 2) + b1 * x + c1; intersect_points = make_pair(Point2D(x, y), Point2D(x, y)); } else { double x1 = (-b - sqrt(D)) / (2 * a); double x2 = (-b + sqrt(D)) / (2 * a); double y1 = a1 * pow(x1, 2) + b1 * x1 + c1; double y2 = a1 * pow(x2, 2) + b1 * x2 + c1; intersect_points = make_pair(Point2D(x1, y1), Point2D(x2, y2)); } } return intersect_points; } 

 Circle::Circle(const Point2D& p1, const Point2D& p2, const Point2D& p3) { Segment2D first_segment(p1, p2); Segment2D second_segment(p2, p3); Line2D first_perpendicular(first_segment.GetCenter(), first_segment.NormalVec()); Line2D second_perpendicular(second_segment.GetCenter(), second_segment.NormalVec()); center_ = first_perpendicular.GetIntersection(second_perpendicular); little_haxis_ = big_haxis_ = center_.l2_distance(p1); }