Order-related issue of MinimumVolumeBox3

[AnsonInTheAstroworld]

Hello, We are trying to use Geometric Tools on our project, but it appears to go wrong when using MiniVolumeBox3 to get the OBB of a point set. However, when I reverse the array of points, it goes well.

    for (size_t i = 0; i < mNumVertices; ++i)
    {
        size_t vLocalMax = vMax;
        ComputeType dLocalMax = dMax;
        size_t const* adjacent = &mAdjacentPool[mVertexAdjacent[vMax]];
        size_t numAdjacent = *adjacent++;
        for (size_t j = 1; j <= numAdjacent; ++j)
        {
            size_t vCandidate = *adjacent++;
            ComputeType dCandidate = Dot(direction, mVertices[vCandidate]);
            if (dCandidate > dLocalMax)
            {
                vLocalMax = vCandidate;
                dLocalMax = dCandidate;
            }
        }
        if (vMax != vLocalMax)
        {
            vMax = vLocalMax;
            dMax = dLocalMax;
        }
        else
        {
            break;
        }
    }

I suppose the problem is in GetExtreme function. It exits early when the local maximum dot product of adjacent points to input direction is equal to dmax, so it can't visit the real extreme vertex. I don't understand why the iteration continues only when the local maximum is greater than dmax, since I don't see it monotonically increasing in the adjacency. I tried to use std queue to handle it, but it takes too long, even much longer than exhausted method. I don't have other solution for now, so do you have any suggestions?

[davideberly]

Are you computing using a floating-point type or a rational type BSNumber? As with many computational geometry problems, floating-point rounding errors can cause misclassifications of signs that lead to incorrect results. The minimum-volume-box code computes the convex hull of the input points, and the hull algorithm is definitely affected by rounding errors.

If you want to provide your dataset that fails, I am happy to look at it and determine whether or not this is a floating-point rounding-error problem. The rational version has been used extensively by other users, and I have not yet received bug reports when using rational.

[AnsonInTheAstroworld]

One of the test case is listed below.

    std::vector<gte::Vector3<double>> pointSet = {
        {242.469421,577.558533,30.638069},
        {242.469421,577.558533,5.895439},
        {192.984161,627.043823,5.895439},
        {192.984161,602.301147,5.895439},
        {217.726807,627.043823,55.380684},
        {192.984161,627.043823,55.380684},
        {242.469421,577.558533,5.895439},
        {217.726807,577.558533,5.895439},
        {217.726807,602.301147,5.895439},
        {192.984161,602.301147,5.895439},
        {192.984161,627.043823,30.638069},
        {192.984161,627.043823,5.895439},
        {217.726807,602.301147,5.895439},
        {217.726807,577.558472,5.895439},
        {192.984161,627.043823,55.380684},
        {192.984161,627.043823,30.638054},
        {217.726807,627.043823,55.380684},
        {217.726807,602.301147,55.380684},
        {242.469421,602.301147,30.638054},
        {242.469421,577.558472,30.638054},
        {242.469421,602.301147,55.380684},
        {217.726807,602.301147,55.380684},
        {242.469421,602.301147,55.380684},
        {242.469421,602.301147,30.638054},
    };

    {
        gte::OrientedBox3<double> obb;
        double vol;
        gte::MinimumVolumeBox3<double, false> solver(8);
        solver(pointSet.size(), pointSet.data(), 8, obb, vol);
        printf("OBB vol: %f\n", vol);
    }

    std::reverse(pointSet.begin(), pointSet.end());

    {
        gte::OrientedBox3<double> obb;
        double vol;
        gte::MinimumVolumeBox3<double, false> solver(8);
        solver(pointSet.size(), pointSet.data(), 8, obb, vol);
        printf("Reversed OBB vol: %f\n", vol);
    }

And the output is

    OBB vol: 0.000000
    Reversed OBB vol: 96525.547185

Using double or rational number doesn't really matter in this case. The result of reversed point array is slightly different for double numbers, but for unreversed one it is still 0. Besides, we are building OBB with voxels, so the points are selected from axis-aligned edges, and it's not always a manifold mesh. Is it has to be a manifold one?

[davideberly]

Thank you for a clearly defined experiment to show that indeed there is a problem with the minimum-box code. It does not matter how you generate the input points, so there are no constraints on them. I will figure out how to fix the code and once fixed, respond in this thread. I appreciate your taking the time to investigate and report this.

