0024624: Lost word in license statement in source files
[occt.git] / src / AppParCurves / AppParCurves_Variational_1.gxx
CommitLineData
b311480e 1// Created on: 1997-09-17
2// Created by: Philippe MANGIN /Igor Feoktistov (1998)
3// Copyright (c) 1997-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 10// by the Free Software Foundation, with special exception defined in the file
11// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12// distribution for complete text of the license and disclaimer of any warranty.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
17#include <SortTools_StraightInsertionSortOfReal.hxx>
18#include <TCollection_CompareOfReal.hxx>
19#include <TColStd_HArray2OfInteger.hxx>
20#include <TColStd_HArray1OfReal.hxx>
21#include <TColStd_Array2OfInteger.hxx>
22#include <FEmTool_Assembly.hxx>
23#include <FEmTool_AssemblyTable.hxx>
24#include <FEmTool_Curve.hxx>
25
26//====================== Private Methodes =============================//
27//=======================================================================
28//function : TheMotor
29//purpose : Smoothing's motor like STRIM routine "MOTLIS"
30//=======================================================================
7fd59977 31void AppParCurves_Variational::TheMotor(
32 Handle(AppParCurves_SmoothCriterion)& J,
33// const Standard_Real WQuadratic,
34 const Standard_Real ,
35 const Standard_Real WQuality,
36 Handle(FEmTool_Curve)& TheCurve,
37 TColStd_Array1OfReal& Ecarts)
38{
39// ...
40
41 const Standard_Real BigValue = 1.e37, SmallValue = 1.e-6, SmallestValue = 1.e-9;
42
43// SortTools_StraightInsertionSortOfReal Sort;
44 TCollection_CompareOfReal CompReal;
45 Handle(TColStd_HArray1OfReal) CurrentTi, NewTi, OldTi;
46 Handle(TColStd_HArray2OfInteger) Dependence;
47 Standard_Boolean lestim, lconst, ToOptim, iscut;
1d47d8d0 48 Standard_Boolean isnear = Standard_False, again = Standard_True;
7fd59977 49 Standard_Integer NbEst, ICDANA, NumPnt, Iter;
50 Standard_Integer MaxNbEst =5;
51 Standard_Real VOCRI[3] = {BigValue, BigValue, BigValue}, EROLD = BigValue,
1d47d8d0 52 VALCRI[3], ERRMAX = BigValue, ERRMOY, ERRQUA;
7fd59977 53 Standard_Real CBLONG, LNOLD;
54 Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
55 Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
56 Handle(FEmTool_Curve) CCurrent, COld, CNew;
1d47d8d0 57 Standard_Real EpsLength = SmallValue;
58 Standard_Real EpsDeg;
7fd59977 59
60 Standard_Real e1, e2, e3;
61 Standard_Real J1min, J2min, J3min;
62 Standard_Integer iprog;
63
64// (0) Init
65
66 J->GetEstimation(e1, e2, e3);
67 J1min = 1.e-8; J2min = J3min = (e1 + 1.e-8) * 1.e-6;
68
69 if(e1 < J1min) e1 = J1min;// Like in
70 if(e2 < J2min) e2 = J2min;// MOTLIS
71 if(e3 < J3min) e3 = J3min;
72
73
74 J->SetEstimation(e1, e2, e3);
75
76 CCurrent = TheCurve;
77 CurrentTi = new TColStd_HArray1OfReal(1, myParameters->Length());
78 CurrentTi->ChangeArray1() = myParameters->Array1();
79 OldTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
80 OldTi->ChangeArray1() = CurrentTi->Array1();
81 COld = CCurrent;
82 LNOLD = CBLONG = J->EstLength();
83 Dependence = J->DependenceTable();
84
85 J->SetCurve(CCurrent);
86 FEmTool_Assembly * TheAssembly =
87 new FEmTool_Assembly (Dependence->Array2(), J->AssemblyTable());
88
89//============ Optimization ============================
90// Standard_Integer inagain = 0;
91 while (again) {
92
93 // (1) Loop Optimization / Estimation
94 lestim = Standard_True;
95 lconst = Standard_True;
96 NbEst = 0;
97
98 J->SetCurve(CCurrent);
99
100 while(lestim) {
101
102 // (1.1) Curve's Optimization.
103 EpsLength = SmallValue * CBLONG / NbrPnt;
104 CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
105 CCurrent->NbElements(), CCurrent->Base(), EpsLength);
106 CNew->Knots() = CCurrent->Knots();
107
108 J->SetParameters(CurrentTi);
109 EpsDeg = Min(WQuality * .1, CBLONG * .001);
110
111 Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
112
113 lconst = Standard_False;
114
115 // (1.2) calcul des criteres de qualites et amelioration
116 // des estimation.
117 ICDANA = J->QualityValues(J1min, J2min, J3min,
118 VALCRI[0], VALCRI[1], VALCRI[2]);
119
120 if(ICDANA > 0) lconst = Standard_True;
121
122 J->ErrorValues(ERRMAX, ERRQUA, ERRMOY);
123
124 isnear = ((Sqrt(ERRQUA / NbrPnt) < 2*WQuality) &&
125 (myNbIterations > 1));
126
127// (1.3) Optimisation des ti par proj orthogonale
128// et calcul de l'erreur aux points.
129
130 if (isnear) {
131 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
132 Project(CNew, CurrentTi->Array1(),
133 NewTi->ChangeArray1(),
134 Ecarts, NumPnt,
135 ERRMAX, ERRQUA, ERRMOY, 2);
136 }
137 else NewTi = CurrentTi;
138
139
140// (1.4) Progression's test
141 iprog = 0;
142 if ((EROLD > WQuality) && (ERRMAX < 0.95*EROLD)) iprog++;
143 if ((EROLD > WQuality) && (ERRMAX < 0.8*EROLD)) iprog++;
144 if ((EROLD > WQuality) && (ERRMAX < WQuality)) iprog++;
145 if ((EROLD > WQuality) && (ERRMAX < 0.99*EROLD)
146 && (ERRMAX < 1.1*WQuality)) iprog++;
147 if ( VALCRI[0] < 0.975 * VOCRI[0]) iprog++;
148 if ( VALCRI[0] < 0.9 * VOCRI[0]) iprog++;
149 if ( VALCRI[1] < 0.95 * VOCRI[1]) iprog++;
150 if ( VALCRI[1] < 0.8 * VOCRI[1]) iprog++;
151 if ( VALCRI[2] < 0.95 * VOCRI[2]) iprog++;
152 if ( VALCRI[2] < 0.8 * VOCRI[2]) iprog++;
153 if ((VOCRI[1]>SmallestValue)&&(VOCRI[2]>SmallestValue)) {
154 if ((VALCRI[1]/VOCRI[1] + 2*VALCRI[2]/VOCRI[2]) < 2.8) iprog++;
155 }
156
157 if (iprog < 2 && NbEst == 0 ) {
158 // (1.5) Invalidation of new knots.
159 VALCRI[0] = VOCRI[0];
160 VALCRI[1] = VOCRI[1];
161 VALCRI[2] = VOCRI[2];
162 ERRMAX = EROLD;
163 CBLONG = LNOLD;
164 CCurrent = COld;
165 CurrentTi = OldTi;
166
167 goto L8000; // exit
168 }
169
170 VOCRI[0] = VALCRI[0];
171 VOCRI[1] = VALCRI[1];
172 VOCRI[2] = VALCRI[2];
173 LNOLD = CBLONG;
174 EROLD = ERRMAX;
175
176 CCurrent = CNew;
177 CurrentTi = NewTi;
178
179 // (1.6) Test if the Estimations seems OK, else repeat
180 NbEst++;
181 lestim = ( (NbEst<MaxNbEst) && (ICDANA == 2)&& (iprog > 0) );
182
183 if (lestim && isnear) {
184 // (1.7) Optimization of ti by ACR.
185// Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
186 SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
187 Standard_Integer Decima = 4;
188
189 CCurrent->Length(0., 1., CBLONG);
190 J->EstLength() = CBLONG;
191
192 ACR(CCurrent, CurrentTi->ChangeArray1(), Decima);
193 lconst = Standard_True;
194
195 }
196 }
197
198
199 // (2) loop of parametric / geometric optimization
200
201 Iter = 1;
202 ToOptim = ((Iter < myNbIterations) && (isnear));
203
204 while(ToOptim) {
205 Iter++;
206 // (2.1) Save curent results
207 VOCRI[0] = VALCRI[0];
208 VOCRI[1] = VALCRI[1];
209 VOCRI[2] = VALCRI[2];
210 EROLD = ERRMAX;
211 LNOLD = CBLONG;
212 COld = CCurrent;
213 OldTi->ChangeArray1() = CurrentTi->Array1();
214
215 // (2.2) Optimization des ti by ACR.
216// Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
217 SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
218 Standard_Integer Decima = 4;
219
220 CCurrent->Length(0., 1., CBLONG);
221 J->EstLength() = CBLONG;
222
223 ACR(CCurrent, CurrentTi->ChangeArray1(), Decima);
224 lconst = Standard_True;
225
226 // (2.3) Optimisation des courbes
227 EpsLength = SmallValue * CBLONG / NbrPnt;
228
229 CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
230 CCurrent->NbElements(), CCurrent->Base(), EpsLength);
231 CNew->Knots() = CCurrent->Knots();
232
233 J->SetParameters(CurrentTi);
234
235 EpsDeg = Min(WQuality * .1, CBLONG * .001);
236 Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
237
238 CCurrent = CNew;
239
240 // (2.4) calcul des criteres de qualites et amelioration
241 // des estimation.
242
243 ICDANA = J->QualityValues(J1min, J2min, J3min, VALCRI[0], VALCRI[1], VALCRI[2]);
244 if(ICDANA > 0) lconst = Standard_True;
245
246 J->GetEstimation(e1, e2, e3);
247 // (2.5) Optimisation des ti par proj orthogonale
248
249 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
250 Project(CCurrent, CurrentTi->Array1(),
251 NewTi->ChangeArray1(),
252 Ecarts, NumPnt,
253 ERRMAX, ERRQUA, ERRMOY, 2);
254
255 // (2.6) Test de non regression
256
257 Standard_Integer iregre = 0;
258 if (NbrConstraint < NbrPnt) {
259 if ( (ERRMAX > WQuality) && (ERRMAX > 1.05*EROLD)) iregre++;
260 if ( (ERRMAX > WQuality) && (ERRMAX > 2*EROLD)) iregre++;
261 if ( (EROLD > WQuality) && (ERRMAX <= 0.5*EROLD)) iregre--;
262 }
263 Standard_Real E1, E2, E3;
264 J->GetEstimation(E1, E2, E3);
265 if ( (VALCRI[0] > E1) && (VALCRI[0] > 1.1*VOCRI[0])) iregre++;
266 if ( (VALCRI[1] > E2) && (VALCRI[1] > 1.1*VOCRI[1])) iregre++;
267 if ( (VALCRI[2] > E3) && (VALCRI[2] > 1.1*VOCRI[2])) iregre++;
268
269
270 if (iregre >= 2) {
271// if (iregre >= 1) {
272 // (2.7) on restaure l'iteration precedente
273 VALCRI[0] = VOCRI[0];
274 VALCRI[1] = VOCRI[1];
275 VALCRI[2] = VOCRI[2];
276 ERRMAX = EROLD;
277 CBLONG = LNOLD;
278 CCurrent = COld;
279 CurrentTi->ChangeArray1() = OldTi->Array1();
280 ToOptim = Standard_False;
281 }
282 else { // Iteration is Ok.
283 CCurrent = CNew;
284 CurrentTi = NewTi;
285 }
286 if (Iter >= myNbIterations) ToOptim = Standard_False;
287 }
288
289 // (3) Decoupe eventuelle
290
291 if ( (CCurrent->NbElements() < myMaxSegment) && myWithCutting ) {
292
293 // (3.1) Sauvgarde de l'etat precedent
294 VOCRI[0] = VALCRI[0];
295 VOCRI[1] = VALCRI[1];
296 VOCRI[2] = VALCRI[2];
297 EROLD = ERRMAX;
298 COld = CCurrent;
299 OldTi->ChangeArray1() = CurrentTi->Array1();
300
301 // (3.2) On arrange les ti : Trie + recadrage sur (0,1)
302 // ---> On trie, afin d'assurer l'ordre par la suite.
303// Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
304 SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
305
306 if ((CurrentTi->Value(1)!= 0.) ||
307 (CurrentTi->Value(NbrPnt)!= 1.)) {
308 Standard_Real t, DelatT =
309 1.0 /(CurrentTi->Value(NbrPnt)-CurrentTi->Value(1));
310 for (Standard_Integer ii=2; ii<NbrPnt; ii++) {
311 t = (CurrentTi->Value(ii)-CurrentTi->Value(1))*DelatT;
312 CurrentTi->SetValue(ii, t);
313 }
314 CurrentTi->SetValue(1, 0.);
315 CurrentTi->SetValue(NbrPnt, 1.);
316 }
317
318 // (3.3) Insert new Knots
319
320 SplitCurve(CCurrent, CurrentTi->Array1(), EpsLength, CNew, iscut);
321 if (!iscut) again = Standard_False;
322 else {
323 CCurrent = CNew;
324 // New Knots => New Assembly.
325 J->SetCurve(CNew);
326 delete TheAssembly;
327 TheAssembly = new FEmTool_Assembly (Dependence->Array2(),
328 J->AssemblyTable());
329 }
330 }
331 else { again = Standard_False;}
332 }
333
334 // ================ Great loop end ===================
335
336 L8000:
337
338 // (4) Compute the best Error.
339 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
340 Project(CCurrent, CurrentTi->Array1(),
341 NewTi->ChangeArray1(),
342 Ecarts, NumPnt,
343 ERRMAX, ERRQUA, ERRMOY, 10);
344
345// (5) field's update
346
347 TheCurve = CCurrent;
348 J->EstLength() = CBLONG;
349 myParameters->ChangeArray1() = NewTi->Array1();
350 myCriterium[0] = ERRQUA;
351 myCriterium[1] = Sqrt(VALCRI[0]);
352 myCriterium[2] = Sqrt(VALCRI[1]);
353 myCriterium[3] = Sqrt(VALCRI[2]);
354 myMaxError = ERRMAX;
355 myMaxErrorIndex = NumPnt;
356 if(NbrPnt > NbrConstraint)
357 myAverageError = ERRMOY / (NbrPnt - NbrConstraint);
358 else
359 myAverageError = ERRMOY / NbrConstraint;
360
361 delete TheAssembly;
362}