dsga : Data Structures for Geometric Algorithms

dsga is a single header-only C++20 library that implements the vectors and matrices from the OpenGL Shading Language 4.6 specification (pdf | html). It is inspired by the spec, but does deviate in some small ways, mostly to make it work well in C++20. It is intended to be used for array programming other than rendering. Our requirements in general are for things like 3D CAD/CAM applications and other geometric and algebraic things. See motivation for more details. This library does not use SIMD instructions or types under the hood, beyond whatever the compiler provides through optimization.

Home

https://github.com/davidbrowne/dsga

Current Version

v2.2.2

Latest Major Changes

v2.2.1
- Replaced home-brew asserts with exceptions. Attempting to make safer through bounds checking and other input checking, enforced by throwing exceptions.
v2.2.0
- Reverted major changes between v2.0.5 and v2.1.4, so we no longer have the view structs that wrap a pointer. Wrapping a pointer turned the data structures from owning to non-owning for the view structs, but we want the vector and matrix structs to be owning. The point of the view experiment was built on lack of insight on the nature of owning vs. non-owning and what this library was trying to achieve.
v2.0.5
- Fixed wrong matrix type (reversed dimensions) being returned from outerProduct().
v2.0.4
- Renamed logicalNot() to compNot(). Deprecated logicalNot().
- Added compAnd() and compOr() functions to complement compNot().
- Added missing scalar versions of non-geometric vector functions.
v2.0.3
- Tolerance checking functions moved to examples/tolerance.hxx.
v2.0.2
- Potentially breaking change: removed an implicit dsga::basic_matrix constructor, now requiring the use of a constructor that is explicit.
v2.0.1
- Added query() function (not in GLSL nor std::valarray) to vector_base. It works like apply(), but expects a boolean predicate, and returns a vector of boolean values instead of element type T.

Tested Compilers

Regularly Tested

Microsoft Visual Studio 2022 v17.12.3
gcc v14.2.0
clang v19.1.5
icx v2024.1.0

Minimum Version

Microsoft Visual Studio 2022 v17.x
gcc v11.4
clang v16.0.6
icx v2023.1.0 - using Compiler Explorer for basic compilation test, but test suite not run

Some Quick Examples

// get a 2D vector that is perpendicular (rotated 90 degrees counter-clockwise)
// to a 2D vector in the plane
template <dsga::floating_point_scalar T>
constexpr auto get_perpendicular1(const dsga::basic_vector<T, 2> &some_vec) noexcept
{
    auto cos90 = 0.0f;
    auto sin90 = 1.0f;

    // rotation matrix -- components in column major order
    return dsga::basic_matrix<T, 2, 2>(cos90, sin90, -sin90, cos90) * some_vec;
}

// same as above, different implementation
template <dsga::floating_point_scalar T>
constexpr auto get_perpendicular2(const dsga::basic_vector<T, 2> &some_vec) noexcept
{
    return dsga::basic_vector<T, 2>(-1, 1) * some_vec.yx;
}

// gives closest projection point from point to a line made from line segment p1 <=> p2
constexpr auto project_to_line1(const dsga::dvec3 &point,
                                const dsga::dvec3 &p1,
                                const dsga::dvec3 &p2) noexcept
{
    auto hyp = point - p1;
    auto v1 = p2 - p1;
    auto t = dsga::dot(hyp, v1) / dsga::dot(v1, v1);

    return p1 + (t * v1);
}

// same as above, different implementation
constexpr auto project_to_line2(const dsga::dvec3 &point,
                                const dsga::dvec3 &p1,
                                const dsga::dvec3 &p2) noexcept
{
    auto hyp = point - p1;
    auto v1 = p2 - p1;
    return p1 + dsga::outerProduct(v1, v1) * hyp / dsga::dot(v1, v1);
}

//
// evaluate a 2D cubic bezier curve at t
//

#if LINEAR_INTERPOLATE

