Skip to content

Instantly share code, notes, and snippets.

@sad1e

sad1e/fmath.h

Last active September 3, 2024 08:56
Show Gist options
  • Save sad1e/a445e464f0d72b6c28c05a6f73e315f3 to your computer and use it in GitHub Desktop.
Save sad1e/a445e464f0d72b6c28c05a6f73e315f3 to your computer and use it in GitHub Desktop.
Download ZIP
parula colormap generator.
Raw
fmath.h
#include <cmath>
#include <concepts>
#define Float double
namespace fmath {
constexpr Float pi = 3.141592653589793;
template <std::floating_point T>
constexpr T clamp(T v, T mi, T ma) {
return v <= mi ? mi : v >= ma ? ma : v;
}
template <std::floating_point T>
constexpr T saturate(T v) {
return clamp(v, 0.0, 1.0);
}
template <std::floating_point T>
constexpr T smoothstep(T a, T b, T x) {
T t = saturate((x - a) / (b - a));
return t * t * (3.0 - (2.0 * t));
}
template <std::floating_point T>
constexpr T mix(T x, T y, T a) {
return x * (1 - a) + y * a;
}
struct vec2 {
vec2() = default;
vec2(const vec2& v) = default;
vec2& operator=(const vec2& v) = default;
vec2(vec2&& v) = default;
vec2& operator=(vec2&& v) = default;
explicit vec2(Float v) : x{v}, y{v} {}
vec2(Float x_, Float y_) : x{x_}, y{y_} {}
vec2 operator+(const vec2& v) const { return vec2{x + v.x, y + v.y}; }
vec2& operator+=(const vec2& v) {
x += v.x;
y += v.y;
return *this;
}
vec2 operator+(Float c) const { return vec2{x + c, y + c}; }
vec2& operator+=(Float c) {
x += c;
y += c;
return *this;
}
vec2 operator-(const vec2& v) const { return vec2{x - v.x, y - v.y}; }
vec2& operator-=(const vec2& v) {
x -= v.x;
y -= v.y;
return *this;
}
vec2 operator-(Float c) const { return vec2{x - c, y - c}; }
vec2& operator-=(Float c) {
x -= c;
y -= c;
return *this;
}
vec2 operator*(const vec2& v) const { return vec2{x * v.x, y * v.y}; }
vec2& operator*=(const vec2& v) {
x *= v.x;
y *= v.y;
return *this;
}
vec2 operator*(Float c) const { return vec2{x * c, y * c}; }
vec2& operator*=(Float c) {
x *= c;
y *= c;
return *this;
}
vec2 operator/(const vec2& v) const { return vec2{x / v.x, y / v.y}; }
vec2& operator/=(const vec2& v) {
x /= v.x;
y /= v.y;
return *this;
}
vec2 operator/(Float c) const { return vec2{x / c, y / c}; }
vec2& operator/=(Float c) {
x /= c;
y /= c;
return *this;
}
vec2 operator-() { return vec2{-x, -y}; }
friend vec2 operator+(Float c, vec2 v);
friend vec2 operator-(Float c, vec2 v);
friend vec2 operator*(Float c, vec2 v);
friend vec2 operator/(Float c, vec2 v);
Float x{0}, y{0};
};
struct vec3 {
vec3() = default;
vec3(const vec3& v) = default;
vec3& operator=(const vec3& v) = default;
vec3(vec3&& v) = default;
vec3& operator=(vec3&& v) = default;
explicit vec3(Float v) : x{v}, y{v}, z{v} {}
vec3(Float x_, Float y_, Float z_) : x{x_}, y{y_}, z{z_} {}
vec3 operator+(const vec3& v) const {
return vec3{x + v.x, y + v.y, z + v.z};
}
vec3& operator+=(const vec3& v) {
x += v.x;
y += v.y;
z += v.z;
return *this;
}
vec3 operator+(Float c) const { return vec3{x + c, y + c, z + c}; }
vec3& operator+=(Float c) {
x += c;
y += c;
z += c;
return *this;
}
vec3 operator-(const vec3& v) const {
return vec3{x - v.x, y - v.y, z - v.z};
}
vec3& operator-=(const vec3& v) {
x -= v.x;
y -= v.y;
z -= v.z;
return *this;
}
vec3 operator-(Float c) const { return vec3{x - c, y - c, z - c}; }
vec3& operator-=(Float c) {
x -= c;
y -= c;
z -= c;
return *this;
}
vec3 operator*(const vec3& v) const {
return vec3{x * v.x, y * v.y, z * v.z};
}
vec3& operator*=(const vec3& v) {
x *= v.x;
y *= v.y;
z *= v.z;
return *this;
}
vec3 operator*(Float c) const { return vec3{x * c, y * c, z * c}; }
vec3& operator*=(Float c) {
x *= c;
y *= c;
z *= c;
return *this;
}
vec3 operator/(const vec3& v) const {
return vec3{x / v.x, y / v.y, z / v.z};
}
vec3& operator/=(const vec3& v) {
x /= v.x;
y /= v.y;
z /= v.z;
return *this;
}
vec3 operator/(Float c) const { return vec3{x / c, y / c, z / c}; }
vec3& operator/=(Float c) {
x /= c;
y /= c;
z /= c;
return *this;
}
vec3 operator-() { return vec3{-x, -y, -z}; }
friend vec3 operator+(Float c, vec3 v);
friend vec3 operator-(Float c, vec3 v);
friend vec3 operator*(Float c, vec3 v);
friend vec3 operator/(Float c, vec3 v);
Float x{0.0}, y{0.0}, z{0.0};
};
vec2 operator+(Float c, vec2 v) { return vec2{c + v.x, c + v.y}; }
vec2 operator-(Float c, vec2 v) { return vec2{c - v.x, c - v.y}; }
vec2 operator*(Float c, vec2 v) { return vec2{c * v.x, c * v.y}; }
vec2 operator/(Float c, vec2 v) { return vec2{c / v.x, c / v.y}; }
vec3 operator+(Float c, vec3 v) { return vec3{c + v.x, c + v.y, c + v.z}; }
vec3 operator-(Float c, vec3 v) { return vec3{c - v.x, c - v.y, c - v.z}; }
vec3 operator*(Float c, vec3 v) { return vec3{c * v.x, c * v.y, c * v.z}; }
vec3 operator/(Float c, vec3 v) { return vec3{c / v.x, c / v.y, c / v.z}; }
Float dot(const vec3& v0, const vec3& v1) {
return v0.x * v1.x + v0.y * v1.y + v0.z * v1.z;
}
Float luminance(const vec3& v) {
return v.x * 0.299 + v.y * 0.587 + v.z * 0.114;
}
Float pow2(Float x) { return x * x; }
Float pow3(Float x) { return x * x * x; }
Float pow1_4(Float x) { return std::sqrt(std::sqrt(x)); }
Float pow4(Float x) { return x * x * x * x; }
Float pow5(Float x) { return x * x * x * x * x; }
vec3 pow(const vec3& b, const vec3& e) {
return vec3{std::pow(b.x, e.x), std::pow(b.y, e.y), std::pow(b.z, e.z)};
}
vec3 sin(const vec3& v) {
return vec3{std::sin(v.x), std::sin(v.y), std::sin(v.z)};
}
vec3 cos(const vec3& v) {
return vec3{std::cos(v.x), std::cos(v.y), std::cos(v.z)};
}
vec3 mix(const vec3& x, const vec3& y, const vec3& a) {
return vec3{
mix(x.x, y.x, a.x),
mix(x.y, y.y, a.y),
mix(x.z, y.z, a.z),
};
}
vec3 cross(const vec3& v0, const vec3& v1) {
vec3 v0_yzx = vec3(v0.y, v0.z, v0.x);
vec3 v1_yzx = vec3(v1.y, v1.z, v1.x);
vec3 v2 = v0 * v1_yzx - v0_yzx * v1;
return vec3{v2.y, v2.z, v2.x};
}
vec3 sqrt(const vec3& v) {
return vec3{std::sqrt(v.x), std::sqrt(v.y), std::sqrt(v.z)};
}
vec3 rsqrt(const vec3& v) { return 1.0 / fmath::sqrt(v); }
vec3 normalize(const vec3& v) { return rsqrt(vec3(dot(v, v))) * v; }
Float rcp(Float x) { return (Float)1.0 / x; }
} // namespace fmath
Raw
generate_parula.cpp
namespace parula {
Float cheaplobe(Float x, Float m, Float w) {
return fmath::smoothstep(m - w, m, x) * fmath::smoothstep(m + w, m, x);
}
fmath::vec3 heatmap(Float c) {
return fmath::pow(
0.5 * (1.0 + fmath::cos(3.1 * c - fmath::vec3(2.7, 1.3, -0.2))),
fmath::vec3(3., 3., 4.));
}
fmath::vec3 heatmap(fmath::vec3 c) {
return fmath::pow(
0.5 * (1.0 + fmath::cos(3.1 * c - fmath::vec3(2.7, 1.3, -0.2))),
fmath::vec3(3., 3., 4.));
}
fmath::vec3 parula(Float c) {
fmath::vec3 d = fmath::vec3(
-0.1 + 0.22 * cheaplobe(c, 0.83, 0.19) - 0.07 * cheaplobe(c, 0.55, 0.16),
0.02 - 0.21 * cheaplobe(c, 0.83, 0.26) -
0.23 * cheaplobe(c, -0.02, 0.14) - 0.03 * cheaplobe(c, 0.59, 0.12),
-0.11 * cheaplobe(c, 0.57, 0.2));
fmath::vec3 w =
fmath::vec3(0.52, 0.5, 0.53) +
fmath::vec3(0.52, 0.52, 0.45) *
fmath::sin(fmath::vec3(2.80, 2.80, 2.04) *
(fmath::vec3(c) - fmath::vec3(0.4, 0.4, -0.63)));
fmath::vec3 p = fmath::vec3(fmath::mix(0.7, 0.8, c), fmath::mix(0.3, 0.7, c),
fmath::mix(2.9, 2.2, c));
return fmath::pow(d + fmath::pow(w, p), fmath::vec3(1.5));
}
fmath::vec3 parula(fmath::vec3 c) {
Float l = luminance(c);
fmath::vec3 d = fmath::vec3(
-0.1 + 0.22 * cheaplobe(l, 0.83, 0.19) - 0.07 * cheaplobe(l, 0.55, 0.16),
0.02 - 0.21 * cheaplobe(l, 0.83, 0.26) -
0.23 * cheaplobe(l, -0.02, 0.14) - 0.03 * cheaplobe(l, 0.59, 0.12),
-0.11 * cheaplobe(l, 0.57, 0.2));
fmath::vec3 w = fmath::vec3(0.52, 0.5, 0.53) +
fmath::vec3(0.52, 0.52, 0.45) *
fmath::sin(fmath::vec3(2.80, 2.80, 2.04) *
(c - fmath::vec3(0.4, 0.4, -0.63)));
fmath::vec3 p =
fmath::vec3(fmath::mix(0.7, 0.8, c.x), fmath::mix(0.3, 0.7, c.y),
fmath::mix(2.9, 2.2, c.z));
return fmath::pow(d + fmath::pow(w, p), fmath::vec3(1.5));
}
} // namespace parula
void GenerateParulaColormap() {
int width = 1024;
int height = 32;
pngwriter png(width, height, 0, "Parula.png");
for (int i = 1; i <= width; ++i) {
double uv_x = (double)i / width;
fmath::vec3 v = parula::parula(uv_x);
for (int j = 1; j <= height; ++j) {
png.plot(i, j, v.x, v.y, v.z);
}
}
png.close();
}
Raw
parula.md
@sad1e
Copy link
Author

sad1e commented Aug 17, 2024

Parula

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment