#pragma once
#include "vde/core/point.h"
#include <functional>
#include <tuple>

namespace vde::sdf {

using core::Point3D;
using core::Vector3D;

/// Compute numerical gradient of a scalar field f at point p
/// Uses central finite differences with step size h
template <typename Fn>
[[nodiscard]] Vector3D gradient(Fn&& f, const Point3D& p, double h = 1e-6) {
    double fx = f(p);
    const double eps_x = h * (1.0 + std::abs(p.x()));
    const double eps_y = h * (1.0 + std::abs(p.y()));
    const double eps_z = h * (1.0 + std::abs(p.z()));

    double dfdx = (f(Point3D(p.x() + eps_x, p.y(), p.z())) -
                   f(Point3D(p.x() - eps_x, p.y(), p.z()))) / (2.0 * eps_x);
    double dfdy = (f(Point3D(p.x(), p.y() + eps_y, p.z())) -
                   f(Point3D(p.x(), p.y() - eps_y, p.z()))) / (2.0 * eps_y);
    double dfdz = (f(Point3D(p.x(), p.y(), p.z() + eps_z)) -
                   f(Point3D(p.x(), p.y(), p.z() - eps_z))) / (2.0 * eps_z);

    return Vector3D(dfdx, dfdy, dfdz);
}

/// Evaluate scalar field and its gradient at point p (joint pass)
/// Returns (value, gradient)
template <typename Fn>
[[nodiscard]] std::tuple<double, Vector3D> evaluate_with_gradient(
        Fn&& f, const Point3D& p, double h = 1e-6) {
    return {f(p), gradient(std::forward<Fn>(f), p, h)};
}

} // namespace vde::sdf