View Issue Details
ID  Project  Category  View Status  Date Submitted  Last Update 

0027112  Community  OCCT:Modeling Algorithms  public  20160125 20:22  20221024 10:42 
Reporter  Roman Lygin  Assigned To  aml  
Priority  normal  Severity  major  
Status  assigned  Resolution  open  
Product Version  6.9.1  
Target Version  7.8.0  
Summary  0027112: GeomAPI_Interpolate produces wrong result  
Description  Apparently GeomAPI_Interpolate tries to compute a curve using cubic spline interpolation. However it seems to incorrectly compute first and last tangents, leading to result which can be significantly far from expected. Enclosed is the screenshot of the interpolation of a curve specified via its points and approximate expected result (in blue). The source data are from a Parasolid file which contains a set of points used to define an interpolation curve and vertices an edge lying on this curve passes through. Edge vertices are shown as cubes. Open CASCADE computes first tangent calling PLib::EvalLagrange() which returns a vector lo (3.64055896 1.24782193 0.351508766). However if to compute a vector lg using formula from Golavanov's book "Geometric modeling" as: q0 = 3/2 (p1p0)/(t1t0)  1/2 q1, where t0 and t1 are parameters (0 and 29.73384) p0 and p1 are first two points q0 and q1 are derivatives at points p0 and p1 respectively, q1=p0*(t1t2)/((t0t1)(t0t2)) + p1 (2t1t0t2)/((t1t0)(t1t2))+ p2(t1t0)/((t2t0)(t2t1)). lg is (0.013074235, 0.993197391, 0.279352588) which is much closer to expected. Perhaps the bug is in broken logic in BuildTangents() in GeomAPI_Interpolate.cxx which computes tangents at first and last points. Severity is set to major as this affects data correctness of the interpolation algorithm which is a core one in OCC. A fix is left for discretion of the OCC team. Perhaps any modification may lead to significant test failures and update of the reference base.  
Steps To Reproduce  pload MODELING interpol c points.txt #compare visually and some property (e.g. length) with the expected result #test bugs moddata_3 bug27112_1 #test bugs moddata_3 bug27112_2 #test bugs moddata_3 bug27112_3 #test bugs moddata_3 bug27112_4 # No shapes needed  public available  
Tags  No tags attached.  
Test case number  

points.txt (776 bytes) 

interpolated.png (27,389 bytes) 

interpolated_expected.png (23,298 bytes) 

interpolated_tangents.png (21,963 bytes) 

interpolated_result_with_explicit_boundary_tangents.png (19,211 bytes) 

If to precompute boundary tangents according to formulas satisfying f"=0, i.e. q0 = 3/2 (p1p0)/(t1t0)  1/2 q1 qn = 3/2 (pn  pn1)/(tn  tn1)  1/2 qn1 then the curve is very much close to expected  see screenshot interpolated_result_with_explicit_boundary_tangents.png This leads to assumption that GeomAPI_Interpolate uses another boundary condition, i.e. not f"=0 (so called natural b.c.) but f'(start)=c1, f'(end)=c2 (so called 'clamped' b.c.). And for c1 and c2 it calculates tangents of a Lagrange polynom. Lagrange interpolation may produce high oscillations and respectively boundary tangents can go significantly from expected result, what is observed in this case. Of course, the user has a workaround to use Load() API to provide explicit tangents but perhaps the default behavior should be changed to better accommodate a general case. 

Branch CR27112 has been created by aml. SHA1: 3abb1ba4040703ec44aec7b31dffe0c382ad19a8 Detailed log of new commits: Author: aml Date: Sun Feb 7 09:31:58 2016 +0300 0027112: GeomAPI_Interpolate produces wrong result Code refactoring. Raw implementation of the second derivative initial condition added into tangent computation. First point tangent computed only with Golovanov formulas Last point tangent computed with help of current occt eval_lagrange. 

Branch CR27112 has been updated by aml. SHA1: c76b4885e4a20386557c2d3489111a18a649cb1a Detailed log of new commits: Author: aml Date: Tue Feb 9 14:37:26 2016 +0300 Geom2dInterploate conversion to the f'' = 0 condition approximation. 

