-
Notifications
You must be signed in to change notification settings - Fork 1
/
nodeGrps_vesSegment.m
189 lines (173 loc) · 5.97 KB
/
nodeGrps_vesSegment.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
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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
function im = nodeGrps_vesSegment( nodePos,nodeEdges )
%%%%%%%%%%
% TO DO:
% correct eLen (segLength) to be in um (non-cubical voxels problem) - easy
% correct segDiam to be in um (non-cubical voxels problem) - more
% complicated - we already have calculated diameters in imView in voxels
%%%%%%%%%%
% if ~isfield(im,'nodeGroupFlag')
% im.nodeGroupFlag = 1;
% end
% if im.nodeGroupFlag==0
% answer = questdlg('Flag is 1. Do you want to run this again?','Run Again?','Yes','No','Yes');
% if ~strcmp(answer,'Yes'),
% return;
% end;
% end
% im.nodeGroupFlag = 0;
hwait = waitbar(0,'Getting Segment and Group Info');
% [nB,im] = nBupdate( im );
nB = zeros(size(nodePos));
nN = size(nodePos,1);
for ii=1:nN
nB(ii) = length(find(nodeEdges(:,1)==ii | nodeEdges(:,2)==ii));
end
nN = size(nodePos,1);
% hxy = im.Hvox(1);
% hz = im.Hvox(3);
hxy = 2;
hz = 2;
% create nodePos_um variable, to have positions in um
nodePos_um = nodePos; %this is a simple fix but we also need to
%update im.nodePos_um whenever we update im.nodePos
nodePos_um(:,1:2) = nodePos(:,1:2).*hxy;
nodePos_um(:,3) = nodePos(:,3).*hz;
lst3p = find(nB>=3 | nB==1);
edgeSegN = zeros(size(nodeEdges,1),1);
nodeSegN = zeros(size(nodePos,1),1);
nodeGrp = zeros(nN,1);
nSeg = 1;
nGrp = 0;
segNedges = [];
segEndNodes = [];
for ii = 1:length(lst3p)
ii
if isequal(rem(ii,100),0)
waitbar(ii/length(lst3p),hwait);
end
iN = lst3p(ii);
if nodeGrp(iN)==0
nGrp = nGrp + 1;
nodeGrp(iN) = nGrp;
end
nGrpN = nodeGrp(iN); %nGrp;
lstE = find(nodeEdges(:,1)==iN | nodeEdges(:,2)==iN);
iNstart = iN;
for jj = 1:length(lstE)
iN = lst3p(ii);
eIdx =lstE(jj);
% correct this to be in um (non-cubical voxels)
eLen = sum(diff( nodePos(squeeze(nodeEdges(eIdx,:)),:), 1, 1).^2).^0.5;
eLen_um = sum(diff( nodePos_um(squeeze(nodeEdges(eIdx,:)),:), 1, 1).^2).^0.5;
iN = setdiff(unique(nodeEdges(eIdx,:)), iN);
if nodeGrp(iN)==0
nodeGrp(iN) = nGrpN;
elseif nodeGrp(iN) < nGrpN
lst = find(nodeGrp==nGrpN);
nodeGrp(lst) = nodeGrp(iN);
nGrpN = nodeGrp(iN);
elseif nodeGrp(iN) > nGrpN
% I think this is fine
lst = find(nodeGrp==nodeGrp(iN));
nodeGrp(lst) = nGrpN;
% if nB(iN)<3
% error('We should never get here');
% end
end
% if nodeSegN(iN)==0 % remove 6/3/09 since it appears below
% nodeSegN(iN) = nSeg; % and should resolve an issue in
% % nodeGrps()
% nodeSegN(iNstart) = nSeg;
% end
if nodeSegN(iN)==0 % added back 8/3/09 because if had nodeSegN=0 one node from bifurcation
nodeSegN(iN) = nSeg;
end
nE = 1;
nLst = [];
if edgeSegN(eIdx)==0
while nB(iN)==2
nLst(end+1) = iN;
edgeSegN(eIdx) = nSeg;
lstE2 = find(nodeEdges(:,1)==iN | nodeEdges(:,2)==iN);
eIdx = setdiff(lstE2,eIdx);
% correct this to be in um (non-cubical voxels)
eLen = eLen + sum(diff( nodePos(squeeze(nodeEdges(eIdx,:)),:), 1, 1).^2).^0.5;
eLen_um = eLen_um + sum(diff( nodePos_um(squeeze(nodeEdges(eIdx,:)),:), 1, 1).^2).^0.5;
iN = setdiff(unique(nodeEdges(eIdx,:)), iN);
if nodeGrp(iN)==0
nodeGrp(iN) = nGrpN;
elseif nodeGrp(iN) < nGrpN
lst = find(nodeGrp==nGrpN);
nodeGrp(lst) = nodeGrp(iN);
nGrpN = nodeGrp(iN);
elseif nodeGrp(iN) > nGrpN
lst = find(nodeGrp==nodeGrp(iN));
nodeGrp(lst) = nGrpN;
end
if nodeSegN(iN)==0
nodeSegN(iN) = nSeg;
end
nE = nE + 1;
end
nodeSegN(iNstart) = nSeg; % added 6/3/09 when deleted above
nodeSegN(iN) = nSeg; % added 6/3/09 to resolve issue with nodeGrps
% i hope this doesn't cause trouble
edgeSegN(eIdx) = nSeg;
segNedges(nSeg) = nE;
segLen(nSeg) = eLen;
segLen_um(nSeg) = eLen_um;
% if ~isempty(nLst)
% % segDiam(nSeg) = median(nodeDiam(nLst));
% segVesType(nSeg) = median(nodeType(nLst));
% else
% % segDiam(nSeg) = mean(nodeDiam(nodeEdges(eIdx,:)));
% segVesType(nSeg) = max(nodeType(nodeEdges(eIdx,:)));
% end
kk = find(lst3p==iN);
if ~isempty(kk)
segEndNodes(end+1,:) = lst3p([ii kk]);
else
error('We should not get here')
end
nSeg = nSeg + 1;
end
end
end
close(hwait)
% remove groups with zero nodes
nGrpN = 0;
for ii=1:nGrp
lst = find(nodeGrp==ii);
if length(lst)>0
nGrpN = nGrpN + 1;
nodeGrp(lst) = nGrpN;
end
end
nGrp = nGrpN;
% save stats
im.nodeGrp = nodeGrp;
im.nodeSegN = nodeSegN;
im.segNedges = segNedges;
im.segLen = segLen;
im.segLen_um = segLen_um;
% im.segDiam = segDiam;
% im.segVesType = segVesType;
im.edgeSegN = edgeSegN;
im.segEndNodes = segEndNodes;
%im.segNodeMap = lst3p;
%im.nB = nB;
im.segPos = squeeze(mean(reshape(nodePos(im.segEndNodes,:),[2 length(segLen) 3]),1));
%
grpNnodes = [];
disp( sprintf('\n\n# Groups = %d\nGrp\t\tE1\t\tE2\t\tE3\t\tE4+\t\tTotal',nGrp) )
for ii=1:nGrp
lst = find(nodeGrp==ii);
grpNnodes(ii) = length(lst);
im.Stats.grpE1234(ii,1) = length(find(nB(lst)==1));
im.Stats.grpE1234(ii,2) = length(find(nB(lst)==2));
im.Stats.grpE1234(ii,3) = length(find(nB(lst)==3));
im.Stats.grpE1234(ii,4) = length(find(nB(lst)>3));
im.Stats.grpE1234(ii,5) = sum(im.Stats.grpE1234(ii,1:4));
disp(sprintf('%d\t\t',[ii im.Stats.grpE1234(ii,:)]))
end
im.Stats.grpNnodes = grpNnodes;