Определение простых циклов связного плоского ориентированного графа

Article
01/22/2010

Решение домашнего задания из предыдущего поста.

Предупреждение: много буков.

Хотя в данном случае решение задачи оказывается проще, чем ее описание. В качестве исходных данных имеется таблица map, импортированная из csv, импортированного из ArcView в посте Как я затаскивал в таблицу карту (Рис.2). В ней содержится куча LINESTRING'ов, т.е. куски границ регионов, образующих все вместе тот самый связный плоский граф, простые циклы в котором (т.е. контуры регионов) нам предстоит выделить. На рисунке куски обозначены линиями разного цвета. Все точки внутри куска имеют степень 2, т.е. ветвление границы внутри куска невозможно. Крайние точки могут быть большей степени (ветвление границы), а могут быть тоже степени 2, т.е. между точками ветвления может протянуться цепочка из нескольких LINESTRING'ов.

Рис.1

Несмотря на то, что на каждом LINESTRING'е имеется ориентация, заданная порядком перечисления точек внутри LINESTRING'а, она исходно ни на что не влияет. Каждый LINESTRING, являющийся внутренним, т.е. кусок границы двух смежных регионов, должен проходиться дважды: один раз для одного региона и второй - для его соседа в обратном порядке. Сделаем для каждого LINESTRING'a его копию, которую развернем в противоположном направлении:

if object_id('RegionsBoundaryFragments', 'U') is not null drop table RegionsBoundaryFragments

select * into RegionsBoundaryFragments from map

insert RegionsBoundaryFragments select -id, dbo.ReverseLinestring(g) from map

alter table RegionsBoundaryFragments add region int, ord int

Скрипт 1

В таблицу RegionsBoundaryFragments перенеслись все исходные LINESTRING'и, а также создались их антиподы. Функция ReverseLinestring была определена в предыдущем посте (Скрипт 6). Она переставляет все точки внутри LINESTRING'а в обратном порядке. Обратный LINESTRING будет иметь тот же id, что и оригинал, но с минусом. Также в таблице были подготовлены поля region, которое будет иметь одинаковое значение у всех LINESTRING'ов, относящихся к одному региону, и ord, задающее порядок прохождения LINESTRING'ов в пределах одного региона. Теперь каждый внутренний LINESTRING будет принадлежать строго одному региону и проходиться лишь однажды. То же справедливо и для внешних LINESTRING'ов, составляющих внешнюю границу, если принять за еще один фиктивный регион плоскость, окружающую внешнюю границу. Эта плоскость будет внутренней областью данного региона.

В силу погрешностей имеющейся карты ребра не всегда сходятся в вершину, а могут заканчиваться где-то поблизости, образуя несколько отдельных вершин, которые, на самом деле, надо считать за одну. То есть ребро, не доходящее до вершины на некоторую величину эпсилон, все равно считается ей инцидентным. Не обязательно вершины при этом должны иметь степень 1. Из N сходящихся в окрестность ребер n1, ..., nm могут сходиться в свои локальные вершинки в этой эпсилон-окрестности. Правильный подбор эпсилон имеет критическое значение для корректной работы алгоритма. Возьмем слишком маленький - потеряем какие-то веточки. Возьмем слишком большой - зацепим лишние.

Хорошо было в предыдущей задаче. Там требовалось сориентировать и упорядочить ребра для создания простого цикла (обхода границы региона). Степень каждой вершины = 2, т.е. из окончания ребра выходит только одно следующее ребро и никаких ветвлений. При допущении простых ограничений на отсутствие циклов экзотических форм (напр., квазивосьмерка) достаточно брать то, что ближе лежит, и никакого эпсилон не нужно. Ближайшее ребро будет продолжением с точностью до ориентации. В нынешней ситуации степень вершины может быть и 3, и 4, и, наверное, больше (см. Рис.1).

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

