A.I.M.S algorithms


meshcc_d.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 
35 
36 #ifndef AIMS_MESH_CONNECTIVITY_MESHCC_D_H
37 #define AIMS_MESH_CONNECTIVITY_MESHCC_D_H
38 
41 #include <cstdlib>
42 #include <map>
43 #include <set>
44 #include <stack>
45 #include <cstdlib>
46 //#include <stdio.h>
47 
48 
49 template<class T>
51  const Texture<T> & inittex, T threshold, int mode )
52 {
53  Texture<T> tex;
54  const std::vector<Point3df> & vert = mesh.vertex();
55  const std::vector< AimsVector<uint,3> > & poly = mesh.polygon();
56  unsigned i, n = vert.size();
57 
58  ASSERT( inittex.nItem() == n );
59  tex.reserve( n );
60 
61  // neighbours map
62 
63  std::map<unsigned, std::set<unsigned> > neighbours;
64  unsigned v1, v2, v3;
65  int nnode = 0;
66 
67  // init texture
68  std::cout << "Thresholding: " << threshold << std::endl;
69  if (mode==1)
70  {
71  for( i=0; i<n; ++i )
72  {
73  if(inittex.item(i)>=threshold) tex.push_back( FORBIDDEN );
74  else
75  {
76  tex.push_back( 0 );
77  nnode++;
78  }
79  }
80  }
81  else if (mode==0)
82  {
83  for( i=0; i<n; ++i )
84  {
85  if(inittex.item(i)<=threshold) tex.push_back( FORBIDDEN );
86  else
87  {
88  tex.push_back( 0 );
89  nnode++;
90  }
91  }
92  }
93  else
94  {
95  std::cerr << "AimsMeshLabelConnectedComponent : thresholding mode unknown" << std::endl;
96  exit(1);
97  }
98  std::cout << nnode << "/" << n << std::flush;
99 
100  //Detect connectivity
101  for( i=0; i<poly.size(); ++i )
102  {
103  v1 = poly[i][0];
104  v2 = poly[i][1];
105  v3 = poly[i][2];
106  if(tex.item(v1)!=FORBIDDEN
107  && tex.item(v2)!=FORBIDDEN)
108  {
109  neighbours[v1].insert( v2 );
110  neighbours[v2].insert( v1 );
111  }
112  if(tex.item(v1)!=FORBIDDEN
113  && tex.item(v3)!=FORBIDDEN)
114  {
115  neighbours[v1].insert( v3 );
116  neighbours[v3].insert( v1 );
117  }
118  if(tex.item(v2)!=FORBIDDEN
119  && tex.item(v3)!=FORBIDDEN)
120  {
121  neighbours[v2].insert( v3 );
122  neighbours[v3].insert( v2 );
123  }
124  }
125 
126  T label = 0;
127  unsigned point;
128  std::stack<unsigned> current;
129  std::set<unsigned>::iterator in, fn;
130 
131 
132  printf("Computing connected component\n");
133  fflush(stdout);
134  for( i=0; i<n; ++i )
135  {
136  if(tex.item(i)==0)
137  {
138  label++;
139  current.push(i);
140  fflush(stdout);
141  while(!current.empty())
142  {
143  point = current.top();
144  current.pop();
145  tex.item(point)=label;
146  for( in=neighbours[point].begin(), fn=neighbours[point].end(); in!=fn; ++in )
147  {
148  if(tex.item(*in)==0) current.push(*in);
149  }
150  }
151  }
152  }
153  std::cout << "nb cc: " << label << std::endl;
154 
155  return tex;
156 }
157 
158 
159 template<class T>
161  const Texture<T> & inittex, const T lab )
162 {
163  Texture<T> tex;
164  const std::vector<Point3df> & vert = mesh.vertex();
165  const std::vector< AimsVector<uint,3> > & poly = mesh.polygon();
166  unsigned i, n = vert.size();
167 
168  ASSERT( inittex.nItem() == n );
169  tex.reserve( n );
170 
171  // init texture
172  for (i=0; i<n; ++i)
173  if (inittex.item(i) == lab)
174  tex.push_back(0);
175  else
176  tex.push_back(FORBIDDEN);
177 
178  //Detect connectivity
179  std::map<unsigned, std::set<unsigned> > neighbours;
180  unsigned v1, v2, v3;
181  for( i=0; i<poly.size(); ++i )
182  {
183  v1 = poly[i][0];
184  v2 = poly[i][1];
185  v3 = poly[i][2];
186  if(tex.item(v1)!=FORBIDDEN
187  && tex.item(v2)!=FORBIDDEN
188  && inittex.item(v1) == inittex.item(v2) )
189 
190  {
191  neighbours[v1].insert( v2 );
192  neighbours[v2].insert( v1 );
193  }
194  if(tex.item(v1)!=FORBIDDEN
195  && tex.item(v3)!=FORBIDDEN
196  && inittex.item(v1) == inittex.item(v3) )
197  {
198  neighbours[v1].insert( v3 );
199  neighbours[v3].insert( v1 );
200  }
201  if(tex.item(v2)!=FORBIDDEN
202  && tex.item(v3)!=FORBIDDEN
203  && inittex.item(v2) == inittex.item(v3) )
204  {
205  neighbours[v2].insert( v3 );
206  neighbours[v3].insert( v2 );
207  }
208  }
209 
210  T label = 0;
211  unsigned inc =0;
212  unsigned point;
213  std::stack<unsigned> current;
214  std::set<unsigned>::iterator in, fn;
215 
216 
217  for( i=0; i<n; ++i )
218  {
219  if(tex.item(i)==0)
220  {
221  label++;
222  ++inc;
223  current.push(i);
224  while(!current.empty())
225  {
226  point = current.top();
227  current.pop();
228  tex.item(point)=label;
229  for( in=neighbours[point].begin(), fn=neighbours[point].end(); in!=fn; ++in )
230  {
231  if(tex.item(*in)==0) current.push(*in);
232  }
233  }
234  }
235  }
236 
237  return inc;
238 }
239 
240 
241 // Give a connected map of inittex.
242 // An initial connex area composed of different labels
243 // is splitted in sub-area of same labels
244 template<class T>
246  const Texture<T> & inittex, T threshold )
247 {
248  Texture<T> tex;
249  const std::vector<Point3df> & vert = mesh.vertex();
250  const std::vector< AimsVector<uint,3> > & poly = mesh.polygon();
251  unsigned i, n = vert.size();
252 
253  ASSERT( inittex.nItem() == n );
254  tex.reserve( n );
255 
256  // neighbours map
257 
258  std::map<unsigned, std::set<unsigned> > neighbours;
259  unsigned v1, v2, v3;
260 
261 
262  // init texture
263  for( i=0; i<n; ++i )
264  {
265  if(inittex.item(i) < threshold)
266  tex.push_back( FORBIDDEN );
267  else
268  tex.push_back( 0 );
269  }
270 
271  //Detect connectivity
272  for( i=0; i<poly.size(); ++i )
273  {
274  v1 = poly[i][0];
275  v2 = poly[i][1];
276  v3 = poly[i][2];
277  if(tex.item(v1)!=FORBIDDEN
278  && tex.item(v2)!=FORBIDDEN
279  && inittex.item(v1) == inittex.item(v2) )
280 
281  {
282  neighbours[v1].insert( v2 );
283  neighbours[v2].insert( v1 );
284  }
285  if(tex.item(v1)!=FORBIDDEN
286  && tex.item(v3)!=FORBIDDEN
287  && inittex.item(v1) == inittex.item(v3) )
288  {
289  neighbours[v1].insert( v3 );
290  neighbours[v3].insert( v1 );
291  }
292  if(tex.item(v2)!=FORBIDDEN
293  && tex.item(v3)!=FORBIDDEN
294  && inittex.item(v2) == inittex.item(v3) )
295  {
296  neighbours[v2].insert( v3 );
297  neighbours[v3].insert( v2 );
298  }
299  }
300 
301  T label = 0;
302  unsigned point;
303  std::stack<unsigned> current;
304  std::set<unsigned>::iterator in, fn;
305 
306 
307  for( i=0; i<n; ++i )
308  if(tex.item(i) == 0)
309  {
310  label++;
311  current.push(i);
312  while(!current.empty())
313  {
314  point = current.top();
315  current.pop();
316  tex.item(point)=label;
317  for( in=neighbours[point].begin(), fn=neighbours[point].end(); in!=fn; ++in )
318  {
319  if(tex.item(*in)==0)
320  {
321  current.push(*in);
322  }
323  }
324  }
325  }
326 
327  std::cout << "Nb of cc : "<< label << std::endl;
328 
329  for (i=0;i<n;++i)
330  if (tex.item(i) == FORBIDDEN)
331  tex.item(i) = 0;
332 
333  return tex;
334 }
335 
336 template<class T>
338  const Texture<T> & inittex, T threshold )
339 {
340  Texture<T> tex;
341  const std::vector<Point3df> & vert = mesh.vertex();
342  const std::vector< AimsVector<uint,3> > & poly = mesh.polygon();
343  unsigned i, n = vert.size();
344 
345  std::map<T,std::set<unsigned> > labels;
346  std::map<unsigned,T> nbLabels;
347 
348  ASSERT( inittex.nItem() == n );
349  tex.reserve( n );
350 
351  // neighbours map
352 
353  std::map<unsigned, std::set<unsigned> > neighbours;
354  unsigned v1, v2, v3;
355 
356 
357  // init texture
358  for( i=0; i<n; ++i )
359  {
360  if(inittex.item(i) != threshold)
361  tex.push_back( FORBIDDEN );
362  else
363  tex.push_back( 0 );
364  }
365 
366  //Detect connectivity
367  for( i=0; i<poly.size(); ++i )
368  {
369  v1 = poly[i][0];
370  v2 = poly[i][1];
371  v3 = poly[i][2];
372  if(tex.item(v1)!=FORBIDDEN
373  && tex.item(v2)!=FORBIDDEN )
374 
375  {
376  neighbours[v1].insert( v2 );
377  neighbours[v2].insert( v1 );
378  }
379  if(tex.item(v1)!=FORBIDDEN
380  && tex.item(v3)!=FORBIDDEN )
381  {
382  neighbours[v1].insert( v3 );
383  neighbours[v3].insert( v1 );
384  }
385  if(tex.item(v2)!=FORBIDDEN
386  && tex.item(v3)!=FORBIDDEN )
387  {
388  neighbours[v2].insert( v3 );
389  neighbours[v3].insert( v2 );
390  }
391  }
392 
393  T label = 0;
394  unsigned point;
395  std::stack<unsigned> current;
396  std::set<unsigned>::iterator in, fn;
397 
398 
399  for( i=0; i<n; ++i )
400  if(tex.item(i) == 0)
401  {
402  label++;
403  current.push(i);
404  while(!current.empty())
405  {
406  point = current.top();
407  current.pop();
408  tex.item(point)=label;
409  labels[label].insert(point);
410  for( in=neighbours[point].begin(), fn=neighbours[point].end(); in!=fn; ++in )
411  {
412  if(tex.item(*in)==0)
413  current.push(*in);
414  }
415  }
416  nbLabels[ labels[label].size() ] = label;
417  std::cout << "The " << label << " th cc (" << threshold << ")has " << labels[label].size() << " points " << std::endl;
418  }
419 
420  std::cout << "Nb of cc : "<< label << std::endl;
421 
422 
423  T maxLab = nbLabels.rbegin()->first;
424  std::cout << "maxlab= " << maxLab << "-> " << nbLabels[maxLab] << std::endl;
425 
426 
427  for (i=0;i<n;++i)
428  {
429  if (tex.item(i) == FORBIDDEN )
430  tex.item(i) = inittex.item(i);
431  else
432  if (tex.item(i) == nbLabels[maxLab])
433  tex.item(i) = threshold;
434  else
435  tex.item(i) = 0;
436  }
437 
438  return tex;
439 }
440 
441 
442 #endif
const T & item(int n) const
const int FORBIDDEN
hum, should not be there as a global, non-namespaced variable...
Definition: meshcc.h:65
void reserve(size_t size)
unsigned AimsMeshLabelNbConnectedComponent(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, const T lab)
Definition: meshcc_d.h:160
const std::vector< Point3df > & vertex() const
const std::vector< AimsVector< uint, D > > & polygon() const
Texture< T > AimsMeshFilterConnectedComponent(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, T threshold)
Definition: meshcc_d.h:337
Texture< T > AimsMeshLabelConnectedComponent2Texture(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, T threshold)
Definition: meshcc_d.h:245
Texture< T > AimsMeshLabelConnectedComponent(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, T threshold, int mode)
Definition: meshcc_d.h:50
#define ASSERT(EX)
size_t nItem() const
void push_back(const T &item)