1 // Created on: 1994-08-18
2 // Created by: Laurent PAINNOT
3 // Copyright (c) 1994-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.
17 // 8-Aug-95 : xab : interpolation uses BSplCLib::Interpolate
19 #include <BSplCLib.hxx>
20 #include <Geom2d_BSplineCurve.hxx>
21 #include <Geom2dAPI_Interpolate.hxx>
22 #include <gp_Vec2d.hxx>
24 #include <Standard_ConstructionError.hxx>
25 #include <StdFail_NotDone.hxx>
26 #include <TColgp_Array1OfPnt2d.hxx>
27 #include <TColStd_Array1OfBoolean.hxx>
28 #include <TColStd_Array1OfInteger.hxx>
29 #include <TColStd_HArray1OfBoolean.hxx>
30 #include <TColStd_HArray1OfReal.hxx>
32 //=======================================================================
33 //function : CheckPoints
35 //=======================================================================
36 static Standard_Boolean CheckPoints(const TColgp_Array1OfPnt2d& PointArray,
37 const Standard_Real Tolerance)
40 Standard_Real tolerance_squared = Tolerance * Tolerance,
42 Standard_Boolean result = Standard_True ;
43 for (ii = PointArray.Lower() ; result && ii < PointArray.Upper() ; ii++) {
45 PointArray.Value(ii).SquareDistance(PointArray.Value(ii+1)) ;
46 result = (distance_squared >= tolerance_squared) ;
51 //=======================================================================
52 //function : CheckTangents
54 //=======================================================================
55 static Standard_Boolean CheckTangents(
56 const TColgp_Array1OfVec2d& Tangents,
57 const TColStd_Array1OfBoolean& TangentFlags,
58 const Standard_Real Tolerance)
62 Standard_Real tolerance_squared = Tolerance * Tolerance,
64 Standard_Boolean result = Standard_True ;
65 index = TangentFlags.Lower() ;
66 for (ii = Tangents.Lower(); result && ii <= Tangents.Upper() ; ii++) {
67 if(TangentFlags.Value(index)) {
69 Tangents.Value(ii).SquareMagnitude() ;
70 result = (distance_squared >= tolerance_squared) ;
77 //=======================================================================
78 //function : CheckParameters
80 //=======================================================================
81 static Standard_Boolean CheckParameters(const
82 TColStd_Array1OfReal& Parameters)
85 Standard_Real distance ;
86 Standard_Boolean result = Standard_True ;
87 for (ii = Parameters.Lower() ; result && ii < Parameters.Upper() ; ii++) {
89 Parameters.Value(ii+1) - Parameters.Value(ii) ;
90 result = (distance >= RealSmall()) ;
94 //=======================================================================
95 //function : BuildParameters
97 //=======================================================================
98 static void BuildParameters(const Standard_Boolean PeriodicFlag,
99 const TColgp_Array1OfPnt2d& PointsArray,
100 Handle(TColStd_HArray1OfReal)& ParametersPtr)
104 Standard_Real distance ;
106 num_parameters = PointsArray.Length() ;
108 num_parameters += 1 ;
111 new TColStd_HArray1OfReal(1,
113 ParametersPtr->SetValue(1,0.0e0) ;
115 for (ii = PointsArray.Lower() ; ii < PointsArray.Upper() ; ii++) {
117 PointsArray.Value(ii).Distance(PointsArray.Value(ii+1)) ;
118 ParametersPtr->SetValue(index,
119 ParametersPtr->Value(ii) + distance) ;
124 PointsArray.Value(PointsArray.Upper()).
125 Distance(PointsArray.Value(PointsArray.Lower())) ;
126 ParametersPtr->SetValue(index,
127 ParametersPtr->Value(ii) + distance) ;
130 //=======================================================================
131 //function : BuildPeriodicTangents
133 //=======================================================================
135 static void BuildPeriodicTangent(
136 const TColgp_Array1OfPnt2d& PointsArray,
137 TColgp_Array1OfVec2d& TangentsArray,
138 TColStd_Array1OfBoolean& TangentFlags,
139 const TColStd_Array1OfReal& ParametersArray)
144 Standard_Real *point_array,
150 if (PointsArray.Length() < 3) {
151 throw Standard_ConstructionError();
154 if (!TangentFlags.Value(1)) {
156 if (PointsArray.Length() == 3) {
159 point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ;
161 (Standard_Real *) &ParametersArray.Value(1) ;
162 TangentFlags.SetValue(1,Standard_True) ;
163 PLib::EvalLagrange(ParametersArray.Value(1),
170 for (ii = 1 ; ii <= 2 ; ii++) {
171 a_vector.SetCoord(ii,eval_result[1][ii-1]) ;
173 TangentsArray.SetValue(1,a_vector) ;
176 //=======================================================================
177 //function : BuildTangents
179 //=======================================================================
181 static void BuildTangents(const TColgp_Array1OfPnt2d& PointsArray,
182 TColgp_Array1OfVec2d& TangentsArray,
183 TColStd_Array1OfBoolean& TangentFlags,
184 const TColStd_Array1OfReal& ParametersArray)
188 Standard_Real *point_array,
196 if ( PointsArray.Length() < 3) {
197 throw Standard_ConstructionError();
199 if (PointsArray.Length() == 3) {
202 if (!TangentFlags.Value(1)) {
203 point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ;
205 (Standard_Real *) &ParametersArray.Value(1) ;
206 TangentFlags.SetValue(1,Standard_True) ;
207 PLib::EvalLagrange(ParametersArray.Value(1),
214 for (ii = 1 ; ii <= 2 ; ii++) {
215 a_vector.SetCoord(ii,eval_result[1][ii-1]) ;
217 TangentsArray.SetValue(1,a_vector) ;
219 if (! TangentFlags.Value(TangentFlags.Upper())) {
221 (Standard_Real *) &PointsArray.Value(PointsArray.Upper() - degree) ;
222 TangentFlags.SetValue(TangentFlags.Upper(),Standard_True) ;
224 (Standard_Real *)&ParametersArray.Value(ParametersArray.Upper() - degree) ;
225 PLib::EvalLagrange(ParametersArray.Value(ParametersArray.Upper()),
232 for (ii = 1 ; ii <= 2 ; ii++) {
233 a_vector.SetCoord(ii,eval_result[1][ii-1]) ;
235 TangentsArray.SetValue(TangentsArray.Upper(),a_vector) ;
238 //=======================================================================
239 //function : BuildTangents
240 //purpose : scale the given tangent so that they have the length of
241 // the size of the derivative of the lagrange interpolation
243 //=======================================================================
244 static void ScaleTangents(const TColgp_Array1OfPnt2d& PointsArray,
245 TColgp_Array1OfVec2d& TangentsArray,
246 const TColStd_Array1OfBoolean& TangentFlags,
247 const TColStd_Array1OfReal& ParametersArray)
255 Standard_Real *point_array,
263 num_points = PointsArray.Length() ;
264 if (num_points == 2) {
267 else if (num_points >= 3) {
271 index = PointsArray.Lower() ;
272 for (ii = TangentFlags.Lower() ; ii <= TangentFlags.Upper() ; ii++) {
273 if (TangentFlags.Value(ii)) {
275 (Standard_Real *) &PointsArray.Value(index) ;
277 (Standard_Real *) &ParametersArray.Value(index) ;
278 PLib::EvalLagrange(ParametersArray.Value(ii),
287 for (jj = 1 ; jj <= 2 ; jj++) {
288 value[0] += Abs(TangentsArray.Value(ii).Coord(jj)) ;
289 value[1] += Abs(eval_result[1][jj-1]) ;
291 ratio = value[1] / value[0] ;
292 for (jj = 1 ; jj <= 2 ; jj++) {
293 a_vector.SetCoord(jj, ratio *
294 TangentsArray.Value(ii).Coord(jj)) ;
296 TangentsArray.SetValue(ii, a_vector) ;
297 if (ii != TangentFlags.Lower()) {
300 if (index > PointsArray.Upper() - degree) {
301 index = PointsArray.Upper() - degree ;
308 //=======================================================================
309 //function : Geom2dAPI_Interpolate
311 //=======================================================================
313 Geom2dAPI_Interpolate::Geom2dAPI_Interpolate
314 (const Handle(TColgp_HArray1OfPnt2d)& PointsPtr,
315 const Standard_Boolean PeriodicFlag,
316 const Standard_Real Tolerance) :
317 myTolerance(Tolerance),
319 myIsDone(Standard_False),
320 myPeriodic(PeriodicFlag),
321 myTangentRequest(Standard_False)
324 Standard_Integer ii ;
325 Standard_Boolean result =
326 CheckPoints(PointsPtr->Array1(),
329 new TColgp_HArray1OfVec2d(myPoints->Lower(),
332 new TColStd_HArray1OfBoolean(myPoints->Lower(),
336 throw Standard_ConstructionError();
338 BuildParameters(PeriodicFlag,
342 for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) {
343 myTangentFlags->SetValue(ii,Standard_False) ;
350 //=======================================================================
351 //function : Geom2dAPI_Interpolate
353 //=======================================================================
355 Geom2dAPI_Interpolate::Geom2dAPI_Interpolate
356 (const Handle(TColgp_HArray1OfPnt2d)& PointsPtr,
357 const Handle(TColStd_HArray1OfReal)& ParametersPtr,
358 const Standard_Boolean PeriodicFlag,
359 const Standard_Real Tolerance) :
360 myTolerance(Tolerance),
362 myIsDone(Standard_False),
363 myParameters(ParametersPtr),
364 myPeriodic(PeriodicFlag),
365 myTangentRequest(Standard_False)
367 Standard_Integer ii ;
370 Standard_Boolean result =
371 CheckPoints(PointsPtr->Array1(),
375 if ((PointsPtr->Length()) + 1 != ParametersPtr->Length()) {
376 throw Standard_ConstructionError();
380 new TColgp_HArray1OfVec2d(myPoints->Lower(),
383 new TColStd_HArray1OfBoolean(myPoints->Lower(),
387 throw Standard_ConstructionError();
391 CheckParameters(ParametersPtr->Array1()) ;
393 throw Standard_ConstructionError();
396 for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) {
397 myTangentFlags->SetValue(ii,Standard_False) ;
401 //=======================================================================
404 //=======================================================================
406 void Geom2dAPI_Interpolate::Load(
407 const TColgp_Array1OfVec2d& Tangents,
408 const Handle(TColStd_HArray1OfBoolean)& TangentFlagsPtr)
411 Standard_Boolean result ;
412 Standard_Integer ii ;
413 myTangentRequest = Standard_True ;
414 myTangentFlags = TangentFlagsPtr ;
415 if (Tangents.Length() != myPoints->Length() ||
416 TangentFlagsPtr->Length() != myPoints->Length()) {
417 throw Standard_ConstructionError();
420 CheckTangents(Tangents,
421 TangentFlagsPtr->Array1(),
425 new TColgp_HArray1OfVec2d(Tangents.Lower(),Tangents.Upper()) ;
426 for (ii = Tangents.Lower() ; ii <= Tangents.Upper() ; ii++ ) {
427 myTangents->SetValue(ii,Tangents.Value(ii)) ;
429 ScaleTangents(myPoints->Array1(),
430 myTangents->ChangeArray1(),
431 TangentFlagsPtr->Array1(),
432 myParameters->Array1()) ;
435 throw Standard_ConstructionError();
441 //=======================================================================
444 //=======================================================================
446 void Geom2dAPI_Interpolate::Load(const gp_Vec2d& InitialTangent,
447 const gp_Vec2d& FinalTangent)
449 Standard_Boolean result ;
450 myTangentRequest = Standard_True ;
451 myTangentFlags->SetValue(1,Standard_True) ;
452 myTangentFlags->SetValue(myPoints->Length(),Standard_True) ;
453 myTangents->SetValue(1,InitialTangent) ;
454 myTangents->SetValue(myPoints->Length(),FinalTangent);
456 CheckTangents(myTangents->Array1(),
457 myTangentFlags->Array1(),
460 throw Standard_ConstructionError();
462 ScaleTangents(myPoints->Array1(),
463 myTangents->ChangeArray1(),
464 myTangentFlags->Array1(),
465 myParameters->Array1()) ;
468 //=======================================================================
471 //=======================================================================
473 void Geom2dAPI_Interpolate::Perform()
479 PerformNonPeriodic() ;
482 //=======================================================================
483 //function : PerformPeriodic
485 //=======================================================================
487 void Geom2dAPI_Interpolate::PerformPeriodic()
489 Standard_Integer degree,
502 Standard_Real period ;
506 num_points = myPoints->Length() ;
507 period = myParameters->Value(myParameters->Upper()) -
508 myParameters->Value(myParameters->Lower()) ;
509 num_poles = num_points + 1 ;
510 if (num_points == 2 && !myTangentRequest) {
512 // build a periodic curve of degree 1
516 TColStd_Array1OfInteger deg1_mults(1,num_poles) ;
517 for (ii = 1 ; ii <= num_poles ; ii++) {
518 deg1_mults.SetValue(ii,1) ;
522 new Geom2d_BSplineCurve(myPoints->Array1(),
523 myParameters->Array1(),
527 myIsDone = Standard_True ;
531 num_distinct_knots = num_points + 1 ;
535 if (myTangentRequest)
536 for (ii = myTangentFlags->Lower() + 1 ;
537 ii <= myTangentFlags->Upper() ; ii++) {
538 if (myTangentFlags->Value(ii)) {
543 TColStd_Array1OfReal parameters(1,num_poles) ;
544 TColStd_Array1OfReal flatknots(1,num_poles + degree + 1) ;
545 TColStd_Array1OfInteger mults(1,num_distinct_knots) ;
546 TColStd_Array1OfInteger contact_order_array(1, num_poles) ;
547 TColgp_Array1OfPnt2d poles(1,num_poles) ;
549 for (ii = 1 ; ii <= half_order ; ii++) {
550 flatknots.SetValue(ii,myParameters->Value(myParameters->Upper() -1) -
552 flatknots.SetValue(ii + half_order,myParameters->
553 Value(myParameters->Lower())) ;
554 flatknots.SetValue(num_poles + ii,
555 myParameters->Value(myParameters->Upper())) ;
556 flatknots.SetValue(num_poles + half_order + ii,
557 myParameters->Value(half_order) + period) ;
559 for (ii = 1 ; ii <= num_poles ; ii++) {
560 contact_order_array.SetValue(ii,0) ;
562 for (ii = 2 ; ii < num_distinct_knots ; ii++) {
563 mults.SetValue(ii,1) ;
565 mults.SetValue(1,half_order) ;
566 mults.SetValue(num_distinct_knots ,half_order) ;
567 if (num_points >= 3) {
570 // only enter here if there are more than 3 points otherwise
571 // it means we have already the tangent
573 BuildPeriodicTangent(myPoints->Array1(),
574 myTangents->ChangeArray1(),
575 myTangentFlags->ChangeArray1(),
576 myParameters->Array1()) ;
578 contact_order_array.SetValue(2,1) ;
579 parameters.SetValue(1,myParameters->Value(1)) ;
580 parameters.SetValue(2,myParameters->Value(1)) ;
581 poles.SetValue(1,myPoints->Value(1)) ;
582 for (jj = 1 ; jj <= 2 ; jj++) {
583 a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
585 poles.SetValue(2,a_point) ;
589 index1 = degree + 2 ;
590 if (myTangentRequest) {
591 for (ii = myTangentFlags->Lower() + 1 ;
592 ii <= myTangentFlags->Upper() ; ii++) {
593 parameters.SetValue(index,myParameters->Value(ii)) ;
594 flatknots.SetValue(index1,myParameters->Value(ii)) ;
595 poles.SetValue(index,myPoints->Value(ii)) ;
598 if (myTangentFlags->Value(ii)) {
599 mults.SetValue(mult_index,mults.Value(mult_index) + 1) ;
600 contact_order_array(index) = 1 ;
602 parameters.SetValue(index,
603 myParameters->Value(ii)) ;
604 flatknots.SetValue(index1,myParameters->Value(ii)) ;
605 for (jj = 1 ; jj <= 2 ; jj++) {
606 a_point.SetCoord(jj,myTangents->Value(ii).Coord(jj)) ;
608 poles.SetValue(index,a_point) ;
618 for(ii = myParameters->Lower() ; ii <= myParameters->Upper() ; ii++) {
619 parameters.SetValue(index1,
620 myParameters->Value(ii)) ;
621 flatknots.SetValue(index,
622 myParameters->Value(ii)) ;
627 for (ii = myPoints->Lower() + 1 ; ii <= myPoints->Upper() ; ii++) {
629 // copy all the given points since the last one will be initialized
630 // below by the first point in the array myPoints
632 poles.SetValue(index,
633 myPoints->Value(ii)) ;
638 contact_order_array.SetValue(num_poles - 1, 1) ;
639 parameters.SetValue(num_poles-1,
640 myParameters->Value(myParameters->Upper())) ;
642 // for the periodic curve ONLY the tangent of the first point
643 // will be used since the curve should close itself at the first
644 // point See BuildPeriodicTangent
646 for (jj = 1 ; jj <= 2 ; jj++) {
647 a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
649 poles.SetValue(num_poles-1,a_point) ;
651 parameters.SetValue(num_poles,
652 myParameters->Value(myParameters->Upper())) ;
654 poles.SetValue(num_poles,
655 myPoints->Value(1)) ;
658 BSplCLib::Interpolate(degree,
664 if (!inversion_problem) {
665 TColgp_Array1OfPnt2d newpoles(poles.Value(1),
669 new Geom2d_BSplineCurve(newpoles,
670 myParameters->Array1(),
674 myIsDone = Standard_True ;
680 //=======================================================================
681 //function : PerformNonPeriodic
683 //=======================================================================
685 void Geom2dAPI_Interpolate::PerformNonPeriodic()
687 Standard_Integer degree,
704 num_poles = myPoints->Length() ;
705 if (num_poles == 2 && !myTangentRequest) {
708 else if (num_poles == 3 && !myTangentRequest) {
710 num_distinct_knots = 2 ;
715 if (myTangentRequest)
716 for (ii = myTangentFlags->Lower() + 1 ;
717 ii < myTangentFlags->Upper() ; ii++) {
718 if (myTangentFlags->Value(ii)) {
725 TColStd_Array1OfReal parameters(1,num_poles) ;
726 TColStd_Array1OfReal flatknots(1,num_poles + degree + 1) ;
727 TColStd_Array1OfInteger mults(1,num_distinct_knots) ;
728 TColStd_Array1OfReal knots(1,num_distinct_knots) ;
729 TColStd_Array1OfInteger contact_order_array(1, num_poles) ;
730 TColgp_Array1OfPnt2d poles(1,num_poles) ;
732 for (ii = 1 ; ii <= degree + 1 ; ii++) {
733 flatknots.SetValue(ii,myParameters->Value(1)) ;
734 flatknots.SetValue(ii + num_poles,
735 myParameters->Value(num_points)) ;
737 for (ii = 1 ; ii <= num_poles ; ii++) {
738 contact_order_array.SetValue(ii,0) ;
740 for (ii = 2 ; ii < num_distinct_knots ; ii++) {
741 mults.SetValue(ii,1) ;
743 mults.SetValue(1,degree + 1) ;
744 mults.SetValue(num_distinct_knots ,degree + 1) ;
748 for (ii = 1 ; ii <= num_poles ; ii++) {
749 poles.SetValue(ii ,myPoints->Value(ii)) ;
752 new Geom2d_BSplineCurve(poles,
753 myParameters->Array1(),
756 myIsDone = Standard_True ;
759 knots.SetValue(1,myParameters->Value(1)) ;
760 knots.SetValue(2,myParameters->Value(3)) ;
761 for (ii = 1 ; ii <= num_poles ; ii++) {
762 poles.SetValue(ii,myPoints->Value(ii)) ;
765 BSplCLib::Interpolate(degree,
767 myParameters->Array1(),
771 if (!inversion_problem) {
773 new Geom2d_BSplineCurve(poles,
777 myIsDone = Standard_True ;
782 // check if the boundary conditions are set
784 if (num_points >= 3) {
786 // cannot build the tangents with degree 3 with only 2 points
787 // if those where not given in advance
789 BuildTangents(myPoints->Array1(),
790 myTangents->ChangeArray1(),
791 myTangentFlags->ChangeArray1(),
792 myParameters->Array1()) ;
794 contact_order_array.SetValue(2,1) ;
795 parameters.SetValue(1,myParameters->Value(1)) ;
796 parameters.SetValue(2,myParameters->Value(1)) ;
797 poles.SetValue(1,myPoints->Value(1)) ;
798 for (jj = 1 ; jj <= 2 ; jj++) {
799 a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
802 poles.SetValue(2,a_point) ;
806 index2 = myPoints->Lower() + 1 ;
807 index3 = degree + 2 ;
808 if (myTangentRequest) {
809 for (ii = myParameters->Lower() + 1 ;
810 ii < myParameters->Upper() ; ii++) {
811 parameters.SetValue(index,myParameters->Value(ii)) ;
812 poles.SetValue(index,myPoints->Value(index2)) ;
813 flatknots.SetValue(index3,myParameters->Value(ii)) ;
816 if (myTangentFlags->Value(index1)) {
818 // set the multiplicities, the order of the contact, the
821 mults.SetValue(mult_index,mults.Value(mult_index) + 1) ;
822 contact_order_array(index) = 1 ;
823 flatknots.SetValue(index3, myParameters->Value(ii)) ;
824 parameters.SetValue(index,
825 myParameters->Value(ii)) ;
826 for (jj = 1 ; jj <= 2 ; jj++) {
827 a_point.SetCoord(jj,myTangents->Value(ii).Coord(jj)) ;
829 poles.SetValue(index,a_point) ;
841 for(ii = myParameters->Lower() ; ii <= myParameters->Upper() ; ii++) {
842 parameters.SetValue(index1,
843 myParameters->Value(ii)) ;
847 for (ii = myPoints->Lower() + 1 ; ii <= myPoints->Upper() - 1 ; ii++) {
848 poles.SetValue(index,
849 myPoints->Value(ii)) ;
855 for(ii = myParameters->Lower() ; ii <= myParameters->Upper() ; ii++) {
856 flatknots.SetValue(index,
857 myParameters->Value(ii)) ;
861 for (jj = 1 ; jj <= 2 ; jj++) {
863 myTangents->Value(num_points).Coord(jj)) ;
865 poles.SetValue(num_poles-1 ,a_point) ;
867 contact_order_array.SetValue(num_poles - 1,1) ;
868 parameters.SetValue(num_poles,
869 myParameters->Value(myParameters->Upper())) ;
870 parameters.SetValue(num_poles -1,
871 myParameters->Value(myParameters->Upper())) ;
873 poles.SetValue(num_poles,
874 myPoints->Value(num_points)) ;
876 BSplCLib::Interpolate(degree,
882 if (!inversion_problem) {
884 new Geom2d_BSplineCurve(poles,
885 myParameters->Array1(),
888 myIsDone = Standard_True ;
894 //=======================================================================
895 //function : Handle(Geom2d_BSplineCurve)&
897 //=======================================================================
899 const Handle(Geom2d_BSplineCurve)& Geom2dAPI_Interpolate::Curve() const
902 throw StdFail_NotDone(" ");
908 //=======================================================================
909 //function : Geom2d_BSplineCurve
911 //=======================================================================
913 Geom2dAPI_Interpolate::operator Handle(Geom2d_BSplineCurve)() const
919 //=======================================================================
922 //=======================================================================
924 Standard_Boolean Geom2dAPI_Interpolate::IsDone() const