src/VTKViewer/VTKViewer_ArcBuilder.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   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.
   7 //
   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.
  12 //
  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
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 //  File   : VTKViewer_ArcBuilder.cxx
  21 //  Author : Roman NIKOLAEV
  22 //  Module : GUI
  23 //
  24 #include "VTKViewer_ArcBuilder.h"
  25
  26 #include <math.h>
  27 #include <float.h>
  28
  29 //VTK includes
  30 #include <vtkMath.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>
  40
  41 #define PRECISION 10e-4
  42 #define ANGLE_PRECISION 0.5
  43 //#define _MY_DEBUG_
  44
  45 #ifdef _MY_DEBUG_
  46 #include <iostream>
  47 #endif
  48
  49
  50 bool CheckAngle(const double compare, const double angle){
  51   if((angle <= compare - ANGLE_PRECISION) || (angle >= compare + ANGLE_PRECISION))
  52     return true;
  53   else
  54     return false;
  55 }
  56
  57 //------------------------------------------------------------------------
  58 /*!
  59  *Class XYZ
  60  *Constructor
  61  */
  62 XYZ::XYZ(double X, double Y, double Z):
  63  x(X),y(Y),z(Z) {}
  64
  65 /*!
  66  * Empty Constructor
  67  */
  68 XYZ::XYZ(){}
  69
  70 /*!
  71  *  Destructor
  72  */
  73 XYZ::~XYZ()
  74 {}
  75
  76 /*
  77  * Calculate modulus
  78  */
  79 double XYZ::Modulus () const {
  80   return sqrt (x * x + y * y + z * z);
  81 }
  82
  83 //------------------------------------------------------------------------
  84 /*!
  85  * Class Pnt
  86  * Constructor
  87  */
  88 Pnt::Pnt(double X,
  89          double Y,
  90          double Z,
  91          double ScalarValue):
  92   coord(X,Y,Z),
  93   scalarValue(ScalarValue)
  94 {
  95 }
  96
  97 Pnt::Pnt()
  98 {
  99 }
 100
 101 /*!
 102  * Destructor
 103  */
 104 Pnt::~Pnt()
 105 {
 106 }
 107
 108 //------------------------------------------------------------------------
 109 /*!
 110  * Class Vec
 111  * Constructor
 112  */
 113 Vec::Vec (const double Xv,
 114           const double Yv,
 115           const double Zv)
 116 {
 117   double D = sqrt (Xv * Xv + Yv * Yv + Zv * Zv);
 118   if(D != 0) {
 119     coord.SetX(Xv / D);
 120     coord.SetY(Yv / D);
 121     coord.SetZ(Zv / D);
 122   }
 123 }
 124
 125 /*!
 126  * Destructor
 127  */
 128 Vec::~Vec(){}
 129
 130
 131 /*!
 132  * Calculate angle between vectors in radians
 133  */
 134 double Vec::AngleBetween(const Vec & Other)
 135 {
 136   double res;
 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 )
 141     d = -1;
 142   else if( d > 1 && d < 1 + PRECISION )
 143     d = 1;
 144   res = acos( d );
 145   return res;
 146 }
 147
 148 /*!
 149  * Calculate angle between vectors in degrees
 150  */
 151 double Vec::AngleBetweenInGrad(const Vec & Other){
 152   return AngleBetween(Other)*vtkMath::DegreesFromRadians(1.);
 153 }
 154
 155 /*
 156  * Calculate vector multiplication
 157 */
 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);
 163   return *aRes;
 164 }
 165
 166 /*---------------------Class Plane --------------------------------*/
 167 /*!
 168  * Constructor
 169  */
 170 Plane::Plane(const Pnt& thePnt1, const Pnt& thePnt2, const Pnt& thePnt3)
 171 {
 172   CalculatePlane(thePnt1,thePnt2,thePnt3);
 173 }
 174
 175 /*!
 176  * Destructor
 177  */
 178 Plane::~Plane()
 179 {}
 180
 181 /*
 182  * Return plane normale
 183  */
 184 Vec Plane::GetNormale() const{
 185   return Vec(myA,myB,myC);
 186 }
 187
 188 /*
 189  * Calculate A,B,C coeefs of plane
 190  */
 191 void Plane::CalculatePlane(const Pnt& thePnt1, const Pnt& thePnt2, const Pnt& thePnt3){
 192
 193   double x1 = thePnt1.GetXYZ().X();
 194   double x2 = thePnt2.GetXYZ().X();
 195   double x3 = thePnt3.GetXYZ().X();
 196
 197   double y1 = thePnt1.GetXYZ().Y();
 198   double y2 = thePnt2.GetXYZ().Y();
 199   double y3 = thePnt3.GetXYZ().Y();
 200
 201   double z1 = thePnt1.GetXYZ().Z();
 202   double z2 = thePnt2.GetXYZ().Z();
 203   double z3 = thePnt3.GetXYZ().Z();
 204
 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);
 208
 209 #ifdef _MY_DEBUG_
 210   cout<<"Plane A: "<<myA<<endl;
 211   cout<<"Plane B: "<<myB<<endl;
 212   cout<<"Plane C: "<<myC<<endl;
 213 #endif
 214
 215 }
 216
 217
 218 /*---------------------Class VTKViewer_ArcBuilder --------------------------------*/
 219 /*!
 220  * Constructor
 221  */
 222 VTKViewer_ArcBuilder::VTKViewer_ArcBuilder(const Pnt& thePnt1,
 223                                            const Pnt& thePnt2,
 224                                            const Pnt& thePnt3,
 225                                            double theAngle):
 226   myStatus(Arc_Error),
 227   myAngle(theAngle)
 228 {
 229   Vec V1(thePnt2.GetXYZ().X()-thePnt1.GetXYZ().X(),
 230          thePnt2.GetXYZ().Y()-thePnt1.GetXYZ().Y(),
 231          thePnt2.GetXYZ().Z()-thePnt1.GetXYZ().Z());
 232
 233   Vec V2(thePnt2.GetXYZ().X()-thePnt3.GetXYZ().X(),
 234          thePnt2.GetXYZ().Y()-thePnt3.GetXYZ().Y(),
 235          thePnt2.GetXYZ().Z()-thePnt3.GetXYZ().Z());
 236
 237   double angle = V1.AngleBetweenInGrad(V2);
 238
 239   //Check that points are not belong one line
 240 #ifdef _MY_DEBUG_
 241   cout<<"Angle for check: "<<angle<<endl;
 242 #endif
 243   if(CheckAngle(180,angle)) {
 244
 245     // Build plane by three points
 246     Plane aPlane(thePnt1,thePnt2,thePnt3);
 247
 248     //Plane normales
 249     Vec aPlaneNormale = aPlane.GetNormale();
 250 #ifdef _MY_DEBUG_
 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;
 254 #endif
 255     //OZ vector
 256     Vec OZ(0,0,1);
 257
 258     //Projection plane normale on XOY
 259     Vec aNormaleProjection (aPlaneNormale.GetXYZ().X(),
 260                             aPlaneNormale.GetXYZ().Y(),
 261                             0);
 262 #ifdef _MY_DEBUG_
 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;
 266 #endif
 267
 268     //Rotation axis
 269     Vec aAxis = aNormaleProjection.VectMultiplication(OZ);
 270 #ifdef _MY_DEBUG_
 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;
 274 #endif
 275     //Rotation angle
 276     double anAngle = OZ.AngleBetweenInGrad(aPlaneNormale);
 277 #ifdef _MY_DEBUG_
 278     std::cout<<"Rotation Angle :"<<anAngle<<endl;
 279 #endif
 280     PntList aInputPnts;
 281     aInputPnts.push_back(thePnt1);
 282     aInputPnts.push_back(thePnt2);
 283     aInputPnts.push_back(thePnt3);
 284
 285     vtkUnstructuredGrid* aGrid = BuildGrid(aInputPnts);
 286
 287     bool needRotation = true;
 288     if(anAngle  == 0 || anAngle == 180)
 289       needRotation = false;
 290
 291     if(aGrid) {
 292       vtkUnstructuredGrid* aTransformedGrid;
 293       if(needRotation) {
 294         aTransformedGrid = TransformGrid(aGrid,aAxis,anAngle);
 295 #ifdef _MY_DEBUG_
 296         cout<<"Need Rotation!!!"<<endl;
 297 #endif
 298       }
 299       else {
 300         aTransformedGrid = aGrid;
 301 #ifdef _MY_DEBUG_
 302         cout<<"Rotation does not need!!!"<<endl;
 303 #endif
 304       }
 305
 306       double coords[3];
 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;
 316       if(needRotation)
 317         anTransArc = TransformGrid(anArc,aAxis,-anAngle);
 318       else
 319         anTransArc = anArc;
 320
 321       myPoints = anTransArc->GetPoints();
 322       myScalarValues = aScalarValues;
 323       myStatus = Arc_Done;
 324     }
 325   }
 326   else{
 327 #ifdef _MY_DEBUG_
 328     std::cout<<"Points lay on the one line !"<<endl;
 329 #endif
 330     PntList aList;
 331     aList.push_back(thePnt1);
 332     aList.push_back(thePnt2);
 333     aList.push_back(thePnt3);
 334     vtkUnstructuredGrid* aGrid = BuildGrid(aList);
 335     myPoints = aGrid->GetPoints();
 336
 337     myScalarValues.clear();
 338     myScalarValues.push_back(thePnt1.GetScalarValue());
 339     myScalarValues.push_back(thePnt2.GetScalarValue());
 340     myScalarValues.push_back(thePnt3.GetScalarValue());
 341     myStatus = Arc_Done;
 342   }
 343 }
 344
 345 /*!
 346  * Destructor
 347  */
 348 VTKViewer_ArcBuilder::~VTKViewer_ArcBuilder()
 349 {}
 350
 351
 352 /*
 353  * Add to the vtkUnstructuredGrid points from input list
 354  */
 355 vtkUnstructuredGrid* VTKViewer_ArcBuilder::BuildGrid(const PntList& theList) const
 356 {
 357   int aListsize = theList.size();
 358   vtkUnstructuredGrid* aGrid = NULL;
 359
 360   if(aListsize != 0) {
 361     aGrid = vtkUnstructuredGrid::New();
 362     vtkPoints* aPoints = vtkPoints::New();
 363     aPoints->SetNumberOfPoints(aListsize);
 364
 365     aGrid->Allocate(aListsize);
 366
 367     PntList::const_iterator it = theList.begin();
 368     int aCounter = 0;
 369     for(;it != theList.end();it++) {
 370       Pnt aPnt =  *it;
 371       aPoints->InsertPoint(aCounter,
 372                            aPnt.GetXYZ().X(),
 373                            aPnt.GetXYZ().Y(),
 374                            aPnt.GetXYZ().Z());
 375       vtkVertex* aVertex = vtkVertex::New();
 376       aVertex->GetPointIds()->SetId(0, aCounter);
 377       aGrid->InsertNextCell(aVertex->GetCellType(), aVertex->GetPointIds());
 378       aCounter++;
 379       aVertex->Delete();
 380     }
 381     aGrid->SetPoints(aPoints);
 382     aPoints->Delete();
 383   }
 384   return aGrid;
 385 }
 386
 387
 388 vtkUnstructuredGrid*
 389 VTKViewer_ArcBuilder::TransformGrid(vtkUnstructuredGrid* theGrid,
 390                                     const Vec& theAxis, const double angle) const
 391 {
 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();
 399 }
 400
 401
 402 void VTKViewer_ArcBuilder::GetAngle(const double theAngle){
 403   myAngle = theAngle;
 404 }
 405
 406 double InterpolateScalarValue(int index, int count, double firstValue, double middleValue, double lastValue)
 407 {
 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;
 414   return value;
 415 }
 416
 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
 421
 422
 423   theScalarValues.clear();
 424
 425   double K1 = 0;
 426   double K2 = 0;
 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;
 432
 433 #ifdef _MY_DEBUG_
 434   std::cout<<"K1 : "<< K1 <<endl;
 435   std::cout<<"K2 : "<< K2 <<endl;
 436 #endif
 437   double xCenter;
 438   if( !okK2 ) //K2 --> infinity
 439     xCenter = (K1*(y1-y3) + (x1+x2))/2.0;
 440
 441   else if( !okK1 ) //K1 --> infinity
 442     xCenter = (K2*(y1-y3) - (x2+x3))/(-2.0);
 443
 444   else
 445     xCenter = (K1*K2*(y1-y3) + K2*(x1+x2) - K1*(x2+x3))/ (2.0*(K2-K1));
 446
 447   double yCenter;
 448
 449   if ( K1 == 0 )
 450     yCenter =  (-1/K2)*(xCenter - (x2+x3)/2.0) + (y2 + y3)/2.0;
 451   else
 452     yCenter =  (-1/K1)*(xCenter - (x1+x2)/2.0) + (y1 + y2)/2.0;
 453
 454 #ifdef _MY_DEBUG_
 455   double zCenter = z;
 456   std::cout<<"xCenter : "<<xCenter<<endl;
 457   std::cout<<"yCenter : "<<yCenter<<endl;
 458   std::cout<<"zCenter : "<<zCenter<<endl;
 459 #endif
 460   double aRadius = sqrt((x1 - xCenter)*(x1 - xCenter) + (y1 - yCenter)*(y1 - yCenter));
 461
 462   double angle1 = GetPointAngleOnCircle(xCenter,yCenter,x1,y1);
 463   double angle2 = GetPointAngleOnCircle(xCenter,yCenter,x2,y2);
 464   double angle3 = GetPointAngleOnCircle(xCenter,yCenter,x3,y3);
 465
 466
 467   double aMaxAngle = vtkMath::RadiansFromDegrees(1.)*myAngle*2;
 468
 469   /*  double aTotalAngle =  fabs(angle3 - angle1);
 470
 471   if (aTotalAngle > vtkMath::Pi())
 472     aTotalAngle = 2*vtkMath::Pi()-aTotalAngle;
 473   */
 474
 475   double aTotalAngle = 0;
 476   IncOrder aOrder = GetArcAngle(angle1,angle2,angle3,&aTotalAngle);
 477
 478   vtkUnstructuredGrid *aC = NULL;
 479
 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;
 486     else
 487       aCurrentAngle+=anIncrementAngle;
 488 #ifdef _MY_DEBUG_
 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;
 493 #endif
 494
 495     PntList aList;
 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);
 503
 504       aList.push_back(aPnt);
 505       theScalarValues.push_back(aPnt.GetScalarValue());
 506       if(aOrder == VTKViewer_ArcBuilder::MINUS)
 507         aCurrentAngle-=anIncrementAngle;
 508       else
 509         aCurrentAngle+=anIncrementAngle;
 510     }
 511     aList.push_back(myPnt3);
 512     theScalarValues.push_back(myPnt3.GetScalarValue());
 513
 514     aC = BuildGrid(aList);
 515 #ifdef _MY_DEBUG_
 516   cout<<"angle1 : "<<angle1*vtkMath::DoubleRadiansToDegrees()<<endl;
 517   cout<<"angle2 : "<<angle2*vtkMath::DoubleRadiansToDegrees()<<endl;
 518   cout<<"angle3 : "<<angle3*vtkMath::DoubleRadiansToDegrees()<<endl;
 519 #endif
 520   }
 521   else{
 522     PntList aList;
 523     aList.push_back(myPnt1);
 524     aList.push_back(myPnt2);
 525     aList.push_back(myPnt3);
 526     aC = BuildGrid(aList);
 527
 528     theScalarValues.push_back(myPnt1.GetScalarValue());
 529     theScalarValues.push_back(myPnt2.GetScalarValue());
 530     theScalarValues.push_back(myPnt3.GetScalarValue());
 531   }
 532   return aC;
 533 }
 534
 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);
 539   if(result < 0 )
 540     result = result+vtkMath::Pi()*2;
 541   return vtkMath::Pi()*2-result;
 542   return result;
 543 }
 544
 545 vtkPoints* VTKViewer_ArcBuilder::GetPoints(){
 546   return myPoints;
 547 }
 548
 549 const std::vector<double>& VTKViewer_ArcBuilder::GetScalarValues()
 550 {
 551   return myScalarValues;
 552 }
 553
 554 VTKViewer_ArcBuilder::IncOrder VTKViewer_ArcBuilder::GetArcAngle( const double& P1, const double& P2, const double& P3,double* Ang){
 555   IncOrder aResult;
 556   if(P1 < P2 && P2 < P3){
 557     *Ang = P3 - P1;
 558     aResult = VTKViewer_ArcBuilder::PLUS;
 559   }
 560   else if((P1 < P3 && P3 < P2) || (P2 < P1 && P1 < P3)){
 561     *Ang = 2*vtkMath::Pi() - P3 + P1;
 562     aResult = VTKViewer_ArcBuilder::MINUS;
 563   }
 564   else if((P2 < P3 && P3 < P1) || (P3 < P1 && P1 < P2)){
 565     *Ang = 2*vtkMath::Pi() - P1 + P3;
 566     aResult = VTKViewer_ArcBuilder::PLUS;
 567   }
 568   else if(P3 < P2 && P2 < P1){
 569     *Ang = P1 - P3;
 570     aResult = VTKViewer_ArcBuilder::MINUS;
 571   }
 572   return aResult;
 573 }
 574
 575 //------------------------------------------------------------------------
 576 Pnt CreatePnt(vtkCell* cell, vtkDataArray* scalars, vtkIdType index)
 577 {
 578   double coord[3];
 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);
 583   return point;
 584 }
 585
 586 //------------------------------------------------------------------------
 587 vtkIdType Build1DArc(vtkIdType cellId, vtkUnstructuredGrid* input,
 588                      vtkPolyData *output,vtkIdType *pts,
 589                      double myMaxArcAngle){
 590
 591   vtkIdType aResult = -1;
 592   vtkIdType *aNewPoints;
 593
 594   vtkDataArray* inputScalars = input->GetPointData()->GetScalars();
 595   vtkDataArray* outputScalars = output->GetPointData()->GetScalars();
 596
 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 );
 602
 603   VTKViewer_ArcBuilder aBuilder(P0,P2,P1,myMaxArcAngle);
 604   if (aBuilder.GetStatus() != VTKViewer_ArcBuilder::Arc_Done) {
 605     return aResult;
 606   }
 607   else{
 608     vtkPoints* aPoints = aBuilder.GetPoints();
 609     std::vector<double> aScalarValues = aBuilder.GetScalarValues();
 610     vtkIdType aNbPts = aPoints->GetNumberOfPoints();
 611     aNewPoints = new vtkIdType[aNbPts];
 612     vtkIdType curID;
 613     vtkIdType aCellType = VTK_POLY_LINE;
 614
 615     aNewPoints[0] = pts[0];
 616     for(vtkIdType idx = 1; idx < aNbPts-1;idx++) {
 617       curID = output->GetPoints()->InsertNextPoint(aPoints->GetPoint(idx));
 618       if( outputScalars )
 619         outputScalars->InsertNextTuple1(aScalarValues[idx]);
 620       aNewPoints[idx] = curID;
 621     }
 622     aNewPoints[aNbPts-1] = pts[1];
 623
 624     aResult = output->InsertNextCell(aCellType,aNbPts,aNewPoints);
 625     return aResult;
 626   }
 627 }
 628
 629 /*!
 630  * Add all points from the input vector theCollection into thePoints.
 631  * Array theIds - it is array with ids of added points.
 632  */
 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,
 637                          vtkIdType* &theIds){
 638   vtkIdType aNbPoints = 0;
 639   vtkIdType anIdCounter = 0;
 640   vtkIdType aNewPntId = 0;
 641
 642   //Compute number of points
 643   std::vector<vtkPoints*>::const_iterator it = theCollection.begin();
 644   for(;it != theCollection.end();it++){
 645     vtkPoints* aPoints = *it;
 646     if(aPoints) {
 647       aNbPoints += aPoints->GetNumberOfPoints()-1; //Because we add all points except last
 648     }
 649   }
 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;
 657
 658     if(aPoints){
 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];
 663         anIdCounter++;
 664       }
 665     }
 666   }
 667   return aNbPoints;
 668 }