Для подбора эпсилон я вытащил все точки, имеющиеся в LINESTRING'ах и положил их в отдельную таблицу points. Поле id_fragment - это id соответствующего LINESTRING'а, поле id содержит первоначальный порядок следования точек внутри каждого LINESTRING'а, поле reverse_id - обратный.

if object_id('points', 'U') is not null drop table points

with cte(id_fragment, id, point) as (

select m.id, f.id, f.p from map m cross apply dbo.LinestringToPoints(g) f

)

select *, row_number() over(partition by id_fragment order by id desc) reverse_id into points from cte

Скрипт 2

Вспомогательная функция LinestringToPoints вытаскивает из LINESTRING все составляющие точки

и складывает их в том же порядке в табличную переменную.

if object_id('LinestringToPoints', 'TF') is not null drop function LinestringToPoints

create function dbo.LinestringToPoints(@l geometry) returns @t table (id int identity, p geometry) as begin

declare @i int = 1

while @i <= @l.STNumPoints() begin

insert @t(p) values (@l.STPointN(@i))

set @i += 1

end

return

end

Скрипт 3

С помощью таблицы points легко посчитать расстояния между соседними точками в каждом LINESTRING'е.

select t1.id_fragment, t1.id, t2.id, t1.point.ToString(), t1.point.STDistance(t2.point) s

from points t1 join points t2 on t1.id_fragment = t2.id_fragment and t1.id + 1 = t2.id

order by s

Скрипт 4

и соотнести с минимальным расстоянием между начальными и конечными точками при попарном сравнении всех LINESTRING'ов.

if object_id('EdgePointDistances', 'U') is not null drop table EdgePointDistances

with EdgePoints as (

select * from points where id = 1

union all

select * from points where reverse_id = 1

)

select t1.id_fragment fragment1, t2.id_fragment fragment2, t1.id id_point1, t2.id id_point2, t1.point.ToString() point1, t2.point.ToString() point2, t1.point.STDistance(t2.point) s

into EdgePointDistances

from EdgePoints t1 join EdgePoints t2 on t1.id_fragment > t2.id_fragment

--00:01:39 -> 923440 row(s)

Скрипт 5

Когда расстояние между LINESTRING'ами 0, здесь без вариантов - один является продолжением другого. Логично предположить, что конец одного LINESTRING'а считается совпадающим с началом следующего, если расстояние между ними порядка среднего между точками внутри LINESTRING'ов. После вдумчивого разглядывания результатов Скрипт 2 вкупе с select * from EdgePointDistances order by s я выбрал значение эпсилон 10000.

Алгоритм выявления простых циклов тоже очень простой . Берем произвольный LINESTRING из таблицы RegionBoundaries (Скрипт 1) из числа еще не отнесенных ни к какому региону. (Выше отмечалось, что каждый ориентированный LINESTRING может относиться к одному и только одному региону). Создаем новый номер региона и присваиваем его этому LINESTRING'у. Смотрим, какие LINESTRING'и растут из его конечной точки, то есть чьи начальные точки находятся в эпсилон-окрестности конечной точки данного LINESTRING'а. Выбираем самый левый. Если он уже имеет номер региона, останавливаемся - LINESTRING не может относиться к двум разным регионам, следовательно, это тот, с которого мы начали, т.е. мы замкнули цикл. Если нет, повторяем. Проще говоря, всегда идем налево. В смысле, в самом положительном направлении обхода контура (против ч.с.) Предположим, мы достигли развилки - существует несколько LINESTRING'ов, продолжающих данный. Чтобы определить самый левый, возьмем последнее звено данного LINESTRING'а, т.е. вектор, идущий от его предпоследней точки к последней. Возьмем первое звено LINESTRING'а-кандидата на продолжение, т.е. вектор, идущий от его начальной точки ко второй. Перемножим первый вектор на второй векторным образом. Поск. вектора двумерны, результатом векторного произведения будет (0, 0, axby - aybx). Длина равна a cos alpha b sin beta - a sin alpha b cos beta = ab sin (beta - alpha), то есть произведение длин векторов на синус угла от первого вектора до второго. Если его поделить на длины векторов получится синус угла поворота. Изо всех кандидатов надо брать тот, у кого синус больше, потому что он круче всех заворачивает налево.Такой подход исключает ситуации непростых циклов, когда мы сначала обойдем регион 1, потом регион 2, а потом их объединение (см. Рис.1).

