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 <cmath>
  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 #include <limits>
  46 #ifdef _MY_DEBUG_
  47 #include <iostream>
  48 #endif
  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 > std::numeric_limits<double>::min() ) {
 119     coord.SetX(Xv / D);
 120     coord.SetY(Yv / D);
 121     coord.SetZ(Zv / D);
 122   }
 123   else {
 124     coord.SetX(0);
 125     coord.SetY(0);
 126     coord.SetZ(0);
 127   }
 128 }
 129
 130 /*!
 131  * Destructor
 132  */
 133 Vec::~Vec(){}
 134
 135
 136 /*!
 137  * Calculate angle between vectors in radians
 138  */
 139 double Vec::AngleBetween(const Vec & Other)
 140 {
 141   double res;
 142   double numerator = GetXYZ().X()*Other.GetXYZ().X()+GetXYZ().Y()*Other.GetXYZ().Y()+GetXYZ().Z()*Other.GetXYZ().Z();
 143   double denumerator = GetXYZ().Modulus()*Other.GetXYZ().Modulus();
 144   double d = numerator/denumerator;
 145   if( d < -1 && d > -1 - PRECISION )
 146     d = -1;
 147   else if( d > 1 && d < 1 + PRECISION )
 148     d = 1;
 149   res = acos( d );
 150   return res;
 151 }
 152
 153 /*!
 154  * Calculate angle between vectors in degrees
 155  */
 156 double Vec::AngleBetweenInGrad(const Vec & Other){
 157   return AngleBetween(Other)*vtkMath::DegreesFromRadians(1.);
 158 }
 159
 160 /*
 161  * Calculate vector multiplication
 162 */
 163 Vec Vec::VectMultiplication(const Vec & Other) const{
 164   double x = GetXYZ().Y()*Other.GetXYZ().Z() - GetXYZ().Z()*Other.GetXYZ().Y();
 165   double y = GetXYZ().Z()*Other.GetXYZ().X() - GetXYZ().X()*Other.GetXYZ().Z();
 166   double z = GetXYZ().X()*Other.GetXYZ().Y() - GetXYZ().Y()*Other.GetXYZ().X();
 167   return Vec(x,y,z);
 168 }
 169
 170 /*---------------------Class Plane --------------------------------*/
 171 /*!
 172  * Constructor
 173  */
 174 Plane::Plane(const Pnt& thePnt1, const Pnt& thePnt2, const Pnt& thePnt3)
 175 {
 176   CalculatePlane(thePnt1,thePnt2,thePnt3);
 177 }
 178
 179 /*!
 180  * Destructor
 181  */
 182 Plane::~Plane()
 183 {}
 184
 185 /*
 186  * Return plane normale
 187  */
 188 Vec Plane::GetNormale() const{
 189   return Vec(myA,myB,myC);
 190 }
 191
 192 /*
 193  * Calculate A,B,C coeefs of plane
 194  */
 195 void Plane::CalculatePlane(const Pnt& thePnt1, const Pnt& thePnt2, const Pnt& thePnt3){
 196
 197   double x1 = thePnt1.GetXYZ().X();
 198   double x2 = thePnt2.GetXYZ().X();
 199   double x3 = thePnt3.GetXYZ().X();
 200
 201   double y1 = thePnt1.GetXYZ().Y();
 202   double y2 = thePnt2.GetXYZ().Y();
 203   double y3 = thePnt3.GetXYZ().Y();
 204
 205   double z1 = thePnt1.GetXYZ().Z();
 206   double z2 = thePnt2.GetXYZ().Z();
 207   double z3 = thePnt3.GetXYZ().Z();
 208
 209   myA = y1*(z2 - z3) + y2*(z3 - z1) + y3*(z1 - z2);
 210   myB = z1*(x2 - x3) + z2*(x3 - x1) + z3*(x1 - x2);
 211   myC = x1*(y2 - y3) + x2*(y3 - y1) + x3*(y1 - y2);
 212
 213 #ifdef _MY_DEBUG_
 214   cout<<"Plane A: "<<myA<<endl;
 215   cout<<"Plane B: "<<myB<<endl;
 216   cout<<"Plane C: "<<myC<<endl;
 217 #endif
 218
 219 }
 220
 221
 222 /*---------------------Class VTKViewer_ArcBuilder --------------------------------*/
 223 /*!
 224  * Constructor
 225  */
 226 VTKViewer_ArcBuilder::VTKViewer_ArcBuilder(const Pnt& thePnt1,
 227                                            const Pnt& thePnt2,
 228                                            const Pnt& thePnt3,
 229                                            double theAngle):
 230   myStatus(Arc_Error),
 231   myAngle(theAngle)
 232 {
 233   Vec V1(thePnt2.GetXYZ().X()-thePnt1.GetXYZ().X(),
 234          thePnt2.GetXYZ().Y()-thePnt1.GetXYZ().Y(),
 235          thePnt2.GetXYZ().Z()-thePnt1.GetXYZ().Z());
 236
 237   Vec V2(thePnt2.GetXYZ().X()-thePnt3.GetXYZ().X(),
 238          thePnt2.GetXYZ().Y()-thePnt3.GetXYZ().Y(),
 239          thePnt2.GetXYZ().Z()-thePnt3.GetXYZ().Z());
 240
 241   double angle = V1.AngleBetweenInGrad(V2);
 242
 243   //Check that points are not belong one line
 244 #ifdef _MY_DEBUG_
 245   cout<<"Angle for check: "<<angle<<endl;
 246 #endif
 247   if(CheckAngle(180,angle)) {
 248
 249     // Build plane by three points
 250     Plane aPlane(thePnt1,thePnt2,thePnt3);
 251
 252     //Plane normales
 253     Vec aPlaneNormale = aPlane.GetNormale();
 254 #ifdef _MY_DEBUG_
 255     std::cout<<"X Normale: "<<aPlaneNormale.GetXYZ().X()<<std::endl;
 256     std::cout<<"Y Normale: "<<aPlaneNormale.GetXYZ().Y()<<std::endl;
 257     std::cout<<"Z Normale: "<<aPlaneNormale.GetXYZ().Z()<<std::endl;
 258 #endif
 259     //OZ vector
 260     Vec OZ(0,0,1);
 261
 262     //Projection plane normale on XOY
 263     Vec aNormaleProjection (aPlaneNormale.GetXYZ().X(),
 264                             aPlaneNormale.GetXYZ().Y(),
 265                             0);
 266 #ifdef _MY_DEBUG_
 267     std::cout<<"X Normale Projection: "<<aNormaleProjection.GetXYZ().X()<<std::endl;
 268     std::cout<<"Y Normale Projection: "<<aNormaleProjection.GetXYZ().Y()<<std::endl;
 269     std::cout<<"Z Normale Projection: "<<aNormaleProjection.GetXYZ().Z()<<std::endl;
 270 #endif
 271
 272     //Rotation axis
 273     Vec aAxis = aNormaleProjection.VectMultiplication(OZ);
 274 #ifdef _MY_DEBUG_
 275     std::cout<<"X Axis: "<<aAxis.GetXYZ().X()<<std::endl;
 276     std::cout<<"Y Axis: "<<aAxis.GetXYZ().Y()<<std::endl;
 277     std::cout<<"Z Axis: "<<aAxis.GetXYZ().Z()<<std::endl;
 278 #endif
 279     //Rotation angle
 280     double anAngle = OZ.AngleBetweenInGrad(aPlaneNormale);
 281 #ifdef _MY_DEBUG_
 282     std::cout<<"Rotation Angle :"<<anAngle<<endl;
 283 #endif
 284     PntList aInputPnts;
 285     aInputPnts.push_back(thePnt1);
 286     aInputPnts.push_back(thePnt2);
 287     aInputPnts.push_back(thePnt3);
 288
 289     vtkSmartPointer<vtkUnstructuredGrid> aGrid = BuildGrid(aInputPnts);
 290     aGrid->Delete();
 291
 292     bool needRotation = true;
 293     if(anAngle  == 0 || anAngle == 180)
 294       needRotation = false;
 295
 296     if(aGrid) {
 297       if(needRotation) {
 298         aGrid = TransformGrid(aGrid,aAxis,anAngle);
 299         aGrid->Delete();
 300 #ifdef _MY_DEBUG_
 301         cout<<"Need Rotation!!!"<<endl;
 302 #endif
 303       }
 304       else {
 305 #ifdef _MY_DEBUG_
 306         cout<<"Rotation does not need!!!"<<endl;
 307 #endif
 308       }
 309
 310       double coords[3];
 311       aGrid->GetPoint(0,coords);
 312       myPnt1 = Pnt(coords[0],coords[1],coords[2], thePnt1.GetScalarValue());
 313       aGrid->GetPoint(1,coords);
 314       myPnt2 = Pnt(coords[0],coords[1],coords[2], thePnt2.GetScalarValue());
 315       aGrid->GetPoint(2,coords);
 316       myPnt3 = Pnt(coords[0],coords[1],coords[2], thePnt3.GetScalarValue());
 317
 318       vtkSmartPointer<vtkUnstructuredGrid> anArc = BuildArc(myScalarValues);
 319       anArc->Delete();
 320       if(needRotation) {
 321         anArc = TransformGrid(anArc,aAxis,-anAngle);
 322         anArc->Delete();
 323       }
 324       myPoints = anArc->GetPoints();
 325       myStatus = Arc_Done;
 326     }
 327   }
 328   else{
 329 #ifdef _MY_DEBUG_
 330     std::cout<<"Points lay on the one line !"<<endl;
 331 #endif
 332     PntList aList;
 333     aList.push_back(thePnt1);
 334     aList.push_back(thePnt2);
 335     aList.push_back(thePnt3);
 336     vtkUnstructuredGrid* aGrid = BuildGrid(aList);
 337     myPoints = aGrid->GetPoints();
 338     aGrid->Delete();
 339
 340     myScalarValues.clear();
 341     myScalarValues.push_back(thePnt1.GetScalarValue());
 342     myScalarValues.push_back(thePnt2.GetScalarValue());
 343     myScalarValues.push_back(thePnt3.GetScalarValue());
 344     myStatus = Arc_Done;
 345   }
 346 }
 347
 348 /*!
 349  * Destructor
 350  */
 351 VTKViewer_ArcBuilder::~VTKViewer_ArcBuilder()
 352 {}
 353
 354
 355 /*
 356  * Add to the vtkUnstructuredGrid points from input list
 357  */
 358 vtkUnstructuredGrid* VTKViewer_ArcBuilder::BuildGrid(const PntList& theList) const
 359 {
 360   int aListsize = theList.size();
 361   vtkUnstructuredGrid* aGrid = NULL;
 362
 363   if(aListsize != 0) {
 364     aGrid = vtkUnstructuredGrid::New();
 365     vtkPoints* aPoints = vtkPoints::New();
 366     aPoints->SetNumberOfPoints(aListsize);
 367
 368     aGrid->Allocate(aListsize);
 369
 370     PntList::const_iterator it = theList.begin();
 371     int aCounter = 0;
 372     for(;it != theList.end();it++) {
 373       Pnt aPnt =  *it;
 374       aPoints->InsertPoint(aCounter,
 375                            aPnt.GetXYZ().X(),
 376                            aPnt.GetXYZ().Y(),
 377                            aPnt.GetXYZ().Z());
 378       vtkVertex* aVertex = vtkVertex::New();
 379       aVertex->GetPointIds()->SetId(0, aCounter);
 380       aGrid->InsertNextCell(aVertex->GetCellType(), aVertex->GetPointIds());
 381       aCounter++;
 382       aVertex->Delete();
 383     }
 384     aGrid->SetPoints(aPoints);
 385     aPoints->Delete();
 386   }
 387   return aGrid;
 388 }
 389
 390
 391 vtkUnstructuredGrid*
 392 VTKViewer_ArcBuilder::TransformGrid(vtkUnstructuredGrid* theGrid,
 393                                     const Vec& theAxis, const double angle) const
 394 {
 395   vtkTransform *aTransform = vtkTransform::New();
 396   aTransform->RotateWXYZ(angle, theAxis.GetXYZ().X(), theAxis.GetXYZ().Y(), theAxis.GetXYZ().Z());
 397   vtkTransformFilter* aTransformFilter  = vtkTransformFilter::New();
 398   aTransformFilter->SetTransform(aTransform);
 399   aTransform->Delete();
 400   aTransformFilter->SetInputData(theGrid);
 401   aTransformFilter->Update();
 402   vtkUnstructuredGrid * aGrid = aTransformFilter->GetUnstructuredGridOutput();
 403   return aGrid;
 404 }
 405
 406
 407 void VTKViewer_ArcBuilder::GetAngle(const double theAngle){
 408   myAngle = theAngle;
 409 }
 410
 411 double InterpolateScalarValue(int index, int count, double firstValue, double middleValue, double lastValue)
 412 {
 413   bool isFirstHalf = index <= count / 2;
 414   double first = isFirstHalf ? firstValue : lastValue;
 415   double last = isFirstHalf ? middleValue : middleValue;
 416   double ratio = (double)index / (double)count;
 417   double position = isFirstHalf ? ratio * 2 : ( 1 - ratio ) * 2;
 418   double value = first + (last - first) * position;
 419   return value;
 420 }
 421
 422 vtkUnstructuredGrid* VTKViewer_ArcBuilder::BuildArc(std::vector<double>& theScalarValues)
 423 {
 424   double x1 = myPnt1.GetXYZ().X(); double x2 = myPnt2.GetXYZ().X(); double x3 = myPnt3.GetXYZ().X();
 425   double y1 = myPnt1.GetXYZ().Y(); double y2 = myPnt2.GetXYZ().Y(); double y3 = myPnt3.GetXYZ().Y();
 426   double z =  myPnt1.GetXYZ().Z();  //Points on plane || XOY
 427
 428
 429   theScalarValues.clear();
 430
 431   double K1 = 0;
 432   double K2 = 0;
 433   bool okK1 = false, okK2 = false;
 434   if ( fabs(x2 - x1) > DBL_MIN )
 435     K1 = (y2 - y1)/(x2 - x1), okK1 = true;
 436   if ( fabs(x3 - x2) > DBL_MIN )
 437     K2 = (y3 - y2)/(x3 - x2), okK2 = true;
 438
 439 #ifdef _MY_DEBUG_
 440   std::cout<<"K1 : "<< K1 <<endl;
 441   std::cout<<"K2 : "<< K2 <<endl;
 442 #endif
 443   double xCenter;
 444   if( !okK2 ) //K2 --> infinity
 445     xCenter = (K1*(y1-y3) + (x1+x2))/2.0;
 446
 447   else if( !okK1 ) //K1 --> infinity
 448     xCenter = (K2*(y1-y3) - (x2+x3))/(-2.0);
 449
 450   else
 451     xCenter = (K1*K2*(y1-y3) + K2*(x1+x2) - K1*(x2+x3))/ (2.0*(K2-K1));
 452
 453   double yCenter;
 454
 455   if ( K1 == 0 )
 456     yCenter =  (-1/K2)*(xCenter - (x2+x3)/2.0) + (y2 + y3)/2.0;
 457   else
 458     yCenter =  (-1/K1)*(xCenter - (x1+x2)/2.0) + (y1 + y2)/2.0;
 459
 460 #ifdef _MY_DEBUG_
 461   double zCenter = z;
 462   std::cout<<"xCenter : "<<xCenter<<endl;
 463   std::cout<<"yCenter : "<<yCenter<<endl;
 464   std::cout<<"zCenter : "<<zCenter<<endl;
 465 #endif
 466   double aRadius = sqrt((x1 - xCenter)*(x1 - xCenter) + (y1 - yCenter)*(y1 - yCenter));
 467
 468   double angle1 = GetPointAngleOnCircle(xCenter,yCenter,x1,y1);
 469   double angle2 = GetPointAngleOnCircle(xCenter,yCenter,x2,y2);
 470   double angle3 = GetPointAngleOnCircle(xCenter,yCenter,x3,y3);
 471
 472
 473   double aMaxAngle = vtkMath::RadiansFromDegrees(1.)*myAngle*2;
 474
 475   /*  double aTotalAngle =  fabs(angle3 - angle1);
 476
 477   if (aTotalAngle > vtkMath::Pi())
 478     aTotalAngle = 2*vtkMath::Pi()-aTotalAngle;
 479   */
 480
 481   double aTotalAngle = 0;
 482   IncOrder aOrder = GetArcAngle(angle1,angle2,angle3,&aTotalAngle);
 483
 484   vtkUnstructuredGrid *aC = NULL;
 485
 486   if(aTotalAngle > aMaxAngle) {
 487     int nbSteps = int(aTotalAngle/aMaxAngle)+1;
 488     double anIncrementAngle = aTotalAngle/nbSteps;
 489     double aCurrentAngle = angle1;
 490     if(aOrder == VTKViewer_ArcBuilder::MINUS)
 491       aCurrentAngle-=anIncrementAngle;
 492     else
 493       aCurrentAngle+=anIncrementAngle;
 494 #ifdef _MY_DEBUG_
 495     cout<<"Total angle :"<<aTotalAngle<<endl;
 496     cout<<"Max Increment Angle :"<<aMaxAngle<<endl;
 497     cout<<"Real Increment angle :"<<anIncrementAngle<<endl;
 498     cout<<"Nb Steps : "<<nbSteps<<endl;
 499 #endif
 500
 501     PntList aList;
 502     aList.push_back(myPnt1);
 503     theScalarValues.push_back(myPnt1.GetScalarValue());
 504     for(int i=1;i<=nbSteps-1;i++){
 505       double x = xCenter + aRadius*cos(aCurrentAngle);
 506       double y = yCenter + aRadius*sin(aCurrentAngle);
 507       double value = InterpolateScalarValue(i, nbSteps, myPnt1.GetScalarValue(), myPnt2.GetScalarValue(), myPnt3.GetScalarValue());
 508       Pnt aPnt(x,y,z,value);
 509
 510       aList.push_back(aPnt);
 511       theScalarValues.push_back(aPnt.GetScalarValue());
 512       if(aOrder == VTKViewer_ArcBuilder::MINUS)
 513         aCurrentAngle-=anIncrementAngle;
 514       else
 515         aCurrentAngle+=anIncrementAngle;
 516     }
 517     aList.push_back(myPnt3);
 518     theScalarValues.push_back(myPnt3.GetScalarValue());
 519
 520     aC = BuildGrid(aList);
 521 #ifdef _MY_DEBUG_
 522   cout<<"angle1 : "<<angle1*vtkMath::DoubleRadiansToDegrees()<<endl;
 523   cout<<"angle2 : "<<angle2*vtkMath::DoubleRadiansToDegrees()<<endl;
 524   cout<<"angle3 : "<<angle3*vtkMath::DoubleRadiansToDegrees()<<endl;
 525 #endif
 526   }
 527   else{
 528     PntList aList;
 529     aList.push_back(myPnt1);
 530     aList.push_back(myPnt2);
 531     aList.push_back(myPnt3);
 532     aC = BuildGrid(aList);
 533
 534     theScalarValues.push_back(myPnt1.GetScalarValue());
 535     theScalarValues.push_back(myPnt2.GetScalarValue());
 536     theScalarValues.push_back(myPnt3.GetScalarValue());
 537   }
 538   return aC;
 539 }
 540
 541 double VTKViewer_ArcBuilder::
 542 GetPointAngleOnCircle(const double theXCenter, const double theYCenter,
 543                       const double theXPoint, const double theYPoint){
 544   double result = atan2(theYCenter - theYPoint, theXPoint - theXCenter);
 545   if(result < 0 )
 546     result = result+vtkMath::Pi()*2;
 547   return vtkMath::Pi()*2-result;
 548   return result;
 549 }
 550
 551 vtkPoints* VTKViewer_ArcBuilder::GetPoints(){
 552   return myPoints;
 553 }
 554
 555 const std::vector<double>& VTKViewer_ArcBuilder::GetScalarValues()
 556 {
 557   return myScalarValues;
 558 }
 559
 560 VTKViewer_ArcBuilder::IncOrder VTKViewer_ArcBuilder::GetArcAngle( const double& P1, const double& P2, const double& P3,double* Ang){
 561   IncOrder aResult;
 562   if(P1 < P2 && P2 < P3){
 563     *Ang = P3 - P1;
 564     aResult = VTKViewer_ArcBuilder::PLUS;
 565   }
 566   else if((P1 < P3 && P3 < P2) || (P2 < P1 && P1 < P3)){
 567     *Ang = 2*vtkMath::Pi() - P3 + P1;
 568     aResult = VTKViewer_ArcBuilder::MINUS;
 569   }
 570   else if((P2 < P3 && P3 < P1) || (P3 < P1 && P1 < P2)){
 571     *Ang = 2*vtkMath::Pi() - P1 + P3;
 572     aResult = VTKViewer_ArcBuilder::PLUS;
 573   }
 574   else if(P3 < P2 && P2 < P1){
 575     *Ang = P1 - P3;
 576     aResult = VTKViewer_ArcBuilder::MINUS;
 577   }
 578   return aResult;
 579 }
 580
 581 //------------------------------------------------------------------------
 582 Pnt CreatePnt(vtkCell* cell, vtkDataArray* scalars, vtkIdType index)
 583 {
 584   double coord[3];
 585   cell->GetPoints()->GetPoint(index, coord);
 586   vtkIdType pointId = cell->GetPointId(index);
 587   double scalarValue = scalars ? scalars->GetTuple1(pointId) : 0;
 588   Pnt point(coord[0], coord[1], coord[2], scalarValue);
 589   return point;
 590 }
 591
 592 //------------------------------------------------------------------------
 593 vtkIdType Build1DArc(vtkIdType cellId, vtkUnstructuredGrid* input,
 594                      vtkPolyData *output,vtkIdType *pts,
 595                      double myMaxArcAngle){
 596
 597   vtkIdType aResult = -1;
 598   vtkIdType *aNewPoints;
 599
 600   vtkDataArray* inputScalars = input->GetPointData()->GetScalars();
 601   vtkDataArray* outputScalars = output->GetPointData()->GetScalars();
 602
 603   vtkCell* aCell = input->GetCell(cellId);
 604   //Get All points from input cell
 605   Pnt P0 = CreatePnt( aCell, inputScalars, 0 );
 606   Pnt P1 = CreatePnt( aCell, inputScalars, 1 );
 607   Pnt P2 = CreatePnt( aCell, inputScalars, 2 );
 608
 609   VTKViewer_ArcBuilder aBuilder(P0,P2,P1,myMaxArcAngle);
 610   if (aBuilder.GetStatus() != VTKViewer_ArcBuilder::Arc_Done) {
 611     return aResult;
 612   }
 613   else{
 614     vtkPoints* aPoints = aBuilder.GetPoints();
 615     std::vector<double> aScalarValues = aBuilder.GetScalarValues();
 616     vtkIdType aNbPts = aPoints->GetNumberOfPoints();
 617     std::vector< vtkIdType > aNewPoints( aNbPts );
 618     vtkIdType curID;
 619     vtkIdType aCellType = VTK_POLY_LINE;
 620
 621     aNewPoints[0] = pts[0];
 622     for(vtkIdType idx = 1; idx < aNbPts-1;idx++) {
 623       curID = output->GetPoints()->InsertNextPoint(aPoints->GetPoint(idx));
 624       if( outputScalars )
 625         outputScalars->InsertNextTuple1(aScalarValues[idx]);
 626       aNewPoints[idx] = curID;
 627     }
 628     aNewPoints[aNbPts-1] = pts[1];
 629
 630     aResult = output->InsertNextCell(aCellType,aNbPts,&aNewPoints[0]);
 631     return aResult;
 632   }
 633 }
 634
 635 /*!
 636  * Add all points from the input vector theCollection into thePoints.
 637  * Array theIds - it is array with ids of added points.
 638  */
 639 vtkIdType MergevtkPoints(const std::vector< vtkSmartPointer< vtkPoints > >& theCollection,
 640                          const std::vector< std::vector<double> >& theScalarCollection,
 641                          vtkPoints* thePoints,
 642                          std::map<int, double>& thePntId2ScalarValue,
 643                          vtkIdType* &theIds){
 644   vtkIdType aNbPoints = 0;
 645   vtkIdType anIdCounter = 0;
 646   vtkIdType aNewPntId = 0;
 647
 648   //Compute number of points
 649   std::vector< vtkSmartPointer< vtkPoints > >::const_iterator it = theCollection.begin();
 650   for(;it != theCollection.end();it++){
 651     vtkPoints* aPoints = *it;
 652     if(aPoints) {
 653       aNbPoints += aPoints->GetNumberOfPoints()-1; //Because we add all points except last
 654     }
 655   }
 656   it = theCollection.begin();
 657   std::vector< std::vector<double> >::const_iterator itScalar = theScalarCollection.begin();
 658   theIds = new vtkIdType[aNbPoints];
 659   // ..and add all points
 660   for(;it != theCollection.end() && itScalar != theScalarCollection.end(); it++, itScalar++){
 661     vtkPoints* aPoints = *it;
 662     std::vector<double> aScalarValues = *itScalar;
 663
 664     if(aPoints){
 665       for(vtkIdType idx = 0;idx < aPoints->GetNumberOfPoints()-1;idx++){
 666         aNewPntId = thePoints->InsertNextPoint(aPoints->GetPoint(idx));
 667         theIds[anIdCounter] = aNewPntId;
 668         thePntId2ScalarValue[ aNewPntId ] = aScalarValues[idx];
 669         anIdCounter++;
 670       }
 671     }
 672   }
 673   return aNbPoints;
 674 }