// STREAMLINES :: https://prideout.net/blog/par_streamlines/
// Simple C library for triangulating wide lines, curves, and streamlines.
//
// Usage example:
//
//   #define PAR_STREAMLINES_IMPLEMENTATION
//   #include "par_streamlines.h"
//
//   parsl_context* ctx = parsl_create_context({ .thickness = 3 });
//   parsl_position vertices[] = { {0, 0}, {2, 1}, {4, 0} };
//   uint16_t spine_lengths[] = { 3 };
//   parsl_mesh* mesh = parsl_mesh_from_lines(ctx, {
//       .num_vertices = sizeof(vertices) / sizeof(parsl_position),
//       .num_spines = sizeof(spine_lengths) / sizeof(uint16_t),
//       .vertices = vertices,
//       .spine_lengths = spine_lengths
//   });
//   ...
//   parsl_destroy_context(ctx);
//
// Distributed under the MIT License, see bottom of file.

#ifndef PAR_STREAMLINES_H
#define PAR_STREAMLINES_H

#ifdef __cplusplus
extern "C" {
#endif

#include <stdint.h>
#include <stdbool.h>

// Configures how the library assigns UV coordinates.
typedef enum {
    PAR_U_MODE_NORMALIZED_DISTANCE, // this is the default
    PAR_U_MODE_DISTANCE,            // non-normalized distance along the curve
    PAR_U_MODE_SEGMENT_INDEX,       // starts at zero for each curve, counts up
    PAR_U_MODE_SEGMENT_FRACTION,    // 0.0, 1.0 / COUNT, 2.0 / COUNT, etc...
} parsl_u_mode;

// Layout for generated vertex attributes.
typedef struct {
    float u_along_curve;   // longitudinal coordinate (see parsl_u_mode)
    float v_across_curve;  // either + or - depending on the side
    float spine_to_edge_x; // normalized vector from spine to edge
    float spine_to_edge_y; // normalized vector from spine to edge
} parsl_annotation;

// Simple two-tuple math type used for mesh and spine vertices.
typedef struct {
    float x;
    float y;
} parsl_position;

// Triangle mesh generated by the library. The vertex data is owned by
// streamlines context and becomes invalid on any subsequent call to the API.
// The annotations, spine_lengths, and random_offsets fields are null unless
// their corresponding flags have been set in parsl_config.
typedef struct {
    uint32_t num_vertices;
    uint32_t num_triangles;
    uint32_t* triangle_indices;
    parsl_position* positions;
    parsl_annotation* annotations;
    float* spine_lengths;
    float* random_offsets;
} parsl_mesh;

// Viewport for streamline seed placement.
typedef struct {
    float left;
    float top;
    float right;
    float bottom;
} parsl_viewport;

#define PARSL_FLAG_WIREFRAME      (1 << 0) // enables 4 indices per triangle
#define PARSL_FLAG_ANNOTATIONS    (1 << 1) // populates mesh.annotations
#define PARSL_FLAG_SPINE_LENGTHS  (1 << 2) // populates mesh.lengths
#define PARSL_FLAG_RANDOM_OFFSETS (1 << 3) // populates mesh.random_offsets
#define PARSL_FLAG_CURVE_GUIDES   (1 << 4) // draws control points

// Immutable configuration for a streamlines context.
typedef struct {
    float thickness;
    uint32_t flags;
    parsl_u_mode u_mode;
    float curves_max_flatness;
    float streamlines_seed_spacing;
    parsl_viewport streamlines_seed_viewport;
    float miter_limit;
} parsl_config;

// Client-owned list of line strips that will be tessellated.
typedef struct {
    uint32_t num_vertices;
    uint16_t num_spines;
    parsl_position* vertices;
    uint16_t* spine_lengths;
    bool closed;
} parsl_spine_list;

// Opaque handle to a streamlines context and its memory arena.
typedef struct parsl_context_s parsl_context;

// Client function that moves a streamline particle by a single time step.
typedef void (*parsl_advection_callback)(parsl_position* point, void* userdata);

parsl_context* parsl_create_context(parsl_config config);

void parsl_destroy_context(parsl_context* ctx);

// Low-level function that simply generates two triangles for each line segment.
parsl_mesh* parsl_mesh_from_lines(parsl_context* ctx, parsl_spine_list spines);

// High-level function that can be used to visualize a vector field.
parsl_mesh* parsl_mesh_from_streamlines(parsl_context* context,
    parsl_advection_callback advect, uint32_t first_tick, uint32_t num_ticks,
    void* userdata);

// High-level function that tessellates a series of curves into triangles,
// where each spine is a series of chained cubic Bézier curves.
//
// The first curve of each spine is defined by an endpoint, followed by two
// control points, followed by an endpoint. Every subsequent curve in the spine
// is defined by a single control point followed by an endpoint. Only one
// control point is required because the first control point is computed via
// reflection over the endpoint.
//
// The number of vertices in each spine should be 4+(n-1)*2 where n is the
// number of piecewise curves.
//
// Each spine is equivalent to an SVG path that looks like M C S S S.
parsl_mesh* parsl_mesh_from_curves_cubic(parsl_context* context,
    parsl_spine_list spines);

// High-level function that tessellates a series of curves into triangles,
// where each spine is a series of chained quadratic Bézier curves.
//
// The first curve of each spine is defined by an endpoint, followed by one
// control point, followed by an endpoint. Every subsequent curve in the spine
// is defined by a single control point followed by an endpoint.
//
// The number of vertices in each spine should be 3+(n-1)*2 where n is the
// number of piecewise curves.
//
// Each spine is equivalent to an SVG path that looks like M Q M Q M Q.
parsl_mesh* parsl_mesh_from_curves_quadratic(parsl_context* context,
    parsl_spine_list spines);

#ifdef __cplusplus
}
#endif

// -----------------------------------------------------------------------------
// END PUBLIC API
// -----------------------------------------------------------------------------
#ifdef PAR_STREAMLINES_IMPLEMENTATION

#include <assert.h>
#include <limits.h>
#include <math.h>
#include <memory.h>
#include <stdlib.h>

