projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0024510: Remove unused local variables
[occt.git] / src / GeomFill / GeomFill_CorrectedFrenet.cxx
CommitLineData
b311480e 1// Created on: 1997-12-19
2// Created by: Roman BORISOV /Philippe MANGIN
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//
973c2be1 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.
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 <stdio.h>
18
19#include <GeomFill_CorrectedFrenet.ixx>
20#include <GeomAbs_CurveType.hxx>
21#include <Adaptor3d_HCurve.hxx>
22#include <gp_Trsf.hxx>
23#include <Precision.hxx>
24#include <TColStd_HArray1OfReal.hxx>
25#include <Law_Interpolate.hxx>
26#include <TColStd_SequenceOfReal.hxx>
27#include <gp_Vec2d.hxx>
28#include <BndLib_Add3dCurve.hxx>
29#include <Bnd_Box.hxx>
30#include <GeomLib.hxx>
31#include <Law_Composite.hxx>
32#include <Law_Constant.hxx>
33#include <Law_BSpFunc.hxx>
34#include <Law_BSpline.hxx>
35#include <GeomFill_SnglrFunc.hxx>
36//Patch
37#include <Geom_Plane.hxx>
38#include <Geom_BezierCurve.hxx>
39#include <Geom_BSplineCurve.hxx>
40#include <TColgp_HArray1OfPnt.hxx>
41
a31abc03 42
7fd59977 43#ifdef DEB
44static Standard_Boolean Affich=0;
45#endif
46
47#ifdef DRAW
48static Standard_Integer CorrNumber = 0;
49#include <Draw_Appli.hxx>
50#include <DrawTrSurf.hxx>
51#include <Draw_Segment2D.hxx>
52//#include <Draw.hxx>
53#include <TColgp_Array1OfPnt.hxx>
54#include <TColStd_Array1OfReal.hxx>
55#include <TColStd_HArray1OfInteger.hxx>
56#endif
57
58#ifdef DRAW
59static void draw(const Handle(Law_Function)& law)
60{
61 Standard_Real Step, u, v, tmin;
62 Standard_Integer NbInt, i, j, jmax;
63 NbInt = law->NbIntervals(GeomAbs_C3);
64 TColStd_Array1OfReal Int(1, NbInt+1);
65 law->Intervals(Int, GeomAbs_C3);
66 gp_Pnt2d old;
67 Handle(Draw_Segment2D) tg2d;
68
69 for(i = 1; i <= NbInt; i++){
70 tmin = Int(i);
71 Step = (Int(i+1)-Int(i))/4;
72 if (i == NbInt) jmax = 4;
73 else jmax = 3;
74 for (j=1; j<=jmax; j++) {
75 u = tmin + (j-1)*Step;
76 v = law->Value(u);
77 gp_Pnt2d point2d(u,v);
78 if ((i>1)||(j>1)) {
79 tg2d = new Draw_Segment2D(old, point2d,Draw_kaki);
80 dout << tg2d;
81 }
82 old = point2d;
83 }
84 }
85 dout.Flush();
86}
87#endif
88
a31abc03 89
90static Standard_Real ComputeTorsion(const Standard_Real Param,
91 const Handle(Adaptor3d_HCurve)& aCurve)
92{
93 Standard_Real Torsion;
94
95 gp_Pnt aPoint;
96 gp_Vec DC1, DC2, DC3;
97 aCurve->D3(Param, aPoint, DC1, DC2, DC3);
98 gp_Vec DC1crossDC2 = DC1 ^ DC2;
99 Standard_Real Norm_DC1crossDC2 = DC1crossDC2.Magnitude();
100
101 Standard_Real DC1DC2DC3 = DC1crossDC2 * DC3 ; //mixed product
102
103 Standard_Real Tol = gp::Resolution();
104 Standard_Real SquareNorm_DC1crossDC2 = Norm_DC1crossDC2 * Norm_DC1crossDC2;
105 if (SquareNorm_DC1crossDC2 <= Tol)
106 Torsion = 0.;
107 else
108 Torsion = DC1DC2DC3 / SquareNorm_DC1crossDC2 ;
109
110 return Torsion;
111}
112
7fd59977 113//===============================================================
114// Function : smoothlaw
115// Purpose : to smooth a law : Reduce the number of knots
116//===============================================================
117static void smoothlaw(Handle(Law_BSpline)& Law,
118 const Handle(TColStd_HArray1OfReal)& Points,
119 const Handle(TColStd_HArray1OfReal)& Param,
120 const Standard_Real Tol)
121{
122 Standard_Real tol, d;
123 Standard_Integer ii, Nbk;
124 Standard_Boolean B, Ok;
125 Handle(Law_BSpline) BS = Law->Copy();
126
127 Nbk = BS->NbKnots();
128 tol = Tol/10;
129 Ok = Standard_False;
130
131 for (ii=Nbk-1; ii>1; ii--) { // Une premiere passe tolerance serres
132 B = BS->RemoveKnot(ii, 0, tol);
133 if (B) Ok = Standard_True;
134 }
135
136 if (Ok) { // controle
137 tol = 0.;
138 for (ii=1; ii<=Param->Length() && Ok; ii++) {
139 d = Abs(BS->Value(Param->Value(ii))-Points->Value(ii));
140 if (d > tol) tol = d;
141 Ok = (tol <= Tol);
142 }
143 if (Ok)
144 tol = (Tol-tol)/2;
145 else {
146#if DEB
147 cout << "smooth law echec" << endl;
148#endif
149 return; // Echec
150 }
151 }
152 else {
153 tol = Tol/2;
154 }
155
156
157 if (Ok) Law = BS;
158
159 Ok = Standard_False; // Une deuxieme passe tolerance desserre
160 Nbk = BS->NbKnots();
161 for (ii=Nbk-1; ii>1; ii--) {
162 B = BS->RemoveKnot(ii, 0, tol);
163 if (B) Ok = Standard_True;
164 }
165
166 if (Ok) { // controle
167 tol = 0.;
168 for (ii=1; ii<=Param->Length() && Ok; ii++) {
169 d = Abs(BS->Value(Param->Value(ii))-Points->Value(ii));
170 if (d > tol) tol = d;
171 Ok = (tol <= Tol);
172 }
173 if (!Ok) {
174#if DEB
175 cout << "smooth law echec" << endl;
176#endif
177 }
178 }
179 if (Ok) Law = BS;
180
181#if DEB
182 if (Affich) {
183 cout << "Knots Law : " << endl;
184 for (ii=1; ii<=BS->NbKnots(); ii++) {
185 cout << ii << " : " << BS->Knot(ii) << endl;
186 }
187 }
188#endif
189}
190
191//===============================================================
192// Function : FindPlane
193// Purpose :
194//===============================================================
195static Standard_Boolean FindPlane ( const Handle(Adaptor3d_HCurve)& c,
196 Handle( Geom_Plane )& P )
197{
198 Standard_Boolean found = Standard_True;
199 Handle(TColgp_HArray1OfPnt) TabP;
200
201 switch (c->GetType()) {
202
203 case GeomAbs_Line:
204 {
205 found = Standard_False;
206 }
207 break;
208
209 case GeomAbs_Circle:
210 P = new Geom_Plane(gp_Ax3(c->Circle().Position()));
211 break;
212
213 case GeomAbs_Ellipse:
214 P = new Geom_Plane(gp_Ax3(c->Ellipse().Position()));
215 break;
216
217 case GeomAbs_Hyperbola:
218 P = new Geom_Plane(gp_Ax3(c->Hyperbola().Position()));
219 break;
220
221 case GeomAbs_Parabola:
222 P = new Geom_Plane(gp_Ax3(c->Parabola().Position()));
223 break;
224
225 case GeomAbs_BezierCurve:
226 {
227 Handle(Geom_BezierCurve) GC = c->Bezier();
228 Standard_Integer nbp = GC->NbPoles();
229 if ( nbp < 2)
230 found = Standard_False;
231 else if ( nbp == 2) {
232 found = Standard_False;
233 }
234 else {
235 TabP = new (TColgp_HArray1OfPnt) (1, nbp);
236 GC->Poles(TabP->ChangeArray1());
237 }
238 }
239 break;
240
241 case GeomAbs_BSplineCurve:
242 {
243 Handle(Geom_BSplineCurve) GC = c->BSpline();
244 Standard_Integer nbp = GC->NbPoles();
245 if ( nbp < 2)
246 found = Standard_False;
247 else if ( nbp == 2) {
248 found = Standard_False;
249 }
250 else {
251 TabP = new (TColgp_HArray1OfPnt) (1, nbp);
252 GC->Poles(TabP->ChangeArray1());
253 }
254 }
255 break;
256
257 default:
258 { // On utilise un echantillonage
259 Standard_Integer nbp = 15 + c->NbIntervals(GeomAbs_C3);
260 Standard_Real f, l, t, inv;
261 Standard_Integer ii;
262 f = c->FirstParameter();
263 l = c->LastParameter();
264 inv = 1./(nbp-1);
265 for (ii=1; ii<=nbp; ii++) {
266 t = ( f*(nbp-ii) + l*(ii-1));
267 t *= inv;
268 TabP->SetValue(ii, c->Value(t));
269 }
270 }
271 }
272
273 if (! TabP.IsNull()) { // Recherche d'un plan moyen et controle
274 Standard_Boolean issingular;
275 gp_Ax2 inertia;
276 GeomLib::AxeOfInertia(TabP->Array1(), inertia, issingular);
277 if (issingular) {
278 found = Standard_False;
279 }
280 else {
281 P = new Geom_Plane(inertia);
282 }
283 if (found)
284 {
285 //control = Controle(TabP->Array1(), P, myTolerance);
286// Standard_Boolean isOnPlane;
287 Standard_Real a,b,c,d, dist;
288 Standard_Integer ii;
289 P->Coefficients(a,b,c,d);
290 for (ii=1; ii<=TabP->Length() && found; ii++) {
291 const gp_XYZ& xyz = TabP->Value(ii).XYZ();
292 dist = a*xyz.X() + b*xyz.Y() + c*xyz.Z() + d;
293 found = (Abs(dist) <= Precision::Confusion());
294 }
295 return found;
296 }
297 }
298
299 return found;
300}
301
302//===============================================================
a31abc03 303// Function : Constructor
7fd59977 304// Purpose :
305//===============================================================
306GeomFill_CorrectedFrenet::GeomFill_CorrectedFrenet()
307 : isFrenet(Standard_False)
308{
309 frenet = new GeomFill_Frenet();
a31abc03 310 myForEvaluation = Standard_False;
311}
312
313//===============================================================
314// Function : Constructor
315// Purpose :
316//===============================================================
317GeomFill_CorrectedFrenet::GeomFill_CorrectedFrenet(const Standard_Boolean ForEvaluation)
318 : isFrenet(Standard_False)
319{
320 frenet = new GeomFill_Frenet();
321 myForEvaluation = ForEvaluation;
7fd59977 322}
323
324Handle(GeomFill_TrihedronLaw) GeomFill_CorrectedFrenet::Copy() const
325{
326 Handle(GeomFill_CorrectedFrenet) copy = new (GeomFill_CorrectedFrenet)();
327 if (!myCurve.IsNull()) copy->SetCurve(myCurve);
328 return copy;
329}
330
331 void GeomFill_CorrectedFrenet::SetCurve(const Handle(Adaptor3d_HCurve)& C)
332{
333
334 GeomFill_TrihedronLaw::SetCurve(C);
335 if (! C.IsNull()) {
336 frenet->SetCurve(C);
337
338 GeomAbs_CurveType type;
339 type = C->GetType();
340 switch (type) {
341 case GeomAbs_Circle:
342 case GeomAbs_Ellipse:
343 case GeomAbs_Hyperbola:
344 case GeomAbs_Parabola:
345 case GeomAbs_Line:
346 {
347 // No probleme isFrenet
348 isFrenet = Standard_True;
a31abc03 349 break;
7fd59977 350 }
351 default :
352 {
353 // We have to search singulaties
354 isFrenet = Standard_True;
355 Init();
356 }
357 }
358 }
359}
360
361
362//===============================================================
363// Function : Init
364// Purpose : Compute angle's law
365//===============================================================
366 void GeomFill_CorrectedFrenet::Init()
367{
368 EvolAroundT = new Law_Composite();
369 Standard_Integer NbI = frenet->NbIntervals(GeomAbs_C0), i;
370 TColStd_Array1OfReal T(1, NbI + 1);
371 frenet->Intervals(T, GeomAbs_C0);
372 Handle(Law_Function) Func;
373 //OCC78
374 TColStd_SequenceOfReal SeqPoles, SeqAngle;
375 TColgp_SequenceOfVec SeqTangent, SeqNormal;
376
377 gp_Vec Tangent, Normal, BN;
378 frenet->D0(myTrimmed->FirstParameter(), Tangent, Normal, BN);
379 Standard_Integer NbStep;
380// Standard_Real StartAng = 0, AvStep, Step, t;
381 Standard_Real StartAng = 0, AvStep, Step;
382
383#if DRAW
384 Standard_Real t;
385
386 if (Affich) { // Display the curve C'^C''(t)
387 GeomFill_SnglrFunc CS(myCurve);
388 NbStep = 99;
389 AvStep = (myTrimmed->LastParameter() -
390 myTrimmed->FirstParameter())/NbStep;
391 TColgp_Array1OfPnt TabP(1, NbStep+1);
392
393 TColStd_Array1OfReal TI(1, NbStep+1);
394 TColStd_Array1OfInteger M(1,NbStep+1);
395 M.Init(1);
396 M(1) = M(NbStep+1) = 2;
397 for (i=1; i<=NbStep+1; i++) {
398 t = (myTrimmed->FirstParameter()+ (i-1)*AvStep);
399 CS.D0(t, TabP(i));
400 TI(i) = t;
401 }
402 char tname[100];
403 Standard_CString name = tname ;
404 sprintf(name,"Binorm_%d", ++CorrNumber);
405 Handle(Geom_BSplineCurve) BS = new
406 (Geom_BSplineCurve) (TabP, TI, M, 1);
407// DrawTrSurf::Set(&name[0], BS);
408 DrawTrSurf::Set(name, BS);
409 }
410#endif
411
412
413 NbStep = 10;
414 AvStep = (myTrimmed->LastParameter() - myTrimmed->FirstParameter())/NbStep;
415 for(i = 1; i <= NbI; i++) {
416 NbStep = Max(Standard_Integer((T(i+1) - T(i))/AvStep), 3);
417 Step = (T(i+1) - T(i))/NbStep;
a31abc03 418 if(!InitInterval(T(i), T(i+1), Step, StartAng, Tangent, Normal, AT, AN, Func,
419 SeqPoles, SeqAngle, SeqTangent, SeqNormal))
420 {
421 if(isFrenet)
422 isFrenet = Standard_False;
423 }
7fd59977 424 Handle(Law_Composite)::DownCast(EvolAroundT)->ChangeLaws().Append(Func);
425 }
426 if(myTrimmed->IsPeriodic())
427 Handle(Law_Composite)::DownCast(EvolAroundT)->SetPeriodic();
428
429 TLaw = EvolAroundT;
430 //OCC78
431 Standard_Integer iEnd = SeqPoles.Length();
432 HArrPoles = new TColStd_HArray1OfReal(1, iEnd);
433 HArrAngle = new TColStd_HArray1OfReal(1, iEnd);
434 HArrTangent = new TColgp_HArray1OfVec(1, iEnd);
435 HArrNormal = new TColgp_HArray1OfVec(1, iEnd);
436 for(i = 1; i <= iEnd; i++){
437 HArrPoles->ChangeValue(i) = SeqPoles(i);
438 HArrAngle->ChangeValue(i) = SeqAngle(i);
439 HArrTangent->ChangeValue(i) = SeqTangent(i);
440 HArrNormal->ChangeValue(i) = SeqNormal(i);
441 };
a31abc03 442
7fd59977 443#if DRAW
444 if (Affich) {
445 draw(EvolAroundT);
446 }
447#endif
448}
449
450//===============================================================
451// Function : InitInterval
452// Purpose : Compute the angle law on a span
453//===============================================================
454 Standard_Boolean GeomFill_CorrectedFrenet::
455 InitInterval(const Standard_Real First, const Standard_Real Last,
456 const Standard_Real Step,
457 Standard_Real& startAng, gp_Vec& prevTangent,
458 gp_Vec& prevNormal, gp_Vec& aT, gp_Vec& aN,
459 Handle(Law_Function)& FuncInt,
460 TColStd_SequenceOfReal& SeqPoles,
461 TColStd_SequenceOfReal& SeqAngle,
462 TColgp_SequenceOfVec& SeqTangent,
463 TColgp_SequenceOfVec& SeqNormal) const
464{
465 Bnd_Box Boite;
466 gp_Vec Tangent, Normal, BN, cross;
467 TColStd_SequenceOfReal parameters;
468 TColStd_SequenceOfReal EvolAT;
96a95605 469 Standard_Real Param = First, L, norm;
7fd59977 470 Standard_Boolean isZero = Standard_True, isConst = Standard_True;
a31abc03 471 const Standard_Real minnorm = 1.e-16;
7fd59977 472 Standard_Integer i;
473 gp_Pnt PonC;
474 gp_Vec D1;
475
476 frenet->SetInterval(First, Last); //To have the rigth evaluation at bounds
477 GeomFill_SnglrFunc CS(myCurve);
478 BndLib_Add3dCurve::Add(CS, First, Last, 1.e-2, Boite);
7fd59977 479
480 aT = gp_Vec(0, 0, 0);
481 aN = gp_Vec(0, 0, 0);
482
1d47d8d0 483 Standard_Real angleAT = 0., currParam, currStep = Step;
7fd59977 484
485 Handle( Geom_Plane ) aPlane;
a31abc03 486 Standard_Boolean isPlanar = Standard_False;
487 if (!myForEvaluation)
488 isPlanar = FindPlane( myCurve, aPlane );
7fd59977 489
490 i = 1;
491 currParam = Param;
492 Standard_Real DLast = Last - Precision::PConfusion();
493
494 while (Param < Last) {
495 if (currParam > DLast) {
496 currStep = DLast - Param;
497 currParam = Last;
498 }
499 if (isPlanar)
500 currParam = Last;
501
502 frenet->D0(currParam, Tangent, Normal, BN);
c6541a0c 503 if (prevTangent.Angle(Tangent) < M_PI/3 || i == 1) {
7fd59977 504 parameters.Append(currParam);
505 //OCC78
506 SeqPoles.Append(Param);
507 SeqAngle.Append(i > 1? EvolAT(i-1) : startAng);
508 SeqTangent.Append(prevTangent);
509 SeqNormal.Append(prevNormal);
510 angleAT = CalcAngleAT(Tangent,Normal,prevTangent,prevNormal);
511
512 if(isConst && i > 1)
513 if(Abs(angleAT) > Precision::PConfusion())
514 isConst = Standard_False;
515
516 angleAT += (i > 1) ? EvolAT(i-1) : startAng;
517 EvolAT.Append(angleAT);
518 prevNormal = Normal;
519
520 if(isZero)
521 if(Abs(angleAT) > Precision::PConfusion())
522 isZero = Standard_False;
523
524 aT += Tangent;
525 cross = Tangent.Crossed(Normal);
526 aN.SetLinearForm(Sin(angleAT), cross,
527 1 - Cos(angleAT), Tangent.Crossed(cross),
528 Normal+aN);
529 prevTangent = Tangent;
530 Param = currParam;
531 i++;
532
533 //Evaluate the Next step
534 CS.D1(Param, PonC, D1);
a31abc03 535
536 L = PonC.XYZ().Modulus()/2;
7fd59977 537 norm = D1.Magnitude();
a31abc03 538 if (norm <= gp::Resolution())
539 {
540 //norm = 2.*gp::Resolution();
541 norm = minnorm;
7fd59977 542 }
543 currStep = L / norm;
a31abc03 544 if (currStep <= gp::Resolution()) //L = 0 => curvature = 0, linear segment
545 currStep = Step;
546 if (currStep < Precision::Confusion()) //too small step
547 currStep = Precision::Confusion();
548 if (currStep > Step) //too big step
549 currStep = Step;//default value
7fd59977 550 }
551 else
552 currStep /= 2; // Step too long !
553
554 currParam = Param + currStep;
555 }
556
557 if (! isPlanar)
558 {
559 aT /= parameters.Length() - 1;
560 aN /= parameters.Length() - 1;
561 }
562 startAng = angleAT;
563
564// Interpolation
565 if (isConst || isPlanar) {
566 FuncInt = new Law_Constant();
567 Handle(Law_Constant)::DownCast(FuncInt)->Set( angleAT, First, Last );
568 }
569
570 else {
571 Standard_Integer Length = parameters.Length();
572 Handle(TColStd_HArray1OfReal) pararr =
573 new TColStd_HArray1OfReal(1, Length);
574 Handle(TColStd_HArray1OfReal) angleATarr =
575 new TColStd_HArray1OfReal(1, Length);
576
577
578 for (i = 1; i <= Length; i++) {
579 pararr->ChangeValue(i) = parameters(i);
580 angleATarr->ChangeValue(i) = EvolAT(i);
581 }
582
583#if DEB
584 if (Affich) {
585 cout<<"NormalEvolution"<<endl;
586 for (i = 1; i <= Length; i++) {
587 cout<<"("<<pararr->Value(i)<<", "<<angleATarr->Value(i)<<")" << endl;
588 }
589 cout<<endl;
590 }
591#endif
592
593 Law_Interpolate lawAT(angleATarr, pararr,
594 Standard_False, Precision::PConfusion());
595 lawAT.Perform();
596 Handle(Law_BSpline) BS = lawAT.Curve();
597 smoothlaw(BS, angleATarr, pararr, 0.1);
598
599 FuncInt = new Law_BSpFunc(BS, First, Last);
600 }
601 return isZero;
602}
603//===============================================================
604// Function : CalcAngleAT (OCC78)
605// Purpose : Calculate angle of rotation of trihedron normal and its derivatives relative
606// at any position on his curve
607//===============================================================
608Standard_Real GeomFill_CorrectedFrenet::CalcAngleAT(const gp_Vec& Tangent, const gp_Vec& Normal,
609 const gp_Vec& prevTangent, const gp_Vec& prevNormal) const
610{
611 Standard_Real angle;
612 gp_Vec Normal_rot, cross;
613 angle = Tangent.Angle(prevTangent);
614 if (Abs(angle) > Precision::Angular()) {
615 cross = Tangent.Crossed(prevTangent).Normalized();
616 Normal_rot = Normal + sin(angle)*cross.Crossed(Normal) +
617 (1 - cos(angle))*cross.Crossed(cross.Crossed(Normal));
618 }
619 else
620 Normal_rot = Normal;
621 Standard_Real angleAT = Normal_rot.Angle(prevNormal);
c6541a0c 622 if(angleAT > Precision::Angular() && M_PI - angleAT > Precision::Angular())
7fd59977 623 if (Normal_rot.Crossed(prevNormal).IsOpposite(prevTangent, Precision::Angular()))
624 angleAT = -angleAT;
625 return angleAT;
626};
627//===============================================================
628// Function : ... (OCC78)
629// Purpose : This family of functions produce conversion of angle utility
630//===============================================================
7fd59977 631static Standard_Real corr2PI_PI(Standard_Real Ang){
c6541a0c 632 return Ang = (Ang < M_PI? Ang: Ang-2*M_PI);
7fd59977 633};
634static Standard_Real diffAng(Standard_Real A, Standard_Real Ao){
c6541a0c 635 Standard_Real dA = (A-Ao) - Floor((A-Ao)/2.0/M_PI)*2.0*M_PI;
7fd59977 636 return dA = dA >= 0? corr2PI_PI(dA): -corr2PI_PI(-dA);
637};
638//===============================================================
639// Function : CalcAngleAT (OCC78)
640// Purpose : Calculate angle of rotation of trihedron normal and its derivatives relative
641// at any position on his curve
642//===============================================================
643Standard_Real GeomFill_CorrectedFrenet::GetAngleAT(const Standard_Real Param) const{
644 // Search index of low margin from poles of TLaw by bisection method
645 Standard_Integer iB = 1, iE = HArrPoles->Length(), iC = (iE+iB)/2;
646 if(Param == HArrPoles->Value(iB)) return TLaw->Value(Param);
647 if(Param > HArrPoles->Value(iE)) iC = iE;
648 if(iC < iE){
649 while(!(HArrPoles->Value(iC) <= Param && Param <= HArrPoles->Value(iC+1))){
650 if(HArrPoles->Value(iC) < Param) iB = iC; else iE = iC;
651 iC = (iE+iB)/2;
652 };
653 if(HArrPoles->Value(iC) == Param || Param == HArrPoles->Value(iC+1)) return TLaw->Value(Param);
654 };
7fd59977 655 // Calculate differenciation between apporoximated and local values of AngleAT
656 Standard_Real AngP = TLaw->Value(Param), AngPo = HArrAngle->Value(iC), dAng = AngP - AngPo;
657 gp_Vec Tangent, Normal, BN;
658 frenet->D0(Param, Tangent, Normal, BN);
659 Standard_Real DAng = CalcAngleAT(Tangent, Normal, HArrTangent->Value(iC), HArrNormal->Value(iC));
660 Standard_Real DA = diffAng(DAng,dAng);
661 // The correction (there is core of OCC78 bug)
c6541a0c 662 if(Abs(DA) > M_PI/2.0){
7fd59977 663 AngP = AngPo + DAng;
664 };
665 return AngP;
666};
667//===============================================================
668// Function : D0
669// Purpose :
670//===============================================================
671 Standard_Boolean GeomFill_CorrectedFrenet::D0(const Standard_Real Param,
672 gp_Vec& Tangent,
673 gp_Vec& Normal,
674 gp_Vec& BiNormal)
675{
676 frenet->D0(Param, Tangent, Normal, BiNormal);
677 if (isFrenet) return Standard_True;
678
679 Standard_Real angleAT;
680 //angleAT = TLaw->Value(Param);
681 angleAT = GetAngleAT(Param); //OCC78
682
683// rotation around Tangent
684 gp_Vec cross;
685 cross = Tangent.Crossed(Normal);
686 Normal.SetLinearForm(Sin(angleAT), cross,
687 (1 - Cos(angleAT)), Tangent.Crossed(cross),
688 Normal);
689 BiNormal = Tangent.Crossed(Normal);
690
691 return Standard_True;
692}
693
694//===============================================================
695// Function : D1
696// Purpose :
697//===============================================================
698
699 Standard_Boolean GeomFill_CorrectedFrenet::D1(const Standard_Real Param,
700 gp_Vec& Tangent,
701 gp_Vec& DTangent,
702 gp_Vec& Normal,
703 gp_Vec& DNormal,
704 gp_Vec& BiNormal,
705 gp_Vec& DBiNormal)
706{
707 frenet->D1(Param, Tangent, DTangent, Normal, DNormal, BiNormal, DBiNormal);
708 if (isFrenet) return Standard_True;
709
710 Standard_Real angleAT, d_angleAT;
711 Standard_Real sina, cosa;
712
713 TLaw->D1(Param, angleAT, d_angleAT);
714 angleAT = GetAngleAT(Param); //OCC78
715
716 gp_Vec cross, dcross, tcross, dtcross, aux;
717 sina = Sin(angleAT);
718 cosa = Cos(angleAT);
719
720 cross = Tangent.Crossed(Normal);
721 dcross.SetLinearForm(1, DTangent.Crossed(Normal),
722 Tangent.Crossed(DNormal));
723
724 tcross = Tangent.Crossed(cross);
725 dtcross.SetLinearForm(1, DTangent.Crossed(cross),
726 Tangent.Crossed(dcross));
727
728 aux.SetLinearForm(sina, dcross,
729 cosa*d_angleAT, cross);
730 aux.SetLinearForm(1 - cosa, dtcross,
731 sina*d_angleAT, tcross,
732 aux);
733 DNormal+=aux;
734
735 Normal.SetLinearForm( sina, cross,
736 (1 - cosa), tcross,
737 Normal);
738
739 BiNormal = Tangent.Crossed(Normal);
740
741 DBiNormal.SetLinearForm(1, DTangent.Crossed(Normal),
742 Tangent.Crossed(DNormal));
743
744// for test
745/* gp_Vec FDN, Tf, Nf, BNf;
746 Standard_Real h;
747 h = 1.0e-8;
748 if (Param + h > myTrimmed->LastParameter()) h = -h;
749 D0(Param + h, Tf, Nf, BNf);
750 FDN = (Nf - Normal)/h;
751 cout<<"Param = "<<Param<<endl;
752 cout<<"DN = ("<<DNormal.X()<<", "<<DNormal.Y()<<", "<<DNormal.Z()<<")"<<endl;
753 cout<<"FDN = ("<<FDN.X()<<", "<<FDN.Y()<<", "<<FDN.Z()<<")"<<endl;
754*/
755
756 return Standard_True;
757}
758
759//===============================================================
760// Function : D2
761// Purpose :
762//===============================================================
763 Standard_Boolean GeomFill_CorrectedFrenet::D2(const Standard_Real Param,
764 gp_Vec& Tangent,
765 gp_Vec& DTangent,
766 gp_Vec& D2Tangent,
767 gp_Vec& Normal,
768 gp_Vec& DNormal,
769 gp_Vec& D2Normal,
770 gp_Vec& BiNormal,
771 gp_Vec& DBiNormal,
772 gp_Vec& D2BiNormal)
773{
774 frenet->D2(Param, Tangent, DTangent, D2Tangent,
775 Normal, DNormal, D2Normal,
776 BiNormal, DBiNormal, D2BiNormal);
777 if (isFrenet) return Standard_True;
778
779 Standard_Real angleAT, d_angleAT, d2_angleAT;
780 Standard_Real sina, cosa;
781 TLaw->D2(Param, angleAT, d_angleAT, d2_angleAT);
782 angleAT = GetAngleAT(Param); //OCC78
783
784 gp_Vec cross, dcross, d2cross, tcross, dtcross, d2tcross, aux;
785 sina = Sin(angleAT);
786 cosa = Cos(angleAT);
787 cross = Tangent.Crossed(Normal);
788 dcross.SetLinearForm(1, DTangent.Crossed(Normal),
789 Tangent.Crossed(DNormal));
790 d2cross.SetLinearForm(1, D2Tangent.Crossed(Normal),
791 2, DTangent.Crossed(DNormal),
792 Tangent.Crossed(D2Normal));
793
794
795 tcross = Tangent.Crossed(cross);
796 dtcross.SetLinearForm(1, DTangent.Crossed(cross),
797 Tangent.Crossed(dcross));
798 d2tcross.SetLinearForm(1, D2Tangent.Crossed(cross),
799 2, DTangent.Crossed(dcross),
800 Tangent.Crossed(d2cross));
801
802
803 aux.SetLinearForm(sina, d2cross,
804 2*cosa*d_angleAT, dcross,
805 cosa*d2_angleAT - sina*d_angleAT*d_angleAT, cross);
806
807 aux.SetLinearForm(1 - cosa, d2tcross,
808 2*sina*d_angleAT, dtcross,
809 cosa*d_angleAT*d_angleAT + sina*d2_angleAT, tcross,
810 aux);
811 D2Normal += aux;
812
813/* D2Normal += sina*(D2Tangent.Crossed(Normal) + 2*DTangent.Crossed(DNormal) + Tangent.Crossed(D2Normal)) +
814 2*cosa*d_angleAT*(DTangent.Crossed(Normal) + Tangent.Crossed(DNormal)) +
815 (cosa*d2_angleAT - sina*d_angleAT*d_angleAT)*Tangent.Crossed(Normal) +
8162*sina*d_angleAT*(DTangent.Crossed(Tangent.Crossed(Normal)) + Tangent.Crossed(DTangent.Crossed(Normal)) + Tangent.Crossed(Tangent.Crossed(DNormal))) +
817(1 - cosa)*(D2Tangent.Crossed(Tangent.Crossed(Normal)) + Tangent.Crossed(D2Tangent.Crossed(Normal)) + Tangent.Crossed(Tangent.Crossed(D2Normal)) + 2*DTangent.Crossed(DTangent.Crossed(Normal)) + 2*DTangent.Crossed(Tangent.Crossed(DNormal)) + 2*Tangent.Crossed(DTangent.Crossed(DNormal)))
818+
819(cosa*d_angleAT*d_angleAT + sina*d2_angleAT)*Tangent.Crossed(Tangent.Crossed(Normal));*/
820
821
822 aux.SetLinearForm(sina, dcross,
823 cosa*d_angleAT, cross);
824 aux.SetLinearForm(1 - cosa, dtcross,
825 sina*d_angleAT, tcross,
826 aux);
827 DNormal+=aux;
828
829
830 Normal.SetLinearForm( sina, cross,
831 (1 - cosa), tcross,
832 Normal);
833
834 BiNormal = Tangent.Crossed(Normal);
835
836 DBiNormal.SetLinearForm(1, DTangent.Crossed(Normal),
837 Tangent.Crossed(DNormal));
838
839 D2BiNormal.SetLinearForm(1, D2Tangent.Crossed(Normal),
840 2, DTangent.Crossed(DNormal),
841 Tangent.Crossed(D2Normal));
842
843// for test
844/* gp_Vec FD2N, FD2T, FD2BN, Tf, DTf, Nf, DNf, BNf, DBNf;
845 Standard_Real h;
846 h = 1.0e-8;
847 if (Param + h > myTrimmed->LastParameter()) h = -h;
848 D1(Param + h, Tf, DTf, Nf, DNf, BNf, DBNf);
849 FD2N = (DNf - DNormal)/h;
850 FD2T = (DTf - DTangent)/h;
851 FD2BN = (DBNf - DBiNormal)/h;
852 cout<<"Param = "<<Param<<endl;
853 cout<<"D2N = ("<<D2Normal.X()<<", "<<D2Normal.Y()<<", "<<D2Normal.Z()<<")"<<endl;
854 cout<<"FD2N = ("<<FD2N.X()<<", "<<FD2N.Y()<<", "<<FD2N.Z()<<")"<<endl<<endl;
855 cout<<"D2T = ("<<D2Tangent.X()<<", "<<D2Tangent.Y()<<", "<<D2Tangent.Z()<<")"<<endl;
856 cout<<"FD2T = ("<<FD2T.X()<<", "<<FD2T.Y()<<", "<<FD2T.Z()<<")"<<endl<<endl;
857 cout<<"D2BN = ("<<D2BiNormal.X()<<", "<<D2BiNormal.Y()<<", "<<D2BiNormal.Z()<<")"<<endl;
858 cout<<"FD2BN = ("<<FD2BN.X()<<", "<<FD2BN.Y()<<", "<<FD2BN.Z()<<")"<<endl<<endl;
859*/
860//
861 return Standard_True;
862}
863
864//===============================================================
865// Function : NbIntervals
866// Purpose :
867//===============================================================
868 Standard_Integer GeomFill_CorrectedFrenet::NbIntervals(const GeomAbs_Shape S) const
869{
870 Standard_Integer NbFrenet, NbLaw;
871 NbFrenet = frenet->NbIntervals(S);
872 if (isFrenet) return NbFrenet;
873
874 NbLaw = EvolAroundT->NbIntervals(S);
875 if (NbFrenet == 1)
876 return NbLaw;
877
878 TColStd_Array1OfReal FrenetInt(1, NbFrenet + 1);
879 TColStd_Array1OfReal LawInt(1, NbLaw + 1);
880 TColStd_SequenceOfReal Fusion;
881
882 frenet->Intervals(FrenetInt, S);
883 EvolAroundT->Intervals(LawInt, S);
884 GeomLib::FuseIntervals(FrenetInt, LawInt, Fusion);
885
886 return Fusion.Length()-1;
887}
888
889//===============================================================
890// Function : Intervals
891// Purpose :
892//===============================================================
893 void GeomFill_CorrectedFrenet::Intervals(TColStd_Array1OfReal& T,
894 const GeomAbs_Shape S) const
895{
896 Standard_Integer NbFrenet, NbLaw;
897 if (isFrenet) {
898 frenet->Intervals(T, S);
899 return;
900 }
901
902 NbFrenet = frenet->NbIntervals(S);
903 if(NbFrenet==1) {
904 EvolAroundT->Intervals(T, S);
905 }
906
907 NbLaw = EvolAroundT->NbIntervals(S);
908
909 TColStd_Array1OfReal FrenetInt(1, NbFrenet + 1);
910 TColStd_Array1OfReal LawInt(1, NbLaw + 1);
911 TColStd_SequenceOfReal Fusion;
912
913 frenet->Intervals(FrenetInt, S);
914 EvolAroundT->Intervals(LawInt, S);
915 GeomLib::FuseIntervals(FrenetInt, LawInt, Fusion);
916
917 for(Standard_Integer i = 1; i <= Fusion.Length(); i++)
918 T.ChangeValue(i) = Fusion.Value(i);
919}
920
921//===============================================================
922// Function : SetInterval
923// Purpose :
924//===============================================================
925 void GeomFill_CorrectedFrenet::SetInterval(const Standard_Real First,
926 const Standard_Real Last)
927{
928 GeomFill_TrihedronLaw::SetInterval(First, Last);
929 frenet->SetInterval(First, Last);
930 if (!isFrenet) TLaw = EvolAroundT->Trim(First, Last,
931 Precision::PConfusion()/2);
932}
933
934//===============================================================
a31abc03 935// Function : EvaluateBestMode
936// Purpose :
937//===============================================================
938GeomFill_Trihedron GeomFill_CorrectedFrenet::EvaluateBestMode()
939{
940 if (EvolAroundT.IsNull())
941 return GeomFill_IsFrenet; //Frenet
942
943 const Standard_Real MaxAngle = 3.*M_PI/4.;
944 const Standard_Real MaxTorsion = 100.;
945
946 Standard_Real Step, u, v, tmin, tmax;
947 Standard_Integer NbInt, i, j, k = 1;
948 NbInt = EvolAroundT->NbIntervals(GeomAbs_CN);
949 TColStd_Array1OfReal Int(1, NbInt+1);
950 EvolAroundT->Intervals(Int, GeomAbs_CN);
951 gp_Pnt2d old;
952 gp_Vec2d aVec, PrevVec;
953
954 Standard_Integer NbSamples = 10;
955 for(i = 1; i <= NbInt; i++){
956 tmin = Int(i);
957 tmax = Int(i+1);
958 Standard_Real Torsion = ComputeTorsion(tmin, myTrimmed);
959 if (Abs(Torsion) > MaxTorsion)
960 return GeomFill_IsDiscreteTrihedron; //DiscreteTrihedron
961
962 Handle(Law_Function) trimmedlaw = EvolAroundT->Trim(tmin, tmax, Precision::PConfusion()/2);
963 Step = (Int(i+1)-Int(i))/NbSamples;
964 for (j = 0; j <= NbSamples; j++) {
965 u = tmin + j*Step;
966 v = trimmedlaw->Value(u);
967 gp_Pnt2d point2d(u,v);
968 if (j != 0)
969 {
970 aVec.SetXY(point2d.XY() - old.XY());
971 if (k > 2)
972 {
973 Standard_Real theAngle = PrevVec.Angle(aVec);
974 if (Abs(theAngle) > MaxAngle)
975 return GeomFill_IsDiscreteTrihedron; //DiscreteTrihedron
976 }
977 PrevVec = aVec;
978 }
979 old = point2d;
980 k++;
981 }
982 }
983
984 return GeomFill_IsCorrectedFrenet; //CorrectedFrenet
985}
986
987//===============================================================
7fd59977 988// Function : GetAverageLaw
989// Purpose :
990//===============================================================
991 void GeomFill_CorrectedFrenet::GetAverageLaw(gp_Vec& ATangent,
992 gp_Vec& ANormal,
993 gp_Vec& ABiNormal)
994{
995 if (isFrenet) frenet->GetAverageLaw(ATangent, ANormal, ABiNormal);
996 else {
997 ATangent = AT;
998 ANormal = AN;
999 ABiNormal = ATangent;
1000 ABiNormal.Cross(ANormal);
1001 }
1002}
1003
1004//===============================================================
1005// Function : IsConstant
1006// Purpose :
1007//===============================================================
1008 Standard_Boolean GeomFill_CorrectedFrenet::IsConstant() const
1009{
1010 return (myCurve->GetType() == GeomAbs_Line);
1011}
1012
1013//===============================================================
1014// Function : IsOnlyBy3dCurve
1015// Purpose :
1016//===============================================================
1017 Standard_Boolean GeomFill_CorrectedFrenet::IsOnlyBy3dCurve() const
1018{
1019 return Standard_True;
1020}
Open CASCADE Technology repository