/*
    n3xd_ocp.sample — a face's UV grid, without the per-sample round trip.

    cad/operations/surface_pattern.py walks an (res+1)² UV grid in nested
    Python loops: one BRepAdaptor_Surface::D1 plus roughly ten accessor calls
    per sample, and the cross product in Python.  The kernel can fill the whole
    grid in one call with the GIL released.

    The output reproduces what the Python loop produces:

      * the grid is BRepTools::UVBounds, sampled the way numpy's linspace
        samples it — start + i*step with the endpoint forced exactly, so the
        sample locations are bit-identical and the fitted surface does not
        move.  The bounds come back with the arrays for the same reason: the
        caller derives its normalised u/v from these numbers rather than
        recomputing them.
      * normals are du × dv normalised, negated for a reversed face.
      * a sample whose cross product is degenerate (a pole or apex) gets the
        zero vector, not a unit vector in an arbitrary direction.  Callers
        test it — `_inward_thickness` skips a sample on `nrm.any()`.
*/

#include "../common/occt_module.h"

#include <nanobind/ndarray.h>

#include <BRepAdaptor_Surface.hxx>
#include <BRepTools.hxx>
#include <TopAbs_Orientation.hxx>
#include <TopoDS_Face.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>

#include <cmath>
#include <stdexcept>

namespace {

template <typename T> nb::capsule owner_of(T *data) {
    return nb::capsule(data, [](void *p) noexcept { delete[] (T *) p; });
}

/// numpy's linspace endpoint handling: the last sample is the bound itself,
/// not start + (n-1)*step, which can miss it by an ulp.
inline double linspace_at(double lo, double hi, int i, int n) {
    if (n <= 1)
        return lo;
    return i == n - 1 ? hi : lo + (double) i * (hi - lo) / (double) (n - 1);
}

constexpr double MIN_NORMAL = 1e-12;

}  // namespace

void register_ext_sample() {
    nb::module_ m = ocp_named_module("n3xd_ocp.sample");

    m.def(
        "face_grid",
        [](const TopoDS_Face &face, int n) {
            if (n < 1)
                throw std::invalid_argument("face_grid: n must be >= 1");

            double umin = 0.0, umax = 0.0, vmin = 0.0, vmax = 0.0;
            BRepTools::UVBounds(face, umin, umax, vmin, vmax);
            const double sign =
                face.Orientation() == TopAbs_REVERSED ? -1.0 : 1.0;

            const size_t count = (size_t) n * (size_t) n;
            double *points = new double[3 * count];
            double *normals = new double[3 * count];

            {
                nb::gil_scoped_release nogil;
                BRepAdaptor_Surface adaptor(face);
                gp_Pnt point;
                gp_Vec du, dv;
                for (int i = 0; i < n; ++i) {
                    const double u = linspace_at(umin, umax, i, n);
                    for (int j = 0; j < n; ++j) {
                        const double v = linspace_at(vmin, vmax, j, n);
                        adaptor.D1(u, v, point, du, dv);

                        const size_t at = 3 * ((size_t) i * (size_t) n + j);
                        points[at + 0] = point.X();
                        points[at + 1] = point.Y();
                        points[at + 2] = point.Z();

                        const double nx = du.Y() * dv.Z() - du.Z() * dv.Y();
                        const double ny = du.Z() * dv.X() - du.X() * dv.Z();
                        const double nz = du.X() * dv.Y() - du.Y() * dv.X();
                        const double mag = std::sqrt(nx * nx + ny * ny + nz * nz);
                        if (mag > MIN_NORMAL) {
                            normals[at + 0] = sign * nx / mag;
                            normals[at + 1] = sign * ny / mag;
                            normals[at + 2] = sign * nz / mag;
                        } else {
                            normals[at + 0] = 0.0;
                            normals[at + 1] = 0.0;
                            normals[at + 2] = 0.0;
                        }
                    }
                }
            }

            size_t shape[3] = {(size_t) n, (size_t) n, 3};
            return nb::make_tuple(
                nb::ndarray<nb::numpy, double, nb::ndim<3>>(points, 3, shape,
                                                            owner_of(points)),
                nb::ndarray<nb::numpy, double, nb::ndim<3>>(normals, 3, shape,
                                                            owner_of(normals)),
                nb::make_tuple(umin, umax, vmin, vmax));
        },
        "face"_a, "n"_a,
        R"doc(Sample *face* on an n x n grid over its UV bounds.

Returns ``(points[n, n, 3], normals[n, n, 3], (umin, umax, vmin, vmax))`` as
float64 arrays, u along axis 0.  Normals are outward (negated for a reversed
face) and unit length, or the zero vector where the surface derivatives are
degenerate.  Sample locations match ``numpy.linspace`` over the returned
bounds exactly.)doc");
}