Cell-Counter/CellCount.m at master · mdaddysman/Cell-Counter · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
%CellCount.m
clearvars;

av1=1; %width of the narrow gaussian
av2=20; %width of the wide gaussian
thresh=1.8; %difference between the gaussians for detection
cellsize = 150; %the minimum size of the cell
cutoff=150;
name='example.png';

%box
x = 70;
y = 25;
xsize = 290;
ysize = 295;

xbkg = 1;
ybkg = 1;
textcolor = 'r';

A = 256*uint16(imread(name));
A_box = A;
A_box(y:y+ysize,x:x+1) = 65536;
A_box(y:y+1,x:x+xsize) = 65536;
A_box(y:y+ysize,x+xsize:x+xsize+1) = 65536;
A_box(y+ysize:y+ysize+1,x:x+xsize+1) = 65536;
A_box(ybkg:ybkg+20,xbkg:xbkg+1) = 65536;
A_box(ybkg:ybkg+1,xbkg:xbkg+20) = 65536;
A_box(ybkg:ybkg+20,xbkg+20:xbkg+21) = 65536;
A_box(ybkg+20:ybkg+21,xbkg:xbkg+21) = 65536;
figure(1); imshow(A_box)

H1=fspecial('gaussian',40,av1);
H2=fspecial('gaussian',40,av2);
a1=imfilter(A,H1,'replicate');
a2=imfilter(A,H2,'replicate');
a=a1-a2;

a_bw=medfilt2(a>(thresh*256));
a_bw2=bwmorph(a_bw,'open',inf);
a_bw3=imfill(a_bw2,'holes');

Amasked1=uint16(a_bw3).*A;
Am_box1 = Amasked1;
Am_box1(y:y+ysize,x:x+1) = 65536;
Am_box1(y:y+1,x:x+xsize) = 65536;
Am_box1(y:y+ysize,x+xsize:x+xsize+1) = 65536;
Am_box1(y+ysize:y+ysize+1,x:x+xsize+1) = 65536;
Am_box1(ybkg:ybkg+20,xbkg:xbkg+1) = 65536;
Am_box1(ybkg:ybkg+1,xbkg:xbkg+20) = 65536;
Am_box1(ybkg:ybkg+20,xbkg+20:xbkg+21) = 65536;
Am_box1(ybkg+20:ybkg+21,xbkg:xbkg+21) = 65536;

figure(2), imshow(Am_box1)

bkg=mean2(A(ybkg:20+ybkg,xbkg:20+xbkg))/256;
meanpix1=sum(sum(Amasked1(y:y+ysize,x:x+xsize)))/sum(sum(a_bw3(y:y+ysize,x:x+xsize)))/256;
meanpix1_bkgcorr=meanpix1-bkg

number_bins=20;
L=bwlabel(a_bw3);
objects=max(max(L));
a_bw4=false(size(a_bw3));
for m=1:objects
    av(m)=sum(sum(uint16(L==m).*A))./sum(sum(uint16(L==m)))/256;
    if av(m)<cutoff
        a_bw4=a_bw4+(L==m);
    end
end

bins=(floor(0.9*min(av))):(ceil(1.1*max(av)));
figure(3)
hist(av,bins), axis tight
avg=mean(av);
stdev=std(av);
av_2stdev=avg+2*stdev;

Amasked2=uint16(a_bw4).*A;
Am_box2 = Amasked2;
Am_box2(y:y+ysize,x:x+1) = 65536;
Am_box2(y:y+1,x:x+xsize) = 65536;
Am_box2(y:y+ysize,x+xsize:x+xsize+1) = 65536;
Am_box2(y+ysize:y+ysize+1,x:x+xsize+1) = 65536;
Am_box2(ybkg:ybkg+20,xbkg:xbkg+1) = 65536;
Am_box2(ybkg:ybkg+1,xbkg:xbkg+20) = 65536;
Am_box2(ybkg:ybkg+20,xbkg+20:xbkg+21) = 65536;
Am_box2(ybkg+20:ybkg+21,xbkg:xbkg+21) = 65536;
figure(4), imshow(Am_box2)
meanpix2=sum(sum(Amasked2(y:y+ysize,x:x+xsize)))/sum(sum(a_bw4(y:y+ysize,x:x+xsize)))/256;
meanpix2_bkgcorr=meanpix2-bkg

a_interest = a_bw4(y:y+ysize,x:x+xsize);

[L cells] = bwlabel(a_interest);
cells
STATS = regionprops(L,'basic');
smallcells = 0;
scidenity = 0;
for m=1:cells
    centroid(m,1) = STATS(m).Centroid(1,1);
    centroid(m,2) = STATS(m).Centroid(1,2);
    area(m) = STATS(m).Area;
    if(area(m) < cellsize)
        smallcells = smallcells + 1;
        scidenity(smallcells) = m;
    end
end
cells_threshold = cells - smallcells
figure(5)
imshow(uint16(a_interest).*A(y:y+ysize,x:x+xsize))
for m=1:cells
    text(centroid(m,1),centroid(m,2),num2str(m),'Color',textcolor)
end
text(0,-10,['Ex. Cells: ' num2str(scidenity)],'Color','k')