How do I characterize these cells by area and diameter?

I have this picture:

The particles of interest are labeled in red. I would like to quantify these labeled images by area and diameter.

For reference, the above image came from this image:

I think there might be a way to take this image and contrast the white to the black and come up with an area and diameter, but I'm not really.

As always,

Thanks for your time,

Respectfully,

Tushar

Answer 1

[John BG]

Hi Tushar

this is John BG <mailto:[e-mail] [e-mail]>

Comment: the start image of this answer is part of the batch including the image posted in the question. This solution works on both non-uniformly illuminated and low contrast samples.

This sample has been taken with spot illumination causing low contrast. Tubular straight or circle fluorescent lighting come really handy to achieve uniform illumination. Attached script.

Anyway let's crack on:

1.

Acquiring image

clear all;clc;close all
A=imread('011.jpg');
figure(1);imshow(A);
% 011

.

Checking possible ways to improve contrast

.

minA=min(min(min(A)))
maxA=max(max(max(A)))
rangeA=abs(maxA-minA)
 minA =
   uint8
    36
 maxA =
   uint8
    134
 rangeA =
   uint8
    98

.

problem 1: poor contrast problem 2: no colour markers to use

% A=10*(A-minA); doesn't work because of problem 3 point illumination

A2=A(:,:,1)+A(:,:,2)+A(:,:,3);    % a bit more contrast
[sz1 sz2 sz3]=size(A2);

minA2=min(min(A2))
maxA2=max(max(A2))
rangeA2=abs(maxA2-minA2)
% minA2 =
%   uint8
%    108
% maxA2 =
%   uint8
%    255
% rangeA2 =
%   uint8
%    147
% contrast improved a bit, but the main problem

A2=A2-minA2;
minA2=min(min(A2))
maxA2=max(max(A2))
rangeA2=abs(maxA2-minA2)
% minA2 =
%   uint8
%    0
% maxA2 =
%   uint8
%    147
% rangeA2 =
%   uint8
%    147
% zeroed

A2=255/rangeA2*A2;
minA2=min(min(A2))
maxA2=max(max(A2))
rangeA2=abs(maxA2-minA2)
figure(2);imshow(uint8(A2))
% 012

.

now contrast reaches full range, but problem 3: poor illumination, non uniform, centre excessive lighting while corner not enough light

. Trying Laplacian

del2A2=del2(double(A2));
figure(3);imshow(del2A2)
% 013

.

% doesn't work

may be var

% varA2sq=var(del2A2,40,3);
% figure(4);h4=surf(varA2sq);h4.EdgeColor='none'
% 
% minvarA2sq=min(min(varA2sq))
% maxvarA2sq=max(max(varA2sq))
% rangevarA2sq=abs(maxvarA2sq-minvarA2sq)
% all nulls

% doesn't work either
% A1=A(:,:,1);
% varA1=var(double(A1),0,3);
% figure(6);h6=surf(varA1);h6.EdgeColor='none

.

applying corrective surface

x_range=linspace(-2,2,863);y_range=linspace(0,3,644)+2;[X Y]=meshgrid(x_range,y_range);
Z=cos(.5*X)+.8*sin(.4*Y);
f4=figure(4);h4=surf(Z);h4.EdgeColor='none'
ax4=gca(f4);ax4.DataAspectRatio=[1 1 1]
% ax4.Position=[0.130000000000000   0.110000000000000   0.775000000000000   0.815000000000000]

taking some measurements on the start picture

%  scaling
Z=Z*128/1.795;

f7=figure(7);h5=surf(Z);h7.EdgeColor='none'
ax7=gca(f7);ax7.DataAspectRatio=[1 1 1]
ax7.CameraPosition=[450 350 5700]

% Z=uint8(Z);
f8=figure(8);h6=surf(Z);h8.EdgeColor='none'
ax8=gca(f8);ax8.DataAspectRatio=[1 1 1]
ax8.CameraPosition=[450 350 5700]

imadjust same problem

A1=A(:,:,1);A3=imadjust(A1);
figure(1);imshow(A1)
% 011-2

.

2.- Counting only across illuminated belts, one at a time, no Gauss filter

