0024624: Lost word in license statement in source files
[occt.git] / src / AppParCurves / AppParCurves_Gradient.gxx
CommitLineData
b311480e 1// Copyright (c) 1995-1999 Matra Datavision
973c2be1 2// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 3//
973c2be1 4// This file is part of Open CASCADE Technology software library.
b311480e 5//
d5f74e42 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
973c2be1 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.
b311480e 11//
973c2be1 12// Alternatively, this file may be used under the terms of Open CASCADE
13// commercial license or contractual agreement.
b311480e 14
7fd59977 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
51AppParCurves_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
213AppParCurves_MultiCurve AppParCurves_Gradient::Value() const {
214 return SCU;
215}
216
217
218Standard_Boolean AppParCurves_Gradient::IsDone() const {
219 return Done;
220}
221
222
223Standard_Real AppParCurves_Gradient::Error(const Standard_Integer Index) const {
224 return ParError(Index);
225}
226
227Standard_Real AppParCurves_Gradient::AverageError() const {
228 return AvError;
229}
230
231Standard_Real AppParCurves_Gradient::MaxError3d() const {
232 return MError3d;
233}
234
235Standard_Real AppParCurves_Gradient::MaxError2d() const {
236 return MError2d;
237}
238
239