detect line in scatter

Question

I have repeated samples from the following process: Most samples will only contain points that are randomly distributed on a 2-dimensional plane. Sometimes, however, the sample will contain not only randomly distributed points but also some points that will be arranged roughly on a line. I say roughly because the points I observe may not all be on the line exactly but be randomly displaced somewhat. Here's the R code for a little toy example to demonstrate what I mean:

line_active <- TRUE

line_n <- 50
noise_n <- 100
perturbation.sd <- 0.01

# background noise
noise_x <- runif(noise_n)
noise_y <- runif(noise_n)

# line coordinates
start <- list(x=0.2,y=0.1)
end <- list(x=0.8,y=0.9)

if(line_active) {
  line_pos <- runif(line_n)
  line_x <- start$x+line_pos*(end$x - start$x)
      line_y <- start$y+line_pos*(end$y - start$y)

  line_x_perturbed <- line_x+rnorm(length(line_x),0,perturbation.sd)
  line_y_perturbed <- line_y+rnorm(length(line_y),0,perturbation.sd)
} else {
  line_x_perturbed <- c()
  line_y_perturbed <- c()
}

points <- data.frame(x=c(noise_x,line_x_perturbed),y=c(noise_y,line_y_perturbed))

# plot!
plot(points$x,points$y)
points(line_x_perturbed,line_y_perturbed,col='red')

My questions are:

1. How would I go about detecting whether or not there is a line in any sample I'm looking at?
2. How would I go about detecting the slope and intercept of the line if it is there?
3. How much more tricky does this problem get if the noise does not come from some known distribution but has to be estimated?

Answer 1

A Hough Transformation should help. It's mainly used for detecting lines is images, but (x,y) pairs can be considered a sparse representation of a bitmap image.

The idea is for each (x,y) point in the input, compute a list of (slope, intercept) pairs which represent "all" lines that pass through that point. "all" is constrained to some ranges of slope and intercept and granularity limit. If you look at or analyze the resulting collection of (slope, intercept) pairs, you'll see a cluster for each line in the input.

Often the (slope, intercept) pairs are represented as an image or scatter plot with transparency, so that clusters or co-incident values appear as dark spots (or bright spots depending on you color settings).

(Addition after comment from @whuber: There are other possible transformations besides slope/intercept. Angle/intercept puts the first coordinate in a nicer range, and angle/distance-to-origin avoids the intercept bounds issue for near-vertical lines.)

Here's a quick R attempt starting with your code.

#############################
# Hough transform:
#   for each point find slopes and intercepts that go through the line
#############################

# first, set up a grid of intercepts to cycle through
dy <- max(points$y) - min(points$y)
intercepts <- seq(min(points$y) - dy, max(points$y) + dy, dy/50) # the intercept grid for the lines to consider

# a function that takes a point and a grid of intercepts and returns a data frame of slope-intercept pairs
compute_slopes_and_intercepts <- function(x,y,intercepts) {
  data.frame(intercept=intercepts,
             slope=(y-intercepts) / x)
}

# apply the function above to all points
all_slopes_and_intercepts.list <- apply(points,1, function(point) compute_slopes_and_intercepts(point['x'],point['y'],intercepts))
# bind together all resulting data frames 
all_slopes_and_intercepts <- do.call(rbind,all_slopes_and_intercepts.list)


# plot the slope-intercept representation
plot(all_slopes_and_intercepts$intercept, all_slopes_and_intercepts$slope, pch=19,col=rgb(50,50,50,2,maxColorValue=255),ylim=c(-5,5))
# circle the true value
slope <- (end$y - start$y) / (end$x - start$x)
intercept <- start$y - start$x * slope
points(intercept, slope, col='red', cex = 4)

This generates the following plot:

In the plot, the actual slope and intercept of the true line is circled. Alternatively, use ggplot2 and stat_bin2d showing count per bin.

We'll do something similar to ggplot above to find a "best guess" estimate:

# Make a best guess. Bin the data according to a fixed grid and count the number of slope-intercept pairs in each bin
slope_intercepts = all_slopes_and_intercepts
bin_width.slope=0.05
bin_width.intercept=0.05
slope_intercepts$slope.cut <- cut(slope_intercepts$slope,seq(-5,5,by=bin_width.slope))
slope_intercepts$intercept.cut <- cut(slope_intercepts$intercept,seq(-5,5,by=bin_width.intercept))
accumulator <- aggregate(slope ~ slope.cut + intercept.cut, data=slope_intercepts, length)
head(accumulator[order(-accumulator$slope),]) # the best guesses
(best.grid_cell <- accumulator[which.max(accumulator$slope),c('slope.cut','intercept.cut')]) # the best guess

# as the best guess take the mean of slope and intercept in the best grid cell
best.slope_intercepts <- slope_intercepts[slope_intercepts$slope.cut == best.grid_cell$slope.cut & slope_intercepts$intercept.cut == best.grid_cell$intercept.cut,]
(best.guess <- colMeans(best.slope_intercepts[,1:2],na.rm = TRUE))
points(best.guess['intercept'], best.guess['slope'], col='blue', cex = 4)

This could be improved in all sorts of ways, e.g. by running a kernel density estimation on the data and picking the likelihood maximum.