Search code examples
matlabimage-processinggrayscalethreshold

Counting colonies on a Petri dish

I have a bunch of Petri dishes full of dots which I'd like to count in Matlab. Can this be done reliably and in batches?

E.g. This plate has 352 colonies

enter image description here

I've tried ImageJ but need to do quite a bit of cropping out of the border and get variable results.

Do you have any suggestions?

Solution

  • My approach to this problem is as follows:

    1. Use Hough transform to identify circles corresponding to the Petri dish.
    2. Global thresholding with Otsu's method, restricted to the dish.
    3. Count colonies as the regional maxima of the original image, which are represented in the segmented image.

    This file exchange toolbox provides us with a working circular Hough transform. Things are pretty straightforward from there:

    function [count,colonies,bw] = colony_count(I)

I = rgb2gray(im2double(I)); %# Color-to-gray conversion.
[m,n] = size(I);

%# Uncomment this if you have might have some images with light background
%# and dark colonies. It will invert any that seem that way.
%#if graythresh(I) < 0.5
%#    I = imcomplement(I);
%#end

bw = I > graythresh(I); %# Otsu's method.
radii = 115:1:130; %# Approx. size of plate, narrower range = faster.
h = circle_hough(bw,radii,'same','normalise'); %# Circular HT.
peaks = circle_houghpeaks(h, radii, 'npeaks', 10); %# Pick top 10 circles.

roi = true(m,n);
for peak = peaks
    [x, y] = circlepoints(peak(3)); %# Points on the circle of this radius.
    x = x + peak(1); %# Translate the circle appropriately.
    y = y + peak(2);
    roi = roi & poly2mask(x,y,m,n); %# Cumulative union of all circles.
end

%# Restrict segmentation to dish. The erosion is to make sure no dish pixels
%# are included in the segmentation.
bw = bw & bwmorph(roi,'erode');

%# Colonies are merged in the segmented image. Observing that colonies are 
%# quite bright, we can find a single point per colony by as the regional
%# maxima (the brightest points in the image) which occur in the segmentation.
colonies = imregionalmax(I) & bw;

%# Component labeling with 4-connectivity to avoid merging adjacent colonies.
bwcc = bwconncomp(colonies,4);
count = bwcc.NumObjects;

    We use this code like this:

    I = imread('https://i.sstatic.net/TiLS3.jpg');
[count,colonies,mask] = colony_count(I);

    I have also uploaded the colony_count function on the file exchange. If you have an image which doesn't work but you think should, leave a comment there.

    The count is 359, which I'd say is pretty close. You can inspect the segmentation (mask) and colony markers (colonies) to see where mistakes are made:

    %# Leave out the changes to mask to just see the colony markers.
%# Then you can see why we are getting some false colonies.
R = I; R(mask) = 255; R(colonies) = 0;
G = I; G(mask) = 0; G(colonies) = 255;
B = I; B(mask) = 0; B(colonies) = 0;
RGB = cat(3,R,G,B);
imshow(RGB);