Warning, /include/opencascade/AppParCurves_BSpGradient.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
0023
0024 #define No_Standard_RangeError
0025 #define No_Standard_OutOfRange
0026
0027 #include <AppParCurves_Constraint.hxx>
0028 #include <AppParCurves_ConstraintCouple.hxx>
0029 #include <math_BFGS.hxx>
0030 #include <StdFail_NotDone.hxx>
0031 #include <AppParCurves_MultiPoint.hxx>
0032 #include <gp_Pnt.hxx>
0033 #include <gp_Pnt2d.hxx>
0034 #include <gp_Vec.hxx>
0035 #include <gp_Vec2d.hxx>
0036 #include <TColgp_Array1OfPnt.hxx>
0037 #include <TColgp_Array1OfPnt2d.hxx>
0038 #include <TColgp_Array1OfVec.hxx>
0039 #include <TColgp_Array1OfVec2d.hxx>
0040
0041 static AppParCurves_Constraint FirstConstraint
0042 (const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
0043 const Standard_Integer FirstPoint)
0044 {
0045 Standard_Integer i, myindex;
0046 Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
0047 AppParCurves_ConstraintCouple mycouple;
0048 AppParCurves_Constraint Cons = AppParCurves_NoConstraint;
0049
0050 for (i = low; i <= upp; i++) {
0051 mycouple = TheConstraints->Value(i);
0052 Cons = mycouple.Constraint();
0053 myindex = mycouple.Index();
0054 if (myindex == FirstPoint) {
0055 break;
0056 }
0057 }
0058 return Cons;
0059 }
0060
0061
0062 static AppParCurves_Constraint LastConstraint
0063 (const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
0064 const Standard_Integer LastPoint)
0065 {
0066 Standard_Integer i, myindex;
0067 Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
0068 AppParCurves_ConstraintCouple mycouple;
0069 AppParCurves_Constraint Cons = AppParCurves_NoConstraint;
0070
0071 for (i = low; i <= upp; i++) {
0072 mycouple = TheConstraints->Value(i);
0073 Cons = mycouple.Constraint();
0074 myindex = mycouple.Index();
0075 if (myindex == LastPoint) {
0076 break;
0077 }
0078 }
0079 return Cons;
0080 }
0081
0082
0083
0084
0085
0086 AppParCurves_BSpGradient::
0087 AppParCurves_BSpGradient(const MultiLine& SSP,
0088 const Standard_Integer FirstPoint,
0089 const Standard_Integer LastPoint,
0090 const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
0091 math_Vector& Parameters,
0092 const TColStd_Array1OfReal& Knots,
0093 const TColStd_Array1OfInteger& Mults,
0094 const Standard_Integer Deg,
0095 const Standard_Real Tol3d,
0096 const Standard_Real Tol2d,
0097 const Standard_Integer NbIterations):
0098 ParError(FirstPoint, LastPoint,0.0),
0099 mylambda1(0.0),
0100 mylambda2(0.0),
0101 myIsLambdaDefined(Standard_False)
0102 {
0103 Perform(SSP, FirstPoint, LastPoint, TheConstraints, Parameters,
0104 Knots, Mults, Deg, Tol3d, Tol2d, NbIterations);
0105 }
0106
0107
0108 AppParCurves_BSpGradient::
0109 AppParCurves_BSpGradient(const MultiLine& SSP,
0110 const Standard_Integer FirstPoint,
0111 const Standard_Integer LastPoint,
0112 const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
0113 math_Vector& Parameters,
0114 const TColStd_Array1OfReal& Knots,
0115 const TColStd_Array1OfInteger& Mults,
0116 const Standard_Integer Deg,
0117 const Standard_Real Tol3d,
0118 const Standard_Real Tol2d,
0119 const Standard_Integer NbIterations,
0120 const Standard_Real lambda1,
0121 const Standard_Real lambda2):
0122 ParError(FirstPoint, LastPoint,0.0),
0123 mylambda1(lambda1),
0124 mylambda2(lambda2),
0125 myIsLambdaDefined(Standard_True)
0126 {
0127 Perform(SSP, FirstPoint, LastPoint, TheConstraints, Parameters,
0128 Knots, Mults, Deg, Tol3d, Tol2d, NbIterations);
0129 }
0130
0131
0132
0133
0134
0135
0136
0137 void AppParCurves_BSpGradient::
0138 Perform(const MultiLine& SSP,
0139 const Standard_Integer FirstPoint,
0140 const Standard_Integer LastPoint,
0141 const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
0142 math_Vector& Parameters,
0143 const TColStd_Array1OfReal& Knots,
0144 const TColStd_Array1OfInteger& Mults,
0145 const Standard_Integer Deg,
0146 const Standard_Real Tol3d,
0147 const Standard_Real Tol2d,
0148 const Standard_Integer NbIterations)
0149 {
0150
0151 // Standard_Boolean grad = Standard_True;
0152 Standard_Integer i, j, k, i2, l;
0153 Standard_Real UF, DU, Fval = 0.0, FU, DFU;
0154 Standard_Integer nbP3d = ToolLine::NbP3d(SSP);
0155 Standard_Integer nbP2d = ToolLine::NbP2d(SSP);
0156 Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
0157 Standard_Integer nbP = nbP3d + nbP2d;
0158 // gp_Pnt Pt, P1, P2;
0159 gp_Pnt Pt;
0160 // gp_Pnt2d Pt2d, P12d, P22d;
0161 gp_Pnt2d Pt2d;
0162 // gp_Vec V1, V2, MyV;
0163 gp_Vec V1, MyV;
0164 // gp_Vec2d V12d, V22d, MyV2d;
0165 gp_Vec2d V12d, MyV2d;
0166 Done = Standard_False;
0167
0168 if (nbP3d == 0) mynbP3d = 1;
0169 if (nbP2d == 0) mynbP2d = 1;
0170 TColgp_Array1OfPnt TabP(1, mynbP3d);
0171 TColgp_Array1OfPnt2d TabP2d(1, mynbP2d);
0172 TColgp_Array1OfVec TabV(1, mynbP3d);
0173 TColgp_Array1OfVec2d TabV2d(1, mynbP2d);
0174
0175 // Calcul de la fonction F= somme(||C(ui)-Ptli||2):
0176 // Appel a une fonction heritant de MultipleVarFunctionWithGradient
0177 // pour calculer F et grad_F.
0178 // ================================================================
0179
0180 Standard_Integer nbpoles = -Deg -1;
0181 for (i = Mults.Lower() ;i <= Mults.Upper(); i++) {
0182 nbpoles += Mults(i);
0183 }
0184
0185 TColgp_Array1OfPnt TabPole(1, nbpoles);
0186 TColgp_Array1OfPnt2d TabPole2d(1, nbpoles);
0187 TColgp_Array1OfPnt ThePoles(1, (nbpoles)*mynbP3d);
0188 TColgp_Array1OfPnt2d ThePoles2d(1, (nbpoles)*mynbP2d);
0189
0190
0191 AppParCurves_Constraint
0192 FirstCons = FirstConstraint(TheConstraints, FirstPoint),
0193 LastCons = LastConstraint(TheConstraints, LastPoint);
0194
0195
0196 AppParCurves_BSpParFunction MyF(SSP, FirstPoint,LastPoint, TheConstraints,
0197 Parameters, Knots, Mults, nbpoles);
0198
0199
0200
0201 if (FirstCons >= AppParCurves_TangencyPoint ||
0202 LastCons >= AppParCurves_TangencyPoint) {
0203
0204 if (!myIsLambdaDefined) {
0205 AppParCurves_BSpParLeastSquare thefitt(SSP, Knots, Mults, FirstPoint,
0206 LastPoint, FirstCons, LastCons,
0207 Parameters, nbpoles);
0208 if (FirstCons >= AppParCurves_TangencyPoint) {
0209 mylambda1 = thefitt.FirstLambda();
0210 MyF.SetFirstLambda(mylambda1);
0211 }
0212 if (LastCons >= AppParCurves_TangencyPoint) {
0213 mylambda2 = thefitt.LastLambda();
0214 MyF.SetLastLambda(mylambda2);
0215 }
0216 }
0217 else {
0218 MyF.SetFirstLambda(mylambda1);
0219 MyF.SetLastLambda(mylambda2);
0220 }
0221 }
0222
0223
0224 MyF.Value(Parameters, Fval);
0225 MError3d = MyF.MaxError3d();
0226 MError2d = MyF.MaxError2d();
0227 SCU = MyF.CurveValue();
0228
0229 if (MError3d > Tol3d || MError2d > Tol2d) {
0230
0231 // Stockage des Poles des courbes pour projeter:
0232 // ============================================
0233 i2 = 0;
0234 for (k = 1; k <= nbP3d; k++) {
0235 SCU.Curve(k, TabPole);
0236 for (j=1; j<=nbpoles; j++) ThePoles(j+i2) = TabPole(j);
0237 i2 += nbpoles;
0238 }
0239 i2 = 0;
0240 for (k = 1; k <= nbP2d; k++) {
0241 SCU.Curve(nbP3d+k, TabPole2d);
0242 for (j=1; j<=nbpoles; j++) ThePoles2d(j+i2) = TabPole2d(j);
0243 i2 += nbpoles;
0244 }
0245
0246 // Une iteration rapide de projection est faite par la methode de
0247 // Rogers & Fog 89, methode equivalente a Hoschek 88 qui ne necessite pas
0248 // le calcul de D2.
0249
0250 const math_Matrix& A = MyF.FunctionMatrix();
0251 const math_Matrix& DA = MyF.DerivativeFunctionMatrix();
0252 const math_IntegerVector& Index = MyF.Index();
0253 Standard_Real aa, da, a, b, c, d , e , f, px, py, pz;
0254 Standard_Integer indexdeb, indexfin;
0255
0256 for (j = FirstPoint+1; j <= LastPoint-1; j++) {
0257
0258 UF = Parameters(j);
0259 if (nbP3d != 0 && nbP2d != 0) ToolLine::Value(SSP, j, TabP, TabP2d);
0260 else if (nbP2d != 0) ToolLine::Value(SSP, j, TabP2d);
0261 else ToolLine::Value(SSP, j, TabP);
0262
0263 FU = 0.0;
0264 DFU = 0.0;
0265 i2 = 0;
0266
0267 indexdeb = Index(j) + 1;
0268 indexfin = indexdeb + Deg;
0269
0270 for (k = 1; k <= nbP3d; k++) {
0271 a = b = c = d = e = f = 0.0;
0272 for (l = indexdeb; l <= indexfin; l++) {
0273 Pt = ThePoles(l+i2);
0274 px = Pt.X(); py = Pt.Y(); pz = Pt.Z();
0275 aa = A(j, l); da = DA(j, l);
0276 a += aa* px; d += da* px;
0277 b += aa* py; e += da* py;
0278 c += aa* pz; f += da* pz;
0279 }
0280 Pt.SetCoord(a, b, c);
0281 V1.SetCoord(d, e, f);
0282 i2 += nbpoles;
0283
0284 MyV = gp_Vec(Pt, TabP(k));
0285 FU += MyV*V1;
0286 DFU += V1.SquareMagnitude();
0287 }
0288 i2 = 0;
0289 for (k = 1; k <= nbP2d; k++) {
0290 a = b = d = e = 0.0;
0291 for (l = indexdeb; l <= indexfin; l++) {
0292 Pt2d = ThePoles2d(l+i2);
0293 px = Pt2d.X(); py = Pt2d.Y();
0294 aa = A(j, l); da = DA(j, l);
0295 a += aa* px; d += da* px;
0296 b += aa* py; e += da* py;
0297 }
0298 Pt2d.SetCoord(a, b);
0299 V12d.SetCoord(d, e);
0300 i2 += nbpoles;
0301
0302 MyV2d = gp_Vec2d(Pt2d, TabP2d(k));
0303 FU += MyV2d*V12d;
0304 DFU += V12d.SquareMagnitude();
0305 }
0306
0307 if (DFU >= RealEpsilon()) {
0308 DU = FU/DFU;
0309 DU = Sign(Min(5.e-02, Abs(DU)), DU);
0310 UF += DU;
0311 Parameters(j) = UF;
0312 }
0313 }
0314
0315 MyF.Value(Parameters, Fval);
0316 MError3d = MyF.MaxError3d();
0317 MError2d = MyF.MaxError2d();
0318 }
0319
0320
0321 if (MError3d<= Tol3d && MError2d <= Tol2d) {
0322 Done = Standard_True;
0323 }
0324 else if (NbIterations != 0) {
0325 // NbIterations de gradient conjugue:
0326 // =================================
0327 Standard_Real Eps = 1.e-07;
0328 AppParCurves_BSpGradient_BFGS FResol(MyF, Parameters, Tol3d,
0329 Tol2d, Eps, NbIterations);
0330 }
0331
0332 SCU = MyF.CurveValue();
0333
0334
0335 AvError = 0.;
0336 for (j = FirstPoint; j <= LastPoint; j++) {
0337 Parameters(j) = MyF.NewParameters()(j);
0338 // Recherche des erreurs maxi et moyenne a un index donne:
0339 for (k = 1; k <= nbP; k++) {
0340 ParError(j) = Max(ParError(j), MyF.Error(j, k));
0341 }
0342 AvError += ParError(j);
0343 }
0344 AvError = AvError/(LastPoint-FirstPoint+1);
0345
0346
0347 MError3d = MyF.MaxError3d();
0348 MError2d = MyF.MaxError2d();
0349 if (MError3d <= Tol3d && MError2d <= Tol2d) {
0350 Done = Standard_True;
0351 }
0352
0353 }
0354
0355
0356
0357 AppParCurves_MultiBSpCurve AppParCurves_BSpGradient::Value() const {
0358 return SCU;
0359 }
0360
0361
0362 Standard_Boolean AppParCurves_BSpGradient::IsDone() const {
0363 return Done;
0364 }
0365
0366
0367 Standard_Real AppParCurves_BSpGradient::Error(const Standard_Integer Index) const {
0368 return ParError(Index);
0369 }
0370
0371 Standard_Real AppParCurves_BSpGradient::AverageError() const {
0372 return AvError;
0373 }
0374
0375 Standard_Real AppParCurves_BSpGradient::MaxError3d() const {
0376 return MError3d;
0377 }
0378
0379 Standard_Real AppParCurves_BSpGradient::MaxError2d() const {
0380 return MError2d;
0381 }
0382
0383
