Moka libraries
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
compute-homology.cc
Go to the documentation of this file.
1 /*
2  * lib-homology : Computation of homology generators.
3  * Copyright (C) 2010, Samuel Peltier, Université de Poitiers, Laboratoire SIC-XLIM
4  * http://www.sic.sp2mi.univ-poitiers.fr/
5  * Copyright (C) 2010, Guillaume Damiand, CNRS, LIRIS,
6  * guillaume.damiand@liris.cnrs.fr, http://liris.cnrs.fr/
7  *
8  * This file is part of lib-homology
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU Lesser General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU Lesser General Public License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public License
21  * along with this program. If not, see <http://www.gnu.org/licenses/>.
22  */
23 
24 //******************************************************************************
25 #include "g-map-vertex.hh"
26 #include "MatricePMQ.hh"
27 #include "compute-homology.hh"
28 #include <cassert>
29 using namespace GMap3d;
30 //******************************************************************************
31 CHomology::CHomology(CGMapVertex* AMap, int AMark) :
32  FMap(AMap),
33  FMark(AMark),
34  FNbVertices(0),
35  FNbEdges(0),
36  FNbFaces(0),
37  FNbVolumes(0),
38  FShowH0(false),
39  FShowH1free(true),
40  FShowH1torsion(true),
41  FShowH2free(false),
42  FShowH2torsion(false),
43  FShowH3(false)
44 {
45  // std::cout<<"Marks: "<<AMark1<<", "<<AMark2<<std::endl;
46  for (int i=0;i<3;++i)
47  FMatrix[i]=NULL;
48 }
49 //******************************************************************************
50 CHomology::~CHomology()
51 {
52  for (int i=0;i<3;++i)
53  delete FMatrix[i];
54 }
55 //******************************************************************************
56 unsigned long CHomology::size() const
57 {
58  unsigned long res = 0;
59  if ( FMatrix[0]!=NULL ) res += FMatrix[0]->size() + FCells[0].size();
60  if ( FMatrix[1]!=NULL ) res += FMatrix[1]->size() + FCells[1].size();
61  if ( FMatrix[2]!=NULL ) res += FMatrix[2]->size() + FCells[2].size();
62  return res;
63 }
64 //******************************************************************************
65 int CHomology::computeIncidenceNumber(CDart* ADart, int ADim,
66  CDart* ADart2, int AIndex)
67 {
68  assert((0 <= ADim) && (ADim < 3));
69 
70  // 1) marquer m1 la i cell orientée incidence à ADart
71  // 2) marquer m2 la i+1 cell orientée incidence à ADart2
72 
73  // 2) Parcourir a01,...,a0i incidente à d2
74 
75 
76  // si brin non traité
77  // si brin marqué m1 alors +1
78  // sinon -1
79  // marquer traité orbite a0,...,i-1
80 
81  int incidenceNumber=0;
82 
83  int m2oriented = FMap->getNewMark();
84  if ( ADim==0 )
85  {
86  FMap->setMark(ADart2,m2oriented);
87  FMap->setMark(ADart2->getAlpha(0),m2oriented);
88  }
89  else
90  {
91  FMap->halfMarkOrbit(ADart2, ORBIT_INF[ADim+1], m2oriented);
92  }
93 
94  int treated = FMap->getNewMark();
95 
96  CCoverage * cov = FMap->getDynamicCoverage(ADart, ORBIT_CELL[ADim]);
97  for (; cov->cont(); ++(*cov))
98  {
99  if ( !FMap->isMarked(**cov,treated) && FMap->isMarked(**cov,m2oriented) )
100  {
101  if ( (ADim>0 && (long int)FMap->getDirectInfo(**cov,AIndex)>0) ||
102  (ADim==0 && **cov==ADart2) )
103  {
104  // std::cout<<"inc incidenceNumber"<<std::endl;
105  ++incidenceNumber;
106  }
107  else
108  {
109  // std::cout<<"dec incidenceNumber"<<std::endl;
110  --incidenceNumber;
111  }
112 
113  FMap-> markOrbit(**cov,ORBIT_INF[ADim],treated);
114  }
115  }
116 
117  for (cov->reinit(); cov->cont(); ++(*cov))
118  {
119  if ( FMap->isMarked(**cov,treated) )
120  {
121  FMap->unmarkOrbit(**cov,ORBIT_INF[ADim],treated);
122  }
123  }
124  delete cov;
125 
126  if ( ADim==0 )
127  {
128  FMap->unsetMark(ADart2,m2oriented);
129  FMap->unsetMark(ADart2->getAlpha(0),m2oriented);
130  }
131  else
132  {
133  FMap->halfUnmarkOrbit(ADart2, ORBIT_INF[ADim+1], m2oriented);
134  }
135 
136  assert(FMap->isWholeMapUnmarked(treated));
137  assert(FMap->isWholeMapUnmarked(m2oriented));
138 
139  // Libérations:
140  FMap->freeMark(treated);
141  FMap->freeMark(m2oriented);
142 
143  return incidenceNumber;
144 }
145 //******************************************************************************
146 void CHomology::computeIncidence(int ADim, bool AComputeNextCells)
147 {
148  assert( 0<=ADim && ADim<3 );
149 
150  long int currentcell = 1;
151  long int number = 0;
152 
153  int treated = FMap->getNewMark();
154  int positive = FMap->getNewMark();
155  int index = FMap->getNewDirectInfo();
156 
157  int treated2 = -1;
158  long int currentcell2 = 1;
159  if ( AComputeNextCells ) treated2 = FMap->getNewMark();
160 
161  // 1) We number each (ADim)-cell of the map
162  CDynamicCoverageAll it(FMap);
163  for (; it.cont(); ++it)
164  {
165  if (!FMap->isMarked(*it, treated))
166  {
167  if ( ADim==1 ) FMap->markOrbit(*it,ORBIT_23,positive);
168  else if ( ADim==2 )
169  { FMap->halfMarkOrbit(*it,ORBIT_01,positive); }
170 
171  FCells[ADim][currentcell-1]=*it;
172  // if ( ADim==0 )
173  // std::cout<<"sommet["<<currentcell-1<<"] : "<<*FMap->findVertex(*it)<<std::endl;
174 
175  CCoverage* it2 = FMap->getDynamicCoverage(*it,ORBIT_CELL[ADim]);
176  for(;it2->cont();++(*it2))
177  {
178  if ( ADim==1 && !FMap->isMarked(**it2,positive) )
179  number = -currentcell;
180  else if (ADim==2 && !FMap->isMarked(**it2,positive) &&
181  !FMap->isMarked((**it2)->getAlpha(3),positive))
182  number = -currentcell;
183  else
184  number = currentcell;
185 
186  //std::cout<<"Number="<<number<<std::endl;
187 
188  FMap->setDirectInfo(**it2,index,(void*)number);
189  FMap->setMark(**it2, treated);
190  }
191 
192  if ( ADim==1 ) FMap->unmarkOrbit(*it,ORBIT_23,positive);
193  else if ( ADim==2 )
194  { FMap->halfUnmarkOrbit(*it,ORBIT_01,positive); }
195 
196  delete it2;
197  ++currentcell;
198  }
199  if ( treated2!=-1 && !FMap->isMarked(*it, treated2) )
200  {
201  FCells[ADim+1][currentcell2-1]=*it;
202  FMap->markOrbit(*it, ORBIT_CELL[ADim+1],treated2);
203  ++currentcell2;
204  }
205  }
206 
207  FMap->negateMaskMark(treated);
208  if ( treated2!=-1 )
209  {
210  FMap->negateMaskMark(treated2);
211  FMap->freeMark(treated2);
212  }
213 
214  // 2) We run through all the (ADim+1)-cell of the map and compute the incidence
215  // number with ADim-cells
216  currentcell = 0;
217  for ( it.reinit(); it.cont(); ++it )
218  {
219  if ( !FMap->isMarked(*it, treated) )
220  {
221  CCoverage* it2 = FMap->getDynamicCoverage(*it,ORBIT_CELL[ADim+1]);
222  for( ;it2->cont();++(*it2) )
223  {
224  int i = (long int)FMap->getDirectInfo(**it2, index);
225  int icorrected = i;
226 
227  if ( i<0 ) icorrected = -i;
228  --icorrected;
229 
230  if ( FMatrix[ADim]->getValPMQ( icorrected, currentcell)==0 )
231  {
232  int val = computeIncidenceNumber(**it2, ADim, *it, index);
233 
234  // std::cout<<"Val="<<val<<", puis ";
235 
236  //if ( i<0 ) { val = -val; }
237 
238  // std::cout<<"("<<i<<","<<currentcell<<")="<<val<<std::endl;
239 
240  // if ( val!=0 ) val = 1;
241  FMatrix[ADim]->setValPMQ( icorrected, currentcell, val);
242  }
243 
244  FMap->setMark(**it2, treated);
245  }
246  delete it2;
247  ++currentcell;
248  // std::cout<<"Next cell : "<<currentcell<<std::endl;
249  }
250  }
251 
252  FMap->negateMaskMark(treated);
253 
254  assert(FMap->isWholeMapUnmarked(treated));
255 
256  FMap->freeMark(treated);
257  FMap->freeDirectInfo(index);
258 }
259 //******************************************************************************
260 bool CHomology::computeSurfacicHomology()
261 {
262  FMap->countCells(-1,&FNbVertices,&FNbEdges,&FNbFaces,NULL,NULL,NULL);
263 
264  delete FMatrix[0]; FMatrix[0] = new MatricePMQ(FNbVertices, FNbEdges);
265  delete FMatrix[1]; FMatrix[1] = new MatricePMQ(FNbEdges, FNbFaces);
266  delete FMatrix[2]; FMatrix[2] = NULL;
267 
268  FCells[0].clear(); FCells[0].reserve(FNbVertices);
269  FCells[1].clear(); FCells[1].reserve(FNbEdges);
270  FCells[2].clear(); FCells[2].reserve(FNbFaces);
271 
272  if (FCells[0].capacity()<FNbVertices ||
273  FCells[1].capacity()<FNbEdges ||
274  FCells[2].capacity()<FNbFaces ||
275  !FMatrix[0]->valid() ||
276  !FMatrix[1]->valid())
277  {
278  FCells[0].clear(); FCells[1].clear(); FCells[2].clear();
279  delete FMatrix[0]; delete FMatrix[1];
280  FMatrix[0]=NULL; FMatrix[1]=NULL;
281  return false;
282  }
283 
284  computeIncidence(0);
285  computeIncidence(1,true);
286 
287  FMatrix[0]->smithForm();
288 
289  FMatrix[1]->getM()->multGauche(FMatrix[0]->getQinv());
290  FMatrix[1]->getP()->setMatrice( FMatrix[0]->getQ());
291  FMatrix[1]->getPinv()->setMatrice(FMatrix[0]->getQinv());
292  FMatrix[1]->smithForm();
293 
294  FNbCycleDim0 = FMatrix[0]->getM()->getnbli();
295  FNbCycleDim1 = FMatrix[0]->getM()->nbCycle();
296  FNbCycleDim2 = FMatrix[1]->getM()->nbCycle();
297  FNbCycleDim3 = 0;
298 
299  FNbBordFaibleDim0 = FMatrix[0]->getM()->getnbcol()-FNbCycleDim1;
300  FNbBordFaibleDim1 = FMatrix[1]->getM()->getnbcol()-FNbCycleDim2;
301  FNbBordFaibleDim2 = 0;
302 
303  FNbGenLibreDim0 = FNbCycleDim0 - FNbBordFaibleDim0;
304  FNbGenLibreDim1 = FNbCycleDim1 - FNbBordFaibleDim1;
305  FNbGenLibreDim2 = FNbCycleDim2;
306  FNbGenLibreDim3 = 0;
307 
308  FNbGenTorsionDim1 = FMatrix[1]->getM()->nbTorsion();
309  FNbGenTorsionDim2 = 0;
310 
311  updateSelectedDarts();
312 
313  return true;
314 }
315 //******************************************************************************
316 bool CHomology::computeVolumicHomology()
317 {
318  FMap->countCells(-1,&FNbVertices,&FNbEdges,&FNbFaces,&FNbVolumes,NULL,NULL);
319 
320  delete FMatrix[0]; FMatrix[0] = new MatricePMQ(FNbVertices, FNbEdges);
321  delete FMatrix[1]; FMatrix[1] = new MatricePMQ(FNbEdges, FNbFaces);
322  delete FMatrix[2]; FMatrix[2] = new MatricePMQ(FNbFaces, FNbVolumes);
323 
324  FCells[0].clear(); FCells[0].reserve(FNbVertices);
325  FCells[1].clear(); FCells[1].reserve(FNbEdges);
326  FCells[2].clear(); FCells[2].reserve(FNbFaces);
327 
328  if (FCells[0].capacity()<FNbVertices ||
329  FCells[1].capacity()<FNbEdges ||
330  FCells[2].capacity()<FNbFaces ||
331  !FMatrix[0]->valid() ||
332  !FMatrix[1]->valid() ||
333  !FMatrix[2]->valid())
334  {
335  FCells[0].clear(); FCells[1].clear(); FCells[2].clear();
336  delete FMatrix[0]; delete FMatrix[1]; delete FMatrix[2];
337  FMatrix[0]=NULL; FMatrix[1]=NULL; FMatrix[2]=NULL;
338  return false;
339  }
340 
341  computeIncidence(0);
342  computeIncidence(1);
343  computeIncidence(2);
344 
345  FMatrix[0]->smithForm();
346 
347 // std::cout<<"FMatrix[1]->getM() avant Smith"<<std::endl;
348 // FMatrix[1]->getM()->affiche();
349 
350  FMatrix[1]->getM()->multGauche(FMatrix[0]->getQinv());
351  FMatrix[1]->getP()->setMatrice( FMatrix[0]->getQ());
352  FMatrix[1]->getPinv()->setMatrice(FMatrix[0]->getQinv());
353  FMatrix[1]->smithForm();
354 
355 // std::cout<<"FMatrix[1]->getM() après Smith"<<std::endl;
356 // FMatrix[1]->getM()->affiche();
357 
358 // std::cout<<"FMatrix[2]->getM() avant Smith"<<std::endl;
359 // FMatrix[2]->getM()->affiche();
360 
361  FMatrix[2]->getM()->multGauche(FMatrix[1]->getQinv());
362  FMatrix[2]->getP()->setMatrice(FMatrix[1]->getQ());
363  FMatrix[2]->getPinv()->setMatrice(FMatrix[1]->getQinv());
364  FMatrix[2]->smithForm();
365 
366 // std::cout<<"FMatrix[2]->getM() après Smith"<<std::endl;
367 // FMatrix[2]->getM()->affiche();
368 
369  FNbCycleDim0 = FMatrix[0]->getM()->getnbli();
370  FNbCycleDim1 = FMatrix[0]->getM()->nbCycle();
371  FNbCycleDim2 = FMatrix[1]->getM()->nbCycle();
372  FNbCycleDim3 = FMatrix[2]->getM()->nbCycle();
373 
374  FNbBordFaibleDim0 = FMatrix[0]->getM()->getnbcol()-FNbCycleDim1;
375  FNbBordFaibleDim1 = FMatrix[1]->getM()->getnbcol()-FNbCycleDim2;
376  FNbBordFaibleDim2 = FMatrix[2]->getM()->getnbcol()-FNbCycleDim3;
377 
378  FNbGenLibreDim0 = FNbCycleDim0 - FNbBordFaibleDim0;
379  FNbGenLibreDim1 = FNbCycleDim1 - FNbBordFaibleDim1;
380  FNbGenLibreDim2 = FNbCycleDim2 - FNbBordFaibleDim2;
381  FNbGenLibreDim3 = FNbCycleDim3;
382 
383  FNbGenTorsionDim1 = FMatrix[1]->getM()->nbTorsion();
384  FNbGenTorsionDim2 = FMatrix[2]->getM()->nbTorsion();
385 
386  updateSelectedDarts();
387 
388  return true;
389 }
390 
391 //******************************************************************************
392 unsigned int CHomology::getH0FreeGenerators()
393 { return FNbGenLibreDim0; }
394 //------------------------------------------------------------------------------
395 unsigned int CHomology::getH1FreeGenerators()
396 { return FNbGenLibreDim1; }
397 //------------------------------------------------------------------------------
398 unsigned int CHomology::getH1TorsionGenerators()
399 { return FNbGenTorsionDim1; }
400 //------------------------------------------------------------------------------
401 unsigned int CHomology::getH2FreeGenerators()
402 { return FNbGenLibreDim2; }
403 //------------------------------------------------------------------------------
404 unsigned int CHomology::getH2TorsionGenerators()
405 { return FNbGenTorsionDim2; }
406 //------------------------------------------------------------------------------
407 unsigned int CHomology::getH3FreeGenerators()
408 { return FNbGenLibreDim3; }
409 //******************************************************************************
410 bool CHomology::getShowH0() const
411 { return FShowH0; }
412 bool CHomology::getShowH1free() const
413 { return FShowH1free; }
414 bool CHomology::getShowH1torsion() const
415 { return FShowH1torsion; }
416 bool CHomology::getShowH2free() const
417 { return FShowH2free; }
418 bool CHomology::getShowH2torsion() const
419 { return FShowH2torsion; }
420 bool CHomology::getShowH3() const
421 { return FShowH3; }
422 //******************************************************************************
423 void CHomology::setShowH0(bool AValue)
424 {
425  if ( AValue!=FShowH0 )
426  {
427  FShowH0=AValue;
428  updateSelectedDarts();
429  }
430 }
431 //------------------------------------------------------------------------------
432 void CHomology::setShowH1free(bool AValue)
433 {
434  if ( AValue!=FShowH1free )
435  {
436  FShowH1free=AValue;
437  updateSelectedDarts();
438  }
439 }
440 //------------------------------------------------------------------------------
441 void CHomology::setShowH1torsion(bool AValue)
442 {
443  if ( AValue!=FShowH1torsion )
444  {
445  FShowH1torsion=AValue;
446  updateSelectedDarts();
447  }
448 }
449 //------------------------------------------------------------------------------
450 void CHomology::setShowH2free(bool AValue)
451 {
452  if ( AValue!=FShowH2free )
453  {
454  FShowH2free=AValue;
455  updateSelectedDarts();
456  }
457 }
458 //------------------------------------------------------------------------------
459 void CHomology::setShowH2torsion(bool AValue)
460 {
461  if ( AValue!=FShowH2torsion )
462  {
463  FShowH2torsion=AValue;
464  updateSelectedDarts();
465  }
466 }
467 //------------------------------------------------------------------------------
468 void CHomology::setShowH3(bool AValue)
469 {
470  if ( AValue!=FShowH3 )
471  {
472  FShowH3=AValue;
473  updateSelectedDarts();
474  }
475 }
476 //******************************************************************************
477 void CHomology::updateSelectedDarts()
478 {
479  if ( FMatrix[0]==NULL || FMark==-1 ) return;
480 
481  // 1) We unmark all the darts
482  for (CDynamicCoverageAll it(FMap); it.cont(); ++it)
483  FMap->unsetMark(*it,FMark);
484 
485  // 2) we mark the required generators.
486  if ( FShowH0 )
487  {
488  //marquage des cellules faisant partie des H0
489  for (int j=FNbBordFaibleDim0;j<FNbBordFaibleDim0+FNbGenLibreDim0;++j)
490  {
491  for (int i=0;i<FMatrix[0]->getP()->getnbli();++i)
492  {
493  if (FMatrix[0]->getP()->getVal(i,j)!=0)
494  {
495  //marquerLibre la d cellule numero i
496  FMap->markOrbit(FCells[0][i],ORBIT_VERTEX,FMark);
497  // std::cout<<"Generateur libre "<<j<<" - Cellule libre: "<<i<<" ("<<FCells[1][i]
498  // <<")"<<std::endl;
499  }
500  }
501  }
502  }
503 
504  if ( FShowH1free )
505  {
506  //marquage des cellules faisant partie des générateurs libres
507  //int deb_libre=nb_bf;
508  //int fin_libre=nb_bf+nb_l;
509  for (int j=FNbBordFaibleDim1;j<FNbBordFaibleDim1+FNbGenLibreDim1;++j)
510  {
511  for (int i=0;i<FMatrix[1]->getP()->getnbli();++i)
512  {
513  if(FMatrix[1]->getP()->getVal(i,j)!=0)
514  {
515  //marquerLibre la d cellule numero i
516  FMap->markOrbit(FCells[1][i],ORBIT_EDGE,FMark);
517  // std::cout<<"Generateur libre "<<j<<" - Cellule libre: "<<i<<" ("<<FCells[1][i]
518  // <<")"<<std::endl;
519  }
520  }
521  }
522  }
523 
524  if ( FShowH1torsion )
525  {
526  //marquage des cellules faisant partie des générateurs de torsion
527  for (int j=0;j<FNbGenTorsionDim1;j++){
528  for (int i=0;i<FMatrix[1]->getP()->getnbli();++i)
529  {
530  if(FMatrix[1]->getP()->getVal(i,j)!=0)
531  {
532  //marquerTorsion la d cellule numero i
533  FMap->markOrbit(FCells[1][i],ORBIT_EDGE,FMark);
534  // std::cout<<"Cellule torsion: "<<i<<" ("<<FCells[1][i]
535  // <<")"<<std::endl;
536  }
537  }
538  }
539  }
540 
541  if ( FShowH2free )
542  {
543  if ( FMatrix[2]!=NULL )
544  {
545  //marquage des cellules faisant partie des H2 libres
546  for (int j=FNbBordFaibleDim2;j<FNbBordFaibleDim2+FNbGenLibreDim2;++j)
547  {
548  for (int i=0;i<FMatrix[2]->getP()->getnbli();++i)
549  {
550  if(FMatrix[2]->getP()->getVal(i,j)!=0)
551  {
552  //marquerLibre la d cellule numero i
553  FMap->markOrbit(FCells[2][i],ORBIT_FACE,FMark);
554  // std::cout<<"Generateur libre "<<j<<" - Cellule libre: "<<i<<" ("<<FCells[1][i]
555  // <<")"<<std::endl;
556  }
557  }
558  }
559  }
560  else
561  {
562  //marquage des cellules faisant partie des H2 libres
563  for (int j=0;j<FNbCycleDim2;++j)
564  {
565  for (int i=0;i<FMatrix[1]->getQ()->getnbli();++i)
566  {
567  if(FMatrix[1]->getQ()->getVal(i,j)!=0)
568  {
569  //marquerLibre la d cellule numero i
570  FMap->markOrbit(FCells[2][i],ORBIT_FACE,FMark);
571  }
572  }
573  }
574  }
575  }
576 
577  if ( FShowH2torsion && FMatrix[2]!=NULL )
578  {
579  //marquage des cellules faisant partie des H2 libres
580  for (int j=0;j<FNbGenTorsionDim2;++j )
581  {
582  for(int i=0;i<FMatrix[2]->getP()->getnbli();++i)
583  {
584  if(FMatrix[2]->getP()->getVal(i,j)!=0)
585  {
586  //marquerLibre la d cellule numero i
587  FMap->markOrbit(FCells[2][i],ORBIT_FACE,FMark);
588  // std::cout<<"Generateur libre "<<j<<" - Cellule libre: "<<i<<" ("<<FCells[1][i]
589  // <<")"<<std::endl;
590  }
591  }
592  }
593  }
594 }