%% How Long Does My Beard Grow in a Week?
% A mind-bending tale of adventure.  A mildly distasteful yarn.
% 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. He is a coordinator and collaborator for biology, biotech, and pharmaceutical applications. Previously, Rob spent five years as a Senior Scientist at Invicro, a Konica Minolta company, where he designed, implemented, and communicated image analysis algorithms for drug discovery and development. Rob earned his PhD at Dartmouth College, where he focused on molecular cancer imaging through multimodal medical image synthesis. Follow more of Rob’s antics on his Twitterand his LinkedIn.
%     I was trimming my beard last weekend.  While my mind was wandering due to the lack of a nearby computational device, I thought "how fast does my beard hair grow?"  Ir was one of those thoughts that I couldn't un-think.  So I started to plan.
%     I figured that I could let it grow for a month and then measure using a ruler, but then I would probably have to go from clean shaven to disheveled, sacrificing what professionalism I retain.
%     Then inspiration struck.  Trim the beard, let it grow for exactly a week, then trim the beard again to the same length.  If we have a picture of the trimmed hairs, and something in the picture that is of known size, we can extract some kind of length information.  Brilliant!
% Table of Contents
% Read and Show the Image
% Use the Penny to get the Pixel Size
% Segment the Hairs
% Try to Keep Unclumped Hairs
% Estimate the Length of all the Hairs
% Calculate the Average Hair Length
% Summary and Discussion
%
%
%% Read and Show the Image

% load the image
imname = 'mySink.jpg';
im = imread(imname);
figure;
imshow(im);
title('Original Image');

%     Note that the field of view has a penny in it.  The reason I threw that in the sink is that the diameter of a penny is well-defined to be 19.05mm.  This serves as a known reference to determine the pixel size later.  It is also what I could afford.
%
%% Use the Penny to get the Pixel Size
% The function imfindcircles can be used to find circular objects in the field of view.  Note that I have a rough estimate of the radius range, and that since the penny is darker than the background, we set the object polarity to be dark.

% crop the image about the penny
pennyImage = im(1400:1800,1000:1500,:);

% filter the penny image a little
pennyImageFilt = medfilt3(pennyImage,[15 15 1]);

% use the imfindcircles function to find the penny
radrange = [100 300]; %
[centers, radii] = imfindcircles(pennyImageFilt,radrange,...
                             'ObjectPolarity','dark');

% visualize all the circles that are found
figure;
imshow(pennyImage)
viscircles(centers, radii,'EdgeColor','b');
title('Penny Detail With Circle Annotation');

% calculate the pixel size from the radius of the penny
% using mm diameter / pixel diameter
pSize = 19.05/(radii(1)*2); % pixel size in mm
fprintf('The pixel size is %d microns\n',round(pSize*1000))

%
%% Segment the Hairs
% Now that we have the pixel size, we can segment the hairs and try to find their length.  To keep things simpler, let's crop out a section with just hair.

justHair = im(500:1800,1500:3000,:);

figure;
subplot(121)
imshow(justHair);
title('Just Hair');
subplot(122)
imshow(justHair(150:500,900:1350,:));
title('Detail on Hair');

% Now we need to segment the hair from the background.  Since the hair is dark and the background is white, this is fairly straightforward.  This can also be done on a grayscale image since color isn't really adding any information to this task.

grayHair = rgb2gray(justHair);

figure;
subplot(121)
imshow(grayHair);
title('My Gray Hair');
subplot(122)
imshow(grayHair(150:500,900:1350))
title('Detail on Gray Hair');

% Hopefully that's the only gray hair I have in the next few years.
% Now we can use a classic Otsu threshold, which automatically separates an image into two classes based on intensity.  This can be done using the command graythresh, which determines the threshold that separates the two classes by intensity, and imbinarize, which uses that threshold to make a binary image.

thresh = graythresh(grayHair);
hairSeg = ~imbinarize(grayHair,thresh);

figure;
subplot(121)
imshow(hairSeg);
title('Hair Segmentation');
subplot(122)
imshow(hairSeg(150:500,900:1350));
title('Detail on Hair Segmentation')

% In this case, imbinarize makes the brighter part "true", so the ~ is to logically flip that to make the hair part of hairSeg, the segmentation of the hair, to be true.
%
%% Try to Keep Unclumped Hairs
% Now we have all the hairs.  Some of them are in clumps, but most of them look like they're standalone.  For a first pass estimate, let's remove the larger clusters, and also the smaller clusters, leaving mostly unique hairs.

% remove any hairs that get cut off at the border
hairSeg = imclearborder(hairSeg); 

% remove any clusters that are too small
hairSeg = bwareaopen(hairSeg,30); 

% remove clusters that are too large
hairSeg = hairSeg & ~ bwareaopen(hairSeg,300);

figure;
subplot(121)
imshow(hairSeg)
title('Hairs Examined for Length');
subplot(122)
imshow(hairSeg(150:500,900:1350))
title('Detail on Hairs Examined for Length');

%
%% Estimate the Length of all the Hairs
% We have a segmentation of the best hairs for us to examine now.  What we can do to estimate their length is first skeletonize (make each hair at most one pixel wide), and then count the pixels in each skeletonized hair.  This is done using the property "area", but since each hair is one pixel wide, the area is roughly equal to the length of the hair.

% skeletonize the hair
hairSkel = bwskel(hairSeg);

% show the skeleton
figure;
subplot(121)
imshow(hairSkel);
title('Hair Skeletonization');
subplot(122)
imshow(hairSkel(150:500,900:1350));
title('Detail on Hair Skeletonization');

% calculate the length of all hairs
stats = regionprops(hairSkel,'Area');

% get the areas out of a struc
hairLengthInPixels = cat(1,stats.Area);

% transform to real world lengths
hairLength = hairLengthInPixels*pSize;

%
%% Calculate the Average Hair Length
% So now we have a list of the lengths of all the hairs.  Thus we can visualize a histogram as well as calculate the average length.

figure;
histogram(hairLength);
xlabel('Hair Length (mm)');
ylabel('Number of Hairs');

meanHairLength = mean(hairLength);

fprintf('Average length of hair: %.1fmm\n',meanHairLength);

%
%% Summary and Discussion
% So that's it!  My beard hair grows about 2.2mm a week, or about 0.31mm per day.  That is very close to Gillette's estimate of 0.27mm per day.
% Of course, this model isn't perfect.  The length calculation could have been more accurate by fitting a curve to each hair.  The segmentation could have been more accurate if we found a way to separate more of the larger clumps of hairs.  But for most practical purposes, this is going to be an excellent estimate.
% This kind of workflow can also be applied to other image processing schemes.  For example, measuring the length of caterpillars or worms in life sciences studies, or measuring shard length for industrial manufacturing quality control.
% Do you have any strange projects you've been working on lately?  Did you know you could use MATLAB for this kind of problem?  Do you think that a high-resolution picture of my beard hair is inappropriate?  Please share your thoughts and opinions in the comments here.
%
% Copyright 2021 The MathWorks, Inc.