reasampler/src/core/instrument/engine/velocity_curve.h

// velocity_curve.h — THE monotone spline, shared by every consumer: the three velocity
// transfer curves (amp gain, pitch offset, filter cutoff offset), evaluated once per note-on,
// and the spline EGs, evaluated per voice per sample through SplineCursor. Editor
// hit-test/inverse-map take an explicit pixel Box rather than a Rect: this module sits below
// sampler_core in the link graph and must not gain a dependency on editor-layout types.

#pragma once

#include <cstddef>
#include <cstdint>
#include <vector>

namespace reasampler::instrument::engine {

// The curve's canonical X span. For the three velocity consumers it IS the MIDI velocity
// domain; a spline EG maps normalized sample time onto the same span, which is what lets one
// implementation serve both without a second X domain to keep in sync.
inline constexpr double kCurveXMin = 0.0;
inline constexpr double kCurveXMax = 127.0;
inline constexpr double kVelMin = kCurveXMin;  // the velocity consumers' spelling of the span
inline constexpr double kVelMax = kCurveXMax;
inline constexpr double kCurveYMax = 1.0;

// Point-count ceiling. A MUSICAL bound, not a performance one: long rhythmic phrases need the
// resolution, and at roughly two points per articulation event 128 is about four bars of 16ths.
// Segment lookup is logarithmic and the editor's node separation is the real density limit, so
// there is nothing to buy by lowering it. DO NOT LOWER.
inline constexpr std::size_t kMaxCurvePoints = 128;

// The curve's Y range. UNIPOLAR [0,1] is a GAIN — the amp's domain, where the do-nothing
// curve is flat at 1. BIPOLAR [-1,1] is a SIGNED modulation shape — the pitch and filter
// domains, where the do-nothing curve is flat at 0 and the sign picks the direction. A
// bipolar curve does not preclude a depth control beside it: the filter has one, and the two
// compose multiplicatively (play_params.h).
enum class CurveDomain { Unipolar, Bipolar };

constexpr double curveYMin(CurveDomain d) { return d == CurveDomain::Bipolar ? -1.0 : 0.0; }

// The value that changes nothing in each domain — unity gain, or zero modulation. THE one home
// for that value: eval()'s own empty-curve fallback reads it directly, and flat()/zero() (what
// fromPoints' sub-2-point fallback constructs) are built from it too, so a corrupt blob always
// loses the shaping rather than inventing one, however the fallback is reached.
constexpr double curveNeutral(CurveDomain d) { return d == CurveDomain::Bipolar ? 0.0 : 1.0; }

// A raw-constructed point is NOT auto-clamped (the mutators own that invariant) — build curves
// through the named constructors / addPoint rather than pushing raw points.
struct VelocityPoint {
    double velocity = 0.0;   // X, over the canonical span
    double value    = 0.0;   // Y, in the owning curve's domain
    // A HARD point does no smoothing on either side: it terminates the monotone sub-curve, so
    // the two adjacent segments meet at their own natural angle instead of a shared derivative.
    // Points are smooth by default; see segmentTangents for the mechanism.
    bool   hard     = false;
};

// Fritsch-Carlson monotone-cubic tangent: a sign change (or flat) neighbour is a local extremum,
// so the tangent pins to 0 to avoid overshoot; otherwise the weighted-harmonic-mean tangent,
// which for collinear knots (dPrev==dNext) reduces exactly to the shared secant — this is what
// makes the spline reproduce a straight line for linear()-style input.
inline double fritschCarlsonTangent(double dPrev, double dNext, double spanPrev, double spanNext) {
    if (dPrev * dNext <= 0.0) return 0.0;
    const double w1 = 2.0 * spanNext + spanPrev;
    const double w2 = spanNext + 2.0 * spanPrev;
    return (w1 + w2) / (w1 / dPrev + w2 / dNext);
}

struct SegmentTangents {
    double mA = 0.0;
    double mB = 0.0;
};

// The Hermite tangents for segment [i, i+1] of an X-ordered point array, where `d` is that
// segment's secant slope and `span` its X width (> 0).
//
// A HARD point is treated exactly as the array's own end is: the tangent there is the segment's
// own secant, so smoothing stops at it. That single rule is the whole hard-point enhancement —
// the contour becomes one or more monotone splines joined at their natural angles, and each
// sub-curve keeps Fritsch-Carlson's no-overshoot guarantee because m == d satisfies its bound.
inline SegmentTangents segmentTangents(const VelocityPoint* p, std::size_t n, std::size_t i,
                                       double d, double span) {
    SegmentTangents t{d, d};
    if (i > 0 && !p[i].hard) {
        const double spanPrev = p[i].velocity - p[i - 1].velocity;
        t.mA = (spanPrev > 0.0)
                   ? fritschCarlsonTangent((p[i].value - p[i - 1].value) / spanPrev, d, spanPrev,
                                           span)
                   : 0.0;
    }
    if (i + 2 < n && !p[i + 1].hard) {
        const double spanNext = p[i + 2].velocity - p[i + 1].velocity;
        t.mB = (spanNext > 0.0)
                   ? fritschCarlsonTangent(d, (p[i + 2].value - p[i + 1].value) / spanNext, span,
                                           spanNext)
                   : 0.0;
    }
    return t;
}

// The cubic Hermite basis evaluated at t in [0,1] across a segment of width `span`.
inline double hermiteAt(double y0, double y1, double span, double mA, double mB, double t) {
    const double t2 = t * t;
    const double t3 = t2 * t;
    return (2.0 * t3 - 3.0 * t2 + 1.0) * y0 + (t3 - 2.0 * t2 + t) * span * mA +
           (-2.0 * t3 + 3.0 * t2) * y1 + (t3 - t2) * span * mB;
}

// Pick radius (px) around a node's drawn point for the editor hit-test.
inline constexpr int kCurveNodeGrabRadius = 6;

// An X-ordered list of control points spanning the canonical X span, evaluated as ONE OR MORE
// monotone cubic Hermite splines (Fritsch-Carlson slope limiting) joined at the hard points — a
// genuine curve, not a polyline, that provably never overshoots any segment's value range. The
// guarantee is PER SEGMENT, so a contour is free to rise and fall. For collinear knots the
// tangents reduce to the secant slope, so the spline reproduces linear()'s straight line to
// within ~1e-15. The two endpoints are load-bearing: they keep eval total over the domain and
// are never deletable.
class VelocityCurve {
public:
    // flat() (endpoints (0,1)/(127,1), every velocity -> unity) is the unipolar default — see
    // velocity_curve in the directory CLAUDE.md for why this isn't bit-identical to the
    // pre-existing linear() response.
    static VelocityCurve flat();
    static VelocityCurve linear();
    // The bipolar default: flat at 0, so velocity modulates nothing until a curve is drawn.
    static VelocityCurve zero();
    // y = 1 - x: the smooth downward slope a freshly created spline EG opens on. Two collinear
    // knots, so it is straight — and straight is smooth. NOT a change to any velocity curve's
    // own default.
    static VelocityCurve rampDown();

