Generating polygonal blocks using growing lines

Question

We first generate a set of locations (i.e., points(x,y)), randomly. Then we assign for each point a number showing the angle of a line being centred at that point. The criteria for the length of lines is that, the length of lines (all together) grows slowly until a line intersects another line (as shown below). Then that end-point of growing line will be frozen (no further growth) for that line. The final result will be a set of closed areas.

We are looking for algorithms and ideas regarding the implementation of a solution for the mentioned problem.

Update:
Here we add the result of our implementation in Python of the R code developed below by whuber, the accepted answer. His mathematical solution is very nice, such that even its Python edition took only 186s (~3m) for 500 points, Wow! This is awesome. Note that in this solution there is no further requirement for refining the resulting lines.

Final Note:
In order to get lines instead of rays you can replace (in R):

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

with something like (in Python):

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

Really Final Note
To mention that to make blocks from touching lines (the product of the post) if you are interested in using of DCEL there is a simple (with some bugs) implementation at {here} for Python lovers. Also you may follow {this post} for other methods in generating blocks (polygons) from touching (intersecting) lines.

Answer 1

This looks like it must be a O(n^2) algorithm for n points (although I have been unable to prove this). That means it will scale poorly and you're doomed to long computation times with more than a few thousand points. But some observations will help:

1. Each "direction" is really two directions, that of a ray and another ray in the opposite direction. It does no harm, and actually simplifies the problem, to generalize it to where the input consists of rays, each one given as its base point together with a velocity vector. (This also allows for the "growth" rates to vary from point point.)

2. During the "growing," nothing happens until one ray meets its intersection with another ray. Let us call these intersections "events" and let us say the "time of occurrence" of an event is the first (shortest) time for any ray to meet the event. It is necessary only to process the events, rather than sequences of densely spaced points along each ray.

3. Whether a growing ray is blocked by another one at an event depends on whether the other one has already been blocked or is still growing itself and has already passed that event. Therefore, the events must be processed in the order (in time) at which they occur.

Observations (2) and (3) are the entire algorithm! The rest is just details. Some of them concern what to do when two rays meet at the same time or when they meet head-on. One solution is to declare them both blocked. Here's a picture of what happens:

The input was four rays originating at the four points shown and moving in the cardinal directions North, Northeast, Southeast, and South (in order from point 1 through point 4). The full rays are plotted in gray and the segments they determine are plotted in color. (If different behavior at mutual collision points is desired, modify the intersect function in the code below.)

When coded in R (which will be around Python's speed and slower than Fortran), 100 points (200 rays) take five seconds to process. 500 points (1000 rays) will therefore take (500/100)^2 * 5 = 125 seconds = 2 minutes, compared to (500/30)^2 * 20 seconds = 1.5 hours for the Python implementation. I would expect to see a substantial speedup in a Fortran port of this algorithm.

It's a simple matter to modify the algorithm to output only the finite line segments it creates. To make them into closed polygons, convert them into a doubly connected edge list (DCEL) and process it in the usual ways.

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))

Answer 2

We put our recent idea here as an answer instead of updates to the question.

Here are our efforts so far as solutions. Pseudo-code of our initial solution:

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

The above idea (pseudo-code) works but takes too long time {in Python} for execution (depending on the growth rate and the resolution of intersection detection). For example, for only 30 points it took about 20s. Note that in this approach the growth rate needs to be very small otherwise the output will be totally incorrect.

---

Note We implemented the idea in Fortran and found it reasonably OK, 500 points in ~1.5 minutes, growth 0.0001 where dimension of the study region is 1x1.

---

The following shows our results: