View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0027112CommunityOCCT:Modeling Algorithmspublic2016-01-25 20:222023-08-14 17:26
ReporterRoman Lygin Assigned Toakaftasev  
PrioritynormalSeveritymajor 
Status assignedResolutionopen 
Product Version6.9.1 
Target VersionUnscheduled 
Summary0027112: GeomAPI_Interpolate produces wrong result
DescriptionApparently 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.
Steps To Reproducepload 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
TagsNo tags attached.
Test case number

Attached Files

  • 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)

Activities

Roman Lygin

2016-01-25 20:22

developer  

points.txt (776 bytes)

Roman Lygin

2016-01-25 20:22

developer  

interpolated.png (27,389 bytes)

Roman Lygin

2016-01-25 20:23

developer  

interpolated_expected.png (23,298 bytes)

Roman Lygin

2016-01-25 20:25

developer  

interpolated_tangents.png (21,963 bytes)

Roman Lygin

2016-01-27 10:32

developer  

interpolated_result_with_explicit_boundary_tangents.png (19,211 bytes)

Roman Lygin

2016-01-27 10:47

developer   ~0050148

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.

git

2016-02-07 10:05

administrator   ~0050455

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.

git

2016-02-09 14:38

administrator   ~0050525

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.

git

2016-02-09 14:43

administrator   ~0050527

Branch CR27112 has been updated forcibly by aml.

SHA-1: 6f137cbb2c6aa8462ff0437ce7905b8258f69dd6

git

2016-02-09 14:52

administrator   ~0050531

Branch CR27112 has been updated forcibly by aml.

SHA-1: 17e9de5edf0127ccb66df33243d12e1045aa010b

aml

2016-02-10 16:02

developer   ~0050582

Last edited: 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

developer   ~0050593

Last edited: 2016-02-10 17:38

 Dear AML,

Another case recently discovered - http://www.opencascade.com/comment/19122.
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.

aml

2016-02-10 17:44

developer   ~0050594

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

developer   ~0050716

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:
    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 = (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;

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;
};

aml

2016-02-15 17:50

developer   ~0050770

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.

3-4
I will take into account these remarks.

abv

2016-02-16 10:11

manager   ~0050794

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.

git

2016-02-19 13:25

administrator   ~0051012

Branch CR27112 has been updated forcibly by aml.

SHA-1: 59499c3b3d46b35145f9b0267c5856fdcc672d0d

git

2016-02-20 07:37

administrator   ~0051033

Branch CR27112 has been updated forcibly by aml.

SHA-1: a67a5b1ab7deb33135d12b19d8befe978ac28c56

git

2016-02-20 10:02

administrator   ~0051034

Branch CR27112 has been updated forcibly by aml.

SHA-1: 94e01d856b0e47768a35c449b2501a116def3aa5

aml

2016-02-20 10:04

developer   ~0051035

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

Roman Lygin

2016-02-20 17:08

developer   ~0051045

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!

msv

2016-02-20 18:01

developer   ~0051046

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 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.

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.

git

2016-02-24 08:09

administrator   ~0051068

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)

git

2016-02-24 08:39

administrator   ~0051069

Branch CR27112_2 has been updated forcibly by aml.

SHA-1: 5c7f01f093929f2e087216aaa03744c067cdb486

git

2016-02-24 11:16

administrator   ~0051072

Branch CR27112_2 has been updated forcibly by aml.

SHA-1: 628580f6b1f87854dad1e8cc52b280c18956249f

aml

2016-02-24 11:21

developer   ~0051074

Last edited: 2016-02-24 14:26

 Dear msv,
Please check updated branch bug27112_2.

msv

2016-02-24 15:32

developer   ~0051084

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

developer   ~0051092

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
PLib_BCType.hxx
enum PLib_BCType {
PLib_BC_CLAMPED
PLib_BC_NATURAL
};
If this no longer applies, please disregard this comment.

git

2016-02-25 08:46

administrator   ~0051094

Branch CR27112_2 has been updated forcibly by aml.

SHA-1: e5976dde58caa53059d33f05f7e72e4dd9ba4908

aml

2016-02-25 08:49

developer   ~0051095

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.

aml

2016-03-23 10:26

developer   ~0051866

It is decided to take into account Roman remarks.

dpasukhi

2023-08-14 17:26

administrator   ~0113976

There were only code remarks and now it can be restored again with adapting to new version of OCCT.

