math::Sphere::OptimalEnclosingSphere() breaks down for torus with sufficient detail

[tksuoran]

Torus that shows a problem:

• major axis = 0.5
• minor axis = 0.125
• major axis steps = 40
• minor axis steps = 32 This gives a bounding sphere that is too large to be optimal (center = (0.619, 0.0, 0.048), radius = 1.247):

Meanwhile, torus with less detail:

• major axis = 0.5
• minor axis = 0.125
• major axis steps = 30
• minor axis steps = 24 This gives the expected bounding sphere (center = (0,0,0), radius = 0.625)

[juj]

Oh, very interesting!

Would you be able to dump the good and bad input vertices into a single .cpp test case file to showcase the error?

[juj]

First suspect might be to play around with the epsilon value at

https://github.com/juj/MathGeoLib/blob/55053da5e3e55a83043af7324944407b174c3724/src/Geometry/Sphere.cpp#L325-L326

which was needed to improve the numerical stability of the algorithm.

[tksuoran]

Indeed, 1e-3f works, even with further more detail. Thanks!

[tksuoran]

I can try to make a repro case some time later.

[juj]

Np.. another possibility to try if numerical stability is giving grief can be to duplicate the function into a version that uses doubles.

[tksuoran]

Repro:

#include <stdio.h>
#include <stdlib.h>

#include "../src/MathGeoLib.h"
#include "../src/Geometry/Sphere.h"
#include "../src/Math/myassert.h"
#include "TestRunner.h"
#include "TestData.h"

MATH_IGNORE_UNUSED_VARS_WARNING

MATH_BEGIN_NAMESPACE

std::vector<math::vec> GenerateTorus(
    double major_radius,
    double minor_radius,
    int    major_axis_steps,
    int    minor_axis_steps
)
{
    std::vector<math::vec> points;
    constexpr double pi = 3.14159265358979323846264338327950288;
    const double R = major_radius;
    const double r = minor_radius;
    for (int major = 0; major < major_axis_steps; ++major)
    {
        const double rel_major = static_cast<double>(major) / static_cast<double>(major_axis_steps);
        for (int minor = 0; minor < minor_axis_steps; ++minor)
        {
            const double rel_minor = static_cast<double>(minor) / static_cast<double>(minor_axis_steps);
            const double theta     = (pi * 2.0 * rel_major);
            const double phi       = (pi * 2.0 * rel_minor);
            const double sin_theta = std::sin(theta);
            const double cos_theta = std::cos(theta);
            const double sin_phi   = std::sin(phi);
            const double cos_phi   = std::cos(phi);
            const double x         = (R + r * cos_phi) * cos_theta;
            const double z         = (R + r * cos_phi) * sin_theta;
            const double y         =      r * sin_phi;
            const float  fx        = static_cast<float>(x);
            const float  fy        = static_cast<float>(y);
            const float  fz        = static_cast<float>(z);
            points.push_back(math::vec{fx, fy, fz});
        }
    }
    return points;
}

UNIQUE_TEST(BoundingSphere)
{
    const auto points = GenerateTorus(0.5, 0.125, 40, 32);

    const math::Sphere sphere = math::Sphere::OptimalEnclosingSphere(
        points.data(),
        static_cast<int>(points.size())
    );
    assert(sphere.r < (0.5 + 0.125 + 0.001));
}

MATH_END_NAMESPACE