[davideberly]

I posted a fix for GetExtreme. Your dataset has coplanar triangles sharing the initial vertex in the extreme-point search, and all vertices adjacent to the initial vertex have the same projection value. As you noticed, the GetExtreme function exits, thinking it has found the extreme point. The hill climbing had to be modified to allow the search to continue through the "flat region" until it finds vertices of larger projection value.

When you run the rational version, you get exactly the same output for the volume whether you use the original points or the reversed points.

[AnsonInTheAstroworld]

It works perfect now. Thank you for your help!👍 👍

[AnsonInTheAstroworld]

The issue appears again. Try this point set. {346.499969, 277.199982, 1068.900024}, { 346.499969, 277.199982, 1045.799927 }, { 346.499969, 277.199982, 1045.799927 }, { 323.399994, 277.199982, 1045.799927 }, { 346.499969, 277.199982, 1068.900024 }, { 323.399994, 277.199982, 1068.900024 }, { 323.399994, 277.199982, 1068.900024 }, { 323.399994, 277.199982, 1045.799927 }

I read the GetExtreme codes again, I think the latest commit just simply skips to the last one of all qualified candidates. So I modified it with queue to make sure it can visit every qualified candidates.

std::vector<uint32_t> visited(mNumVertices, 0);
size_t vMax = 0;
dMax = Dot(direction, mVertices[vMax]);
std::vector<size_t> vQue;
vQue.emplace_back(vMax);
while(!vQue.empty())
{
    std::vector<size_t> nxtQue;
    for(const auto vCur : vQue)
    {
        for (size_t i = 0; i < mNumVertices; ++i)
        {
            size_t const* adjacent = &mAdjacentPool[mVertexAdjacent[vCur]];
            size_t numAdjacent = *adjacent++;
            for (size_t j = 1; j <= numAdjacent; ++j)
            {
                size_t vCandidate = *adjacent++;
                if (visited[vCandidate] == 0)
                {
                    ComputeType dCandidate = Dot(direction, mVertices[vCandidate]);
                    if (dCandidate >= dMax)
                    {
                        vMax = vCandidate;
                        dMax = dCandidate;
                        nxtQue.emplace_back(vCandidate);
                    }
                    visited[vCandidate] = 1;
                }
            }
        }
    }
    vQue = nxtQue;
}
return vMax;

It works fine in my test cases, and slightly puts on time cost. Maybe you can try more cases to test its robustness.

[davideberly]

I can certainly analyze GetExtreme for your dataset, but before doing so, would you mind clarifying what you are expecting. The y-value of all your points are the same, so the points are coplanar and necessarily the box has a volume of 0. The output extent[2] is 0. The OBBs are the same whether I use the points in your specifed order or in reverse order.

[davideberly]

Regarding your test dataset, I should have mentioned that the points are coplanar which leads to a 2D convex hull. The execution path in MinimumVolumeBox3::operator() has 'dimension 2', in which case there is a call to MinimumAreaBox2. The function GetExtreme is never called.

[AnsonInTheAstroworld]

