Сглаживание полигонов на контурной карте?

52

Вот контурная карта, для которой доступны все полигоны уровней.

Давайте спросим, как сгладить полигоны, сохранив все вершины в их точных местоположениях?

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

На самом деле я ищу кусок кода (желательно на Python ), который может выполнять сглаживание 2D-полигонов (любого типа: выпуклый, вогнутый, самопересекающийся и т. Д.), Достаточно безболезненный (забудьте страницы кодов) и точный.

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

введите описание изображения здесь


1)

Scipy.interpolate:

введите описание изображения здесь

Как видите, полученные сплайны (красные) не являются удовлетворительными!

2)

Вот результат, используя код, приведенный здесь . Это не работает хорошо!

введите описание изображения здесь

3)

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

введите описание изображения здесь

Удовлетворяя условию, что сплайн проходит точки:

введите описание изображения здесь

4)

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

введите описание изображения здесь

K = 2 - это катастрофа и поэтому для k> = 4.

5)

Я удалил одну точку в проблемном месте, и полученный сплайн теперь идентичен сплайну. Но это все еще вопрос, почему метод не работает для всех случаев?

введите описание изображения здесь

6)

Хорошее сглаживание для данных whuber может быть следующим (нарисованным программным обеспечением векторной графики), в котором плавно добавляется дополнительная точка (сравните с обновлением)

4):

введите описание изображения здесь

7)

Смотрите результат из Python-версии кода whuber для некоторых иконических фигур:

введите описание изображения здесь
Обратите внимание, что этот метод не работает для полилиний. Для угла полилинии (контур) зеленый это то, что я хочу, но получил красный. Это необходимо учитывать, поскольку контурные карты всегда являются полилиниями, хотя замкнутые полилинии можно рассматривать как многоугольники, как в моих примерах. Также не то, что проблема, возникшая в обновлении 4, еще не решена.

8) [мой последний]

Вот окончательное решение (не идеальное!):

введите описание изображения здесь

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

разработчик
источник
как вы генерируете контуры "многоугольника"? разве они не всегда будут линиями, поскольку контур, пересекающий край матрицы высот, никогда не закроется сам по себе?
фисташковый
Я использовал функцию v.generalize в GRASS, чтобы сгладить контурные линии с приличными результатами, хотя это может занять некоторое время для карт с очень плотными контурами.
фисташковый
@pistachionut Уровни контуров можно считать полилиниями. Я ищу чистый код на первом этапе. Если нет в наличии, то легкий пакет для Python.
Разработчик
Возможно, посмотрите на scipy.org/Cookbook/Interpolation, потому что это звучит так, как будто вы хотите сплайн
PolyGeo
1
Кривая @Pablo Bezier в вашей ссылке хорошо работает для полилиний. Whuber работает почти хорошо для полигонов. Чтобы они вместе могли решить вопрос. Большое спасибо за то, что поделились своими знаниями бесплатно.
Разработчик

Ответы:

37

Большинство методов сплайнования последовательностей чисел сплайн-полигонов. Хитрость заключается в том, чтобы сплайны плавно «закрывались» на конечных точках. Для этого «оберните» вершины вокруг концов. Затем сплайн x- и y-координаты отдельно.

Вот рабочий пример в R. Он использует кубическую splineпроцедуру по умолчанию, доступную в пакете базовой статистики. Для большего контроля, замените почти любую процедуру, которую вы предпочитаете: просто убедитесь, что она проходит через числа (то есть интерполирует их), а не просто используя их в качестве «контрольных точек».

