1 // Created on: 1997-09-17
2 // Created by: Philippe MANGIN /Igor Feoktistov (1998)
3 // Copyright (c) 1997-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
23 #include <SortTools_StraightInsertionSortOfReal.hxx>
24 #include <TCollection_CompareOfReal.hxx>
25 #include <TColStd_HArray2OfInteger.hxx>
26 #include <TColStd_HArray1OfReal.hxx>
27 #include <TColStd_Array2OfInteger.hxx>
28 #include <FEmTool_Assembly.hxx>
29 #include <FEmTool_AssemblyTable.hxx>
30 #include <FEmTool_Curve.hxx>
32 //====================== Private Methodes =============================//
33 //=======================================================================
35 //purpose : Smoothing's motor like STRIM routine "MOTLIS"
36 //=======================================================================
37 void AppParCurves_Variational::TheMotor(
38 Handle(AppParCurves_SmoothCriterion)& J,
39 // const Standard_Real WQuadratic,
41 const Standard_Real WQuality,
42 Handle(FEmTool_Curve)& TheCurve,
43 TColStd_Array1OfReal& Ecarts)
47 const Standard_Real BigValue = 1.e37, SmallValue = 1.e-6, SmallestValue = 1.e-9;
49 // SortTools_StraightInsertionSortOfReal Sort;
50 TCollection_CompareOfReal CompReal;
51 Handle(TColStd_HArray1OfReal) CurrentTi, NewTi, OldTi;
52 Handle(TColStd_HArray2OfInteger) Dependence;
53 Standard_Boolean lestim, lconst, ToOptim, iscut;
54 Standard_Boolean isnear, again = Standard_True;
55 Standard_Integer NbEst, ICDANA, NumPnt, Iter;
56 Standard_Integer MaxNbEst =5;
57 Standard_Real VOCRI[3] = {BigValue, BigValue, BigValue}, EROLD = BigValue,
58 VALCRI[3], ERRMAX, ERRMOY, ERRQUA;
59 Standard_Real CBLONG, LNOLD;
60 Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
61 Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
62 Handle(FEmTool_Curve) CCurrent, COld, CNew;
63 Standard_Real EpsLength;// = 1.e-6 * CBLONG / NbrPnt;
64 Standard_Real EpsDeg; // = Min(WQuality * .1, CBLONG * .001);
66 Standard_Real e1, e2, e3;
67 Standard_Real J1min, J2min, J3min;
68 Standard_Integer iprog;
72 J->GetEstimation(e1, e2, e3);
73 J1min = 1.e-8; J2min = J3min = (e1 + 1.e-8) * 1.e-6;
75 if(e1 < J1min) e1 = J1min;// Like in
76 if(e2 < J2min) e2 = J2min;// MOTLIS
77 if(e3 < J3min) e3 = J3min;
80 J->SetEstimation(e1, e2, e3);
83 CurrentTi = new TColStd_HArray1OfReal(1, myParameters->Length());
84 CurrentTi->ChangeArray1() = myParameters->Array1();
85 OldTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
86 OldTi->ChangeArray1() = CurrentTi->Array1();
88 LNOLD = CBLONG = J->EstLength();
89 Dependence = J->DependenceTable();
91 J->SetCurve(CCurrent);
92 FEmTool_Assembly * TheAssembly =
93 new FEmTool_Assembly (Dependence->Array2(), J->AssemblyTable());
95 //============ Optimization ============================
96 // Standard_Integer inagain = 0;
99 // (1) Loop Optimization / Estimation
100 lestim = Standard_True;
101 lconst = Standard_True;
104 J->SetCurve(CCurrent);
108 // (1.1) Curve's Optimization.
109 EpsLength = SmallValue * CBLONG / NbrPnt;
110 CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
111 CCurrent->NbElements(), CCurrent->Base(), EpsLength);
112 CNew->Knots() = CCurrent->Knots();
114 J->SetParameters(CurrentTi);
115 EpsDeg = Min(WQuality * .1, CBLONG * .001);
117 Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
119 lconst = Standard_False;
121 // (1.2) calcul des criteres de qualites et amelioration
123 ICDANA = J->QualityValues(J1min, J2min, J3min,
124 VALCRI[0], VALCRI[1], VALCRI[2]);
126 if(ICDANA > 0) lconst = Standard_True;
128 J->ErrorValues(ERRMAX, ERRQUA, ERRMOY);
130 isnear = ((Sqrt(ERRQUA / NbrPnt) < 2*WQuality) &&
131 (myNbIterations > 1));
133 // (1.3) Optimisation des ti par proj orthogonale
134 // et calcul de l'erreur aux points.
137 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
138 Project(CNew, CurrentTi->Array1(),
139 NewTi->ChangeArray1(),
141 ERRMAX, ERRQUA, ERRMOY, 2);
143 else NewTi = CurrentTi;
146 // (1.4) Progression's test
148 if ((EROLD > WQuality) && (ERRMAX < 0.95*EROLD)) iprog++;
149 if ((EROLD > WQuality) && (ERRMAX < 0.8*EROLD)) iprog++;
150 if ((EROLD > WQuality) && (ERRMAX < WQuality)) iprog++;
151 if ((EROLD > WQuality) && (ERRMAX < 0.99*EROLD)
152 && (ERRMAX < 1.1*WQuality)) iprog++;
153 if ( VALCRI[0] < 0.975 * VOCRI[0]) iprog++;
154 if ( VALCRI[0] < 0.9 * VOCRI[0]) iprog++;
155 if ( VALCRI[1] < 0.95 * VOCRI[1]) iprog++;
156 if ( VALCRI[1] < 0.8 * VOCRI[1]) iprog++;
157 if ( VALCRI[2] < 0.95 * VOCRI[2]) iprog++;
158 if ( VALCRI[2] < 0.8 * VOCRI[2]) iprog++;
159 if ((VOCRI[1]>SmallestValue)&&(VOCRI[2]>SmallestValue)) {
160 if ((VALCRI[1]/VOCRI[1] + 2*VALCRI[2]/VOCRI[2]) < 2.8) iprog++;
163 if (iprog < 2 && NbEst == 0 ) {
164 // (1.5) Invalidation of new knots.
165 VALCRI[0] = VOCRI[0];
166 VALCRI[1] = VOCRI[1];
167 VALCRI[2] = VOCRI[2];
176 VOCRI[0] = VALCRI[0];
177 VOCRI[1] = VALCRI[1];
178 VOCRI[2] = VALCRI[2];
185 // (1.6) Test if the Estimations seems OK, else repeat
187 lestim = ( (NbEst<MaxNbEst) && (ICDANA == 2)&& (iprog > 0) );
189 if (lestim && isnear) {
190 // (1.7) Optimization of ti by ACR.
191 // Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
192 SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
193 Standard_Integer Decima = 4;
195 CCurrent->Length(0., 1., CBLONG);
196 J->EstLength() = CBLONG;
198 ACR(CCurrent, CurrentTi->ChangeArray1(), Decima);
199 lconst = Standard_True;
205 // (2) loop of parametric / geometric optimization
208 ToOptim = ((Iter < myNbIterations) && (isnear));
212 // (2.1) Save curent results
213 VOCRI[0] = VALCRI[0];
214 VOCRI[1] = VALCRI[1];
215 VOCRI[2] = VALCRI[2];
219 OldTi->ChangeArray1() = CurrentTi->Array1();
221 // (2.2) Optimization des ti by ACR.
222 // Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
223 SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
224 Standard_Integer Decima = 4;
226 CCurrent->Length(0., 1., CBLONG);
227 J->EstLength() = CBLONG;
229 ACR(CCurrent, CurrentTi->ChangeArray1(), Decima);
230 lconst = Standard_True;
232 // (2.3) Optimisation des courbes
233 EpsLength = SmallValue * CBLONG / NbrPnt;
235 CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
236 CCurrent->NbElements(), CCurrent->Base(), EpsLength);
237 CNew->Knots() = CCurrent->Knots();
239 J->SetParameters(CurrentTi);
241 EpsDeg = Min(WQuality * .1, CBLONG * .001);
242 Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
246 // (2.4) calcul des criteres de qualites et amelioration
249 ICDANA = J->QualityValues(J1min, J2min, J3min, VALCRI[0], VALCRI[1], VALCRI[2]);
250 if(ICDANA > 0) lconst = Standard_True;
252 J->GetEstimation(e1, e2, e3);
253 // (2.5) Optimisation des ti par proj orthogonale
255 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
256 Project(CCurrent, CurrentTi->Array1(),
257 NewTi->ChangeArray1(),
259 ERRMAX, ERRQUA, ERRMOY, 2);
261 // (2.6) Test de non regression
263 Standard_Integer iregre = 0;
264 if (NbrConstraint < NbrPnt) {
265 if ( (ERRMAX > WQuality) && (ERRMAX > 1.05*EROLD)) iregre++;
266 if ( (ERRMAX > WQuality) && (ERRMAX > 2*EROLD)) iregre++;
267 if ( (EROLD > WQuality) && (ERRMAX <= 0.5*EROLD)) iregre--;
269 Standard_Real E1, E2, E3;
270 J->GetEstimation(E1, E2, E3);
271 if ( (VALCRI[0] > E1) && (VALCRI[0] > 1.1*VOCRI[0])) iregre++;
272 if ( (VALCRI[1] > E2) && (VALCRI[1] > 1.1*VOCRI[1])) iregre++;
273 if ( (VALCRI[2] > E3) && (VALCRI[2] > 1.1*VOCRI[2])) iregre++;
277 // if (iregre >= 1) {
278 // (2.7) on restaure l'iteration precedente
279 VALCRI[0] = VOCRI[0];
280 VALCRI[1] = VOCRI[1];
281 VALCRI[2] = VOCRI[2];
285 CurrentTi->ChangeArray1() = OldTi->Array1();
286 ToOptim = Standard_False;
288 else { // Iteration is Ok.
292 if (Iter >= myNbIterations) ToOptim = Standard_False;
295 // (3) Decoupe eventuelle
297 if ( (CCurrent->NbElements() < myMaxSegment) && myWithCutting ) {
299 // (3.1) Sauvgarde de l'etat precedent
300 VOCRI[0] = VALCRI[0];
301 VOCRI[1] = VALCRI[1];
302 VOCRI[2] = VALCRI[2];
305 OldTi->ChangeArray1() = CurrentTi->Array1();
307 // (3.2) On arrange les ti : Trie + recadrage sur (0,1)
308 // ---> On trie, afin d'assurer l'ordre par la suite.
309 // Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
310 SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
312 if ((CurrentTi->Value(1)!= 0.) ||
313 (CurrentTi->Value(NbrPnt)!= 1.)) {
314 Standard_Real t, DelatT =
315 1.0 /(CurrentTi->Value(NbrPnt)-CurrentTi->Value(1));
316 for (Standard_Integer ii=2; ii<NbrPnt; ii++) {
317 t = (CurrentTi->Value(ii)-CurrentTi->Value(1))*DelatT;
318 CurrentTi->SetValue(ii, t);
320 CurrentTi->SetValue(1, 0.);
321 CurrentTi->SetValue(NbrPnt, 1.);
324 // (3.3) Insert new Knots
326 SplitCurve(CCurrent, CurrentTi->Array1(), EpsLength, CNew, iscut);
327 if (!iscut) again = Standard_False;
330 // New Knots => New Assembly.
333 TheAssembly = new FEmTool_Assembly (Dependence->Array2(),
337 else { again = Standard_False;}
340 // ================ Great loop end ===================
344 // (4) Compute the best Error.
345 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
346 Project(CCurrent, CurrentTi->Array1(),
347 NewTi->ChangeArray1(),
349 ERRMAX, ERRQUA, ERRMOY, 10);
351 // (5) field's update
354 J->EstLength() = CBLONG;
355 myParameters->ChangeArray1() = NewTi->Array1();
356 myCriterium[0] = ERRQUA;
357 myCriterium[1] = Sqrt(VALCRI[0]);
358 myCriterium[2] = Sqrt(VALCRI[1]);
359 myCriterium[3] = Sqrt(VALCRI[2]);
361 myMaxErrorIndex = NumPnt;
362 if(NbrPnt > NbrConstraint)
363 myAverageError = ERRMOY / (NbrPnt - NbrConstraint);
365 myAverageError = ERRMOY / NbrConstraint;