1 // Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE
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/ or email : webmaster.salome@opencascade.com
20 // File : VTKViewer_ArcBuilder.cxx
21 // Author : Roman NIKOLAEV
24 #include "VTKViewer_ArcBuilder.h"
31 #include <vtkUnstructuredGrid.h>
32 #include <vtkTransformFilter.h>
33 #include <vtkTransform.h>
34 #include <vtkPoints.h>
35 #include <vtkVertex.h>
36 #include <vtkCellArray.h>
37 #include <vtkTriangle.h>
38 #include <vtkPolyData.h>
39 #include <vtkPointData.h>
41 #define PRECISION 10e-4
42 #define ANGLE_PRECISION 0.5
50 bool CheckAngle(const double compare, const double angle){
51 if((angle <= compare - ANGLE_PRECISION) || (angle >= compare + ANGLE_PRECISION))
57 //------------------------------------------------------------------------
62 XYZ::XYZ(double X, double Y, double Z):
79 double XYZ::Modulus () const {
80 return sqrt (x * x + y * y + z * z);
83 //------------------------------------------------------------------------
93 scalarValue(ScalarValue)
108 //------------------------------------------------------------------------
113 Vec::Vec (const double Xv,
117 double D = sqrt (Xv * Xv + Yv * Yv + Zv * Zv);
132 * Calculate angle between vectors in radians
134 double Vec::AngleBetween(const Vec & Other)
137 double numerator = GetXYZ().X()*Other.GetXYZ().X()+GetXYZ().Y()*Other.GetXYZ().Y()+GetXYZ().Z()*Other.GetXYZ().Z();
138 double denumerator = GetXYZ().Modulus()*Other.GetXYZ().Modulus();
139 double d = numerator/denumerator;
140 if( d < -1 && d > -1 - PRECISION )
142 else if( d > 1 && d < 1 + PRECISION )
149 * Calculate angle between vectors in degrees
151 double Vec::AngleBetweenInGrad(const Vec & Other){
152 return AngleBetween(Other)*vtkMath::DegreesFromRadians(1.);
156 * Calculate vector multiplication
158 Vec Vec::VectMultiplication(const Vec & Other) const{
159 double x = GetXYZ().Y()*Other.GetXYZ().Z() - GetXYZ().Z()*Other.GetXYZ().Y();
160 double y = GetXYZ().Z()*Other.GetXYZ().X() - GetXYZ().X()*Other.GetXYZ().Z();
161 double z = GetXYZ().X()*Other.GetXYZ().Y() - GetXYZ().Y()*Other.GetXYZ().X();
162 Vec *aRes = new Vec(x,y,z);
166 /*---------------------Class Plane --------------------------------*/
170 Plane::Plane(const Pnt& thePnt1, const Pnt& thePnt2, const Pnt& thePnt3)
172 CalculatePlane(thePnt1,thePnt2,thePnt3);
182 * Return plane normale
184 Vec Plane::GetNormale() const{
185 return Vec(myA,myB,myC);
189 * Calculate A,B,C coeefs of plane
191 void Plane::CalculatePlane(const Pnt& thePnt1, const Pnt& thePnt2, const Pnt& thePnt3){
193 double x1 = thePnt1.GetXYZ().X();
194 double x2 = thePnt2.GetXYZ().X();
195 double x3 = thePnt3.GetXYZ().X();
197 double y1 = thePnt1.GetXYZ().Y();
198 double y2 = thePnt2.GetXYZ().Y();
199 double y3 = thePnt3.GetXYZ().Y();
201 double z1 = thePnt1.GetXYZ().Z();
202 double z2 = thePnt2.GetXYZ().Z();
203 double z3 = thePnt3.GetXYZ().Z();
205 myA = y1*(z2 - z3) + y2*(z3 - z1) + y3*(z1 - z2);
206 myB = z1*(x2 - x3) + z2*(x3 - x1) + z3*(x1 - x2);
207 myC = x1*(y2 - y3) + x2*(y3 - y1) + x3*(y1 - y2);
210 cout<<"Plane A: "<<myA<<endl;
211 cout<<"Plane B: "<<myB<<endl;
212 cout<<"Plane C: "<<myC<<endl;
218 /*---------------------Class VTKViewer_ArcBuilder --------------------------------*/
222 VTKViewer_ArcBuilder::VTKViewer_ArcBuilder(const Pnt& thePnt1,
229 Vec V1(thePnt2.GetXYZ().X()-thePnt1.GetXYZ().X(),
230 thePnt2.GetXYZ().Y()-thePnt1.GetXYZ().Y(),
231 thePnt2.GetXYZ().Z()-thePnt1.GetXYZ().Z());
233 Vec V2(thePnt2.GetXYZ().X()-thePnt3.GetXYZ().X(),
234 thePnt2.GetXYZ().Y()-thePnt3.GetXYZ().Y(),
235 thePnt2.GetXYZ().Z()-thePnt3.GetXYZ().Z());
237 double angle = V1.AngleBetweenInGrad(V2);
239 //Check that points are not belong one line
241 cout<<"Angle for check: "<<angle<<endl;
243 if(CheckAngle(180,angle)) {
245 // Build plane by three points
246 Plane aPlane(thePnt1,thePnt2,thePnt3);
249 Vec aPlaneNormale = aPlane.GetNormale();
251 std::cout<<"X Normale: "<<aPlaneNormale.GetXYZ().X()<<std::endl;
252 std::cout<<"Y Normale: "<<aPlaneNormale.GetXYZ().Y()<<std::endl;
253 std::cout<<"Z Normale: "<<aPlaneNormale.GetXYZ().Z()<<std::endl;
258 //Projection plane normale on XOY
259 Vec aNormaleProjection (aPlaneNormale.GetXYZ().X(),
260 aPlaneNormale.GetXYZ().Y(),
263 std::cout<<"X Normale Projection: "<<aNormaleProjection.GetXYZ().X()<<std::endl;
264 std::cout<<"Y Normale Projection: "<<aNormaleProjection.GetXYZ().Y()<<std::endl;
265 std::cout<<"Z Normale Projection: "<<aNormaleProjection.GetXYZ().Z()<<std::endl;
269 Vec aAxis = aNormaleProjection.VectMultiplication(OZ);
271 std::cout<<"X Axis: "<<aAxis.GetXYZ().X()<<std::endl;
272 std::cout<<"Y Axis: "<<aAxis.GetXYZ().Y()<<std::endl;
273 std::cout<<"Z Axis: "<<aAxis.GetXYZ().Z()<<std::endl;
276 double anAngle = OZ.AngleBetweenInGrad(aPlaneNormale);
278 std::cout<<"Rotation Angle :"<<anAngle<<endl;
281 aInputPnts.push_back(thePnt1);
282 aInputPnts.push_back(thePnt2);
283 aInputPnts.push_back(thePnt3);
285 vtkUnstructuredGrid* aGrid = BuildGrid(aInputPnts);
287 bool needRotation = true;
288 if(anAngle == 0 || anAngle == 180)
289 needRotation = false;
292 vtkUnstructuredGrid* aTransformedGrid;
294 aTransformedGrid = TransformGrid(aGrid,aAxis,anAngle);
296 cout<<"Need Rotation!!!"<<endl;
300 aTransformedGrid = aGrid;
302 cout<<"Rotation does not need!!!"<<endl;
307 aTransformedGrid->GetPoint(0,coords);
308 myPnt1 = Pnt(coords[0],coords[1],coords[2], thePnt1.GetScalarValue());
309 aTransformedGrid->GetPoint(1,coords);
310 myPnt2 = Pnt(coords[0],coords[1],coords[2], thePnt2.GetScalarValue());
311 aTransformedGrid->GetPoint(2,coords);
312 myPnt3 = Pnt(coords[0],coords[1],coords[2], thePnt3.GetScalarValue());
313 std::vector<double> aScalarValues;
314 vtkUnstructuredGrid* anArc = BuildArc(aScalarValues);
315 vtkUnstructuredGrid* anTransArc;
317 anTransArc = TransformGrid(anArc,aAxis,-anAngle);
321 myPoints = anTransArc->GetPoints();
322 myScalarValues = aScalarValues;
328 std::cout<<"Points lay on the one line !"<<endl;
331 aList.push_back(thePnt1);
332 aList.push_back(thePnt2);
333 aList.push_back(thePnt3);
334 vtkUnstructuredGrid* aGrid = BuildGrid(aList);
335 myPoints = aGrid->GetPoints();
337 myScalarValues.clear();
338 myScalarValues.push_back(thePnt1.GetScalarValue());
339 myScalarValues.push_back(thePnt2.GetScalarValue());
340 myScalarValues.push_back(thePnt3.GetScalarValue());
348 VTKViewer_ArcBuilder::~VTKViewer_ArcBuilder()
353 * Add to the vtkUnstructuredGrid points from input list
355 vtkUnstructuredGrid* VTKViewer_ArcBuilder::BuildGrid(const PntList& theList) const
357 int aListsize = theList.size();
358 vtkUnstructuredGrid* aGrid = NULL;
361 aGrid = vtkUnstructuredGrid::New();
362 vtkPoints* aPoints = vtkPoints::New();
363 aPoints->SetNumberOfPoints(aListsize);
365 aGrid->Allocate(aListsize);
367 PntList::const_iterator it = theList.begin();
369 for(;it != theList.end();it++) {
371 aPoints->InsertPoint(aCounter,
375 vtkVertex* aVertex = vtkVertex::New();
376 aVertex->GetPointIds()->SetId(0, aCounter);
377 aGrid->InsertNextCell(aVertex->GetCellType(), aVertex->GetPointIds());
381 aGrid->SetPoints(aPoints);
389 VTKViewer_ArcBuilder::TransformGrid(vtkUnstructuredGrid* theGrid,
390 const Vec& theAxis, const double angle) const
392 vtkTransform *aTransform = vtkTransform::New();
393 aTransform->RotateWXYZ(angle, theAxis.GetXYZ().X(), theAxis.GetXYZ().Y(), theAxis.GetXYZ().Z());
394 vtkTransformFilter* aTransformFilter = vtkTransformFilter::New();
395 aTransformFilter->SetTransform(aTransform);
396 aTransformFilter->SetInputData(theGrid);
397 aTransform->Delete();
398 return aTransformFilter->GetUnstructuredGridOutput();
402 void VTKViewer_ArcBuilder::GetAngle(const double theAngle){
406 double InterpolateScalarValue(int index, int count, double firstValue, double middleValue, double lastValue)
408 bool isFirstHalf = index <= count / 2;
409 double first = isFirstHalf ? firstValue : lastValue;
410 double last = isFirstHalf ? middleValue : middleValue;
411 double ratio = (double)index / (double)count;
412 double position = isFirstHalf ? ratio * 2 : ( 1 - ratio ) * 2;
413 double value = first + (last - first) * position;
417 vtkUnstructuredGrid* VTKViewer_ArcBuilder::BuildArc(std::vector<double>& theScalarValues){
418 double x1 = myPnt1.GetXYZ().X(); double x2 = myPnt2.GetXYZ().X(); double x3 = myPnt3.GetXYZ().X();
419 double y1 = myPnt1.GetXYZ().Y(); double y2 = myPnt2.GetXYZ().Y(); double y3 = myPnt3.GetXYZ().Y();
420 double z = myPnt1.GetXYZ().Z(); //Points on plane || XOY
423 theScalarValues.clear();
427 bool okK1 = false, okK2 = false;
428 if ( fabs(x2 - x1) > DBL_MIN )
429 K1 = (y2 - y1)/(x2 - x1), okK1 = true;
430 if ( fabs(x3 - x2) > DBL_MIN )
431 K2 = (y3 - y2)/(x3 - x2), okK2 = true;
434 std::cout<<"K1 : "<< K1 <<endl;
435 std::cout<<"K2 : "<< K2 <<endl;
438 if( !okK2 ) //K2 --> infinity
439 xCenter = (K1*(y1-y3) + (x1+x2))/2.0;
441 else if( !okK1 ) //K1 --> infinity
442 xCenter = (K2*(y1-y3) - (x2+x3))/(-2.0);
445 xCenter = (K1*K2*(y1-y3) + K2*(x1+x2) - K1*(x2+x3))/ (2.0*(K2-K1));
450 yCenter = (-1/K2)*(xCenter - (x2+x3)/2.0) + (y2 + y3)/2.0;
452 yCenter = (-1/K1)*(xCenter - (x1+x2)/2.0) + (y1 + y2)/2.0;
456 std::cout<<"xCenter : "<<xCenter<<endl;
457 std::cout<<"yCenter : "<<yCenter<<endl;
458 std::cout<<"zCenter : "<<zCenter<<endl;
460 double aRadius = sqrt((x1 - xCenter)*(x1 - xCenter) + (y1 - yCenter)*(y1 - yCenter));
462 double angle1 = GetPointAngleOnCircle(xCenter,yCenter,x1,y1);
463 double angle2 = GetPointAngleOnCircle(xCenter,yCenter,x2,y2);
464 double angle3 = GetPointAngleOnCircle(xCenter,yCenter,x3,y3);
467 double aMaxAngle = vtkMath::RadiansFromDegrees(1.)*myAngle*2;
469 /* double aTotalAngle = fabs(angle3 - angle1);
471 if (aTotalAngle > vtkMath::Pi())
472 aTotalAngle = 2*vtkMath::Pi()-aTotalAngle;
475 double aTotalAngle = 0;
476 IncOrder aOrder = GetArcAngle(angle1,angle2,angle3,&aTotalAngle);
478 vtkUnstructuredGrid *aC = NULL;
480 if(aTotalAngle > aMaxAngle) {
481 int nbSteps = int(aTotalAngle/aMaxAngle)+1;
482 double anIncrementAngle = aTotalAngle/nbSteps;
483 double aCurrentAngle = angle1;
484 if(aOrder == VTKViewer_ArcBuilder::MINUS)
485 aCurrentAngle-=anIncrementAngle;
487 aCurrentAngle+=anIncrementAngle;
489 cout<<"Total angle :"<<aTotalAngle<<endl;
490 cout<<"Max Increment Angle :"<<aMaxAngle<<endl;
491 cout<<"Real Increment angle :"<<anIncrementAngle<<endl;
492 cout<<"Nb Steps : "<<nbSteps<<endl;
496 aList.push_back(myPnt1);
497 theScalarValues.push_back(myPnt1.GetScalarValue());
498 for(int i=1;i<=nbSteps-1;i++){
499 double x = xCenter + aRadius*cos(aCurrentAngle);
500 double y = yCenter + aRadius*sin(aCurrentAngle);
501 double value = InterpolateScalarValue(i, nbSteps, myPnt1.GetScalarValue(), myPnt2.GetScalarValue(), myPnt3.GetScalarValue());
502 Pnt aPnt(x,y,z,value);
504 aList.push_back(aPnt);
505 theScalarValues.push_back(aPnt.GetScalarValue());
506 if(aOrder == VTKViewer_ArcBuilder::MINUS)
507 aCurrentAngle-=anIncrementAngle;
509 aCurrentAngle+=anIncrementAngle;
511 aList.push_back(myPnt3);
512 theScalarValues.push_back(myPnt3.GetScalarValue());
514 aC = BuildGrid(aList);
516 cout<<"angle1 : "<<angle1*vtkMath::DoubleRadiansToDegrees()<<endl;
517 cout<<"angle2 : "<<angle2*vtkMath::DoubleRadiansToDegrees()<<endl;
518 cout<<"angle3 : "<<angle3*vtkMath::DoubleRadiansToDegrees()<<endl;
523 aList.push_back(myPnt1);
524 aList.push_back(myPnt2);
525 aList.push_back(myPnt3);
526 aC = BuildGrid(aList);
528 theScalarValues.push_back(myPnt1.GetScalarValue());
529 theScalarValues.push_back(myPnt2.GetScalarValue());
530 theScalarValues.push_back(myPnt3.GetScalarValue());
535 double VTKViewer_ArcBuilder::
536 GetPointAngleOnCircle(const double theXCenter, const double theYCenter,
537 const double theXPoint, const double theYPoint){
538 double result = atan2(theYCenter - theYPoint, theXPoint - theXCenter);
540 result = result+vtkMath::Pi()*2;
541 return vtkMath::Pi()*2-result;
545 vtkPoints* VTKViewer_ArcBuilder::GetPoints(){
549 const std::vector<double>& VTKViewer_ArcBuilder::GetScalarValues()
551 return myScalarValues;
554 VTKViewer_ArcBuilder::IncOrder VTKViewer_ArcBuilder::GetArcAngle( const double& P1, const double& P2, const double& P3,double* Ang){
556 if(P1 < P2 && P2 < P3){
558 aResult = VTKViewer_ArcBuilder::PLUS;
560 else if((P1 < P3 && P3 < P2) || (P2 < P1 && P1 < P3)){
561 *Ang = 2*vtkMath::Pi() - P3 + P1;
562 aResult = VTKViewer_ArcBuilder::MINUS;
564 else if((P2 < P3 && P3 < P1) || (P3 < P1 && P1 < P2)){
565 *Ang = 2*vtkMath::Pi() - P1 + P3;
566 aResult = VTKViewer_ArcBuilder::PLUS;
568 else if(P3 < P2 && P2 < P1){
570 aResult = VTKViewer_ArcBuilder::MINUS;
575 //------------------------------------------------------------------------
576 Pnt CreatePnt(vtkCell* cell, vtkDataArray* scalars, vtkIdType index)
579 cell->GetPoints()->GetPoint(index, coord);
580 vtkIdType pointId = cell->GetPointId(index);
581 double scalarValue = scalars ? scalars->GetTuple1(pointId) : 0;
582 Pnt point(coord[0], coord[1], coord[2], scalarValue);
586 //------------------------------------------------------------------------
587 vtkIdType Build1DArc(vtkIdType cellId, vtkUnstructuredGrid* input,
588 vtkPolyData *output,vtkIdType *pts,
589 double myMaxArcAngle){
591 vtkIdType aResult = -1;
592 vtkIdType *aNewPoints;
594 vtkDataArray* inputScalars = input->GetPointData()->GetScalars();
595 vtkDataArray* outputScalars = output->GetPointData()->GetScalars();
597 vtkCell* aCell = input->GetCell(cellId);
598 //Get All points from input cell
599 Pnt P0 = CreatePnt( aCell, inputScalars, 0 );
600 Pnt P1 = CreatePnt( aCell, inputScalars, 1 );
601 Pnt P2 = CreatePnt( aCell, inputScalars, 2 );
603 VTKViewer_ArcBuilder aBuilder(P0,P2,P1,myMaxArcAngle);
604 if (aBuilder.GetStatus() != VTKViewer_ArcBuilder::Arc_Done) {
608 vtkPoints* aPoints = aBuilder.GetPoints();
609 std::vector<double> aScalarValues = aBuilder.GetScalarValues();
610 vtkIdType aNbPts = aPoints->GetNumberOfPoints();
611 aNewPoints = new vtkIdType[aNbPts];
613 vtkIdType aCellType = VTK_POLY_LINE;
615 aNewPoints[0] = pts[0];
616 for(vtkIdType idx = 1; idx < aNbPts-1;idx++) {
617 curID = output->GetPoints()->InsertNextPoint(aPoints->GetPoint(idx));
619 outputScalars->InsertNextTuple1(aScalarValues[idx]);
620 aNewPoints[idx] = curID;
622 aNewPoints[aNbPts-1] = pts[1];
624 aResult = output->InsertNextCell(aCellType,aNbPts,aNewPoints);
630 * Add all points from the input vector theCollection into thePoints.
631 * Array theIds - it is array with ids of added points.
633 vtkIdType MergevtkPoints(const std::vector<vtkPoints*>& theCollection,
634 const std::vector< std::vector<double> >& theScalarCollection,
635 vtkPoints* thePoints,
636 std::map<int, double>& thePntId2ScalarValue,
638 vtkIdType aNbPoints = 0;
639 vtkIdType anIdCounter = 0;
640 vtkIdType aNewPntId = 0;
642 //Compute number of points
643 std::vector<vtkPoints*>::const_iterator it = theCollection.begin();
644 for(;it != theCollection.end();it++){
645 vtkPoints* aPoints = *it;
647 aNbPoints += aPoints->GetNumberOfPoints()-1; //Because we add all points except last
650 it = theCollection.begin();
651 std::vector< std::vector<double> >::const_iterator itScalar = theScalarCollection.begin();
652 theIds = new vtkIdType[aNbPoints];
653 // ..and add all points
654 for(;it != theCollection.end() && itScalar != theScalarCollection.end(); it++, itScalar++){
655 vtkPoints* aPoints = *it;
656 std::vector<double> aScalarValues = *itScalar;
659 for(vtkIdType idx = 0;idx < aPoints->GetNumberOfPoints()-1;idx++){
660 aNewPntId = thePoints->InsertNextPoint(aPoints->GetPoint(idx));
661 theIds[anIdCounter] = aNewPntId;
662 thePntId2ScalarValue[ aNewPntId ] = aScalarValues[idx];