static float parsl__dot(parsl_position a, parsl_position b) {
    return a.x * b.x + a.y * b.y;
}

static parsl_position parsl__sub(parsl_position a, parsl_position b) {
    return (parsl_position) { a.x - b.x, a.y - b.y };
}

static parsl_position parsl__add(parsl_position a, parsl_position b) {
    return (parsl_position) { a.x + b.x, a.y + b.y };
}

static parsl_position parsl_mul(parsl_position v, float s) {
    return (parsl_position) { v.x * s, v.y * s };
}

#define PARSL_MAX_RECURSION 16

#ifndef PAR_PI
#define PAR_PI (3.14159265359)
#define PAR_MIN(a, b) (a > b ? b : a)
#define PAR_MAX(a, b) (a > b ? a : b)
#define PAR_CLAMP(v, lo, hi) PAR_MAX(lo, PAR_MIN(hi, v))
#define PAR_SWAP(T, A, B) { T tmp = B; B = A; A = tmp; }
#define PAR_SQR(a) ((a) * (a))
#endif

#ifndef PAR_MALLOC
#define PAR_MALLOC(T, N) ((T*) malloc(N * sizeof(T)))
#define PAR_CALLOC(T, N) ((T*) calloc(N * sizeof(T), 1))
#define PAR_REALLOC(T, BUF, N) ((T*) realloc(BUF, sizeof(T) * (N)))
#define PAR_FREE(BUF) free(BUF)
#endif

#ifndef PAR_ARRAY
#define PAR_ARRAY
#define pa_free(a) ((a) ? PAR_FREE(pa___raw(a)), 0 : 0)
#define pa_push(a, v) (pa___maybegrow(a, (int) 1), (a)[pa___n(a)++] = (v))
#define pa_count(a) ((a) ? pa___n(a) : 0)
#define pa_add(a, n) (pa___maybegrow(a, (int) n), pa___n(a) += (n))
#define pa_last(a) ((a)[pa___n(a) - 1])
#define pa_end(a) (a + pa_count(a))
#define pa_clear(arr) if (arr) pa___n(arr) = 0

#define pa___raw(a) ((int*) (a) -2)
#define pa___m(a) pa___raw(a)[0]
#define pa___n(a) pa___raw(a)[1]
#define pa___needgrow(a, n) ((a) == 0 || pa___n(a) + ((int) n) >= pa___m(a))
#define pa___maybegrow(a, n) (pa___needgrow(a, (n)) ? pa___grow(a, n) : 0)
#define pa___grow(a, n) (*((void**)& (a)) = pa___growf((void*) (a), (n), \
    sizeof(*(a))))

// ptr[-2] is capacity, ptr[-1] is size.
static void* pa___growf(void* arr, int increment, int itemsize)
{
    int dbl_cur = arr ? 2 * pa___m(arr) : 0;
    int min_needed = pa_count(arr) + increment;
    int m = dbl_cur > min_needed ? dbl_cur : min_needed;
    int* p = (int *) PAR_REALLOC(uint8_t, arr ? pa___raw(arr) : 0,
        itemsize * m + sizeof(int) * 2);
    if (p) {
        if (!arr) {
            p[1] = 0;
        }
        p[0] = m;
        return p + 2;
    }
    return (void*) (2 * sizeof(int));
}

#endif

struct parsl_context_s {
    parsl_config config;
    parsl_mesh result;
    parsl_position* streamline_seeds;
    parsl_position* streamline_points;
    parsl_spine_list streamline_spines;
    parsl_spine_list curve_spines;
    uint16_t guideline_start;
};

parsl_context* parsl_create_context(parsl_config config)
{
    parsl_context* context = PAR_CALLOC(parsl_context, 1);
    context->config = config;
    return context;
}

void parsl_destroy_context(parsl_context* context)
{
    pa_free(context->result.triangle_indices);
    pa_free(context->result.spine_lengths);
    pa_free(context->result.annotations);
    pa_free(context->result.positions);
    pa_free(context->result.random_offsets);
    pa_free(context->streamline_seeds);
    pa_free(context->streamline_points);
    pa_free(context->streamline_spines.spine_lengths);
    pa_free(context->streamline_spines.vertices);
    pa_free(context->curve_spines.spine_lengths);
    pa_free(context->curve_spines.vertices);
    PAR_FREE(context);
}

