src/VTKViewer/VTKViewer_ArcBuilder.cxx

   1 // Copyright (C) 2007-2012  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         aTransformedGrid->Update();
 296 #ifdef _MY_DEBUG_
 297         cout<<"Need Rotation!!!"<<endl;
 298 #endif
 299       }
 300       else {
 301         aTransformedGrid = aGrid;
 302 #ifdef _MY_DEBUG_
 303         cout<<"Rotation does not need!!!"<<endl;
 304 #endif
 305       }
 306
 307       double coords[3];
 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;
 317       if(needRotation) {
 318         anTransArc = TransformGrid(anArc,aAxis,-anAngle);
 319         anTransArc->Update();
 320       }
 321       else
 322         anTransArc = anArc;
 323
 324       myPoints = anTransArc->GetPoints();
 325       myScalarValues = aScalarValues;
 326       myStatus = Arc_Done;
 327     }
 328   }
 329   else{
 330 #ifdef _MY_DEBUG_
 331     std::cout<<"Points lay on the one line !"<<endl;
 332 #endif
 333     PntList aList;
 334     aList.push_back(thePnt1);
 335     aList.push_back(thePnt2);
 336     aList.push_back(thePnt3);
 337     vtkUnstructuredGrid* aGrid = BuildGrid(aList);
 338     aGrid->Update();
 339     myPoints = aGrid->GetPoints();
 340
 341     myScalarValues.clear();
 342     myScalarValues.push_back(thePnt1.GetScalarValue());
 343     myScalarValues.push_back(thePnt2.GetScalarValue());
 344     myScalarValues.push_back(thePnt3.GetScalarValue());
 345     myStatus = Arc_Done;
 346   }
 347 }
 348
 349 /*!
 350  * Destructor
 351  */
 352 VTKViewer_ArcBuilder::~VTKViewer_ArcBuilder()
 353 {}
 354
 355
 356 /*
 357  * Add to the vtkUnstructuredGrid points from input list
 358  */
 359 vtkUnstructuredGrid* VTKViewer_ArcBuilder::BuildGrid(const PntList& theList) const
 360 {
 361   int aListsize = theList.size();
 362   vtkUnstructuredGrid* aGrid = NULL;
 363
 364   if(aListsize != 0) {
 365     aGrid = vtkUnstructuredGrid::New();
 366     vtkPoints* aPoints = vtkPoints::New();
 367     aPoints->SetNumberOfPoints(aListsize);
 368
 369     aGrid->Allocate(aListsize);
 370
 371     PntList::const_iterator it = theList.begin();
 372     int aCounter = 0;
 373     for(;it != theList.end();it++) {
 374       Pnt aPnt =  *it;
 375       aPoints->InsertPoint(aCounter,
 376                            aPnt.GetXYZ().X(),
 377                            aPnt.GetXYZ().Y(),
 378                            aPnt.GetXYZ().Z());
 379       vtkVertex* aVertex = vtkVertex::New();
 380       aVertex->GetPointIds()->SetId(0, aCounter);
 381       aGrid->InsertNextCell(aVertex->GetCellType(), aVertex->GetPointIds());
 382       aCounter++;
 383       aVertex->Delete();
 384     }
 385     aGrid->SetPoints(aPoints);
 386     aPoints->Delete();
 387   }
 388   return aGrid;
 389 }
 390
 391
 392 vtkUnstructuredGrid*
 393 VTKViewer_ArcBuilder::TransformGrid(vtkUnstructuredGrid* theGrid,
 394                                     const Vec& theAxis, const double angle) const
 395 {
 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();
 403 }
 404
 405
 406 void VTKViewer_ArcBuilder::GetAngle(const double theAngle){
 407   myAngle = theAngle;
 408 }
 409
 410 double InterpolateScalarValue(int index, int count, double firstValue, double middleValue, double lastValue)
 411 {
 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;
 418   return value;
 419 }
 420
 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
 425
 426
 427   theScalarValues.clear();
 428
 429   double K1 = 0;
 430   double K2 = 0;
 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;
 436
 437 #ifdef _MY_DEBUG_
 438   std::cout<<"K1 : "<< K1 <<endl;
 439   std::cout<<"K2 : "<< K2 <<endl;
 440 #endif
 441   double xCenter;
 442   if( !okK2 ) //K2 --> infinity
 443     xCenter = (K1*(y1-y3) + (x1+x2))/2.0;
 444
 445   else if( !okK1 ) //K1 --> infinity
 446     xCenter = (K2*(y1-y3) - (x2+x3))/(-2.0);
 447
 448   else
 449     xCenter = (K1*K2*(y1-y3) + K2*(x1+x2) - K1*(x2+x3))/ (2.0*(K2-K1));
 450
 451   double yCenter;
 452
 453   if ( K1 == 0 )
 454     yCenter =  (-1/K2)*(xCenter - (x2+x3)/2.0) + (y2 + y3)/2.0;
 455   else
 456     yCenter =  (-1/K1)*(xCenter - (x1+x2)/2.0) + (y1 + y2)/2.0;
 457
 458 #ifdef _MY_DEBUG_
 459   double zCenter = z;
 460   std::cout<<"xCenter : "<<xCenter<<endl;
 461   std::cout<<"yCenter : "<<yCenter<<endl;
 462   std::cout<<"zCenter : "<<zCenter<<endl;
 463 #endif
 464   double aRadius = sqrt((x1 - xCenter)*(x1 - xCenter) + (y1 - yCenter)*(y1 - yCenter));
 465
 466   double angle1 = GetPointAngleOnCircle(xCenter,yCenter,x1,y1);
 467   double angle2 = GetPointAngleOnCircle(xCenter,yCenter,x2,y2);
 468   double angle3 = GetPointAngleOnCircle(xCenter,yCenter,x3,y3);
 469
 470
 471   double aMaxAngle = vtkMath::RadiansFromDegrees(1.)*myAngle*2;
 472
 473   /*  double aTotalAngle =  fabs(angle3 - angle1);
 474
 475   if (aTotalAngle > vtkMath::DoublePi())
 476     aTotalAngle = 2*vtkMath::DoublePi()-aTotalAngle;
 477   */
 478
 479   double aTotalAngle = 0;
 480   IncOrder aOrder = GetArcAngle(angle1,angle2,angle3,&aTotalAngle);
 481
 482   vtkUnstructuredGrid *aC = NULL;
 483
 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;
 490     else
 491       aCurrentAngle+=anIncrementAngle;
 492 #ifdef _MY_DEBUG_
 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;
 497 #endif
 498
 499     PntList aList;
 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);
 507
 508       aList.push_back(aPnt);
 509       theScalarValues.push_back(aPnt.GetScalarValue());
 510       if(aOrder == VTKViewer_ArcBuilder::MINUS)
 511         aCurrentAngle-=anIncrementAngle;
 512       else
 513         aCurrentAngle+=anIncrementAngle;
 514     }
 515     aList.push_back(myPnt3);
 516     theScalarValues.push_back(myPnt3.GetScalarValue());
 517
 518     aC = BuildGrid(aList);
 519 #ifdef _MY_DEBUG_
 520   cout<<"angle1 : "<<angle1*vtkMath::DoubleRadiansToDegrees()<<endl;
 521   cout<<"angle2 : "<<angle2*vtkMath::DoubleRadiansToDegrees()<<endl;
 522   cout<<"angle3 : "<<angle3*vtkMath::DoubleRadiansToDegrees()<<endl;
 523 #endif
 524   }
 525   else{
 526     PntList aList;
 527     aList.push_back(myPnt1);
 528     aList.push_back(myPnt2);
 529     aList.push_back(myPnt3);
 530     aC = BuildGrid(aList);
 531
 532     theScalarValues.push_back(myPnt1.GetScalarValue());
 533     theScalarValues.push_back(myPnt2.GetScalarValue());
 534     theScalarValues.push_back(myPnt3.GetScalarValue());
 535   }
 536   return aC;
 537 }
 538
 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);
 543   if(result < 0 )
 544     result = result+vtkMath::DoublePi()*2;
 545   return vtkMath::DoublePi()*2-result;
 546   return result;
 547 }
 548
 549 vtkPoints* VTKViewer_ArcBuilder::GetPoints(){
 550   return myPoints;
 551 }
 552
 553 const std::vector<double>& VTKViewer_ArcBuilder::GetScalarValues()
 554 {
 555   return myScalarValues;
 556 }
 557
 558 VTKViewer_ArcBuilder::IncOrder VTKViewer_ArcBuilder::GetArcAngle( const double& P1, const double& P2, const double& P3,double* Ang){
 559   IncOrder aResult;
 560   if(P1 < P2 && P2 < P3){
 561     *Ang = P3 - P1;
 562     aResult = VTKViewer_ArcBuilder::PLUS;
 563   }
 564   else if((P1 < P3 && P3 < P2) || (P2 < P1 && P1 < P3)){
 565     *Ang = 2*vtkMath::DoublePi() - P3 + P1;
 566     aResult = VTKViewer_ArcBuilder::MINUS;
 567   }
 568   else if((P2 < P3 && P3 < P1) || (P3 < P1 && P1 < P2)){
 569     *Ang = 2*vtkMath::DoublePi() - P1 + P3;
 570     aResult = VTKViewer_ArcBuilder::PLUS;
 571   }
 572   else if(P3 < P2 && P2 < P1){
 573     *Ang = P1 - P3;
 574     aResult = VTKViewer_ArcBuilder::MINUS;
 575   }
 576   return aResult;
 577 }
 578
 579 //------------------------------------------------------------------------
 580 Pnt CreatePnt(vtkCell* cell, vtkDataArray* scalars, vtkIdType index)
 581 {
 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);
 587   return point;
 588 }
 589
 590 //------------------------------------------------------------------------
 591 vtkIdType Build1DArc(vtkIdType cellId, vtkUnstructuredGrid* input,
 592                      vtkPolyData *output,vtkIdType *pts,
 593                      vtkFloatingPointType myMaxArcAngle){
 594
 595   vtkIdType aResult = -1;
 596   vtkIdType *aNewPoints;
 597
 598   vtkDataArray* inputScalars = input->GetPointData()->GetScalars();
 599   vtkDataArray* outputScalars = output->GetPointData()->GetScalars();
 600
 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 );
 606
 607   VTKViewer_ArcBuilder aBuilder(P0,P2,P1,myMaxArcAngle);
 608   if (aBuilder.GetStatus() != VTKViewer_ArcBuilder::Arc_Done) {
 609     return aResult;
 610   }
 611   else{
 612     vtkPoints* aPoints = aBuilder.GetPoints();
 613     std::vector<double> aScalarValues = aBuilder.GetScalarValues();
 614     vtkIdType aNbPts = aPoints->GetNumberOfPoints();
 615     aNewPoints = new vtkIdType[aNbPts];
 616     vtkIdType curID;
 617     vtkIdType aCellType = VTK_POLY_LINE;
 618
 619     aNewPoints[0] = pts[0];
 620     for(vtkIdType idx = 1; idx < aNbPts-1;idx++) {
 621       curID = output->GetPoints()->InsertNextPoint(aPoints->GetPoint(idx));
 622       if( outputScalars )
 623         outputScalars->InsertNextTuple1(aScalarValues[idx]);
 624       aNewPoints[idx] = curID;
 625     }
 626     aNewPoints[aNbPts-1] = pts[1];
 627
 628     aResult = output->InsertNextCell(aCellType,aNbPts,aNewPoints);
 629     return aResult;
 630   }
 631 }
 632
 633 /*!
 634  * Add all points from the input vector theCollection into thePoints.
 635  * Array theIds - it is array with ids of added points.
 636  */
 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,
 641                          vtkIdType* &theIds){
 642   vtkIdType aNbPoints = 0;
 643   vtkIdType anIdCounter = 0;
 644   vtkIdType aNewPntId = 0;
 645
 646   //Compute number of points
 647   std::vector<vtkPoints*>::const_iterator it = theCollection.begin();
 648   for(;it != theCollection.end();it++){
 649     vtkPoints* aPoints = *it;
 650     if(aPoints) {
 651       aNbPoints += aPoints->GetNumberOfPoints()-1; //Because we add all points except last
 652     }
 653   }
 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;
 661
 662     if(aPoints){
 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];
 667         anIdCounter++;
 668       }
 669     }
 670   }
 671   return aNbPoints;
 672 }