-// Copyright (C) 2005 OPEN CASCADE, CEA/DEN, EDF R&D, PRINCIPIA R&D
-//
+// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
-//
-// This library is distributed in the hope that it will be useful
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
-// See http://www.salome-platform.org/
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
-// This library is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
-//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-//
-// See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
#include "VTKViewer_ConvexTool.h"
+#include "VTKViewer_GeometryFilter.h"
#include <set>
#include <map>
+#include <algorithm>
+#include <iterator>
#include <vtkUnstructuredGrid.h>
#include <vtkGeometryFilter.h>
#include <vtkCell.h>
#include <vtkPlane.h>
#include <vtkMath.h>
+#include <vtkCellArray.h>
+#include <vtkTriangle.h>
+#include <vtkOrderedTriangulator.h>
+
+#ifdef _DEBUG_
+static int DEBUG_TRIA_EXECUTE = 0;
+#else
+static int DEBUG_TRIA_EXECUTE = 0;
+#endif
namespace
{
struct TPolygon
{
TConnectivities myConnectivities;
- vtkFloatingPointType myOrigin[3];
- vtkFloatingPointType myNormal[3];
+ double myOrigin[3];
+ double myNormal[3];
TPolygon(const TConnectivities& theConnectivities,
- vtkFloatingPointType theOrigin[3],
- vtkFloatingPointType theNormal[3]):
+ double theOrigin[3],
+ double theNormal[3]):
myConnectivities(theConnectivities)
{
myOrigin[0] = theOrigin[0];
typedef std::vector<TPolygon> TPolygons;
}
-/*!
- Constructor
-*/
+
+//----------------------------------------------------------------------------
VTKViewer_Triangulator
::VTKViewer_Triangulator():
- myInput(NULL),
- myCellId(-1),
- myShowInside(-1),
- myAllVisible(-1),
- myCellsVisibility(NULL),
- myCellIds(vtkIdList::New())
+ myCellIds(vtkIdList::New()),
+ myPointIds(NULL),
+ myFaceIds(vtkIdList::New()),
+ myPoints(vtkPoints::New())
{}
-/*!
- Destructor
-*/
+//----------------------------------------------------------------------------
VTKViewer_Triangulator
::~VTKViewer_Triangulator()
{
myCellIds->Delete();
+ myFaceIds->Delete();
+ myPoints->Delete();
}
+//----------------------------------------------------------------------------
+vtkPoints*
+VTKViewer_Triangulator
+::InitPoints(vtkUnstructuredGrid *theInput,
+ vtkIdType theCellId)
+{
+ myPoints->Reset();
+ myPoints->Modified(); // the VTK bug
+
+ vtkIdType aNumPts;
+ theInput->GetCellPoints(theCellId, aNumPts, myPointIds);
+ if ( aNumPts > 0 ) {
+ double anAbsoluteCoord[3];
+ myPoints->SetNumberOfPoints(aNumPts);
+ vtkPoints *anInputPoints = theInput->GetPoints();
+ for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
+ anInputPoints->GetPoint(myPointIds[aPntId], anAbsoluteCoord);
+ myPoints->SetPoint(aPntId, anAbsoluteCoord);
+ }
+ }
+
+ return myPoints;
+}
+
+
+//----------------------------------------------------------------------------
+vtkIdType
+VTKViewer_Triangulator
+::GetNbOfPoints()
+{
+ return myPoints->GetNumberOfPoints();
+}
+
+
+//----------------------------------------------------------------------------
+vtkIdType
+VTKViewer_Triangulator
+::GetPointId(vtkIdType thePointId)
+{
+ return thePointId;
+}
+
+
+//----------------------------------------------------------------------------
+double
+VTKViewer_Triangulator
+::GetCellLength()
+{
+ double aBounds[6];
+ myPoints->GetBounds(aBounds);
+
+ double aCoordDiff[3];
+ aCoordDiff[0] = (aBounds[1] - aBounds[0]);
+ aCoordDiff[1] = (aBounds[3] - aBounds[2]);
+ aCoordDiff[2] = (aBounds[5] - aBounds[4]);
+
+ return sqrt(aCoordDiff[0]*aCoordDiff[0] +
+ aCoordDiff[1]*aCoordDiff[1] +
+ aCoordDiff[2]*aCoordDiff[2]);
+}
+
+//----------------------------------------------------------------------------
+void
+VTKViewer_Triangulator
+::GetCellNeighbors(vtkUnstructuredGrid *theInput,
+ vtkIdType theCellId,
+ vtkCell* theFace,
+ vtkIdList* theCellIds)
+{
+ myFaceIds->Reset();
+ vtkIdList *anIdList = theFace->PointIds;
+ myFaceIds->InsertNextId(myPointIds[anIdList->GetId(0)]);
+ myFaceIds->InsertNextId(myPointIds[anIdList->GetId(1)]);
+ myFaceIds->InsertNextId(myPointIds[anIdList->GetId(2)]);
+
+ theInput->GetCellNeighbors(theCellId, myFaceIds, theCellIds);
+}
+
+
+//----------------------------------------------------------------------------
+vtkIdType
+VTKViewer_Triangulator
+::GetConnectivity(vtkIdType thePntId)
+{
+ return myPointIds[thePntId];
+}
+
+
+//----------------------------------------------------------------------------
bool
VTKViewer_Triangulator
::Execute(vtkUnstructuredGrid *theInput,
- vtkCellData* thInputCD,
- vtkIdType theCellId,
- int theShowInside,
- int theAllVisible,
- const char* theCellsVisibility,
- vtkPolyData *theOutput,
- vtkCellData* theOutputCD,
- int theStoreMapping,
- std::vector<vtkIdType>& theVTK2ObjIds,
- bool theIsCheckConvex)
+ vtkCellData* thInputCD,
+ vtkIdType theCellId,
+ int theShowInside,
+ int theAllVisible,
+ int theAppendCoincident3D,
+ const char* theCellsVisibility,
+ vtkPolyData *theOutput,
+ vtkCellData* theOutputCD,
+ int theStoreMapping,
+ std::vector<vtkIdType>& theVTK2ObjIds,
+ std::vector< std::vector<vtkIdType> >& theDimension2VTK2ObjIds,
+ bool theIsCheckConvex)
{
- myInput = theInput;
- myCellId = theCellId;
- myShowInside = theShowInside;
- myAllVisible = theAllVisible;
- myCellsVisibility = theCellsVisibility;
-
- vtkPoints *aPoints = InitPoints();
+ vtkPoints *aPoints = InitPoints(theInput, theCellId);
vtkIdType aNumPts = GetNbOfPoints();
-
//cout<<"Triangulator - aNumPts = "<<aNumPts<<"\n";
+
if(DEBUG_TRIA_EXECUTE) cout<<"Triangulator - aNumPts = "<<aNumPts<<"\n";
if(aNumPts == 0)
return true;
// To calculate the bary center of the cell
- vtkFloatingPointType aCellCenter[3] = {0.0, 0.0, 0.0};
+ double aCellCenter[3] = {0.0, 0.0, 0.0};
{
- vtkFloatingPointType aPntCoord[3];
+ double aPntCoord[3];
for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
aPoints->GetPoint(GetPointId(aPntId),aPntCoord);
-
//cout<<"\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}\n";
+
if(DEBUG_TRIA_EXECUTE) cout<<"\taPntId = "<<GetPointId(aPntId)<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}\n";
aCellCenter[0] += aPntCoord[0];
aCellCenter[1] += aPntCoord[1];
aCellCenter[2] += aPntCoord[2];
aCellCenter[2] /= aNumPts;
}
- vtkFloatingPointType aCellLength = GetCellLength();
+ double aCellLength = GetCellLength();
int aNumFaces = GetNumFaces();
- static vtkFloatingPointType EPS = 1.0E-2;
- vtkFloatingPointType aDistEps = aCellLength * EPS;
-
//cout<<"\taCellLength = "<<aCellLength<<"; aDistEps = "<<aDistEps<<"\n";
+ static double EPS = 1.0E-2;
+ double aDistEps = aCellLength/3.0 * EPS;
+
if(DEBUG_TRIA_EXECUTE) cout<<"\taNumFaces = "<<aNumFaces<<"; aCellLength = "<<aCellLength<<"; aDistEps = "<<aDistEps<<"\n";
// To initialize set of points that belong to the cell
typedef std::set<vtkIdType> TPointIds;
for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) {
vtkCell* aFace = GetFace(aFaceId);
- GetCellNeighbors(theCellId, aFace, myCellIds);
- if((!myAllVisible && !myCellsVisibility[myCellIds->GetId(0)]) ||
- myCellIds->GetNumberOfIds() <= 0 || myShowInside)
+ GetCellNeighbors(theInput, theCellId, aFace, myCellIds);
+ bool process = myCellIds->GetNumberOfIds() <= 0 ? true : theAppendCoincident3D;
+ if((!theAllVisible && !theCellsVisibility[myCellIds->GetId(0)]) ||
+ myCellIds->GetNumberOfIds() <= 0 || theShowInside || process)
{
TPointIds aPointIds;
vtkIdList *anIdList = aFace->PointIds;
vtkIdList *anIdList = aFace->PointIds;
vtkIdType aNewPts[3] = {anIdList->GetId(0), anIdList->GetId(1), anIdList->GetId(2)};
-
+
// To initialize set of points for the plane where the trinangle face belong to
TPointIds aPointIds;
aPointIds.insert(aNewPts[0]);
// To get know, if the points of the trinagle were already observed
bool anIsObserved = aFace2PointIds.find(aPointIds) == aFace2PointIds.end();
- //cout<<"\taFaceId = "<<aFaceId<<"; anIsObserved = "<<anIsObserved;
- //cout<<"; aNewPts = {"<<aNewPts[0]<<", "<<aNewPts[1]<<", "<<aNewPts[2]<<"}\n";
+ if(DEBUG_TRIA_EXECUTE) {
+ cout<<"\taFaceId = "<<aFaceId<<"; anIsObserved = "<<anIsObserved;
+ cout<<"; aNewPts = {"<<aNewPts[0]<<", "<<aNewPts[1]<<", "<<aNewPts[2]<<"}\n";
+ }
if(!anIsObserved){
// To get coordinates of the points of the traingle face
- vtkFloatingPointType aCoord[3][3];
+ double aCoord[3][3];
aPoints->GetPoint(aNewPts[0],aCoord[0]);
aPoints->GetPoint(aNewPts[1],aCoord[1]);
aPoints->GetPoint(aNewPts[2],aCoord[2]);
- // To calculate plane normal
- vtkFloatingPointType aVector01[3] = { aCoord[1][0] - aCoord[0][0],
- aCoord[1][1] - aCoord[0][1],
- aCoord[1][2] - aCoord[0][2] };
+ /* To calculate plane normal for face (aFace)
+
+
+ ^ aNormal
+ |
+ | ^ aVector01
+ | /
+ /_________> aVector02
+
+
+ */
+ double aVector01[3] = { aCoord[1][0] - aCoord[0][0],
+ aCoord[1][1] - aCoord[0][1],
+ aCoord[1][2] - aCoord[0][2] };
- vtkFloatingPointType aVector02[3] = { aCoord[2][0] - aCoord[0][0],
- aCoord[2][1] - aCoord[0][1],
- aCoord[2][2] - aCoord[0][2] };
+ double aVector02[3] = { aCoord[2][0] - aCoord[0][0],
+ aCoord[2][1] - aCoord[0][1],
+ aCoord[2][2] - aCoord[0][2] };
+ vtkMath::Normalize(aVector01);
+ vtkMath::Normalize(aVector02);
+
// To calculate the normal for the triangle
- vtkFloatingPointType aNormal[3];
+ double aNormal[3];
vtkMath::Cross(aVector02,aVector01,aNormal);
vtkMath::Normalize(aNormal);
// To calculate what points belong to the plane
// To calculate bounds of the point set
- vtkFloatingPointType aCenter[3] = {0.0, 0.0, 0.0};
+ double aCenter[3] = {0.0, 0.0, 0.0};
{
- TPointIds::const_iterator anIter = anInitialPointIds.begin();
- TPointIds::const_iterator anEndIter = anInitialPointIds.end();
- for(; anIter != anEndIter; anIter++){
- vtkFloatingPointType aPntCoord[3];
- vtkIdType aPntId = *anIter;
- aPoints->GetPoint(aPntId,aPntCoord);
- vtkFloatingPointType aDist = vtkPlane::DistanceToPlane(aPntCoord,aNormal,aCoord[0]);
- //cout<<"\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}; aDist = "<<aDist<<"\n";
- if(fabs(aDist) < aDistEps){
- aPointIds.insert(aPntId);
- aCenter[0] += aPntCoord[0];
- aCenter[1] += aPntCoord[1];
- aCenter[2] += aPntCoord[2];
- }
- }
- int aNbPoints = aPointIds.size();
- aCenter[0] /= aNbPoints;
- aCenter[1] /= aNbPoints;
- aCenter[2] /= aNbPoints;
+ TPointIds::const_iterator anIter = anInitialPointIds.begin();
+ TPointIds::const_iterator anEndIter = anInitialPointIds.end();
+ for(; anIter != anEndIter; anIter++){
+ double aPntCoord[3];
+ vtkIdType aPntId = *anIter;
+ aPoints->GetPoint(aPntId,aPntCoord);
+
+ double aVector0Pnt[3] = { aPntCoord[0] - aCoord[0][0],
+ aPntCoord[1] - aCoord[0][1],
+ aPntCoord[2] - aCoord[0][2] };
+
+
+ vtkMath::Normalize(aVector0Pnt);
+
+ double aNormalPnt[3];
+ // calculate aNormalPnt
+ {
+ double aCosPnt01 = vtkMath::Dot(aVector0Pnt,aVector01);
+ double aCosPnt02 = vtkMath::Dot(aVector0Pnt,aVector02);
+ if(aCosPnt01<-1)
+ aCosPnt01 = -1;
+ if(aCosPnt01>1)
+ aCosPnt01 = 1;
+ if(aCosPnt02<-1)
+ aCosPnt02 = -1;
+ if(aCosPnt02>1)
+ aCosPnt02 = 1;
+
+ double aDist01,aDist02;// deflection from Pi/3 angle (equilateral triangle)
+ double aAngPnt01 = fabs(acos(aCosPnt01));
+ double aAngPnt02 = fabs(acos(aCosPnt02));
+
+ /* check that triangle similar to equilateral triangle
+ AOC or COB ?
+ aVector0Pnt = (OC)
+ aVector01 = (OB)
+ aVector02 = (OA)
+
+ B
+ ^ aVector01 C
+ | ^ aVector0Pnt
+ | _____/
+ | ___/
+ |/________> aVector02
+ O A
+ */
+ aDist01 = fabs(aAngPnt01-(vtkMath::Pi())/3.0);
+ aDist02 = fabs(aAngPnt02-(vtkMath::Pi())/3.0);
+
+ // caculate a normal for best triangle
+ if(aDist01 <= aDist02)
+ vtkMath::Cross(aVector0Pnt,aVector01,aNormalPnt);
+ else
+ vtkMath::Cross(aVector0Pnt,aVector02,aNormalPnt);
+
+ }
+
+ vtkMath::Normalize(aNormalPnt);
+
+ if(DEBUG_TRIA_EXECUTE)
+ cout<<"\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"};";
+
+ double aDist = vtkPlane::DistanceToPlane(aPntCoord,aNormal,aCoord[0]);
+ if(DEBUG_TRIA_EXECUTE) cout<<": aDist = "<<aDist;
+ if(fabs(aDist) < aDistEps){
+ aPointIds.insert(aPntId);
+ aCenter[0] += aPntCoord[0];
+ aCenter[1] += aPntCoord[1];
+ aCenter[2] += aPntCoord[2];
+ if(DEBUG_TRIA_EXECUTE) cout << "; Added = TRUE" << endl;
+ } else {
+ if(DEBUG_TRIA_EXECUTE) cout << "; Added = FALSE" << endl;
+ }
+ }
+ int aNbPoints = aPointIds.size();
+ aCenter[0] /= aNbPoints;
+ aCenter[1] /= aNbPoints;
+ aCenter[2] /= aNbPoints;
}
//To sinchronize orientation of the cell and its face
- vtkFloatingPointType aVectorC[3] = { aCenter[0] - aCellCenter[0],
- aCenter[1] - aCellCenter[1],
- aCenter[2] - aCellCenter[2] };
+ double aVectorC[3] = { aCenter[0] - aCellCenter[0],
+ aCenter[1] - aCellCenter[1],
+ aCenter[2] - aCellCenter[2] };
vtkMath::Normalize(aVectorC);
- vtkFloatingPointType aDot = vtkMath::Dot(aNormal,aVectorC);
- //cout<<"\t\taNormal = {"<<aNormal[0]<<", "<<aNormal[1]<<", "<<aNormal[2]<<"}";
- //cout<<"; aVectorC = {"<<aVectorC[0]<<", "<<aVectorC[1]<<", "<<aVectorC[2]<<"}\n";
- //cout<<"\t\taDot = "<<aDot<<"\n";
+ double aDot = vtkMath::Dot(aNormal,aVectorC);
+ if(DEBUG_TRIA_EXECUTE) {
+ cout<<"\t\taNormal = {"<<aNormal[0]<<", "<<aNormal[1]<<", "<<aNormal[2]<<"}";
+ cout<<"; aVectorC = {"<<aVectorC[0]<<", "<<aVectorC[1]<<", "<<aVectorC[2]<<"}\n";
+ cout<<"\t\taDot = "<<aDot<<"\n";
+ }
if(aDot > 0){
- aNormal[0] = -aNormal[0];
- aNormal[1] = -aNormal[1];
- aNormal[2] = -aNormal[2];
+ aNormal[0] = -aNormal[0];
+ aNormal[1] = -aNormal[1];
+ aNormal[2] = -aNormal[2];
}
// To calculate the primary direction for point set
- vtkFloatingPointType aVector0[3] = { aCoord[0][0] - aCenter[0],
- aCoord[0][1] - aCenter[1],
- aCoord[0][2] - aCenter[2] };
+ double aVector0[3] = { aCoord[0][0] - aCenter[0],
+ aCoord[0][1] - aCenter[1],
+ aCoord[0][2] - aCenter[2] };
vtkMath::Normalize(aVector0);
- //cout<<"\t\taCenter = {"<<aCenter[0]<<", "<<aCenter[1]<<", "<<aCenter[2]<<"}";
- //cout<<"; aVector0 = {"<<aVector0[0]<<", "<<aVector0[1]<<", "<<aVector0[2]<<"}\n";
+ if(DEBUG_TRIA_EXECUTE) {
+ cout<<"\t\taCenter = {"<<aCenter[0]<<", "<<aCenter[1]<<", "<<aCenter[2]<<"}";
+ cout<<"; aVector0 = {"<<aVector0[0]<<", "<<aVector0[1]<<", "<<aVector0[2]<<"}\n";
+ }
// To calculate the set of points by face those that belong to the plane
TFace2PointIds aRemoveFace2PointIds;
{
- TFace2PointIds::const_iterator anIter = aFace2PointIds.begin();
- TFace2PointIds::const_iterator anEndIter = aFace2PointIds.end();
- for(; anIter != anEndIter; anIter++){
- const TPointIds& anIds = *anIter;
- TPointIds anIntersection;
- std::set_intersection(aPointIds.begin(),aPointIds.end(),
- anIds.begin(),anIds.end(),
- std::inserter(anIntersection,anIntersection.begin()));
-
- if(anIntersection == anIds){
- aRemoveFace2PointIds.insert(anIds);
- }
- }
+ TFace2PointIds::const_iterator anIter = aFace2PointIds.begin();
+ TFace2PointIds::const_iterator anEndIter = aFace2PointIds.end();
+ for(; anIter != anEndIter; anIter++){
+ const TPointIds& anIds = *anIter;
+ TPointIds anIntersection;
+ std::set_intersection(aPointIds.begin(),aPointIds.end(),
+ anIds.begin(),anIds.end(),
+ std::inserter(anIntersection,anIntersection.begin()));
+
+
+ if(DEBUG_TRIA_EXECUTE) {
+ cout << "anIntersection:";
+ TPointIds::iterator aII = anIntersection.begin();
+ for(;aII!=anIntersection.end();aII++)
+ cout << *aII << ",";
+ cout << endl;
+ cout << "anIds :";
+ TPointIds::const_iterator aIIds = anIds.begin();
+ for(;aIIds!=anIds.end();aIIds++)
+ cout << *aIIds << ",";
+ cout << endl;
+ }
+ if(anIntersection == anIds){
+ aRemoveFace2PointIds.insert(anIds);
+ }
+ }
}
// To remove from the set of points by face those that belong to the plane
{
- TFace2PointIds::const_iterator anIter = aRemoveFace2PointIds.begin();
- TFace2PointIds::const_iterator anEndIter = aRemoveFace2PointIds.end();
- for(; anIter != anEndIter; anIter++){
- const TPointIds& anIds = *anIter;
- aFace2PointIds.erase(anIds);
- }
+ TFace2PointIds::const_iterator anIter = aRemoveFace2PointIds.begin();
+ TFace2PointIds::const_iterator anEndIter = aRemoveFace2PointIds.end();
+ for(; anIter != anEndIter; anIter++){
+ const TPointIds& anIds = *anIter;
+ aFace2PointIds.erase(anIds);
+ }
}
// To sort the planar set of the points accrding to the angle
{
- typedef std::map<vtkFloatingPointType,vtkIdType> TSortedPointIds;
- TSortedPointIds aSortedPointIds;
-
- TPointIds::const_iterator anIter = aPointIds.begin();
- TPointIds::const_iterator anEndIter = aPointIds.end();
- for(; anIter != anEndIter; anIter++){
- vtkFloatingPointType aPntCoord[3];
- vtkIdType aPntId = *anIter;
- aPoints->GetPoint(aPntId,aPntCoord);
- vtkFloatingPointType aVector[3] = { aPntCoord[0] - aCenter[0],
- aPntCoord[1] - aCenter[1],
- aPntCoord[2] - aCenter[2] };
- vtkMath::Normalize(aVector);
-
- vtkFloatingPointType aCross[3];
- vtkMath::Cross(aVector,aVector0,aCross);
- bool aGreaterThanPi = vtkMath::Dot(aCross,aNormal) < 0;
- vtkFloatingPointType aCosinus = vtkMath::Dot(aVector,aVector0);
- if(aCosinus > 1.0)
- aCosinus = 1.0;
- if(aCosinus < -1.0)
- aCosinus = -1.0;
- static vtkFloatingPointType a2Pi = 2.0 * vtkMath::Pi();
- vtkFloatingPointType anAngle = acos(aCosinus);
- //cout<<"\t\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}";
- //cout<<"; aGreaterThanPi = "<<aGreaterThanPi<<"; aCosinus = "<<aCosinus<<"; anAngle = "<<anAngle<<"\n";
- if(aGreaterThanPi){
- anAngle = a2Pi - anAngle;
- //cout<<"\t\t\t\tanAngle = "<<anAngle<<"\n";
- }
- aSortedPointIds[anAngle] = aPntId;
- }
-
- if(!aSortedPointIds.empty()){
- int aNumFacePts = aSortedPointIds.size();
- ::TConnectivities aConnectivities(aNumFacePts);
- TSortedPointIds::const_iterator anIter = aSortedPointIds.begin();
- TSortedPointIds::const_iterator anEndIter = aSortedPointIds.end();
- for(vtkIdType anId = 0; anIter != anEndIter; anIter++, anId++){
- vtkIdType aPntId = anIter->second;
- aConnectivities[anId] = GetConnectivity(aPntId);
- }
- aPolygons.push_back(::TPolygon(aConnectivities,aCenter,aNormal));
- }
+ typedef std::map<double,vtkIdType> TSortedPointIds;
+ TSortedPointIds aSortedPointIds;
+
+ TPointIds::const_iterator anIter = aPointIds.begin();
+ TPointIds::const_iterator anEndIter = aPointIds.end();
+ for(; anIter != anEndIter; anIter++){
+ double aPntCoord[3];
+ vtkIdType aPntId = *anIter;
+ aPoints->GetPoint(aPntId,aPntCoord);
+ double aVector[3] = { aPntCoord[0] - aCenter[0],
+ aPntCoord[1] - aCenter[1],
+ aPntCoord[2] - aCenter[2] };
+ vtkMath::Normalize(aVector);
+
+ double aCross[3];
+ vtkMath::Cross(aVector,aVector0,aCross);
+ double aCr = vtkMath::Dot(aCross,aNormal);
+ bool aGreaterThanPi = aCr < 0;
+ double aCosinus = vtkMath::Dot(aVector,aVector0);
+ double anAngle = 0.0;
+ if(aCosinus >= 1.0){
+ aCosinus = 1.0;
+ } else if (aCosinus <= -1.0){
+ aCosinus = -1.0;
+ anAngle = vtkMath::Pi();
+ } else {
+ anAngle = acos(aCosinus);
+ if(aGreaterThanPi)
+ anAngle = 2*vtkMath::Pi() - anAngle;
+ }
+
+ if(DEBUG_TRIA_EXECUTE) {
+ cout << "\t\t\t vtkMath::Dot(aCross,aNormal)="<<aCr<<endl;
+ cout<<"\t\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}";
+ cout<<"; aGreaterThanPi = "<<aGreaterThanPi<<"; aCosinus = "<<aCosinus<<"; anAngle = "<<anAngle<<"\n";
+ }
+ aSortedPointIds[anAngle] = aPntId;
+ }
+
+ if(!aSortedPointIds.empty()){
+ int aNumFacePts = aSortedPointIds.size();
+ ::TConnectivities aConnectivities(aNumFacePts);
+ TSortedPointIds::const_iterator anIter = aSortedPointIds.begin();
+ TSortedPointIds::const_iterator anEndIter = aSortedPointIds.end();
+ if(DEBUG_TRIA_EXECUTE) cout << "Polygon:";
+ for(vtkIdType anId = 0; anIter != anEndIter; anIter++, anId++){
+ vtkIdType aPntId = anIter->second;
+ aConnectivities[anId] = GetConnectivity(aPntId);
+ if(DEBUG_TRIA_EXECUTE) cout << aPntId << ",";
+ }
+ if(DEBUG_TRIA_EXECUTE) cout << endl;
+ aPolygons.push_back(::TPolygon(aConnectivities,aCenter,aNormal));
+ }
}
}
}
-
if(aPolygons.empty())
return true;
int aNbPolygons = aPolygons.size();
for (int aPolygonId = 0; aPolygonId < aNbPolygons; aPolygonId++) {
::TPolygon& aPolygon = aPolygons[aPolygonId];
- vtkFloatingPointType* aNormal = aPolygon.myNormal;
- vtkFloatingPointType* anOrigin = aPolygon.myOrigin;
- //cout<<"\taPolygonId = "<<aPolygonId<<"\n";
- //cout<<"\t\taNormal = {"<<aNormal[0]<<", "<<aNormal[1]<<", "<<aNormal[2]<<"}";
- //cout<<"; anOrigin = {"<<anOrigin[0]<<", "<<anOrigin[1]<<", "<<anOrigin[2]<<"}\n";
+ double* aNormal = aPolygon.myNormal;
+ double* anOrigin = aPolygon.myOrigin;
+ if(DEBUG_TRIA_EXECUTE) {
+ cout<<"\taPolygonId = "<<aPolygonId<<"\n";
+ cout<<"\t\taNormal = {"<<aNormal[0]<<", "<<aNormal[1]<<", "<<aNormal[2]<<"}";
+ cout<<"; anOrigin = {"<<anOrigin[0]<<", "<<anOrigin[1]<<", "<<anOrigin[2]<<"}\n";
+ }
for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++){
- vtkFloatingPointType aPntCoord[3];
- vtkIdType anId = GetPointId(aPntId);
- aPoints->GetPoint(anId,aPntCoord);
- vtkFloatingPointType aDist = vtkPlane::Evaluate(aNormal,anOrigin,aPntCoord);
-
//cout<<"\t\taPntId = "<<anId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}; aDist = "<<aDist<<"\n";
- if(aDist < -aDistEps)
- return false;
+ double aPntCoord[3];
+ vtkIdType anId = GetPointId(aPntId);
+ aPoints->GetPoint(anId,aPntCoord);
+ double aDist = vtkPlane::Evaluate(aNormal,anOrigin,aPntCoord);
+
if(DEBUG_TRIA_EXECUTE) cout<<"\t\taPntId = "<<anId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}; aDist = "<<aDist<<"\n";
+ if(aDist < -aDistEps)
+ return false;
}
}
}
int aNbPolygons = aPolygons.size();
for (int aPolygonId = 0; aPolygonId < aNbPolygons; aPolygonId++) {
::TPolygon& aPolygon = aPolygons[aPolygonId];
+ if(DEBUG_TRIA_EXECUTE) cout << "PoilygonId="<<aPolygonId<<" | ";
TConnectivities& aConnectivities = aPolygon.myConnectivities;
+ if(DEBUG_TRIA_EXECUTE) {
+ for(size_t i=0;i<aConnectivities.size();i++)
+ cout << aConnectivities[i] << ",";
+ cout << endl;
+ }
int aNbPoints = aConnectivities.size();
vtkIdType aNewCellId = theOutput->InsertNextCell(VTK_POLYGON,aNbPoints,&aConnectivities[0]);
if(theStoreMapping)
- theVTK2ObjIds.push_back(theCellId);
+ VTKViewer_GeometryFilter::InsertId( theCellId, VTK_POLYGON, theVTK2ObjIds, theDimension2VTK2ObjIds );
theOutputCD->CopyData(thInputCD,theCellId,aNewCellId);
}
}
- //cout<<"\tTriangulator - Ok\n";
+ if(DEBUG_TRIA_EXECUTE) cout<<"\tTriangulator - Ok\n";
+
return true;
}
-/*!
- Constructor
-*/
+
+//----------------------------------------------------------------------------
VTKViewer_OrderedTriangulator
::VTKViewer_OrderedTriangulator():
- myCell(vtkGenericCell::New())
-{}
+ myTriangulator(vtkOrderedTriangulator::New()),
+ myBoundaryTris(vtkCellArray::New()),
+ myTriangle(vtkTriangle::New())
+{
+ myBoundaryTris->Allocate(VTK_CELL_SIZE);
+ myTriangulator->PreSortedOff();
+}
-/*!
- Destructor
-*/
+
+//----------------------------------------------------------------------------
VTKViewer_OrderedTriangulator
::~VTKViewer_OrderedTriangulator()
{
- myCell->Delete();
+ myTriangle->Delete();
+ myBoundaryTris->Delete();
+ myTriangulator->Delete();
}
+
+//----------------------------------------------------------------------------
vtkPoints*
VTKViewer_OrderedTriangulator
-::InitPoints()
+::InitPoints(vtkUnstructuredGrid *theInput,
+ vtkIdType theCellId)
{
- myInput->GetCell(myCellId,myCell);
- return myInput->GetPoints();
-}
+ myBoundaryTris->Reset();
+
+ vtkPoints* aPoints = VTKViewer_Triangulator::InitPoints(theInput, theCellId);
+ vtkIdType aNumPts = myPoints->GetNumberOfPoints();
+ if ( aNumPts > 0 ) {
+ myTriangulator->InitTriangulation(0.0, 1.0, 0.0, 1.0, 0.0, 1.0, aNumPts);
+
+ double aBounds[6];
+ myPoints->GetBounds(aBounds);
+ double xSize, ySize, zSize;
+ xSize = aBounds[1] - aBounds[0];
+ ySize = aBounds[3] - aBounds[2];
+ zSize = aBounds[5] - aBounds[4];
+ double anAbsoluteCoord[3];
+ double aParamentrucCoord[3];
+ for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
+ myPoints->GetPoint(aPntId, anAbsoluteCoord);
+ aParamentrucCoord[0] = xSize==0. ? 0. : ((anAbsoluteCoord[0] - aBounds[0]) / xSize);
+ aParamentrucCoord[1] = ySize==0. ? 0. : ((anAbsoluteCoord[1] - aBounds[2]) / ySize);
+ aParamentrucCoord[2] = zSize==0. ? 0. : ((anAbsoluteCoord[2] - aBounds[4]) / zSize);
+ myTriangulator->InsertPoint(aPntId, anAbsoluteCoord, aParamentrucCoord, 0);
+ }
-vtkIdType
-VTKViewer_OrderedTriangulator
-::GetNbOfPoints()
-{
- return myCell->GetNumberOfPoints();
-}
+ myTriangulator->Triangulate();
+ myTriangulator->AddTriangles(myBoundaryTris);
+ }
-vtkIdType
-VTKViewer_OrderedTriangulator
-::GetPointId(vtkIdType thePointId)
-{
- return myCell->GetPointId(thePointId);
+ return aPoints;
}
-vtkFloatingPointType
-VTKViewer_OrderedTriangulator
-::GetCellLength()
-{
- return sqrt(myCell->GetLength2());
-}
+//----------------------------------------------------------------------------
vtkIdType
VTKViewer_OrderedTriangulator
::GetNumFaces()
{
- return myCell->GetNumberOfFaces();
+ return myBoundaryTris->GetNumberOfCells();
}
+
+//----------------------------------------------------------------------------
vtkCell*
VTKViewer_OrderedTriangulator
::GetFace(vtkIdType theFaceId)
{
- return myCell->GetFace(theFaceId);
-}
+ vtkIdType aNumCells = myBoundaryTris->GetNumberOfCells();
+ if ( theFaceId < 0 || theFaceId >= aNumCells )
+ return NULL;
-void
-VTKViewer_OrderedTriangulator
-::GetCellNeighbors(vtkIdType theCellId,
- vtkCell* theFace,
- vtkIdList* theCellIds)
-{
- vtkIdList *anIdList = theFace->PointIds;
- myInput->GetCellNeighbors(theCellId, anIdList, theCellIds);
-}
+ vtkIdType *aCells = myBoundaryTris->GetPointer();
-vtkIdType
-VTKViewer_OrderedTriangulator
-::GetConnectivity(vtkIdType thePntId)
-{
- return thePntId;
+ // Each triangle has three points plus number of points
+ vtkIdType *aCellPtr = aCells + 4*theFaceId;
+
+ myTriangle->PointIds->SetId(0, aCellPtr[1]);
+ myTriangle->Points->SetPoint(0, myPoints->GetPoint(aCellPtr[1]));
+
+ myTriangle->PointIds->SetId(1, aCellPtr[2]);
+ myTriangle->Points->SetPoint(1, myPoints->GetPoint(aCellPtr[2]));
+
+ myTriangle->PointIds->SetId(2, aCellPtr[3]);
+ myTriangle->Points->SetPoint(2, myPoints->GetPoint(aCellPtr[3]));
+
+ return myTriangle;
}
-/*!
- Constructor
-*/
+
+//----------------------------------------------------------------------------
VTKViewer_DelaunayTriangulator
::VTKViewer_DelaunayTriangulator():
myUnstructuredGrid(vtkUnstructuredGrid::New()),
myGeometryFilter(vtkGeometryFilter::New()),
myDelaunay3D(vtkDelaunay3D::New()),
- myFaceIds(vtkIdList::New()),
- myPoints(vtkPoints::New()),
- myPolyData(NULL),
- myPointIds(NULL)
+ myPolyData(NULL)
{
- myDelaunay3D->SetInput(myUnstructuredGrid);
- myGeometryFilter->SetInput(myDelaunay3D->GetOutput());
-}
+ myUnstructuredGrid->Initialize();
+ myUnstructuredGrid->Allocate();
+ myUnstructuredGrid->SetPoints(myPoints);
+ myDelaunay3D->SetInputData(myUnstructuredGrid);
+ myGeometryFilter->SetInputConnection(myDelaunay3D->GetOutputPort());
+ myPolyData = myGeometryFilter->GetOutput();
+}
-/*!
- Destructor
-*/
+//----------------------------------------------------------------------------
VTKViewer_DelaunayTriangulator
::~VTKViewer_DelaunayTriangulator()
{
myUnstructuredGrid->Delete();
myGeometryFilter->Delete();
myDelaunay3D->Delete();
- myFaceIds->Delete();
- myPoints->Delete();
}
+//----------------------------------------------------------------------------
vtkPoints*
VTKViewer_DelaunayTriangulator
-::InitPoints()
+::InitPoints(vtkUnstructuredGrid *theInput,
+ vtkIdType theCellId)
{
- myUnstructuredGrid->Initialize();
- myUnstructuredGrid->Allocate();
- myUnstructuredGrid->SetPoints(myPoints);
-
- vtkIdType aNumPts;
- myInput->GetCellPoints(myCellId,aNumPts,myPointIds);
- {
- vtkFloatingPointType aPntCoord[3];
- myPoints->SetNumberOfPoints(aNumPts);
- vtkPoints *anInputPoints = myInput->GetPoints();
- for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
- anInputPoints->GetPoint(myPointIds[aPntId],aPntCoord);
- myPoints->SetPoint(aPntId,aPntCoord);
- }
- }
+ vtkPoints* aPoints = VTKViewer_Triangulator::InitPoints(theInput, theCellId);
+ myPoints->Modified();
+ myUnstructuredGrid->Modified();
myGeometryFilter->Update();
- myPolyData = myGeometryFilter->GetOutput();
-
- return myPoints;
-}
-
-vtkIdType
-VTKViewer_DelaunayTriangulator
-::GetNbOfPoints()
-{
- return myPoints->GetNumberOfPoints();
-}
-
-vtkIdType
-VTKViewer_DelaunayTriangulator
-::GetPointId(vtkIdType thePointId)
-{
- return thePointId;
+
+ return aPoints;
}
-vtkFloatingPointType
-VTKViewer_DelaunayTriangulator
-::GetCellLength()
-{
- return myPolyData->GetLength();
-}
+//----------------------------------------------------------------------------
vtkIdType
VTKViewer_DelaunayTriangulator
::GetNumFaces()
return myPolyData->GetNumberOfCells();
}
+
+//----------------------------------------------------------------------------
vtkCell*
VTKViewer_DelaunayTriangulator
::GetFace(vtkIdType theFaceId)
{
return myPolyData->GetCell(theFaceId);
}
-
-void
-VTKViewer_DelaunayTriangulator
-::GetCellNeighbors(vtkIdType theCellId,
- vtkCell* theFace,
- vtkIdList* theCellIds)
-{
- myFaceIds->Reset();
- vtkIdList *anIdList = theFace->PointIds;
- myFaceIds->InsertNextId(myPointIds[anIdList->GetId(0)]);
- myFaceIds->InsertNextId(myPointIds[anIdList->GetId(1)]);
- myFaceIds->InsertNextId(myPointIds[anIdList->GetId(2)]);
-
- myInput->GetCellNeighbors(theCellId, myFaceIds, theCellIds);
-}
-
-
-vtkIdType
-VTKViewer_DelaunayTriangulator
-::GetConnectivity(vtkIdType thePntId)
-{
- return myPointIds[thePntId];
-}
