|
Notes |
|
|
|
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. |
|
|
|
(0050455)
|
|
git
|
|
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.
|
|
|
|
(0050525)
|
|
git
|
|
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.
|
|
|
|
(0050527)
|
|
git
|
|
2016-02-09 14:43
|
|
Branch CR27112 has been updated forcibly by aml.
SHA-1: 6f137cbb2c6aa8462ff0437ce7905b8258f69dd6 |
|
|
|
(0050531)
|
|
git
|
|
2016-02-09 14:52
|
|
Branch CR27112 has been updated forcibly by aml.
SHA-1: 17e9de5edf0127ccb66df33243d12e1045aa010b |
|
|
|
(0050582)
|
|
aml
|
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?
|
|
|
|
(0050593)
|
|
Roman Lygin
|
2016-02-10 17:36
(edited on: 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.
|
|
|
|
(0050594)
|
|
aml
|
|
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. |
|
|
|
|
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;
};
|
|
|
|
(0050770)
|
|
aml
|
|
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.
3-4
I will take into account these remarks. |
|
|
|
(0050794)
|
|
abv
|
|
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 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. |
|
|
|
(0051012)
|
|
git
|
|
2016-02-19 13:25
|
|
Branch CR27112 has been updated forcibly by aml.
SHA-1: 59499c3b3d46b35145f9b0267c5856fdcc672d0d |
|
|
|
(0051033)
|
|
git
|
|
2016-02-20 07:37
|
|
Branch CR27112 has been updated forcibly by aml.
SHA-1: a67a5b1ab7deb33135d12b19d8befe978ac28c56 |
|
|
|
(0051034)
|
|
git
|
|
2016-02-20 10:02
|
|
Branch CR27112 has been updated forcibly by aml.
SHA-1: 94e01d856b0e47768a35c449b2501a116def3aa5 |
|
|
|
(0051035)
|
|
aml
|
|
2016-02-20 10:04
|
|
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 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! |
|
|
|
(0051046)
|
|
msv
|
|
2016-02-20 18:01
|
|
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. |
|
|
|
(0051068)
|
|
git
|
|
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)
|
|
|
|
(0051069)
|
|
git
|
|
2016-02-24 08:39
|
|
Branch CR27112_2 has been updated forcibly by aml.
SHA-1: 5c7f01f093929f2e087216aaa03744c067cdb486 |
|
|
|
(0051072)
|
|
git
|
|
2016-02-24 11:16
|
|
Branch CR27112_2 has been updated forcibly by aml.
SHA-1: 628580f6b1f87854dad1e8cc52b280c18956249f |
|
|
|
(0051074)
|
|
aml
|
2016-02-24 11:21
(edited on: 2016-02-24 14:26) |
|
Dear msv,
Please check updated branch bug27112_2.
|
|
|
|
(0051084)
|
|
msv
|
|
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". |
|
|
|
|
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. |
|
|
|
(0051094)
|
|
git
|
|
2016-02-25 08:46
|
|
Branch CR27112_2 has been updated forcibly by aml.
SHA-1: e5976dde58caa53059d33f05f7e72e4dd9ba4908 |
|
|
|
(0051095)
|
|
aml
|
|
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. |
|
|
|
(0051866)
|
|
aml
|
|
2016-03-23 10:26
|
|
|
It is decided to take into account Roman remarks. |
|
points.txt (776) 2016-01-25 20:22
interpolated.png (27,389) 2016-01-25 20:22