src/AppCont/AppCont_LeastSquare.cxx

   1 // Created on: 1995-03-14
   2 // Created by: Modelistation
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #ifndef OCCT_DEBUG
  18 #define No_Standard_OutOfRange
  19 #define No_Standard_RangeError
  20 #endif
  21 #include <AppCont_LeastSquare.hxx>
  22
  23 #include <math.hxx>
  24 #include <AppParCurves_MultiPoint.hxx>
  25 #include <AppCont_ContMatrices.hxx>
  26 #include <PLib.hxx>
  27
  28
  29 //=======================================================================
  30 //function : AppCont_LeastSquare
  31 //purpose  :
  32 //=======================================================================
  33 void AppCont_LeastSquare::FixSingleBorderPoint(const AppCont_Function&       theSSP,
  34                                                const Standard_Real           theU,
  35                                                const Standard_Real           theU0,
  36                                                const Standard_Real           theU1,
  37                                                NCollection_Array1<gp_Pnt2d>& theFix2d,
  38                                                NCollection_Array1<gp_Pnt>&   theFix)
  39 {
  40   Standard_Real aMaxIter = 15.0;
  41   Standard_Integer j, i2;
  42   NCollection_Array1<gp_Pnt>   aTabP(1, Max (myNbP, 1)),     aPrevP(1, Max (myNbP, 1));
  43   NCollection_Array1<gp_Pnt2d> aTabP2d(1,  Max (myNbP2d, 1)), aPrevP2d(1,  Max (myNbP2d, 1));
  44   Standard_Real aMult = ((theU - theU0) > (theU1 - theU)) ? 1.0: -1.0;
  45   Standard_Real aStartParam = (theU0 + theU1) / 2.0,
  46                 aCurrParam, aPrevDist = 1.0, aCurrDist = 1.0;
  47
  48   for (Standard_Real anIter = 1.0; anIter < aMaxIter; anIter += 1.0)
  49   {
  50     aCurrParam = aStartParam + aMult * (1 - pow(10, -anIter)) * (theU1 - theU0) / 2.0;
  51     theSSP.Value(aCurrParam, aTabP2d, aTabP);
  52
  53     // from second iteration
  54     if (anIter > 1.5)
  55     {
  56       aCurrDist = 0.0;
  57
  58       i2 = 1;
  59       for (j = 1; j <= myNbP; j++)
  60       {
  61         aCurrDist += aTabP(j).Distance(aPrevP(j));
  62         i2 += 3;
  63       }
  64       for (j = 1; j <= myNbP2d; j++)
  65       {
  66         aCurrDist += aTabP2d(j).Distance(aPrevP2d(j));
  67         i2 += 2;
  68       }
  69
  70       // from the third iteration
  71       if (anIter > 2.5 && aCurrDist / aPrevDist > 10.0)
  72         break;
  73     }
  74     aPrevP = aTabP;
  75     aPrevP2d = aTabP2d;
  76     aPrevDist = aCurrDist;
  77   }
  78   theFix2d = aPrevP2d;
  79   theFix   = aPrevP;
  80 }
  81
  82
  83 //=======================================================================
  84 //function : AppCont_LeastSquare
  85 //purpose  :
  86 //=======================================================================
  87
  88 AppCont_LeastSquare::AppCont_LeastSquare(const AppCont_Function&       SSP,
  89                                          const Standard_Real           U0,
  90                                          const Standard_Real           U1,
  91                                          const AppParCurves_Constraint FirstCons,
  92                                          const AppParCurves_Constraint LastCons,
  93                                          const Standard_Integer        Deg,
  94                                          const Standard_Integer        myNbPoints)
  95 : mySCU(Deg+1),
  96   myPoints(1, myNbPoints, 1, 3 * SSP.GetNbOf3dPoints() + 2 * SSP.GetNbOf2dPoints()),
  97   myPoles(1, Deg + 1, 1, 3 * SSP.GetNbOf3dPoints() + 2 * SSP.GetNbOf2dPoints(), 0.0),
  98   myParam(1, myNbPoints),
  99   myVB(1, Deg+1, 1, myNbPoints),
 100   myPerInfo(1, 3 * SSP.GetNbOf3dPoints() + 2 * SSP.GetNbOf2dPoints() )
 101 {
 102   myDone = Standard_False;
 103   myDegre = Deg;
 104   math_Matrix InvM(1, Deg+1, 1, Deg + 1);
 105   Standard_Integer i, j, k, c, i2;
 106   Standard_Integer classe = Deg + 1, cl1 = Deg;
 107   Standard_Real U, dU, Coeff, Coeff2;
 108   Standard_Real IBij, IBPij;
 109
 110   Standard_Integer FirstP = 1, LastP = myNbPoints;
 111   Standard_Integer nbcol = 3 * SSP.GetNbOf3dPoints() + 2 * SSP.GetNbOf2dPoints();
 112   math_Matrix B(1, classe, 1, nbcol, 0.0);
 113   Standard_Integer bdeb = 1, bfin = classe;
 114   AppParCurves_Constraint myFirstC = FirstCons, myLastC = LastCons;
 115   SSP.GetNumberOfPoints(myNbP, myNbP2d);
 116
 117   Standard_Integer i2plus1, i2plus2;
 118   myNbdiscret = myNbPoints;
 119   NCollection_Array1<gp_Pnt>   aTabP(1, Max (myNbP, 1));
 120   NCollection_Array1<gp_Pnt2d> aTabP2d(1, Max (myNbP2d, 1));
 121   NCollection_Array1<gp_Vec>   aTabV(1, Max (myNbP, 1));
 122   NCollection_Array1<gp_Vec2d> aTabV2d(1, Max (myNbP2d, 1));
 123
 124   for(Standard_Integer aDimIdx = 1; aDimIdx <= myNbP * 3 + myNbP2d * 2; aDimIdx++)
 125   {
 126     SSP.PeriodInformation(aDimIdx,
 127                           myPerInfo(aDimIdx).isPeriodic,
 128                           myPerInfo(aDimIdx).myPeriod);
 129   }
 130
 131   Standard_Boolean Ok;
 132   if (myFirstC == AppParCurves_TangencyPoint)
 133   {
 134     Ok = SSP.D1(U0, aTabV2d, aTabV);
 135     if (!Ok) myFirstC = AppParCurves_PassPoint;
 136   }
 137
 138   if (myLastC == AppParCurves_TangencyPoint)
 139   {
 140     Ok = SSP.D1(U1, aTabV2d, aTabV);
 141     if (!Ok) myLastC = AppParCurves_PassPoint;
 142   }
 143
 144   // Compute control points params on which approximation will be built.
 145   math_Vector GaussP(1, myNbPoints), GaussW(1, myNbPoints);
 146   math::GaussPoints(myNbPoints, GaussP);
 147   math::GaussWeights(myNbPoints, GaussW);
 148   math_Vector TheWeights(1, myNbPoints), VBParam(1, myNbPoints);
 149   dU = 0.5*(U1-U0);
 150   for (i = FirstP; i <= LastP; i++)
 151   {
 152     U  = 0.5 * (U1 + U0) + dU * GaussP(i);
 153     if (i <=  (myNbPoints+1)/2)
 154     {
 155       myParam(LastP - i + 1)  = U;
 156       VBParam(LastP - i + 1)  = 0.5 * (1 + GaussP(i));
 157       TheWeights(LastP - i + 1) = 0.5 * GaussW(i);
 158     }
 159     else
 160     {
 161       VBParam(i - (myNbPoints + 1) / 2)  = 0.5*(1 + GaussP(i));
 162       myParam(i - (myNbPoints + 1) / 2) = U;
 163       TheWeights(i - (myNbPoints+ 1) / 2) = 0.5 * GaussW(i);
 164     }
 165   }
 166
 167   // Compute control points.
 168   for (i = FirstP; i <= LastP; i++)
 169   {
 170     U = myParam(i);
 171     SSP.Value(U, aTabP2d, aTabP);
 172
 173     i2 = 1;
 174     for (j = 1; j <= myNbP; j++)
 175     {
 176       (aTabP(j)).Coord(myPoints(i, i2), myPoints(i, i2+1), myPoints(i, i2+2));
 177       i2 += 3;
 178     }
 179     for (j = 1; j <= myNbP2d; j++)
 180     {
 181       (aTabP2d(j)).Coord(myPoints(i, i2), myPoints(i, i2+1));
 182       i2 += 2;
 183     }
 184   }
 185
 186   // Fix possible "period jump".
 187   Standard_Integer aMaxDim = 3 * myNbP + 2 * myNbP2d;
 188   for(Standard_Integer aDimIdx = 1; aDimIdx <= aMaxDim; aDimIdx++)
 189   {
 190     if (myPerInfo(aDimIdx).isPeriodic &&
 191         Abs (myPoints(1, aDimIdx) - myPoints(2, aDimIdx)) > myPerInfo(aDimIdx).myPeriod / 2.01 &&
 192         Abs (myPoints(2, aDimIdx) - myPoints(3, aDimIdx)) < myPerInfo(aDimIdx).myPeriod / 2.01)
 193     {
 194       Standard_Real aPeriodMult = (myPoints(1, aDimIdx) < myPoints(2, aDimIdx)) ? 1.0 : -1.0;
 195       Standard_Real aNewParam = myPoints(1, aDimIdx) + aPeriodMult * myPerInfo(aDimIdx).myPeriod;
 196       myPoints(1, aDimIdx) = aNewParam;
 197     }
 198   }
 199   for (Standard_Integer aPntIdx = 1; aPntIdx < myNbPoints; aPntIdx++)
 200   {
 201     for(Standard_Integer aDimIdx = 1; aDimIdx <= aMaxDim; aDimIdx++)
 202     {
 203       if (myPerInfo(aDimIdx).isPeriodic &&
 204         Abs ( myPoints(aPntIdx, aDimIdx) - myPoints(aPntIdx + 1, aDimIdx) ) > myPerInfo(aDimIdx).myPeriod / 2.01)
 205       {
 206         Standard_Real aPeriodMult = (myPoints(aPntIdx, aDimIdx) > myPoints(aPntIdx + 1, aDimIdx)) ? 1.0 : -1.0;
 207         Standard_Real aNewParam = myPoints(aPntIdx + 1, aDimIdx) + aPeriodMult * myPerInfo(aDimIdx).myPeriod;
 208         myPoints(aPntIdx + 1, aDimIdx) = aNewParam;
 209       }
 210     }
 211   }
 212
 213   VBernstein(classe, myNbPoints, myVB);
 214
 215   // Traitement du second membre:
 216   NCollection_Array1<Standard_Real> tmppoints(1, nbcol);
 217
 218   for (c = 1; c <= classe; c++)
 219   {
 220     tmppoints.Init(0.0);
 221     for (i = 1; i <= myNbPoints; i++)
 222     {
 223       Coeff = TheWeights(i) * myVB(c, i);
 224       for (j = 1; j <= nbcol; j++)
 225       {
 226         tmppoints(j) += myPoints(i, j)*Coeff;
 227       }
 228     }
 229     for (k = 1; k <= nbcol; k++)
 230     {
 231       B(c, k) += tmppoints(k);
 232     }
 233   }
 234
 235   if (myFirstC == AppParCurves_NoConstraint &&
 236       myLastC  == AppParCurves_NoConstraint) {
 237
 238     math_Matrix InvM(1, classe, 1, classe);
 239     InvMMatrix(classe, InvM);
 240     // Calcul direct des poles:
 241
 242     for (i = 1; i <= classe; i++) {
 243       for (j = 1; j <= classe; j++) {
 244         IBij = InvM(i, j);
 245         for (k = 1; k <= nbcol; k++) {
 246           myPoles(i, k)   += IBij * B(j, k);
 247         }
 248       }
 249     }
 250   }
 251
 252
 253   else
 254   {
 255     math_Matrix M(1, classe, 1, classe);
 256     MMatrix(classe, M);
 257     NCollection_Array1<gp_Pnt2d> aFixP2d(1, Max (myNbP2d, 1));
 258     NCollection_Array1<gp_Pnt>   aFixP(1, Max (myNbP, 1));
 259
 260     if (myFirstC == AppParCurves_PassPoint ||
 261         myFirstC == AppParCurves_TangencyPoint)
 262     {
 263       SSP.Value(U0, aTabP2d, aTabP);
 264       FixSingleBorderPoint(SSP, U0, U0, U1, aFixP2d, aFixP);
 265
 266       i2 = 1;
 267       for (k = 1; k<= myNbP; k++)
 268       {
 269         if (aFixP(k).Distance(aTabP(k)) > 0.1)
 270           (aFixP(k)).Coord(myPoles(1, i2), myPoles(1, i2 + 1), myPoles(1, i2 + 2));
 271         else
 272           (aTabP(k)).Coord(myPoles(1, i2), myPoles(1, i2 + 1), myPoles(1, i2 + 2));
 273         i2 += 3;
 274       }
 275       for (k = 1; k<= myNbP2d; k++)
 276       {
 277         if (aFixP2d(k).Distance(aTabP2d(k)) > 0.1)
 278           (aFixP2d(k)).Coord(myPoles(1, i2), myPoles(1, i2 + 1));
 279         else
 280           (aTabP2d(k)).Coord(myPoles(1, i2), myPoles(1, i2 + 1));
 281         i2 += 2;
 282       }
 283
 284       for (Standard_Integer aDimIdx = 1; aDimIdx <= aMaxDim; aDimIdx++)
 285       {
 286         if (myPerInfo(aDimIdx).isPeriodic &&
 287             Abs ( myPoles(1, aDimIdx) - myPoints(1, aDimIdx) ) > myPerInfo(aDimIdx).myPeriod / 2.01 )
 288         {
 289           Standard_Real aMult = myPoles(1, aDimIdx) < myPoints(1, aDimIdx)? 1.0: -1.0;
 290           myPoles(1,aDimIdx) += aMult * myPerInfo(aDimIdx).myPeriod;
 291         }
 292       }
 293     }
 294
 295     if (myLastC == AppParCurves_PassPoint ||
 296         myLastC == AppParCurves_TangencyPoint)
 297     {
 298       SSP.Value(U1, aTabP2d, aTabP);
 299       FixSingleBorderPoint(SSP, U1, U0, U1, aFixP2d, aFixP);
 300
 301       i2 = 1;
 302       for (k = 1; k<= myNbP; k++)
 303       {
 304         if (aFixP(k).Distance(aTabP(k)) > 0.1)
 305           (aFixP(k)).Coord(myPoles(classe, i2), myPoles(classe, i2 + 1), myPoles(classe, i2 + 2));
 306         else
 307           (aTabP(k)).Coord(myPoles(classe, i2), myPoles(classe, i2 + 1), myPoles(classe, i2 + 2));
 308         i2 += 3;
 309       }
 310       for (k = 1; k<= myNbP2d; k++)
 311       {
 312         if (aFixP2d(k).Distance(aTabP2d(k)) > 0.1)
 313           (aFixP2d(k)).Coord(myPoles(classe, i2), myPoles(classe, i2 + 1));
 314         else
 315           (aTabP2d(k)).Coord(myPoles(classe, i2), myPoles(classe, i2 + 1));
 316         i2 += 2;
 317       }
 318
 319
 320       for (Standard_Integer aDimIdx = 1; aDimIdx <= 2; aDimIdx++)
 321       {
 322         if (myPerInfo(aDimIdx).isPeriodic &&
 323           Abs ( myPoles(classe, aDimIdx) - myPoints(myNbPoints, aDimIdx) ) > myPerInfo(aDimIdx).myPeriod / 2.01 )
 324         {
 325           Standard_Real aMult = myPoles(classe, aDimIdx) < myPoints(myNbPoints, aDimIdx)? 1.0: -1.0;
 326           myPoles(classe,aDimIdx) += aMult * myPerInfo(aDimIdx).myPeriod;
 327         }
 328       }
 329     }
 330
 331     if (myFirstC == AppParCurves_PassPoint) {
 332       bdeb = 2;
 333       // mise a jour du second membre:
 334       for (i = 1; i <= classe; i++) {
 335         Coeff = M(i, 1);
 336         for (k = 1; k <= nbcol; k++) {
 337           B(i, k) -= myPoles(1, k)*Coeff;
 338         }
 339       }
 340     }
 341
 342     if (myLastC == AppParCurves_PassPoint) {
 343       bfin = cl1;
 344       for (i = 1; i <= classe; i++) {
 345         Coeff = M(i, classe);
 346         for (k = 1; k <= nbcol; k++) {
 347           B(i, k) -= myPoles(classe, k)*Coeff;
 348         }
 349       }
 350     }
 351
 352     if (myFirstC == AppParCurves_TangencyPoint) {
 353       // On fixe le second pole::
 354       bdeb = 3;
 355       SSP.D1(U0, aTabV2d, aTabV);
 356
 357       i2 = 1;
 358       Coeff = (U1-U0)/myDegre;
 359       for (k = 1; k<= myNbP; k++) {
 360         i2plus1 = i2+1; i2plus2 = i2+2;
 361         myPoles(2, i2)      = myPoles(1, i2)      + aTabV(k).X()*Coeff;
 362         myPoles(2, i2plus1) = myPoles(1, i2plus1) + aTabV(k).Y()*Coeff;
 363         myPoles(2, i2plus2) = myPoles(1, i2plus2) + aTabV(k).Z()*Coeff;
 364         i2 += 3;
 365       }
 366       for (k = 1; k<= myNbP2d; k++) {
 367         i2plus1 = i2+1;
 368         myPoles(2, i2)      = myPoles(1, i2)      + aTabV2d(k).X()*Coeff;
 369         myPoles(2, i2plus1) = myPoles(1, i2plus1) + aTabV2d(k).Y()*Coeff;
 370         i2 += 2;
 371       }
 372
 373       for (i = 1; i <= classe; i++) {
 374         Coeff = M(i, 1); Coeff2 = M(i, 2);
 375         for (k = 1; k <= nbcol; k++) {
 376           B(i, k) -= myPoles(1, k)*Coeff+myPoles(2, k)*Coeff2;
 377         }
 378       }
 379     }
 380
 381     if (myLastC == AppParCurves_TangencyPoint) {
 382       bfin = classe-2;
 383       SSP.D1(U1, aTabV2d, aTabV);
 384       i2 = 1;
 385       Coeff = (U1-U0)/myDegre;
 386       for (k = 1; k<= myNbP; k++) {
 387         i2plus1 = i2+1; i2plus2 = i2+2;
 388         myPoles(cl1,i2)      = myPoles(classe, i2)      - aTabV(k).X()*Coeff;
 389         myPoles(cl1,i2plus1) = myPoles(classe, i2plus1) - aTabV(k).Y()*Coeff;
 390         myPoles(cl1,i2plus2) = myPoles(classe, i2plus2) - aTabV(k).Z()*Coeff;
 391         i2 += 3;
 392       }
 393       for (k = 1; k<= myNbP2d; k++) {
 394         i2plus1 = i2+1;
 395         myPoles(cl1,i2) = myPoles(classe, i2) - aTabV2d(k).X()*Coeff;
 396         myPoles(cl1,i2plus1) = myPoles(classe, i2plus1) - aTabV2d(k).Y()*Coeff;
 397         i2 += 2;
 398       }
 399
 400       for (i = 1; i <= classe; i++) {
 401         Coeff = M(i, classe); Coeff2 = M(i, cl1);
 402         for (k = 1; k <= nbcol; k++) {
 403           B(i, k) -= myPoles(classe, k)*Coeff + myPoles(cl1, k)*Coeff2;
 404         }
 405       }
 406     }
 407
 408
 409     if (bdeb <= bfin) {
 410       math_Matrix B2(bdeb, bfin, 1, B.UpperCol(), 0.0);
 411
 412       for (i = bdeb; i <= bfin; i++) {
 413         for (j = 1; j <= classe; j++) {
 414           Coeff = M(i, j);
 415           for (k = 1; k <= nbcol; k++) {
 416             B2(i, k) += B(j, k)*Coeff;
 417           }
 418         }
 419       }
 420
 421       // Resolution:
 422       // ===========
 423       math_Matrix IBP(bdeb, bfin, bdeb, bfin);
 424
 425       // dans IBPMatrix at IBTMatrix ne sont stockees que les resultats pour
 426       // une classe inferieure ou egale a 26 (pour l instant du moins.)
 427
 428       if (bdeb == 2 && bfin == classe-1 && classe <= 26) {
 429         IBPMatrix(classe, IBP);
 430       }
 431       else if (bdeb == 3 && bfin == classe-2 && classe <= 26) {
 432         IBTMatrix(classe, IBP);
 433       }
 434       else {
 435         math_Matrix MP(1, classe, bdeb, bfin);
 436         for (i = 1; i <= classe; i++) {
 437           for (j = bdeb; j <= bfin; j++) {
 438             MP(i, j) = M(i, j);
 439           }
 440         }
 441         math_Matrix IBP1(bdeb, bfin, bdeb, bfin);
 442         IBP1 = MP.Transposed()*MP;
 443         IBP = IBP1.Inverse();
 444       }
 445
 446       myDone = Standard_True;
 447       for (i = bdeb; i <= bfin; i++) {
 448         for (j = bdeb; j <= bfin; j++) {
 449           IBPij = IBP(i, j);;
 450           for (k = 1; k<= nbcol; k++) {
 451             myPoles(i, k)   += IBPij * B2(j, k);
 452           }
 453         }
 454       }
 455     }
 456   }
 457 }
 458
 459 //=======================================================================
 460 //function : Value
 461 //purpose  :
 462 //=======================================================================
 463
 464 const AppParCurves_MultiCurve& AppCont_LeastSquare::Value()
 465 {
 466
 467   Standard_Integer i, j, j2;
 468   gp_Pnt Pt;
 469   gp_Pnt2d Pt2d;
 470   Standard_Integer ideb = 1, ifin = myDegre+1;
 471
 472   // On met le resultat dans les curves correspondantes
 473   for (i = ideb; i <= ifin; i++) {
 474     j2 = 1;
 475     AppParCurves_MultiPoint MPole(myNbP, myNbP2d);
 476     for (j = 1; j <= myNbP; j++) {
 477       Pt.SetCoord(myPoles(i, j2), myPoles(i, j2+1), myPoles(i,j2+2));
 478       MPole.SetPoint(j, Pt);
 479       j2 += 3;
 480     }
 481     for (j = myNbP+1;j <= myNbP+myNbP2d; j++) {
 482       Pt2d.SetCoord(myPoles(i, j2), myPoles(i, j2+1));
 483       MPole.SetPoint2d(j, Pt2d);
 484       j2 += 2;
 485     }
 486     mySCU.SetValue(i, MPole);
 487   }
 488   return mySCU;
 489 }
 490
 491
 492
 493 //=======================================================================
 494 //function : Error
 495 //purpose  :
 496 //=======================================================================
 497
 498 void AppCont_LeastSquare::Error(Standard_Real& F,
 499                                 Standard_Real& MaxE3d,
 500                                 Standard_Real& MaxE2d) const
 501 {
 502   Standard_Integer i, j, k, c, i2, classe = myDegre + 1;
 503   Standard_Real Coeff, err3d = 0.0, err2d = 0.0;
 504   Standard_Integer ncol = myPoints.UpperCol() - myPoints.LowerCol() + 1;
 505
 506   math_Matrix MyPoints(1, myNbdiscret, 1, ncol);
 507   MyPoints = myPoints;
 508
 509   MaxE3d = MaxE2d = F = 0.0;
 510
 511   NCollection_Array1<Standard_Real> tmppoles(1, ncol);
 512
 513   for (c = 1; c <= classe; c++)
 514   {
 515     for (k = 1; k <= ncol; k++)
 516     {
 517       tmppoles(k) = myPoles(c, k);
 518     }
 519     for (i = 1; i <= myNbdiscret; i++)
 520     {
 521       Coeff = myVB(c, i);
 522       for (j = 1; j <= ncol; j++)
 523       {
 524         MyPoints(i, j) -= tmppoles(j) * Coeff;
 525       }
 526     }
 527   }
 528
 529   Standard_Real e1, e2, e3;
 530   for (i = 1; i <= myNbdiscret; i++)
 531   {
 532     i2 = 1;
 533     for (j = 1; j<= myNbP; j++) {
 534       e1 = MyPoints(i, i2);
 535       e2 = MyPoints(i, i2+1);
 536       e3 = MyPoints(i, i2+2);
 537       err3d = e1*e1+e2*e2+e3*e3;
 538       MaxE3d = Max(MaxE3d, err3d);
 539       F += err3d;
 540       i2 += 3;
 541     }
 542     for (j = 1; j<= myNbP2d; j++) {
 543       e1 = MyPoints(i, i2);
 544       e2 = MyPoints(i, i2+1);
 545       err2d = e1*e1+e2*e2;
 546       MaxE2d = Max(MaxE2d, err2d);
 547       F += err2d;
 548       i2 += 2;
 549     }
 550   }
 551
 552   MaxE3d = Sqrt(MaxE3d);
 553   MaxE2d = Sqrt(MaxE2d);
 554
 555 }
 556
 557
 558 //=======================================================================
 559 //function : IsDone
 560 //purpose  :
 561 //=======================================================================
 562
 563 Standard_Boolean AppCont_LeastSquare::IsDone() const
 564 {
 565   return myDone;
 566 }