Issue History

Date Modified Username Field Change
2016-01-25 20:22 Roman Lygin New Issue
2016-01-25 20:22 Roman Lygin Assigned To => msv
2016-01-25 20:22 Roman Lygin File Added: points.txt
2016-01-25 20:22 Roman Lygin File Added: interpolated.png
2016-01-25 20:23 Roman Lygin File Added: interpolated_expected.png
2016-01-25 20:25 Roman Lygin File Added: interpolated_tangents.png
2016-01-27 10:32 Roman Lygin File Added: interpolated_result_with_explicit_boundary_tangents.png
2016-01-27 10:47 Roman Lygin Note Added: 0050148
2016-01-28 14:01 msv Assigned To msv => aml
2016-01-28 14:01 msv Status new => assigned
2016-02-07 10:05 git Note Added: 0050455
2016-02-09 14:38 git Note Added: 0050525
2016-02-09 14:43 git Note Added: 0050527
2016-02-09 14:52 git Note Added: 0050531
2016-02-10 16:02 aml Note Added: 0050582
2016-02-10 16:02 aml Assigned To aml => Roman Lygin
2016-02-10 16:02 aml Status assigned => feedback
2016-02-10 16:03 aml Note Edited: 0050582
2016-02-10 17:36 Roman Lygin Note Added: 0050593
2016-02-10 17:38 Roman Lygin Note Edited: 0050593
2016-02-10 17:44 aml Note Added: 0050594
2016-02-13 09:02 Roman Lygin Note Added: 0050716
2016-02-13 09:02 Roman Lygin Assigned To Roman Lygin => aml
2016-02-13 09:02 Roman Lygin Status feedback => assigned
2016-02-15 17:50 aml Note Added: 0050770
2016-02-16 10:11 abv Note Added: 0050794
2016-02-19 13:25 git Note Added: 0051012
2016-02-20 07:37 git Note Added: 0051033
2016-02-20 10:02 git Note Added: 0051034
2016-02-20 10:04 aml Note Added: 0051035
2016-02-20 10:04 aml Assigned To aml => msv
2016-02-20 10:04 aml Status assigned => resolved
2016-02-20 10:04 aml Steps to Reproduce Updated
2016-02-20 17:08 Roman Lygin Note Added: 0051045
2016-02-20 18:01 msv Note Added: 0051046
2016-02-20 18:01 msv Assigned To msv => aml
2016-02-20 18:01 msv Status resolved => assigned
2016-02-24 08:09 git Note Added: 0051068
2016-02-24 08:39 git Note Added: 0051069
2016-02-24 11:16 git Note Added: 0051072
2016-02-24 11:21 aml Note Added: 0051074
2016-02-24 11:21 aml Assigned To aml => msv
2016-02-24 11:21 aml Status assigned => resolved
2016-02-24 11:21 aml Steps to Reproduce Updated
2016-02-24 14:26 aml Note Edited: 0051074
2016-02-24 15:32 msv Note Added: 0051084
2016-02-24 15:32 msv Assigned To msv => aml
2016-02-24 15:32 msv Status resolved => assigned
2016-02-24 20:50 Roman Lygin Note Added: 0051092
2016-02-25 08:46 git Note Added: 0051094
2016-02-25 08:49 aml Note Added: 0051095
2016-02-25 08:49 aml Assigned To aml => msv
2016-02-25 08:49 aml Status assigned => resolved
2016-03-23 10:26 aml Note Added: 0051866
2016-03-23 10:26 aml Assigned To msv => aml
2016-03-23 10:26 aml Status resolved => assigned
2016-10-28 16:32 msv Target Version 7.1.0 => 7.2.0
2017-07-21 11:34 msv Target Version 7.2.0 => 7.3.0
2017-12-05 17:09 msv Target Version 7.3.0 => 7.4.0
2019-08-12 16:44 msv Target Version 7.4.0 => 7.5.0
2020-09-14 22:55 msv Target Version 7.5.0 => 7.6.0
2021-08-29 18:52 msv Target Version 7.6.0 => 7.7.0
2022-10-24 10:42 szy Target Version 7.7.0 => 7.8.0
2023-08-01 15:08 dpasukhi Target Version 7.8.0 => Unscheduled
2023-08-14 17:25 dpasukhi Assigned To aml => akaftasev
2023-08-14 17:26 dpasukhi Note Added: 0113976