Диаграмма Вороного и её применения

Доброго всем времени суток, уважаемые посетители сайта Хабрахабр. В данной статье я бы хотел рассказать вам о том, что такое диаграмма Вороного (изображена на картинке ниже), о различных алгоритмах её построения (за , — пересечение полуплоскостей, — алгоритм Форчуна) и некоторых тонкостях реализации (на языке C++).

Также будет рассмотрено много интересных применений диаграммы и несколько любопытных фактов о ней. Будет интересно!

План статьи:

— Необходимые понятия и определения
— Алгоритмы построения
— Применения
— Интересные факты
— Список литературы

Стоит отметить, что в данной статье будут рассматриваться только алгоритмы построения диаграммы Вороного на плоскости. Попутно будут рассмотрены некоторые другие алгоритмы, необходимые для построения диаграммы — алгоритм определения точки пересечения двух отрезков, алгоритм О`Рурка пересечения двух выпуклых многоугольников, алгоритм построения «береговой линии».

Необходимые понятия и определения

Сразу скажу, что всё, что будет дальше происходить, находится на плоскости.

Что ж, перед тем, как начать разбираться, что это такое — диаграмма Вороного, напомню некоторые понятия нужных нам геометрических объектов (предполагается, однако, что вы уже знакомы с определениями точки, прямой, луча, отрезка, многоугольника, вершины и ребра многоугольника, вектора и интуитивного понятия разбиения плоскости):

Простой многоугольник — это многоугольник без самопересечений. Мы будем рассматривать только простые многоугольники.

Невыпуклый многоугольник — это многоугольник, в котором найдутся такие две вершины, что через них проводится прямая, пересекающая данный многоугольник где-либо ещё, кроме ребра, соединяющего эти вершины (см. картинку),

Выпуклый многоугольник — это многоугольник, у которого продолжения сторон не пересекают других его сторон (см. картинку). С другими вариантами определений можно ознакомиться на википедии.

Именно из выпуклых многоугольников и будет состоять диаграмма. Почему именно из выпуклых? Потому что они являются ничем иным, как пересечением полуплоскостей (как мы увидим чуть позже), которые являются выпуклыми фигурами, а вот почему пересечение выпуклых фигур — выпуклая фигура, предложу вам узнать самим (доказательство есть, например, в книге [2]).

Раз уж я начал говорить про полуплоскости, то можно плавно перейти и к самой диаграмме — она состоит из так называемых локусов — областей, в которых присутствуют все точки, которые находятся ближе к данной точке, чем ко всем остальным. В диаграмме Вороного локусы являются выпуклыми многоугольниками.

Как строить локус? По определению он будет строиться так: пусть дано множество из n точек, для которого мы строим диаграмму. Возьмём конкретную точку p, для которой строим локус, и ещё одну точку из данного нам множества — q (не равную p). Проведём отрезок, соединяющий эти две точки, и проведём прямую, которая будет являться серединным перпендикуляром данного отрезка. Эта прямая делит плоскость на две полуплоскости — в одной лежит точка p, в другой лежит точка q. В данном случае локусами этих двух точек являются полученные полуплоскости. То есть для того, чтобы построить локус точки p, нужно получить пересечение всех таких полуплоскостей (то есть на месте q побывают все точки данного множества, кроме p).

Точку, для которой строится локус, называют сайтом (site). На следующей картинке локусы помечены разными цветами.

Алгоритмы построения диаграммы как раз и есть не что иное, как алгоритмы построения этих самых локусов для всех точек из заданного набора. Локусы в данной задаче также называют многоугольниками Вороного или ячейками Вороного.

Наконец, сформулируем определение диаграммы Вороного n точек на плоскости (n — натуральное) — это разбиение плоскости, состоящее из n локусов (для каждой точки по локусу). Опять же, другой вариант определения можно найти на википедии.

Кстати, вот сайт с интерактивным цветным визуализатором диаграммы Вороного, можно самому нажимать на область и залипать видеть, как строится диаграмма.

Если интересно, как появилась диаграмма и почему она носит фамилию Вороного, то вам стоит взглянуть под спойлер ниже.

Алгоритмы построения

Научимся строить диаграмму Вороного. Мы рассмотрим 4 алгоритма, 2 из которых подробно (1 с реализацией, полной реализации алгоритма Форчуна будет посвящена отдельная статья), ещё 2 кратко (и без реализации):

1. Алгоритм построения диаграммы Вороного «в лоб». Сложность: ;
2. Алгоритм построения диаграммы Вороного путём пересечения полуплоскостей. Сложность: ;
3. Алгоритм Форчуна построения диаграммы Вороного на плоскости. Сложность: ;
4. Рекурсивный алгоритм построения диаграммы Вороного. Сложность: .

После описания некоторых алгоритмов будет приведена его реализация на языке C++. Для реализации использовалась написанная нами библиотека SplashGeom © — ссылка на github, в которой есть всё необходимое для реализации многих алгоритмов вычислительной геометрии на плоскости и некоторых в пространстве. Прошу не судить строго данную библиотеку, она ещё находится в стадии активной разработки и улучшения, однако все замечания будут услышаны.

Если же вам интересны другие реализации, то вот ещё некоторые:

• boost.polygon
• Сайт Стива Форчуна
• aewllin/openvoronoi
• Реализация на java
• Реализация на JavaScript (там ещё и визуализатор есть)

Теперь перейдём к непосредственному рассмотрению алгоритмов:

Алгоритм построения диаграммы Вороного «в лоб» за

Здесь идея в том, чтобы пересекать не полуплоскости, а именно серединные перпендикуляры отрезков (потому что это проще, согласитесь), соединяющих данную точку со всеми другими точками. То есть мы, следуя определению ячейки Вороного, будем строить локус для точки p так:

1. Получаем n-1 прямую (серединные перпендикуляры), так как мы провели серединные перпендикуляры всех отрезков, соединяющих данную точку p с остальными;
2. Пересекаем попарно все прямые, получаем точек пересечения (потому что каждая прямая может пересечь все другие, в «худшем случае»);
3. Проверяем все эти точек на принадлежность каждой из n-1 полуплоскостей, то есть получаем уже асимптотику . Соответственно те точки, которые принадлежат всем полуплоскостям, и будут вершинами ячейки Вороного точки p;
4. Проделываем первые три шага для всех n точек, получаем итоговую асимптотику .

С алгоритмом также можно ознакомиться на e-maxx.ru.

При желании Вы самостоятельно можете реализовать этот алгоритм с помощью SplashGeom ©. В данной статье его реализация не приводится, потому как данный алгоритм на практике сильно уступает хотя бы следующему…

Алгоритм построения диаграммы Вороного путём пересечения полуплоскостей за

Этот алгоритм уже можно использовать на практике, так как он не имеет столь большой вычислительной сложности. Для него нам потребуется: уметь пересекать отрезки и прямые, уметь пересекать выпуклые многоугольники, уметь пересекать полуплоскости, уметь объединять полученные локусы в диаграмму.

Алгоритм

1. Получаем n-1 прямую для текущего сайта (как в предыдущем алгоритме — серединные перпендикуляры). Это будут «образующие» полуплоскостей;
2. Теперь мы имеем n-1 полуплоскость. Каждая из этих полуплоскостей задаётся какой-либо прямой из пред. пункта и ориентацией, то есть с какой стороны от прямой она расположена. Ориентацию можно определить по текущему сайту, для которого строим локус — он лежит в искомой полуплоскости, а значит и его локус должен лежать в ней;
3. Пересекаем все полуплоскости — мы сможем делать это за — получаем локус для текущего сайта;
4. Проделываем первые три шага для всех n точек, получаем итоговую асимптотику * = .

Реализация

Основная загвоздка здесь, я считаю, реализовать нормальное пересечение выпуклых многоугольников, потому как там есть неприятные вырожденные случаи (совпадение вершин и/или сторон многоугольников; в случае пересечения многоугольника с самим собой надо неплохо подумать, чтобы адаптировать алгоритм О`Рурка для корректной работы в этом случае).