Процедура BuildRegionBoundary принимает на входе некоторый кусок границы (LINESTRING) и, отталкиваясь от него, очерчивает полностью границу региона, к которому он относится, согласно вышеприведенному алгоритму.

if object_id('BuildRegionBoundary', 'P') is not null drop proc BuildRegionBoundary

create proc BuildRegionBoundary @fragment int, @eps float as

if (select region from RegionsBoundaryFragments where id = @fragment) is not null return --проверяем, если данный кусок границы уже отнесен к какому-то региону, сразу выходим

declare @region int = isnull((select max(region) from RegionsBoundaryFragments), 0) + 1 --заводим новый код региона, который сейчас будем очерчивать

while 1 = 1 begin

update RegionsBoundaryFragments set region = @region where id = @fragment --присвоили куску границы код региона

update RegionsBoundaryFragments set ord = isnull((select max(ord) from RegionsBoundaryFragments where region = @region), 0) + 1 where id = @fragment --нарастили порядковый номер в пределах данного региона

--Получили первый вектор, от к-го будем мерить угол; предпоследняя -> последняя точка фрагмента @fragment,

--если он находится в своей исходной ориентации, и вторая -> первая, если в противоположной

--противоположная ориентация означает, что @fragment < 0

declare @start bit = case when @fragment > 0 then 0 else 1 end, @direction bit = case when @fragment > 0 then 1 else 0 end

declare @v1 geometry = dbo.GetLinkVector(abs(@fragment), @start, @direction)

--Строим вектора потенциальных продолжений и смотрим синус угла между вектором окончания текущего куска и потенциальным продолжением.

--Из всех возможных ветвей сворачиваем туда, где синус больше.

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

--Если > 1, значит, это последняя точка куска и он будет проходиться в обратном порядке. Вектор продолжения для него - последняя точка -> предпоследняя.

--Временная таблица с кандидатами на продолжение. Все поля, кроме fragment2, для нужд отладки.

declare @t table (point1 varchar(50), fragment2 int, point2 varchar(50), v2 geometry, sinfi float)

delete @t;

with cte as (select *, dbo.GetLinkVector(fragment2, case when fragment2 > 0 then 1 else 0 end, case when fragment2 > 0 then 1 else 0 end) v2 from dbo.FindOutgoingFragments(@fragment, @eps))

--Если след.кусок границы будет проходиться в первоначальном направлении, его вектор направлен от 1-й точке ко 2-й, если в обратном - от последней к предпоследней.

insert @t select *, dbo.SinAngleBetween2Vectors(@v1, v2) from cte

--Отладочная выдача

select @fragment fragment1, point1, @v1.ToString() v1, fragment2, point2, v2.ToString() v2, sinfi from @t order by sinfi desc

select top 1 @fragment = fragment2 from @t order by sinfi desc --в кач-ве след.итерации берем из возможных кусков-продолжений круче всех сворачивающий налево

if (select region from RegionsBoundaryFragments where id = @fragment) is not null break; --если вышли на уже проходившийся кусок, цикл замкнулся - регион готов

end

Скрипт 6

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

--Функция возвращает первое или последнее звено куска границы @fragment в зависимости от @start

--и ориентирует его в зависимости от @direction.

--Если @direction = 1, порядок считается исходным, eсли 0, противоположным.

--Кусок у нас LINESTRING, след., содержит как минимум 2 точки.