#
# Splining a polygon.
#
#   The rows of 'xy' give coordinates of the boundary vertices, in order.
#   'vertices' is the number of spline vertices to create.
#              (Not all are used: some are clipped from the ends.)
#   'k' is the number of points to wrap around the ends to obtain
#       a smooth periodic spline.
#
#   Returns an array of points. 
# 
spline.poly <- function(xy, vertices, k=3, ...) {
    # Assert: xy is an n by 2 matrix with n >= k.

    # Wrap k vertices around each end.
    n <- dim(xy)[1]
    if (k >= 1) {
        data <- rbind(xy[(n-k+1):n,], xy, xy[1:k, ])
    } else {
        data <- xy
    }

    # Spline the x and y coordinates.
    data.spline <- spline(1:(n+2*k), data[,1], n=vertices, ...)
    x <- data.spline$x
    x1 <- data.spline$y
    x2 <- spline(1:(n+2*k), data[,2], n=vertices, ...)$y

    # Retain only the middle part.
    cbind(x1, x2)[k < x & x <= n+k, ]
}

Чтобы проиллюстрировать его использование, давайте создадим небольшой (но сложный) многоугольник.

#
# Example polygon, randomly generated.
#
set.seed(17)
n.vertices <- 10
theta <- (runif(n.vertices) + 1:n.vertices - 1) * 2 * pi / n.vertices
r <- rgamma(n.vertices, shape=3)
xy <- cbind(cos(theta) * r, sin(theta) * r)

Сплайн это, используя предыдущий код. Чтобы сделать сплайн более гладким, увеличьте количество вершин со 100; чтобы сделать его менее плавным, уменьшите количество вершин.

s <- spline.poly(xy, 100, k=3)

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

plot(s, type="l", lwd=2, col="Gray")
lines(xy, col="Red", lty=2, lwd=2)
points(xy, col="Red", pch=19)
points(s, cex=0.8)
points(s[c(1,dim(s)[1]),], col="Blue", pch=19)

Шлицевый многоугольник

Whuber
источник
5
Хороший ответ. Есть ли способ гарантировать, что контуры не пересекаются в результате сглаживания?
Кирк Куйкендалл
Это хороший вопрос, @Kirk. Я не знаю ни одного метода, гарантирующего не пересечение этой формы сглаживания. (На самом деле, я даже не вижу, как гарантировать, что сглаженная ломаная линия не является самопересекающейся. Хотя это не является большой проблемой для большинства контуров.) Для этого вам нужно будет вернуться к исходному DEM и вместо этого использовать лучший метод для вычисления контуров в первую очередь. (Там являются лучшими методами - они известны уже давно. - но AFAIK некоторые из наиболее популярного GISes не использовать их)
whuber
Во-первых, я все еще работаю над реализацией вашего ответа на Python, но пока не увенчался успехом. Во-вторых, что получится, если вы примените свой метод к квадрату? Вы можете обратиться к тем, кого я нарисовал в этом вопросе.
Разработчик
1
Я принял это как ответ, так как это дает хорошее решение. Несмотря на то, что он не идеален, но он дал мне некоторые идеи для работы, надеюсь, я найду решение, которое удовлетворяет пунктам, которые я упомянул выше в своем вопросе и комментариях. Вы также можете рассмотреть комментарии uuber к вопросу [QC], там есть хорошие приемы. Наконец, я должен сказать, что перевод на python является почти простым с установленным прекрасным пакетом Scipy. Также рассмотрим комментарий Пабло в КК как возможное решение для полилиний, т. Е. Кривых Безье. Всем удачи.
Разработчик
1
видя ваши ответы, я сожалею, что не позаботился о моей математике !!!
Винаян
2

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

Я не рассматриваю это как упражнение по подгонке 2D кривой, а скорее как 3D. Рассматривая данные как трехмерные, мы можем гарантировать, что кривые никогда не пересекаются, и можем использовать информацию из других контуров, чтобы улучшить нашу оценку для текущего.

В следующем извлечении iPython используется кубическая интерполяция, предоставляемая SciPy. Обратите внимание, что значения z, которые я вычерчивал, не важны, поскольку все контуры равноудалены по высоте.

In [1]: %pylab inline
        pylab.rcParams['figure.figsize'] = (10, 10)
        Populating the interactive namespace from numpy and matplotlib