Сразу скажу, что мы ограничиваем всю область действий ограничивающим прямоугольником — полуплоскости будут отсекать от него части, которые мы и пересечём друг с другом. Это решает проблему бесконечных ячеек в диаграмме.

Пересечение прямых и отрезков

Причём нам нужны именно точки пересечения, а не просто определение его наличия. В SplashGeom © это есть — вот, к примеру, реализация пересечения прямых и отрезков:

// пересечение прямых
// 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;
}

Пересечение выпуклых многоугольников

Теперь реализуем пересечение многоугольников. Можно, конечно, делать это за , где n и m — количество вершин первого и второго многоугольника соответственно — пересекать каждую сторону первого многоугольника с каждой стороной из второго, записывая точки пересечения, а также проверять точки на принадлежность другому многоугольнику, но мы, дабы достичь лучшей скорости, будем пересекать по алгоритму О`Рурка (оригинальное описание алгоритма — [3]).

Также с одним из вариантов его описания можно ознакомиться на algolist.ru. У нас же реализация основана на описании алгоритма в книге [1] (стр. 334), с некоторыми дополняющими идеями. Сразу отмечу, что данная реализация не учитывает случаи, когда у многоугольников есть общие стороны (как отмечено в книге [1], этот случай требует отдельного рассмотрения), однако случаи с общими вершинами работают корректно.

Под спойлером ниже можно увидеть общее описание алгоритма.

Ну а в SplashGeom © это выглядит так:

// для получения многоугольников, которые полуплоскости отсекают от ограничивающего прямоугольника
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;
}

Так, мы научились строить диаграмму Вороного за , она имеет вид вектора (списка) локусов. Недостаток данного решения — информации о соседях не получить (возможно, этот недостаток можно ликвидировать улучшенной реализацией).

Куда больше информации можно получить из РСДС (DCEL). Эта структура будет использована в алгоритме Форчуна.

Далее будет описан алгоритм Форчуна с использованием заметающей прямой и «береговой линии» (готовьте купальники и плавки — идём на пляж). По моему мнению, это самый приемлемый вариант, если вы хотите реализовать построение диаграммы Вороного и строить её со «скоростью» света .

Алгоритм Форчуна построения диаграммы Вороного за

В 1987 году Стив Форчун (Steve Fortune) предложил алгоритм построения диаграммы за . Конечно, он является не единственным алгоритмом построения с такой асимптотикой, но он достаточно понятен и не очень сложен в реализации (да и, к тому же, очень красив и математичен!), поэтому я выбрал именно его.

Материалы по алгоритму Форчуна можно найти тут, здесь, вот здесь и ещё вот здесь.

Кстати, рассмотрению данного алгоритма уже была посвящена статья на Хабрахабре.

Итак, основная идея алгоритма — это так называемая заметающая прямая (ЗП) (sweep line). Она применяется во многих алгоритмах вычислительной геометрии, потому что позволяет удобно моделировать движение прямой по некоторому множеству объектов (например, в алгоритме пересечения n отрезков тоже используется sweep line).

Перед тем, как начать говорить про то, как и что мы будем делать, давайте посмотрим, как движется sweep line (взято отсюда):

Красиво, не правда ли? В реализации всё примерно так же, только ЗП обычно движется сверху вниз, а не слева направо, и на самом деле всё не так плавно, а происходит от события к событию (см. ниже), то есть дискретно.

Суть алгоритма

Есть n сайтов (точек на плоскости). Есть заметающая прямая, которая двигается (например) «сверху вниз», то есть от сайта с наибольшей ординатой к сайту с меньшей (от события к событию, если быть точным). Сразу стоит отметить, что влияние на построение диаграммы оказывают только те сайты, которые находятся выше или на заметающей прямой.

Когда ЗП попадает на очередной сайт (происходит событие точки (point event)), создаётся новая парабола (arch), фокусом которой является данный сайт, а директрисой — заметающая прямая (про параболу на википедии). Эта парабола делит плоскость на две части — «внутренняя» область параболы соответствует точкам, которые сейчас ближе к сайту, а «внешняя» область — точкам, которые ближе к sweep line, ну а точки, лежащие на параболе — равноудалены от сайта и ЗП. Парабола будет меняться в зависимости от положения ЗП к сайту — чем дальше ЗП уходит от сайта вниз, тем больше расширяется парабола, однако в самом начале она вообще является отрезком («направленным» вверх).

Далее парабола расширяется, у неё появляются две контрольные точки (break points) — точки её пересечения с остальными параболами («береговой линией»). В «береговой линии» мы храним дуги парабол от одной точки пересечения их друг с другом до другой, так и получается beach line. По сути, в этом алгоритме мы моделируем движение этой «береговой линии». потому как эти самые break point`ы движутся аккурат по рёбрам ячеек Вороного (ведь получается, что контрольные точки равноудалены от обоих сайтов, которым соответствуют эти параболы, да ещё и от ЗП).

