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