// cubic bezier linear interpolation, one ordinate at a time, e.g., x, y, z, or w
// very slow implementation (de Casteljau algorithm), but illustrates the library
constexpr auto single_ordinate_cubic_bezier_eval(const dsga::vec4 &cubic_control_points, float t) noexcept
{
    auto quadratic_control_points = dsga::mix(cubic_control_points.xyz, cubic_control_points.yzw, t);
    auto linear_control_points = dsga::mix(quadratic_control_points.xy, quadratic_control_points.yz, t);
    return dsga::mix(linear_control_points.x, linear_control_points.y, t);
}

#else

// ~10-25x faster - Bernstein polynomials
constexpr auto single_ordinate_cubic_bezier_eval(const dsga::vec4 &cubic_control_points, T t) noexcept
{
    auto t_complement = T(1) - t;
    return
        t_complement * t_complement * t_complement * cubic_control_points[0] +
        T(3) * t * t_complement * t_complement * cubic_control_points[1] +
        T(3) * t * t * t_complement * cubic_control_points[2] +
        t * t * t * cubic_control_points[3];
}

#endif

// main cubic bezier eval function -- takes 2D control points with float values.
// returns the 2D point on the curve at t
constexpr auto simple_cubic_bezier_eval(dsga::vec2 p0, dsga::vec2 p1, dsga::vec2 p2, dsga::vec2 p3, float t) noexcept
{
    // each control point is a column of the matrix.
    // the rows represent x coords and y coords.
    auto AoS = dsga::mat4x2(p0, p1, p2, p3);

    // lambda pack wrapper -- would be better solution if vector size was generic
    return [&]<std::size_t ...Is>(std::index_sequence<Is...>) noexcept
    {
        return dsga::vec2(single_ordinate_cubic_bezier_eval(AoS.row(Is), t)...);
    }(std::make_index_sequence<2>{});
}

//
// find the minimum positive angle between 2 vectors and/or indexed vectors (swizzles).
// Uses base class for vector types to be inclusive to both types.
// 2D or 3D only.
//

template <bool W1, dsga::floating_point_scalar T, std::size_t C, class D1, bool W2, class D2>
requires ((C > 1) && (C < 4))
auto angle_between(const dsga::vector_base<W1, T, C, D1> &v1,
                   const dsga::vector_base<W2, T, C, D2> &v2)
{
    auto a = v1 * dsga::length(v2);
    auto b = v2 * dsga::length(v1);
    auto numerator = dsga::length(a - b);
    auto denominator = dsga::length(a + b);

    if (numerator == T(0))
        return T(0);
    else if (denominator == T(0))
        return std::numbers::pi_v<T>;

    return T(2) * std::atan(numerator / denominator);
}

//
// STL file format read/write helpers
//

// make sure data has no infinities or NaNs
constexpr bool definite_coordinate_triple(const dsga::vec3 &data) noexcept
{
    return !(dsga::any(dsga::isinf(data)) || dsga::any(dsga::isnan(data)));
}

// make sure normal vector has no infinities or NaNs and is not the zero-vector { 0, 0, 0 }
constexpr bool valid_normal_vector(const dsga::vec3 &normal) noexcept
{
    return definite_coordinate_triple(normal) && dsga::any(dsga::notEqual(normal, dsga::vec3(0)));
}

// not checking for positive-only first octant data -- we are allowing zeros and negative values
constexpr bool valid_vertex_relaxed(const dsga::vec3 &vertex) noexcept
{
    return definite_coordinate_triple(vertex);
}

// strict version where all vertex coordinates must be positive-definite
constexpr bool valid_vertex_strict(const dsga::vec3 &vertex) noexcept
{
    return definite_coordinate_triple(vertex) && dsga::all(dsga::greaterThan(vertex, dsga::vec3(0)));
}

// right-handed unit normal vector for a triangle facet,
// inputs are triangle vertices in counter-clockwise order
constexpr dsga::vec3 right_handed_normal(const dsga::vec3 &v1, const dsga::vec3 &v2, const dsga::vec3 &v3) noexcept
{
    return dsga::normalize(dsga::cross(v2 - v1, v3 - v1));
}

//
// cross product
//

