src/VTKViewer/VTKViewer_ArcBuilder.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF, 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, or (at your option) any later version.
   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   myAngle(theAngle),
 231   myStatus(Arc_Error)
 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 = (int)theList.size(); //!< TODO: conversion from size_t to int
 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   aGrid->Register(0);
 404   aTransformFilter->Delete();
 405   return aGrid;
 406 }
 407
 408
 409 void VTKViewer_ArcBuilder::GetAngle(const double theAngle){
 410   myAngle = theAngle;
 411 }
 412
 413 double InterpolateScalarValue(int index, int count, double firstValue, double middleValue, double lastValue)
 414 {
 415   bool isFirstHalf = index <= count / 2;
 416   double first = isFirstHalf ? firstValue : lastValue;
 417   double last = isFirstHalf ? middleValue : middleValue;
 418   double ratio = (double)index / (double)count;
 419   double position = isFirstHalf ? ratio * 2 : ( 1 - ratio ) * 2;
 420   double value = first + (last - first) * position;
 421   return value;
 422 }
 423
 424 vtkUnstructuredGrid* VTKViewer_ArcBuilder::BuildArc(std::vector<double>& theScalarValues)
 425 {
 426   double x1 = myPnt1.GetXYZ().X(); double x2 = myPnt2.GetXYZ().X(); double x3 = myPnt3.GetXYZ().X();
 427   double y1 = myPnt1.GetXYZ().Y(); double y2 = myPnt2.GetXYZ().Y(); double y3 = myPnt3.GetXYZ().Y();
 428   double z =  myPnt1.GetXYZ().Z();  //Points on plane || XOY
 429
 430
 431   theScalarValues.clear();
 432
 433   double K1 = 0;
 434   double K2 = 0;
 435   bool okK1 = false, okK2 = false;
 436   if ( fabs(x2 - x1) > DBL_MIN )
 437     K1 = (y2 - y1)/(x2 - x1), okK1 = true;
 438   if ( fabs(x3 - x2) > DBL_MIN )
 439     K2 = (y3 - y2)/(x3 - x2), okK2 = true;
 440
 441 #ifdef _MY_DEBUG_
 442   std::cout<<"K1 : "<< K1 <<endl;
 443   std::cout<<"K2 : "<< K2 <<endl;
 444 #endif
 445   double xCenter;
 446   if( !okK2 ) //K2 --> infinity
 447     xCenter = (K1*(y1-y3) + (x1+x2))/2.0;
 448
 449   else if( !okK1 ) //K1 --> infinity
 450     xCenter = (K2*(y1-y3) - (x2+x3))/(-2.0);
 451
 452   else
 453     xCenter = (K1*K2*(y1-y3) + K2*(x1+x2) - K1*(x2+x3))/ (2.0*(K2-K1));
 454
 455   double yCenter;
 456
 457   if ( K1 == 0 )
 458     yCenter =  (-1/K2)*(xCenter - (x2+x3)/2.0) + (y2 + y3)/2.0;
 459   else
 460     yCenter =  (-1/K1)*(xCenter - (x1+x2)/2.0) + (y1 + y2)/2.0;
 461
 462 #ifdef _MY_DEBUG_
 463   double zCenter = z;
 464   std::cout<<"xCenter : "<<xCenter<<endl;
 465   std::cout<<"yCenter : "<<yCenter<<endl;
 466   std::cout<<"zCenter : "<<zCenter<<endl;
 467 #endif
 468   double aRadius = sqrt((x1 - xCenter)*(x1 - xCenter) + (y1 - yCenter)*(y1 - yCenter));
 469
 470   double angle1 = GetPointAngleOnCircle(xCenter,yCenter,x1,y1);
 471   double angle2 = GetPointAngleOnCircle(xCenter,yCenter,x2,y2);
 472   double angle3 = GetPointAngleOnCircle(xCenter,yCenter,x3,y3);
 473
 474
 475   double aMaxAngle = vtkMath::RadiansFromDegrees(1.)*myAngle*2;
 476
 477   /*  double aTotalAngle =  fabs(angle3 - angle1);
 478
 479   if (aTotalAngle > vtkMath::Pi())
 480     aTotalAngle = 2*vtkMath::Pi()-aTotalAngle;
 481   */
 482
 483   double aTotalAngle = 0;
 484   IncOrder aOrder = GetArcAngle(angle1,angle2,angle3,&aTotalAngle);
 485
 486   vtkUnstructuredGrid *aC = NULL;
 487
 488   if(aTotalAngle > aMaxAngle) {
 489     int nbSteps = int(aTotalAngle/aMaxAngle)+1;
 490     double anIncrementAngle = aTotalAngle/nbSteps;
 491     double aCurrentAngle = angle1;
 492     if(aOrder == VTKViewer_ArcBuilder::MINUS)
 493       aCurrentAngle-=anIncrementAngle;
 494     else
 495       aCurrentAngle+=anIncrementAngle;
 496 #ifdef _MY_DEBUG_
 497     cout<<"Total angle :"<<aTotalAngle<<endl;
 498     cout<<"Max Increment Angle :"<<aMaxAngle<<endl;
 499     cout<<"Real Increment angle :"<<anIncrementAngle<<endl;
 500     cout<<"Nb Steps : "<<nbSteps<<endl;
 501 #endif
 502
 503     PntList aList;
 504     aList.push_back(myPnt1);
 505     theScalarValues.push_back(myPnt1.GetScalarValue());
 506     for(int i=1;i<=nbSteps-1;i++){
 507       double x = xCenter + aRadius*cos(aCurrentAngle);
 508       double y = yCenter + aRadius*sin(aCurrentAngle);
 509       double value = InterpolateScalarValue(i, nbSteps, myPnt1.GetScalarValue(), myPnt2.GetScalarValue(), myPnt3.GetScalarValue());
 510       Pnt aPnt(x,y,z,value);
 511
 512       aList.push_back(aPnt);
 513       theScalarValues.push_back(aPnt.GetScalarValue());
 514       if(aOrder == VTKViewer_ArcBuilder::MINUS)
 515         aCurrentAngle-=anIncrementAngle;
 516       else
 517         aCurrentAngle+=anIncrementAngle;
 518     }
 519     aList.push_back(myPnt3);
 520     theScalarValues.push_back(myPnt3.GetScalarValue());
 521
 522     aC = BuildGrid(aList);
 523 #ifdef _MY_DEBUG_
 524   cout<<"angle1 : "<<angle1*vtkMath::DoubleRadiansToDegrees()<<endl;
 525   cout<<"angle2 : "<<angle2*vtkMath::DoubleRadiansToDegrees()<<endl;
 526   cout<<"angle3 : "<<angle3*vtkMath::DoubleRadiansToDegrees()<<endl;
 527 #endif
 528   }
 529   else{
 530     PntList aList;
 531     aList.push_back(myPnt1);
 532     aList.push_back(myPnt2);
 533     aList.push_back(myPnt3);
 534     aC = BuildGrid(aList);
 535
 536     theScalarValues.push_back(myPnt1.GetScalarValue());
 537     theScalarValues.push_back(myPnt2.GetScalarValue());
 538     theScalarValues.push_back(myPnt3.GetScalarValue());
 539   }
 540   return aC;
 541 }
 542
 543 double VTKViewer_ArcBuilder::
 544 GetPointAngleOnCircle(const double theXCenter, const double theYCenter,
 545                       const double theXPoint, const double theYPoint){
 546   double result = atan2(theYCenter - theYPoint, theXPoint - theXCenter);
 547   if(result < 0 )
 548     result = result+vtkMath::Pi()*2;
 549   return vtkMath::Pi()*2-result;
 550   return result;
 551 }
 552
 553 vtkPoints* VTKViewer_ArcBuilder::GetPoints(){
 554   return myPoints;
 555 }
 556
 557 const std::vector<double>& VTKViewer_ArcBuilder::GetScalarValues()
 558 {
 559   return myScalarValues;
 560 }
 561
 562 VTKViewer_ArcBuilder::IncOrder VTKViewer_ArcBuilder::GetArcAngle( const double& P1, const double& P2, const double& P3,double* Ang){
 563   IncOrder aResult = MINUS;
 564   if(P1 < P2 && P2 < P3){
 565     *Ang = P3 - P1;
 566     aResult = PLUS;
 567   }
 568   else if((P1 < P3 && P3 < P2) || (P2 < P1 && P1 < P3)){
 569     *Ang = 2*vtkMath::Pi() - P3 + P1;
 570     aResult = MINUS;
 571   }
 572   else if((P2 < P3 && P3 < P1) || (P3 < P1 && P1 < P2)){
 573     *Ang = 2*vtkMath::Pi() - P1 + P3;
 574     aResult = PLUS;
 575   }
 576   else if(P3 < P2 && P2 < P1){
 577     *Ang = P1 - P3;
 578     aResult = MINUS;
 579   }
 580   return aResult;
 581 }
 582
 583 //------------------------------------------------------------------------
 584 Pnt CreatePnt(vtkCell* cell, vtkDataArray* scalars, vtkIdType index)
 585 {
 586   double coord[3];
 587   cell->GetPoints()->GetPoint(index, coord);
 588   vtkIdType pointId = cell->GetPointId(index);
 589   double scalarValue = scalars ? scalars->GetTuple1(pointId) : 0;
 590   Pnt point(coord[0], coord[1], coord[2], scalarValue);
 591   return point;
 592 }
 593
 594 //------------------------------------------------------------------------
 595 vtkIdType Build1DArc(vtkIdType cellId, vtkUnstructuredGrid* input,
 596                      vtkPolyData *output,vtkIdType *pts,
 597                      double myMaxArcAngle){
 598
 599   vtkIdType aResult = -1;
 600
 601   vtkDataArray* inputScalars = input->GetPointData()->GetScalars();
 602   vtkDataArray* outputScalars = output->GetPointData()->GetScalars();
 603
 604   vtkCell* aCell = input->GetCell(cellId);
 605   //Get All points from input cell
 606   Pnt P0 = CreatePnt( aCell, inputScalars, 0 );
 607   Pnt P1 = CreatePnt( aCell, inputScalars, 1 );
 608   Pnt P2 = CreatePnt( aCell, inputScalars, 2 );
 609
 610   VTKViewer_ArcBuilder aBuilder(P0,P2,P1,myMaxArcAngle);
 611   if (aBuilder.GetStatus() != VTKViewer_ArcBuilder::Arc_Done) {
 612     return aResult;
 613   }
 614   else{
 615     vtkPoints* aPoints = aBuilder.GetPoints();
 616     std::vector<double> aScalarValues = aBuilder.GetScalarValues();
 617     vtkIdType aNbPts = aPoints->GetNumberOfPoints();
 618     std::vector< vtkIdType > aNewPoints( aNbPts );
 619     vtkIdType curID;
 620     vtkIdType aCellType = VTK_POLY_LINE;
 621
 622     aNewPoints[0] = pts[0];
 623     for(vtkIdType idx = 1; idx < aNbPts-1;idx++) {
 624       curID = output->GetPoints()->InsertNextPoint(aPoints->GetPoint(idx));
 625       if( outputScalars )
 626         outputScalars->InsertNextTuple1(aScalarValues[idx]);
 627       aNewPoints[idx] = curID;
 628     }
 629     aNewPoints[aNbPts-1] = pts[1];
 630
 631     aResult = output->InsertNextCell(aCellType,aNbPts,&aNewPoints[0]);
 632     return aResult;
 633   }
 634 }
 635
 636 /*!
 637  * Add all points from the input vector theCollection into thePoints.
 638  * Array theIds - it is array with ids of added points.
 639  */
 640 vtkIdType MergevtkPoints(const std::vector< vtkSmartPointer< vtkPoints > >& theCollection,
 641                          const std::vector< std::vector<double> >& theScalarCollection,
 642                          vtkPoints* thePoints,
 643                          std::map<vtkIdType, double>& thePntId2ScalarValue,
 644                          vtkIdType* &theIds){
 645   vtkIdType aNbPoints = 0;
 646   vtkIdType anIdCounter = 0;
 647   vtkIdType aNewPntId = 0;
 648
 649   //Compute number of points
 650   std::vector< vtkSmartPointer< vtkPoints > >::const_iterator it = theCollection.begin();
 651   for(;it != theCollection.end();it++){
 652     vtkPoints* aPoints = *it;
 653     if(aPoints) {
 654       aNbPoints += aPoints->GetNumberOfPoints()-1; //Because we add all points except last
 655     }
 656   }
 657   it = theCollection.begin();
 658   std::vector< std::vector<double> >::const_iterator itScalar = theScalarCollection.begin();
 659   theIds = new vtkIdType[aNbPoints];
 660   // ..and add all points
 661   for(;it != theCollection.end() && itScalar != theScalarCollection.end(); it++, itScalar++){
 662     vtkPoints* aPoints = *it;
 663     std::vector<double> aScalarValues = *itScalar;
 664
 665     if(aPoints){
 666       for(vtkIdType idx = 0;idx < aPoints->GetNumberOfPoints()-1;idx++){
 667         aNewPntId = thePoints->InsertNextPoint(aPoints->GetPoint(idx));
 668         theIds[anIdCounter] = aNewPntId;
 669         thePntId2ScalarValue[ aNewPntId ] = aScalarValues[idx];
 670         anIdCounter++;
 671       }
 672     }
 673   }
 674   return aNbPoints;
 675 }