How can I detect bubbles on this image using Matlab Image Processing Toolbox?

Hello,

I've been trying to detect the bubbles on this image:

I have tried several different methods, both taken from the Matlab help and from this same forum. Below is an example of what I have used:

clearvars;
clc;
imtool close all;
close all;
%% Import image
I = imread('DSC_0438.jpg'); % picture with bubbles
% From previous experience, it is futile to manipulate the image when it
% has been binarized; most of the information should be discriminated while
% it is still a graysccale image
%% Convert image to grayscale
I = rgb2gray(I);
I2 = adapthisteq(I);
%% Remove background from image
background1 = imopen(I,strel('diamond',20));
I3 = I2 - background1;
%% Binarize image
I4 = imbinarize(I3,'adaptive');
I4 = bwmorph(I4,'majority');
se = strel('disk',6);
I4 = imclose(I4,se);
cc = bwconncomp(I4);
%% Filter out small pixel clusters by size and remove non-circular shapes
% Taken from the Tips section of the regionprops help
% The discriminating area ranges (100 and 1000) are equivalent to a
% bubble diameter of 44 and 139 microns, respectively. The maximum
% allowable eccentricity is 0.5 (0 is a circle and 1 is a line)
min_area = 100;
max_area = 10000;
max_eccentricity_global = 0.5;
stats = regionprops(cc,'Area','BoundingBox','Eccentricity','Centroid',...
    'MajorAxisLength','MinorAxisLength');
areas = [stats.Area]';
eccentricities = [stats.Eccentricity]';
idx = find(areas > min_area & areas < max_area & ...
    eccentricities < max_eccentricity_global);
I5 = ismember(labelmatrix(cc),idx);

%% Find the detected bubbles' centroids, diameters and radii
% Taken from the Examples section of the regionprops help
cc2 = bwconncomp(I5);
stats2 = regionprops('table',cc2,'Centroid','MajorAxisLength',...
    'MinorAxisLength');
centers = stats2.Centroid;
diameters = mean([stats2.MajorAxisLength stats2.MinorAxisLength],2);
radii = diameters/2;
diameters_um = diameters .* 3.9; %3.9 um is the pixel size of the D7200
mean_diameter_um = mean(diameters_um);
%% Draw circles around the detected bubbles and overlap with original image
figure
subplot(2,2,1)
imshow(I)
title('Original image')
subplot(2,2,2)
imshow(I)
hold on
viscircles(centers,radii);
hold off
title('Bubble detection')
subplot(2,2,3)
imshowpair(I4,I5)
title('Comparison of discriminated areas')

Basically, what I get is this:

I tried playing with the eccentricity and the thresholds I define for the areas and for the binarization, as well with the bwmorph function to eliminate noise, but I still get a lot of false positives.

I know that the bubbles on the binary photo look like (portions of) rings (see below)

I was wondering if I could define a neighborhood matrix for the strel function to keep those and eliminate the rest, but I am not sure how to implement it. Any help would be appreciated.

Regards,

Julio

[Rik]

It might be easier to try to remove the background. Even as a human I have trouble spotting the bubbles.

[Julio Cesar Garcia Navarro]

I have been trying that but, since a lot of bubbles have the same optical properties as the background, they ger removed whenever I remove the background.

[Rik]

How have you tried to remove the background? I would generate a prediction for the background and then look at the squared difference between predicted background and image.

[Julio Cesar Garcia Navarro]

What I did was use imopen and a manually-adjusted strel, but what I get is a lower saturation on the background, which still matches the saturation of the smallest bubbles.

The portion of the code I use is this one:

%%Remove background from image
background1 = imopen(I,strel('diamond',20));
I3 = I2 - background1;

I also tried to binarize both images (oroginal and background) before substracting them but I remove some bubbles as well.

[Rik]

The background is essentially a tiled image. I would try to find a suitable tile size, and then average all tiles to find the average background tile, replicate back to full image size and then subtract that from the original.

Maybe some playing around with a 2D Fourier transform would be easier.

[Julio Cesar Garcia Navarro]

Could you please give me a sample code for that? I get the idea but I do not know how to implement it.

[Rik]

Use imshow and zooming to find coordinates.

%I'll assume you have converted the corner coordinates to [number of rows, number of columns]

%IM=rand(4000,6000);
IM=imread('DSC_0438_modified.jpg');
IM=double(IM)/255;

ext=[425 553];
%from 615 to 5710 there are 12 repetitions, so ext_row=425
%from 290 to 3610 there are  6 repetitions, so ext_col=553

%replicate the first rows and cols of IM so it can be easily be split into tiles
row_idx=[1:size(IM,1) 1:size(IM,1)];
col_idx=[1:size(IM,2) 1:size(IM,2)];
row_idx=row_idx(1:ceil(size(IM,1)/ext(1))*ext(1));
col_idx=col_idx(1:ceil(size(IM,2)/ext(2))*ext(2));
IM2=IM(row_idx,col_idx);

n_tiles=prod(size(IM2)./ext);
tile=zeros(ext(1),ext(2),n_tiles);
tilecounter=0;
for cur_row=1:(size(IM2,1)/ext(1))
    for cur_col=1:(size(IM2,2)/ext(2))
        tilecounter=tilecounter+1;
        r=ext(1).*[cur_row-1 cur_row]+[1 0];
        c=ext(2).*[cur_col-1 cur_col]+[1 0];
        r=r(1):r(2);
        c=c(1):c(2);
        tile(:,:,tilecounter)=IM2(r,c);
    end
