0

I am currently doing a project on morphology of filamentous fungi during batch fermentation (Yes, I am not a software engineer.. Biotech). Where I am taken pictures of the morphology in a petri dish. I am developing a "fast" method to describe the pellets (small aggregates of fungi) that occurs during the fermentation. To do this I am writing a code in MatLab. Depending on the color of the pellets (light or dark) the pictures are taken on differen backgrounds, black or white. I am inverting the picture if the mean gray value is below 70 to distinguish between backgrounds.

Pictures: White background Dark background

I have several problems:

  1. Detecting the edge of the petri dish so it won't be regarded as an object (Currently done with the edge('log',) function). The edge is detected, but i miss some parts, think because of the lower light in top.
  2. Proper thresholding inside the dish
  3. Detection of pellets - right now it is done by a combination of running through each color channel, but might be done with some blob detection?

Does anybody have some inputs?

My code is as following:

close all
clear all
clc

%Empty arrays to hold data
metricD=[];
areaD=[];
perimeterD=[];

% Specify the folder where the files live.
myFolder = pwd;
% Check to make sure that folder actually exists.  Warn user if it doesn't.
if ~isdir(myFolder)
  errorMessage = sprintf('Error: The following folder does not exist:\n%s', myFolder);
  uiwait(warndlg(errorMessage));
  return;
end
% Get a list of all files in the folder with the desired file name pattern.
filePattern = fullfile(myFolder, '*.jpg'); % Change to whatever pattern you need.
theFiles = dir(filePattern);

% Show debugging plots
plotFig = 0;

% parameters that can be tuned
% how many colors channels we minimum want to see a spore in
% e.g. set to 1 for image "P. f Def C.tif"
labelcutOff = 1;

% remove areas larger than
removeLargerthan = 500000;

for k = 1 : length(theFiles)
    baseFileName = theFiles(k).name;
    fullFileName = fullfile(myFolder, baseFileName);

    %% reading as an image array with im
    I = imread(fullFileName);
    % convert to grayscale
    Ig = rgb2gray(I);
    if plotFig
        figure;imagesc(I)
        figure;imagesc(Ig)
    end

    mm=mean(mean(Ig)); 
    if mm < 70
        I=imcomplement(I);
        Ig = imcomplement(Ig);
    end


    % BLOB DETCTION
%     h = fspecial('log', [15 15], 2);
%     imLOG = imfilter(Ig, h);
%     figure;imagesc(imLOG)

    %% find petridish by edges and binary operations
    % HACK - NOT HOW IT SHOULD BE DONE
    Ig = wiener2(Ig,[5 5]);
    imEdge = edge(Ig,'log');
    circle = bwareaopen(imEdge,50);
    circle = imclose(circle,strel('disk',30));
    circle = bwareaopen(circle,8000);
%     circle = imfill(circle,'holes');
    circle = bwconvhull(circle);
    circle = imerode(circle,strel('disk',150));
    if plotFig
        figure;imagesc(circle)
    end

    %% Get thresholds inside dish using otsu on each channel
    imR = double(I(:,:,1)) .* circle;
    imG = double(I(:,:,2)) .* circle;
    imB = double(I(:,:,3)) .* circle;

    thresR = graythresh(uint8(imR(circle))) *max(imR(circle));
    thresG = graythresh(uint8(imG(circle))) *max(imG(circle));
    thresB = graythresh(uint8(imB(circle))) *max(imB(circle));

    if plotFig
    figure;imagesc(imR)
    figure;imagesc(imG)
    figure;imagesc(imB)
    end

    %% classify inside dish
    % check if it should be smaller or larger than
    if sum(imR(circle) < thresR) >  sum(imR(circle) > thresR)
        labelR = imR > thresR;
    else
        labelR = imR < thresR;
    end

    if sum(imG(circle) < thresG) >  sum(imG(circle) > thresG)
        labelG = imG > thresG;
    else
        labelG = imG < thresG;
    end

    if sum(imB(circle) < thresB) >  sum(imB(circle) > thresB)
        labelB = imB > thresB;
    else
        labelB = imB < thresB;
    end

    if plotFig
        figure;imagesc(labelR)
        figure;imagesc(labelG)
        figure;imagesc(labelB)
    end

    labels = (labelR + labelG + labelB) .* circle;
    labels(labels < labelcutOff) = 0;
    labels = imfill(labels,'holes');
    labels = bwareaopen(labels,30);
    if plotFig
    figure;imagesc(labels)
    end

    %% clean up labels
    labelBig = bwareaopen(labels,removeLargerthan);
    labels = labels - labelBig;
    if plotFig
    figure;imagesc(labels)
    end


    BN = labels;


    %% old script


    stats = regionprops(BN,'Basic');
    obj2 = numel(stats);
    [B,L] = bwboundaries(BN,'holes');

    figure
%     imshow(label2rgb(L, @jet, [.5 .5 .5]))
    imshow(I)
    hold on
    title(baseFileName)

    for j = 1:length(B)
      boundary = B{j};

      plot(boundary(:,2), boundary(:,1),'w','LineWidth',2)
    end
    %region stats
    stats = regionprops(L,'Area','Centroid');
    %Threshold for printing in end
    threshold = 0.2;
    %Conversion factor pixel to cm
    conversionF=9/2125;

    % loop over the boundaries
    for j = 1:length(B)

      % obtain (X,Y) boundary coordinates corresponding to label 'j'
      boundary = B{j};

      % compute a simple estimate of the object's perimeter
      delta_sq = diff(boundary).^2;
      perimeter = sum(sqrt(sum(delta_sq,2)));
      perimeterD(j,k)=perimeter*conversionF;

      % obtain the area calculation corresponding to label 'k'
      area = stats(j).Area;
      areaD(j,k)=area*conversionF^2;

      % compute the roundness metric
      metric = 4*pi*area/perimeter^2;
      metricD(j,k)=metric;

      % display the results
      metric_string = sprintf('%d. %2.2f', j,metric);

      text(boundary(1,2)-50,boundary(1,1)+23,metric_string,'Color','k',...
           'FontSize',14,'FontWeight','bold');
    end

      drawnow; % Force display to update immediately.
end
%Calculating stats
areaM=mean(areaD);
pM=mean(perimeterD);
metricD(metricD==Inf)=0;
mM=mean(metricD);
Martin Agerbæk
  • 1
  • 2

1 Answers1

0

Hint:

A morphological top-hat filter fllowed by binarization (with a constant threshold ?) can be a good start. And filtering on the blob size will do a reasonable cleanup.

For the edges, try circular Hough.

enter image description here