1 // Created on: 1998-12-21
2 // Created by: Igor FEOKTISTOV
3 // Copyright (c) 1998-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
8 // This library is free software; you can redistribute it and / or modify it
9 // under the terms of the GNU Lesser General Public version 2.1 as published
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559
19 #include <PLib_Base.hxx>
20 #include <PLib_JacobiPolynomial.hxx>
21 #include <PLib_HermitJacobi.hxx>
22 #include <FEmTool_Curve.hxx>
23 #include <math_Vector.hxx>
24 #include <TColStd_Array1OfReal.hxx>
26 //=======================================================================
27 //function : InitSmoothCriterion
28 //purpose : Initializes the SmoothCriterion
29 //=======================================================================
30 void AppParCurves_Variational::InitSmoothCriterion()
33 const Standard_Real Eps2 = 1.e-6, Eps3 = 1.e-9;
34 // const Standard_Real J1 = .01, J2 = .001, J3 = .001;
40 InitParameters(Length);
42 mySmoothCriterion->SetParameters(myParameters);
44 Standard_Real E1, E2, E3;
46 InitCriterionEstimations(Length, E1, E2, E3);
48 J1 = 1.e-8; J2 = J3 = (E1 + 1.e-8) * 1.e-6;
54 mySmoothCriterion->EstLength() = Length;
55 mySmoothCriterion->SetEstimation(E1, E2, E3);
57 Standard_Real WQuadratic, WQuality;
59 if(!myWithMinMax && myTolerance != 0.)
60 WQuality = myTolerance;
61 else if (myTolerance == 0.)
64 WQuality = Max(myTolerance, Eps2 * Length);
66 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
67 WQuadratic = Sqrt((Standard_Real)(myNbPoints - NbConstr)) * WQuality;
68 if(WQuadratic > Eps3) WQuadratic = 1./ WQuadratic;
70 if(WQuadratic == 0.) WQuadratic = Max(Sqrt(E1), 1.);
73 mySmoothCriterion->SetWeight(WQuadratic, WQuality,
74 myPercent[0], myPercent[1], myPercent[2]);
77 Handle(PLib_Base) TheBase = new PLib_HermitJacobi(myMaxDegree, myContinuity);
78 Handle(FEmTool_Curve) TheCurve;
79 Standard_Integer NbElem;
80 Standard_Real CurvTol = Eps2 * Length / myNbPoints;
82 // Decoupe de l'intervalle en fonction des contraintes
83 if(myWithCutting == Standard_True && NbConstr != 0) {
85 InitCutting(TheBase, CurvTol, TheCurve);
91 TheCurve = new FEmTool_Curve(myDimension, NbElem, TheBase, CurvTol);
92 TheCurve->Knots().SetValue(TheCurve->Knots().Lower(),
93 myParameters->Value(myFirstPoint));
94 TheCurve->Knots().SetValue(TheCurve->Knots().Upper(),
95 myParameters->Value(myLastPoint));
99 mySmoothCriterion->SetCurve(TheCurve);
106 //=======================================================================
107 //function : InitParameters
108 //purpose : Calculation of initial estimation of the multicurve's length
109 // and parameters for multipoints.
110 //=======================================================================
112 void AppParCurves_Variational::InitParameters(Standard_Real& Length)
115 const Standard_Real Eps1 = Precision::Confusion() * .01;
117 Standard_Real aux, dist;
118 Standard_Integer i, i0, i1 = 0, ipoint;
122 myParameters->SetValue(myFirstPoint, Length);
124 for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint; ipoint++) {
125 i0 = i1; i1 += myDimension;
127 for(i = 1; i <= myDimension; i++) {
128 aux = myTabPoints->Value(i1 + i) - myTabPoints->Value(i0 + i);
131 Length += Sqrt(dist);
132 myParameters->SetValue(ipoint, Length);
137 Standard_ConstructionError::Raise("AppParCurves_Variational::InitParameters");
140 for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint - 1; ipoint++)
141 myParameters->SetValue(ipoint, myParameters->Value(ipoint) / Length);
143 myParameters->SetValue(myLastPoint, 1.);
145 // Avec peu de point il y a sans doute sous estimation ...
146 if(myNbPoints < 10) Length *= (1. + 0.1 / (myNbPoints - 1));
150 //=======================================================================
151 //function : InitCriterionEstimations
152 //purpose : Calculation of initial estimation of the linear criteria
154 //=======================================================================
156 void AppParCurves_Variational::InitCriterionEstimations(const Standard_Real Length,
159 Standard_Real& E3) const
161 E1 = Length * Length;
163 const Standard_Real Eps1 = Precision::Confusion() * .01;
165 math_Vector VTang1(1, myDimension), VTang2(1, myDimension), VTang3(1, myDimension),
166 VScnd1(1, myDimension), VScnd2(1, myDimension), VScnd3(1, myDimension);
168 // ========== Treatment of first point =================
170 Standard_Integer ipnt = myFirstPoint;
172 EstTangent(ipnt, VTang1);
174 EstTangent(ipnt, VTang2);
176 EstTangent(ipnt, VTang3);
179 EstSecnd(ipnt, VTang1, VTang2, Length, VScnd1);
181 EstSecnd(ipnt, VTang1, VTang3, Length, VScnd2);
183 // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 Begin
184 // Standard_Real Delta = .5 * (myParameters->Value(ipnt) - myParameters->Value(--ipnt));
185 Standard_Integer anInd = ipnt;
186 Standard_Real Delta = .5 * (myParameters->Value(anInd) - myParameters->Value(--ipnt));
187 // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 End
189 if(Delta <= Eps1) Delta = 1.;
191 E2 = VScnd1.Norm2() * Delta;
193 E3 = (Delta > Eps1) ? VScnd2.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
194 // ========== Treatment of internal points =================
196 Standard_Integer CurrPoint = 2;
198 for(ipnt = myFirstPoint + 1; ipnt < myLastPoint; ipnt++) {
200 Delta = .5 * (myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1));
203 if(ipnt + 1 != myLastPoint) {
204 EstTangent(ipnt + 2, VTang3);
205 EstSecnd(ipnt + 1, VTang1, VTang3, Length, VScnd2);
208 EstSecnd(ipnt + 1, VTang1, VTang2, Length, VScnd2);
210 E2 += VScnd1.Norm2() * Delta;
211 E3 += (Delta > Eps1) ? VScnd2.Subtracted(VScnd3).Norm2() / (4. * Delta) : 0.;
214 else if(CurrPoint == 2) {
215 if(ipnt + 1 != myLastPoint) {
216 EstTangent(ipnt + 2, VTang1);
217 EstSecnd(ipnt + 1, VTang2, VTang1, Length, VScnd3);
220 EstSecnd(ipnt + 1, VTang2, VTang3, Length, VScnd3);
222 E2 += VScnd2.Norm2() * Delta;
223 E3 += (Delta > Eps1) ? VScnd3.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
227 if(ipnt + 1 != myLastPoint) {
228 EstTangent(ipnt + 2, VTang2);
229 EstSecnd(ipnt + 1, VTang3, VTang2, Length, VScnd1);
232 EstSecnd(ipnt + 1, VTang3, VTang1, Length, VScnd1);
234 E2 += VScnd3.Norm2() * Delta;
235 E3 += (Delta > Eps1) ? VScnd1.Subtracted(VScnd2).Norm2() / (4. * Delta) : 0.;
239 CurrPoint++; if(CurrPoint == 4) CurrPoint = 1;
242 // ========== Treatment of last point =================
244 Delta = .5 * (myParameters->Value(myLastPoint) - myParameters->Value(myLastPoint - 1));
245 if(Delta <= Eps1) Delta = 1.;
251 E2 += VScnd1.Norm2() * Delta;
252 aux = VScnd1.Subtracted(VScnd3).Norm2();
253 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
256 else if(CurrPoint == 2) {
258 E2 += VScnd2.Norm2() * Delta;
259 aux = VScnd2.Subtracted(VScnd1).Norm2();
260 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
265 E2 += VScnd3.Norm2() * Delta;
266 aux = VScnd3.Subtracted(VScnd2).Norm2();
267 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
271 aux = Length * Length;
273 E2 *= aux; E3 *= aux;
279 //=======================================================================
280 //function : EstTangent
281 //purpose : Calculation of estimation of the Tangent
282 // (see fortran routine MMLIPRI)
283 //=======================================================================
286 void AppParCurves_Variational::EstTangent(const Standard_Integer ipnt,
287 math_Vector& VTang) const
291 const Standard_Real Eps1 = Precision::Confusion() * .01;
292 const Standard_Real EpsNorm = 1.e-9;
294 Standard_Real Wpnt = 1.;
297 if(ipnt == myFirstPoint) {
298 // Estimation at first point
303 Standard_Integer adr1 = 1, adr2 = adr1 + myDimension,
304 adr3 = adr2 + myDimension;
306 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
307 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
308 math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
310 // Parabolic interpolation
311 // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
312 // first derivative for t=0 is A1 = ((d2-1)*P1 + P2 - d2*P3)/(d*(1-d))
313 // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*d
314 Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
315 Standard_Real V2 = 0.;
316 if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
318 Standard_Real d = V1 / (V1 + V2), d1;
319 d1 = 1. / (d * (1 - d)); d *= d;
320 VTang = ((d - 1.) * Pnt1 + Pnt2 - d * Pnt3) * d1;
323 // Simple 2-point estimation
329 else if (ipnt == myLastPoint) {
330 // Estimation at last point
335 Standard_Integer adr1 = (myLastPoint - 3) * myDimension + 1,
336 adr2 = adr1 + myDimension,
337 adr3 = adr2 + myDimension;
339 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
340 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
341 math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
343 // Parabolic interpolation
344 // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
345 // first derivative for t=1 is 2*A2 + A1 = ((d2+1)*P1 - P2 - d2*P3)/(d*(1-d))
346 // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*(d-2)
347 Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
348 Standard_Real V2 = 0.;
349 if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
351 Standard_Real d = V1 / (V1 + V2), d1;
352 d1 = 1. / (d * (1 - d)); d *= d - 2;
353 VTang = ((d + 1.) * Pnt1 - Pnt2 - d * Pnt3) * d1;
356 // Simple 2-point estimation
364 Standard_Integer adr1 = (ipnt - myFirstPoint - 1) * myDimension + 1,
365 adr2 = adr1 + 2 * myDimension;
367 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
368 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
374 Standard_Real Vnorm = VTang.Norm();
381 // Estimation with constraints
383 Standard_Real Wcnt = 0.;
384 Standard_Integer IdCnt = 1;
386 // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
388 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
389 math_Vector VCnt(1, myDimension, 0.);
393 while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
394 IdCnt <= NbConstr) IdCnt++;
395 if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
396 (myTypConstraints->Value(2 * IdCnt) >= 1)) {
398 Standard_Integer i0 = 2 * myDimension * (IdCnt - 1), k = 0;
399 for( i = 1; i <= myNbP3d; i++) {
400 for(Standard_Integer j = 1; j <= 3; j++)
401 VCnt(++k) = myTabConstraints->Value(++i0);
404 for(i = 1; i <= myNbP2d; i++) {
405 for(Standard_Integer j = 1; j <= 2; j++)
406 VCnt(++k) = myTabConstraints->Value(++i0);
412 // Averaging of estimation
414 Standard_Real Denom = Wpnt + Wcnt;
415 if(Denom == 0.) Denom = 1.;
416 else Denom = 1. / Denom;
418 VTang = (Wpnt * VTang + Wcnt * VCnt) * Denom;
420 Vnorm = VTang.Norm();
431 //=======================================================================
432 //function : EstSecnd
433 //purpose : Calculation of estimation of the second derivative
434 // (see fortran routine MLIMSCN)
435 //=======================================================================
437 void AppParCurves_Variational::EstSecnd(const Standard_Integer ipnt,
438 const math_Vector& VTang1,
439 const math_Vector& VTang2,
440 const Standard_Real Length,
441 math_Vector& VScnd) const
445 const Standard_Real Eps = 1.e-9;
447 Standard_Real Wpnt = 1.;
451 if(ipnt == myFirstPoint)
452 aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt);
453 else if(ipnt == myLastPoint)
454 aux = myParameters->Value(ipnt) - myParameters->Value(ipnt - 1);
456 aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1);
463 VScnd = (VTang2 - VTang1) * aux;
465 // Estimation with constraints
467 Standard_Real Wcnt = 0.;
468 Standard_Integer IdCnt = 1;
470 // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
472 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
473 math_Vector VCnt(1, myDimension, 0.);
477 while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
478 IdCnt <= NbConstr) IdCnt++;
480 if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
481 (myTypConstraints->Value(2 * IdCnt) >= 2)) {
483 Standard_Integer i0 = 2 * myDimension * (IdCnt - 1) + 3, k = 0;
484 for( i = 1; i <= myNbP3d; i++) {
485 for(Standard_Integer j = 1; j <= 3; j++)
486 VCnt(++k) = myTabConstraints->Value(++i0);
490 for(i = 1; i <= myNbP2d; i++) {
491 for(Standard_Integer j = 1; j <= 2; j++)
492 VCnt(++k) = myTabConstraints->Value(++i0);
498 // Averaging of estimation
500 Standard_Real Denom = Wpnt + Wcnt;
501 if(Denom == 0.) Denom = 1.;
502 else Denom = 1. / Denom;
504 VScnd = (Wpnt * VScnd + (Wcnt * Length) * VCnt) * Denom;
510 //=======================================================================
511 //function : InitCutting
512 //purpose : Realisation of curve's cutting a priory accordingly to
513 // constraints (see fortran routine MLICUT)
514 //=======================================================================
516 void AppParCurves_Variational::InitCutting(const Handle(PLib_Base)& aBase,
517 const Standard_Real CurvTol,
518 Handle(FEmTool_Curve)& aCurve) const
521 // Definition of number of elements
522 Standard_Integer ORCMx = -1, NCont = 0, i, kk, NbElem;
523 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
525 for(i = 1; i <= NbConstr; i++) {
526 kk = Abs(myTypConstraints->Value(2 * i)) + 1;
527 ORCMx = Max(ORCMx, kk);
531 if(ORCMx > myMaxDegree - myNivCont)
532 Standard_ConstructionError::Raise("AppParCurves_Variational::InitCutting");
534 Standard_Integer NLibre = Max(myMaxDegree - myNivCont - (myMaxDegree + 1) / 4,
537 NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
539 while((NbElem > myMaxSegment) && (NLibre < myMaxDegree - myNivCont)) {
542 NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
547 if(NbElem > myMaxSegment)
548 Standard_ConstructionError::Raise("AppParCurves_Variational::InitCutting");
551 aCurve = new FEmTool_Curve(myDimension, NbElem, aBase, CurvTol);
553 Standard_Integer NCnt = (NCont - 1) / NbElem + 1;
554 Standard_Integer NPlus = NbElem - (NCnt * NbElem - NCont);
556 TColStd_Array1OfReal& Knot = aCurve->Knots();
558 Standard_Integer IDeb = 0, IFin = NbConstr + 1,
560 IndEl = Knot.Lower(), IUpper = Knot.Upper(),
564 Knot(IndEl) = myParameters->Value(myFirstPoint);
565 Knot(IUpper) = myParameters->Value(myLastPoint);
567 while(NbElem - NbEl > 1) {
570 if(NPlus == 0) NCnt--;
572 while(NDeb < NCnt && IDeb < IFin) {
574 NDeb += Abs(myTypConstraints->Value(2 * IDeb)) + 1;
580 myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) > Knot(IndEl - 1))
582 Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
584 Knot(IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) +
585 myParameters->Value(myTypConstraints->Value(2 * IDeb + 1))) / 2;
590 Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
594 if(NPlus == 0) NCnt--;
596 if(NbElem - NbEl == 1) break;
600 while(NFin < NCnt && IDeb < IFin) {
602 NFin += Abs(myTypConstraints->Value(2 * IFin)) + 1;
607 Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
608 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
612 if(myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) < Knot(IUpper - IndEl + 1))
613 Knot(IUpper + 1 - IndEl) = myParameters->Value(myTypConstraints->Value(2 * IFin - 1));
615 Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
616 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;