// arguments are of the vector_base class type, and this function will be used if any passed argument is of type indexed_vector
template <bool W1, dsga::floating_point_scalar T1, typename D1, bool W2, dsga::floating_point_scalar T2, typename D2>
[[nodiscard]] constexpr auto cross(const dsga::vector_base<W1, T1, 3, D1> &a,
                                   const dsga::vector_base<W2, T2, 3, D2> &b) noexcept
{
    // CTAD gets us the type and size for the vector
    return dsga::basic_vector((a[1] * b[2]) - (b[1] * a[2]),
                              (a[2] * b[0]) - (b[2] * a[0]),
                              (a[0] * b[1]) - (b[0] * a[1]));
}

// arguments are of type basic_vector, and there is a compact swizzled implementation
template <dsga::floating_point_scalar T1, dsga::floating_point_scalar T2>
[[nodiscard]] constexpr auto cross(const dsga::basic_vector<T1, 3> &a,
                                   const dsga::basic_vector<T2, 3> &b) noexcept
{
    return (a.yzx * b.zxy) - (a.zxy * b.yzx);
}

Relevant GLSL Overview

Our programming environment is C++20, not a GLSL shader program, so the entire GLSL Shading language specification is a super-set of what we are trying to achieve. We really just want the vector and matrix data structures (and their corresponding functions and behavior) to be usable in a C++20 environment. Another term for this type of programming is array programming.

The following links to the shading specification should help with understanding what we are trying to implement with this header-only library.

Variables and Types
- Basic Types: we have added support for vectors to also hold values of type std::size_t, unsigned long long (which is what std::size_t really is for x64), and signed long long.
- Vectors: GLSL does not have 1-dimensional vectors, but we do, which we have using directives to give them names that describe them as scalars and not as vectors, e.g., dsga::iscal, dsga::dscal, dsga::bscal. We support 1-dimensional vectors because GLSL does something special with the fundamental types, allowing them to be swizzled. We use the 1-dimensional vectors to mimic that ability.
```
// glsl
double value = 10.0;
dvec3 swizzled_value = value.xxx; 

// dsga
// dscal is an alias for dsga::basic_vector<double, 1>
dsga::dscal value = 10.0;
dsga::dvec3 swizzled_value = value.xxx; 
```
  1-dimensional vectors types are also the return type for single component swizzles, e.g., val.x, val.y, val.z, val.w. They are designed to be easily convertible to the underlying type of the vector elements.
- Matrices
Operators and Expressions
- Vector and Matrix Constructors
- Vector and Scalar Components and Length: we only allow swizzling with the { x, y, z, w } component names. Support for { r, g, b, a } and { s, t, p , q } has not been implemented.
  
  In addition, you cannot swizzle a swizzle. I am currently unclear if this is a constraint of the specification, but it is a constraint of dsga's implementation:
```
auto my_vec = dsga::vec3(10, 20, 30);
auto double_swiz = my_vec.zxy.x;           // error: no such data member x
auto swiz = my_vec.zxy;                    // swizzle type is not dsga::vec3
auto swiz_again = swiz.x;                  // error: no such data member x
auto try_swiz_again = dsga::vec3(swiz).x;  // wrapping with dsga::vec3 works
dsga::vec3 swiz_reborn = my_vec.zxy;       // dsga::vec3 constructor from swizzle
auto and_swiz_again = swiz_reborn.x;       // works
```
- Matrix Components
- Assignments
- Expressions
- Vector and Matrix Operations
- Out-of-Bounds Accesses: we have asserts for operator[] vectors and matrices, which help with debug runtimes and constexpr variable validity. These asserts (and others in the library) can be disabled by adding the following prior to including the header dsga.hxx:
```
#define DSGA_DISABLE_ASSERTS
#include "dsga.hxx"
```
Built-In Functions: we support the additional types std::size_t, unsigned long long, and signed long long in the functions where appropriate. We also added bit conversion functions between these 64-bit integral types and double.

