From bd1675247d4845dfcfbd108c651afc6b10bc37d4 Mon Sep 17 00:00:00 2001 From: apo Date: Mon, 20 Mar 2006 09:18:21 +0000 Subject: [PATCH] Fix for Bug IPAL11800 Fatal SIGSEGV on create presentation on cells for external file from PAL11736 Bug IPAL11667 SIGSEGV after trying to display mesh from "polyedres.med" in Mesh and Post-Pro modules. --- src/VTKViewer/VTKViewer_ConvexTool.cxx | 1284 ++++++++------------ src/VTKViewer/VTKViewer_ConvexTool.h | 196 ++- src/VTKViewer/VTKViewer_GeometryFilter.cxx | 344 +----- 3 files changed, 691 insertions(+), 1133 deletions(-) diff --git a/src/VTKViewer/VTKViewer_ConvexTool.cxx b/src/VTKViewer/VTKViewer_ConvexTool.cxx index 679e48428..464b27ce0 100644 --- a/src/VTKViewer/VTKViewer_ConvexTool.cxx +++ b/src/VTKViewer/VTKViewer_ConvexTool.cxx @@ -29,862 +29,552 @@ #include "VTKViewer_ConvexTool.h" -#include -#include -#include -#include -#include -#include - #include -#include -#include -#include - -typedef vtkUnstructuredGrid TInput; +#include -typedef std::set TUIDS; // unique ids -typedef std::map TPTOIDS; // id points -> unique ids - -namespace CONVEX_TOOL -{ - // all pairs - typedef std::pair TPair; - typedef std::set TSet; - - void - WriteToFile(vtkUnstructuredGrid* theDataSet, const std::string& theFileName) - { - vtkUnstructuredGridWriter* aWriter = vtkUnstructuredGridWriter::New(); - aWriter->SetFileName(theFileName.c_str()); - aWriter->SetInput(theDataSet); - aWriter->Write(); - aWriter->Delete(); - } +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include -static float FACE_ANGLE_TOLERANCE=1.5; -#define EPS 1.0e-38 -#define EPS_T 1.0e-3 - -#ifdef _DEBUG_ - static int MYDEBUG = 0; - static int MYDEBUG_REMOVE = 0; -#else - static int MYDEBUG = 0; - static int MYDEBUG_REMOVE = 0; -#endif - -/*! \fn static void GetCenter(TInput* theGrid,TCell theptIds,float *center) - * \brief Calculation of geometry center. - * \param theGrid - TInput cell. - * \param theptIds - point ids. - * \retval center - output array[3] with coordinates of geometry center. - */ -static void GetCenter(vtkPoints* thePoints,float center[3]) +//---------------------------------------------------------------------------- +namespace { - float p[3]; - center[0] = center[1] = center[2] = 0.0; + typedef std::vector TConnectivities; - int numPts = thePoints->GetNumberOfPoints(); - if (numPts == 0) return; - - // get the center of the cell - for (int i = 0; i < numPts; i++) + struct TPolygon { - thePoints->GetPoint(i, p); - for (int j = 0; j < 3; j++) + TConnectivities myConnectivities; + float myOrigin[3]; + float myNormal[3]; + TPolygon(const TConnectivities& theConnectivities, + float theOrigin[3], + float theNormal[3]): + myConnectivities(theConnectivities) { - center[j] += p[j]; + myOrigin[0] = theOrigin[0]; + myOrigin[1] = theOrigin[1]; + myOrigin[2] = theOrigin[2]; + + myNormal[0] = theNormal[0]; + myNormal[1] = theNormal[1]; + myNormal[2] = theNormal[2]; } - } - for (int j = 0; j < 3; j++) - { - center[j] /= numPts; - } + }; + + typedef std::vector TPolygons; } -/*! \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) +//---------------------------------------------------------------------------- +VTKViewer_Triangulator +::VTKViewer_Triangulator(): + myInput(NULL), + myCellId(-1), + myShowInside(-1), + myAllVisible(-1), + myCellsVisibility(NULL), + myCellIds(vtkIdList::New()) +{} + + +//---------------------------------------------------------------------------- +VTKViewer_Triangulator +::~VTKViewer_Triangulator() { - 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; - } + myCellIds->Delete(); } -/*! \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) -{ - TCellArray::const_iterator f2pIter = f2points.begin(); - 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 ; i& theVTK2ObjIds, + bool theIsCheckConvex) +{ + myInput = theInput; + myCellId = theCellId; + myShowInside = theShowInside; + myAllVisible = theAllVisible; + myCellsVisibility = theCellsVisibility; - id2 = face_points[i]; + vtkPoints *aPoints = InitPoints(); + vtkIdType aNumPts = GetNbOfPoints(); + //cout<<"Triangulator - aNumPts = "<second.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(aNumPts == 0) + return true; - if(MYDEBUG){ - cout << "fId[" << faceId <<"]: "; - std::copy(output_faces.begin(), output_faces.end(), std::ostream_iterator(cout, " ")); - cout << endl; + // To calculate the bary center of the cell + float aCellCenter[3] = {0.0, 0.0, 0.0}; + { + float aPntCoord[3]; + for (int aPntId = 0; aPntId < aNumPts; aPntId++) { + aPoints->GetPoint(GetPointId(aPntId),aPntCoord); + //cout<<"\taPntId = "<"; - 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); + float aCellLength = GetCellLength(); + int aNumFaces = GetNumFaces(); - float b1 = vtkMath::Norm(vec_b1); + static float EPS = 1.0E-2; + float aDistEps = aCellLength * EPS; + //cout<<"\taCellLength = "<=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] << ") "; - } - if(status) cout << "angle="< TPointIds; + TPointIds anInitialPointIds; + for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++) + anInitialPointIds.insert(GetPointId(aPntId)); - return status; -} - -/*! \fn void GetAllFacesOnOnePlane( TPTOIDS theFaces, vtkIdType faceId,TUIDS &new_faces, TCell &new_faces_v2 ) - * \brief Calculate faces which on one plane. - * \param theFaces - - * \param faceId - - * \param new_faces - - * \param new_faces_v2 - - */ -void GetAllFacesOnOnePlane( TPTOIDS theFaces, vtkIdType faceId, - TUIDS &new_faces, TCell &new_faces_v2 ) -{ - if (new_faces.find(faceId) != new_faces.end()) return; + // To initialize set of points by face that belong to the cell and backward + typedef std::set TFace2Visibility; + TFace2Visibility aFace2Visibility; - 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 - } - } - return; -} + typedef std::set TFace2PointIds; + TFace2PointIds aFace2PointIds; -/*! \fn void GetSumm(TCell v1,TCell v2,TCell &output) - * \brief Gluing two faces (gluing points ids) - * \param v1 - first face - * \param v2 - second face - * \param output - output face. - */ -void GetSumm(TCell v1,TCell v2,TCell &output) -{ - output.clear(); - - if(MYDEBUG) cout << "========================================="<(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; - - TCell::iterator v1_iter = v1.begin(); - - for(;v1_iter!=v1.end();v1_iter++){ - - vtkIdType curr_id = *v1_iter; + for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) { + vtkCell* aFace = GetFace(aFaceId); - output.push_back(curr_id); - - 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(theCellId, aFace, myCellIds); + if((!myAllVisible && !myCellsVisibility[myCellIds->GetId(0)]) || + myCellIds->GetNumberOfIds() <= 0 || myShowInside) + { + 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; -} - -static void GetAndRemoveIdsOnOneLine(vtkPoints* points, - TUIDS input_points_ids, - TUIDS input_two_points_ids, - TUIDS& out_two_points_ids, - TUIDS& removed_points_ids){ - if (MYDEBUG_REMOVE) cout << EPS <GetPoint(current_points_ids[0],P[0]); - points->GetPoint(current_points_ids[1],P[1]); - TUIDS::iterator aPointsIter = input_points_ids.begin(); - for(;aPointsIter!=input_points_ids.end();aPointsIter++){ - if(iscurrent_points_changed){ - points->GetPoint(current_points_ids[0],P[0]); - points->GetPoint(current_points_ids[1],P[1]); - iscurrent_points_changed = false; - if (MYDEBUG_REMOVE) - cout << " " << current_points_ids[0] << " " << current_points_ids[1] << endl; - } - // check: is point on line input_two_points_ids - points->GetPoint(*aPointsIter,P[2]); - if (MYDEBUG_REMOVE) { - cout << "\t" << current_points_ids[0] << ":"< coeff[0][0] = (x-x1), coeff[0][1] = x2-x1 - // y-y1=(y2-y1)*t -> coeff[1][0] = (y-y1), coeff[1][1] = y2-y1 - // z-z1=(z2-z1)*t -> coeff[2][0] = (z-z1), coeff[2][1] = z2-z1 - float coeff[3][2]; - for(int i=0;i<3;i++){ - coeff[i][0] = P[2][i]-P[0][i]; - coeff[i][1] = P[1][i]-P[0][i]; - } - bool isok_coord[3]; - bool isok = true; - float t[3]; - for(int i=0;i<3;i++){ - isok_coord[i] = false; - if( fabs(coeff[i][0]) <= EPS && fabs(coeff[i][1]) <= EPS) { - isok_coord[i] = true; - continue; - } - if( fabs(coeff[i][1]) <= EPS && fabs(coeff[i][0]) > EPS) {isok = false;t[i]=1.0/EPS;break;} - t[i] = (coeff[i][0])/(coeff[i][1]); - } - for(int i=0;i<3;i++) - if (MYDEBUG_REMOVE) - cout << __LINE__ << " " - << coeff[i][0] << ","<t[0],t[1],t[2] - // anilize bounds of line - for(int i=0;i<3;i++){ - for(int j=0;j<3;j++) - if(!isok_coord[j]) param[i] = (P[i][j]-P[0][j])/(P[1][j]-P[0][j]); - } - if (MYDEBUG_REMOVE) cout << "Params: " << param[0] << "," << param[1] << "," << param[2] << endl; - vtkIdType imax,imin; - float min,max; - for(vtkIdType i=0;i<3;i++) - if(!isok_coord[i]){ - min = param[0];imin=0; - max = param[0];imax=0; - break; - } - for(vtkIdType i=0;i<3;i++){ - if(min > param[i]) {min = param[i]; imin=i;} - if(max < param[i]) {max = param[i]; imax=i;} - } - if (MYDEBUG_REMOVE) - cout << "\t min="<"< RemoveAllUnneededPoints(vtkConvexPointSet* convex){ - vtkSmartPointer out = vtkConvexPointSet::New(); - - TUIDS two_points,input_points,out_two_points_ids,removed_points_ids,loc_removed_points_ids; - vtkIdList* aPointIds = convex->GetPointIds(); - int numIds = aPointIds->GetNumberOfIds(); - if (numIds<2) return out; - TSet good_point_ids; - TSet aLists[numIds-2]; - for(int i=0;iGetId(aFirId) << "," << aSecId <<":"<GetId(aSecId)<< " --- "; - for(vtkIdType k=aSecId+1;kGetId(k) << ","; - } - if (MYDEBUG_REMOVE) { - cout << endl; - cout << "\t"; - for(TUIDS::iterator aDelIter = loc_removed_points_ids.begin();aDelIter!=loc_removed_points_ids.end();aDelIter++) - cout << *aDelIter<<","; - cout << endl; - } - GetAndRemoveIdsOnOneLine(convex->Points, - input_points, - two_points, - out_two_points_ids, - loc_removed_points_ids); - TUIDS::iterator aOutIter = out_two_points_ids.begin(); - vtkIdType aFirst=*aOutIter;aOutIter++;vtkIdType aSecond=*aOutIter; - TPair aPair(aFirst,aSecond); - good_point_ids.insert(aPair); - if (MYDEBUG_REMOVE){ - cout << "Output: "; - TUIDS::iterator aIter = out_two_points_ids.begin(); - for(;aIter!=out_two_points_ids.end();aIter++) - cout << *aIter << ","; - cout << " --- "; - } - TUIDS::iterator aDelIter = loc_removed_points_ids.begin(); - for(;aDelIter!=loc_removed_points_ids.end();aDelIter++){ - removed_points_ids.insert(*aDelIter); - if (MYDEBUG_REMOVE) cout << *aDelIter << ","; - } - if (MYDEBUG_REMOVE) cout << endl; - } - } - if (MYDEBUG_REMOVE) { - cout << "============ Resultat ================" <Points->SetNumberOfPoints(result_ids.size()); - out->PointIds->SetNumberOfIds(result_ids.size()); - int aId=0; - if(MYDEBUG_REMOVE) cout << "Result out:"; - for(TUIDS::iterator aIter=result_ids.begin();aIter!=result_ids.end();aIter++,aId++){ - float P[3]; - convex->Points->GetPoint(*aIter,P); - out->Points->SetPoint(aId,P); - out->PointIds->SetId(aId,aPointIds->GetId(*aIter)); - if (MYDEBUG_REMOVE) cout << *aIter << ":" << aPointIds->GetId(*aIter) << " , "; - } - if(MYDEBUG_REMOVE) cout << endl; - out->Modified(); - out->Initialize(); - - return out; -} -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_in = static_cast(theGrid->GetCell(cellId))){ - vtkSmartPointer convex = RemoveAllUnneededPoints(convex_in); - TCellArray f2points; - float convex_center[3]; // convex center point coorinat - int aNbFaces = convex->GetNumberOfFaces(); - int numPts = convex->GetNumberOfPoints(); - if(MYDEBUG_REMOVE) cout << "aNbFaces="<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]; - } + ::TPolygons aPolygons; - if (vtkMath::Determinant3x3(v_a,v_b,v_convex2face) < 0){ - ReverseIds(aIds); - } + for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) { + if(aFace2Visibility.find(aFaceId) == aFace2Visibility.end()) + continue; - for(i=0;i<(int)aIds.size();i++){ - TUIDS &acell = p2faces[aIds[i]]; - acell.insert(faceId); - } - - f2points[faceId] = aIds; + vtkCell* aFace = GetFace(aFaceId); - } - - TPTOIDS face2face; - GetFriends(p2faces,f2points,face2face); + 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(); + //cout<<"\taFaceId = "<GetPoint(aNewPts[0],aCoord[0]); + aPoints->GetPoint(aNewPts[1],aCoord[1]); + aPoints->GetPoint(aNewPts[2],aCoord[2]); - TPTOIDS face2points; + // To calculate plane normal + float aVector01[3] = { aCoord[1][0] - aCoord[0][0], + aCoord[1][1] - aCoord[0][1], + aCoord[1][2] - aCoord[0][2] }; - // copy TCellArray::f2points to TPTOIDS::face2points - for(TCellArray::iterator f2points_iter=f2points.begin(); - f2points_iter!=f2points.end(); - f2points_iter++){ - - TUIDS tmp; - for(TCell::iterator points_iter=(f2points_iter->second).begin(); - points_iter!=(f2points_iter->second).end(); - points_iter++) - tmp.insert(*points_iter); - - face2points[f2points_iter->first] = tmp; - } // end copy - + float aVector02[3] = { aCoord[2][0] - aCoord[0][0], + aCoord[2][1] - aCoord[0][1], + aCoord[2][2] - aCoord[0][2] }; - TPTOIDS new_face2faces; // which connected and in one plane + // To calculate the normal for the triangle + float aNormal[3]; + vtkMath::Cross(aVector02,aVector01,aNormal); - 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 + vtkMath::Normalize(aNormal); - if(MYDEBUG) + // To calculate what points belong to the plane + // To calculate bounds of the point set + float aCenter[3] = {0.0, 0.0, 0.0}; { - 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; + TPointIds::const_iterator anIter = anInitialPointIds.begin(); + TPointIds::const_iterator anEndIter = anInitialPointIds.end(); + for(; anIter != anEndIter; anIter++){ + float aPntCoord[3]; + vtkIdType aPntId = *anIter; + aPoints->GetPoint(aPntId,aPntCoord); + float aDist = vtkPlane::DistanceToPlane(aPntCoord,aNormal,aCoord[0]); + //cout<<"\t\taPntId = "< 0){ + aNormal[0] = -aNormal[0]; + aNormal[1] = -aNormal[1]; + aNormal[2] = -aNormal[2]; } - TPTOIDS output_newid2face; - TCellArray output_newid2face_v2; + // To calculate the primary direction for point set + float aVector0[3] = { aCoord[0][0] - aCenter[0], + aCoord[0][1] - aCenter[1], + aCoord[0][2] - aCenter[2] }; + vtkMath::Normalize(aVector0); + + //cout<<"\t\taCenter = {"<first) != already_in.end()) - continue; - if(new_face2face_iter->second.size() > 1) - continue; + 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())); - 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); - } + if(anIntersection == anIds){ + aRemoveFace2PointIds.insert(anIds); } - k++; } } - 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++){ + typedef std::map TSortedPointIds; + TSortedPointIds aSortedPointIds; + + TPointIds::const_iterator anIter = aPointIds.begin(); + TPointIds::const_iterator anEndIter = aPointIds.end(); + for(; anIter != anEndIter; anIter++){ + float aPntCoord[3]; + vtkIdType aPntId = *anIter; + aPoints->GetPoint(aPntId,aPntCoord); + float aVector[3] = { aPntCoord[0] - aCenter[0], + aPntCoord[1] - aCenter[1], + aPntCoord[2] - aCenter[2] }; + vtkMath::Normalize(aVector); - 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; + float aCross[3]; + vtkMath::Cross(aVector,aVector0,aCross); + bool aGreaterThanPi = vtkMath::Dot(aCross,aNormal) < 0; + float aCosinus = vtkMath::Dot(aVector,aVector0); + if(aCosinus > 1.0) + aCosinus = 1.0; + if(aCosinus < -1.0) + aCosinus = -1.0; + static float a2Pi = 2.0 * vtkMath::Pi(); + float anAngle = acos(aCosinus); + //cout<<"\t\t\taPntId = "<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; + aSortedPointIds[anAngle] = aPntId; + } - } // end new_faces_iter - - } // new_face2face_iter + 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)); + } } - - if(MYDEBUG) { - cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<GetPoint(anId,aPntCoord); + float aDist = vtkPlane::Evaluate(aNormal,anOrigin,aPntCoord); + //cout<<"\t\taPntId = "<InsertNextCell(VTK_POLYGON,aNbPoints,&aConnectivities[0]); + if(theStoreMapping) + theVTK2ObjIds.push_back(theCellId); + theOutputCD->CopyData(thInputCD,theCellId,aNewCellId); + } + } + + //cout<<"\tTriangulator - Ok\n"; + return true; +} + +//---------------------------------------------------------------------------- +VTKViewer_OrderedTriangulator +::VTKViewer_OrderedTriangulator(): + myCell(vtkGenericCell::New()) +{} + + +//---------------------------------------------------------------------------- +VTKViewer_OrderedTriangulator +::~VTKViewer_OrderedTriangulator() +{ + myCell->Delete(); +} + +//---------------------------------------------------------------------------- +vtkPoints* +VTKViewer_OrderedTriangulator +::InitPoints() +{ + myInput->GetCell(myCellId,myCell); + return myInput->GetPoints(); +} + +vtkIdType +VTKViewer_OrderedTriangulator +::GetNbOfPoints() +{ + return myCell->GetNumberOfPoints(); +} + +vtkIdType +VTKViewer_OrderedTriangulator +::GetPointId(vtkIdType thePointId) +{ + return myCell->GetPointId(thePointId); } + +float +VTKViewer_OrderedTriangulator +::GetCellLength() +{ + return sqrt(myCell->GetLength2()); +} + +vtkIdType +VTKViewer_OrderedTriangulator +::GetNumFaces() +{ + return myCell->GetNumberOfFaces(); +} + +vtkCell* +VTKViewer_OrderedTriangulator +::GetFace(vtkIdType theFaceId) +{ + return myCell->GetFace(theFaceId); +} + +void +VTKViewer_OrderedTriangulator +::GetCellNeighbors(vtkIdType theCellId, + vtkCell* theFace, + vtkIdList* theCellIds) +{ + vtkIdList *anIdList = theFace->PointIds; + myInput->GetCellNeighbors(theCellId, anIdList, theCellIds); +} + +vtkIdType +VTKViewer_OrderedTriangulator +::GetConnectivity(vtkIdType thePntId) +{ + return thePntId; +} + +//---------------------------------------------------------------------------- +VTKViewer_DelaunayTriangulator +::VTKViewer_DelaunayTriangulator(): + myUnstructuredGrid(vtkUnstructuredGrid::New()), + myGeometryFilter(vtkGeometryFilter::New()), + myDelaunay3D(vtkDelaunay3D::New()), + myFaceIds(vtkIdList::New()), + myPoints(vtkPoints::New()), + myPolyData(NULL), + myPointIds(NULL) +{ + myDelaunay3D->SetInput(myUnstructuredGrid); + myGeometryFilter->SetInput(myDelaunay3D->GetOutput()); +} + + +//---------------------------------------------------------------------------- +VTKViewer_DelaunayTriangulator +::~VTKViewer_DelaunayTriangulator() +{ + myUnstructuredGrid->Delete(); + myGeometryFilter->Delete(); + myDelaunay3D->Delete(); + myFaceIds->Delete(); + myPoints->Delete(); +} + + +//---------------------------------------------------------------------------- +vtkPoints* +VTKViewer_DelaunayTriangulator +::InitPoints() +{ + myUnstructuredGrid->Initialize(); + myUnstructuredGrid->Allocate(); + myUnstructuredGrid->SetPoints(myPoints); + + vtkIdType aNumPts; + myInput->GetCellPoints(myCellId,aNumPts,myPointIds); + { + float 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); + } + } + + myGeometryFilter->Update(); + myPolyData = myGeometryFilter->GetOutput(); + + return myPoints; +} + +vtkIdType +VTKViewer_DelaunayTriangulator +::GetNbOfPoints() +{ + return myPoints->GetNumberOfPoints(); +} + +vtkIdType +VTKViewer_DelaunayTriangulator +::GetPointId(vtkIdType thePointId) +{ + return thePointId; +} + +float +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]; } diff --git a/src/VTKViewer/VTKViewer_ConvexTool.h b/src/VTKViewer/VTKViewer_ConvexTool.h index ed303af03..f56a502b4 100644 --- a/src/VTKViewer/VTKViewer_ConvexTool.h +++ b/src/VTKViewer/VTKViewer_ConvexTool.h @@ -21,28 +21,186 @@ #ifndef _VTKViewer_ConvexTool_H #define _VTKViewer_ConvexTool_H -#include +#include "VTKViewer.h" + #include -#include -typedef std::vector TCell; // ptsIds -typedef std::map TCellArray; // CellId, TCell +#include + +class vtkUnstructuredGrid; +class vtkGeometryFilter; +class vtkGenericCell; +class vtkDelaunay3D; +class vtkPolyData; +class vtkCellData; +class vtkPoints; +class vtkIdList; +class vtkCell; + +//---------------------------------------------------------------------------- +class VTKVIEWER_EXPORT VTKViewer_Triangulator +{ + public: + VTKViewer_Triangulator(); + + ~VTKViewer_Triangulator(); + + bool + Execute(vtkUnstructuredGrid *theInput, + vtkCellData* thInputCD, + vtkIdType theCellId, + int theShowInside, + int theAllVisible, + const char* theCellsVisibility, + vtkPolyData *theOutput, + vtkCellData* theOutputCD, + int theStoreMapping, + std::vector& theVTK2ObjIds, + bool theIsCheckConvex); + + protected: + vtkIdList* myCellIds; + + vtkUnstructuredGrid *myInput; + vtkIdType myCellId; + int myShowInside; + int myAllVisible; + const char* myCellsVisibility; + + virtual + vtkPoints* + InitPoints() = 0; + + virtual + vtkIdType + GetNbOfPoints() = 0; + + virtual + vtkIdType + GetPointId(vtkIdType thePointId) = 0; -/*! This package \namespace CONVEX_TOOL used for: - * calculation of VTK_POLYGON cell array from VTK_TRIANGLE (triangulation) - * of VTK_CONVEX_POINT_SET cell type. - */ -namespace CONVEX_TOOL + virtual + float + GetCellLength() = 0; + + virtual + vtkIdType + GetNumFaces() = 0; + + virtual + vtkCell* + GetFace(vtkIdType theFaceId) = 0; + + virtual + void + GetCellNeighbors(vtkIdType theCellId, + vtkCell* theFace, + vtkIdList* theCellIds) = 0; + + virtual + vtkIdType + GetConnectivity(vtkIdType thePntId) = 0; +}; + + +//---------------------------------------------------------------------------- +class VTKVIEWER_EXPORT VTKViewer_OrderedTriangulator : public VTKViewer_Triangulator { - /*! \fn void CONVEX_TOOL::GetPolygonalFaces(vtkUnstructuredGrid* theCell,int cellId,TCellArray &outputCellArray) - * \brief Main function. - * \param theCell - vtkUnstructuredGrid cell pointer - * \param cellId - id of cell type VTK_CONVEX_POINT_SET - * \retval outputCellArray - output array with new cells types VTK_POLYGON - */ - void - WriteToFile(vtkUnstructuredGrid* theDataSet, const std::string& theFileName); - void GetPolygonalFaces(vtkUnstructuredGrid* theCell,int cellId,TCellArray &outputCellArray); -} + public: + + VTKViewer_OrderedTriangulator(); + + ~VTKViewer_OrderedTriangulator(); + + protected: + vtkGenericCell *myCell; + + virtual + vtkPoints* + InitPoints(); + + virtual + vtkIdType + GetNbOfPoints(); + + vtkIdType + GetPointId(vtkIdType thePointId); + + virtual + float + GetCellLength(); + + virtual + vtkIdType + GetNumFaces(); + + virtual + vtkCell* + GetFace(vtkIdType theFaceId); + + virtual + void + GetCellNeighbors(vtkIdType theCellId, + vtkCell* theFace, + vtkIdList* theCellIds); + + virtual + vtkIdType + GetConnectivity(vtkIdType thePntId); +}; + + +//---------------------------------------------------------------------------- +class VTKVIEWER_EXPORT VTKViewer_DelaunayTriangulator : public VTKViewer_Triangulator +{ + public: + + VTKViewer_DelaunayTriangulator(); + + ~VTKViewer_DelaunayTriangulator(); + + protected: + vtkUnstructuredGrid* myUnstructuredGrid; + vtkGeometryFilter* myGeometryFilter; + vtkDelaunay3D* myDelaunay3D; + vtkPolyData* myPolyData; + vtkIdType *myPointIds; + vtkIdList* myFaceIds; + vtkPoints* myPoints; + + virtual + vtkPoints* + InitPoints(); + + virtual + vtkIdType + GetNbOfPoints(); + + vtkIdType + GetPointId(vtkIdType thePointId); + + virtual + float + GetCellLength(); + + virtual + vtkIdType + GetNumFaces(); + + virtual + vtkCell* + GetFace(vtkIdType theFaceId); + + virtual + void + GetCellNeighbors(vtkIdType theCellId, + vtkCell* theFace, + vtkIdList* theCellIds); + + virtual + vtkIdType + GetConnectivity(vtkIdType thePntId); +}; + #endif // _VTKViewer_ConvexTool_H diff --git a/src/VTKViewer/VTKViewer_GeometryFilter.cxx b/src/VTKViewer/VTKViewer_GeometryFilter.cxx index df096cf7a..83d6093ce 100755 --- a/src/VTKViewer/VTKViewer_GeometryFilter.cxx +++ b/src/VTKViewer/VTKViewer_GeometryFilter.cxx @@ -27,6 +27,7 @@ // $Header$ #include "VTKViewer_GeometryFilter.h" +#include "VTKViewer_ConvexTool.h" #include #include @@ -45,11 +46,6 @@ #include #include -#include -#include -#include -#include - #include #include #include @@ -62,7 +58,7 @@ #endif #endif -#define USE_ROBUST_TRIANGULATION +//#define USE_ROBUST_TRIANGULATION //---------------------------------------------------------------------------- vtkCxxRevisionMacro(VTKViewer_GeometryFilter, "$Revision$"); @@ -128,21 +124,12 @@ VTKViewer_GeometryFilter vtkPolyData *output = this->GetOutput(); vtkPointData *outputPD = output->GetPointData(); -#ifdef USE_ROBUST_TRIANGULATION - vtkUnstructuredGrid* anUnstructuredGrid = vtkUnstructuredGrid::New(); - vtkPoints* aDelaunayPoints = vtkPoints::New(); - - vtkDelaunay3D* aDelaunay3D = vtkDelaunay3D::New(); - aDelaunay3D->SetInput(anUnstructuredGrid); - - vtkGeometryFilter* aGeometryFilter = vtkGeometryFilter::New(); - aGeometryFilter->SetInput(aDelaunay3D->GetOutput()); -#endif + VTKViewer_OrderedTriangulator anOrderedTriangulator; + VTKViewer_DelaunayTriangulator aDelaunayTriangulator; vtkCellData *outputCD = output->GetCellData(); vtkGenericCell *cell = vtkGenericCell::New(); - vtkIdList *cellIds = vtkIdList::New(); vtkIdList *faceIds = vtkIdList::New(); @@ -306,299 +293,30 @@ VTKViewer_GeometryFilter break; case VTK_CONVEX_POINT_SET: { - //cout<<"cellId = "<Initialize(); - anUnstructuredGrid->Allocate(); - anUnstructuredGrid->SetPoints(aDelaunayPoints); - - vtkIdType *aPts; - input->GetCellPoints(cellId,aNumPts,aPts); - { - float aPntCoord[3]; - aDelaunayPoints->SetNumberOfPoints(aNumPts); - vtkPoints *anInputPoints = input->GetPoints(); - for (int aPntId = 0; aPntId < aNumPts; aPntId++) { - anInputPoints->GetPoint(aPts[aPntId],aPntCoord); - aDelaunayPoints->SetPoint(aPntId,aPntCoord); - } - } -#else - input->GetCell(cellId,cell); - aPoints = input->GetPoints(); - aNumPts = cell->GetNumberOfPoints(); -#endif - // To calculate the bary center of the cell - float aCellCenter[3] = {0.0, 0.0, 0.0}; - { - float aPntCoord[3]; - for (int aPntId = 0; aPntId < aNumPts; aPntId++) { -#ifdef USE_ROBUST_TRIANGULATION - aPoints->GetPoint(aPntId,aPntCoord); -#else - aPoints->GetPoint(cell->GetPointId(aPntId),aPntCoord); -#endif - //cout<<"\t\taPntId = "<Update(); - vtkPolyData* aPolyData = aGeometryFilter->GetOutput(); - - float aCellLength = aPolyData->GetLength(); - int aNumFaces = aPolyData->GetNumberOfCells(); -#else - float aCellLength = sqrt(cell->GetLength2()); - int aNumFaces = cell->GetNumberOfFaces(); -#endif - - static float EPS = 1.0E-5; - float aDistEps = aCellLength * EPS; - - // To initialize set of points that belong to the cell - typedef std::set TPointIds; - TPointIds anInitialPointIds; - for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++){ -#ifdef USE_ROBUST_TRIANGULATION - anInitialPointIds.insert(aPntId); -#else - anInitialPointIds.insert(cell->GetPointId(aPntId)); -#endif - } - - // To initialize set of points by face that belong to the cell and backward - typedef std::set TFace2Visibility; - TFace2Visibility aFace2Visibility; - - typedef std::set TFace2PointIds; - TFace2PointIds aFace2PointIds; - - for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) { -#ifdef USE_ROBUST_TRIANGULATION - vtkCell* aFace = aPolyData->GetCell(aFaceId); -#else - vtkCell* aFace = cell->GetFace(aFaceId); -#endif - vtkIdList *anIdList = aFace->PointIds; - aNewPts[0] = anIdList->GetId(0); - aNewPts[1] = anIdList->GetId(1); - aNewPts[2] = anIdList->GetId(2); - -#ifdef USE_ROBUST_TRIANGULATION - faceIds->Reset(); - faceIds->InsertNextId(aPts[aNewPts[0]]); - faceIds->InsertNextId(aPts[aNewPts[1]]); - faceIds->InsertNextId(aPts[aNewPts[2]]); - input->GetCellNeighbors(cellId, faceIds, cellIds); -#else - input->GetCellNeighbors(cellId, anIdList, cellIds); -#endif - if((!allVisible && !cellVis[cellIds->GetId(0)]) || - cellIds->GetNumberOfIds() <= 0 || - myShowInside) - { - TPointIds aPointIds; - aPointIds.insert(aNewPts[0]); - aPointIds.insert(aNewPts[1]); - aPointIds.insert(aNewPts[2]); - - aFace2PointIds.insert(aPointIds); - aFace2Visibility.insert(aFaceId); - } - } - - for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) { - if(aFace2Visibility.find(aFaceId) == aFace2Visibility.end()) - continue; - -#ifdef USE_ROBUST_TRIANGULATION - vtkCell* aFace = aPolyData->GetCell(aFaceId); -#else - vtkCell* aFace = cell->GetFace(aFaceId); -#endif - vtkIdList *anIdList = aFace->PointIds; - aNewPts[0] = anIdList->GetId(0); - aNewPts[1] = anIdList->GetId(1); - aNewPts[2] = 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(); - //cout<<"\taFaceId = "<GetPoint(aNewPts[0],aCoord[0]); - aPoints->GetPoint(aNewPts[1],aCoord[1]); - aPoints->GetPoint(aNewPts[2],aCoord[2]); - - // To calculate plane normal - float aVector01[3] = { aCoord[1][0] - aCoord[0][0], - aCoord[1][1] - aCoord[0][1], - aCoord[1][2] - aCoord[0][2] }; - - float aVector02[3] = { aCoord[2][0] - aCoord[0][0], - aCoord[2][1] - aCoord[0][1], - aCoord[2][2] - aCoord[0][2] }; + bool anIsOk = anOrderedTriangulator.Execute(input, + cd, + cellId, + myShowInside, + allVisible, + cellVis, + output, + outputCD, + myStoreMapping, + myVTK2ObjIds, + true); + if(!anIsOk) + aDelaunayTriangulator.Execute(input, + cd, + cellId, + myShowInside, + allVisible, + cellVis, + output, + outputCD, + myStoreMapping, + myVTK2ObjIds, + false); - float aCross21[3]; - vtkMath::Cross(aVector02,aVector01,aCross21); - - vtkMath::Normalize(aCross21); - - // To calculate what points belong to the plane - // To calculate bounds of the point set - float aCenter[3] = {0.0, 0.0, 0.0}; - { - TPointIds::const_iterator anIter = anInitialPointIds.begin(); - TPointIds::const_iterator anEndIter = anInitialPointIds.end(); - for(; anIter != anEndIter; anIter++){ - float aPntCoord[3]; - vtkIdType aPntId = *anIter; - aPoints->GetPoint(aPntId,aPntCoord); - float aDist = vtkPlane::DistanceToPlane(aPntCoord,aCross21,aCoord[0]); - //cout<<"\t\taPntId = "< 0){ - aCross21[0] = -aCross21[0]; - aCross21[1] = -aCross21[1]; - aCross21[2] = -aCross21[2]; - } - - vtkMath::Normalize(aVector0); - - //cout<<"\t\taCenter = {"< TSortedPointIds; - TSortedPointIds aSortedPointIds; - - TPointIds::const_iterator anIter = aPointIds.begin(); - TPointIds::const_iterator anEndIter = aPointIds.end(); - for(; anIter != anEndIter; anIter++){ - float aPntCoord[3]; - vtkIdType aPntId = *anIter; - aPoints->GetPoint(aPntId,aPntCoord); - float aVector[3] = { aPntCoord[0] - aCenter[0], - aPntCoord[1] - aCenter[1], - aPntCoord[2] - aCenter[2] }; - vtkMath::Normalize(aVector); - - float aCross[3]; - vtkMath::Cross(aVector,aVector0,aCross); - bool aGreaterThanPi = vtkMath::Dot(aCross,aCross21) < 0; - float aCosinus = vtkMath::Dot(aVector,aVector0); - if(aCosinus > 1.0) - aCosinus = 1.0; - if(aCosinus < -1.0) - aCosinus = -1.0; - static float a2Pi = 2.0 * vtkMath::Pi(); - float anAngle = acos(aCosinus); - //cout<<"\t\taPntId = "< aConnectivities(numFacePts); - TSortedPointIds::const_iterator anIter = aSortedPointIds.begin(); - TSortedPointIds::const_iterator anEndIter = aSortedPointIds.end(); - for(vtkIdType anId = 0; anIter != anEndIter; anIter++, anId++){ - vtkIdType aPntId = anIter->second; -#ifdef USE_ROBUST_TRIANGULATION - aConnectivities[anId] = aPts[aPntId]; -#else - aConnectivities[anId] = aPntId; -#endif - } - newCellId = output->InsertNextCell(aCellType,numFacePts,&aConnectivities[0]); - if(myStoreMapping) - myVTK2ObjIds.push_back(cellId); - outputCD->CopyData(cd,cellId,newCellId); - } - } - } - } - break; } case VTK_TETRA: { @@ -1029,14 +747,6 @@ VTKViewer_GeometryFilter vtkDebugMacro(<<"Extracted " << input->GetNumberOfPoints() << " points," << output->GetNumberOfCells() << " cells."); -#ifdef USE_ROBUST_TRIANGULATION - anUnstructuredGrid->Delete(); - aDelaunayPoints->Delete(); - - aDelaunay3D->Delete(); - aGeometryFilter->Delete(); -#endif - cell->Delete(); cellIds->Delete(); -- 2.39.2