{"id":468,"date":"2022-10-13T10:24:39","date_gmt":"2022-10-13T14:24:39","guid":{"rendered":"https:\/\/blogs.mathworks.com\/matlab\/?p=468"},"modified":"2022-10-18T08:33:12","modified_gmt":"2022-10-18T12:33:12","slug":"how-well-can-i-draw-circle","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/matlab\/2022\/10\/13\/how-well-can-i-draw-circle\/","title":{"rendered":"How Well Can I Draw a Circle? A MATLAB Adventure"},"content":{"rendered":"
Today's guest blogger is Rob Holt, who works at MathWorks in Natick, Massachusetts. Rob currently serves as the Manager for Biological Sciences at MathWorks, where he is a coordinator and collaborator for biology, biotech, and pharmaceutical applications. Follow more of Rob's antics on <\/span>his Twitter<\/span><\/a>and <\/span>his LinkedIn<\/span><\/a>.<\/span><\/div>

Am I better than the rest of my team?<\/span><\/h2>
When I was looking out my office window the other day, I was thinking about how I can compare my skills to those of my office mates in the most inane, inconsequential fashion. As I constantly strive for perfection (an admittedly <\/span>dangerous<\/span> goal) and elegance, I was looking for some physical feat that would be straight forward to analyze quantitatively. Then it hit me: draw a circle.<\/span><\/div>
The circle has been a symbol of perfection (and I guess unity or whatever) for centuries, but perfect circles remain an elusive goal. Squaring the circle, calculating pi, forming the ideal cookie: all demonstrably impossible. Sure, one can use a compass, but what about on a white board? Surely this becomes the quintessential metric by which to judge human skill.<\/span><\/div>
But how to quantify? In this example, we'll review how to segment a hand-drawn shape on a whiteboard. We will also calculate how round it is. Then we'll compare the results from my teammates to see whose circle reigns supreme. We'll run the example of my own circle, and then apply to the others.<\/span><\/div>

Read and Show an Image<\/span><\/h2>
First one draws a shape on a whiteboard (or in MS Paint or whatever), takes a picture, and shows it. See my original file at <\/span>Rob.jpg<\/span><\/a><\/div>
% read an image and display it<\/span><\/span><\/div><\/div>
imname = fullfile(<\/span>'images'<\/span>,<\/span>'Rob.jpg'<\/span>); <\/span>% images are named after the artist<\/span><\/span><\/div><\/div>
[~,person,~] = fileparts(imname); <\/span>% extract the person's name<\/span><\/span><\/div><\/div>
 <\/div><\/div>
im = imread(imname);<\/span><\/span><\/div><\/div>
imshow(im)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>

Segment the Drawing from the White Space<\/span><\/h2>
This segmentation task is relatively straightforward, since we are trying to segment a dark shape from a white background. But since the background lighting in the office is nonuniform and the white board is somewhat reflective, we'll be sure to use an <\/span>adaptive threshold.<\/span> That is, we'll use a different local threshold value in different parts of the image to separate foreground from background and overcome the effect of the nonuniform lighting and reflections of the lights in the image. This can be done using <\/span>imbinarize<\/span><\/a> with the method set to <\/span>'adaptive'<\/span><\/a>,<\/span><\/div>
% We don't really need color information, so let's make that image gray<\/span><\/span><\/div><\/div>
im = rgb2gray(im);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% The images are super large from phones, and might differ from person to<\/span><\/span><\/div><\/div>
% person, so let's make that image smaller and a little more uniform to make<\/span><\/span><\/div><\/div>
% processing easier<\/span><\/span><\/div><\/div>
im = imresize(im,[500 nan]);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% this imbinarize command uses adaptive thresholding to segment foreground<\/span><\/span><\/div><\/div>
% from background<\/span><\/span><\/div><\/div>
seg = imbinarize(im,<\/span>'adaptive'<\/span>,<\/span>'ForegroundPolarity'<\/span>,<\/span>'dark'<\/span>);<\/span><\/span><\/div><\/div>
seg = ~seg; <\/span>% now the dark part is bright<\/span><\/span><\/div><\/div>
 <\/div><\/div>
%% clean up that segmentation by keeping only the largest connected element<\/span><\/span><\/div><\/div>
seg = bwareafilt(seg,1);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
figure<\/span><\/span><\/div><\/div>
imshow(seg)<\/span><\/span><\/div><\/div>
title(<\/span>'Image Segmentation'<\/span>)<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>

See How Round That Thing Is<\/span><\/h2>
There are many different ways of calculating how good I am at things. For drawing circles, we're going to take a look at two different methods. First, we calculate the <\/span>eccentricity<\/span><\/a> of the shape. The eccentricity of an ideal circle is zero, and higher is worse. <\/span><\/div>
Then we calculate the <\/span>circularity<\/span><\/a>. The circularity of an ideal circle is 1, and lower is worse.<\/span><\/div>
Since circularity relies on having a well-defined perimeter, we will close the shape into a convex hull before running this calculation.<\/span><\/div>
Both eccentricity and circularity can be calculated directly using <\/span>regionprops<\/span><\/a>.<\/span><\/div>
% make a convex hull of the shape<\/span><\/span><\/div><\/div>
segfill = bwconvhull(seg);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% calculate the area and perimeter of this shape<\/span><\/span><\/div><\/div>
convprops = regionprops(segfill,<\/span>'Eccentricity'<\/span>,<\/span>'circularity'<\/span>);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
fprintf(<\/span>'\\tEccentricity is %.3f\\n'<\/span>,convprops.Eccentricity);<\/span><\/span><\/div>
\tEccentricity is 0.296<\/div><\/div><\/div><\/div>
fprintf(<\/span>'\\tCircularity is %.3f\\n'<\/span>,convprops.Circularity);<\/span><\/span><\/div>
\tCircularity is 0.995<\/div><\/div><\/div><\/div><\/div>

See How I Compare to my Teammates<\/span><\/h2>
Now it's time to see how I stack up to my teammates. We've encapsulated most of this processing in a function called <\/span>circleProcess<\/span>.<\/span><\/div>
imdir = <\/span>'images'<\/span>; <\/span>% where the images are stored<\/span><\/span><\/div><\/div>
imnames = dir(fullfile(imdir,<\/span>'*.*g'<\/span>)); <\/span>% the images are png or jpg, but all end in g<\/span><\/span><\/div><\/div>
numpeople = numel(imnames);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% set up some space for us to store the results<\/span><\/span><\/div><\/div>
people = cell(numpeople,1);<\/span><\/span><\/div><\/div>
ecc = zeros(numpeople,1);<\/span><\/span><\/div><\/div>
cir = zeros(numpeople,1);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% loop over the names and get some information<\/span><\/span><\/div><\/div>
for <\/span>somecounter = 1:numpeople<\/span><\/span><\/div><\/div>
<\/span>% process the images<\/span><\/span><\/div><\/div>
imname = fullfile(imdir,imnames(somecounter).name);<\/span><\/span><\/div><\/div>
[~,person,~] = fileparts(imname);<\/span><\/span><\/div><\/div>
[convprops,segfill,im] = circleProcess(imname);<\/span><\/span><\/div><\/div>
<\/span><\/span><\/div><\/div>
<\/span>% store the results<\/span><\/span><\/div><\/div>
people{somecounter} = person;<\/span><\/span><\/div><\/div>
cir(somecounter) = convprops.Circularity;<\/span><\/span><\/div><\/div>
ecc(somecounter) = convprops.Eccentricity;<\/span><\/span><\/div><\/div>
<\/span><\/span><\/div><\/div>
<\/span>% take a look at the segmentations for some quality control<\/span><\/span><\/div><\/div>
subplot(numpeople,2,(somecounter-1)*2+1)<\/span><\/span><\/div><\/div>
imshow(im)<\/span><\/span><\/div><\/div>
title(sprintf(<\/span>'%s''s Drawing'<\/span>,person))<\/span><\/span><\/div><\/div>
subplot(numpeople,2,(somecounter-1)*2+2)<\/span><\/span><\/div><\/div>
imshow(segfill)<\/span><\/span><\/div><\/div>
title(sprintf(<\/span>'%s''s Segmentation'<\/span>,person))<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Now plot the quantitative results so we can compare<\/span><\/div>
figure;<\/span><\/span><\/div><\/div>
subplot(211)<\/span><\/span><\/div><\/div>
% make a bar graph for eccentricity<\/span><\/span><\/div><\/div>
b1 = bar(ecc);<\/span><\/span><\/div><\/div>
% set the x axis to be the people<\/span><\/span><\/div><\/div>
set(gca,<\/span>'xticklabel'<\/span>,people,<\/span>'ylim'<\/span>,[0 0.7])<\/span><\/span><\/div><\/div>
% label the tops of the bars<\/span><\/span><\/div><\/div>
xtips1 = b1(1).XEndPoints;<\/span><\/span><\/div><\/div>
ytips1 = b1(1).YEndPoints;<\/span><\/span><\/div><\/div>
labels1 = string(b1(1).YData);<\/span><\/span><\/div><\/div>
text(xtips1,ytips1,labels1,<\/span>'HorizontalAlignment'<\/span>,<\/span>'center'<\/span>,<\/span>...<\/span><\/span><\/div><\/div>
<\/span>'VerticalAlignment'<\/span>,<\/span>'bottom'<\/span>)<\/span><\/span><\/div><\/div>
title(<\/span>'Eccentricity (Smaller is Better)'<\/span>)<\/span><\/span><\/div><\/div>
 <\/div><\/div>
subplot(212)<\/span><\/span><\/div><\/div>
% do the same thing for circularity<\/span><\/span><\/div><\/div>
b2 = bar(cir);<\/span><\/span><\/div><\/div>
set(gca,<\/span>'xticklabel'<\/span>,people,<\/span>'ylim'<\/span>,[0.8 1.1])<\/span><\/span><\/div><\/div>
xtips2 = b2(1).XEndPoints;<\/span><\/span><\/div><\/div>
ytips2 = b2(1).YEndPoints;<\/span><\/span><\/div><\/div>
labels2 = string(b2(1).YData);<\/span><\/span><\/div><\/div>
text(xtips2,ytips2,labels2,<\/span>'HorizontalAlignment'<\/span>,<\/span>'center'<\/span>,<\/span>...<\/span><\/span><\/div><\/div>
<\/span>'VerticalAlignment'<\/span>,<\/span>'bottom'<\/span>)<\/span><\/span><\/div><\/div>
title(<\/span>'Circularity (Larger is Better)'<\/span>);<\/span><\/span><\/div>
<\/div>
<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
Well . . . looks like on eccentricity I, Rob Holt, was beaten by our summer intern Xin. That's ok. I still have my pride. <\/span><\/div>
Oh, and Temo beat me on circularity. And apparently beat math, because the circularity of a perfect circle is 1 (this is a relic of the digitization of the image, and I suspect the <\/span>coastline paradox<\/span><\/a>). Anyways, my circle is nearly the least circular.<\/span><\/div>
So . . . ummm . . . I . . . no longer wish I did this study.<\/span><\/div>
Hope you learned something. Now on to squares?<\/span><\/div>
Helper functions<\/span><\/div>
function <\/span>[convprops,segfill,im] = circleProcess(imname)<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% read, grayscale, and resice<\/span><\/span><\/div><\/div>
im = imread(imname);<\/span><\/span><\/div><\/div>
im = rgb2gray(im);<\/span><\/span><\/div><\/div>
im = imresize(im,[500 nan]);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% perform the segmentation<\/span><\/span><\/div><\/div>
seg = imbinarize(im,<\/span>'adaptive'<\/span>,<\/span>'ForegroundPolarity'<\/span>,<\/span>'dark'<\/span>);<\/span><\/span><\/div><\/div>
seg = ~seg; <\/span>% now the dark par is bright<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% area filter<\/span><\/span><\/div><\/div>
seg = bwareafilt(seg,1);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% make a convex hull of the shape<\/span><\/span><\/div><\/div>
segfill = bwconvhull(seg);<\/span><\/span><\/div><\/div>
 <\/div><\/div>
% calculate the area and perimeter of this shape<\/span><\/span><\/div><\/div>
convprops = regionprops(segfill,<\/span>'Eccentricity'<\/span>,<\/span>'circularity'<\/span>);<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div><\/div><\/div>\r\n<\/div>