Computing the circularity of a binary image

问题:

I am trying to compute the circularity of a given binary image. After some research its clear to me that the formula for circularity is

4π*area/perimeter^2

Which should range from 0 to 1, 1 being most circular.

Given the binary matrix im

Computing the area is trivial

area = sum(im)

I am computing the perimeter following this rule: 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)

Then i compute circularity like this:

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

However, I get values larger than 1 sometimes and things dont add up. For example:

this image gives me circ=1.155119

and this image gives me circ=1.148728

Any idea whats going on? Shouldn't the values be more like 0.95 and 0.7

回答1:

Your definition for the "binary perimeter" is not a good approximation of the smooth perimeter.

# 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

One worrying feature is that this perimeter is not rotation-invariant: when you rotate a square of side 1 (with sides parallel to the axes) by 45 degrees, its area remains the same, but its perimeter is divided by 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

Here is a slightly better approximation of the perimeter. For each point on the perimeter, look at which points in its 8-neighbourhood are also in the perimeter; if they form a vertical or horizontal segment, the contribution of the pair to the perimeter is 1, if they are in diagonal, the contribution is 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

回答2:

Let me suggest a different algorithm.

You should be able to get the area of the blob with a high degree of accuracy, just by counting the pixels. You can also find the center by taking an average of the coordinates of each inside pixel. Now you can find the radius with sqrt(area/pi). With the radius and the center you can draw a perfect circle that has nearly the same area - count the number of pixels that are part of both the blob and the perfect circle, and divide by the area calculated earlier.

回答3:

If I am right, a perfect circle encloses maximum area for a given perimeter. So if the ratio is 1.0 for a circle, then according to your formula, it will be greater than 1.0 for any other shape.