murrayreadccdc commented 1 year ago

We are using MinimumVolumeBox3 around shapes with high symmetry centred near the origin. We find that sometimes the box returned does not include all of volume needed to include the points. Sometimes the box is terminated at the origin. The behaviour seems to be sensitive to the order of the input points. I have managed to reproduce this problem with the following test code:

#include <iostream>
#include <Mathematics/MinimumVolumeBox3.h>
#include <random>

void test_minimum_volume_box()
{
    std::vector< gte::Vector3<double> > vertices =
    {
        {-2.9523419,11.413461,-16.111856},
        {-12.314037,11.413461,10.801515},
        {-12.314037,-11.413461,10.801515},
        {-2.9523419,-11.413461,-16.111856},
        {12.314037,11.413461,-10.801515},
        {2.9523419,11.413461,16.111856},
        {2.9523419,-11.413461,16.111856},
        {12.314037,-11.413461,-10.801515},
        {-3.7572584,-17.136135,10.801515},
        {3.7572584,-17.136135,-10.801515},
        {1.9792732,-15.362445,-13.31444},
        {-3.7572584,17.136135,10.801515},
        {3.7572584,17.136135,-10.801515},
        {1.9792732,15.362445,-13.31444},
        {-1.9792732,-15.362445,13.31444},
        {-1.9792732,15.362445,13.31444},
        {-2.9523419,-11.413461,16.111856},
        {-2.9523419,11.413461,16.111856},
        {2.9523419,11.413461,-16.111856},
        {2.9523419,-11.413461,-16.111856}
    };

    std::mt19937 mte;
    for (int i = 0; i < 20; i++)
    {
        std::shuffle(vertices.begin(), vertices.end(), mte);

        MinimumVolumeBox3<double, true > gte_minimum_volume_box_calculator;

        size_t lgMaxSample = 5;
        OrientedBox3<double> oriented_box;
        double volume = 0.0;
        gte_minimum_volume_box_calculator(static_cast<int>(vertices.size()), &vertices[0], lgMaxSample, oriented_box, volume);

        std::cout << volume << std::endl;
    }
}

This is producing the following results when I run it:

Since the test finds the problem with random reordering of the input, it suggests a workaround I can use which would be to retry with random input ordering until I get a valid result.

murrayreadccdc commented 1 year ago

This is a smaller set of points describing a hexagonal prism that produce very unstable results. The y values close to but not quite zero seem to cause the problems here.

    {4.7581165, -8.2412995, 15.005264},
    {4.7581165, -8.2412995, -15.005264},
    {-4.7581165, -8.2412995, -15.005264},
    {-4.7581165, -8.2412995, 15.005264},
    {4.7581165, 8.2412995, 15.005264},
    {4.7581165, 8.2412995, -15.005264},
    {-4.7581165, 8.2412995, -15.005264},
    {-4.7581165, 8.2412995, 15.005264},
    {-9.516233, -9.9321042e-13, -15.005264},
    {-9.516233, -9.8865609e-13, 15.005264},
    {9.516233, 1.9292248e-12, 15.005264},
    {9.516233, 1.9264354e-12, -15.005264},

volumes:

davideberly commented 1 year ago

Hello Murray.

The minimum-volume box algorithm is a typical computational geometry problem where you need exact arithmetic to obtain the correct results. You have the ComputeDouble template parameter set to ‘true’, so the computations are using floating-point arithmetic.

Have you tried the exact arithmetic version? Set the ComputeDouble parameter to ‘false’. Yes, it runs slower than with double-precision arithmetic, but that is a trade-off common to such geometry problems.

Regardless, I will investigate with the data you sent and verify for myself that the problem is what I think it is.

-- Dave Eberly https://www.geometrictools.com

