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
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
|
#include "filter_cluster.h"
#include "utils.h" // for on_grid(...)
typedef STD_list<unsigned int> Cluster;
typedef STD_list<Cluster> ClusterList;
typedef STD_vector<TinyVector<int,3> > NextNeighbours;
typedef STD_list<TinyVector<int,3> > ToDoList;
/////////////////////////////////////////////////////////////////////
void print_clusterlist(const ClusterList& clustelist) {
Log<Filter> odinlog("","print_clusterlist");
unsigned int i=0;
for(ClusterList::const_iterator it=clustelist.begin(); it!=clustelist.end(); ++it) {
ODINLOG(odinlog,normalDebug) << i << " : ";
for(Cluster::const_iterator it2=it->begin(); it2!=it->end(); ++it2) {
ODINLOG(odinlog,normalDebug) << *it2 << " ";
}
ODINLOG(odinlog,normalDebug) << STD_endl;
i++;
}
}
/////////////////////////////////////////////////////////////////////
bool FilterCluster::process(Data<float,4>& data, Protocol& prot) const {
Log<Filter> odinlog(c_label(),"process");
Range all=Range::all();
TinyVector<int,4> shape=data.shape();
TinyVector<int,3> spatshape(shape(sliceDim), shape(phaseDim), shape(readDim));
ODINLOG(odinlog,normalDebug) << "shape=" << shape << STD_endl;
Data<unsigned int, 3> indexarray(spatshape); // linear index in every voxel
for(unsigned int i=0; i<indexarray.size(); i++) {
TinyVector<int,3> index=indexarray.create_index(i);
indexarray(index)=i;
}
NextNeighbours next_neighb;
int n=1;
for(int k=-n; k<=n; k++) {
for(int j=-n; j<=n; j++) {
for(int i=-n; i<=n; i++) {
int abssum=abs(i)+abs(j)+abs(k);
if(abssum<=1) { // only adjacent voxels (with surface-to-surface contact)
next_neighb.push_back(TinyVector<int,3>( k, j, i));
}
}
}
}
unsigned int nn=next_neighb.size();
Data<float,4> outdata(shape);
outdata=0.0;
for(int irep=0; irep<shape(timeDim); irep++) {
ClusterList clusterlist;
Data<float,3> voldata(spatshape);
voldata(all,all,all)=data(irep,all,all,all);
Data<bool,3> done(spatshape);
done=false;
ToDoList todo;
TinyVector<int,3> seedindex;
TinyVector<int,3> centerindex;
TinyVector<int,3> neighbindex;
for(unsigned int i=0; i<voldata.size(); i++) {
seedindex=voldata.create_index(i);
if(voldata(seedindex)>0.0 && !done(seedindex)) {
todo.push_back(seedindex);
Cluster indexlist;
while(todo.size()) {
ToDoList nexttodo;
for(ToDoList::iterator it=todo.begin(); it!=todo.end(); ++it) {
centerindex=(*it);
indexlist.push_back(indexarray(centerindex));
// find next neighbours
for(unsigned int in=0; in<nn; in++) {
neighbindex=centerindex+next_neighb[in];
if(on_grid<3>(spatshape, neighbindex)) {
if(voldata(neighbindex)>0.0 && !done(neighbindex)) {
nexttodo.push_back(neighbindex);
done(neighbindex)=true;
}
}
}
}
todo=nexttodo;
ODINLOG(odinlog,normalDebug) << "todo.size=" << todo.size() << STD_endl;
}
clusterlist.push_back(indexlist); // add new cluster
}
}
// print_clusterlist(clusterlist);
bool cont=true;
int nclusters=clusterlist.size();
while(cont) {
ODINLOG(odinlog,normalDebug) << nclusters << " clusters" << STD_endl;
cont=false;
for(ClusterList::iterator it=clusterlist.begin(); it!=clusterlist.end(); ++it) {
// sort and make unique before going through the elements
it->sort();
it->unique();
ODINLOG(odinlog,normalDebug) << "clustersize=" << it->size() << STD_endl;
for(Cluster::const_iterator it2=it->begin(); it2!=it->end(); ++it2) {
unsigned int i=*it2; // find this index in following clusters
ODINLOG(odinlog,normalDebug) << "finding index " << i << STD_endl;
ClusterList::iterator nextcluster=it; ++nextcluster;
for(ClusterList::iterator itf=nextcluster; itf!=clusterlist.end(); ++itf) {
if(itf->size()) { // skip empty clusters
bool found_in_cluster=false;
for(Cluster::const_iterator it3=itf->begin(); it3!=itf->end(); ++it3) {
if(i==(*it3)) {
found_in_cluster=true;
break;
}
}
if(found_in_cluster) {
it->merge(*itf); // sort into list, itf will be emptied
it->unique();
cont=true; // reiterating clustering while no mergers left
nclusters--;
ODINLOG(odinlog,normalDebug) << "found " << i << " ..." << STD_endl;
} else {
ODINLOG(odinlog,normalDebug) << "not found " << i << " ..." << STD_endl;
}
}
}
}
}
}
ClusterList clusterlist_sort;
for(ClusterList::iterator it=clusterlist.begin(); it!=clusterlist.end(); ++it) {
if(it->size()) {
Cluster indexlist;
indexlist.push_back(it->size()); // store size of cluster as first element for sorting
for(Cluster::const_iterator it2=it->begin(); it2!=it->end(); ++it2) {
indexlist.push_back(*it2);
}
clusterlist_sort.push_back(indexlist);
}
}
clusterlist_sort.sort();
// print_clusterlist(clusterlist_sort);
unsigned int ncluster=clusterlist_sort.size();
unsigned int iclust=0;
for(ClusterList::iterator it=clusterlist_sort.begin(); it!=clusterlist_sort.end(); ++it) {
for(Cluster::const_iterator it2=it->begin(); it2!=it->end(); ++it2) {
if(it2!=it->begin()) { // omit 1st element which contains list size
TinyVector<int,3> spatindex=voldata.create_index(*it2);
outdata(irep,spatindex(0),spatindex(1),spatindex(2))=ncluster-iclust;
}
}
iclust++;
}
}
data.reference(outdata);
return true;
}
|
