Issues with cone fitting using ApprCone3

[ruellm]

I am getting a very wrong result when using GTE for fitting a cone to the points from a mesh i created I uploaded the mesh in question and I also provided a text file that contains only the points here: https://drive.google.com/file/d/1Lhbuk2AkdTGUXaOuDZ03tlF4UC_-zWlJ/view?usp=sharing

Which looks like this (a simple cone)

I used the following code adapted from the sample here https://github.com/davideberly/GeometricTools/blob/master/GTE/Samples/Mathematics/FitCone/FitConeWindow3.cpp

The input are ALL the vertices of the STL in question.

`` size_t const maxIterations = 4; double const updateLengthTolerance = (double) 1e-04; double const errorDifferenceTolerance = (double) 1e-08; bool useConeInputAsInitialGuess = false;

double const lambdaFactor = (double) 0.001;
double const lambdaAdjust = (double) 10.0;
size_t const maxAdjustments = 8;

fitter(numPoints, conePoints, maxIterations, updateLengthTolerance, 
    errorDifferenceTolerance, lambdaFactor, lambdaAdjust, maxAdjustments, useConeInputAsInitialGuess,
    coneVertex, coneAxis, coneAngle);

the thing is, the resulting coneVertex is way off or wrong I rendered a point/circle on it just to demonstrate

my other test cases are OK except for this one, why is this? and how should I handle the issue. I used the latest version from github.

[davideberly]

I was able to reproduce the problem. The Levenberg-Marquardt minimizer cannot reduce the least-squares error for the initial cone parameters that are computed internally (you have useConeInputAsInitialGuess = false). The minimizer iterations keep multiplying lambdaFactor by lambdaAdjust, but each new factor does not reduce the error.

After some investigation, I believe there is a flaw in the derivation of the initialization of the cone parameters. The initial cone axis direction is estimated by my current code to be approximately (0,0,1), but in fact for your single-sided cone it needs to be (0,0,-1). I modified the direction when debugging to be initially (0,0,-1), and the minimizer then finds a good estimate. This argues I need to modify the initialization code to produce the axis direction that points in the direction of opening of the cone. I speculate that your other cases worked correctly because the initialization code produced an appropriate axis direction (and not its negative).

I have a fix locally that appears to work. I need to integrate this into the current ApprCone3.h file and unit test.

Finally, your data set consists of the vertex and a collection of points on a circle. If you know that you are going to fit cones to triangulated data obtained from a cone (in Blender), there are simpler methods to estimate the parameters. My "fix" to the initialization is actually one of them (a small amount of code), and it does not require you to use ApprCone3.h.

[ruellm]

Thank you for looking into it David, really appreciated, I will wait once your update is avaiable. If its not much to ask, do you have an estimate when will the update be available? :)

Can you also share the simple method you just described to estimate the parameters without using ApprCone3.h? Maybe I can use that as the input for the initial guess for ApprCone3 as well. s well.

I also have another case with another cone below, I think its the same issue related to axis result. (I uploaded in the same folder as before as cone6.stl or cone6-points.txt)

when extracted, the resulting conevertex is like below when using all points.

Another behaviour is, if I just use few points from cone6, it looks like this (with cone drawing and conevertex illustration) side view

This is how I did to generate the other parameters of the cone to generate my mesh drawing

    Real hmin = std::numeric_limits<Real>::max();
Real hmax = 0.0;
for (size_t i = 0; i < numPoints; ++i)
{
    Vector3<Real> diff = conePoints[i] - coneVertex;
    Real h = Dot(coneAxis, diff);
    if (h > hmax)
    {
        hmax = h;
    }
    if (h < hmin)
    {
        hmin = h;
    }
}

// clamp and avoid negative
hmin = std::max<Real>(hmin, 0);

// compute the radius based on computed height
topRadius = hmin * tan(coneAngle);
bottomRadius = hmax * tan(coneAngle);

// top and bottom position
top = vertex + (axis * hmin);
bottom = vertex + (axis * hmax);
height = (bottom - top).Length();

the points used are uploaded in the same folder as cone-6-limited.txt Will this be related to your fix or is this a different issue?

Hope you can enlighten me what is happening here, Thank you in advance!

[davideberly]

I will try to post the modification today or tomorrow. I tried the same Google drive link for the new dataset, but all I get is the old one Cone.stl. Perhaps you can give me new links or email the files to me.

