0023024: Update headers of OCCT files
[occt.git] / src / AppParCurves / AppParCurves_Variational_6.gxx
CommitLineData
b311480e 1// Created on: 1998-12-21
2// Created by: Igor FEOKTISTOV
3// Copyright (c) 1998-1999 Matra Datavision
4// Copyright (c) 1999-2012 OPEN CASCADE SAS
5//
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.
10//
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.
13//
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.
20
7fd59977 21
22
23// Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559
24
25#include <PLib_Base.hxx>
26#include <PLib_JacobiPolynomial.hxx>
27#include <PLib_HermitJacobi.hxx>
28#include <FEmTool_Curve.hxx>
29#include <math_Vector.hxx>
30#include <TColStd_Array1OfReal.hxx>
31
7fd59977 32//=======================================================================
33//function : InitSmoothCriterion
34//purpose : Initializes the SmoothCriterion
7fd59977 35//=======================================================================
7fd59977 36void AppParCurves_Variational::InitSmoothCriterion()
37{
38
39 const Standard_Real Eps2 = 1.e-6, Eps3 = 1.e-9;
40// const Standard_Real J1 = .01, J2 = .001, J3 = .001;
41
42
43
44 Standard_Real Length;
45
46 InitParameters(Length);
47
48 mySmoothCriterion->SetParameters(myParameters);
49
50 Standard_Real E1, E2, E3;
51
52 InitCriterionEstimations(Length, E1, E2, E3);
53/*
54 J1 = 1.e-8; J2 = J3 = (E1 + 1.e-8) * 1.e-6;
55
56 if(E1 < J1) E1 = J1;
57 if(E2 < J2) E2 = J2;
58 if(E3 < J3) E3 = J3;
59*/
60 mySmoothCriterion->EstLength() = Length;
61 mySmoothCriterion->SetEstimation(E1, E2, E3);
62
63 Standard_Real WQuadratic, WQuality;
64
65 if(!myWithMinMax && myTolerance != 0.)
66 WQuality = myTolerance;
67 else if (myTolerance == 0.)
68 WQuality = 1.;
69 else
70 WQuality = Max(myTolerance, Eps2 * Length);
71
72 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
73 WQuadratic = Sqrt((Standard_Real)(myNbPoints - NbConstr)) * WQuality;
74 if(WQuadratic > Eps3) WQuadratic = 1./ WQuadratic;
75
76 if(WQuadratic == 0.) WQuadratic = Max(Sqrt(E1), 1.);
77
78
79 mySmoothCriterion->SetWeight(WQuadratic, WQuality,
80 myPercent[0], myPercent[1], myPercent[2]);
81
82
83 Handle(PLib_Base) TheBase = new PLib_HermitJacobi(myMaxDegree, myContinuity);
84 Handle(FEmTool_Curve) TheCurve;
85 Standard_Integer NbElem;
86 Standard_Real CurvTol = Eps2 * Length / myNbPoints;
87
88 // Decoupe de l'intervalle en fonction des contraintes
89 if(myWithCutting == Standard_True && NbConstr != 0) {
90
91 InitCutting(TheBase, CurvTol, TheCurve);
92
93 }
94 else {
95
96 NbElem = 1;
97 TheCurve = new FEmTool_Curve(myDimension, NbElem, TheBase, CurvTol);
98 TheCurve->Knots().SetValue(TheCurve->Knots().Lower(),
99 myParameters->Value(myFirstPoint));
100 TheCurve->Knots().SetValue(TheCurve->Knots().Upper(),
101 myParameters->Value(myLastPoint));
102
103 }
104
105 mySmoothCriterion->SetCurve(TheCurve);
106
107 return;
108
109}
110
111//
112//=======================================================================
113//function : InitParameters
114//purpose : Calculation of initial estimation of the multicurve's length
115// and parameters for multipoints.
116//=======================================================================
117//
118void AppParCurves_Variational::InitParameters(Standard_Real& Length)
119{
120
121 const Standard_Real Eps1 = Precision::Confusion() * .01;
122
123 Standard_Real aux, dist;
124 Standard_Integer i, i0, i1 = 0, ipoint;
125
126
127 Length = 0.;
128 myParameters->SetValue(myFirstPoint, Length);
129
130 for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint; ipoint++) {
131 i0 = i1; i1 += myDimension;
132 dist = 0;
133 for(i = 1; i <= myDimension; i++) {
134 aux = myTabPoints->Value(i1 + i) - myTabPoints->Value(i0 + i);
135 dist += aux * aux;
136 }
137 Length += Sqrt(dist);
138 myParameters->SetValue(ipoint, Length);
139 }
140
141
142 if(Length <= Eps1)
143 Standard_ConstructionError::Raise("AppParCurves_Variational::InitParameters");
144
145
146 for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint - 1; ipoint++)
147 myParameters->SetValue(ipoint, myParameters->Value(ipoint) / Length);
148
149 myParameters->SetValue(myLastPoint, 1.);
150
151 // Avec peu de point il y a sans doute sous estimation ...
152 if(myNbPoints < 10) Length *= (1. + 0.1 / (myNbPoints - 1));
153}
154
155//
156//=======================================================================
157//function : InitCriterionEstimations
158//purpose : Calculation of initial estimation of the linear criteria
159//
160//=======================================================================
161//
162void AppParCurves_Variational::InitCriterionEstimations(const Standard_Real Length,
163 Standard_Real& E1,
164 Standard_Real& E2,
165 Standard_Real& E3) const
166{
167 E1 = Length * Length;
168
169 const Standard_Real Eps1 = Precision::Confusion() * .01;
170
171 math_Vector VTang1(1, myDimension), VTang2(1, myDimension), VTang3(1, myDimension),
172 VScnd1(1, myDimension), VScnd2(1, myDimension), VScnd3(1, myDimension);
173
174 // ========== Treatment of first point =================
175
176 Standard_Integer ipnt = myFirstPoint;
177
178 EstTangent(ipnt, VTang1);
179 ipnt++;
180 EstTangent(ipnt, VTang2);
181 ipnt++;
182 EstTangent(ipnt, VTang3);
183
184 ipnt = myFirstPoint;
185 EstSecnd(ipnt, VTang1, VTang2, Length, VScnd1);
186 ipnt++;
187 EstSecnd(ipnt, VTang1, VTang3, Length, VScnd2);
188
189// Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 Begin
190// Standard_Real Delta = .5 * (myParameters->Value(ipnt) - myParameters->Value(--ipnt));
191 Standard_Integer anInd = ipnt;
192 Standard_Real Delta = .5 * (myParameters->Value(anInd) - myParameters->Value(--ipnt));
193// Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 End
194
195 if(Delta <= Eps1) Delta = 1.;
196
197 E2 = VScnd1.Norm2() * Delta;
198
199 E3 = (Delta > Eps1) ? VScnd2.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
200 // ========== Treatment of internal points =================
201
202 Standard_Integer CurrPoint = 2;
203
204 for(ipnt = myFirstPoint + 1; ipnt < myLastPoint; ipnt++) {
205
206 Delta = .5 * (myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1));
207
208 if(CurrPoint == 1) {
209 if(ipnt + 1 != myLastPoint) {
210 EstTangent(ipnt + 2, VTang3);
211 EstSecnd(ipnt + 1, VTang1, VTang3, Length, VScnd2);
212 }
213 else
214 EstSecnd(ipnt + 1, VTang1, VTang2, Length, VScnd2);
215
216 E2 += VScnd1.Norm2() * Delta;
217 E3 += (Delta > Eps1) ? VScnd2.Subtracted(VScnd3).Norm2() / (4. * Delta) : 0.;
218
219 }
220 else if(CurrPoint == 2) {
221 if(ipnt + 1 != myLastPoint) {
222 EstTangent(ipnt + 2, VTang1);
223 EstSecnd(ipnt + 1, VTang2, VTang1, Length, VScnd3);
224 }
225 else
226 EstSecnd(ipnt + 1, VTang2, VTang3, Length, VScnd3);
227
228 E2 += VScnd2.Norm2() * Delta;
229 E3 += (Delta > Eps1) ? VScnd3.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
230
231 }
232 else {
233 if(ipnt + 1 != myLastPoint) {
234 EstTangent(ipnt + 2, VTang2);
235 EstSecnd(ipnt + 1, VTang3, VTang2, Length, VScnd1);
236 }
237 else
238 EstSecnd(ipnt + 1, VTang3, VTang1, Length, VScnd1);
239
240 E2 += VScnd3.Norm2() * Delta;
241 E3 += (Delta > Eps1) ? VScnd1.Subtracted(VScnd2).Norm2() / (4. * Delta) : 0.;
242
243 }
244
245 CurrPoint++; if(CurrPoint == 4) CurrPoint = 1;
246 }
247
248 // ========== Treatment of last point =================
249
250 Delta = .5 * (myParameters->Value(myLastPoint) - myParameters->Value(myLastPoint - 1));
251 if(Delta <= Eps1) Delta = 1.;
252
253 Standard_Real aux;
254
255 if(CurrPoint == 1) {
256
257 E2 += VScnd1.Norm2() * Delta;
258 aux = VScnd1.Subtracted(VScnd3).Norm2();
259 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
260
261 }
262 else if(CurrPoint == 2) {
263
264 E2 += VScnd2.Norm2() * Delta;
265 aux = VScnd2.Subtracted(VScnd1).Norm2();
266 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
267
268 }
269 else {
270
271 E2 += VScnd3.Norm2() * Delta;
272 aux = VScnd3.Subtracted(VScnd2).Norm2();
273 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
274
275 }
276
277 aux = Length * Length;
278
279 E2 *= aux; E3 *= aux;
280
281
282}
283
284//
285//=======================================================================
286//function : EstTangent
287//purpose : Calculation of estimation of the Tangent
288// (see fortran routine MMLIPRI)
289//=======================================================================
290//
291
292void AppParCurves_Variational::EstTangent(const Standard_Integer ipnt,
293 math_Vector& VTang) const
294
295{
296 Standard_Integer i ;
297 const Standard_Real Eps1 = Precision::Confusion() * .01;
298 const Standard_Real EpsNorm = 1.e-9;
299
300 Standard_Real Wpnt = 1.;
301
302
303 if(ipnt == myFirstPoint) {
304 // Estimation at first point
305 if(myNbPoints < 3)
306 Wpnt = 0.;
307 else {
308
309 Standard_Integer adr1 = 1, adr2 = adr1 + myDimension,
310 adr3 = adr2 + myDimension;
311
312 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
313 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
314 math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
315
316 // Parabolic interpolation
317 // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
318 // first derivative for t=0 is A1 = ((d2-1)*P1 + P2 - d2*P3)/(d*(1-d))
319 // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*d
320 Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
321 Standard_Real V2 = 0.;
322 if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
323 if(V2 > Eps1) {
324 Standard_Real d = V1 / (V1 + V2), d1;
325 d1 = 1. / (d * (1 - d)); d *= d;
326 VTang = ((d - 1.) * Pnt1 + Pnt2 - d * Pnt3) * d1;
327 }
328 else {
329 // Simple 2-point estimation
330
331 VTang = Pnt2 - Pnt1;
332 }
333 }
334 }
335 else if (ipnt == myLastPoint) {
336 // Estimation at last point
337 if(myNbPoints < 3)
338 Wpnt = 0.;
339 else {
340
341 Standard_Integer adr1 = (myLastPoint - 3) * myDimension + 1,
342 adr2 = adr1 + myDimension,
343 adr3 = adr2 + myDimension;
344
345 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
346 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
347 math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
348
349 // Parabolic interpolation
350 // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
351 // first derivative for t=1 is 2*A2 + A1 = ((d2+1)*P1 - P2 - d2*P3)/(d*(1-d))
352 // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*(d-2)
353 Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
354 Standard_Real V2 = 0.;
355 if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
356 if(V2 > Eps1) {
357 Standard_Real d = V1 / (V1 + V2), d1;
358 d1 = 1. / (d * (1 - d)); d *= d - 2;
359 VTang = ((d + 1.) * Pnt1 - Pnt2 - d * Pnt3) * d1;
360 }
361 else {
362 // Simple 2-point estimation
363
364 VTang = Pnt3 - Pnt2;
365 }
366 }
367 }
368 else {
369
370 Standard_Integer adr1 = (ipnt - myFirstPoint - 1) * myDimension + 1,
371 adr2 = adr1 + 2 * myDimension;
372
373 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
374 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
375
376 VTang = Pnt2 - Pnt1;
377
378 }
379
380 Standard_Real Vnorm = VTang.Norm();
381
382 if(Vnorm <= EpsNorm)
383 VTang.Init(0.);
384 else
385 VTang /= Vnorm;
386
387 // Estimation with constraints
388
389 Standard_Real Wcnt = 0.;
390 Standard_Integer IdCnt = 1;
391
392// Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
393
394 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
395 math_Vector VCnt(1, myDimension, 0.);
396
397 if(NbConstr > 0) {
398
399 while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
400 IdCnt <= NbConstr) IdCnt++;
401 if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
402 (myTypConstraints->Value(2 * IdCnt) >= 1)) {
403 Wcnt = 1.;
404 Standard_Integer i0 = 2 * myDimension * (IdCnt - 1), k = 0;
405 for( i = 1; i <= myNbP3d; i++) {
406 for(Standard_Integer j = 1; j <= 3; j++)
407 VCnt(++k) = myTabConstraints->Value(++i0);
408 i0 += 3;
409 }
410 for(i = 1; i <= myNbP2d; i++) {
411 for(Standard_Integer j = 1; j <= 2; j++)
412 VCnt(++k) = myTabConstraints->Value(++i0);
413 i0 += 2;
414 }
415 }
416 }
417
418 // Averaging of estimation
419
420 Standard_Real Denom = Wpnt + Wcnt;
421 if(Denom == 0.) Denom = 1.;
422 else Denom = 1. / Denom;
423
424 VTang = (Wpnt * VTang + Wcnt * VCnt) * Denom;
425
426 Vnorm = VTang.Norm();
427
428 if(Vnorm <= EpsNorm)
429 VTang.Init(0.);
430 else
431 VTang /= Vnorm;
432
433
434}
435
436//
437//=======================================================================
438//function : EstSecnd
439//purpose : Calculation of estimation of the second derivative
440// (see fortran routine MLIMSCN)
441//=======================================================================
442//
443void AppParCurves_Variational::EstSecnd(const Standard_Integer ipnt,
444 const math_Vector& VTang1,
445 const math_Vector& VTang2,
446 const Standard_Real Length,
447 math_Vector& VScnd) const
448{
449 Standard_Integer i ;
450
451 const Standard_Real Eps = 1.e-9;
452
453 Standard_Real Wpnt = 1.;
454
455 Standard_Real aux;
456
457 if(ipnt == myFirstPoint)
458 aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt);
459 else if(ipnt == myLastPoint)
460 aux = myParameters->Value(ipnt) - myParameters->Value(ipnt - 1);
461 else
462 aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1);
463
464 if(aux <= Eps)
465 aux = 1.;
466 else
467 aux = 1. / aux;
468
469 VScnd = (VTang2 - VTang1) * aux;
470
471 // Estimation with constraints
472
473 Standard_Real Wcnt = 0.;
474 Standard_Integer IdCnt = 1;
475
476// Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
477
478 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
479 math_Vector VCnt(1, myDimension, 0.);
480
481 if(NbConstr > 0) {
482
483 while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
484 IdCnt <= NbConstr) IdCnt++;
485
486 if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
487 (myTypConstraints->Value(2 * IdCnt) >= 2)) {
488 Wcnt = 1.;
489 Standard_Integer i0 = 2 * myDimension * (IdCnt - 1) + 3, k = 0;
490 for( i = 1; i <= myNbP3d; i++) {
491 for(Standard_Integer j = 1; j <= 3; j++)
492 VCnt(++k) = myTabConstraints->Value(++i0);
493 i0 += 3;
494 }
495 i0--;
496 for(i = 1; i <= myNbP2d; i++) {
497 for(Standard_Integer j = 1; j <= 2; j++)
498 VCnt(++k) = myTabConstraints->Value(++i0);
499 i0 += 2;
500 }
501 }
502 }
503
504 // Averaging of estimation
505
506 Standard_Real Denom = Wpnt + Wcnt;
507 if(Denom == 0.) Denom = 1.;
508 else Denom = 1. / Denom;
509
510 VScnd = (Wpnt * VScnd + (Wcnt * Length) * VCnt) * Denom;
511
512}
513
514
515//
516//=======================================================================
517//function : InitCutting
518//purpose : Realisation of curve's cutting a priory accordingly to
519// constraints (see fortran routine MLICUT)
520//=======================================================================
521//
522void AppParCurves_Variational::InitCutting(const Handle(PLib_Base)& aBase,
523 const Standard_Real CurvTol,
524 Handle(FEmTool_Curve)& aCurve) const
525{
526
527 // Definition of number of elements
528 Standard_Integer ORCMx = -1, NCont = 0, i, kk, NbElem;
529 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
530
531 for(i = 1; i <= NbConstr; i++) {
532 kk = Abs(myTypConstraints->Value(2 * i)) + 1;
533 ORCMx = Max(ORCMx, kk);
534 NCont += kk;
535 }
536
537 if(ORCMx > myMaxDegree - myNivCont)
538 Standard_ConstructionError::Raise("AppParCurves_Variational::InitCutting");
539
540 Standard_Integer NLibre = Max(myMaxDegree - myNivCont - (myMaxDegree + 1) / 4,
541 myNivCont + 1);
542
543 NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
544
545 while((NbElem > myMaxSegment) && (NLibre < myMaxDegree - myNivCont)) {
546
547 NLibre++;
548 NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
549
550 }
551
552
553 if(NbElem > myMaxSegment)
554 Standard_ConstructionError::Raise("AppParCurves_Variational::InitCutting");
555
556
557 aCurve = new FEmTool_Curve(myDimension, NbElem, aBase, CurvTol);
558
559 Standard_Integer NCnt = (NCont - 1) / NbElem + 1;
560 Standard_Integer NPlus = NbElem - (NCnt * NbElem - NCont);
561
562 TColStd_Array1OfReal& Knot = aCurve->Knots();
563
564 Standard_Integer IDeb = 0, IFin = NbConstr + 1,
565 NDeb = 0, NFin = 0,
566 IndEl = Knot.Lower(), IUpper = Knot.Upper(),
567 NbEl = 0;
568
569
570 Knot(IndEl) = myParameters->Value(myFirstPoint);
571 Knot(IUpper) = myParameters->Value(myLastPoint);
572
573 while(NbElem - NbEl > 1) {
574
575 IndEl++; NbEl++;
576 if(NPlus == 0) NCnt--;
577
578 while(NDeb < NCnt && IDeb < IFin) {
579 IDeb++;
580 NDeb += Abs(myTypConstraints->Value(2 * IDeb)) + 1;
581 }
582
583 if(NDeb == NCnt) {
584 NDeb = 0;
585 if(NPlus == 1 &&
586 myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) > Knot(IndEl - 1))
587
588 Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
589 else
590 Knot(IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) +
591 myParameters->Value(myTypConstraints->Value(2 * IDeb + 1))) / 2;
592
593 }
594 else {
595 NDeb -= NCnt;
596 Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
597 }
598
599 NPlus--;
600 if(NPlus == 0) NCnt--;
601
602 if(NbElem - NbEl == 1) break;
603
604 NbEl++;
605
606 while(NFin < NCnt && IDeb < IFin) {
607 IFin--;
608 NFin += Abs(myTypConstraints->Value(2 * IFin)) + 1;
609 }
610
611 if(NFin == NCnt) {
612 NFin = 0;
613 Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
614 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
615 }
616 else {
617 NFin -= NCnt;
618 if(myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) < Knot(IUpper - IndEl + 1))
619 Knot(IUpper + 1 - IndEl) = myParameters->Value(myTypConstraints->Value(2 * IFin - 1));
620 else
621 Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
622 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
623 }
624
625 }
626
627
628}