src/AppParCurves/AppParCurves_Variational.gxx

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