Generación de bloques poligonales usando líneas de crecimiento

Question

Primero vamos a generar un conjunto de ubicaciones (es decir, los puntos(x,y)), de forma aleatoria. Entonces le asignamos para cada punto de un número que muestra el ángulo de la línea que se está centrada en ese punto. Los criterios para la longitud de las líneas es que, la longitud de las líneas (todos juntos) crece lentamente hasta que la línea se cruza con otra línea (como se muestra a continuación). Luego de que el punto final de la línea de crecimiento serán congelados (no más de crecimiento) para esa línea. El resultado final será un conjunto de áreas cerradas.

Estamos buscando los algoritmos y las ideas con respecto a la implementación de una solución para el problema.

Actualización:
Aquí se añade el resultado de nuestra implementación en Python de la R código desarrollado a continuación por whuber, la aceptación de la respuesta. Su solución matemática es muy agradable, de tal forma que incluso sus Python edición tuvo sólo 186s (~3m) 500 puntos, Wow! Esto es impresionante. Tenga en cuenta que en esta solución no hay ningún requisito adicional para la refinación de las líneas resultantes.

Nota Final:
Con el fin de obtener líneas en lugar de rayos puede reemplazar (en R):

rays <- array(as.vector(cbind(pts,pts+limits*vec)),dim=c(n.points, 2, 2))

con algo como esto (en Python):

pts += limits*vec
n = len(pts)//2
lines = hstack(pts[:n],pts[n:])

Realmente Se Nota Final
Mencionar que para hacer blocks de touching lines (el producto de la post) si usted está interesado en el uso de DCEL no es un simple (con algunos errores) aplicación a {aquí} para Python amantes. También puede seguir {este post} por otros métodos en la generación de bloques (polígonos) de contacto (intersección) líneas.

Answer 1

Ponemos nuestra reciente idea aquí como una respuesta en lugar de las actualizaciones a la pregunta.

Aquí están nuestros esfuerzos hasta ahora como soluciones. Pseudo-código de nuestra solución inicial:

while any(endpoint can grow):
    for all lines:
        for both endpoints of a line:
            b0 <== check if endpoint can grow
            if not b0: next
            grow line towards this endpoint (very small length ~ growth rate)
            b1 <== check if new endpoint is still inside the region
            b2 <== check if grown line is crossing other lines
            if (not b1) | b2:
                reject the growth
                freeze the endpoint
        update line

La idea (pseudo-código) funciona, pero tarda demasiado tiempo {en Python} para la ejecución (dependiendo de la tasa de crecimiento y la resolución de intersección de detección). Por ejemplo, para sólo 30 puntos que tomó cerca de 20 años. Tenga en cuenta que en este enfoque de la tasa de crecimiento debe ser muy pequeño lo contrario, el resultado será totalmente incorrecta.

---

Nota de implementar la idea en Fortran y se encontró que razonablemente ACEPTAR, 500 puntos en ~1.5 minutos, crecimiento 0.0001 donde la dimensión de la región de estudio es 1x1.

---

La siguiente muestra los resultados:

Answer 2

Esto parece que debe ser un O(n^2) algoritmo para n puntos (aunque he sido incapaz de probar esto). Que significa una escala muy mal, y que estamos condenados a largos tiempos de computación con más de un par de miles de puntos. Pero algunas observaciones ayudará a:

1. Cada "dirección" es en realidad dos direcciones, la de un rayo y otro rayo en la dirección opuesta. No le hace ningún daño, y de hecho simplifica el problema, para generalizar a que la entrada consiste de rayos, cada uno según su punto de base, junto con un vector de velocidad. (Esto también permite que el "crecimiento" a tasas que varían de punto a punto).

2. Durante la "creciente" nada sucede hasta que uno de los rayos cumple su intersección con otro rayo. Llamemos a estas intersecciones "eventos" y digamos que el "tiempo de presencia" de un evento es la primera (la más corta) tiempo para cualquier rayo para cubrir el evento. Es necesario sólo para el proceso de los acontecimientos, en lugar de secuencias de la densidad de puntos espaciados a lo largo de cada rayo.

3. Si una creciente ray es bloqueada por otro, en un evento depende de si el otro ya ha sido bloqueado o todavía está creciendo y ha pasado ya ese evento. Por lo tanto, los eventos deben ser procesadas en el orden (en el tiempo) en el que se producen.

Observaciones (2) y (3) son la totalidad del algoritmo! El resto es sólo detalles. Algunos de ellos se refieren a ¿qué hacer cuando dos rayos de satisfacer al mismo tiempo o cuando se encuentran en la cabeza -. Una solución es declarar bloqueados. Aquí una imagen de lo que sucede:

La entrada era de cuatro rayos que se originan en los cuatro puntos de la muestra y el movimiento en las direcciones cardinales Norte, Noreste, Sureste y el Sur (en orden desde el punto 1 hasta el punto 4). El total de los rayos se representan en gris y los segmentos que determinan están representadas en color. (Si es diferente comportamiento en la colisión mutua puntos que se desee, modificar el intersect función en el código de abajo.)