Here is the simple test code for your Cone.stl file.

    STLBinaryFile stl{};
    stl.Load("Cone.stl");

    std::vector<Vector3<double>> points(3 * stl.triangles.size());
    points.clear();
    for (auto const& tri : stl.triangles)
    {
        for (size_t i = 0; i < 3; ++i)
        {
            auto const& P = tri.vertex[i];
            points.push_back({ (double)P[0], (double)P[1], (double)P[2] });
        }
    }
    // 62 triangles, 3*62 = 186 points of which 33 are unique.
    // 1 point is the vertex, the other 32 points are on a circle.

    // The following code should work well when the points are approximately
    // uniformly distributed on the cone surface.

    // Fit the points with a Gaussian distribution. 2 of the eigenvalues are
    // nearly equal, choose the eigenvector corresponding to the other
    // eigenvalue as the cone axis direction. The box.extent stores the
    // eigenvalues and box.extent[1] and box.extent[2] are nearly equal. The
    // cone axis direction is chosen to be box.axis[0].
    ApprGaussian3<double> gfitter{};
    gfitter.Fit(points);
    auto const& box = gfitter.GetParameters();
    // box.center = {2.5155299396672451e-09, -9.8058829109195312e-09, -0.65591397849462374}
    // box.axis[0] = {0.98078546828497160, 0.19508937746845470, 1.1830806899135137e-09}
    // box.axis[1] = {0.19508937746845462, -0.98078546828497126, -8.5132267664456776e-10}
    // box.axis[2] = {5.6712562731620983e-09, -2.2447173221576151e-08, -0.99999999999999967}
    // box.extent = {0.40860214688365271, 0.41935483688848924, 0.56977685281535617}

    // To illustrate the cone axis direction issue, I am choosing the
    // initial estimate to point away from the single-sided cone.
    Vector3<double> coneAxis = -box.axis[2];
    // coneAxis = {-5.6712562731620983e-09, 2.2447173221576151e-08, 0.99999999999999967}

    // The box.center is generally not the cone vertex, but we will use the
    // center and modify it to obtain a cone vertex estimate.

    // It is possible that coneAxis is directed away from the single-sided
    // cone. The next code attempts to determine that. Compute the minimum
    // and maximum projections of the points onto the line
    // box.center + t * coneAxis.
    double hmin = std::numeric_limits<double>::max(), hmax = -hmin;
    size_t imin = 0, imax = 0;
    for (size_t i = 0; i < points.size(); ++i)
    {
        double h = Dot(coneAxis, points[i] - box.center);
        if (h < hmin)
        {
            hmin = h;
            imin = i;
        }
        if (h > hmax)
        {
            hmax = h;
            imax = i;
        }
    }

    // hmin = -0.34408604462764081, imin = 44
    // hmax = 1.6559139784946235, imax = 1

    // This is the critical step. points[imin] projects to hmin and has
    // distance rmin to the cone axis. points[imax] projects to hmax and
    // has distance rmax to the cone axis. The names rmin and rmax are
    // misleading, because it is possible that rmin > rmax. The subnames
    // 'min' and 'max' refer to the ordering of hmin < hmax.
    //
    // If coneAxis is directed towards the cone opening, then rmin < rmax.
    // But if the coneAxis is directed away from the single-sided cone, then
    // rmin > rmax. In this case we need to use -coneAxis. The measurements
    // need to be swapped: rmin with rmax, hmin with hmax, points[imin] with
    // points[imax].
    //
    // NOTE: The code below uses single points for computing rmin and rmax.
    // The actual fix will compute (h,r) values and fit them with a line.
    // The slope of the line determines whether or not coneAxis must be
    // negated.
    Vector3<double> diffMin = points[imin] - box.center;
    diffMin -= Dot(coneAxis, diffMin) * coneAxis;
    double rmin = Length(diffMin);
    // rmin = 1.0000000086458720
    Vector3<double> diffMax = points[imax] - box.center;
    diffMax -= Dot(coneAxis, diffMax) * coneAxis;
    double rmax = Length(diffMax);
    // rmax = 2.8215256795713097e-08 
    if (rmax < rmin)
    {
        coneAxis = -coneAxis;
        std::swap(hmin, hmax);
        hmin = -hmin;
        hmax = -hmax;
        std::swap(rmin, rmax);
        std::swap(diffMin, diffMax);
    }

    // Create a vector diffExtreme so that coneAxis and diffExtreme lie in
    // the same plane. This allows us to compute the cone angle using a dot
    // product.
    Vector3<double> maxExtreme = hmax * coneAxis + diffMax;
    Vector3<double> minExtreme = hmin * coneAxis + (rmin / rmax) * diffMax;
    Vector3<double> diffExtreme = maxExtreme - minExtreme;
    Normalize(diffExtreme);
    double cosConeAngle = Dot(coneAxis, diffExtreme);  // cosConeAngle = 0.89442719656871184
    double coneAngle = std::acos(cosConeAngle);  // coneAngle = 0.46364759654859966

    // Once you have the cone angle and cone axis direction, the cone
    // vertex can be computed. This uses trigonometry and similar right
    // triangles.
    double offset = rmax / std::tan(coneAngle) - std::fabs(hmax);
    Vector3<double> coneVertex = box.center - offset * coneAxis;
    // coneVertex = {-6.8755829187191039e-09, 2.7364706271084139e-08, 1.0000000564305136}
    // coneAxis = {5.6712562731620983e-09, -2.2447173221576151e-08, -0.99999999999999967}
    // coneAngle = 0.46364759654859966

    ApprCone3<double> fitter{};
    size_t const maxIterations = 4;
    double const updateLengthTolerance = 1e-04;
    double const errorDifferenceTolerance = 1e-08;
    double const lambdaFactor = 0.001;
    double const lambdaAdjust = 10.0;
    size_t const maxAdjustments = 8;
    // Use the cone parameter estimates computed previously.
    bool useConeInputAsInitialGuess = true;

    int32_t numPoints = static_cast<int32_t>(points.size());
    auto const* conePoints = points.data();
    auto result = fitter(numPoints, conePoints, maxIterations, updateLengthTolerance,
        errorDifferenceTolerance, lambdaFactor, lambdaAdjust, maxAdjustments, useConeInputAsInitialGuess,
        coneVertex, coneAxis, coneAngle);
    // result.minLocation[0] = -8.5847082075293818e-09
    // result.minLocation[1] = 3.1845565499156142e-08
    // result.minLocation[2] = 1.0000000564208569
    // result.minLocation[3] = 2.5189507820619284e-09
    // result.minLocation[4] = -1.5077215317894200e-08
    // result.minLocation[5] = -1.1180339817025293
    // result.minError = 1.9785953518769534e-13
    // result.minErrorDifference = 2.5897897369161533e-13
    // result.minUpdateLength = 1.1747439061655418e-08
    // result.numIterations = 1
    // result.numAdjustments = 0
    // result.converged = true
    //
    // coneVertex = {-8.5847082075293818e-09, 3.1845565499156142e-08, 1.0000000564208569}
    // coneAxis = {2.2530180864682656e-09, -1.3485471429888732e-08, -1.0000000000000000}
    // coneAngle = 0.46364759639409514

