MantisBT - Community
View Issue Details
0027112Community[OCCT] OCCT:Modeling Algorithmspublic2016-01-25 20:222020-09-14 22:55
Roman Lygin 
[OCCT] 6.9.1 
[OCCT] 7.6.0* 
0027112: GeomAPI_Interpolate produces wrong result
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 (p1-p0)/(t1-t0) - 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*(t1-t2)/((t0-t1)(t0-t2)) + p1 (2t1-t0-t2)/((t1-t0)(t1-t2))+ p2(t1-t0)/((t2-t0)(t2-t1)).

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.
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
No tags attached.
txt points.txt (776) 2016-01-25 20:22
png interpolated.png (27,389) 2016-01-25 20:22
png interpolated_expected.png (23,298) 2016-01-25 20:23
png interpolated_tangents.png (21,963) 2016-01-25 20:25
png interpolated_result_with_explicit_boundary_tangents.png (19,211) 2016-01-27 10:32
Issue History
2016-01-25 20:22Roman LyginNew Issue
2016-01-25 20:22Roman LyginAssigned To => msv
2016-01-25 20:22Roman LyginFile Added: points.txt
2016-01-25 20:22Roman LyginFile Added: interpolated.png
2016-01-25 20:23Roman LyginFile Added: interpolated_expected.png
2016-01-25 20:25Roman LyginFile Added: interpolated_tangents.png
2016-01-27 10:32Roman LyginFile Added: interpolated_result_with_explicit_boundary_tangents.png
2016-01-27 10:47Roman LyginNote Added: 0050148
2016-01-28 14:01msvAssigned Tomsv => aml
2016-01-28 14:01msvStatusnew => assigned
2016-02-07 10:05gitNote Added: 0050455
2016-02-09 14:38gitNote Added: 0050525
2016-02-09 14:43gitNote Added: 0050527
2016-02-09 14:52gitNote Added: 0050531
2016-02-10 16:02amlNote Added: 0050582
2016-02-10 16:02amlAssigned Toaml => Roman Lygin
2016-02-10 16:02amlStatusassigned => feedback
2016-02-10 16:03amlNote Edited: 0050582bug_revision_view_page.php?bugnote_id=50582#r12902
2016-02-10 17:36Roman LyginNote Added: 0050593
2016-02-10 17:38Roman LyginNote Edited: 0050593bug_revision_view_page.php?bugnote_id=50593#r12906
2016-02-10 17:44amlNote Added: 0050594
2016-02-13 09:02Roman LyginNote Added: 0050716
2016-02-13 09:02Roman LyginAssigned ToRoman Lygin => aml
2016-02-13 09:02Roman LyginStatusfeedback => assigned
2016-02-15 17:50amlNote Added: 0050770
2016-02-16 10:11abvNote Added: 0050794
2016-02-19 13:25gitNote Added: 0051012
2016-02-20 07:37gitNote Added: 0051033
2016-02-20 10:02gitNote Added: 0051034
2016-02-20 10:04amlNote Added: 0051035
2016-02-20 10:04amlAssigned Toaml => msv
2016-02-20 10:04amlStatusassigned => resolved
2016-02-20 10:04amlSteps to Reproduce Updatedbug_revision_view_page.php?rev_id=13036#r13036
2016-02-20 17:08Roman LyginNote Added: 0051045
2016-02-20 18:01msvNote Added: 0051046
2016-02-20 18:01msvAssigned Tomsv => aml
2016-02-20 18:01msvStatusresolved => assigned
2016-02-24 08:09gitNote Added: 0051068
2016-02-24 08:39gitNote Added: 0051069
2016-02-24 11:16gitNote Added: 0051072
2016-02-24 11:21amlNote Added: 0051074
2016-02-24 11:21amlAssigned Toaml => msv
2016-02-24 11:21amlStatusassigned => resolved
2016-02-24 11:21amlSteps to Reproduce Updatedbug_revision_view_page.php?rev_id=13047#r13047
2016-02-24 14:26amlNote Edited: 0051074bug_revision_view_page.php?bugnote_id=51074#r13051
2016-02-24 15:32msvNote Added: 0051084
2016-02-24 15:32msvAssigned Tomsv => aml
2016-02-24 15:32msvStatusresolved => assigned
2016-02-24 20:50Roman LyginNote Added: 0051092
2016-02-25 08:46gitNote Added: 0051094
2016-02-25 08:49amlNote Added: 0051095
2016-02-25 08:49amlAssigned Toaml => msv
2016-02-25 08:49amlStatusassigned => resolved
2016-03-23 10:26amlNote Added: 0051866
2016-03-23 10:26amlAssigned Tomsv => aml
2016-03-23 10:26amlStatusresolved => assigned
2016-10-28 16:32msvTarget Version7.1.0 => 7.2.0
2017-07-21 11:34msvTarget Version7.2.0 => 7.3.0
2017-12-05 17:09msvTarget Version7.3.0 => 7.4.0
2019-08-12 16:44msvTarget Version7.4.0 => 7.5.0
2020-09-14 22:55msvTarget Version7.5.0 => 7.6.0*

