Computing the circularity of a binary image

我试图计算给定的二进制图像的圆度。经过一番研究，它清楚地知道，对于圆公式为

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