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

조회 수: 62 (최근 30일)
Eric Chadwick
Eric Chadwick 2019년 1월 30일
편집: 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

채택된 답변

John D'Errico
John D'Errico 2019년 1월 30일
편집: John D'Errico 2019년 1월 30일
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.
  댓글 수: 8
Eric Chadwick
Eric Chadwick 2020년 1월 14일
Nevermind I do not need to do this. Thank you again for your help.
Steve Eddins
Steve Eddins 2023년 3월 22일
See also this updated answer regarding circularity values greater than 1.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Steve Eddins
Steve Eddins 2023년 3월 22일
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
Nikoo Moradpour
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

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Bounding Regions에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by