Renderive/render_2D/plottable/Curve_Sampling.h

#pragma once

#include "../axis/Axis.h"
#include "Heatmap_Utils.h"

#include <algorithm>
#include <cmath>
#include <span>
#include <vector>

namespace renderive::detail {
namespace curve_sampling {

struct Sample {
    double coordinate{};
    double value{};
};

inline double power_domain_lerp(double first, double second, double ratio) {
    const double first_power = std::pow(10.0, std::clamp(first, -3'000.0, 3'000.0) / 10.0);
    const double second_power = std::pow(10.0, std::clamp(second, -3'000.0, 3'000.0) / 10.0);
    const double power = first_power + (second_power - first_power) * ratio;
    return 10.0 * std::log10(std::max(power, 1e-300));
}

inline double cubic_value(double previous, double first, double second, double next,
                          double ratio) {
    const double ratio2 = ratio * ratio;
    const double ratio3 = ratio2 * ratio;
    return 0.5 * ((2.0 * first) + (-previous + second) * ratio +
                  (2.0 * previous - 5.0 * first + 4.0 * second - next) * ratio2 +
                  (-previous + 3.0 * first - 3.0 * second + next) * ratio3);
}

inline std::vector<Sample> interpolate(std::span<const double> values, Range domain,
                                       Line_Interpolation_Mode mode) {
    std::vector<Sample> result;
    if (values.empty())
        return result;
    if (values.size() == 1) {
        result.push_back({domain.origin, values.front()});
        return result;
    }

    constexpr int smooth_subdivisions = 4;
    const std::size_t multiplier =
        mode == Line_Interpolation_Mode::Linear_Power_Domain ||
                mode == Line_Interpolation_Mode::Cubic_Value
            ? smooth_subdivisions
            : mode == Line_Interpolation_Mode::Linear_Value ? 1 : 3;
    result.reserve(1 + (values.size() - 1) * multiplier);

    const double denominator = static_cast<double>(values.size() - 1);
    const auto coordinate_at = [domain, denominator](std::size_t index) {
        return domain.origin + domain.length() * static_cast<double>(index) / denominator;
    };
    result.push_back({coordinate_at(0), values.front()});
    for (std::size_t index = 0; index + 1 < values.size(); ++index) {
        const double first_coordinate = coordinate_at(index);
        const double second_coordinate = coordinate_at(index + 1);
        const double first_value = values[index];
        const double second_value = values[index + 1];
        switch (mode) {
        case Line_Interpolation_Mode::Nearest_Sample: {
            const double middle = (first_coordinate + second_coordinate) * 0.5;
            result.push_back({middle, first_value});
            result.push_back({middle, second_value});
            result.push_back({second_coordinate, second_value});
            break;
        }
        case Line_Interpolation_Mode::Linear_Value:
            result.push_back({second_coordinate, second_value});
            break;
        case Line_Interpolation_Mode::Linear_Power_Domain:
            for (int part = 1; part <= smooth_subdivisions; ++part) {
                const double ratio = static_cast<double>(part) / smooth_subdivisions;
                result.push_back({first_coordinate +
                                      (second_coordinate - first_coordinate) * ratio,
                                  power_domain_lerp(first_value, second_value, ratio)});
            }
            break;
        case Line_Interpolation_Mode::Step_Left:
            result.push_back({second_coordinate, first_value});
            result.push_back({second_coordinate, second_value});
            break;
        case Line_Interpolation_Mode::Step_Right:
            result.push_back({first_coordinate, second_value});
            result.push_back({second_coordinate, second_value});
            break;
        case Line_Interpolation_Mode::Cubic_Value: {
            const double previous = values[index == 0 ? 0 : index - 1];
            const double next = values[std::min(index + 2, values.size() - 1)];
            for (int part = 1; part <= smooth_subdivisions; ++part) {
                const double ratio = static_cast<double>(part) / smooth_subdivisions;
                result.push_back({first_coordinate +
                                      (second_coordinate - first_coordinate) * ratio,
                                  cubic_value(previous, first_value, second_value, next, ratio)});
            }
            break;
        }
        }
    }
    return result;
}

inline std::vector<Sample> visible_samples(std::vector<Sample> samples,
                                           Range visible_range) {
    if (samples.size() < 2)
        return visible_range.contains(samples.empty() ? 0.0 : samples.front().coordinate)
                   ? std::move(samples)
                   : std::vector<Sample>{};
    const auto [visible_low, visible_high] =
        std::minmax(visible_range.origin, visible_range.target);
    std::size_t first = samples.size();
    std::size_t last{};
    for (std::size_t index = 0; index + 1 < samples.size(); ++index) {
        const auto [segment_low, segment_high] =
            std::minmax(samples[index].coordinate, samples[index + 1].coordinate);
        if (segment_high < visible_low || segment_low > visible_high)
            continue;
        first = std::min(first, index);
        last = std::max(last, index + 1);
    }
    if (first == samples.size())
        return {};
    return {samples.begin() + static_cast<std::ptrdiff_t>(first),
            samples.begin() + static_cast<std::ptrdiff_t>(last + 1)};
}

} // namespace curve_sampling

inline std::vector<PointF> curve_points(std::span<const double> values,
                                        Range domain,
                                        const Axis_Transform& x_axis,
                                        const Axis_Transform& y_axis,
                                        bool visible_only,
                                        Line_Interpolation_Mode mode) {
    auto samples = curve_sampling::interpolate(values, domain, mode);
    if (visible_only)
        samples = curve_sampling::visible_samples(std::move(samples), x_axis.coordinate_range);
    std::vector<PointF> points;
    points.reserve(samples.size());
    for (const auto& sample : samples) {
        if (std::isfinite(sample.coordinate) && std::isfinite(sample.value))
            points.push_back(mapped_point(x_axis, sample.coordinate, y_axis, sample.value));
    }
    return points;
}

} // namespace renderive::detail