Sorry I gave a wrong one. You can try this again: {{-438.900024,-138.600006,75.599960,}, {-462.000000,-138.600006,75.599960,}, {-462.000031,-138.600006,98.699966,}, {-462.000031,-138.600006,75.599960,}, {-485.100037,-115.500000,6.299975,}, {-508.200012,-115.500000,6.299975,}, {-438.899994,-138.600006,121.800003,}, {-438.899994,-138.600006,98.699997,}, {-438.899994,-69.299995,52.499989,}, {-438.899994,-92.400002,52.499989,}, {-438.899994,-138.600006,75.599991,}, {-438.899994,-138.600006,52.499989,}, {-438.899994,-115.500000,121.800003,}, {-438.899994,-138.600006,121.800003,}, {-438.900024,-138.600006,98.699966,}, {-462.000000,-138.600006,98.699966,}, {-531.299988,-92.399994,29.399975,}, {-531.299988,-92.399994,6.299975,}, {-438.899994,-115.500000,52.499989,}, {-438.899994,-138.600006,52.499989,}, {-438.900024,-115.500000,52.499989,}, {-462.000000,-115.500000,52.499989,}, {-462.000031,-115.500000,52.499950,}, {-462.000031,-115.500000,29.399952,}, {-531.299988,-115.500000,29.399975,}, {-531.299988,-115.500000,6.299975,}, {-438.900024,-115.500000,121.800003,}, {-462.000000,-115.500000,121.800003,}, {-462.000031,-115.500000,29.399975,}, {-462.000031,-115.500000,6.299975,}, {-462.000031,-92.399994,52.499950,}, {-462.000031,-92.399994,29.399952,}, {-462.000031,-92.399994,121.800003,}, {-462.000031,-115.500000,121.800003,}, {-438.900024,-92.399994,52.499989,}, {-462.000000,-92.399994,52.499989,}, {-462.000031,-115.500000,6.299975,}, {-485.100006,-115.500000,6.299975,}, {-438.900024,-92.399994,121.800003,}, {-462.000000,-92.399994,121.800003,}, {-438.899994,-69.299995,121.800003,}, {-438.899994,-92.400002,121.800003,}, {-508.200012,-115.500000,6.299975,}, {-531.299988,-115.500000,6.299975,}, {-531.299988,-115.500000,98.699966,}, {-531.299988,-115.500000,75.599960,}, {-438.899994,-69.299995,121.800003,}, {-438.899994,-69.299995,98.699997,}, {-438.900024,-69.299995,98.699966,}, {-462.000000,-69.299995,98.699966,}, {-531.299988,-115.500000,52.499950,}, {-531.299988,-115.500000,29.399952,}, {-462.000031,-69.299995,98.699966,}, {-462.000031,-69.299995,75.599960,}, {-438.899994,-69.299995,75.599991,}, {-438.899994,-69.299995,52.499989,}, {-531.299988,-115.500000,75.599991,}, {-531.299988,-115.500000,52.499989,}, {-438.900024,-69.299995,75.599960,}, {-462.000000,-69.299995,75.599960,}, {-531.299988,-92.399994,98.699966,}, {-531.299988,-115.500000,98.699966,}, {-531.299988,-92.399994,98.699966,}, {-531.299988,-92.399994,75.599960,}, {-531.299988,-92.399994,75.599991,}, {-531.299988,-92.399994,52.499989,}, {-531.299988,-92.399994,52.499950,}, {-531.299988,-92.399994,29.399952,}, {-508.200012,-92.399994,6.299975,}, {-531.299988,-92.399994,6.299975,}, {-462.000031,-92.399994,29.399975,}, {-462.000031,-92.399994,6.299975,}, {-462.000031,-92.399994,6.299975,}, {-485.100006,-92.399994,6.299975,}, {-485.100037,-92.399994,6.299975,}, {-508.200012,-92.399994,6.299975,}} This point set should generate an OBB degraded to 2d. However, if change the element order as below, it would generate a correct OBB. {{-438.900024,-138.600006,75.599960}, {-462.000000,-138.600006,75.599960}, {-462.000031,-138.600006,98.699966}, {-462.000031,-138.600006,75.599960}, {-438.899994,-138.600006,121.800003}, {-438.899994,-138.600006,98.699997}, {-485.100037,-115.500000,6.299975}, {-508.200012,-115.500000,6.299975}, {-438.899994,-138.600006,75.599991}, {-438.899994,-138.600006,52.499989}, {-438.899994,-69.299995,52.499989}, {-438.899994,-92.400002,52.499989}, {-438.899994,-115.500000,121.800003}, {-438.899994,-138.600006,121.800003}, {-438.900024,-138.600006,98.699966}, {-462.000000,-138.600006,98.699966}, {-438.899994,-115.500000,52.499989}, {-438.899994,-138.600006,52.499989}, {-531.299988,-92.399994,29.399975}, {-531.299988,-92.399994,6.299975}, {-438.900024,-115.500000,52.499989}, {-462.000000,-115.500000,52.499989}, {-462.000031,-115.500000,52.499950}, {-462.000031,-115.500000,29.399952}, {-438.900024,-115.500000,121.800003}, {-462.000000,-115.500000,121.800003}, {-531.299988,-115.500000,29.399975}, {-531.299988,-115.500000,6.299975}, {-462.000031,-115.500000,29.399975}, {-462.000031,-115.500000,6.299975}, {-462.000031,-92.399994,121.800003}, {-462.000031,-115.500000,121.800003}, {-462.000031,-92.399994,52.499950}, {-462.000031,-92.399994,29.399952}, {-462.000031,-115.500000,6.299975}, {-485.100006,-115.500000,6.299975}, {-438.900024,-92.399994,52.499989}, {-462.000000,-92.399994,52.499989}, {-438.900024,-92.399994,121.800003}, {-462.000000,-92.399994,121.800003}, {-438.899994,-69.299995,121.800003}, {-438.899994,-92.400002,121.800003}, {-508.200012,-115.500000,6.299975}, {-531.299988,-115.500000,6.299975}, {-438.899994,-69.299995,121.800003}, {-438.899994,-69.299995,98.699997}, {-531.299988,-115.500000,98.699966}, {-531.299988,-115.500000,75.599960}, {-438.900024,-69.299995,98.699966}, {-462.000000,-69.299995,98.699966}, {-531.299988,-115.500000,52.499950}, {-531.299988,-115.500000,29.399952}, {-462.000031,-69.299995,98.699966}, {-462.000031,-69.299995,75.599960}, {-531.299988,-115.500000,75.599991}, {-531.299988,-115.500000,52.499989}, {-438.899994,-69.299995,75.599991}, {-438.899994,-69.299995,52.499989}, {-438.900024,-69.299995,75.599960}, {-462.000000,-69.299995,75.599960}, {-531.299988,-92.399994,98.699966}, {-531.299988,-115.500000,98.699966}, {-531.299988,-92.399994,98.699966}, {-531.299988,-92.399994,75.599960}, {-531.299988,-92.399994,75.599991}, {-531.299988,-92.399994,52.499989}, {-531.299988,-92.399994,52.499950}, {-531.299988,-92.399994,29.399952}, {-462.000031,-92.399994,29.399975}, {-462.000031,-92.399994,6.299975}, {-508.200012,-92.399994,6.299975}, {-531.299988,-92.399994,6.299975}, {-462.000031,-92.399994,6.299975}, {-485.100006,-92.399994,6.299975}, {-485.100037,-92.399994,6.299975}, {-508.200012,-92.399994,6.299975}} If use the GetExtreme function I modified above, two data will generate same OBB.

