-
Notifications
You must be signed in to change notification settings - Fork 0
/
countSpotsNearNuclei.m
44 lines (29 loc) · 1.05 KB
/
countSpotsNearNuclei.m
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
function spotTable = countSpotsNearNuclei(spotIm, dapiIm)
minNucleusSize = 1000;
donutRadius = 60;
% Deal with spots
h = -fspecial('log',20,2);
filt = imfilter(spotIm,h,'replicate');
irm = imregionalmax(filt);
% Get nuclei
T = adaptthresh(dapiIm,0.3,'ForegroundPolarity','bright');
dp = imbinarize(dapiIm,T);
CC = bwconncomp(dp);
rp = regionprops(CC);
area = [rp.Area];
idx = area > minNucleusSize; % Get rid of small stuff
CC2 = CC;
CC2.NumObjects = sum(idx);
CC2.PixelIdxList = CC2.PixelIdxList(idx);
lab = labelmatrix(CC2);
allSpots = [];
cellNumber = [];
for i = 1:CC2.NumObjects
tempbw = lab == i; % select the ith cell
tempbw = imdilate(tempbw,strel('disk',donutRadius)) - tempbw; % donut around ith cell
temprm = irm.*tempbw; % keep only spots in the donut
tempSpots = filt(temprm==1)'; % keep signal intensities for the spots in the donut
allSpots = [allSpots tempSpots];
cellNumber = [cellNumber i*ones(1,numel(tempSpots))];
end
spotTable = table(allSpots',cellNumber','VariableNames',["spotIntensities","cellNumber"]);