parsl_mesh* parsl_mesh_from_lines(parsl_context* context,
    parsl_spine_list spines)
{
    typedef parsl_position Position;
    typedef parsl_annotation Annotation;

    parsl_mesh* mesh = &context->result;
    const bool closed = spines.closed;
    const bool wireframe = context->config.flags & PARSL_FLAG_WIREFRAME;
    const bool has_annotations = context->config.flags & PARSL_FLAG_ANNOTATIONS;
    const bool has_lengths = context->config.flags & PARSL_FLAG_SPINE_LENGTHS;
    const float miter_limit = context->config.miter_limit ?
        context->config.miter_limit : (context->config.thickness * 2);
    const float miter_acos_max = +1.0;
    const float miter_acos_min = -1.0;
    const uint32_t ind_per_tri = wireframe ? 4 : 3;

    mesh->num_vertices = 0;
    mesh->num_triangles = 0;

    for (uint32_t spine = 0; spine < spines.num_spines; spine++) {
        assert(spines.spine_lengths[spine] > 1);
        mesh->num_vertices += 2 * spines.spine_lengths[spine];
        mesh->num_triangles += 2 * (spines.spine_lengths[spine] - 1);
        if (closed) {
            mesh->num_vertices += 2;
            mesh->num_triangles += 2;
        }
    }

    pa_clear(mesh->spine_lengths);
    pa_clear(mesh->annotations);
    pa_clear(mesh->positions);
    pa_clear(mesh->triangle_indices);

    if (has_lengths) {
        pa_add(mesh->spine_lengths, mesh->num_vertices);
    }
    if (has_annotations) {
        pa_add(mesh->annotations, mesh->num_vertices);
    }

    pa_add(mesh->positions, mesh->num_vertices);
    pa_add(mesh->triangle_indices, ind_per_tri * mesh->num_triangles);

    float* dst_lengths = mesh->spine_lengths;
    Annotation* dst_annotations = mesh->annotations;
    Position* dst_positions = mesh->positions;
    uint32_t* dst_indices = mesh->triangle_indices;

    const Position* src_position = spines.vertices;
    uint32_t base_index = 0;

    for (uint16_t spine = 0; spine < spines.num_spines; spine++) {
        const bool thin = context->guideline_start > 0 &&
            spine >= context->guideline_start;
        const float thickness = thin ? 1.0f : context->config.thickness;
        const uint16_t spine_length = spines.spine_lengths[spine];
        float dx = src_position[1].x - src_position[0].x;
        float dy = src_position[1].y - src_position[0].y;
        float segment_length = sqrtf(dx * dx + dy * dy);
        float invlen = segment_length ? 1.0f / segment_length : 0.0f;
        const float nx = -dy * invlen;
        const float ny = dx * invlen;

        const Position first_src_position = src_position[0];
        const Position last_src_position = src_position[spine_length - 1];

        float ex = nx * thickness / 2;
        float ey = ny * thickness / 2;

        if (closed) {
            const float dx = src_position[0].x - last_src_position.x;
            const float dy = src_position[0].y - last_src_position.y;
            const float segment_length = sqrtf(dx * dx + dy * dy);
            float invlen = segment_length ? 1.0f / segment_length : 0.0f;
            const float pnx = -dy * invlen;
            const float pny = dx * invlen;

            // NOTE: sin(pi / 2 - acos(X) / 2) == sqrt(1 + X) / sqrt(2)
            float extent = 0.5 * thickness;
            const float dotp = (pnx * nx + pny * ny);
            if (dotp < miter_acos_max && dotp > miter_acos_min) {
                const float phi = acos(dotp) / 2;
                const float theta = PAR_PI / 2 - phi;
                extent = PAR_CLAMP(extent / sin(theta), -miter_limit,
                    miter_limit);
            }

            ex = pnx + nx;
            ey = pny + ny;
            const float len = sqrtf(ex * ex + ey * ey);
            invlen = len == 0.0 ? 0.0 : (1.0f / len);
            ex *= invlen * extent;
            ey *= invlen * extent;
        }

        dst_positions[0].x = src_position[0].x + ex;
        dst_positions[0].y = src_position[0].y + ey;
        dst_positions[1].x = src_position[0].x - ex;
        dst_positions[1].y = src_position[0].y - ey;

        float pnx = nx;
        float pny = ny;

        const Position first_dst_positions[2] = {
            dst_positions[0],
            dst_positions[1]
        };

        src_position++;
        dst_positions += 2;

        if (has_annotations) {
            dst_annotations[0].u_along_curve = 0;
            dst_annotations[1].u_along_curve = 0;
            dst_annotations[0].v_across_curve = 1;
            dst_annotations[1].v_across_curve = -1;
            dst_annotations[0].spine_to_edge_x = ex;
            dst_annotations[1].spine_to_edge_x = -ex;
            dst_annotations[0].spine_to_edge_y = ey;
            dst_annotations[1].spine_to_edge_y = -ey;
            dst_annotations += 2;
        }

        float distance_along_spine = segment_length;

        uint16_t segment_index = 1;
        for (; segment_index < spine_length - 1; segment_index++) {

            const float dx = src_position[1].x - src_position[0].x;
            const float dy = src_position[1].y - src_position[0].y;
            const float segment_length = sqrtf(dx * dx + dy * dy);
            float invlen = segment_length ? 1.0f / segment_length : 0.0f;
            const float nx = -dy * invlen;
            const float ny = dx * invlen;

            // NOTE: sin(pi / 2 - acos(X) / 2) == sqrt(1 + X) / sqrt(2)
            float extent = 0.5 * thickness;
            const float dotp = (pnx * nx + pny * ny);
            if (dotp < miter_acos_max && dotp > miter_acos_min) {
                const float phi = acos(dotp) / 2;
                const float theta = PAR_PI / 2 - phi;
                extent = PAR_CLAMP(extent / sin(theta), -miter_limit,
                    miter_limit);
            }

            float ex = pnx + nx;
            float ey = pny + ny;
            const float len = sqrtf(ex * ex + ey * ey);
            invlen = len == 0.0 ? 0.0 : (1.0f / len);
            ex *= invlen * extent;
            ey *= invlen * extent;

            dst_positions[0].x = src_position[0].x + ex;
            dst_positions[0].y = src_position[0].y + ey;
            dst_positions[1].x = src_position[0].x - ex;
            dst_positions[1].y = src_position[0].y - ey;
            src_position++;
            dst_positions += 2;

            pnx = nx;
            pny = ny;

            if (has_annotations) {
                dst_annotations[0].u_along_curve = distance_along_spine;
                dst_annotations[1].u_along_curve = distance_along_spine;
                dst_annotations[0].v_across_curve = 1;
                dst_annotations[1].v_across_curve = -1;
                dst_annotations[0].spine_to_edge_x = ex;
                dst_annotations[1].spine_to_edge_x = -ex;
                dst_annotations[0].spine_to_edge_y = ey;
                dst_annotations[1].spine_to_edge_y = -ey;
                dst_annotations += 2;
            }
            distance_along_spine += segment_length;

            if (wireframe) {
                dst_indices[0] = base_index + (segment_index - 1) * 2;
                dst_indices[1] = base_index + (segment_index - 1) * 2 + 1;
                dst_indices[2] = base_index + (segment_index - 0) * 2;
                dst_indices[3] = base_index + (segment_index - 1) * 2;

                dst_indices[4] = base_index + (segment_index - 0) * 2;
                dst_indices[5] = base_index + (segment_index - 1) * 2 + 1;
                dst_indices[6] = base_index + (segment_index - 0) * 2 + 1;
                dst_indices[7] = base_index + (segment_index - 0) * 2;
                dst_indices += 8;
            } else {
                dst_indices[0] = base_index + (segment_index - 1) * 2;
                dst_indices[1] = base_index + (segment_index - 1) * 2 + 1;
                dst_indices[2] = base_index + (segment_index - 0) * 2;

                dst_indices[3] = base_index + (segment_index - 0) * 2;
                dst_indices[4] = base_index + (segment_index - 1) * 2 + 1;
                dst_indices[5] = base_index + (segment_index - 0) * 2 + 1;
                dst_indices += 6;
            }
        }

        ex = pnx * thickness / 2;
        ey = pny * thickness / 2;

        if (closed) {
            const float dx = first_src_position.x - src_position[0].x;
            const float dy = first_src_position.y - src_position[0].y;
            segment_length = sqrtf(dx * dx + dy * dy);
            float invlen = segment_length ? 1.0f / segment_length : 0.0f;
            const float nx = -dy * invlen;
            const float ny = dx * invlen;

            // NOTE: sin(pi / 2 - acos(X) / 2) == sqrt(1 + X) / sqrt(2)
            float extent = 0.5 * thickness;
            const float dotp = (pnx * nx + pny * ny);
            if (dotp < miter_acos_max && dotp > miter_acos_min) {
                const float phi = acos(dotp) / 2;
                const float theta = PAR_PI / 2 - phi;
                extent = PAR_CLAMP(extent / sin(theta), -miter_limit,
                    miter_limit);
            }

            ex = pnx + nx;
            ey = pny + ny;
            const float len = sqrtf(ex * ex + ey * ey);
            invlen = len == 0.0 ? 0.0 : (1.0f / len);
            ex *= invlen * extent;
            ey *= invlen * extent;
        }

        dst_positions[0].x = src_position[0].x + ex;
        dst_positions[0].y = src_position[0].y + ey;
        dst_positions[1].x = src_position[0].x - ex;
        dst_positions[1].y = src_position[0].y - ey;
        src_position++;
        dst_positions += 2;

        pnx = nx;
        pny = ny;

        if (has_annotations) {
            dst_annotations[0].u_along_curve = distance_along_spine;
            dst_annotations[1].u_along_curve = distance_along_spine;
            dst_annotations[0].v_across_curve = 1;
            dst_annotations[1].v_across_curve = -1;
            dst_annotations[0].spine_to_edge_x = ex;
            dst_annotations[1].spine_to_edge_x = -ex;
            dst_annotations[0].spine_to_edge_y = ey;
            dst_annotations[1].spine_to_edge_y = -ey;
            dst_annotations += 2;
        }

        if (wireframe) {
            dst_indices[0] = base_index + (segment_index - 1) * 2;
            dst_indices[1] = base_index + (segment_index - 1) * 2 + 1;
            dst_indices[2] = base_index + (segment_index - 0) * 2;
            dst_indices[3] = base_index + (segment_index - 1) * 2;

            dst_indices[4] = base_index + (segment_index - 0) * 2;
            dst_indices[5] = base_index + (segment_index - 1) * 2 + 1;
            dst_indices[6] = base_index + (segment_index - 0) * 2 + 1;
            dst_indices[7] = base_index + (segment_index - 0) * 2;
            dst_indices += 8;
        } else {
            dst_indices[0] = base_index + (segment_index - 1) * 2;
            dst_indices[1] = base_index + (segment_index - 1) * 2 + 1;
            dst_indices[2] = base_index + (segment_index - 0) * 2;

            dst_indices[3] = base_index + (segment_index - 0) * 2;
            dst_indices[4] = base_index + (segment_index - 1) * 2 + 1;
            dst_indices[5] = base_index + (segment_index - 0) * 2 + 1;
            dst_indices += 6;
        }

        if (closed) {
            segment_index++;
            distance_along_spine += segment_length;

            dst_positions[0] = first_dst_positions[0];
            dst_positions[1] = first_dst_positions[1];
            dst_positions += 2;

            if (has_annotations) {
                dst_annotations[0].u_along_curve = distance_along_spine;
                dst_annotations[1].u_along_curve = distance_along_spine;
                dst_annotations[0].v_across_curve = 1;
                dst_annotations[1].v_across_curve = -1;
                dst_annotations[0].spine_to_edge_x = ex;
                dst_annotations[1].spine_to_edge_x = -ex;
                dst_annotations[0].spine_to_edge_y = ey;
                dst_annotations[1].spine_to_edge_y = -ey;
                dst_annotations += 2;
            }

            if (wireframe) {
                dst_indices[0] = base_index + (segment_index - 1) * 2;
                dst_indices[1] = base_index + (segment_index - 1) * 2 + 1;
                dst_indices[2] = base_index + (segment_index - 0) * 2;
                dst_indices[3] = base_index + (segment_index - 1) * 2;

                dst_indices[4] = base_index + (segment_index - 0) * 2;
                dst_indices[5] = base_index + (segment_index - 1) * 2 + 1;
                dst_indices[6] = base_index + (segment_index - 0) * 2 + 1;
                dst_indices[7] = base_index + (segment_index - 0) * 2;
                dst_indices += 8;
            } else {
                dst_indices[0] = base_index + (segment_index - 1) * 2;
                dst_indices[1] = base_index + (segment_index - 1) * 2 + 1;
                dst_indices[2] = base_index + (segment_index - 0) * 2;

                dst_indices[3] = base_index + (segment_index - 0) * 2;
                dst_indices[4] = base_index + (segment_index - 1) * 2 + 1;
                dst_indices[5] = base_index + (segment_index - 0) * 2 + 1;
                dst_indices += 6;
            }
        }

        base_index += spine_length * 2 + (closed ? 2  : 0);

        const uint16_t nverts = spine_length + (closed ? 1 : 0);

        if (has_lengths) {
            for (uint16_t i = 0; i < nverts; i++) {
                dst_lengths[0] = distance_along_spine;
                dst_lengths[1] = distance_along_spine;
                dst_lengths += 2;
            }
        }

        // Go back through the curve and fix up the U coordinates.
        if (has_annotations) {
            const float invlength = 1.0f / distance_along_spine;
            const float invcount = 1.0f / spine_length;
            switch (context->config.u_mode) {
            case PAR_U_MODE_DISTANCE:
                break;
            case PAR_U_MODE_NORMALIZED_DISTANCE:
                dst_annotations -= nverts * 2;
                for (uint16_t i = 0; i < nverts; i++) {
                    dst_annotations[0].u_along_curve *= invlength;
                    dst_annotations[1].u_along_curve *= invlength;
                    dst_annotations += 2;
                }
                break;
            case PAR_U_MODE_SEGMENT_INDEX:
                dst_annotations -= nverts * 2;
                for (uint16_t i = 0; i < nverts; i++) {
                    dst_annotations[0].u_along_curve = i;
                    dst_annotations[1].u_along_curve = i;
                    dst_annotations += 2;
                }
                break;
            case PAR_U_MODE_SEGMENT_FRACTION:
                dst_annotations -= nverts * 2;
                for (uint16_t i = 0; i < nverts; i++) {
                    dst_annotations[0].u_along_curve = invcount * i;
                    dst_annotations[1].u_along_curve = invcount * i;
                    dst_annotations += 2;
                }
                break;
            }
        }
    }

    assert(src_position - spines.vertices == spines.num_vertices);
    assert(dst_positions - mesh->positions == mesh->num_vertices);
    assert(dst_indices - mesh->triangle_indices ==
        mesh->num_triangles * ind_per_tri);

    if (context->config.flags & PARSL_FLAG_RANDOM_OFFSETS) {
        pa_clear(mesh->random_offsets);
        pa_add(mesh->random_offsets, mesh->num_vertices);
        float* pvertex = mesh->random_offsets;
        for (uint16_t spine = 0; spine < spines.num_spines; spine++) {
            const uint16_t num_segments = spines.spine_lengths[spine];
            const float r = (float) rand() / RAND_MAX;
            for (uint16_t segment = 0; segment < num_segments; segment++) {
                *pvertex++ = r;
            }
        }
    }

    return mesh;
}

