src/VTKViewer/VTKViewer_ConvexTool.cxx

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