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

Публикация № 164095

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

41
Нередко встает необходимость ведения аналитики по территориальному признаку. Вручную, человеческими руками определить к какому региону относится контрагент, адрес доставки, розничная точка или ещё что-либо проблем не составляет. А если этих точек необходимо обрабатывать 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);

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

 

 

41

Скачать файлы

Наименование Файл Версия Размер
Обработка проверки вхождения точки в полигон
.epf 9,37Kb
28.11.12
38
.epf 9,37Kb 38 Скачать
Теория
.pdf 241,15Kb
28.11.12
7
.pdf 241,15Kb 7 Скачать

См. также

Специальные предложения

Комментарии
Избранное Подписка Сортировка: Древо
1. Serj1C 475 28.11.12 14:02 Сейчас в теме
(0) > но стало интересно в какие сроки этот алгоритм будет работать на реальных данных
Какие выводы? быстро ли? правильно?
2. logos 140 28.11.12 14:29 Сейчас в теме
(1) Serj1C, Ну, конечно не быстро, не идеально, но: обсчет 400 точек на вхождение в 50 зон по 5-87 точек в каждой проходит за 30 секунд. Тут нужно смотреть, что будет более ресурсоемко: работать с накладными расходами на транслятор языка 1С или с накладными расходами на COM или взаимодействие между процессами в каком-либо другом виде. Считаю, что для больших объемов данных будет более целесообразно сделать внешнюю компоненту, которая будет брать на входе 2 массива: точек и полигонов и возвращать массив результатов. В моём случае такой необходимости не было, не те объемы данных.
19. kuzz 02.07.19 17:19 Сейчас в теме
(2) Спасибо за публикацию! Отлично работает!
3. Serj1C 475 02.12.12 22:07 Сейчас в теме
Тема заинтересовала. Посидел на выходных, в теории не особо разобрался, но удалось оптимизировать код. Обошел транслятор и интерпретатор 1С тем, что вывел все вычисления в запрос. Вычисления выполняются на стороне "сервера". Вот код и пример обработки:
//Полигон	ТаблицаЗначений: НомерСтроки, Х, У
//Точки		ТаблицаЗначений: НомерСтроки, Х, У
Функция ПроверитьТочкиЗапросом(Полигон, Точки) Экспорт
	Если Полигон.Количество() < 3 Тогда
		Возврат Ложь;
	КонецЕсли;
	
	Запрос = Новый Запрос;
	Запрос.УстановитьПараметр("Полигон", Полигон);
	Запрос.УстановитьПараметр("РазмерПолигона", Полигон.Количество());
	Запрос.УстановитьПараметр("Точки", Точки);
	Запрос.Текст =
		"ВЫБРАТЬ
		|	Полигон.НомерСтроки КАК НомерТочкиПолигона,
		|	Полигон.Х,
		|	Полигон.У
		|ПОМЕСТИТЬ Полигон
		|ИЗ
		|	&Полигон КАК Полигон
		|;
		|
		|////////////////////////////////////////////////////////////­////////////////////
		|ВЫБРАТЬ
		|	Точки.НомерСтроки КАК НаименованиеТочки,
		|	Точки.Х КАК Х,
		|	Точки.У КАК У
		|ПОМЕСТИТЬ Точки
		|ИЗ
		|	&Точки КАК Точки
		|;
		|
		|////////////////////////////////////////////////////////////­////////////////////
		|ВЫБРАТЬ
		|	Точки.НаименованиеТочки,
		|	Точки.Х КАК Точка_Х,
		|	Точки.У КАК Точка_У,
		|	Полигон.НомерТочкиПолигона,
		|	Полигон.Х КАК Полигон_Х,
		|	Полигон.У КАК Полигон_У,
		|	Полигон.Х - Точки.Х КАК Разница_Х,
		|	Полигон.У - Точки.У КАК Разница_У,
		|	ВЫБОР
		|		КОГДА Полигон.Х - Точки.Х < 0
		|				И Полигон.У - Точки.У < 0
		|			ТОГДА 2
		|		КОГДА Полигон.Х - Точки.Х < 0
		|				И Полигон.У - Точки.У >= 0
		|			ТОГДА 1
		|		КОГДА Полигон.Х - Точки.Х >= 0
		|				И Полигон.У - Точки.У < 0
		|			ТОГДА 3
		|		КОГДА Полигон.Х - Точки.Х >= 0
		|				И Полигон.У - Точки.У >= 0
		|			ТОГДА 0
		|	КОНЕЦ КАК К
		|ПОМЕСТИТЬ Таб
		|ИЗ
		|	Точки КАК Точки,
		|	Полигон КАК Полигон
		|
		|ИНДЕКСИРОВАТЬ ПО
		|	Точки.НаименованиеТочки,
		|	Полигон.НомерТочкиПолигона
		|;
		|
		|////////////////////////////////////////////////////////////­////////////////////
		|ВЫБРАТЬ
		|	Таб.НаименованиеТочки,
		|	Таб.Точка_Х,
		|	Таб.Точка_У,
		|	Таб.НомерТочкиПолигона,
		|	Таб.Полигон_Х,
		|	Таб.Полигон_У,
		|	Таб.К КАК К,
		|	ТабПред.К КАК ПредК,
		|	Таб.К - ТабПред.К КАК ДельтаКу,
		|	ВЫБОР
		|		КОГДА Таб.К - ТабПред.К = -3
		|			ТОГДА 1
		|		КОГДА Таб.К - ТабПред.К = 3
		|			ТОГДА -1
		|		КОГДА Таб.К - ТабПред.К = -2
		|				И ТабПред.Разница_Х * Таб.Разница_У >= ТабПред.Разница_У * Таб.Разница_Х
		|			ТОГДА 1
		|		КОГДА Таб.К - ТабПред.К = 2
		|				И НЕ ТабПред.Разница_Х * Таб.Разница_У >= ТабПред.Разница_У * Таб.Разница_Х
		|			ТОГДА -1
		|		ИНАЧЕ 0
		|	КОНЕЦ КАК Результат
		|ПОМЕСТИТЬ Развернуто
		|ИЗ
		|	Таб КАК Таб
		|		ВНУТРЕННЕЕ СОЕДИНЕНИЕ Таб КАК ТабПред
		|		ПО Таб.НаименованиеТочки = ТабПред.НаименованиеТочки
		|			И (Таб.НомерТочкиПолигона = ТабПред.НомерТочкиПолигона + 1
		|				ИЛИ Таб.НомерТочкиПолигона = 1
		|					И ТабПред.НомерТочкиПолигона = &РазмерПолигона)
		|;
		|
		|////////////////////////////////////////////////////////////­////////////////////
		|ВЫБРАТЬ
		|	Развернуто.НаименованиеТочки,
		|	Развернуто.Точка_Х КАК Х,
		|	Развернуто.Точка_У КАК У,
		|	ВЫБОР
		|		КОГДА СУММА(Развернуто.Результат) = 0
		|			ТОГДА ЛОЖЬ
		|		ИНАЧЕ ИСТИНА
		|	КОНЕЦ КАК Результат
		|ИЗ
		|	Развернуто КАК Развернуто
		|
		|СГРУППИРОВАТЬ ПО
		|	Развернуто.НаименованиеТочки,
		|	Развернуто.Точка_Х,
		|	Развернуто.Точка_У";
	Возврат Запрос.Выполнить().Выгрузить();