Branch CR27112 has been updated forcibly by aml. SHA1: 6f137cbb2c6aa8462ff0437ce7905b8258f69dd6 

Branch CR27112 has been updated forcibly by aml. SHA1: 17e9de5edf0127ccb66df33243d12e1045aa010b 

Dear Roman, We've checked natural conditions and decided to make them available via flag. We cannot accept your idea as default behavior due to significant number of regressions. Can you provide more 3D test cases with improvement like in bug description? 

Dear AML, Another case recently discovered  http://www.opencascade.com/comment/19122. On our side, we now workaround the issue by explicit specifying tangents in all points and thus create an Hermite spline of degree 3 (which is C1continuous) instead of a cubic spline (which is C2continuous). There are various strategies for computing tangents. 

Dear Roman, I've implemented solution in the same manner  only approximation of this natural condition as tangents on boundaries. Our vision is to add such approach as flag. Correct implementation on "equations" level seems to be too difficult. 

Dear AML, I have looked at CR27112 branch. A few comments: 1. AFAIU, you computed q1 using Golovanov's 2011 formula and qn1 using Eval::Lagrange(). Is there a reason for this inconsistency ? 2. Could you avoid code duplication between computing first and last point and between 3D and 2D cases ? All the four versions could be merged into a single templatebased implementation. Below is an excerpt from our code to give you an idea how this might look like. 3. In general (as a separate effort though), GeomAPI_Interpolate and Geom2dAPI_Interpolate could be reimplemented using single templatebased implementation. This would help avoid inconsistencies which exist now  e.g. support of Scaling option in GeomAPI and missed in Geom2dAPI. 4. As you say, computation of the derivatives can be user parameterbased to preserve current behavior. So you might want to introduce some way (e.g. enumerationbased) to precompute q1 and qn1 and boundary conditions. As I said, we had to resort to complete computing tangents on our own and just calling GeomAPI_Interpolate to generate a C1continous Hermite spline. For various strategies/computational schemes you might want refer to comments in the below code. /*! Computes a tangent for a medium point \a pi, i.e. not on a boundary using a finite differences scheme. There are multiple alternative ways to compute tangents for an Hermite spline  see for instance here: https://en.wikipedia.org/wiki/Cubic_Hermite_spline. */ template<int D> static typename ModelData_gpTraits<D>::VecType InterimTangentByFiniteDiff ( const typename ModelData_gpTraits<D>::PointType& pi_1, const typename ModelData_gpTraits<D>::PointType& pi, const typename ModelData_gpTraits<D>::PointType& pi_p1, double ti_1, double ti, double ti_p1) { auto qi = 0.5 * ((pi_p1.Coord()  pi .Coord()) / (ti_p1  ti ) + (pi .Coord()  pi_1.Coord()) / (ti  ti_1)); return qi; } /*! Computes a boundary tangent to satisfy a requirement f" = 0. See Golovanov's book, the formula can be proved by computing second derivative for the Hermite spline r(w) = a0(w)pi + a1(w)pi+1 + b0(w)(ti+1  ti)qi + b1(w)(ti+1  ti)qi+1, w = (tti)/(ti+1  ti). There can be other approaches to compute a boundary tangent, for instance: \li using finite difference: q0 = (p1  p0) / (t1  t0). \li ensuring third derivative to be null: q0 = 2 (p1  p0) / (t1  t0)  see Golovanov's book (edition 2004, page 84). The derivative \a q1 at the neighbor point \a p1 must already be precomputed using any of the computational scheme (e.g. finite differences). */ template<int D> static typename ModelData_gpTraits<D>::VecType BoundaryTangentForNullSecondDerivative ( const typename ModelData_gpTraits<D>::PointType& p0, const typename ModelData_gpTraits<D>::PointType& p1, double t0, double t1, const typename ModelData_gpTraits<D>::VecType& q1) { auto q0 = 1.5 * typename ModelData_gpTraits<D>::VecType (p1.Coord()  p0.Coord()) / (t1  t0)  0.5 * q1; return q0; } template<int D> static typename ModelData_gpTraits<D>::VecType FirstTangentForNullSecondDerivative ( const typename ModelData_gpTraits<D>::PointType& p0, const typename ModelData_gpTraits<D>::PointType& p1, double t0, double t1, const typename ModelData_gpTraits<D>::VecType& q1) { return BoundaryTangentForNullSecondDerivative<D> (p0, p1, t0, t1, q1); } /*! Computes last tangent. Refer to FirstTangentForNullSecondDerivative() for alternatives. */ template<int D> static typename ModelData_gpTraits<D>::VecType LastTangentForNullSecondDerivative ( const typename ModelData_gpTraits<D>::PointType& pn_1, const typename ModelData_gpTraits<D>::PointType& pn, double tn_1, double tn, const typename ModelData_gpTraits<D>::VecType& qn_1) { //use reversed order return BoundaryTangentForNullSecondDerivative<D> (pn, pn_1, tn, tn_1, qn_1); } where : template<int D> struct ModelData_gpTraits; template<> struct ModelData_gpTraits<2> { typedef gp_XY CoordType; typedef gp_Pnt2d PointType; typedef gp_Dir2d DirType; typedef gp_Vec2d VecType; typedef gp_Lin2d LineType; typedef gp_Circ2d CircleType; typedef gp_Elips2d EllipseType; typedef gp_Hypr2d HyperbolaType; typedef gp_Parab2d ParabolaType; typedef TColgp_Array1OfPnt2d PointArrayType; typedef TColgp_HArray1OfPnt2d PointHArrayType; typedef Handle(TColgp_HArray1OfPnt2d) HPointArrayType; }; template<> struct ModelData_gpTraits<3> { typedef gp_XYZ CoordType; typedef gp_Pnt PointType; typedef gp_Dir DirType; typedef gp_Vec VecType; typedef gp_Lin LineType; typedef gp_Circ CircleType; typedef gp_Elips EllipseType; typedef gp_Hypr HyperbolaType; typedef gp_Parab ParabolaType; typedef TColgp_Array1OfPnt PointArrayType; typedef TColgp_HArray1OfPnt PointHArrayType; typedef Handle(TColgp_HArray1OfPnt) HPointArrayType; }; 

