Fitting a line to an image

Question

I am trying to fit a line to an image based on the intensity (or color) of the pixels. The figure below shows a typical test image in panel 1 with a line manually drawn in panel 2. The test image (matrix) can be downloaded here: .RData from dropbox .

I would like to use a regression analysis to produce something similar to the manually drawn line in panel 2. However, I can not use a simple linear regression because, as with all images, there are errors in both the x and y axes.

I am open to algorithm descriptions with relevant equations, links, etc... and not necessarily code that I can copy and paste.

METHODS I WANT TO AVOID

1. Correlating a series synthetic binary images of pixels drawn at various slopes with the actual data image. For example the correlation of the two images below would be quite good, but again, I want to avoid this method.

2. Using a skeletonization algorithm to reduce the image such that a simple linear regression can be used.

Answer 1

Seismologists, interestingly enough, deal with similar problems where they correct reflection data based on the distance between a seismic source and a receiver with a process known as normal move out (Normal Moveout). I used a similar process.

The general algorithm is:

1. load in the image
2. define a series of slopes to investigate
3. define a window length that is < number of image columns
4. loop over the series of slopes and...
   • define index locations (x,y) over the image based on the slope and the size of the window (gray points in row one of image below).
   • build a matrix from those original matrix indexed at the x,y locations from above (plots in row two of image below).
   • sum the matrix then normalize the sum by dividing by the length of the summed matrix.
   • save the each sum (there will be 1 sum for every velocity you loop over)
5. The velocity vector corresponding to the max (or min) index of the sum vector is the best slope/velocity of the image at that current pixel column (row three in image below).
6. Perform the above steps along the columns of the image.

The algorithm is visually described in the image below.

The code to perform the above procedure is on one column of the test data given in the question is:

load('test.RData')

## INPUTS ## 
img=test

vel.min=1 ## minimum velocity (or slope) to test
vel.max=20 ## max velocity to test
vel.number=100 ## how many velocities to test
win=10 ## size of window to investigate 

## define a time index
ti=nrow(img)/2

## set up a vector to hold the velocity correlation values
vel.corrs <- rep(NA,vel.number)

## define the set of velocities to search over
vels <- seq(vel.min,vel.max,length.out=vel.number)

## define a velocity index
vi=1

while(vi<=length(vels)) {

    ## build a binary matrix with corresponding to the window and velocity
    bin.mat <- matrix(0,ncol=ncol(img),nrow=nrow(img))
    slope.line <- seq(0,ncol(bin.mat)/vels[vi],length.out=ncol(bin.mat))
    bin.mat[(ti-win/2):(ti+win/2),]=1


    ## define the offeset
    offset <- rep(slope.line,each=win+1)

    ## define the indices of array points according to velocity and window
    win.vel.ind <- cbind(which(bin.mat==1,arr.ind=TRUE)[,1]+offset,which(bin.mat==1,arr.ind=TRUE)[,2])

    ## limit the points to the dimensions of the image
    if(any(floor(win.vel.ind[,1]) > nrow(img))){
        win.vel.ind[(which(floor(win.vel.ind[,1])>nrow(img))),]=NA
        ##win.vel.ind <- win.vel.ind[-(which(floor(win.vel.ind[,1])>nrow(img))),]
    }

    ## pluck the values of the image associated with those non-NA indices
    slice <- img[win.vel.ind]

    ## build a matrix of the slice vector with nrow=win+1
    slice.mat <- matrix(slice,nrow=win+1,ncol=ncol(img),byrow=FALSE)

    ## apply a hamming window
    ##ham.mat <- matrix(hamming(win+1),ncol=ncol(slice.mat),nrow=nrow(slice.mat))
    ##slice.ham <- slice.mat*ham.mat

    ## sum this 'slice' and normalize and store
    vel.corrs[vi] <- sum(slice,na.rm=TRUE)/length(na.omit(slice))

    vi=vi+1
}