    // Rebuilds from a deserialized point list, repairing the invariant defensively: box-clamps
    // each point into `domain`, stable-sorts by velocity, forces both endpoints present
    // (synthesized if missing), falls back to the domain's neutral curve if fewer than 2 usable
    // points remain. A corrupt/truncated blob yields a well-formed curve, never an
    // invariant-violating one.
    static VelocityCurve fromPoints(std::vector<VelocityPoint> pts, CurveDomain domain);

    CurveDomain domain() const { return domain_; }
    const std::vector<VelocityPoint>& points() const { return points_; }
    std::size_t size() const { return points_.size(); }

    // Degenerate cases (shouldn't occur post-construction): empty curve returns the domain's
    // neutral; a one-point curve returns that point's value.
    double eval(double velocity) const;

    // Inserted at a velocity duplicating an existing point lands immediately after it, so a
    // subsequent move can separate them. Returns the inserted index, or -1 when the curve is
    // already at kMaxCurvePoints — a refusal leaves the contour bit-identical.
    int addPoint(double velocity, double value);

    // Flips a point between hard and smooth. Out-of-range index is a no-op returning false.
    // Permitted on the endpoints, where it changes nothing evaluable: an endpoint's outward
    // tangent is already its own secant, which is what hard means.
    bool toggleHard(std::size_t index);
    bool setHard(std::size_t index, bool hard);

    // Box-clamped and X-clamped between immediate neighbours (monotonic-X grammar). The two
    // endpoints are pinned in X (only their value moves); out-of-range index is a no-op.
    VelocityPoint movePoint(std::size_t index, double velocity, double value);

    // Endpoints (index 0 and last) are not deletable; that or an out-of-range index is a no-op
    // returning false.
    bool deletePoint(std::size_t index);

    // The drawn box, in pixels: X = velocity across the width, Y = value UP the height (the
    // domain's max at top). Passed explicitly rather than a Rect — see header preamble.
    struct Box {
        int left = 0;
        int top = 0;
        int width = 0;
        int height = 0;
    };

    // Index of the first point within the pick radius on both axes, or -1 for a miss. First-match
    // in point order for determinism.
    int pointAtPixel(const Box& box, int x, int y) const;

    // The one point->pixel mapping, exposed so drawing and hit-testing can never drift apart.
    struct CurvePixel {
        int x = 0;
        int y = 0;
    };
    CurvePixel pixelFromPoint(const Box& box, const VelocityPoint& p) const;