Dear Roman, 1. AFAIU, you computed q1 using Golovanov's 2011 formula and qn1 using Eval::Lagrange(). Is there a reason for this inconsistency ? Yes, but it used for testing. In result version there will be only Golovanov's approach for f''= 0 condition. The main reason for it  minimal number of points. For Eval::Lagrange() will require at least 4 points and Golovanov formulas will require only 3 points. 2. Could you avoid code duplication between computing first and last point and between 3D and 2D cases ? All the four versions could be merged into a single templatebased implementation. Below is an excerpt from our code to give you an idea how this might look like. Eval::Lagrange approach is not symmetric. So it is not easy to merge tangent computations. 34 I will take into account these remarks. 

Hello Roman, Just to make it clear: we have tried the proposed change on OCCT tests and observed only multiple regressions in DE. This means that current behavior looks better on our cases. The possibility to add an option to allow the user selecting the method to compute tangents will be considered but with low priority. I believe that if you use interpolator directly, e.g. for modeling purposes, you will definitely find many more points to improve, and we will be glad to see your improvements submitted to OCCT. Regarding the case discussed at http://www.opencascade.com/comment/19122 I guess that Rhino may take into account the fact that start and end points are the same, to produce symmetric curve. 

Branch CR27112 has been updated forcibly by aml. SHA1: 59499c3b3d46b35145f9b0267c5856fdcc672d0d 

Branch CR27112 has been updated forcibly by aml. SHA1: a67a5b1ab7deb33135d12b19d8befe978ac28c56 

Branch CR27112 has been updated forcibly by aml. SHA1: 94e01d856b0e47768a35c449b2501a116def3aa5 

Dear msv, Please check current state of the branch CR27112. 

Hello Andrey, Thank you for the comment. Are the regressions in tolerance and/or checkshape ? These regressions might be due to obvious difference in resulting geometry (before and after the fix). The fix does help to mitigate the situations of uneven distribution of sample points (like in the reproducer). AFAIR, in DE (ShapeConstruct_ProjectCurveOnSurface) the number of points is large enough (23 + 22*k or something like that). Therefore even for unevenly parametrized BSplines this number is sufficient to generate a more dense set of sample points, so interpolation is not vulnerable to cases like in the reproducer. I understand your hesitation to change the default behavior. This is OK with us as we found a workaround on our side to pregenerate tangents in all points. So please prioritize at your discretion. Still hope the above suggestions make sense to align codes of 3D and 2D cases, and whenever you proceed to reworking implementation they can be considered. Good luck! 