end
figure(888),clf(888)%random figure
subplot(2,1,1)
imshow(mean(tile,3),[])
subplot(2,1,2)
montage(permute(tile,[1 2 4 3]))

From this result you can easily see that the image is not as straight as I hoped, which means you might need to keystone-correct the image first. Or you could average the 4-neighborhood of each tile to get a local average, because the shift is not large from tile to tile. (the 4-neighborhood is the tile itself, one left, one right, one up and one down)

[Julio Cesar Garcia Navarro]

Great, thanks mate. I can make progress now from that.

Answer 1

[Julio Cesar Garcia Navarro]

UPDATE: I used rgb2hsv and got rid of the h component and now I have something like this:

I can say now that I am looking for the blue blobs that are surrounded by green. How do I program that in MATLAB?

[Image Analyst]

How can you do that? The image, even if RGB, is essentially gray scale. So the H and S components would be zero, and what you're showing must be the V component, which is essentially what you'd get by converting your image to gray scale with rgb2gray(). It looks like you just have a gray scale image (the v component) with some colormap applied so essentially you're back to dealing with a monochrome/grayscale thresholding situation.

[Julio Cesar Garcia Navarro]

Im sorry; the original image is this one:

What I did was this:

I2 = rgb2hsv(I);
h = I2(:,:,1);
s = I2(:,:,2);
v = I2(:,:,3);
a = zeros(size(I, 1), size(I, 2));
sv = cat(3, a, s, v);

Thus, s is green and v is blue.

[Image Analyst]

It doesn't look like the bubbles are a different color so doing color classification doesn't make sense or help.

You could find particles that are surrounded by black, but it looks like there are some that are touching the grid. So I'd recommend creating a template by imaging a grid known to have no bubbles are all. Then align the image and subtract the actual from the template to show only what's different in the two images, which is presumably only the bubbles.

[Julio Cesar Garcia Navarro]

Momentarily taking more pictures is not an option. How do I find particles surrounded by black? Those are the ones I am actually interested in.

[Julio Cesar Garcia Navarro]

I happened to find a picture with almost no bubbles, so I used imregister on the red channel of both images. What I end up with is the following:

I have tried to use regionprops to discriminate bubbles from lines using circularity and diameter (averaging the Major- and MinorAxisLength), but I end up with many false positives.

Is there a way to remove the gray lines without touching bubbles?

[Image Analyst]

I'm not really sure what are bubbles. Are the big white C-shaped things bubbles? What about the clusters of 3 or 4 small points inside a grid cell? What threshold did you use? And please indicate what diameter and circularity thresholds you used, and what blobs were detected that should not have been.

[Julio Cesar Garcia Navarro]

I took a portion of the image from my previous post:

The red circles denote the bubbles I am trying to detect, while the blue circles denote the clusters of 3 or 4 small points you mention, which are actually reflections of light on larger bubbles. What I am after are the red-marked bubbles.

So far, the most effective algorithm I have come up with is calling regionprops after I do opening-closing by reconstruction:

I = imread('DSC_0438_nobgrnd.jpg'); %The image form my previous post
se = strel('disk',10);
Io = imopen(I, se);
Ie = imerode(I, se);
Iobr = imreconstruct(Ie, I);
Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);

This portion of the algorithm produes something like this:

This image gets rid (form my perspective) of most of the background that was not removed using imregister; the downside is a loss of the smallest bubbles (I can live with that). From here I use the following:

I2 = imbinarize(Iobrcbr);

I2 = bwareaopen(I2,30);

I2 = imclose(I2,se);

I2 = imfill(I2,'holes');

cc = bwconncomp(I2);

min_circ = 0.6; %Minimum allowable circularity
max_radius = 70; %Maximum allowable radius

stats = regionprops('table',cc,'Area','Perimeter','MajorAxisLength',...
    'MinorAxisLength');

areas = stats.Area;
perims = stats.Perimeter;

circs = 4.*pi.*areas./perims.^2;
diameters = mean([stats.MajorAxisLength stats.MinorAxisLength],2);
radii = diameters/2;

idx = find(circs > min_circ & radii < max_radius);

I3 = ismember(labelmatrix(cc),idx);

cc2 = bwconncomp(I3);
stats2 = regionprops('table',cc2,'Centroid','MajorAxisLength',...
    'MinorAxisLength');
centers = stats2.Centroid;
diameters2 = mean([stats2.MajorAxisLength stats2.MinorAxisLength],2);
radii2 = diameters2/2;
diameters_um = diameters2 .* 3.9; %3.9 um is the pixel size of the D7200
mean_diameter_um = mean(diameters_um);
figure,imshow(I)
viscircles(centers,radii2)

I call regionprops twice to use the first one as discrimination and the second one to find and measure the bubbles:

And here is the complete picture:

I got rid of most of the false positives but the detection turnout is quite low: on the third image of this post there are 4 red circles; the one on the bottom-left is not really a small bubble but part of a larger one.

I hand-adjusted the values of min_circ and max_radius. Could you think of a different algorithm that is more successful with the smaller bubbles?

Thank you