diff options
Diffstat (limited to 'thirdparty/misc')
| -rw-r--r-- | thirdparty/misc/ok_color.h | 688 | ||||
| -rw-r--r-- | thirdparty/misc/ok_color_shader.h | 663 | ||||
| -rw-r--r-- | thirdparty/misc/open-simplex-noise-LICENSE | 25 | ||||
| -rw-r--r-- | thirdparty/misc/open-simplex-noise-no-allocate.patch | 133 | ||||
| -rw-r--r-- | thirdparty/misc/open-simplex-noise.c | 2255 | ||||
| -rw-r--r-- | thirdparty/misc/open-simplex-noise.h | 58 | ||||
| -rw-r--r-- | thirdparty/misc/patches/polypartition-godot-types.patch | 87 | ||||
| -rw-r--r-- | thirdparty/misc/polypartition.cpp | 16 | ||||
| -rw-r--r-- | thirdparty/misc/polypartition.h | 6 | ||||
| -rw-r--r-- | thirdparty/misc/stb_rect_pack.h | 35 |
10 files changed, 1420 insertions, 2546 deletions
diff --git a/thirdparty/misc/ok_color.h b/thirdparty/misc/ok_color.h new file mode 100644 index 0000000000..dbc7dafc36 --- /dev/null +++ b/thirdparty/misc/ok_color.h @@ -0,0 +1,688 @@ +// Copyright(c) 2021 Björn Ottosson +// +// 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. + +#ifndef OK_COLOR_H +#define OK_COLOR_H + +#include <cmath> +#include <cfloat> + +class ok_color +{ +public: + +struct Lab { float L; float a; float b; }; +struct RGB { float r; float g; float b; }; +struct HSV { float h; float s; float v; }; +struct HSL { float h; float s; float l; }; +struct LC { float L; float C; }; + +// Alternative representation of (L_cusp, C_cusp) +// Encoded so S = C_cusp/L_cusp and T = C_cusp/(1-L_cusp) +// The maximum value for C in the triangle is then found as fmin(S*L, T*(1-L)), for a given L +struct ST { float S; float T; }; + +static constexpr float pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062f; + +float clamp(float x, float min, float max) +{ + if (x < min) + return min; + if (x > max) + return max; + + return x; +} + +float sgn(float x) +{ + return (float)(0.f < x) - (float)(x < 0.f); +} + +float srgb_transfer_function(float a) +{ + return .0031308f >= a ? 12.92f * a : 1.055f * powf(a, .4166666666666667f) - .055f; +} + +float srgb_transfer_function_inv(float a) +{ + return .04045f < a ? powf((a + .055f) / 1.055f, 2.4f) : a / 12.92f; +} + +Lab linear_srgb_to_oklab(RGB c) +{ + float l = 0.4122214708f * c.r + 0.5363325363f * c.g + 0.0514459929f * c.b; + float m = 0.2119034982f * c.r + 0.6806995451f * c.g + 0.1073969566f * c.b; + float s = 0.0883024619f * c.r + 0.2817188376f * c.g + 0.6299787005f * c.b; + + float l_ = cbrtf(l); + float m_ = cbrtf(m); + float s_ = cbrtf(s); + + return { + 0.2104542553f * l_ + 0.7936177850f * m_ - 0.0040720468f * s_, + 1.9779984951f * l_ - 2.4285922050f * m_ + 0.4505937099f * s_, + 0.0259040371f * l_ + 0.7827717662f * m_ - 0.8086757660f * s_, + }; +} + +RGB oklab_to_linear_srgb(Lab c) +{ + float l_ = c.L + 0.3963377774f * c.a + 0.2158037573f * c.b; + float m_ = c.L - 0.1055613458f * c.a - 0.0638541728f * c.b; + float s_ = c.L - 0.0894841775f * c.a - 1.2914855480f * c.b; + + float l = l_ * l_ * l_; + float m = m_ * m_ * m_; + float s = s_ * s_ * s_; + + return { + +4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s, + -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s, + -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s, + }; +} + +// Finds the maximum saturation possible for a given hue that fits in sRGB +// Saturation here is defined as S = C/L +// a and b must be normalized so a^2 + b^2 == 1 +float compute_max_saturation(float a, float b) +{ + // Max saturation will be when one of r, g or b goes below zero. + + // Select different coefficients depending on which component goes below zero first + float k0, k1, k2, k3, k4, wl, wm, ws; + + if (-1.88170328f * a - 0.80936493f * b > 1) + { + // Red component + k0 = +1.19086277f; k1 = +1.76576728f; k2 = +0.59662641f; k3 = +0.75515197f; k4 = +0.56771245f; + wl = +4.0767416621f; wm = -3.3077115913f; ws = +0.2309699292f; + } + else if (1.81444104f * a - 1.19445276f * b > 1) + { + // Green component + k0 = +0.73956515f; k1 = -0.45954404f; k2 = +0.08285427f; k3 = +0.12541070f; k4 = +0.14503204f; + wl = -1.2684380046f; wm = +2.6097574011f; ws = -0.3413193965f; + } + else + { + // Blue component + k0 = +1.35733652f; k1 = -0.00915799f; k2 = -1.15130210f; k3 = -0.50559606f; k4 = +0.00692167f; + wl = -0.0041960863f; wm = -0.7034186147f; ws = +1.7076147010f; + } + + // Approximate max saturation using a polynomial: + float S = k0 + k1 * a + k2 * b + k3 * a * a + k4 * a * b; + + // Do one step Halley's method to get closer + // this gives an error less than 10e6, except for some blue hues where the dS/dh is close to infinite + // this should be sufficient for most applications, otherwise do two/three steps + + float k_l = +0.3963377774f * a + 0.2158037573f * b; + float k_m = -0.1055613458f * a - 0.0638541728f * b; + float k_s = -0.0894841775f * a - 1.2914855480f * b; + + { + float l_ = 1.f + S * k_l; + float m_ = 1.f + S * k_m; + float s_ = 1.f + S * k_s; + + float l = l_ * l_ * l_; + float m = m_ * m_ * m_; + float s = s_ * s_ * s_; + + float l_dS = 3.f * k_l * l_ * l_; + float m_dS = 3.f * k_m * m_ * m_; + float s_dS = 3.f * k_s * s_ * s_; + + float l_dS2 = 6.f * k_l * k_l * l_; + float m_dS2 = 6.f * k_m * k_m * m_; + float s_dS2 = 6.f * k_s * k_s * s_; + + float f = wl * l + wm * m + ws * s; + float f1 = wl * l_dS + wm * m_dS + ws * s_dS; + float f2 = wl * l_dS2 + wm * m_dS2 + ws * s_dS2; + + S = S - f * f1 / (f1 * f1 - 0.5f * f * f2); + } + + return S; +} + +// finds L_cusp and C_cusp for a given hue +// a and b must be normalized so a^2 + b^2 == 1 +LC find_cusp(float a, float b) +{ + // First, find the maximum saturation (saturation S = C/L) + float S_cusp = compute_max_saturation(a, b); + + // Convert to linear sRGB to find the first point where at least one of r,g or b >= 1: + RGB rgb_at_max = oklab_to_linear_srgb({ 1, S_cusp * a, S_cusp * b }); + float L_cusp = cbrtf(1.f / fmax(fmax(rgb_at_max.r, rgb_at_max.g), rgb_at_max.b)); + float C_cusp = L_cusp * S_cusp; + + return { L_cusp , C_cusp }; +} + +// Finds intersection of the line defined by +// L = L0 * (1 - t) + t * L1; +// C = t * C1; +// a and b must be normalized so a^2 + b^2 == 1 +float find_gamut_intersection(float a, float b, float L1, float C1, float L0, LC cusp) +{ + // Find the intersection for upper and lower half seprately + float t; + if (((L1 - L0) * cusp.C - (cusp.L - L0) * C1) <= 0.f) + { + // Lower half + + t = cusp.C * L0 / (C1 * cusp.L + cusp.C * (L0 - L1)); + } + else + { + // Upper half + + // First intersect with triangle + t = cusp.C * (L0 - 1.f) / (C1 * (cusp.L - 1.f) + cusp.C * (L0 - L1)); + + // Then one step Halley's method + { + float dL = L1 - L0; + float dC = C1; + + float k_l = +0.3963377774f * a + 0.2158037573f * b; + float k_m = -0.1055613458f * a - 0.0638541728f * b; + float k_s = -0.0894841775f * a - 1.2914855480f * b; + + float l_dt = dL + dC * k_l; + float m_dt = dL + dC * k_m; + float s_dt = dL + dC * k_s; + + + // If higher accuracy is required, 2 or 3 iterations of the following block can be used: + { + float L = L0 * (1.f - t) + t * L1; + float C = t * C1; + + float l_ = L + C * k_l; + float m_ = L + C * k_m; + float s_ = L + C * k_s; + + float l = l_ * l_ * l_; + float m = m_ * m_ * m_; + float s = s_ * s_ * s_; + + float ldt = 3 * l_dt * l_ * l_; + float mdt = 3 * m_dt * m_ * m_; + float sdt = 3 * s_dt * s_ * s_; + + float ldt2 = 6 * l_dt * l_dt * l_; + float mdt2 = 6 * m_dt * m_dt * m_; + float sdt2 = 6 * s_dt * s_dt * s_; + + float r = 4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s - 1; + float r1 = 4.0767416621f * ldt - 3.3077115913f * mdt + 0.2309699292f * sdt; + float r2 = 4.0767416621f * ldt2 - 3.3077115913f * mdt2 + 0.2309699292f * sdt2; + + float u_r = r1 / (r1 * r1 - 0.5f * r * r2); + float t_r = -r * u_r; + + float g = -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s - 1; + float g1 = -1.2684380046f * ldt + 2.6097574011f * mdt - 0.3413193965f * sdt; + float g2 = -1.2684380046f * ldt2 + 2.6097574011f * mdt2 - 0.3413193965f * sdt2; + + float u_g = g1 / (g1 * g1 - 0.5f * g * g2); + float t_g = -g * u_g; + + b = -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s - 1; + float b1 = -0.0041960863f * ldt - 0.7034186147f * mdt + 1.7076147010f * sdt; + float b2 = -0.0041960863f * ldt2 - 0.7034186147f * mdt2 + 1.7076147010f * sdt2; + + float u_b = b1 / (b1 * b1 - 0.5f * b * b2); + float t_b = -b * u_b; + + t_r = u_r >= 0.f ? t_r : FLT_MAX; + t_g = u_g >= 0.f ? t_g : FLT_MAX; + t_b = u_b >= 0.f ? t_b : FLT_MAX; + + t += fmin(t_r, fmin(t_g, t_b)); + } + } + } + + return t; +} + +float find_gamut_intersection(float a, float b, float L1, float C1, float L0) +{ + // Find the cusp of the gamut triangle + LC cusp = find_cusp(a, b); + + return find_gamut_intersection(a, b, L1, C1, L0, cusp); +} + +RGB gamut_clip_preserve_chroma(RGB rgb) +{ + if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0) + return rgb; + + Lab lab = linear_srgb_to_oklab(rgb); + + float L = lab.L; + float eps = 0.00001f; + float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b)); + float a_ = lab.a / C; + float b_ = lab.b / C; + + float L0 = clamp(L, 0, 1); + + float t = find_gamut_intersection(a_, b_, L, C, L0); + float L_clipped = L0 * (1 - t) + t * L; + float C_clipped = t * C; + + return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ }); +} + +RGB gamut_clip_project_to_0_5(RGB rgb) +{ + if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0) + return rgb; + + Lab lab = linear_srgb_to_oklab(rgb); + + float L = lab.L; + float eps = 0.00001f; + float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b)); + float a_ = lab.a / C; + float b_ = lab.b / C; + + float L0 = 0.5; + + float t = find_gamut_intersection(a_, b_, L, C, L0); + float L_clipped = L0 * (1 - t) + t * L; + float C_clipped = t * C; + + return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ }); +} + +RGB gamut_clip_project_to_L_cusp(RGB rgb) +{ + if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0) + return rgb; + + Lab lab = linear_srgb_to_oklab(rgb); + + float L = lab.L; + float eps = 0.00001f; + float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b)); + float a_ = lab.a / C; + float b_ = lab.b / C; + + // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once. + LC cusp = find_cusp(a_, b_); + + float L0 = cusp.L; + + float t = find_gamut_intersection(a_, b_, L, C, L0); + + float L_clipped = L0 * (1 - t) + t * L; + float C_clipped = t * C; + + return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ }); +} + +RGB gamut_clip_adaptive_L0_0_5(RGB rgb, float alpha = 0.05f) +{ + if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0) + return rgb; + + Lab lab = linear_srgb_to_oklab(rgb); + + float L = lab.L; + float eps = 0.00001f; + float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b)); + float a_ = lab.a / C; + float b_ = lab.b / C; + + float Ld = L - 0.5f; + float e1 = 0.5f + fabs(Ld) + alpha * C; + float L0 = 0.5f * (1.f + sgn(Ld) * (e1 - sqrtf(e1 * e1 - 2.f * fabs(Ld)))); + + float t = find_gamut_intersection(a_, b_, L, C, L0); + float L_clipped = L0 * (1.f - t) + t * L; + float C_clipped = t * C; + + return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ }); +} + +RGB gamut_clip_adaptive_L0_L_cusp(RGB rgb, float alpha = 0.05f) +{ + if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0) + return rgb; + + Lab lab = linear_srgb_to_oklab(rgb); + + float L = lab.L; + float eps = 0.00001f; + float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b)); + float a_ = lab.a / C; + float b_ = lab.b / C; + + // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once. + LC cusp = find_cusp(a_, b_); + + float Ld = L - cusp.L; + float k = 2.f * (Ld > 0 ? 1.f - cusp.L : cusp.L); + + float e1 = 0.5f * k + fabs(Ld) + alpha * C / k; + float L0 = cusp.L + 0.5f * (sgn(Ld) * (e1 - sqrtf(e1 * e1 - 2.f * k * fabs(Ld)))); + + float t = find_gamut_intersection(a_, b_, L, C, L0); + float L_clipped = L0 * (1.f - t) + t * L; + float C_clipped = t * C; + + return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ }); +} + +float toe(float x) +{ + constexpr float k_1 = 0.206f; + constexpr float k_2 = 0.03f; + constexpr float k_3 = (1.f + k_1) / (1.f + k_2); + return 0.5f * (k_3 * x - k_1 + sqrtf((k_3 * x - k_1) * (k_3 * x - k_1) + 4 * k_2 * k_3 * x)); +} + +float toe_inv(float x) +{ + constexpr float k_1 = 0.206f; + constexpr float k_2 = 0.03f; + constexpr float k_3 = (1.f + k_1) / (1.f + k_2); + return (x * x + k_1 * x) / (k_3 * (x + k_2)); +} + +ST to_ST(LC cusp) +{ + float L = cusp.L; + float C = cusp.C; + return { C / L, C / (1 - L) }; +} + +// Returns a smooth approximation of the location of the cusp +// This polynomial was created by an optimization process +// It has been designed so that S_mid < S_max and T_mid < T_max +ST get_ST_mid(float a_, float b_) +{ + float S = 0.11516993f + 1.f / ( + +7.44778970f + 4.15901240f * b_ + + a_ * (-2.19557347f + 1.75198401f * b_ + + a_ * (-2.13704948f - 10.02301043f * b_ + + a_ * (-4.24894561f + 5.38770819f * b_ + 4.69891013f * a_ + ))) + ); + + float T = 0.11239642f + 1.f / ( + +1.61320320f - 0.68124379f * b_ + + a_ * (+0.40370612f + 0.90148123f * b_ + + a_ * (-0.27087943f + 0.61223990f * b_ + + a_ * (+0.00299215f - 0.45399568f * b_ - 0.14661872f * a_ + ))) + ); + + return { S, T }; +} + +struct Cs { float C_0; float C_mid; float C_max; }; +Cs get_Cs(float L, float a_, float b_) +{ + LC cusp = find_cusp(a_, b_); + + float C_max = find_gamut_intersection(a_, b_, L, 1, L, cusp); + ST ST_max = to_ST(cusp); + + // Scale factor to compensate for the curved part of gamut shape: + float k = C_max / fmin((L * ST_max.S), (1 - L) * ST_max.T); + + float C_mid; + { + ST ST_mid = get_ST_mid(a_, b_); + + // Use a soft minimum function, instead of a sharp triangle shape to get a smooth value for chroma. + float C_a = L * ST_mid.S; + float C_b = (1.f - L) * ST_mid.T; + C_mid = 0.9f * k * sqrtf(sqrtf(1.f / (1.f / (C_a * C_a * C_a * C_a) + 1.f / (C_b * C_b * C_b * C_b)))); + } + + float C_0; + { + // for C_0, the shape is independent of hue, so ST are constant. Values picked to roughly be the average values of ST. + float C_a = L * 0.4f; + float C_b = (1.f - L) * 0.8f; + + // Use a soft minimum function, instead of a sharp triangle shape to get a smooth value for chroma. + C_0 = sqrtf(1.f / (1.f / (C_a * C_a) + 1.f / (C_b * C_b))); + } + + return { C_0, C_mid, C_max }; +} + +RGB okhsl_to_srgb(HSL hsl) +{ + float h = hsl.h; + float s = hsl.s; + float l = hsl.l; + + if (l == 1.0f) + { + return { 1.f, 1.f, 1.f }; + } + + else if (l == 0.f) + { + return { 0.f, 0.f, 0.f }; + } + + float a_ = cosf(2.f * pi * h); + float b_ = sinf(2.f * pi * h); + float L = toe_inv(l); + + Cs cs = get_Cs(L, a_, b_); + float C_0 = cs.C_0; + float C_mid = cs.C_mid; + float C_max = cs.C_max; + + float mid = 0.8f; + float mid_inv = 1.25f; + + float C, t, k_0, k_1, k_2; + + if (s < mid) + { + t = mid_inv * s; + + k_1 = mid * C_0; + k_2 = (1.f - k_1 / C_mid); + + C = t * k_1 / (1.f - k_2 * t); + } + else + { + t = (s - mid)/ (1 - mid); + + k_0 = C_mid; + k_1 = (1.f - mid) * C_mid * C_mid * mid_inv * mid_inv / C_0; + k_2 = (1.f - (k_1) / (C_max - C_mid)); + + C = k_0 + t * k_1 / (1.f - k_2 * t); + } + + RGB rgb = oklab_to_linear_srgb({ L, C * a_, C * b_ }); + return { + srgb_transfer_function(rgb.r), + srgb_transfer_function(rgb.g), + srgb_transfer_function(rgb.b), + }; +} + +HSL srgb_to_okhsl(RGB rgb) +{ + Lab lab = linear_srgb_to_oklab({ + srgb_transfer_function_inv(rgb.r), + srgb_transfer_function_inv(rgb.g), + srgb_transfer_function_inv(rgb.b) + }); + + float C = sqrtf(lab.a * lab.a + lab.b * lab.b); + float a_ = lab.a / C; + float b_ = lab.b / C; + + float L = lab.L; + float h = 0.5f + 0.5f * atan2f(-lab.b, -lab.a) / pi; + + Cs cs = get_Cs(L, a_, b_); + float C_0 = cs.C_0; + float C_mid = cs.C_mid; + float C_max = cs.C_max; + + // Inverse of the interpolation in okhsl_to_srgb: + + float mid = 0.8f; + float mid_inv = 1.25f; + + float s; + if (C < C_mid) + { + float k_1 = mid * C_0; + float k_2 = (1.f - k_1 / C_mid); + + float t = C / (k_1 + k_2 * C); + s = t * mid; + } + else + { + float k_0 = C_mid; + float k_1 = (1.f - mid) * C_mid * C_mid * mid_inv * mid_inv / C_0; + float k_2 = (1.f - (k_1) / (C_max - C_mid)); + + float t = (C - k_0) / (k_1 + k_2 * (C - k_0)); + s = mid + (1.f - mid) * t; + } + + float l = toe(L); + return { h, s, l }; +} + + +RGB okhsv_to_srgb(HSV hsv) +{ + float h = hsv.h; + float s = hsv.s; + float v = hsv.v; + + float a_ = cosf(2.f * pi * h); + float b_ = sinf(2.f * pi * h); + + LC cusp = find_cusp(a_, b_); + ST ST_max = to_ST(cusp); + float S_max = ST_max.S; + float T_max = ST_max.T; + float S_0 = 0.5f; + float k = 1 - S_0 / S_max; + + // first we compute L and V as if the gamut is a perfect triangle: + + // L, C when v==1: + float L_v = 1 - s * S_0 / (S_0 + T_max - T_max * k * s); + float C_v = s * T_max * S_0 / (S_0 + T_max - T_max * k * s); + + float L = v * L_v; + float C = v * C_v; + + // then we compensate for both toe and the curved top part of the triangle: + float L_vt = toe_inv(L_v); + float C_vt = C_v * L_vt / L_v; + + float L_new = toe_inv(L); + C = C * L_new / L; + L = L_new; + + RGB rgb_scale = oklab_to_linear_srgb({ L_vt, a_ * C_vt, b_ * C_vt }); + float scale_L = cbrtf(1.f / fmax(fmax(rgb_scale.r, rgb_scale.g), fmax(rgb_scale.b, 0.f))); + + L = L * scale_L; + C = C * scale_L; + + RGB rgb = oklab_to_linear_srgb({ L, C * a_, C * b_ }); + return { + srgb_transfer_function(rgb.r), + srgb_transfer_function(rgb.g), + srgb_transfer_function(rgb.b), + }; +} + +HSV srgb_to_okhsv(RGB rgb) +{ + Lab lab = linear_srgb_to_oklab({ + srgb_transfer_function_inv(rgb.r), + srgb_transfer_function_inv(rgb.g), + srgb_transfer_function_inv(rgb.b) + }); + + float C = sqrtf(lab.a * lab.a + lab.b * lab.b); + float a_ = lab.a / C; + float b_ = lab.b / C; + + float L = lab.L; + float h = 0.5f + 0.5f * atan2f(-lab.b, -lab.a) / pi; + + LC cusp = find_cusp(a_, b_); + ST ST_max = to_ST(cusp); + float S_max = ST_max.S; + float T_max = ST_max.T; + float S_0 = 0.5f; + float k = 1 - S_0 / S_max; + + // first we find L_v, C_v, L_vt and C_vt + + float t = T_max / (C + L * T_max); + float L_v = t * L; + float C_v = t * C; + + float L_vt = toe_inv(L_v); + float C_vt = C_v * L_vt / L_v; + + // we can then use these to invert the step that compensates for the toe and the curved top part of the triangle: + RGB rgb_scale = oklab_to_linear_srgb({ L_vt, a_ * C_vt, b_ * C_vt }); + float scale_L = cbrtf(1.f / fmax(fmax(rgb_scale.r, rgb_scale.g), fmax(rgb_scale.b, 0.f))); + + L = L / scale_L; + C = C / scale_L; + + C = C * toe(L) / L; + L = toe(L); + + // we can now compute v and s: + + float v = L / L_v; + float s = (S_0 + T_max) * C_v / ((T_max * S_0) + T_max * k * C_v); + + return { h, s, v }; +} + +}; +#endif // OK_COLOR_H diff --git a/thirdparty/misc/ok_color_shader.h b/thirdparty/misc/ok_color_shader.h new file mode 100644 index 0000000000..40d83366ee --- /dev/null +++ b/thirdparty/misc/ok_color_shader.h @@ -0,0 +1,663 @@ +// Copyright(c) 2021 Björn Ottosson +// +// 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. + +#ifndef OK_COLOR_SHADER_H +#define OK_COLOR_SHADER_H + +#include "core/string/ustring.h" + +static String OK_COLOR_SHADER = R"(shader_type canvas_item; + +const float M_PI = 3.1415926535897932384626433832795; + +float cbrt( float x ) +{ + return sign(x)*pow(abs(x),1.0f/3.0f); +} + +float srgb_transfer_function(float a) +{ + return .0031308f >= a ? 12.92f * a : 1.055f * pow(a, .4166666666666667f) - .055f; +} + +float srgb_transfer_function_inv(float a) +{ + return .04045f < a ? pow((a + .055f) / 1.055f, 2.4f) : a / 12.92f; +} + +vec3 linear_srgb_to_oklab(vec3 c) +{ + float l = 0.4122214708f * c.r + 0.5363325363f * c.g + 0.0514459929f * c.b; + float m = 0.2119034982f * c.r + 0.6806995451f * c.g + 0.1073969566f * c.b; + float s = 0.0883024619f * c.r + 0.2817188376f * c.g + 0.6299787005f * c.b; + + float l_ = cbrt(l); + float m_ = cbrt(m); + float s_ = cbrt(s); + + return vec3( + 0.2104542553f * l_ + 0.7936177850f * m_ - 0.0040720468f * s_, + 1.9779984951f * l_ - 2.4285922050f * m_ + 0.4505937099f * s_, + 0.0259040371f * l_ + 0.7827717662f * m_ - 0.8086757660f * s_ + ); +} + +vec3 oklab_to_linear_srgb(vec3 c) +{ + float l_ = c.x + 0.3963377774f * c.y + 0.2158037573f * c.z; + float m_ = c.x - 0.1055613458f * c.y - 0.0638541728f * c.z; + float s_ = c.x - 0.0894841775f * c.y - 1.2914855480f * c.z; + + float l = l_ * l_ * l_; + float m = m_ * m_ * m_; + float s = s_ * s_ * s_; + + return vec3( + +4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s, + -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s, + -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s + ); +} + +// Finds the maximum saturation possible for a given hue that fits in sRGB +// Saturation here is defined as S = C/L +// a and b must be normalized so a^2 + b^2 == 1 +float compute_max_saturation(float a, float b) +{ + // Max saturation will be when one of r, g or b goes below zero. + + // Select different coefficients depending on which component goes below zero first + float k0, k1, k2, k3, k4, wl, wm, ws; + + if (-1.88170328f * a - 0.80936493f * b > 1.f) + { + // Red component + k0 = +1.19086277f; k1 = +1.76576728f; k2 = +0.59662641f; k3 = +0.75515197f; k4 = +0.56771245f; + wl = +4.0767416621f; wm = -3.3077115913f; ws = +0.2309699292f; + } + else if (1.81444104f * a - 1.19445276f * b > 1.f) + { + // Green component + k0 = +0.73956515f; k1 = -0.45954404f; k2 = +0.08285427f; k3 = +0.12541070f; k4 = +0.14503204f; + wl = -1.2684380046f; wm = +2.6097574011f; ws = -0.3413193965f; + } + else + { + // Blue component + k0 = +1.35733652f; k1 = -0.00915799f; k2 = -1.15130210f; k3 = -0.50559606f; k4 = +0.00692167f; + wl = -0.0041960863f; wm = -0.7034186147f; ws = +1.7076147010f; + } + + // Approximate max saturation using a polynomial: + float S = k0 + k1 * a + k2 * b + k3 * a * a + k4 * a * b; + + // Do one step Halley's method to get closer + // this gives an error less than 10e6, except for some blue hues where the dS/dh is close to infinite + // this should be sufficient for most applications, otherwise do two/three steps + + float k_l = +0.3963377774f * a + 0.2158037573f * b; + float k_m = -0.1055613458f * a - 0.0638541728f * b; + float k_s = -0.0894841775f * a - 1.2914855480f * b; + + { + float l_ = 1.f + S * k_l; + float m_ = 1.f + S * k_m; + float s_ = 1.f + S * k_s; + + float l = l_ * l_ * l_; + float m = m_ * m_ * m_; + float s = s_ * s_ * s_; + + float l_dS = 3.f * k_l * l_ * l_; + float m_dS = 3.f * k_m * m_ * m_; + float s_dS = 3.f * k_s * s_ * s_; + + float l_dS2 = 6.f * k_l * k_l * l_; + float m_dS2 = 6.f * k_m * k_m * m_; + float s_dS2 = 6.f * k_s * k_s * s_; + + float f = wl * l + wm * m + ws * s; + float f1 = wl * l_dS + wm * m_dS + ws * s_dS; + float f2 = wl * l_dS2 + wm * m_dS2 + ws * s_dS2; + + S = S - f * f1 / (f1 * f1 - 0.5f * f * f2); + } + + return S; +} + +// finds L_cusp and C_cusp for a given hue +// a and b must be normalized so a^2 + b^2 == 1 +vec2 find_cusp(float a, float b) +{ + // First, find the maximum saturation (saturation S = C/L) + float S_cusp = compute_max_saturation(a, b); + + // Convert to linear sRGB to find the first point where at least one of r,g or b >= 1: + vec3 rgb_at_max = oklab_to_linear_srgb(vec3( 1, S_cusp * a, S_cusp * b )); + float L_cusp = cbrt(1.f / max(max(rgb_at_max.r, rgb_at_max.g), rgb_at_max.b)); + float C_cusp = L_cusp * S_cusp; + + return vec2( L_cusp , C_cusp ); +} )" +R"(// Finds intersection of the line defined by +// L = L0 * (1 - t) + t * L1; +// C = t * C1; +// a and b must be normalized so a^2 + b^2 == 1 +float find_gamut_intersection(float a, float b, float L1, float C1, float L0, vec2 cusp) +{ + // Find the intersection for upper and lower half seprately + float t; + if (((L1 - L0) * cusp.y - (cusp.x - L0) * C1) <= 0.f) + { + // Lower half + + t = cusp.y * L0 / (C1 * cusp.x + cusp.y * (L0 - L1)); + } + else + { + // Upper half + + // First intersect with triangle + t = cusp.y * (L0 - 1.f) / (C1 * (cusp.x - 1.f) + cusp.y * (L0 - L1)); + + // Then one step Halley's method + { + float dL = L1 - L0; + float dC = C1; + + float k_l = +0.3963377774f * a + 0.2158037573f * b; + float k_m = -0.1055613458f * a - 0.0638541728f * b; + float k_s = -0.0894841775f * a - 1.2914855480f * b; + + float l_dt = dL + dC * k_l; + float m_dt = dL + dC * k_m; + float s_dt = dL + dC * k_s; + + + // If higher accuracy is required, 2 or 3 iterations of the following block can be used: + { + float L = L0 * (1.f - t) + t * L1; + float C = t * C1; + + float l_ = L + C * k_l; + float m_ = L + C * k_m; + float s_ = L + C * k_s; + + float l = l_ * l_ * l_; + float m = m_ * m_ * m_; + float s = s_ * s_ * s_; + + float ldt = 3.f * l_dt * l_ * l_; + float mdt = 3.f * m_dt * m_ * m_; + float sdt = 3.f * s_dt * s_ * s_; + + float ldt2 = 6.f * l_dt * l_dt * l_; + float mdt2 = 6.f * m_dt * m_dt * m_; + float sdt2 = 6.f * s_dt * s_dt * s_; + + float r = 4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s - 1.f; + float r1 = 4.0767416621f * ldt - 3.3077115913f * mdt + 0.2309699292f * sdt; + float r2 = 4.0767416621f * ldt2 - 3.3077115913f * mdt2 + 0.2309699292f * sdt2; + + float u_r = r1 / (r1 * r1 - 0.5f * r * r2); + float t_r = -r * u_r; + + float g = -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s - 1.f; + float g1 = -1.2684380046f * ldt + 2.6097574011f * mdt - 0.3413193965f * sdt; + float g2 = -1.2684380046f * ldt2 + 2.6097574011f * mdt2 - 0.3413193965f * sdt2; + + float u_g = g1 / (g1 * g1 - 0.5f * g * g2); + float t_g = -g * u_g; + + float b = -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s - 1.f; + float b1 = -0.0041960863f * ldt - 0.7034186147f * mdt + 1.7076147010f * sdt; + float b2 = -0.0041960863f * ldt2 - 0.7034186147f * mdt2 + 1.7076147010f * sdt2; + + float u_b = b1 / (b1 * b1 - 0.5f * b * b2); + float t_b = -b * u_b; + + t_r = u_r >= 0.f ? t_r : 10000.f; + t_g = u_g >= 0.f ? t_g : 10000.f; + t_b = u_b >= 0.f ? t_b : 10000.f; + + t += min(t_r, min(t_g, t_b)); + } + } + } + + return t; +} + +float find_gamut_intersection_5(float a, float b, float L1, float C1, float L0) +{ + // Find the cusp of the gamut triangle + vec2 cusp = find_cusp(a, b); + + return find_gamut_intersection(a, b, L1, C1, L0, cusp); +})" +R"( + +vec3 gamut_clip_preserve_chroma(vec3 rgb) +{ + if (rgb.r < 1.f && rgb.g < 1.f && rgb.b < 1.f && rgb.r > 0.f && rgb.g > 0.f && rgb.b > 0.f) + return rgb; + + vec3 lab = linear_srgb_to_oklab(rgb); + + float L = lab.x; + float eps = 0.00001f; + float C = max(eps, sqrt(lab.y * lab.y + lab.z * lab.z)); + float a_ = lab.y / C; + float b_ = lab.z / C; + + float L0 = clamp(L, 0.f, 1.f); + + float t = find_gamut_intersection_5(a_, b_, L, C, L0); + float L_clipped = L0 * (1.f - t) + t * L; + float C_clipped = t * C; + + return oklab_to_linear_srgb(vec3( L_clipped, C_clipped * a_, C_clipped * b_ )); +} + +vec3 gamut_clip_project_to_0_5(vec3 rgb) +{ + if (rgb.r < 1.f && rgb.g < 1.f && rgb.b < 1.f && rgb.r > 0.f && rgb.g > 0.f && rgb.b > 0.f) + return rgb; + + vec3 lab = linear_srgb_to_oklab(rgb); + + float L = lab.x; + float eps = 0.00001f; + float C = max(eps, sqrt(lab.y * lab.y + lab.z * lab.z)); + float a_ = lab.y / C; + float b_ = lab.z / C; + + float L0 = 0.5; + + float t = find_gamut_intersection_5(a_, b_, L, C, L0); + float L_clipped = L0 * (1.f - t) + t * L; + float C_clipped = t * C; + + return oklab_to_linear_srgb(vec3( L_clipped, C_clipped * a_, C_clipped * b_ )); +} + +vec3 gamut_clip_project_to_L_cusp(vec3 rgb) +{ + if (rgb.r < 1.f && rgb.g < 1.f && rgb.b < 1.f && rgb.r > 0.f && rgb.g > 0.f && rgb.b > 0.f) + return rgb; + + vec3 lab = linear_srgb_to_oklab(rgb); + + float L = lab.x; + float eps = 0.00001f; + float C = max(eps, sqrt(lab.y * lab.y + lab.z * lab.z)); + float a_ = lab.y / C; + float b_ = lab.z / C; + + // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once. + vec2 cusp = find_cusp(a_, b_); + + float L0 = cusp.x; + + float t = find_gamut_intersection_5(a_, b_, L, C, L0); + + float L_clipped = L0 * (1.f - t) + t * L; + float C_clipped = t * C; + + return oklab_to_linear_srgb(vec3( L_clipped, C_clipped * a_, C_clipped * b_ )); +} + +vec3 gamut_clip_adaptive_L0_0_5(vec3 rgb, float alpha) +{ + if (rgb.r < 1.f && rgb.g < 1.f && rgb.b < 1.f && rgb.r > 0.f && rgb.g > 0.f && rgb.b > 0.f) + return rgb; + + vec3 lab = linear_srgb_to_oklab(rgb); + + float L = lab.x; + float eps = 0.00001f; + float C = max(eps, sqrt(lab.y * lab.y + lab.z * lab.z)); + float a_ = lab.y / C; + float b_ = lab.z / C; + + float Ld = L - 0.5f; + float e1 = 0.5f + abs(Ld) + alpha * C; + float L0 = 0.5f * (1.f + sign(Ld) * (e1 - sqrt(e1 * e1 - 2.f * abs(Ld)))); + + float t = find_gamut_intersection_5(a_, b_, L, C, L0); + float L_clipped = L0 * (1.f - t) + t * L; + float C_clipped = t * C; + + return oklab_to_linear_srgb(vec3( L_clipped, C_clipped * a_, C_clipped * b_ )); +} + +vec3 gamut_clip_adaptive_L0_L_cusp(vec3 rgb, float alpha) +{ + if (rgb.r < 1.f && rgb.g < 1.f && rgb.b < 1.f && rgb.r > 0.f && rgb.g > 0.f && rgb.b > 0.f) + return rgb; + + vec3 lab = linear_srgb_to_oklab(rgb); + + float L = lab.x; + float eps = 0.00001f; + float C = max(eps, sqrt(lab.y * lab.y + lab.z * lab.z)); + float a_ = lab.y / C; + float b_ = lab.z / C; + + // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once. + vec2 cusp = find_cusp(a_, b_); + + float Ld = L - cusp.x; + float k = 2.f * (Ld > 0.f ? 1.f - cusp.x : cusp.x); + + float e1 = 0.5f * k + abs(Ld) + alpha * C / k; + float L0 = cusp.x + 0.5f * (sign(Ld) * (e1 - sqrt(e1 * e1 - 2.f * k * abs(Ld)))); + + float t = find_gamut_intersection_5(a_, b_, L, C, L0); + float L_clipped = L0 * (1.f - t) + t * L; + float C_clipped = t * C; + + return oklab_to_linear_srgb(vec3( L_clipped, C_clipped * a_, C_clipped * b_ )); +} + +float toe(float x) +{ + float k_1 = 0.206f; + float k_2 = 0.03f; + float k_3 = (1.f + k_1) / (1.f + k_2); + return 0.5f * (k_3 * x - k_1 + sqrt((k_3 * x - k_1) * (k_3 * x - k_1) + 4.f * k_2 * k_3 * x)); +} + +float toe_inv(float x) +{ + float k_1 = 0.206f; + float k_2 = 0.03f; + float k_3 = (1.f + k_1) / (1.f + k_2); + return (x * x + k_1 * x) / (k_3 * (x + k_2)); +} +)" +R"(vec2 to_ST(vec2 cusp) +{ + float L = cusp.x; + float C = cusp.y; + return vec2( C / L, C / (1.f - L) ); +} + +// Returns a smooth approximation of the location of the cusp +// This polynomial was created by an optimization process +// It has been designed so that S_mid < S_max and T_mid < T_max +vec2 get_ST_mid(float a_, float b_) +{ + float S = 0.11516993f + 1.f / ( + +7.44778970f + 4.15901240f * b_ + + a_ * (-2.19557347f + 1.75198401f * b_ + + a_ * (-2.13704948f - 10.02301043f * b_ + + a_ * (-4.24894561f + 5.38770819f * b_ + 4.69891013f * a_ + ))) + ); + + float T = 0.11239642f + 1.f / ( + +1.61320320f - 0.68124379f * b_ + + a_ * (+0.40370612f + 0.90148123f * b_ + + a_ * (-0.27087943f + 0.61223990f * b_ + + a_ * (+0.00299215f - 0.45399568f * b_ - 0.14661872f * a_ + ))) + ); + + return vec2( S, T ); +} + +vec3 get_Cs(float L, float a_, float b_) +{ + vec2 cusp = find_cusp(a_, b_); + + float C_max = find_gamut_intersection(a_, b_, L, 1.f, L, cusp); + vec2 ST_max = to_ST(cusp); + + // Scale factor to compensate for the curved part of gamut shape: + float k = C_max / min((L * ST_max.x), (1.f - L) * ST_max.y); + + float C_mid; + { + vec2 ST_mid = get_ST_mid(a_, b_); + + // Use a soft minimum function, instead of a sharp triangle shape to get a smooth value for chroma. + float C_a = L * ST_mid.x; + float C_b = (1.f - L) * ST_mid.y; + C_mid = 0.9f * k * sqrt(sqrt(1.f / (1.f / (C_a * C_a * C_a * C_a) + 1.f / (C_b * C_b * C_b * C_b)))); + } + + float C_0; + { + // for C_0, the shape is independent of hue, so vec2 are constant. Values picked to roughly be the average values of vec2. + float C_a = L * 0.4f; + float C_b = (1.f - L) * 0.8f; + + // Use a soft minimum function, instead of a sharp triangle shape to get a smooth value for chroma. + C_0 = sqrt(1.f / (1.f / (C_a * C_a) + 1.f / (C_b * C_b))); + } + + return vec3( C_0, C_mid, C_max ); +} + +vec3 okhsl_to_srgb(vec3 hsl) +{ + float h = hsl.x; + float s = hsl.y; + float l = hsl.z; + + if (l == 1.0f) + { + return vec3( 1.f, 1.f, 1.f ); + } + + else if (l == 0.f) + { + return vec3( 0.f, 0.f, 0.f ); + } + + float a_ = cos(2.f * M_PI * h); + float b_ = sin(2.f * M_PI * h); + float L = toe_inv(l); + + vec3 cs = get_Cs(L, a_, b_); + float C_0 = cs.x; + float C_mid = cs.y; + float C_max = cs.z; + + float mid = 0.8f; + float mid_inv = 1.25f; + + float C, t, k_0, k_1, k_2; + + if (s < mid) + { + t = mid_inv * s; + + k_1 = mid * C_0; + k_2 = (1.f - k_1 / C_mid); + + C = t * k_1 / (1.f - k_2 * t); + } + else + { + t = (s - mid)/ (1.f - mid); + + k_0 = C_mid; + k_1 = (1.f - mid) * C_mid * C_mid * mid_inv * mid_inv / C_0; + k_2 = (1.f - (k_1) / (C_max - C_mid)); + + C = k_0 + t * k_1 / (1.f - k_2 * t); + } + + vec3 rgb = oklab_to_linear_srgb(vec3( L, C * a_, C * b_ )); + return vec3( + srgb_transfer_function(rgb.r), + srgb_transfer_function(rgb.g), + srgb_transfer_function(rgb.b) + ); +} + +vec3 srgb_to_okhsl(vec3 rgb) +{ + vec3 lab = linear_srgb_to_oklab(vec3( + srgb_transfer_function_inv(rgb.r), + srgb_transfer_function_inv(rgb.g), + srgb_transfer_function_inv(rgb.b) + )); + + float C = sqrt(lab.y * lab.y + lab.z * lab.z); + float a_ = lab.y / C; + float b_ = lab.z / C; + + float L = lab.x; + float h = 0.5f + 0.5f * atan(-lab.z, -lab.y) / M_PI; + + vec3 cs = get_Cs(L, a_, b_); + float C_0 = cs.x; + float C_mid = cs.y; + float C_max = cs.z; + + // Inverse of the interpolation in okhsl_to_srgb: + + float mid = 0.8f; + float mid_inv = 1.25f; + + float s; + if (C < C_mid) + { + float k_1 = mid * C_0; + float k_2 = (1.f - k_1 / C_mid); + + float t = C / (k_1 + k_2 * C); + s = t * mid; + } + else + { + float k_0 = C_mid; + float k_1 = (1.f - mid) * C_mid * C_mid * mid_inv * mid_inv / C_0; + float k_2 = (1.f - (k_1) / (C_max - C_mid)); + + float t = (C - k_0) / (k_1 + k_2 * (C - k_0)); + s = mid + (1.f - mid) * t; + } + + float l = toe(L); + return vec3( h, s, l ); +} + + +vec3 okhsv_to_srgb(vec3 hsv) +{ + float h = hsv.x; + float s = hsv.y; + float v = hsv.z; + + float a_ = cos(2.f * M_PI * h); + float b_ = sin(2.f * M_PI * h); + + vec2 cusp = find_cusp(a_, b_); + vec2 ST_max = to_ST(cusp); + float S_max = ST_max.x; + float T_max = ST_max.y; + float S_0 = 0.5f; + float k = 1.f- S_0 / S_max; + + // first we compute L and V as if the gamut is a perfect triangle: + + // L, C when v==1: + float L_v = 1.f - s * S_0 / (S_0 + T_max - T_max * k * s); + float C_v = s * T_max * S_0 / (S_0 + T_max - T_max * k * s); + + float L = v * L_v; + float C = v * C_v; + + // then we compensate for both toe and the curved top part of the triangle: + float L_vt = toe_inv(L_v); + float C_vt = C_v * L_vt / L_v; + + float L_new = toe_inv(L); + C = C * L_new / L; + L = L_new; + + vec3 rgb_scale = oklab_to_linear_srgb(vec3( L_vt, a_ * C_vt, b_ * C_vt )); + float scale_L = cbrt(1.f / max(max(rgb_scale.r, rgb_scale.g), max(rgb_scale.b, 0.f))); + + L = L * scale_L; + C = C * scale_L; + + vec3 rgb = oklab_to_linear_srgb(vec3( L, C * a_, C * b_ )); + return vec3( + srgb_transfer_function(rgb.r), + srgb_transfer_function(rgb.g), + srgb_transfer_function(rgb.b) + ); +} +)" +R"( +vec3 srgb_to_okhsv(vec3 rgb) +{ + vec3 lab = linear_srgb_to_oklab(vec3( + srgb_transfer_function_inv(rgb.r), + srgb_transfer_function_inv(rgb.g), + srgb_transfer_function_inv(rgb.b) + )); + + float C = sqrt(lab.y * lab.y + lab.z * lab.z); + float a_ = lab.y / C; + float b_ = lab.z / C; + + float L = lab.x; + float h = 0.5f + 0.5f * atan(-lab.z, -lab.y) / M_PI; + + vec2 cusp = find_cusp(a_, b_); + vec2 ST_max = to_ST(cusp); + float S_max = ST_max.x; + float T_max = ST_max.y; + float S_0 = 0.5f; + float k = 1.f - S_0 / S_max; + + // first we find L_v, C_v, L_vt and C_vt + + float t = T_max / (C + L * T_max); + float L_v = t * L; + float C_v = t * C; + + float L_vt = toe_inv(L_v); + float C_vt = C_v * L_vt / L_v; + + // we can then use these to invert the step that compensates for the toe and the curved top part of the triangle: + vec3 rgb_scale = oklab_to_linear_srgb(vec3( L_vt, a_ * C_vt, b_ * C_vt )); + float scale_L = cbrt(1.f / max(max(rgb_scale.r, rgb_scale.g), max(rgb_scale.b, 0.f))); + + L = L / scale_L; + C = C / scale_L; + + C = C * toe(L) / L; + L = toe(L); + + // we can now compute v and s: + + float v = L / L_v; + float s = (S_0 + T_max) * C_v / ((T_max * S_0) + T_max * k * C_v); + + return vec3 (h, s, v ); +})"; + +#endif diff --git a/thirdparty/misc/open-simplex-noise-LICENSE b/thirdparty/misc/open-simplex-noise-LICENSE deleted file mode 100644 index a84c395662..0000000000 --- a/thirdparty/misc/open-simplex-noise-LICENSE +++ /dev/null @@ -1,25 +0,0 @@ -This is free and unencumbered software released into the public domain. - -Anyone is free to copy, modify, publish, use, compile, sell, or -distribute this software, either in source code form or as a compiled -binary, for any purpose, commercial or non-commercial, and by any -means. - -In jurisdictions that recognize copyright laws, the author or authors -of this software dedicate any and all copyright interest in the -software to the public domain. We make this dedication for the benefit -of the public at large and to the detriment of our heirs and -successors. We intend this dedication to be an overt act of -relinquishment in perpetuity of all present and future rights to this -software under copyright law. - -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 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. - -For more information, please refer to <http://unlicense.org> - diff --git a/thirdparty/misc/open-simplex-noise-no-allocate.patch b/thirdparty/misc/open-simplex-noise-no-allocate.patch deleted file mode 100644 index fc3abe7d00..0000000000 --- a/thirdparty/misc/open-simplex-noise-no-allocate.patch +++ /dev/null @@ -1,133 +0,0 @@ -diff -u orig/open-simplex-noise.c misc/open-simplex-noise.c ---- orig/open-simplex-noise.c 2018-09-14 11:11:40.049810000 +0200 -+++ misc/open-simplex-noise.c 2018-09-14 11:09:39.726457000 +0200 -@@ -13,6 +13,11 @@ - * of any particular randomization library, so results - * will be the same when ported to other languages. - */ -+ -+// -- GODOT start -- -+// Modified to work without allocating memory, also removed some unused function. -+// -- GODOT end -- -+ - #include <math.h> - #include <stdlib.h> - #include <stdint.h> -@@ -34,11 +39,12 @@ - - #define DEFAULT_SEED (0LL) - --struct osn_context { -+// -- GODOT start -- -+/*struct osn_context { - int16_t *perm; - int16_t *permGradIndex3D; --}; -- -+};*/ -+// -- GODOT end -- - #define ARRAYSIZE(x) (sizeof((x)) / sizeof((x)[0])) - - /* -@@ -126,7 +132,9 @@ - int xi = (int) x; - return x < xi ? xi - 1 : xi; - } -- -+ -+// -- GODOT start -- -+/* - static int allocate_perm(struct osn_context *ctx, int nperm, int ngrad) - { - if (ctx->perm) -@@ -154,18 +162,21 @@ - memcpy(ctx->perm, p, sizeof(*ctx->perm) * nelements); - - for (i = 0; i < 256; i++) { -- /* Since 3D has 24 gradients, simple bitmask won't work, so precompute modulo array. */ -+ // Since 3D has 24 gradients, simple bitmask won't work, so precompute modulo array. - ctx->permGradIndex3D[i] = (int16_t)((ctx->perm[i] % (ARRAYSIZE(gradients3D) / 3)) * 3); - } - return 0; - } -+*/ -+// -- GODOT end -- - - /* - * Initializes using a permutation array generated from a 64-bit seed. - * Generates a proper permutation (i.e. doesn't merely perform N successive pair - * swaps on a base array). Uses a simple 64-bit LCG. - */ --int open_simplex_noise(int64_t seed, struct osn_context **ctx) -+// -- GODOT start -- -+int open_simplex_noise(int64_t seed, struct osn_context *ctx) - { - int rc; - int16_t source[256]; -@@ -174,20 +185,9 @@ - int16_t *permGradIndex3D; - int r; - -- *ctx = (struct osn_context *) malloc(sizeof(**ctx)); -- if (!(*ctx)) -- return -ENOMEM; -- (*ctx)->perm = NULL; -- (*ctx)->permGradIndex3D = NULL; -- -- rc = allocate_perm(*ctx, 256, 256); -- if (rc) { -- free(*ctx); -- return rc; -- } -- -- perm = (*ctx)->perm; -- permGradIndex3D = (*ctx)->permGradIndex3D; -+ perm = ctx->perm; -+ permGradIndex3D = ctx->permGradIndex3D; -+// -- GODOT end -- - - for (i = 0; i < 256; i++) - source[i] = (int16_t) i; -@@ -206,6 +206,8 @@ - return 0; - } - -+// -- GODOT start -- -+/* - void open_simplex_noise_free(struct osn_context *ctx) - { - if (!ctx) -@@ -220,6 +222,8 @@ - } - free(ctx); - } -+*/ -+// -- GODOT end -- - - /* 2D OpenSimplex (Simplectic) Noise. */ - double open_simplex_noise2(struct osn_context *ctx, double x, double y) -diff -u orig/open-simplex-noise.h misc/open-simplex-noise.h ---- orig/open-simplex-noise.h 2018-09-14 11:11:19.659807000 +0200 -+++ misc/open-simplex-noise.h 2018-09-14 11:10:05.006460000 +0200 -@@ -35,11 +35,18 @@ - extern "C" { - #endif - --struct osn_context; -+// -- GODOT start -- -+// Modified to work without allocating memory, also removed some unused function. - --int open_simplex_noise(int64_t seed, struct osn_context **ctx); -+struct osn_context { -+ int16_t perm[256]; -+ int16_t permGradIndex3D[256]; -+}; -+ -+int open_simplex_noise(int64_t seed, struct osn_context *ctx); -+//int open_simplex_noise_init_perm(struct osn_context *ctx, int16_t p[], int nelements); -+// -- GODOT end -- - void open_simplex_noise_free(struct osn_context *ctx); --int open_simplex_noise_init_perm(struct osn_context *ctx, int16_t p[], int nelements); - double open_simplex_noise2(struct osn_context *ctx, double x, double y); - double open_simplex_noise3(struct osn_context *ctx, double x, double y, double z); - double open_simplex_noise4(struct osn_context *ctx, double x, double y, double z, double w); diff --git a/thirdparty/misc/open-simplex-noise.c b/thirdparty/misc/open-simplex-noise.c deleted file mode 100644 index 44a072cad1..0000000000 --- a/thirdparty/misc/open-simplex-noise.c +++ /dev/null @@ -1,2255 +0,0 @@ -/* - * OpenSimplex (Simplectic) Noise in C. - * Ported by Stephen M. Cameron from Kurt Spencer's java implementation - * - * v1.1 (October 5, 2014) - * - Added 2D and 4D implementations. - * - Proper gradient sets for all dimensions, from a - * dimensionally-generalizable scheme with an actual - * rhyme and reason behind it. - * - Removed default permutation array in favor of - * default seed. - * - Changed seed-based constructor to be independent - * of any particular randomization library, so results - * will be the same when ported to other languages. - */ - -// -- GODOT start -- -// Modified to work without allocating memory, also removed some unused function. -// -- GODOT end -- - -#include <math.h> -#include <stdlib.h> -#include <stdint.h> -#include <string.h> -#include <errno.h> - -#include "open-simplex-noise.h" - -#define STRETCH_CONSTANT_2D (-0.211324865405187) /* (1 / sqrt(2 + 1) - 1 ) / 2; */ -#define SQUISH_CONSTANT_2D (0.366025403784439) /* (sqrt(2 + 1) -1) / 2; */ -#define STRETCH_CONSTANT_3D (-1.0 / 6.0) /* (1 / sqrt(3 + 1) - 1) / 3; */ -#define SQUISH_CONSTANT_3D (1.0 / 3.0) /* (sqrt(3+1)-1)/3; */ -#define STRETCH_CONSTANT_4D (-0.138196601125011) /* (1 / sqrt(4 + 1) - 1) / 4; */ -#define SQUISH_CONSTANT_4D (0.309016994374947) /* (sqrt(4 + 1) - 1) / 4; */ - -#define NORM_CONSTANT_2D (47.0) -#define NORM_CONSTANT_3D (103.0) -#define NORM_CONSTANT_4D (30.0) - -#define DEFAULT_SEED (0LL) - -// -- GODOT start -- -/*struct osn_context { - int16_t *perm; - int16_t *permGradIndex3D; -};*/ -// -- GODOT end -- -#define ARRAYSIZE(x) (sizeof((x)) / sizeof((x)[0])) - -/* - * Gradients for 2D. They approximate the directions to the - * vertices of an octagon from the center. - */ -static const int8_t gradients2D[] = { - 5, 2, 2, 5, - -5, 2, -2, 5, - 5, -2, 2, -5, - -5, -2, -2, -5, -}; - -/* - * Gradients for 3D. They approximate the directions to the - * vertices of a rhombicuboctahedron from the center, skewed so - * that the triangular and square facets can be inscribed inside - * circles of the same radius. - */ -static const signed char gradients3D[] = { - -11, 4, 4, -4, 11, 4, -4, 4, 11, - 11, 4, 4, 4, 11, 4, 4, 4, 11, - -11, -4, 4, -4, -11, 4, -4, -4, 11, - 11, -4, 4, 4, -11, 4, 4, -4, 11, - -11, 4, -4, -4, 11, -4, -4, 4, -11, - 11, 4, -4, 4, 11, -4, 4, 4, -11, - -11, -4, -4, -4, -11, -4, -4, -4, -11, - 11, -4, -4, 4, -11, -4, 4, -4, -11, -}; - -/* - * Gradients for 4D. They approximate the directions to the - * vertices of a disprismatotesseractihexadecachoron from the center, - * skewed so that the tetrahedral and cubic facets can be inscribed inside - * spheres of the same radius. - */ -static const signed char gradients4D[] = { - 3, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3, 1, 1, 1, 1, 3, - -3, 1, 1, 1, -1, 3, 1, 1, -1, 1, 3, 1, -1, 1, 1, 3, - 3, -1, 1, 1, 1, -3, 1, 1, 1, -1, 3, 1, 1, -1, 1, 3, - -3, -1, 1, 1, -1, -3, 1, 1, -1, -1, 3, 1, -1, -1, 1, 3, - 3, 1, -1, 1, 1, 3, -1, 1, 1, 1, -3, 1, 1, 1, -1, 3, - -3, 1, -1, 1, -1, 3, -1, 1, -1, 1, -3, 1, -1, 1, -1, 3, - 3, -1, -1, 1, 1, -3, -1, 1, 1, -1, -3, 1, 1, -1, -1, 3, - -3, -1, -1, 1, -1, -3, -1, 1, -1, -1, -3, 1, -1, -1, -1, 3, - 3, 1, 1, -1, 1, 3, 1, -1, 1, 1, 3, -1, 1, 1, 1, -3, - -3, 1, 1, -1, -1, 3, 1, -1, -1, 1, 3, -1, -1, 1, 1, -3, - 3, -1, 1, -1, 1, -3, 1, -1, 1, -1, 3, -1, 1, -1, 1, -3, - -3, -1, 1, -1, -1, -3, 1, -1, -1, -1, 3, -1, -1, -1, 1, -3, - 3, 1, -1, -1, 1, 3, -1, -1, 1, 1, -3, -1, 1, 1, -1, -3, - -3, 1, -1, -1, -1, 3, -1, -1, -1, 1, -3, -1, -1, 1, -1, -3, - 3, -1, -1, -1, 1, -3, -1, -1, 1, -1, -3, -1, 1, -1, -1, -3, - -3, -1, -1, -1, -1, -3, -1, -1, -1, -1, -3, -1, -1, -1, -1, -3, -}; - -static double extrapolate2(const struct osn_context *ctx, int xsb, int ysb, double dx, double dy) -{ - const int16_t *perm = ctx->perm; - int index = perm[(perm[xsb & 0xFF] + ysb) & 0xFF] & 0x0E; - return gradients2D[index] * dx - + gradients2D[index + 1] * dy; -} - -static double extrapolate3(const struct osn_context *ctx, int xsb, int ysb, int zsb, double dx, double dy, double dz) -{ - const int16_t *perm = ctx->perm; - const int16_t *permGradIndex3D = ctx->permGradIndex3D; - int index = permGradIndex3D[(perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF]; - return gradients3D[index] * dx - + gradients3D[index + 1] * dy - + gradients3D[index + 2] * dz; -} - -static double extrapolate4(const struct osn_context *ctx, int xsb, int ysb, int zsb, int wsb, double dx, double dy, double dz, double dw) -{ - const int16_t *perm = ctx->perm; - int index = perm[(perm[(perm[(perm[xsb & 0xFF] + ysb) & 0xFF] + zsb) & 0xFF] + wsb) & 0xFF] & 0xFC; - return gradients4D[index] * dx - + gradients4D[index + 1] * dy - + gradients4D[index + 2] * dz - + gradients4D[index + 3] * dw; -} - -static INLINE int fastFloor(double x) { - int xi = (int) x; - return x < xi ? xi - 1 : xi; -} - -// -- GODOT start -- -/* -static int allocate_perm(struct osn_context *ctx, int nperm, int ngrad) -{ - if (ctx->perm) - free(ctx->perm); - if (ctx->permGradIndex3D) - free(ctx->permGradIndex3D); - ctx->perm = (int16_t *) malloc(sizeof(*ctx->perm) * nperm); - if (!ctx->perm) - return -ENOMEM; - ctx->permGradIndex3D = (int16_t *) malloc(sizeof(*ctx->permGradIndex3D) * ngrad); - if (!ctx->permGradIndex3D) { - free(ctx->perm); - return -ENOMEM; - } - return 0; -} - -int open_simplex_noise_init_perm(struct osn_context *ctx, int16_t p[], int nelements) -{ - int i, rc; - - rc = allocate_perm(ctx, nelements, 256); - if (rc) - return rc; - memcpy(ctx->perm, p, sizeof(*ctx->perm) * nelements); - - for (i = 0; i < 256; i++) { - // Since 3D has 24 gradients, simple bitmask won't work, so precompute modulo array. - ctx->permGradIndex3D[i] = (int16_t)((ctx->perm[i] % (ARRAYSIZE(gradients3D) / 3)) * 3); - } - return 0; -} -*/ -// -- GODOT end -- - -/* - * Initializes using a permutation array generated from a 64-bit seed. - * Generates a proper permutation (i.e. doesn't merely perform N successive pair - * swaps on a base array). Uses a simple 64-bit LCG. - */ -// -- GODOT start -- -int open_simplex_noise(int64_t seed, struct osn_context *ctx) -{ - int rc; - int16_t source[256]; - int i; - int16_t *perm; - int16_t *permGradIndex3D; - int r; - - perm = ctx->perm; - permGradIndex3D = ctx->permGradIndex3D; -// -- GODOT end -- - - uint64_t seedU = seed; - for (i = 0; i < 256; i++) - source[i] = (int16_t) i; - seedU = seedU * 6364136223846793005ULL + 1442695040888963407ULL; - seedU = seedU * 6364136223846793005ULL + 1442695040888963407ULL; - seedU = seedU * 6364136223846793005ULL + 1442695040888963407ULL; - for (i = 255; i >= 0; i--) { - seedU = seedU * 6364136223846793005ULL + 1442695040888963407ULL; - r = (int)((seedU + 31) % (i + 1)); - if (r < 0) - r += (i + 1); - perm[i] = source[r]; - permGradIndex3D[i] = (short)((perm[i] % (ARRAYSIZE(gradients3D) / 3)) * 3); - source[r] = source[i]; - } - return 0; -} - -// -- GODOT start -- -/* -void open_simplex_noise_free(struct osn_context *ctx) -{ - if (!ctx) - return; - if (ctx->perm) { - free(ctx->perm); - ctx->perm = NULL; - } - if (ctx->permGradIndex3D) { - free(ctx->permGradIndex3D); - ctx->permGradIndex3D = NULL; - } - free(ctx); -} -*/ -// -- GODOT end -- - -/* 2D OpenSimplex (Simplectic) Noise. */ -double open_simplex_noise2(const struct osn_context *ctx, double x, double y) -{ - - /* Place input coordinates onto grid. */ - double stretchOffset = (x + y) * STRETCH_CONSTANT_2D; - double xs = x + stretchOffset; - double ys = y + stretchOffset; - - /* Floor to get grid coordinates of rhombus (stretched square) super-cell origin. */ - int xsb = fastFloor(xs); - int ysb = fastFloor(ys); - - /* Skew out to get actual coordinates of rhombus origin. We'll need these later. */ - double squishOffset = (xsb + ysb) * SQUISH_CONSTANT_2D; - double xb = xsb + squishOffset; - double yb = ysb + squishOffset; - - /* Compute grid coordinates relative to rhombus origin. */ - double xins = xs - xsb; - double yins = ys - ysb; - - /* Sum those together to get a value that determines which region we're in. */ - double inSum = xins + yins; - - /* Positions relative to origin point. */ - double dx0 = x - xb; - double dy0 = y - yb; - - /* We'll be defining these inside the next block and using them afterwards. */ - double dx_ext, dy_ext; - int xsv_ext, ysv_ext; - - double dx1; - double dy1; - double attn1; - double dx2; - double dy2; - double attn2; - double zins; - double attn0; - double attn_ext; - - double value = 0; - - /* Contribution (1,0) */ - dx1 = dx0 - 1 - SQUISH_CONSTANT_2D; - dy1 = dy0 - 0 - SQUISH_CONSTANT_2D; - attn1 = 2 - dx1 * dx1 - dy1 * dy1; - if (attn1 > 0) { - attn1 *= attn1; - value += attn1 * attn1 * extrapolate2(ctx, xsb + 1, ysb + 0, dx1, dy1); - } - - /* Contribution (0,1) */ - dx2 = dx0 - 0 - SQUISH_CONSTANT_2D; - dy2 = dy0 - 1 - SQUISH_CONSTANT_2D; - attn2 = 2 - dx2 * dx2 - dy2 * dy2; - if (attn2 > 0) { - attn2 *= attn2; - value += attn2 * attn2 * extrapolate2(ctx, xsb + 0, ysb + 1, dx2, dy2); - } - - if (inSum <= 1) { /* We're inside the triangle (2-Simplex) at (0,0) */ - zins = 1 - inSum; - if (zins > xins || zins > yins) { /* (0,0) is one of the closest two triangular vertices */ - if (xins > yins) { - xsv_ext = xsb + 1; - ysv_ext = ysb - 1; - dx_ext = dx0 - 1; - dy_ext = dy0 + 1; - } else { - xsv_ext = xsb - 1; - ysv_ext = ysb + 1; - dx_ext = dx0 + 1; - dy_ext = dy0 - 1; - } - } else { /* (1,0) and (0,1) are the closest two vertices. */ - xsv_ext = xsb + 1; - ysv_ext = ysb + 1; - dx_ext = dx0 - 1 - 2 * SQUISH_CONSTANT_2D; - dy_ext = dy0 - 1 - 2 * SQUISH_CONSTANT_2D; - } - } else { /* We're inside the triangle (2-Simplex) at (1,1) */ - zins = 2 - inSum; - if (zins < xins || zins < yins) { /* (0,0) is one of the closest two triangular vertices */ - if (xins > yins) { - xsv_ext = xsb + 2; - ysv_ext = ysb + 0; - dx_ext = dx0 - 2 - 2 * SQUISH_CONSTANT_2D; - dy_ext = dy0 + 0 - 2 * SQUISH_CONSTANT_2D; - } else { - xsv_ext = xsb + 0; - ysv_ext = ysb + 2; - dx_ext = dx0 + 0 - 2 * SQUISH_CONSTANT_2D; - dy_ext = dy0 - 2 - 2 * SQUISH_CONSTANT_2D; - } - } else { /* (1,0) and (0,1) are the closest two vertices. */ - dx_ext = dx0; - dy_ext = dy0; - xsv_ext = xsb; - ysv_ext = ysb; - } - xsb += 1; - ysb += 1; - dx0 = dx0 - 1 - 2 * SQUISH_CONSTANT_2D; - dy0 = dy0 - 1 - 2 * SQUISH_CONSTANT_2D; - } - - /* Contribution (0,0) or (1,1) */ - attn0 = 2 - dx0 * dx0 - dy0 * dy0; - if (attn0 > 0) { - attn0 *= attn0; - value += attn0 * attn0 * extrapolate2(ctx, xsb, ysb, dx0, dy0); - } - - /* Extra Vertex */ - attn_ext = 2 - dx_ext * dx_ext - dy_ext * dy_ext; - if (attn_ext > 0) { - attn_ext *= attn_ext; - value += attn_ext * attn_ext * extrapolate2(ctx, xsv_ext, ysv_ext, dx_ext, dy_ext); - } - - return value / NORM_CONSTANT_2D; -} - -/* - * 3D OpenSimplex (Simplectic) Noise - */ -double open_simplex_noise3(const struct osn_context *ctx, double x, double y, double z) -{ - - /* Place input coordinates on simplectic honeycomb. */ - double stretchOffset = (x + y + z) * STRETCH_CONSTANT_3D; - double xs = x + stretchOffset; - double ys = y + stretchOffset; - double zs = z + stretchOffset; - - /* Floor to get simplectic honeycomb coordinates of rhombohedron (stretched cube) super-cell origin. */ - int xsb = fastFloor(xs); - int ysb = fastFloor(ys); - int zsb = fastFloor(zs); - - /* Skew out to get actual coordinates of rhombohedron origin. We'll need these later. */ - double squishOffset = (xsb + ysb + zsb) * SQUISH_CONSTANT_3D; - double xb = xsb + squishOffset; - double yb = ysb + squishOffset; - double zb = zsb + squishOffset; - - /* Compute simplectic honeycomb coordinates relative to rhombohedral origin. */ - double xins = xs - xsb; - double yins = ys - ysb; - double zins = zs - zsb; - - /* Sum those together to get a value that determines which region we're in. */ - double inSum = xins + yins + zins; - - /* Positions relative to origin point. */ - double dx0 = x - xb; - double dy0 = y - yb; - double dz0 = z - zb; - - /* We'll be defining these inside the next block and using them afterwards. */ - double dx_ext0, dy_ext0, dz_ext0; - double dx_ext1, dy_ext1, dz_ext1; - int xsv_ext0, ysv_ext0, zsv_ext0; - int xsv_ext1, ysv_ext1, zsv_ext1; - - double wins; - int8_t c, c1, c2; - int8_t aPoint, bPoint; - double aScore, bScore; - int aIsFurtherSide; - int bIsFurtherSide; - double p1, p2, p3; - double score; - double attn0, attn1, attn2, attn3, attn4, attn5, attn6; - double dx1, dy1, dz1; - double dx2, dy2, dz2; - double dx3, dy3, dz3; - double dx4, dy4, dz4; - double dx5, dy5, dz5; - double dx6, dy6, dz6; - double attn_ext0, attn_ext1; - - double value = 0; - if (inSum <= 1) { /* We're inside the tetrahedron (3-Simplex) at (0,0,0) */ - - /* Determine which two of (0,0,1), (0,1,0), (1,0,0) are closest. */ - aPoint = 0x01; - aScore = xins; - bPoint = 0x02; - bScore = yins; - if (aScore >= bScore && zins > bScore) { - bScore = zins; - bPoint = 0x04; - } else if (aScore < bScore && zins > aScore) { - aScore = zins; - aPoint = 0x04; - } - - /* Now we determine the two lattice points not part of the tetrahedron that may contribute. - This depends on the closest two tetrahedral vertices, including (0,0,0) */ - wins = 1 - inSum; - if (wins > aScore || wins > bScore) { /* (0,0,0) is one of the closest two tetrahedral vertices. */ - c = (bScore > aScore ? bPoint : aPoint); /* Our other closest vertex is the closest out of a and b. */ - - if ((c & 0x01) == 0) { - xsv_ext0 = xsb - 1; - xsv_ext1 = xsb; - dx_ext0 = dx0 + 1; - dx_ext1 = dx0; - } else { - xsv_ext0 = xsv_ext1 = xsb + 1; - dx_ext0 = dx_ext1 = dx0 - 1; - } - - if ((c & 0x02) == 0) { - ysv_ext0 = ysv_ext1 = ysb; - dy_ext0 = dy_ext1 = dy0; - if ((c & 0x01) == 0) { - ysv_ext1 -= 1; - dy_ext1 += 1; - } else { - ysv_ext0 -= 1; - dy_ext0 += 1; - } - } else { - ysv_ext0 = ysv_ext1 = ysb + 1; - dy_ext0 = dy_ext1 = dy0 - 1; - } - - if ((c & 0x04) == 0) { - zsv_ext0 = zsb; - zsv_ext1 = zsb - 1; - dz_ext0 = dz0; - dz_ext1 = dz0 + 1; - } else { - zsv_ext0 = zsv_ext1 = zsb + 1; - dz_ext0 = dz_ext1 = dz0 - 1; - } - } else { /* (0,0,0) is not one of the closest two tetrahedral vertices. */ - c = (int8_t)(aPoint | bPoint); /* Our two extra vertices are determined by the closest two. */ - - if ((c & 0x01) == 0) { - xsv_ext0 = xsb; - xsv_ext1 = xsb - 1; - dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_3D; - dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_3D; - } else { - xsv_ext0 = xsv_ext1 = xsb + 1; - dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; - dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; - } - - if ((c & 0x02) == 0) { - ysv_ext0 = ysb; - ysv_ext1 = ysb - 1; - dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_3D; - dy_ext1 = dy0 + 1 - SQUISH_CONSTANT_3D; - } else { - ysv_ext0 = ysv_ext1 = ysb + 1; - dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; - dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; - } - - if ((c & 0x04) == 0) { - zsv_ext0 = zsb; - zsv_ext1 = zsb - 1; - dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_3D; - dz_ext1 = dz0 + 1 - SQUISH_CONSTANT_3D; - } else { - zsv_ext0 = zsv_ext1 = zsb + 1; - dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; - dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; - } - } - - /* Contribution (0,0,0) */ - attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0; - if (attn0 > 0) { - attn0 *= attn0; - value += attn0 * attn0 * extrapolate3(ctx, xsb + 0, ysb + 0, zsb + 0, dx0, dy0, dz0); - } - - /* Contribution (1,0,0) */ - dx1 = dx0 - 1 - SQUISH_CONSTANT_3D; - dy1 = dy0 - 0 - SQUISH_CONSTANT_3D; - dz1 = dz0 - 0 - SQUISH_CONSTANT_3D; - attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; - if (attn1 > 0) { - attn1 *= attn1; - value += attn1 * attn1 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 0, dx1, dy1, dz1); - } - - /* Contribution (0,1,0) */ - dx2 = dx0 - 0 - SQUISH_CONSTANT_3D; - dy2 = dy0 - 1 - SQUISH_CONSTANT_3D; - dz2 = dz1; - attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; - if (attn2 > 0) { - attn2 *= attn2; - value += attn2 * attn2 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 0, dx2, dy2, dz2); - } - - /* Contribution (0,0,1) */ - dx3 = dx2; - dy3 = dy1; - dz3 = dz0 - 1 - SQUISH_CONSTANT_3D; - attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; - if (attn3 > 0) { - attn3 *= attn3; - value += attn3 * attn3 * extrapolate3(ctx, xsb + 0, ysb + 0, zsb + 1, dx3, dy3, dz3); - } - } else if (inSum >= 2) { /* We're inside the tetrahedron (3-Simplex) at (1,1,1) */ - - /* Determine which two tetrahedral vertices are the closest, out of (1,1,0), (1,0,1), (0,1,1) but not (1,1,1). */ - aPoint = 0x06; - aScore = xins; - bPoint = 0x05; - bScore = yins; - if (aScore <= bScore && zins < bScore) { - bScore = zins; - bPoint = 0x03; - } else if (aScore > bScore && zins < aScore) { - aScore = zins; - aPoint = 0x03; - } - - /* Now we determine the two lattice points not part of the tetrahedron that may contribute. - This depends on the closest two tetrahedral vertices, including (1,1,1) */ - wins = 3 - inSum; - if (wins < aScore || wins < bScore) { /* (1,1,1) is one of the closest two tetrahedral vertices. */ - c = (bScore < aScore ? bPoint : aPoint); /* Our other closest vertex is the closest out of a and b. */ - - if ((c & 0x01) != 0) { - xsv_ext0 = xsb + 2; - xsv_ext1 = xsb + 1; - dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_3D; - dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; - } else { - xsv_ext0 = xsv_ext1 = xsb; - dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_3D; - } - - if ((c & 0x02) != 0) { - ysv_ext0 = ysv_ext1 = ysb + 1; - dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; - if ((c & 0x01) != 0) { - ysv_ext1 += 1; - dy_ext1 -= 1; - } else { - ysv_ext0 += 1; - dy_ext0 -= 1; - } - } else { - ysv_ext0 = ysv_ext1 = ysb; - dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_3D; - } - - if ((c & 0x04) != 0) { - zsv_ext0 = zsb + 1; - zsv_ext1 = zsb + 2; - dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; - dz_ext1 = dz0 - 2 - 3 * SQUISH_CONSTANT_3D; - } else { - zsv_ext0 = zsv_ext1 = zsb; - dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_3D; - } - } else { /* (1,1,1) is not one of the closest two tetrahedral vertices. */ - c = (int8_t)(aPoint & bPoint); /* Our two extra vertices are determined by the closest two. */ - - if ((c & 0x01) != 0) { - xsv_ext0 = xsb + 1; - xsv_ext1 = xsb + 2; - dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; - dx_ext1 = dx0 - 2 - 2 * SQUISH_CONSTANT_3D; - } else { - xsv_ext0 = xsv_ext1 = xsb; - dx_ext0 = dx0 - SQUISH_CONSTANT_3D; - dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; - } - - if ((c & 0x02) != 0) { - ysv_ext0 = ysb + 1; - ysv_ext1 = ysb + 2; - dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; - dy_ext1 = dy0 - 2 - 2 * SQUISH_CONSTANT_3D; - } else { - ysv_ext0 = ysv_ext1 = ysb; - dy_ext0 = dy0 - SQUISH_CONSTANT_3D; - dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; - } - - if ((c & 0x04) != 0) { - zsv_ext0 = zsb + 1; - zsv_ext1 = zsb + 2; - dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; - dz_ext1 = dz0 - 2 - 2 * SQUISH_CONSTANT_3D; - } else { - zsv_ext0 = zsv_ext1 = zsb; - dz_ext0 = dz0 - SQUISH_CONSTANT_3D; - dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; - } - } - - /* Contribution (1,1,0) */ - dx3 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; - dy3 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; - dz3 = dz0 - 0 - 2 * SQUISH_CONSTANT_3D; - attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; - if (attn3 > 0) { - attn3 *= attn3; - value += attn3 * attn3 * extrapolate3(ctx, xsb + 1, ysb + 1, zsb + 0, dx3, dy3, dz3); - } - - /* Contribution (1,0,1) */ - dx2 = dx3; - dy2 = dy0 - 0 - 2 * SQUISH_CONSTANT_3D; - dz2 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; - attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; - if (attn2 > 0) { - attn2 *= attn2; - value += attn2 * attn2 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 1, dx2, dy2, dz2); - } - - /* Contribution (0,1,1) */ - dx1 = dx0 - 0 - 2 * SQUISH_CONSTANT_3D; - dy1 = dy3; - dz1 = dz2; - attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; - if (attn1 > 0) { - attn1 *= attn1; - value += attn1 * attn1 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 1, dx1, dy1, dz1); - } - - /* Contribution (1,1,1) */ - dx0 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; - dy0 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; - dz0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; - attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0; - if (attn0 > 0) { - attn0 *= attn0; - value += attn0 * attn0 * extrapolate3(ctx, xsb + 1, ysb + 1, zsb + 1, dx0, dy0, dz0); - } - } else { /* We're inside the octahedron (Rectified 3-Simplex) in between. - Decide between point (0,0,1) and (1,1,0) as closest */ - p1 = xins + yins; - if (p1 > 1) { - aScore = p1 - 1; - aPoint = 0x03; - aIsFurtherSide = 1; - } else { - aScore = 1 - p1; - aPoint = 0x04; - aIsFurtherSide = 0; - } - - /* Decide between point (0,1,0) and (1,0,1) as closest */ - p2 = xins + zins; - if (p2 > 1) { - bScore = p2 - 1; - bPoint = 0x05; - bIsFurtherSide = 1; - } else { - bScore = 1 - p2; - bPoint = 0x02; - bIsFurtherSide = 0; - } - - /* The closest out of the two (1,0,0) and (0,1,1) will replace the furthest out of the two decided above, if closer. */ - p3 = yins + zins; - if (p3 > 1) { - score = p3 - 1; - if (aScore <= bScore && aScore < score) { - aScore = score; - aPoint = 0x06; - aIsFurtherSide = 1; - } else if (aScore > bScore && bScore < score) { - bScore = score; - bPoint = 0x06; - bIsFurtherSide = 1; - } - } else { - score = 1 - p3; - if (aScore <= bScore && aScore < score) { - aScore = score; - aPoint = 0x01; - aIsFurtherSide = 0; - } else if (aScore > bScore && bScore < score) { - bScore = score; - bPoint = 0x01; - bIsFurtherSide = 0; - } - } - - /* Where each of the two closest points are determines how the extra two vertices are calculated. */ - if (aIsFurtherSide == bIsFurtherSide) { - if (aIsFurtherSide) { /* Both closest points on (1,1,1) side */ - - /* One of the two extra points is (1,1,1) */ - dx_ext0 = dx0 - 1 - 3 * SQUISH_CONSTANT_3D; - dy_ext0 = dy0 - 1 - 3 * SQUISH_CONSTANT_3D; - dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_3D; - xsv_ext0 = xsb + 1; - ysv_ext0 = ysb + 1; - zsv_ext0 = zsb + 1; - - /* Other extra point is based on the shared axis. */ - c = (int8_t)(aPoint & bPoint); - if ((c & 0x01) != 0) { - dx_ext1 = dx0 - 2 - 2 * SQUISH_CONSTANT_3D; - dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; - dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; - xsv_ext1 = xsb + 2; - ysv_ext1 = ysb; - zsv_ext1 = zsb; - } else if ((c & 0x02) != 0) { - dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; - dy_ext1 = dy0 - 2 - 2 * SQUISH_CONSTANT_3D; - dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; - xsv_ext1 = xsb; - ysv_ext1 = ysb + 2; - zsv_ext1 = zsb; - } else { - dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; - dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; - dz_ext1 = dz0 - 2 - 2 * SQUISH_CONSTANT_3D; - xsv_ext1 = xsb; - ysv_ext1 = ysb; - zsv_ext1 = zsb + 2; - } - } else { /* Both closest points on (0,0,0) side */ - - /* One of the two extra points is (0,0,0) */ - dx_ext0 = dx0; - dy_ext0 = dy0; - dz_ext0 = dz0; - xsv_ext0 = xsb; - ysv_ext0 = ysb; - zsv_ext0 = zsb; - - /* Other extra point is based on the omitted axis. */ - c = (int8_t)(aPoint | bPoint); - if ((c & 0x01) == 0) { - dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_3D; - dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; - dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; - xsv_ext1 = xsb - 1; - ysv_ext1 = ysb + 1; - zsv_ext1 = zsb + 1; - } else if ((c & 0x02) == 0) { - dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; - dy_ext1 = dy0 + 1 - SQUISH_CONSTANT_3D; - dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_3D; - xsv_ext1 = xsb + 1; - ysv_ext1 = ysb - 1; - zsv_ext1 = zsb + 1; - } else { - dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_3D; - dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_3D; - dz_ext1 = dz0 + 1 - SQUISH_CONSTANT_3D; - xsv_ext1 = xsb + 1; - ysv_ext1 = ysb + 1; - zsv_ext1 = zsb - 1; - } - } - } else { /* One point on (0,0,0) side, one point on (1,1,1) side */ - if (aIsFurtherSide) { - c1 = aPoint; - c2 = bPoint; - } else { - c1 = bPoint; - c2 = aPoint; - } - - /* One contribution is a permutation of (1,1,-1) */ - if ((c1 & 0x01) == 0) { - dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_3D; - dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; - dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; - xsv_ext0 = xsb - 1; - ysv_ext0 = ysb + 1; - zsv_ext0 = zsb + 1; - } else if ((c1 & 0x02) == 0) { - dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; - dy_ext0 = dy0 + 1 - SQUISH_CONSTANT_3D; - dz_ext0 = dz0 - 1 - SQUISH_CONSTANT_3D; - xsv_ext0 = xsb + 1; - ysv_ext0 = ysb - 1; - zsv_ext0 = zsb + 1; - } else { - dx_ext0 = dx0 - 1 - SQUISH_CONSTANT_3D; - dy_ext0 = dy0 - 1 - SQUISH_CONSTANT_3D; - dz_ext0 = dz0 + 1 - SQUISH_CONSTANT_3D; - xsv_ext0 = xsb + 1; - ysv_ext0 = ysb + 1; - zsv_ext0 = zsb - 1; - } - - /* One contribution is a permutation of (0,0,2) */ - dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_3D; - dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_3D; - dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_3D; - xsv_ext1 = xsb; - ysv_ext1 = ysb; - zsv_ext1 = zsb; - if ((c2 & 0x01) != 0) { - dx_ext1 -= 2; - xsv_ext1 += 2; - } else if ((c2 & 0x02) != 0) { - dy_ext1 -= 2; - ysv_ext1 += 2; - } else { - dz_ext1 -= 2; - zsv_ext1 += 2; - } - } - - /* Contribution (1,0,0) */ - dx1 = dx0 - 1 - SQUISH_CONSTANT_3D; - dy1 = dy0 - 0 - SQUISH_CONSTANT_3D; - dz1 = dz0 - 0 - SQUISH_CONSTANT_3D; - attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1; - if (attn1 > 0) { - attn1 *= attn1; - value += attn1 * attn1 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 0, dx1, dy1, dz1); - } - - /* Contribution (0,1,0) */ - dx2 = dx0 - 0 - SQUISH_CONSTANT_3D; - dy2 = dy0 - 1 - SQUISH_CONSTANT_3D; - dz2 = dz1; - attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2; - if (attn2 > 0) { - attn2 *= attn2; - value += attn2 * attn2 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 0, dx2, dy2, dz2); - } - - /* Contribution (0,0,1) */ - dx3 = dx2; - dy3 = dy1; - dz3 = dz0 - 1 - SQUISH_CONSTANT_3D; - attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3; - if (attn3 > 0) { - attn3 *= attn3; - value += attn3 * attn3 * extrapolate3(ctx, xsb + 0, ysb + 0, zsb + 1, dx3, dy3, dz3); - } - - /* Contribution (1,1,0) */ - dx4 = dx0 - 1 - 2 * SQUISH_CONSTANT_3D; - dy4 = dy0 - 1 - 2 * SQUISH_CONSTANT_3D; - dz4 = dz0 - 0 - 2 * SQUISH_CONSTANT_3D; - attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4; - if (attn4 > 0) { - attn4 *= attn4; - value += attn4 * attn4 * extrapolate3(ctx, xsb + 1, ysb + 1, zsb + 0, dx4, dy4, dz4); - } - - /* Contribution (1,0,1) */ - dx5 = dx4; - dy5 = dy0 - 0 - 2 * SQUISH_CONSTANT_3D; - dz5 = dz0 - 1 - 2 * SQUISH_CONSTANT_3D; - attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5; - if (attn5 > 0) { - attn5 *= attn5; - value += attn5 * attn5 * extrapolate3(ctx, xsb + 1, ysb + 0, zsb + 1, dx5, dy5, dz5); - } - - /* Contribution (0,1,1) */ - dx6 = dx0 - 0 - 2 * SQUISH_CONSTANT_3D; - dy6 = dy4; - dz6 = dz5; - attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6; - if (attn6 > 0) { - attn6 *= attn6; - value += attn6 * attn6 * extrapolate3(ctx, xsb + 0, ysb + 1, zsb + 1, dx6, dy6, dz6); - } - } - - /* First extra vertex */ - attn_ext0 = 2 - dx_ext0 * dx_ext0 - dy_ext0 * dy_ext0 - dz_ext0 * dz_ext0; - if (attn_ext0 > 0) - { - attn_ext0 *= attn_ext0; - value += attn_ext0 * attn_ext0 * extrapolate3(ctx, xsv_ext0, ysv_ext0, zsv_ext0, dx_ext0, dy_ext0, dz_ext0); - } - - /* Second extra vertex */ - attn_ext1 = 2 - dx_ext1 * dx_ext1 - dy_ext1 * dy_ext1 - dz_ext1 * dz_ext1; - if (attn_ext1 > 0) - { - attn_ext1 *= attn_ext1; - value += attn_ext1 * attn_ext1 * extrapolate3(ctx, xsv_ext1, ysv_ext1, zsv_ext1, dx_ext1, dy_ext1, dz_ext1); - } - - return value / NORM_CONSTANT_3D; -} - -/* - * 4D OpenSimplex (Simplectic) Noise. - */ -double open_simplex_noise4(const struct osn_context *ctx, double x, double y, double z, double w) -{ - double uins; - double dx1, dy1, dz1, dw1; - double dx2, dy2, dz2, dw2; - double dx3, dy3, dz3, dw3; - double dx4, dy4, dz4, dw4; - double dx5, dy5, dz5, dw5; - double dx6, dy6, dz6, dw6; - double dx7, dy7, dz7, dw7; - double dx8, dy8, dz8, dw8; - double dx9, dy9, dz9, dw9; - double dx10, dy10, dz10, dw10; - double attn0, attn1, attn2, attn3, attn4; - double attn5, attn6, attn7, attn8, attn9, attn10; - double attn_ext0, attn_ext1, attn_ext2; - int8_t c, c1, c2; - int8_t aPoint, bPoint; - double aScore, bScore; - int aIsBiggerSide; - int bIsBiggerSide; - double p1, p2, p3, p4; - double score; - - /* Place input coordinates on simplectic honeycomb. */ - double stretchOffset = (x + y + z + w) * STRETCH_CONSTANT_4D; - double xs = x + stretchOffset; - double ys = y + stretchOffset; - double zs = z + stretchOffset; - double ws = w + stretchOffset; - - /* Floor to get simplectic honeycomb coordinates of rhombo-hypercube super-cell origin. */ - int xsb = fastFloor(xs); - int ysb = fastFloor(ys); - int zsb = fastFloor(zs); - int wsb = fastFloor(ws); - - /* Skew out to get actual coordinates of stretched rhombo-hypercube origin. We'll need these later. */ - double squishOffset = (xsb + ysb + zsb + wsb) * SQUISH_CONSTANT_4D; - double xb = xsb + squishOffset; - double yb = ysb + squishOffset; - double zb = zsb + squishOffset; - double wb = wsb + squishOffset; - - /* Compute simplectic honeycomb coordinates relative to rhombo-hypercube origin. */ - double xins = xs - xsb; - double yins = ys - ysb; - double zins = zs - zsb; - double wins = ws - wsb; - - /* Sum those together to get a value that determines which region we're in. */ - double inSum = xins + yins + zins + wins; - - /* Positions relative to origin point. */ - double dx0 = x - xb; - double dy0 = y - yb; - double dz0 = z - zb; - double dw0 = w - wb; - - /* We'll be defining these inside the next block and using them afterwards. */ - double dx_ext0, dy_ext0, dz_ext0, dw_ext0; - double dx_ext1, dy_ext1, dz_ext1, dw_ext1; - double dx_ext2, dy_ext2, dz_ext2, dw_ext2; - int xsv_ext0, ysv_ext0, zsv_ext0, wsv_ext0; - int xsv_ext1, ysv_ext1, zsv_ext1, wsv_ext1; - int xsv_ext2, ysv_ext2, zsv_ext2, wsv_ext2; - - double value = 0; - if (inSum <= 1) { /* We're inside the pentachoron (4-Simplex) at (0,0,0,0) */ - - /* Determine which two of (0,0,0,1), (0,0,1,0), (0,1,0,0), (1,0,0,0) are closest. */ - aPoint = 0x01; - aScore = xins; - bPoint = 0x02; - bScore = yins; - if (aScore >= bScore && zins > bScore) { - bScore = zins; - bPoint = 0x04; - } else if (aScore < bScore && zins > aScore) { - aScore = zins; - aPoint = 0x04; - } - if (aScore >= bScore && wins > bScore) { - bScore = wins; - bPoint = 0x08; - } else if (aScore < bScore && wins > aScore) { - aScore = wins; - aPoint = 0x08; - } - - /* Now we determine the three lattice points not part of the pentachoron that may contribute. - This depends on the closest two pentachoron vertices, including (0,0,0,0) */ - uins = 1 - inSum; - if (uins > aScore || uins > bScore) { /* (0,0,0,0) is one of the closest two pentachoron vertices. */ - c = (bScore > aScore ? bPoint : aPoint); /* Our other closest vertex is the closest out of a and b. */ - if ((c & 0x01) == 0) { - xsv_ext0 = xsb - 1; - xsv_ext1 = xsv_ext2 = xsb; - dx_ext0 = dx0 + 1; - dx_ext1 = dx_ext2 = dx0; - } else { - xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb + 1; - dx_ext0 = dx_ext1 = dx_ext2 = dx0 - 1; - } - - if ((c & 0x02) == 0) { - ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; - dy_ext0 = dy_ext1 = dy_ext2 = dy0; - if ((c & 0x01) == 0x01) { - ysv_ext0 -= 1; - dy_ext0 += 1; - } else { - ysv_ext1 -= 1; - dy_ext1 += 1; - } - } else { - ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; - dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 1; - } - - if ((c & 0x04) == 0) { - zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; - dz_ext0 = dz_ext1 = dz_ext2 = dz0; - if ((c & 0x03) != 0) { - if ((c & 0x03) == 0x03) { - zsv_ext0 -= 1; - dz_ext0 += 1; - } else { - zsv_ext1 -= 1; - dz_ext1 += 1; - } - } else { - zsv_ext2 -= 1; - dz_ext2 += 1; - } - } else { - zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; - dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 1; - } - - if ((c & 0x08) == 0) { - wsv_ext0 = wsv_ext1 = wsb; - wsv_ext2 = wsb - 1; - dw_ext0 = dw_ext1 = dw0; - dw_ext2 = dw0 + 1; - } else { - wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb + 1; - dw_ext0 = dw_ext1 = dw_ext2 = dw0 - 1; - } - } else { /* (0,0,0,0) is not one of the closest two pentachoron vertices. */ - c = (int8_t)(aPoint | bPoint); /* Our three extra vertices are determined by the closest two. */ - - if ((c & 0x01) == 0) { - xsv_ext0 = xsv_ext2 = xsb; - xsv_ext1 = xsb - 1; - dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_4D; - dx_ext1 = dx0 + 1 - SQUISH_CONSTANT_4D; - dx_ext2 = dx0 - SQUISH_CONSTANT_4D; - } else { - xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb + 1; - dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - dx_ext1 = dx_ext2 = dx0 - 1 - SQUISH_CONSTANT_4D; - } - - if ((c & 0x02) == 0) { - ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; - dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_4D; - dy_ext1 = dy_ext2 = dy0 - SQUISH_CONSTANT_4D; - if ((c & 0x01) == 0x01) { - ysv_ext1 -= 1; - dy_ext1 += 1; - } else { - ysv_ext2 -= 1; - dy_ext2 += 1; - } - } else { - ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; - dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - dy_ext1 = dy_ext2 = dy0 - 1 - SQUISH_CONSTANT_4D; - } - - if ((c & 0x04) == 0) { - zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; - dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_4D; - dz_ext1 = dz_ext2 = dz0 - SQUISH_CONSTANT_4D; - if ((c & 0x03) == 0x03) { - zsv_ext1 -= 1; - dz_ext1 += 1; - } else { - zsv_ext2 -= 1; - dz_ext2 += 1; - } - } else { - zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; - dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - dz_ext1 = dz_ext2 = dz0 - 1 - SQUISH_CONSTANT_4D; - } - - if ((c & 0x08) == 0) { - wsv_ext0 = wsv_ext1 = wsb; - wsv_ext2 = wsb - 1; - dw_ext0 = dw0 - 2 * SQUISH_CONSTANT_4D; - dw_ext1 = dw0 - SQUISH_CONSTANT_4D; - dw_ext2 = dw0 + 1 - SQUISH_CONSTANT_4D; - } else { - wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb + 1; - dw_ext0 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - dw_ext1 = dw_ext2 = dw0 - 1 - SQUISH_CONSTANT_4D; - } - } - - /* Contribution (0,0,0,0) */ - attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0 - dw0 * dw0; - if (attn0 > 0) { - attn0 *= attn0; - value += attn0 * attn0 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 0, wsb + 0, dx0, dy0, dz0, dw0); - } - - /* Contribution (1,0,0,0) */ - dx1 = dx0 - 1 - SQUISH_CONSTANT_4D; - dy1 = dy0 - 0 - SQUISH_CONSTANT_4D; - dz1 = dz0 - 0 - SQUISH_CONSTANT_4D; - dw1 = dw0 - 0 - SQUISH_CONSTANT_4D; - attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; - if (attn1 > 0) { - attn1 *= attn1; - value += attn1 * attn1 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 0, wsb + 0, dx1, dy1, dz1, dw1); - } - - /* Contribution (0,1,0,0) */ - dx2 = dx0 - 0 - SQUISH_CONSTANT_4D; - dy2 = dy0 - 1 - SQUISH_CONSTANT_4D; - dz2 = dz1; - dw2 = dw1; - attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; - if (attn2 > 0) { - attn2 *= attn2; - value += attn2 * attn2 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 0, wsb + 0, dx2, dy2, dz2, dw2); - } - - /* Contribution (0,0,1,0) */ - dx3 = dx2; - dy3 = dy1; - dz3 = dz0 - 1 - SQUISH_CONSTANT_4D; - dw3 = dw1; - attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; - if (attn3 > 0) { - attn3 *= attn3; - value += attn3 * attn3 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 1, wsb + 0, dx3, dy3, dz3, dw3); - } - - /* Contribution (0,0,0,1) */ - dx4 = dx2; - dy4 = dy1; - dz4 = dz1; - dw4 = dw0 - 1 - SQUISH_CONSTANT_4D; - attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; - if (attn4 > 0) { - attn4 *= attn4; - value += attn4 * attn4 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 0, wsb + 1, dx4, dy4, dz4, dw4); - } - } else if (inSum >= 3) { /* We're inside the pentachoron (4-Simplex) at (1,1,1,1) - Determine which two of (1,1,1,0), (1,1,0,1), (1,0,1,1), (0,1,1,1) are closest. */ - aPoint = 0x0E; - aScore = xins; - bPoint = 0x0D; - bScore = yins; - if (aScore <= bScore && zins < bScore) { - bScore = zins; - bPoint = 0x0B; - } else if (aScore > bScore && zins < aScore) { - aScore = zins; - aPoint = 0x0B; - } - if (aScore <= bScore && wins < bScore) { - bScore = wins; - bPoint = 0x07; - } else if (aScore > bScore && wins < aScore) { - aScore = wins; - aPoint = 0x07; - } - - /* Now we determine the three lattice points not part of the pentachoron that may contribute. - This depends on the closest two pentachoron vertices, including (0,0,0,0) */ - uins = 4 - inSum; - if (uins < aScore || uins < bScore) { /* (1,1,1,1) is one of the closest two pentachoron vertices. */ - c = (bScore < aScore ? bPoint : aPoint); /* Our other closest vertex is the closest out of a and b. */ - - if ((c & 0x01) != 0) { - xsv_ext0 = xsb + 2; - xsv_ext1 = xsv_ext2 = xsb + 1; - dx_ext0 = dx0 - 2 - 4 * SQUISH_CONSTANT_4D; - dx_ext1 = dx_ext2 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; - } else { - xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb; - dx_ext0 = dx_ext1 = dx_ext2 = dx0 - 4 * SQUISH_CONSTANT_4D; - } - - if ((c & 0x02) != 0) { - ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; - dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; - if ((c & 0x01) != 0) { - ysv_ext1 += 1; - dy_ext1 -= 1; - } else { - ysv_ext0 += 1; - dy_ext0 -= 1; - } - } else { - ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; - dy_ext0 = dy_ext1 = dy_ext2 = dy0 - 4 * SQUISH_CONSTANT_4D; - } - - if ((c & 0x04) != 0) { - zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; - dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; - if ((c & 0x03) != 0x03) { - if ((c & 0x03) == 0) { - zsv_ext0 += 1; - dz_ext0 -= 1; - } else { - zsv_ext1 += 1; - dz_ext1 -= 1; - } - } else { - zsv_ext2 += 1; - dz_ext2 -= 1; - } - } else { - zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; - dz_ext0 = dz_ext1 = dz_ext2 = dz0 - 4 * SQUISH_CONSTANT_4D; - } - - if ((c & 0x08) != 0) { - wsv_ext0 = wsv_ext1 = wsb + 1; - wsv_ext2 = wsb + 2; - dw_ext0 = dw_ext1 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; - dw_ext2 = dw0 - 2 - 4 * SQUISH_CONSTANT_4D; - } else { - wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb; - dw_ext0 = dw_ext1 = dw_ext2 = dw0 - 4 * SQUISH_CONSTANT_4D; - } - } else { /* (1,1,1,1) is not one of the closest two pentachoron vertices. */ - c = (int8_t)(aPoint & bPoint); /* Our three extra vertices are determined by the closest two. */ - - if ((c & 0x01) != 0) { - xsv_ext0 = xsv_ext2 = xsb + 1; - xsv_ext1 = xsb + 2; - dx_ext0 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - dx_ext1 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; - dx_ext2 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; - } else { - xsv_ext0 = xsv_ext1 = xsv_ext2 = xsb; - dx_ext0 = dx0 - 2 * SQUISH_CONSTANT_4D; - dx_ext1 = dx_ext2 = dx0 - 3 * SQUISH_CONSTANT_4D; - } - - if ((c & 0x02) != 0) { - ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb + 1; - dy_ext0 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - dy_ext1 = dy_ext2 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; - if ((c & 0x01) != 0) { - ysv_ext2 += 1; - dy_ext2 -= 1; - } else { - ysv_ext1 += 1; - dy_ext1 -= 1; - } - } else { - ysv_ext0 = ysv_ext1 = ysv_ext2 = ysb; - dy_ext0 = dy0 - 2 * SQUISH_CONSTANT_4D; - dy_ext1 = dy_ext2 = dy0 - 3 * SQUISH_CONSTANT_4D; - } - - if ((c & 0x04) != 0) { - zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb + 1; - dz_ext0 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - dz_ext1 = dz_ext2 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; - if ((c & 0x03) != 0) { - zsv_ext2 += 1; - dz_ext2 -= 1; - } else { - zsv_ext1 += 1; - dz_ext1 -= 1; - } - } else { - zsv_ext0 = zsv_ext1 = zsv_ext2 = zsb; - dz_ext0 = dz0 - 2 * SQUISH_CONSTANT_4D; - dz_ext1 = dz_ext2 = dz0 - 3 * SQUISH_CONSTANT_4D; - } - - if ((c & 0x08) != 0) { - wsv_ext0 = wsv_ext1 = wsb + 1; - wsv_ext2 = wsb + 2; - dw_ext0 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - dw_ext1 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; - dw_ext2 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; - } else { - wsv_ext0 = wsv_ext1 = wsv_ext2 = wsb; - dw_ext0 = dw0 - 2 * SQUISH_CONSTANT_4D; - dw_ext1 = dw_ext2 = dw0 - 3 * SQUISH_CONSTANT_4D; - } - } - - /* Contribution (1,1,1,0) */ - dx4 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; - dy4 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; - dz4 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; - dw4 = dw0 - 3 * SQUISH_CONSTANT_4D; - attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; - if (attn4 > 0) { - attn4 *= attn4; - value += attn4 * attn4 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 1, wsb + 0, dx4, dy4, dz4, dw4); - } - - /* Contribution (1,1,0,1) */ - dx3 = dx4; - dy3 = dy4; - dz3 = dz0 - 3 * SQUISH_CONSTANT_4D; - dw3 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; - attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; - if (attn3 > 0) { - attn3 *= attn3; - value += attn3 * attn3 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 0, wsb + 1, dx3, dy3, dz3, dw3); - } - - /* Contribution (1,0,1,1) */ - dx2 = dx4; - dy2 = dy0 - 3 * SQUISH_CONSTANT_4D; - dz2 = dz4; - dw2 = dw3; - attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; - if (attn2 > 0) { - attn2 *= attn2; - value += attn2 * attn2 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 1, wsb + 1, dx2, dy2, dz2, dw2); - } - - /* Contribution (0,1,1,1) */ - dx1 = dx0 - 3 * SQUISH_CONSTANT_4D; - dz1 = dz4; - dy1 = dy4; - dw1 = dw3; - attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; - if (attn1 > 0) { - attn1 *= attn1; - value += attn1 * attn1 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 1, wsb + 1, dx1, dy1, dz1, dw1); - } - - /* Contribution (1,1,1,1) */ - dx0 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; - dy0 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; - dz0 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; - dw0 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; - attn0 = 2 - dx0 * dx0 - dy0 * dy0 - dz0 * dz0 - dw0 * dw0; - if (attn0 > 0) { - attn0 *= attn0; - value += attn0 * attn0 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 1, wsb + 1, dx0, dy0, dz0, dw0); - } - } else if (inSum <= 2) { /* We're inside the first dispentachoron (Rectified 4-Simplex) */ - aIsBiggerSide = 1; - bIsBiggerSide = 1; - - /* Decide between (1,1,0,0) and (0,0,1,1) */ - if (xins + yins > zins + wins) { - aScore = xins + yins; - aPoint = 0x03; - } else { - aScore = zins + wins; - aPoint = 0x0C; - } - - /* Decide between (1,0,1,0) and (0,1,0,1) */ - if (xins + zins > yins + wins) { - bScore = xins + zins; - bPoint = 0x05; - } else { - bScore = yins + wins; - bPoint = 0x0A; - } - - /* Closer between (1,0,0,1) and (0,1,1,0) will replace the further of a and b, if closer. */ - if (xins + wins > yins + zins) { - score = xins + wins; - if (aScore >= bScore && score > bScore) { - bScore = score; - bPoint = 0x09; - } else if (aScore < bScore && score > aScore) { - aScore = score; - aPoint = 0x09; - } - } else { - score = yins + zins; - if (aScore >= bScore && score > bScore) { - bScore = score; - bPoint = 0x06; - } else if (aScore < bScore && score > aScore) { - aScore = score; - aPoint = 0x06; - } - } - - /* Decide if (1,0,0,0) is closer. */ - p1 = 2 - inSum + xins; - if (aScore >= bScore && p1 > bScore) { - bScore = p1; - bPoint = 0x01; - bIsBiggerSide = 0; - } else if (aScore < bScore && p1 > aScore) { - aScore = p1; - aPoint = 0x01; - aIsBiggerSide = 0; - } - - /* Decide if (0,1,0,0) is closer. */ - p2 = 2 - inSum + yins; - if (aScore >= bScore && p2 > bScore) { - bScore = p2; - bPoint = 0x02; - bIsBiggerSide = 0; - } else if (aScore < bScore && p2 > aScore) { - aScore = p2; - aPoint = 0x02; - aIsBiggerSide = 0; - } - - /* Decide if (0,0,1,0) is closer. */ - p3 = 2 - inSum + zins; - if (aScore >= bScore && p3 > bScore) { - bScore = p3; - bPoint = 0x04; - bIsBiggerSide = 0; - } else if (aScore < bScore && p3 > aScore) { - aScore = p3; - aPoint = 0x04; - aIsBiggerSide = 0; - } - - /* Decide if (0,0,0,1) is closer. */ - p4 = 2 - inSum + wins; - if (aScore >= bScore && p4 > bScore) { - bScore = p4; - bPoint = 0x08; - bIsBiggerSide = 0; - } else if (aScore < bScore && p4 > aScore) { - aScore = p4; - aPoint = 0x08; - aIsBiggerSide = 0; - } - - /* Where each of the two closest points are determines how the extra three vertices are calculated. */ - if (aIsBiggerSide == bIsBiggerSide) { - if (aIsBiggerSide) { /* Both closest points on the bigger side */ - c1 = (int8_t)(aPoint | bPoint); - c2 = (int8_t)(aPoint & bPoint); - if ((c1 & 0x01) == 0) { - xsv_ext0 = xsb; - xsv_ext1 = xsb - 1; - dx_ext0 = dx0 - 3 * SQUISH_CONSTANT_4D; - dx_ext1 = dx0 + 1 - 2 * SQUISH_CONSTANT_4D; - } else { - xsv_ext0 = xsv_ext1 = xsb + 1; - dx_ext0 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; - dx_ext1 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - } - - if ((c1 & 0x02) == 0) { - ysv_ext0 = ysb; - ysv_ext1 = ysb - 1; - dy_ext0 = dy0 - 3 * SQUISH_CONSTANT_4D; - dy_ext1 = dy0 + 1 - 2 * SQUISH_CONSTANT_4D; - } else { - ysv_ext0 = ysv_ext1 = ysb + 1; - dy_ext0 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; - dy_ext1 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - } - - if ((c1 & 0x04) == 0) { - zsv_ext0 = zsb; - zsv_ext1 = zsb - 1; - dz_ext0 = dz0 - 3 * SQUISH_CONSTANT_4D; - dz_ext1 = dz0 + 1 - 2 * SQUISH_CONSTANT_4D; - } else { - zsv_ext0 = zsv_ext1 = zsb + 1; - dz_ext0 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; - dz_ext1 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - } - - if ((c1 & 0x08) == 0) { - wsv_ext0 = wsb; - wsv_ext1 = wsb - 1; - dw_ext0 = dw0 - 3 * SQUISH_CONSTANT_4D; - dw_ext1 = dw0 + 1 - 2 * SQUISH_CONSTANT_4D; - } else { - wsv_ext0 = wsv_ext1 = wsb + 1; - dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; - dw_ext1 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - } - - /* One combination is a permutation of (0,0,0,2) based on c2 */ - xsv_ext2 = xsb; - ysv_ext2 = ysb; - zsv_ext2 = zsb; - wsv_ext2 = wsb; - dx_ext2 = dx0 - 2 * SQUISH_CONSTANT_4D; - dy_ext2 = dy0 - 2 * SQUISH_CONSTANT_4D; - dz_ext2 = dz0 - 2 * SQUISH_CONSTANT_4D; - dw_ext2 = dw0 - 2 * SQUISH_CONSTANT_4D; - if ((c2 & 0x01) != 0) { - xsv_ext2 += 2; - dx_ext2 -= 2; - } else if ((c2 & 0x02) != 0) { - ysv_ext2 += 2; - dy_ext2 -= 2; - } else if ((c2 & 0x04) != 0) { - zsv_ext2 += 2; - dz_ext2 -= 2; - } else { - wsv_ext2 += 2; - dw_ext2 -= 2; - } - - } else { /* Both closest points on the smaller side */ - /* One of the two extra points is (0,0,0,0) */ - xsv_ext2 = xsb; - ysv_ext2 = ysb; - zsv_ext2 = zsb; - wsv_ext2 = wsb; - dx_ext2 = dx0; - dy_ext2 = dy0; - dz_ext2 = dz0; - dw_ext2 = dw0; - - /* Other two points are based on the omitted axes. */ - c = (int8_t)(aPoint | bPoint); - - if ((c & 0x01) == 0) { - xsv_ext0 = xsb - 1; - xsv_ext1 = xsb; - dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_4D; - dx_ext1 = dx0 - SQUISH_CONSTANT_4D; - } else { - xsv_ext0 = xsv_ext1 = xsb + 1; - dx_ext0 = dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_4D; - } - - if ((c & 0x02) == 0) { - ysv_ext0 = ysv_ext1 = ysb; - dy_ext0 = dy_ext1 = dy0 - SQUISH_CONSTANT_4D; - if ((c & 0x01) == 0x01) - { - ysv_ext0 -= 1; - dy_ext0 += 1; - } else { - ysv_ext1 -= 1; - dy_ext1 += 1; - } - } else { - ysv_ext0 = ysv_ext1 = ysb + 1; - dy_ext0 = dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_4D; - } - - if ((c & 0x04) == 0) { - zsv_ext0 = zsv_ext1 = zsb; - dz_ext0 = dz_ext1 = dz0 - SQUISH_CONSTANT_4D; - if ((c & 0x03) == 0x03) - { - zsv_ext0 -= 1; - dz_ext0 += 1; - } else { - zsv_ext1 -= 1; - dz_ext1 += 1; - } - } else { - zsv_ext0 = zsv_ext1 = zsb + 1; - dz_ext0 = dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_4D; - } - - if ((c & 0x08) == 0) - { - wsv_ext0 = wsb; - wsv_ext1 = wsb - 1; - dw_ext0 = dw0 - SQUISH_CONSTANT_4D; - dw_ext1 = dw0 + 1 - SQUISH_CONSTANT_4D; - } else { - wsv_ext0 = wsv_ext1 = wsb + 1; - dw_ext0 = dw_ext1 = dw0 - 1 - SQUISH_CONSTANT_4D; - } - - } - } else { /* One point on each "side" */ - if (aIsBiggerSide) { - c1 = aPoint; - c2 = bPoint; - } else { - c1 = bPoint; - c2 = aPoint; - } - - /* Two contributions are the bigger-sided point with each 0 replaced with -1. */ - if ((c1 & 0x01) == 0) { - xsv_ext0 = xsb - 1; - xsv_ext1 = xsb; - dx_ext0 = dx0 + 1 - SQUISH_CONSTANT_4D; - dx_ext1 = dx0 - SQUISH_CONSTANT_4D; - } else { - xsv_ext0 = xsv_ext1 = xsb + 1; - dx_ext0 = dx_ext1 = dx0 - 1 - SQUISH_CONSTANT_4D; - } - - if ((c1 & 0x02) == 0) { - ysv_ext0 = ysv_ext1 = ysb; - dy_ext0 = dy_ext1 = dy0 - SQUISH_CONSTANT_4D; - if ((c1 & 0x01) == 0x01) { - ysv_ext0 -= 1; - dy_ext0 += 1; - } else { - ysv_ext1 -= 1; - dy_ext1 += 1; - } - } else { - ysv_ext0 = ysv_ext1 = ysb + 1; - dy_ext0 = dy_ext1 = dy0 - 1 - SQUISH_CONSTANT_4D; - } - - if ((c1 & 0x04) == 0) { - zsv_ext0 = zsv_ext1 = zsb; - dz_ext0 = dz_ext1 = dz0 - SQUISH_CONSTANT_4D; - if ((c1 & 0x03) == 0x03) { - zsv_ext0 -= 1; - dz_ext0 += 1; - } else { - zsv_ext1 -= 1; - dz_ext1 += 1; - } - } else { - zsv_ext0 = zsv_ext1 = zsb + 1; - dz_ext0 = dz_ext1 = dz0 - 1 - SQUISH_CONSTANT_4D; - } - - if ((c1 & 0x08) == 0) { - wsv_ext0 = wsb; - wsv_ext1 = wsb - 1; - dw_ext0 = dw0 - SQUISH_CONSTANT_4D; - dw_ext1 = dw0 + 1 - SQUISH_CONSTANT_4D; - } else { - wsv_ext0 = wsv_ext1 = wsb + 1; - dw_ext0 = dw_ext1 = dw0 - 1 - SQUISH_CONSTANT_4D; - } - - /* One contribution is a permutation of (0,0,0,2) based on the smaller-sided point */ - xsv_ext2 = xsb; - ysv_ext2 = ysb; - zsv_ext2 = zsb; - wsv_ext2 = wsb; - dx_ext2 = dx0 - 2 * SQUISH_CONSTANT_4D; - dy_ext2 = dy0 - 2 * SQUISH_CONSTANT_4D; - dz_ext2 = dz0 - 2 * SQUISH_CONSTANT_4D; - dw_ext2 = dw0 - 2 * SQUISH_CONSTANT_4D; - if ((c2 & 0x01) != 0) { - xsv_ext2 += 2; - dx_ext2 -= 2; - } else if ((c2 & 0x02) != 0) { - ysv_ext2 += 2; - dy_ext2 -= 2; - } else if ((c2 & 0x04) != 0) { - zsv_ext2 += 2; - dz_ext2 -= 2; - } else { - wsv_ext2 += 2; - dw_ext2 -= 2; - } - } - - /* Contribution (1,0,0,0) */ - dx1 = dx0 - 1 - SQUISH_CONSTANT_4D; - dy1 = dy0 - 0 - SQUISH_CONSTANT_4D; - dz1 = dz0 - 0 - SQUISH_CONSTANT_4D; - dw1 = dw0 - 0 - SQUISH_CONSTANT_4D; - attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; - if (attn1 > 0) { - attn1 *= attn1; - value += attn1 * attn1 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 0, wsb + 0, dx1, dy1, dz1, dw1); - } - - /* Contribution (0,1,0,0) */ - dx2 = dx0 - 0 - SQUISH_CONSTANT_4D; - dy2 = dy0 - 1 - SQUISH_CONSTANT_4D; - dz2 = dz1; - dw2 = dw1; - attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; - if (attn2 > 0) { - attn2 *= attn2; - value += attn2 * attn2 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 0, wsb + 0, dx2, dy2, dz2, dw2); - } - - /* Contribution (0,0,1,0) */ - dx3 = dx2; - dy3 = dy1; - dz3 = dz0 - 1 - SQUISH_CONSTANT_4D; - dw3 = dw1; - attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; - if (attn3 > 0) { - attn3 *= attn3; - value += attn3 * attn3 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 1, wsb + 0, dx3, dy3, dz3, dw3); - } - - /* Contribution (0,0,0,1) */ - dx4 = dx2; - dy4 = dy1; - dz4 = dz1; - dw4 = dw0 - 1 - SQUISH_CONSTANT_4D; - attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; - if (attn4 > 0) { - attn4 *= attn4; - value += attn4 * attn4 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 0, wsb + 1, dx4, dy4, dz4, dw4); - } - - /* Contribution (1,1,0,0) */ - dx5 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - dy5 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - dz5 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; - dw5 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; - attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5 - dw5 * dw5; - if (attn5 > 0) { - attn5 *= attn5; - value += attn5 * attn5 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 0, wsb + 0, dx5, dy5, dz5, dw5); - } - - /* Contribution (1,0,1,0) */ - dx6 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - dy6 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; - dz6 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - dw6 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; - attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6 - dw6 * dw6; - if (attn6 > 0) { - attn6 *= attn6; - value += attn6 * attn6 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 1, wsb + 0, dx6, dy6, dz6, dw6); - } - - /* Contribution (1,0,0,1) */ - dx7 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - dy7 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; - dz7 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; - dw7 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - attn7 = 2 - dx7 * dx7 - dy7 * dy7 - dz7 * dz7 - dw7 * dw7; - if (attn7 > 0) { - attn7 *= attn7; - value += attn7 * attn7 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 0, wsb + 1, dx7, dy7, dz7, dw7); - } - - /* Contribution (0,1,1,0) */ - dx8 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; - dy8 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - dz8 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - dw8 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; - attn8 = 2 - dx8 * dx8 - dy8 * dy8 - dz8 * dz8 - dw8 * dw8; - if (attn8 > 0) { - attn8 *= attn8; - value += attn8 * attn8 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 1, wsb + 0, dx8, dy8, dz8, dw8); - } - - /* Contribution (0,1,0,1) */ - dx9 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; - dy9 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - dz9 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; - dw9 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - attn9 = 2 - dx9 * dx9 - dy9 * dy9 - dz9 * dz9 - dw9 * dw9; - if (attn9 > 0) { - attn9 *= attn9; - value += attn9 * attn9 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 0, wsb + 1, dx9, dy9, dz9, dw9); - } - - /* Contribution (0,0,1,1) */ - dx10 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; - dy10 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; - dz10 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - dw10 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - attn10 = 2 - dx10 * dx10 - dy10 * dy10 - dz10 * dz10 - dw10 * dw10; - if (attn10 > 0) { - attn10 *= attn10; - value += attn10 * attn10 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 1, wsb + 1, dx10, dy10, dz10, dw10); - } - } else { /* We're inside the second dispentachoron (Rectified 4-Simplex) */ - aIsBiggerSide = 1; - bIsBiggerSide = 1; - - /* Decide between (0,0,1,1) and (1,1,0,0) */ - if (xins + yins < zins + wins) { - aScore = xins + yins; - aPoint = 0x0C; - } else { - aScore = zins + wins; - aPoint = 0x03; - } - - /* Decide between (0,1,0,1) and (1,0,1,0) */ - if (xins + zins < yins + wins) { - bScore = xins + zins; - bPoint = 0x0A; - } else { - bScore = yins + wins; - bPoint = 0x05; - } - - /* Closer between (0,1,1,0) and (1,0,0,1) will replace the further of a and b, if closer. */ - if (xins + wins < yins + zins) { - score = xins + wins; - if (aScore <= bScore && score < bScore) { - bScore = score; - bPoint = 0x06; - } else if (aScore > bScore && score < aScore) { - aScore = score; - aPoint = 0x06; - } - } else { - score = yins + zins; - if (aScore <= bScore && score < bScore) { - bScore = score; - bPoint = 0x09; - } else if (aScore > bScore && score < aScore) { - aScore = score; - aPoint = 0x09; - } - } - - /* Decide if (0,1,1,1) is closer. */ - p1 = 3 - inSum + xins; - if (aScore <= bScore && p1 < bScore) { - bScore = p1; - bPoint = 0x0E; - bIsBiggerSide = 0; - } else if (aScore > bScore && p1 < aScore) { - aScore = p1; - aPoint = 0x0E; - aIsBiggerSide = 0; - } - - /* Decide if (1,0,1,1) is closer. */ - p2 = 3 - inSum + yins; - if (aScore <= bScore && p2 < bScore) { - bScore = p2; - bPoint = 0x0D; - bIsBiggerSide = 0; - } else if (aScore > bScore && p2 < aScore) { - aScore = p2; - aPoint = 0x0D; - aIsBiggerSide = 0; - } - - /* Decide if (1,1,0,1) is closer. */ - p3 = 3 - inSum + zins; - if (aScore <= bScore && p3 < bScore) { - bScore = p3; - bPoint = 0x0B; - bIsBiggerSide = 0; - } else if (aScore > bScore && p3 < aScore) { - aScore = p3; - aPoint = 0x0B; - aIsBiggerSide = 0; - } - - /* Decide if (1,1,1,0) is closer. */ - p4 = 3 - inSum + wins; - if (aScore <= bScore && p4 < bScore) { - bScore = p4; - bPoint = 0x07; - bIsBiggerSide = 0; - } else if (aScore > bScore && p4 < aScore) { - aScore = p4; - aPoint = 0x07; - aIsBiggerSide = 0; - } - - /* Where each of the two closest points are determines how the extra three vertices are calculated. */ - if (aIsBiggerSide == bIsBiggerSide) { - if (aIsBiggerSide) { /* Both closest points on the bigger side */ - c1 = (int8_t)(aPoint & bPoint); - c2 = (int8_t)(aPoint | bPoint); - - /* Two contributions are permutations of (0,0,0,1) and (0,0,0,2) based on c1 */ - xsv_ext0 = xsv_ext1 = xsb; - ysv_ext0 = ysv_ext1 = ysb; - zsv_ext0 = zsv_ext1 = zsb; - wsv_ext0 = wsv_ext1 = wsb; - dx_ext0 = dx0 - SQUISH_CONSTANT_4D; - dy_ext0 = dy0 - SQUISH_CONSTANT_4D; - dz_ext0 = dz0 - SQUISH_CONSTANT_4D; - dw_ext0 = dw0 - SQUISH_CONSTANT_4D; - dx_ext1 = dx0 - 2 * SQUISH_CONSTANT_4D; - dy_ext1 = dy0 - 2 * SQUISH_CONSTANT_4D; - dz_ext1 = dz0 - 2 * SQUISH_CONSTANT_4D; - dw_ext1 = dw0 - 2 * SQUISH_CONSTANT_4D; - if ((c1 & 0x01) != 0) { - xsv_ext0 += 1; - dx_ext0 -= 1; - xsv_ext1 += 2; - dx_ext1 -= 2; - } else if ((c1 & 0x02) != 0) { - ysv_ext0 += 1; - dy_ext0 -= 1; - ysv_ext1 += 2; - dy_ext1 -= 2; - } else if ((c1 & 0x04) != 0) { - zsv_ext0 += 1; - dz_ext0 -= 1; - zsv_ext1 += 2; - dz_ext1 -= 2; - } else { - wsv_ext0 += 1; - dw_ext0 -= 1; - wsv_ext1 += 2; - dw_ext1 -= 2; - } - - /* One contribution is a permutation of (1,1,1,-1) based on c2 */ - xsv_ext2 = xsb + 1; - ysv_ext2 = ysb + 1; - zsv_ext2 = zsb + 1; - wsv_ext2 = wsb + 1; - dx_ext2 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - dy_ext2 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - dz_ext2 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - dw_ext2 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - if ((c2 & 0x01) == 0) { - xsv_ext2 -= 2; - dx_ext2 += 2; - } else if ((c2 & 0x02) == 0) { - ysv_ext2 -= 2; - dy_ext2 += 2; - } else if ((c2 & 0x04) == 0) { - zsv_ext2 -= 2; - dz_ext2 += 2; - } else { - wsv_ext2 -= 2; - dw_ext2 += 2; - } - } else { /* Both closest points on the smaller side */ - /* One of the two extra points is (1,1,1,1) */ - xsv_ext2 = xsb + 1; - ysv_ext2 = ysb + 1; - zsv_ext2 = zsb + 1; - wsv_ext2 = wsb + 1; - dx_ext2 = dx0 - 1 - 4 * SQUISH_CONSTANT_4D; - dy_ext2 = dy0 - 1 - 4 * SQUISH_CONSTANT_4D; - dz_ext2 = dz0 - 1 - 4 * SQUISH_CONSTANT_4D; - dw_ext2 = dw0 - 1 - 4 * SQUISH_CONSTANT_4D; - - /* Other two points are based on the shared axes. */ - c = (int8_t)(aPoint & bPoint); - - if ((c & 0x01) != 0) { - xsv_ext0 = xsb + 2; - xsv_ext1 = xsb + 1; - dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; - dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; - } else { - xsv_ext0 = xsv_ext1 = xsb; - dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_4D; - } - - if ((c & 0x02) != 0) { - ysv_ext0 = ysv_ext1 = ysb + 1; - dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; - if ((c & 0x01) == 0) - { - ysv_ext0 += 1; - dy_ext0 -= 1; - } else { - ysv_ext1 += 1; - dy_ext1 -= 1; - } - } else { - ysv_ext0 = ysv_ext1 = ysb; - dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_4D; - } - - if ((c & 0x04) != 0) { - zsv_ext0 = zsv_ext1 = zsb + 1; - dz_ext0 = dz_ext1 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; - if ((c & 0x03) == 0) - { - zsv_ext0 += 1; - dz_ext0 -= 1; - } else { - zsv_ext1 += 1; - dz_ext1 -= 1; - } - } else { - zsv_ext0 = zsv_ext1 = zsb; - dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_4D; - } - - if ((c & 0x08) != 0) - { - wsv_ext0 = wsb + 1; - wsv_ext1 = wsb + 2; - dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; - dw_ext1 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; - } else { - wsv_ext0 = wsv_ext1 = wsb; - dw_ext0 = dw_ext1 = dw0 - 3 * SQUISH_CONSTANT_4D; - } - } - } else { /* One point on each "side" */ - if (aIsBiggerSide) { - c1 = aPoint; - c2 = bPoint; - } else { - c1 = bPoint; - c2 = aPoint; - } - - /* Two contributions are the bigger-sided point with each 1 replaced with 2. */ - if ((c1 & 0x01) != 0) { - xsv_ext0 = xsb + 2; - xsv_ext1 = xsb + 1; - dx_ext0 = dx0 - 2 - 3 * SQUISH_CONSTANT_4D; - dx_ext1 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; - } else { - xsv_ext0 = xsv_ext1 = xsb; - dx_ext0 = dx_ext1 = dx0 - 3 * SQUISH_CONSTANT_4D; - } - - if ((c1 & 0x02) != 0) { - ysv_ext0 = ysv_ext1 = ysb + 1; - dy_ext0 = dy_ext1 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; - if ((c1 & 0x01) == 0) { - ysv_ext0 += 1; - dy_ext0 -= 1; - } else { - ysv_ext1 += 1; - dy_ext1 -= 1; - } - } else { - ysv_ext0 = ysv_ext1 = ysb; - dy_ext0 = dy_ext1 = dy0 - 3 * SQUISH_CONSTANT_4D; - } - - if ((c1 & 0x04) != 0) { - zsv_ext0 = zsv_ext1 = zsb + 1; - dz_ext0 = dz_ext1 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; - if ((c1 & 0x03) == 0) { - zsv_ext0 += 1; - dz_ext0 -= 1; - } else { - zsv_ext1 += 1; - dz_ext1 -= 1; - } - } else { - zsv_ext0 = zsv_ext1 = zsb; - dz_ext0 = dz_ext1 = dz0 - 3 * SQUISH_CONSTANT_4D; - } - - if ((c1 & 0x08) != 0) { - wsv_ext0 = wsb + 1; - wsv_ext1 = wsb + 2; - dw_ext0 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; - dw_ext1 = dw0 - 2 - 3 * SQUISH_CONSTANT_4D; - } else { - wsv_ext0 = wsv_ext1 = wsb; - dw_ext0 = dw_ext1 = dw0 - 3 * SQUISH_CONSTANT_4D; - } - - /* One contribution is a permutation of (1,1,1,-1) based on the smaller-sided point */ - xsv_ext2 = xsb + 1; - ysv_ext2 = ysb + 1; - zsv_ext2 = zsb + 1; - wsv_ext2 = wsb + 1; - dx_ext2 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - dy_ext2 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - dz_ext2 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - dw_ext2 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - if ((c2 & 0x01) == 0) { - xsv_ext2 -= 2; - dx_ext2 += 2; - } else if ((c2 & 0x02) == 0) { - ysv_ext2 -= 2; - dy_ext2 += 2; - } else if ((c2 & 0x04) == 0) { - zsv_ext2 -= 2; - dz_ext2 += 2; - } else { - wsv_ext2 -= 2; - dw_ext2 += 2; - } - } - - /* Contribution (1,1,1,0) */ - dx4 = dx0 - 1 - 3 * SQUISH_CONSTANT_4D; - dy4 = dy0 - 1 - 3 * SQUISH_CONSTANT_4D; - dz4 = dz0 - 1 - 3 * SQUISH_CONSTANT_4D; - dw4 = dw0 - 3 * SQUISH_CONSTANT_4D; - attn4 = 2 - dx4 * dx4 - dy4 * dy4 - dz4 * dz4 - dw4 * dw4; - if (attn4 > 0) { - attn4 *= attn4; - value += attn4 * attn4 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 1, wsb + 0, dx4, dy4, dz4, dw4); - } - - /* Contribution (1,1,0,1) */ - dx3 = dx4; - dy3 = dy4; - dz3 = dz0 - 3 * SQUISH_CONSTANT_4D; - dw3 = dw0 - 1 - 3 * SQUISH_CONSTANT_4D; - attn3 = 2 - dx3 * dx3 - dy3 * dy3 - dz3 * dz3 - dw3 * dw3; - if (attn3 > 0) { - attn3 *= attn3; - value += attn3 * attn3 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 0, wsb + 1, dx3, dy3, dz3, dw3); - } - - /* Contribution (1,0,1,1) */ - dx2 = dx4; - dy2 = dy0 - 3 * SQUISH_CONSTANT_4D; - dz2 = dz4; - dw2 = dw3; - attn2 = 2 - dx2 * dx2 - dy2 * dy2 - dz2 * dz2 - dw2 * dw2; - if (attn2 > 0) { - attn2 *= attn2; - value += attn2 * attn2 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 1, wsb + 1, dx2, dy2, dz2, dw2); - } - - /* Contribution (0,1,1,1) */ - dx1 = dx0 - 3 * SQUISH_CONSTANT_4D; - dz1 = dz4; - dy1 = dy4; - dw1 = dw3; - attn1 = 2 - dx1 * dx1 - dy1 * dy1 - dz1 * dz1 - dw1 * dw1; - if (attn1 > 0) { - attn1 *= attn1; - value += attn1 * attn1 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 1, wsb + 1, dx1, dy1, dz1, dw1); - } - - /* Contribution (1,1,0,0) */ - dx5 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - dy5 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - dz5 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; - dw5 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; - attn5 = 2 - dx5 * dx5 - dy5 * dy5 - dz5 * dz5 - dw5 * dw5; - if (attn5 > 0) { - attn5 *= attn5; - value += attn5 * attn5 * extrapolate4(ctx, xsb + 1, ysb + 1, zsb + 0, wsb + 0, dx5, dy5, dz5, dw5); - } - - /* Contribution (1,0,1,0) */ - dx6 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - dy6 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; - dz6 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - dw6 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; - attn6 = 2 - dx6 * dx6 - dy6 * dy6 - dz6 * dz6 - dw6 * dw6; - if (attn6 > 0) { - attn6 *= attn6; - value += attn6 * attn6 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 1, wsb + 0, dx6, dy6, dz6, dw6); - } - - /* Contribution (1,0,0,1) */ - dx7 = dx0 - 1 - 2 * SQUISH_CONSTANT_4D; - dy7 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; - dz7 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; - dw7 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - attn7 = 2 - dx7 * dx7 - dy7 * dy7 - dz7 * dz7 - dw7 * dw7; - if (attn7 > 0) { - attn7 *= attn7; - value += attn7 * attn7 * extrapolate4(ctx, xsb + 1, ysb + 0, zsb + 0, wsb + 1, dx7, dy7, dz7, dw7); - } - - /* Contribution (0,1,1,0) */ - dx8 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; - dy8 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - dz8 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - dw8 = dw0 - 0 - 2 * SQUISH_CONSTANT_4D; - attn8 = 2 - dx8 * dx8 - dy8 * dy8 - dz8 * dz8 - dw8 * dw8; - if (attn8 > 0) { - attn8 *= attn8; - value += attn8 * attn8 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 1, wsb + 0, dx8, dy8, dz8, dw8); - } - - /* Contribution (0,1,0,1) */ - dx9 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; - dy9 = dy0 - 1 - 2 * SQUISH_CONSTANT_4D; - dz9 = dz0 - 0 - 2 * SQUISH_CONSTANT_4D; - dw9 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - attn9 = 2 - dx9 * dx9 - dy9 * dy9 - dz9 * dz9 - dw9 * dw9; - if (attn9 > 0) { - attn9 *= attn9; - value += attn9 * attn9 * extrapolate4(ctx, xsb + 0, ysb + 1, zsb + 0, wsb + 1, dx9, dy9, dz9, dw9); - } - - /* Contribution (0,0,1,1) */ - dx10 = dx0 - 0 - 2 * SQUISH_CONSTANT_4D; - dy10 = dy0 - 0 - 2 * SQUISH_CONSTANT_4D; - dz10 = dz0 - 1 - 2 * SQUISH_CONSTANT_4D; - dw10 = dw0 - 1 - 2 * SQUISH_CONSTANT_4D; - attn10 = 2 - dx10 * dx10 - dy10 * dy10 - dz10 * dz10 - dw10 * dw10; - if (attn10 > 0) { - attn10 *= attn10; - value += attn10 * attn10 * extrapolate4(ctx, xsb + 0, ysb + 0, zsb + 1, wsb + 1, dx10, dy10, dz10, dw10); - } - } - - /* First extra vertex */ - attn_ext0 = 2 - dx_ext0 * dx_ext0 - dy_ext0 * dy_ext0 - dz_ext0 * dz_ext0 - dw_ext0 * dw_ext0; - if (attn_ext0 > 0) - { - attn_ext0 *= attn_ext0; - value += attn_ext0 * attn_ext0 * extrapolate4(ctx, xsv_ext0, ysv_ext0, zsv_ext0, wsv_ext0, dx_ext0, dy_ext0, dz_ext0, dw_ext0); - } - - /* Second extra vertex */ - attn_ext1 = 2 - dx_ext1 * dx_ext1 - dy_ext1 * dy_ext1 - dz_ext1 * dz_ext1 - dw_ext1 * dw_ext1; - if (attn_ext1 > 0) - { - attn_ext1 *= attn_ext1; - value += attn_ext1 * attn_ext1 * extrapolate4(ctx, xsv_ext1, ysv_ext1, zsv_ext1, wsv_ext1, dx_ext1, dy_ext1, dz_ext1, dw_ext1); - } - - /* Third extra vertex */ - attn_ext2 = 2 - dx_ext2 * dx_ext2 - dy_ext2 * dy_ext2 - dz_ext2 * dz_ext2 - dw_ext2 * dw_ext2; - if (attn_ext2 > 0) - { - attn_ext2 *= attn_ext2; - value += attn_ext2 * attn_ext2 * extrapolate4(ctx, xsv_ext2, ysv_ext2, zsv_ext2, wsv_ext2, dx_ext2, dy_ext2, dz_ext2, dw_ext2); - } - - return value / NORM_CONSTANT_4D; -} - diff --git a/thirdparty/misc/open-simplex-noise.h b/thirdparty/misc/open-simplex-noise.h deleted file mode 100644 index fd9248c3a1..0000000000 --- a/thirdparty/misc/open-simplex-noise.h +++ /dev/null @@ -1,58 +0,0 @@ -#ifndef OPEN_SIMPLEX_NOISE_H__ -#define OPEN_SIMPLEX_NOISE_H__ - -/* - * OpenSimplex (Simplectic) Noise in C. - * Ported to C from Kurt Spencer's java implementation by Stephen M. Cameron - * - * v1.1 (October 6, 2014) - * - Ported to C - * - * v1.1 (October 5, 2014) - * - Added 2D and 4D implementations. - * - Proper gradient sets for all dimensions, from a - * dimensionally-generalizable scheme with an actual - * rhyme and reason behind it. - * - Removed default permutation array in favor of - * default seed. - * - Changed seed-based constructor to be independent - * of any particular randomization library, so results - * will be the same when ported to other languages. - */ - -#if ((__GNUC_STDC_INLINE__) || (__STDC_VERSION__ >= 199901L)) - #include <stdint.h> - #define INLINE inline -#elif (defined (_MSC_VER) || defined (__GNUC_GNU_INLINE__)) - #include <stdint.h> - #define INLINE __inline -#else - /* ANSI C doesn't have inline or stdint.h. */ - #define INLINE -#endif - -#ifdef __cplusplus - extern "C" { -#endif - -// -- GODOT start -- -// Modified to work without allocating memory, also removed some unused function. - -struct osn_context { - int16_t perm[256]; - int16_t permGradIndex3D[256]; -}; - -int open_simplex_noise(int64_t seed, struct osn_context *ctx); -//int open_simplex_noise_init_perm(struct osn_context *ctx, int16_t p[], int nelements); -// -- GODOT end -- -void open_simplex_noise_free(struct osn_context *ctx); -double open_simplex_noise2(const struct osn_context *ctx, double x, double y); -double open_simplex_noise3(const struct osn_context *ctx, double x, double y, double z); -double open_simplex_noise4(const struct osn_context *ctx, double x, double y, double z, double w); - -#ifdef __cplusplus - } -#endif - -#endif diff --git a/thirdparty/misc/patches/polypartition-godot-types.patch b/thirdparty/misc/patches/polypartition-godot-types.patch index 782f02e8dc..5d8aba3437 100644 --- a/thirdparty/misc/patches/polypartition-godot-types.patch +++ b/thirdparty/misc/patches/polypartition-godot-types.patch @@ -1,19 +1,16 @@ diff --git a/thirdparty/misc/polypartition.cpp b/thirdparty/misc/polypartition.cpp -index 3a8a6efa83..5e94793b79 100644 +index 3a8a6efa83..8c5409bf24 100644 --- a/thirdparty/misc/polypartition.cpp +++ b/thirdparty/misc/polypartition.cpp -@@ -23,10 +23,7 @@ - - #include "polypartition.h" - --#include <math.h> --#include <string.h> +@@ -26,7 +26,6 @@ + #include <math.h> + #include <string.h> #include <algorithm> -#include <vector> TPPLPoly::TPPLPoly() { hole = false; -@@ -186,7 +183,7 @@ int TPPLPartition::Intersects(TPPLPoint &p11, TPPLPoint &p12, TPPLPoint &p21, TP +@@ -186,7 +185,7 @@ int TPPLPartition::Intersects(TPPLPoint &p11, TPPLPoint &p12, TPPLPoint &p21, TP // Removes holes from inpolys by merging them with non-holes. int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { TPPLPolyList polys; @@ -22,7 +19,7 @@ index 3a8a6efa83..5e94793b79 100644 long i, i2, holepointindex, polypointindex; TPPLPoint holepoint, polypoint, bestpolypoint; TPPLPoint linep1, linep2; -@@ -198,15 +195,15 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { +@@ -198,15 +197,15 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { // Check for the trivial case of no holes. hasholes = false; @@ -42,7 +39,7 @@ index 3a8a6efa83..5e94793b79 100644 } return 1; } -@@ -216,8 +213,8 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { +@@ -216,8 +215,8 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { while (1) { // Find the hole point with the largest x. hasholes = false; @@ -53,7 +50,7 @@ index 3a8a6efa83..5e94793b79 100644 continue; } -@@ -227,8 +224,8 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { +@@ -227,8 +226,8 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { holepointindex = 0; } @@ -64,7 +61,7 @@ index 3a8a6efa83..5e94793b79 100644 holeiter = iter; holepointindex = i; } -@@ -237,24 +234,24 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { +@@ -237,24 +236,24 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { if (!hasholes) { break; } @@ -98,13 +95,13 @@ index 3a8a6efa83..5e94793b79 100644 if (pointfound) { v1 = Normalize(polypoint - holepoint); v2 = Normalize(bestpolypoint - holepoint); -@@ -263,13 +260,13 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { +@@ -263,13 +262,13 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { } } pointvisible = true; - for (iter2 = polys.begin(); iter2 != polys.end(); iter2++) { - if (iter2->IsHole()) { -+ for (iter2 = polys.front(); iter2; iter2->next()) { ++ for (iter2 = polys.front(); iter2; iter2 = iter2->next()) { + if (iter2->get().IsHole()) { continue; } @@ -117,7 +114,7 @@ index 3a8a6efa83..5e94793b79 100644 if (Intersects(holepoint, polypoint, linep1, linep2)) { pointvisible = false; break; -@@ -292,18 +289,18 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { +@@ -292,18 +291,18 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { return 0; } @@ -142,7 +139,7 @@ index 3a8a6efa83..5e94793b79 100644 i2++; } -@@ -312,8 +309,8 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { +@@ -312,8 +311,8 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { polys.push_back(newpoly); } @@ -153,7 +150,7 @@ index 3a8a6efa83..5e94793b79 100644 } return 1; -@@ -524,13 +521,13 @@ int TPPLPartition::Triangulate_EC(TPPLPoly *poly, TPPLPolyList *triangles) { +@@ -524,13 +523,13 @@ int TPPLPartition::Triangulate_EC(TPPLPoly *poly, TPPLPolyList *triangles) { int TPPLPartition::Triangulate_EC(TPPLPolyList *inpolys, TPPLPolyList *triangles) { TPPLPolyList outpolys; @@ -170,7 +167,7 @@ index 3a8a6efa83..5e94793b79 100644 return 0; } } -@@ -543,7 +540,7 @@ int TPPLPartition::ConvexPartition_HM(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -543,7 +542,7 @@ int TPPLPartition::ConvexPartition_HM(TPPLPoly *poly, TPPLPolyList *parts) { } TPPLPolyList triangles; @@ -179,7 +176,7 @@ index 3a8a6efa83..5e94793b79 100644 TPPLPoly *poly1 = NULL, *poly2 = NULL; TPPLPoly newpoly; TPPLPoint d1, d2, p1, p2, p3; -@@ -578,19 +575,19 @@ int TPPLPartition::ConvexPartition_HM(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -578,19 +577,19 @@ int TPPLPartition::ConvexPartition_HM(TPPLPoly *poly, TPPLPolyList *parts) { return 0; } @@ -203,7 +200,7 @@ index 3a8a6efa83..5e94793b79 100644 for (i21 = 0; i21 < poly2->GetNumPoints(); i21++) { if ((d2.x != poly2->GetPoint(i21).x) || (d2.y != poly2->GetPoint(i21).y)) { -@@ -660,16 +657,16 @@ int TPPLPartition::ConvexPartition_HM(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -660,16 +659,16 @@ int TPPLPartition::ConvexPartition_HM(TPPLPoly *poly, TPPLPolyList *parts) { } triangles.erase(iter2); @@ -224,7 +221,7 @@ index 3a8a6efa83..5e94793b79 100644 } return 1; -@@ -677,13 +674,13 @@ int TPPLPartition::ConvexPartition_HM(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -677,13 +676,13 @@ int TPPLPartition::ConvexPartition_HM(TPPLPoly *poly, TPPLPolyList *parts) { int TPPLPartition::ConvexPartition_HM(TPPLPolyList *inpolys, TPPLPolyList *parts) { TPPLPolyList outpolys; @@ -241,7 +238,7 @@ index 3a8a6efa83..5e94793b79 100644 return 0; } } -@@ -824,8 +821,8 @@ int TPPLPartition::Triangulate_OPT(TPPLPoly *poly, TPPLPolyList *triangles) { +@@ -824,8 +823,8 @@ int TPPLPartition::Triangulate_OPT(TPPLPoly *poly, TPPLPolyList *triangles) { newdiagonal.index1 = 0; newdiagonal.index2 = n - 1; diagonals.push_back(newdiagonal); @@ -252,7 +249,7 @@ index 3a8a6efa83..5e94793b79 100644 diagonals.pop_front(); bestvertex = dpstates[diagonal.index2][diagonal.index1].bestvertex; if (bestvertex == -1) { -@@ -873,10 +870,10 @@ void TPPLPartition::UpdateState(long a, long b, long w, long i, long j, DPState2 +@@ -873,10 +872,10 @@ void TPPLPartition::UpdateState(long a, long b, long w, long i, long j, DPState2 pairs->push_front(newdiagonal); dpstates[a][b].weight = w; } else { @@ -265,7 +262,7 @@ index 3a8a6efa83..5e94793b79 100644 pairs->pop_front(); } pairs->push_front(newdiagonal); -@@ -885,7 +882,7 @@ void TPPLPartition::UpdateState(long a, long b, long w, long i, long j, DPState2 +@@ -885,7 +884,7 @@ void TPPLPartition::UpdateState(long a, long b, long w, long i, long j, DPState2 void TPPLPartition::TypeA(long i, long j, long k, PartitionVertex *vertices, DPState2 **dpstates) { DiagonalList *pairs = NULL; @@ -274,7 +271,7 @@ index 3a8a6efa83..5e94793b79 100644 long top; long w; -@@ -902,23 +899,23 @@ void TPPLPartition::TypeA(long i, long j, long k, PartitionVertex *vertices, DPS +@@ -902,23 +901,23 @@ void TPPLPartition::TypeA(long i, long j, long k, PartitionVertex *vertices, DPS } if (j - i > 1) { pairs = &(dpstates[i][j].pairs); @@ -305,7 +302,7 @@ index 3a8a6efa83..5e94793b79 100644 } } } -@@ -927,7 +924,7 @@ void TPPLPartition::TypeA(long i, long j, long k, PartitionVertex *vertices, DPS +@@ -927,7 +926,7 @@ void TPPLPartition::TypeA(long i, long j, long k, PartitionVertex *vertices, DPS void TPPLPartition::TypeB(long i, long j, long k, PartitionVertex *vertices, DPState2 **dpstates) { DiagonalList *pairs = NULL; @@ -314,7 +311,7 @@ index 3a8a6efa83..5e94793b79 100644 long top; long w; -@@ -946,21 +943,21 @@ void TPPLPartition::TypeB(long i, long j, long k, PartitionVertex *vertices, DPS +@@ -946,21 +945,21 @@ void TPPLPartition::TypeB(long i, long j, long k, PartitionVertex *vertices, DPS if (k - j > 1) { pairs = &(dpstates[j][k].pairs); @@ -343,7 +340,7 @@ index 3a8a6efa83..5e94793b79 100644 } } else { w++; -@@ -981,11 +978,11 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -981,11 +980,11 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { DiagonalList diagonals, diagonals2; Diagonal diagonal, newdiagonal; DiagonalList *pairs = NULL, *pairs2 = NULL; @@ -358,7 +355,7 @@ index 3a8a6efa83..5e94793b79 100644 bool ijreal, jkreal; n = poly->GetNumPoints(); -@@ -1110,35 +1107,35 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -1110,35 +1109,35 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { newdiagonal.index1 = 0; newdiagonal.index2 = n - 1; diagonals.push_front(newdiagonal); @@ -403,7 +400,7 @@ index 3a8a6efa83..5e94793b79 100644 pairs2->pop_back(); } else { break; -@@ -1153,21 +1150,21 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -1153,21 +1152,21 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { diagonals.push_front(newdiagonal); } } else { @@ -431,7 +428,7 @@ index 3a8a6efa83..5e94793b79 100644 pairs2->pop_front(); } else { break; -@@ -1197,8 +1194,8 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -1197,8 +1196,8 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { newdiagonal.index1 = 0; newdiagonal.index2 = n - 1; diagonals.push_front(newdiagonal); @@ -442,7 +439,7 @@ index 3a8a6efa83..5e94793b79 100644 diagonals.pop_front(); if ((diagonal.index2 - diagonal.index1) <= 1) { continue; -@@ -1210,8 +1207,8 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -1210,8 +1209,8 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { indices.push_back(diagonal.index2); diagonals2.push_front(diagonal); @@ -453,7 +450,7 @@ index 3a8a6efa83..5e94793b79 100644 diagonals2.pop_front(); if ((diagonal.index2 - diagonal.index1) <= 1) { continue; -@@ -1220,16 +1217,16 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -1220,16 +1219,16 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { jkreal = true; pairs = &(dpstates[diagonal.index1][diagonal.index2].pairs); if (!vertices[diagonal.index1].isConvex) { @@ -476,7 +473,7 @@ index 3a8a6efa83..5e94793b79 100644 jkreal = false; } } -@@ -1253,11 +1250,12 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -1253,11 +1252,12 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { indices.push_back(j); } @@ -492,7 +489,7 @@ index 3a8a6efa83..5e94793b79 100644 k++; } parts->push_back(newpoly); -@@ -1281,7 +1279,7 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { +@@ -1281,7 +1281,7 @@ int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) { // "Computational Geometry: Algorithms and Applications" // by Mark de Berg, Otfried Cheong, Marc van Kreveld, and Mark Overmars. int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monotonePolys) { @@ -501,7 +498,7 @@ index 3a8a6efa83..5e94793b79 100644 MonotoneVertex *vertices = NULL; long i, numvertices, vindex, vindex2, newnumvertices, maxnumvertices; long polystartindex, polyendindex; -@@ -1291,11 +1289,8 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto +@@ -1291,11 +1291,8 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto bool error = false; numvertices = 0; @@ -515,7 +512,7 @@ index 3a8a6efa83..5e94793b79 100644 } maxnumvertices = numvertices * 3; -@@ -1303,8 +1298,8 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto +@@ -1303,8 +1300,8 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto newnumvertices = numvertices; polystartindex = 0; @@ -526,7 +523,7 @@ index 3a8a6efa83..5e94793b79 100644 polyendindex = polystartindex + poly->GetNumPoints() - 1; for (i = 0; i < poly->GetNumPoints(); i++) { vertices[i + polystartindex].p = poly->GetPoint(i); -@@ -1360,14 +1355,14 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto +@@ -1360,14 +1357,14 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto // Note that while set doesn't actually have to be implemented as // a tree, complexity requirements for operations are the same as // for the balanced binary search tree. @@ -546,7 +543,7 @@ index 3a8a6efa83..5e94793b79 100644 } // For each vertex. -@@ -1387,13 +1382,14 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto +@@ -1387,13 +1384,14 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto newedge.p1 = v->p; newedge.p2 = vertices[v->next].p; newedge.index = vindex; @@ -564,7 +561,7 @@ index 3a8a6efa83..5e94793b79 100644 error = true; break; } -@@ -1412,29 +1408,30 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto +@@ -1412,29 +1410,30 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto newedge.p1 = v->p; newedge.p2 = v->p; edgeIter = edgeTree.lower_bound(newedge); @@ -601,7 +598,7 @@ index 3a8a6efa83..5e94793b79 100644 error = true; break; } -@@ -1452,25 +1449,25 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto +@@ -1452,25 +1451,25 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto newedge.p1 = v->p; newedge.p2 = v->p; edgeIter = edgeTree.lower_bound(newedge); @@ -632,7 +629,7 @@ index 3a8a6efa83..5e94793b79 100644 error = true; break; } -@@ -1488,27 +1485,28 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto +@@ -1488,27 +1487,28 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto newedge.p1 = v2->p; newedge.p2 = vertices[v2->next].p; newedge.index = vindex2; @@ -668,7 +665,7 @@ index 3a8a6efa83..5e94793b79 100644 } break; } -@@ -1569,8 +1567,8 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto +@@ -1569,8 +1569,8 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto // Adds a diagonal to the doubly-connected list of vertices. void TPPLPartition::AddDiagonal(MonotoneVertex *vertices, long *numvertices, long index1, long index2, @@ -679,7 +676,7 @@ index 3a8a6efa83..5e94793b79 100644 long newindex1, newindex2; newindex1 = *numvertices; -@@ -1597,14 +1595,14 @@ void TPPLPartition::AddDiagonal(MonotoneVertex *vertices, long *numvertices, lon +@@ -1597,14 +1597,14 @@ void TPPLPartition::AddDiagonal(MonotoneVertex *vertices, long *numvertices, lon vertextypes[newindex1] = vertextypes[index1]; edgeTreeIterators[newindex1] = edgeTreeIterators[index1]; helpers[newindex1] = helpers[index1]; @@ -698,7 +695,7 @@ index 3a8a6efa83..5e94793b79 100644 } } -@@ -1830,13 +1828,13 @@ int TPPLPartition::TriangulateMonotone(TPPLPoly *inPoly, TPPLPolyList *triangles +@@ -1830,13 +1830,13 @@ int TPPLPartition::TriangulateMonotone(TPPLPoly *inPoly, TPPLPolyList *triangles int TPPLPartition::Triangulate_MONO(TPPLPolyList *inpolys, TPPLPolyList *triangles) { TPPLPolyList monotone; diff --git a/thirdparty/misc/polypartition.cpp b/thirdparty/misc/polypartition.cpp index 5e94793b79..a725125ed0 100644 --- a/thirdparty/misc/polypartition.cpp +++ b/thirdparty/misc/polypartition.cpp @@ -23,6 +23,8 @@ #include "polypartition.h" +#include <math.h> +#include <string.h> #include <algorithm> TPPLPoly::TPPLPoly() { @@ -260,7 +262,7 @@ int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) { } } pointvisible = true; - for (iter2 = polys.front(); iter2; iter2->next()) { + for (iter2 = polys.front(); iter2; iter2 = iter2->next()) { if (iter2->get().IsHole()) { continue; } @@ -1355,12 +1357,12 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto // Note that while set doesn't actually have to be implemented as // a tree, complexity requirements for operations are the same as // for the balanced binary search tree. - Set<ScanLineEdge> edgeTree; + RBSet<ScanLineEdge> edgeTree; // Store iterators to the edge tree elements. // This makes deleting existing edges much faster. - Set<ScanLineEdge>::Element **edgeTreeIterators, *edgeIter; - edgeTreeIterators = new Set<ScanLineEdge>::Element *[maxnumvertices]; - //Pair<Set<ScanLineEdge>::iterator, bool> edgeTreeRet; + RBSet<ScanLineEdge>::Element **edgeTreeIterators, *edgeIter; + edgeTreeIterators = new RBSet<ScanLineEdge>::Element *[maxnumvertices]; + //Pair<RBSet<ScanLineEdge>::iterator, bool> edgeTreeRet; for (i = 0; i < numvertices; i++) { edgeTreeIterators[i] = nullptr; } @@ -1567,8 +1569,8 @@ int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monoto // Adds a diagonal to the doubly-connected list of vertices. void TPPLPartition::AddDiagonal(MonotoneVertex *vertices, long *numvertices, long index1, long index2, - TPPLVertexType *vertextypes, Set<ScanLineEdge>::Element **edgeTreeIterators, - Set<ScanLineEdge> *edgeTree, long *helpers) { + TPPLVertexType *vertextypes, RBSet<ScanLineEdge>::Element **edgeTreeIterators, + RBSet<ScanLineEdge> *edgeTree, long *helpers) { long newindex1, newindex2; newindex1 = *numvertices; diff --git a/thirdparty/misc/polypartition.h b/thirdparty/misc/polypartition.h index b2d905a3ef..fae7909079 100644 --- a/thirdparty/misc/polypartition.h +++ b/thirdparty/misc/polypartition.h @@ -26,7 +26,7 @@ #include "core/math/vector2.h" #include "core/templates/list.h" -#include "core/templates/set.h" +#include "core/templates/rb_set.h" typedef double tppl_float; @@ -224,8 +224,8 @@ public: // Helper functions for MonotonePartition. bool Below(TPPLPoint &p1, TPPLPoint &p2); void AddDiagonal(MonotoneVertex *vertices, long *numvertices, long index1, long index2, - TPPLVertexType *vertextypes, Set<ScanLineEdge>::Element **edgeTreeIterators, - Set<ScanLineEdge> *edgeTree, long *helpers); + TPPLVertexType *vertextypes, RBSet<ScanLineEdge>::Element **edgeTreeIterators, + RBSet<ScanLineEdge> *edgeTree, long *helpers); // Triangulates a monotone polygon, used in Triangulate_MONO. int TriangulateMonotone(TPPLPoly *inPoly, TPPLPolyList *triangles); diff --git a/thirdparty/misc/stb_rect_pack.h b/thirdparty/misc/stb_rect_pack.h index 5c848de0e7..6a633ce666 100644 --- a/thirdparty/misc/stb_rect_pack.h +++ b/thirdparty/misc/stb_rect_pack.h @@ -1,9 +1,15 @@ -// stb_rect_pack.h - v1.00 - public domain - rectangle packing +// stb_rect_pack.h - v1.01 - public domain - rectangle packing // Sean Barrett 2014 // // Useful for e.g. packing rectangular textures into an atlas. // Does not do rotation. // +// Before #including, +// +// #define STB_RECT_PACK_IMPLEMENTATION +// +// in the file that you want to have the implementation. +// // Not necessarily the awesomest packing method, but better than // the totally naive one in stb_truetype (which is primarily what // this is meant to replace). @@ -35,6 +41,7 @@ // // Version history: // +// 1.01 (2021-07-11) always use large rect mode, expose STBRP__MAXVAL in public section // 1.00 (2019-02-25) avoid small space waste; gracefully fail too-wide rectangles // 0.99 (2019-02-07) warning fixes // 0.11 (2017-03-03) return packing success/fail result @@ -75,11 +82,10 @@ typedef struct stbrp_context stbrp_context; typedef struct stbrp_node stbrp_node; typedef struct stbrp_rect stbrp_rect; -#ifdef STBRP_LARGE_RECTS typedef int stbrp_coord; -#else -typedef unsigned short stbrp_coord; -#endif + +#define STBRP__MAXVAL 0x7fffffff +// Mostly for internal use, but this is the maximum supported coordinate value. STBRP_DEF int stbrp_pack_rects (stbrp_context *context, stbrp_rect *rects, int num_rects); // Assign packed locations to rectangles. The rectangles are of type @@ -209,8 +215,10 @@ struct stbrp_context #ifdef _MSC_VER #define STBRP__NOTUSED(v) (void)(v) +#define STBRP__CDECL __cdecl #else #define STBRP__NOTUSED(v) (void)sizeof(v) +#define STBRP__CDECL #endif enum @@ -253,9 +261,6 @@ STBRP_DEF void stbrp_setup_allow_out_of_mem(stbrp_context *context, int allow_ou STBRP_DEF void stbrp_init_target(stbrp_context *context, int width, int height, stbrp_node *nodes, int num_nodes) { int i; -#ifndef STBRP_LARGE_RECTS - STBRP_ASSERT(width <= 0xffff && height <= 0xffff); -#endif for (i=0; i < num_nodes-1; ++i) nodes[i].next = &nodes[i+1]; @@ -274,11 +279,7 @@ STBRP_DEF void stbrp_init_target(stbrp_context *context, int width, int height, context->extra[0].y = 0; context->extra[0].next = &context->extra[1]; context->extra[1].x = (stbrp_coord) width; -#ifdef STBRP_LARGE_RECTS context->extra[1].y = (1<<30); -#else - context->extra[1].y = 65535; -#endif context->extra[1].next = NULL; } @@ -520,7 +521,7 @@ static stbrp__findresult stbrp__skyline_pack_rectangle(stbrp_context *context, i return res; } -static int rect_height_compare(const void *a, const void *b) +static int STBRP__CDECL rect_height_compare(const void *a, const void *b) { const stbrp_rect *p = (const stbrp_rect *) a; const stbrp_rect *q = (const stbrp_rect *) b; @@ -531,19 +532,13 @@ static int rect_height_compare(const void *a, const void *b) return (p->w > q->w) ? -1 : (p->w < q->w); } -static int rect_original_order(const void *a, const void *b) +static int STBRP__CDECL rect_original_order(const void *a, const void *b) { const stbrp_rect *p = (const stbrp_rect *) a; const stbrp_rect *q = (const stbrp_rect *) b; return (p->was_packed < q->was_packed) ? -1 : (p->was_packed > q->was_packed); } -#ifdef STBRP_LARGE_RECTS -#define STBRP__MAXVAL 0xffffffff -#else -#define STBRP__MAXVAL 0xffff -#endif - STBRP_DEF int stbrp_pack_rects(stbrp_context *context, stbrp_rect *rects, int num_rects) { int i, all_rects_packed = 1; |