[ruellm]

I have sent you the test files thru your email and also made the folder public, folder link is in the email as well. the test cases can also be downloaded here https://file.io/DjngfTC8ro6L (expires in 2 weeks)

[davideberly]

Thanks for the files. The file cone6.stl does not have the "wrong direction" cone axis problem, but the Levenberg-Marquardt minimizer does not appear to update the initial estimate the way I expected. Part of the problem are all the parameters you need to pass to the cone fitter, and the parameters control the minimizer. It is not clear how to choose the parameters. Numerical minimizers are more of an art than science.

I have been reading the cone-fitting section in my least-squares PDF. The mathematical equations involving the integrals are correct--I computed them with Mathematica previously but re-verified them today. A statement I made in the PDF is that several matrices and vectors, when estimated accurately, should produce a good fit. It is not clear whether this is happening. I need to verify additional equations and framework to see whether there is a mathematical error or implementation error. Finally, if what I have looks good mathematically, floating-point rounding errors could be causing problems. I might have to reformulate the fitting framework and try something different.

Fitting quadric surfaces to point sets is a challenging problem. Years ago I had a contract to fit ellipsoids to data (arose in 3D magnetic compass calibration). The points lived on a small patch of the surface (not uniformly distributed on an entire ellipsoid). Very small changes in the fitted parameters did not change the errors much near the patch, but the corresponding ellipsoids were of significantly different sizes.

The bottom line is that this changes my estimate for when a "fix" is ready. I do not have an estimate at the moment. Sorry.

[davideberly]

I found an error in the mathematical derivation of Section 8.1.3 Attempt to Reconstruct the Cone Vertex (of LeastSquaresFitting.pdf). This is most likely the cause of the incorrect vertex in your cone6.stl dataset.

[owai1980]

Congratulations!

Le ven. 9 juin 2023 à 21:24, David Eberly @.***> a écrit :

I found an error in the mathematical derivation of Section 8.1.3 Attempt to Reconstruct the Cone Vertex (of LeastSquaresFitting.pdf). This is most likely the cause of the incorrect vertex in your cone6.stl dataset.

— Reply to this email directly, view it on GitHub https://github.com/davideberly/GeometricTools/issues/63#issuecomment-1585036285, or unsubscribe https://github.com/notifications/unsubscribe-auth/AG3BAQLKZYO3OXHU64GS66TXKNZ5PANCNFSM6AAAAAAY37ZSMI . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[davideberly]

I have posted an update to ApprCone3.h and to LeastSquaresFitting.pdf, the latter file to my website. The PDF had an algorithm for estimating the cone vertex and cone angle, which required computing roots of a degree-4 polynomial, but the mathematical derivation was wrong. The correct derivation requires computing roots to a degree-10 polynomial.

I removed the description from the PDF and replaced it with a description similar to what I mentioned previously. The new algorithm is simpler and fixes the problem with the cone direction selected as U but needed to be -U.

Your cone.stl and cone6.stl files now have reasonable cone fits. The file where you limited the number of points (to 35) is not fit properly. As I had mentioned, if you have points not approximately uniformly distributed on the cone frustum, it is a very difficult problem to fit the points. Your limited set has points on approximately half a frustum. You might want to consider ways to choose your own initial cone given any information you have about the generation of your point set.

There is nothing I can do about the limited-point case for now. My cylinder fitting code reduces the minimizer to a complicated function of a unit-length vector. I sample unit-length vectors on a hemisphere and evaluate the error function, choosing the vector that produces the smallest error. The same idea can be applied to the cone fitting, but my investigation of this convinces me this will take a lot of development time. I do not have that time right now.

[ruellm]

Thank you for your help @davideberly really appreciated!