61 views (last 30 days)

Hi all,

I am using the output of bwconncomp as an input to regionprops to calculate some values for blobs in a binary image. One of the values that I calculate is "solidity" which is a direct output of regionprops. However, I am also interested in calculating circularity, which I define as 4*pi*Area/Perimeter^2. However, I keep getting values that are greater than one, which is unexpected. From running bwconncomp and regionprops on some artificial binary images with some pixels arbitrarily set to 1 to represent particles I can tell that regionprops does not calculate perimeter in the way that I would normally think of it. Does anybody know if the perimeter calculation fails for very small clusters of "on" pixels, or is there some rationale behind what regionprops is doing, and I just don't understand it?

I'll post the code I use to troubleshoot this:

%make an empty matrix

a = zeros(5,5)

%put a small cross shaped particle onto it

a(2, 2:4) = 1

a(1:3, 3) = 1

%call bwconncomp

cc = bwconncomp(a)

%call regionprops

blobFacts = regionprops(cc, 'Perimeter', 'Solidity', 'Area')

%report values

p = blobFacts.Perimeter

s = blobFacts.Solidity

a = blobFacts.Area

%manually calculate circularity

c = (4*pi*a)/p^2

This gives a circularity of 1.9865, which is higher than 1 (what I intuitively think of as the max possible circularity). The area of the particle is 5, which is expected, but the perimeter is 5.6240, and I can't seem to figure out how regionprops is getting that value. Thanks in advance,

Clay

nsunjaya
on 25 Sep 2017

The 'perimeter' calculation in "regionprops" is done based on a formula explained in the following paper:

A.M. Vossepoel, and A.W.M. Smeulders, "Vector code probability and metrication error in the representation of straight lines of finite length", Computer Graphics and Image Processing, Volume 20, Issue 4, December 1982, Pages 347–364.

You can also refer to the following blog post which gives a general overview of the algorithm used (please note that this is not an official documentation):

You can see the perimeter is computed using a local function called "computePerimeterFromBoundary" in the "regionprops" function which corresponds to the formula discussed in the paper and blog post above:

perimeter = sum(isEven)*0.980 + sum(~isEven)*1.406 - sum(isCorner)*0.091;

Matt J
on 23 Sep 2017

Edited: Matt J
on 23 Sep 2017

I think the example of a 2x2 square below is easier for discussion purposes. I assume you expect the calculated perimeter to be 8, since that is the theoretical perimeter of an ideal square of area 4. That would be the calculated result if regionprops traversed the boundary pixel edges. For most shapes, that will give a very poor approximation to the perimeter of the imaged object, so I think it's clear why regionprops doesn't do that.

Based on the documentation for regionprops, I would expect the calculated perimeter to be 4. In the doc, it sounds like the calculation traverses the boundary pixels center-to-center, accumulating the inter-pixel distance. Why it comes out to be 3.5560, I just don't know. However, even a computed perimeter of 4 will yield circularity=pi>1 in your formula. That is, unless you modify the area calculation so that only the area of the region bounded by the pixel centers (in this case 1) is computed. You could do that with polyarea().

Bottom line: the way regionprops computes area is not consistent with the way it computes perimeter, and so your formula for circularity will indeed break down for small pixel clusters.

>> a

a =

0 0 0 0 0

0 1 1 0 0

0 1 1 0 0

0 0 0 0 0

0 0 0 0 0

>> regionprops(a,'Area','Perimeter')

ans =

struct with fields:

Area: 4

Perimeter: 3.5560

Image Analyst
on 24 Sep 2017

If you edit regionprops() you can see how it computes perimeter:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function perimeter = computePerimeterFromBoundaryOld(B)

delta = diff(B).^2;

perimeter = sum(sqrt(sum(delta,2)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function perimeter = computePerimeterFromBoundary(B)

delta = diff(B).^2;

if(size(delta,1) > 1)

isCorner = any(diff([delta;delta(1,:)]),2); % Count corners.

isEven = any(~delta,2);

perimeter = sum(isEven)*0.980 + sum(~isEven)*1.406 - sum(isCorner)*0.091;

else

perimeter = 0; % if the number of pixels is 1 or less.

end

Image Analyst
on 24 Sep 2017

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.