И как раз-таки в тот момент, когда две контрольные точки — по одной из разных парабол — «встречаются», то есть как бы превращаются в одну, эта точка и становится вершиной ячейки Вороного (происходит событие круга (circle event)), причём в это время та дуга, которая находилась между этими двумя точками — «схлопывается» и удаляется из «береговой линии». Далее мы просто соединяем эту точку с предыдущей соответствующей ей и получаем ребро ячейки Вороного.

Алгоритм

Итак, при движении sweep line вниз мы встречаемся с двумя типами событий:

Событие точки (point event)

Событие точки — это попадание ЗП на один из сайтов, поэтому мы создаём новую параболу, соответствующую данному сайту, а также добавляются две контрольные точки (break points) (на самом деле сначала — одна, а при расширении арки уже две) — точки пересечения этой параболы с береговой линией (то есть с фронтом уже существующих парабол). Стоит отметить, что в данном алгоритме парабола (а точнее её часть, принадлежащая «береговой лини» — арка) «вставляется в береговую линию» только в случае события точки, то есть новая арка может появиться только при обработке события точки.

Кстати, на следующих картинках видно, почему объединения таких «кусков парабол» называют «береговой линией».

Событие круга (circle event)

Событие круга — это возникновение новой вершины ячейки Вороного вместе с удалением одной арки, потому что возникновение новой вершины диаграммы здесь означает, что было три арки, левая, средняя и правая, средняя «схлопывается» вследствие сближения левой и правой точек арок и получается новая вершина диаграммы Вороного. Стоит отметить, что в данном алгоритме парабола (арка) удаляется из «береговой линии» только в случае события круга, то есть арка может удалиться только при обработке события круга.