[davideberly]

The GetExtreme function had another bug. Before using your queue approach, I am hoping the performance is better by fixing this bug. If you encounter another bug, please re-open the issue.

Your illustrative example, when run single-threaded, showed the failure on the third call to GetExtreme. The dMax value is computed for mVertices[0] and is 0. All vertices adjacent to mVertices[0] have d-values of 0. The loop then processes mVertices[1], and all its adjacent vertices have d-values of 0. The last of those vertices is mVertices[0]. This causes vMax to be set to 0, which invalidates the search. It is necessary before the loop to mark mVertices[0] as visited.

I have posted the fix. Once again, thank you for reporting the problem and providing the dataset!

[AnsonInTheAstroworld]

The bug appears again. Try the point set below. {-88.704002,205.128021, 546.839966}, {-99.792000,205.128021, 546.839966}, {-77.616005,205.128021, 546.839966}, {-88.704002,205.128021, 546.839966}, {-99.792007,216.216019, 546.839966}, {-99.792007,205.128006, 546.839966}, {-77.616005,216.216019, 546.839966}, {-77.616005,205.128006, 546.839966}, {-99.792007,216.216019, 546.839966}, {-110.880005,216.216019, 546.839966}, {-110.880005,216.216019, 569.015991}, {-110.880005,205.128006, 569.015991}, {-110.880005,216.216019, 557.927979}, {-110.880005,216.216019, 546.839966}, {-66.528008,216.216019, 546.839966}, {-77.616005,216.216019, 546.839966}, {-66.528008,216.216019, 557.927979}, {-66.528008,216.216019, 546.839966}, {-110.880005,216.216019, 557.927979}, {-110.880005,205.128006, 557.927979}, {-66.528008,216.216019, 557.927979}, {-66.528008,205.128006, 557.927979}, {-110.880005,205.128021, 569.015991}, {-110.880005,205.128021, 557.927979}, {-66.528008,205.128021, 569.015991}, {-66.528008,205.128021, 557.927979}, {-66.528008,216.216019, 569.015991}, {-66.528008,205.128006, 569.015991}, {-110.880005,216.216019, 591.192017}, {-110.880005,216.216019, 580.104004}, {-110.880005,216.216019, 580.104004}, {-110.880005,216.216019, 569.015991}, {-66.528008,216.216019, 580.104004}, {-66.528008,216.216019, 569.015991}, {-66.528008,216.216019, 591.192017}, {-66.528008,216.216019, 580.104004}, {-99.792007,216.216019, 591.192017}, {-110.880005,216.216019, 591.192017}, {-66.528008,216.216019, 591.192017}, {-77.616005,216.216019, 591.192017}, {-88.704002,216.216019, 591.192017}, {-99.792000,216.216019, 591.192017}, {-77.616005,216.216019, 591.192017}, {-88.704002,216.216019, 591.192017} There is 0 axis extent in result obb, while it is correct if change its order as below. {-88.704002,205.128021, 546.839966}, {-99.792000,205.128021, 546.839966}, {-77.616005,205.128021, 546.839966}, {-88.704002,205.128021, 546.839966}, {-99.792007,216.216019, 546.839966}, {-99.792007,205.128006, 546.839966}, {-99.792007,216.216019, 546.839966}, {-110.880005,216.216019, 546.839966}, {-77.616005,216.216019, 546.839966}, {-77.616005,205.128006, 546.839966}, {-110.880005,216.216019, 557.927979}, {-110.880005,216.216019, 546.839966}, {-110.880005,216.216019, 569.015991}, {-110.880005,205.128006, 569.015991}, {-66.528008,216.216019, 546.839966}, {-77.616005,216.216019, 546.839966}, {-66.528008,216.216019, 557.927979}, {-66.528008,216.216019, 546.839966}, {-110.880005,216.216019, 557.927979}, {-110.880005,205.128006, 557.927979}, {-66.528008,216.216019, 557.927979}, {-66.528008,205.128006, 557.927979}, {-110.880005,205.128021, 569.015991}, {-110.880005,205.128021, 557.927979}, {-66.528008,205.128021, 569.015991}, {-66.528008,205.128021, 557.927979}, {-66.528008,216.216019, 569.015991}, {-66.528008,205.128006, 569.015991}, {-110.880005,216.216019, 580.104004}, {-110.880005,216.216019, 569.015991}, {-110.880005,216.216019, 591.192017}, {-110.880005,216.216019, 580.104004}, {-66.528008,216.216019, 580.104004}, {-66.528008,216.216019, 569.015991}, {-66.528008,216.216019, 591.192017}, {-66.528008,216.216019, 580.104004}, {-99.792007,216.216019, 591.192017}, {-110.880005,216.216019, 591.192017}, {-66.528008,216.216019, 591.192017}, {-77.616005,216.216019, 591.192017}, {-88.704002,216.216019, 591.192017}, {-99.792000,216.216019, 591.192017}, {-77.616005,216.216019, 591.192017}, {-88.704002,216.216019, 591.192017} Still the queue method works fine for both.

