Moka libraries
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
image-3d-extraction.icc
Go to the documentation of this file.
1 /*
2  * lib-extraction-images : Extraction de cartes à partir d'images 2D et 3D.
3  * Copyright (C) 2004, Moka Team, Université de Poitiers, Laboratoire SIC
4  * http://www.sic.sp2mi.univ-poitiers.fr/
5  * Copyright (C) 2009, Guillaume Damiand, CNRS, LIRIS,
6  * guillaume.damiand@liris.cnrs.fr, http://liris.cnrs.fr/
7  *
8  * This file is part of lib-extraction-images
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 
26 // La position des brins last, up est behind est donnée dans la thèse,
27 // (p 129) sauf qu'ici on est en G-carte donc le brin est celui du coté
28 // de la pointe des flèches
29 //
30 // Étant donné les trois brins last, up et behind, par rotation ils deviennent :
31 // last -> FMap->alpha01(up) -> FMap->alpha01(behind)
32 // up -> behind -> FMap->alpha10(last)
33 // behind -> FMap->alpha10(last) -> up
34 // (exprimés partout en last-up-behind d'origine)
35 
36 //******************************************************************************
37 #include "image-3d-magick.hh"
38 #include "image-3d-cea.hh"
39 #include "geometry.hh"
40 #include "pixel-step.hh"
41 
42 namespace GMap3d
43 {
44 
45 INLINE
46 CDartVertex* CExtractionImage::createFaceForExtract3d()
47 {
48  CDartVertex* prev;
49  CDartVertex* actu;
50  CDartVertex* first = (CDartVertex *) FMap->addMapDart();
51  FMap->linkAlpha0(first, FMap->addMapDart());
52 
53  prev = first;
54 
55  actu = (CDartVertex *) FMap->addMapDart();
56  FMap->linkAlpha0(actu, FMap->addMapDart());
57  FMap->linkAlpha1(actu,FMap->alpha0(prev));
58  prev = actu;
59 
60  actu = (CDartVertex *) FMap->addMapDart();
61  FMap->linkAlpha0(actu, FMap->addMapDart());
62  FMap->linkAlpha1(actu,FMap->alpha0(prev));
63  prev = actu;
64 
65  actu = (CDartVertex *) FMap->addMapDart();
66  FMap->linkAlpha0(actu, FMap->addMapDart());
67  FMap->linkAlpha1(actu,FMap->alpha0(prev));
68 
69  FMap->linkAlpha1(FMap->alpha0(actu),first);
70 
71  return first;
72 }
73 
74 //******************************************************************************
75 INLINE
76 CDartVertex* CExtractionImage::createFaceForExtract3d( const CVertex & AVertex )
77 {
78  CDartVertex* first = createFaceForExtract3d();
79 
80  FMap->addAttribute(first, ORBIT_VERTEX,
81  new CAttributeVertex(AVertex.getX()*STEP3D_X,
82  AVertex.getY()*STEP3D_Y,
83  AVertex.getZ()*STEP3D_Z));
84 
85  return first;
86 }
87 
88 //******************************************************************************
89 INLINE
90 CDartVertex* CExtractionImage::createCubeForExtract3d( const CVertex & AVertex )
91 {
92  CDartVertex* d1 = createFaceForExtract3d(); // face de gauche
93  CDartVertex* d2 = createFaceForExtract3d(); // face de derrière
94  CDartVertex* d3 = createFaceForExtract3d(); // face du bas
95  CDartVertex* d4 = createFaceForExtract3d(); // face du haut
96  CDartVertex* d5 = createFaceForExtract3d(); // face de devant
97  CDartVertex* last = createFaceForExtract3d(AVertex); // face de droite
98 
99  FMap->linkAlpha2(FMap->alpha101(d1),d2); FMap->linkAlpha2(FMap->alpha1010(d1), FMap->alpha0(d2));
100  FMap->linkAlpha2(FMap->alpha1(d1),FMap->alpha01(d3)); FMap->linkAlpha2(FMap->alpha10(d1),FMap->alpha010(d3));
101  FMap->linkAlpha2(FMap->alpha1(d2),FMap->alpha1010(d3)); FMap->linkAlpha2(FMap->alpha10(d2),FMap->alpha101(d3));
102  FMap->linkAlpha2(d4,FMap->alpha01(d1)); FMap->linkAlpha2(FMap->alpha0(d4),FMap->alpha010(d1));
103  FMap->linkAlpha2(FMap->alpha01(d4),FMap->alpha01(d2)); FMap->linkAlpha2(FMap->alpha010(d4),FMap->alpha010(d2));
104  FMap->linkAlpha2(FMap->alpha1(d5),d3); FMap->linkAlpha2(FMap->alpha10(d5),FMap->alpha0(d3));
105  FMap->linkAlpha2(FMap->alpha01(d5),FMap->alpha10(d4)); FMap->linkAlpha2(FMap->alpha010(d5),FMap->alpha1(d4));
106  FMap->linkAlpha2(FMap->alpha101(d5),d1); FMap->linkAlpha2(FMap->alpha1010(d5),FMap->alpha0(d1));
107 
108  FMap->linkAlpha2(last, d5); FMap->linkAlpha2(FMap->alpha0(last), FMap->alpha0(d5));
109  FMap->linkAlpha2(FMap->alpha01(last),FMap->alpha101(d4)); FMap->linkAlpha2(FMap->alpha010(last), FMap->alpha1010(d4));
110  FMap->linkAlpha2(FMap->alpha101(last), FMap->alpha101(d2)); FMap->linkAlpha2(FMap->alpha1010(last), FMap->alpha1010(d2));
111  FMap->linkAlpha2(FMap->alpha1(last), FMap->alpha1(d3)); FMap->linkAlpha2(FMap->alpha10(last), FMap->alpha10(d3));
112 
113  return last; // brin de devant en bas de la face de droite
114 }
115 
116 //******************************************************************************
117 INLINE
118 CDart* CExtractionImage::createBorder3d(unsigned int columns,unsigned int rows)
119 {
120  CDart* prev;
121  CDart* prev2;
122  CDart* actu;
123 
124  // Construction de la premiere face, FMap->alpha gauche du premier voxel
125  CDart* first = createFaceForExtract3d(CVertex(0,1,1));
126 
127  prev = FMap->alpha101(first);
128 
129  // Construction des faces derrieres les voxels de la premiere ligne
130  for (unsigned int i=0;i<=columns;++i)
131  {
132  actu = createFaceForExtract3d(CVertex(i,0,1));
133  FMap->linkAlpha2(prev,actu); FMap->linkAlpha2(FMap->alpha0(prev),FMap->alpha0(actu));
134 
135  prev = FMap->alpha101(actu);
136  }
137 
138  // Couture pour le repliement
139  FMap->linkAlpha2(prev, first); FMap->linkAlpha2(FMap->alpha0(prev), FMap->alpha0(first));
140 
141  // Construction des faces au dessus des voxels de la premiere plaque
142  prev = FMap->alpha010(first);
143  prev2 = FMap->alpha101201(first);
144  for (unsigned int j=0;j<=rows;++j)
145  {
146  for (unsigned int i=0;i<=columns;++i)
147  {
148  actu = createFaceForExtract3d(CVertex(i,j,0));
149 
150  FMap->linkAlpha2(actu,prev); FMap->linkAlpha2(FMap->alpha0(actu),FMap->alpha0(prev));
151 
152  FMap->linkAlpha2(FMap->alpha1(actu),prev2); FMap->linkAlpha2(FMap->alpha10(actu),FMap->alpha0(prev2));
153 
154  prev = FMap->alpha101(actu);
155  prev2 = FMap->alpha0121(prev2);
156  }
157 
158  if (j==0) prev2 = FMap->alpha21(prev2); // correctif pour la premiere ligne
159  }
160 
161  // Couture pour le repliement
162  FMap->linkAlpha2(prev,FMap->alpha10(first)); FMap->linkAlpha2(FMap->alpha0(prev),FMap->alpha1(first));
163 
164  // Coutures des bord de devant de la plaque sup avec le bord inf de la ligne derriere
165  prev = FMap->alpha10121(first);
166  for (unsigned int i=0;i<=columns;++i)
167  {
168  FMap->linkAlpha2(prev,prev2); FMap->linkAlpha2(FMap->alpha0(prev),FMap->alpha0(prev2));
169 
170  prev = FMap->alpha0121(prev);
171  prev2 = FMap->alpha0121(prev2);
172  }
173 
174  return FMap->alpha101(first); // Le brin derriere et en bas de la face FMap->alpha gauche du premier voxel
175 }
176 
177 //******************************************************************************
178 INLINE
179 CDart* CExtractionImage::computeUpFromLast3d(CDart* ADart)
180 {
181  CDart* res = FMap->alpha012(ADart);
182 
183  while(!FMap->isFree3(res))
184  res = FMap->alpha32(res);
185 
186  return FMap->alpha010(res);
187 }
188 
189 //******************************************************************************
190 INLINE
191 CDart* CExtractionImage::computeBehindFromLast3d(CDart* ADart)
192 {
193  CDart* res = FMap->alpha2(ADart);
194 
195  while(!FMap->isFree3(res))
196  res = FMap->alpha32(res);
197 
198  return FMap->alpha101(res);
199 }
200 //******************************************************************************
201 INLINE
202 void CExtractionImage::mergeImage3dVolumes(CDart* ADart)
203 {
204  CDart* tmp;
205  CDart* tmp2;
206  int i;
207 
208  for (i=0;i<4;++i)
209  {
210  if (FMap->alpha3(ADart)!=FMap->alpha2(ADart))
211  {
212  tmp = FMap->alpha32(ADart); tmp2=FMap->alpha2(ADart);
213  FMap->unsew2(ADart); FMap->unsew2(FMap->alpha3(ADart));
214  FMap->sew2(tmp, tmp2);
215  }
216 
217  ADart = FMap->alpha01(ADart);
218  }
219 
220  for (i=0;i<4;++i)
221  {
222  tmp = ADart;
223  ADart = FMap->alpha01(ADart);
224 
225  FMap->delMapDart(FMap->alpha30(tmp));
226  FMap->delMapDart(FMap->alpha3(tmp));
227  FMap->delMapDart(FMap->alpha0(tmp));
228  FMap->delMapDart(tmp);
229  }
230 }
231 
232 //******************************************************************************
233 INLINE
234 void CExtractionImage::mergeImage3dFaces(CDart* ADart)
235 {
236  CDart *tmp, *tmp2;
237 
238  if (FMap->alpha2(ADart)!=FMap->alpha1(ADart))
239  {
240  tmp=FMap->alpha1(ADart); tmp2=FMap->alpha21(ADart);
241 
242  FMap->unsew1(ADart); FMap->unsew1(FMap->alpha2(ADart));
243 
244  FMap->sew1(tmp,tmp2);
245  }
246 
247  if (FMap->alpha02(ADart)!=FMap->alpha01(ADart))
248  {
249  tmp=FMap->alpha01(ADart); tmp2=FMap->alpha021(ADart);
250 
251  FMap->unsew1(FMap->alpha0(ADart));
252  FMap->unsew1(FMap->alpha02(ADart));
253 
254  FMap->sew1(tmp,tmp2);
255  }
256 
257  FMap->delMapDart(FMap->alpha302(ADart));
258  FMap->delMapDart(FMap->alpha32(ADart));
259  FMap->delMapDart(FMap->alpha30(ADart));
260  FMap->delMapDart(FMap->alpha3(ADart));
261  FMap->delMapDart(FMap->alpha02(ADart));
262  FMap->delMapDart(FMap->alpha2(ADart));
263  FMap->delMapDart(FMap->alpha0(ADart));
264  FMap->delMapDart(ADart);
265 }
266 
267 //******************************************************************************
268 INLINE
269 void CExtractionImage::mergeImage3dEdges(CDart* ADart)
270  // Attention : cette méthode plante si le sommet est de degré un.
271 {
272  CDart *tmp, *tmp2;
273 
274  tmp = FMap->alpha0(ADart);
275  tmp2 = FMap->alpha10(ADart);
276 
277  FMap->topoUnsew0(ADart);
278  FMap->topoUnsew0(FMap->alpha1(ADart));
279 
280  FMap->topoSew0(tmp,tmp2);
281 
282  CStaticCoverage123 it(FMap,ADart);
283  while (it.cont())
284  FMap->delMapDart(it++);
285 }
286 
287 //******************************************************************************
288 INLINE
289 bool CExtractionImage::isDegre2Vertex3d(CDart* ADart)
290 {
291  CDynamicCoverage23 it(FMap,ADart);
292  while (it.cont())
293  {
294  if (FMap->alpha2(*it)==FMap->alpha1(*it)) return false;
295  if (FMap->alpha12(ADart)!=FMap->alpha21(ADart)) return false;
296  ++it;
297  }
298  return true;
299 }
300 
301 //******************************************************************************
302 INLINE
303 CDart * CExtractionImage::isRealDegre2Vertex3d(CDart* ADart, int AMark)
304 {
305  CDart * other = NULL;
306  int mark = FMap->getNewMark();
307  int nb=0;
308 
309  CDynamicCoverage123 it(FMap,ADart);
310  while (nb<3 && it.cont())
311  {
312  if (!FMap->isMarked(*it,mark))
313  {
314  if (!FMap->isMarked(*it,AMark))
315  {
316  other = *it;
317  ++nb;
318  }
319 
320  CDynamicCoverage023 it2(FMap,*it);
321  while(it2.cont())
322  FMap->setMark(it2++,mark);
323  }
324  ++it;
325  }
326 
327  FMap->unmarkAll(mark);
328  FMap->freeMark(mark);
329 
330  if (nb!=2) other=NULL;
331 
332  return other;
333 }
334 
335 //******************************************************************************
336 INLINE
337 bool CExtractionImage::areFacesColinear(CDart* ADart)
338 {
339  CVertex vector1 = FMap->faceNormalVector(ADart);
340  CVertex vector2 = FMap->faceNormalVector(FMap->alpha2(ADart));
341 
342  return CGeometry::areColinear(vector1, vector2);
343 }
344 
345 //******************************************************************************
346 INLINE
347 bool CExtractionImage::isDegre2Edge3d(CDart* ADart)
348 {
349  return (FMap->alpha23(ADart)==FMap->alpha32(ADart));
350 }
351 
352 //******************************************************************************
353 INLINE
354 void CExtractionImage::sew3Cube(CDart* ADart, CDart* TheLastDart)
355 {
356  FMap->topoSew3(FMap->alpha012(ADart),FMap->alpha1(computeUpFromLast3d(TheLastDart)));
357 
358  FMap->topoSew3(FMap->alpha1012(ADart),computeBehindFromLast3d(TheLastDart));
359 
360  FMap->topoSew3(FMap->alpha21012(ADart),FMap->alpha101(TheLastDart));
361 }
362 
363 //******************************************************************************
364 INLINE
365 bool CExtractionImage::faceMergingWillDisconnect(CDart* ADart)
366 {
367  if (FMap->alpha1(ADart)==FMap->alpha2(ADart) ||
368  FMap->alpha01(ADart)==FMap->alpha02(ADart)) return false;
369 
370  CDynamicCoverage01 it(FMap,ADart);
371  while (it.cont())
372  {
373  if (it++==FMap->alpha2(ADart)) return true;
374  }
375 
376  return false;
377 }
378 //******************************************************************************
379 INLINE
380 bool CExtractionImage::extract3dImage(CImage3d * AImage, int ALevel,
381  bool ADestroyBorder,
382  bool AKeepFictiveEdges,
383  bool AShiftFictiveEdges)
384 {
385  assert(AImage != NULL);
386 
387  if (! AImage->isOk())
388  return false;
389 
390 #ifdef DEBUG_EXTRACT_IMAGE
391  cout<<"Avant de créer le bord supérieur : "<<endl;
392 #endif
393 
394  CDart * last = createBorder3d(AImage->getWidth(),AImage->getHeight());
395 
396 #ifdef DEBUG_EXTRACT_IMAGE
397  cout<<"Ok on FMap->alpha créé le bord supérieur."<<endl;
398 #endif
399 
400  if (ALevel>=0)
401  {
402  CDart * tmp;
403  unsigned int k=0;
404 
405  // Brins utilisés pour la fusion de volumes
406  CDart * leftFace;
407  CDart * behindFace;
408  CDart * upFace;
409 
410  // Brins utilisés pour la fusion de faces
411  CDart * leftFace2;
412  CDart * behindFace2;
413  CDart * upFace2;
414 
415  // Brin utilisé pour la fusion d'arête
416  CDart * vertex;
417  CDart * vertex1;
418  CDart * vertex2;
419  CDart * vertex3;
420 
421  // Marque pour les arêtes fictives
422  int fictiveMark = FMap->getNewMark();
423 
424  // Les coordonnées du pointel courant
425  // (multipliés ensuite par STEP3D_X, STEP3D_Y et STEP3D_Z)
426  CVertex CurrentPointel(1,1,1);
427 
428  do // On boucle tant qu'on peut lire un prochain plan
429  {
430  for (unsigned int j=0;j<=AImage->getHeight();++j)
431  {
432  for (unsigned int i=0;i<=AImage->getWidth();++i)
433  {
434 #ifdef DEBUG_EXTRACT_IMAGE
435  cout<<"On traite le voxel ("<<i<<", "<<j<<", "<<k<<")"<<endl;
436 #endif
437 
438  tmp = createCubeForExtract3d(CurrentPointel);
439  sew3Cube(tmp,last);
440 
441  if (ALevel>=1)
442  {
443  leftFace = FMap->alpha0(last);
444  upFace = FMap->alpha3012(last);
445  behindFace = FMap->alpha3201(last);
446 
447  vertex1 = FMap->alpha23(leftFace);
448  vertex2 = FMap->alpha323(upFace);
449  vertex3 = FMap->alpha323(behindFace);
450 
451  // Les trois fusions de volumes possibles
452  if (AImage->sameVoxelActuLeft(i,j))
453  {
454  if (vertex2==FMap->alpha13(leftFace)) vertex2=NULL;
455 
456  mergeImage3dVolumes(leftFace);
457  leftFace = NULL;
458  }
459 
460  if (AImage->sameVoxelActuBehind(i,j))
461  {
462  if (vertex1==FMap->alpha1(behindFace)) vertex1=NULL;
463 
464  mergeImage3dVolumes(behindFace);
465  behindFace = NULL;
466  }
467 
468  if (AImage->sameVoxelActuUp(i,j))
469  {
470  if (vertex3==FMap->alpha1(upFace)) vertex3=NULL;
471 
472  mergeImage3dVolumes(upFace);
473  upFace = NULL;
474  }
475 
476  if (ALevel>=2)
477  {
478  // Les six fusions de faces coplanaires possibles
479  if (leftFace!=NULL)
480  {
481  leftFace2 = FMap->alpha1(leftFace);
482 
483  if ( isDegre2Edge3d(leftFace) &&
484  areFacesColinear(leftFace) )
485  {
486  mergeImage3dFaces(leftFace);
487  leftFace = NULL;
488  vertex1=NULL;
489  }
490 
491  if ( isDegre2Edge3d(leftFace2) &&
492  areFacesColinear(leftFace2) )
493  {
494  if ( AKeepFictiveEdges &&
495  leftFace == NULL &&
496  FMap->isMarked(FMap->alpha1(leftFace2),fictiveMark) )
497  {
498  FMap->markOrbit(leftFace2,ORBIT_023,fictiveMark);
499  mergeImage3dEdges(leftFace2);
500  leftFace2 = NULL;
501  vertex2 = NULL;
502  }
503  else
504  if ( leftFace != NULL ||
505  !AKeepFictiveEdges ||
506  !faceMergingWillDisconnect(leftFace2) )
507  {
508  mergeImage3dFaces(leftFace2);
509  leftFace2 = NULL;
510  vertex2 = NULL;
511  }
512  else
513  {
514  if (AKeepFictiveEdges)
515  FMap->markOrbit(leftFace2,ORBIT_023,fictiveMark);
516  }
517  }
518  }
519 
520  if (behindFace!=NULL)
521  {
522  behindFace2 = FMap->alpha1(behindFace);
523 
524  if ( isDegre2Edge3d(behindFace) &&
525  areFacesColinear(behindFace) )
526  {
527  mergeImage3dFaces(behindFace);
528  behindFace = NULL;
529  vertex3 = NULL;
530  }
531 
532  if ( isDegre2Edge3d(behindFace2) &&
533  areFacesColinear(behindFace2) )
534  {
535  if ( behindFace == NULL &&
536  AKeepFictiveEdges &&
537  FMap->isMarked(FMap->alpha1(behindFace2),fictiveMark) )
538  {
539  FMap->markOrbit(behindFace2,ORBIT_023,fictiveMark);
540  mergeImage3dEdges(behindFace2);
541  behindFace2 = NULL;
542  vertex1 = NULL;
543  }
544  else
545  if ( behindFace != NULL ||
546  !AKeepFictiveEdges ||
547  !faceMergingWillDisconnect(behindFace2) )
548  {
549  mergeImage3dFaces(behindFace2);
550  behindFace2 = NULL;
551  vertex1 = NULL;
552  }
553  else
554  {
555  if (AKeepFictiveEdges)
556  FMap->markOrbit(behindFace2,ORBIT_023,fictiveMark);
557  }
558  }
559  }
560 
561  if (upFace!=NULL)
562  {
563  upFace2 = FMap->alpha1(upFace);
564 
565  if ( isDegre2Edge3d(upFace) &&
566  areFacesColinear(upFace) )
567  {
568  mergeImage3dFaces(upFace);
569  upFace = NULL;
570  vertex2 = NULL;
571  }
572 
573  if ( isDegre2Edge3d(upFace2) &&
574  areFacesColinear(upFace2) )
575  {
576  if ( upFace == NULL &&
577  AKeepFictiveEdges &&
578  FMap->isMarked(FMap->alpha1(upFace2),fictiveMark) )
579  {
580  FMap->markOrbit(upFace2,ORBIT_023,fictiveMark);
581  mergeImage3dEdges(upFace2);
582  upFace2 = NULL;
583  vertex3 = NULL;
584  }
585  else
586  if ( upFace != NULL ||
587  !AKeepFictiveEdges ||
588  !faceMergingWillDisconnect(upFace2) )
589  {
590  mergeImage3dFaces(upFace2);
591  upFace2 = NULL;
592  vertex3 = NULL;
593  }
594  else
595  {
596  if (AKeepFictiveEdges)
597  FMap->markOrbit(upFace2,ORBIT_023,fictiveMark);
598  }
599  }
600  }
601 
602  if (ALevel>=3)
603  {
604  vertex = NULL;
605 
606  if (vertex1!=NULL && vertex2==NULL && vertex3==NULL)
607  vertex = vertex1;
608  else if (vertex2!=NULL && vertex1==NULL && vertex3==NULL)
609  vertex = vertex2;
610  else if (vertex3!=NULL && vertex1==NULL && vertex2==NULL)
611  vertex = vertex3;
612 
613  if ( vertex!=NULL )
614  {
615  if (AKeepFictiveEdges && AShiftFictiveEdges)
616  {
617  vertex2=isRealDegre2Vertex3d(vertex,fictiveMark);
618  if ( vertex2!=NULL &&
619  areEdgesAlign(vertex,vertex2) )
620  {
621  FMap->shiftAllFictiveEdges(vertex,fictiveMark);
622  mergeImage3dEdges(vertex);
623  }
624  }
625  else
626  {
627  if ( isDegre2Vertex3d(vertex) &&
628  areEdgesAlign(vertex,FMap->alpha1(vertex)) )
629  {
630  mergeImage3dEdges(vertex);
631  }
632  }
633  }
634  }
635  }
636  }
637 
638  last = FMap->alpha101(tmp);
639 
640  // Calcule la prochaine position du pointel courant
641  // Équivalent à faire :
642  // CurrentPointel.setX((i+1)%(columns+1));
643  // CurrentPointel.setY((j+1+((i+1)/(columns+1)))%(rows+1));
644  // CurrentPointel.setZ(k+1+((j+1+((i+1)/(columns+1)))/(rows+1)));
645  // Mais ici c'est optimisé car moins de formule math :)
646  if ( i == AImage->getWidth()-1 )
647  {
648  CurrentPointel.setX(0);
649 
650  if ( j == AImage->getHeight()-1 )
651  {
652  CurrentPointel.setY(0);
653  CurrentPointel.setZ(CurrentPointel.getZ()+1);
654  }
655  else
656  CurrentPointel.setY(CurrentPointel.getY()+1);
657  }
658  else
659  CurrentPointel.setX(CurrentPointel.getX()+1);
660  }
661  }
662 
663  std::cout << "Ok pour l'extraction de la plaque num "
664  << AImage->getPlaneActu() << std::endl;
665  ++k;
666  }
667  while(AImage->readNextPlane());
668 
669  if (ADestroyBorder && ALevel > 0)
670  {
671  CStaticCoverage012 it(FMap,last);
672  while(it.cont())
673  FMap->delMapDart(it++);
674  }
675 
676  FMap->unmarkAll(fictiveMark);
677  FMap->freeMark(fictiveMark);
678  }
679 
680  return true;
681 }
682 //******************************************************************************
683 INLINE
684 bool CExtractionImage::extract3dImage( const std::string & AFilename,
685  int FirstPlane,
686  int NbPlaneToRead,
687  int Level,
688  int Lg,
689  bool DestroyBorder,
690  bool KeepFictiveEdges,
691  bool ShiftFictiveEdges )
692 {
693 #ifdef MODULE_EXTRACTION_IMAGE_WITHOUT_MAGICK
694  return false;
695 #else
696  CImage3dMagick image(AFilename,FirstPlane,NbPlaneToRead,Lg);
697  return extract3dImage(& image, Level,
698  DestroyBorder, KeepFictiveEdges, ShiftFictiveEdges);
699 #endif // MODULE_EXTRACTION_IMAGE_WITHOUT_MAGICK
700 }
701 //******************************************************************************
702 INLINE
703 bool CExtractionImage::extract3dImageCEA( const std::string & AFilename,
704  int FirstPlane,
705  int NbPlaneToRead,
706  int Level,
707  bool DestroyBorder,
708  bool KeepFictiveEdges,
709  bool ShiftFictiveEdges )
710 {
711  CImage3dCEA image(AFilename,FirstPlane,NbPlaneToRead);
712 
713  return extract3dImage(& image, Level,
714  DestroyBorder, KeepFictiveEdges, ShiftFictiveEdges);
715 }
716 //******************************************************************************
717 INLINE
718 bool CExtractionImage::extractOneRegionVoxels(CImage3d * AImage,
719  int ARed, int AGreen,
720  int ABlue, int AAlpha )
721 {
722  assert(AImage != NULL);
723 
724  if (! AImage->isOk())
725  return false;
726 
727 #ifdef DEBUG_EXTRACT_IMAGE
728  std::cout<<"Début extractOneRegionVoxels:"<<std::flush;
729 #endif
730 
731  int toKeep = FMap->getNewMark();
732 
733  // Pour garder les brins déjà existant.
734  FMap->negateMaskMark(toKeep);
735 
736 #ifdef DEBUG_EXTRACT_IMAGE
737  cout<<"\nAvant de créer le bord supérieur : "<<endl;
738 #endif
739 
740  CDart * last = createBorder3d(AImage->getWidth(),AImage->getHeight());
741 
742 #ifdef DEBUG_EXTRACT_IMAGE
743  std::cout<<"Ok on a créé le bord supérieur."<<std::endl;
744 #endif
745 
746 
747 #ifdef DEBUG_EXTRACT_IMAGE
748  cout<<"Ok on FMap->alpha créé le bord supérieur."<<endl;
749 #endif
750 
751  CDart * tmp;
752  unsigned int k=0;
753  unsigned int nbpoints = 0;
754 
755  // Brins utilisés pour la fusion de volumes
756  CDart * leftFace;
757  CDart * behindFace;
758  CDart * upFace;
759 
760  // Les coordonnées du pointel courant
761  // (multipliés ensuite par STEP3D_X, STEP3D_Y et STEP3D_Z)
762  CVertex CurrentPointel(1,1,1);
763 
764  do // On boucle tant qu'on peut lire un prochain plan
765  {
766  for (unsigned int j=0;j<=AImage->getHeight();++j)
767  {
768  for (unsigned int i=0;i<=AImage->getWidth();++i)
769  {
770 #ifdef DEBUG_EXTRACT_IMAGE
771  cout<<"On traite le voxel ("<<i<<", "<<j<<", "<<k<<")"<<endl;
772 #endif
773 
774  tmp = createCubeForExtract3d(CurrentPointel);
775  sew3Cube(tmp,last);
776 
777  // Si le voxel est dans la région AnId, on le marque.
778  if ( i<AImage->getWidth() && j<AImage->getHeight() &&
779  !AImage->lastPlane() &&
780  AImage->sameVoxel(i,j,true,ARed, AGreen, ABlue, AAlpha) )
781  {
782  FMap->markOrbit(tmp, ORBIT_VOLUME, toKeep);
783  ++nbpoints;
784  }
785  else
786  {
787  // Sinon on fait les fusions avec les autres voxels.
788  leftFace = FMap->alpha0(last);
789  upFace = FMap->alpha3012(last);
790  behindFace = FMap->alpha3201(last);
791 
792  // Les trois fusions de volumes possibles
793  if (AImage->sameVoxelActuLeft(i,j))
794  {
795  mergeImage3dVolumes(leftFace);
796  leftFace = NULL;
797  }
798 
799  if (AImage->sameVoxelActuBehind(i,j))
800  {
801  mergeImage3dVolumes(behindFace);
802  behindFace = NULL;
803  }
804 
805  if (AImage->sameVoxelActuUp(i,j))
806  {
807  mergeImage3dVolumes(upFace);
808  upFace = NULL;
809  }
810 
811  }
812 
813  last = FMap->alpha101(tmp);
814  // Calcule la prochaine position du pointel courant
815  // Équivalent à faire :
816  // CurrentPointel.setX((i+1)%(columns+1));
817  // CurrentPointel.setY((j+1+((i+1)/(columns+1)))%(rows+1));
818  // CurrentPointel.setZ(k+1+((j+1+((i+1)/(columns+1)))/(rows+1)));
819  // Mais ici c'est optimisé car moins de formule math :)
820  if ( i == AImage->getWidth()-1 )
821  {
822  CurrentPointel.setX(0);
823 
824  if ( j == AImage->getHeight()-1 )
825  {
826  CurrentPointel.setY(0);
827  CurrentPointel.setZ(CurrentPointel.getZ()+1);
828  }
829  else
830  CurrentPointel.setY(CurrentPointel.getY()+1);
831  }
832  else
833  CurrentPointel.setX(CurrentPointel.getX()+1);
834  }
835  }
836 
837 #ifdef DEBUG_EXTRACT_IMAGE
838  std::cout << "Ok pour l'extraction de la plaque num "
839  << AImage->getPlaneActu() << std::endl;
840 #endif
841 
842  ++k;
843  }
844  while(AImage->readNextPlane());
845 
846  // On détruit le bord
847  CStaticCoverage012 it(FMap,last);
848  while(it.cont())
849  FMap->delMapDart(it++);
850 
851 #ifdef DEBUG_EXTRACT_IMAGE
852  std::cout << "Ok pour la destruction du bord." << std::endl;
853 #endif
854 
855  // Maintenant on détruit tout les brins qui ne sont pas dans la région
856  // à garder, en 3-décousant préalablement les volumes.
857  CDynamicCoverageAll itAll(FMap);
858  while( itAll.cont() )
859  {
860  tmp = itAll++;
861  if ( !FMap->isMarked(tmp, toKeep) )
862  {
863  if ( !FMap->isFree3(tmp) )
864  {
865  if ( FMap->isMarked(FMap->alpha3(tmp), toKeep) )
866  {
867  FMap->unsew3(tmp);
868  }
869  }
870  }
871  }
872 
873 #ifdef DEBUG_EXTRACT_IMAGE
874  std::cout << "Ok pour la décousure par alpha3." << std::endl;
875 #endif
876 
877  itAll.reinit();
878  while( itAll.cont() )
879  {
880  tmp = itAll++;
881 
882  if ( !FMap->isMarked(tmp, toKeep) )
883  FMap->delMapDart(tmp);
884  }
885 
886  FMap->negateMaskMark(toKeep);
887 
888  assert( FMap->isWholeMapUnmarked(toKeep) );
889 
890  FMap->freeMark(toKeep);
891 
892 
893 #ifdef DEBUG_EXTRACT_IMAGE
894  std::cout << "Ok pour l'extraction d'une région : terminée. ";
895  std::cout<<nbpoints<<" voxels."<<std::endl;
896 #endif
897 
898  return true;
899 }
900 //******************************************************************************
901 INLINE
902 bool CExtractionImage::extractOneRegionVoxels( const std::string & AFilename,
903  int FirstPlane,
904  int NbPlaneToRead,
905  int Lg,
906  int ARed, int AGreen,
907  int ABlue, int AAlpha )
908 {
909 #ifdef MODULE_EXTRACTION_IMAGE_WITHOUT_MAGICK
910  return false;
911 #else
912  CImage3dMagick image(AFilename,FirstPlane,NbPlaneToRead,Lg);
913  return extractOneRegionVoxels(& image, ARed, AGreen, ABlue, AAlpha);
914 #endif // MODULE_EXTRACTION_IMAGE_WITHOUT_MAGICK
915 }
916 //******************************************************************************
917 } // namespace GMap3d
918 //******************************************************************************