为图像拟合一条线
Fitting a line to an image
我正在尝试根据像素的强度(或颜色)为图像拟合一条线。下图显示了面板 1 中的典型测试图像,在面板 2 中手动绘制了一条线 。测试图像(矩阵)可以在此处下载:.RData from dropbox 。
我想使用回归分析来生成类似于面板 2 中手动绘制的线的结果。但是,我不能使用简单 线性回归,因为与所有图像一样,x 轴和 y 轴都有误差。
我愿意接受带有相关方程式、链接等的算法描述...但不一定是我可以复制和粘贴的代码。
我想避免的方法
- 将以各种斜率绘制的像素的一系列合成二值图像与实际数据图像相关联。例如下面两张图片的相关性会很好,但同样,我想避免这种方法。
- 使用骨架化算法缩小图像,以便可以使用简单的线性回归。
有趣的是,地震学家处理类似的问题,他们根据震源和接收器之间的距离校正反射数据,过程称为正常移出 (Normal Moveout)。我使用了类似的过程。
一般算法为:
- 载入图片
- 定义一系列斜率进行调查
- 定义一个 window 长度 < 图像列数
- 在一系列斜坡上循环并...
- 根据 window 的斜率和大小定义图像上的索引位置 (x,y)(下图第一行中的灰色点)。
- 根据在上方的 x,y 位置索引的那些原始矩阵构建一个矩阵(下图第二行中的绘图)。
- 对矩阵求和,然后通过除以求和矩阵的长度对总和进行归一化。
- 保存每个总和(你循环的每个速度都会有 1 个总和)
- 对应于和矢量的最大(或最小)索引的速度矢量是当前像素列(下图中的第三行)图像中最好的slope/velocity。
- 沿着图像的列执行上述步骤。
下图直观地描述了该算法。
执行上述过程的代码在问题中给出的测试数据的一列上是:
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
}
我正在尝试根据像素的强度(或颜色)为图像拟合一条线。下图显示了面板 1 中的典型测试图像,在面板 2 中手动绘制了一条线 。测试图像(矩阵)可以在此处下载:.RData from dropbox
我想使用回归分析来生成类似于面板 2 中手动绘制的线的结果。但是,我不能使用简单 线性回归,因为与所有图像一样,x 轴和 y 轴都有误差。
我愿意接受带有相关方程式、链接等的算法描述...但不一定是我可以复制和粘贴的代码。
我想避免的方法
- 将以各种斜率绘制的像素的一系列合成二值图像与实际数据图像相关联。例如下面两张图片的相关性会很好,但同样,我想避免这种方法。
- 使用骨架化算法缩小图像,以便可以使用简单的线性回归。
有趣的是,地震学家处理类似的问题,他们根据震源和接收器之间的距离校正反射数据,过程称为正常移出 (Normal Moveout)。我使用了类似的过程。
一般算法为:
- 载入图片
- 定义一系列斜率进行调查
- 定义一个 window 长度 < 图像列数
- 在一系列斜坡上循环并...
- 根据 window 的斜率和大小定义图像上的索引位置 (x,y)(下图第一行中的灰色点)。
- 根据在上方的 x,y 位置索引的那些原始矩阵构建一个矩阵(下图第二行中的绘图)。
- 对矩阵求和,然后通过除以求和矩阵的长度对总和进行归一化。
- 保存每个总和(你循环的每个速度都会有 1 个总和)
- 对应于和矢量的最大(或最小)索引的速度矢量是当前像素列(下图中的第三行)图像中最好的slope/velocity。
- 沿着图像的列执行上述步骤。
下图直观地描述了该算法。
执行上述过程的代码在问题中给出的测试数据的一列上是:
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
}