dc-terrain / libdc / include / dc_qef.h

#ifndef DC_QEF_H
#define DC_QEF_H

// QEF (Quadratic Error Function) solver for Dual Contouring.
// Pure C++ — no Godot dependencies. Uses plain float arrays so it
// compiles standalone with g++ for testing.
//
// Given N intersection planes (point_i, normal_i), finds the position x
// that minimises  E(x) = sum_i (normal_i . (x - point_i))^2
//
// Internally accumulates the normal equations  (A^T A) x = A^T b  and
// solves via pseudoinverse (Jacobi SVD on the 3x3 symmetric A^T A).

struct QEFData {
    // Upper triangle of the symmetric 3x3 matrix A^T A.
    //   [0]=xx  [1]=xy  [2]=xz  [3]=yy  [4]=yz  [5]=zz
    float ata[6];

    // Right-hand side vector A^T b  (length 3).
    float atb[3];

    // Scalar b^T b — used for error computation.
    float btb;

    // Running average of all intersection points added so far.
    float masspoint[3];

    // Number of intersection points accumulated.
    int point_count;

    // Zero-initialise everything.
    void clear();

    // Add one intersection plane (point, normal).  normal need not be
    // unit-length but results are best when it is.
    void add(const float point[3], const float normal[3]);

    // Merge another QEFData into this one (for octree LOD simplification).
    void accumulate(const QEFData& other);

    // Solve (A^T A)^{-1} A^T b via SVD pseudoinverse.
    // Writes the optimal position into result[3].
    // For degenerate cases the masspoint supplies the missing degrees of
    // freedom — the solver projects the masspoint onto the solution
    // subspace so the result is always well-defined.
    void solve(float result[3]) const;

    // Solve with Tikhonov regularization (Boris the Brave center bias).
    // Adds bias_weight * ||x - center||^2 to the QEF, pulling the solution
    // toward cell center when the system is ill-conditioned (e.g., shallow
    // slopes where all normals are nearly colinear). This prevents the
    // solution from flying outside the cell and avoids stair-stepping
    // artifacts from hard clamping.
    void solve_biased(float result[3], const float center[3], float bias_weight) const;

    // Evaluate E(x) = x^T ATA x - 2 x^T ATB + BTB  at a given point.
    float error(const float point[3]) const;
};

#endif // DC_QEF_H