d2n_evaluator.cpp

/* ***********************************************************
 *
 * Description:
 *  d2n_evaluator.h/cpp computes evaluation metrics between an estimated
 *  normal map and the ground truth normal map.
 *
 *  Metrics are computed only on valid pixels (valid_mask != 0).
 *  
 * ----------------------------------------------------------
 *
 *  Copyright (C) 2026 Claudio Z. (cloudofoz)
 *
 *  Permission is hereby granted, free of charge, to any person obtaining a copy
 *  of this software and associated documentation files (the "Software"), to deal
 *  in the Software without restriction, including without limitation the rights
 *  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 *  copies of the Software, and to permit persons to whom the Software is
 *  furnished to do so, subject to the following conditions:
 *
 *  The above copyright notice and this permission notice shall be included in all
 *  copies or substantial portions of the Software.
 *
 *  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 *  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 *  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 *  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 *  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 *  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 *  SOFTWARE.
 *
 * **********************************************************/

#include "d2n_evaluator.h"

#include <algorithm>
#include <cmath>
#include <limits>
#include <raylib.h>

 //---------------------------------------------------------
 // Utility

static float clampf(float x, float a, float b)
{
    return std::max(a, std::min(b, x));
}

// Euclidean norm for a 3D vector (Vec3f).
static float compute_norm(const cv::Vec3f& v)
{
    return std::sqrt(v.dot(v));
}

// Note: modifies vector order (std::nth_element performs a partial reorder).
// Returns the "upper median" when v.size() is even
static float compute_median(std::vector<float>& v)
{
    const auto middleIndex = v.size() / 2;
    std::nth_element(v.begin(), v.begin() + middleIndex, v.end());
    return v[middleIndex];
}

// Variance via E[x^2] - (E[x])^2, where E is the mean over valid samples.
// Uses max(var, 0) for numerical safety (floating point errors).
static double compute_variance_from_moments(double sum, double sum2, int count)
{
    if (count <= 0)
        return 0.0;

    const double mean = sum / static_cast<double>(count);
    double var = (sum2 / static_cast<double>(count)) - mean * mean;
    return std::max(var, 0.0);
}

//---------------------------------------------------------

