]> SALOME platform Git repositories - modules/gui.git/blob - src/VTKViewer/VTKViewer_ArcBuilder.cxx
Salome HOME
8cc110c7c5421b5861c1634c8560f559cc227e7d
[modules/gui.git] / src / VTKViewer / VTKViewer_ArcBuilder.cxx
1 // Copyright (C) 2007-2015  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, 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   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   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;
564   if(P1 < P2 && P2 < P3){
565     *Ang = P3 - P1;
566     aResult = VTKViewer_ArcBuilder::PLUS;
567   }
568   else if((P1 < P3 && P3 < P2) || (P2 < P1 && P1 < P3)){
569     *Ang = 2*vtkMath::Pi() - P3 + P1;
570     aResult = VTKViewer_ArcBuilder::MINUS;
571   }
572   else if((P2 < P3 && P3 < P1) || (P3 < P1 && P1 < P2)){
573     *Ang = 2*vtkMath::Pi() - P1 + P3;
574     aResult = VTKViewer_ArcBuilder::PLUS;
575   }
576   else if(P3 < P2 && P2 < P1){
577     *Ang = P1 - P3;
578     aResult = VTKViewer_ArcBuilder::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   vtkIdType *aNewPoints;
601
602   vtkDataArray* inputScalars = input->GetPointData()->GetScalars();
603   vtkDataArray* outputScalars = output->GetPointData()->GetScalars();
604
605   vtkCell* aCell = input->GetCell(cellId);
606   //Get All points from input cell
607   Pnt P0 = CreatePnt( aCell, inputScalars, 0 );
608   Pnt P1 = CreatePnt( aCell, inputScalars, 1 );
609   Pnt P2 = CreatePnt( aCell, inputScalars, 2 );
610
611   VTKViewer_ArcBuilder aBuilder(P0,P2,P1,myMaxArcAngle);
612   if (aBuilder.GetStatus() != VTKViewer_ArcBuilder::Arc_Done) {
613     return aResult;
614   }
615   else{
616     vtkPoints* aPoints = aBuilder.GetPoints();
617     std::vector<double> aScalarValues = aBuilder.GetScalarValues();
618     vtkIdType aNbPts = aPoints->GetNumberOfPoints();
619     std::vector< vtkIdType > aNewPoints( aNbPts );
620     vtkIdType curID;
621     vtkIdType aCellType = VTK_POLY_LINE;
622
623     aNewPoints[0] = pts[0];
624     for(vtkIdType idx = 1; idx < aNbPts-1;idx++) {
625       curID = output->GetPoints()->InsertNextPoint(aPoints->GetPoint(idx));
626       if( outputScalars )
627         outputScalars->InsertNextTuple1(aScalarValues[idx]);
628       aNewPoints[idx] = curID;
629     }
630     aNewPoints[aNbPts-1] = pts[1];
631
632     aResult = output->InsertNextCell(aCellType,aNbPts,&aNewPoints[0]);
633     return aResult;
634   }
635 }
636
637 /*!
638  * Add all points from the input vector theCollection into thePoints. 
639  * Array theIds - it is array with ids of added points.
640  */
641 vtkIdType MergevtkPoints(const std::vector< vtkSmartPointer< vtkPoints > >& theCollection,
642                          const std::vector< std::vector<double> >& theScalarCollection,
643                          vtkPoints* thePoints,
644                          std::map<int, double>& thePntId2ScalarValue,
645                          vtkIdType* &theIds){
646   vtkIdType aNbPoints = 0;
647   vtkIdType anIdCounter = 0;
648   vtkIdType aNewPntId = 0;
649
650   //Compute number of points
651   std::vector< vtkSmartPointer< vtkPoints > >::const_iterator it = theCollection.begin();
652   for(;it != theCollection.end();it++){
653     vtkPoints* aPoints = *it;
654     if(aPoints) { 
655       aNbPoints += aPoints->GetNumberOfPoints()-1; //Because we add all points except last
656     }
657   }
658   it = theCollection.begin();
659   std::vector< std::vector<double> >::const_iterator itScalar = theScalarCollection.begin();
660   theIds = new vtkIdType[aNbPoints];
661   // ..and add all points
662   for(;it != theCollection.end() && itScalar != theScalarCollection.end(); it++, itScalar++){
663     vtkPoints* aPoints = *it;
664     std::vector<double> aScalarValues = *itScalar;
665     
666     if(aPoints){
667       for(vtkIdType idx = 0;idx < aPoints->GetNumberOfPoints()-1;idx++){
668         aNewPntId = thePoints->InsertNextPoint(aPoints->GetPoint(idx));
669         theIds[anIdCounter] = aNewPntId;
670         thePntId2ScalarValue[ aNewPntId ] = aScalarValues[idx];
671         anIdCounter++;
672       }
673     }
674   }
675   return aNbPoints;
676 }