src/VTKViewer/VTKViewer_ConvexTool.cxx

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