    // Exact inverse of pixelFromPoint (within the one-pixel quantum) — where an empty-space
    // click lands as a new point. Degenerate box: zero-width reads velocity 0; height <= 1
    // reads the domain's max (the top row is what a collapsed box draws).
    VelocityPoint pointFromPixel(const Box& box, int x, int y) const;

    // `grabCurve` is the curve as of mouse-down (shell snapshots it so the delta is absolute).
    // Maps the pixel delta to velocity/value over the box, then applies movePoint's clamp. Zero
    // width/height box or out-of-range index returns grabCurve unchanged.
    static VelocityCurve resolvePointDrag(const VelocityCurve& grabCurve, std::size_t index,
                                          const Box& box, int dxPixels, int dyPixels);

    bool equals(const VelocityCurve& other, double eps = 1e-9) const;

private:
    // Private: an implicit-default curve is empty (no endpoints) and Unipolar, so a stray
    // default-construction wouldn't fail loudly — it would eval() to unity gain everywhere,
    // or a full +/-1 (a full-scale transpose / wide-open filter) if ever read as bipolar. Build
    // through flat()/linear()/zero()/fromPoints(), all of which establish the endpoint invariant.
    VelocityCurve() = default;

    // Always X-ordered with an endpoint at 0 and 127; constructors + deserialize establish the
    // invariant, mutators preserve it.
    std::vector<VelocityPoint> points_;
    CurveDomain domain_ = CurveDomain::Unipolar;
};

// The RT read head over a contour: an indexed segment search plus one Hermite evaluation, with
// the segment and its two tangents cached across samples so a monotone read costs one compare.
// Header-inline, branch-only, NO allocation and NO virtual dispatch — it runs per voice per
// sample. A jump (a loop wrap, a fresh note) falls back to a binary search, <= 7 steps at the
// 128-point ceiling.
//
// Holds a RAW POINTER into the bound curve's point array: the caller guarantees the curve
// outlives the cursor. The voice binds against its SampleData, which has exactly that lifetime.
class SplineCursor {
public:
    // Binds `c` if it has an evaluable segment; a shorter curve leaves the cursor inactive so
    // the caller's `if (active())` skips the whole spline path.
    void bind(const VelocityCurve& c) {
        const std::vector<VelocityPoint>& pts = c.points();
        if (pts.size() < 2) { clear(); return; }
        pts_ = pts.data();
        n_ = pts.size();
        select(0);
    }
    void clear() { pts_ = nullptr; n_ = 0; }
    bool active() const { return n_ >= 2; }

    // `phase` is normalized position over the contour's whole span, [0,1]; out-of-range clamps
    // to the terminal values (a note past its span holds the contour's last level).
    double eval(double phase) {
        const double x = (phase <= 0.0)   ? kCurveXMin
                         : (phase >= 1.0) ? kCurveXMax
                                          : kCurveXMin + phase * (kCurveXMax - kCurveXMin);
        if (x <= x0_ && seg_ == 0) return y0_;
        if (x >= x1_ && seg_ + 2 == n_) return y1_;
        if (x < x0_ || x > x1_) locate(x);
        if (span_ <= 0.0) return y1_;  // coincident-X knots: a step, no interior to blend
        return hermiteAt(y0_, y1_, span_, mA_, mB_, (x - x0_) / span_);
    }

private:
    // The common case is the next segment (a monotone read walking forward); anything else is a
    // binary search over the X-ordered array.
    void locate(double x) {
        if (x > x1_ && seg_ + 2 < n_ && x <= pts_[seg_ + 2].velocity) { select(seg_ + 1); return; }
        std::size_t lo = 0, hi = n_ - 2;
        while (lo < hi) {
            const std::size_t mid = lo + (hi - lo + 1) / 2;
            if (pts_[mid].velocity <= x) lo = mid; else hi = mid - 1;
        }
        select(lo);
    }

    void select(std::size_t i) {
        seg_ = i;
        x0_ = pts_[i].velocity;
        x1_ = pts_[i + 1].velocity;
        y0_ = pts_[i].value;
        y1_ = pts_[i + 1].value;
        span_ = x1_ - x0_;
        const SegmentTangents t =
            segmentTangents(pts_, n_, i, span_ > 0.0 ? (y1_ - y0_) / span_ : 0.0, span_);
        mA_ = t.mA;
        mB_ = t.mB;
    }

    const VelocityPoint* pts_ = nullptr;
    std::size_t n_ = 0;
    std::size_t seg_ = 0;
    double x0_ = 0.0, x1_ = 0.0, y0_ = 0.0, y1_ = 0.0, span_ = 0.0, mA_ = 0.0, mB_ = 0.0;
};

}  // namespace reasampler::instrument::engine