// This function is designed to be called in two passes. In the first pass, the
// points pointer is null, so this simply determines the number of required
// points to fulfill the flatness criterion. On the second pass, points is
// non-null it actually writes out the point positions.
static void parsl__tesselate_cubic(
    parsl_position* points, uint32_t* num_points,
    float x0, float y0, float x1, float y1,
    float x2, float y2, float x3, float y3,
    float max_flatness_squared, int recursion_depth)
{
    float dx0 = x1-x0;
    float dy0 = y1-y0;
    float dx1 = x2-x1;
    float dy1 = y2-y1;
    float dx2 = x3-x2;
    float dy2 = y3-y2;
    float dx = x3-x0;
    float dy = y3-y0;
    float longlen = (float) (sqrt(dx0*dx0 + dy0*dy0) + sqrt(dx1*dx1 + dy1*dy1) +
        sqrt(dx2*dx2 + dy2*dy2));
    float shortlen = (float) sqrt(dx*dx+dy*dy);
    float flatness_squared = longlen*longlen - shortlen*shortlen;

    if (recursion_depth > PARSL_MAX_RECURSION) {
        return;
    }

    if (flatness_squared > max_flatness_squared) {
        const float x01 = (x0+x1) / 2;
        const float y01 = (y0+y1) / 2;
        const float x12 = (x1+x2) / 2;
        const float y12 = (y1+y2) / 2;
        const float x23 = (x2+x3) / 2;
        const float y23 = (y2+y3) / 2;
        const float xa = (x01+x12) / 2;
        const float ya = (y01+y12) / 2;
        const float xb = (x12+x23) / 2;
        const float yb = (y12+y23) / 2;
        const float mx = (xa+xb) / 2;
        const float my = (ya+yb) / 2;
        parsl__tesselate_cubic(points, num_points, x0,y0, x01,y01, xa,ya, mx,my,
            max_flatness_squared, recursion_depth + 1);
        parsl__tesselate_cubic(points, num_points, mx,my, xb,yb, x23,y23, x3,y3,
            max_flatness_squared, recursion_depth + 1);
        return;
    }

    int n = *num_points;
    if (points) {
        points[n].x = x3;
        points[n].y = y3;
    }
    *num_points = n + 1;
}

