graduate student from 01.01.2018 to 01.01.2024
Russian Federation
The problem of global route planning for a mobile robot between two given points in a known area with static obstacles is considered. To solve the problem of constructing a route in an area with a large number of obstacles of complex shape, an integrated approach based on graph theory methods is proposed. It includes the Voronoi diagram, visibility graph and the Dijkstra's algorithm. At the first stage, the study area is represented as a polygonal object, the space outside the object is considered as obstacles. Next, to get a safe distance from obstacles, an internal buffer of the polygonal object is built using the Minkowski difference. Then the vertices of the polygon are compacted, and the Voronoi polygons are constructed from the resulting vertices. The median axis of the polygon is calculated from the Voronoi polygons. Then the Dijkstra's algorithm is applied to calculate the shortest path. The resulting path is used to construct a visibility graph, and the Dijkstra's algorithm is reapplied to the resulting graph. The proposed approach allows to build a route that is optimal in terms of length and distance to obstacles. It significantly reduces the computational complexity of constructing a visibility graph. The approach was implemented in the freely distributed QGIS geographic information system for planning the route of a mobile robot in an aquatic environment. The results of the experiment showed that the Voronoi diagram reduced the number of vertices required to construct the visibility graph by 8.3 times, while the visibility graph improved the path obtained from the Voronoi diagram by 8%. The proposed approach can be used for global planning of routes for mobile robots in various environments.
mobile robots, route planning, territory, obstacle, shortest path, the Voronoi diagram, visibility graph, the Dijkstra's algorithm
Введение
В настоящее время мобильные роботы используются практически повсеместно – от бытовой сферы и промышленности до применения в научных исследованиях. Ключевой особенностью таких роботов является отсутствие привязки к конкретной точке пространства и способность перемещаться по неограниченной территории. Существует большое количество разновидностей мобильных роботов. По среде передвижения их можно разделить на наземные транспортные средства, воздушные аппараты, надводные и подводные мобильные аппараты, а также гибридные роботы, которые могут перемещаться сразу в нескольких средах. По способу управления мобильные роботы можно разделить на управляемые оператором, полуавтономные и полностью автономные (беспилотные аппараты), которые способны перемещаться в среде и выполнять задачи самостоятельно, без участия человека.
Несмотря на разнообразие мобильных роботов, все они для выполнения своей цели должны уметь планировать маршрут от точки запуска до точки назначения. Планирование маршрута можно разделить на глобальное планирование, т. е. планирование движения по известной местности со статическими препятствиями, и локальное планирование, т. е. планирование в динамической среде непосредственно на местности. Как правило, предварительно составленный в ходе глобального планирования маршрут далее корректируется роботом в ходе локального планирования в режиме реального времени. При этом задача составления оптимального глобального маршрута мобильного робота остается не менее актуальной, чем локальное планирование, особенно в условиях перемещения робота по большой территории со значительным количеством препятствий.
Методы планирования пути
К методам планирования пути относятся методы на основе графов, метод клеточной декомпозиции, метод потенциальных полей, методы нейронных сетей [1, 2]. Методы на основе графов включают в себя граф видимости, диаграмму Вороного, метод вероятностной дорожной карты и метод быстро исследующих случайных деревьев.
В случае графа видимости препятствия представляются в виде полигональных объектов. Граф строится по вершинам препятствий, вершины должны быть видимы. Построение графа видимости рассматривается в работах [3, 4]. Следует отметить, что с помощью графа видимости можно получить кратчайший путь, но при увеличении количества вершин значительно возрастает сложность метода, поэтому граф видимости чаще всего строят на объектах с малым числом препятствий.
В случае диаграммы Вороного плоскость с препятствиями разбивается многоугольниками таким образом, что точки каждого многоугольника будут ближе к точке препятствия, чем к любой другой. Сложность алгоритма ниже, чем графа видимости. Метод поиска пути с помощью диаграммы Вороного и итеративного сглаживания рассматривается в работе [5]. Полученный путь является гладким и безопасным, но не оптимальным с точки зрения длины.
Метод потенциальных полей основан на предположении, что конечная точка притягивает объект, а препятствия его отталкивают. Метод вычислительно эффективен, но подвержен проблеме локальных минимумов, в результате чего робот может зациклиться. В работе [6] используется подход, сочетающий в себе граф видимости, диаграмму Вороного и метод потенциальных полей. Диаграмма Вороного используется для получения безопасной от препятствий области. Далее в ней итеративно ищется кратчайший путь с помощью метода потенциальных полей и графа видимости. Следует отметить, что метод достаточно сложен в реализации.
В методе клеточной декомпозиции территория разбивается на клетки, сетка из клеток рассматривается как связный граф. Результат метода сильно зависит от масштаба сетки. С повышением точности сетки сложность метода возрастает, при низкой точности путь может получиться неоптимальным. Поиск пути при помощи метода клеточной декомпозиции рассматривается в работе [7]. Для вычисления оптимального пути используется волновой метод с поиском в ширину. В работе [8] применяются методы клеточной декомпозиции и графа видимости для построения маршрутов в помещениях по заранее известному плану, но не рассматривается возможность построения маршрута на открытой территории с большим количеством препятствий.
Авторы [9] рассматривают возможность локального планирования маршрута робота с помощью сверхточных нейронных сетей. Была обучена нейронная сеть YOLOv5 для перемещения мобильного робота в искусственном бассейне в рамках соревнования. Для применения в естественной среде метод потребует доработки и дополнительного обучения.
Несмотря на разнообразие методов поиска пути они имеют недостатки в случае большого количества и сложной формы препятствий: высокая вычислительная сложность алгоритма, получение неоптимального пути на выходе, зацикливание в области локального минимума или необходимость повторного обучения. Таким образом, существует потребность в разработке подхода к эффективному планированию пути мобильного робота на большой территории с большим количеством препятствий.
Подход к глобальному планированию маршрута мобильного робота на основе графовых методов
Графовые методы являются одними из самых точных методов построения кратчайшего пути, но исследователи редко используют их для построения маршрутов на территории со значительным количеством препятствий, обосновывая это высокой вычислительной сложностью. Особенно это касается графа видимости, т. к. сложность алгоритма построения графа видимости может описываться функцией
ƒ(n) = O(n2), (1)
где n – количество вершин; О – нотация для описания сложности алгоритма, и, следовательно, квадратично возрастать с увеличением количества вершин. При этом существует возможность эффективно использовать граф видимости для построения оптимального пути мобильного робота в среде с большим количеством препятствий.
Целью данного исследования является разработка подхода к глобальному планированию маршрута мобильного робота в среде с большим количеством препятствий на основе графовых методов, а именно на основе диаграммы Вороного и графа видимости. Критериями оптимальности в данном случае служат минимальная длина маршрута и безопасное расстояние от препятствий.
Подход основан на предположении о том, что в случае протяженной территории с большим количеством препятствий сложной формы для нахождения маршрута между двумя заданными точками не требуется строить граф видимости на всей территории. Достаточно построить диаграмму Вороного, которая обладает значительно меньшей сложностью, описываемой функцией ƒ(n) = O(nlog(n)), где n – количество вершин, используемых для построения диаграммы, затем найти кратчайший путь в ней, построить по вершинам пути граф видимости и найти кратчайший путь в графе видимости. Такой подход позволит объединить сильные стороны методов: для диаграммы Вороного – скорость построения и удаленность пути от препятствий, для графа видимости – возможность строить кратчайшие траектории.
Основными этапами предлагаемого подхода являются:
– подготовка исходных данных;
– обеспечение безопасного расстояния от препятствий;
– построение диаграммы Вороного и нахождение кратчайшего пути в ней с помощью алгоритма Дейкстры на основе очереди с минимальным приоритетом;
– построение графа видимости и повторное нахождение кратчайшего пути с помощью алгоритма Дейкстры на основе очереди с минимальным приоритетом.
Далее рассмотрим подробнее каждый этап.
Подготовка исходных данных
В первую очередь территорию, на которой возможно перемещение робота, необходимо представить в виде многоугольника. Пусть M – многоугольник, {(Xi, Yi)}, i = 1, 2, …, n – количество вершин n-угольника (количество вершин, используемых для дальнейшего построения диаграммы Вороного).
Пространство за пределами многоугольника будет считаться препятствиями на пути робота. Характер территории может быть любым, подходящим для перемещения робота (рис. 1).
Рис. 1. Представление разрешенной территории в виде многоугольника
Fig. 1. Representation of the permitted territory in the form of a polygon
Пусть даны две точки s, t ϵ М, между которыми требуется построить маршрут. Возможные траектории движения в пределах многоугольника можно представить в виде графа G = (V, E), где V – множество вершин графа G, V Î {(X, Y)}È{s, t}; E – множество ребер графа G. Ребра представляют собой все возможные соединения между парами вершин графа, которые не пересекают границ препятствий. Необходимо найти кратчайший путь по графу между точками s и t.
Для нахождения кратчайшего пути в многоугольнике только с помощью графа видимости потребовалось бы создать ребра между всеми видимыми вершинами, и вычислительная сложность метода описывалась бы формулой (1). Кроме того, путь между начальной и конечной точками зачастую пролегает лишь по небольшой части территории, вследствие чего построение графа видимости на всей территории будет избыточным. Поэтому перед построением графа необходимо выполнить предварительный этап, который позволит сократить сложность его построения.
Обеспечение безопасного расстояния от препятствий
Как правило, самый короткий путь по территории с препятствиями проходит по границам препятствий, что может привести к столкновениям мобильного робота. Для исключения этого, а также для учета размеров робота и пространства на повороты, необходимо установить минимальное безопасное расстояние от препятствий. Это также значительно облегчает планирование, т. к. робота при планировании можно будет представить в виде точки в пространстве.
Так как в данной работе многоугольником представлены не препятствия, а разрешенная для перемещения территория, то для обеспечения безопасного расстояния ее необходимо уменьшить. Уменьшение территории многоугольника можно осуществить путем применения к нему операции буферизации с отрицательным значением буферного расстояния. С точки зрения математических методов в основе буферизации лежит вычисление суммы или разности Минковского между исходным объектом и кругом с радиусом, равным буферному расстоянию (рис. 2).
а б
Рис. 2. Буферизация многоугольника: а – внешняя буферизация; б – внутренняя буферизация: А, В – множества
Fig. 2. Buffering of the polygon: a – external buffering; б – internal buffering; A, B – sets
Суммой Минковского двух подмножеств A и B является множество C, состоящее из сумм всех векторов из A и B [10, 11]:
С = А + В = {a + b | a ϵ A, b ϵ B}.
Разностью Минковского двух подмножеств A и B является максимальное множество C, которое в сумме с B даст подмножество A [11]. Таким образом
А – В = {С | С + В = А}.
Сумма Минковского позволяет расширить границы объекта (внешняя буферизация, см. рис. 2, а), разность – уменьшить (внутренняя буферизация, см. рис. 2, б).
Операция буферизации является стандартной в геоинформационных системах, таких как QGIS, и не требует дополнительной реализации. Буферное расстояние задается пользователем.
Построение диаграммы Вороного
Этап включает уплотнение вершин многоугольника, построение диаграммы Вороного, выделение из нее срединной оси и поиск кратчайшего пути с помощью алгоритма Дейкстры.
Первым шагом является уплотнение вершин, т. к. вершины в многоугольнике располагаются на разном расстоянии друг от друга, в местах наиболее разреженного их размещения будет минимальное количество полигонов Вороного, что может снизить точность полученной срединной оси. Поэтому необходимо выполнить уплотнение вершин следующим образом. Пусть препятствия представляют собой конечный набор точек {p1, … pn}
в эвклидовой плоскости (рис. 3, а).
Рис. 3. Построение диаграммы Вороного:
а – исходные точки; б – диаграмма Вороного; в – добавление начальной и конечной точки; г – кратчайший путь
Fig. 3. Construction of the Voronoi diagram:
a – starting points; б – the Voronoi diagram; в – adding the starting and ending points; г – the shortest path
Пара точек p1 = (x1, y1) и p2 = (x2, y2) лежит на прямой m, заданной уравнением
Расстояние между p1 и p2 определяется по формуле эвклидова расстояния
(2)
где D – эвклидово расстояние.
Для уплотнения отрезка [p1, p2] требуется найти такую точку p3 = (x3, y3), которая находится на расстоянии d от точки p1:
где d – расстояние уплотнения; d < D.
Следующим шагом является построение диаграммы Вороного. Для этого вокруг каждой точки-препятствия pk необходимо построить многоугольник Rk таким образом, чтобы расстояние от него до точки pk было меньше или равно минимальному расстоянию до любой другой точки pj [12]:
,
где Х – множество всех точек евклидова пространства; d – евклидово расстояние между точками x и pk, определяемое по формуле (2).
Диаграмма Вороного представлена на рис. 3, б.
Следующим шагом является получение срединной оси из диаграммы Вороного. Пуcть диаграмма Вороного состоит из множества отрезков L, разрешенная для перемещения территория представлена многоугольником M. Сначала необходимо удалить отрезки за пределами многоугольника M (l ∉ M). Далее необходимо удалить ребра диаграммы Вороного, которые ведут к препятствиям, т. е. к границам многоугольника. Для каждого отрезка диаграммы Вороного l = [(x1, y1), (x2, y2)] проверяем пересечение с границей препятствия h =
= [(x3, y3), (x4, y4)]. Отрезок l подлежит удалению, если l ∩ h. Отрезки l и h пересекаются в точке
(x0, y0), если существуют такие параметры s0, t0, для которых существует решение системы уравнений
при 0 ≤ s0, t0 ≤ 1.
В результате вышеописанных действий был получен взвешенный неориентированный граф Gd = (V, E, w), где вес ребра w равен евклидову расстоянию между узлами по формуле (2). Граф содержит серединную ось многоугольника.
Далее необходимо добавить начальную (s) и конечную (t) точки пути в граф. Предполагая, что граф содержит n вершин, выбираем два узла s' и t' ∊ V, которые имеют минимальное эвклидово расстояние (2) до точек s и t соответственно. Добавляем вершины в граф и добавляем ребра между s и s', t и t'. Результат представлен на рис. 3, в.
Следующим шагом является поиск кратчайшего пути в графе. Для этого используется классический алгоритм жадного поиска – алгоритм Дейкстры. Алгоритм Дейкстры работает для графов с положительным весом ребер. Поиск пути осуществляется последовательным перебором всех вершин и расчетом расстояния от начальной точки до каждой вершины. Блок-схема используемой в подходе реализации алгоритма Дейкстры представлена на рис. 4.
Рис. 4. Блок-схема алгоритма Дейкстры
Fig. 4. Block diagram of the Dijkstra's algorithm
Итоговым результатом этого этапа является кратчайший путь, полученный из диаграммы Вороного (см. рис. 3, г). Он максимально удален от границ препятствий, но не оптимален с точки зрения длины и содержит избыточные повороты.
Построение графа видимости
В общем случае при построении графа видимости вершины препятствий становятся узлами графа, соединения между видимыми вершинами – ребрами. Общий принцип построения графа видимости демонстрирует рис. 5.
Рис. 5. Граф видимости:
а – исходные точки и препятствия; б – полный граф; в – граф видимости; г – кратчайший путь
Fig. 5. Visibility graph:
a – starting points and obstacles; б – full graph; в – visibility graph; г – shortest path
В данном исследовании многоугольником представлена разрешенная для движения робота территория.
Представим взвешенный неориентированный граф Gv:
Gv = (V, E, w),
где V – количество вершин графа; Е – конечное множество ребер графа, соединяющих подмножество вершин; w – вес ребра графа, равный евклидову расстоянию между вершинами.
Вершины графа Gv равны вершинам кратчайшего пути, полученного из диаграммы Вороного. Количество ребер графа равно числу сочетаний из q объектов по 2:
,
где q – количество вершин в кратчайшем пути, полученном из диаграммы Вороного.
С учетом вышеизложенного для графа справедливы следующие условия:
Вторым шагом является сокращение построенного графа. Ребра графа должны удовлетворять условию видимости вершин. Две вершины (xa, ya) и (xb, yb) считаются видимыми, если между ними можно провести прямую таким образом, что эта прямая не проходит через вершину препятствия (xc, yc). Таким образом, должно выполняться условие
Если это условие не выполняется, то ребра пересекают границу буферизованного многоугольника и подлежат удалению. Итоговый результат представлен на рис. 5, в.
Третьим шагом данного этапа является поиск кратчайшего пути в графе видимости путем повторного применения алгоритма Дейкстры. Подход аналогичен представленному на рис. 4, результат поиска пути представлен на рис. 5, г.
Результаты исследования
Для исследования был взят векторный полигон (рис. 6, а), представляющий собой один из притоков реки Северная Двина на территории Архангельской области.
Исследуемый объект имеет площадь 1,7 км2 и представлен 4 297 вершинами. На рис. 6, а также отмечены начальная и конечная точки маршрута.
Для выбранной области был реализован вышеописанный подход в геоинформационной системе QGIS.
Для полигона, представленного на рис. 6, а, был построен внутренний буфер в 5 м, в результате чего была получена уменьшенная копия многоугольника. В результате буферизации и уплотнения многоугольника был получен многоугольник, содержащий 5 320 вершин. На рис. 6, б представлена построенная диаграмма Вороного. На рис. 6, в–д продемонстрированы результаты формирования скелета многоугольника из диаграммы Вороного. Найденный с помощью алгоритма Дейкстры путь представлен на рис. 6, е, ж.
Рис. 6. Результат построения кратчайшего пути из диаграммы Вороного:
а – исходная область; б – полигоны Вороного; в, г – сокращенная диаграмма Вороного; д – скелет многоугольника; е, ж – кратчайший путь из диаграммы Вороного
Fig. 6. Result of constructing the shortest path from the Voronoi diagram:
a – research area; б – the Voronoi polygons; в, г – reduced Voronoi diagram; д – the skeleton of the polygon; е, ж – the shortest path from the Voronoi diagram
Путь является безопасным, но неоптимальным с точки зрения длины и содержит избыточные повороты, поэтому далее был построен граф видимости, в котором найден кратчайший путь. На рис. 7 представлен итоговый результат практического применения подхода – кратчайший путь, полученный с применением комплексного подхода (диаграмма Вороного и граф видимости).
Рис. 7. Кратчайшие пути. Сплошная линия – путь из диаграммы Вороного,
пунктирная линия – путь из графа видимости
Fig. 7. The shortest paths. Solid line – path based on the Voronoi diagram,
dashed line – the path from the visibility graph
Для сравнения на том же изображении представлен только кратчайший путь из диаграммы Вороного.
Путь из диаграммы Вороного имеет длину 6 144 м, в то время как путь из графа видимости – 5 669 м. Таким образом, длина пути была сокращена примерно на 8 %. Кроме того, путь из графа видимости получился более плавным, чем путь из диаграммы Вороного, даже без применения дополнительных методов сглаживания. При этом он удовлетворяет условию безопасности (удален на минимальное расстояние от препятствий) и является вычислительно эффективным в отличие от полного графа видимости. Количество вершин, используемых для построения графа видимости, удалось сократить в 8,3 раза, что также положительно повлияло на время работы алгоритма графа видимости. Наибольшую часть времени в предложенном подходе в этом и других экспериментах заняло построение диаграммы Вороного, что позволяет предположить, что сложность предложенного подхода можно выразить той же функцией, что и сложность диаграммы Вороного, т. е. O(nlog(n)), где n – количество вершин многоугольника, по которым строится диаграмма Вороного.
С целью автоматизации предложенного подхода был разработан модуль для геоинформационной системы QGIS на языке программирования Python. Внешний вид интерфейса модуля представлен на рис. 8, результат работы модуля представлен на рис. 9.
Рис. 8. Интерфейс модуля QGIS для автоматизации предлагаемого подхода
Fig. 8. Interface of QGIS module to automate the proposed approach
Рис. 9. Результат работы модуля – кратчайший маршрут из точки s в точку t
Fig. 9. QGIS module result – the shortest route from point s to point t
Для работы с векторными объектами была использована библиотека qgis.core, для работы с графами – библиотека NetworkX, а также ряд вспомогательных библиотек (math, queue и т. д.). Для запуска модуля требуется указать полигон с разрешенной для движения робота территорией, начальную и конечную точки маршрута, используемое разрешение (точность в метрах) и величину безопасного расстояния от границы полигона (в метрах). Пример работы модуля представлен на рис. 9. В данном случае разрешенная для передвижения робота территория представлена полигоном площадью 411 км2, кратчайший путь на основе описанного подхода был найден за 35 с, длина маршрута составила 54,7 км.
Таким образом, модуль выполняется в конечное время в пределах от 1 с для простого полигона до 35 с для полигона сложной формы и обеспечивает нахождение кратчайших маршрутов между заданными точками на заданном пользователем расстоянии от препятствий.
Заключение
В статье рассмотрен новый подход к глобальному планированию маршрутов мобильных роботов. Предлагаемый подход предназначен для построения маршрута робота между двумя заданными точками по известной карте и основан на применении методов теории графов, а именно на комплексном применении методов диаграммы Вороного, графа видимости и алгоритма Дейкстры. Результатом подхода является построение маршрута мобильного робота, оптимального с точки зрения длины и безопасного расстояния от препятствий. Подход сочетает сильные стороны обоих методов, при этом устраняет такие их недостатки, как низкая вычислительная эффективность графа видимости, избыточная длина и избыточные повороты маршрута графа Вороного. Применение графа Вороного позволяет в разы сократить количество вершин, требуемых для построения графа видимости, а граф видимости улучшает характеристики пути графа Вороного на 8 %. При этом полученный путь продолжает оставаться безопасным, предотвращая столкновение мобильного робота с препятствиями.
Для автоматизации предложенного подхода был разработан модуль для геоинформационной системы QGIS. Дальнейшая работа будет направлена на улучшение характеристик предложенного подхода, а также проведение полевых испытаний с реальным роботом.
1. Liu V. Metody planirovaniia puti v srede s prepiatstviiami (obzor) [Path planning methods in an obstacle envi-ronment (overview)]. Matematika i matematicheskoe modelirovanie, 2018, no. 1, pp. 15-58.
2. Raja P., Pugazhenthi S. Optimal path planning of mobile robots: A review. International Journal of the Physical Sciences, 2012, vol. 7 (9), pp. 1314-1320.
3. Zaeva K. A., Semenov A. B. Sistema poiska mini-mal'nogo puti v srede s poligonal'nymi prepiatstviiami [A system for finding the minimum path in an environment with polygonal obstacles]. Graphicon 2014 Proceedings. 2014. Pp. 163-166. Available at: https://studylib.ru/doc/2682224/sistema-poiska-minimal._nogo-puti-v-srede-s (accessed:13.01.2024).
4. Han-Pang Huang, Shu-Yun Chung. Dynamic visibility graph for path planning. International Conference on Intelligent Robots and Systems (IROS) (IEEE Cat. No.04CH37566). Sendai, Japan, 2004. Vol. 3. Pp. 2813-2818.
5. Bhattacharya Priyadarshi, Gavrilova M. Roadmap-Based Path Planning - Using the Voronoi Diagram for a Clear-ance-Based Shortest Path. Robotics & Automation Magazine, IEEE, 2008, vol. 15, pp. 58-66.
6. Masehian Ellips, Amin-Naseri M. R. A Voronoi diagram-visibility graph-potential field compound algorithm for robot path planning. J. Field Robotics, 2004, vol. 21, pp. 275-300.
7. Jan Gene, Chang Ki Yin, Parberry Ian. Optimal Path Planning for Mobile Robot Navigation. IEEE/ASME Trans-actions on Mechatronics, 2008, vol. 13 (4), pp. 451-460.
8. Andrienko Iu. A., Donetskii S. Iu. Razrabotka algo-ritma postroeniia optimal'nogo puti na osnove plana po-meshcheniia dlia izmeritel'nogo robota [Development of an algorithm for constructing an optimal path based on a room plan for a measuring robot]. International Journal of Open Information Technologies, 2023, no. 8, pp. 73-77.
9. Plotnikov V. A., Starykh G. K., Davtian A. A. Primenenie neironnoi seti YOLOv5 dlia raspoznavaniia ob"ektov na sorevnovaniiakh po avtonomnoi podvodnoi robototekhnike Singapore AUV Challenge [Application of the YOLOv5 Neural Network for object Recognition at the Singapore AUV Challenge Autonomous Underwater Robotics Competition]. Vestnik MGTU «Stankin», 2023, no. 1 (64), pp. 60-66.
10. Paniukov A. V. Predstavlenie summy Minkovskogo dlia dvukh poliedrov sistemoi lineinykh neravenstv [Rep-resentation of the Minkowski sum for two polyhedra by a system of linear inequalities]. Vestnik IuUrGU. Seriia: Matematicheskoe modelirovanie i programmirovanie, 2012, no. 40 (299), pp. 108-119.
11. Berg M. T., Kreveld M. J., Overmars M. H. Computational Geometry: Algorithms and Applications. Springer-Verlag, 2008. 367 p.
12. Burrough P. A., McDonnell R. A., Lloyd C. D. Nearest neighbours: Thiessen (Dirichlet/Voroni) polygons. Oxford University Press, 2015. Pp. 160-168.
