Алгоритм определения точек перегиба для полилинии

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

Каковы были бы шаги высокого уровня для автоматизации генерации этих очков? У меня есть рабочий стол ArcGIS, и я очень удобен с ArcObjects.

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

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

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

Вот рабочий Rкод, чтобы проиллюстрировать эти идеи. Давайте начнем со строки строки, выраженной в виде последовательности координат:

xy <- matrix(c(5,20, 3,18, 2,19, 1.5,16, 5.5,9, 4.5,8, 3.5,12, 2.5,11, 3.5,3, 
               2,3, 2,6, 0,6, 2.5,-4, 4,-5, 6.5,-2, 7.5,-2.5, 7.7,-3.5, 6.5,-8), ncol=2, byrow=TRUE)

Сплайн координаты х и у отдельно для достижения параметризации кривой. (Параметр будет вызван time.)

n <- dim(xy)[1]
fx <- splinefun(1:n, xy[,1], method="natural")
fy <- splinefun(1:n, xy[,2], method="natural")

Интерполируйте сплайны для построения и вычисления:

time <- seq(1,n,length.out=511)
uv <- sapply(time, function(t) c(fx(t), fy(t)))

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

curvature <- function(t, fx, fy) {
  # t is an argument to spline functions fx and fy.
  xp <- fx(t,1); yp <- fy(t,1)            # First derivatives
  xpp <- fx(t,2); ypp <- fy(t,2)          # Second derivatives
  v <- sqrt(xp^2 + yp^2)                  # Speed
  (xp*ypp - yp*xpp) / v^3                 # (Signed) curvature
  # (Left turns have positive curvature; right turns, negative.)
}

kappa <- abs(curvature(time, fx, fy))     # Absolute curvature of the data

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

curvature.zero <- 2*pi / max(range(xy[,1]), range(xy[,2])) # A small threshold
i.col <- 1 + floor(127 * curvature.zero/(curvature.zero + kappa)) 
palette(terrain.colors(max(i.col)))                        # Colors

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

plot(xy, asp=1, xlab="x",ylab="y", type="n")
tmp <- sapply(2:length(kappa), function(i) lines(rbind(uv[,i-1],uv[,i]), lwd=2, col=i.col[i]))
points(t(sapply(time[diff(kappa < curvature.zero/2) != 0], 
       function(t) c(fx(t), fy(t)))), pch=19, col="Black")
points(xy)

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

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

На этой картинке мы видим три шага:

1- Увеличьте линию, используя угол. Я использовал 10 градусов в качестве параметра, и мы использовали сплитлайн. На картинке изогнутая линия находится в начальной фазе.

arcpy.Densify_edit("line" , "ANGLE" , "","",10)
arcpy.SplitLine_management("line" , "line_split")

2- Выберите сегменты, где shape_length не является избыточным. Как мы видим из таблицы, я не выбрал эти избыточные длины. Затем я выбираю их в новый класс пространственных объектов.

arcpy.Select_analysis("line_split" , "line_split_selected")

3- Мы извлекли вершины, расположенные по краям линий, которые являются точками перегиба.

arcpy.FeatureVerticesToPoints_management("line_split_selected" , "line_split_pnt" , "DANGLE")

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

Если мы назовем исходную строку «line_cur», а обобщенную - «line_gen», то мы можем обрезать «line_cur» по «line_gen». Результатом будет прямой сегмент "line_cur". Затем мы можем очистить некоторые очень короткие сегменты, удалив их с помощью SQL-запроса, который выбирает Shape_length больше минимальной длины дороги.