霍夫变换 R PET 包检测线

Question

我有如下图像，我想使用 PET 包的霍夫变换检测图像中的线条。我需要帮助来了解如何从该图像中获取线条。

library("PET", lib.loc="~/R/win-library/3.1")
library("raster", lib.loc="~/R/win-library/3.1")

p=matrix(diag(100), 100)
library(raster)
r <- raster(p)
plot(r)
abc=hough(p)

viewData(list(p, abc$hData), list("Phantom", "Hough transformed phantom"))

我如上所述应用了霍夫变换。我运行最后一行后得到的原始图像和图像如下

关于如何获取线的坐标（来自原始图像）的任何输入？我知道第二个图像右侧窗格中的白点代表这条线。该线是使用极坐标系绘制的。但我不知道如何使用第二张图像来获取原线的坐标

我查看了 PET 包的文档，但发现很难理解 :( 我运行了他们的示例代码，但我不明白

================================================ ================================

我按照用户 NiceE 在 cmets 中给出的建议更新了我的代码，如下所示

library("PET", lib.loc="~/R/win-library/3.1")
library("raster", lib.loc="~/R/win-library/3.1")

#p=matrix(diag(1000), 1000)
p=matrix(rep(0,10000), 100, 100)
# for (i in 1:100)
# {p[i,100-i+1]=1
# }
for (i in 1:100)
{p[i,50]=1
}


# library(raster)
# r <- raster(p)
# plot(r)
abc=hough(p)

maxPoint<-which(abc$hData==max(abc$hData),arr.ind=T)
library(pracma)
a<-cot(maxPoint[1,"row"]*pi/180)
b<-maxPoint[1,"col"]/sin(maxPoint[1,"row"]*pi/180)
a
b
par(pty="s")
par(mfrow=c(1,2))
#image(r, main="org")
image(p,main="original")
image(abc$hData, main="Houghmatrix")

a 和 b 的新值是否正确？我觉得 b 应该是 50（原线与 (0,0) 的垂直距离）。我做错了什么？

我也想知道为什么 abc$hData 有 181 行和 143 列。我可以想象 181 行与 PI 弧度是 180 度有关。但我对 143 列一无所知...

================================================ ========================= 更新 2 如果我更新我的原始矩阵，因为我觉得我得到了奇怪的答案。我得到 a=-57.6 和 b=1786.12。

p=matrix(rep(0,10000), 100, 100)

for (i in 1:100)
{p[80,i]=1
}

Answer 1

一旦你对数据进行了霍夫变换，找到矩阵最大值的索引（假设你只有一行）：

maxPoint<-which(abc$hData==max(abc$hData),arr.ind=T)

如果你知道只有一行，你也可以取所有 maxPoitns 的平均值。在你的情况下，你会得到这个：

    row col
[1,] 137  72

你还需要函数的默认参数

houghParam

这些为您提供 RhoMin、ThethaMin、DeltaMin 和 DeltaRho，这两个变量的增量。有了这些，你可以从矩阵坐标中得到 rho 和 theta。

theta=(maxPoint[1,"row"]-1)*houghParam["DeltaXY1"]+houghParam["XYmin1"]
rho=(maxPoint[1,"col"]-1)*houghParam["DeltaXY2"]+houghParam["XYmin2"]

如果线的方程是y=ax+b，你可以得到a和b使用：

library(pracma)
a<--cot(theta)
b<-(rho)/sin(theta)

另外，在hough 的手册页中，他们声明他们将图像的中心作为 (0,0) 点。

有关数学解释，请查看Wikipedia page of Hough Transform ...

编辑：更改了公式并删除了一些错误信息

总代码为：

library(PET)
library(pracma)

a=matrix(rep(0,10000), 100, 100)
for (i in 1:100)
{a[i,60]=1
}

d=matrix(rep(0,10000), 100, 100)
for (i in 1:100)
{d[60,i]=1
}

e=matrix(diag(100), 100)

getLineHough<-function(p){
abc=hough(p)
#get the brightest point in the hough tranform
maxPoint<-which(abc$hData==max(abc$hData),arr.ind=T)
#if there is only one line, can average the results in case there are several brightest points
maxPoint<-apply(maxPoint,2,mean)

houghParam<-unlist(abc$Header)

theta=(maxPoint[1]-1)*houghParam["DeltaXY1"]+houghParam["XYmin1"]
rho=(maxPoint[2]-1)*houghParam["DeltaXY2"]+houghParam["XYmin2"]

a<--cot(theta)
b<-rho/sin(theta)


par(mfrow=c(1,2))
image(p,main="original")
#add the predicted lines, also have to change the slope and intercept because
#the origin of the plot function is not the center of the image the bottom left corner
if(theta==0){
        abline(v=(rho+50)/100)
} else{
        abline((b+50-a*50)/100,a)
}
image(abc$hData, main="Houghmatrix")
}

getLineHough(a)
getLineHough(d)
getLineHough(e)

Edit2：文档并没有真正说明矩阵第一行的值是什么。由于有 181 行，它应该从 0 而不是 1*houghParam["DeltaXY1"] 开始。相应地更改了代码

Edit3：将代码制作成函数并将预测线添加到绘图中