How Well Can I Draw a Circle? A MATLAB Adventure
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 his Twitterand his LinkedIn.
Am I better than the rest of my team?
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 dangerous 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.
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.
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.
Read and Show an Image
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 Rob.jpg
% read an image and display it
imname = fullfile('images','Rob.jpg'); % images are named after the artist
[~,person,~] = fileparts(imname); % extract the person's name
im = imread(imname);
imshow(im)
Segment the Drawing from the White Space
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 adaptive threshold. 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 imbinarize with the method set to 'adaptive',
% We don't really need color information, so let's make that image gray
im = rgb2gray(im);
% The images are super large from phones, and might differ from person to
% person, so let's make that image smaller and a little more uniform to make
% processing easier
im = imresize(im,[500 nan]);
% this imbinarize command uses adaptive thresholding to segment foreground
% from background
seg = imbinarize(im,'adaptive','ForegroundPolarity','dark');
seg = ~seg; % now the dark part is bright
%% clean up that segmentation by keeping only the largest connected element
seg = bwareafilt(seg,1);
figure
imshow(seg)
title('Image Segmentation')
See How Round That Thing Is
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 eccentricity of the shape. The eccentricity of an ideal circle is zero, and higher is worse.
Since circularity relies on having a well-defined perimeter, we will close the shape into a convex hull before running this calculation.
% make a convex hull of the shape
segfill = bwconvhull(seg);
% calculate the area and perimeter of this shape
convprops = regionprops(segfill,'Eccentricity','circularity');
fprintf('\tEccentricity is %.3f\n',convprops.Eccentricity);
fprintf('\tCircularity is %.3f\n',convprops.Circularity);
See How I Compare to my Teammates
Now it's time to see how I stack up to my teammates. We've encapsulated most of this processing in a function called circleProcess.
imdir = 'images'; % where the images are stored
imnames = dir(fullfile(imdir,'*.*g')); % the images are png or jpg, but all end in g
numpeople = numel(imnames);
% set up some space for us to store the results
people = cell(numpeople,1);
ecc = zeros(numpeople,1);
cir = zeros(numpeople,1);
% loop over the names and get some information
for somecounter = 1:numpeople
% process the images
imname = fullfile(imdir,imnames(somecounter).name);
[~,person,~] = fileparts(imname);
[convprops,segfill,im] = circleProcess(imname);
% store the results
people{somecounter} = person;
cir(somecounter) = convprops.Circularity;
ecc(somecounter) = convprops.Eccentricity;
% take a look at the segmentations for some quality control
subplot(numpeople,2,(somecounter-1)*2+1)
imshow(im)
title(sprintf('%s''s Drawing',person))
subplot(numpeople,2,(somecounter-1)*2+2)
imshow(segfill)
title(sprintf('%s''s Segmentation',person))
end
Now plot the quantitative results so we can compare
figure;
subplot(211)
% make a bar graph for eccentricity
b1 = bar(ecc);
% set the x axis to be the people
set(gca,'xticklabel',people,'ylim',[0 0.7])
% label the tops of the bars
xtips1 = b1(1).XEndPoints;
ytips1 = b1(1).YEndPoints;
labels1 = string(b1(1).YData);
text(xtips1,ytips1,labels1,'HorizontalAlignment','center',...
'VerticalAlignment','bottom')
title('Eccentricity (Smaller is Better)')
subplot(212)
% do the same thing for circularity
b2 = bar(cir);
set(gca,'xticklabel',people,'ylim',[0.8 1.1])
xtips2 = b2(1).XEndPoints;
ytips2 = b2(1).YEndPoints;
labels2 = string(b2(1).YData);
text(xtips2,ytips2,labels2,'HorizontalAlignment','center',...
'VerticalAlignment','bottom')
title('Circularity (Larger is Better)');
Well . . . looks like on eccentricity I, Rob Holt, was beaten by our summer intern Xin. That's ok. I still have my pride.
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 coastline paradox). Anyways, my circle is nearly the least circular.
So . . . ummm . . . I . . . no longer wish I did this study.
Hope you learned something. Now on to squares?
Helper functions
function [convprops,segfill,im] = circleProcess(imname)
% read, grayscale, and resice
im = imread(imname);
im = rgb2gray(im);
im = imresize(im,[500 nan]);
% perform the segmentation
seg = imbinarize(im,'adaptive','ForegroundPolarity','dark');
seg = ~seg; % now the dark par is bright
% area filter
seg = bwareafilt(seg,1);
% make a convex hull of the shape
segfill = bwconvhull(seg);
% calculate the area and perimeter of this shape
convprops = regionprops(segfill,'Eccentricity','circularity');
end
- Category:
- Guest posts,
- image processing
See Also
-
A New View of Our Logo
Blogs
-
Finding the Channels
Blogs
Comments
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.