[AnsonInTheAstroworld]

The GetExtreme function had another bug. Before using your queue approach, I am hoping the performance is better by fixing this bug. If you encounter another bug, please re-open the issue.

Your illustrative example, when run single-threaded, showed the failure on the third call to GetExtreme. The dMax value is computed for mVertices[0] and is 0. All vertices adjacent to mVertices[0] have d-values of 0. The loop then processes mVertices[1], and all its adjacent vertices have d-values of 0. The last of those vertices is mVertices[0]. This causes vMax to be set to 0, which invalidates the search. It is necessary before the loop to mark mVertices[0] as visited.

I have posted the fix. Once again, thank you for reporting the problem and providing the dataset!

I think I can't reopen an issue closed by you in GitHub. I'll open a new one if you can't see my comments.

[davideberly]

I must be missing something. The hill-climbing seems like it should work without the queue. Regardless, I will read through your implementation to make sure I understand it and verify that it is a correct algorithm.

[davideberly]

I replaced the code by a simple iteration over all vertices. Intel's oneAPI VTune profiler shows that this performs slightly better than the queue-based approach. The check-in comments have the profiling details. Thank you yet one more time for being persistent about having me fix my errors.

As a bonus, the performance was awful because of UIntegerFP32 spending an enormous amount of time initializing the large mBits arrays (because of mBits{} in the initializer list). I removed those initializations, which led to a large increase in performance.

[davideberly]

Closing inactive issue. Without feedback, assuming the problem has been corrected.