КонецФункции
Показать


Мой слабенький нетбук на исходных данных считает алгоритмически 93 секунды, а запросом - 27.
frkbvfnjh; nomadon; Дмитрий Рудаков; vadimlp77; smit1c; logos; +6 Ответить
4. Serj1C 475 02.12.12 22:09 Сейчас в теме
(3) Serj1C, Файл с обработкой не могу приложить - в комментариях кнопочки нет(
5. logos 140 03.12.12 13:34 Сейчас в теме
(4) Serj1C, Интересно, хотя это всё полумеры и по делу нужно выносить вычисления в библиотеку.
15. DanilaDru 249 01.11.17 10:52 Сейчас в теме
(3)
| КОГДА СУММА(Развернуто.Результат) = 0


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

Столкнулся с тем что в выборку попадают противоположные относительно друг друга полигоны. Так как у одного из них результат = -1, а у другого 1.
user779781; +1 Ответить
6. CagoBHuK 31 06.12.12 10:01 Сейчас в теме
Все работает для проскости. Но если полигон располагается на северном полюсе и представляет собой плоскость 200км на 200км? Необходимо добавить расчет кривизны земной поверхности.
7. logos 140 06.12.12 10:21 Сейчас в теме
(6) CagoBHuK, Да, безусловно, Вы правы. Но в этом случае для описания поверхности потребуется другая модель определения попадания точки в область. Необходимо будет определять множество полигонов, видимо треугольных, как минимально допустимая плоскость, это множество должно будет описывать все неровности рельефа. Вы же не забыли, что точка может находиться выше или ниже полигона? В этом случае задача вырождается в определение, попадает ли точка в треугольник или нет для каждого треугольника из множества описывающего поверхность, являющуюся нашей областью.
vadimlp77; +1 Ответить
21. kuzz 02.07.19 17:24 Сейчас в теме
(7) Скажите пожалуйста, а каково Ваше мнение по поводу перевода координат на плоскость? Нужно ли этим заниматься если задача стоит в определении принадлежности точек к областям среднестатистического города нашей страны? С точки зрения увеличения точности
8. CagoBHuK 31 06.12.12 14:08 Сейчас в теме
Если говорить о гуглокартах (а это, несомненно, они изображены на спойлере), то вполне хватит простого перевода изгиба поверхности планеты в плоскость через тригонометрическую функцию sin(широта). Почему я клоню в эту сторону? Да потому, что уже решал подобные задачи, а выложить на Инфостарт не додумался. Поэтому я направляю Вас в правильном направлении, чтобы не наступать на те же грабли, что и я.
10. tav13 5 07.08.14 11:59 Сейчас в теме
(8) CagoBHuK, почитав вас сделал обработку перевода координат на плоскость http://infostart.ru/public/295388/
9. SunShinne 615 09.01.14 22:56 Сейчас в теме
Автору огромное спасибо! Однозначно +. Единственное есть два нюанса - во первых в коде закралась ошибка (путаница между латинскими X Y и кириллическими Х У), во вторых алгоритм работает на вхождение именно внутри полигона, но вот если точка на границе полигона то алгоритм выдает отрицательный результат. Но это вопрос спорный, однако отметить стоит. Еще раз спасибо - помогло!
11. frkbvfnjh 556 05.03.15 11:41 Сейчас в теме
По крайней мере с координатами Google карты не работает. Спер алгоритм на JS, там все вхождения четко определил а этот какую то ерунду выдает, я так понял, что координаты еще предварительно нужно обработать обработкой http://infostart.ru/public/295388/ что ли, что бы этот алгоритм заработал?
20. kuzz 02.07.19 17:20 Сейчас в теме
(11) У нас все работает на ура
12. frkbvfnjh 556 05.03.15 11:48 Сейчас в теме
я так понимаю широта это X долгота это Y, или нет?
13. ildarovich 6672 05.03.15 12:18 Сейчас в теме
Я бы попробовал применить в этой задаче готовый объект "ГеографическаяСхема". Там есть метод "ПроверитьПоРасположению". По описанию делает как раз то, что надо.
16. Жолтокнижниг 247 14.11.18 17:31 Сейчас в теме
(13) К сожалению этот вариант работает не корректно
17. vs84 13.01.19 10:50 Сейчас в теме
(16) Проверяли, заработало, но некорректно определяет вхождение? Планировали использовать этот метод для аналогичной задачи (предполагаем, что он определяет с учетом сферичности поверхности), но смутил ваш комментарий.
18. Жолтокнижниг 247 16.01.19 15:51 Сейчас в теме
(17) Сейчас не помню в чем конкретно были проблемы, но нам не подошел этот метод.
Возможно вообще плохо определял попадания, точно не скажу, потратив день на разбор пришли к выводу что не пригодно.

В итоге использую алгоритм статьи и http://infostart.ru/public/295388/ с небольшой шлифовкой. Работает на ура.
22. frkbvfnjh 556 15.07.19 10:24 Сейчас в теме
(18) А с какой небольшой шлифовкой? Может сделаете статью уже с окончательным вариантом алгоритма, который 100% работает, а то сколько не читаю статей в сети, везде потом в конце пишут, что нужно допилить
23. Жолтокнижниг 247 15.07.19 22:32 Сейчас в теме
(22) К сожалению не могу предоставить код, вот примерный алгоритм

Вводная
1. Используется проекция UTM.
2. Используются небольшие геозоны.
Алгоритм
1. Для перевода координат необходимо определить UTM зону (реквизит Zome из обработки https://infostart.ru/public/295388/)
2. Зона определялась, кажется, по медианной долготе (код определения зоны по координатам есть в обработке)
3. Используя алгоритм обработки и зону переводим геозону в полигон
4. По алгоритму из этой статьи определяем вхождение

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

ПС: работу с геоданными поддерживают СУБД(MS и PostgreSQL) при больших объемах данных лучше использовать их.
25. frkbvfnjh 556 16.07.19 06:13 Сейчас в теме
14. Sergius79 03.09.15 07:38 Сейчас в теме
Спасибо, за решение. В благодарность прикладываю эту же функцию для 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);
КонецФункции 
Показать
24. acanta 64 15.07.19 22:56 Сейчас в теме
Можно подробнее про геоданные в СУБД ?
Никогда не слышала об этом.
Оставьте свое сообщение