How to detect circularity more accurately than 4*pi*A/P^2?

Hello,

I have noticed when calculating circularity using the above formula I frequently get values greater than 1. I have done some digging to find out why and I came across this paper that talks about the shortcomings of this formula when used on images of circular objects. I have read it, but aside from a few specific examples they give at the end, I am not really sure how to impliment their strategy. Does anyone know how to do this or know of any other methods of calculating circularity.

Note: I have tried regionprop's eccentricity and it is not adequate either.

Thanks in advanced,

Eric

[Steve Eddins]

See this updated answer regarding circularity values greater than 1.

Answer 1

[John D'Errico]

I remember discussing this exact question recently. In fact, the paper you reference answers the question as to why things fail. Let me see how well I can explain things.

Note the expectation is that all objects will give a smaller ratio than 1, because a circle has the largest area for a given perimeter. The problem is when you work in terms of pixels, the perimeter is a fractally quantized thing. No matter how close you get to a circle, no matter how small are the pixels, the perimeter is still measured in terms of a city block measure, not a smooth curve.

So your formula will fail, it MUST fail, as long as you compute the perimeter in terms of the edges of all those perimeter pixels.

You can do better by making a better measure of the perimeter. For example, if you connect the centers of the perimeter pixels (as suggested in the paper) then using the length of those diagonal lines to measure the perimeter as a polygon, now you will get a better measure of the perimeter. As such, your nice little formula will now be a better measure of the circularity of the object.

So what happens with a unit square? The area of a unit square is 1. The perimeter is 4.

4*pi*1/4^2
ans =
       0.7854

So lets see what happens when I take a circle of radius 4.5, but convert it to a pixelated curve.

A = zeros(15);
[I,J] = meshgrid(1:15,1:15);
k = sqrt((I - 8).^2 + (J - 8).^2) <= 4.5;
A(k) = 1;
A
A =
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     1     1     1     1     0     0     0     0     0
     0     0     0     0     1     1     1     1     1     1     1     0     0     0     0
     0     0     0     1     1     1     1     1     1     1     1     1     0     0     0
     0     0     0     1     1     1     1     1     1     1     1     1     0     0     0
     0     0     0     1     1     1     1     1     1     1     1     1     0     0     0
     0     0     0     1     1     1     1     1     1     1     1     1     0     0     0
     0     0     0     1     1     1     1     1     1     1     1     1     0     0     0
     0     0     0     0     1     1     1     1     1     1     1     0     0     0     0
     0     0     0     0     0     1     1     1     1     1     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0

Now, what is the area, as measured by the pixels we chose?

sum(A(:))
ans =
    69
    
pi*4.5^2
ans =
    63.617

So the pixelated region is only an approximation to the circle I tried to draw. Anyway, lets find the border pixels, and sum the lengths thereof.

fractalperim = sum(sum(abs(conv2(A,[1 -1],'valid')))) + sum(sum(abs(conv2(A,[1;-1],'valid'))))
fractalperim =
    36

So the pixelated fractal perimeter is 36.

4*pi*sum(A(:))/fractalperim^2
ans =
      0.66904

This pixellated approximation to a circle actually has a smaller ratio than I got from a square region.

But, remember that this perimeter really is a fractal thing. The length of that curve is NOT the length of a circle. In fact, no matter how finely you refine that curve, it still has a fractal dimension greater than 1. As such, even a very finely pixelated region will fail, because the fractal length of that perimeter curve is far larger than the corresponding circular arc would have been.

2*pi*r
ans =
       28.274

So the pixelated perimeter overestimated the circular perimeter by over 27%.

36/28.274 - 1
ans =
      0.27325

Yet if we used that more accurate estimate for the perimeter, we do pretty well. (The area is still only an approximation, so we don't get 1 here.)

4*pi*sum(A(:))/(2*pi*r)^2
ans =
       1.0846

Hmm. Can we get a better estimate of the perimeter? Perhaps. Here is my attempt:

So perhaps if I use the coordiantes of the perimeter pixels, then compute a convex hull.

perimpix = abs(conv2(A,[1 -1],'same')) | abs(conv2(A,[1;-1],'same'));
[Ip,Jp] = find(perimpix);

% use a convexhull to generate a polygonal perimeter
edgelist = convhull(Ip,Jp);
polyperim = sum(sqrt(diff(Ip(edgelist)).^2 + diff(Jp(edgelist)).^2))
polyperim =
       30.728

4*pi*sum(A(:))/polyperim.^2
ans =
      0.91832

So this comes much closer to a unit area. By avoiding the fractal perimeter we computed before, things get MUCH better.

What happens if I used a MUCH larger radius circle, so much smaller pixels? Again, I'll use a non-integer radius to avoid edge effects.

A = zeros(200);
[I,J] = meshgrid(1:200,1:200);
r = 89.5;
k = sqrt((I - 100).^2 + (J - 100).^2) <= r;
A(k) = 1;

Areaest = sum(A(:));
fractalperim = sum(sum(abs(conv2(A,[1 -1],'valid')))) + sum(sum(abs(conv2(A,[1;-1],'valid'))));

perimpix = abs(conv2(A,[1 -1],'same')) | abs(conv2(A,[1;-1],'same'));
[Ip,Jp] = find(perimpix);

% use a convexhull to generate a polygonal perimeter
edgelist = convhull(Ip,Jp);
polyperim = sum(sqrt(diff(Ip(edgelist)).^2 + diff(Jp(edgelist)).^2));

4*pi*Areaest/fractalperim^2
ans =
      0.61734

4*pi*Areaest/polyperim^2
ans =
       0.9908

So the much larger circle with the very fine pixels is now coming quite close to 1 for a ratio, yet the fractal perimeter still failed terribly.

But that depends on the convexity of the region. It would fail of course, if the region is not well represented by a convex hull. So non-convex regions will be problematic. But that too can be fixed with some effort. At least you should understand the issue now.

[Eric Chadwick]

Thank you so much for the quick response! That seems to work and I understand it. I tried with a star shape to see what you mean about the convexity and got a circularity of 0.88, so yeah I will have to find a way to deal with that.

I was thinking I could compare the area of the convex hull to the original area and if it increases too much then it is not very circular. The problem there is what to do with that information? For my purposes I am still very much interested in images with low circularity as it has implications for my research, so I want to include all outputs - circular or not.

[John D'Errico]

I admit, I was not impressed with the paper's metric. Mathematical syntax that they never fully seem to define.

You could certainly compare the convex hull area to that of the pixels identified. Significant deviations would certainly indicate the region is not convex, and therefore a convex hull computed perimeter would fail to be accurate.

Perhaps we can come up with a better solution. How about using an alpha shape? They are in MATLAB now. Start with the large circular region.

A = zeros(200);
[I,J] = meshgrid(1:200,1:200);
r = 89.5;
k = sqrt((I - 100).^2 + (J - 100).^2) <= r;
A(k) = 1;

[Ip,Jp] = find(A);
sum(A(:))
ans =
       25185

S = alphaShape(Ip,Jp,2)
S = 
  alphaShape with properties:

             Points: [25185×2 double]
              Alpha: 2
      HoleThreshold: 0
    RegionThreshold: 0

4*pi*area(S)/perimeter(S)^2
ans =
      0.97079

So it tells me this region is pretty close to circular. The nice thing about an alpha shape is it can be robust to non-convex regions. And the nice thing is it is very easy to compute. An alpha of 2 makes complete sense to apply to pixelated regions too.

A = zeros(200);
[I,J] = meshgrid(1:200,1:200);
k = inpolygon(I,J,[50 150 75 150 50 50],[50 50 100 150 150 50]);
A(k) = 1;
pcolor(A)

So a reasonably non-convex domain. Clearly non-circular.

[Ip,Jp] = find(A);
S = alphaShape(Ip,Jp,2);
4*pi*area(S)/perimeter(S)^2
ans =
       0.3372

How about a square region? We expect a ratio of pi/4=0.7854... for a square.

A = zeros(200);
[I,J] = meshgrid(1:200,1:200);
k = inpolygon(I,J,[50 150 150 50 50],[50 50 150 150 50]);
A(k) = 1;

[Ip,Jp] = find(A);
S = alphaShape(Ip,Jp,2);
4*pi*area(S)/perimeter(S)^2
ans =
       0.7854

So the alpha shape seems to work well enough. And it is easy to compute.

I'd go this way if I were you.

[Eric Chadwick]

This is incredible! Thank you very much! I will use definitely use this. Now to learn more about alpha shapes so I can scientifically defend this choice :)

Cheers,

Eric

[Eric Chadwick]

Hi John,

Thanks again for your help last year. I was wondering if you might be able to help me again on a similar problem. I have noticed that alpha shapes "enclose" all points, but they don't surround all points. Is there a way to modify the generation of an alpha shape in order to truel enclose all pixels given?

What I have tried with some success is simply dilating the image that I am extracting the alpha shape from. For my purposes I am using a 2D circle so this works quite well. However, I worry about the robustness if I was ever to move to more complex shapes.

Any ideas? I have attached photos to describe what I have done and what I want.

Cheers,

Eric

[John D'Errico]

As soon as you dilate the result, then by definition it will no longer enclose the points that lie on the perimeter. It will be smoother, but I thought you wanted something that described the points you have.

So there are two solutions. First, you can use an alpha shape with a larger alpha radius. I assume the name defaultalpha means you are allowing alphaShape to choose the value it thinks might be correct.

If you read the help for alphaShape, you will find it allows another argument, an alpha radius of YOUR choosing.

alphaShape(rand(10,1),rand(10,1))
ans = 
  alphaShape with properties:

             Points: [10×2 double]
              Alpha: 0.36111
      HoleThreshold: 0
    RegionThreshold: 0
alphaShape(rand(10,1),rand(10,1),.5)
ans = 
  alphaShape with properties:

             Points: [10×2 double]
              Alpha: 0.5
      HoleThreshold: 0
    RegionThreshold: 0

In the first case, it chose 0.36111. In the second case, I forced it to use a larger value.

As the alpha radius grows to be very large, the alpha shape will approach a convex hull. I don't know if your data ever has concavities in it, that an alpha shape can detect. If not, then you really just wanted a convex hull after all.

An alpha shape can and will result in a non-convex result. In fact, that is what it is designed to allow. A convex hull will ALWAYS be convex. It is the smallest polygonal region that encloses all of your points.

x = rand(10,1);
y = rand(10,1);
S = alphaShape(x,y)
S = 
  alphaShape with properties:

             Points: [10×2 double]
              Alpha: 0.32455
      HoleThreshold: 0
    RegionThreshold: 0

alphaSpectrum(S)
ans =
       4.0277
      0.43071
      0.40049
      0.32455
      0.30012
       0.2948
      0.25527
      0.22377
      0.21632
      0.20858
      0.17767
      0.17281
      0.13751
       0.1286

Sch = alphaShape(x,y,4.03);
plot(Sch,'facecolor','r')
hold on
plot(S,'facecolor','y')

So when I used a value for the alpha radius that exceeds the maximum value from alphaSpectrum(S), I got the convex hull.

In fact, it is often a bad idea to use the default alpha shape. You can also change the value of alpha used in an alpha shape like this:

S.Alpha = 1;

And you may sometimes find it useful to adjust the RegionThreshold or HoldThreshold parameters. But since I don't know what your data normally looks like, that is hard to know.

[Eric Chadwick]

Hi John,

Thanks for the quick reply. My data looks like the image in circle_pix.jpg. I think I was not clear enough in my first query. The default file is with an alpha shape of 2 which produces the closest thing to a circle. Any value greater than 2 yields identical results, so increase the alpha radius doesn't help me enclose the pixels.

When I dilated the image, the idea was to use the apha shape of the dilated image as the alpha shape for the original image because in that case every pixel was contained int he alpha shape as desired. I was just wondering if there is away to achieve that result without dilating the original image.

Cheers,

Eric

[Eric Chadwick]

Nevermind I do not need to do this. Thank you again for your help.

[Steve Eddins]

See also this updated answer regarding circularity values greater than 1.

Answer 2

[Steve Eddins]

R2023a Update: Correcting regionprops Circularity Measurements That Are Greater Than 1.0

Image Processing Toolbox R2023a or Later

In releases R2023a or later, the Circularity computation in regionprops has been revised to correct a bias towards higher values, and it no longer returns values greater than 1.0.

See this 21-Mar-2023 blog post for a detailed explanation.

if verLessThan("images","11.7")
    error("Use this code with releases R2023a or later.")
end
A = imread("text.png");
props = regionprops("table",A,"Circularity");

Image Processing Toolbox R2019a to R2022b

In earlier releases, you can apply the same correction that was introduced in R2023a by adapting the following example. Note that the call to regionprops below uses the form that returns a table because it makes the correction steps easier.

if verLessThan("images","10.4") || ~verLessThan("images","11.7")
    error("Use this code with releases R2019a through R2022b.")
end
A = imread("text.png");
props = regionprops("table",A,["Circularity" "Perimeter"]);
c = props.Circularity;
p = props.Perimeter;
r_eq = p/(2*pi) + 0.5;
correction = (1 - (0.5./r_eq)).^2;
props.Circularity = min(c .* correction, 1);

Image Processing Toolbox R2018b or Earlier

In these older releases, the function regionprops did not compute Circularity. You can compute it by adapting the following example.

if ~verLessThan("images","10.4")
    error("Use this code with releases R2018b or earlier.")
end
A = imread("text.png");
props = regionprops("table",A,["Area" "Perimeter"]);
a = props.Area;
p = props.Perimeter;
c = 4*pi*a ./ (p.^2);
r_eq = p/(2*pi) + 0.5;
correction = (1 - (0.5./r_eq)).^2;
props.Circularity = min(c .* correction, 1);  

[Nikoo Moradpour]

Hello,

Thanks for your comment I have used that for my analysis. Just one question I have is that in the corrected codes you have used perimetere to find the equivalent radius but in this post, it has been mentioned to use the area to find the equivalent radius of an object and then correct the circularity. Would you please correct me if I am wrong?

Thanks