john/prosperon
1
0
Fork 0
Code Issues Pull requests Actions Releases Activity
prosperon/source/engine/thirdparty/par/par_streamlines.h
John Alanbrook 9a325543ae par and mipmaps
2023-05-13 03:21:11 +00:00

1260 lines
47 KiB
C
Raw Blame History

// 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.
Reference in a new issue View git blame Copy permalink