Remarks: 1) I think the methods to compute tangents are better to be moved out of Geom package. The package Geom contains only geometric structures. And auxiliary methods and calculations are placed somewhere else, e.g. GeomLib or PLib. src\Geom2dAPI\Geom2dAPI_Interpolate.hxx src\GeomAPI\GeomAPI_Interpolate.hxx 2) The method to set boundary condition type is called GetBoundaryCondition. src\Geom2dAPI\Geom2dAPI_Interpolate.cxx 3) The lines 180181 are extra. 4) In the original code of BuildPeriodicTangent(), when PointsArray.Length() == 3 then degree = 2. In the new code, degree is always 3. Is it a mistake or intention? 4) In BuildTangents(), there was check of condition PointsArray.Length() < 3, but new code is not protected against it, that is unsafe. 5) In BuildTangents(), dimension == 3 is passed in EvalLagrange, must be 2. And below iteration on ii is to be performed from 1 to 2. 6) The code in lines 263271 is better to remove, and set corresponding values directly where they are computed. src\Geom2dAPI\Geom2dAPI_Interpolate.cxx 7) In BuildPeriodicTangent, the condition PointsArray.Length() < 3 is not checked in case of GeomInterpolate_NATURAL, that is unsafe. 8) The same as 4) 9) The same as 6) src\GeometryTest\GeometryTest_ConstraintCommands.cxx 10) In help text, change interpol [3d/2d] [Natural/Clamped] [x y [z]] to interpol {3d/2d} {Natural/Clamped} {x y [z] ...} as [] usually mean optional parameter. 11) There are 3 test cases with 2D dimension, it is worth also to create a case with 2D points. 

Branch CR27112_2 has been created by aml. SHA1: 4d3688d3ce01e3a8f886bc6398f61dbd81cbd614 Detailed log of new commits: Author: aml Date: Sun Feb 7 09:31:58 2016 +0300 0027112: GeomAPI_Interpolate produces wrong result Support of the approximation of natural boundary condition (f''= 0) 

Branch CR27112_2 has been updated forcibly by aml. SHA1: 5c7f01f093929f2e087216aaa03744c067cdb486 

Branch CR27112_2 has been updated forcibly by aml. SHA1: 628580f6b1f87854dad1e8cc52b280c18956249f 

Dear msv, Please check updated branch bug27112_2. 

More remarks: 1) Please move new static methods to PLib. 2) Check 2d curve length in the 4th test case using the command "length". 

Remarks per CR27112_2 branch: 1. Code duplication is excessive. This undermines readability and creates risks of discrepancies. Please consider use of templates as suggested above and provide a single implementation. 1a. Don't make distinct names  ComputeLagrangeTangent2d() and ComputeLagrangeTangent(). Instead rely on C++ overload and name equally ComputeLagrangeTangent(). This will enable you to use templates. 1b. Provide implementation of ComputeLagrangeTangent() for 3D and 2D case using templates, avoid code duplication. Refer to provided excerpts to see how this can be easily done. 2. In ComputeLagrangeTangent() does use of direct formulas match computing by PLib::EvalLagrange()? If yes (and this must be yes!), then you should rather use PLib::EvalLagrange() because that is its purpose, again avoid code duplication  don't do manual Lagrange computation when you already have that code elsewhere. You can certainly keep formulas in the comments for easier comprehension. If using formulas provide different results from PLib::EvalLagrange() then please understand the rootcause and keep single implementation anyway. 3. Naming of enumeration. By conventions OCC enums should (or at least used to) be placed in a separate header file. Names should also contain an acronym for the type (this is an old technique to avoid name conflict, and perhaps has to be maintained until you drop preC++11 compilers). With that, it should be something like PLib_BCType.hxx enum PLib_BCType { PLib_BC_CLAMPED PLib_BC_NATURAL }; If this no longer applies, please disregard this comment. 

Branch CR27112_2 has been updated forcibly by aml. SHA1: e5976dde58caa53059d33f05f7e72e4dd9ba4908 

