X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FVTKViewer%2FVTKViewer_ConvexTool.cxx;h=f15527b73d83c9aae6ab574adfecafc8b3c4cf2a;hb=refs%2Fheads%2Fngr%2Fpython3_dev_pv5.4;hp=1c4a8a710bb95594d759932fbaaa1931458db599;hpb=84e739898389c91887f9f6325dacd12348efb35e;p=modules%2Fgui.git diff --git a/src/VTKViewer/VTKViewer_ConvexTool.cxx b/src/VTKViewer/VTKViewer_ConvexTool.cxx index 1c4a8a710..f15527b73 100644 --- a/src/VTKViewer/VTKViewer_ConvexTool.cxx +++ b/src/VTKViewer/VTKViewer_ConvexTool.cxx @@ -1,601 +1,748 @@ -// 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 +// 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, 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 +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// #include "VTKViewer_ConvexTool.h" - -#include -#include -#include -#include +#include "VTKViewer_GeometryFilter.h" #include -#include +#include #include -#include - -typedef std::set TUIDS; // unique ids -typedef std::map TPTOIDS; // id points -> unique ids - -namespace CONVEX_TOOL -{ +#include -static float FACE_ANGLE_TOLERANCE=1.5; +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #ifdef _DEBUG_ - static int MYDEBUG = 0; +static int DEBUG_TRIA_EXECUTE = 0; #else - static int MYDEBUG = 0; +static int DEBUG_TRIA_EXECUTE = 0; #endif -/*! \fn static void GetCenter(vtkUnstructuredGrid* theGrid,TCell theptIds,float *center) - * \brief Calculation of geometry center. - * \param theGrid - vtkUnstructuredGrid cell. - * \param theptIds - point ids. - * \retval center - output array[3] with coordinates of geometry center. - */ -static void GetCenter(vtkUnstructuredGrid* theGrid,TCell theptIds,float *center) +namespace { - float *p; - center[0] = center[1] = center[2] = 0.0; - - int numIds=theptIds.size(); - - // get the center of the cell - for (int i=0; i < numIds; i++) - { - p = theGrid->GetPoint(theptIds[i]); - for (int j=0; j < 3; j++) - { - center[j] += p[j]; - } - } - for (int j=0; j<3; j++) + typedef std::vector TConnectivities; + + struct TPolygon + { + TConnectivities myConnectivities; + double myOrigin[3]; + double myNormal[3]; + TPolygon(const TConnectivities& theConnectivities, + double theOrigin[3], + double theNormal[3]): + myConnectivities(theConnectivities) { - center[j] /= numIds; + myOrigin[0] = theOrigin[0]; + myOrigin[1] = theOrigin[1]; + myOrigin[2] = theOrigin[2]; + + myNormal[0] = theNormal[0]; + myNormal[1] = theNormal[1]; + myNormal[2] = theNormal[2]; } + }; + + typedef std::vector TPolygons; +} + + +//---------------------------------------------------------------------------- +VTKViewer_Triangulator +::VTKViewer_Triangulator(): + myCellIds(vtkIdList::New()), + myPointIds(NULL), + myFaceIds(vtkIdList::New()), + myPoints(vtkPoints::New()) +{} + + +//---------------------------------------------------------------------------- +VTKViewer_Triangulator +::~VTKViewer_Triangulator() +{ + myCellIds->Delete(); + myFaceIds->Delete(); + myPoints->Delete(); } -/*! \fn static void ReverseIds(TCell &theIds) - * \brief Reverse ids. - * \param theIds - points ids. - * \retval theIds - example input:(1,2,3,4) -> output:(4,3,2,1) - */ -static void ReverseIds(TCell &theIds) + +//---------------------------------------------------------------------------- +vtkPoints* +VTKViewer_Triangulator +::InitPoints(vtkUnstructuredGrid *theInput, + vtkIdType theCellId) { - int i; - vtkIdType tmp; - vtkIdType npts=theIds.size(); - - for(i=0;i<(npts/2);i++){ - tmp = theIds[i]; - theIds[i] = theIds[npts-i-1]; - theIds[npts-i-1] = tmp; + 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; } -/*! \fn void GetFriends(const TPTOIDS p2faces,const TCellArray f2points,TPTOIDS& face2face_output) - * \brief Caclulation of connected faces (faceId -> (faceId1,faceId2, ...)) - * \param p2faces - point to faces ids map. - * \param f2points - faces to points ids map. - * \retval face2face_output - faces to faces ids map. - */ -void GetFriends(const TPTOIDS p2faces,const TCellArray f2points,TPTOIDS& face2face_output) + +//---------------------------------------------------------------------------- +vtkIdType +VTKViewer_Triangulator +::GetNbOfPoints() { - TCellArray::const_iterator f2pIter = f2points.begin(); + return myPoints->GetNumberOfPoints(); +} - for( ; f2pIter!=f2points.end() ; f2pIter++ ){ - vtkIdType faceId = f2pIter->first; - TCell face_points = f2pIter->second; - int nb_face_points = face_points.size(); - - vtkIdType id1; - vtkIdType id2; - TPTOIDS::const_iterator faces1; - TPTOIDS::const_iterator faces2; - - id1 = face_points[0]; - faces1 = p2faces.find(id1); - - TUIDS output_faces; - - for(int i=1 ; isecond.begin(), faces1->second.end(), faces2->second.begin(), faces2->second.end(), - std::inserter(output_faces,output_faces.begin())); - - id1 = id2; - faces1 = faces2; - } - id1 = face_points[0]; - faces1 = p2faces.find(id1); - std::set_intersection(faces1->second.begin(), faces1->second.end(), faces2->second.begin(), faces2->second.end(), - std::inserter(output_faces,output_faces.begin())); - - output_faces.erase(faceId); // erase the face id for which we found friends - if(MYDEBUG){ - cout << "fId[" << faceId <<"]: "; - std::copy(output_faces.begin(), output_faces.end(), std::ostream_iterator(cout, " ")); - cout << endl; - } - - face2face_output[faceId] = output_faces; - } +//---------------------------------------------------------------------------- +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]); } -/*! \fn bool IsConnectedFacesOnOnePlane( vtkUnstructuredGrid* theGrid,vtkIdType theFId1, vtkIdType theFId2,TUIDS FpIds1, TUIDS FpIds2 ) - * \brief Check is connected faces on one plane. - * \param theGrid - vtkUnstructuredGrid - * \param theFId1 - id of first face - * \param theFId2 - id of second face - * \param FpIds1 - first face points ids. - * \param FpIds2 - second face points ids. - * \return TRUE if two faces on one plane, else FALSE. - */ -bool IsConnectedFacesOnOnePlane( vtkUnstructuredGrid* theGrid, - vtkIdType theFId1, vtkIdType theFId2, - TUIDS FpIds1, TUIDS FpIds2 ) + +//---------------------------------------------------------------------------- +void +VTKViewer_Triangulator +::GetCellNeighbors(vtkUnstructuredGrid *theInput, + vtkIdType theCellId, + vtkCell* theFace, + vtkIdList* theCellIds) { - bool status = false; - TUIDS common_ids; - std::set_intersection(FpIds1.begin(), FpIds1.end(), FpIds2.begin(), FpIds2.end(), - std::inserter(common_ids,common_ids.begin())); - - /* Number of common ids = 2 (A1,A2) - - - _ _ _ _ _ _ _ _ _ vectors: - | \ / | v1 {A2,A1} - \ / v2 {A1,B1} - | | A2 | v3 {A1,C1} - | - | | | - | - | | A1 | - / \ - |_ _ _ _ _/ \_ _ _ _ _| - B2 B1 C1 C2 - - */ - TUIDS::iterator common_iter = common_ids.begin(); - if(common_ids.size() == 2){ - TUIDS::iterator loc_id1_0 = FpIds1.find(*(common_iter)); - common_iter++; - TUIDS::iterator loc_id1_1 = FpIds1.find(*(common_iter)); - - TUIDS::iterator loc_id2_0 = FpIds1.begin(); - TUIDS::iterator loc_id2_1 = FpIds2.begin(); - - vtkIdType A1 = *loc_id1_0; - vtkIdType A2 = *loc_id1_1; - vtkIdType B1; - vtkIdType C1; - - for(;loc_id2_0!=FpIds1.end();loc_id2_0++) - if(*loc_id2_0 != A1 && *loc_id2_0!= A2){ - B1 = *loc_id2_0; - break; - } - for(;loc_id2_1!=FpIds2.end();loc_id2_1++) - if(*loc_id2_1 != A1 && *loc_id2_1!= A2){ - C1 = *loc_id2_1; - break; - } - if(MYDEBUG) cout <"; - float *p[4]; - float v1[3],v2[3],v3[3]; - p[0] = theGrid->GetPoint(A1); - p[1] = theGrid->GetPoint(A2); - p[2] = theGrid->GetPoint(B1); - p[3] = theGrid->GetPoint(C1); - - for(int i=0;i<3;i++){ - v1[i] = p[1][i] - p[0][i]; - v2[i] = p[2][i] - p[0][i]; - v3[i] = p[3][i] - p[0][i]; - } - - - float vec_b1[3]; - vtkMath::Cross(v2,v1,vec_b1); - float vec_b2[3]; - vtkMath::Cross(v1,v3,vec_b2); + myFaceIds->Reset(); + vtkIdList *anIdList = theFace->PointIds; + myFaceIds->InsertNextId(myPointIds[anIdList->GetId(0)]); + myFaceIds->InsertNextId(myPointIds[anIdList->GetId(1)]); + myFaceIds->InsertNextId(myPointIds[anIdList->GetId(2)]); - float b1 = vtkMath::Norm(vec_b1); + theInput->GetCellNeighbors(theCellId, myFaceIds, theCellIds); +} - float b2 = vtkMath::Norm(vec_b2); - float aCos = vtkMath::Dot(vec_b1,vec_b2)/(b1*b2); - - float angle=90.0; - angle = aCos>=1.0 ? 0.0 : 180*acosf(aCos)/vtkMath::Pi(); - - if( angle <= FACE_ANGLE_TOLERANCE) - status = true; - if (MYDEBUG){ - for(int k=0;k<4;k++){ - cout << " ("; - for(int j=0;j<2;j++){ - cout << p[k][j] << ","; - } - cout << p[k][2] << ") "; - } - cout << "angle="<& theVTK2ObjIds, + std::vector< std::vector >& theDimension2VTK2ObjIds, + bool theIsCheckConvex) { - if (new_faces.find(faceId) != new_faces.end()) return; - - new_faces.insert(new_faces.begin(),faceId); - new_faces_v2.push_back(faceId); - - TPTOIDS::const_iterator aIter1 = theFaces.find(faceId); - if(aIter1!=theFaces.end()){ - TUIDS::const_iterator aIter2 = (aIter1->second).begin(); - for(;aIter2!=(aIter1->second).end();aIter2++){ - if (new_faces.find(*aIter2) != new_faces.end()) continue; - GetAllFacesOnOnePlane(theFaces,*aIter2, - new_faces,new_faces_v2); // recurvise + vtkPoints *aPoints = InitPoints(theInput, theCellId); + vtkIdType aNumPts = GetNbOfPoints(); + if(DEBUG_TRIA_EXECUTE) cout<<"Triangulator - aNumPts = "<GetPoint(GetPointId(aPntId),aPntCoord); + if(DEBUG_TRIA_EXECUTE) cout<<"\taPntId = "<(cout, " ")); - if(MYDEBUG) cout << "\tv2:"; - if(MYDEBUG) std::copy(v2.begin(), v2.end(), std::ostream_iterator(cout, " ")); - if(MYDEBUG) cout << endl; - - TUIDS v1_set; - std::copy(v1.begin(), v1.end(), std::inserter(v1_set,v1_set.begin())); - TUIDS v2_set; - std::copy(v2.begin(), v2.end(), std::inserter(v2_set,v2_set.begin())); - TUIDS tmpIntersection; - std::set_intersection(v1_set.begin(),v1_set.end(),v2_set.begin(),v2_set.end(), std::inserter(tmpIntersection,tmpIntersection.begin())); - if(MYDEBUG) std::copy(tmpIntersection.begin(),tmpIntersection.end(), std::ostream_iterator(cout, " ")); - if(MYDEBUG) cout << endl; - - if(tmpIntersection.size() < 2) - if(MYDEBUG) cout << __FILE__ << "[" << __LINE__ << "]: Warning ! Wrong ids" << endl; + double aCellLength = GetCellLength(); + int aNumFaces = GetNumFaces(); + + static double EPS = 1.0E-2; + double aDistEps = aCellLength/3.0 * EPS; + if(DEBUG_TRIA_EXECUTE) cout<<"\taNumFaces = "< TPointIds; + TPointIds anInitialPointIds; + for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++) + anInitialPointIds.insert(GetPointId(aPntId)); - TCell::iterator v1_iter = v1.begin(); + // To initialize set of points by face that belong to the cell and backward + typedef std::set TFace2Visibility; + TFace2Visibility aFace2Visibility; - for(;v1_iter!=v1.end();v1_iter++){ - - vtkIdType curr_id = *v1_iter; - - output.push_back(curr_id); + typedef std::set TFace2PointIds; + TFace2PointIds aFace2PointIds; + + for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) { + vtkCell* aFace = GetFace(aFaceId); - if(tmpIntersection.find(curr_id) != tmpIntersection.end()){ - TCell::iterator v1_iter_tmp; - v1_iter_tmp = v1_iter; - v1_iter++; - - if(v1_iter==v1.end()) v1_iter=v1.begin(); - - curr_id = *v1_iter; - - if(tmpIntersection.find(curr_id) != tmpIntersection.end()){ - TCell::iterator v2_iter = v2.begin(); - for(;v2_iter!=v2.end();v2_iter++){ - vtkIdType v2_id = *v2_iter; - if(tmpIntersection.find(v2_id) == tmpIntersection.end()) - output.push_back(v2_id); - } - } - - v1_iter = v1_iter_tmp; - curr_id = *v1_iter; + 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; + aPointIds.insert(anIdList->GetId(0)); + aPointIds.insert(anIdList->GetId(1)); + aPointIds.insert(anIdList->GetId(2)); + aFace2PointIds.insert(aPointIds); + aFace2Visibility.insert(aFaceId); } } - if(MYDEBUG) cout << "Result: " ; - if(MYDEBUG) std::copy(output.begin(),output.end(),std::ostream_iterator(cout, " ")); - if(MYDEBUG) cout << endl; -} -void GetPolygonalFaces(vtkUnstructuredGrid* theGrid,int cellId,TCellArray &outputCellArray) -{ - if (theGrid->GetCellType(cellId) == VTK_CONVEX_POINT_SET){ - // get vtkCell type = VTK_CONVEX_POINT_SET - if(vtkConvexPointSet* convex = static_cast(theGrid->GetCell(cellId))){ - TCellArray f2points; - float convex_center[3]; // convex center point coorinat - int aNbFaces = convex->GetNumberOfFaces(); - int numPts = convex->GetNumberOfPoints(); - TCell convex_ids; - TPTOIDS p2faces; // key=pointId , value facesIds set + ::TPolygons aPolygons; + + for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) { + if(aFace2Visibility.find(aFaceId) == aFace2Visibility.end()) + continue; + + vtkCell* aFace = GetFace(aFaceId); + + 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]); + aPointIds.insert(aNewPts[1]); + aPointIds.insert(aNewPts[2]); + + // To get know, if the points of the trinagle were already observed + bool anIsObserved = aFace2PointIds.find(aPointIds) == aFace2PointIds.end(); + if(DEBUG_TRIA_EXECUTE) { + cout<<"\taFaceId = "<GetPoint(aNewPts[0],aCoord[0]); + aPoints->GetPoint(aNewPts[1],aCoord[1]); + aPoints->GetPoint(aNewPts[2],aCoord[2]); - for(int i=0;iGetPointId(i)); + /* 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] }; - GetCenter(theGrid,convex_ids,convex_center); - - for (vtkIdType faceId=0; faceId < aNbFaces; faceId++){ - vtkCell *aFace = convex->GetFace(faceId); - int numFacePts = aFace->GetNumberOfPoints(); - TCell aIds; - - int i = 0; - for(i=0;iGetPointId(i)); - - float v_a[3],v_b[3],v_convex2face[3]; // vectors - float *id_0,*id_1,*id_n; - /*============================================= - ,+- - - - _ - _ / id_n | v_b {id_0,id_n} - v_b / _ - / | v_a {id_0,id_1} - / - / | - + id_0 - \ - _ \ | - v_a \ - \ id_1 | - "+- - - - - - ============================================*/ - id_0 = theGrid->GetPoint(aIds[0]); - id_1 = theGrid->GetPoint(aIds[1]); - id_n = theGrid->GetPoint(aIds[numFacePts-1]); - - for(i=0;i<3;i++){ - v_a[i] = id_1[i] - id_0[i]; - v_b[i] = id_n[i] - id_0[i]; - v_convex2face[i] = id_0[i] - convex_center[i]; - } - - if (vtkMath::Determinant3x3(v_a,v_b,v_convex2face) < 0){ - ReverseIds(aIds); - } - - for(i=0;i<(int)aIds.size();i++){ - TUIDS &acell = p2faces[aIds[i]]; - acell.insert(faceId); - } - - f2points[faceId] = aIds; + 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 + double aNormal[3]; + vtkMath::Cross(aVector02,aVector01,aNormal); - TPTOIDS face2face; - GetFriends(p2faces,f2points,face2face); + vtkMath::Normalize(aNormal); - TPTOIDS face2points; + // To calculate what points belong to the plane + // To calculate bounds of the point set + 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++){ + 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 = "<second).begin(); - points_iter!=(f2points_iter->second).end(); - points_iter++) - tmp.insert(*points_iter); - - face2points[f2points_iter->first] = tmp; - } // end copy - + //To sinchronize orientation of the cell and its face + double aVectorC[3] = { aCenter[0] - aCellCenter[0], + aCenter[1] - aCellCenter[1], + aCenter[2] - aCellCenter[2] }; + vtkMath::Normalize(aVectorC); - TPTOIDS new_face2faces; // which connected and in one plane + double aDot = vtkMath::Dot(aNormal,aVectorC); + if(DEBUG_TRIA_EXECUTE) { + cout<<"\t\taNormal = {"< 0){ + aNormal[0] = -aNormal[0]; + aNormal[1] = -aNormal[1]; + aNormal[2] = -aNormal[2]; + } - TPTOIDS::const_iterator aF2FIter = face2face.begin(); - for(;aF2FIter!=face2face.end();aF2FIter++){ - vtkIdType f_key = aF2FIter->first; - TUIDS &faces = new_face2faces[f_key]; - //faces.insert(f_key); - TUIDS f_friends = aF2FIter->second; - TUIDS::const_iterator a_friends_iter = f_friends.begin(); - for(;a_friends_iter!=f_friends.end();a_friends_iter++){ - vtkIdType friend_id = *a_friends_iter; - if( IsConnectedFacesOnOnePlane(theGrid,f_key,friend_id, - (face2points.find(f_key))->second, - (face2points.find(friend_id))->second)){ - faces.insert(friend_id); - } // end if - - } // end a_friends_iter - } // end aF2FIter + // To calculate the primary direction for point set + double aVector0[3] = { aCoord[0][0] - aCenter[0], + aCoord[0][1] - aCenter[1], + aCoord[0][2] - aCenter[2] }; + vtkMath::Normalize(aVector0); - if(MYDEBUG) - { - TPTOIDS::const_iterator new_face2face_iter = new_face2faces.begin(); - cout << "Connected faces and on plane:" << endl; - for(;new_face2face_iter!=new_face2faces.end();new_face2face_iter++){ - cout << "Group ["<first<<"] :"; - TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin(); - for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++) - cout << " " << *new_faces_iter ; - cout << endl; - } + if(DEBUG_TRIA_EXECUTE) { + cout<<"\t\taCenter = {"<first) != already_in.end()) - continue; - if(new_face2face_iter->second.size() > 1) - continue; - - TCell &tmp_v2 = output_newid2face_v2[k]; - tmp_v2.push_back(new_face2face_iter->first); - already_in.insert(new_face2face_iter->first); - - TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin(); - for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){ - if(already_in.find(*new_faces_iter) != already_in.end()) continue; - already_in.insert(*new_faces_iter); - - already_in_tmp.clear(); - already_in_tmp.insert(new_face2face_iter->first); - - TUIDS &tmp = output_newid2face[k]; - GetAllFacesOnOnePlane(new_face2faces,*new_faces_iter, - already_in_tmp,tmp_v2); - - for(TUIDS::const_iterator aIter=already_in_tmp.begin(); - aIter!=already_in_tmp.end(); - aIter++) - { - already_in.insert(*aIter); - tmp.insert(*aIter); - } - } - k++; - } + 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); + } + } } - if(MYDEBUG) { - cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<first<<"] :"; - TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin(); - for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++) - cout << " " << *new_faces_iter ; - cout << endl; - } - } - cout << "++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++"<first<<"] :"; - TCell::const_iterator new_faces_iter = (new_face2face_iter->second).begin(); - for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++) - cout << " " << *new_faces_iter ; - cout << endl; - } - } + // 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); + } } - TCellArray output_new_face2ids; -// { -// TPTOIDS::const_iterator new_face2face_iter = output_newid2face.begin(); -// for(;new_face2face_iter!=output_newid2face.end();new_face2face_iter++){ - -// vtkIdType new_faceId = new_face2face_iter->first; -// TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin(); -// vtkIdType fId0 = *new_faces_iter; -// TCellArray::const_iterator pIds0_iter = f2points.find(fId0); -// TCell pIds0 = pIds0_iter->second; -// TCell &output = output_new_face2ids[new_faceId]; -// new_faces_iter++; -// if(new_face2face_iter->second.size() > 2 ){} -// for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){ - -// vtkIdType faceId = *new_faces_iter; -// // find how much nodes in face (f2points) -// TCellArray::const_iterator pIds_iter = f2points.find(faceId); -// TCell pIds = pIds_iter->second; - -// GetSumm(pIds0,pIds,output); -// pIds0 = output; - -// } // end new_faces_iter - -// } // new_face2face_iter -// } + // To sort the planar set of the points accrding to the angle { - TCellArray::const_iterator new_face2face_iter = output_newid2face_v2.begin(); - for(;new_face2face_iter!=output_newid2face_v2.end();new_face2face_iter++){ - - vtkIdType new_faceId = new_face2face_iter->first; - TCell::const_iterator new_faces_iter = (new_face2face_iter->second).begin(); - vtkIdType fId0 = *new_faces_iter; - TCellArray::const_iterator pIds0_iter = f2points.find(fId0); - TCell pIds0 = pIds0_iter->second; - TCell &output = output_new_face2ids[new_faceId]; - new_faces_iter++; - if(new_face2face_iter->second.size() == 1 ){ - TCellArray::const_iterator pIds_iter = f2points.find(fId0); - TCell pIds = pIds_iter->second; - output = pIds; - continue; - } - for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){ - - vtkIdType faceId = *new_faces_iter; - // find how much nodes in face (f2points) - TCellArray::const_iterator pIds_iter = f2points.find(faceId); - TCell pIds = pIds_iter->second; - - GetSumm(pIds0,pIds,output); - pIds0 = output; - - } // end new_faces_iter - - } // new_face2face_iter + typedef std::map 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)="<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(MYDEBUG) { - cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<GetPoint(anId,aPntCoord); + double aDist = vtkPlane::Evaluate(aNormal,anOrigin,aPntCoord); + if(DEBUG_TRIA_EXECUTE) cout<<"\t\taPntId = "<InsertNextCell(VTK_POLYGON,aNbPoints,&aConnectivities[0]); + if(theStoreMapping) + VTKViewer_GeometryFilter::InsertId( theCellId, VTK_POLYGON, theVTK2ObjIds, theDimension2VTK2ObjIds ); + theOutputCD->CopyData(thInputCD,theCellId,aNewCellId); + } + } + + if(DEBUG_TRIA_EXECUTE) cout<<"\tTriangulator - Ok\n"; + + return true; +} + + +//---------------------------------------------------------------------------- +VTKViewer_OrderedTriangulator +::VTKViewer_OrderedTriangulator(): + myTriangulator(vtkOrderedTriangulator::New()), + myBoundaryTris(vtkCellArray::New()), + myTriangle(vtkTriangle::New()) +{ + myBoundaryTris->Allocate(VTK_CELL_SIZE); + myTriangulator->PreSortedOff(); +} + + +//---------------------------------------------------------------------------- +VTKViewer_OrderedTriangulator +::~VTKViewer_OrderedTriangulator() +{ + myTriangle->Delete(); + myBoundaryTris->Delete(); + myTriangulator->Delete(); } + + +//---------------------------------------------------------------------------- +vtkPoints* +VTKViewer_OrderedTriangulator +::InitPoints(vtkUnstructuredGrid *theInput, + vtkIdType theCellId) +{ + 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); + } + + myTriangulator->Triangulate(); + myTriangulator->AddTriangles(myBoundaryTris); + } + + return aPoints; +} + + +//---------------------------------------------------------------------------- +vtkIdType +VTKViewer_OrderedTriangulator +::GetNumFaces() +{ + return myBoundaryTris->GetNumberOfCells(); +} + + +//---------------------------------------------------------------------------- +vtkCell* +VTKViewer_OrderedTriangulator +::GetFace(vtkIdType theFaceId) +{ + vtkIdType aNumCells = myBoundaryTris->GetNumberOfCells(); + if ( theFaceId < 0 || theFaceId >= aNumCells ) + return NULL; + + vtkIdType *aCells = myBoundaryTris->GetPointer(); + + // 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; +} + + +//---------------------------------------------------------------------------- +VTKViewer_DelaunayTriangulator +::VTKViewer_DelaunayTriangulator(): + myUnstructuredGrid(vtkUnstructuredGrid::New()), + myGeometryFilter(vtkGeometryFilter::New()), + myDelaunay3D(vtkDelaunay3D::New()), + myPolyData(NULL) +{ + myUnstructuredGrid->Initialize(); + myUnstructuredGrid->Allocate(); + myUnstructuredGrid->SetPoints(myPoints); + + myDelaunay3D->SetInputData(myUnstructuredGrid); + myGeometryFilter->SetInputConnection(myDelaunay3D->GetOutputPort()); + myPolyData = myGeometryFilter->GetOutput(); +} + + +//---------------------------------------------------------------------------- +VTKViewer_DelaunayTriangulator +::~VTKViewer_DelaunayTriangulator() +{ + myUnstructuredGrid->Delete(); + myGeometryFilter->Delete(); + myDelaunay3D->Delete(); +} + + +//---------------------------------------------------------------------------- +vtkPoints* +VTKViewer_DelaunayTriangulator +::InitPoints(vtkUnstructuredGrid *theInput, + vtkIdType theCellId) +{ + vtkPoints* aPoints = VTKViewer_Triangulator::InitPoints(theInput, theCellId); + + myPoints->Modified(); + myUnstructuredGrid->Modified(); + myGeometryFilter->Update(); + + return aPoints; +} + + +//---------------------------------------------------------------------------- +vtkIdType +VTKViewer_DelaunayTriangulator +::GetNumFaces() +{ + return myPolyData->GetNumberOfCells(); +} + + +//---------------------------------------------------------------------------- +vtkCell* +VTKViewer_DelaunayTriangulator +::GetFace(vtkIdType theFaceId) +{ + return myPolyData->GetCell(theFaceId); }