src/VTKViewer/VTKViewer_GeometryFilter.cxx

   1 //  SALOME OBJECT : kernel of SALOME component
   2 //
   3 //  Copyright (C) 2003  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 //  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  21 //
  22 //
  23 //
  24 //  File   : VTKViewer_GeometryFilter.cxx
  25 //  Author : Michael ZORIN
  26 //  Module : SALOME
  27 //  $Header$
  28
  29 #include "VTKViewer_GeometryFilter.h"
  30
  31 #include <vtkSmartPointer.h>
  32 #include <vtkCellArray.h>
  33 #include <vtkCellData.h>
  34 #include <vtkGenericCell.h>
  35 #include <vtkHexahedron.h>
  36 #include <vtkMergePoints.h>
  37 #include <vtkObjectFactory.h>
  38 #include <vtkPointData.h>
  39 #include <vtkPolyData.h>
  40 #include <vtkPyramid.h>
  41 #include <vtkStructuredGrid.h>
  42 #include <vtkTetra.h>
  43 #include <vtkUnsignedCharArray.h>
  44 #include <vtkUnstructuredGrid.h>
  45 #include <vtkVoxel.h>
  46 #include <vtkWedge.h>
  47
  48 #include <vtkMath.h>
  49 #include <vtkPlane.h>
  50 #include <vtkDelaunay3D.h>
  51 #include <vtkGeometryFilter.h>
  52
  53 #include <algorithm>
  54 #include <iterator>
  55 #include <vector>
  56 #include <map>
  57 #include <set>
  58
  59 #if defined __GNUC__
  60   #if __GNUC__ == 2
  61     #define __GNUC_2__
  62   #endif
  63 #endif
  64
  65 #define USE_ROBUST_TRIANGULATION
  66
  67 //----------------------------------------------------------------------------
  68 vtkCxxRevisionMacro(VTKViewer_GeometryFilter, "$Revision$");
  69 vtkStandardNewMacro(VTKViewer_GeometryFilter);
  70
  71 VTKViewer_GeometryFilter
  72 ::VTKViewer_GeometryFilter():
  73   myShowInside(0),
  74   myStoreMapping(0),
  75   myIsWireframeMode(0)
  76 {}
  77
  78
  79 VTKViewer_GeometryFilter
  80 ::~VTKViewer_GeometryFilter()
  81 {}
  82
  83
  84 //----------------------------------------------------------------------------
  85 void
  86 VTKViewer_GeometryFilter
  87 ::Execute()
  88 {
  89   vtkDataSet *input= this->GetInput();
  90   vtkIdType numCells=input->GetNumberOfCells();
  91
  92   if (numCells == 0)
  93     {
  94       return;
  95     }
  96
  97   if (input->GetDataObjectType() == VTK_UNSTRUCTURED_GRID){
  98     this->UnstructuredGridExecute();
  99     return;
 100   }else
 101     vtkGeometryFilter::Execute();
 102 }
 103
 104
 105 //----------------------------------------------------------------------------
 106 void
 107 VTKViewer_GeometryFilter
 108 ::UnstructuredGridExecute()
 109 {
 110   vtkUnstructuredGrid *input= (vtkUnstructuredGrid *)this->GetInput();
 111   vtkCellArray *Connectivity = input->GetCells();
 112   // Check input
 113   if ( Connectivity == NULL )
 114     {
 115     vtkDebugMacro(<<"Nothing to extract");
 116     return;
 117     }
 118
 119   vtkIdType cellId;
 120   int i;
 121   int allVisible;
 122   vtkIdType npts = 0;
 123   vtkIdType *pts = 0;
 124   vtkPoints *p = input->GetPoints();
 125   vtkIdType numCells=input->GetNumberOfCells();
 126   vtkPointData *pd = input->GetPointData();
 127   vtkCellData *cd = input->GetCellData();
 128   vtkPolyData *output = this->GetOutput();
 129   vtkPointData *outputPD = output->GetPointData();
 130
 131 #ifdef USE_ROBUST_TRIANGULATION
 132   vtkUnstructuredGrid* anUnstructuredGrid = vtkUnstructuredGrid::New();
 133   vtkPoints* aDelaunayPoints = vtkPoints::New();
 134
 135   vtkDelaunay3D* aDelaunay3D = vtkDelaunay3D::New();
 136   aDelaunay3D->SetInput(anUnstructuredGrid);
 137
 138   vtkGeometryFilter* aGeometryFilter = vtkGeometryFilter::New();
 139   aGeometryFilter->SetInput(aDelaunay3D->GetOutput());
 140 #endif
 141
 142   vtkCellData *outputCD = output->GetCellData();
 143   vtkGenericCell *cell = vtkGenericCell::New();
 144
 145
 146   vtkIdList *cellIds = vtkIdList::New();
 147   vtkIdList *faceIds = vtkIdList::New();
 148
 149   char *cellVis;
 150   vtkIdType newCellId;
 151   int faceId, *faceVerts, numFacePts;
 152   float *x;
 153   int PixelConvert[4], aNewPts[VTK_CELL_SIZE];
 154   // ghost cell stuff
 155   unsigned char  updateLevel = (unsigned char)(output->GetUpdateGhostLevel());
 156   unsigned char  *cellGhostLevels = 0;
 157
 158   PixelConvert[0] = 0;
 159   PixelConvert[1] = 1;
 160   PixelConvert[2] = 3;
 161   PixelConvert[3] = 2;
 162
 163   vtkDebugMacro(<<"Executing geometry filter for unstructured grid input");
 164
 165   vtkDataArray* temp = 0;
 166   if (cd)
 167     {
 168     temp = cd->GetArray("vtkGhostLevels");
 169     }
 170   if ( (!temp) || (temp->GetDataType() != VTK_UNSIGNED_CHAR)
 171     || (temp->GetNumberOfComponents() != 1))
 172     {
 173     vtkDebugMacro("No appropriate ghost levels field available.");
 174     }
 175   else
 176     {
 177     cellGhostLevels = ((vtkUnsignedCharArray*)temp)->GetPointer(0);
 178     }
 179
 180   // Determine nature of what we have to do
 181   if ( (!this->CellClipping) && (!this->PointClipping) &&
 182        (!this->ExtentClipping) )
 183     {
 184     allVisible = 1;
 185     cellVis = NULL;
 186     }
 187   else
 188     {
 189     allVisible = 0;
 190     cellVis = new char[numCells];
 191     }
 192
 193   // Just pass points through, never merge
 194   output->SetPoints(input->GetPoints());
 195   outputPD->PassData(pd);
 196
 197   outputCD->CopyAllocate(cd,numCells,numCells/2);
 198
 199   output->Allocate(numCells/4+1,numCells);
 200
 201   // Loop over the cells determining what's visible
 202   if (!allVisible)
 203     {
 204     for (cellId=0, Connectivity->InitTraversal();
 205          Connectivity->GetNextCell(npts,pts);
 206          cellId++)
 207       {
 208       cellVis[cellId] = 1;
 209       if ( this->CellClipping && cellId < this->CellMinimum ||
 210            cellId > this->CellMaximum )
 211         {
 212         cellVis[cellId] = 0;
 213         }
 214       else
 215         {
 216         for (i=0; i < npts; i++)
 217           {
 218           x = p->GetPoint(pts[i]);
 219           if ( (this->PointClipping && (pts[i] < this->PointMinimum ||
 220                                         pts[i] > this->PointMaximum) ) ||
 221                (this->ExtentClipping &&
 222                 (x[0] < this->Extent[0] || x[0] > this->Extent[1] ||
 223                  x[1] < this->Extent[2] || x[1] > this->Extent[3] ||
 224                  x[2] < this->Extent[4] || x[2] > this->Extent[5] )) )
 225             {
 226             cellVis[cellId] = 0;
 227             break;
 228             }//point/extent clipping
 229           }//for each point
 230         }//if point clipping needs checking
 231       }//for all cells
 232     }//if not all visible
 233
 234   // Loop over all cells now that visibility is known
 235   // (Have to compute visibility first for 3D cell boundarys)
 236   int progressInterval = numCells/20 + 1;
 237   if(myStoreMapping){
 238     myVTK2ObjIds.clear();
 239     myVTK2ObjIds.reserve(numCells);
 240   }
 241   for (cellId=0, Connectivity->InitTraversal();
 242        Connectivity->GetNextCell(npts,pts);
 243        cellId++)
 244     {
 245     //Progress and abort method support
 246     if ( !(cellId % progressInterval) )
 247       {
 248       vtkDebugMacro(<<"Process cell #" << cellId);
 249       this->UpdateProgress ((float)cellId/numCells);
 250       }
 251
 252     // Handle ghost cells here.  Another option was used cellVis array.
 253     if (cellGhostLevels && cellGhostLevels[cellId] > updateLevel)
 254       { // Do not create surfaces in outer ghost cells.
 255       continue;
 256       }
 257
 258     if (allVisible || cellVis[cellId])  //now if visible extract geometry
 259       {
 260       //special code for nonlinear cells - rarely occurs, so right now it
 261       //is slow.
 262       vtkIdType aCellType = input->GetCellType(cellId);
 263       switch (aCellType)
 264         {
 265         case VTK_EMPTY_CELL:
 266           break;
 267
 268         case VTK_VERTEX:
 269         case VTK_POLY_VERTEX:
 270           newCellId = output->InsertNextCell(aCellType,npts,pts);
 271           if(myStoreMapping){
 272             myVTK2ObjIds.push_back(cellId); //apo
 273           }
 274           outputCD->CopyData(cd,cellId,newCellId);
 275           break;
 276
 277         case VTK_LINE:
 278         case VTK_POLY_LINE:
 279           newCellId = output->InsertNextCell(aCellType,npts,pts);
 280           if(myStoreMapping)
 281             myVTK2ObjIds.push_back(cellId);
 282           outputCD->CopyData(cd,cellId,newCellId);
 283           break;
 284
 285         case VTK_TRIANGLE:
 286         case VTK_QUAD:
 287         case VTK_POLYGON:
 288           newCellId = output->InsertNextCell(aCellType,npts,pts);
 289           if(myStoreMapping)
 290             myVTK2ObjIds.push_back(cellId);
 291           outputCD->CopyData(cd,cellId,newCellId);
 292           break;
 293
 294         case VTK_TRIANGLE_STRIP:
 295           newCellId = output->InsertNextCell(aCellType,npts,pts);
 296           if(myStoreMapping)
 297             myVTK2ObjIds.push_back(cellId);
 298           outputCD->CopyData(cd,cellId,newCellId);
 299           break;
 300
 301         case VTK_PIXEL:
 302           newCellId = output->InsertNextCell(aCellType,npts,pts);
 303           if(myStoreMapping)
 304             myVTK2ObjIds.push_back(cellId);
 305           outputCD->CopyData(cd,cellId,newCellId);
 306           break;
 307
 308         case VTK_CONVEX_POINT_SET: {
 309           //cout<<"cellId = "<<cellId<<"\n";
 310
 311           vtkIdType aNumPts;
 312           vtkPoints *aPoints;
 313 #ifdef USE_ROBUST_TRIANGULATION
 314           aPoints = aDelaunayPoints;
 315           anUnstructuredGrid->Initialize();
 316           anUnstructuredGrid->Allocate();
 317           anUnstructuredGrid->SetPoints(aDelaunayPoints);
 318
 319           vtkIdType *aPts;
 320           input->GetCellPoints(cellId,aNumPts,aPts);
 321           {
 322             float aPntCoord[3];
 323             aDelaunayPoints->SetNumberOfPoints(aNumPts);
 324             vtkPoints *anInputPoints = input->GetPoints();
 325             for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
 326               anInputPoints->GetPoint(aPts[aPntId],aPntCoord);
 327               aDelaunayPoints->SetPoint(aPntId,aPntCoord);
 328             }
 329           }
 330 #else
 331           input->GetCell(cellId,cell);
 332           aPoints = input->GetPoints();
 333           aNumPts = cell->GetNumberOfPoints();
 334 #endif
 335           // To calculate the bary center of the cell
 336           float aCellCenter[3] = {0.0, 0.0, 0.0};
 337           {
 338             float aPntCoord[3];
 339             for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
 340 #ifdef USE_ROBUST_TRIANGULATION
 341               aPoints->GetPoint(aPntId,aPntCoord);
 342 #else
 343               aPoints->GetPoint(cell->GetPointId(aPntId),aPntCoord);
 344 #endif
 345               //cout<<"\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}\n";
 346               aCellCenter[0] += aPntCoord[0];
 347               aCellCenter[1] += aPntCoord[1];
 348               aCellCenter[2] += aPntCoord[2];
 349             }
 350             aCellCenter[0] /= aNumPts;
 351             aCellCenter[1] /= aNumPts;
 352             aCellCenter[2] /= aNumPts;
 353           }
 354
 355 #ifdef USE_ROBUST_TRIANGULATION
 356           aGeometryFilter->Update();
 357           vtkPolyData* aPolyData = aGeometryFilter->GetOutput();
 358
 359           float aCellLength = aPolyData->GetLength();
 360           int aNumFaces = aPolyData->GetNumberOfCells();
 361 #else
 362           float aCellLength = sqrt(cell->GetLength2());
 363           int aNumFaces = cell->GetNumberOfFaces();
 364 #endif
 365
 366           static float EPS = 1.0E-5;
 367           float aDistEps = aCellLength * EPS;
 368
 369           // To initialize set of points that belong to the cell
 370           typedef std::set<vtkIdType> TPointIds;
 371           TPointIds anInitialPointIds;
 372           for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++){
 373 #ifdef USE_ROBUST_TRIANGULATION
 374             anInitialPointIds.insert(aPntId);
 375 #else
 376             anInitialPointIds.insert(cell->GetPointId(aPntId));
 377 #endif
 378           }
 379
 380           // To initialize set of points by face that belong to the cell and backward
 381           typedef std::set<vtkIdType> TFace2Visibility;
 382           TFace2Visibility aFace2Visibility;
 383
 384           typedef std::set<TPointIds> TFace2PointIds;
 385           TFace2PointIds aFace2PointIds;
 386
 387           for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) {
 388 #ifdef USE_ROBUST_TRIANGULATION
 389             vtkCell* aFace = aPolyData->GetCell(aFaceId);
 390 #else
 391             vtkCell* aFace = cell->GetFace(aFaceId);
 392 #endif
 393             vtkIdList *anIdList = aFace->PointIds;
 394             aNewPts[0] = anIdList->GetId(0);
 395             aNewPts[1] = anIdList->GetId(1);
 396             aNewPts[2] = anIdList->GetId(2);
 397
 398 #ifdef USE_ROBUST_TRIANGULATION
 399             faceIds->Reset();
 400             faceIds->InsertNextId(aPts[aNewPts[0]]);
 401             faceIds->InsertNextId(aPts[aNewPts[1]]);
 402             faceIds->InsertNextId(aPts[aNewPts[2]]);
 403             input->GetCellNeighbors(cellId, faceIds, cellIds);
 404 #else
 405             input->GetCellNeighbors(cellId, anIdList, cellIds);
 406 #endif
 407             if((!allVisible && !cellVis[cellIds->GetId(0)]) ||
 408                cellIds->GetNumberOfIds() <= 0 ||
 409                myShowInside)
 410             {
 411               TPointIds aPointIds;
 412               aPointIds.insert(aNewPts[0]);
 413               aPointIds.insert(aNewPts[1]);
 414               aPointIds.insert(aNewPts[2]);
 415
 416               aFace2PointIds.insert(aPointIds);
 417               aFace2Visibility.insert(aFaceId);
 418             }
 419           }
 420
 421           for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) {
 422             if(aFace2Visibility.find(aFaceId) == aFace2Visibility.end())
 423               continue;
 424
 425 #ifdef USE_ROBUST_TRIANGULATION
 426             vtkCell* aFace = aPolyData->GetCell(aFaceId);
 427 #else
 428             vtkCell* aFace = cell->GetFace(aFaceId);
 429 #endif
 430             vtkIdList *anIdList = aFace->PointIds;
 431             aNewPts[0] = anIdList->GetId(0);
 432             aNewPts[1] = anIdList->GetId(1);
 433             aNewPts[2] = anIdList->GetId(2);
 434
 435             // To initialize set of points for the plane where the trinangle face belong to
 436             TPointIds aPointIds;
 437             aPointIds.insert(aNewPts[0]);
 438             aPointIds.insert(aNewPts[1]);
 439             aPointIds.insert(aNewPts[2]);
 440
 441             // To get know, if the points of the trinagle were already observed
 442             bool anIsObserved = aFace2PointIds.find(aPointIds) == aFace2PointIds.end();
 443             //cout<<"\taFaceId = "<<aFaceId<<"; anIsObserved = "<<anIsObserved;
 444             //cout<<"; aNewPts = {"<<aNewPts[0]<<", "<<aNewPts[1]<<", "<<aNewPts[2]<<"}\n";
 445
 446             if(!anIsObserved){
 447               // To get coordinates of the points of the traingle face
 448               float aCoord[3][3];
 449               aPoints->GetPoint(aNewPts[0],aCoord[0]);
 450               aPoints->GetPoint(aNewPts[1],aCoord[1]);
 451               aPoints->GetPoint(aNewPts[2],aCoord[2]);
 452
 453               // To calculate plane normal
 454               float aVector01[3] = { aCoord[1][0] - aCoord[0][0],
 455                                      aCoord[1][1] - aCoord[0][1],
 456                                      aCoord[1][2] - aCoord[0][2] };
 457
 458               float aVector02[3] = { aCoord[2][0] - aCoord[0][0],
 459                                      aCoord[2][1] - aCoord[0][1],
 460                                      aCoord[2][2] - aCoord[0][2] };
 461
 462               float aCross21[3];
 463               vtkMath::Cross(aVector02,aVector01,aCross21);
 464
 465               vtkMath::Normalize(aCross21);
 466
 467               // To calculate what points belong to the plane
 468               // To calculate bounds of the point set
 469               float aCenter[3] = {0.0, 0.0, 0.0};
 470               {
 471                 TPointIds::const_iterator anIter = anInitialPointIds.begin();
 472                 TPointIds::const_iterator anEndIter = anInitialPointIds.end();
 473                 for(; anIter != anEndIter; anIter++){
 474                   float aPntCoord[3];
 475                   vtkIdType aPntId = *anIter;
 476                   aPoints->GetPoint(aPntId,aPntCoord);
 477                   float aDist = vtkPlane::DistanceToPlane(aPntCoord,aCross21,aCoord[0]);
 478                   //cout<<"\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}; aDist = "<<aDist<<"\n";
 479                   if(fabs(aDist) < aDistEps){
 480                     aPointIds.insert(aPntId);
 481                     aCenter[0] += aPntCoord[0];
 482                     aCenter[1] += aPntCoord[1];
 483                     aCenter[2] += aPntCoord[2];
 484                   }
 485                 }
 486                 int aNbPoints = aPointIds.size();
 487                 aCenter[0] /= aNbPoints;
 488                 aCenter[1] /= aNbPoints;
 489                 aCenter[2] /= aNbPoints;
 490               }
 491
 492               // To calculate the primary direction for point set
 493               float aVector0[3] = { aCoord[0][0] - aCenter[0],
 494                                     aCoord[0][1] - aCenter[1],
 495                                     aCoord[0][2] - aCenter[2] };
 496
 497               //To sinchronize orientation of the cell and its face
 498               float aVectorC[3] = { aCenter[0] - aCellCenter[0],
 499                                     aCenter[1] - aCellCenter[1],
 500                                     aCenter[2] - aCellCenter[2] };
 501               vtkMath::Normalize(aVectorC);
 502
 503               float aDot = vtkMath::Dot(aCross21,aVectorC);
 504               //cout<<"\t\taCross21 = {"<<aCross21[0]<<", "<<aCross21[1]<<", "<<aCross21[2]<<"}";
 505               //cout<<"; aVectorC = {"<<aVectorC[0]<<", "<<aVectorC[1]<<", "<<aVectorC[2]<<"}\n";
 506               //cout<<"\t\taDot = "<<aDot<<"\n";
 507               if(aDot > 0){
 508                 aCross21[0] = -aCross21[0];
 509                 aCross21[1] = -aCross21[1];
 510                 aCross21[2] = -aCross21[2];
 511               }
 512
 513               vtkMath::Normalize(aVector0);
 514
 515               //cout<<"\t\taCenter = {"<<aCenter[0]<<", "<<aCenter[1]<<", "<<aCenter[2]<<"}";
 516               //cout<<"; aVector0 = {"<<aVector0[0]<<", "<<aVector0[1]<<", "<<aVector0[2]<<"}\n";
 517
 518               // To calculate the set of points by face those that belong to the plane
 519               TFace2PointIds aRemoveFace2PointIds;
 520               {
 521                 TFace2PointIds::const_iterator anIter = aFace2PointIds.begin();
 522                 TFace2PointIds::const_iterator anEndIter = aFace2PointIds.end();
 523                 for(; anIter != anEndIter; anIter++){
 524                   const TPointIds& anIds = *anIter;
 525                   TPointIds anIntersection;
 526                   std::set_intersection(aPointIds.begin(),aPointIds.end(),
 527                                         anIds.begin(),anIds.end(),
 528                                         std::inserter(anIntersection,anIntersection.begin()));
 529
 530                   if(anIntersection == anIds){
 531                     aRemoveFace2PointIds.insert(anIds);
 532                   }
 533                 }
 534               }
 535
 536               // To remove from the set of points by face those that belong to the plane
 537               {
 538                 TFace2PointIds::const_iterator anIter = aRemoveFace2PointIds.begin();
 539                 TFace2PointIds::const_iterator anEndIter = aRemoveFace2PointIds.end();
 540                 for(; anIter != anEndIter; anIter++){
 541                   const TPointIds& anIds = *anIter;
 542                   aFace2PointIds.erase(anIds);
 543                 }
 544               }
 545
 546               // To sort the planar set of the points accrding to the angle
 547               {
 548                 typedef std::map<float,vtkIdType> TSortedPointIds;
 549                 TSortedPointIds aSortedPointIds;
 550
 551                 TPointIds::const_iterator anIter = aPointIds.begin();
 552                 TPointIds::const_iterator anEndIter = aPointIds.end();
 553                 for(; anIter != anEndIter; anIter++){
 554                   float aPntCoord[3];
 555                   vtkIdType aPntId = *anIter;
 556                   aPoints->GetPoint(aPntId,aPntCoord);
 557                   float aVector[3] = { aPntCoord[0] - aCenter[0],
 558                                        aPntCoord[1] - aCenter[1],
 559                                        aPntCoord[2] - aCenter[2] };
 560                   vtkMath::Normalize(aVector);
 561
 562                   float aCross[3];
 563                   vtkMath::Cross(aVector,aVector0,aCross);
 564                   bool aGreaterThanPi = vtkMath::Dot(aCross,aCross21) < 0;
 565                   float aCosinus = vtkMath::Dot(aVector,aVector0);
 566                   if(aCosinus > 1.0)
 567                     aCosinus = 1.0;
 568                   if(aCosinus < -1.0)
 569                     aCosinus = -1.0;
 570                   static float a2Pi = 2.0 * vtkMath::Pi();
 571                   float anAngle = acos(aCosinus);
 572                   //cout<<"\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}";
 573                   //cout<<"; aGreaterThanPi = "<<aGreaterThanPi<<"; aCosinus = "<<aCosinus<<"; anAngle = "<<anAngle<<"\n";
 574                   if(aGreaterThanPi)
 575                     anAngle = a2Pi - anAngle;
 576                   aSortedPointIds[anAngle] = aPntId;
 577                   //cout<<"\t\t\tanAngle = "<<anAngle<<"\n";
 578                 }
 579                 if(!aSortedPointIds.empty()){
 580                   aCellType = VTK_POLYGON;
 581                   int numFacePts = aSortedPointIds.size();
 582                   std::vector<vtkIdType> aConnectivities(numFacePts);
 583                   TSortedPointIds::const_iterator anIter = aSortedPointIds.begin();
 584                   TSortedPointIds::const_iterator anEndIter = aSortedPointIds.end();
 585                   for(vtkIdType anId = 0; anIter != anEndIter; anIter++, anId++){
 586                     vtkIdType aPntId = anIter->second;
 587 #ifdef USE_ROBUST_TRIANGULATION
 588                     aConnectivities[anId] = aPts[aPntId];
 589 #else
 590                     aConnectivities[anId] = aPntId;
 591 #endif
 592                   }
 593                   newCellId = output->InsertNextCell(aCellType,numFacePts,&aConnectivities[0]);
 594                   if(myStoreMapping)
 595                     myVTK2ObjIds.push_back(cellId);
 596                   outputCD->CopyData(cd,cellId,newCellId);
 597                 }
 598               }
 599             }
 600           }
 601
 602           break;
 603         }
 604         case VTK_TETRA: {
 605           for (faceId = 0; faceId < 4; faceId++)
 606             {
 607             faceIds->Reset();
 608             faceVerts = vtkTetra::GetFaceArray(faceId);
 609             faceIds->InsertNextId(pts[faceVerts[0]]);
 610             faceIds->InsertNextId(pts[faceVerts[1]]);
 611             faceIds->InsertNextId(pts[faceVerts[2]]);
 612             aCellType = VTK_TRIANGLE;
 613             numFacePts = 3;
 614             input->GetCellNeighbors(cellId, faceIds, cellIds);
 615             if ( cellIds->GetNumberOfIds() <= 0 || myShowInside ||
 616                  (!allVisible && !cellVis[cellIds->GetId(0)]) )
 617               {
 618               for ( i=0; i < numFacePts; i++)
 619                 aNewPts[i] = pts[faceVerts[i]];
 620               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 621               if(myStoreMapping)
 622                 myVTK2ObjIds.push_back(cellId);
 623               outputCD->CopyData(cd,cellId,newCellId);
 624               }
 625             }
 626           break;
 627         }
 628         case VTK_VOXEL: {
 629           for (faceId = 0; faceId < 6; faceId++)
 630             {
 631             faceIds->Reset();
 632             faceVerts = vtkVoxel::GetFaceArray(faceId);
 633             faceIds->InsertNextId(pts[faceVerts[0]]);
 634             faceIds->InsertNextId(pts[faceVerts[1]]);
 635             faceIds->InsertNextId(pts[faceVerts[2]]);
 636             faceIds->InsertNextId(pts[faceVerts[3]]);
 637             aCellType = VTK_QUAD;
 638             numFacePts = 4;
 639             input->GetCellNeighbors(cellId, faceIds, cellIds);
 640             if ( cellIds->GetNumberOfIds() <= 0 || myShowInside ||
 641                  (!allVisible && !cellVis[cellIds->GetId(0)]) )
 642               {
 643               for ( i=0; i < numFacePts; i++)
 644                 aNewPts[i] = pts[faceVerts[PixelConvert[i]]];
 645               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 646               if(myStoreMapping)
 647                 myVTK2ObjIds.push_back(cellId);
 648               outputCD->CopyData(cd,cellId,newCellId);
 649               }
 650             }
 651           break;
 652         }
 653         case VTK_HEXAHEDRON: {
 654           for (faceId = 0; faceId < 6; faceId++)
 655             {
 656             faceIds->Reset();
 657             faceVerts = vtkHexahedron::GetFaceArray(faceId);
 658             faceIds->InsertNextId(pts[faceVerts[0]]);
 659             faceIds->InsertNextId(pts[faceVerts[1]]);
 660             faceIds->InsertNextId(pts[faceVerts[2]]);
 661             faceIds->InsertNextId(pts[faceVerts[3]]);
 662             aCellType = VTK_QUAD;
 663             numFacePts = 4;
 664             input->GetCellNeighbors(cellId, faceIds, cellIds);
 665             if ( cellIds->GetNumberOfIds() <= 0 || myShowInside ||
 666                  (!allVisible && !cellVis[cellIds->GetId(0)]) )
 667               {
 668               for ( i=0; i < numFacePts; i++)
 669                 aNewPts[i] = pts[faceVerts[i]];
 670               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 671               if(myStoreMapping)
 672                 myVTK2ObjIds.push_back(cellId);
 673               outputCD->CopyData(cd,cellId,newCellId);
 674               }
 675             }
 676           break;
 677         }
 678         case VTK_WEDGE: {
 679           for (faceId = 0; faceId < 5; faceId++)
 680             {
 681             faceIds->Reset();
 682             faceVerts = vtkWedge::GetFaceArray(faceId);
 683             faceIds->InsertNextId(pts[faceVerts[0]]);
 684             faceIds->InsertNextId(pts[faceVerts[1]]);
 685             faceIds->InsertNextId(pts[faceVerts[2]]);
 686             aCellType = VTK_TRIANGLE;
 687             numFacePts = 3;
 688             if (faceVerts[3] >= 0)
 689               {
 690               faceIds->InsertNextId(pts[faceVerts[3]]);
 691               aCellType = VTK_QUAD;
 692               numFacePts = 4;
 693               }
 694             input->GetCellNeighbors(cellId, faceIds, cellIds);
 695             if ( cellIds->GetNumberOfIds() <= 0 || myShowInside ||
 696                  (!allVisible && !cellVis[cellIds->GetId(0)]) )
 697               {
 698               for ( i=0; i < numFacePts; i++)
 699                 aNewPts[i] = pts[faceVerts[i]];
 700               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 701               if(myStoreMapping)
 702                 myVTK2ObjIds.push_back(cellId);
 703               outputCD->CopyData(cd,cellId,newCellId);
 704               }
 705             }
 706           break;
 707         }
 708         case VTK_PYRAMID: {
 709           for (faceId = 0; faceId < 5; faceId++)
 710             {
 711             faceIds->Reset();
 712             faceVerts = vtkPyramid::GetFaceArray(faceId);
 713             faceIds->InsertNextId(pts[faceVerts[0]]);
 714             faceIds->InsertNextId(pts[faceVerts[1]]);
 715             faceIds->InsertNextId(pts[faceVerts[2]]);
 716             aCellType = VTK_TRIANGLE;
 717             numFacePts = 3;
 718             if (faceVerts[3] >= 0)
 719               {
 720               faceIds->InsertNextId(pts[faceVerts[3]]);
 721               aCellType = VTK_QUAD;
 722               numFacePts = 4;
 723               }
 724             input->GetCellNeighbors(cellId, faceIds, cellIds);
 725             if ( cellIds->GetNumberOfIds() <= 0 || myShowInside ||
 726                  (!allVisible && !cellVis[cellIds->GetId(0)]) )
 727               {
 728               for ( i=0; i < numFacePts; i++)
 729                 aNewPts[i] = pts[faceVerts[i]];
 730               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 731               if(myStoreMapping)
 732                 myVTK2ObjIds.push_back(cellId);
 733               outputCD->CopyData(cd,cellId,newCellId);
 734               }
 735             }
 736           break;
 737         }
 738         //Quadratic cells
 739         case VTK_QUADRATIC_EDGE:
 740         case VTK_QUADRATIC_TRIANGLE:
 741         case VTK_QUADRATIC_QUAD:
 742         case VTK_QUADRATIC_TETRA:
 743         case VTK_QUADRATIC_HEXAHEDRON:
 744           if(!myIsWireframeMode){
 745             input->GetCell(cellId,cell);
 746             vtkIdList *pts = vtkIdList::New();
 747             vtkPoints *coords = vtkPoints::New();
 748             vtkIdList *cellIds = vtkIdList::New();
 749             vtkIdType newCellId;
 750
 751             if ( cell->GetCellDimension() == 1 ) {
 752               aCellType = VTK_LINE;
 753               numFacePts = 2;
 754               cell->Triangulate(0,pts,coords);
 755               for (i=0; i < pts->GetNumberOfIds(); i+=2) {
 756                 aNewPts[0] = pts->GetId(i);
 757                 aNewPts[1] = pts->GetId(i+1);
 758                 newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 759                 if(myStoreMapping)
 760                   myVTK2ObjIds.push_back(cellId);
 761                 outputCD->CopyData(cd,cellId,newCellId);
 762               }
 763             }
 764             else if ( cell->GetCellDimension() == 2 ) {
 765               aCellType = VTK_TRIANGLE;
 766               numFacePts = 3;
 767               cell->Triangulate(0,pts,coords);
 768               for (i=0; i < pts->GetNumberOfIds(); i+=3) {
 769                 aNewPts[0] = pts->GetId(i);
 770                 aNewPts[1] = pts->GetId(i+1);
 771                 aNewPts[2] = pts->GetId(i+2);
 772                 newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 773                 if(myStoreMapping)
 774                   myVTK2ObjIds.push_back(cellId);
 775                 outputCD->CopyData(cd,cellId,newCellId);
 776               }
 777             }
 778             else //3D nonlinear cell
 779             {
 780               aCellType = VTK_TRIANGLE;
 781               numFacePts = 3;
 782               for (int j=0; j < cell->GetNumberOfFaces(); j++){
 783                 vtkCell *face = cell->GetFace(j);
 784                 input->GetCellNeighbors(cellId, face->PointIds, cellIds);
 785                 if ( cellIds->GetNumberOfIds() <= 0 || myShowInside ) {
 786                   face->Triangulate(0,pts,coords);
 787                   for (i=0; i < pts->GetNumberOfIds(); i+=3) {
 788                     aNewPts[0] = pts->GetId(i);
 789                     aNewPts[1] = pts->GetId(i+1);
 790                     aNewPts[2] = pts->GetId(i+2);
 791                     newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 792                     if(myStoreMapping)
 793                       myVTK2ObjIds.push_back(cellId);
 794                     outputCD->CopyData(cd,cellId,newCellId);
 795                   }
 796                 }
 797               }
 798             } //3d cell
 799             cellIds->Delete();
 800             coords->Delete();
 801             pts->Delete();
 802             break;
 803           }else{
 804             switch(aCellType){
 805             case VTK_QUADRATIC_EDGE: {
 806               aCellType = VTK_POLY_LINE;
 807               numFacePts = 3;
 808
 809               aNewPts[0] = pts[0];
 810               aNewPts[2] = pts[1];
 811               aNewPts[1] = pts[2];
 812
 813               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 814               if(myStoreMapping)
 815                 myVTK2ObjIds.push_back(cellId);
 816
 817               outputCD->CopyData(cd,cellId,newCellId);
 818               break;
 819             }
 820             case VTK_QUADRATIC_TRIANGLE: {
 821               aCellType = VTK_POLYGON;
 822               numFacePts = 6;
 823
 824               aNewPts[0] = pts[0];
 825               aNewPts[1] = pts[3];
 826               aNewPts[2] = pts[1];
 827               aNewPts[3] = pts[4];
 828               aNewPts[4] = pts[2];
 829               aNewPts[5] = pts[5];
 830
 831               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 832               if(myStoreMapping)
 833                 myVTK2ObjIds.push_back(cellId);
 834
 835               outputCD->CopyData(cd,cellId,newCellId);
 836               break;
 837             }
 838             case VTK_QUADRATIC_QUAD: {
 839               aCellType = VTK_POLYGON;
 840               numFacePts = 8;
 841
 842               aNewPts[0] = pts[0];
 843               aNewPts[1] = pts[4];
 844               aNewPts[2] = pts[1];
 845               aNewPts[3] = pts[5];
 846               aNewPts[4] = pts[2];
 847               aNewPts[5] = pts[6];
 848               aNewPts[6] = pts[3];
 849               aNewPts[7] = pts[7];
 850
 851               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 852               if(myStoreMapping)
 853                 myVTK2ObjIds.push_back(cellId);
 854
 855               outputCD->CopyData(cd,cellId,newCellId);
 856               break;
 857             }
 858             case VTK_QUADRATIC_TETRA: {
 859               aCellType = VTK_POLYGON;
 860               numFacePts = 6;
 861
 862               //---------------------------------------------------------------
 863               aNewPts[0] = pts[0];
 864               aNewPts[1] = pts[4];
 865               aNewPts[2] = pts[1];
 866               aNewPts[3] = pts[5];
 867               aNewPts[4] = pts[2];
 868               aNewPts[5] = pts[6];
 869
 870               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 871               if(myStoreMapping)
 872                 myVTK2ObjIds.push_back(cellId);
 873
 874               outputCD->CopyData(cd,cellId,newCellId);
 875
 876               //---------------------------------------------------------------
 877               aNewPts[0] = pts[0];
 878               aNewPts[1] = pts[7];
 879               aNewPts[2] = pts[3];
 880               aNewPts[3] = pts[8];
 881               aNewPts[4] = pts[1];
 882               aNewPts[5] = pts[4];
 883
 884               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 885               if(myStoreMapping)
 886                 myVTK2ObjIds.push_back(cellId);
 887
 888               outputCD->CopyData(cd,cellId,newCellId);
 889
 890               //---------------------------------------------------------------
 891               aNewPts[0] = pts[1];
 892               aNewPts[1] = pts[8];
 893               aNewPts[2] = pts[3];
 894               aNewPts[3] = pts[9];
 895               aNewPts[4] = pts[2];
 896               aNewPts[5] = pts[5];
 897
 898               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 899               if(myStoreMapping)
 900                 myVTK2ObjIds.push_back(cellId);
 901
 902               outputCD->CopyData(cd,cellId,newCellId);
 903
 904               //---------------------------------------------------------------
 905               aNewPts[0] = pts[2];
 906               aNewPts[1] = pts[9];
 907               aNewPts[2] = pts[3];
 908               aNewPts[3] = pts[7];
 909               aNewPts[4] = pts[0];
 910               aNewPts[5] = pts[6];
 911
 912               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 913               if(myStoreMapping)
 914                 myVTK2ObjIds.push_back(cellId);
 915
 916               outputCD->CopyData(cd,cellId,newCellId);
 917
 918               break;
 919             }
 920             case VTK_QUADRATIC_HEXAHEDRON: {
 921               aCellType = VTK_POLYGON;
 922               numFacePts = 8;
 923
 924               //---------------------------------------------------------------
 925               aNewPts[0] = pts[0];
 926               aNewPts[1] = pts[8];
 927               aNewPts[2] = pts[1];
 928               aNewPts[3] = pts[17];
 929               aNewPts[4] = pts[5];
 930               aNewPts[5] = pts[12];
 931               aNewPts[6] = pts[4];
 932               aNewPts[7] = pts[16];
 933
 934               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 935               if(myStoreMapping)
 936                 myVTK2ObjIds.push_back(cellId);
 937
 938               outputCD->CopyData(cd,cellId,newCellId);
 939
 940               //---------------------------------------------------------------
 941               aNewPts[0] = pts[1];
 942               aNewPts[1] = pts[9];
 943               aNewPts[2] = pts[2];
 944               aNewPts[3] = pts[18];
 945               aNewPts[4] = pts[6];
 946               aNewPts[5] = pts[13];
 947               aNewPts[6] = pts[5];
 948               aNewPts[7] = pts[17];
 949
 950               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 951               if(myStoreMapping)
 952                 myVTK2ObjIds.push_back(cellId);
 953
 954               outputCD->CopyData(cd,cellId,newCellId);
 955
 956               //---------------------------------------------------------------
 957               aNewPts[0] = pts[2];
 958               aNewPts[1] = pts[10];
 959               aNewPts[2] = pts[3];
 960               aNewPts[3] = pts[19];
 961               aNewPts[4] = pts[7];
 962               aNewPts[5] = pts[14];
 963               aNewPts[6] = pts[6];
 964               aNewPts[7] = pts[18];
 965
 966               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 967               if(myStoreMapping)
 968                 myVTK2ObjIds.push_back(cellId);
 969
 970               outputCD->CopyData(cd,cellId,newCellId);
 971
 972               //---------------------------------------------------------------
 973               aNewPts[0] = pts[3];
 974               aNewPts[1] = pts[11];
 975               aNewPts[2] = pts[0];
 976               aNewPts[3] = pts[16];
 977               aNewPts[4] = pts[4];
 978               aNewPts[5] = pts[15];
 979               aNewPts[6] = pts[7];
 980               aNewPts[7] = pts[19];
 981
 982               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 983               if(myStoreMapping)
 984                 myVTK2ObjIds.push_back(cellId);
 985
 986               outputCD->CopyData(cd,cellId,newCellId);
 987
 988               //---------------------------------------------------------------
 989               aNewPts[0] = pts[0];
 990               aNewPts[1] = pts[8];
 991               aNewPts[2] = pts[1];
 992               aNewPts[3] = pts[9];
 993               aNewPts[4] = pts[2];
 994               aNewPts[5] = pts[10];
 995               aNewPts[6] = pts[3];
 996               aNewPts[7] = pts[11];
 997
 998               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
 999               if(myStoreMapping)
1000                 myVTK2ObjIds.push_back(cellId);
1001
1002               outputCD->CopyData(cd,cellId,newCellId);
1003
1004               //---------------------------------------------------------------
1005               aNewPts[0] = pts[4];
1006               aNewPts[1] = pts[12];
1007               aNewPts[2] = pts[5];
1008               aNewPts[3] = pts[13];
1009               aNewPts[4] = pts[6];
1010               aNewPts[5] = pts[14];
1011               aNewPts[6] = pts[7];
1012               aNewPts[7] = pts[15];
1013
1014               newCellId = output->InsertNextCell(aCellType,numFacePts,aNewPts);
1015               if(myStoreMapping)
1016                 myVTK2ObjIds.push_back(cellId);
1017
1018               outputCD->CopyData(cd,cellId,newCellId);
1019
1020               break;
1021             }}
1022           }
1023         } //switch
1024       } //if visible
1025     } //for all cells
1026
1027   output->Squeeze();
1028
1029   vtkDebugMacro(<<"Extracted " << input->GetNumberOfPoints() << " points,"
1030   << output->GetNumberOfCells() << " cells.");
1031
1032 #ifdef USE_ROBUST_TRIANGULATION
1033   anUnstructuredGrid->Delete();
1034   aDelaunayPoints->Delete();
1035
1036   aDelaunay3D->Delete();
1037   aGeometryFilter->Delete();
1038 #endif
1039
1040   cell->Delete();
1041
1042   cellIds->Delete();
1043   faceIds->Delete();
1044
1045   if ( cellVis )
1046     {
1047     delete [] cellVis;
1048     }
1049 }
1050
1051
1052 //----------------------------------------------------------------------------
1053 void
1054 VTKViewer_GeometryFilter
1055 ::SetInside(int theShowInside)
1056 {
1057   if(myShowInside == theShowInside)
1058     return;
1059
1060   myShowInside = theShowInside;
1061   this->Modified();
1062 }
1063
1064 int
1065 VTKViewer_GeometryFilter
1066 ::GetInside()
1067 {
1068   return myShowInside;
1069 }
1070
1071
1072 //----------------------------------------------------------------------------
1073 void
1074 VTKViewer_GeometryFilter
1075 ::SetWireframeMode(int theIsWireframeMode)
1076 {
1077   if(myIsWireframeMode == theIsWireframeMode)
1078     return;
1079
1080   myIsWireframeMode = theIsWireframeMode;
1081   this->Modified();
1082 }
1083
1084 int
1085 VTKViewer_GeometryFilter
1086 ::GetWireframeMode()
1087 {
1088   return myIsWireframeMode;
1089 }
1090
1091
1092 //----------------------------------------------------------------------------
1093 void
1094 VTKViewer_GeometryFilter
1095 ::SetStoreMapping(int theStoreMapping)
1096 {
1097   if(myStoreMapping == theStoreMapping)
1098     return;
1099
1100   myStoreMapping = theStoreMapping;
1101   this->Modified();
1102 }
1103
1104 int
1105 VTKViewer_GeometryFilter
1106 ::GetStoreMapping()
1107 {
1108   return myStoreMapping;
1109 }
1110
1111
1112 //----------------------------------------------------------------------------
1113 vtkIdType VTKViewer_GeometryFilter::GetElemObjId(int theVtkID){
1114   if(myVTK2ObjIds.empty() || theVtkID > myVTK2ObjIds.size()) return -1;
1115 #if defined __GNUC_2__
1116   return myVTK2ObjIds[theVtkID];
1117 #else
1118   return myVTK2ObjIds.at(theVtkID);
1119 #endif
1120 }