Hi Tushar
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