if object_id('GetLinkVector', 'FN') is not null drop function GetLinkVector

create function GetLinkVector(@fragment int, @start bit, @direction bit) returns geometry as begin

--Таблица points содержит все точки из всех кусков границы

--Поле id_fragment - кусок границы (LINESTRING)

--id - порядковый номер точки внутри данного LINESTRING

--reverse_id - обратный порядок

set @fragment = abs(@fragment)

declare @p0 geometry, @p1 geometry

if @start = 1 begin

select @p0 = point from points where id = 1 and id_fragment = @fragment

select @p1 = point from points where id = 2 and id_fragment = @fragment

end

else begin

select @p0 = point from points where reverse_id = 2 and id_fragment = @fragment

select @p1 = point from points where reverse_id = 1 and id_fragment = @fragment

end

if @direction = 0 begin

declare @p geometry = @p0; set @p0 = @p1; set @p1 = @p

end

declare @s varchar(max) = 'LINESTRING('

set @s += CAST(cast(@p0.STX as decimal(20, 9)) as varchar(20)) + ' ' + CAST(cast(@p0.STY as decimal(20, 9)) as varchar(20))

set @s += ','

set @s += CAST(cast(@p1.STX as decimal(20, 9)) as varchar(20)) + ' ' + CAST(cast(@p1.STY as decimal(20, 9)) as varchar(20))

set @s += ')'

return geometry::STLineFromText(@s, 0)

end

Скрипт 7

--Функция определяет, крайние точки каких фрагментов попадают в эпсилон-окрестность конечной точки данного куска

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

--По существу имеет значение fragment2, остальные поля - для отладочных нужд

if object_id('FindOutgoingFragments', 'TF') is not null drop function FindOutgoingFragments

create function dbo.FindOutgoingFragments(@fragment int, @eps float)

returns @t table (point1 varchar(50), fragment2 int, point2 varchar(50))

as begin

if @fragment > 0

with cte as (select

case fragment1 when abs(@fragment) then point1 else point2 end point1,

case fragment1 when abs(@fragment) then fragment2 else fragment1 end fragment2,

case fragment1 when abs(@fragment) then id_point2 else id_point1 end id_point2,

case fragment1 when abs(@fragment) then point2 else point1 end point2

from EdgePointDistances where (fragment1 = abs(@fragment) and id_point1 > 1 or fragment2 = abs(@fragment) and id_point2 > 1) and s < @eps)

--abs - потому что если это противоположный фрагмент, у него будет отрицательный id, хотя все расстояния будут такие же

--or и case - потому что таблица EdgePointDistances сделана треугольной для экономии

--point > 1 в EdgePointDistances, потому что при @fragment > 0 кусок проходится в прямом направлении, след., это конечная точка

insert @t select point1,

case id_point2 when 1 then 1 else -1 end * fragment2, --если ближайшая точка потенциально след.куска имеет порядковый номер 1, следующий кусок границы будет проходиться в прямом порядке и номер фрагмента берем положительным. Если > 1, значит, это последняя точка куска и он будет проходиться в обратном порядке. Это будет его антипод с отриц.номером

point2 from cte

else

with cte as (select

case fragment1 when abs(@fragment) then point1 else point2 end point1,

case fragment1 when abs(@fragment) then fragment2 else fragment1 end fragment2,

case fragment1 when abs(@fragment) then id_point2 else id_point1 end id_point2,

case fragment1 when abs(@fragment) then point2 else point1 end point2

from EdgePointDistances where (fragment1 = abs(@fragment) and id_point1 = 1 or fragment2 = abs(@fragment) and id_point2 = 1) and s < @eps)

--id_point = 1 в EdgePointDistances, потому что это кусок проходится в противоположном направлении (@fragment < 0), след., это начальная точка

insert @t select point1,

case id_point2 when 1 then 1 else -1 end * fragment2,

point2 from cte