From: Murray Read @.> Sent: Monday, September 11, 2023 3:18 AM To: davideberly/GeometricTools @.> Cc: Subscribed @.***> Subject: [davideberly/GeometricTools] MinimumVolumeBox3 sometimes misses large parts of the required volume (Issue #72)

We are using MinimumVolumeBox3 around shapes with high symmetry centred near the origin. We find that sometimes the box returned does not include all of volume needed to include the points. Sometimes the box is terminated at the origin. The behaviour seems to be sensitive to the order of the input points. I have managed to reproduce this problem with the following test code:

include

include <Mathematics/MinimumVolumeBox3.h>

include

void test_minimum_volume_box()

{

std::vector< gte::Vector3<double> > vertices =

{

    {-2.9523419,11.413461,-16.111856},

    {-12.314037,11.413461,10.801515},

    {-12.314037,-11.413461,10.801515},

    {-2.9523419,-11.413461,-16.111856},

    {12.314037,11.413461,-10.801515},

    {2.9523419,11.413461,16.111856},

    {2.9523419,-11.413461,16.111856},

    {12.314037,-11.413461,-10.801515},

    {-3.7572584,-17.136135,10.801515},

    {3.7572584,-17.136135,-10.801515},

    {1.9792732,-15.362445,-13.31444},

    {-3.7572584,17.136135,10.801515},

    {3.7572584,17.136135,-10.801515},

    {1.9792732,15.362445,-13.31444},

    {-1.9792732,-15.362445,13.31444},

    {-1.9792732,15.362445,13.31444},

    {-2.9523419,-11.413461,16.111856},

    {-2.9523419,11.413461,16.111856},

    {2.9523419,11.413461,-16.111856},

    {2.9523419,-11.413461,-16.111856}

};

std::mt19937 mte;

for (int i = 0; i < 20; i++)

{

    std::shuffle(vertices.begin(), vertices.end(), mte);

    MinimumVolumeBox3<double, true > gte_minimum_volume_box_calculator;

    size_t lgMaxSample = 5;

    OrientedBox3<double> oriented_box;

    double volume = 0.0;

    gte_minimum_volume_box_calculator(static_cast<int>(vertices.size()), &vertices[0], lgMaxSample, oriented_box, volume);

    std::cout << volume << std::endl;

}

}

This is producing the following results when I run it:

17934.5

8967.26

17934.5

Since the test finds the problem with random reordering of the input, it suggests a workaround I can use which would be to retry with random input ordering until I get a valid result.

— Reply to this email directly, view it on GitHubhttps://github.com/davideberly/GeometricTools/issues/72, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AQ5SOLJKYDOYZEFIS36YSTDXZ3QNDANCNFSM6AAAAAA4TCSHUQ. You are receiving this because you are subscribed to this thread.Message ID: @.**@.>>

davideberly commented 1 year ago

I added the rational version to your code.

    std::mt19937 mte;
    for (int i = 0; i < 20; i++)
    {
        std::shuffle(vertices.begin(), vertices.end(), mte);

        MinimumVolumeBox3<double, true > gte_minimum_volume_box_calculator;

        size_t lgMaxSample = 5;
        OrientedBox3<double> oriented_box;
        double volume = 0.0;
        gte_minimum_volume_box_calculator(static_cast<int>(vertices.size()), &vertices[0], lgMaxSample, oriented_box, volume);
        std::cout << "double volume = " << volume << "; ";

        size_t constexpr numThreads = 8;
        MinimumVolumeBox3<double, false> mvb(numThreads);
        mvb(static_cast<int>(vertices.size()), vertices.data(), lgMaxSample, oriented_box, volume);
        std::cout << "rational volume = " << volume << std::endl;
    }

The volumes computed using rational arithmetic are the same regardless of the shuffle: 17934.5 for your first dataset and 9414.44 for your second dataset. The MinimumVolumeBox3 constructor allows you to specify number of threads for computing. In my modification to your code, I chose 8 (running on an Intel i9-10900 that has 10 cores [20 logical processors]).

So indeed the problem is that the output using double-precision is not always correct because of floating-point rounding errors.

For the record, in my GTL development track, I am removing the option to select floating-point for those geometric algorithms where you cannot ensure the correct results when using floating-point.

murrayreadccdc commented 1 year ago

Hi Dave,

Thanks for looking into this. I will experiment with the exact arithmetic setting to see if I can make it work with my code.

Rounding errors are hard to deal with. My work has a number of coordinate optimisation problems which have a habit of falling into different local minima depending on the order of calculation as rounding errors are magnified by iteration. But are these minimum volume box differences rounding errors? I didn't make it clear in the test code above, but some of the boxes I've seen returned only cover a fraction, maybe half, of the space covered by the input points. The boxes occasionally have faces that cut the origin rather than being symmetric around the origin, which is a property of the input points used.

To demonstrate the issue, this version of the test code tests the furthest that any of the input points lie beyond the extent of the box. You would expect the value to be 1, IE no points are outside the box, but some are on the edges of the box.

std::mt19937 mte;
for (int r = 0; r < 20; r++)
{
    std::shuffle(vertices.begin(), vertices.end(), mte);

    MinimumVolumeBox3<double, true > gte_minimum_volume_box_calculator;

    size_t lgMaxSample = 5;
    OrientedBox3<double> oriented_box;
    double volume = 0.0;
    gte_minimum_volume_box_calculator(static_cast<int>(vertices.size()), &vertices[0], lgMaxSample, oriented_box, volume);

    // Get axis transpose
    std::array< gte::Vector3<double>, 3> axis_tp;

    for (auto i = 0; i != 3; ++i)
    {
        for (auto j = 0; j != 3; ++j)
        {
            axis_tp[j][i] = oriented_box.axis[i][j];
        }
    }
    gte::Vector3<double> ct = (oriented_box.center[0] * axis_tp[0]) + (oriented_box.center[1] * axis_tp[1]) + (oriented_box.center[2] * axis_tp[2]);

    double max_outside = 0.0;
    for (const gte::Vector3<double>& v : vertices)
    {
        // convert the vertex to the box local coordinate system
        gte::Vector3<double> local = (v[0] * axis_tp[0]) + (v[1] * axis_tp[1]) + (v[2] * axis_tp[2]) - ct;
        for (int d = 0; d < 3; d++)
        {
            double outside = std::abs(local[d]) / oriented_box.extent[d];
            if (outside > max_outside)
                max_outside = outside;
        }
    }

    std::cout << volume << " " << max_outside << std::endl;
}

On running this with the test second test set it seems that some of the input points are up to 65% beyond the extent of the box.

9406.11 1.03437
7498.01 1.58377
9389.86 1.05908
9389.86 1.05908
9414.44 1
9414.44 1
9414.44 1
7498.01 1.58377
9381.15 1.06879
9414.44 1
9406.11 1.03437
9414.44 1
9405.08 1.03642
7498.01 1.58377
9407.45 1.03146
9405.08 1.03642
9389.86 1.05908
7132.4 1.65315
9414.44 1
9405.08 1.03642

davideberly commented 1 year ago

"But are these minimum volume box differences rounding errors?" Absolutely! In your declaration for MinimumVolumeBox3<double, true>, change the second template parameter to 'false'. Run your program. The reported volumes are all the same (9414.44). The second column of numbers are all the same (1).

My favorite example of how bad the errors can be is in my "Robust Geometric Computing" book, Section 1.1.3 A Nonrobust Floating-Point Implementation. If you do not have access to a copy, let me know and I can summarize the example here.

There are closed-form equations for the roots of a cubic polynomial. There is a number computed from the coefficients that is used for classifying whether there are complex roots. When it is positive, there are such roots. Using 'float', the number is 2.74877972e+11. When using rational arithmetic, the number is -8.807716e-08. (The root finder was in the DirectX math library. I have not checked recently to see whether it has been replaced by robust code.)

This is why in GTL I am preventing folks from selecting compute type as 'float' or 'double' when rational arithmetic is required for correct output.

murrayreadccdc commented 1 year ago

Thanks again Dave, I wasn't appreciating the subtleties of robust geometric computation. My application involves calculations on scientific data which is necessarily imprecise, so we are typically wanting valid maybe approximate calculation results rather than exact maybe invalid results. I will see what I can do with the rational mode calculation, though the warning on stack size have put me off so far. Thanks for everything you're doing!

davideberly commented 1 year ago

Perhaps you can view the exact rational fit to inputs with errors as “valid maybe approximate”?

I forgot to mention about stack size. Rationals use a lot of memory, stack and/or heap. If using Visual Studio, I set the “stack reserve size” to 256MB.

murrayreadccdc commented 1 year ago

I've been experimenting with the rational calculation, and fortunately when running under Python there does seem to be enough stack space to use this. This gives me an effective solution, which is to try the real number calculation with repetition to find a valid box, otherwise use the rational calculation. I'd welcome it if you ever decide to investigate real number based algorithms that always give valid answers. There is some interesting test data in this ticket. Otherwise I'm happy for this to be closed.

davideberly commented 1 year ago

I do not understand your reluctance to use the rational algorithm. Selecting random permutations of the input hoping to find the minimum volume box using floating-point and then falling back to the rational computations if you do not find one seems like more work than simply using the rational computations.

I already spent time over several years looking at floating-point-only algorithms. This is fruitless. You want accurate results? You have to use exact predicates involving rational arithmetic, or for smaller amortized cost: a mixture of interval arithmetic and rational arithmetic.

Sorry I am unable to help you with your problem. I wish you luck with your approach.

davideberly / GeometricTools

MinimumVolumeBox3 sometimes misses large parts of the required volume #72

include

include <Mathematics/MinimumVolumeBox3.h>

include