.

clear all;clc;close all

th1=[140:-5:45];

A=imread('011.jpg');
hf1=figure(1);imshow(A);hold all

A1=A(:,:,1);
[sz1 sz2 sz3]=size(A1);

A5=uint8(zeros(sz1,sz2,numel(th1)));

for k=1:1:numel(th1)
    A10=A1;
    A10(A10>th1(k))=255;
    A10(A10<=th1(k))=0;
     A5(:,:,k)=A10;
    figure(2);imshow(A10);
    drawnow
end

Radii_cell={};Metric_cell={};Centers_cell={};
for k=1:1:numel(th1)

      A21=A5(:,:,k);
  %     A21=-(A21-255);
      hf=figure(k+10);imshow(A21);hold on 

%     A30=A30(1,:,:);

%     A30(A30>125)=130;A30(A30<=125)=255;A30(A30==130)=0;  % invert
%     figure(k+10);imshow(A30)

      [centers, radii, metric] = imfindcircles(A21,[15 45],'Sensitivity',1,'EdgeThreshold',.85,'ObjectPolarity','dark','Method','PhaseCode')

      if ~isempty(centers)
          [sz1cent sz2cent]=size(centers)
          centersStrong5 = centers(1:min(sz1cent,20),:);
          radiiStrong5 = radii(1:min(sz1cent,20));
          metricStrong5 = metric(1:min(sz1cent,20));
          viscircles(gca(hf1),centersStrong5, radiiStrong5,'EdgeColor','r');hold all

          viscircles(gca(hf),centersStrong5, radiiStrong5,'EdgeColor',[.8 .3 .3],'LineWidth',1);hold all

          % print circles on start image
          viscircles(gca(hf),centersStrong5, radiiStrong5,'EdgeColor',[.8 .3 .3],'LineWidth',1);hold all 
          Radii_cell{k}=radiiStrong5;
          Metric_cell{k}=metricStrong5;
          Centers_cell{k}=centersStrong5;
      end

end

% 020: radii range [4 45]

% 020-2: radii range [15 45]

.

yes but not quite Karlsberg, forget Carling

3.- With Gauss filter

clear all;clc;close all

    th1=[140:-5:45];
    PSF = fspecial('gaussian',5,5);

    A=imread('011.jpg');
    hf1=figure(1);imshow(A);hold all

    A1=A(:,:,1);
    [sz1 sz2 sz3]=size(A1);

    A5=uint8(zeros(sz1,sz2,numel(th1)));

    for k=1:1:numel(th1)
        A10=A1;
        A10(A10>th1(k))=255;
        A10(A10<=th1(k))=0;
         A5(:,:,k)=A10;
        figure(2);imshow(A10);
        drawnow
    end

    % 1 loop explained

    Boundaries_cell={};
    for k=1:1:numel(th1)

        A21=A5(:,:,k);
    %     A21=-(A21-255);
        hf=figure(k+10);imshow(A21);hold on 
     % 040

         A3 = deconvlucy(uint8(A21),PSF,5);
         figure(k+10);imshow(A3);hold on
         % 041

         A3(A3>125)=130;A3(A3<=125)=255;A3(A3==130)=0;  % invert
         figure(k+10);imshow(A3);hold on
         % 042

         [B,L] = bwboundaries(A3,'noholes');

    %      if ~isempty(B)  % remove boundary exceedingly long
    %         LG=[];
    %         for s=1:1:length(B) 
    %             L1=B{s};
    %             LG=[LG length(L1)]; 
    %         end % measuring lengths of each boundary
    %         if max(LG)>50 % only remove long shore boundaries
    %             n_remove=find(LG==max(LG)); % find longest boundary
    %             n0=[1:1:length(B)];
    % 
    %             n0(n_remove)=[];
    %             B=B(n0);
    %         end
    %     end

        for d = 1:length(B)
            boundary = B{d};
            plot(gca(hf),boundary(:,2), boundary(:,1), 'b', 'LineWidth', 1);hold on;      
        end

         for d = 1:length(B)  % repeat cell skin marking on start image
            boundary = B{d};
            plot(gca(hf1),boundary(:,2), boundary(:,1), 'b', 'LineWidth', 1);hold on;      
         end

        % log boundaries for each area measured in each for loop round
        Boundaries_cell{k}=B;

    end

    % 044 iso marks are counted as boundaries

