tei/MapLibreParticleLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

core: Update Hillshade Algorithms and add Color-Relief Layer support (#3965)

Browse Source
This commit is contained in:
Andrew Calcutt
2026-02-25 20:50:06 -05:00
committed by GitHub
parent 9682e09ae4
commit 369af13ead
187 changed files with 8804 additions and 2352 deletions
Download Patch File Download Diff File Expand all files Collapse all files

77
shaders/color_relief.fragment.glsl Normal file
View File

@@ -0,0 +1,77 @@
layout(std140) uniform ColorReliefDrawableUBO {
highp mat4 u_matrix;
};
layout(std140) uniform ColorReliefTilePropsUBO {
highp vec4 u_unpack;
highp vec2 u_dimension;
highp int u_color_ramp_size;
highp float pad_tile0;
};
layout(std140) uniform ColorReliefEvaluatedPropsUBO {
highp float u_opacity;
highp float pad_eval0;
highp float pad_eval1;
highp float pad_eval2;
};
uniform sampler2D u_image;
uniform sampler2D u_elevation_stops;
uniform sampler2D u_color_stops;
in vec2 v_pos;
float getElevation(vec2 coord) {
// Convert encoded elevation value to meters
vec4 data = texture(u_image, coord) * 255.0;
data.a = -1.0;
return dot(data, u_unpack);
}
float getElevationStop(int stop) {
// Elevation stops are plain float values, not terrain-RGB encoded
float x = (float(stop) + 0.5) / float(u_color_ramp_size);
return texture(u_elevation_stops, vec2(x, 0.5)).r;
}
vec4 getColorStop(int stop) {
float x = (float(stop) + 0.5) / float(u_color_ramp_size);
return texture(u_color_stops, vec2(x, 0.5));
}
void main() {
float el = getElevation(v_pos);
// Binary search for color stops
int r = (u_color_ramp_size - 1);
int l = 0;
while (r - l > 1) {
int m = (r + l) / 2;
float el_m = getElevationStop(m);
if (el < el_m) {
r = m;
} else {
l = m;
}
}
// Get elevation values for interpolation
float el_l = getElevationStop(l);
float el_r = getElevationStop(r);
// Get colors for interpolation
vec4 color_l = getColorStop(l);
vec4 color_r = getColorStop(r);
// Interpolate between the two colors
// Guard against division by zero when el_r == el_l
float denom = el_r - el_l;
float t = (abs(denom) < 0.0001) ? 0.0 : clamp((el - el_l) / denom, 0.0, 1.0);
fragColor = u_opacity * mix(color_l, color_r, t);
#ifdef OVERDRAW_INSPECTOR
fragColor = vec4(1.0);
#endif
}

26
shaders/color_relief.vertex.glsl Normal file
View File

@@ -0,0 +1,26 @@
layout(std140) uniform ColorReliefDrawableUBO {
highp mat4 u_matrix;
};
layout(std140) uniform ColorReliefTilePropsUBO {
highp vec4 u_unpack;
highp vec2 u_dimension;
highp int u_color_ramp_size;
highp float pad_tile0;
};
layout (location = 0) in vec2 a_pos;
layout (location = 1) in vec2 a_texture_pos;
out vec2 v_pos;
void main() {
gl_Position = u_matrix * vec4(a_pos, 0, 1);
highp vec2 epsilon = 1.0 / u_dimension;
float scale = (u_dimension.x - 2.0) / u_dimension.x;
v_pos = (a_texture_pos / 8192.0) * scale + epsilon;
// Handle poles
if (a_pos.y < -32767.5) v_pos.y = 0.0;
if (a_pos.y > 32766.5) v_pos.y = 1.0;
}

2
shaders/heatmap.vertex.glsl
View File

@@ -46,7 +46,7 @@ void main(void) {
// weight * u_intensity * GAUSS_COEF * exp(-0.5 * 3.0^2 * S^2) == ZERO
// Which solves to:
// S = sqrt(-2.0 * log(ZERO / (weight * u_intensity * GAUSS_COEF))) / 3.0
float S = sqrt(-2.0 * log(ZERO / weight / u_intensity / GAUSS_COEF)) / 3.0;
float S = sqrt(-2.0 * log(ZERO / (max(weight, ZERO) * max(u_intensity, ZERO) * GAUSS_COEF))) / 3.0;
// Pass the varying in units of radius
v_extrude = S * unscaled_extrude;

202
shaders/hillshade.fragment.glsl
View File

@@ -1,55 +1,185 @@
in vec2 v_pos;
uniform sampler2D u_image;
layout (std140) uniform HillshadeTilePropsUBO {
layout(std140) uniform HillshadeTilePropsUBO {
highp vec2 u_latrange;
highp vec2 u_light;
highp float u_exaggeration;
highp int u_method; // Hillshade method (0: STANDARD, 1: COMBINED, 2: IGOR, 3: MULTIDIRECTIONAL, 4: BASIC)
highp int u_num_lights; // Number of light sources (1-4)
highp float u_pad0;
highp float u_pad1;
highp float u_pad2;
};
layout (std140) uniform HillshadeEvaluatedPropsUBO {
highp vec4 u_highlight;
highp vec4 u_shadow;
layout(std140) uniform HillshadeEvaluatedPropsUBO {
highp vec4 u_accent;
highp vec4 u_altitudes; // Up to 4 altitude values (in radians)
highp vec4 u_azimuths; // Up to 4 azimuth values (in radians)
highp vec4 u_shadows[4]; // Shadow colors (up to 4 lights)
highp vec4 u_highlights[4]; // Highlight colors (up to 4 lights)
};
#define PI 3.141592653589793
#define STANDARD 0
#define COMBINED 1
#define IGOR 2
#define MULTIDIRECTIONAL 3
#define BASIC 4
float get_aspect(vec2 deriv) {
return deriv.x != 0.0 ?
atan(deriv.y, -deriv.x) : PI / 2.0 * (deriv.y > 0.0 ? 1.0 : -1.0);
}
// MapLibre's legacy hillshade algorithm (Method 0)
void standard_hillshade(vec2 deriv) {
float azimuth = u_azimuths.x + PI;
float slope = atan(0.625 * length(deriv));
float aspect = get_aspect(deriv);
// Note: This implementation uses u_exaggeration as intensity, though it typically should be derived from u_altitudes.x.
float intensity = u_exaggeration;
// Scale the slope exponentially based on intensity
float base = 1.875 - intensity * 1.75;
float maxValue = 0.5 * PI;
float scaledSlope = abs(intensity - 0.5) > 1e-6 ? ((pow(base, slope) - 1.0) / (pow(base, maxValue) - 1.0)) * maxValue : slope;
float accent = cos(scaledSlope);
vec4 accent_color = (1.0 - accent) * u_accent * clamp(intensity * 2.0, 0.0, 1.0);
// Calculate shadow/highlight based on aspect
float shade = abs(mod((aspect + azimuth) / PI + 0.5, 2.0) - 1.0);
vec4 shade_color = mix(u_shadows[0], u_highlights[0], shade) * sin(scaledSlope) * clamp(intensity * 2.0, 0.0, 1.0);
fragColor = accent_color * (1.0 - shade_color.a) + shade_color;
}
// Basic directional hillshade (Method 4)
void basic_hillshade(vec2 deriv) {
deriv = deriv * u_exaggeration * 2.0; // Exaggerate the slope derivative
float azimuth = u_azimuths.x + PI;
float cos_az = cos(azimuth);
float sin_az = sin(azimuth);
float cos_alt = cos(u_altitudes.x);
float sin_alt = sin(u_altitudes.x);
// Calculate the cosine of the angle between the light vector and the surface normal
float cang = (sin_alt - (deriv.y * cos_az * cos_alt - deriv.x * sin_az * cos_alt)) / sqrt(1.0 + dot(deriv, deriv));
float shade = clamp(cang, 0.0, 1.0); // cang is the hillshade intensity [0, 1]
// Blend shadow and highlight based on intensity
if (shade > 0.5) {
fragColor = u_highlights[0] * (2.0 * shade - 1.0); // Highlight strength [0, 1]
} else {
fragColor = u_shadows[0] * (1.0 - 2.0 * shade); // Shadow strength [0, 1]
}
}
// Multidirectional hillshade (Method 3)
void multidirectional_hillshade(vec2 deriv) {
deriv = deriv * u_exaggeration * 2.0; // Exaggerate the slope derivative
fragColor = vec4(0, 0, 0, 0);
// Iterate through all light sources (up to u_num_lights, max 4)
for (int i = 0; i < 4; i++) {
if (i >= u_num_lights) break;
// Access altitude and azimuth from vec4 UBOs
float altitude = (i == 0) ? u_altitudes.x : (i == 1) ? u_altitudes.y : (i == 2) ? u_altitudes.z : u_altitudes.w;
float azimuth = (i == 0) ? u_azimuths.x : (i == 1) ? u_azimuths.y : (i == 2) ? u_azimuths.z : u_azimuths.w;
float cos_alt = cos(altitude);
float sin_alt = sin(altitude);
// Negate cos/sin azimuth for correct light direction in the normal calculation
float cos_az = -cos(azimuth);
float sin_az = -sin(azimuth);
// Calculate the cosine of the angle between the light vector and the surface normal
float cang = (sin_alt - (deriv.y * cos_az * cos_alt - deriv.x * sin_az * cos_alt)) /
sqrt(1.0 + dot(deriv, deriv));
float shade = clamp(cang, 0.0, 1.0); // cang is the hillshade intensity [0, 1]
// Accumulate shadow/highlight contribution from each light
if (shade > 0.5) {
fragColor += u_highlights[i] * (2.0 * shade - 1.0) / float(u_num_lights);
} else {
fragColor += u_shadows[i] * (1.0 - 2.0 * shade) / float(u_num_lights);
}
}
}
// Combined shadow and highlight method (Method 1)
void combined_hillshade(vec2 deriv) {
// Only supports one light source (index 0)
deriv = deriv * u_exaggeration * 2.0;
float azimuth = u_azimuths.x + PI;
float cos_az = cos(azimuth);
float sin_az = sin(azimuth);
float cos_alt = cos(u_altitudes.x);
float sin_alt = sin(u_altitudes.x);
// Calculate the angle between the light vector and the surface normal
float cang = acos((sin_alt - (deriv.y * cos_az * cos_alt - deriv.x * sin_az * cos_alt)) /
sqrt(1.0 + dot(deriv, deriv)));
cang = clamp(cang, 0.0, PI / 2.0); // Clamp angle to 90 degrees (half hemisphere)
// The "shade" and "highlight" components are calculated from cang (angle) and the magnitude of the slope
float shade = cang * atan(length(deriv)) * 4.0 / PI / PI;
float highlight = (PI / 2.0 - cang) * atan(length(deriv)) * 4.0 / PI / PI;
fragColor = u_shadows[0] * shade + u_highlights[0] * highlight;
}
// Igor's shadow/highlight method (Method 2)
void igor_hillshade(vec2 deriv) {
// Only supports one light source (index 0)
deriv = deriv * u_exaggeration * 2.0;
float aspect = get_aspect(deriv);
float azimuth = u_azimuths.x + PI;
// Slope strength is magnitude of slope vector, normalized to [0, 1]
float slope_strength = atan(length(deriv)) * 2.0 / PI;
// Aspect strength is difference between aspect and light azimuth, normalized to [0, 1]
float aspect_strength = 1.0 - abs(mod((aspect + azimuth) / PI + 0.5, 2.0) - 1.0);
float shadow_strength = slope_strength * aspect_strength;
float highlight_strength = slope_strength * (1.0 - aspect_strength);
fragColor = u_shadows[0] * shadow_strength + u_highlights[0] * highlight_strength;
}
void main() {
vec4 pixel = texture(u_image, v_pos);
vec2 deriv = ((pixel.rg * 2.0) - 1.0);
// We divide the slope by a scale factor based on the cosin of the pixel's approximate latitude
// to account for mercator projection distortion. see #4807 for details
// Scale the derivative based on the mercator distortion at this latitude (v_pos.y)
float scaleFactor = cos(radians((u_latrange[0] - u_latrange[1]) * (1.0 - v_pos.y) + u_latrange[1]));
// We also multiply the slope by an arbitrary z-factor of 1.25
float slope = atan(1.25 * length(deriv) / scaleFactor);
float aspect = deriv.x != 0.0 ? atan(deriv.y, -deriv.x) : PI / 2.0 * (deriv.y > 0.0 ? 1.0 : -1.0);
float intensity = u_light.x;
// We add PI to make this property match the global light object, which adds PI/2 to the light's azimuthal
// position property to account for 0deg corresponding to north/the top of the viewport in the style spec
// and the original shader was written to accept (-illuminationDirection - 90) as the azimuthal.
float azimuth = u_light.y + PI;
// The derivative is scaled back from [0, 1] texture range to world-space slope
// Texture range [0, 1] corresponds to slope range [-4, 4] (8.0 * 0.5 * deriv)
vec2 deriv = ((pixel.rg * 8.0) - 4.0) / scaleFactor;
// We scale the slope exponentially based on intensity, using a calculation similar to
// the exponential interpolation function in the style spec:
// https://github.com/mapbox/mapbox-gl-js/blob/master/src/style-spec/expression/definitions/interpolate.js#L217-L228
// so that higher intensity values create more opaque hillshading.
float base = 1.875 - intensity * 1.75;
float maxValue = 0.5 * PI;
float scaledSlope = intensity != 0.5 ? ((pow(base, slope) - 1.0) / (pow(base, maxValue) - 1.0)) * maxValue : slope;
// The accent color is calculated with the cosine of the slope while the shade color is calculated with the sine
// so that the accent color's rate of change eases in while the shade color's eases out.
float accent = cos(scaledSlope);
// We multiply both the accent and shade color by a clamped intensity value
// so that intensities >= 0.5 do not additionally affect the color values
// while intensity values < 0.5 make the overall color more transparent.
vec4 accent_color = (1.0 - accent) * u_accent * clamp(intensity * 2.0, 0.0, 1.0);
float shade = abs(mod((aspect + azimuth) / PI + 0.5, 2.0) - 1.0);
vec4 shade_color = mix(u_shadow, u_highlight, shade) * sin(scaledSlope) * clamp(intensity * 2.0, 0.0, 1.0);
fragColor = accent_color * (1.0 - shade_color.a) + shade_color;
// Dispatch to the selected hillshade method
switch (u_method) {
case BASIC:
basic_hillshade(deriv);
break;
case COMBINED:
combined_hillshade(deriv);
break;
case IGOR:
igor_hillshade(deriv);
break;
case MULTIDIRECTIONAL:
multidirectional_hillshade(deriv);
break;
case STANDARD:
default:
standard_hillshade(deriv);
break;
}
#ifdef OVERDRAW_INSPECTOR
fragColor = vec4(1.0);

34
shaders/hillshade_prepare.fragment.glsl
View File

@@ -4,7 +4,6 @@ precision highp float;
in vec2 v_pos;
uniform sampler2D u_image;
layout (std140) uniform HillshadePrepareTilePropsUBO {
highp vec4 u_unpack;
highp vec2 u_dimension;
@@ -16,13 +15,14 @@ float getElevation(vec2 coord, float bias) {
// Convert encoded elevation value to meters
vec4 data = texture(u_image, coord) * 255.0;
data.a = -1.0;
return dot(data, u_unpack) / 4.0;
return dot(data, u_unpack);
}
void main() {
vec2 epsilon = 1.0 / u_dimension;
float tileSize = u_dimension.x - 2.0;
// queried pixels:
// queried pixels (using Sobel operator kernel):
// +-----------+
// | | | |
// | a | b | c |
@@ -47,30 +47,28 @@ void main() {
float h = getElevation(v_pos + vec2(0, epsilon.y), 0.0);
float i = getElevation(v_pos + vec2(epsilon.x, epsilon.y), 0.0);
// here we divide the x and y slopes by 8 * pixel size
// where pixel size (aka meters/pixel) is:
// circumference of the world / (pixels per tile * number of tiles)
// which is equivalent to: 8 * 40075016.6855785 / (512 * pow(2, u_zoom))
// which can be reduced to: pow(2, 19.25619978527 - u_zoom)
// we want to vertically exaggerate the hillshading though, because otherwise
// it is barely noticeable at low zooms. to do this, we multiply this by some
// scale factor pow(2, (u_zoom - u_maxzoom) * a) where a is an arbitrary value
// Here we use a=0.3 which works out to the expression below. see
// nickidlugash's awesome breakdown for more info
// Convert the raw pixel-space derivative (slope) into world-space slope.
// The conversion factor is: tileSize / (8 * meters_per_pixel).
// meters_per_pixel is calculated as pow(2.0, 28.2562 - u_zoom).
// The exaggeration factor is applied to scale the effect at lower zooms.
// See nickidlugash's awesome breakdown for more info
// https://github.com/mapbox/mapbox-gl-js/pull/5286#discussion_r148419556
float exaggeration = u_zoom < 2.0 ? 0.4 : u_zoom < 4.5 ? 0.35 : 0.3;
float exaggerationFactor = u_zoom < 2.0 ? 0.4 : u_zoom < 4.5 ? 0.35 : 0.3;
float exaggeration = u_zoom < 15.0 ? (u_zoom - 15.0) * exaggerationFactor : 0.0;
vec2 deriv = vec2(
(c + f + f + i) - (a + d + d + g),
(g + h + h + i) - (a + b + b + c)
) / pow(2.0, (u_zoom - u_maxzoom) * exaggeration + 19.2562 - u_zoom);
) * tileSize / pow(2.0, exaggeration + (28.2562 - u_zoom));
// Encode the derivative into the color channels (r and g)
// The derivative is scaled from world-space slope to the range [0, 1] for texture storage.
// The maximum possible world-space derivative is assumed to be 4 (hence division by 8.0).
fragColor = clamp(vec4(
deriv.x / 2.0 + 0.5,
deriv.y / 2.0 + 0.5,
deriv.x / 8.0 + 0.5,
deriv.y / 8.0 + 0.5,
1.0,
1.0), 0.0, 1.0);
#ifdef OVERDRAW_INSPECTOR
fragColor = vec4(1.0);
#endif

7
shaders/manifest.json
View File

@@ -137,6 +137,13 @@
"glsl_frag": "hillshade_prepare.fragment.glsl",
"uses_ubos": true
},
{
"name": "ColorReliefShader",
"header": "color_relief",
"glsl_vert": "color_relief.vertex.glsl",
"glsl_frag": "color_relief.fragment.glsl",
"uses_ubos": true
},
{
"name": "HillshadeShader",
"header": "hillshade",