// Evaluates estimated normals against ground truth normals.
//
// - For each valid pixel, it compares the estimated normal with the GT normal.
// - The comparison is the angular distance (in degrees) between the two vectors.
// - It then aggregates those per-pixel errors into summary metrics (mean, median, std, thresholds).
//
// Note:
// - Only pixels where est.valid_mask != 0 are considered.
// - Optionally, normals can be normalized here.
// - Optionally, the angle can be "unsigned" (n and -n as equivalent).
// - Optionally, a per-pixel heatmap of angular error can be produced.
D2NEvaluatorMetrics d2n_evaluate(const D2NEstimatorOutput& est,
    const cv::Mat& gt_normals,
    const D2NEvaluatorParams& opt,
    cv::Mat* outAngleErrorHeatmap)
{
    // Sanity checks: normal maps (3 floats per pixel) and a valid mask (1 byte per pixel).
    CV_Assert(est.estimated_normals.type() == CV_32FC3);
    CV_Assert(est.valid_mask.type() == CV_8UC1);
    CV_Assert(gt_normals.type() == CV_32FC3);
    CV_Assert(est.estimated_normals.size() == gt_normals.size());
    CV_Assert(est.estimated_normals.size() == est.valid_mask.size());

    const int rows = est.estimated_normals.rows;
    const int cols = est.estimated_normals.cols;

    // Optional output: per-pixel angular error heatmap (in degrees).
    // Useful for debugging and qualitative analysis (where are the errors located).
    if (opt.output_error_heatmap)
    {
        // A valid pointer has to be provided.
        CV_Assert(outAngleErrorHeatmap != nullptr);

        outAngleErrorHeatmap->create(rows, cols, CV_32FC1);
        outAngleErrorHeatmap->setTo(0.0f); // keep invalid pixels at 0
    }

    // Store all per-pixel angles in a vector because the median requires sorting.
    // Reserve an estimate of the needed size to avoid repeated reallocations.
    std::vector<float> angles_deg;
    angles_deg.reserve(est.valid_pixels > 0 ? est.valid_pixels : rows * cols);

    // Accumulate sums for mean and standard deviation.
    // Store sum and sum of squares so we do not need a second pass.
    double sum = 0.0; // sum(angle)
    double sum2 = 0.0; // sum(angle^2)

    // Track min/max angular error.
    double minA = +std::numeric_limits<double>::infinity();
    double maxA = -std::numeric_limits<double>::infinity();

    // Simple "how many pixels are below threshold" counters.
    int count5 = 0, count10 = 0, count20 = 0;
    int valid = 0;

    for (int v = 0; v < rows; ++v)
    {
        // Row pointers for performance: evaluation runs per pixel and can be large.
        const cv::Vec3f* estPtr = est.estimated_normals.ptr<cv::Vec3f>(v);
        const cv::Vec3f* gtPtr = gt_normals.ptr<cv::Vec3f>(v);
        const uchar* maskPtr = est.valid_mask.ptr<uchar>(v);

        // If the heatmap is enabled, get the row pointer; otherwise keep it null.
        float* hmPtr = (opt.output_error_heatmap) ? outAngleErrorHeatmap->ptr<float>(v) : nullptr;

        for (int u = 0; u < cols; ++u)
        {
            // Skip invalid pixels.
            if (maskPtr[u] == 0)
                continue;

            cv::Vec3f a = estPtr[u];
            cv::Vec3f b = gtPtr[u];

            // Optional normalization:
            // Normal maps should already be unit vectors, this optionally increases robustness
            if (opt.normalize_estimated)
            {
                const float length = compute_norm(a);
                if (length > FLT_EPSILON) a /= length;
                else continue; // -> skip
            }

            if (opt.normalize_gt)
            {
                const float length = compute_norm(b);
                if (length > FLT_EPSILON) b /= length;
                else continue; // -> skip
            }

            // Angle between two unit vectors is computed from the dot product:
            // dot = cos(theta), so theta = acos(dot).
            //
            // Forcing [-1, +1] (floating points potential issues).
            float dot = clampf(a.dot(b), -1.0f, 1.0f);

            // Unsigned angle:
            // n and -n represent the same surface orientation
            // Taking abs(dot) makes theta in [0, 90] instead of [0, 180].
            if (opt.use_unsigned_angle)
                dot = std::fabs(dot);

            // Convert from radians to degrees to make the results easier to read.
            const float angRad = std::acos(dot);
            const float angDeg = angRad * RAD2DEG;

            // Track extremes.
            minA = std::min(minA, static_cast<double>(angDeg));
            maxA = std::max(maxA, static_cast<double>(angDeg));

            // Store optional per-pixel error.
            // Note: invalid pixels remain at 0 because the heatmap was initialized to 0.
            if (hmPtr)
                hmPtr[u] = angDeg;

            // Update statistics.
            sum += angDeg;
            sum2 += static_cast<double>(angDeg) * static_cast<double>(angDeg);

            // Threshold counters (strictly "below").
            if (angDeg < 5.0f)  ++count5;
            if (angDeg < 10.0f) ++count10;
            if (angDeg < 20.0f) ++count20;

            // Keep per-pixel values for median computation.
            angles_deg.push_back(angDeg);
            ++valid;
        }
    }

    D2NEvaluatorMetrics M;

    // If no valid pixels is valid, return the default metrics struct.
    if (valid == 0)
        return M;

    // Mean absolute angular error (in degrees).
    // (The mean of per-pixel angular errors.)
    M.mae_deg = sum / static_cast<double>(valid);

    // Standard deviation of the angular error distribution.
    // We compute it from the first and second moments (sum, sum2).
    const double var = compute_variance_from_moments(sum, sum2, valid);
    M.std_deg = std::sqrt(var);

    // Median is robust to outliers ( requires collecting all angles ).
    M.median_deg = compute_median(angles_deg);

    // Percent of pixels below common thresholds.
    M.pct_below_5deg = 100.0 * static_cast<double>(count5) / static_cast<double>(valid);
    M.pct_below_10deg = 100.0 * static_cast<double>(count10) / static_cast<double>(valid);
    M.pct_below_20deg = 100.0 * static_cast<double>(count20) / static_cast<double>(valid);

    M.min_deg = minA;
    M.max_deg = maxA;

    return M;
}