0023024: Update headers of OCCT files
[occt.git] / src / AppParCurves / AppParCurves_Gradient.gxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
3 //
4 // The content of this file is subject to the Open CASCADE Technology Public
5 // License Version 6.5 (the "License"). You may not use the content of this file
6 // except in compliance with the License. Please obtain a copy of the License
7 // at http://www.opencascade.org and read it completely before using this file.
8 //
9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
11 //
12 // The Original Code and all software distributed under the License is
13 // distributed on an "AS IS" basis, without warranty of any kind, and the
14 // Initial Developer hereby disclaims all such warranties, including without
15 // limitation, any warranties of merchantability, fitness for a particular
16 // purpose or non-infringement. Please see the License for the specific terms
17 // and conditions governing the rights and limitations under the License.
18
19 // lpa, le 11/09/91
20
21
22 // Application de la methode du gradient corrige pour minimiser 
23 // F  = somme(||C(ui, Poles(ui)) - ptli||2.
24 // La methode de gradient conjugue est programmee dans la bibliotheque 
25 // mathematique: math_BFGS.
26 // cet algorithme doit etre appele uniquement lorsque on a affaire a un set 
27 // de points contraints (ailleurs qu aux extremites). En effet, l appel de la 
28 // fonction F a minimiser implique un appel a ParLeastSquare et ResConstraint.
29 // Si ce n est pas le cas, l appel a ResConstraint est equivalent a une 
30 // seconde resolution par les moindres carres donc beaucoup de temps perdu.
31
32
33 #define No_Standard_RangeError
34 #define No_Standard_OutOfRange
35
36 #include <AppParCurves_Constraint.hxx>
37 #include <math_BFGS.hxx>
38 #include <StdFail_NotDone.hxx>
39 #include <AppParCurves_MultiPoint.hxx>
40 #include <gp_Pnt.hxx>
41 #include <gp_Pnt2d.hxx>
42 #include <gp_Vec.hxx>
43 #include <gp_Vec2d.hxx>
44 #include <TColgp_Array1OfPnt.hxx>
45 #include <TColgp_Array1OfPnt2d.hxx>
46 #include <TColgp_Array1OfVec.hxx>
47 #include <TColgp_Array1OfVec2d.hxx>
48 #include <BSplCLib.hxx>
49 #include <PLib.hxx>
50
51 // #define AppParCurves_Gradient_BFGS BFGS_/**/AppParCurves_Gradient
52
53
54
55 AppParCurves_Gradient::
56    AppParCurves_Gradient(const MultiLine& SSP,
57          const Standard_Integer FirstPoint,
58          const Standard_Integer LastPoint,
59          const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
60          math_Vector& Parameters,
61          const Standard_Integer Deg,
62          const Standard_Real Tol3d,
63          const Standard_Real Tol2d,
64          const Standard_Integer NbIterations):
65          ParError(FirstPoint, LastPoint,0.0) {
66
67 //  Standard_Boolean grad = Standard_True;
68   Standard_Integer j, k, i2, l;
69   Standard_Real UF, DU, Fval = 0.0, FU, DFU;
70   Standard_Integer nbP3d = ToolLine::NbP3d(SSP);
71   Standard_Integer nbP2d = ToolLine::NbP2d(SSP);
72   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
73   Standard_Integer nbP = nbP3d + nbP2d;
74 //  gp_Pnt Pt, P1, P2;
75   gp_Pnt Pt;
76 //  gp_Pnt2d Pt2d, P12d, P22d;
77   gp_Pnt2d Pt2d;
78 //  gp_Vec V1, V2, MyV;
79   gp_Vec V1, MyV;
80 //  gp_Vec2d V12d, V22d, MyV2d;
81   gp_Vec2d V12d, MyV2d;
82   Done = Standard_False;
83   
84   if (nbP3d == 0) mynbP3d = 1;
85   if (nbP2d == 0) mynbP2d = 1;
86   TColgp_Array1OfPnt TabP(1, mynbP3d);
87   TColgp_Array1OfPnt2d TabP2d(1, mynbP2d);
88   TColgp_Array1OfVec TabV(1, mynbP3d);
89   TColgp_Array1OfVec2d TabV2d(1, mynbP2d);
90
91   // Calcul de la fonction F= somme(||C(ui)-Ptli||2):
92   // Appel a une fonction heritant de MultipleVarFunctionWithGradient
93   // pour calculer F et grad_F.
94   // ================================================================
95
96   AppParCurves_ParFunction MyF(SSP, FirstPoint,LastPoint, TheConstraints, Parameters, Deg);
97
98
99   if (!MyF.Value(Parameters, Fval)) {
100     Done = Standard_False;
101     return;
102   }
103
104   SCU = MyF.CurveValue();
105   Standard_Integer deg = SCU.NbPoles()-1;
106   TColgp_Array1OfPnt   TabPole(1, deg+1),   TabCoef(1, deg+1);
107   TColgp_Array1OfPnt2d TabPole2d(1, deg+1), TabCoef2d(1, deg+1);
108   TColgp_Array1OfPnt    TheCoef(1, (deg+1)*mynbP3d);
109   TColgp_Array1OfPnt2d  TheCoef2d(1, (deg+1)*mynbP2d);
110
111   
112   // Stockage des Poles des courbes pour projeter:
113   // ============================================
114   i2 = 0;
115   for (k = 1; k <= nbP3d; k++) {
116     SCU.Curve(k, TabPole);
117     BSplCLib::PolesCoefficients(TabPole, PLib::NoWeights(),
118                                 TabCoef, PLib::NoWeights());
119     for (j=1; j<=deg+1; j++) TheCoef(j+i2) = TabCoef(j);
120     i2 += deg+1;
121   }
122   i2 = 0;
123   for (k = 1; k <= nbP2d; k++) {
124     SCU.Curve(nbP3d+k, TabPole2d);
125     BSplCLib::PolesCoefficients(TabPole2d, PLib::NoWeights(),
126                                 TabCoef2d, PLib::NoWeights());
127     for (j=1; j<=deg+1; j++) TheCoef2d(j+i2) = TabCoef2d(j);
128     i2 += deg+1;
129   }
130
131   //  Une iteration rapide de projection est faite par la methode de 
132   //  Rogers & Fog 89, methode equivalente a Hoschek 88 qui ne necessite pas
133   //  le calcul de D2.
134
135
136   // Iteration de Projection:
137   // =======================
138   for (j = FirstPoint+1; j <= LastPoint-1; j++) {
139     UF = Parameters(j);
140     if (nbP != 0 && nbP2d != 0) ToolLine::Value(SSP, j, TabP, TabP2d);
141     else if (nbP2d != 0)        ToolLine::Value(SSP, j, TabP2d);
142     else                        ToolLine::Value(SSP, j, TabP);
143     
144     FU = 0.0;
145     DFU = 0.0;
146     i2 = 0;
147     for (k = 1; k <= nbP3d; k++) {
148       for (l=1; l<=deg+1; l++) TabCoef(l) = TheCoef(l+i2);
149       i2 += deg+1;
150       BSplCLib::CoefsD1(UF, TabCoef, PLib::NoWeights(), Pt, V1);
151       MyV = gp_Vec(Pt, TabP(k));
152       FU += MyV*V1;
153       DFU += V1.SquareMagnitude();
154     }
155     i2 = 0;
156     for (k = 1; k <= nbP2d; k++) {
157       for (l=1; l<=deg+1; l++) TabCoef2d(l) = TheCoef2d(l+i2);
158       i2 += deg+1;
159       BSplCLib::CoefsD1(UF, TabCoef2d, PLib::NoWeights(), Pt2d, V12d);
160       MyV2d = gp_Vec2d(Pt2d, TabP2d(k));
161       FU += MyV2d*V12d;
162       DFU += V12d.SquareMagnitude();
163     }
164     
165     if (DFU >= RealEpsilon()) {
166       DU = FU/DFU;
167       DU = Sign(Min(5.e-02, Abs(DU)), DU);
168       UF += DU;
169       Parameters(j) = UF;
170     }
171   }
172
173   
174   if (!MyF.Value(Parameters, Fval)) {
175     SCU = AppParCurves_MultiCurve();
176     Done = Standard_False;
177     return;
178   }
179   MError3d = MyF.MaxError3d();
180   MError2d = MyF.MaxError2d();
181   
182   if (MError3d<= Tol3d && MError2d <= Tol2d) {
183     Done = Standard_True;
184     SCU = MyF.CurveValue();
185   }
186   else if (NbIterations != 0) {
187     // NbIterations de gradient conjugue:
188     // =================================
189     Standard_Real Eps = 1.e-07;
190     AppParCurves_Gradient_BFGS FResol(MyF, Parameters, Tol3d, Tol2d, Eps, NbIterations);
191     Parameters = MyF.NewParameters();
192     SCU = MyF.CurveValue();
193   }
194
195     
196   AvError = 0.;
197   for (j = FirstPoint; j <= LastPoint; j++) {  
198     // Recherche des erreurs maxi et moyenne a un index donne:
199     for (k = 1; k <= nbP; k++) {
200       ParError(j) = Max(ParError(j), MyF.Error(j, k));
201     }
202     AvError += ParError(j);
203   }
204   AvError = AvError/(LastPoint-FirstPoint+1);
205
206
207   MError3d = MyF.MaxError3d();
208   MError2d = MyF.MaxError2d();
209   if (MError3d <= Tol3d && MError2d <= Tol2d) {
210     Done = Standard_True;
211   }
212
213 }
214
215
216
217 AppParCurves_MultiCurve AppParCurves_Gradient::Value() const {
218   return SCU;
219 }
220
221
222 Standard_Boolean AppParCurves_Gradient::IsDone() const {
223   return Done;
224 }
225
226
227 Standard_Real AppParCurves_Gradient::Error(const Standard_Integer Index) const {
228   return ParError(Index);
229 }
230
231 Standard_Real AppParCurves_Gradient::AverageError() const {
232   return AvError;
233 }
234
235 Standard_Real AppParCurves_Gradient::MaxError3d() const {
236   return MError3d;
237 }
238
239 Standard_Real AppParCurves_Gradient::MaxError2d() const {
240   return MError2d;
241 }
242
243