将线条拟合到图像

Question

我正在尝试根据像素的强度（或颜色）为图像拟合一条线。下图显示了面板 1 中的典型测试图像，其中在面板 2 中手动绘制了一条线。测试图像（矩阵）可以在这里下载：.RData from dropbox。

我想使用 回归分析 来生成类似于面板 2 中手动绘制的线的东西。但是，我不能使用 简单 线性回归，因为，与所有图像一样，x 轴和 y 轴都有误差。

我愿意接受包含相关方程式、链接等的算法描述...而不一定是我可以复制和粘贴的代码。

我想避免的方法

1. 将在不同斜率处绘制的像素的一系列合成二值图像与实际数据图像相关联。例如下面两张图片的相关性会很好，但我想再次避免这种方法。

1. 使用骨架化算法来缩小图像，以便可以使用简单的线性回归。

Answer 1

有趣的是，地震学家处理类似的问题，他们根据震源和接收器之间的距离校正反射数据，过程称为正常移出 (Normal Moveout)。我使用了类似的过程。

一般算法是：

1. 加载图片
2. 定义一系列斜坡进行调查
3. 定义一个的窗口长度
4. 在一系列斜坡上循环，然后...
   • 根据斜率和窗口大小（下图第一行中的灰点）在图像上定义索引位置 (x,y)。
   • 从上面的 x,y 位置索引的原始矩阵构建一个矩阵（下图第二行中的绘图）。
   • 对矩阵求和，然后通过除以求和矩阵的长度来归一化总和。
   • 保存每个总和（循环的每个速度都会有 1 个总和）
5. 与和矢量的最大（或最小）索引对应的速度矢量是图像在当前像素列（下图中的第三行）处的最佳斜率/速度。
6. 沿图像的列执行上述步骤。

下图直观地描述了该算法。

执行上述程序的代码是问题中给出的一列测试数据：

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
}