// This function is designed to be called in two passes. In the first pass, the
// points pointer is null, so this simply determines the number of required
// points to fulfill the flatness criterion. On the second pass, points is
// non-null it actually writes out the point positions.
static void parsl__tesselate_quadratic(
    parsl_position* points, uint32_t* num_points,
    float x0, float y0, float x1, float y1, float x2, float y2,
    float max_flatness_squared, int recursion_depth)
{
    const float mx = (x0 + 2 * x1 + x2) / 4;
    const float my = (y0 + 2 * y1 + y2) / 4;
    const float dx = (x0 + x2) / 2 - mx;
    const float dy = (y0 + y2) / 2 - my;
    const float flatness_squared = dx * dx + dy * dy;

    if (recursion_depth++ > PARSL_MAX_RECURSION) {
        return;
    }

    if (flatness_squared > max_flatness_squared) {
        parsl__tesselate_quadratic(points, num_points, x0,y0,
            (x0 + x1) / 2.0f, (y0 + y1) / 2.0f,
            mx, my,
            max_flatness_squared, recursion_depth);
        parsl__tesselate_quadratic(points, num_points, mx,my,
            (x1 + x2) / 2.0f, (y1 + y2) / 2.0f,
            x2, y2,
            max_flatness_squared, recursion_depth);
        return;
    }

    int n = *num_points;
    if (points) {
        points[n].x = x2;
        points[n].y = y2;
    }
    *num_points = n + 1;
}