return

end

Скрипт 8

--Считает синус угла между двумя векторами

if object_id('SinAngleBetween2Vectors', 'FN') is not null drop function SinAngleBetween2Vectors

create function SinAngleBetween2Vectors(@v0 geometry, @v1 geometry) returns float as begin

declare @sinfi float = ((@v0.STPointN(2).STX - @v0.STPointN(1).STX) * (@v1.STPointN(2).STY - @v1.STPointN(1).STY) - (@v0.STPointN(2).STY - @v0.STPointN(1).STY) * ((@v1.STPointN(2).STX - @v1.STPointN(1).STX))) / @v0.STLength() / @v1.STLength()

return @sinfi

end

Скрипт 9

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

update dbo.RegionsBoundaryFragments set region = null

declare @fragment int = (select top 1 id from RegionsBoundaryFragments where region is null) --выбрали кусок произвольно

select @fragment

exec BuildRegionBoundary @fragment = @fragment, @eps = 10000

Скрипт 10

Это оказался фрагмент с id = -1 где-то на Чукотке. Все работает нормально. Чукотка обрисовалась:

Рис.2

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

В прошлом посте была написана функция OrderedBoundaryToPolygon, превращающая упорядоченную последовательность LINESTRING'ов в POLYGON (с внутренней областью). Она брала по порядку LINESTRING'и в соответствии с полем ord и доставала из каждого поочередно его точки в определение многоугольника. В данном посте все точки уже вытащены и сложены в таблицу points (см. Скрипт 2), соответственно, можно заменить вызовы метода STPointN() на чтение этой таблицы.

if object_id('OrderedBoundaryToPolygon', 'FN') is not null drop function OrderedBoundaryToPolygon

create function dbo.OrderedBoundaryToPolygon(@region int) returns geometry as begin

declare @id_fragment int, @s varchar(max) = 'POLYGON((', @ord int = 0, @p geometry

while 1 = 1 begin

select top 1 @id_fragment = id, @ord = ord from RegionsBoundaryFragments where region = @region and ord > @ord order by ord

if @@ROWCOUNT = 0 break

declare @cur cursor

if @id_fragment > 0 set @cur = cursor fast_forward for select point from points where id_fragment = abs(@id_fragment) order by id

else set @cur = cursor fast_forward for select point from points where id_fragment = abs(@id_fragment) order by id desc

open @cur

while 1 = 1 begin

fetch next from @cur into @p

if @@fetch_status <> 0 break

set @s += cast(cast(@p.STX as decimal(20, 9)) as varchar(20)) + ' ' + cast(cast(@p.STY as decimal(20, 9)) as varchar(20)) + ', '

end

close @cur; deallocate @cur

end

select top 1 @id_fragment = id from RegionsBoundaryFragments where region = @region order by ord

if @id_fragment > 0 select top 1 @p = point from points where @id_fragment = abs(id_fragment) order by id

else select top 1 @p = point from points where @id_fragment = abs(id_fragment) order by id desc

set @s += cast(cast(@p.STX as decimal(20, 9)) as varchar(20)) + ' ' + cast(cast(@p.STY as decimal(20, 9)) as varchar(20)) + '))'

return geometry::STGeomFromText(@s, 0)

end

Скрипт 11

Вот искомый многоугольник для Чукотского региона:

Рис.3

Если тот же фрагмент внешней границы теперь пройти в противоположном направлении, будет очерчен фиктивный регион, внутренней областью которого является вся внешняя граница, а внутренней областью - окружающая плоскость:

Рис.4

Извините, это не совсем то, что я ожидал. Надо разобраться.

Позднейшая вставка.

Самую малость не хватило, дабы триумфально завершить сей доблестный труд. Причина ошибки и окончательное решение приводятся здесь.

Алексей Шуленин

Partager via

Определение простых циклов связного плоского ориентированного графа

Ressources supplémentaires