计算二进制图像的圆度

Question

我试图计算给定二进制图像的圆度.经过一些研究后,我明白循环的公式是

4?*area/perimeter^2

其范围应为0到1,1更多为圆形.

给定二进制矩阵 im

计算该区域是微不足道的

area = sum(im)

我按照这条规则计算周长: A pixel is part of the perimeter if it is nonzero and it is connected to at least one zero-valued pixel

per = matrix(0, nrow(im), ncol(im))
for(i in 2:(nrow(im)-1)){
  for(j in 2:(ncol(im)-1)){
    if(im[i,j] != 0){
      x=c(im[i-1,j],im[i+1,j],im[i,j-1], im[i,j+1])
      if(0 %in% x) per[i,j] = 1
    }
  }
}
perimeter = sum(per)

然后我像这样计算圆度:

circ = (4*pi*area)/(perimiter^2)

但是,我有时会得到大于1的值,而事情并没有加起来.例如:

这张图片给了我 circ=1.155119

这个形象给了我 circ=1.148728

有什么想法吗？价值不应该更像0.95和0.7

Answer 1

您对"二进制周长"的定义并不是光滑周长的良好近似值.

# Sample data
n <- 100
im <- matrix(0, 3*n, 3*n+1)
x <- ( col(im) - 1.5*n ) / n
y <- ( row(im) - 1.5*n ) / n
im[ x^2 + y^2 <= 1 ] <- 1
image(im)

# Shift the image in one direction
s1 <- function(z) cbind(rep(0,nrow(z)), z[,-ncol(z)] )
s2 <- function(z) cbind(z[,-1], rep(0,nrow(z)) )
s3 <- function(z) rbind(rep(0,ncol(z)), z[-nrow(z),] )
s4 <- function(z) rbind(z[-1,], rep(0,ncol(z)) )

# Area, perimeter and circularity
area <- function(z) sum(z)
perimeter <- function(z) sum( z != 0 & s1(z)*s2(z)*s3(z)*s4(z) == 0)
circularity <- function(z) 4*pi*area(z) / perimeter(z)^2

circularity(im)
# [1] 1.241127

area(im)
# [1] 31417
n^2*pi
# [1] 31415.93

perimeter(im)
# [1] 564
2*pi*n
# [1] 628.3185

一个令人担忧的特征是这个周长不是旋转不变的:当你将第1边的正方形(边与轴平行)旋转45度时,它的面积保持不变,但是它的周长除以sqrt(2). ..

square1 <- -1 <= x & x <= 1 & -1 <= y & y <= 1
c( perimeter(square1), area(square1) )
# [1]   800 40401

square2 <- abs(x) + abs(y) <= sqrt(2)
c( perimeter(square2), area(square2) )
# [1]   564 40045

这是一个稍微好一些的周长近似值.对于周边的每个点,查看其8个邻域中的哪些点也在周边; 如果它们形成垂直或水平段,则该对对周长的贡献为1,如果它们是对角线,则贡献为sqrt(2).

edge <- function(z) z & !(s1(z)&s2(z)&s3(z)&s4(z))
perimeter <- function(z) {
  e <- edge(z)
  ( 
    # horizontal and vertical segments
    sum( e & s1(e) ) + sum( e & s2(e) ) + sum( e & s3(e) ) + sum( e & s4(e) ) + 
    # diagonal segments
    sqrt(2)*( sum(e & s1(s3(e))) + sum(e & s1(s4(e))) + sum(e & s2(s3(e))) + sum(e & s2(s4(e))) )
  ) / 2  # Each segment was counted twice, once for each end
}

perimeter(im)
# [1] 661.7544
c( perimeter(square1), area(square1) )
# [1]   805.6569 40401.0000
c( perimeter(square2), area(square2) )
# [1]   797.6164 40045.0000

circularity(im)
# [1] 0.9015315
circularity(square1)
# [1] 0.7821711
circularity(square2)
# [1] 0.7909881

Answer 2

让我建议一个不同的算法.

您应该能够通过计算像素来获得高精度的斑点区域.您还可以通过获取每个内部像素的平均值来找到中心.现在你可以找到半径sqrt(area/pi).使用半径和中心,您可以绘制一个具有几乎相同面积的完美圆 - 计算斑点和完美圆的一部分像素数,并除以之前计算的面积.