.

4.- added the following, it removes watermark long lines, but not all of them

.

   if ~isempty(B)  % remove boundary exceedingly long
        LG=[];
        for s=1:1:length(B) 
            L1=B{s};
            LG=[LG length(L1)]; 
        end % measuring lengths of each boundary
        if max(LG)>50 % only remove long shore boundaries
            n_remove=find(LG==max(LG)); % find longest boundary
            n0=[1:1:length(B)];

            n0(n_remove)=[];
            B=B(n0);
        end
   end

.

5.- result approximation

sum(count_Boundaries)

   ans =
           1760

.

Comment: this count is rough since there may be repeated cells due to the low contrast in the illumination of the sample

.

6.- correction attempt: remove double counts

.

% out of simple visual inspection, way too many, there are repeated boundaries
match_pairs=[0 0];
for k=1:1:numel(Boundaries_cell)-1
    L3=Boundaries_cell{k};
    L4=Boundaries_cell{k+1};
    % build combinations with repetition
    %     N1=combinator(numel(L3),numel(L4),'c','r'); when k=10 amounts exceed
    % combinator capacity
    N1=[0 0];  
    for t1=1:1:numel(L3)
        for t2=1:1:numel(L4)
            N1=[N1;t1 t2];
        end
    end
    N1(1,:)=[];

      for s=1:1:numel(N1)
          if isequal(L3{N1(k,1)},L4{N1(k,2)})
              match_pairs=[match_pairs ;N1(k,:)];
          end
      end

end

.

7.- recap

.

hf30=figure(30);imshow(A);hold all
for k1=1:1:numel(Boundaries_cell)
    B0=Boundaries_cell{k1}
    for k2 = 1:length(B0)  % repeat cell skin marking on start image
        boundary = B0{k2};
        plot(gca(hf30),boundary(:,2), boundary(:,1), 'b', 'LineWidth', 1);hold all;      
    end
end

hf40=figure(40);imshow(A1);

.

save these to to retrieve them after clean slate .

    saveas(hf30,'watermarks01.jpg')
    save data1.mat 'Boundaries_cell' 'count_Boundaries'

.

8.- start again with the already refined image

.

clear all;close all;clc
load data1.mat
A=imread('watermarks01.jpg');
hf1=figure(1);imshow(A);

<</matlabcentral/answers/uploaded_files/90608/055.jpg>>

A1=A(:,:,1);   % on blue only
th10=45;      % threshold again, have to try different thresholds manually until you see all clear enough to count
A1(A1>th10)=255;
A1(A1<=th10)=0;
figure(2);imshow(A1);
% z01

.

9.- crop a bit, only a littel bit

.

% look at z02: 
% z02

% on far right only about 29 cells, eye count, and on far
% left site about 6 only, so, crop and add 25 to end count
[A2,crop_rectangle]=imcrop(A1)
% figure(3);imshow(A2)

% 10.- binarise
A3=imbinarize(A2);
A3=~A3;
figure(4);imshow(A3);
%z03

.

11.- fill up hollow cells

.

A4=imfill(A3,'holes');
figure(5);imshow(A4);
%z04

% 12.- remove small specks that are clearly not cells
area_threshold=10;
A5 = bwareaopen(A4,area_threshold);
hf6=figure(6);imshow(A5);

% attempting to fit circles doesn't work, not really many circle shapes around
% [centers, radii, metric] = imfindcircles(A5,[5 45],'Sensitivity',.8,'EdgeThreshold',1,'ObjectPolarity','dark','Method','PhaseCode')
% [sz1cent sz2cent]=size(centers)
% centersStrong5 = centers(1:min(sz1cent,20),:);
% radiiStrong5 = radii(1:min(sz1cent,20));
% metricStrong5 = metric(1:min(sz1cent,20));
% viscircles(gca(hf6),centersStrong5, radiiStrong5,'EdgeColor','r');hold all
% 
% viscircles(gca(hf1),centersStrong5, radiiStrong5,'EdgeColor','r');hold all
% viscircles(gca(hf),centersStrong5, radiiStrong5,'EdgeColor',[.8 .3 .3],'LineWidth',1);hold all