Dear msv, Please check updated version of the CR27112_2. I've updated test cases with "checklength" command to check length of the resulting curves. 

It is decided to take into account Roman remarks. 
Date Modified  Username  Field  Change 

20160125 20:22  Roman Lygin  New Issue  
20160125 20:22  Roman Lygin  Assigned To  => msv 
20160125 20:22  Roman Lygin  File Added: points.txt  
20160125 20:22  Roman Lygin  File Added: interpolated.png  
20160125 20:23  Roman Lygin  File Added: interpolated_expected.png  
20160125 20:25  Roman Lygin  File Added: interpolated_tangents.png  
20160127 10:32  Roman Lygin  File Added: interpolated_result_with_explicit_boundary_tangents.png  
20160127 10:47  Roman Lygin  Note Added: 0050148  
20160128 14:01  msv  Assigned To  msv => aml 
20160128 14:01  msv  Status  new => assigned 
20160207 10:05  git  Note Added: 0050455  
20160209 14:38  git  Note Added: 0050525  
20160209 14:43  git  Note Added: 0050527  
20160209 14:52  git  Note Added: 0050531  
20160210 16:02  aml  Note Added: 0050582  
20160210 16:02  aml  Assigned To  aml => Roman Lygin 
20160210 16:02  aml  Status  assigned => feedback 
20160210 16:03  aml  Note Edited: 0050582  
20160210 17:36  Roman Lygin  Note Added: 0050593  
20160210 17:38  Roman Lygin  Note Edited: 0050593  
20160210 17:44  aml  Note Added: 0050594  
20160213 09:02  Roman Lygin  Note Added: 0050716  
20160213 09:02  Roman Lygin  Assigned To  Roman Lygin => aml 
20160213 09:02  Roman Lygin  Status  feedback => assigned 
20160215 17:50  aml  Note Added: 0050770  
20160216 10:11  abv  Note Added: 0050794  
20160219 13:25  git  Note Added: 0051012  
20160220 07:37  git  Note Added: 0051033  
20160220 10:02  git  Note Added: 0051034  
20160220 10:04  aml  Note Added: 0051035  
20160220 10:04  aml  Assigned To  aml => msv 
20160220 10:04  aml  Status  assigned => resolved 
20160220 10:04  aml  Steps to Reproduce Updated  
20160220 17:08  Roman Lygin  Note Added: 0051045  
20160220 18:01  msv  Note Added: 0051046  
20160220 18:01  msv  Assigned To  msv => aml 
20160220 18:01  msv  Status  resolved => assigned 
20160224 08:09  git  Note Added: 0051068  
20160224 08:39  git  Note Added: 0051069  
20160224 11:16  git  Note Added: 0051072  
20160224 11:21  aml  Note Added: 0051074  
20160224 11:21  aml  Assigned To  aml => msv 
20160224 11:21  aml  Status  assigned => resolved 
20160224 11:21  aml  Steps to Reproduce Updated  
20160224 14:26  aml  Note Edited: 0051074  
20160224 15:32  msv  Note Added: 0051084  
20160224 15:32  msv  Assigned To  msv => aml 
20160224 15:32  msv  Status  resolved => assigned 
20160224 20:50  Roman Lygin  Note Added: 0051092  
20160225 08:46  git  Note Added: 0051094  
20160225 08:49  aml  Note Added: 0051095  
20160225 08:49  aml  Assigned To  aml => msv 
20160225 08:49  aml  Status  assigned => resolved 
20160323 10:26  aml  Note Added: 0051866  
20160323 10:26  aml  Assigned To  msv => aml 
20160323 10:26  aml  Status  resolved => assigned 
20161028 16:32  msv  Target Version  7.1.0 => 7.2.0 
20170721 11:34  msv  Target Version  7.2.0 => 7.3.0 
20171205 17:09  msv  Target Version  7.3.0 => 7.4.0 
20190812 16:44  msv  Target Version  7.4.0 => 7.5.0 
20200914 22:55  msv  Target Version  7.5.0 => 7.6.0 
20210829 18:52  msv  Target Version  7.6.0 => 7.7.0 
20221024 10:42  szy  Target Version  7.7.0 => 7.8.0 