Prawdopodobnie wymaga to trochę skryptów na dowolnej platformie GIS.
Najbardziej efektywną metodą (asymptotycznie) jest przeciąganie linii pionowej: wymaga sortowania krawędzi według ich minimalnych współrzędnych y, a następnie przetwarzania krawędzi od dołu (minimum y) do góry (maksimum y), dla O (e * log ( e)) algorytm, gdy zaangażowane są e zbocza.
Procedura, choć prosta, jest zaskakująco trudna do wykonania we wszystkich przypadkach. Wieloboki mogą być paskudne: mogą mieć zwisające, odłamki, dziury, być rozłączane, mieć zduplikowane wierzchołki, biegi wierzchołków wzdłuż prostych linii i mieć nierozpuszczone granice między dwoma sąsiednimi składnikami. Oto przykład pokazujący wiele z tych cech (i więcej):
Będziemy w szczególności poszukiwać segmentów poziomych o maksymalnej długości leżących całkowicie w obrębie zamknięcia wielokąta. Na przykład eliminuje to zwis pomiędzy x = 20 a x = 40 emanujący z otworu między x = 10 a x = 25. Łatwo jest wtedy pokazać, że co najmniej jeden z poziomych odcinków o maksymalnej długości przecina co najmniej jeden wierzchołek. (Jeśli istnieją rozwiązania nie przecinające się wierzchołki, będą leżeć we wnętrzu pewnego równoległoboku ograniczonego na górze i na dole przez rozwiązań zrobienia przecinają co najmniej jeden wierzchołek. To daje nam środki, aby znaleźć wszystkie rozwiązania).
W związku z tym przeciąganie linii musi zaczynać się od najniższych wierzchołków, a następnie przesuwać się w górę (to znaczy w kierunku wyższych wartości y), aby zatrzymać się na każdym wierzchołku. Na każdym przystanku znajdujemy wszelkie nowe krawędzie wychodzące z tej wysokości; wyeliminuj wszelkie krawędzie kończące się od dołu na tej wysokości (jest to jeden z kluczowych pomysłów: upraszcza algorytm i eliminuje połowę potencjalnego przetwarzania); i ostrożnie przetwarzaj wszystkie krawędzie leżące całkowicie na stałym poziomie (krawędzie poziome).
Weźmy na przykład stan, gdy osiągnięty zostanie poziom y = 10. Od lewej do prawej znajdują się następujące krawędzie:
x.min x.max y.min y.max
[1,] 10 0 0 30
[2,] 10 24 10 20
[3,] 20 24 10 20
[4,] 20 40 10 10
[5,] 40 20 10 10
[6,] 60 0 5 30
[7,] 60 60 5 30
[8,] 60 70 5 20
[9,] 60 70 5 15
[10,] 90 100 10 40
W tej tabeli (x.min, y.min) są współrzędnymi dolnego punktu końcowego krawędzi, a (x.max, y.max) są współrzędnymi jego górnego punktu końcowego. Na tym poziomie (y = 10) pierwsza krawędź jest przechwytywana w jej wnętrzu, druga krawędź jest przechwytywana na jej dole i tak dalej. Niektóre krawędzie kończące się na tym poziomie, takie jak od (10,0) do (10,10), nie są uwzględnione na liście.
Aby określić, gdzie znajdują się punkty wewnętrzne i zewnętrzne, wyobraź sobie, że zaczynasz od skrajnej lewej - oczywiście poza wielokątem - i przesuwasz się poziomo w prawo. Za każdym razem, gdy przekraczamy krawędź, która nie jest pozioma , na przemian przełączamy się z zewnątrz na wnętrze i z tyłu. (To kolejny kluczowy pomysł.) Jednak wszystkie punkty w obrębie dowolnej poziomej krawędzi są określone jako znajdujące się wewnątrz wielokąta, bez względu na wszystko. (Zamknięcie wielokąta zawsze obejmuje jego krawędzie.)
Kontynuując przykład, oto posortowana lista współrzędnych x, gdzie krawędzie inne niż poziome zaczynają się od linii y = 10 lub przecinają ją:
x.array 6.7 10 20 48 60 63.3 65 90
interior 1 0 1 0 1 0 1 0
(Zauważ, że x = 40 nie ma na tej liście.) Wartości interiortablicy oznaczają lewe punkty końcowe segmentów wewnętrznych: 1 oznacza interwał wewnętrzny, 0 interwał zewnętrzny. Zatem pierwsza 1 wskazuje, że przedział od x = 6,7 do x = 10 znajduje się wewnątrz wielokąta. Następne 0 wskazuje, że przedział od x = 10 do x = 20 znajduje się poza wielokątem. I tak dalej: tablica identyfikuje cztery oddzielne przedziały jako wewnątrz wielokąta.
Niektóre z tych przedziałów, takie jak od x = 60 do x = 63,3, nie przecinają żadnych wierzchołków: szybkie sprawdzenie współrzędnych x wszystkich wierzchołków z y = 10 eliminuje takie przedziały.
Podczas skanowania możemy monitorować długości tych przedziałów, zachowując dane dotyczące znalezionych dotychczas przedziałów maksymalnej długości.
Zwróć uwagę na niektóre konsekwencje tego podejścia. Kiedy napotkamy wierzchołek w kształcie litery „v”, powstają dwie krawędzie. Dlatego podczas przełączania występują dwa przełączniki. Te przełączniki się anulują. Każde odwrócone „v” nie jest nawet przetwarzane, ponieważ obie jego krawędzie są eliminowane przed rozpoczęciem skanowania od lewej do prawej. W obu przypadkach taki wierzchołek nie blokuje odcinka poziomego.
Więcej niż dwie krawędzie mogą współdzielić wierzchołek: ilustruje to (10,0), (60,5), (25, 20) i - choć trudno powiedzieć - w (20,10) i (40 , 10). (To dlatego, że dynda idzie (20,10) -> (40,10) -> (40,0) -> (40, -50) -> (40, 10) -> (20, 10) Zauważ, że wierzchołek przy (40,0) znajduje się również wewnątrz innej krawędzi ... to paskudne.) Ten algorytm radzi sobie z tymi sytuacjami w porządku.
Trudna sytuacja została zilustrowana na samym dole: współrzędne x segmentów innych niż poziome
30, 50
Powoduje to, że wszystko po lewej stronie x = 30 jest uważane za zewnętrzne, wszystko między 30 a 50 jest wewnętrzne, a wszystko po 50 ponownie jest zewnętrzne. Wierzchołek przy x = 40 nigdy nie jest nawet brany pod uwagę w tym algorytmie.
Oto jak wygląda wielokąt na końcu skanu. Pokazuję wszystkie wewnętrzne odstępy zawierające wierzchołki w kolorze ciemnoszarym, wszelkie odstępy o maksymalnej długości w kolorze czerwonym i koloruję wierzchołki zgodnie z ich współrzędnymi y. Maksymalny odstęp wynosi 64 jednostki.
Jedynymi wymaganymi obliczeniami geometrycznymi są obliczenia, w których krawędzie przecinają linie poziome: jest to prosta interpolacja liniowa. Obliczenia są również potrzebne do określenia, które zawierają segmenty wewnętrzne wierzchołki: są betweenness ustalenia, łatwo obliczone z kilkoma nierówności. Ta prostota sprawia, że algorytm jest solidny i odpowiedni zarówno dla reprezentacji współrzędnych całkowitych, jak i zmiennoprzecinkowych.
Jeśli współrzędne są geograficzne , wówczas linie poziome są rzeczywiście na okręgach szerokości geograficznej. Ich długości nie są trudne do obliczenia: wystarczy pomnożyć ich długości euklidesowe przez cosinus ich szerokości geograficznej (w modelu sferycznym). Dlatego ten algorytm ładnie dostosowuje się do współrzędnych geograficznych. (Aby poradzić sobie z zawijaniem się wokół studni południka + -180, może być konieczne znalezienie krzywej od bieguna południowego do bieguna północnego, która nie przechodzi przez wielokąt. Po ponownym wyrażeniu wszystkich współrzędnych x jako przemieszczenia poziome względem tego krzywej, ten algorytm poprawnie znajdzie maksymalny segment poziomy).
Poniżej znajduje się Rkod zaimplementowany do wykonywania obliczeń i tworzenia ilustracji.
#
# Plotting functions.
#
points.polygon <- function(p, ...) {
points(p$v, ...)
}
plot.polygon <- function(p, ...) {
apply(p$e, 1, function(e) lines(matrix(e[c("x.min", "x.max", "y.min", "y.max")], ncol=2), ...))
}
expand <- function(bb, e=1) {
a <- matrix(c(e, 0, 0, e), ncol=2)
origin <- apply(bb, 2, mean)
delta <- origin %*% a - origin
t(apply(bb %*% a, 1, function(x) x - delta))
}
#
# Convert polygon to a better data structure.
#
# A polygon class has three attributes:
# v is an array of vertex coordinates "x" and "y" sorted by increasing y;
# e is an array of edges from (x.min, y.min) to (x.max, y.max) with y.max >= y.min, sorted by y.min;
# bb is its rectangular extent (x0,y0), (x1,y1).
#
as.polygon <- function(p) {
#
# p is a list of linestrings, each represented as a sequence of 2-vectors
# with coordinates in columns "x" and "y".
#
f <- function(p) {
g <- function(i) {
v <- p[(i-1):i, ]
v[order(v[, "y"]), ]
}
sapply(2:nrow(p), g)
}
vertices <- do.call(rbind, p)
edges <- t(do.call(cbind, lapply(p, f)))
colnames(edges) <- c("x.min", "x.max", "y.min", "y.max")
#
# Sort by y.min.
#
vertices <- vertices[order(vertices[, "y"]), ]
vertices <- vertices[!duplicated(vertices), ]
edges <- edges[order(edges[, "y.min"]), ]
# Maintaining an extent is useful.
bb <- apply(vertices <- vertices[, c("x","y")], 2, function(z) c(min(z), max(z)))
# Package the output.
l <- list(v=vertices, e=edges, bb=bb); class(l) <- "polygon"
l
}
#
# Compute the maximal horizontal interior segments of a polygon.
#
fetch.x <- function(p) {
#
# Update moves the line from the previous level to a new, higher level, changing the
# state to represent all edges originating or strictly passing through level `y`.
#
update <- function(y) {
if (y > state$level) {
state$level <<- y
#
# Remove edges below the new level from state$current.
#
current <- state$current
current <- current[current[, "y.max"] > y, ]
#
# Adjoin edges at this level.
#
i <- state$i
while (i <= nrow(p$e) && p$e[i, "y.min"] <= y) {
current <- rbind(current, p$e[i, ])
i <- i+1
}
state$i <<- i
#
# Sort the current edges by x-coordinate.
#
x.coord <- function(e, y) {
if (e["y.max"] > e["y.min"]) {
((y - e["y.min"]) * e["x.max"] + (e["y.max"] - y) * e["x.min"]) / (e["y.max"] - e["y.min"])
} else {
min(e["x.min"], e["x.max"])
}
}
if (length(current) > 0) {
x.array <- apply(current, 1, function(e) x.coord(e, y))
i.x <- order(x.array)
current <- current[i.x, ]
x.array <- x.array[i.x]
#
# Scan and mark each interval as interior or exterior.
#
status <- FALSE
interior <- numeric(length(x.array))
for (i in 1:length(x.array)) {
if (current[i, "y.max"] == y) {
interior[i] <- TRUE
} else {
status <- !status
interior[i] <- status
}
}
#
# Simplify the data structure by retaining the last value of `interior`
# within each group of common values of `x.array`.
#
interior <- sapply(split(interior, x.array), function(i) rev(i)[1])
x.array <- sapply(split(x.array, x.array), function(i) i[1])
print(y)
print(current)
print(rbind(x.array, interior))
markers <- c(1, diff(interior))
intervals <- x.array[markers != 0]
#
# Break into a list structure.
#
if (length(intervals) > 1) {
if (length(intervals) %% 2 == 1)
intervals <- intervals[-length(intervals)]
blocks <- 1:length(intervals) - 1
blocks <- blocks - (blocks %% 2)
intervals <- split(intervals, blocks)
} else {
intervals <- list()
}
} else {
intervals <- list()
}
#
# Update the state.
#
state$current <<- current
}
list(y=y, x=intervals)
} # Update()
process <- function(intervals, x, y) {
# intervals is a list of 2-vectors. Each represents the endpoints of
# an interior interval of a polygon.
# x is an array of x-coordinates of vertices.
#
# Retains only the intervals containing at least one vertex.
between <- function(i) {
1 == max(mapply(function(a,b) a && b, i[1] <= x, x <= i[2]))
}
is.good <- lapply(intervals$x, between)
list(y=y, x=intervals$x[unlist(is.good)])
#intervals
}
#
# Group the vertices by common y-coordinate.
#
vertices.x <- split(p$v[, "x"], p$v[, "y"])
vertices.y <- lapply(split(p$v[, "y"], p$v[, "y"]), max)
#
# The "state" is a collection of segments and an index into edges.
# It will updated during the vertical line sweep.
#
state <- list(level=-Inf, current=c(), i=1, x=c(), interior=c())
#
# Sweep vertically from bottom to top, processing the intersection
# as we go.
#
mapply(function(x,y) process(update(y), x, y), vertices.x, vertices.y)
}
scale <- 10
p.raw = list(scale * cbind(x=c(0:10,7,6,0), y=c(3,0,0,-1,-1,-1,0,-0.5,0.75,1,4,1.5,0.5,3)),
scale *cbind(x=c(1,1,2.4,2,4,4,4,4,2,1), y=c(0,1,2,1,1,0,-0.5,1,1,0)),
scale *cbind(x=c(6,7,6,6), y=c(.5,2,3,.5)))
#p.raw = list(cbind(x=c(0,2,1,1/2,0), y=c(0,0,2,1,0)))
#p.raw = list(cbind(x=c(0, 35, 100, 65, 0), y=c(0, 50, 100, 50, 0)))
p <- as.polygon(p.raw)
results <- fetch.x(p)
#
# Find the longest.
#
dx <- matrix(unlist(results["x", ]), nrow=2)
length.max <- max(dx[2,] - dx[1,])
#
# Draw pictures.
#
segment.plot <- function(s, length.max, colors, ...) {
lapply(s$x, function(x) {
col <- ifelse (diff(x) >= length.max, colors[1], colors[2])
lines(x, rep(s$y,2), col=col, ...)
})
}
gray <- "#f0f0f0"
grayer <- "#d0d0d0"
plot(expand(p$bb, 1.1), type="n", xlab="x", ylab="y", main="After the Scan")
sapply(1:length(p.raw), function(i) polygon(p.raw[[i]], col=c(gray, "White", grayer)[i]))
apply(results, 2, function(s) segment.plot(s, length.max, colors=c("Red", "#b8b8a8"), lwd=4))
plot(p, col="Black", lty=3)
points(p, pch=19, col=round(2 + 2*p$v[, "y"]/scale, 0))
points(p, cex=1.25)