Есть теорема, в которой говорится, что вершина диаграммы Вороного всегда лежит на пересечении ровно трёх рёбер диаграммы, и есть следствие из этой теоремы, которое гласит, что вершина диаграммы является центром окружности, проходящей через три сайта и расстояние от этой точки до заметающей прямой тоже равно радиусу этой окружности (это свойство точек, лежащих на «береговой линии»). Это — ключевой момент, потому что именно когда самая нижняя точка окружности, проходящей через три сайта, лежит ниже или на заметающей прямой, мы пушим в очередь событий событие круга с этой самой нижней точкой, потому что когда прямая на неё попадёт, мы получим вершину диаграммы Вороного.

Важно, что с любым событием (точки или круга) связана одна конкретная арка, и наоборот. Это пригодится при обработке событий. Также надо не забыть, что нужно вовремя добавлять рёбра в РСДС (DCEL) (пункт 1 в структурах, см. ниже), так что надо понимать связь арок с рёбрами.

Таким образом, движение прямой дискретно — прямая в любой момент времени либо на сайте, либо в нижней точке окружности, проходящей через три сайта, центр которой — новая вершина диаграммы Вороного. Прекрасно.

Общий алгоритм:

1. Создаём очередь (с приоритетом) событий, изначально инициализируя событиями точки — данным множеством сайтов (ведь каждому сайту соответствует событие точки);
2. Пока очередь не пуста:

   а). Берём из неё событие;
   б). Если это — событие точки, то обрабатываем событие точки;
   в). Если это — событие круга, то обрабатываем событие круга;