parsl_mesh* parsl_mesh_from_curves_cubic(parsl_context* context,
    parsl_spine_list source_spines)
{
    float max_flatness = context->config.curves_max_flatness;
    if (max_flatness == 0) {
        max_flatness = 1.0f;
    }
    const float max_flatness_squared = max_flatness * max_flatness;
    parsl_spine_list* target_spines = &context->curve_spines;
    const bool has_guides = context->config.flags & PARSL_FLAG_CURVE_GUIDES;

    // Determine the number of spines in the target list.
    target_spines->num_spines = source_spines.num_spines;
    if (has_guides) {
        for (uint32_t spine = 0; spine < source_spines.num_spines; spine++) {
            uint32_t spine_length = source_spines.spine_lengths[spine];
            uint32_t num_piecewise = 1 + (spine_length - 4) / 2;
            target_spines->num_spines += num_piecewise * 2;
        }
    }
    pa_clear(target_spines->spine_lengths);
    pa_add(target_spines->spine_lengths, target_spines->num_spines);

    // First pass: determine the number of required vertices.
    uint32_t total_required_spine_points = 0;
    const parsl_position* psource = source_spines.vertices;
    for (uint32_t spine = 0; spine < source_spines.num_spines; spine++) {

        // Source vertices look like: P1 C1 C2 P2 [C2 P2]*
        uint32_t spine_length = source_spines.spine_lengths[spine];
        assert(spine_length >= 4);
        assert((spine_length % 2) == 0);
        uint32_t num_piecewise = 1 + (spine_length - 4) / 2;

        // First piecewise curve.
        uint32_t num_required_spine_points = 1;
        parsl__tesselate_cubic(NULL, &num_required_spine_points,
            psource[0].x, psource[0].y, psource[1].x, psource[1].y,
            psource[2].x, psource[2].y, psource[3].x, psource[3].y,
            max_flatness_squared, 0);
        psource += 4;

        // Subsequent piecewise curves.
        for (uint32_t piecewise = 1; piecewise < num_piecewise; piecewise++) {
            parsl_position p1 = psource[-1];
            parsl_position previous_c2 = psource[-2];
            parsl_position c1 = parsl__sub(p1, parsl__sub(previous_c2, p1));
            parsl_position c2 = psource[0];
            parsl_position p2 = psource[1];
            parsl__tesselate_cubic(NULL, &num_required_spine_points,
                p1.x, p1.y, c1.x, c1.y, c2.x, c2.y, p2.x, p2.y,
                max_flatness_squared, 0);
            psource += 2;
        }

        target_spines->spine_lengths[spine] = num_required_spine_points;
        total_required_spine_points += num_required_spine_points;
    }

    if (has_guides) {
        uint32_t nsrcspines = source_spines.num_spines;
        uint16_t* guide_lengths = &target_spines->spine_lengths[nsrcspines];
        for (uint32_t spine = 0; spine < nsrcspines; spine++) {
            uint32_t spine_length = source_spines.spine_lengths[spine];
            uint32_t num_piecewise = 1 + (spine_length - 4) / 2;
            for (uint32_t pw = 0; pw < num_piecewise; pw++) {
                guide_lengths[0] = 2;
                guide_lengths[1] = 2;
                guide_lengths += 2;
                total_required_spine_points += 4;
            }
        }
    }

    // Allocate memory.
    target_spines->num_vertices = total_required_spine_points;
    pa_clear(target_spines->vertices);
    pa_add(target_spines->vertices, total_required_spine_points);

    // Second pass: write out the data.
    psource = source_spines.vertices;
    parsl_position* ptarget = target_spines->vertices;
    for (uint32_t spine = 0; spine < source_spines.num_spines; spine++) {

        // Source vertices look like: P1 C1 C2 P2 [C2 P2]*
        uint32_t spine_length = source_spines.spine_lengths[spine];
        uint32_t num_piecewise = 1 + (spine_length - 4) / 2;

        __attribute__((unused))
        parsl_position* target_spine_start = ptarget;

        // First piecewise curve.
        ptarget[0].x = psource[0].x;
        ptarget[0].y = psource[0].y;
        ptarget++;
        uint32_t num_written_points = 0;
        parsl__tesselate_cubic(ptarget, &num_written_points,
            psource[0].x, psource[0].y, psource[1].x, psource[1].y,
            psource[2].x, psource[2].y, psource[3].x, psource[3].y,
            max_flatness_squared, 0);
        psource += 4;
        ptarget += num_written_points;

        // Subsequent piecewise curves.
        for (uint32_t piecewise = 1; piecewise < num_piecewise; piecewise++) {
            parsl_position p1 = psource[-1];
            parsl_position previous_c2 = psource[-2];
            parsl_position c1 = parsl__sub(p1, parsl__sub(previous_c2, p1));
            parsl_position c2 = psource[0];
            parsl_position p2 = psource[1];
            num_written_points = 0;
            parsl__tesselate_cubic(ptarget, &num_written_points,
                p1.x, p1.y, c1.x, c1.y, c2.x, c2.y, p2.x, p2.y,
                max_flatness_squared, 0);
            psource += 2;
            ptarget += num_written_points;
        }

        __attribute__((unused))
        uint32_t num_written = ptarget - target_spine_start;
        assert(num_written == (uint32_t) target_spines->spine_lengths[spine]);
    }

    // Source vertices look like: P1 C1 C2 P2 [C2 P2]*
    if (has_guides) {
        uint32_t nsrcspines = source_spines.num_spines;
        context->guideline_start = nsrcspines;
        psource = source_spines.vertices;
        for (uint32_t spine = 0; spine < nsrcspines; spine++) {
            uint32_t spine_length = source_spines.spine_lengths[spine];
            uint32_t num_piecewise = 1 + (spine_length - 4) / 2;
            *ptarget++ = psource[0];
            *ptarget++ = psource[1];
            *ptarget++ = psource[2];
            *ptarget++ = psource[3];
            psource += 4;
            for (uint32_t pw = 1; pw < num_piecewise; pw++) {
                parsl_position p1 = psource[-1];
                parsl_position previous_c2 = psource[-2];
                parsl_position c1 = parsl__sub(p1, parsl__sub(previous_c2, p1));
                parsl_position c2 = psource[0];
                parsl_position p2 = psource[1];
                *ptarget++ = p1;
                *ptarget++ = c1;
                *ptarget++ = p2;
                *ptarget++ = c2;
                psource += 2;
            }
        }
    }

    assert(ptarget - target_spines->vertices == total_required_spine_points);
    parsl_mesh_from_lines(context, context->curve_spines);
    context->guideline_start = 0;
    return &context->result;
}

