Warning, /include/opencascade/AppParCurves_Gradient.gxx is written in an unsupported language. File is not indexed.
0001 // Copyright (c) 1995-1999 Matra Datavision
0002 // Copyright (c) 1999-2014 OPEN CASCADE SAS
0003 //
0004 // This file is part of Open CASCADE Technology software library.
0005 //
0006 // This library is free software; you can redistribute it and/or modify it under
0007 // the terms of the GNU Lesser General Public License version 2.1 as published
0008 // by the Free Software Foundation, with special exception defined in the file
0009 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0010 // distribution for complete text of the license and disclaimer of any warranty.
0011 //
0012 // Alternatively, this file may be used under the terms of Open CASCADE
0013 // commercial license or contractual agreement.
0014
0015 // lpa, le 11/09/91
0016
0017
0018 // Application de la methode du gradient corrige pour minimiser
0019 // F = somme(||C(ui, Poles(ui)) - ptli||2.
0020 // La methode de gradient conjugue est programmee dans la bibliotheque
0021 // mathematique: math_BFGS.
0022 // cet algorithme doit etre appele uniquement lorsque on a affaire a un set
0023 // de points contraints (ailleurs qu aux extremites). En effet, l appel de la
0024 // fonction F a minimiser implique un appel a ParLeastSquare et ResConstraint.
0025 // Si ce n est pas le cas, l appel a ResConstraint est equivalent a une
0026 // seconde resolution par les moindres carres donc beaucoup de temps perdu.
0027
0028
0029 #define No_Standard_RangeError
0030 #define No_Standard_OutOfRange
0031
0032 #include <AppParCurves_Constraint.hxx>
0033 #include <math_BFGS.hxx>
0034 #include <StdFail_NotDone.hxx>
0035 #include <AppParCurves_MultiPoint.hxx>
0036 #include <gp_Pnt.hxx>
0037 #include <gp_Pnt2d.hxx>
0038 #include <gp_Vec.hxx>
0039 #include <gp_Vec2d.hxx>
0040 #include <TColgp_Array1OfPnt.hxx>
0041 #include <TColgp_Array1OfPnt2d.hxx>
0042 #include <TColgp_Array1OfVec.hxx>
0043 #include <TColgp_Array1OfVec2d.hxx>
0044 #include <BSplCLib.hxx>
0045 #include <PLib.hxx>
0046
0047 // #define AppParCurves_Gradient_BFGS BFGS_/**/AppParCurves_Gradient
0048
0049
0050
0051 AppParCurves_Gradient::
0052 AppParCurves_Gradient(const MultiLine& SSP,
0053 const Standard_Integer FirstPoint,
0054 const Standard_Integer LastPoint,
0055 const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
0056 math_Vector& Parameters,
0057 const Standard_Integer Deg,
0058 const Standard_Real Tol3d,
0059 const Standard_Real Tol2d,
0060 const Standard_Integer NbIterations):
0061 ParError(FirstPoint, LastPoint,0.0),
0062 AvError(0.0),
0063 MError3d(0.0),
0064 MError2d(0.0)
0065 {
0066
0067 // Standard_Boolean grad = Standard_True;
0068 Standard_Integer j, k, i2, l;
0069 Standard_Real UF, DU, Fval = 0.0, FU, DFU;
0070 Standard_Integer nbP3d = ToolLine::NbP3d(SSP);
0071 Standard_Integer nbP2d = ToolLine::NbP2d(SSP);
0072 Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
0073 Standard_Integer nbP = nbP3d + nbP2d;
0074 // gp_Pnt Pt, P1, P2;
0075 gp_Pnt Pt;
0076 // gp_Pnt2d Pt2d, P12d, P22d;
0077 gp_Pnt2d Pt2d;
0078 // gp_Vec V1, V2, MyV;
0079 gp_Vec V1, MyV;
0080 // gp_Vec2d V12d, V22d, MyV2d;
0081 gp_Vec2d V12d, MyV2d;
0082 Done = Standard_False;
0083
0084 if (nbP3d == 0) mynbP3d = 1;
0085 if (nbP2d == 0) mynbP2d = 1;
0086 TColgp_Array1OfPnt TabP(1, mynbP3d);
0087 TColgp_Array1OfPnt2d TabP2d(1, mynbP2d);
0088 TColgp_Array1OfVec TabV(1, mynbP3d);
0089 TColgp_Array1OfVec2d TabV2d(1, mynbP2d);
0090
0091 // Calcul de la fonction F= somme(||C(ui)-Ptli||2):
0092 // Appel a une fonction heritant de MultipleVarFunctionWithGradient
0093 // pour calculer F et grad_F.
0094 // ================================================================
0095
0096 AppParCurves_ParFunction MyF(SSP, FirstPoint,LastPoint, TheConstraints, Parameters, Deg);
0097
0098
0099 if (!MyF.Value(Parameters, Fval)) {
0100 Done = Standard_False;
0101 return;
0102 }
0103
0104 SCU = MyF.CurveValue();
0105 Standard_Integer deg = SCU.NbPoles()-1;
0106 TColgp_Array1OfPnt TabPole(1, deg+1), TabCoef(1, deg+1);
0107 TColgp_Array1OfPnt2d TabPole2d(1, deg+1), TabCoef2d(1, deg+1);
0108 TColgp_Array1OfPnt TheCoef(1, (deg+1)*mynbP3d);
0109 TColgp_Array1OfPnt2d TheCoef2d(1, (deg+1)*mynbP2d);
0110
0111
0112 // Stockage des Poles des courbes pour projeter:
0113 // ============================================
0114 i2 = 0;
0115 for (k = 1; k <= nbP3d; k++) {
0116 SCU.Curve(k, TabPole);
0117 BSplCLib::PolesCoefficients(TabPole, PLib::NoWeights(),
0118 TabCoef, PLib::NoWeights());
0119 for (j=1; j<=deg+1; j++) TheCoef(j+i2) = TabCoef(j);
0120 i2 += deg+1;
0121 }
0122 i2 = 0;
0123 for (k = 1; k <= nbP2d; k++) {
0124 SCU.Curve(nbP3d+k, TabPole2d);
0125 BSplCLib::PolesCoefficients(TabPole2d, PLib::NoWeights(),
0126 TabCoef2d, PLib::NoWeights());
0127 for (j=1; j<=deg+1; j++) TheCoef2d(j+i2) = TabCoef2d(j);
0128 i2 += deg+1;
0129 }
0130
0131 // Une iteration rapide de projection est faite par la methode de
0132 // Rogers & Fog 89, methode equivalente a Hoschek 88 qui ne necessite pas
0133 // le calcul de D2.
0134
0135
0136 // Iteration de Projection:
0137 // =======================
0138 for (j = FirstPoint+1; j <= LastPoint-1; j++) {
0139 UF = Parameters(j);
0140 if (nbP != 0 && nbP2d != 0) ToolLine::Value(SSP, j, TabP, TabP2d);
0141 else if (nbP2d != 0) ToolLine::Value(SSP, j, TabP2d);
0142 else ToolLine::Value(SSP, j, TabP);
0143
0144 FU = 0.0;
0145 DFU = 0.0;
0146 i2 = 0;
0147 for (k = 1; k <= nbP3d; k++) {
0148 for (l=1; l<=deg+1; l++) TabCoef(l) = TheCoef(l+i2);
0149 i2 += deg+1;
0150 BSplCLib::CoefsD1(UF, TabCoef, BSplCLib::NoWeights(), Pt, V1);
0151 MyV = gp_Vec(Pt, TabP(k));
0152 FU += MyV*V1;
0153 DFU += V1.SquareMagnitude();
0154 }
0155 i2 = 0;
0156 for (k = 1; k <= nbP2d; k++) {
0157 for (l=1; l<=deg+1; l++) TabCoef2d(l) = TheCoef2d(l+i2);
0158 i2 += deg+1;
0159 BSplCLib::CoefsD1(UF, TabCoef2d, BSplCLib::NoWeights(), Pt2d, V12d);
0160 MyV2d = gp_Vec2d(Pt2d, TabP2d(k));
0161 FU += MyV2d*V12d;
0162 DFU += V12d.SquareMagnitude();
0163 }
0164
0165 if (DFU >= RealEpsilon()) {
0166 DU = FU/DFU;
0167 DU = Sign(Min(5.e-02, Abs(DU)), DU);
0168 UF += DU;
0169 Parameters(j) = UF;
0170 }
0171 }
0172
0173
0174 if (!MyF.Value(Parameters, Fval)) {
0175 SCU = AppParCurves_MultiCurve();
0176 Done = Standard_False;
0177 return;
0178 }
0179 MError3d = MyF.MaxError3d();
0180 MError2d = MyF.MaxError2d();
0181
0182 if (MError3d<= Tol3d && MError2d <= Tol2d) {
0183 Done = Standard_True;
0184 SCU = MyF.CurveValue();
0185 }
0186 else if (NbIterations != 0) {
0187 // NbIterations de gradient conjugue:
0188 // =================================
0189 Standard_Real Eps = 1.e-07;
0190 AppParCurves_Gradient_BFGS FResol(MyF, Parameters, Tol3d, Tol2d, Eps, NbIterations);
0191 Parameters = MyF.NewParameters();
0192 SCU = MyF.CurveValue();
0193 }
0194
0195
0196 AvError = 0.;
0197 for (j = FirstPoint; j <= LastPoint; j++) {
0198 // Recherche des erreurs maxi et moyenne a un index donne:
0199 for (k = 1; k <= nbP; k++) {
0200 ParError(j) = Max(ParError(j), MyF.Error(j, k));
0201 }
0202 AvError += ParError(j);
0203 }
0204 AvError = AvError/(LastPoint-FirstPoint+1);
0205
0206
0207 MError3d = MyF.MaxError3d();
0208 MError2d = MyF.MaxError2d();
0209 if (MError3d <= Tol3d && MError2d <= Tol2d) {
0210 Done = Standard_True;
0211 }
0212
0213 }
0214
0215
0216
0217 AppParCurves_MultiCurve AppParCurves_Gradient::Value() const {
0218 return SCU;
0219 }
0220
0221
0222 Standard_Boolean AppParCurves_Gradient::IsDone() const {
0223 return Done;
0224 }
0225
0226
0227 Standard_Real AppParCurves_Gradient::Error(const Standard_Integer Index) const {
0228 return ParError(Index);
0229 }
0230
0231 Standard_Real AppParCurves_Gradient::AverageError() const {
0232 return AvError;
0233 }
0234
0235 Standard_Real AppParCurves_Gradient::MaxError3d() const {
0236 return MError3d;
0237 }
0238
0239 Standard_Real AppParCurves_Gradient::MaxError2d() const {
0240 return MError2d;
0241 }
0242
0243