In [2]: import scipy.interpolate as si

        xs = np.array([0.0, 0.0, 4.5, 4.5,
                       0.3, 1.5, 2.3, 3.8, 3.7, 2.3,
                       1.5, 2.2, 2.8, 2.2,
                       2.1, 2.2, 2.3])
        ys = np.array([0.0, 3.0, 3.0, 0.0,
                       1.1, 2.3, 2.5, 2.3, 1.1, 0.5,
                       1.1, 2.1, 1.1, 0.8,
                       1.1, 1.3, 1.1])
        zs = np.array([0,   0,   0,   0,
                       1,   1,   1,   1,   1,   1,
                       2,   2,   2,   2,
                       3,   3,   3])
        pts = np.array([xs, ys]).transpose()

        # set up a grid for us to resample onto
        nx, ny = (100, 100)
        xrange = np.linspace(np.min(xs[zs!=0])-0.1, np.max(xs[zs!=0])+0.1, nx)
        yrange = np.linspace(np.min(ys[zs!=0])-0.1, np.max(ys[zs!=0])+0.1, ny)
        xv, yv = np.meshgrid(xrange, yrange)
        ptv = np.array([xv, yv]).transpose()

        # interpolate over the grid
        out = si.griddata(pts, zs, ptv, method='cubic').transpose()

        def close(vals):
            return np.concatenate((vals, [vals[0]]))

        # plot the results
        levels = [1, 2, 3]
        plt.plot(close(xs[zs==1]), close(ys[zs==1]))
        plt.plot(close(xs[zs==2]), close(ys[zs==2]))
        plt.plot(close(xs[zs==3]), close(ys[zs==3]))
        plt.contour(xrange, yrange, out, levels)
        plt.show()

Кубический интерполированный результат

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

Джилли
источник
Установленные плавные кривые должны оставаться как можно ближе к исходному многоугольнику / полилинии.
Разработчик
1

Я написал почти тот пакет, который вы ищете ... но он был на Perl и был более десяти лет назад: GD :: Polyline . Он использовал 2D кубические кривые Безье и «сгладил» произвольный многоугольник или «полилинию» (тогда меня называли тем, что сейчас обычно называют «LineString»).

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

Вторая часть проста; первая часть была немного искусства. Здесь было прозрение: рассмотреть «сегмент управления» вершинного N: vN. Сегмент управления было три коллинеарных точки: [cNa, vN, cNb]. Центральной точкой была вершина. Наклон этого контрольного сегмента был равен наклону от вершины N-1 до вершины N + 1. Длина левой части этого сегмента составляла 1/3 длины от вершины N-1 до вершины N, а длина правой части этого сегмента составляла 1/3 длины от вершины N до вершины N + 1.

Если исходная кривая была четыре вершины: [v1, v2, v3, v4]то каждая вершина теперь получить сегмент управления в виде: [c2a, v2, c2b]. Соедините их вместе вот так: [v1, c1b, c2a, v2, c2b, c3a, v3, c3b, c4a, v4]и жуйте их по четыре за четыре точки Безье: [v1, c1b, c2a, v2]затем [v2, c2b, c3a, v3], и так далее. Поскольку [c2a, v2, c2b]были коллинеарны, результирующая кривая будет гладкой в ​​каждой вершине.

Таким образом, это также соответствует вашему требованию для параметризации «плотности» кривой: используйте меньшее значение, чем 1/3, для «более жесткой» кривой, большее значение для «более плотной» подгонки. В любом случае результирующая кривая всегда проходит через исходные заданные точки.

Это привело к плавной кривой, которая «описала» исходный полигон. У меня также был какой-то способ «вписать» плавную кривую ... но я не вижу этого в коде CPAN.

Во всяком случае, в настоящее время у меня нет версии, доступной на Python, и у меня нет цифр. НО ... если / когда я перенесу это на Python, я обязательно опубликую здесь.

Дэн Х
источник
Невозможно оценить код Perl, добавьте графику, чтобы продемонстрировать, как он работает, если это возможно.
Разработчик