// cr_sweep_bounds.h
// Header-only reference implementation (C++11) of per-segment bounding for
// Centripedal Catmull–Rom (α=1/2) sphere sweeps using CR→Bezier conversion
// and the cubic Bezier convex-hull property. Inflate by r_max via scalar Bézier
// control points. Intended as a drop-in helper you can adapt to POV-Ray sources.
//
// Sources:
//  * Yuksel, Schaefer, Keyser: Catmull–Rom parameterization & Bézier control
//    point formulation (see http://www.cemyuksel.com/research/catmullrom_param/)
//  * Bezier convex-hull property (e.g., Stanford CS148 notes)
//
// License: Public domain / CC0

#ifndef CR_SWEEP_BOUNDS_H
#define CR_SWEEP_BOUNDS_H

#include <algorithm>
#include <array>
#include <cmath>
#include <cfloat>

namespace crbound {

struct Vec3 {
    double x, y, z;
    Vec3() : x(0), y(0), z(0) {}
    Vec3(double X,double Y,double Z):x(X),y(Y),z(Z){}
    Vec3 operator+(const Vec3& b) const { return Vec3(x+b.x,y+b.y,z+b.z); }
    Vec3 operator-(const Vec3& b) const { return Vec3(x-b.x,y-b.y,z-b.z); }
    Vec3 operator*(double s) const { return Vec3(x*s,y*s,z*s); }
};

inline double length(const Vec3& v) { return std::sqrt(v.x*v.x + v.y*v.y + v.z*v.z); }
inline Vec3 vmin(const Vec3& a, const Vec3& b) { return Vec3(std::min(a.x,b.x), std::min(a.y,b.y), std::min(a.z,b.z)); }
inline Vec3 vmax(const Vec3& a, const Vec3& b) { return Vec3(std::max(a.x,b.x), std::max(a.y,b.y), std::max(a.z,b.z)); }

struct AABB { Vec3 min, max; };

// Compute centripetal CR knots t0..t3 for four points
inline void centripetalKnots(const Vec3& P0, const Vec3& P1, const Vec3& P2, const Vec3& P3,
                             double alpha, std::array<double,4>& T)
{
    const double eps = 1e-12;
    T[0] = 0.0;
    double d01 = std::pow(std::max(length(P1-P0), eps), alpha);
    T[1] = T[0] + d01;
    double d12 = std::pow(std::max(length(P2-P1), eps), alpha);
    T[2] = T[1] + d12;
    double d23 = std::pow(std::max(length(P3-P2), eps), alpha);
    T[3] = T[2] + d23;
}

// Convert CR segment (P1->P2 with neighbors P0,P3) to cubic Bezier control points
inline void crSegmentToBezier(const Vec3& P0, const Vec3& P1, const Vec3& P2, const Vec3& P3,
                              const std::array<double,4>& T, std::array<Vec3,4>& B)
{
    const double eps = 1e-12;
    const double t0=T[0], t1=T[1], t2=T[2], t3=T[3];
    // endpoint derivatives w.r.t original parameter t
    Vec3 m1 = (P2 - P0) * ((t2 - t1) / std::max(t2 - t0, eps));
    Vec3 m2 = (P3 - P1) * ((t2 - t1) / std::max(t3 - t1, eps));
    // reparameterize to s in [0,1]
    double k = (t2 - t1);
    Vec3 m1s = m1 * k;
    Vec3 m2s = m2 * k;
    B[0] = P1;
    B[1] = P1 + (m1s * (1.0/3.0));
    B[2] = P2 - (m2s * (1.0/3.0));
    B[3] = P2;
}

// Scalar (radius) version: returns 4 scalar Bezier controls for radius
inline void crScalarToBezier(double R0,double R1,double R2,double R3,
                             const std::array<double,4>& T, std::array<double,4>& BR)
{
    const double eps = 1e-12;
    const double t0=T[0], t1=T[1], t2=T[2], t3=T[3];
    double m1 = (R2 - R0) * ((t2 - t1) / std::max(t2 - t0, eps));
    double m2 = (R3 - R1) * ((t2 - t1) / std::max(t3 - t1, eps));
    double k  = (t2 - t1);
    double m1s = m1 * k;
    double m2s = m2 * k;
    BR[0] = R1;
    BR[1] = R1 + (m1s / 3.0);
    BR[2] = R2 - (m2s / 3.0);
    BR[3] = R2;
}

// Per-segment AABB from Bezier control points, inflated by r_max (from scalar Bezier controls)
inline AABB segmentAABB(const std::array<Vec3,4>& B, const std::array<double,4>& BR)
{
    AABB aabb; aabb.min = aabb.max = B[0];
    for (int i=1;i<4;++i) { aabb.min = vmin(aabb.min, B[i]); aabb.max = vmax(aabb.max, B[i]); }
    double rmax = std::max(std::max(BR[0],BR[1]), std::max(BR[2],BR[3]));
    Vec3 e(rmax,rmax,rmax);
    aabb.min = aabb.min - e;
    aabb.max = aabb.max + e;
    return aabb;
}

// Union of two AABBs
inline AABB unite(const AABB& a, const AABB& b)
{
    AABB r; r.min = vmin(a.min,b.min); r.max = vmax(a.max,b.max); return r;
}

// Compute global AABB for a CR sphere-sweep from arrays of points/radii.
// Expects at least 4 points; each segment is (P_{i}..P_{i+1}) for i=1..n-3.
inline AABB sweepAABB(const Vec3* pts, const double* rad, int n, double alpha = 0.5)
{
    AABB global; global.min = Vec3(DBL_MAX,DBL_MAX,DBL_MAX); global.max = Vec3(-DBL_MAX,-DBL_MAX,-DBL_MAX);
    if (n < 4) return global;
    for (int i=1; i <= n-3; ++i)
    {
        const Vec3& P0 = pts[i-1];
        const Vec3& P1 = pts[i  ];
        const Vec3& P2 = pts[i+1];
        const Vec3& P3 = pts[i+2];
        double R0 = rad[i-1], R1 = rad[i], R2 = rad[i+1], R3 = rad[i+2];
        std::array<double,4> T; centripetalKnots(P0,P1,P2,P3, alpha, T);
        std::array<Vec3,4>   B; crSegmentToBezier(P0,P1,P2,P3, T, B);
        std::array<double,4> BR; crScalarToBezier(R0,R1,R2,R3, T, BR);
        AABB seg = segmentAABB(B, BR);
        if (global.min.x == DBL_MAX) global = seg; else global = unite(global, seg);
    }
    return global;
}

} // namespace crbound

#endif // CR_SWEEP_BOUNDS_H