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