r - 计算二进制图像的圆度

Question

我正在尝试计算给定二进制图像的圆度。经过一番研究，我很清楚圆度的公式是

4π*area/perimeter^2

范围从0到1，最圆形的为1。

给定二进制矩阵im
计算面积是微不足道的

面积=总和(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