5
$\begingroup$

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?
$\endgroup$
0
7
$\begingroup$

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:

enter image description here

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

enter image description here

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.

$\endgroup$
8
  • $\begingroup$ Thanks, that looks good! Any pointers on how to do this with ready-made packages in R or Python? $\endgroup$ – RoyalTS Jul 29 '14 at 15:34
  • $\begingroup$ I'm not familiar with any, but any you find are probably focused in image processing, anyway. The starting point should be straightforward to do manually: for each point(for a range of plausible intercepts(append (computed slope for line through that point, intercept))). $\endgroup$ – xan Jul 29 '14 at 16:59
  • $\begingroup$ For anyone who finds this, I found this Python tutorial useful: scikit-image.org/docs/0.9.x/auto_examples/… $\endgroup$ – RoyalTS Jul 29 '14 at 21:28
  • $\begingroup$ Man, that's far more help than I was expecting. Thank you so very much! I've got an additional bit of code that will actually spit out a best guess slope and intercept but don't want to edit your answer. $\endgroup$ – RoyalTS Jul 31 '14 at 14:22
  • 2
    $\begingroup$ The Hough transform will work better and exhibit less bias when you measure slope as an angle rather than as a slope (which is the tangent of the angle). The plot of an individual point in the Hough transform should be sinusoidal, as illustrated here, rather than linear. $\endgroup$ – whuber Aug 18 '14 at 14:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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