3. Закончить все оставшиеся рёбра (поработать с border_box).

Реализация

Реализация алгоритма Форчуна будет рассмотрена подробно в отдельной статье, однако здесь приводятся некоторые наработки, которые могут помочь в его понимании.

Необходимые структуры

Для реализации данного алгоритма нам понадобятся несколько структур (классов), а именно:

• РСДС (DCEL) — список для хранения уже найденный рёбер диаграммы Вороного;
• Приоритетная очередь с событиями;
• Двоичное дерево (у нас это BeachSearchTree) — для хранения «береговой линии» — текущего положения парабол и точек. Стоит отметить, что это дерево является сбалансированным — у узла либо ровно два сына, либо ноль (является листом). Подробнее об этой структуре можно прочитать, например, в этой статье на Хабре (у нас она представлена немного по-другому).

Имея такие структуры данных, можем написать реализацию общего алгоритма:

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

Вся логика и сложность в HandlePointEvent() и HandleCircleEvent(), им и будет посвящена отдельная статья, далее же приведу некоторые вспомогательные функции, которые потом помогут в реализации.

Вспомогательные функции

Пересечение парабол (арок)

Нам нужно уметь получать пересечение двух парабол (арок) в зависимости от положения ЗП. Уравнение параболы с фокусом в точке x' и y' и директрисой, положение которой по оси y равно l, задаётся следующим уравнением (его можно вывести):

Кстати, дробь перед скобкой — это фокальный параметр данной параболы. Отсюда мы можем «вытащить» соответствующие коэффициенты уравнения параболы и решить систему из двух нелинейных уравнений простым способом — вычтем одно из другого, а потом подставим найденные корни в первое, получим две точки. Нас интересует та точка, которая ниже (то есть с меньшей ординатой), потому что высокая точка лежит за «береговой линией». Описанные действия отражаются в следующем коде:

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

Обработка события точки

Событие точки — это когда мы вытащили из очереди PointEvent. Что в нём есть? В нём есть только сайт, с которым ассоциируется это событие.

Что мы делаем при его обработке? Мы добавляем новую арку в «береговую линию», настраивая в дереве все связи «как надо», и проверяем, не появилось ли событие круга в одном из трёх возможных случаев — нам нужно проверить все случаи, в которых может участвовать новый сайт.

Что мы делаем при добавлении арки? Мы ищем для неё место в нашем двоичном дереве «береговой линии» (по координате x), затем вставляем её.

Что мы делаем при вставке? Мы нашли указатель на арку, с которой новая имеет пересечение в двух точках (случай с попаданием точки пересечения ровно на одну из контрольных точек рассматривается отдельно — там пересечение будет с двумя параболами — слева и справа).