parsl_mesh* parsl_mesh_from_curves_quadratic(parsl_context* context,
    parsl_spine_list source_spines)
{
    float max_flatness = context->config.curves_max_flatness;
    if (max_flatness == 0) {
        max_flatness = 1.0f;
    }
    const float max_flatness_squared = max_flatness * max_flatness;
    parsl_spine_list* target_spines = &context->curve_spines;
    const bool has_guides = context->config.flags & PARSL_FLAG_CURVE_GUIDES;

    // Determine the number of spines in the target list.
    target_spines->num_spines = source_spines.num_spines;
    if (has_guides) {
        target_spines->num_spines += source_spines.num_spines;
    }
    pa_clear(target_spines->spine_lengths);
    pa_add(target_spines->spine_lengths, target_spines->num_spines);

    // First pass: determine the number of required vertices.
    uint32_t total_required_spine_points = 0;
    const parsl_position* psource = source_spines.vertices;
    for (uint32_t spine = 0; spine < source_spines.num_spines; spine++) {

        // Source vertices look like: PT C PT [C PT]*
        uint32_t spine_length = source_spines.spine_lengths[spine];
        assert(spine_length >= 3);
        assert((spine_length % 2) == 1);
        uint32_t num_piecewise = 1 + (spine_length - 3) / 2;

        // First piecewise curve.
        uint32_t num_required_spine_points = 1;
        parsl__tesselate_quadratic(NULL, &num_required_spine_points,
            psource[0].x, psource[0].y, psource[1].x, psource[1].y,
            psource[2].x, psource[2].y,  max_flatness_squared, 0);
        psource += 3;

        // Subsequent piecewise curves.
        for (uint32_t piecewise = 1; piecewise < num_piecewise; piecewise++) {
            parsl_position p1 = psource[-1];
            parsl_position c1 = psource[0];
            parsl_position p2 = psource[1];
            parsl__tesselate_quadratic(NULL, &num_required_spine_points,
                p1.x, p1.y, c1.x, c1.y, p2.x, p2.y,  max_flatness_squared, 0);
            psource += 2;
        }

        target_spines->spine_lengths[spine] = num_required_spine_points;
        total_required_spine_points += num_required_spine_points;
    }

    if (has_guides) {
        uint32_t nsrcspines = source_spines.num_spines;
        uint16_t* guide_lengths = &target_spines->spine_lengths[nsrcspines];
        for (uint32_t spine = 0; spine < nsrcspines; spine++) {
            uint32_t spine_length = source_spines.spine_lengths[spine];
            uint32_t num_piecewise = 1 + (spine_length - 3) / 2;
            guide_lengths[0] = 3 + (num_piecewise - 1) * 2;
            total_required_spine_points += guide_lengths[0];
            guide_lengths++;
        }
    }

    // Allocate memory.
    target_spines->num_vertices = total_required_spine_points;
    pa_clear(target_spines->vertices);
    pa_add(target_spines->vertices, total_required_spine_points);

    // Second pass: write out the data.
    psource = source_spines.vertices;
    parsl_position* ptarget = target_spines->vertices;
    for (uint32_t spine = 0; spine < source_spines.num_spines; spine++) {

        // Source vertices look like: PT C PT [C PT]*
        uint32_t spine_length = source_spines.spine_lengths[spine];
        uint32_t num_piecewise = 1 + (spine_length - 3) / 2;

        __attribute__((unused))
        parsl_position* target_spine_start = ptarget;

        // First piecewise curve.
        ptarget[0].x = psource[0].x;
        ptarget[0].y = psource[0].y;
        ptarget++;
        uint32_t num_written_points = 0;
        parsl__tesselate_quadratic(ptarget, &num_written_points,
            psource[0].x, psource[0].y, psource[1].x, psource[1].y,
            psource[2].x, psource[2].y, max_flatness_squared, 0);
        psource += 3;
        ptarget += num_written_points;

        // Subsequent piecewise curves.
        for (uint32_t piecewise = 1; piecewise < num_piecewise; piecewise++) {
            parsl_position p1 = psource[-1];
            parsl_position c1 = psource[0];
            parsl_position p2 = psource[1];
            num_written_points = 0;
            parsl__tesselate_quadratic(ptarget, &num_written_points,
                p1.x, p1.y, c1.x, c1.y, p2.x, p2.y, max_flatness_squared, 0);
            psource += 2;
            ptarget += num_written_points;
        }

        __attribute__((unused))
        uint32_t num_written = ptarget - target_spine_start;
        assert(num_written == (uint32_t) target_spines->spine_lengths[spine]);
    }

    // Source vertices look like: PT C PT [C PT]*
    if (has_guides) {
        uint32_t nsrcspines = source_spines.num_spines;
        context->guideline_start = nsrcspines;
        psource = source_spines.vertices;
        for (uint32_t spine = 0; spine < nsrcspines; spine++) {
            uint32_t spine_length = source_spines.spine_lengths[spine];
            uint32_t num_piecewise = 1 + (spine_length - 3) / 2;
            *ptarget++ = psource[0];
            *ptarget++ = psource[1];
            *ptarget++ = psource[2];
            psource += 3;
            for (uint32_t pw = 1; pw < num_piecewise; pw++) {
                *ptarget++ = psource[0];
                *ptarget++ = psource[1];
                psource += 2;
            }
        }
    }

    assert(ptarget - target_spines->vertices == total_required_spine_points);
    parsl_mesh_from_lines(context, context->curve_spines);
    context->guideline_start = 0;
    return &context->result;
}

static unsigned int par__randhash(unsigned int seed) {
    unsigned int i = (seed ^ 12345391u) * 2654435769u;
    i ^= (i << 6) ^ (i >> 26);
    i *= 2654435769u;
    i += (i << 5) ^ (i >> 12);
    return i;
}

static float par__randhashf(unsigned int seed, float a, float b) {
    return (b - a) * par__randhash(seed) / (float) UINT_MAX + a;
}

static parsl_position par__sample_annulus(float radius, parsl_position center,
    int* seedptr) {
    unsigned int seed = *seedptr;
    parsl_position r;
    float rscale = 1.0f / UINT_MAX;
    while (1) {
        r.x = 4 * rscale * par__randhash(seed++) - 2;
        r.y = 4 * rscale * par__randhash(seed++) - 2;
        float r2 = parsl__dot(r, r);
        if (r2 > 1 && r2 <= 4) {
            break;
        }
    }
    *seedptr = seed;
    return (parsl_position) {
               r.x * radius + center.x,
                   r.y * radius + center.y
    };
}

