Back to index

salome-gui  6.5.0
VTKViewer_ConvexTool.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 #include "VTKViewer_ConvexTool.h"
00024 #include "VTKViewer_GeometryFilter.h"
00025 
00026 #include <set>
00027 #include <map>
00028 #include <algorithm>
00029 
00030 #include <vtkUnstructuredGrid.h>
00031 #include <vtkGeometryFilter.h>
00032 #include <vtkDelaunay3D.h>
00033 #include <vtkGenericCell.h>
00034 #include <vtkPolyData.h>
00035 #include <vtkCellData.h>
00036 #include <vtkPoints.h>
00037 #include <vtkIdList.h>
00038 #include <vtkCell.h>
00039 #include <vtkPlane.h>
00040 #include <vtkMath.h>
00041 #include <vtkCellArray.h>
00042 #include <vtkTriangle.h>
00043 #include <vtkOrderedTriangulator.h>
00044 
00045 #ifdef _DEBUG_
00046 static int DEBUG_TRIA_EXECUTE = 0;
00047 #else
00048 static int DEBUG_TRIA_EXECUTE = 0;
00049 #endif
00050 
00051 namespace
00052 {
00053   typedef std::vector<vtkIdType> TConnectivities;
00054 
00055   struct TPolygon
00056   {
00057     TConnectivities myConnectivities;
00058     vtkFloatingPointType myOrigin[3];
00059     vtkFloatingPointType myNormal[3];
00060     TPolygon(const TConnectivities& theConnectivities,
00061              vtkFloatingPointType theOrigin[3],
00062              vtkFloatingPointType theNormal[3]):
00063       myConnectivities(theConnectivities)
00064     {
00065       myOrigin[0] = theOrigin[0];
00066       myOrigin[1] = theOrigin[1];
00067       myOrigin[2] = theOrigin[2];
00068 
00069       myNormal[0] = theNormal[0];
00070       myNormal[1] = theNormal[1];
00071       myNormal[2] = theNormal[2];
00072     }
00073   };
00074 
00075   typedef std::vector<TPolygon> TPolygons;
00076 }
00077 
00078 
00079 //----------------------------------------------------------------------------
00080 VTKViewer_Triangulator
00081 ::VTKViewer_Triangulator():
00082   myCellIds(vtkIdList::New()),
00083   myFaceIds(vtkIdList::New()),
00084   myPoints(vtkPoints::New()),
00085   myPointIds(NULL)
00086 {}
00087 
00088 
00089 //----------------------------------------------------------------------------
00090 VTKViewer_Triangulator
00091 ::~VTKViewer_Triangulator()
00092 {
00093   myCellIds->Delete();
00094   myFaceIds->Delete();
00095   myPoints->Delete();
00096 }
00097 
00098 
00099 //----------------------------------------------------------------------------
00100 vtkPoints*
00101 VTKViewer_Triangulator
00102 ::InitPoints(vtkUnstructuredGrid *theInput,
00103              vtkIdType theCellId)
00104 {
00105   myPoints->Reset();
00106   myPoints->Modified(); // the VTK bug
00107 
00108   vtkIdType aNumPts;
00109   theInput->GetCellPoints(theCellId, aNumPts, myPointIds); 
00110   if ( aNumPts > 0 ) {
00111     vtkFloatingPointType anAbsoluteCoord[3];
00112     myPoints->SetNumberOfPoints(aNumPts);
00113     vtkPoints *anInputPoints = theInput->GetPoints();
00114     for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
00115       anInputPoints->GetPoint(myPointIds[aPntId], anAbsoluteCoord);
00116       myPoints->SetPoint(aPntId, anAbsoluteCoord);
00117     }
00118   }
00119 
00120   return myPoints;
00121 }
00122 
00123 
00124 //----------------------------------------------------------------------------
00125 vtkIdType 
00126 VTKViewer_Triangulator
00127 ::GetNbOfPoints()
00128 {
00129   return myPoints->GetNumberOfPoints();
00130 }
00131 
00132 
00133 //----------------------------------------------------------------------------
00134 vtkIdType 
00135 VTKViewer_Triangulator
00136 ::GetPointId(vtkIdType thePointId)
00137 {
00138   return thePointId;
00139 }
00140 
00141 
00142 //----------------------------------------------------------------------------
00143 vtkFloatingPointType 
00144 VTKViewer_Triangulator
00145 ::GetCellLength()
00146 {
00147   vtkFloatingPointType aBounds[6];
00148   myPoints->GetBounds(aBounds);
00149 
00150   vtkFloatingPointType aCoordDiff[3];
00151   aCoordDiff[0] = (aBounds[1] - aBounds[0]);
00152   aCoordDiff[1] = (aBounds[3] - aBounds[2]);
00153   aCoordDiff[2] = (aBounds[5] - aBounds[4]);
00154 
00155   return sqrt(aCoordDiff[0]*aCoordDiff[0] + 
00156               aCoordDiff[1]*aCoordDiff[1] + 
00157               aCoordDiff[2]*aCoordDiff[2]);
00158 }
00159 
00160 
00161 //----------------------------------------------------------------------------
00162 void 
00163 VTKViewer_Triangulator
00164 ::GetCellNeighbors(vtkUnstructuredGrid *theInput,
00165                    vtkIdType theCellId,
00166                    vtkCell* theFace,
00167                    vtkIdList* theCellIds)
00168 {
00169   myFaceIds->Reset();
00170   vtkIdList *anIdList = theFace->PointIds;  
00171   myFaceIds->InsertNextId(myPointIds[anIdList->GetId(0)]);
00172   myFaceIds->InsertNextId(myPointIds[anIdList->GetId(1)]);
00173   myFaceIds->InsertNextId(myPointIds[anIdList->GetId(2)]);
00174 
00175   theInput->GetCellNeighbors(theCellId, myFaceIds, theCellIds);
00176 }
00177 
00178 
00179 //----------------------------------------------------------------------------
00180 vtkIdType 
00181 VTKViewer_Triangulator
00182 ::GetConnectivity(vtkIdType thePntId)
00183 {
00184   return myPointIds[thePntId];
00185 }
00186 
00187 
00188 //----------------------------------------------------------------------------
00189 bool 
00190 VTKViewer_Triangulator
00191 ::Execute(vtkUnstructuredGrid *theInput,
00192           vtkCellData* thInputCD,
00193           vtkIdType theCellId,
00194           int theShowInside,
00195           int theAllVisible,
00196           int theAppendCoincident3D,
00197           const char* theCellsVisibility,
00198           vtkPolyData *theOutput,
00199           vtkCellData* theOutputCD,
00200           int theStoreMapping,
00201           std::vector<vtkIdType>& theVTK2ObjIds,
00202           std::map< vtkIdType, std::vector<vtkIdType> >& theDimension2VTK2ObjIds,
00203           bool theIsCheckConvex)
00204 {
00205   vtkPoints *aPoints = InitPoints(theInput, theCellId);
00206   vtkIdType aNumPts = GetNbOfPoints();
00207   if(DEBUG_TRIA_EXECUTE) cout<<"Triangulator - aNumPts = "<<aNumPts<<"\n";
00208 
00209   if(aNumPts == 0)
00210     return true;
00211 
00212   // To calculate the bary center of the cell
00213   vtkFloatingPointType aCellCenter[3] = {0.0, 0.0, 0.0};
00214   {
00215     vtkFloatingPointType aPntCoord[3];
00216     for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
00217       aPoints->GetPoint(GetPointId(aPntId),aPntCoord);
00218       if(DEBUG_TRIA_EXECUTE) cout<<"\taPntId = "<<GetPointId(aPntId)<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}\n";
00219       aCellCenter[0] += aPntCoord[0];
00220       aCellCenter[1] += aPntCoord[1];
00221       aCellCenter[2] += aPntCoord[2];
00222     }
00223     aCellCenter[0] /= aNumPts;
00224     aCellCenter[1] /= aNumPts;
00225     aCellCenter[2] /= aNumPts;
00226   }
00227 
00228   vtkFloatingPointType aCellLength = GetCellLength();
00229   int aNumFaces = GetNumFaces();
00230 
00231   static vtkFloatingPointType EPS = 1.0E-2;
00232   vtkFloatingPointType aDistEps = aCellLength/3.0 * EPS;
00233   if(DEBUG_TRIA_EXECUTE) cout<<"\taNumFaces = "<<aNumFaces<<"; aCellLength = "<<aCellLength<<"; aDistEps = "<<aDistEps<<"\n";
00234 
00235   // To initialize set of points that belong to the cell
00236   typedef std::set<vtkIdType> TPointIds;
00237   TPointIds anInitialPointIds;
00238   for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++)
00239     anInitialPointIds.insert(GetPointId(aPntId));
00240   
00241   // To initialize set of points by face that belong to the cell and backward
00242   typedef std::set<vtkIdType> TFace2Visibility;
00243   TFace2Visibility aFace2Visibility;
00244   
00245   typedef std::set<TPointIds> TFace2PointIds;
00246   TFace2PointIds aFace2PointIds;
00247 
00248   for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) {
00249     vtkCell* aFace = GetFace(aFaceId);
00250     
00251     GetCellNeighbors(theInput, theCellId, aFace, myCellIds);
00252     bool process = myCellIds->GetNumberOfIds() <= 0 ? true : theAppendCoincident3D;
00253     if((!theAllVisible && !theCellsVisibility[myCellIds->GetId(0)]) || 
00254        myCellIds->GetNumberOfIds() <= 0 || theShowInside || process)
00255     {
00256       TPointIds aPointIds;
00257       vtkIdList *anIdList = aFace->PointIds;  
00258       aPointIds.insert(anIdList->GetId(0));
00259       aPointIds.insert(anIdList->GetId(1));
00260       aPointIds.insert(anIdList->GetId(2));
00261       
00262       aFace2PointIds.insert(aPointIds);
00263       aFace2Visibility.insert(aFaceId);
00264     }
00265   }
00266 
00267 
00268   ::TPolygons aPolygons;
00269 
00270   for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) {
00271     if(aFace2Visibility.find(aFaceId) == aFace2Visibility.end())
00272       continue;
00273 
00274     vtkCell* aFace = GetFace(aFaceId);
00275 
00276     vtkIdList *anIdList = aFace->PointIds;
00277     vtkIdType aNewPts[3] = {anIdList->GetId(0), anIdList->GetId(1), anIdList->GetId(2)};
00278             
00279     // To initialize set of points for the plane where the trinangle face belong to
00280     TPointIds aPointIds;
00281     aPointIds.insert(aNewPts[0]);
00282     aPointIds.insert(aNewPts[1]);
00283     aPointIds.insert(aNewPts[2]);
00284 
00285     // To get know, if the points of the trinagle were already observed
00286     bool anIsObserved = aFace2PointIds.find(aPointIds) == aFace2PointIds.end();
00287     if(DEBUG_TRIA_EXECUTE) {
00288       cout<<"\taFaceId = "<<aFaceId<<"; anIsObserved = "<<anIsObserved;
00289       cout<<"; aNewPts = {"<<aNewPts[0]<<", "<<aNewPts[1]<<", "<<aNewPts[2]<<"}\n";
00290     }
00291     
00292     if(!anIsObserved){
00293       // To get coordinates of the points of the traingle face
00294       vtkFloatingPointType aCoord[3][3];
00295       aPoints->GetPoint(aNewPts[0],aCoord[0]);
00296       aPoints->GetPoint(aNewPts[1],aCoord[1]);
00297       aPoints->GetPoint(aNewPts[2],aCoord[2]);
00298       
00299       /* To calculate plane normal for face (aFace)
00300 
00301 
00302         ^ aNormal
00303         |     
00304         |   ^ aVector01
00305         | /
00306         /_________> aVector02
00307        
00308       
00309       */
00310       vtkFloatingPointType aVector01[3] = { aCoord[1][0] - aCoord[0][0],
00311                                             aCoord[1][1] - aCoord[0][1],
00312                                             aCoord[1][2] - aCoord[0][2] };
00313       
00314       vtkFloatingPointType aVector02[3] = { aCoord[2][0] - aCoord[0][0],
00315                                             aCoord[2][1] - aCoord[0][1],
00316                                             aCoord[2][2] - aCoord[0][2] };
00317       
00318       vtkMath::Normalize(aVector01);
00319       vtkMath::Normalize(aVector02);
00320 
00321       // To calculate the normal for the triangle
00322       vtkFloatingPointType aNormal[3];
00323       vtkMath::Cross(aVector02,aVector01,aNormal);
00324       
00325       vtkMath::Normalize(aNormal);
00326       
00327       // To calculate what points belong to the plane
00328       // To calculate bounds of the point set
00329       vtkFloatingPointType aCenter[3] = {0.0, 0.0, 0.0};
00330       {
00331         TPointIds::const_iterator anIter = anInitialPointIds.begin();
00332         TPointIds::const_iterator anEndIter = anInitialPointIds.end();
00333         for(; anIter != anEndIter; anIter++){
00334           vtkFloatingPointType aPntCoord[3];
00335           vtkIdType aPntId = *anIter;
00336           aPoints->GetPoint(aPntId,aPntCoord);
00337           
00338           vtkFloatingPointType aVector0Pnt[3] = { aPntCoord[0] - aCoord[0][0],
00339                                                   aPntCoord[1] - aCoord[0][1],
00340                                                   aPntCoord[2] - aCoord[0][2] };
00341 
00342           
00343           vtkMath::Normalize(aVector0Pnt);
00344           
00345           vtkFloatingPointType aNormalPnt[3];
00346           // calculate aNormalPnt
00347           {
00348             vtkFloatingPointType aCosPnt01 = vtkMath::Dot(aVector0Pnt,aVector01);
00349             vtkFloatingPointType aCosPnt02 = vtkMath::Dot(aVector0Pnt,aVector02);
00350             if(aCosPnt01<-1)
00351               aCosPnt01 = -1;
00352             if(aCosPnt01>1)
00353               aCosPnt01 = 1;
00354             if(aCosPnt02<-1)
00355               aCosPnt02 = -1;
00356             if(aCosPnt02>1)
00357               aCosPnt02 = 1;
00358 
00359             vtkFloatingPointType aDist01,aDist02;// deflection from Pi/3 angle (equilateral triangle)
00360             vtkFloatingPointType aAngPnt01 = fabs(acos(aCosPnt01));
00361             vtkFloatingPointType aAngPnt02 = fabs(acos(aCosPnt02));
00362 
00363             /*  check that triangle similar to equilateral triangle
00364                 AOC or COB ?
00365                 aVector0Pnt = (OC)
00366                 aVector01   = (OB)
00367                 aVector02   = (OA)
00368             
00369             B
00370             ^ aVector01  C     
00371             |           ^ aVector0Pnt  
00372             |     _____/ 
00373             | ___/
00374             |/________> aVector02
00375             O          A
00376             */
00377             aDist01 = fabs(aAngPnt01-(vtkMath::Pi())/3.0); 
00378             aDist02 = fabs(aAngPnt02-(vtkMath::Pi())/3.0);
00379             
00380             // caculate a normal for best triangle
00381             if(aDist01 <= aDist02)
00382               vtkMath::Cross(aVector0Pnt,aVector01,aNormalPnt);
00383             else
00384               vtkMath::Cross(aVector0Pnt,aVector02,aNormalPnt);
00385 
00386           }
00387           
00388           vtkMath::Normalize(aNormalPnt);
00389           
00390           if(DEBUG_TRIA_EXECUTE)
00391             cout<<"\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"};";
00392           
00393           vtkFloatingPointType aDist = vtkPlane::DistanceToPlane(aPntCoord,aNormal,aCoord[0]);
00394           if(DEBUG_TRIA_EXECUTE) cout<<": aDist = "<<aDist;
00395           if(fabs(aDist) < aDistEps){
00396             aPointIds.insert(aPntId);
00397             aCenter[0] += aPntCoord[0];
00398             aCenter[1] += aPntCoord[1];
00399             aCenter[2] += aPntCoord[2];
00400             if(DEBUG_TRIA_EXECUTE) cout  << "; Added = TRUE" << endl;
00401           } else {
00402             if(DEBUG_TRIA_EXECUTE) cout  << "; Added = FALSE" << endl;
00403           }
00404         }
00405         int aNbPoints = aPointIds.size();
00406         aCenter[0] /= aNbPoints;
00407         aCenter[1] /= aNbPoints;
00408         aCenter[2] /= aNbPoints;
00409       }
00410       
00411       //To sinchronize orientation of the cell and its face
00412       vtkFloatingPointType aVectorC[3] = { aCenter[0] - aCellCenter[0],
00413                                            aCenter[1] - aCellCenter[1],
00414                                            aCenter[2] - aCellCenter[2] };
00415       vtkMath::Normalize(aVectorC);
00416       
00417       vtkFloatingPointType aDot = vtkMath::Dot(aNormal,aVectorC);
00418       if(DEBUG_TRIA_EXECUTE) {
00419         cout<<"\t\taNormal = {"<<aNormal[0]<<", "<<aNormal[1]<<", "<<aNormal[2]<<"}";
00420         cout<<"; aVectorC = {"<<aVectorC[0]<<", "<<aVectorC[1]<<", "<<aVectorC[2]<<"}\n";
00421         cout<<"\t\taDot = "<<aDot<<"\n";
00422       }
00423       if(aDot > 0){
00424         aNormal[0] = -aNormal[0];
00425         aNormal[1] = -aNormal[1];
00426         aNormal[2] = -aNormal[2];
00427       }
00428       
00429       // To calculate the primary direction for point set
00430       vtkFloatingPointType aVector0[3] = { aCoord[0][0] - aCenter[0],
00431                                            aCoord[0][1] - aCenter[1],
00432                                            aCoord[0][2] - aCenter[2] };
00433       vtkMath::Normalize(aVector0);
00434       
00435       if(DEBUG_TRIA_EXECUTE) {
00436         cout<<"\t\taCenter = {"<<aCenter[0]<<", "<<aCenter[1]<<", "<<aCenter[2]<<"}";
00437         cout<<"; aVector0 = {"<<aVector0[0]<<", "<<aVector0[1]<<", "<<aVector0[2]<<"}\n";
00438       }
00439       
00440       // To calculate the set of points by face those that belong to the plane
00441       TFace2PointIds aRemoveFace2PointIds;
00442       {
00443         TFace2PointIds::const_iterator anIter = aFace2PointIds.begin();
00444         TFace2PointIds::const_iterator anEndIter = aFace2PointIds.end();
00445         for(; anIter != anEndIter; anIter++){
00446           const TPointIds& anIds = *anIter;
00447           TPointIds anIntersection;
00448           std::set_intersection(aPointIds.begin(),aPointIds.end(),
00449                                 anIds.begin(),anIds.end(),
00450                                 std::inserter(anIntersection,anIntersection.begin()));
00451           
00452 
00453           if(DEBUG_TRIA_EXECUTE) {
00454             cout << "anIntersection:";
00455             TPointIds::iterator aII = anIntersection.begin();
00456             for(;aII!=anIntersection.end();aII++)
00457               cout << *aII << ",";
00458             cout << endl;
00459             cout << "anIds         :";
00460             TPointIds::const_iterator aIIds = anIds.begin();
00461             for(;aIIds!=anIds.end();aIIds++)
00462               cout << *aIIds << ",";
00463             cout << endl;
00464           }
00465           if(anIntersection == anIds){
00466             aRemoveFace2PointIds.insert(anIds);
00467           }
00468         }
00469       }
00470       
00471       // To remove from the set of points by face those that belong to the plane
00472       {
00473         TFace2PointIds::const_iterator anIter = aRemoveFace2PointIds.begin();
00474         TFace2PointIds::const_iterator anEndIter = aRemoveFace2PointIds.end();
00475         for(; anIter != anEndIter; anIter++){
00476           const TPointIds& anIds = *anIter;
00477           aFace2PointIds.erase(anIds);
00478         }
00479       }
00480       
00481       // To sort the planar set of the points accrding to the angle
00482       {
00483         typedef std::map<vtkFloatingPointType,vtkIdType> TSortedPointIds;
00484         TSortedPointIds aSortedPointIds;
00485         
00486         TPointIds::const_iterator anIter = aPointIds.begin();
00487         TPointIds::const_iterator anEndIter = aPointIds.end();
00488         for(; anIter != anEndIter; anIter++){
00489           vtkFloatingPointType aPntCoord[3];
00490           vtkIdType aPntId = *anIter;
00491           aPoints->GetPoint(aPntId,aPntCoord);
00492           vtkFloatingPointType aVector[3] = { aPntCoord[0] - aCenter[0],
00493                                               aPntCoord[1] - aCenter[1],
00494                                               aPntCoord[2] - aCenter[2] };
00495           vtkMath::Normalize(aVector);
00496           
00497           vtkFloatingPointType aCross[3];
00498           vtkMath::Cross(aVector,aVector0,aCross);
00499           vtkFloatingPointType aCr = vtkMath::Dot(aCross,aNormal);
00500           bool aGreaterThanPi = aCr < 0;
00501           vtkFloatingPointType aCosinus = vtkMath::Dot(aVector,aVector0);
00502           vtkFloatingPointType anAngle = 0.0;
00503           if(aCosinus >= 1.0){
00504             aCosinus = 1.0;
00505           } else if (aCosinus <= -1.0){
00506             aCosinus = -1.0;
00507             anAngle = vtkMath::Pi();
00508           } else {
00509             anAngle = acos(aCosinus);
00510             if(aGreaterThanPi)
00511               anAngle = 2*vtkMath::Pi() - anAngle;
00512           }
00513           
00514           if(DEBUG_TRIA_EXECUTE) {
00515             cout << "\t\t\t vtkMath::Dot(aCross,aNormal)="<<aCr<<endl;
00516             cout<<"\t\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}";
00517             cout<<"; aGreaterThanPi = "<<aGreaterThanPi<<"; aCosinus = "<<aCosinus<<"; anAngle = "<<anAngle<<"\n";
00518           }
00519           aSortedPointIds[anAngle] = aPntId;
00520         }
00521 
00522         if(!aSortedPointIds.empty()){
00523           int aNumFacePts = aSortedPointIds.size();
00524           ::TConnectivities aConnectivities(aNumFacePts);
00525           TSortedPointIds::const_iterator anIter = aSortedPointIds.begin();
00526           TSortedPointIds::const_iterator anEndIter = aSortedPointIds.end();
00527           if(DEBUG_TRIA_EXECUTE) cout << "Polygon:";
00528           for(vtkIdType anId = 0; anIter != anEndIter; anIter++, anId++){
00529             vtkIdType aPntId = anIter->second;
00530             aConnectivities[anId] = GetConnectivity(aPntId);
00531             if(DEBUG_TRIA_EXECUTE) cout << aPntId << ",";
00532           }
00533           if(DEBUG_TRIA_EXECUTE) cout << endl;
00534           aPolygons.push_back(::TPolygon(aConnectivities,aCenter,aNormal));
00535         }
00536       }
00537     }
00538   }
00539   if(aPolygons.empty())
00540     return true;
00541 
00542   // To check, whether the polygons give a convex polyhedron or not
00543   if(theIsCheckConvex){
00544     int aNbPolygons = aPolygons.size();
00545     for (int aPolygonId = 0; aPolygonId < aNbPolygons; aPolygonId++) {
00546       ::TPolygon& aPolygon = aPolygons[aPolygonId];
00547       vtkFloatingPointType* aNormal = aPolygon.myNormal;
00548       vtkFloatingPointType* anOrigin = aPolygon.myOrigin;
00549       if(DEBUG_TRIA_EXECUTE) {
00550         cout<<"\taPolygonId = "<<aPolygonId<<"\n";
00551         cout<<"\t\taNormal = {"<<aNormal[0]<<", "<<aNormal[1]<<", "<<aNormal[2]<<"}";
00552         cout<<"; anOrigin = {"<<anOrigin[0]<<", "<<anOrigin[1]<<", "<<anOrigin[2]<<"}\n";
00553       }
00554       for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++){
00555         vtkFloatingPointType aPntCoord[3];
00556         vtkIdType anId = GetPointId(aPntId);
00557         aPoints->GetPoint(anId,aPntCoord);
00558         vtkFloatingPointType aDist = vtkPlane::Evaluate(aNormal,anOrigin,aPntCoord);
00559         if(DEBUG_TRIA_EXECUTE) cout<<"\t\taPntId = "<<anId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}; aDist = "<<aDist<<"\n";
00560         if(aDist < -aDistEps)
00561           return false;
00562       }
00563     }
00564   }
00565 
00566 
00567   // To pass resulting set of the polygons to the output
00568   {
00569     int aNbPolygons = aPolygons.size();
00570     for (int aPolygonId = 0; aPolygonId < aNbPolygons; aPolygonId++) {
00571       ::TPolygon& aPolygon = aPolygons[aPolygonId];
00572       if(DEBUG_TRIA_EXECUTE) cout << "PoilygonId="<<aPolygonId<<" | ";
00573       TConnectivities& aConnectivities = aPolygon.myConnectivities;
00574       if(DEBUG_TRIA_EXECUTE) {
00575         for(int i=0;i<aConnectivities.size();i++)
00576           cout << aConnectivities[i] << ",";
00577         cout << endl;
00578       }
00579       int aNbPoints = aConnectivities.size();
00580       vtkIdType aNewCellId = theOutput->InsertNextCell(VTK_POLYGON,aNbPoints,&aConnectivities[0]);
00581       if(theStoreMapping)
00582         VTKViewer_GeometryFilter::InsertId( theCellId, VTK_POLYGON, theVTK2ObjIds, theDimension2VTK2ObjIds );
00583       theOutputCD->CopyData(thInputCD,theCellId,aNewCellId);
00584     }
00585   }
00586 
00587   if(DEBUG_TRIA_EXECUTE) cout<<"\tTriangulator - Ok\n";
00588   
00589   return true;
00590 }
00591 
00592 
00593 //----------------------------------------------------------------------------
00594 VTKViewer_OrderedTriangulator
00595 ::VTKViewer_OrderedTriangulator():
00596   myTriangulator(vtkOrderedTriangulator::New()),
00597   myBoundaryTris(vtkCellArray::New()),
00598   myTriangle(vtkTriangle::New())
00599 {
00600   myBoundaryTris->Allocate(VTK_CELL_SIZE);
00601   myTriangulator->PreSortedOff();
00602 }
00603 
00604 
00605 //----------------------------------------------------------------------------
00606 VTKViewer_OrderedTriangulator
00607 ::~VTKViewer_OrderedTriangulator()
00608 {
00609   myTriangle->Delete();
00610   myBoundaryTris->Delete();
00611   myTriangulator->Delete();
00612 }
00613 
00614 
00615 //----------------------------------------------------------------------------
00616 vtkPoints*
00617 VTKViewer_OrderedTriangulator
00618 ::InitPoints(vtkUnstructuredGrid *theInput,
00619              vtkIdType theCellId)
00620 {
00621   myBoundaryTris->Reset();
00622 
00623   vtkPoints* aPoints = VTKViewer_Triangulator::InitPoints(theInput, theCellId);
00624   vtkIdType aNumPts = myPoints->GetNumberOfPoints();
00625   if ( aNumPts > 0 ) {
00626     myTriangulator->InitTriangulation(0.0, 1.0, 0.0, 1.0, 0.0, 1.0, aNumPts);
00627 
00628     vtkFloatingPointType aBounds[6];
00629     myPoints->GetBounds(aBounds);
00630     vtkFloatingPointType xSize, ySize, zSize;
00631     xSize = aBounds[1] - aBounds[0];
00632     ySize = aBounds[3] - aBounds[2];
00633     zSize = aBounds[5] - aBounds[4];
00634     vtkFloatingPointType anAbsoluteCoord[3];
00635     vtkFloatingPointType aParamentrucCoord[3];
00636     for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
00637       myPoints->GetPoint(aPntId, anAbsoluteCoord);
00638       aParamentrucCoord[0] = xSize==0. ? 0. : ((anAbsoluteCoord[0] - aBounds[0]) / xSize);
00639       aParamentrucCoord[1] = ySize==0. ? 0. : ((anAbsoluteCoord[1] - aBounds[2]) / ySize);
00640       aParamentrucCoord[2] = zSize==0. ? 0. : ((anAbsoluteCoord[2] - aBounds[4]) / zSize);
00641       myTriangulator->InsertPoint(aPntId, anAbsoluteCoord, aParamentrucCoord, 0);
00642     }
00643 
00644     myTriangulator->Triangulate();
00645     myTriangulator->AddTriangles(myBoundaryTris);
00646   }
00647 
00648   return aPoints;
00649 }
00650 
00651 
00652 //----------------------------------------------------------------------------
00653 vtkIdType 
00654 VTKViewer_OrderedTriangulator
00655 ::GetNumFaces()
00656 {
00657   return myBoundaryTris->GetNumberOfCells();
00658 }
00659 
00660 
00661 //----------------------------------------------------------------------------
00662 vtkCell*
00663 VTKViewer_OrderedTriangulator
00664 ::GetFace(vtkIdType theFaceId)
00665 {
00666   vtkIdType aNumCells = myBoundaryTris->GetNumberOfCells();
00667   if ( theFaceId < 0 || theFaceId >= aNumCells ) 
00668     return NULL;
00669 
00670   vtkIdType *aCells = myBoundaryTris->GetPointer();
00671 
00672   // Each triangle has three points plus number of points
00673   vtkIdType *aCellPtr = aCells + 4*theFaceId;
00674   
00675   myTriangle->PointIds->SetId(0, aCellPtr[1]);
00676   myTriangle->Points->SetPoint(0, myPoints->GetPoint(aCellPtr[1]));
00677 
00678   myTriangle->PointIds->SetId(1, aCellPtr[2]);
00679   myTriangle->Points->SetPoint(1, myPoints->GetPoint(aCellPtr[2]));
00680 
00681   myTriangle->PointIds->SetId(2, aCellPtr[3]);
00682   myTriangle->Points->SetPoint(2, myPoints->GetPoint(aCellPtr[3]));
00683 
00684   return myTriangle;
00685 }
00686 
00687 
00688 //----------------------------------------------------------------------------
00689 VTKViewer_DelaunayTriangulator
00690 ::VTKViewer_DelaunayTriangulator():
00691   myUnstructuredGrid(vtkUnstructuredGrid::New()),
00692   myGeometryFilter(vtkGeometryFilter::New()),
00693   myDelaunay3D(vtkDelaunay3D::New()),
00694   myPolyData(NULL)
00695 {
00696   myUnstructuredGrid->Initialize();
00697   myUnstructuredGrid->Allocate();
00698   myUnstructuredGrid->SetPoints(myPoints);
00699 
00700   myDelaunay3D->SetInput(myUnstructuredGrid);
00701   myGeometryFilter->SetInput(myDelaunay3D->GetOutput());
00702   myPolyData = myGeometryFilter->GetOutput();
00703 }
00704 
00705 
00706 //----------------------------------------------------------------------------
00707 VTKViewer_DelaunayTriangulator
00708 ::~VTKViewer_DelaunayTriangulator()
00709 {
00710   myUnstructuredGrid->Delete();
00711   myGeometryFilter->Delete();
00712   myDelaunay3D->Delete();
00713 }
00714 
00715 
00716 //----------------------------------------------------------------------------
00717 vtkPoints* 
00718 VTKViewer_DelaunayTriangulator
00719 ::InitPoints(vtkUnstructuredGrid *theInput,
00720              vtkIdType theCellId)
00721 {
00722   vtkPoints* aPoints = VTKViewer_Triangulator::InitPoints(theInput, theCellId);
00723 
00724   myPoints->Modified();
00725   myUnstructuredGrid->Modified();
00726   myGeometryFilter->Update();
00727   
00728   return aPoints;
00729 }
00730 
00731 
00732 //----------------------------------------------------------------------------
00733 vtkIdType 
00734 VTKViewer_DelaunayTriangulator
00735 ::GetNumFaces()
00736 {
00737   return myPolyData->GetNumberOfCells();
00738 }
00739 
00740 
00741 //----------------------------------------------------------------------------
00742 vtkCell*
00743 VTKViewer_DelaunayTriangulator
00744 ::GetFace(vtkIdType theFaceId)
00745 {
00746   return myPolyData->GetCell(theFaceId);
00747 }