То есть мы берём эту арку, которая «разбивается», и вместо неё вставляем аж 5 узлов (был 1, стало 5, да) — arch1, bp1, arch2, bp2, arch3. Arch1 — это левый кусок арки, которую пересекает новая, то есть это кусок слева от левого break point`а, bp1 — левый break point (левое пересечение новой арки), arch2 — это новая арка собственной персоной, bp2 — правый break point (правое пересечение новой арки), arch3 — это правый кусок арки, которую пересекает новая.

Стоит отметить, что событие точки даёт начало новому ребру (или рёбрам, есть случаи) диаграммы Вороного.

Один из возможных вариантов добавления новой арки в «береговую линию» представлен в этой статье на Хабре, наш же вариант описан выше, и больше похож на тот, который предложен здесь.

Обработка события круга

Событие круга — это когда мы вытащили из очереди CircleEvent.

Что в нём есть? В нём есть точка — это самая нижняя точка некоторой окружности, которая проходит через три каких-то сайта, и арка, которую следует удалить. В дереве есть две её контрольные точки и она сама, контрольные точки в итоге превратятся в одну, а арку нужно аккуратно удалить из дерева, перестроив все «связи родителей-детей». По сути, обработка этого события заменяет в дереве три узла на один (два break point`а и арку на один break point).

Стоит отметить, что событие круга завершает два ребра диаграммы Вороного, то есть при обработке этого события рёбра будут завершаться.

Также скажу, что важной частью обработки событий является то, как мы следим за ростом рёбер, добавлением их в список и завершением рёбер, то есть за тем, чтобы в итоге они все были «закончены», и либо были конечны, либо упирались в ограничивающий прямоугольник.

Небольшой анализ результата

Так как мы получим на выходе уже РСДС (DCEL), то сможем получить достаточно информации — для каждого ребра мы знаем соответствующий сайт, можем получить «близнеца» этого ребра, узнать сайт для него и вуаля — мы узнали соседа (так, пройдя по всему списку, можем сформировать списки соседей для всех сайтов, а это уже достижение, так как всё произошло в итоге за + = ).

Причём если у ребра «нулевой» сайт, то мы оказались на «граничном сайте» — сайте с «бесконечным» локусом. Хмм, но ведь это позволит нам построить выпуклую оболочку начального множества точек за в итоге, здорово.

Ну и вообще, РСДС — это, по сути, граф, так что можно продолжать работать с этим списком во многих алгоритмах на графах.

Рекурсивный алгоритм построения диаграммы Вороного за

Данный алгоритм приводится в книге [1] (стр. 260), здесь я приведу только сам алгоритм построения, поскольку данный вариант мы не реализовывали, хотя он является хорошей аналогией алгоритму Форчуна.

Алгоритм

1. Делим всё множество сайтов S на две примерно равные части (может быть нечётное количество точек) S₁ и S₂;
2. Рекурсивно строим диаграммы Вороного для S₁ и S₂;
3. Объединить полученные диаграммы и получить диаграмму для S.

Общее описание алгоритма не сложное, однако в алгоритме есть свои тонкости, с которыми Вы можете ознакомиться в указанной книге.

Применения

Исчерпывающий список всех применений диаграммы Вороного находится здесь, я же упомяну некоторые из них, показавшиеся мне наиболее интересными (многие сведения взяты из [1]):

В программировании, разработке игр и картографии

В вычислительной геометрии диаграмма Вороного нужна прежде всего для решения задачи близости точек, а точнее, особый выигрыш диаграмма даёт в решении задачи ВСЕ ближайшие соседи (не те, которые громко включают музыку, а хотя..), потому как аналогичные ей способы не так просты ([1]). Используя диаграмму Вороного можно построить выпуклую оболочку за (смотреть на «лучевые» рёбра, находить сайты, к которым они принадлежат, и включать их в оболочку).

Также существует важная связь диаграммы Вороного с триангуляцией Делоне, которая позволяет строить одно по другому за , потому как они являются двойственными друг к другу (соединяем рёбрами соседние сайты, в итоге получим триангуляцию Делоне — викиконспекты).

Пример использования диаграммы Вороного в геймдеве можно найти, например, в этой статье — здесь система навигации в игровом движке основана на диаграмме.

