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.

enter image description here

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.

enter image description here

2 Answers 2


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:

Figure 1

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.

Figure 2, 500 points

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)}
  # 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))
# Create data.
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.
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))
  • 1
    Thanks for such a great solution. Took sometime to get a working edition in Python. But the resulting is just amazing. Two quick notes about the above R code: 1) events[,"i"] != events[,"j"] is not necessary as i and j are always different for a row. 2) the second is.finite(events[,"i.t"]) should be corrected to is.finite(events[,"j.t"])
    – Developer
    Commented Apr 30, 2013 at 11:06
  • and one more thing to be said that when you do such a great job solving problems in a such nice mathematical/programming way, we wonder how much you enjoy doing research when you interested in. Lucky you! We are still amazed with the answer!
    – Developer
    Commented Apr 30, 2013 at 11:11
  • Re the comments: I was aware of (1) but kept it, despite its redundancy, in case previous logical tests were changed. Thanks for catching (2). Fortunately it makes no difference to the outcome. Incidentally, I now believe there may be an algorithm that has expected asymptotic performance better than O(n^2). It would proceed by adding one point (randomly) at a time to a growing configuration of segments and rays maintained in an appropriate data structure. The coding would be rather tricky to carry out, but would pay off if you needed to apply this algorithm to thousands of points.
    – whuber
    Commented Apr 30, 2013 at 14:26
  • Two more comments: 1) about DCEL, we were unable to find a neat code in Python. :( , 2) one post-processing stage is also required after all to merge the two rays as one line (or modifying the code to do so).
    – Developer
    Commented May 2, 2013 at 13:00
  • Creating a DCEL is a standard first exercise in computational geometry--see the Preparata & Shamos text for example-so the lack of an implementation in your favorite platform isn't much of a hurdle. I don't see any reason to merge rays into lines (and some advantages to not doing that), but if it appears necessary, it's an easy post-processing step. The hard part would be implementing the algorithm; it's not straightforward to insert a new ray into an existing configuration.
    – whuber
    Commented May 2, 2013 at 13:19

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:

enter image description here

  • Note that this method is more flexible than mathematical solution. For example, you can have predefined lines as obstacles to control the growth of rays. Each ray can have different speed of growth. Or a limit can be put on the length of growing lines. etc.
    – Developer
    Commented May 5, 2013 at 4:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.