Алгоритм определения вхождения точки в полигон

16.05.2017 25 Comments

Практика программирования, Разработка

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

Недавно возник вопрос, как можно определить входит ли точка в произвольный полигон. Различных решений в интернетах море, но, как ни странно, никто не реализовывал эту задачу под 1С. Многие решения были достаточно громоздки, использовали тригонометрические функции, что приводило к довольно длительным вычислениям на больших массивах данных. Я выбрал алгоритм из http://habrahabr.ru/post/125356/ как отрабатывающий за время кратное количеству вершин в полигоне и требующий на каждую вершину лишь: 2 сложения, 4 вычитания, 1 умножение, 1 деление и одну операция взятия знака числа. Конечно, интерпретатор 1С далёк от лаконичности, но стало интересно в какие сроки этот алгоритм будет работать на реальных данных. Расчет реализован с помощью функции, принимающей на входе 2 параметра: массив структур, содержащих координаты вершин полигона и структуру, содержащую координаты точек. На выходе функция возвращает «Истина» если точка входит в полигон и «Ложь» в ином случае.

С теорией по работе алгоритма можно ознакомиться по следующей ссылке: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.88.5498&rep=rep1&type=pdf или из pdf во вложенном файле.

Код функции следующий:

//Полигон — Массив структур{Х,У} //Точка — Структура{Х,У}

Функция ПроверитьТочку(Знач Полигон,Знач Точка) Экспорт  

Перем Результат;    

Если Полигон.Количество() < 3 Тогда   

Возврат Ложь;  

КонецЕсли;  

Результат = 0;  

Шаблон = Новый Массив;  

Строка1 = Новый Массив;  

Строка1.Добавить(0);  

Строка1.Добавить(1);  

Шаблон.Добавить(Строка1);  

Строка2 = Новый Массив;  

Строка2.Добавить(3);  

Строка2.Добавить(2);  

Шаблон.Добавить(Строка2);    

Размер = Полигон.Количество()-1;  

Конец = Полигон[Размер];  

ПредыдущаяТочка = Новый Структура(«Х,У»,Полигон[Размер].Х,Полигон[Размер].У);  

ПредыдущаяТочка.Х = ПредыдущаяТочка.Х — Точка.Х;  

ПредыдущаяТочка.У = ПредыдущаяТочка.У — Точка.У;  

Если ПредыдущаяТочка.У < 0 Тогда   

Х = 1;  

Иначе   

Х = 0;  

КонецЕсли;    

Если ПредыдущаяТочка.Х < 0 Тогда   

У = 1;  

Иначе   

У = 0;  

КонецЕсли;  

Пред_КУ = Шаблон[Х][У];  

Для Сч = 0 По Размер Цикл   

ТекТочка = Новый Структура(«Х,У»,Полигон[Сч].Х,Полигон[Сч].У);   

ТекТочка.Х = ТекТочка.Х — Точка.Х;   

ТекТочка.У = ТекТочка.У — Точка.У;   

Если ТекТочка.У < 0 Тогда    

Х = 1;   

Иначе    

Х = 0;   

КонецЕсли;      

Если ТекТочка.Х < 0 Тогда    

У = 1;   

Иначе    

У = 0;   

КонецЕсли;   

Ку = Шаблон[Х][У];   

ДельтаКу = Ку — Пред_Ку;   

Если ДельтаКу = -3 Тогда    

Результат = Результат + 1;   

ИначеЕсли ДельтаКу = 3 Тогда    

Результат = Результат — 1;   

ИначеЕсли ДельтаКу = -2 Тогда    

Если ПредыдущаяТочка.Х*ТекТочка.У >= ПредыдущаяТочка.У*ТекТочка.Х Тогда     

Результат = Результат + 1;    

КонецЕсли;   

ИначеЕсли ДельтаКу = 2 Тогда    

Если НЕ (ПредыдущаяТочка.Х*ТекТочка.У >= ПредыдущаяТочка.У*ТекТочка.Х) Тогда     

Результат = Результат — 1;    

КонецЕсли;   

КонецЕсли;   

ПредыдущаяТочка = ТекТочка;   

Пред_КУ = Ку;  

КонецЦикла;  

Возврат  НЕ(Результат = 0);

КонецФункции

 

 

Related Posts

25 Comments

  1. Serj1C

    (0) > но стало интересно в какие сроки этот алгоритм будет работать на реальных данных

    Какие выводы? быстро ли? правильно?

    Reply ↓
  2. logos

    (1) Serj1C, Ну, конечно не быстро, не идеально, но: обсчет 400 точек на вхождение в 50 зон по 5-87 точек в каждой проходит за 30 секунд. Тут нужно смотреть, что будет более ресурсоемко: работать с накладными расходами на транслятор языка 1С или с накладными расходами на COM или взаимодействие между процессами в каком-либо другом виде. Считаю, что для больших объемов данных будет более целесообразно сделать внешнюю компоненту, которая будет брать на входе 2 массива: точек и полигонов и возвращать массив результатов. В моём случае такой необходимости не было, не те объемы данных.

    Reply ↓
  3. Serj1C
    Reply ↓
  4. Serj1C

    (3) Serj1C, Файл с обработкой не могу приложить — в комментариях кнопочки нет(

    Reply ↓
  5. logos

    (4) Serj1C, Интересно, хотя это всё полумеры и по делу нужно выносить вычисления в библиотеку.

    Reply ↓
  6. CagoBHuK

    Все работает для проскости. Но если полигон располагается на северном полюсе и представляет собой плоскость 200км на 200км? Необходимо добавить расчет кривизны земной поверхности.

    Reply ↓
  7. logos

    (6) CagoBHuK, Да, безусловно, Вы правы. Но в этом случае для описания поверхности потребуется другая модель определения попадания точки в область. Необходимо будет определять множество полигонов, видимо треугольных, как минимально допустимая плоскость, это множество должно будет описывать все неровности рельефа. Вы же не забыли, что точка может находиться выше или ниже полигона? В этом случае задача вырождается в определение, попадает ли точка в треугольник или нет для каждого треугольника из множества описывающего поверхность, являющуюся нашей областью.

    Reply ↓
  8. CagoBHuK

    Если говорить о гуглокартах (а это, несомненно, они изображены на спойлере), то вполне хватит простого перевода изгиба поверхности планеты в плоскость через тригонометрическую функцию sin(широта). Почему я клоню в эту сторону? Да потому, что уже решал подобные задачи, а выложить на Инфостарт не додумался. Поэтому я направляю Вас в правильном направлении, чтобы не наступать на те же грабли, что и я.

    Reply ↓
  9. SunShinne

    Автору огромное спасибо! Однозначно +. Единственное есть два нюанса — во первых в коде закралась ошибка (путаница между латинскими X Y и кириллическими Х У), во вторых алгоритм работает на вхождение именно внутри полигона, но вот если точка на границе полигона то алгоритм выдает отрицательный результат. Но это вопрос спорный, однако отметить стоит. Еще раз спасибо — помогло!

    Reply ↓
  10. tav13

    (8) CagoBHuK, почитав вас сделал обработку перевода координат на плоскость http://infostart.ru/public/295388/

    Reply ↓
  11. frkbvfnjh

    По крайней мере с координатами Google карты не работает. Спер алгоритм на JS, там все вхождения четко определил а этот какую то ерунду выдает, я так понял, что координаты еще предварительно нужно обработать обработкой http://infostart.ru/public/295388/ что ли, что бы этот алгоритм заработал?

    Reply ↓
  12. frkbvfnjh

    я так понимаю широта это X долгота это Y, или нет?

    Reply ↓
  13. ildarovich

    Я бы попробовал применить в этой задаче готовый объект «ГеографическаяСхема». Там есть метод «ПроверитьПоРасположению». По описанию делает как раз то, что надо.

    Reply ↓
  14. Sergius79

    Спасибо, за решение. В благодарность прикладываю эту же функцию для 1С 7.7

    // ======================================
Функция ГлПроверитьТочку(СписПолигонов,ТочкаХ,ТочкаУ) Экспорт
Перем Результат, ШаблонКоординат;

Результат = 0;
Для Инд=1 по СписПолигонов.РазмерСписка() Цикл
Результат = 0;
Полигон=СписПолигонов.ПолучитьЗначение(Инд);
Размер=Полигон.КоличествоСтрок();
Если  Размер< 3 Тогда
Возврат 0;
КонецЕсли;
ШаблонКоординат=СоздатьОбъект(«ТаблицаЗначений»);
ШаблонКоординат.НоваяКолонка(«У»,»Число»);
ШаблонКоординат.НоваяКолонка(«У1″,»Число»);

ШаблонКоординат.НоваяСтрока();
ШаблонКоординат.У = 0;
ШаблонКоординат.У1 = 1;
ШаблонКоординат.НоваяСтрока();
ШаблонКоординат.У = 3;
ШаблонКоординат.У1 = 2;


Полигон.ПолучитьСтрокуПоНомеру(Размер);
Конец = Полигон;

ПредыдущаяТочкаХ = Полигон.Х;
ПредыдущаяТочкаУ = Полигон.У;

ПредыдущаяТочкаХ = ПредыдущаяТочкаХ — ТочкаХ;

ПредыдущаяТочкаУ = ПредыдущаяТочкаУ — ТочкаУ;

Если ПредыдущаяТочкаУ < 0 Тогда
Х = 1;
Иначе
Х = 0;
КонецЕсли;

Если ПредыдущаяТочкаХ < 0 Тогда
У = 1;
Иначе
У = 0;
КонецЕсли;

Пред_КУ = ШаблонКоординат.ПолучитьЗначение(Х+1,У+1);

Полигон.ВыбратьСтроки();
Пока Полигон.ПолучитьСтроку()=1 Цикл
ТекТочкаХ = Полигон.Х;
ТекТочкаУ = Полигон.У;

ТекТочкаХ = ТекТочкаХ — ТочкаХ;
ТекТочкаУ = ТекТочкаУ — ТочкаУ;

Если ТекТочкаУ < 0 Тогда
Х = 1;
Иначе
Х = 0;
КонецЕсли;

Если ТекТочкаХ < 0 Тогда
У = 1;
Иначе
У = 0;
КонецЕсли;

Ку = ШаблонКоординат.ПолучитьЗначение(Х+1,У+1);

ДельтаКу = Ку — Пред_Ку;

Если ДельтаКу = -3 Тогда
Результат = Результат + 1;
ИначеЕсли ДельтаКу = 3 Тогда
Результат = Результат — 1;
ИначеЕсли ДельтаКу = -2 Тогда
Если ПредыдущаяТочкаХ*ТекТочкаУ >= ПредыдущаяТочкаУ*ТекТочкаХ Тогда
Результат = Результат + 1;
КонецЕсли;
ИначеЕсли ДельтаКу = 2 Тогда
Если НЕ (ПредыдущаяТочкаХ*ТекТочкаУ >= ПредыдущаяТочкаУ*ТекТочкаХ) Тогда
Результат = Результат — 1;
КонецЕсли;
КонецЕсли;

ПредыдущаяТочкаХ = ТекТочкаХ;
ПредыдущаяТочкаУ = ТекТочкаУ;
Пред_КУ = Ку;
КонецЦикла;
Если Результат <>  0 Тогда
Возврат 1;
КонецЕсли;
КонецЦикла;
Возврат  ?(Результат = 0,0,1);
КонецФункции

    Показать

    Reply ↓
  15. DanilaDru

    (3)

    | КОГДА СУММА(Развернуто.Результат) = 0

    Со временем использования не возникало вопроса что сумму надо проверять на > 0 для Истина ?

    Столкнулся с тем что в выборку попадают противоположные относительно друг друга полигоны. Так как у одного из них результат = -1, а у другого 1.

    Reply ↓
  16. Жолтокнижниг

    (13) К сожалению этот вариант работает не корректно

    Reply ↓
  17. vs84

    (16) Проверяли, заработало, но некорректно определяет вхождение? Планировали использовать этот метод для аналогичной задачи (предполагаем, что он определяет с учетом сферичности поверхности), но смутил ваш комментарий.

    Reply ↓
  18. Жолтокнижниг

    (17) Сейчас не помню в чем конкретно были проблемы, но нам не подошел этот метод.

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

    В итоге использую алгоритм статьи и http://infostart.ru/public/295388/ с небольшой шлифовкой. Работает на ура.

    Reply ↓
  19. kuzz

    (2) Спасибо за публикацию! Отлично работает!

    Reply ↓
  20. kuzz

    (11) У нас все работает на ура

    Reply ↓
  21. kuzz

    (7) Скажите пожалуйста, а каково Ваше мнение по поводу перевода координат на плоскость? Нужно ли этим заниматься если задача стоит в определении принадлежности точек к областям среднестатистического города нашей страны? С точки зрения увеличения точности

    Reply ↓
  22. frkbvfnjh

    (18) А с какой небольшой шлифовкой? Может сделаете статью уже с окончательным вариантом алгоритма, который 100% работает, а то сколько не читаю статей в сети, везде потом в конце пишут, что нужно допилить

    Reply ↓
  23. Жолтокнижниг

    (22) К сожалению не могу предоставить код, вот примерный алгоритм

    Вводная

    1. Используется проекция UTM.

    2. Используются небольшие геозоны.

    Алгоритм

    1. Для перевода координат необходимо определить UTM зону (реквизит Zome из обработки https://infostart.ru/public/295388/)

    2. Зона определялась, кажется, по медианной долготе (код определения зоны по координатам есть в обработке)

    3. Используя алгоритм обработки и зону переводим геозону в полигон

    4. По алгоритму из этой статьи определяем вхождение

    Поток геоданных небольшой, поэтому оптимизаций не делали.

    ПС: работу с геоданными поддерживают СУБД(MS и PostgreSQL) при больших объемах данных лучше использовать их.

    Reply ↓
  24. acanta

    Можно подробнее про геоданные в СУБД ?

    Никогда не слышала об этом.

    Reply ↓
  25. frkbvfnjh

    (23) Спасибо

    Reply ↓

Leave a Comment

Ваш адрес email не будет опубликован. Обязательные поля помечены *