src/VTKViewer/VTKViewer_ConvexTool.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  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.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "VTKViewer_ConvexTool.h"
  24 #include "VTKViewer_GeometryFilter.h"
  25
  26 #include <set>
  27 #include <map>
  28 #include <algorithm>
  29
  30 #include <vtkUnstructuredGrid.h>
  31 #include <vtkGeometryFilter.h>
  32 #include <vtkDelaunay3D.h>
  33 #include <vtkGenericCell.h>
  34 #include <vtkPolyData.h>
  35 #include <vtkCellData.h>
  36 #include <vtkPoints.h>
  37 #include <vtkIdList.h>
  38 #include <vtkCell.h>
  39 #include <vtkPlane.h>
  40 #include <vtkMath.h>
  41 #include <vtkCellArray.h>
  42 #include <vtkTriangle.h>
  43 #include <vtkOrderedTriangulator.h>
  44
  45 #ifdef _DEBUG_
  46 static int DEBUG_TRIA_EXECUTE = 0;
  47 #else
  48 static int DEBUG_TRIA_EXECUTE = 0;
  49 #endif
  50
  51 namespace
  52 {
  53   typedef std::vector<vtkIdType> TConnectivities;
  54
  55   struct TPolygon
  56   {
  57     TConnectivities myConnectivities;
  58     vtkFloatingPointType myOrigin[3];
  59     vtkFloatingPointType myNormal[3];
  60     TPolygon(const TConnectivities& theConnectivities,
  61              vtkFloatingPointType theOrigin[3],
  62              vtkFloatingPointType theNormal[3]):
  63       myConnectivities(theConnectivities)
  64     {
  65       myOrigin[0] = theOrigin[0];
  66       myOrigin[1] = theOrigin[1];
  67       myOrigin[2] = theOrigin[2];
  68
  69       myNormal[0] = theNormal[0];
  70       myNormal[1] = theNormal[1];
  71       myNormal[2] = theNormal[2];
  72     }
  73   };
  74
  75   typedef std::vector<TPolygon> TPolygons;
  76 }
  77
  78
  79 //----------------------------------------------------------------------------
  80 VTKViewer_Triangulator
  81 ::VTKViewer_Triangulator():
  82   myCellIds(vtkIdList::New()),
  83   myFaceIds(vtkIdList::New()),
  84   myPoints(vtkPoints::New()),
  85   myPointIds(NULL)
  86 {}
  87
  88
  89 //----------------------------------------------------------------------------
  90 VTKViewer_Triangulator
  91 ::~VTKViewer_Triangulator()
  92 {
  93   myCellIds->Delete();
  94   myFaceIds->Delete();
  95   myPoints->Delete();
  96 }
  97
  98
  99 //----------------------------------------------------------------------------
 100 vtkPoints*
 101 VTKViewer_Triangulator
 102 ::InitPoints(vtkUnstructuredGrid *theInput,
 103              vtkIdType theCellId)
 104 {
 105   myPoints->Reset();
 106   myPoints->Modified(); // the VTK bug
 107
 108   vtkIdType aNumPts;
 109   theInput->GetCellPoints(theCellId, aNumPts, myPointIds);
 110   if ( aNumPts > 0 ) {
 111     vtkFloatingPointType anAbsoluteCoord[3];
 112     myPoints->SetNumberOfPoints(aNumPts);
 113     vtkPoints *anInputPoints = theInput->GetPoints();
 114     for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
 115       anInputPoints->GetPoint(myPointIds[aPntId], anAbsoluteCoord);
 116       myPoints->SetPoint(aPntId, anAbsoluteCoord);
 117     }
 118   }
 119
 120   return myPoints;
 121 }
 122
 123
 124 //----------------------------------------------------------------------------
 125 vtkIdType
 126 VTKViewer_Triangulator
 127 ::GetNbOfPoints()
 128 {
 129   return myPoints->GetNumberOfPoints();
 130 }
 131
 132
 133 //----------------------------------------------------------------------------
 134 vtkIdType
 135 VTKViewer_Triangulator
 136 ::GetPointId(vtkIdType thePointId)
 137 {
 138   return thePointId;
 139 }
 140
 141
 142 //----------------------------------------------------------------------------
 143 vtkFloatingPointType
 144 VTKViewer_Triangulator
 145 ::GetCellLength()
 146 {
 147   vtkFloatingPointType aBounds[6];
 148   myPoints->GetBounds(aBounds);
 149
 150   vtkFloatingPointType aCoordDiff[3];
 151   aCoordDiff[0] = (aBounds[1] - aBounds[0]);
 152   aCoordDiff[1] = (aBounds[3] - aBounds[2]);
 153   aCoordDiff[2] = (aBounds[5] - aBounds[4]);
 154
 155   return sqrt(aCoordDiff[0]*aCoordDiff[0] +
 156               aCoordDiff[1]*aCoordDiff[1] +
 157               aCoordDiff[2]*aCoordDiff[2]);
 158 }
 159
 160
 161 //----------------------------------------------------------------------------
 162 void
 163 VTKViewer_Triangulator
 164 ::GetCellNeighbors(vtkUnstructuredGrid *theInput,
 165                    vtkIdType theCellId,
 166                    vtkCell* theFace,
 167                    vtkIdList* theCellIds)
 168 {
 169   myFaceIds->Reset();
 170   vtkIdList *anIdList = theFace->PointIds;
 171   myFaceIds->InsertNextId(myPointIds[anIdList->GetId(0)]);
 172   myFaceIds->InsertNextId(myPointIds[anIdList->GetId(1)]);
 173   myFaceIds->InsertNextId(myPointIds[anIdList->GetId(2)]);
 174
 175   theInput->GetCellNeighbors(theCellId, myFaceIds, theCellIds);
 176 }
 177
 178
 179 //----------------------------------------------------------------------------
 180 vtkIdType
 181 VTKViewer_Triangulator
 182 ::GetConnectivity(vtkIdType thePntId)
 183 {
 184   return myPointIds[thePntId];
 185 }
 186
 187
 188 //----------------------------------------------------------------------------
 189 bool
 190 VTKViewer_Triangulator
 191 ::Execute(vtkUnstructuredGrid *theInput,
 192           vtkCellData* thInputCD,
 193           vtkIdType theCellId,
 194           int theShowInside,
 195           int theAllVisible,
 196           int theAppendCoincident3D,
 197           const char* theCellsVisibility,
 198           vtkPolyData *theOutput,
 199           vtkCellData* theOutputCD,
 200           int theStoreMapping,
 201           std::vector<vtkIdType>& theVTK2ObjIds,
 202           std::map< vtkIdType, std::vector<vtkIdType> >& theDimension2VTK2ObjIds,
 203           bool theIsCheckConvex)
 204 {
 205   vtkPoints *aPoints = InitPoints(theInput, theCellId);
 206   vtkIdType aNumPts = GetNbOfPoints();
 207   if(DEBUG_TRIA_EXECUTE) cout<<"Triangulator - aNumPts = "<<aNumPts<<"\n";
 208
 209   if(aNumPts == 0)
 210     return true;
 211
 212   // To calculate the bary center of the cell
 213   vtkFloatingPointType aCellCenter[3] = {0.0, 0.0, 0.0};
 214   {
 215     vtkFloatingPointType aPntCoord[3];
 216     for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
 217       aPoints->GetPoint(GetPointId(aPntId),aPntCoord);
 218       if(DEBUG_TRIA_EXECUTE) cout<<"\taPntId = "<<GetPointId(aPntId)<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}\n";
 219       aCellCenter[0] += aPntCoord[0];
 220       aCellCenter[1] += aPntCoord[1];
 221       aCellCenter[2] += aPntCoord[2];
 222     }
 223     aCellCenter[0] /= aNumPts;
 224     aCellCenter[1] /= aNumPts;
 225     aCellCenter[2] /= aNumPts;
 226   }
 227
 228   vtkFloatingPointType aCellLength = GetCellLength();
 229   int aNumFaces = GetNumFaces();
 230
 231   static vtkFloatingPointType EPS = 1.0E-2;
 232   vtkFloatingPointType aDistEps = aCellLength/3.0 * EPS;
 233   if(DEBUG_TRIA_EXECUTE) cout<<"\taNumFaces = "<<aNumFaces<<"; aCellLength = "<<aCellLength<<"; aDistEps = "<<aDistEps<<"\n";
 234
 235   // To initialize set of points that belong to the cell
 236   typedef std::set<vtkIdType> TPointIds;
 237   TPointIds anInitialPointIds;
 238   for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++)
 239     anInitialPointIds.insert(GetPointId(aPntId));
 240
 241   // To initialize set of points by face that belong to the cell and backward
 242   typedef std::set<vtkIdType> TFace2Visibility;
 243   TFace2Visibility aFace2Visibility;
 244
 245   typedef std::set<TPointIds> TFace2PointIds;
 246   TFace2PointIds aFace2PointIds;
 247
 248   for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) {
 249     vtkCell* aFace = GetFace(aFaceId);
 250
 251     GetCellNeighbors(theInput, theCellId, aFace, myCellIds);
 252     bool process = myCellIds->GetNumberOfIds() <= 0 ? true : theAppendCoincident3D;
 253     if((!theAllVisible && !theCellsVisibility[myCellIds->GetId(0)]) ||
 254        myCellIds->GetNumberOfIds() <= 0 || theShowInside || process)
 255     {
 256       TPointIds aPointIds;
 257       vtkIdList *anIdList = aFace->PointIds;
 258       aPointIds.insert(anIdList->GetId(0));
 259       aPointIds.insert(anIdList->GetId(1));
 260       aPointIds.insert(anIdList->GetId(2));
 261
 262       aFace2PointIds.insert(aPointIds);
 263       aFace2Visibility.insert(aFaceId);
 264     }
 265   }
 266
 267
 268   ::TPolygons aPolygons;
 269
 270   for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) {
 271     if(aFace2Visibility.find(aFaceId) == aFace2Visibility.end())
 272       continue;
 273
 274     vtkCell* aFace = GetFace(aFaceId);
 275
 276     vtkIdList *anIdList = aFace->PointIds;
 277     vtkIdType aNewPts[3] = {anIdList->GetId(0), anIdList->GetId(1), anIdList->GetId(2)};
 278
 279     // To initialize set of points for the plane where the trinangle face belong to
 280     TPointIds aPointIds;
 281     aPointIds.insert(aNewPts[0]);
 282     aPointIds.insert(aNewPts[1]);
 283     aPointIds.insert(aNewPts[2]);
 284
 285     // To get know, if the points of the trinagle were already observed
 286     bool anIsObserved = aFace2PointIds.find(aPointIds) == aFace2PointIds.end();
 287     if(DEBUG_TRIA_EXECUTE) {
 288       cout<<"\taFaceId = "<<aFaceId<<"; anIsObserved = "<<anIsObserved;
 289       cout<<"; aNewPts = {"<<aNewPts[0]<<", "<<aNewPts[1]<<", "<<aNewPts[2]<<"}\n";
 290     }
 291
 292     if(!anIsObserved){
 293       // To get coordinates of the points of the traingle face
 294       vtkFloatingPointType aCoord[3][3];
 295       aPoints->GetPoint(aNewPts[0],aCoord[0]);
 296       aPoints->GetPoint(aNewPts[1],aCoord[1]);
 297       aPoints->GetPoint(aNewPts[2],aCoord[2]);
 298
 299       /* To calculate plane normal for face (aFace)
 300
 301
 302         ^ aNormal
 303         |
 304         |   ^ aVector01
 305         | /
 306         /_________> aVector02
 307
 308
 309       */
 310       vtkFloatingPointType aVector01[3] = { aCoord[1][0] - aCoord[0][0],
 311                                             aCoord[1][1] - aCoord[0][1],
 312                                             aCoord[1][2] - aCoord[0][2] };
 313
 314       vtkFloatingPointType aVector02[3] = { aCoord[2][0] - aCoord[0][0],
 315                                             aCoord[2][1] - aCoord[0][1],
 316                                             aCoord[2][2] - aCoord[0][2] };
 317
 318       vtkMath::Normalize(aVector01);
 319       vtkMath::Normalize(aVector02);
 320
 321       // To calculate the normal for the triangle
 322       vtkFloatingPointType aNormal[3];
 323       vtkMath::Cross(aVector02,aVector01,aNormal);
 324
 325       vtkMath::Normalize(aNormal);
 326
 327       // To calculate what points belong to the plane
 328       // To calculate bounds of the point set
 329       vtkFloatingPointType aCenter[3] = {0.0, 0.0, 0.0};
 330       {
 331         TPointIds::const_iterator anIter = anInitialPointIds.begin();
 332         TPointIds::const_iterator anEndIter = anInitialPointIds.end();
 333         for(; anIter != anEndIter; anIter++){
 334           vtkFloatingPointType aPntCoord[3];
 335           vtkIdType aPntId = *anIter;
 336           aPoints->GetPoint(aPntId,aPntCoord);
 337
 338           vtkFloatingPointType aVector0Pnt[3] = { aPntCoord[0] - aCoord[0][0],
 339                                                   aPntCoord[1] - aCoord[0][1],
 340                                                   aPntCoord[2] - aCoord[0][2] };
 341
 342
 343           vtkMath::Normalize(aVector0Pnt);
 344
 345           vtkFloatingPointType aNormalPnt[3];
 346           // calculate aNormalPnt
 347           {
 348             vtkFloatingPointType aCosPnt01 = vtkMath::Dot(aVector0Pnt,aVector01);
 349             vtkFloatingPointType aCosPnt02 = vtkMath::Dot(aVector0Pnt,aVector02);
 350             if(aCosPnt01<-1)
 351               aCosPnt01 = -1;
 352             if(aCosPnt01>1)
 353               aCosPnt01 = 1;
 354             if(aCosPnt02<-1)
 355               aCosPnt02 = -1;
 356             if(aCosPnt02>1)
 357               aCosPnt02 = 1;
 358
 359             vtkFloatingPointType aDist01,aDist02;// deflection from Pi/3 angle (equilateral triangle)
 360             vtkFloatingPointType aAngPnt01 = fabs(acos(aCosPnt01));
 361             vtkFloatingPointType aAngPnt02 = fabs(acos(aCosPnt02));
 362
 363             /*  check that triangle similar to equilateral triangle
 364                 AOC or COB ?
 365                 aVector0Pnt = (OC)
 366                 aVector01   = (OB)
 367                 aVector02   = (OA)
 368
 369             B
 370             ^ aVector01  C
 371             |           ^ aVector0Pnt
 372             |     _____/
 373             | ___/
 374             |/________> aVector02
 375             O          A
 376             */
 377             aDist01 = fabs(aAngPnt01-(vtkMath::Pi())/3.0);
 378             aDist02 = fabs(aAngPnt02-(vtkMath::Pi())/3.0);
 379
 380             // caculate a normal for best triangle
 381             if(aDist01 <= aDist02)
 382               vtkMath::Cross(aVector0Pnt,aVector01,aNormalPnt);
 383             else
 384               vtkMath::Cross(aVector0Pnt,aVector02,aNormalPnt);
 385
 386           }
 387
 388           vtkMath::Normalize(aNormalPnt);
 389
 390           if(DEBUG_TRIA_EXECUTE)
 391             cout<<"\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"};";
 392
 393           vtkFloatingPointType aDist = vtkPlane::DistanceToPlane(aPntCoord,aNormal,aCoord[0]);
 394           if(DEBUG_TRIA_EXECUTE) cout<<": aDist = "<<aDist;
 395           if(fabs(aDist) < aDistEps){
 396             aPointIds.insert(aPntId);
 397             aCenter[0] += aPntCoord[0];
 398             aCenter[1] += aPntCoord[1];
 399             aCenter[2] += aPntCoord[2];
 400             if(DEBUG_TRIA_EXECUTE) cout  << "; Added = TRUE" << endl;
 401           } else {
 402             if(DEBUG_TRIA_EXECUTE) cout  << "; Added = FALSE" << endl;
 403           }
 404         }
 405         int aNbPoints = aPointIds.size();
 406         aCenter[0] /= aNbPoints;
 407         aCenter[1] /= aNbPoints;
 408         aCenter[2] /= aNbPoints;
 409       }
 410
 411       //To sinchronize orientation of the cell and its face
 412       vtkFloatingPointType aVectorC[3] = { aCenter[0] - aCellCenter[0],
 413                                            aCenter[1] - aCellCenter[1],
 414                                            aCenter[2] - aCellCenter[2] };
 415       vtkMath::Normalize(aVectorC);
 416
 417       vtkFloatingPointType aDot = vtkMath::Dot(aNormal,aVectorC);
 418       if(DEBUG_TRIA_EXECUTE) {
 419         cout<<"\t\taNormal = {"<<aNormal[0]<<", "<<aNormal[1]<<", "<<aNormal[2]<<"}";
 420         cout<<"; aVectorC = {"<<aVectorC[0]<<", "<<aVectorC[1]<<", "<<aVectorC[2]<<"}\n";
 421         cout<<"\t\taDot = "<<aDot<<"\n";
 422       }
 423       if(aDot > 0){
 424         aNormal[0] = -aNormal[0];
 425         aNormal[1] = -aNormal[1];
 426         aNormal[2] = -aNormal[2];
 427       }
 428
 429       // To calculate the primary direction for point set
 430       vtkFloatingPointType aVector0[3] = { aCoord[0][0] - aCenter[0],
 431                                            aCoord[0][1] - aCenter[1],
 432                                            aCoord[0][2] - aCenter[2] };
 433       vtkMath::Normalize(aVector0);
 434
 435       if(DEBUG_TRIA_EXECUTE) {
 436         cout<<"\t\taCenter = {"<<aCenter[0]<<", "<<aCenter[1]<<", "<<aCenter[2]<<"}";
 437         cout<<"; aVector0 = {"<<aVector0[0]<<", "<<aVector0[1]<<", "<<aVector0[2]<<"}\n";
 438       }
 439
 440       // To calculate the set of points by face those that belong to the plane
 441       TFace2PointIds aRemoveFace2PointIds;
 442       {
 443         TFace2PointIds::const_iterator anIter = aFace2PointIds.begin();
 444         TFace2PointIds::const_iterator anEndIter = aFace2PointIds.end();
 445         for(; anIter != anEndIter; anIter++){
 446           const TPointIds& anIds = *anIter;
 447           TPointIds anIntersection;
 448           std::set_intersection(aPointIds.begin(),aPointIds.end(),
 449                                 anIds.begin(),anIds.end(),
 450                                 std::inserter(anIntersection,anIntersection.begin()));
 451
 452
 453           if(DEBUG_TRIA_EXECUTE) {
 454             cout << "anIntersection:";
 455             TPointIds::iterator aII = anIntersection.begin();
 456             for(;aII!=anIntersection.end();aII++)
 457               cout << *aII << ",";
 458             cout << endl;
 459             cout << "anIds         :";
 460             TPointIds::const_iterator aIIds = anIds.begin();
 461             for(;aIIds!=anIds.end();aIIds++)
 462               cout << *aIIds << ",";
 463             cout << endl;
 464           }
 465           if(anIntersection == anIds){
 466             aRemoveFace2PointIds.insert(anIds);
 467           }
 468         }
 469       }
 470
 471       // To remove from the set of points by face those that belong to the plane
 472       {
 473         TFace2PointIds::const_iterator anIter = aRemoveFace2PointIds.begin();
 474         TFace2PointIds::const_iterator anEndIter = aRemoveFace2PointIds.end();
 475         for(; anIter != anEndIter; anIter++){
 476           const TPointIds& anIds = *anIter;
 477           aFace2PointIds.erase(anIds);
 478         }
 479       }
 480
 481       // To sort the planar set of the points accrding to the angle
 482       {
 483         typedef std::map<vtkFloatingPointType,vtkIdType> TSortedPointIds;
 484         TSortedPointIds aSortedPointIds;
 485
 486         TPointIds::const_iterator anIter = aPointIds.begin();
 487         TPointIds::const_iterator anEndIter = aPointIds.end();
 488         for(; anIter != anEndIter; anIter++){
 489           vtkFloatingPointType aPntCoord[3];
 490           vtkIdType aPntId = *anIter;
 491           aPoints->GetPoint(aPntId,aPntCoord);
 492           vtkFloatingPointType aVector[3] = { aPntCoord[0] - aCenter[0],
 493                                               aPntCoord[1] - aCenter[1],
 494                                               aPntCoord[2] - aCenter[2] };
 495           vtkMath::Normalize(aVector);
 496
 497           vtkFloatingPointType aCross[3];
 498           vtkMath::Cross(aVector,aVector0,aCross);
 499           vtkFloatingPointType aCr = vtkMath::Dot(aCross,aNormal);
 500           bool aGreaterThanPi = aCr < 0;
 501           vtkFloatingPointType aCosinus = vtkMath::Dot(aVector,aVector0);
 502           vtkFloatingPointType anAngle = 0.0;
 503           if(aCosinus >= 1.0){
 504             aCosinus = 1.0;
 505           } else if (aCosinus <= -1.0){
 506             aCosinus = -1.0;
 507             anAngle = vtkMath::Pi();
 508           } else {
 509             anAngle = acos(aCosinus);
 510             if(aGreaterThanPi)
 511               anAngle = 2*vtkMath::Pi() - anAngle;
 512           }
 513
 514           if(DEBUG_TRIA_EXECUTE) {
 515             cout << "\t\t\t vtkMath::Dot(aCross,aNormal)="<<aCr<<endl;
 516             cout<<"\t\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}";
 517             cout<<"; aGreaterThanPi = "<<aGreaterThanPi<<"; aCosinus = "<<aCosinus<<"; anAngle = "<<anAngle<<"\n";
 518           }
 519           aSortedPointIds[anAngle] = aPntId;
 520         }
 521
 522         if(!aSortedPointIds.empty()){
 523           int aNumFacePts = aSortedPointIds.size();
 524           ::TConnectivities aConnectivities(aNumFacePts);
 525           TSortedPointIds::const_iterator anIter = aSortedPointIds.begin();
 526           TSortedPointIds::const_iterator anEndIter = aSortedPointIds.end();
 527           if(DEBUG_TRIA_EXECUTE) cout << "Polygon:";
 528           for(vtkIdType anId = 0; anIter != anEndIter; anIter++, anId++){
 529             vtkIdType aPntId = anIter->second;
 530             aConnectivities[anId] = GetConnectivity(aPntId);
 531             if(DEBUG_TRIA_EXECUTE) cout << aPntId << ",";
 532           }
 533           if(DEBUG_TRIA_EXECUTE) cout << endl;
 534           aPolygons.push_back(::TPolygon(aConnectivities,aCenter,aNormal));
 535         }
 536       }
 537     }
 538   }
 539   if(aPolygons.empty())
 540     return true;
 541
 542   // To check, whether the polygons give a convex polyhedron or not
 543   if(theIsCheckConvex){
 544     int aNbPolygons = aPolygons.size();
 545     for (int aPolygonId = 0; aPolygonId < aNbPolygons; aPolygonId++) {
 546       ::TPolygon& aPolygon = aPolygons[aPolygonId];
 547       vtkFloatingPointType* aNormal = aPolygon.myNormal;
 548       vtkFloatingPointType* anOrigin = aPolygon.myOrigin;
 549       if(DEBUG_TRIA_EXECUTE) {
 550         cout<<"\taPolygonId = "<<aPolygonId<<"\n";
 551         cout<<"\t\taNormal = {"<<aNormal[0]<<", "<<aNormal[1]<<", "<<aNormal[2]<<"}";
 552         cout<<"; anOrigin = {"<<anOrigin[0]<<", "<<anOrigin[1]<<", "<<anOrigin[2]<<"}\n";
 553       }
 554       for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++){
 555         vtkFloatingPointType aPntCoord[3];
 556         vtkIdType anId = GetPointId(aPntId);
 557         aPoints->GetPoint(anId,aPntCoord);
 558         vtkFloatingPointType aDist = vtkPlane::Evaluate(aNormal,anOrigin,aPntCoord);
 559         if(DEBUG_TRIA_EXECUTE) cout<<"\t\taPntId = "<<anId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}; aDist = "<<aDist<<"\n";
 560         if(aDist < -aDistEps)
 561           return false;
 562       }
 563     }
 564   }
 565
 566
 567   // To pass resulting set of the polygons to the output
 568   {
 569     int aNbPolygons = aPolygons.size();
 570     for (int aPolygonId = 0; aPolygonId < aNbPolygons; aPolygonId++) {
 571       ::TPolygon& aPolygon = aPolygons[aPolygonId];
 572       if(DEBUG_TRIA_EXECUTE) cout << "PoilygonId="<<aPolygonId<<" | ";
 573       TConnectivities& aConnectivities = aPolygon.myConnectivities;
 574       if(DEBUG_TRIA_EXECUTE) {
 575         for(int i=0;i<aConnectivities.size();i++)
 576           cout << aConnectivities[i] << ",";
 577         cout << endl;
 578       }
 579       int aNbPoints = aConnectivities.size();
 580       vtkIdType aNewCellId = theOutput->InsertNextCell(VTK_POLYGON,aNbPoints,&aConnectivities[0]);
 581       if(theStoreMapping)
 582         VTKViewer_GeometryFilter::InsertId( theCellId, VTK_POLYGON, theVTK2ObjIds, theDimension2VTK2ObjIds );
 583       theOutputCD->CopyData(thInputCD,theCellId,aNewCellId);
 584     }
 585   }
 586
 587   if(DEBUG_TRIA_EXECUTE) cout<<"\tTriangulator - Ok\n";
 588
 589   return true;
 590 }
 591
 592
 593 //----------------------------------------------------------------------------
 594 VTKViewer_OrderedTriangulator
 595 ::VTKViewer_OrderedTriangulator():
 596   myTriangulator(vtkOrderedTriangulator::New()),
 597   myBoundaryTris(vtkCellArray::New()),
 598   myTriangle(vtkTriangle::New())
 599 {
 600   myBoundaryTris->Allocate(VTK_CELL_SIZE);
 601   myTriangulator->PreSortedOff();
 602 }
 603
 604
 605 //----------------------------------------------------------------------------
 606 VTKViewer_OrderedTriangulator
 607 ::~VTKViewer_OrderedTriangulator()
 608 {
 609   myTriangle->Delete();
 610   myBoundaryTris->Delete();
 611   myTriangulator->Delete();
 612 }
 613
 614
 615 //----------------------------------------------------------------------------
 616 vtkPoints*
 617 VTKViewer_OrderedTriangulator
 618 ::InitPoints(vtkUnstructuredGrid *theInput,
 619              vtkIdType theCellId)
 620 {
 621   myBoundaryTris->Reset();
 622
 623   vtkPoints* aPoints = VTKViewer_Triangulator::InitPoints(theInput, theCellId);
 624   vtkIdType aNumPts = myPoints->GetNumberOfPoints();
 625   if ( aNumPts > 0 ) {
 626     myTriangulator->InitTriangulation(0.0, 1.0, 0.0, 1.0, 0.0, 1.0, aNumPts);
 627
 628     vtkFloatingPointType aBounds[6];
 629     myPoints->GetBounds(aBounds);
 630     vtkFloatingPointType xSize, ySize, zSize;
 631     xSize = aBounds[1] - aBounds[0];
 632     ySize = aBounds[3] - aBounds[2];
 633     zSize = aBounds[5] - aBounds[4];
 634     vtkFloatingPointType anAbsoluteCoord[3];
 635     vtkFloatingPointType aParamentrucCoord[3];
 636     for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
 637       myPoints->GetPoint(aPntId, anAbsoluteCoord);
 638       aParamentrucCoord[0] = xSize==0. ? 0. : ((anAbsoluteCoord[0] - aBounds[0]) / xSize);
 639       aParamentrucCoord[1] = ySize==0. ? 0. : ((anAbsoluteCoord[1] - aBounds[2]) / ySize);
 640       aParamentrucCoord[2] = zSize==0. ? 0. : ((anAbsoluteCoord[2] - aBounds[4]) / zSize);
 641       myTriangulator->InsertPoint(aPntId, anAbsoluteCoord, aParamentrucCoord, 0);
 642     }
 643
 644     myTriangulator->Triangulate();
 645     myTriangulator->AddTriangles(myBoundaryTris);
 646   }
 647
 648   return aPoints;
 649 }
 650
 651
 652 //----------------------------------------------------------------------------
 653 vtkIdType
 654 VTKViewer_OrderedTriangulator
 655 ::GetNumFaces()
 656 {
 657   return myBoundaryTris->GetNumberOfCells();
 658 }
 659
 660
 661 //----------------------------------------------------------------------------
 662 vtkCell*
 663 VTKViewer_OrderedTriangulator
 664 ::GetFace(vtkIdType theFaceId)
 665 {
 666   vtkIdType aNumCells = myBoundaryTris->GetNumberOfCells();
 667   if ( theFaceId < 0 || theFaceId >= aNumCells )
 668     return NULL;
 669
 670   vtkIdType *aCells = myBoundaryTris->GetPointer();
 671
 672   // Each triangle has three points plus number of points
 673   vtkIdType *aCellPtr = aCells + 4*theFaceId;
 674
 675   myTriangle->PointIds->SetId(0, aCellPtr[1]);
 676   myTriangle->Points->SetPoint(0, myPoints->GetPoint(aCellPtr[1]));
 677
 678   myTriangle->PointIds->SetId(1, aCellPtr[2]);
 679   myTriangle->Points->SetPoint(1, myPoints->GetPoint(aCellPtr[2]));
 680
 681   myTriangle->PointIds->SetId(2, aCellPtr[3]);
 682   myTriangle->Points->SetPoint(2, myPoints->GetPoint(aCellPtr[3]));
 683
 684   return myTriangle;
 685 }
 686
 687
 688 //----------------------------------------------------------------------------
 689 VTKViewer_DelaunayTriangulator
 690 ::VTKViewer_DelaunayTriangulator():
 691   myUnstructuredGrid(vtkUnstructuredGrid::New()),
 692   myGeometryFilter(vtkGeometryFilter::New()),
 693   myDelaunay3D(vtkDelaunay3D::New()),
 694   myPolyData(NULL)
 695 {
 696   myUnstructuredGrid->Initialize();
 697   myUnstructuredGrid->Allocate();
 698   myUnstructuredGrid->SetPoints(myPoints);
 699
 700   myDelaunay3D->SetInput(myUnstructuredGrid);
 701   myGeometryFilter->SetInput(myDelaunay3D->GetOutput());
 702   myPolyData = myGeometryFilter->GetOutput();
 703 }
 704
 705
 706 //----------------------------------------------------------------------------
 707 VTKViewer_DelaunayTriangulator
 708 ::~VTKViewer_DelaunayTriangulator()
 709 {
 710   myUnstructuredGrid->Delete();
 711   myGeometryFilter->Delete();
 712   myDelaunay3D->Delete();
 713 }
 714
 715
 716 //----------------------------------------------------------------------------
 717 vtkPoints*
 718 VTKViewer_DelaunayTriangulator
 719 ::InitPoints(vtkUnstructuredGrid *theInput,
 720              vtkIdType theCellId)
 721 {
 722   vtkPoints* aPoints = VTKViewer_Triangulator::InitPoints(theInput, theCellId);
 723
 724   myPoints->Modified();
 725   myUnstructuredGrid->Modified();
 726   myGeometryFilter->Update();
 727
 728   return aPoints;
 729 }
 730
 731
 732 //----------------------------------------------------------------------------
 733 vtkIdType
 734 VTKViewer_DelaunayTriangulator
 735 ::GetNumFaces()
 736 {
 737   return myPolyData->GetNumberOfCells();
 738 }
 739
 740
 741 //----------------------------------------------------------------------------
 742 vtkCell*
 743 VTKViewer_DelaunayTriangulator
 744 ::GetFace(vtkIdType theFaceId)
 745 {
 746   return myPolyData->GetCell(theFaceId);
 747 }