% 13.- the closing step, count objects
[B,L] = bwboundaries(A5,'noholes');
map=zeros(length(B),3);cmap=colormap(map);
hf7=figure(7);imshow(label2rgb(L,cmap, [.5 .5 .5]))
% z06

hold on
for k = 1:length(B)
   boundary = B{k};
   plot(gca(hf7),boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1)
end
% z07

.

.

14.- Result

When counting directly B

   numel(B)
   ans =
      556

Result=numel(B)+29+6-9
  =
    582

+29 recover cropped far right +6 recover cropped far left -9 do not take into account the few long watermarks still left behind

Comment: splitting the image into for instance 6 different squares may help process over more uniform illumination areas and therefore avoid highlighted problems.

Now the areas and diameters:

From image z06

A1=A(:,:,1);imshow(A1);

A2=imbinarize(A1);
imshow(A2)

[C, R, metric] = imfindcircles(A2,[2 100],'ObjectPolarity','dark','Method','TwoStage');
viscircles(C,R,'EdgeColor','b');

.

.

the centres are listed in C, radii in R, and the sought diameters in D

D=2*R;

note that the majority of cells have a circumscribed circle diameter between 4 and 10

histogram(D)

I have also tried command regionprops, the areas as requested, note that regionprops applies a coarse quantizing of the measured areas. I do not agree with the following, I would write a function that measures, for each cell, departing from each cell centre pixel, the actual amount of pixels, but that would take more time, and using each circumscribing circle as count limit to avoid propagating count to adjacent cells and wrongly measure cells as larger than actually are. Anyway, to close the comment on regionprops:

Table_stats_regionprops_search_Area = regionprops('table',A2,'Area')
 9×1 table
      Area
     ______
     497365
       3380
       1742
          5
        366
          3
          2
          1
        174

Perhaps also of interest:

Table_stats_regionprops_boundingbox = regionprops('table',A2,'BoundingBox')
format short
Table_stats_regionprops_boundingbox{:,:}
=
    0.5000    0.5000  798.0000  664.0000
    0.5000  530.5000   95.0000  134.0000
    0.5000  574.5000   58.0000   90.0000
    0.5000  627.5000    2.0000    4.0000
    0.5000  632.5000   24.0000   32.0000
  400.5000  409.5000    2.0000    2.0000
  454.5000  663.5000    2.0000    1.0000
  696.5000  663.5000    1.0000    1.0000
  793.5000  180.5000    5.0000   53.0000

.

And the broader sides of the rectangles found may be used as ellipses centres :

.

    Table_stats_regionprops_search_centroid = regionprops('table',A2,'Centroid', 'MajorAxisLength','MinorAxisLength')
    Table_stats_regionprops_search_centroid{:,1}
     =
      400.6168  332.0258
       40.4891  609.2843
       23.1045  631.4202
        1.2000  629.8000
        8.1366  650.9945
      401.6667  410.3333
      455.5000  664.0000
      697.0000  664.0000
    796.7644  208.3391

.

approximated ellipses major axis:

.

   Table_stats_regionprops_search_centroid{:,2}
  =
    913.8992
    166.3650
     96.6134
      4.8819
     33.6375
      2.5820
      2.3094
      1.1547
     60.7072

approximated major axis:

.

Table_stats_regionprops_search_centroid{:,3}
 =
  756.5742
   30.4791
   27.0545
    1.7759
   17.7974
    1.7638
    1.1547
    1.1547
    4.3422

.

.

If you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

Answer 2

[Jyotish Robin]

Hi Tushar!

From your problem statement, I understand that you are looking for a way in which you can identify the area and diameter of the circular approximations of the white blobs (as shown in the second figure).

I would suggest you to have a look at the following ML Answers post as it addresses a similar query: https://www.mathworks.com/matlabcentral/answers/94507-how-can-i-visually-isolate-and-measure-a-round-object-within-my-matlab-image

Hope this helps.

Regards,

Jyotish