src/AppParCurves/AppParCurves_Function.gxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 // Lpa, le 20/09/91
  16
  17
  18 // Calcul de la valeur de F et grad_F, connaissant le parametrage.
  19 // Cette fonction, appelee par le gradient conjugue, calcul F et
  20 // DF(ui, Poles(ui)) ce qui implique un calcul des nouveaux poles
  21 //  a chaque appel.
  22
  23 #define No_Standard_RangeError
  24 #define No_Standard_OutOfRange
  25
  26
  27
  28 #include <AppParCurves_MultiCurve.hxx>
  29 #include <AppParCurves_MultiPoint.hxx>
  30 #include <TColStd_HArray1OfInteger.hxx>
  31 #include <gp_Pnt.hxx>
  32 #include <gp_Pnt2d.hxx>
  33 #include <gp_Vec.hxx>
  34 #include <gp_Vec2d.hxx>
  35 #include <TColgp_Array1OfPnt.hxx>
  36 #include <TColgp_Array1OfPnt2d.hxx>
  37 #include <AppParCurves_ConstraintCouple.hxx>
  38
  39 AppParCurves_Function::
  40   AppParCurves_Function(const MultiLine& SSP,
  41          const Standard_Integer FirstPoint,
  42          const Standard_Integer LastPoint,
  43          const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
  44          const math_Vector& Parameters,
  45          const Standard_Integer Deg) :
  46          MyMultiLine(SSP),
  47          MyMultiCurve(Deg+1),
  48          myParameters(Parameters.Lower(), Parameters.Upper()),
  49          ValGrad_F(FirstPoint, LastPoint),
  50          MyF(FirstPoint, LastPoint,
  51              1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
  52          PTLX(FirstPoint, LastPoint,
  53              1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
  54          PTLY(FirstPoint, LastPoint,
  55              1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
  56          PTLZ(FirstPoint, LastPoint,
  57              1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
  58          A(FirstPoint, LastPoint, 1, Deg+1),
  59          DA(FirstPoint, LastPoint, 1, Deg+1),
  60          MyLeastSquare(SSP, FirstPoint, LastPoint,
  61                        FirstConstraint(TheConstraints, FirstPoint),
  62                        LastConstraint(TheConstraints, LastPoint), Deg+1)
  63 {
  64   Standard_Integer i;
  65   for (i=Parameters.Lower(); i<=Parameters.Upper();i++)
  66     myParameters(i)=Parameters(i);
  67   FirstP = FirstPoint;
  68   LastP = LastPoint;
  69   myConstraints = TheConstraints;
  70   NbP = LastP-FirstP+1;
  71   Adeb = FirstP;
  72   Afin = LastP;
  73   Degre = Deg;
  74   Contraintes = Standard_False;
  75   Standard_Integer myindex;
  76   Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
  77   AppParCurves_ConstraintCouple mycouple;
  78   AppParCurves_Constraint Cons;
  79
  80   for (i = low; i <= upp; i++) {
  81     mycouple = TheConstraints->Value(i);
  82     Cons = mycouple.Constraint();
  83     myindex = mycouple.Index();
  84     if (myindex == FirstP) {
  85       if (Cons >= 1) Adeb = Adeb+1;
  86     }
  87     else if (myindex == LastP) {
  88       if (Cons >= 1) Afin = Afin-1;
  89     }
  90     else {
  91       if (Cons >= 1) Contraintes = Standard_True;
  92     }
  93   }
  94
  95   Standard_Integer nb3d = ToolLine::NbP3d(SSP);
  96   Standard_Integer nb2d = ToolLine::NbP2d(SSP);
  97   Standard_Integer mynb3d= nb3d, mynb2d=nb2d;
  98   if (nb3d == 0) mynb3d = 1;
  99   if (nb2d == 0) mynb2d = 1;
 100
 101   NbCu = nb3d+nb2d;
 102   tabdim = new TColStd_HArray1OfInteger(0, NbCu-1);
 103
 104   if (Contraintes) {
 105     for (i = 1; i <= NbCu; i++) {
 106       if (i <= nb3d) tabdim->SetValue(i-1, 3);
 107       else tabdim->SetValue(i-1, 2);
 108     }
 109
 110     TColgp_Array1OfPnt TabP(1, mynb3d);
 111     TColgp_Array1OfPnt2d TabP2d(1, mynb2d);
 112
 113     for ( i = FirstP; i <= LastP; i++) {
 114       if (nb3d != 0 && nb2d != 0) ToolLine::Value(SSP, i, TabP, TabP2d);
 115       else if (nb3d != 0)         ToolLine::Value(SSP, i, TabP);
 116       else                        ToolLine::Value(SSP, i, TabP2d);
 117       for (Standard_Integer j = 1; j <= NbCu; j++) {
 118         if (tabdim->Value(j-1) == 3) {
 119           TabP(j).Coord(PTLX(i, j), PTLY(i, j),PTLZ(i, j));
 120         }
 121         else {
 122           TabP2d(j).Coord(PTLX(i, j), PTLY(i, j));
 123         }
 124       }
 125     }
 126   }
 127 }
 128
 129
 130 AppParCurves_Constraint AppParCurves_Function::FirstConstraint
 131   (const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
 132    const Standard_Integer FirstPoint) const
 133 {
 134   Standard_Integer i, myindex;
 135   Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
 136   AppParCurves_ConstraintCouple mycouple;
 137   AppParCurves_Constraint Cons = AppParCurves_NoConstraint;
 138
 139   for (i = low; i <= upp; i++) {
 140     mycouple = TheConstraints->Value(i);
 141     Cons = mycouple.Constraint();
 142     myindex = mycouple.Index();
 143     if (myindex == FirstPoint) {
 144       break;
 145     }
 146   }
 147   return Cons;
 148 }
 149
 150
 151 AppParCurves_Constraint AppParCurves_Function::LastConstraint
 152   (const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
 153    const Standard_Integer LastPoint) const
 154 {
 155   Standard_Integer i, myindex;
 156   Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
 157   AppParCurves_ConstraintCouple mycouple;
 158   AppParCurves_Constraint Cons = AppParCurves_NoConstraint;
 159
 160   for (i = low; i <= upp; i++) {
 161     mycouple = TheConstraints->Value(i);
 162     Cons = mycouple.Constraint();
 163     myindex = mycouple.Index();
 164     if (myindex == LastPoint) {
 165       break;
 166     }
 167   }
 168   return Cons;
 169 }
 170
 171
 172
 173
 174 Standard_Boolean AppParCurves_Function::Value (const math_Vector& X,
 175                                                Standard_Real& F) {
 176
 177   myParameters = X;
 178
 179   // Resolution moindres carres:
 180   // ===========================
 181   MyLeastSquare.Perform(myParameters);
 182   if (!(MyLeastSquare.IsDone())) {
 183     Done = Standard_False;
 184     return Standard_False;
 185   }
 186   if (!Contraintes) {
 187     MyLeastSquare.Error(FVal, ERR3d, ERR2d);
 188     F = FVal;
 189   }
 190
 191   // Resolution avec contraintes:
 192   // ============================
 193   else {
 194     Standard_Integer Npol = Degre+1;
 195 //    Standard_Boolean Ext = Standard_True;
 196     Standard_Integer Ci, i, j, dimen;
 197     Standard_Real AA, BB, CC, AIJ, FX, FY, FZ, Fi;
 198     math_Vector PTCXCI(1, Npol), PTCYCI(1, Npol), PTCZCI(1, Npol);
 199     ERR3d = ERR2d = 0.0;
 200
 201     MyMultiCurve = MyLeastSquare.BezierValue();
 202
 203     A = MyLeastSquare.FunctionMatrix();
 204     ResolCons Resol(MyMultiLine, MyMultiCurve, FirstP, LastP, myConstraints,
 205                     A, MyLeastSquare.DerivativeFunctionMatrix());
 206     if (!Resol.IsDone()) {
 207       Done = Standard_False;
 208       return Standard_False;
 209     }
 210
 211     // Calcul de F = Sum||C(ui)-Ptli||2  sur toutes les courbes :
 212     // ========================================================================
 213     FVal = 0.0;
 214
 215     for (Ci = 1; Ci <= NbCu; Ci++) {
 216       dimen = tabdim->Value(Ci-1);
 217       for (j = 1; j <= Npol; j++) {
 218         if (dimen == 3){
 219           MyMultiCurve.Value(j).Point(Ci).Coord(PTCXCI(j),PTCYCI(j),PTCZCI(j));
 220         }
 221         else{
 222           MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCXCI(j), PTCYCI(j));
 223         }
 224       }
 225
 226       // Calcul de F:
 227       // ============
 228       for (i = Adeb; i <= Afin; i++) {
 229         AA = 0.0; BB = 0.0; CC = 0.0;
 230         for (j = 1; j <= Npol; j++) {
 231           AIJ = A(i, j);
 232           AA += AIJ*PTCXCI(j);
 233           BB += AIJ*PTCYCI(j);
 234           if (dimen == 3) {
 235             CC += AIJ*PTCZCI(j);
 236           }
 237         }
 238         FX = AA-PTLX(i, Ci);
 239         FY = BB-PTLY(i, Ci);
 240         MyF(i,Ci) = FX*FX + FY*FY;
 241         if (dimen == 3) {
 242           FZ = CC-PTLZ(i,Ci);
 243           MyF(i, Ci) += FZ*FZ;
 244           Fi = MyF(i, Ci);
 245           if (Sqrt(Fi) > ERR3d) ERR3d = Sqrt(Fi);
 246         }
 247         else {
 248           Fi = MyF(i, Ci);
 249           if (Sqrt(Fi) > ERR2d) ERR2d = Sqrt(Fi);
 250         }
 251         FVal += Fi;
 252       }
 253     }
 254     F = FVal;
 255   }
 256   return Standard_True;
 257 }
 258
 259
 260
 261
 262 void AppParCurves_Function::Perform(const math_Vector& X) {
 263   Standard_Integer j;
 264
 265   myParameters = X;
 266   // Resolution moindres carres:
 267   // ===========================
 268   MyLeastSquare.Perform(myParameters);
 269
 270   if (!(MyLeastSquare.IsDone())) {
 271     Done = Standard_False;
 272     return;
 273   }
 274
 275   for(j = myParameters.Lower(); j <= myParameters.Upper(); j++) {
 276     ValGrad_F(j) = 0.0;
 277   }
 278
 279   if (!Contraintes) {
 280     MyLeastSquare.ErrorGradient(ValGrad_F, FVal, ERR3d, ERR2d);
 281   }
 282   else {
 283     Standard_Integer Pi, Ci, i, k, dimen;
 284     Standard_Integer  Npol = Degre+1;
 285     Standard_Real Scal, AA, BB, CC, DAA, DBB, DCC;
 286     Standard_Real FX, FY, FZ, AIJ, DAIJ, px, py, pz, Fi;
 287     AppParCurves_Constraint Cons=AppParCurves_NoConstraint;
 288     math_Matrix Grad_F(FirstP, LastP, 1, NbCu, 0.0);
 289     math_Vector PTCXCI(1, Npol), PTCYCI(1, Npol), PTCZCI(1, Npol);
 290     math_Vector PTCOXCI(1, Npol), PTCOYCI(1, Npol), PTCOZCI(1, Npol);
 291 //    Standard_Boolean Ext = Standard_True;
 292     ERR3d = ERR2d = 0.0;
 293
 294     math_Matrix PTCOX(1, Npol, 1, NbCu), PTCOY(1, Npol, 1, NbCu),
 295                 PTCOZ(1, Npol,1, NbCu);
 296     math_Matrix PTCX(1, Npol, 1, NbCu), PTCY(1, Npol, 1, NbCu),
 297                 PTCZ(1, Npol,1, NbCu);
 298     Standard_Integer Inc;
 299
 300     MyMultiCurve = MyLeastSquare.BezierValue();
 301
 302     for (Ci =1; Ci <= NbCu; Ci++) {
 303       dimen = tabdim->Value(Ci-1);
 304       for (j = 1; j <= Npol; j++) {
 305         if (dimen == 3){
 306           MyMultiCurve.Value(j).Point(Ci).Coord(PTCOX(j, Ci),
 307                                                 PTCOY(j, Ci),
 308                                                 PTCOZ(j, Ci));
 309         }
 310         else{
 311           MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCOX(j, Ci), PTCOY(j, Ci));
 312           PTCOZ(j, Ci) = 0.0;
 313         }
 314       }
 315     }
 316
 317     A = MyLeastSquare.FunctionMatrix();
 318     DA = MyLeastSquare.DerivativeFunctionMatrix();
 319
 320     // Resolution avec contraintes:
 321     // ============================
 322
 323     ResolCons Resol(MyMultiLine, MyMultiCurve, FirstP, LastP,
 324                     myConstraints, A, DA);
 325     if (!Resol.IsDone()) {
 326       Done = Standard_False;
 327       return;
 328     }
 329
 330
 331     // Calcul de F = Sum||C(ui)-Ptli||2 et du gradient non contraint de F pour
 332     // chaque point PointIndex.
 333     // ========================================================================
 334     FVal = 0.0;
 335     for(j = FirstP; j <= LastP; j++) {
 336       ValGrad_F(j) = 0.0;
 337     }
 338
 339     math_Matrix TrA(A.LowerCol(), A.UpperCol(), A.LowerRow(), A.UpperRow());
 340     math_Matrix TrDA(DA.LowerCol(), DA.UpperCol(), DA.LowerRow(), DA.UpperRow());
 341     math_Matrix RESTM(A.LowerCol(), A.UpperCol(), A.LowerCol(), A.UpperCol());
 342
 343     const math_Matrix& K = Resol.ConstraintMatrix();
 344     const math_Matrix& DK = Resol.ConstraintDerivative(MyMultiLine, X, Degre, DA);
 345     math_Matrix TK(K.LowerCol(), K.UpperCol(), K.LowerRow(), K.UpperRow());
 346     TK = K.Transposed();
 347     const math_Vector& Vardua = Resol.Duale();
 348     math_Matrix KK(K.LowerCol(), K.UpperCol(), Vardua.Lower(), Vardua.Upper());
 349     KK = (K.Transposed())*(Resol.InverseMatrix());
 350     math_Matrix DTK(DK.LowerCol(), DK.UpperCol(), DK.LowerRow(), DK.UpperRow());
 351     DTK = DK.Transposed();
 352     TrA = A.Transposed();
 353     TrDA = DA.Transposed();
 354     RESTM = ((A.Transposed()*A).Inverse());
 355
 356     math_Vector DPTCO(1, K.ColNumber());
 357     math_Matrix DPTCO1(FirstP, LastP, 1, K.ColNumber());
 358     math_Vector DKPTC(1, K.RowNumber());
 359
 360
 361
 362
 363     FVal = 0.0;
 364     for (Ci = 1; Ci <= NbCu; Ci++) {
 365       dimen = tabdim->Value(Ci-1);
 366       for (j = 1; j <= Npol; j++) {
 367         if (dimen == 3){
 368           MyMultiCurve.Value(j).Point(Ci).Coord(PTCX(j, Ci),
 369                                                 PTCY(j, Ci),
 370                                                 PTCZ(j, Ci));
 371         }
 372         else{
 373           MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCX(j, Ci), PTCY(j,Ci));
 374           PTCZ(j, Ci) = 0.0;
 375         }
 376       }
 377     }
 378
 379
 380     // Calcul du gradient sans contraintes:
 381     // ====================================
 382
 383     for (Ci = 1; Ci <= NbCu; Ci++) {
 384       dimen = tabdim->Value(Ci-1);
 385       for (i = Adeb; i <= Afin; i++) {
 386         AA = 0.0; BB = 0.0; CC = 0.0; DAA = 0.0; DBB = 0.0; DCC = 0.0;
 387         for (j = 1; j <= Npol; j++) {
 388           AIJ = A(i, j); DAIJ = DA(i, j);
 389           px = PTCX(j, Ci); py = PTCY(j, Ci);
 390           AA += AIJ*px;  BB += AIJ*py;
 391           DAA += DAIJ*px;  DBB += DAIJ*py;
 392           if (dimen == 3) {
 393             pz = PTCZ(j, Ci);
 394             CC += AIJ*pz;  DCC += DAIJ*pz;
 395           }
 396         }
 397         FX = AA-PTLX(i, Ci);
 398         FY = BB-PTLY(i, Ci);
 399         MyF(i,Ci) = FX*FX + FY*FY;
 400         Grad_F(i, Ci) = 2.0*(DAA*FX + DBB*FY);
 401         if (dimen == 3) {
 402           FZ = CC-PTLZ(i,Ci);
 403           MyF(i, Ci) += FZ*FZ;
 404           Grad_F(i, Ci) += 2.0*DCC*FZ;
 405           Fi = MyF(i, Ci);
 406           if (Sqrt(Fi) > ERR3d) ERR3d = Sqrt(Fi);
 407         }
 408         else {
 409           Fi = MyF(i, Ci);
 410           if (Sqrt(Fi) > ERR2d) ERR2d = Sqrt(Fi);
 411         }
 412         FVal += Fi;
 413         ValGrad_F(i) += Grad_F(i, Ci);
 414       }
 415     }
 416
 417
 418     // Calcul de DK*PTC:
 419     // =================
 420     for (i = 1; i <= K.RowNumber(); i++) {
 421       Inc = 0;
 422       for (Ci = 1; Ci <= NbCu; Ci++) {
 423         dimen = tabdim->Value(Ci-1);
 424         DKPTC(i) = 0.0;
 425         for (j = 1; j <= Npol; j++) {
 426           DKPTC(i) += DK(i, j+Inc)*PTCX(j, Ci)+ DK(i, j+Inc+Npol)*PTCY(j, Ci);
 427           if (dimen == 3) {
 428             DKPTC(i) += DK(i, j+Inc+2*Npol)*PTCZ(j, Ci);
 429           }
 430         }
 431         if (dimen == 3) Inc = Inc +3*Npol;
 432         else Inc = Inc +2*Npol;
 433       }
 434     }
 435
 436     math_Vector DERR(DTK.LowerRow(), DTK.UpperRow());
 437     DERR = (DTK)*Vardua-KK* ((DKPTC) + K*(DTK)*Vardua);
 438
 439     // rajout du gradient avec contraintes:
 440     // ====================================
 441     // dPTCO1/duk = [d(TA)/duk*[A*PTCO-PTL] + TA*dA/duk*PTCO]
 442
 443
 444     Inc = 0;
 445
 446     math_Vector Errx(A.LowerRow(), A.UpperRow());
 447     math_Vector Erry(A.LowerRow(), A.UpperRow());
 448     math_Vector Errz(A.LowerRow(), A.UpperRow());
 449     math_Vector Scalx(DA.LowerRow(), DA.UpperRow());
 450     math_Vector Scaly(DA.LowerRow(), DA.UpperRow());
 451     math_Vector Scalz(DA.LowerRow(), DA.UpperRow());
 452     math_Vector Erruzax(PTCXCI.Lower(), PTCXCI.Upper());
 453     math_Vector Erruzay(PTCYCI.Lower(), PTCYCI.Upper());
 454     math_Vector Erruzaz(PTCZCI.Lower(), PTCZCI.Upper());
 455     math_Vector TrDAPI(TrDA.LowerRow(), TrDA.UpperRow());
 456     math_Vector TrAPI(TrA.LowerRow(), TrA.UpperRow());
 457
 458     for (Ci = 1; Ci <= NbCu; Ci++) {
 459       dimen = tabdim->Value(Ci-1);
 460       PTCOXCI = PTCOX.Col(Ci);
 461       PTCOYCI = PTCOY.Col(Ci);
 462       PTCOZCI = PTCOZ.Col(Ci);
 463       PTCXCI = PTCX.Col(Ci);
 464       PTCYCI = PTCY.Col(Ci);
 465       PTCZCI = PTCZ.Col(Ci);
 466
 467
 468       Errx = (A*PTCOXCI - PTLX.Col(Ci));
 469       Erry = (A*PTCOYCI - PTLY.Col(Ci));
 470       Errz = (A*PTCOZCI - PTLZ.Col(Ci));
 471       Scalx = (DA*PTCOXCI);   // Scal = DA * PTCO
 472       Scaly = (DA*PTCOYCI);
 473       Scalz = (DA*PTCOZCI);
 474       Erruzax = (PTCXCI - PTCOXCI);
 475       Erruzay = (PTCYCI - PTCOYCI);
 476       Erruzaz = (PTCZCI - PTCOZCI);
 477
 478       for (Pi = FirstP; Pi <= LastP; Pi++) {
 479         TrDAPI = (TrDA.Col(Pi));
 480         TrAPI = (TrA.Col(Pi));
 481         Standard_Real Taa = TrAPI*A.Row(Pi);
 482         Scal = 0.0;
 483         for (j = 1; j <= Npol; j++) {
 484           DPTCO1(Pi, j + Inc) = (TrDAPI*Errx(Pi)+TrAPI*Scalx(Pi))(j);
 485           DPTCO1(Pi, j + Inc+ Npol) = (TrDAPI*Erry(Pi)+TrAPI*Scaly(Pi))(j);
 486           Scal += DPTCO1(Pi, j+Inc)* Taa*Erruzax(j) + DPTCO1(Pi, j+Inc+Npol)*Taa*Erruzay(j);
 487           if (dimen == 3) {
 488             DPTCO1(Pi, j + Inc+ 2*Npol) = (TrDAPI*Errz(Pi)+TrAPI*Scalz(Pi))(j);
 489             Scal += DPTCO1(Pi, j+Inc+2*Npol)*Taa*Erruzaz(j);
 490           }
 491         }
 492         ValGrad_F(Pi) = ValGrad_F(Pi) - 2*Scal;
 493       }
 494       if (dimen == 3) Inc = Inc + 3*Npol;
 495       else Inc = Inc +2*Npol;
 496     }
 497
 498
 499     // on calcule DPTCO = - RESTM * DPTCO1:
 500
 501     // Calcul de DPTCO/duk:
 502     // dPTCO/duk = -Inv(T(A)*A)*[d(TA)/duk*[A*PTCO-PTL] + TA*dA/duk*PTCO]
 503
 504     Standard_Integer low=myConstraints->Lower(), upp=myConstraints->Upper();
 505     Inc = 0;
 506     for (Pi = FirstP; Pi <= LastP; Pi++) {
 507       for (i = low; i <= upp; i++) {
 508         if (myConstraints->Value(i).Index() == Pi) {
 509           Cons = myConstraints->Value(i).Constraint();
 510           break;
 511         }
 512       }
 513       if (Cons >= 1) {
 514         Inc = 0;
 515         for (Ci = 1; Ci <= NbCu; Ci++) {
 516           dimen = tabdim->Value(Ci-1);
 517           for (j = 1; j <= Npol; j++) {
 518             DPTCO(j+Inc) = 0.0;
 519             DPTCO(j+Inc+Npol) = 0.0;
 520             if (dimen == 3) DPTCO(j+Inc+2*Npol) = 0.0;
 521             for (k = 1; k <= Npol; k++) {
 522               DPTCO(j+Inc) = DPTCO(j+Inc) -RESTM(j, k) * DPTCO1(Pi, j+Inc);
 523               DPTCO(j+Inc+Npol)=DPTCO(j+Inc+Npol)-RESTM(j, k)*DPTCO1(Pi,j+Inc+Npol);
 524               if (dimen == 3) {
 525                 DPTCO(j+Inc+2*Npol) = DPTCO(j+Inc+2*Npol)
 526                   -RESTM(j, k) * DPTCO1(Pi, j+Inc+2*Npol);
 527               }
 528             }
 529           }
 530           if (dimen == 3) Inc += 3*Npol;
 531           else Inc += 2*Npol;
 532         }
 533
 534         DERR = DERR-KK*K*DPTCO;
 535
 536         Inc = 0;
 537         for (Ci = 1; Ci <= NbCu; Ci++) {
 538           dimen = tabdim->Value(Ci-1);
 539           PTCOXCI = PTCOX.Col(Ci);
 540           PTCOYCI = PTCOY.Col(Ci);
 541           PTCOZCI = PTCOZ.Col(Ci);
 542           PTCXCI = PTCX.Col(Ci);
 543           PTCYCI = PTCY.Col(Ci);
 544           PTCZCI = PTCZ.Col(Ci);
 545           Erruzax = (PTCXCI - PTCOXCI);
 546           Erruzay = (PTCYCI - PTCOYCI);
 547           Erruzaz = (PTCZCI - PTCOZCI);
 548           Scal = 0.0;
 549
 550           for (j = 1; j <= Npol ; j++) {
 551             Scal = (A(Pi, j)*Erruzax(j)) * (A(Pi, j)*DERR(j+Inc)) +
 552               (A(Pi, j)*Erruzay(j)) * (A(Pi, j)*DERR(j+Inc+Npol));
 553             if (dimen == 3) {
 554               Scal += (A(Pi, j)*Erruzax(j)) * (A(Pi, j)*DERR(j+Inc+2*Npol));
 555             }
 556           }
 557
 558           ValGrad_F(Pi) = ValGrad_F(Pi) + 2*Scal;
 559           if (dimen == 3) Inc = Inc +3*Npol;
 560           else Inc = Inc + 2*Npol;
 561         }
 562       }
 563     }
 564
 565   }
 566 }
 567
 568
 569 Standard_Integer AppParCurves_Function::NbVariables() const{
 570   return NbP;
 571 }
 572
 573
 574 Standard_Boolean AppParCurves_Function::Gradient (const math_Vector& X,
 575                                                   math_Vector& G) {
 576
 577   Perform(X);
 578   G = ValGrad_F;
 579
 580   return Standard_True;
 581 }
 582
 583
 584 Standard_Boolean AppParCurves_Function::Values (const math_Vector& X,
 585                                                 Standard_Real& F,
 586                                                 math_Vector& G) {
 587
 588
 589   Perform(X);
 590   F = FVal;
 591   G = ValGrad_F;
 592   return Standard_True;
 593 }
 594
 595
 596 const AppParCurves_MultiCurve& AppParCurves_Function::CurveValue() {
 597   if (!Contraintes)  MyMultiCurve = MyLeastSquare.BezierValue();
 598   return MyMultiCurve;
 599 }
 600
 601
 602 Standard_Real AppParCurves_Function::Error(const Standard_Integer IPoint,
 603                                      const Standard_Integer CurveIndex) const {
 604   return Sqrt(MyF(IPoint, CurveIndex));
 605 }
 606
 607 Standard_Real AppParCurves_Function::MaxError3d() const
 608 {
 609   return ERR3d;
 610 }
 611
 612 Standard_Real AppParCurves_Function::MaxError2d() const
 613 {
 614   return ERR2d;
 615 }
 616
 617
 618
 619 const math_Vector& AppParCurves_Function::NewParameters() const
 620 {
 621   return myParameters;
 622 }