dcc3f096eebf361e093204599b06b7bcc0d0d92f
[occt.git] / src / AppParCurves / AppParCurves_Variational.gxx
1 // File AppParCurves_Variational.gxx
2 // Jeannine PANCIATICI le 06/06/96
3 // Igor FEOKTISTOV 14/12/98 - correction of Approximate() and Init().
4 // Approximation d une MultiLine de points decrite par le tool MLineTool.
5 // avec criteres variationnels
6
7
8 #define No_Standard_RangeError
9 #define No_Standard_OutOfRange
10 #define No_Standard_DimensionError
11 #define No_Standard_ConstructionError
12
13 #include <Standard_Stream.hxx>
14
15 #include <AppParCurves.hxx>
16 #include <AppParCurves_Constraint.hxx>
17 #include <AppParCurves_HArray1OfConstraintCouple.hxx>
18 #include <AppParCurves_Array1OfMultiPoint.hxx>
19 #include <AppParCurves_MultiPoint.hxx>
20 #include <AppParCurves_MultiBSpCurve.hxx>
21 #include <Convert_CompPolynomialToPoles.hxx>
22 #include <gp_Pnt.hxx>
23 #include <gp_Pnt2d.hxx>
24 #include <gp_Vec.hxx>
25 #include <gp_Vec2d.hxx>
26 #include <TColgp_Array1OfPnt.hxx>
27 #include <TColgp_Array1OfPnt2d.hxx>
28 #include <TColgp_Array1OfVec.hxx>
29 #include <TColgp_Array1OfVec2d.hxx>
30 #include <TColStd_Array1OfInteger.hxx>
31 #include <TColStd_HArray1OfInteger.hxx>
32 #include <TColStd_Array1OfReal.hxx>
33 #include <TColStd_HArray1OfReal.hxx>
34 #include <TColStd_HArray2OfReal.hxx>
35 #include <StdFail_NotDone.hxx>
36 #include <Standard_SStream.hxx>
37 #include <Standard_NoSuchObject.hxx>
38 #include <Precision.hxx>
39 //#include <Smoothing.h>
40
41 #if defined(WNT)
42 # include <stdio.h>
43 # include <stdarg.h>
44 #endif  /* WNT */
45
46
47 //
48 //=======================================================================
49 //function : AppParCurves_Variational
50 //purpose  : Initialization of   the   fields.
51 //=======================================================================
52 //
53  AppParCurves_Variational::AppParCurves_Variational(const MultiLine& SSP,
54                                                     const Standard_Integer FirstPoint, 
55                                                     const Standard_Integer LastPoint,
56                                                     const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
57                                                     const Standard_Integer MaxDegree,
58                                                     const Standard_Integer MaxSegment,
59                                                     const GeomAbs_Shape Continuity, 
60                                                     const Standard_Boolean WithMinMax, 
61                                                     const Standard_Boolean WithCutting, 
62                                                     const Standard_Real Tolerance, 
63                                                     const Standard_Integer NbIterations):
64                                         mySSP(SSP),
65                                         myFirstPoint(FirstPoint),
66                                         myLastPoint(LastPoint),
67                                         myConstraints(TheConstraints),
68                                         myMaxDegree(MaxDegree),
69                                         myMaxSegment(MaxSegment),
70                                         myNbIterations(NbIterations),
71                                         myTolerance(Tolerance),
72                                         myContinuity(Continuity),
73                                         myWithMinMax(WithMinMax),
74                                         myWithCutting(WithCutting)
75
76
77 {
78 // Verifications:
79     if (myMaxDegree < 1) Standard_DomainError::Raise();
80     myMaxDegree = Min (30, myMaxDegree);
81 //
82     if (myMaxSegment < 1) Standard_DomainError::Raise();
83 //
84     if (myWithMinMax != 0 && myWithMinMax !=1 ) Standard_DomainError::Raise();
85     if (myWithCutting != 0 && myWithCutting !=1 ) Standard_DomainError::Raise();
86 //
87     myIsOverConstr = Standard_False;
88     myIsCreated    = Standard_False;
89     myIsDone       = Standard_False;  
90     switch (myContinuity) {
91     case GeomAbs_C0:
92        myNivCont=0;
93        break ;
94     case GeomAbs_C1: 
95        myNivCont=1;
96        break ;
97     case GeomAbs_C2:  
98        myNivCont=2;
99        break ;
100     default:
101     Standard_ConstructionError::Raise();
102   }
103 //
104     myNbP2d = ToolLine::NbP2d(SSP); 
105     myNbP3d = ToolLine::NbP3d(SSP);
106     myDimension = 2 * myNbP2d + 3* myNbP3d ;
107 //
108     myPercent[0]=0.4;
109     myPercent[1]=0.2;
110     myPercent[2]=0.4;
111     myKnots= new TColStd_HArray1OfReal(1,2);
112     myKnots->SetValue(1,0.);
113     myKnots->SetValue(2,1.);
114     
115 //  Declaration
116 //  
117     mySmoothCriterion = new AppParCurves_MyCriterion(mySSP, myFirstPoint, myLastPoint);
118     myParameters = new TColStd_HArray1OfReal(myFirstPoint, myLastPoint);
119     myNbPoints=myLastPoint-myFirstPoint+1;
120     if (myNbPoints <= 0) Standard_ConstructionError::Raise();
121 // 
122     myTabPoints= new TColStd_HArray1OfReal(1,myDimension*myNbPoints);
123 //
124 //  Table of Points initialization
125 //
126     Standard_Integer ipoint,jp2d,jp3d,index;
127     TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
128     TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));    
129     gp_Pnt2d P2d;
130     gp_Pnt P3d;
131     index=1;
132     
133     for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
134       {
135
136         if(myNbP2d !=0 && myNbP3d ==0 ) {
137                           ToolLine::Value(mySSP,ipoint,TabP2d);
138
139                           for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)
140                               {  P2d = TabP2d.Value(jp2d);
141                                  
142                                  myTabPoints->SetValue(index++,P2d.X());
143                                  myTabPoints->SetValue(index++,P2d.Y());
144                               }
145                          } 
146         if(myNbP3d !=0 && myNbP2d == 0 ) {
147                           ToolLine::Value(mySSP,ipoint,TabP3d);
148                           
149                           for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
150
151                               {  P3d=TabP3d.Value(jp3d);
152                                  
153                                  myTabPoints->SetValue(index++,P3d.X());
154                                  myTabPoints->SetValue(index++,P3d.Y()); 
155                                  myTabPoints->SetValue(index++,P3d.Z());
156
157                               }
158                          }
159         if(myNbP3d !=0 && myNbP2d != 0 ) {
160                           ToolLine::Value(mySSP,ipoint,TabP3d,TabP2d);
161                           
162                           for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
163
164                               {  P3d=TabP3d.Value(jp3d);
165                                  
166                                  myTabPoints->SetValue(index++,P3d.X());
167                                  myTabPoints->SetValue(index++,P3d.Y()); 
168                                  myTabPoints->SetValue(index++,P3d.Z());
169
170                               }
171                           for ( jp2d = 1 ; jp2d <= myNbP2d ; jp2d++)
172
173                               {  P2d=TabP2d.Value(jp2d);
174                                  
175                                  myTabPoints->SetValue(index++,P2d.X());
176                                  myTabPoints->SetValue(index++,P2d.Y()); 
177                               }
178                          }
179       }
180     Init();
181 }
182 //
183 //=======================================================================
184 //function : Init
185 //purpose  : Initializes the tables of constraints
186 //           and verifies if the problem is not over-constrained.
187 //           This method is used in the Create and the method SetConstraint. 
188 //=======================================================================
189 //
190 void AppParCurves_Variational::Init()
191 {
192                               
193     Standard_Integer ipoint,jp2d,jp3d,index,jndex;
194     Standard_Integer CurMultyPoint;
195     TColgp_Array1OfVec TabV3d(1, Max(1,myNbP3d));
196     TColgp_Array1OfVec2d TabV2d(1, Max(1,myNbP2d));
197     TColgp_Array1OfVec TabV3dcurv(1, Max(1,myNbP3d));
198     TColgp_Array1OfVec2d TabV2dcurv(1, Max(1,myNbP2d));
199
200     gp_Vec Vt3d, Vc3d;
201     gp_Vec2d Vt2d, Vc2d;
202
203     myNbConstraints=myConstraints->Length();
204     if (myNbConstraints < 0) Standard_ConstructionError::Raise(); 
205
206     myTypConstraints = new TColStd_HArray1OfInteger(1,Max(1,2*myNbConstraints));
207     myTabConstraints = new TColStd_HArray1OfReal(1,Max(1,2*myDimension*myNbConstraints));
208     myTtheta = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints));
209     myTfthet = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints));
210
211         
212 //
213 // Table of types initialization
214       Standard_Integer iconstr;
215       index=1;
216       jndex=1;
217       CurMultyPoint = 1;
218       myNbPassPoints=0;
219       myNbTangPoints=0;
220       myNbCurvPoints=0;
221       AppParCurves_Constraint valcontr;
222
223       for ( iconstr = myConstraints->Lower() ; iconstr <= myConstraints->Upper() ; iconstr++)
224         {
225         ipoint=(myConstraints->Value(iconstr)).Index();
226
227           valcontr=(myConstraints->Value(iconstr)).Constraint();
228           switch (valcontr) {
229                              case AppParCurves_NoConstraint:
230                                   CurMultyPoint -= myNbP3d * 6 + myNbP2d * 2;
231                                   break ;
232                              case AppParCurves_PassPoint: 
233                                   myTypConstraints->SetValue(index++,ipoint);                                    
234                                   myTypConstraints->SetValue(index++,0);
235                                   myNbPassPoints++;
236                                   if(myNbP2d !=0 ) jndex=jndex+4*myNbP2d;
237                                   if(myNbP3d !=0 ) jndex=jndex+6*myNbP3d;
238                                   break ;
239                              case AppParCurves_TangencyPoint:
240                                   myTypConstraints->SetValue(index++,ipoint);
241                                   myTypConstraints->SetValue(index++,1); 
242                                   myNbTangPoints++;
243                                   if(myNbP2d !=0 && myNbP3d == 0 ) 
244                                     {
245                                      if (ToolLine::Tangency(mySSP,ipoint,TabV2d) == Standard_False)
246                                           Standard_ConstructionError::Raise();
247                                      for (jp2d=1;jp2d<=myNbP2d;jp2d++)
248                                          {  
249                                           Vt2d=TabV2d.Value(jp2d);
250                                           Vt2d.Normalize();
251                                           myTabConstraints->SetValue(jndex++,Vt2d.X());
252                                           myTabConstraints->SetValue(jndex++,Vt2d.Y());
253                                           jndex=jndex+2;
254                                           InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
255                                          }
256                                      } 
257                                    if(myNbP3d !=0 && myNbP2d == 0) 
258                                      {
259                                       if (ToolLine::Tangency(mySSP,ipoint,TabV3d) == Standard_False)
260                                           Standard_ConstructionError::Raise();
261                                       for (jp3d=1;jp3d<=myNbP3d;jp3d++)
262                                           {  
263                                            Vt3d=TabV3d.Value(jp3d);
264                                            Vt3d.Normalize();    
265                                            myTabConstraints->SetValue(jndex++,Vt3d.X());
266                                                  
267                                            myTabConstraints->SetValue(jndex++,Vt3d.Y()); 
268                                                  
269                                            myTabConstraints->SetValue(jndex++,Vt3d.Z());
270                                            jndex=jndex+3;
271                                            InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
272                                           }
273                                       }
274                                    if(myNbP3d !=0 && myNbP2d != 0) 
275                                      {
276                                       if (ToolLine::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False)
277                                           Standard_ConstructionError::Raise();
278                                       for (jp3d=1;jp3d<=myNbP3d;jp3d++)
279                                           {  
280                                            Vt3d=TabV3d.Value(jp3d);
281                                            Vt3d.Normalize();    
282                                            myTabConstraints->SetValue(jndex++,Vt3d.X());
283                                            myTabConstraints->SetValue(jndex++,Vt3d.Y()); 
284                                            myTabConstraints->SetValue(jndex++,Vt3d.Z());
285                                            jndex=jndex+3;
286                                            InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
287                                           }
288
289                                      for (jp2d=1;jp2d<=myNbP2d;jp2d++)
290                                          {  
291                                           Vt2d=TabV2d.Value(jp2d);
292                                           Vt2d.Normalize();
293                                           myTabConstraints->SetValue(jndex++,Vt2d.X());
294                                           myTabConstraints->SetValue(jndex++,Vt2d.Y());
295                                           jndex=jndex+2;
296                                           InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
297                                          }
298                                      }
299
300                                        
301                                        break ;
302
303                              case AppParCurves_CurvaturePoint:
304                                    myTypConstraints->SetValue(index++,ipoint);
305                                    myTypConstraints->SetValue(index++,2);
306                                    myNbCurvPoints++;
307                                    if(myNbP2d !=0 && myNbP3d == 0) 
308                                        {
309                                        if (ToolLine::Tangency(mySSP,ipoint,TabV2d) == Standard_False )
310                                            Standard_ConstructionError::Raise();
311                                        if (ToolLine::Curvature(mySSP,ipoint,TabV2dcurv) == Standard_False)
312                                            Standard_ConstructionError::Raise();
313                                        for (jp2d=1;jp2d<=myNbP2d;jp2d++)
314                                            {  
315                                             Vt2d=TabV2d.Value(jp2d);
316                                             Vt2d.Normalize();
317                                             Vc2d=TabV2dcurv.Value(jp2d);
318                                             if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
319                                                  Standard_ConstructionError::Raise();
320                                             myTabConstraints->SetValue(jndex++,Vt2d.X());
321                                             myTabConstraints->SetValue(jndex++,Vt2d.Y());
322                                             myTabConstraints->SetValue(jndex++,Vc2d.X());
323                                             myTabConstraints->SetValue(jndex++,Vc2d.Y());
324                                             InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
325                                            }
326                                        } 
327
328                                     if(myNbP3d !=0 && myNbP2d == 0 ) 
329                                        {
330                                        if (ToolLine::Tangency(mySSP,ipoint,TabV3d) == Standard_False )
331                                            Standard_ConstructionError::Raise();
332                                        if (ToolLine::Curvature(mySSP,ipoint,TabV3dcurv) == Standard_False)
333                                            Standard_ConstructionError::Raise();
334                                        for (jp3d=1;jp3d<=myNbP3d;jp3d++)
335                                            {  
336                                             Vt3d=TabV3d.Value(jp3d);
337                                             Vt3d.Normalize();
338                                             Vc3d=TabV3dcurv.Value(jp3d);
339                                             if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False) 
340                                                  Standard_ConstructionError::Raise();
341                                             myTabConstraints->SetValue(jndex++,Vt3d.X());
342                                             myTabConstraints->SetValue(jndex++,Vt3d.Y()); 
343                                             myTabConstraints->SetValue(jndex++,Vt3d.Z());
344                                             myTabConstraints->SetValue(jndex++,Vc3d.X());
345                                             myTabConstraints->SetValue(jndex++,Vc3d.Y()); 
346                                             myTabConstraints->SetValue(jndex++,Vc3d.Z());
347                                             InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
348                                            }
349                                        }
350                                     if(myNbP3d !=0 && myNbP2d != 0 ) 
351                                        {
352                                        if (ToolLine::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False )
353                                            Standard_ConstructionError::Raise();
354                                        if (ToolLine::Curvature(mySSP,ipoint,TabV3dcurv,TabV2dcurv) == Standard_False)
355                                            Standard_ConstructionError::Raise();
356                                        for (jp3d=1;jp3d<=myNbP3d;jp3d++)
357                                            {  
358                                             Vt3d=TabV3d.Value(jp3d);
359                                             Vt3d.Normalize();
360                                             Vc3d=TabV3dcurv.Value(jp3d);
361                                             if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False) 
362                                                  Standard_ConstructionError::Raise();
363                                             myTabConstraints->SetValue(jndex++,Vt3d.X());
364                                             myTabConstraints->SetValue(jndex++,Vt3d.Y()); 
365                                             myTabConstraints->SetValue(jndex++,Vt3d.Z());
366                                             myTabConstraints->SetValue(jndex++,Vc3d.X());
367                                             myTabConstraints->SetValue(jndex++,Vc3d.Y()); 
368                                             myTabConstraints->SetValue(jndex++,Vc3d.Z());
369                                             InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
370                                            }
371                                        for (jp2d=1;jp2d<=myNbP2d;jp2d++)
372                                            {  
373                                             Vt2d=TabV2d.Value(jp2d);
374                                             Vt2d.Normalize();
375                                             Vc2d=TabV2dcurv.Value(jp2d);
376                                             if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
377                                                  Standard_ConstructionError::Raise();
378                                             myTabConstraints->SetValue(jndex++,Vt2d.X());
379                                             myTabConstraints->SetValue(jndex++,Vt2d.Y());
380                                             myTabConstraints->SetValue(jndex++,Vc2d.X());
381                                             myTabConstraints->SetValue(jndex++,Vc2d.Y());
382                                             InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
383                                            }
384
385                                       }
386                                       break ; 
387                                  default:
388                                       Standard_ConstructionError::Raise();
389                                 }
390           CurMultyPoint += myNbP3d * 6 + myNbP2d * 2;
391            }
392 // OverConstraint Detection
393         Standard_Integer MaxSeg;
394         if(myWithCutting == Standard_True) MaxSeg = myMaxSegment ;
395         else MaxSeg = 1;
396         if (((myMaxDegree-myNivCont)*MaxSeg-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
397              {
398                myIsOverConstr =Standard_True;
399                myIsCreated = Standard_False;
400              }
401         else {
402           InitSmoothCriterion();
403           myIsCreated = Standard_True;
404         }
405               
406           
407 }
408 //
409 //=======================================================================
410 //function : Approximate
411 //purpose  : Makes the approximation with the current fields.
412 //=======================================================================
413 //
414 void AppParCurves_Variational::Approximate()
415
416 {
417   if (myIsCreated == Standard_False )  StdFail_NotDone:: Raise();
418
419
420   Standard_Real WQuadratic, WQuality;
421
422   TColStd_Array1OfReal Ecarts(myFirstPoint, myLastPoint); 
423
424   mySmoothCriterion->GetWeight(WQuadratic, WQuality);
425   
426   Handle(FEmTool_Curve) TheCurve;  
427
428   mySmoothCriterion->GetCurve(TheCurve);
429
430 //---------------------------------------------------------------------
431
432   TheMotor(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
433
434
435   if(myWithMinMax && myTolerance < myMaxError)
436     Adjusting(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
437     
438 //---------------------------------------------------------------------
439
440   Standard_Integer jp2d,jp3d,index,ipole, 
441                    NbElem = TheCurve->NbElements();
442
443   TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
444   TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
445   Standard_Real debfin[2] = {-1., 1};
446   
447   gp_Pnt2d P2d;
448   gp_Pnt P3d;
449   
450   index=0;
451
452   {
453     Handle(TColStd_HArray2OfReal) PolynomialIntervalsPtr =
454       new TColStd_HArray2OfReal(1,NbElem,1,2) ;
455
456     Handle(TColStd_HArray1OfInteger) NbCoeffPtr = 
457       new TColStd_HArray1OfInteger(1, myMaxSegment);
458
459     Standard_Integer size = myMaxSegment * (myMaxDegree + 1) * myDimension ;
460     Handle(TColStd_HArray1OfReal) CoeffPtr = new TColStd_HArray1OfReal(1,size);
461
462     CoeffPtr->Init(0.);
463
464     Handle(TColStd_HArray1OfReal) IntervallesPtr = 
465       new TColStd_HArray1OfReal(1, NbElem + 1);
466
467     IntervallesPtr->ChangeArray1() = TheCurve->Knots();
468
469     TheCurve->GetPolynom(CoeffPtr->ChangeArray1());
470
471     Standard_Integer ii;
472     
473     for(ii = 1; ii <= NbElem; ii++)  
474       NbCoeffPtr->SetValue(ii, TheCurve->Degree(ii)+1);
475       
476
477     for (ii = PolynomialIntervalsPtr->LowerRow() ;
478          ii <= PolynomialIntervalsPtr->UpperRow() ;ii++) 
479       {
480         PolynomialIntervalsPtr->SetValue(ii,1,debfin[0]) ;
481         PolynomialIntervalsPtr->SetValue(ii,2,debfin[1]) ;
482       }
483 /*   
484         printf("\n =========== Parameters for converting\n");
485         printf("nb_courbes : %d \n", NbElem);
486         printf("tab_option[4] : %d \n", myNivCont);
487         printf("myDimension : %d \n", myDimension);
488         printf("myMaxDegree : %d \n", myMaxDegree);
489         printf("\n NbcoeffPtr :\n");
490         for(ii = 1; ii <= NbElem; ii++)  printf("NbCoeffPtr(%d) = %d \n", ii, NbCoeffPtr->Value(ii));
491         size = NbElem*(myMaxDegree + 1) * myDimension;
492         printf("\n CoeffPtr :\n");
493         for(ii = 1; ii <= size; ii++)  printf("CoeffPtr(%d) = %.8e \n", ii, CoeffPtr->Value(ii));
494         printf("\n PolinimialIntervalsPtr :\n");
495         for (ii = PolynomialIntervalsPtr->LowerRow() ;
496              ii <= PolynomialIntervalsPtr->UpperRow() ;ii++) 
497             {
498              printf(" %d : [%f, %f] \n", ii, PolynomialIntervalsPtr->Value(ii,1),
499                     PolynomialIntervalsPtr->Value(ii,2)) ;
500             }
501         printf("\n IntervallesPtr :\n");
502         for (ii = IntervallesPtr->Lower();
503              ii <= IntervallesPtr->Upper() - 1; ii++) 
504             {
505              printf(" %d : [%f, %f] \n", ii, IntervallesPtr->Value(ii),
506                     IntervallesPtr->Value(ii+1)) ;
507             }
508 */  
509     Convert_CompPolynomialToPoles AConverter(NbElem, 
510                                              myNivCont,
511                                              myDimension,
512                                              myMaxDegree,
513                                              NbCoeffPtr,
514                                              CoeffPtr,
515                                              PolynomialIntervalsPtr,
516                                              IntervallesPtr) ; 
517     if (AConverter.IsDone()) 
518       {
519         Handle(TColStd_HArray2OfReal) PolesPtr ;
520         Handle(TColStd_HArray1OfInteger) Mults;
521         Standard_Integer NbPoles=AConverter.NbPoles();
522 //      Standard_Integer Deg=AConverter.Degree();
523         AppParCurves_Array1OfMultiPoint TabMU(1,NbPoles);
524         AConverter.Poles(PolesPtr) ;
525         AConverter.Knots(myKnots) ;
526         AConverter.Multiplicities(Mults) ;
527         
528         for (ipole=PolesPtr->LowerRow();ipole<=PolesPtr->UpperRow();ipole++)
529           {
530             index=PolesPtr->LowerCol();
531 /*          if(myNbP2d !=0 ) 
532               {
533                 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
534                   {
535                     P2d.SetX(PolesPtr->Value(ipole,index++));
536                     P2d.SetY(PolesPtr->Value(ipole,index++));  
537                     TabP2d.SetValue(jp2d,P2d);
538                   }
539               }*/
540             if(myNbP3d !=0 ) 
541               {
542                 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
543                   {
544                     //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
545                     P3d.SetX(PolesPtr->Value(ipole,index++));
546                     //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
547                     P3d.SetY(PolesPtr->Value(ipole,index++));
548                     //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
549                     P3d.SetZ(PolesPtr->Value(ipole,index++)); 
550                     TabP3d.SetValue(jp3d,P3d);
551                   } 
552               }
553             if(myNbP2d !=0 ) 
554               {
555                 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
556                   {
557                     P2d.SetX(PolesPtr->Value(ipole,index++));
558                     P2d.SetY(PolesPtr->Value(ipole,index++));  
559                     TabP2d.SetValue(jp2d,P2d);
560                   }
561               }
562             if(myNbP2d !=0 && myNbP3d !=0) 
563               {
564                 AppParCurves_MultiPoint aMultiPoint(TabP3d,TabP2d);
565                 TabMU.SetValue(ipole,aMultiPoint);                     
566               }
567             else if (myNbP2d !=0)
568               {
569                 AppParCurves_MultiPoint aMultiPoint(TabP2d);
570                 TabMU.SetValue(ipole,aMultiPoint);                     
571               }
572             else 
573               {
574                 
575                 AppParCurves_MultiPoint aMultiPoint(TabP3d);
576                 TabMU.SetValue(ipole,aMultiPoint);                     
577               }
578             
579           }
580         AppParCurves_MultiBSpCurve aCurve(TabMU,myKnots->Array1(),Mults->Array1());
581         myMBSpCurve=aCurve; 
582         myIsDone = Standard_True;
583       }
584   }
585
586 }
587 //
588 //=======================================================================
589 //function : IsCreated
590 //purpose  : returns True if the creation is done 
591 //=======================================================================
592 //
593 Standard_Boolean AppParCurves_Variational::IsCreated() const 
594 {
595   return myIsCreated;
596 }
597 //
598 //=======================================================================
599 //function : IsDone
600 //purpose  : returns True if the  approximation is ok
601 //=======================================================================
602 //
603 Standard_Boolean AppParCurves_Variational::IsDone() const 
604 {
605   return myIsDone;
606 }
607 //
608 //=======================================================================
609 //function : IsOverConstrained
610 //purpose  : returns True if the problem is overconstrained
611 //           in this case, approximation cannot be done.
612 //=======================================================================
613 //
614 Standard_Boolean AppParCurves_Variational::IsOverConstrained() const 
615 {
616   return myIsOverConstr;
617 }
618 //
619 //=======================================================================
620 //function : Value
621 //purpose  : returns all the BSpline curves approximating the
622 //           MultiLine SSP after minimization of the parameter.
623
624 //=======================================================================
625 //
626 AppParCurves_MultiBSpCurve AppParCurves_Variational::Value() const 
627
628   if (myIsDone == Standard_False)  StdFail_NotDone::Raise(); 
629   return myMBSpCurve;
630
631 }
632 //
633 //=======================================================================
634 //function : MaxError
635 //purpose  : returns the maximum of the distances between 
636 //           the points of the multiline and the approximation 
637 //           curves.
638 //=======================================================================
639 //
640 Standard_Real AppParCurves_Variational::MaxError() const 
641 {
642   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
643   return myMaxError;
644 }
645 //
646 //=======================================================================
647 //function : MaxErrorIndex
648 //purpose  : returns the index of the MultiPoint of ErrorMax
649 //=======================================================================
650 //
651 Standard_Integer AppParCurves_Variational::MaxErrorIndex() const 
652 {
653   if (myIsDone == Standard_False)  StdFail_NotDone::Raise(); 
654   return myMaxErrorIndex;
655 }
656 //
657 //=======================================================================
658 //function : QuadraticError
659 //purpose  :  returns the quadratic average of the distances between 
660 //            the points of the multiline and the approximation 
661 //            curves.
662 //=======================================================================
663 //
664 Standard_Real AppParCurves_Variational::QuadraticError() const 
665 {
666   if (myIsDone == Standard_False)  StdFail_NotDone::Raise(); 
667   return myCriterium[0];
668 }
669 //
670 //=======================================================================
671 //function : Distance
672 //purpose  : returns the distances between the points of the 
673 //           multiline and the approximation curves.
674 //=======================================================================
675 //
676 void AppParCurves_Variational::Distance(math_Matrix& mat)
677
678 {
679       if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
680       Standard_Integer ipoint,jp2d,jp3d,index;
681       TColgp_Array1OfPnt TabP3d(1,Max(1,myNbP3d));
682       TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
683       Standard_Integer j0 = mat.LowerCol() - myFirstPoint;
684
685       gp_Pnt2d P2d;
686       gp_Pnt P3d;
687
688
689       gp_Pnt Pt3d;
690       gp_Pnt2d Pt2d;
691
692       for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
693         {
694           index=1;
695           if(myNbP3d !=0 ) {
696                             ToolLine::Value(mySSP,ipoint,TabP3d);
697
698                             for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
699
700                                       {  P3d=TabP3d.Value(jp3d);
701                                          myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt3d);
702                                          mat(index++, j0 + ipoint)=P3d.Distance(Pt3d);                                 
703
704                                       }
705                             }
706           if(myNbP2d !=0 ) {
707                             if(myNbP3d == 0 ) ToolLine::Value(mySSP,ipoint,TabP2d);
708                             else ToolLine::Value(mySSP,ipoint,TabP3d,TabP2d);
709                             for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)
710
711                                       {  P2d = TabP2d.Value(jp2d);
712                                          myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt2d);
713                                          mat(index++, j0 + ipoint)=P2d.Distance(Pt2d);
714                                       }
715                            } 
716       }
717
718 }
719 //
720 //=======================================================================
721 //function : AverageError
722 //purpose  : returns the average error between          
723 //           the MultiLine and the approximation.
724 //=======================================================================
725 //
726 Standard_Real AppParCurves_Variational::AverageError() const 
727 {
728   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
729   return myAverageError;
730 }
731 //
732 //=======================================================================
733 //function : Parameters
734 //purpose  : returns the parameters uses to the approximations
735 //=======================================================================
736 //
737 const Handle(TColStd_HArray1OfReal)& AppParCurves_Variational::Parameters() const 
738 {
739   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
740   return myParameters;
741 }
742 //
743 //=======================================================================
744 //function : Knots
745 //purpose  : returns the knots uses to the approximations
746 //=======================================================================
747 //
748 const Handle(TColStd_HArray1OfReal)& AppParCurves_Variational::Knots() const 
749 {
750   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
751   return myKnots;
752 }
753 //
754 //=======================================================================
755 //function : Criterium
756 //purpose  : returns the values of the quality criterium.
757 //=======================================================================
758 //
759 void AppParCurves_Variational::Criterium(Standard_Real& VFirstOrder, Standard_Real& VSecondOrder, Standard_Real& VThirdOrder) const 
760 {     
761   if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
762   VFirstOrder=myCriterium[1] ;
763   VSecondOrder=myCriterium[2];
764   VThirdOrder=myCriterium[3];
765 }
766 //
767 //=======================================================================
768 //function : CriteriumWeight
769 //purpose  : returns the Weights (as percent) associed  to the criterium used in
770 //           the  optimization.
771 //=======================================================================
772 //
773 void AppParCurves_Variational::CriteriumWeight(Standard_Real& Percent1, Standard_Real& Percent2, Standard_Real& Percent3) const 
774 {
775   Percent1 = myPercent[0];
776   Percent2 = myPercent[1];
777   Percent3 = myPercent[2] ;
778 }
779 //
780 //=======================================================================
781 //function : MaxDegree
782 //purpose  : returns the Maximum Degree used in the approximation 
783 //=======================================================================
784 //
785 Standard_Integer AppParCurves_Variational::MaxDegree() const 
786 {
787    return myMaxDegree;
788 }
789 //
790 //=======================================================================
791 //function : MaxSegment
792 //purpose  : returns the Maximum of segment used in the approximation
793 //=======================================================================
794 //
795 Standard_Integer AppParCurves_Variational::MaxSegment() const 
796 {
797     return myMaxSegment;
798  }
799 //
800 //=======================================================================
801 //function : Continuity
802 //purpose  : returns the Continuity used in the approximation
803 //=======================================================================
804 //
805 GeomAbs_Shape AppParCurves_Variational::Continuity() const 
806 {
807     return myContinuity;
808 }
809 //
810 //=======================================================================
811 //function : WithMinMax
812 //purpose  : returns if the  approximation  search to  minimize the
813 //           maximum Error or not.
814 //=======================================================================
815 //
816 Standard_Boolean AppParCurves_Variational::WithMinMax() const 
817 {
818     return myWithMinMax;
819 }
820 //
821 //=======================================================================
822 //function : WithCutting
823 //purpose  :  returns if the  approximation can insert new Knots or not.
824 //=======================================================================
825 //
826 Standard_Boolean AppParCurves_Variational::WithCutting() const 
827 {
828    return myWithCutting;
829 }
830 //
831 //=======================================================================
832 //function : Tolerance
833 //purpose  : returns the tolerance used in the approximation.
834 //=======================================================================
835 //
836 Standard_Real AppParCurves_Variational::Tolerance() const 
837 {
838    return myTolerance;
839 }
840 //
841 //=======================================================================
842 //function : NbIterations
843 //purpose  : returns the number of iterations used in the approximation.
844 //=======================================================================
845 //
846 Standard_Integer AppParCurves_Variational::NbIterations() const 
847 {
848    return myNbIterations;
849 }
850 //
851 //=======================================================================
852 //function : Dump
853 //purpose  : Prints on the stream o information on the current state 
854 //           of the object.
855 //=======================================================================
856 //
857 void AppParCurves_Variational::Dump(Standard_OStream& o) const 
858 {
859   o << " \nVariational Smoothing " << endl;
860   o << " Number of multipoints                   "  << myNbPoints << endl;
861   o << " Number of 2d par multipoint "  << myNbP2d << endl;
862   o << " Nombre of 3d par multipoint "  << myNbP3d << endl;
863   o << " Number of PassagePoint      "  << myNbPassPoints << endl;
864   o << " Number of TangencyPoints    "  << myNbTangPoints << endl;
865   o << " Number of CurvaturePoints   "  << myNbCurvPoints << endl;
866   o << " \nTolerance " << o.setf(ios::scientific) << setprecision(3) << setw(9) << myTolerance;
867   if ( WithMinMax()) { o << "  as Max Error." << endl;}
868   else { o << "  as size Error." << endl;}
869   o << "CriteriumWeights : " << myPercent[0] << " , " 
870     << myPercent[1] << " , " << myPercent[2] << endl;
871
872   if (myIsDone ) {
873     o << " MaxError             " << setprecision(3) << setw(9) << myMaxError << endl;
874     o << " Index of  MaxError   " << myMaxErrorIndex << endl;
875     o << " Average Error        " << setprecision(3) << setw(9) << myAverageError << endl;
876     o << " Quadratic Error      " << setprecision(3) << setw(9) << myCriterium[0] << endl;
877     o << " Tension Criterium    " << setprecision(3) << setw(9) << myCriterium[1] << endl;
878     o << " Flexion  Criterium   " << setprecision(3) << setw(9) << myCriterium[2] << endl;
879     o << " Jerk  Criterium      " << setprecision(3) << setw(9) << myCriterium[3] << endl;
880     o << " NbSegments           "  << myKnots->Length()-1 << endl;
881   }
882   else
883    { if (myIsOverConstr) o << "The probleme is overconstraint " << endl;
884      else o << " Erreur dans l''approximation" << endl;
885    }   
886 }
887 //
888 //=======================================================================
889 //function : SetConstraints
890 //purpose  : Define the constraints to approximate
891 //           If this value is incompatible with the others fields
892 //           this method modify nothing and returns false 
893 //=======================================================================
894 //
895 Standard_Boolean AppParCurves_Variational::SetConstraints(const Handle(AppParCurves_HArray1OfConstraintCouple)& aConstraint)
896
897 {
898    myConstraints=aConstraint;
899    Init();
900    if (myIsOverConstr ) return Standard_False;
901    else return Standard_True;
902 }
903 //
904 //=======================================================================
905 //function : SetParameters
906 //purpose  : Defines the parameters used by the approximations.
907 //=======================================================================
908 //
909 void AppParCurves_Variational::SetParameters(const Handle(TColStd_HArray1OfReal)& param)
910 {
911   myParameters->ChangeArray1() = param->Array1();
912 }
913 //
914 //=======================================================================
915 //function : SetKnots
916 //purpose  : Defines the knots used by the approximations
917 //    --          If this value is incompatible with the others fields
918 //    --          this method modify nothing and returns false
919 //=======================================================================
920 //
921 Standard_Boolean AppParCurves_Variational::SetKnots(const Handle(TColStd_HArray1OfReal)& knots)
922 {
923   myKnots->ChangeArray1() = knots->Array1();
924   return Standard_True;
925 }
926 //
927 //=======================================================================
928 //function : SetMaxDegree
929 //purpose  : Define the Maximum Degree used in the approximation
930 //           If this value is incompatible with the others fields
931 //           this method modify nothing and returns false 
932 //=======================================================================
933 //
934 Standard_Boolean AppParCurves_Variational::SetMaxDegree(const Standard_Integer Degree)
935 {
936   if (((Degree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
937        return Standard_False; 
938   else
939      {
940           myMaxDegree=Degree;
941
942           InitSmoothCriterion();
943
944           return Standard_True;
945         }
946
947 }
948 //
949 //=======================================================================
950 //function : SetMaxSegment
951 //purpose  : Define the maximum number of segments used in the approximation
952 //           If this value is incompatible with the others fields
953 //           this method modify nothing and returns false 
954 //=======================================================================
955 //
956 Standard_Boolean AppParCurves_Variational::SetMaxSegment(const Standard_Integer NbSegment)
957 {
958   if ( myWithCutting == Standard_True && 
959       ((myMaxDegree-myNivCont)*NbSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
960        return Standard_False;
961   else
962     {
963        myMaxSegment=NbSegment;
964        return Standard_True;
965     }   
966 }
967 //
968 //=======================================================================
969 //function : SetContinuity
970 //purpose  : Define the Continuity used in the approximation 
971 //           If this value is incompatible with the others fields
972 //           this method modify nothing and returns false 
973 //=======================================================================
974 //
975 Standard_Boolean AppParCurves_Variational::SetContinuity(const GeomAbs_Shape C)
976 {
977     Standard_Integer NivCont=0;
978     switch (C) {
979     case GeomAbs_C0:
980        NivCont=0;
981        break ;
982     case GeomAbs_C1: 
983        NivCont=1;
984        break ;
985     case GeomAbs_C2:  
986        NivCont=2;
987        break ;
988     default:
989     Standard_ConstructionError::Raise();
990   }
991   if (((myMaxDegree-NivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
992        return Standard_False; 
993   else
994     {
995          myContinuity= C;
996          myNivCont=NivCont;
997
998          InitSmoothCriterion();
999        return Standard_True;
1000     }
1001 }
1002 //
1003 //=======================================================================
1004 //function : SetWithMinMax
1005 //purpose  : Define if the  approximation  search to  minimize the
1006 //           maximum Error or not.
1007 //=======================================================================
1008 //
1009 void AppParCurves_Variational::SetWithMinMax(const Standard_Boolean MinMax)
1010 {
1011   myWithMinMax=MinMax;
1012
1013   InitSmoothCriterion();
1014 }
1015 //
1016 //=======================================================================
1017 //function : SetWithCutting
1018 //purpose  : Define if the  approximation can insert new Knots or not.
1019 //           If this value is incompatible with the others fields
1020 //           this method modify nothing and returns false 
1021 //=======================================================================
1022 //
1023 Standard_Boolean AppParCurves_Variational::SetWithCutting(const Standard_Boolean Cutting)
1024 {
1025   if (Cutting == Standard_False) 
1026       {
1027        if (((myMaxDegree-myNivCont)*myKnots->Length()-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1028        return Standard_False;
1029
1030        else
1031           {
1032         myWithCutting=Cutting;
1033         InitSmoothCriterion();
1034         return Standard_True;
1035       }
1036       }
1037   else
1038     {
1039       if (((myMaxDegree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1040        return Standard_False;
1041
1042        else
1043       {
1044         myWithCutting=Cutting;
1045         InitSmoothCriterion();
1046         return Standard_True;
1047       }
1048     }
1049 }
1050 //
1051 //=======================================================================
1052 //function : SetCriteriumWeight
1053 //purpose  : define the Weights (as percent) associed to the criterium used in
1054 //           the  optimization.
1055 //=======================================================================
1056 //
1057 void AppParCurves_Variational::SetCriteriumWeight(const Standard_Real Percent1, const Standard_Real Percent2, const Standard_Real Percent3)
1058 {
1059   if (Percent1 < 0 || Percent2 < 0 || Percent3 < 0 ) Standard_DomainError::Raise();
1060   Standard_Real Total = Percent1 + Percent2 + Percent3;
1061   myPercent[0] = Percent1/Total;
1062   myPercent[1] = Percent2/Total;
1063   myPercent[2] = Percent3/Total;
1064
1065   InitSmoothCriterion();
1066 }
1067 //
1068 //=======================================================================
1069 //function : SetCriteriumWeight
1070 //purpose  :  define the  Weight   (as  percent)  associed  to   the
1071 //            criterium   Order used in   the optimization  : Others
1072 //            weights are updated.
1073 //=======================================================================
1074 //
1075 void AppParCurves_Variational::SetCriteriumWeight(const Standard_Integer Order, const Standard_Real Percent)
1076 {
1077   if ( Percent < 0 ) Standard_DomainError::Raise();
1078   if ( Order < 1 || Order > 3 ) Standard_ConstructionError::Raise();
1079   myPercent[Order-1] = Percent;
1080   Standard_Real Total = myPercent[0] + myPercent[1] + myPercent[2];
1081   myPercent[0] = myPercent[0] / Total;
1082   myPercent[1] = myPercent[1] / Total;
1083   myPercent[2] = myPercent[2] / Total;
1084
1085   InitSmoothCriterion();
1086   
1087 }
1088 //
1089 //=======================================================================
1090 //function : SetTolerance
1091 //purpose  : define the tolerance used in the approximation.
1092 //=======================================================================
1093 //
1094 void AppParCurves_Variational::SetTolerance(const Standard_Real Tol)
1095 {
1096   myTolerance=Tol;
1097   InitSmoothCriterion();
1098 }
1099 //
1100 //=======================================================================
1101 //function : SetNbIterations
1102 //purpose  : define the number of iterations used in the approximation.
1103 //=======================================================================
1104 //
1105 void AppParCurves_Variational::SetNbIterations(const Standard_Integer Iter)
1106 {
1107   myNbIterations=Iter;
1108 }
1109
1110
1111 // Private methods
1112 #include <AppParCurves_Variational_1.gxx> // TheMotor
1113 #include <AppParCurves_Variational_2.gxx> // Optimization
1114 #include <AppParCurves_Variational_3.gxx> // Project
1115 #include <AppParCurves_Variational_4.gxx> // ACR
1116 #include <AppParCurves_Variational_5.gxx> // SplitCurve
1117 #include <AppParCurves_Variational_6.gxx> // InitSmoothCriterion and other Init... methods
1118 #include <AppParCurves_Variational_7.gxx> // Adjusting
1119 #include <AppParCurves_Variational_8.gxx> // AssemblingConstraints
1120 #include <AppParCurves_Variational_9.gxx> // InitTtheta and InitTfthet methods