Automating colour segmentation script

Using the code below, my aim is to get the coordinates of ants that are painted with green blobs.

To do this, I manually draw a shape around the desired paint spot to use as a colour reference to detect similar coloured pixels. I want to be able to do this for one image, and then use the same colour reference for all images in a folder without having to redraw the shape every time.

The code below " ************ " can be run again after reading in a different image, using the first image as a reference. The issue is, I have no idea how to automate this process - I need to automatically read in a new image, and run all the code below *******, saving two files - a text file with the coordinates, and the mask image.

I've tried to create a function to do this, but haven't had any luck yet. Any pointers or advice would be massively appreciated.

% read image, convert to LAB
ant = imread("IMG_4335.GREEN.JPG");
ant = imgaussfilt(ant,2); % introduce gaussian filter, to avoid jagged edges later
antLAB = rgb2lab(ant); % convert RGB to LAB colour space
antL = antLAB(:,:,1); % luminosity
antA = antLAB(:,:,2); % red-green axis
antB = antLAB(:,:,3); % blue-yellow axis

%% Draw shape filled with a paint spot - want to run this for ONE image in a folder,
%  to use the values as a reference for ALL images
imshow(ant)
roi = drawpolygon % draw shape to capture the colour spot you want, then press enter

% use the shape above (roi), to get pixel values of drawn area
BWant = createMask(roi); % create a mask using the drawn shape
aROI = antA(BWant); %get the pixel values for A channel
bROI = antB(BWant); % get the pixel values for B channel
meanROI = [mean(aROI) mean(bROI)]; % take mean values for each channel
aMean = meanROI(1); % index above, channel a
bMean = meanROI(2); %index above, channel b

% ************

% read a second image (remove %s), convert to LAB
% ant = imread("SECOND.IMAGE.JPG");
% ant = imgaussfilt(ant,2); % introduce gaussian filter, to avoid jagged edges later
% antLAB = rgb2lab(ant); %convert RGB to LAB colour space
% antL = antLAB(:,:,1); % luminosity
% antA = antLAB(:,:,2); % red-green axis
% antB = antLAB(:,:,3); % blue-yellow axis

%% create distance matrix - diff between avg colour and every pixel
distLab = sqrt((antA - aMean).^2 + (antB - bMean).^2);

% treshold the distance matrix - which pixels are similar to sample colour
mask = distLab < 13; % smaller - stricter threshold

% dilate image, if two dots close, should merge them. Be careful with busy images
se = strel("disk",6); % create disk shaped area to use to dilate image next
dilatedant = imdilate(mask,se); % dilate iamge - to join blotches that should be one blob
maskfilledholes = imfill(dilatedant,"holes"); % fill in holes
maskcleaned = bwareaopen(maskfilledholes,10); % remove pixels smaller than 10

% find coordinates of centroids, write into text file
antmeasurements = regionprops(maskcleaned, "all"); 
centroids = cat(1,antmeasurements.Centroid); % turn centroids into two columns (x and y coordinates)
dlmwrite("antcoordinates.txt",centroids);

% save output image to filename
imwrite(maskcleaned,"outputantimage.jpg");

[KALYAN ACHARJYA]

It is quite uneasy to read long text, please be specific

[BC]

Sorry, I've tried to make it a bit shorter!

[KALYAN ACHARJYA]

Is your task to count such ants (with green color) or do you have to segments all such ants?

[BC]

My aim is to get a coordinate for each ant with a green paint mark. I don't need to count the number of green ants, and I can ignore any ants that aren't painted.

Answer 1

[BC]

Ok so I think I've written a function that helps somewhat - slightly different from above, same principle. I now need to work out how to read the next image in a directory, and apply the local function to this image, and output the files as the same as the image name in the folder. Updating in case this code is useful for anyone else.

%% read REFERENCE image to get colour
ant = imread("IMG_4335.GREEN.JPG");
antLAB = rgb2lab(ant); % convert RGB to LAB colour space
antL = antLAB(:,:,1); % luminosity
antA = antLAB(:,:,2); % red-green axis
antB = antLAB(:,:,3); % blue-yellow axis

% Draw shape filled with a paint spot - run this at start of new folder
% to get baseline reference colour
imshow(ant)
roi = drawpolygon % draw shape to capture the colour spot you want, then press enter

% use the shape above (roi), to get pixel values of area
BWant = createMask(roi); % create a mask using the drawn shape
aROI = antA(BWant); %get the pixel values for A channel
bROI = antB(BWant); % get the pixel values for B channel
meanROI = [mean(aROI) mean(bROI)]; % take mean values for each channel
aMean = meanROI(1); % index above, channel a
bMean = meanROI(2); %index above, channel b

%% use local function to get output for reference image
[finalantimage centroids] = createantmask(ant, antA, aMean, bMean);
% function in the format: outputs = functionname(inputs)

%% read a different image, use createantmask function again without
%  having to redraw shape
%  using pixel values from first image as a reference
ant = imread("Y2.JPG");
[finalantimage centroids] = createantmask(ant, antA, aMean, bMean);

%% plot centroids on image, check looks ok
imshow(finalantimage)
hold on
plot(centroids(:,1),centroids(:,2),"b*");
hold off

%% display images side by side, check not missing any blobs etc.
fig = figure();
ax(1) = axes('Units','normalized','Position', [ .1 .1 .4 .8]);
ax(2) = axes('Units','normalized','Position', [ .5 .1 .4 .8]);
imshow(finalantimage, 'Parent', ax(1)) %display mask
imshow(ant, 'Parent', ax(2)) %display original image
linkaxes(ax)


% create local function below, to apply to image named "ant"
% have to run the first section first to get the reference colour pixels
% and required inputs

%%
function [finalantimage centroids] = createantmask(ant, antA, aMean, bMean)
ant = imgaussfilt(ant,2); % introduce gaussian filter, to avoid jagged edges later
antLAB = rgb2lab(ant); % convert RGB to LAB colour space
antL = antLAB(:,:,1); % luminosity
antA = antLAB(:,:,2); % red-green axis
antB = antLAB(:,:,3); % blue-yellow axis
distLab = sqrt((antA - aMean).^2 + (antB - bMean).^2); % distance matrix. 
% difference between average colour and every pixel

% treshold the distance matrix - which pixels are similar to sample colour.
% smaller = stricter treshold, but risk cutting off some blobs
mask = distLab < 10; 
% value not much of an issue when just detecting largest blob only.

% dilate image, if two dots close, should merge them.
se = strel("disk",7); % create disk shaped area to use to dilate image next
dilatedant = imdilate(mask,se); % dilate image - to join blotches that should be one blob
maskfilledholes = imfill(dilatedant,"holes"); % fill in holes
maskcleaned = bwareaopen(maskfilledholes,10); % remove leftover pixels smaller than 10

% filter only largest blob - only for the queen coordinate
% DON'T USE FOR NON-QUEEN IMAGES
% finalantimage = bwpropfilt(maskcleaned,"area",1,"largest");

% find coordinates of centroids, write into text file
antmeasurements = regionprops(finalantimage, "all"); 
centroids = cat(1,antmeasurements.Centroid); % turn centroids into two columns
dlmwrite("antcoordinateS.txt",centroids); % (x and y coordinates)

% save jpg output of the final mask
imwrite(maskcleaned,"outputantimage.jpg");
end