#define GRIDF(vec) \
    grid [(int) (vec.x * invcell) + ncols * (int) (vec.y * invcell)]

#define GRIDI(vec) grid [(int) vec.y * ncols + (int) vec.x]

static parsl_position* par__generate_pts(float width, float height,
    float radius, int seed, parsl_position* result) {

    int maxattempts = 30;
    float rscale = 1.0f / UINT_MAX;
    parsl_position rvec;
    rvec.x = rvec.y = radius;
    float r2 = radius * radius;

    // Acceleration grid.
    float cellsize = radius / sqrtf(2);
    float invcell = 1.0f / cellsize;
    int ncols = ceil(width * invcell);
    int nrows = ceil(height * invcell);
    int maxcol = ncols - 1;
    int maxrow = nrows - 1;
    int ncells = ncols * nrows;
    int* grid = (int*) PAR_MALLOC(int, ncells);
    for (int i = 0; i < ncells; i++) {
        grid[i] = -1;
    }

    // Active list and resulting sample list.
    int* actives = (int*) PAR_MALLOC(int, ncells);
    int nactives = 0;

    pa_clear(result);
    pa_add(result, ncells);

    parsl_position* samples = result;
    int nsamples = 0;

    // First sample.
    parsl_position pt;
    pt.x = width * par__randhash(seed++) * rscale;
    pt.y = height * par__randhash(seed++) * rscale;
    GRIDF(pt) = actives[nactives++] = nsamples;
    samples[nsamples++] = pt;

    while (nsamples < ncells) {
        int aindex = PAR_MIN(par__randhashf(seed++, 0, nactives),
            nactives - 1.0f);
        int sindex = actives[aindex];
        int found = 0;
        parsl_position j, minj, maxj, delta;
        int attempt;
        for (attempt = 0; attempt < maxattempts && !found; attempt++) {
            pt = par__sample_annulus(radius, samples[sindex], &seed);

            // Check that this sample is within bounds.
            if (pt.x < 0 || pt.x >= width || pt.y < 0 || pt.y >= height) {
                continue;
            }

            // Test proximity to nearby samples.
            maxj = parsl_mul(parsl__add(pt, rvec), invcell);
            minj = parsl_mul(parsl__sub(pt, rvec), invcell);

            minj.x = PAR_CLAMP((int) minj.x, 0, maxcol);
            minj.y = PAR_CLAMP((int) minj.y, 0, maxrow);
            maxj.x = PAR_CLAMP((int) maxj.x, 0, maxcol);
            maxj.y = PAR_CLAMP((int) maxj.y, 0, maxrow);
            int reject = 0;
            for (j.y = minj.y; j.y <= maxj.y && !reject; j.y++) {
                for (j.x = minj.x; j.x <= maxj.x && !reject; j.x++) {
                    int entry = GRIDI(j);
                    if (entry > -1 && entry != sindex) {
                        delta = parsl__sub(samples[entry], pt);
                        if (parsl__dot(delta, delta) < r2) {
                            reject = 1;
                        }
                    }
                }
            }
            if (reject) {
                continue;
            }
            found = 1;
        }
        if (found) {
            GRIDF(pt) = actives[nactives++] = nsamples;
            samples[nsamples++] = pt;
        } else {
            if (--nactives <= 0) {
                break;
            }
            actives[aindex] = actives[nactives];
        }
    }

    pa___n(result) = nsamples * 2;

    PAR_FREE(grid);
    PAR_FREE(actives);
    return result;
}

#undef GRIDF
#undef GRIDI

parsl_mesh* parsl_mesh_from_streamlines(parsl_context* context,
    parsl_advection_callback advect, uint32_t first_tick, uint32_t num_ticks,
    void* userdata)
{
    const int seed = 42;
    const parsl_viewport vp = context->config.streamlines_seed_viewport;
    const float radius = context->config.streamlines_seed_spacing;
    float width = vp.right - vp.left;
    float height = vp.bottom - vp.top;

    if (context->streamline_seeds == NULL) {
        context->streamline_seeds = par__generate_pts(width, height, radius,
            seed, context->streamline_seeds);
        uint32_t num_points = pa_count(context->streamline_seeds);
        for (uint32_t p = 0; p < num_points; ++p) {
            context->streamline_seeds[p].x += vp.left;
            context->streamline_seeds[p].y += vp.top;
        }
    }

    uint32_t num_points = pa_count(context->streamline_seeds);
    pa_clear(context->streamline_points);
    pa_add(context->streamline_points, num_points);

    parsl_position* points = context->streamline_points;
    memcpy(points, context->streamline_seeds,
        num_points * sizeof(parsl_position));

    context->streamline_spines.num_spines = num_points;
    pa_clear(context->streamline_spines.spine_lengths);
    pa_add(context->streamline_spines.spine_lengths, num_points);
    uint16_t* lengths = context->streamline_spines.spine_lengths;
    for (uint32_t i = 0; i < num_points; i++) {
        lengths[i] = num_ticks;
    }

    context->streamline_spines.num_vertices = num_points * num_ticks;
    pa_clear(context->streamline_spines.vertices);
    pa_add(context->streamline_spines.vertices, num_points * num_ticks);
    parsl_position* vertices = context->streamline_spines.vertices;

    for (uint32_t tick = 0; tick < first_tick; tick++) {
        for (uint32_t i = 0; i < num_points; i++) {
            advect(&points[i], userdata);
        }
    }

    parsl_position* pvertices = vertices;
    for (uint32_t i = 0; i < num_points; i++) {
        for (uint32_t tick = 0; tick < num_ticks; ++tick) {
            advect(&points[i], userdata);
            *pvertices++ = points[i];
        }
    }

    parsl_mesh_from_lines(context, context->streamline_spines);
    return &context->result;
}

#endif // PAR_STREAMLINES_IMPLEMENTATION
#endif // PAR_STREAMLINES_H

// par_streamlines is distributed under the MIT license:
//
// Copyright (c) 2019 Philip Rideout
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.