src/VTKViewer/VTKViewer_ConvexTool.cxx

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