/** * DOC_TBA * * @name czm_infinity * @glslConstant */ const float czm_infinity = 5906376272000.0; // Distance from the Sun to Pluto in meters. TODO: What is best given lowp, mediump, and highp? /** * DOC_TBA * * @name czm_epsilon1 * @glslConstant */ const float czm_epsilon1 = 0.1; /** * DOC_TBA * * @name czm_epsilon2 * @glslConstant */ const float czm_epsilon2 = 0.01; /** * DOC_TBA * * @name czm_epsilon3 * @glslConstant */ const float czm_epsilon3 = 0.001; /** * DOC_TBA * * @name czm_epsilon4 * @glslConstant */ const float czm_epsilon4 = 0.0001; /** * DOC_TBA * * @name czm_epsilon5 * @glslConstant */ const float czm_epsilon5 = 0.00001; /** * DOC_TBA * * @name czm_epsilon6 * @glslConstant */ const float czm_epsilon6 = 0.000001; /** * DOC_TBA * * @name czm_epsilon7 * @glslConstant */ const float czm_epsilon7 = 0.0000001; /** * Compares <code>left</code> and <code>right</code> componentwise. Returns <code>true</code> * if they are within <code>epsilon</code> and <code>false</code> otherwise. The inputs * <code>left</code> and <code>right</code> can be <code>float</code>s, <code>vec2</code>s, * <code>vec3</code>s, or <code>vec4</code>s. * * @name czm_equalsEpsilon * @glslFunction * * @param {} left The first vector. * @param {} right The second vector. * @param {float} epsilon The epsilon to use for equality testing. * @returns {bool} <code>true</code> if the components are within <code>epsilon</code> and <code>false</code> otherwise. * * @example * // GLSL declarations * bool czm_equalsEpsilon(float left, float right, float epsilon); * bool czm_equalsEpsilon(vec2 left, vec2 right, float epsilon); * bool czm_equalsEpsilon(vec3 left, vec3 right, float epsilon); * bool czm_equalsEpsilon(vec4 left, vec4 right, float epsilon); */ bool czm_equalsEpsilon(vec4 left, vec4 right, float epsilon) { return all(lessThanEqual(abs(left - right), vec4(epsilon))); } bool czm_equalsEpsilon(vec3 left, vec3 right, float epsilon) { return all(lessThanEqual(abs(left - right), vec3(epsilon))); } bool czm_equalsEpsilon(vec2 left, vec2 right, float epsilon) { return all(lessThanEqual(abs(left - right), vec2(epsilon))); } bool czm_equalsEpsilon(float left, float right, float epsilon) { return (abs(left - right) <= epsilon); } bool czm_equalsEpsilon(float left, float right) { // Workaround bug in Opera Next 12. Do not delegate to the other czm_equalsEpsilon. return (abs(left - right) <= czm_epsilon7); } /////////////////////////////////////////////////////////////////////////////// /** * Returns the transpose of the matrix. The input <code>matrix</code> can be * a <code>mat2</code>, <code>mat3</code>, or <code>mat4</code>. * * @name czm_transpose * @glslFunction * * @param {} matrix The matrix to transpose. * * @returns {} The transposed matrix. * * @example * // GLSL declarations * mat2 czm_transpose(mat2 matrix); * mat3 czm_transpose(mat3 matrix); * mat4 czm_transpose(mat4 matrix); * * // Tranpose a 3x3 rotation matrix to find its inverse. * mat3 eastNorthUpToEye = czm_eastNorthUpToEyeCoordinates( * positionMC, normalEC); * mat3 eyeToEastNorthUp = czm_transpose(eastNorthUpToEye); */ mat2 czm_transpose(mat2 matrix) { return mat2( matrix[0][0], matrix[1][0], matrix[0][1], matrix[1][1]); } mat3 czm_transpose(mat3 matrix) { return mat3( matrix[0][0], matrix[1][0], matrix[2][0], matrix[0][1], matrix[1][1], matrix[2][1], matrix[0][2], matrix[1][2], matrix[2][2]); } mat4 czm_transpose(mat4 matrix) { return mat4( matrix[0][0], matrix[1][0], matrix[2][0], matrix[3][0], matrix[0][1], matrix[1][1], matrix[2][1], matrix[3][1], matrix[0][2], matrix[1][2], matrix[2][2], matrix[3][2], matrix[0][3], matrix[1][3], matrix[2][3], matrix[3][3]); } /////////////////////////////////////////////////////////////////////////////// /** * Transforms a position from model to window coordinates. The transformation * from model to clip coordinates is done using {@link czm_modelViewProjection}. * The transform from normalized device coordinates to window coordinates is * done using {@link czm_viewportTransformation}, which assumes a depth range * of <code>near = 0</code> and <code>far = 1</code>. * <br /><br /> * This transform is useful when there is a need to manipulate window coordinates * in a vertex shader as done by {@link BillboardCollection}. * <br /><br /> * This function should not be confused with {@link czm_viewportOrthographic}, * which is an orthographic projection matrix that transforms from window * coordinates to clip coordinates. * * @name czm_modelToWindowCoordinates * @glslFunction * * @param {vec4} position The position in model coordinates to transform. * * @returns {vec4} The transformed position in window coordinates. * * @see czm_eyeToWindowCoordinates * @see czm_modelViewProjection * @see czm_viewportTransformation * @see czm_viewportOrthographic * @see BillboardCollection * * @example * vec4 positionWC = czm_modelToWindowCoordinates(positionMC); */ vec4 czm_modelToWindowCoordinates(vec4 position) { vec4 q = czm_modelViewProjection * position; // clip coordinates q.xyz /= q.w; // normalized device coordinates q.xyz = (czm_viewportTransformation * vec4(q.xyz, 1.0)).xyz; // window coordinates return q; } /** * Transforms a position from eye to window coordinates. The transformation * from eye to clip coordinates is done using {@link czm_projection}. * The transform from normalized device coordinates to window coordinates is * done using {@link czm_viewportTransformation}, which assumes a depth range * of <code>near = 0</code> and <code>far = 1</code>. * <br /><br /> * This transform is useful when there is a need to manipulate window coordinates * in a vertex shader as done by {@link BillboardCollection}. * * @name czm_eyeToWindowCoordinates * @glslFunction * * @param {vec4} position The position in eye coordinates to transform. * * @returns {vec4} The transformed position in window coordinates. * * @see czm_modelToWindowCoordinates * @see czm_projection * @see czm_viewportTransformation * @see BillboardCollection * * @example * vec4 positionWC = czm_eyeToWindowCoordinates(positionEC); */ vec4 czm_eyeToWindowCoordinates(vec4 positionEC) { vec4 q = czm_projection * positionEC; // clip coordinates q.xyz /= q.w; // normalized device coordinates q.xyz = (czm_viewportTransformation * vec4(q.xyz, 1.0)).xyz; // window coordinates return q; } /** * Transforms a position from window to eye coordinates. * The transform from window to normalized device coordinates is done using components * of (@link czm_viewport} and {@link czm_viewportTransformation} instead of calculating * the inverse of <code>czm_viewportTransformation</code>. The transformation from * normalized device coordinates to clip coordinates is done using <code>positionWC.w</code>, * which is expected to be the scalar used in the perspective divide. The transformation * from clip to eye coordinates is done using {@link czm_inverseProjection}. * * @name czm_windowToEyeCoordinates * @glslFunction * * @param {vec4} fragmentCoordinate The position in window coordinates to transform. * * @returns {vec4} The transformed position in eye coordinates. * * @see czm_modelToWindowCoordinates * @see czm_eyeToWindowCoordinates * @see czm_inverseProjection * @see czm_viewport * @see czm_viewportTransformation * * @example * vec4 positionEC = czm_windowToEyeCoordinates(gl_FragCoord); */ vec4 czm_windowToEyeCoordinates(vec4 fragmentCoordinate) { float x = 2.0 * (fragmentCoordinate.x - czm_viewport.x) / czm_viewport.z - 1.0; float y = 2.0 * (fragmentCoordinate.y - czm_viewport.y) / czm_viewport.w - 1.0; float z = (fragmentCoordinate.z - czm_viewportTransformation[3][2]) / czm_viewportTransformation[2][2]; vec4 q = vec4(x, y, z, 1.0); q /= fragmentCoordinate.w; q = czm_inverseProjection * q; return q; } /////////////////////////////////////////////////////////////////////////////// /** * DOC_TBA * * @name czm_eyeOffset * @glslFunction * * @param {vec4} positionEC DOC_TBA. * @param {vec3} eyeOffset DOC_TBA. * * @returns {vec4} DOC_TBA. */ vec4 czm_eyeOffset(vec4 positionEC, vec3 eyeOffset) { // This equation is approximate in x and y. vec4 p = positionEC; vec4 zEyeOffset = normalize(p) * eyeOffset.z; p.xy += eyeOffset.xy + zEyeOffset.xy; p.z += zEyeOffset.z; return p; } /////////////////////////////////////////////////////////////////////////////// /** * DOC_TBA * * @name czm_geodeticSurfaceNormal * @glslFunction * * @param {vec3} positionOnEllipsoid DOC_TBA * @param {vec3} ellipsoidCenter DOC_TBA * @param {vec3} oneOverEllipsoidRadiiSquared DOC_TBA * * @returns {vec3} DOC_TBA. */ vec3 czm_geodeticSurfaceNormal(vec3 positionOnEllipsoid, vec3 ellipsoidCenter, vec3 oneOverEllipsoidRadiiSquared) { return normalize((positionOnEllipsoid - ellipsoidCenter) * oneOverEllipsoidRadiiSquared); } /** * DOC_TBA * * @name czm_ellipsoidWgs84TextureCoordinates * @glslFunction */ vec2 czm_ellipsoidWgs84TextureCoordinates(vec3 normal) { return vec2(atan(normal.y, normal.x) * czm_oneOverTwoPi + 0.5, asin(normal.z) * czm_oneOverPi + 0.5); } /** * Computes a 3x3 rotation matrix that transforms vectors from an ellipsoid's east-north-up coordinate system * to eye coordinates. In east-north-up coordinates, x points east, y points north, and z points along the * surface normal. East-north-up can be used as an ellipsoid's tangent space for operations such as bump mapping. * <br /><br /> * The ellipsoid is assumed to be centered at the model coordinate's origin. * * @name czm_eastNorthUpToEyeCoordinates * @glslFunction * * @param {vec3} positionMC The position on the ellipsoid in model coordinates. * @param {vec3} normalEC The normalized ellipsoid surface normal, at <code>positionMC</code>, in eye coordinates. * * @returns {mat3} A 3x3 rotation matrix that transforms vectors from the east-north-up coordinate system to eye coordinates. * * @example * // Transform a vector defined in the east-north-up coordinate * // system, (0, 0, 1) which is the surface normal, to eye * // coordinates. * mat3 m = czm_eastNorthUpToEyeCoordinates(positionMC, normalEC); * vec3 normalEC = m * vec3(0.0, 0.0, 1.0); */ mat3 czm_eastNorthUpToEyeCoordinates(vec3 positionMC, vec3 normalEC) { vec3 tangentMC = normalize(vec3(-positionMC.y, positionMC.x, 0.0)); // normalized surface tangent in model coordinates vec3 tangentEC = normalize(czm_normal3D * tangentMC); // normalized surface tangent in eye coordiantes vec3 bitangentEC = normalize(cross(normalEC, tangentEC)); // normalized surface bitangent in eye coordinates return mat3( tangentEC.x, tangentEC.y, tangentEC.z, bitangentEC.x, bitangentEC.y, bitangentEC.z, normalEC.x, normalEC.y, normalEC.z); } /////////////////////////////////////////////////////////////////////////////// /** * Used as input to every material's czm_getMaterial function. * * @name czm_materialInput * @glslStruct * * @property {float} s 1D texture coordinates. * @property {vec2} st 2D texture coordinates. * @property {vec3} str 3D texture coordinates. * @property {vec3} normalEC Unperturbed surface normal in eye coordinates. * @property {mat3} tangentToEyeMatrix Matrix for converting a tangent space normal to eye space. * @property {vec3} positionToEyeEC Vector from the fragment to the eye in eye coordinates. The magnitude is the distance in meters from the fragment to the eye. */ struct czm_materialInput { float s; vec2 st; vec3 str; vec3 normalEC; mat3 tangentToEyeMatrix; vec3 positionToEyeEC; }; /** * Holds material information that can be used for lighting. Returned by all czm_getMaterial functions. * * @name czm_material * @glslStruct * * @property {vec3} diffuse Incoming light that scatters evenly in all directions. * @property {float} specular Intensity of incoming light reflecting in a single direction. * @property {float} shininess The sharpness of the specular reflection. Higher values create a smaller, more focused specular highlight. * @property {vec3} normal Surface's normal in eye coordinates. It is used for effects such as normal mapping. The default is the surface's unmodified normal. * @property {vec3} emission Light emitted by the material equally in all directions. The default is vec3(0.0), which emits no light. * @property {float} alpha Opacity of this material. 0.0 is completely transparent; 1.0 is completely opaque. */ struct czm_material { vec3 diffuse; float specular; float shininess; vec3 normal; vec3 emission; float alpha; }; /** * An czm_material with default values. Every material's czm_getMaterial * should use this default material as a base for the material it returns. * The default normal value is given by materialInput.normalEC. * * @name czm_getDefaultMaterial * @glslFunction * * @param {czm_materialInput} input The input used to construct the default material. * * @returns {czm_material} The default material. * * @see czm_materialInput * @see czm_material * @see czm_getMaterial */ czm_material czm_getDefaultMaterial(czm_materialInput materialInput) { czm_material material; material.diffuse = vec3(0.0); material.specular = 0.0; material.shininess = 1.0; material.normal = materialInput.normalEC; material.emission = vec3(0.0); material.alpha = 1.0; return material; } float getLambertDiffuse(vec3 lightDirectionEC, vec3 normalEC) { return max(dot(lightDirectionEC, normalEC), 0.0); } float getLambertDiffuseOfMaterial(vec3 lightDirectionEC, czm_material material) { return getLambertDiffuse(lightDirectionEC, material.normal); } float getSpecular(vec3 lightDirectionEC, vec3 toEyeEC, vec3 normalEC, float shininess) { vec3 toReflectedLight = reflect(-lightDirectionEC, normalEC); float specular = max(dot(toReflectedLight, toEyeEC), 0.0); return pow(specular, shininess); } float getSpecularOfMaterial(vec3 lightDirectionEC, vec3 toEyeEC, czm_material material) { return getSpecular(lightDirectionEC, toEyeEC, material.normal, material.shininess); } /** * Computes a color using the Phong lighting model. * * @name czm_phong * @glslFunction * * @param {vec3} toEye A normalized vector from the fragment to the eye in eye coordinates. * @param {czm_material} material The fragment's material. * * @returns {vec4} The computed color. * * @example * vec3 positionToEyeEC = // ... * czm_material material = // ... * gl_FragColor = czm_phong(normalize(positionToEyeEC), material); * * @see czm_getMaterial */ vec4 czm_phong(vec3 toEye, czm_material material) { // Diffuse from directional light sources at eye (for top-down and horizon views) float diffuse = getLambertDiffuseOfMaterial(vec3(0.0, 0.0, 1.0), material) + getLambertDiffuseOfMaterial(vec3(0.0, 1.0, 0.0), material); // Specular from sun and pseudo-moon float specular = getSpecularOfMaterial(czm_sunDirectionEC, toEye, material) + getSpecularOfMaterial(czm_moonDirectionEC, toEye, material); vec3 ambient = vec3(0.0); vec3 color = ambient + material.emission; color += material.diffuse * diffuse; color += material.specular * specular; return vec4(color, material.alpha); } /** * Computes the luminance of a color. * * @name czm_luminance * @glslFunction * * @param {vec3} rgb The color. * * @returns {float} The luminance. * * @example * float light = czm_luminance(vec3(0.0)); // 0.0 * float dark = czm_luminance(vec3(1.0)); // ~1.0 */ float czm_luminance(vec3 rgb) { // Algorithm from Chapter 10 of Graphics Shaders. const vec3 W = vec3(0.2125, 0.7154, 0.0721); return dot(rgb, W); } /** * Adjusts the hue of a color. * * @name czm_hue * @glslFunction * * @param {vec3} rgb The color. * @param {float} adjustment The amount to adjust the hue of the color in radians. * * @returns {float} The color with the hue adjusted. * * @example * vec3 adjustHue = czm_hue(color, czm_pi); // The same as czm_hue(color, -czm_pi) */ vec3 czm_hue(vec3 rgb, float adjustment) { const mat3 toYIQ = mat3(0.299, 0.587, 0.114, 0.595716, -0.274453, -0.321263, 0.211456, -0.522591, 0.311135); const mat3 toRGB = mat3(1.0, 0.9563, 0.6210, 1.0, -0.2721, -0.6474, 1.0, -1.107, 1.7046); vec3 yiq = toYIQ * rgb; float hue = atan(yiq.z, yiq.y) + adjustment; float chroma = sqrt(yiq.z * yiq.z + yiq.y * yiq.y); vec3 color = vec3(yiq.x, chroma * cos(hue), chroma * sin(hue)); return toRGB * color; } /** * Adjusts the saturation of a color. * * @name czm_saturation * @glslFunction * * @param {vec3} rgb The color. * @param {float} adjustment The amount to adjust the saturation of the color. * * @returns {float} The color with the saturation adjusted. * * @example * vec3 greyScale = czm_saturation(color, 0.0); * vec3 doubleSaturation = czm_saturation(color, 2.0); */ vec3 czm_saturation(vec3 rgb, float adjustment) { // Algorithm from Chapter 16 of OpenGL Shading Language vec3 intensity = vec3(czm_luminance(rgb)); return mix(intensity, rgb, adjustment); } /** * DOC_TBA * * @name czm_multiplyWithColorBalance * @glslFunction */ vec3 czm_multiplyWithColorBalance(vec3 left, vec3 right) { // Algorithm from Chapter 10 of Graphics Shaders. const vec3 W = vec3(0.2125, 0.7154, 0.0721); vec3 target = left * right; float leftLuminance = dot(left, W); float rightLuminance = dot(right, W); float targetLuminance = dot(target, W); return ((leftLuminance + rightLuminance) / (2.0 * targetLuminance)) * target; } /** * Procedural anti-aliasing by blurring two colors that meet at a sharp edge. * * @name czm_antialias * @glslFunction * * @param {vec4} color1 The color on one side of the edge. * @param {vec4} color2 The color on the other side of the edge. * @param {vec4} currentcolor The current color, either <code>color1</code> or <code>color2</code>. * @param {float} dist The distance to the edge in texture coordinates. * @param {float} [fuzzFactor=0.1] Controls the blurriness between the two colors. * @returns {vec4} The anti-aliased color. * * @example * // GLSL declarations * vec4 czm_antialias(vec4 color1, vec4 color2, vec4 currentColor, float dist, float fuzzFactor); * vec4 czm_antialias(vec4 color1, vec4 color2, vec4 currentColor, float dist); * * // get the color for a material that has a sharp edge at the line y = 0.5 in texture space * float dist = abs(textureCoordinates.t - 0.5); * vec4 currentColor = mix(bottomColor, topColor, step(0.5, textureCoordinates.t)); * vec4 color = czm_antialias(bottomColor, topColor, currentColor, dist, 0.1); */ vec4 czm_antialias(vec4 color1, vec4 color2, vec4 currentColor, float dist, float fuzzFactor) { float val1 = clamp(dist / fuzzFactor, 0.0, 1.0); float val2 = clamp((dist - 0.5) / fuzzFactor, 0.0, 1.0); val1 = val1 * (1.0 - val2); val1 = val1 * val1 * (3.0 - (2.0 * val1)); val1 = pow(val1, 0.5); //makes the transition nicer vec4 midColor = (color1 + color2) * 0.5; return mix(midColor, currentColor, val1); } vec4 czm_antialias(vec4 color1, vec4 color2, vec4 currentColor, float dist) { return czm_antialias(color1, color2, currentColor, dist, 0.1); } /////////////////////////////////////////////////////////////////////////////// /** * The maximum latitude, in radians, both North and South, supported by a Web Mercator * (EPSG:3857) projection. Technically, the Mercator projection is defined * for any latitude up to (but not including) 90 degrees, but it makes sense * to cut it off sooner because it grows exponentially with increasing latitude. * The logic behind this particular cutoff value, which is the one used by * Google Maps, Bing Maps, and Esri, is that it makes the projection * square. That is, the extent is equal in the X and Y directions. * * The constant value is computed as follows: * czm_pi * 0.5 - (2.0 * atan(exp(-czm_pi))) * * @name czm_webMercatorMaxLatitude * @glslConstant */ const float czm_webMercatorMaxLatitude = 1.4844222297453323669610967939; /** * Specifies a flat, 2D map. * * @name czm_scene2D * @glslConstant */ const int czm_scene2D = 0; /** * Specifies 2.D Columbus View. * * @name czm_columbusView * @glslConstant */ const int czm_columbusView = 1; /** * Specifies a 3D globe. * * @name czm_scene3D * @glslConstant */ const int czm_scene3D = 2; /** * Specifies that the scene is morphing between modes. * * @name czm_morphing * @glslConstant */ const int czm_morphing = 3; /** * DOC_TBA * * @name czm_columbusViewMorph * @glslFunction */ vec4 czm_columbusViewMorph(vec3 position2D, vec3 position3D, float time) { // Just linear for now. vec3 p = mix(position2D, position3D, time); return vec4(p, 1.0); } /////////////////////////////////////////////////////////////////////////////// /** * DOC_TBA * * @name czm_ray * @glslStruct */ struct czm_ray { vec3 origin; vec3 direction; }; /** * Computes the point along a ray at the given time. <code>time</code> can be positive, negative, or zero. * * @name czm_pointAlongRay * @glslFunction * * @param {czm_ray} ray The ray to compute the point along. * @param {float} time The time along the ray. * * @returns {vec3} The point along the ray at the given time. * * @example * czm_ray ray = czm_ray(vec3(0.0), vec3(1.0, 0.0, 0.0)); // origin, direction * vec3 v = czm_pointAlongRay(ray, 2.0); // (2.0, 0.0, 0.0) */ vec3 czm_pointAlongRay(czm_ray ray, float time) { return ray.origin + (time * ray.direction); } /////////////////////////////////////////////////////////////////////////////// /** * DOC_TBA * * @name czm_raySegment * @glslStruct */ struct czm_raySegment { float start; float stop; }; /** * DOC_TBA * * @name czm_emptyRaySegment * @glslConstant */ const czm_raySegment czm_emptyRaySegment = czm_raySegment(-czm_infinity, -czm_infinity); /** * DOC_TBA * * @name czm_fullRaySegment * @glslConstant */ const czm_raySegment czm_fullRaySegment = czm_raySegment(0.0, czm_infinity); /** * Determines if a time interval is empty. * * @name czm_isEmpty * @glslFunction * * @param {czm_raySegment} interval The interval to test. * * @returns {bool} <code>true</code> if the time interval is empty; otherwise, <code>false</code>. * * @example * bool b0 = czm_isEmpty(czm_emptyRaySegment); // true * bool b1 = czm_isEmpty(czm_raySegment(0.0, 1.0)); // false * bool b2 = czm_isEmpty(czm_raySegment(1.0, 1.0)); // false, contains 1.0. */ bool czm_isEmpty(czm_raySegment interval) { return (interval.stop < 0.0); } /** * Determines if a time interval is empty. * * @name czm_isFull * @glslFunction * * @param {czm_raySegment} interval The interval to test. * * @returns {bool} <code>true</code> if the time interval is empty; otherwise, <code>false</code>. * * @example * bool b0 = czm_isEmpty(czm_emptyRaySegment); // true * bool b1 = czm_isEmpty(czm_raySegment(0.0, 1.0)); // false * bool b2 = czm_isEmpty(czm_raySegment(1.0, 1.0)); // false, contains 1.0. */ bool czm_isFull(czm_raySegment interval) { return (interval.start == 0.0 && interval.stop == czm_infinity); } /////////////////////////////////////////////////////////////////////////////// /** * DOC_TBA * * @name czm_ellipsoid * @glslStruct */ struct czm_ellipsoid { vec3 center; vec3 radii; vec3 inverseRadii; vec3 inverseRadiiSquared; }; /** * DOC_TBA * * @name czm_ellipsoidNew * @glslFunction * */ czm_ellipsoid czm_ellipsoidNew(vec3 center, vec3 radii) { vec3 inverseRadii = vec3(1.0 / radii.x, 1.0 / radii.y, 1.0 / radii.z); vec3 inverseRadiiSquared = inverseRadii * inverseRadii; czm_ellipsoid temp = czm_ellipsoid(center, radii, inverseRadii, inverseRadiiSquared); return temp; } /** * DOC_TBA * * @name czm_ellipsoidContainsPoint * @glslFunction * */ bool czm_ellipsoidContainsPoint(czm_ellipsoid ellipsoid, vec3 point) { vec3 scaled = ellipsoid.inverseRadii * (czm_inverseModelView * vec4(point, 1.0)).xyz; return (dot(scaled, scaled) <= 1.0); } /** * DOC_TBA * * * @name czm_rayEllipsoidIntersectionInterval * @glslFunction */ czm_raySegment czm_rayEllipsoidIntersectionInterval(czm_ray ray, czm_ellipsoid ellipsoid) { // ray and ellipsoid center in eye coordinates. radii in model coordinates. vec3 q = ellipsoid.inverseRadii * (czm_inverseModelView * vec4(ray.origin, 1.0)).xyz; vec3 w = ellipsoid.inverseRadii * (czm_inverseModelView * vec4(ray.direction, 0.0)).xyz; q = q - ellipsoid.inverseRadii * (czm_inverseModelView * vec4(ellipsoid.center, 1.0)).xyz; float q2 = dot(q, q); float qw = dot(q, w); if (q2 > 1.0) // Outside ellipsoid. { if (qw >= 0.0) // Looking outward or tangent (0 intersections). { return czm_emptyRaySegment; } else // qw < 0.0. { float qw2 = qw * qw; float difference = q2 - 1.0; // Positively valued. float w2 = dot(w, w); float product = w2 * difference; if (qw2 < product) // Imaginary roots (0 intersections). { return czm_emptyRaySegment; } else if (qw2 > product) // Distinct roots (2 intersections). { float discriminant = qw * qw - product; float temp = -qw + sqrt(discriminant); // Avoid cancellation. float root0 = temp / w2; float root1 = difference / temp; if (root0 < root1) { czm_raySegment i = czm_raySegment(root0, root1); return i; } else { czm_raySegment i = czm_raySegment(root1, root0); return i; } } else // qw2 == product. Repeated roots (2 intersections). { float root = sqrt(difference / w2); czm_raySegment i = czm_raySegment(root, root); return i; } } } else if (q2 < 1.0) // Inside ellipsoid (2 intersections). { float difference = q2 - 1.0; // Negatively valued. float w2 = dot(w, w); float product = w2 * difference; // Negatively valued. float discriminant = qw * qw - product; float temp = -qw + sqrt(discriminant); // Positively valued. czm_raySegment i = czm_raySegment(0.0, temp / w2); return i; } else // q2 == 1.0. On ellipsoid. { if (qw < 0.0) // Looking inward. { float w2 = dot(w, w); czm_raySegment i = czm_raySegment(0.0, -qw / w2); return i; } else // qw >= 0.0. Looking outward or tangent. { return czm_emptyRaySegment; } } } /** * Returns the WGS84 ellipsoid, with its center at the origin of world coordinates, in eye coordinates. * * @name czm_getWgs84EllipsoidEC * @glslFunction * * @returns {czm_ellipsoid} The WGS84 ellipsoid, with its center at the origin of world coordinates, in eye coordinates. * * @see Ellipsoid.getWgs84 * * @example * czm_ellipsoid ellipsoid = czm_getWgs84EllipsoidEC(); */ czm_ellipsoid czm_getWgs84EllipsoidEC() { return czm_ellipsoidNew( czm_view[3].xyz, // center vec3(6378137.0, 6378137.0, 6356752.314245)); // radii } /** * Computes the fraction of a Web Wercator extent at which a given geodetic latitude is located. * * @name czm_latitudeToWebMercatorFraction * @glslFunction * * @param {float} The geodetic latitude, in radians. * @param {float} The low portion of the Web Mercator coordinate of the southern boundary of the extent. * @param {float} The high portion of the Web Mercator coordinate of the southern boundary of the extent. * @param {float} The total height of the extent in Web Mercator coordinates. * * @returns {float} The fraction of the extent at which the latitude occurs. If the latitude is the southern * boundary of the extent, the return value will be zero. If it is the northern boundary, the return * value will be 1.0. Latitudes in between are mapped according to the Web Mercator projection. */ float czm_latitudeToWebMercatorFraction(float latitude, float southMercatorYLow, float southMercatorYHigh, float oneOverMercatorHeight) { float sinLatitude = sin(latitude); float mercatorY = 0.5 * log((1.0 + sinLatitude) / (1.0 - sinLatitude)); // mercatorY - southMercatorY in simulated double precision. float t1 = 0.0 - southMercatorYLow; float e = t1 - 0.0; float t2 = ((-southMercatorYLow - e) + (0.0 - (t1 - e))) + mercatorY - southMercatorYHigh; float highDifference = t1 + t2; float lowDifference = t2 - (highDifference - t1); return highDifference * oneOverMercatorHeight + lowDifference * oneOverMercatorHeight; } /** * Translates a position (or any <code>vec3</code>) that was encoded with {@link EncodedCartesian3}, * and then provided to the shader as separate <code>high</code> and <code>low</code> bits to * be relative to the eye. As shown in the example, the position can then be transformed in eye * or clip coordinates using {@link czm_modelViewRelativeToEye} or {@link czm_modelViewProjectionRelativeToEye}, * respectively. * <p> * This technique, called GPU RTE, eliminates jittering artifacts when using large coordinates as * described in <a href="http://blogs.agi.com/insight3d/index.php/2008/09/03/precisions-precisions/">Precisions, Precisions</a>. * </p> * * @name czm_translateRelativeToEye * @glslFunction * * @param {vec3} high The position's high bits. * @param {vec3} low The position's low bits. * @returns {vec3} The position translated to be relative to the camera's position. * * @example * attribute vec3 positionHigh; * attribute vec3 positionLow; * * void main() * { * vec3 p = czm_translateRelativeToEye(positionHigh, positionLow); * gl_Position = czm_modelViewProjectionRelativeToEye * vec4(p, 1.0); * } * * @see czm_modelViewRelativeToEye * @see czm_modelViewProjectionRelativeToEye * @see EncodedCartesian3 */ vec3 czm_translateRelativeToEye(vec3 high, vec3 low) { vec3 highDifference = high - czm_encodedCameraPositionMCHigh; vec3 lowDifference = low - czm_encodedCameraPositionMCLow; return highDifference + lowDifference; } /** * @private */ vec4 czm_getWaterNoise(sampler2D normalMap, vec2 uv, float time, float angleInRadians) { float cosAngle = cos(angleInRadians); float sinAngle = sin(angleInRadians); // time dependent sampling directions vec2 s0 = vec2(1.0/17.0, 0.0); vec2 s1 = vec2(-1.0/29.0, 0.0); vec2 s2 = vec2(1.0/101.0, 1.0/59.0); vec2 s3 = vec2(-1.0/109.0, -1.0/57.0); // rotate sampling direction by specified angle s0 = vec2((cosAngle * s0.x) - (sinAngle * s0.y), (sinAngle * s0.x) + (cosAngle * s0.y)); s1 = vec2((cosAngle * s1.x) - (sinAngle * s1.y), (sinAngle * s1.x) + (cosAngle * s1.y)); s2 = vec2((cosAngle * s2.x) - (sinAngle * s2.y), (sinAngle * s2.x) + (cosAngle * s2.y)); s3 = vec2((cosAngle * s3.x) - (sinAngle * s3.y), (sinAngle * s3.x) + (cosAngle * s3.y)); vec2 uv0 = (uv/103.0) + (time * s0); vec2 uv1 = uv/107.0 + (time * s1) + vec2(0.23); vec2 uv2 = uv/vec2(897.0, 983.0) + (time * s2) + vec2(0.51); vec2 uv3 = uv/vec2(991.0, 877.0) + (time * s3) + vec2(0.71); uv0 = fract(uv0); uv1 = fract(uv1); uv2 = fract(uv2); uv3 = fract(uv3); vec4 noise = (texture2D(normalMap, uv0)) + (texture2D(normalMap, uv1)) + (texture2D(normalMap, uv2)) + (texture2D(normalMap, uv3)); // average and scale to between -1 and 1 return ((noise / 4.0) - 0.5) * 2.0; } /** * @license * Description : Array and textureless GLSL 2D/3D/4D simplex * noise functions. * Author : Ian McEwan, Ashima Arts. * Maintainer : ijm * Lastmod : 20110822 (ijm) * License : Copyright (C) 2011 Ashima Arts. All rights reserved. * Distributed under the MIT License. See LICENSE file. * https://github.com/ashima/webgl-noise */ vec4 _czm_mod289(vec4 x) { return x - floor(x * (1.0 / 289.0)) * 289.0; } vec3 _czm_mod289(vec3 x) { return x - floor(x * (1.0 / 289.0)) * 289.0; } vec2 _czm_mod289(vec2 x) { return x - floor(x * (1.0 / 289.0)) * 289.0; } float _czm_mod289(float x) { return x - floor(x * (1.0 / 289.0)) * 289.0; } vec4 _czm_permute(vec4 x) { return _czm_mod289(((x*34.0)+1.0)*x); } vec3 _czm_permute(vec3 x) { return _czm_mod289(((x*34.0)+1.0)*x); } float _czm_permute(float x) { return _czm_mod289(((x*34.0)+1.0)*x); } vec4 _czm_taylorInvSqrt(vec4 r) { return 1.79284291400159 - 0.85373472095314 * r; } float _czm_taylorInvSqrt(float r) { return 1.79284291400159 - 0.85373472095314 * r; } vec4 _czm_grad4(float j, vec4 ip) { const vec4 ones = vec4(1.0, 1.0, 1.0, -1.0); vec4 p,s; p.xyz = floor( fract (vec3(j) * ip.xyz) * 7.0) * ip.z - 1.0; p.w = 1.5 - dot(abs(p.xyz), ones.xyz); s = vec4(lessThan(p, vec4(0.0))); p.xyz = p.xyz + (s.xyz*2.0 - 1.0) * s.www; return p; } /** * DOC_TBA * * Implemented by Ian McEwan, Ashima Arts, and distributed under the MIT License. {@link https://github.com/ashima/webgl-noise} * * @name czm_snoise * @glslFunction * * @see <a href="https://github.com/ashima/webgl-noise">https://github.com/ashima/webgl-noise</a> * @see Stefan Gustavson's paper <a href="http://www.itn.liu.se/~stegu/simplexnoise/simplexnoise.pdf">Simplex noise demystified</a> */ float czm_snoise(vec2 v) { const vec4 C = vec4(0.211324865405187, // (3.0-sqrt(3.0))/6.0 0.366025403784439, // 0.5*(sqrt(3.0)-1.0) -0.577350269189626, // -1.0 + 2.0 * C.x 0.024390243902439); // 1.0 / 41.0 // First corner vec2 i = floor(v + dot(v, C.yy) ); vec2 x0 = v - i + dot(i, C.xx); // Other corners vec2 i1; //i1.x = step( x0.y, x0.x ); // x0.x > x0.y ? 1.0 : 0.0 //i1.y = 1.0 - i1.x; i1 = (x0.x > x0.y) ? vec2(1.0, 0.0) : vec2(0.0, 1.0); // x0 = x0 - 0.0 + 0.0 * C.xx ; // x1 = x0 - i1 + 1.0 * C.xx ; // x2 = x0 - 1.0 + 2.0 * C.xx ; vec4 x12 = x0.xyxy + C.xxzz; x12.xy -= i1; // Permutations i = _czm_mod289(i); // Avoid truncation effects in permutation vec3 p = _czm_permute( _czm_permute( i.y + vec3(0.0, i1.y, 1.0 )) + i.x + vec3(0.0, i1.x, 1.0 )); vec3 m = max(0.5 - vec3(dot(x0,x0), dot(x12.xy,x12.xy), dot(x12.zw,x12.zw)), 0.0); m = m*m ; m = m*m ; // Gradients: 41 points uniformly over a line, mapped onto a diamond. // The ring size 17*17 = 289 is close to a multiple of 41 (41*7 = 287) vec3 x = 2.0 * fract(p * C.www) - 1.0; vec3 h = abs(x) - 0.5; vec3 ox = floor(x + 0.5); vec3 a0 = x - ox; // Normalise gradients implicitly by scaling m // Approximation of: m *= inversesqrt( a0*a0 + h*h ); m *= 1.79284291400159 - 0.85373472095314 * ( a0*a0 + h*h ); // Compute final noise value at P vec3 g; g.x = a0.x * x0.x + h.x * x0.y; g.yz = a0.yz * x12.xz + h.yz * x12.yw; return 130.0 * dot(m, g); } float czm_snoise(vec3 v) { const vec2 C = vec2(1.0/6.0, 1.0/3.0) ; const vec4 D = vec4(0.0, 0.5, 1.0, 2.0); // First corner vec3 i = floor(v + dot(v, C.yyy) ); vec3 x0 = v - i + dot(i, C.xxx) ; // Other corners vec3 g = step(x0.yzx, x0.xyz); vec3 l = 1.0 - g; vec3 i1 = min( g.xyz, l.zxy ); vec3 i2 = max( g.xyz, l.zxy ); // x0 = x0 - 0.0 + 0.0 * C.xxx; // x1 = x0 - i1 + 1.0 * C.xxx; // x2 = x0 - i2 + 2.0 * C.xxx; // x3 = x0 - 1.0 + 3.0 * C.xxx; vec3 x1 = x0 - i1 + C.xxx; vec3 x2 = x0 - i2 + C.yyy; // 2.0*C.x = 1/3 = C.y vec3 x3 = x0 - D.yyy; // -1.0+3.0*C.x = -0.5 = -D.y // Permutations i = _czm_mod289(i); vec4 p = _czm_permute( _czm_permute( _czm_permute( i.z + vec4(0.0, i1.z, i2.z, 1.0 )) + i.y + vec4(0.0, i1.y, i2.y, 1.0 )) + i.x + vec4(0.0, i1.x, i2.x, 1.0 )); // Gradients: 7x7 points over a square, mapped onto an octahedron. // The ring size 17*17 = 289 is close to a multiple of 49 (49*6 = 294) float n_ = 0.142857142857; // 1.0/7.0 vec3 ns = n_ * D.wyz - D.xzx; vec4 j = p - 49.0 * floor(p * ns.z * ns.z); // mod(p,7*7) vec4 x_ = floor(j * ns.z); vec4 y_ = floor(j - 7.0 * x_ ); // mod(j,N) vec4 x = x_ *ns.x + ns.yyyy; vec4 y = y_ *ns.x + ns.yyyy; vec4 h = 1.0 - abs(x) - abs(y); vec4 b0 = vec4( x.xy, y.xy ); vec4 b1 = vec4( x.zw, y.zw ); //vec4 s0 = vec4(lessThan(b0,0.0))*2.0 - 1.0; //vec4 s1 = vec4(lessThan(b1,0.0))*2.0 - 1.0; vec4 s0 = floor(b0)*2.0 + 1.0; vec4 s1 = floor(b1)*2.0 + 1.0; vec4 sh = -step(h, vec4(0.0)); vec4 a0 = b0.xzyw + s0.xzyw*sh.xxyy ; vec4 a1 = b1.xzyw + s1.xzyw*sh.zzww ; vec3 p0 = vec3(a0.xy,h.x); vec3 p1 = vec3(a0.zw,h.y); vec3 p2 = vec3(a1.xy,h.z); vec3 p3 = vec3(a1.zw,h.w); //Normalise gradients vec4 norm = _czm_taylorInvSqrt(vec4(dot(p0,p0), dot(p1,p1), dot(p2, p2), dot(p3,p3))); p0 *= norm.x; p1 *= norm.y; p2 *= norm.z; p3 *= norm.w; // Mix final noise value vec4 m = max(0.6 - vec4(dot(x0,x0), dot(x1,x1), dot(x2,x2), dot(x3,x3)), 0.0); m = m * m; return 42.0 * dot( m*m, vec4( dot(p0,x0), dot(p1,x1), dot(p2,x2), dot(p3,x3) ) ); } float czm_snoise(vec4 v) { const vec4 C = vec4( 0.138196601125011, // (5 - sqrt(5))/20 G4 0.276393202250021, // 2 * G4 0.414589803375032, // 3 * G4 -0.447213595499958); // -1 + 4 * G4 // (sqrt(5) - 1)/4 = F4, used once below #define F4 0.309016994374947451 // First corner vec4 i = floor(v + dot(v, vec4(F4)) ); vec4 x0 = v - i + dot(i, C.xxxx); // Other corners // Rank sorting originally contributed by Bill Licea-Kane, AMD (formerly ATI) vec4 i0; vec3 isX = step( x0.yzw, x0.xxx ); vec3 isYZ = step( x0.zww, x0.yyz ); // i0.x = dot( isX, vec3( 1.0 ) ); i0.x = isX.x + isX.y + isX.z; i0.yzw = 1.0 - isX; // i0.y += dot( isYZ.xy, vec2( 1.0 ) ); i0.y += isYZ.x + isYZ.y; i0.zw += 1.0 - isYZ.xy; i0.z += isYZ.z; i0.w += 1.0 - isYZ.z; // i0 now contains the unique values 0,1,2,3 in each channel vec4 i3 = clamp( i0, 0.0, 1.0 ); vec4 i2 = clamp( i0-1.0, 0.0, 1.0 ); vec4 i1 = clamp( i0-2.0, 0.0, 1.0 ); // x0 = x0 - 0.0 + 0.0 * C.xxxx // x1 = x0 - i1 + 1.0 * C.xxxx // x2 = x0 - i2 + 2.0 * C.xxxx // x3 = x0 - i3 + 3.0 * C.xxxx // x4 = x0 - 1.0 + 4.0 * C.xxxx vec4 x1 = x0 - i1 + C.xxxx; vec4 x2 = x0 - i2 + C.yyyy; vec4 x3 = x0 - i3 + C.zzzz; vec4 x4 = x0 + C.wwww; // Permutations i = _czm_mod289(i); float j0 = _czm_permute( _czm_permute( _czm_permute( _czm_permute(i.w) + i.z) + i.y) + i.x); vec4 j1 = _czm_permute( _czm_permute( _czm_permute( _czm_permute ( i.w + vec4(i1.w, i2.w, i3.w, 1.0 )) + i.z + vec4(i1.z, i2.z, i3.z, 1.0 )) + i.y + vec4(i1.y, i2.y, i3.y, 1.0 )) + i.x + vec4(i1.x, i2.x, i3.x, 1.0 )); // Gradients: 7x7x6 points over a cube, mapped onto a 4-cross polytope // 7*7*6 = 294, which is close to the ring size 17*17 = 289. vec4 ip = vec4(1.0/294.0, 1.0/49.0, 1.0/7.0, 0.0) ; vec4 p0 = _czm_grad4(j0, ip); vec4 p1 = _czm_grad4(j1.x, ip); vec4 p2 = _czm_grad4(j1.y, ip); vec4 p3 = _czm_grad4(j1.z, ip); vec4 p4 = _czm_grad4(j1.w, ip); // Normalise gradients vec4 norm = _czm_taylorInvSqrt(vec4(dot(p0,p0), dot(p1,p1), dot(p2, p2), dot(p3,p3))); p0 *= norm.x; p1 *= norm.y; p2 *= norm.z; p3 *= norm.w; p4 *= _czm_taylorInvSqrt(dot(p4,p4)); // Mix contributions from the five corners vec3 m0 = max(0.6 - vec3(dot(x0,x0), dot(x1,x1), dot(x2,x2)), 0.0); vec2 m1 = max(0.6 - vec2(dot(x3,x3), dot(x4,x4) ), 0.0); m0 = m0 * m0; m1 = m1 * m1; return 49.0 * ( dot(m0*m0, vec3( dot( p0, x0 ), dot( p1, x1 ), dot( p2, x2 ))) + dot(m1*m1, vec2( dot( p3, x3 ), dot( p4, x4 ) ) ) ) ; } /////////////////////////////////////////////////////////////////////////////// /** * @license * Cellular noise ("Worley noise") in 2D in GLSL. * Copyright (c) Stefan Gustavson 2011-04-19. All rights reserved. * This code is released under the conditions of the MIT license. * See LICENSE file for details. */ //#ifdef GL_OES_standard_derivatives // #extension GL_OES_standard_derivatives : enable //#endif // //float aastep (float threshold , float value) //{ // float afwidth = 0.7 * length ( vec2 ( dFdx ( value ), dFdy ( value ))); // return smoothstep ( threshold - afwidth , threshold + afwidth , value ); //} // Permutation polynomial: (34x^2 + x) mod 289 vec3 _czm_permute289(vec3 x) { return mod((34.0 * x + 1.0) * x, 289.0); } /** * DOC_TBA * * Implemented by Stefan Gustavson, and distributed under the MIT License. {@link http://openglinsights.git.sourceforge.net/git/gitweb.cgi?p=openglinsights/openglinsights;a=tree;f=proceduraltextures} * * @name czm_cellular * @glslFunction * * @see Stefan Gustavson's chapter, <i>Procedural Textures in GLSL</i>, in <a href="http://www.openglinsights.com/">OpenGL Insights</a>. */ vec2 czm_cellular(vec2 P) // Cellular noise, returning F1 and F2 in a vec2. // Standard 3x3 search window for good F1 and F2 values { #define K 0.142857142857 // 1/7 #define Ko 0.428571428571 // 3/7 #define jitter 1.0 // Less gives more regular pattern vec2 Pi = mod(floor(P), 289.0); vec2 Pf = fract(P); vec3 oi = vec3(-1.0, 0.0, 1.0); vec3 of = vec3(-0.5, 0.5, 1.5); vec3 px = _czm_permute289(Pi.x + oi); vec3 p = _czm_permute289(px.x + Pi.y + oi); // p11, p12, p13 vec3 ox = fract(p*K) - Ko; vec3 oy = mod(floor(p*K),7.0)*K - Ko; vec3 dx = Pf.x + 0.5 + jitter*ox; vec3 dy = Pf.y - of + jitter*oy; vec3 d1 = dx * dx + dy * dy; // d11, d12 and d13, squared p = _czm_permute289(px.y + Pi.y + oi); // p21, p22, p23 ox = fract(p*K) - Ko; oy = mod(floor(p*K),7.0)*K - Ko; dx = Pf.x - 0.5 + jitter*ox; dy = Pf.y - of + jitter*oy; vec3 d2 = dx * dx + dy * dy; // d21, d22 and d23, squared p = _czm_permute289(px.z + Pi.y + oi); // p31, p32, p33 ox = fract(p*K) - Ko; oy = mod(floor(p*K),7.0)*K - Ko; dx = Pf.x - 1.5 + jitter*ox; dy = Pf.y - of + jitter*oy; vec3 d3 = dx * dx + dy * dy; // d31, d32 and d33, squared // Sort out the two smallest distances (F1, F2) vec3 d1a = min(d1, d2); d2 = max(d1, d2); // Swap to keep candidates for F2 d2 = min(d2, d3); // neither F1 nor F2 are now in d3 d1 = min(d1a, d2); // F1 is now in d1 d2 = max(d1a, d2); // Swap to keep candidates for F2 d1.xy = (d1.x < d1.y) ? d1.xy : d1.yx; // Swap if smaller d1.xz = (d1.x < d1.z) ? d1.xz : d1.zx; // F1 is in d1.x d1.yz = min(d1.yz, d2.yz); // F2 is now not in d2.yz d1.y = min(d1.y, d1.z); // nor in d1.z d1.y = min(d1.y, d2.x); // F2 is in d1.y, we're done. return sqrt(d1.xy); } /* // Cellular noise, returning F1 and F2 in a vec2 and the // 2D vectors to each of the two closest points in a vec4. // Standard 3x3 search window for good F1 and F2 values. void czm_cellular(in vec2 P, out vec2 F, out vec4 d1d2) { #define K 0.142857142857 // 1/7 #define Ko 0.428571428571 // 3/7 #define jitter 1.0 // Less gives more regular pattern vec2 Pi = mod(floor(P), 289.0); vec2 Pf = fract(P); vec3 oi = vec3(-1.0, 0.0, 1.0); vec3 of = vec3(-0.5, 0.5, 1.5); vec3 px = _czm_permute289(Pi.x + oi); vec3 p = _czm_permute289(px.x + Pi.y + oi); // p11, p12, p13 vec3 ox = fract(p*K) - Ko; vec3 oy = mod(floor(p*K),7.0)*K - Ko; vec3 d1x = Pf.x + 0.5 + jitter*ox; vec3 d1y = Pf.y - of + jitter*oy; vec3 d1 = d1x * d1x + d1y * d1y; // d11, d12 and d13, squared p = _czm_permute289(px.y + Pi.y + oi); // p21, p22, p23 ox = fract(p*K) - Ko; oy = mod(floor(p*K),7.0)*K - Ko; vec3 d2x = Pf.x - 0.5 + jitter*ox; vec3 d2y = Pf.y - of + jitter*oy; vec3 d2 = d2x * d2x + d2y * d2y; // d21, d22 and d23, squared p = _czm_permute289(px.z + Pi.y + oi); // p31, p32, p33 ox = fract(p*K) - Ko; oy = mod(floor(p*K),7.0)*K - Ko; vec3 d3x = Pf.x - 1.5 + jitter*ox; vec3 d3y = Pf.y - of + jitter*oy; vec3 d3 = d3x * d3x + d3y * d3y; // d31, d32 and d33, squared // Sort out the two smallest distances (F1, F2) // While also swapping dx and dy accordingly vec3 comp3 = step(d2, d1); vec3 d1a = mix(d1, d2, comp3); vec3 d1xa = mix(d1x, d2x, comp3); vec3 d1ya = mix(d1y, d2y, comp3); d2 = mix(d2, d1, comp3); // Swap to keep candidates for F2 d2x = mix(d2x, d1x, comp3); d2y = mix(d2y, d1y, comp3); comp3 = step(d3, d2); d2 = mix(d2, d3, comp3); // neither F1 nor F2 are now in d3 d2x = mix(d2x, d3x, comp3); d2y = mix(d2y, d3y, comp3); comp3 = step(d2, d1a); d1 = mix(d1a, d2, comp3); // F1 is now in d1 d1x = mix(d1xa, d2x, comp3); d1y = mix(d1ya, d2y, comp3); d2 = mix(d2, d1a, comp3); // Swap to keep candidates for F2 d2x = mix(d2x, d1xa, comp3); d2y = mix(d2y, d1ya, comp3); float comp1 = step(d1.y, d1.x); d1.xy = mix(d1.xy, d1.yx, comp1); // Swap if smaller d1x.xy = mix(d1x.xy, d1x.yx, comp1); d1y.xy = mix(d1y.xy, d1y.yx, comp1); comp1 = step(d1.z, d1.x); d1.xz = mix(d1.xz, d1.zx, comp1); // F1 is in d1.x d1x.xz = mix(d1x.xz, d1x.zx, comp1); d1y.xz = mix(d1y.xz, d1y.zx, comp1); vec2 comp2 = step(d2.yz, d1.yz); d1.yz = mix(d1.yz, d2.yz, comp2); // F2 is now not in d2.yz d1x.yz = mix(d1x.yz, d2x.yz, comp2); d1y.yz = mix(d1y.yz, d2y.yz, comp2); comp1 = step(d1.z, d1.y); d1.y = mix(d1.y, d1.z, comp1); // nor in d1.z d1x.y = mix(d1x.y, d1x.z, comp1); d1y.y = mix(d1y.y, d1y.z, comp1); comp1 = step(d2.x, d1.y); d1.y = mix(d1.y, d2.x, comp1); // F2 is in d1.y, we're done. d1x.y = mix(d1x.y, d2x.x, comp1); d1y.y = mix(d1y.y, d2y.x, comp1); F = sqrt(d1.xy); d1d2 = vec4(d1x.x, d1y.x, d1x.y, d1y.y); } */