src/AppParCurves/AppParCurves_Gradient.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
   7 // under the terms of the GNU Lesser General Public version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 // lpa, le 11/09/91
  16
  17
  18 // Application de la methode du gradient corrige pour minimiser
  19 // F  = somme(||C(ui, Poles(ui)) - ptli||2.
  20 // La methode de gradient conjugue est programmee dans la bibliotheque
  21 // mathematique: math_BFGS.
  22 // cet algorithme doit etre appele uniquement lorsque on a affaire a un set
  23 // de points contraints (ailleurs qu aux extremites). En effet, l appel de la
  24 // fonction F a minimiser implique un appel a ParLeastSquare et ResConstraint.
  25 // Si ce n est pas le cas, l appel a ResConstraint est equivalent a une
  26 // seconde resolution par les moindres carres donc beaucoup de temps perdu.
  27
  28
  29 #define No_Standard_RangeError
  30 #define No_Standard_OutOfRange
  31
  32 #include <AppParCurves_Constraint.hxx>
  33 #include <math_BFGS.hxx>
  34 #include <StdFail_NotDone.hxx>
  35 #include <AppParCurves_MultiPoint.hxx>
  36 #include <gp_Pnt.hxx>
  37 #include <gp_Pnt2d.hxx>
  38 #include <gp_Vec.hxx>
  39 #include <gp_Vec2d.hxx>
  40 #include <TColgp_Array1OfPnt.hxx>
  41 #include <TColgp_Array1OfPnt2d.hxx>
  42 #include <TColgp_Array1OfVec.hxx>
  43 #include <TColgp_Array1OfVec2d.hxx>
  44 #include <BSplCLib.hxx>
  45 #include <PLib.hxx>
  46
  47 // #define AppParCurves_Gradient_BFGS BFGS_/**/AppParCurves_Gradient
  48
  49
  50
  51 AppParCurves_Gradient::
  52    AppParCurves_Gradient(const MultiLine& SSP,
  53          const Standard_Integer FirstPoint,
  54          const Standard_Integer LastPoint,
  55          const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
  56          math_Vector& Parameters,
  57          const Standard_Integer Deg,
  58          const Standard_Real Tol3d,
  59          const Standard_Real Tol2d,
  60          const Standard_Integer NbIterations):
  61          ParError(FirstPoint, LastPoint,0.0) {
  62
  63 //  Standard_Boolean grad = Standard_True;
  64   Standard_Integer j, k, i2, l;
  65   Standard_Real UF, DU, Fval = 0.0, FU, DFU;
  66   Standard_Integer nbP3d = ToolLine::NbP3d(SSP);
  67   Standard_Integer nbP2d = ToolLine::NbP2d(SSP);
  68   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
  69   Standard_Integer nbP = nbP3d + nbP2d;
  70 //  gp_Pnt Pt, P1, P2;
  71   gp_Pnt Pt;
  72 //  gp_Pnt2d Pt2d, P12d, P22d;
  73   gp_Pnt2d Pt2d;
  74 //  gp_Vec V1, V2, MyV;
  75   gp_Vec V1, MyV;
  76 //  gp_Vec2d V12d, V22d, MyV2d;
  77   gp_Vec2d V12d, MyV2d;
  78   Done = Standard_False;
  79
  80   if (nbP3d == 0) mynbP3d = 1;
  81   if (nbP2d == 0) mynbP2d = 1;
  82   TColgp_Array1OfPnt TabP(1, mynbP3d);
  83   TColgp_Array1OfPnt2d TabP2d(1, mynbP2d);
  84   TColgp_Array1OfVec TabV(1, mynbP3d);
  85   TColgp_Array1OfVec2d TabV2d(1, mynbP2d);
  86
  87   // Calcul de la fonction F= somme(||C(ui)-Ptli||2):
  88   // Appel a une fonction heritant de MultipleVarFunctionWithGradient
  89   // pour calculer F et grad_F.
  90   // ================================================================
  91
  92   AppParCurves_ParFunction MyF(SSP, FirstPoint,LastPoint, TheConstraints, Parameters, Deg);
  93
  94
  95   if (!MyF.Value(Parameters, Fval)) {
  96     Done = Standard_False;
  97     return;
  98   }
  99
 100   SCU = MyF.CurveValue();
 101   Standard_Integer deg = SCU.NbPoles()-1;
 102   TColgp_Array1OfPnt   TabPole(1, deg+1),   TabCoef(1, deg+1);
 103   TColgp_Array1OfPnt2d TabPole2d(1, deg+1), TabCoef2d(1, deg+1);
 104   TColgp_Array1OfPnt    TheCoef(1, (deg+1)*mynbP3d);
 105   TColgp_Array1OfPnt2d  TheCoef2d(1, (deg+1)*mynbP2d);
 106
 107
 108   // Stockage des Poles des courbes pour projeter:
 109   // ============================================
 110   i2 = 0;
 111   for (k = 1; k <= nbP3d; k++) {
 112     SCU.Curve(k, TabPole);
 113     BSplCLib::PolesCoefficients(TabPole, PLib::NoWeights(),
 114                                 TabCoef, PLib::NoWeights());
 115     for (j=1; j<=deg+1; j++) TheCoef(j+i2) = TabCoef(j);
 116     i2 += deg+1;
 117   }
 118   i2 = 0;
 119   for (k = 1; k <= nbP2d; k++) {
 120     SCU.Curve(nbP3d+k, TabPole2d);
 121     BSplCLib::PolesCoefficients(TabPole2d, PLib::NoWeights(),
 122                                 TabCoef2d, PLib::NoWeights());
 123     for (j=1; j<=deg+1; j++) TheCoef2d(j+i2) = TabCoef2d(j);
 124     i2 += deg+1;
 125   }
 126
 127   //  Une iteration rapide de projection est faite par la methode de
 128   //  Rogers & Fog 89, methode equivalente a Hoschek 88 qui ne necessite pas
 129   //  le calcul de D2.
 130
 131
 132   // Iteration de Projection:
 133   // =======================
 134   for (j = FirstPoint+1; j <= LastPoint-1; j++) {
 135     UF = Parameters(j);
 136     if (nbP != 0 && nbP2d != 0) ToolLine::Value(SSP, j, TabP, TabP2d);
 137     else if (nbP2d != 0)        ToolLine::Value(SSP, j, TabP2d);
 138     else                        ToolLine::Value(SSP, j, TabP);
 139
 140     FU = 0.0;
 141     DFU = 0.0;
 142     i2 = 0;
 143     for (k = 1; k <= nbP3d; k++) {
 144       for (l=1; l<=deg+1; l++) TabCoef(l) = TheCoef(l+i2);
 145       i2 += deg+1;
 146       BSplCLib::CoefsD1(UF, TabCoef, PLib::NoWeights(), Pt, V1);
 147       MyV = gp_Vec(Pt, TabP(k));
 148       FU += MyV*V1;
 149       DFU += V1.SquareMagnitude();
 150     }
 151     i2 = 0;
 152     for (k = 1; k <= nbP2d; k++) {
 153       for (l=1; l<=deg+1; l++) TabCoef2d(l) = TheCoef2d(l+i2);
 154       i2 += deg+1;
 155       BSplCLib::CoefsD1(UF, TabCoef2d, PLib::NoWeights(), Pt2d, V12d);
 156       MyV2d = gp_Vec2d(Pt2d, TabP2d(k));
 157       FU += MyV2d*V12d;
 158       DFU += V12d.SquareMagnitude();
 159     }
 160
 161     if (DFU >= RealEpsilon()) {
 162       DU = FU/DFU;
 163       DU = Sign(Min(5.e-02, Abs(DU)), DU);
 164       UF += DU;
 165       Parameters(j) = UF;
 166     }
 167   }
 168
 169
 170   if (!MyF.Value(Parameters, Fval)) {
 171     SCU = AppParCurves_MultiCurve();
 172     Done = Standard_False;
 173     return;
 174   }
 175   MError3d = MyF.MaxError3d();
 176   MError2d = MyF.MaxError2d();
 177
 178   if (MError3d<= Tol3d && MError2d <= Tol2d) {
 179     Done = Standard_True;
 180     SCU = MyF.CurveValue();
 181   }
 182   else if (NbIterations != 0) {
 183     // NbIterations de gradient conjugue:
 184     // =================================
 185     Standard_Real Eps = 1.e-07;
 186     AppParCurves_Gradient_BFGS FResol(MyF, Parameters, Tol3d, Tol2d, Eps, NbIterations);
 187     Parameters = MyF.NewParameters();
 188     SCU = MyF.CurveValue();
 189   }
 190
 191
 192   AvError = 0.;
 193   for (j = FirstPoint; j <= LastPoint; j++) {
 194     // Recherche des erreurs maxi et moyenne a un index donne:
 195     for (k = 1; k <= nbP; k++) {
 196       ParError(j) = Max(ParError(j), MyF.Error(j, k));
 197     }
 198     AvError += ParError(j);
 199   }
 200   AvError = AvError/(LastPoint-FirstPoint+1);
 201
 202
 203   MError3d = MyF.MaxError3d();
 204   MError2d = MyF.MaxError2d();
 205   if (MError3d <= Tol3d && MError2d <= Tol2d) {
 206     Done = Standard_True;
 207   }
 208
 209 }
 210
 211
 212
 213 AppParCurves_MultiCurve AppParCurves_Gradient::Value() const {
 214   return SCU;
 215 }
 216
 217
 218 Standard_Boolean AppParCurves_Gradient::IsDone() const {
 219   return Done;
 220 }
 221
 222
 223 Standard_Real AppParCurves_Gradient::Error(const Standard_Integer Index) const {
 224   return ParError(Index);
 225 }
 226
 227 Standard_Real AppParCurves_Gradient::AverageError() const {
 228   return AvError;
 229 }
 230
 231 Standard_Real AppParCurves_Gradient::MaxError3d() const {
 232   return MError3d;
 233 }
 234
 235 Standard_Real AppParCurves_Gradient::MaxError2d() const {
 236   return MError2d;
 237 }
 238
 239