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 <Geom2dAPI_Interpolate.ixx>
20 #include <Standard_ConstructionError.hxx>
22 #include <BSplCLib.hxx>
23 #include <gp_Vec2d.hxx>
24 #include <TColgp_Array1OfPnt2d.hxx>
25 #include <TColStd_HArray1OfReal.hxx>
26 #include <TColStd_Array1OfBoolean.hxx>
27 #include <TColStd_Array1OfInteger.hxx>
28 #include <TColStd_HArray1OfBoolean.hxx>
30 //=======================================================================
31 //function : CheckPoints
33 //=======================================================================
35 static Standard_Boolean CheckPoints(const TColgp_Array1OfPnt2d& PointArray,
36 const Standard_Real Tolerance)
39 Standard_Real tolerance_squared = Tolerance * Tolerance,
41 Standard_Boolean result = Standard_True ;
42 for (ii = PointArray.Lower() ; result && ii < PointArray.Upper() ; ii++) {
44 PointArray.Value(ii).SquareDistance(PointArray.Value(ii+1)) ;
45 result = (distance_squared >= tolerance_squared) ;
50 //=======================================================================
51 //function : CheckTangents
53 //=======================================================================
54 static Standard_Boolean CheckTangents(
55 const TColgp_Array1OfVec2d& Tangents,
56 const TColStd_Array1OfBoolean& TangentFlags,
57 const Standard_Real Tolerance)
61 Standard_Real tolerance_squared = Tolerance * Tolerance,
63 Standard_Boolean result = Standard_True ;
64 index = TangentFlags.Lower() ;
65 for (ii = Tangents.Lower(); result && ii <= Tangents.Upper() ; ii++) {
66 if(TangentFlags.Value(index)) {
68 Tangents.Value(ii).SquareMagnitude() ;
69 result = (distance_squared >= tolerance_squared) ;
76 //=======================================================================
77 //function : CheckParameters
79 //=======================================================================
80 static Standard_Boolean CheckParameters(const
81 TColStd_Array1OfReal& Parameters)
84 Standard_Real distance ;
85 Standard_Boolean result = Standard_True ;
86 for (ii = Parameters.Lower() ; result && ii < Parameters.Upper() ; ii++) {
88 Parameters.Value(ii+1) - Parameters.Value(ii) ;
89 result = (distance >= RealSmall()) ;
93 //=======================================================================
94 //function : BuildParameters
96 //=======================================================================
97 static void BuildParameters(const Standard_Boolean PeriodicFlag,
98 const TColgp_Array1OfPnt2d& PointsArray,
99 Handle(TColStd_HArray1OfReal)& ParametersPtr)
103 Standard_Real distance ;
105 num_parameters = PointsArray.Length() ;
107 num_parameters += 1 ;
110 new TColStd_HArray1OfReal(1,
112 ParametersPtr->SetValue(1,0.0e0) ;
114 for (ii = PointsArray.Lower() ; ii < PointsArray.Upper() ; ii++) {
116 PointsArray.Value(ii).Distance(PointsArray.Value(ii+1)) ;
117 ParametersPtr->SetValue(index,
118 ParametersPtr->Value(ii) + distance) ;
123 PointsArray.Value(PointsArray.Upper()).
124 Distance(PointsArray.Value(PointsArray.Lower())) ;
125 ParametersPtr->SetValue(index,
126 ParametersPtr->Value(ii) + distance) ;
129 //=======================================================================
130 //function : BuildPeriodicTangents
132 //=======================================================================
134 static void BuildPeriodicTangent(
135 const TColgp_Array1OfPnt2d& PointsArray,
136 TColgp_Array1OfVec2d& TangentsArray,
137 TColStd_Array1OfBoolean& TangentFlags,
138 const TColStd_Array1OfReal& ParametersArray)
143 Standard_Real *point_array,
149 if (PointsArray.Length() < 3) {
150 Standard_ConstructionError::Raise();
153 if (!TangentFlags.Value(1)) {
155 if (PointsArray.Length() == 3) {
158 point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ;
160 (Standard_Real *) &ParametersArray.Value(1) ;
161 TangentFlags.SetValue(1,Standard_True) ;
162 PLib::EvalLagrange(ParametersArray.Value(1),
169 for (ii = 1 ; ii <= 2 ; ii++) {
170 a_vector.SetCoord(ii,eval_result[1][ii-1]) ;
172 TangentsArray.SetValue(1,a_vector) ;
175 //=======================================================================
176 //function : BuildTangents
178 //=======================================================================
180 static void BuildTangents(const TColgp_Array1OfPnt2d& PointsArray,
181 TColgp_Array1OfVec2d& TangentsArray,
182 TColStd_Array1OfBoolean& TangentFlags,
183 const TColStd_Array1OfReal& ParametersArray)
187 Standard_Real *point_array,
195 if ( PointsArray.Length() < 3) {
196 Standard_ConstructionError::Raise();
198 if (PointsArray.Length() == 3) {
201 if (!TangentFlags.Value(1)) {
202 point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ;
204 (Standard_Real *) &ParametersArray.Value(1) ;
205 TangentFlags.SetValue(1,Standard_True) ;
206 PLib::EvalLagrange(ParametersArray.Value(1),
213 for (ii = 1 ; ii <= 2 ; ii++) {
214 a_vector.SetCoord(ii,eval_result[1][ii-1]) ;
216 TangentsArray.SetValue(1,a_vector) ;
218 if (! TangentFlags.Value(TangentFlags.Upper())) {
220 (Standard_Real *) &PointsArray.Value(PointsArray.Upper() - degree) ;
221 TangentFlags.SetValue(TangentFlags.Upper(),Standard_True) ;
223 (Standard_Real *)&ParametersArray.Value(ParametersArray.Upper() - degree) ;
224 PLib::EvalLagrange(ParametersArray.Value(ParametersArray.Upper()),
231 for (ii = 1 ; ii <= 2 ; ii++) {
232 a_vector.SetCoord(ii,eval_result[1][ii-1]) ;
234 TangentsArray.SetValue(TangentsArray.Upper(),a_vector) ;
237 //=======================================================================
238 //function : BuildTangents
239 //purpose : scale the given tangent so that they have the length of
240 // the size of the derivative of the lagrange interpolation
242 //=======================================================================
243 static void ScaleTangents(const TColgp_Array1OfPnt2d& PointsArray,
244 TColgp_Array1OfVec2d& TangentsArray,
245 const TColStd_Array1OfBoolean& TangentFlags,
246 const TColStd_Array1OfReal& ParametersArray)
254 Standard_Real *point_array,
262 num_points = PointsArray.Length() ;
263 if (num_points == 2) {
266 else if (num_points >= 3) {
270 index = PointsArray.Lower() ;
271 for (ii = TangentFlags.Lower() ; ii <= TangentFlags.Upper() ; ii++) {
272 if (TangentFlags.Value(ii)) {
274 (Standard_Real *) &PointsArray.Value(index) ;
276 (Standard_Real *) &ParametersArray.Value(index) ;
277 PLib::EvalLagrange(ParametersArray.Value(ii),
286 for (jj = 1 ; jj <= 2 ; jj++) {
287 value[0] += Abs(TangentsArray.Value(ii).Coord(jj)) ;
288 value[1] += Abs(eval_result[1][jj-1]) ;
290 ratio = value[1] / value[0] ;
291 for (jj = 1 ; jj <= 2 ; jj++) {
292 a_vector.SetCoord(jj, ratio *
293 TangentsArray.Value(ii).Coord(jj)) ;
295 TangentsArray.SetValue(ii, a_vector) ;
296 if (ii != TangentFlags.Lower()) {
299 if (index > PointsArray.Upper() - degree) {
300 index = PointsArray.Upper() - degree ;
307 //=======================================================================
308 //function : Geom2dAPI_Interpolate
310 //=======================================================================
312 Geom2dAPI_Interpolate::Geom2dAPI_Interpolate
313 (const Handle(TColgp_HArray1OfPnt2d)& PointsPtr,
314 const Standard_Boolean PeriodicFlag,
315 const Standard_Real Tolerance) :
316 myTolerance(Tolerance),
318 myIsDone(Standard_False),
319 myPeriodic(PeriodicFlag),
320 myTangentRequest(Standard_False)
323 Standard_Integer ii ;
324 Standard_Boolean result =
325 CheckPoints(PointsPtr->Array1(),
328 new TColgp_HArray1OfVec2d(myPoints->Lower(),
331 new TColStd_HArray1OfBoolean(myPoints->Lower(),
335 Standard_ConstructionError::Raise();
337 BuildParameters(PeriodicFlag,
341 for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) {
342 myTangentFlags->SetValue(ii,Standard_False) ;
349 //=======================================================================
350 //function : Geom2dAPI_Interpolate
352 //=======================================================================
354 Geom2dAPI_Interpolate::Geom2dAPI_Interpolate
355 (const Handle(TColgp_HArray1OfPnt2d)& PointsPtr,
356 const Handle(TColStd_HArray1OfReal)& ParametersPtr,
357 const Standard_Boolean PeriodicFlag,
358 const Standard_Real Tolerance) :
359 myTolerance(Tolerance),
361 myIsDone(Standard_False),
362 myParameters(ParametersPtr),
363 myPeriodic(PeriodicFlag),
364 myTangentRequest(Standard_False)
366 Standard_Integer ii ;
369 Standard_Boolean result =
370 CheckPoints(PointsPtr->Array1(),
374 if ((PointsPtr->Length()) + 1 != ParametersPtr->Length()) {
375 Standard_ConstructionError::Raise();
379 new TColgp_HArray1OfVec2d(myPoints->Lower(),
382 new TColStd_HArray1OfBoolean(myPoints->Lower(),
386 Standard_ConstructionError::Raise();
390 CheckParameters(ParametersPtr->Array1()) ;
392 Standard_ConstructionError::Raise();
395 for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) {
396 myTangentFlags->SetValue(ii,Standard_False) ;
400 //=======================================================================
403 //=======================================================================
405 void Geom2dAPI_Interpolate::Load(
406 const TColgp_Array1OfVec2d& Tangents,
407 const Handle(TColStd_HArray1OfBoolean)& TangentFlagsPtr)
410 Standard_Boolean result ;
411 Standard_Integer ii ;
412 myTangentRequest = Standard_True ;
413 myTangentFlags = TangentFlagsPtr ;
414 if (Tangents.Length() != myPoints->Length() ||
415 TangentFlagsPtr->Length() != myPoints->Length()) {
416 Standard_ConstructionError::Raise();
419 CheckTangents(Tangents,
420 TangentFlagsPtr->Array1(),
424 new TColgp_HArray1OfVec2d(Tangents.Lower(),Tangents.Upper()) ;
425 for (ii = Tangents.Lower() ; ii <= Tangents.Upper() ; ii++ ) {
426 myTangents->SetValue(ii,Tangents.Value(ii)) ;
428 ScaleTangents(myPoints->Array1(),
429 myTangents->ChangeArray1(),
430 TangentFlagsPtr->Array1(),
431 myParameters->Array1()) ;
434 Standard_ConstructionError::Raise();
440 //=======================================================================
443 //=======================================================================
445 void Geom2dAPI_Interpolate::Load(const gp_Vec2d& InitialTangent,
446 const gp_Vec2d& FinalTangent)
448 Standard_Boolean result ;
449 myTangentRequest = Standard_True ;
450 myTangentFlags->SetValue(1,Standard_True) ;
451 myTangentFlags->SetValue(myPoints->Length(),Standard_True) ;
452 myTangents->SetValue(1,InitialTangent) ;
453 myTangents->SetValue(myPoints->Length(),FinalTangent);
455 CheckTangents(myTangents->Array1(),
456 myTangentFlags->Array1(),
459 Standard_ConstructionError::Raise();
461 ScaleTangents(myPoints->Array1(),
462 myTangents->ChangeArray1(),
463 myTangentFlags->Array1(),
464 myParameters->Array1()) ;
467 //=======================================================================
470 //=======================================================================
472 void Geom2dAPI_Interpolate::Perform()
478 PerformNonPeriodic() ;
481 //=======================================================================
482 //function : PerformPeriodic
484 //=======================================================================
486 void Geom2dAPI_Interpolate::PerformPeriodic()
488 Standard_Integer degree,
501 Standard_Real period ;
505 num_points = myPoints->Length() ;
506 period = myParameters->Value(myParameters->Upper()) -
507 myParameters->Value(myParameters->Lower()) ;
508 num_poles = num_points + 1 ;
509 if (num_points == 2 && !myTangentRequest) {
511 // build a periodic curve of degree 1
515 TColStd_Array1OfInteger deg1_mults(1,num_poles) ;
516 for (ii = 1 ; ii <= num_poles ; ii++) {
517 deg1_mults.SetValue(ii,1) ;
521 new Geom2d_BSplineCurve(myPoints->Array1(),
522 myParameters->Array1(),
526 myIsDone = Standard_True ;
530 num_distinct_knots = num_points + 1 ;
534 if (myTangentRequest)
535 for (ii = myTangentFlags->Lower() + 1 ;
536 ii <= myTangentFlags->Upper() ; ii++) {
537 if (myTangentFlags->Value(ii)) {
542 TColStd_Array1OfReal parameters(1,num_poles) ;
543 TColStd_Array1OfReal flatknots(1,num_poles + degree + 1) ;
544 TColStd_Array1OfInteger mults(1,num_distinct_knots) ;
545 TColStd_Array1OfInteger contact_order_array(1, num_poles) ;
546 TColgp_Array1OfPnt2d poles(1,num_poles) ;
548 for (ii = 1 ; ii <= half_order ; ii++) {
549 flatknots.SetValue(ii,myParameters->Value(myParameters->Upper() -1) -
551 flatknots.SetValue(ii + half_order,myParameters->
552 Value(myParameters->Lower())) ;
553 flatknots.SetValue(num_poles + ii,
554 myParameters->Value(myParameters->Upper())) ;
555 flatknots.SetValue(num_poles + half_order + ii,
556 myParameters->Value(half_order) + period) ;
558 for (ii = 1 ; ii <= num_poles ; ii++) {
559 contact_order_array.SetValue(ii,0) ;
561 for (ii = 2 ; ii < num_distinct_knots ; ii++) {
562 mults.SetValue(ii,1) ;
564 mults.SetValue(1,half_order) ;
565 mults.SetValue(num_distinct_knots ,half_order) ;
566 if (num_points >= 3) {
569 // only enter here if there are more than 3 points otherwise
570 // it means we have already the tangent
572 BuildPeriodicTangent(myPoints->Array1(),
573 myTangents->ChangeArray1(),
574 myTangentFlags->ChangeArray1(),
575 myParameters->Array1()) ;
577 contact_order_array.SetValue(2,1) ;
578 parameters.SetValue(1,myParameters->Value(1)) ;
579 parameters.SetValue(2,myParameters->Value(1)) ;
580 poles.SetValue(1,myPoints->Value(1)) ;
581 for (jj = 1 ; jj <= 2 ; jj++) {
582 a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
584 poles.SetValue(2,a_point) ;
588 index1 = degree + 2 ;
589 if (myTangentRequest) {
590 for (ii = myTangentFlags->Lower() + 1 ;
591 ii <= myTangentFlags->Upper() ; ii++) {
592 parameters.SetValue(index,myParameters->Value(ii)) ;
593 flatknots.SetValue(index1,myParameters->Value(ii)) ;
594 poles.SetValue(index,myPoints->Value(ii)) ;
597 if (myTangentFlags->Value(ii)) {
598 mults.SetValue(mult_index,mults.Value(mult_index) + 1) ;
599 contact_order_array(index) = 1 ;
601 parameters.SetValue(index,
602 myParameters->Value(ii)) ;
603 flatknots.SetValue(index1,myParameters->Value(ii)) ;
604 for (jj = 1 ; jj <= 2 ; jj++) {
605 a_point.SetCoord(jj,myTangents->Value(ii).Coord(jj)) ;
607 poles.SetValue(index,a_point) ;
617 for(ii = myParameters->Lower() ; ii <= myParameters->Upper() ; ii++) {
618 parameters.SetValue(index1,
619 myParameters->Value(ii)) ;
620 flatknots.SetValue(index,
621 myParameters->Value(ii)) ;
626 for (ii = myPoints->Lower() + 1 ; ii <= myPoints->Upper() ; ii++) {
628 // copy all the given points since the last one will be initialized
629 // below by the first point in the array myPoints
631 poles.SetValue(index,
632 myPoints->Value(ii)) ;
637 contact_order_array.SetValue(num_poles - 1, 1) ;
638 parameters.SetValue(num_poles-1,
639 myParameters->Value(myParameters->Upper())) ;
641 // for the periodic curve ONLY the tangent of the first point
642 // will be used since the curve should close itself at the first
643 // point See BuildPeriodicTangent
645 for (jj = 1 ; jj <= 2 ; jj++) {
646 a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
648 poles.SetValue(num_poles-1,a_point) ;
650 parameters.SetValue(num_poles,
651 myParameters->Value(myParameters->Upper())) ;
653 poles.SetValue(num_poles,
654 myPoints->Value(1)) ;
657 BSplCLib::Interpolate(degree,
663 if (!inversion_problem) {
664 TColgp_Array1OfPnt2d newpoles(poles.Value(1),
668 new Geom2d_BSplineCurve(newpoles,
669 myParameters->Array1(),
673 myIsDone = Standard_True ;
679 //=======================================================================
680 //function : PerformNonPeriodic
682 //=======================================================================
684 void Geom2dAPI_Interpolate::PerformNonPeriodic()
686 Standard_Integer degree,
703 num_poles = myPoints->Length() ;
704 if (num_poles == 2 && !myTangentRequest) {
707 else if (num_poles == 3 && !myTangentRequest) {
709 num_distinct_knots = 2 ;
714 if (myTangentRequest)
715 for (ii = myTangentFlags->Lower() + 1 ;
716 ii < myTangentFlags->Upper() ; ii++) {
717 if (myTangentFlags->Value(ii)) {
724 TColStd_Array1OfReal parameters(1,num_poles) ;
725 TColStd_Array1OfReal flatknots(1,num_poles + degree + 1) ;
726 TColStd_Array1OfInteger mults(1,num_distinct_knots) ;
727 TColStd_Array1OfReal knots(1,num_distinct_knots) ;
728 TColStd_Array1OfInteger contact_order_array(1, num_poles) ;
729 TColgp_Array1OfPnt2d poles(1,num_poles) ;
731 for (ii = 1 ; ii <= degree + 1 ; ii++) {
732 flatknots.SetValue(ii,myParameters->Value(1)) ;
733 flatknots.SetValue(ii + num_poles,
734 myParameters->Value(num_points)) ;
736 for (ii = 1 ; ii <= num_poles ; ii++) {
737 contact_order_array.SetValue(ii,0) ;
739 for (ii = 2 ; ii < num_distinct_knots ; ii++) {
740 mults.SetValue(ii,1) ;
742 mults.SetValue(1,degree + 1) ;
743 mults.SetValue(num_distinct_knots ,degree + 1) ;
747 for (ii = 1 ; ii <= num_poles ; ii++) {
748 poles.SetValue(ii ,myPoints->Value(ii)) ;
751 new Geom2d_BSplineCurve(poles,
752 myParameters->Array1(),
755 myIsDone = Standard_True ;
758 knots.SetValue(1,myParameters->Value(1)) ;
759 knots.SetValue(2,myParameters->Value(3)) ;
760 for (ii = 1 ; ii <= num_poles ; ii++) {
761 poles.SetValue(ii,myPoints->Value(ii)) ;
764 BSplCLib::Interpolate(degree,
766 myParameters->Array1(),
770 if (!inversion_problem) {
772 new Geom2d_BSplineCurve(poles,
776 myIsDone = Standard_True ;
781 // check if the boundary conditions are set
783 if (num_points >= 3) {
785 // cannot build the tangents with degree 3 with only 2 points
786 // if those where not given in advance
788 BuildTangents(myPoints->Array1(),
789 myTangents->ChangeArray1(),
790 myTangentFlags->ChangeArray1(),
791 myParameters->Array1()) ;
793 contact_order_array.SetValue(2,1) ;
794 parameters.SetValue(1,myParameters->Value(1)) ;
795 parameters.SetValue(2,myParameters->Value(1)) ;
796 poles.SetValue(1,myPoints->Value(1)) ;
797 for (jj = 1 ; jj <= 2 ; jj++) {
798 a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
801 poles.SetValue(2,a_point) ;
805 index2 = myPoints->Lower() + 1 ;
806 index3 = degree + 2 ;
807 if (myTangentRequest) {
808 for (ii = myParameters->Lower() + 1 ;
809 ii < myParameters->Upper() ; ii++) {
810 parameters.SetValue(index,myParameters->Value(ii)) ;
811 poles.SetValue(index,myPoints->Value(index2)) ;
812 flatknots.SetValue(index3,myParameters->Value(ii)) ;
815 if (myTangentFlags->Value(index1)) {
817 // set the multiplicities, the order of the contact, the
820 mults.SetValue(mult_index,mults.Value(mult_index) + 1) ;
821 contact_order_array(index) = 1 ;
822 flatknots.SetValue(index3, myParameters->Value(ii)) ;
823 parameters.SetValue(index,
824 myParameters->Value(ii)) ;
825 for (jj = 1 ; jj <= 2 ; jj++) {
826 a_point.SetCoord(jj,myTangents->Value(ii).Coord(jj)) ;
828 poles.SetValue(index,a_point) ;
840 for(ii = myParameters->Lower() ; ii <= myParameters->Upper() ; ii++) {
841 parameters.SetValue(index1,
842 myParameters->Value(ii)) ;
846 for (ii = myPoints->Lower() + 1 ; ii <= myPoints->Upper() - 1 ; ii++) {
847 poles.SetValue(index,
848 myPoints->Value(ii)) ;
854 for(ii = myParameters->Lower() ; ii <= myParameters->Upper() ; ii++) {
855 flatknots.SetValue(index,
856 myParameters->Value(ii)) ;
860 for (jj = 1 ; jj <= 2 ; jj++) {
862 myTangents->Value(num_points).Coord(jj)) ;
864 poles.SetValue(num_poles-1 ,a_point) ;
866 contact_order_array.SetValue(num_poles - 1,1) ;
867 parameters.SetValue(num_poles,
868 myParameters->Value(myParameters->Upper())) ;
869 parameters.SetValue(num_poles -1,
870 myParameters->Value(myParameters->Upper())) ;
872 poles.SetValue(num_poles,
873 myPoints->Value(num_points)) ;
875 BSplCLib::Interpolate(degree,
881 if (!inversion_problem) {
883 new Geom2d_BSplineCurve(poles,
884 myParameters->Array1(),
887 myIsDone = Standard_True ;
893 //=======================================================================
894 //function : Handle(Geom2d_BSplineCurve)&
896 //=======================================================================
898 const Handle(Geom2d_BSplineCurve)& Geom2dAPI_Interpolate::Curve() const
901 StdFail_NotDone::Raise(" ");
907 //=======================================================================
908 //function : Geom2d_BSplineCurve
910 //=======================================================================
912 Geom2dAPI_Interpolate::operator Handle(Geom2d_BSplineCurve)() const
918 //=======================================================================
921 //=======================================================================
923 Standard_Boolean Geom2dAPI_Interpolate::IsDone() const