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);
295 aTransformedGrid->Update();
297 cout<<"Need Rotation!!!"<<endl;
301 aTransformedGrid = aGrid;
303 cout<<"Rotation does not need!!!"<<endl;
308 aTransformedGrid->GetPoint(0,coords);
309 myPnt1 = Pnt(coords[0],coords[1],coords[2], thePnt1.GetScalarValue());
310 aTransformedGrid->GetPoint(1,coords);
311 myPnt2 = Pnt(coords[0],coords[1],coords[2], thePnt2.GetScalarValue());
312 aTransformedGrid->GetPoint(2,coords);
313 myPnt3 = Pnt(coords[0],coords[1],coords[2], thePnt3.GetScalarValue());
314 std::vector<double> aScalarValues;
315 vtkUnstructuredGrid* anArc = BuildArc(aScalarValues);
316 vtkUnstructuredGrid* anTransArc;
318 anTransArc = TransformGrid(anArc,aAxis,-anAngle);
319 anTransArc->Update();
324 myPoints = anTransArc->GetPoints();
325 myScalarValues = aScalarValues;
331 std::cout<<"Points lay on the one line !"<<endl;
334 aList.push_back(thePnt1);
335 aList.push_back(thePnt2);
336 aList.push_back(thePnt3);
337 vtkUnstructuredGrid* aGrid = BuildGrid(aList);
339 myPoints = aGrid->GetPoints();
341 myScalarValues.clear();
342 myScalarValues.push_back(thePnt1.GetScalarValue());
343 myScalarValues.push_back(thePnt2.GetScalarValue());
344 myScalarValues.push_back(thePnt3.GetScalarValue());
352 VTKViewer_ArcBuilder::~VTKViewer_ArcBuilder()
357 * Add to the vtkUnstructuredGrid points from input list
359 vtkUnstructuredGrid* VTKViewer_ArcBuilder::BuildGrid(const PntList& theList) const
361 int aListsize = theList.size();
362 vtkUnstructuredGrid* aGrid = NULL;
365 aGrid = vtkUnstructuredGrid::New();
366 vtkPoints* aPoints = vtkPoints::New();
367 aPoints->SetNumberOfPoints(aListsize);
369 aGrid->Allocate(aListsize);
371 PntList::const_iterator it = theList.begin();
373 for(;it != theList.end();it++) {
375 aPoints->InsertPoint(aCounter,
379 vtkVertex* aVertex = vtkVertex::New();
380 aVertex->GetPointIds()->SetId(0, aCounter);
381 aGrid->InsertNextCell(aVertex->GetCellType(), aVertex->GetPointIds());
385 aGrid->SetPoints(aPoints);
393 VTKViewer_ArcBuilder::TransformGrid(vtkUnstructuredGrid* theGrid,
394 const Vec& theAxis, const double angle) const
396 vtkTransform *aTransform = vtkTransform::New();
397 aTransform->RotateWXYZ(angle, theAxis.GetXYZ().X(), theAxis.GetXYZ().Y(), theAxis.GetXYZ().Z());
398 vtkTransformFilter* aTransformFilter = vtkTransformFilter::New();
399 aTransformFilter->SetTransform(aTransform);
400 aTransformFilter->SetInput(theGrid);
401 aTransform->Delete();
402 return aTransformFilter->GetUnstructuredGridOutput();
406 void VTKViewer_ArcBuilder::GetAngle(const double theAngle){
410 double InterpolateScalarValue(int index, int count, double firstValue, double middleValue, double lastValue)
412 bool isFirstHalf = index <= count / 2;
413 double first = isFirstHalf ? firstValue : lastValue;
414 double last = isFirstHalf ? middleValue : middleValue;
415 double ratio = (double)index / (double)count;
416 double position = isFirstHalf ? ratio * 2 : ( 1 - ratio ) * 2;
417 double value = first + (last - first) * position;
421 vtkUnstructuredGrid* VTKViewer_ArcBuilder::BuildArc(std::vector<double>& theScalarValues){
422 double x1 = myPnt1.GetXYZ().X(); double x2 = myPnt2.GetXYZ().X(); double x3 = myPnt3.GetXYZ().X();
423 double y1 = myPnt1.GetXYZ().Y(); double y2 = myPnt2.GetXYZ().Y(); double y3 = myPnt3.GetXYZ().Y();
424 double z = myPnt1.GetXYZ().Z(); //Points on plane || XOY
427 theScalarValues.clear();
431 bool okK1 = false, okK2 = false;
432 if ( fabs(x2 - x1) > DBL_MIN )
433 K1 = (y2 - y1)/(x2 - x1), okK1 = true;
434 if ( fabs(x3 - x2) > DBL_MIN )
435 K2 = (y3 - y2)/(x3 - x2), okK2 = true;
438 std::cout<<"K1 : "<< K1 <<endl;
439 std::cout<<"K2 : "<< K2 <<endl;
442 if( !okK2 ) //K2 --> infinity
443 xCenter = (K1*(y1-y3) + (x1+x2))/2.0;
445 else if( !okK1 ) //K1 --> infinity
446 xCenter = (K2*(y1-y3) - (x2+x3))/(-2.0);
449 xCenter = (K1*K2*(y1-y3) + K2*(x1+x2) - K1*(x2+x3))/ (2.0*(K2-K1));
454 yCenter = (-1/K2)*(xCenter - (x2+x3)/2.0) + (y2 + y3)/2.0;
456 yCenter = (-1/K1)*(xCenter - (x1+x2)/2.0) + (y1 + y2)/2.0;
460 std::cout<<"xCenter : "<<xCenter<<endl;
461 std::cout<<"yCenter : "<<yCenter<<endl;
462 std::cout<<"zCenter : "<<zCenter<<endl;
464 double aRadius = sqrt((x1 - xCenter)*(x1 - xCenter) + (y1 - yCenter)*(y1 - yCenter));
466 double angle1 = GetPointAngleOnCircle(xCenter,yCenter,x1,y1);
467 double angle2 = GetPointAngleOnCircle(xCenter,yCenter,x2,y2);
468 double angle3 = GetPointAngleOnCircle(xCenter,yCenter,x3,y3);
471 double aMaxAngle = vtkMath::RadiansFromDegrees(1.)*myAngle*2;
473 /* double aTotalAngle = fabs(angle3 - angle1);
475 if (aTotalAngle > vtkMath::DoublePi())
476 aTotalAngle = 2*vtkMath::DoublePi()-aTotalAngle;
479 double aTotalAngle = 0;
480 IncOrder aOrder = GetArcAngle(angle1,angle2,angle3,&aTotalAngle);
482 vtkUnstructuredGrid *aC = NULL;
484 if(aTotalAngle > aMaxAngle) {
485 int nbSteps = int(aTotalAngle/aMaxAngle)+1;
486 double anIncrementAngle = aTotalAngle/nbSteps;
487 double aCurrentAngle = angle1;
488 if(aOrder == VTKViewer_ArcBuilder::MINUS)
489 aCurrentAngle-=anIncrementAngle;
491 aCurrentAngle+=anIncrementAngle;
493 cout<<"Total angle :"<<aTotalAngle<<endl;
494 cout<<"Max Increment Angle :"<<aMaxAngle<<endl;
495 cout<<"Real Increment angle :"<<anIncrementAngle<<endl;
496 cout<<"Nb Steps : "<<nbSteps<<endl;
500 aList.push_back(myPnt1);
501 theScalarValues.push_back(myPnt1.GetScalarValue());
502 for(int i=1;i<=nbSteps-1;i++){
503 double x = xCenter + aRadius*cos(aCurrentAngle);
504 double y = yCenter + aRadius*sin(aCurrentAngle);
505 double value = InterpolateScalarValue(i, nbSteps, myPnt1.GetScalarValue(), myPnt2.GetScalarValue(), myPnt3.GetScalarValue());
506 Pnt aPnt(x,y,z,value);
508 aList.push_back(aPnt);
509 theScalarValues.push_back(aPnt.GetScalarValue());
510 if(aOrder == VTKViewer_ArcBuilder::MINUS)
511 aCurrentAngle-=anIncrementAngle;
513 aCurrentAngle+=anIncrementAngle;
515 aList.push_back(myPnt3);
516 theScalarValues.push_back(myPnt3.GetScalarValue());
518 aC = BuildGrid(aList);
520 cout<<"angle1 : "<<angle1*vtkMath::DoubleRadiansToDegrees()<<endl;
521 cout<<"angle2 : "<<angle2*vtkMath::DoubleRadiansToDegrees()<<endl;
522 cout<<"angle3 : "<<angle3*vtkMath::DoubleRadiansToDegrees()<<endl;
527 aList.push_back(myPnt1);
528 aList.push_back(myPnt2);
529 aList.push_back(myPnt3);
530 aC = BuildGrid(aList);
532 theScalarValues.push_back(myPnt1.GetScalarValue());
533 theScalarValues.push_back(myPnt2.GetScalarValue());
534 theScalarValues.push_back(myPnt3.GetScalarValue());
539 double VTKViewer_ArcBuilder::
540 GetPointAngleOnCircle(const double theXCenter, const double theYCenter,
541 const double theXPoint, const double theYPoint){
542 double result = atan2(theYCenter - theYPoint, theXPoint - theXCenter);
544 result = result+vtkMath::DoublePi()*2;
545 return vtkMath::DoublePi()*2-result;
549 vtkPoints* VTKViewer_ArcBuilder::GetPoints(){
553 const std::vector<double>& VTKViewer_ArcBuilder::GetScalarValues()
555 return myScalarValues;
558 VTKViewer_ArcBuilder::IncOrder VTKViewer_ArcBuilder::GetArcAngle( const double& P1, const double& P2, const double& P3,double* Ang){
560 if(P1 < P2 && P2 < P3){
562 aResult = VTKViewer_ArcBuilder::PLUS;
564 else if((P1 < P3 && P3 < P2) || (P2 < P1 && P1 < P3)){
565 *Ang = 2*vtkMath::DoublePi() - P3 + P1;
566 aResult = VTKViewer_ArcBuilder::MINUS;
568 else if((P2 < P3 && P3 < P1) || (P3 < P1 && P1 < P2)){
569 *Ang = 2*vtkMath::DoublePi() - P1 + P3;
570 aResult = VTKViewer_ArcBuilder::PLUS;
572 else if(P3 < P2 && P2 < P1){
574 aResult = VTKViewer_ArcBuilder::MINUS;
579 //------------------------------------------------------------------------
580 Pnt CreatePnt(vtkCell* cell, vtkDataArray* scalars, vtkIdType index)
582 vtkFloatingPointType coord[3];
583 cell->GetPoints()->GetPoint(index, coord);
584 vtkIdType pointId = cell->GetPointId(index);
585 double scalarValue = scalars ? scalars->GetTuple1(pointId) : 0;
586 Pnt point(coord[0], coord[1], coord[2], scalarValue);
590 //------------------------------------------------------------------------
591 vtkIdType Build1DArc(vtkIdType cellId, vtkUnstructuredGrid* input,
592 vtkPolyData *output,vtkIdType *pts,
593 vtkFloatingPointType myMaxArcAngle){
595 vtkIdType aResult = -1;
596 vtkIdType *aNewPoints;
598 vtkDataArray* inputScalars = input->GetPointData()->GetScalars();
599 vtkDataArray* outputScalars = output->GetPointData()->GetScalars();
601 vtkCell* aCell = input->GetCell(cellId);
602 //Get All points from input cell
603 Pnt P0 = CreatePnt( aCell, inputScalars, 0 );
604 Pnt P1 = CreatePnt( aCell, inputScalars, 1 );
605 Pnt P2 = CreatePnt( aCell, inputScalars, 2 );
607 VTKViewer_ArcBuilder aBuilder(P0,P2,P1,myMaxArcAngle);
608 if (aBuilder.GetStatus() != VTKViewer_ArcBuilder::Arc_Done) {
612 vtkPoints* aPoints = aBuilder.GetPoints();
613 std::vector<double> aScalarValues = aBuilder.GetScalarValues();
614 vtkIdType aNbPts = aPoints->GetNumberOfPoints();
615 aNewPoints = new vtkIdType[aNbPts];
617 vtkIdType aCellType = VTK_POLY_LINE;
619 aNewPoints[0] = pts[0];
620 for(vtkIdType idx = 1; idx < aNbPts-1;idx++) {
621 curID = output->GetPoints()->InsertNextPoint(aPoints->GetPoint(idx));
623 outputScalars->InsertNextTuple1(aScalarValues[idx]);
624 aNewPoints[idx] = curID;
626 aNewPoints[aNbPts-1] = pts[1];
628 aResult = output->InsertNextCell(aCellType,aNbPts,aNewPoints);
634 * Add all points from the input vector theCollection into thePoints.
635 * Array theIds - it is array with ids of added points.
637 vtkIdType MergevtkPoints(const std::vector<vtkPoints*>& theCollection,
638 const std::vector< std::vector<double> >& theScalarCollection,
639 vtkPoints* thePoints,
640 std::map<int, double>& thePntId2ScalarValue,
642 vtkIdType aNbPoints = 0;
643 vtkIdType anIdCounter = 0;
644 vtkIdType aNewPntId = 0;
646 //Compute number of points
647 std::vector<vtkPoints*>::const_iterator it = theCollection.begin();
648 for(;it != theCollection.end();it++){
649 vtkPoints* aPoints = *it;
651 aNbPoints += aPoints->GetNumberOfPoints()-1; //Because we add all points except last
654 it = theCollection.begin();
655 std::vector< std::vector<double> >::const_iterator itScalar = theScalarCollection.begin();
656 theIds = new vtkIdType[aNbPoints];
657 // ..and add all points
658 for(;it != theCollection.end() && itScalar != theScalarCollection.end(); it++, itScalar++){
659 vtkPoints* aPoints = *it;
660 std::vector<double> aScalarValues = *itScalar;
663 for(vtkIdType idx = 0;idx < aPoints->GetNumberOfPoints()-1;idx++){
664 aNewPntId = thePoints->InsertNextPoint(aPoints->GetPoint(idx));
665 theIds[anIdCounter] = aNewPntId;
666 thePntId2ScalarValue[ aNewPntId ] = aScalarValues[idx];