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