estimate_normals throws when receiving point inputs belonging to the three origin planes

[cakarsubasi]

I ran into this specific problem when trying to generate sampled mathematical surfaces for testing. Quite possibly the simplest mathematical surface that I can think of is the plane, so I tried generating points that satisfy the equation z = 0 (i.e. arbitrary points). This input crashes rotation system reconstruction during the normal generation step because the covariance matrix given to smallestEigenVector is singular. It also crashes for y = 0 and x = 0 but not for other planes as far as I can tell.

Although it seems extremely unlikely that this case is reached in any "real input", it is also the simplest input that came to my mind.

An example input triggering the crash is provided below:

math_points.obj

[janba]

Ok. This must be fixed, and it should be relatively easy to fix. I think the correct behaviors is that the normal is computed as the cross product of the two largest eigenvectors and the smallest is ignored.

[cakarsubasi]

Seemingly, the current behavior is that the smallest eigenvector is returned. Although the problem is apparently that the covariance matrix is singular?

template <std::ranges::input_range Range>
CGLA::Vec3d estimateNormal(Range&& neighbors)
{
    CGLA::Vec3d centroid(0.0f);
    for (const auto& point : neighbors) {
        centroid += point;
    }
    centroid /= static_cast<double>(neighbors.size());

    CGLA::Mat3x3d covariance(0.0f);
    for (const auto& point : neighbors) {
        CGLA::Vec3d diff = point - centroid;
        covariance += outer_product(diff, diff);
    }
    covariance /= static_cast<double>(neighbors.size());

    // inlined smallestEigenVector function
    Mat3x3d inverseMatrix = invert(covariance);
    Mat3x3d eigenvectors(0);
    Mat3x3d eigenvalues(0);

    power_eigensolution<Mat3x3d>(inverseMatrix, eigenvectors, eigenvalues, 1);
    // end inlined smallestEigenVector function

    const CGLA::Vec3d normal = eigenvectors[0];


    return normalize(normal); // Normalize to ensure unit length
}

[janba]

Right. This code inverts the covariance matrix which preserves eigenvectors but inverts eigenvalues. With this approach the smallest eigenvalue becomes the greatest, and the one returned by power_eigensolution. This works well, but if the matrix is singular the inversion fails. However, even for a singular matrix, power eigensolution should work for the un-inverted matrix as long as you ask it only for the two largest eigensolutions. That should fail only if there are zero or one solutions, but then we do not have a surface, so that is ok.

[cakarsubasi]

I tried this modification:

    CGLA::Mat3x3d eigenvectors(0);
    CGLA::Mat3x3d eigenvalues(0);

    auto sols = CGLA::power_eigensolution<CGLA::Mat3x3d>(covariance, eigenvectors, eigenvalues, 2);
    // GEL_ASSERT(sols > 1);

    auto normal = CGLA::cross(eigenvectors[0], eigenvectors[1]);

This one works for the flat surface, but almost immediately fails for the smooth Stanford bunny.