Не исключено, и даже вероятно, что различный геолокационный софт использует диаграммы Вороного. Геолокационные рекомендательные системы могут использовать диаграмму Вороного для определения, например, ближайшего к вам продуктового магазина, для различного поиска и анализа местоположения.

Здесь же можно упомянуть и применение диаграммы в картографии — для очерчивания границ регионов и дальнейшего анализа на их основе. Да и вообще, любые географические диаграммы, показывающие распределение чего либо, можно наглядно проиллюстрировать с помощью раскрашенных диаграмм Вороного, и там будет виден переход нужного нам показателя (например, температуры):

Ну и, конечно, можно делать различные фильтры-обработчики фото с помощью диаграммы Вороного, получая некую «мозаику».

Но это только начало её применения.

В архитектуре и дизайне

Весьма логично, что людям в голову пришла идея использовать диаграмму Вороного в архитектуре и дизайне, поскольку она сама по себе является красивым рисунком, своего рода «геометрической паутиной», так что есть много случаев применения её в качестве одного из основных элементов композиции или даже каркаса всего творения. Примеры:

В археологии

Из [1]:

В археологии многоугольники Вороного используются для нанесения на карту ареала применения орудия труда в древних культурах и для изучения влияния соперничающих центров торговли.
В экологии возможности организма на выживание зависят от числа соседей, с которыми он должен бороться за пищу и свет.

— что вполне логично, ведь обычно за любое «выживание» бьются именно соседствующие регионы.

В моделировании и распознавании

В этой статье 3D-диаграмма Вороного не рассматривается, однако она имеет многие приложения в физике и 3D-моделировании объектов. Разного рода сетки (и скелеты) объектов в пространстве можно построить с помощью диаграммы Вороного (однако чаще с помощью триангуляции Делоне).

3D-сканирование (и компьютерное зрение (computer vision)) различных объектов тоже может использовать диаграмму Вороного и триангуляцию Делоне, также это тесно связано с робототехникой — движение робота с учётом препятствий на пути.

В биологии и химии

Из [1]:

Совместное влияние электрических и близкодействующих сил, для изучения которых строятся сложные диаграммы Вороного, помогает определять структуру молекул.

Вот ещё одна интересная статья про применение диаграммы Вороного.

Интересные факты

Природа — удивительная вещь, ведь оказывается, что окрас жирафа фактически имеет вид диаграммы Вороного. Это видно невооружённым глазом:

Также примечательно и следующее наблюдение, которое показывает, что диаграмму можно увидеть даже на листьях деревьев:

Кстати, чуть меньше, чем год назад, с диаграммой Вороного также был связан один инцидент — здесь она использовалась в логотипе проекта «Новая Москва».

И — напоследок — видео, в котором видно, как диаграмма появляется при росте окружностей с центрами в сайтах:

Здесь можно провести аналогию с распространением пожара, имеющего несколько очагов возгорания.

Ну, вот и подошло к концу наше небольшое путешествие в мир, где все знают, кто находится к ним ближе всего. Надеюсь, вам понравилась эта статья и вы действительно узнали что-то новое, полезное и интересное.

Спасибо за внимание!

Список литературы

[1] Препарата Ф., Шеймос М. Вычислительная геометрия. Введение (1989)
[2] Александров А. Д., Вернер А. Л., Рыжик В. И. Стереометрия. Геометрия в пространстве
[3] Joseph O'Rourke. Computational Geometry in C

» Статью подготовил студент 1-го курса ФИВТ МФТИ Захаркин Илья.
» В написании библиотеки также помогали студентки 1-го курса ФИВТ МФТИ Кириленко Елена и Касимова Надежда.
» В тестировании библиотеки помогал Ярослав Спирин.
» Отдельная благодарность менторам Гадельшину Ильнуру и Гафарову Рустаму.

Автор: QuantZero

Источник