1 // Copyright (C) 2005 OPEN CASCADE, CEA/DEN, EDF R&D, PRINCIPIA R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
8 // This library is distributed in the hope that it will be useful
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/
19 // This library is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22 // Lesser General Public License for more details.
24 // You should have received a copy of the GNU Lesser General Public
25 // License along with this library; if not, write to the Free Software
26 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28 // See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
30 #include "VTKViewer_ConvexTool.h"
35 #include <vtkUnstructuredGrid.h>
36 #include <vtkGeometryFilter.h>
37 #include <vtkDelaunay3D.h>
38 #include <vtkGenericCell.h>
39 #include <vtkPolyData.h>
40 #include <vtkCellData.h>
41 #include <vtkPoints.h>
42 #include <vtkIdList.h>
49 typedef std::vector<vtkIdType> TConnectivities;
53 TConnectivities myConnectivities;
54 vtkFloatingPointType myOrigin[3];
55 vtkFloatingPointType myNormal[3];
56 TPolygon(const TConnectivities& theConnectivities,
57 vtkFloatingPointType theOrigin[3],
58 vtkFloatingPointType theNormal[3]):
59 myConnectivities(theConnectivities)
61 myOrigin[0] = theOrigin[0];
62 myOrigin[1] = theOrigin[1];
63 myOrigin[2] = theOrigin[2];
65 myNormal[0] = theNormal[0];
66 myNormal[1] = theNormal[1];
67 myNormal[2] = theNormal[2];
71 typedef std::vector<TPolygon> TPolygons;
77 VTKViewer_Triangulator
78 ::VTKViewer_Triangulator():
83 myCellsVisibility(NULL),
84 myCellIds(vtkIdList::New())
91 VTKViewer_Triangulator
92 ::~VTKViewer_Triangulator()
100 VTKViewer_Triangulator
101 ::Execute(vtkUnstructuredGrid *theInput,
102 vtkCellData* thInputCD,
106 const char* theCellsVisibility,
107 vtkPolyData *theOutput,
108 vtkCellData* theOutputCD,
110 std::vector<vtkIdType>& theVTK2ObjIds,
111 bool theIsCheckConvex)
114 myCellId = theCellId;
115 myShowInside = theShowInside;
116 myAllVisible = theAllVisible;
117 myCellsVisibility = theCellsVisibility;
119 vtkPoints *aPoints = InitPoints();
120 vtkIdType aNumPts = GetNbOfPoints();
121 //cout<<"Triangulator - aNumPts = "<<aNumPts<<"\n";
126 // To calculate the bary center of the cell
127 vtkFloatingPointType aCellCenter[3] = {0.0, 0.0, 0.0};
129 vtkFloatingPointType aPntCoord[3];
130 for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
131 aPoints->GetPoint(GetPointId(aPntId),aPntCoord);
132 //cout<<"\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}\n";
133 aCellCenter[0] += aPntCoord[0];
134 aCellCenter[1] += aPntCoord[1];
135 aCellCenter[2] += aPntCoord[2];
137 aCellCenter[0] /= aNumPts;
138 aCellCenter[1] /= aNumPts;
139 aCellCenter[2] /= aNumPts;
142 vtkFloatingPointType aCellLength = GetCellLength();
143 int aNumFaces = GetNumFaces();
145 static vtkFloatingPointType EPS = 1.0E-2;
146 vtkFloatingPointType aDistEps = aCellLength * EPS;
147 //cout<<"\taCellLength = "<<aCellLength<<"; aDistEps = "<<aDistEps<<"\n";
149 // To initialize set of points that belong to the cell
150 typedef std::set<vtkIdType> TPointIds;
151 TPointIds anInitialPointIds;
152 for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++)
153 anInitialPointIds.insert(GetPointId(aPntId));
155 // To initialize set of points by face that belong to the cell and backward
156 typedef std::set<vtkIdType> TFace2Visibility;
157 TFace2Visibility aFace2Visibility;
159 typedef std::set<TPointIds> TFace2PointIds;
160 TFace2PointIds aFace2PointIds;
162 for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) {
163 vtkCell* aFace = GetFace(aFaceId);
165 GetCellNeighbors(theCellId, aFace, myCellIds);
166 if((!myAllVisible && !myCellsVisibility[myCellIds->GetId(0)]) ||
167 myCellIds->GetNumberOfIds() <= 0 || myShowInside)
170 vtkIdList *anIdList = aFace->PointIds;
171 aPointIds.insert(anIdList->GetId(0));
172 aPointIds.insert(anIdList->GetId(1));
173 aPointIds.insert(anIdList->GetId(2));
175 aFace2PointIds.insert(aPointIds);
176 aFace2Visibility.insert(aFaceId);
181 ::TPolygons aPolygons;
183 for (int aFaceId = 0; aFaceId < aNumFaces; aFaceId++) {
184 if(aFace2Visibility.find(aFaceId) == aFace2Visibility.end())
187 vtkCell* aFace = GetFace(aFaceId);
189 vtkIdList *anIdList = aFace->PointIds;
190 vtkIdType aNewPts[3] = {anIdList->GetId(0), anIdList->GetId(1), anIdList->GetId(2)};
192 // To initialize set of points for the plane where the trinangle face belong to
194 aPointIds.insert(aNewPts[0]);
195 aPointIds.insert(aNewPts[1]);
196 aPointIds.insert(aNewPts[2]);
198 // To get know, if the points of the trinagle were already observed
199 bool anIsObserved = aFace2PointIds.find(aPointIds) == aFace2PointIds.end();
200 //cout<<"\taFaceId = "<<aFaceId<<"; anIsObserved = "<<anIsObserved;
201 //cout<<"; aNewPts = {"<<aNewPts[0]<<", "<<aNewPts[1]<<", "<<aNewPts[2]<<"}\n";
204 // To get coordinates of the points of the traingle face
205 vtkFloatingPointType aCoord[3][3];
206 aPoints->GetPoint(aNewPts[0],aCoord[0]);
207 aPoints->GetPoint(aNewPts[1],aCoord[1]);
208 aPoints->GetPoint(aNewPts[2],aCoord[2]);
210 // To calculate plane normal
211 vtkFloatingPointType aVector01[3] = { aCoord[1][0] - aCoord[0][0],
212 aCoord[1][1] - aCoord[0][1],
213 aCoord[1][2] - aCoord[0][2] };
215 vtkFloatingPointType aVector02[3] = { aCoord[2][0] - aCoord[0][0],
216 aCoord[2][1] - aCoord[0][1],
217 aCoord[2][2] - aCoord[0][2] };
219 // To calculate the normal for the triangle
220 vtkFloatingPointType aNormal[3];
221 vtkMath::Cross(aVector02,aVector01,aNormal);
223 vtkMath::Normalize(aNormal);
225 // To calculate what points belong to the plane
226 // To calculate bounds of the point set
227 vtkFloatingPointType aCenter[3] = {0.0, 0.0, 0.0};
229 TPointIds::const_iterator anIter = anInitialPointIds.begin();
230 TPointIds::const_iterator anEndIter = anInitialPointIds.end();
231 for(; anIter != anEndIter; anIter++){
232 vtkFloatingPointType aPntCoord[3];
233 vtkIdType aPntId = *anIter;
234 aPoints->GetPoint(aPntId,aPntCoord);
235 vtkFloatingPointType aDist = vtkPlane::DistanceToPlane(aPntCoord,aNormal,aCoord[0]);
236 //cout<<"\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}; aDist = "<<aDist<<"\n";
237 if(fabs(aDist) < aDistEps){
238 aPointIds.insert(aPntId);
239 aCenter[0] += aPntCoord[0];
240 aCenter[1] += aPntCoord[1];
241 aCenter[2] += aPntCoord[2];
244 int aNbPoints = aPointIds.size();
245 aCenter[0] /= aNbPoints;
246 aCenter[1] /= aNbPoints;
247 aCenter[2] /= aNbPoints;
250 //To sinchronize orientation of the cell and its face
251 vtkFloatingPointType aVectorC[3] = { aCenter[0] - aCellCenter[0],
252 aCenter[1] - aCellCenter[1],
253 aCenter[2] - aCellCenter[2] };
254 vtkMath::Normalize(aVectorC);
256 vtkFloatingPointType aDot = vtkMath::Dot(aNormal,aVectorC);
257 //cout<<"\t\taNormal = {"<<aNormal[0]<<", "<<aNormal[1]<<", "<<aNormal[2]<<"}";
258 //cout<<"; aVectorC = {"<<aVectorC[0]<<", "<<aVectorC[1]<<", "<<aVectorC[2]<<"}\n";
259 //cout<<"\t\taDot = "<<aDot<<"\n";
261 aNormal[0] = -aNormal[0];
262 aNormal[1] = -aNormal[1];
263 aNormal[2] = -aNormal[2];
266 // To calculate the primary direction for point set
267 vtkFloatingPointType aVector0[3] = { aCoord[0][0] - aCenter[0],
268 aCoord[0][1] - aCenter[1],
269 aCoord[0][2] - aCenter[2] };
270 vtkMath::Normalize(aVector0);
272 //cout<<"\t\taCenter = {"<<aCenter[0]<<", "<<aCenter[1]<<", "<<aCenter[2]<<"}";
273 //cout<<"; aVector0 = {"<<aVector0[0]<<", "<<aVector0[1]<<", "<<aVector0[2]<<"}\n";
275 // To calculate the set of points by face those that belong to the plane
276 TFace2PointIds aRemoveFace2PointIds;
278 TFace2PointIds::const_iterator anIter = aFace2PointIds.begin();
279 TFace2PointIds::const_iterator anEndIter = aFace2PointIds.end();
280 for(; anIter != anEndIter; anIter++){
281 const TPointIds& anIds = *anIter;
282 TPointIds anIntersection;
283 std::set_intersection(aPointIds.begin(),aPointIds.end(),
284 anIds.begin(),anIds.end(),
285 std::inserter(anIntersection,anIntersection.begin()));
287 if(anIntersection == anIds){
288 aRemoveFace2PointIds.insert(anIds);
293 // To remove from the set of points by face those that belong to the plane
295 TFace2PointIds::const_iterator anIter = aRemoveFace2PointIds.begin();
296 TFace2PointIds::const_iterator anEndIter = aRemoveFace2PointIds.end();
297 for(; anIter != anEndIter; anIter++){
298 const TPointIds& anIds = *anIter;
299 aFace2PointIds.erase(anIds);
303 // To sort the planar set of the points accrding to the angle
305 typedef std::map<vtkFloatingPointType,vtkIdType> TSortedPointIds;
306 TSortedPointIds aSortedPointIds;
308 TPointIds::const_iterator anIter = aPointIds.begin();
309 TPointIds::const_iterator anEndIter = aPointIds.end();
310 for(; anIter != anEndIter; anIter++){
311 vtkFloatingPointType aPntCoord[3];
312 vtkIdType aPntId = *anIter;
313 aPoints->GetPoint(aPntId,aPntCoord);
314 vtkFloatingPointType aVector[3] = { aPntCoord[0] - aCenter[0],
315 aPntCoord[1] - aCenter[1],
316 aPntCoord[2] - aCenter[2] };
317 vtkMath::Normalize(aVector);
319 vtkFloatingPointType aCross[3];
320 vtkMath::Cross(aVector,aVector0,aCross);
321 bool aGreaterThanPi = vtkMath::Dot(aCross,aNormal) < 0;
322 vtkFloatingPointType aCosinus = vtkMath::Dot(aVector,aVector0);
327 static vtkFloatingPointType a2Pi = 2.0 * vtkMath::Pi();
328 vtkFloatingPointType anAngle = acos(aCosinus);
329 //cout<<"\t\t\taPntId = "<<aPntId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}";
330 //cout<<"; aGreaterThanPi = "<<aGreaterThanPi<<"; aCosinus = "<<aCosinus<<"; anAngle = "<<anAngle<<"\n";
332 anAngle = a2Pi - anAngle;
333 //cout<<"\t\t\t\tanAngle = "<<anAngle<<"\n";
335 aSortedPointIds[anAngle] = aPntId;
338 if(!aSortedPointIds.empty()){
339 int aNumFacePts = aSortedPointIds.size();
340 ::TConnectivities aConnectivities(aNumFacePts);
341 TSortedPointIds::const_iterator anIter = aSortedPointIds.begin();
342 TSortedPointIds::const_iterator anEndIter = aSortedPointIds.end();
343 for(vtkIdType anId = 0; anIter != anEndIter; anIter++, anId++){
344 vtkIdType aPntId = anIter->second;
345 aConnectivities[anId] = GetConnectivity(aPntId);
347 aPolygons.push_back(::TPolygon(aConnectivities,aCenter,aNormal));
353 if(aPolygons.empty())
356 // To check, whether the polygons give a convex polyhedron or not
357 if(theIsCheckConvex){
358 int aNbPolygons = aPolygons.size();
359 for (int aPolygonId = 0; aPolygonId < aNbPolygons; aPolygonId++) {
360 ::TPolygon& aPolygon = aPolygons[aPolygonId];
361 vtkFloatingPointType* aNormal = aPolygon.myNormal;
362 vtkFloatingPointType* anOrigin = aPolygon.myOrigin;
363 //cout<<"\taPolygonId = "<<aPolygonId<<"\n";
364 //cout<<"\t\taNormal = {"<<aNormal[0]<<", "<<aNormal[1]<<", "<<aNormal[2]<<"}";
365 //cout<<"; anOrigin = {"<<anOrigin[0]<<", "<<anOrigin[1]<<", "<<anOrigin[2]<<"}\n";
366 for(vtkIdType aPntId = 0; aPntId < aNumPts; aPntId++){
367 vtkFloatingPointType aPntCoord[3];
368 vtkIdType anId = GetPointId(aPntId);
369 aPoints->GetPoint(anId,aPntCoord);
370 vtkFloatingPointType aDist = vtkPlane::Evaluate(aNormal,anOrigin,aPntCoord);
371 //cout<<"\t\taPntId = "<<anId<<" {"<<aPntCoord[0]<<", "<<aPntCoord[1]<<", "<<aPntCoord[2]<<"}; aDist = "<<aDist<<"\n";
372 if(aDist < -aDistEps)
379 // To pass resulting set of the polygons to the output
381 int aNbPolygons = aPolygons.size();
382 for (int aPolygonId = 0; aPolygonId < aNbPolygons; aPolygonId++) {
383 ::TPolygon& aPolygon = aPolygons[aPolygonId];
384 TConnectivities& aConnectivities = aPolygon.myConnectivities;
385 int aNbPoints = aConnectivities.size();
386 vtkIdType aNewCellId = theOutput->InsertNextCell(VTK_POLYGON,aNbPoints,&aConnectivities[0]);
388 theVTK2ObjIds.push_back(theCellId);
389 theOutputCD->CopyData(thInputCD,theCellId,aNewCellId);
393 //cout<<"\tTriangulator - Ok\n";
400 VTKViewer_OrderedTriangulator
401 ::VTKViewer_OrderedTriangulator():
402 myCell(vtkGenericCell::New())
408 VTKViewer_OrderedTriangulator
409 ::~VTKViewer_OrderedTriangulator()
415 VTKViewer_OrderedTriangulator
418 myInput->GetCell(myCellId,myCell);
419 return myInput->GetPoints();
423 VTKViewer_OrderedTriangulator
426 return myCell->GetNumberOfPoints();
430 VTKViewer_OrderedTriangulator
431 ::GetPointId(vtkIdType thePointId)
433 return myCell->GetPointId(thePointId);
437 VTKViewer_OrderedTriangulator
440 return sqrt(myCell->GetLength2());
444 VTKViewer_OrderedTriangulator
447 return myCell->GetNumberOfFaces();
451 VTKViewer_OrderedTriangulator
452 ::GetFace(vtkIdType theFaceId)
454 return myCell->GetFace(theFaceId);
458 VTKViewer_OrderedTriangulator
459 ::GetCellNeighbors(vtkIdType theCellId,
461 vtkIdList* theCellIds)
463 vtkIdList *anIdList = theFace->PointIds;
464 myInput->GetCellNeighbors(theCellId, anIdList, theCellIds);
468 VTKViewer_OrderedTriangulator
469 ::GetConnectivity(vtkIdType thePntId)
477 VTKViewer_DelaunayTriangulator
478 ::VTKViewer_DelaunayTriangulator():
479 myUnstructuredGrid(vtkUnstructuredGrid::New()),
480 myGeometryFilter(vtkGeometryFilter::New()),
481 myDelaunay3D(vtkDelaunay3D::New()),
482 myFaceIds(vtkIdList::New()),
483 myPoints(vtkPoints::New()),
487 myDelaunay3D->SetInput(myUnstructuredGrid);
488 myGeometryFilter->SetInput(myDelaunay3D->GetOutput());
496 VTKViewer_DelaunayTriangulator
497 ::~VTKViewer_DelaunayTriangulator()
499 myUnstructuredGrid->Delete();
500 myGeometryFilter->Delete();
501 myDelaunay3D->Delete();
508 VTKViewer_DelaunayTriangulator
511 myUnstructuredGrid->Initialize();
512 myUnstructuredGrid->Allocate();
513 myUnstructuredGrid->SetPoints(myPoints);
516 myInput->GetCellPoints(myCellId,aNumPts,myPointIds);
518 vtkFloatingPointType aPntCoord[3];
519 myPoints->SetNumberOfPoints(aNumPts);
520 vtkPoints *anInputPoints = myInput->GetPoints();
521 for (int aPntId = 0; aPntId < aNumPts; aPntId++) {
522 anInputPoints->GetPoint(myPointIds[aPntId],aPntCoord);
523 myPoints->SetPoint(aPntId,aPntCoord);
527 myGeometryFilter->Update();
528 myPolyData = myGeometryFilter->GetOutput();
534 VTKViewer_DelaunayTriangulator
537 return myPoints->GetNumberOfPoints();
541 VTKViewer_DelaunayTriangulator
542 ::GetPointId(vtkIdType thePointId)
548 VTKViewer_DelaunayTriangulator
551 return myPolyData->GetLength();
555 VTKViewer_DelaunayTriangulator
558 return myPolyData->GetNumberOfCells();
562 VTKViewer_DelaunayTriangulator
563 ::GetFace(vtkIdType theFaceId)
565 return myPolyData->GetCell(theFaceId);
569 VTKViewer_DelaunayTriangulator
570 ::GetCellNeighbors(vtkIdType theCellId,
572 vtkIdList* theCellIds)
575 vtkIdList *anIdList = theFace->PointIds;
576 myFaceIds->InsertNextId(myPointIds[anIdList->GetId(0)]);
577 myFaceIds->InsertNextId(myPointIds[anIdList->GetId(1)]);
578 myFaceIds->InsertNextId(myPointIds[anIdList->GetId(2)]);
580 myInput->GetCellNeighbors(theCellId, myFaceIds, theCellIds);
585 VTKViewer_DelaunayTriangulator
586 ::GetConnectivity(vtkIdType thePntId)
588 return myPointIds[thePntId];