We also support using double for all the functions where float is supported, with the exception of the bit conversion functions for float with 32-bit integral types, and double with the 64-bit integral types.
- Angle and Trigonometry Functions: there are also scalar versions of these functions, but where c++ does the same thing, it might be easier to use the std:: version instead of the dsga:: version.
- Exponential Functions: there are also scalar versions of these functions, but where c++ does the same thing, it might be easier to use the std:: version instead of the dsga:: version.
- Common Functions: there are also scalar versions of these functions, but where c++ does the same thing, it might be easier to use the std:: version instead of the dsga:: version.
- Geometric Functions: ftransform() is not implemented as it is only for GLSL vertex shader programs.
- Matrix Functions
- Vector Relational Functions: GLSL has a vector function not(), but not is a c++ keyword. Instead of naming this function not(), we name it logicalNot().
  
  In addition, we have added the non-GLSL convenience function none(), which returns !any().

Implemented Interfaces

To make the vectors and matrices as useful as possible in a C++ context, various C++ customization points were implemented or interfaces partially emulated, e.g., std::valarray<>. There are many options for data access. For dsga vectors and matrices, we have:

Swizzle access like GLSL (vector only)
- Only from the set of { x, y, z, w }, e.g., foo.wyxz
std::tuple protocol, structured bindings
- get
- tuple_size
- tuple_element
Iterator access, ranges, range-for loop
- begin
- cbegin
- rbegin
- crbegin
- end
- cend
- rend
- crend
Index access (logical)
- operator []
- size
- length
Pointer access (physical), std::span (for contiguous range types dsga::basic_vector and dsga::basic_matrix)
- data
  - vector - pointer to scalars of concept type dsga::dimensional_scalar
  - matrix - pointer to column vectors whose scalars are of concept type dsga::floating_point_scalar
- size
- vector only - these ordering facilities allow logical use of data
  - offsets
  - sequence
Type Conversions
- to_vector - from both std::array and C style arrays
- to_matrix - from both std::array and C style arrays
- to_array - from both dsga::basic_matrix and dsga::vector_base to std::array
- std::span example
Text output
- std::ostream example
- std::formatter example
std::valarray API (vector only)
- apply
- query - not in std::valarray nor GLSL - like apply() but for boolean predicates
- shift
- cshift
- min
- max
- sum

Installation

This is a single header library, where you just need the file dsga.hxx. Things are defined in the dsga namespace. The types provided by this library can be seen summarized in the documentation, using directives.

Under the hood, we depend on the cxcm project for constexpr versions of some cmath functions. cxcm has been brought into dsga.hxx, converted to a nested namespace cxcm under namespace dsga, so we don't need to also include the files from cxcm.

There are asserts in the codebase that can be disabled by defining the macro DSGA_DISABLE_ASSERTS.

This may be a single header library, but if Visual Studio is being used, we recommend to also get the dsga.natvis file for debugging and inspecting vectors and matrices in the IDE. While debugging this on Linux (WSL2: Windows Subsystem for Linux) with gcc in Visual Studio Code, we created a .natvis file for that too.

This is a c++20 library, so that needs to be the minimum standard that you tell the compiler to use.

Status

Current version: v2.2.2

Everything major has some tests, but code coverage is not 100%.
Last Released: v2.0.0
Change Log

The next steps

Refining API documentation.
Working on better cmake support.
Add more tests.

Usage

Use it more or less like you would use vectors and matrices in a shader program, but not necessarily for shading. We hope to be able to use it for rapid development of geometric algorithms. See the examples directory.

The documentation explains more about how the vector and matrix classes work, and describes the API.

More in depth explanation can be found in the details.

Testing

This project uses doctest for testing. We occasionally use nanobench for understanding implementation tradeoffs.

All tests are currently 100% PASSING on all the testing platforms and compilers.

The tests have been most recently run on:

Windows 11 Native

MSVC 2022 - v17.12.3

[doctest] doctest version is "2.4.11"
[doctest] run with "--help" for options
===============================================================================
[doctest] test cases:  110 |  110 passed | 0 failed | 0 skipped
[doctest] assertions: 2165 | 2165 passed | 0 failed |
[doctest] Status: SUCCESS!

gcc 14.2.0 on Windows, MSYS2 distribution:

[doctest] doctest version is "2.4.11"
[doctest] run with "--help" for options
===============================================================================
[doctest] test cases:  110 |  110 passed | 0 failed | 0 skipped
[doctest] assertions: 2165 | 2165 passed | 0 failed |
[doctest] Status: SUCCESS!

clang 19.1.5 on Windows, official binaries: