1 // Created on: 1998-05-12
2 // Created by: Roman BORISOV
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 under
9 // the terms of the GNU Lesser General Public License 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.
18 #include <Adaptor2d_HCurve2d.hxx>
19 #include <Adaptor3d_Curve.hxx>
20 #include <Adaptor3d_CurveOnSurface.hxx>
21 #include <Adaptor3d_HCurve.hxx>
22 #include <Adaptor3d_HCurveOnSurface.hxx>
23 #include <Adaptor3d_HSurface.hxx>
24 #include <Approx_CurvlinFunc.hxx>
25 #include <GCPnts_AbscissaPoint.hxx>
26 #include <GeomLib.hxx>
27 #include <Precision.hxx>
28 #include <Standard_ConstructionError.hxx>
29 #include <Standard_DomainError.hxx>
30 #include <Standard_OutOfRange.hxx>
31 #include <Standard_Type.hxx>
32 #include <TColStd_SequenceOfReal.hxx>
34 IMPLEMENT_STANDARD_RTTIEXT(Approx_CurvlinFunc,Standard_Transient)
36 #ifdef OCCT_DEBUG_CHRONO
37 #include <OSD_Timer.hxx>
38 static OSD_Chronometer chr_uparam;
39 Standard_EXPORT Standard_Integer uparam_count;
40 Standard_EXPORT Standard_Real t_uparam;
42 //Standard_IMPORT extern void InitChron(OSD_Chronometer& ch);
43 Standard_IMPORT void InitChron(OSD_Chronometer& ch);
44 //Standard_IMPORT extern void ResultChron( OSD_Chronometer & ch, Standard_Real & time);
45 Standard_IMPORT void ResultChron( OSD_Chronometer & ch, Standard_Real & time);
48 static Standard_Real cubic(const Standard_Real X, const Standard_Real *Xi, const Standard_Real *Yi)
50 Standard_Real I1, I2, I3, I21, I22, I31, Result;
52 I1 = (Yi[0] - Yi[1])/(Xi[0] - Xi[1]);
53 I2 = (Yi[1] - Yi[2])/(Xi[1] - Xi[2]);
54 I3 = (Yi[2] - Yi[3])/(Xi[2] - Xi[3]);
56 I21 = (I1 - I2)/(Xi[0] - Xi[2]);
57 I22 = (I2 - I3)/(Xi[1] - Xi[3]);
59 I31 = (I21 - I22)/(Xi[0] - Xi[3]);
61 Result = Yi[0] + (X - Xi[0])*(I1 + (X - Xi[1])*(I21 + (X - Xi[2])*I31));
66 //static void findfourpoints(const Standard_Real S,
67 static void findfourpoints(const Standard_Real ,
68 Standard_Integer NInterval,
69 const Handle(TColStd_HArray1OfReal)& Si,
70 Handle(TColStd_HArray1OfReal)& Ui,
71 const Standard_Real prevS,
72 const Standard_Real prevU, Standard_Real *Xi,
75 Standard_Integer i, j;
76 Standard_Integer NbInt = Si->Length() - 1;
77 if (NbInt < 3) throw Standard_ConstructionError("Approx_CurvlinFunc::GetUParameter");
79 if(NInterval < 1) NInterval = 1;
80 else if(NInterval > NbInt - 2) NInterval = NbInt - 2;
82 for(i = 0; i < 4; i++) {
83 Xi[i] = Si->Value(NInterval - 1 + i);
84 Yi[i] = Ui->Value(NInterval - 1 + i);
86 // try to insert (S, U)
87 for(i = 0; i < 3; i++) {
88 if(Xi[i] < prevS && prevS < Xi[i+1]) {
89 for(j = 0; j < i; j++) {
100 /*static Standard_Real curvature(const Standard_Real U, const Adaptor3d_Curve& C)
102 Standard_Real k, tau, mod1, mod2, OMEGA;
105 C.D3(U, P, D1, D2, D3);
106 mod1 = D1.Magnitude();
107 mod2 = D1.Crossed(D2).Magnitude();
108 k = mod2/(mod1*mod1*mod1);
109 tau = D1.Dot(D2.Crossed(D3));
111 OMEGA = Sqrt(k*k + tau*tau);
117 Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor3d_HCurve)& C, const Standard_Real Tol) : myC3D(C),
128 Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor2d_HCurve2d)& C2D, const Handle(Adaptor3d_HSurface)& S, const Standard_Real Tol) :
141 Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor2d_HCurve2d)& C2D1, const Handle(Adaptor2d_HCurve2d)& C2D2, const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_HSurface)& S2, const Standard_Real Tol) :
156 void Approx_CurvlinFunc::Init()
158 Adaptor3d_CurveOnSurface CurOnSur;
162 Init(myC3D->GetCurve(), mySi_1, myUi_1);
163 myFirstU1 = myC3D->FirstParameter();
164 myLastU1 = myC3D->LastParameter();
165 myFirstU2 = myLastU2 = 0;
168 CurOnSur.Load(myC2D1);
169 CurOnSur.Load(mySurf1);
170 Init(CurOnSur, mySi_1, myUi_1);
171 myFirstU1 = CurOnSur.FirstParameter();
172 myLastU1 = CurOnSur.LastParameter();
173 myFirstU2 = myLastU2 = 0;
176 CurOnSur.Load(myC2D1);
177 CurOnSur.Load(mySurf1);
178 Init(CurOnSur, mySi_1, myUi_1);
179 myFirstU1 = CurOnSur.FirstParameter();
180 myLastU1 = CurOnSur.LastParameter();
181 CurOnSur.Load(myC2D2);
182 CurOnSur.Load(mySurf2);
183 Init(CurOnSur, mySi_2, myUi_2);
184 myFirstU2 = CurOnSur.FirstParameter();
185 myLastU2 = CurOnSur.LastParameter();
192 //=======================================================================
194 //purpose : Init the values
195 //history : 23/10/1998 PMN : Cut at curve's discontinuities
196 //=======================================================================
197 void Approx_CurvlinFunc::Init(Adaptor3d_Curve& C, Handle(TColStd_HArray1OfReal)& Si,
198 Handle(TColStd_HArray1OfReal)& Ui) const
200 Standard_Real Step, FirstU, LastU;
201 Standard_Integer i, j, k, NbInt, NbIntC3;
202 FirstU = C.FirstParameter();
203 LastU = C.LastParameter();
206 NbIntC3 = C.NbIntervals(GeomAbs_C3);
207 TColStd_Array1OfReal Disc(1, NbIntC3+1);
210 C.Intervals(Disc, GeomAbs_C3);
217 Ui = new TColStd_HArray1OfReal (0,NbIntC3*NbInt);
218 Si = new TColStd_HArray1OfReal (0,NbIntC3*NbInt);
220 Ui->SetValue(0, FirstU);
223 for(j = 1, i=1; j<=NbIntC3; j++) {
224 Step = (Disc(j+1) - Disc(j))/NbInt;
225 for(k = 1; k <= NbInt; k++, i++) {
226 Ui->ChangeValue(i) = Ui->Value(i-1) + Step;
227 Si->ChangeValue(i) = Si->Value(i-1) + Length(C, Ui->Value(i-1), Ui->Value(i));
231 Standard_Real Len = Si->Value(Si->Upper());
232 for(i = Si->Lower(); i<= Si->Upper(); i++)
233 Si->ChangeValue(i) /= Len;
235 // TODO - fields should be mutable
236 const_cast<Approx_CurvlinFunc*>(this)->myPrevS = myFirstS;
237 const_cast<Approx_CurvlinFunc*>(this)->myPrevU = FirstU;
240 void Approx_CurvlinFunc::SetTol(const Standard_Real Tol)
245 Standard_Real Approx_CurvlinFunc::FirstParameter() const
250 Standard_Real Approx_CurvlinFunc::LastParameter() const
255 Standard_Integer Approx_CurvlinFunc::NbIntervals(const GeomAbs_Shape S) const
257 Adaptor3d_CurveOnSurface CurOnSur;
261 return myC3D->NbIntervals(S);
263 CurOnSur.Load(myC2D1);
264 CurOnSur.Load(mySurf1);
265 return CurOnSur.NbIntervals(S);
267 Standard_Integer NbInt;
268 CurOnSur.Load(myC2D1);
269 CurOnSur.Load(mySurf1);
270 NbInt = CurOnSur.NbIntervals(S);
271 TColStd_Array1OfReal T1(1, NbInt+1);
272 CurOnSur.Intervals(T1, S);
273 CurOnSur.Load(myC2D2);
274 CurOnSur.Load(mySurf2);
275 NbInt = CurOnSur.NbIntervals(S);
276 TColStd_Array1OfReal T2(1, NbInt+1);
277 CurOnSur.Intervals(T2, S);
279 TColStd_SequenceOfReal Fusion;
280 GeomLib::FuseIntervals(T1, T2, Fusion);
281 return Fusion.Length() - 1;
288 void Approx_CurvlinFunc::Intervals(TColStd_Array1OfReal& T, const GeomAbs_Shape S) const
290 Adaptor3d_CurveOnSurface CurOnSur;
295 myC3D->Intervals(T, S);
298 CurOnSur.Load(myC2D1);
299 CurOnSur.Load(mySurf1);
300 CurOnSur.Intervals(T, S);
303 Standard_Integer NbInt;
304 CurOnSur.Load(myC2D1);
305 CurOnSur.Load(mySurf1);
306 NbInt = CurOnSur.NbIntervals(S);
307 TColStd_Array1OfReal T1(1, NbInt+1);
308 CurOnSur.Intervals(T1, S);
309 CurOnSur.Load(myC2D2);
310 CurOnSur.Load(mySurf2);
311 NbInt = CurOnSur.NbIntervals(S);
312 TColStd_Array1OfReal T2(1, NbInt+1);
313 CurOnSur.Intervals(T2, S);
315 TColStd_SequenceOfReal Fusion;
316 GeomLib::FuseIntervals(T1, T2, Fusion);
318 for (i = 1; i <= Fusion.Length(); i++)
319 T.ChangeValue(i) = Fusion.Value(i);
322 for(i = 1; i <= T.Length(); i++)
323 T.ChangeValue(i) = GetSParameter(T.Value(i));
326 void Approx_CurvlinFunc::Trim(const Standard_Real First, const Standard_Real Last, const Standard_Real Tol)
328 if (First < 0 || Last >1) throw Standard_OutOfRange("Approx_CurvlinFunc::Trim");
329 if ((Last - First) < Tol) return;
331 Standard_Real FirstU, LastU;
332 Adaptor3d_CurveOnSurface CurOnSur;
333 Handle(Adaptor3d_HCurve) HCurOnSur;
337 myC3D = myC3D->Trim(myFirstU1, myLastU1, Tol);
338 FirstU = GetUParameter(myC3D->GetCurve(), First, 1);
339 LastU = GetUParameter(myC3D->GetCurve(), Last, 1);
340 myC3D = myC3D->Trim(FirstU, LastU, Tol);
343 CurOnSur.Load(myC2D2);
344 CurOnSur.Load(mySurf2);
345 HCurOnSur = CurOnSur.Trim(myFirstU2, myLastU2, Tol);
346 myC2D2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
347 mySurf2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
348 CurOnSur.Load(myC2D2);
349 CurOnSur.Load(mySurf2);
351 FirstU = GetUParameter(CurOnSur, First, 1);
352 LastU = GetUParameter(CurOnSur, Last, 1);
353 HCurOnSur = CurOnSur.Trim(FirstU, LastU, Tol);
354 myC2D2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
355 mySurf2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
359 CurOnSur.Load(myC2D1);
360 CurOnSur.Load(mySurf1);
361 HCurOnSur = CurOnSur.Trim(myFirstU1, myLastU1, Tol);
362 myC2D1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
363 mySurf1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
364 CurOnSur.Load(myC2D1);
365 CurOnSur.Load(mySurf1);
367 FirstU = GetUParameter(CurOnSur, First, 1);
368 LastU = GetUParameter(CurOnSur, Last, 1);
369 HCurOnSur = CurOnSur.Trim(FirstU, LastU, Tol);
370 myC2D1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
371 mySurf1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
377 void Approx_CurvlinFunc::Length()
379 Adaptor3d_CurveOnSurface CurOnSur;
380 Standard_Real FirstU, LastU;
384 FirstU = myC3D->FirstParameter();
385 LastU = myC3D->LastParameter();
386 myLength = Length(myC3D->GetCurve(), FirstU, LastU);
387 myLength1 = myLength2 = 0;
390 CurOnSur.Load(myC2D1);
391 CurOnSur.Load(mySurf1);
392 FirstU = CurOnSur.FirstParameter();
393 LastU = CurOnSur.LastParameter();
394 myLength = Length(CurOnSur, FirstU, LastU);
395 myLength1 = myLength2 = 0;
398 CurOnSur.Load(myC2D1);
399 CurOnSur.Load(mySurf1);
400 FirstU = CurOnSur.FirstParameter();
401 LastU = CurOnSur.LastParameter();
402 myLength1 = Length(CurOnSur, FirstU, LastU);
403 CurOnSur.Load(myC2D2);
404 CurOnSur.Load(mySurf2);
405 FirstU = CurOnSur.FirstParameter();
406 LastU = CurOnSur.LastParameter();
407 myLength2 = Length(CurOnSur, FirstU, LastU);
408 myLength = (myLength1 + myLength2)/2;
413 Standard_Real Approx_CurvlinFunc::Length(Adaptor3d_Curve& C, const Standard_Real FirstU, const Standard_Real LastU) const
415 Standard_Real Length;
417 Length = GCPnts_AbscissaPoint::Length(C, FirstU, LastU, myTolLen);
422 Standard_Real Approx_CurvlinFunc::GetLength() const
427 Standard_Real Approx_CurvlinFunc::GetSParameter(const Standard_Real U) const
429 Standard_Real S=0, S1, S2;
430 Adaptor3d_CurveOnSurface CurOnSur;
434 S = GetSParameter(myC3D->GetCurve(), U, myLength);
437 CurOnSur.Load(myC2D1);
438 CurOnSur.Load(mySurf1);
439 S = GetSParameter(CurOnSur, U, myLength);
442 CurOnSur.Load(myC2D1);
443 CurOnSur.Load(mySurf1);
444 S1 = GetSParameter(CurOnSur, U, myLength1);
445 CurOnSur.Load(myC2D2);
446 CurOnSur.Load(mySurf2);
447 S2 = GetSParameter(CurOnSur, U, myLength2);
455 Standard_Real Approx_CurvlinFunc::GetUParameter(Adaptor3d_Curve& C,
456 const Standard_Real S,
457 const Standard_Integer NumberOfCurve) const
459 Standard_Real deltaS, base, U, Length;
460 Standard_Integer NbInt, NInterval, i;
461 Handle(TColStd_HArray1OfReal) InitUArray, InitSArray;
462 #ifdef OCCT_DEBUG_CHRONO
463 InitChron(chr_uparam);
465 if(S < 0 || S > 1) throw Standard_ConstructionError("Approx_CurvlinFunc::GetUParameter");
467 if(NumberOfCurve == 1) {
481 NbInt = InitUArray->Length() - 1;
483 if(S == 1) NInterval = NbInt - 1;
485 for(i = 0; i < NbInt; i++) {
486 if((InitSArray->Value(i) <= S && S < InitSArray->Value(i+1)))
491 if(S==InitSArray->Value(NInterval)) {
492 return InitUArray->Value(NInterval);
494 if(S==InitSArray->Value(NInterval+1)) {
495 return InitUArray->Value(NInterval+1);
498 base = InitUArray->Value(NInterval);
499 deltaS = (S - InitSArray->Value(NInterval))*Length;
501 // to find an initial point
502 Standard_Real Xi[4], Yi[4], UGuess;
503 findfourpoints(S, NInterval, InitSArray, InitUArray, myPrevS, myPrevU, Xi, Yi);
504 UGuess = cubic(S , Xi, Yi);
506 U = GCPnts_AbscissaPoint(C, deltaS, base, UGuess, myTolLen).Parameter();
508 // TODO - fields should be mutable
509 const_cast<Approx_CurvlinFunc*>(this)->myPrevS = S;
510 const_cast<Approx_CurvlinFunc*>(this)->myPrevU = U;
512 #ifdef OCCT_DEBUG_CHRONO
513 ResultChron(chr_uparam, t_uparam);
520 Standard_Real Approx_CurvlinFunc::GetSParameter(Adaptor3d_Curve& C, const Standard_Real U, const Standard_Real Len) const
522 Standard_Real S, Origin;
524 Origin = C.FirstParameter();
525 S = myFirstS + Length(C, Origin, U)/Len;
529 Standard_Boolean Approx_CurvlinFunc::EvalCase1(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) const
531 if(myCase != 1) throw Standard_ConstructionError("Approx_CurvlinFunc::EvalCase1");
534 gp_Vec dC_dU, dC_dS, d2C_dU2, d2C_dS2;
535 Standard_Real U, Mag, dU_dS, d2U_dS2;
537 U = GetUParameter(myC3D->GetCurve(), S, 1);
550 myC3D->D1(U, C, dC_dU);
551 Mag = dC_dU.Magnitude();
552 dU_dS = myLength/Mag;
555 Result(0) = dC_dS.X();
556 Result(1) = dC_dS.Y();
557 Result(2) = dC_dS.Z();
561 myC3D->D2(U, C, dC_dU, d2C_dU2);
562 Mag = dC_dU.Magnitude();
563 dU_dS = myLength/Mag;
564 d2U_dS2 = -myLength*dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag*Mag);
565 d2C_dS2 = d2C_dU2*dU_dS*dU_dS + dC_dU*d2U_dS2;
567 Result(0) = d2C_dS2.X();
568 Result(1) = d2C_dS2.Y();
569 Result(2) = d2C_dS2.Z();
572 default: Result(0) = Result(1) = Result(2) = 0;
573 return Standard_False;
575 return Standard_True;
578 Standard_Boolean Approx_CurvlinFunc::EvalCase2(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) const
580 if(myCase != 2) throw Standard_ConstructionError("Approx_CurvlinFunc::EvalCase2");
582 Standard_Boolean Done;
584 Done = EvalCurOnSur(S, Order, Result, 1);
589 Standard_Boolean Approx_CurvlinFunc::EvalCase3(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result)
591 if(myCase != 3) throw Standard_ConstructionError("Approx_CurvlinFunc::EvalCase3");
593 TColStd_Array1OfReal tmpRes1(0, 4), tmpRes2(0, 4);
594 Standard_Boolean Done;
596 Done = EvalCurOnSur(S, Order, tmpRes1, 1);
598 Done = EvalCurOnSur(S, Order, tmpRes2, 2) && Done;
600 Result(0) = tmpRes1(0);
601 Result(1) = tmpRes1(1);
602 Result(2) = tmpRes2(0);
603 Result(3) = tmpRes2(1);
604 Result(4) = 0.5*(tmpRes1(2) + tmpRes2(2));
605 Result(5) = 0.5*(tmpRes1(3) + tmpRes2(3));
606 Result(6) = 0.5*(tmpRes1(4) + tmpRes2(4));
611 Standard_Boolean Approx_CurvlinFunc::EvalCurOnSur(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result, const Standard_Integer NumberOfCurve) const
613 Handle(Adaptor2d_HCurve2d) Cur2D;
614 Handle(Adaptor3d_HSurface) Surf;
615 Standard_Real U=0, Length=0;
617 if (NumberOfCurve == 1) {
620 Adaptor3d_CurveOnSurface CurOnSur(myC2D1, mySurf1);
621 U = GetUParameter(CurOnSur, S, 1);
622 if(myCase == 3) Length = myLength1;
623 else Length = myLength;
625 else if (NumberOfCurve == 2) {
628 Adaptor3d_CurveOnSurface CurOnSur(myC2D2, mySurf2);
629 U = GetUParameter(CurOnSur, S, 2);
633 throw Standard_ConstructionError("Approx_CurvlinFunc::EvalCurOnSur");
635 Standard_Real Mag, dU_dS, d2U_dS2, dV_dU, dW_dU, dV_dS, dW_dS, d2V_dS2, d2W_dS2, d2V_dU2, d2W_dU2;
638 gp_Vec2d dC2D_dU, d2C2D_dU2;
639 gp_Vec dC_dU, d2C_dU2, dC_dS, d2C_dS2, dS_dV, dS_dW, d2S_dV2, d2S_dW2, d2S_dVdW;
644 Surf->D0(C2D.X(), C2D.Y(), C);
654 Cur2D->D1(U, C2D, dC2D_dU);
657 Surf->D1(C2D.X(), C2D.Y(), C, dS_dV, dS_dW);
658 dC_dU = dS_dV*dV_dU + dS_dW*dW_dU;
659 Mag = dC_dU.Magnitude();
668 Result(2) = dC_dS.X();
669 Result(3) = dC_dS.Y();
670 Result(4) = dC_dS.Z();
674 Cur2D->D2(U, C2D, dC2D_dU, d2C2D_dU2);
677 d2V_dU2 = d2C2D_dU2.X();
678 d2W_dU2 = d2C2D_dU2.Y();
679 Surf->D2(C2D.X(), C2D.Y(), C, dS_dV, dS_dW, d2S_dV2, d2S_dW2, d2S_dVdW);
680 dC_dU = dS_dV*dV_dU + dS_dW*dW_dU;
681 d2C_dU2 = (d2S_dV2*dV_dU + d2S_dVdW*dW_dU)*dV_dU + dS_dV*d2V_dU2 +
682 (d2S_dVdW*dV_dU + d2S_dW2*dW_dU)*dW_dU + dS_dW*d2W_dU2;
683 Mag = dC_dU.Magnitude();
685 d2U_dS2 = -Length*dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag*Mag);
687 dV_dS = dV_dU * dU_dS;
688 dW_dS = dW_dU * dU_dS;
689 d2V_dS2 = d2V_dU2*dU_dS*dU_dS + dV_dU*d2U_dS2;
690 d2W_dS2 = d2W_dU2*dU_dS*dU_dS + dW_dU*d2U_dS2;
692 d2U_dS2 = -dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag);
693 d2C_dS2 = (d2S_dV2 * dV_dS + d2S_dVdW * dW_dS) * dV_dS + dS_dV * d2V_dS2 +
694 (d2S_dW2 * dW_dS + d2S_dVdW * dV_dS) * dW_dS + dS_dW * d2W_dS2;
698 Result(2) = d2C_dS2.X();
699 Result(3) = d2C_dS2.Y();
700 Result(4) = d2C_dS2.Z();
703 default: Result(0) = Result(1) = Result(2) = Result(3) = Result(4) = 0;
704 return Standard_False;
706 return Standard_True;