How to detect circularity more accurately than 4piA/P^2?

Question

Eric Chadwick 2019-1-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/442374-how-to-detect-circularity-more-accurately-than-4-pi-a-p-2

编辑： Nikoo Moradpour 2023-8-31

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Steve Eddins 2023-3-22

See this updated answer regarding circularity values greater than 1.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2019-1-30

4
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/442374-how-to-detect-circularity-more-accurately-than-4-pi-a-p-2#answer_358796

编辑：John D'Errico 2019-1-30

在 MATLAB Online 中打开

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     1     1     1     1     1     0     0     0     0     0
   0     0     0     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     1     1     1     1     1     1     1     1     1     0     0     0
   0     0     1     1     1     1     1     1     1     1     1     0     0     0
   0     0     1     1     1     1     1     1     1     1     1     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     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

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.

8 个评论
显示 6更早的评论隐藏 6更早的评论

John D'Errico 2019-1-31

在 MATLAB Online 中打开

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.

John D'Errico 2020-1-14

在 MATLAB Online 中打开

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 2020-1-14

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

Steve Eddins 2023-3-22

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

请先登录，再进行评论。

Answer 2

Steve Eddins 2023-3-22

3
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/442374-how-to-detect-circularity-more-accurately-than-4-pi-a-p-2#answer_1198565

在 MATLAB Online 中打开

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);  

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Nikoo Moradpour 2023-8-31

编辑：Nikoo Moradpour 2023-8-31

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

请先登录，再进行评论。

How to detect circularity more accurately than 4piA/P^2?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

8 个评论
显示 6更早的评论隐藏 6更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

8 个评论 显示 6更早的评论隐藏 6更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

How to detect circularity more accurately than 4piA/P^2?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论