Roman Lygin   
2016-01-27 10:47   
If to precompute boundary tangents according to formulas satisfying f"=0, i.e.

q0 = 3/2 (p1-p0)/(t1-t0) - 1/2 q1
qn = 3/2 (pn - pn-1)/(tn - tn-1) - 1/2 qn-1

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 work-around to use Load() API to provide explicit tangents but perhaps the default behavior should be changed to better accommodate a general case.
2016-02-07 10:05   
Branch CR27112 has been created by aml.

SHA-1: 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.
2016-02-09 14:38   
Branch CR27112 has been updated by aml.

SHA-1: 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.

2016-02-09 14:43   
Branch CR27112 has been updated forcibly by aml.

SHA-1: 6f137cbb2c6aa8462ff0437ce7905b8258f69dd6
2016-02-09 14:52   
Branch CR27112 has been updated forcibly by aml.

SHA-1: 17e9de5edf0127ccb66df33243d12e1045aa010b
2016-02-10 16:02   
(edited on: 2016-02-10 16:03)
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?

Roman Lygin   
2016-02-10 17:36   
(edited on: 2016-02-10 17:38)
Dear AML,

Another case recently discovered - [^]
On our side, we now work-around the issue by explicit specifying tangents in all points and thus create an Hermite spline of degree 3 (which is C1-continuous) instead of a cubic spline (which is C2-continuous). There are various strategies for computing tangents.

2016-02-10 17:44   
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.
Roman Lygin   
2016-02-13 09:02   
Dear AML,

I have looked at CR27112 branch. A few comments:

1. AFAIU, you computed q1 using Golovanov's 2011 formula and qn-1 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 template-based 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 template-based 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 parameter-based to preserve current behavior. So you might want to introduce some way (e.g. enumeration-based) to precompute q1 and qn-1 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 C1-continous 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: [^]
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 = (t-ti)/(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;

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;

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;

2016-02-15 17:50   
Dear Roman,

1. AFAIU, you computed q1 using Golovanov's 2011 formula and qn-1 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 template-based 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.

I will take into account these remarks.
2016-02-16 10:11   
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 [^] I guess that Rhino may take into account the fact that start and end points are the same, to produce symmetric curve.
2016-02-19 13:25   
Branch CR27112 has been updated forcibly by aml.

SHA-1: 59499c3b3d46b35145f9b0267c5856fdcc672d0d
2016-02-20 07:37   
Branch CR27112 has been updated forcibly by aml.

SHA-1: a67a5b1ab7deb33135d12b19d8befe978ac28c56
2016-02-20 10:02   
Branch CR27112 has been updated forcibly by aml.

SHA-1: 94e01d856b0e47768a35c449b2501a116def3aa5
2016-02-20 10:04   
Dear msv,
Please check current state of the branch CR27112.
Roman Lygin   
2016-02-20 17:08   
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 B-Splines 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 work-around on our side to pre-generate 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!
2016-02-20 18:01   

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.


2) The method to set boundary condition type is called GetBoundaryCondition.


3) The lines 180-181 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 263-271 is better to remove, and set corresponding values directly where they are computed.


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)


10) In help text, change
interpol [3d/2d] [Natural/Clamped] [x y [z]]
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.
2016-02-24 08:09   
Branch CR27112_2 has been created by aml.

SHA-1: 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)
2016-02-24 08:39   
Branch CR27112_2 has been updated forcibly by aml.

SHA-1: 5c7f01f093929f2e087216aaa03744c067cdb486
2016-02-24 11:16   
Branch CR27112_2 has been updated forcibly by aml.

SHA-1: 628580f6b1f87854dad1e8cc52b280c18956249f
2016-02-24 11:21   
(edited on: 2016-02-24 14:26)
Dear msv,
Please check updated branch bug27112_2.

2016-02-24 15:32   
More remarks:

1) Please move new static methods to PLib.
2) Check 2d curve length in the 4th test case using the command "length".
Roman Lygin   
2016-02-24 20:50   
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 root-cause 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 pre-C++11 compilers).
With that, it should be something like
enum PLib_BCType {
If this no longer applies, please disregard this comment.
2016-02-25 08:46   
Branch CR27112_2 has been updated forcibly by aml.

SHA-1: e5976dde58caa53059d33f05f7e72e4dd9ba4908
2016-02-25 08:49   
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.
2016-03-23 10:26   
It is decided to take into account Roman remarks.