Cuando codificado en R (que será de alrededor de Python, velocidad y más lento que el Fortran), 100 puntos (200 rayos) tiene cinco segundos para el proceso. 500 puntos (1000 rayos) será, por tanto, tomar (500/100)^2 * 5 = 125 segundos = 2 minutos, en comparación con (500/30)^2 * 20 segundos = 1.5 horas para la implementación de Python. Yo esperaría a ver una sustancial aceleración en un Fortran puerto de este algoritmo.

Es una simple cuestión de modificar el algoritmo para la salida de la finita de segmentos de línea que se crea. Para hacer de ellos cerrado polígonos, convertirlos en un doble conectado borde de la lista (DCEL) y el proceso en la forma habitual.

grow <- function(pts, vec) {
  n.points <- dim(pts)[1]

  intersect <- function(x, x.dir, y, y.dir) {
    # Return the intersection of two rays given in point-bearing form.
    # Returns the *times* needed to reach the intersection point.
    if (isTRUE(all.equal(x,y))) return(c(Inf, Inf)) # Same origin -> no intersection
    p <- try(solve(cbind(x.dir, -y.dir), y-x), silent=TRUE)
    if (class(p) == "try-error") { # Collinear directions
      v <- c(-x.dir[2], x.dir[1])
      xy <- sqrt(sum((x-y)*(x-y)))
      if (abs(sum(v * (x-y))) < 1.e-10 * xy && sum(x.dir * y.dir) < 0) {
        # Direct collision
        p <- c(xy, xy) / (sqrt(sum(x.dir*x.dir)) + sqrt(sum(y.dir*y.dir)))
      }
      else {p <- c(Inf, Inf)}
    }
    return(p)
  }
  #
  # Compute the events.
  # Store them as tuples (i, j, i.t, j.t) where i and j index the rays and
  # i.t and j.t are the times taken to reach the event for i and j respectively.
  #
  events <- sapply(1:(n.points-1), function(i) sapply((i+1):n.points, 
            function(j) c(i,j, intersect(pts[i,], vec[i,], 
                                         pts[j,], vec[j,]))))
  events <- matrix(unlist(events), ncol=4, byrow=TRUE)
  colnames(events) <- c("i", "j", "i.t", "j.t")
  #
  # Restrict to actual ray intersections: the times must both be non-negative
  # (and finite) and the two rays should not be coincident.
  #
  events <- events[events[,"i.t"] >= 0 & events[,"j.t"] >= 0 & 
                  (is.finite(events[,"i.t"]) | is.finite(events[,"i.t"])) &
                   events[,"i"] != events[,"j"], ]
  #
  # Sort by the first time to reach the intersection.
  #
  events <- events[order(apply(events[ ,3:4], 1, min)), ]
  #
  # The algorithm:
  # If the "winner" is still growing, stop the loser.
  #
  limits <- rep(Inf, n.points) # The times to hit the first event
  apply(events, 1, function(x) {
    if (x["i.t"] <= x["j.t"] && x["i.t"] <= limits[x["i"]]) {
      limits[x["j"]] <<- min(x["j.t"], limits[x["j"]])
    }
    if (x["j.t"] <= x["i.t"] && x["j.t"] <= limits[x["j"]]) {
      limits[x["i"]] <<- min(x["i.t"], limits[x["i"]])
    }
  })
  #
  # The results are rays (or segments) extending from each base point
  # until the time given in `limits`.  (Because the plotting functions
  # will not know how to handle the rays, truncate them to the longest
  # elapsed time.)
  #
  limits[is.infinite(limits)] <- max(limits[!is.infinite(limits)])
  rays <- array(as.vector(cbind(pts, pts + limits * vec)), 
                dim=c(n.points, 2, 2))
  return(rays)
}
#
# Create data.
# 
set.seed(17)
n.points <- 30
pts <- matrix(runif(n.points*2), ncol=2)
pts <- pts[order(pts[,1]),]
directions <- runif(n.points, 0, 2 * pi)
vec <- t(sapply(directions, function(a) c(cos(a), sin(a))))
pts <- rbind(pts, pts)
vec <- rbind(vec, -vec)
#
# Time the algorithm.
#
system.time(rays <- grow(pts, vec))
#
# Auxiliary plotting functions.
#
vector <- function(pts, vec, length=1) {
  x <- sapply(1:n.points, function(i) {
    c(pts[i, ], pts[i, ] + length * vec[i, ], c(NA, NA))
  })
  return(matrix(x, ncol=2, byrow=TRUE))
}
#
# Plot the data and results.
#
par(mfrow=c(1,1))
c.minmax <- c(-0.5, 1.5)
plot(c.minmax, c.minmax, type="n", asp=1, xlab="x", ylab="y")
lines(vector(pts, vec, 13), col=gray(.9))
colors <- palette(rainbow(n.points))
tmp <- sapply(1:(dim(pts)[1]),  
         function(i) lines(t(rays[i,,]), col=colors[1+(i-1)%%n.points], lwd=2))
points(pts, pch=19, cex=1/2, col=gray(.3))