31:T5f4,md`# Gaussian curvature in triangulated 2-Manifolds. In this notebook I explore the gaussian curvature for some two dimensional surfaces embedded in a 3D space. The surfaces are parametrised as function ${tex`r(u, v)`} of coordinates ${tex`u`} and ${tex`v`}: ${tex`r(u, v) \equiv (x(u, v), y(u, v), z(u, v))`}. The exact analytic calculation of the Gaussian curvature may be done with the coefficients of the first and second fundamental forms, determining the metric properties of a surface (see for example Gray A. in *Modern Differential Geometry of Curves and Surfaces with Mathematica*). Hence we need the gradient and hessian of function ${tex`r`}. It would be interesting to numerically calculate the first and second derivatives of ${tex`r(u, v)`} in each vertex using the hyper-dual numbers [proposed](http://adl.stanford.edu/hyperdual/) by Jeffrey A. Fike and Juan J. Alonso. Instead, in this notebook, I implemented the *discrete* Gaussian curvature as defined by Meyer et al. 2003 in *Discrete Differential-Geometry Operators for Triangulated 2-Manifolds*. The same concept is described by S. L. Devadoss and J. O'Rourke in *Discrete and Computational Geometry*. For the implementation, please have look at the vertex shader in function \`drawSurface\` below. Note that the implementation calculates the position of the neighbourhood vertices through i) the grid mesh created by the function \`createMesh\`, which sets the values of coordinates ${tex`u`} and ${tex`v`}, and ii) the function ${tex`r(u, v)`}. `32:T89a,viewof controls = { const canvas = createCanvas(); canvas.value = {eye: [0, 0, 1], viewingVolumeFactor: 4}; let radius = 1; let theta = Math.PI / 3; // Angle with respect the y axis. let phi = Math.PI / 3; // Angle between the x axis and the projection in the (x,z) plane. let velocity = 0.03; // Speed of the rotation. let viewingVolumeFactor = 15; // Scaling factor used to scale the viewing volume, zooming in/out the scene. let minViewingVolumeFactor = 1.; let maxViewingVolumeFactor = 20.; let previuosTheta = theta; let previousPhi = phi; let isMoving = false; let initialPositions = [0, 0]; // Calculate cartesian coordinates from spherical coordinates. const cartesian = (radius, theta, phi) => [ radius * Math.sin(theta) * Math.cos(phi), radius * Math.cos(theta), radius * Math.sin(theta) * Math.sin(phi) ]; const updateValue = () => { canvas.value.eye = cartesian(radius, theta, phi); canvas.value.viewingVolumeFactor = viewingVolumeFactor; // Emit input event to update variable `controls`. // See https://observablehq.com/@observablehq/introduction-to-views canvas.dispatchEvent(new CustomEvent("input")); } canvas.ontouchstart = canvas.onmousedown = event => { isMoving = true; initialPositions = [event.offsetX, event.offsetY]; }; window.ontouchend = window.onmouseup = () => { isMoving = false; previuosTheta = theta; previousPhi = phi; }; canvas.ontouchmove = canvas.onmousemove = event => { if (isMoving) { theta = clamp( previuosTheta - (event.offsetY - initialPositions[1]) / (Math.PI / velocity), Math.PI / 256, Math.PI / 256 * 255 ); phi = previousPhi - (event.offsetX - initialPositions[0]) / (Math.PI * 2 / velocity); updateValue(); } } canvas.onwheel = event => { event.preventDefault(); viewingVolumeFactor = clamp( viewingVolumeFactor + event.deltaY * 0.1, minViewingVolumeFactor, maxViewingVolumeFactor ); updateValue(); } updateValue(); // Update eye to the initial values of theta and phi. return canvas; }33:T460,kleinBottle = ({ r: ` vec3 r (vec2 uv) { float u = uv.x; float v = uv.y; // Parametrisation described in "Modern Differential Geometry of Curves and Surfaces with Mathematica" // by Gray, A. float a = 2.; float x = (a + cos(u/2.) * sin(v) - sin(u/2.) * sin(2.*v)) * cos(u); float z = (a + cos(u/2.) * sin(v) - sin(u/2.) * sin(2.*v)) * sin(u); float y = sin(u/2.) * sin(v) + cos(u/2.) * sin(2.*v); /* // From wikipedia float x = -2./15. * cos(u) * (3. * cos(v) - 30. * sin(u) + 90. * pow(cos(u), 4.) * sin(u) - 60. * pow(cos(u), 6.) * sin(u) + 5. * cos(u) * cos(v) * sin(u)); float z = -1./15. * sin(u) * (3. * cos(v) - 3. * pow(cos(u), 2.) * cos(v) - 48. * pow(cos(u), 4.) * cos(v) + 48. * pow(cos(u), 6.) * cos(v) - 60. * sin(u) + 5. * cos(u) * cos(v) * sin(u) - 5. * pow(cos(u), 3.) * cos(v) * sin(u) - 80. * pow(cos(u), 5.) * cos(v) * sin(u) + 80. * pow(cos(u), 7.) * cos(v) * sin(u)); float y = 2./15. * (3. + 5. * cos(u) * sin(u)) * sin(v); */ return vec3(x, y, z); } `, uDomain: [-Math.PI, Math.PI], vDomain: [-Math.PI, Math.PI] })34:T1824,drawSurface = gl.regl({ frag:` #extension GL_OES_standard_derivatives : enable precision highp float; varying vec4 v; varying float curvature; vec3 normal(vec3 pos) { vec3 fdx = dFdx(pos); vec3 fdy = dFdy(pos); return normalize(cross(fdx, fdy)); } // https://github.com/rreusser/glsl-solid-wireframe float grid (vec2 parameter, float width, float feather) { float w1 = width - feather * .5; vec2 d = fwidth(parameter); vec2 looped = 0.5 - abs(mod(parameter, 1.0) - 0.5); vec2 a2 = smoothstep(d * w1, d * (w1 + feather), looped); return min(a2.x, a2.y); } vec3 multiplyBlendMode (vec3 b, vec3 s) { return b * s; } vec3 screenBlendMode (vec3 b, vec3 s) { return b + s - (b * s); } vec3 complement(vec3 c) { return vec3(1.) - c; } vec4 generalCompositingFormula (vec4 b, vec4 s, float Fa, float Fb) { return vec4( s.a * Fa * s.rgb + b.a * Fb * b.rgb, s.a * Fa + b.a * Fb ); } vec4 sourceOverCompositing (vec4 b, vec4 s) { return generalCompositingFormula(b, s, 1., 1. - s.a); } vec4 simpleAlphaCompositing (vec4 b, vec4 s) { return sourceOverCompositing(b, s); } vec4 destinationOverCompositing (vec4 b, vec4 s) { return generalCompositingFormula(b, s, 1. - b.a, 1.); } float smootherstep(float edge0, float edge1, float x) { x = clamp((x - edge0) / (edge1 - edge0), 0.0, 1.0); return x * x * x * (x * (x * 6. - 15.) + 10.); } void main() { float e = smoothstep(0., 2., ${e.toFixed(2)}); vec4 b = vec4(mix( vec3(1., .0, .0), vec3(0., 0., 1.), smootherstep(-e, e, curvature) ), 1.); vec4 s = vec4(vec3(grid(v.xz * 5.0, .5, 2.)), .3); // Apply blending in place s.rgb = screenBlendMode(b.rgb, s.rgb); // Apply compositing vec4 color = destinationOverCompositing(s, b); // Use normals instead of curvature. //color = vec4(normal(v.xyz), 1.); gl_FragColor = color; }`, vert:` precision highp float; uniform mat4 projection, view, model; uniform vec2 step; attribute vec2 uv; varying vec4 v; varying float curvature; ${surface.r} float gaussianCurvature (vec2 uv) { /* Analytic evaluation of the gaussian curvature at point r(u,v) = (u, u**2 - v**2, v). References: Gray, A. in "Modern Differential Geometry of Curves and Surfaces with Mathematica". Please note that the dirivatives below depends on the surface r(u,v) and they must be calculated by hand. */ // To calculate by hand. // In comments, vec3 r is function of vec2 uv. vec3 ru = vec3(1., 2. * uv.x, 0.); // ru := ∂r/∂u vec3 rv = vec3(0., -2. * uv.y, 1.); // rv := ∂r/∂v vec3 ruu = vec3(0., 2., 0.); // ruu := ∂ru/∂u vec3 rvv = vec3(0., -2., 0.); // rvv := ∂rv/∂v vec3 ruv = vec3(0.); // ruv := ∂rv/∂u // End calculation by hand vec3 n = normalize(cross(ru, rv)); // normal vector to the surface r float E = dot(ru, ru); float F = dot(ru, rv); float G = dot(rv, rv); float e = dot(ruu, n); float f = dot(rvv, n); float g = dot(ruv, n); return (e * g - f * f) / (E * G - F * F) ; } float angle (vec3 p, vec3 q) { // Angle between vectors p and q. return acos(dot(normalize(p), normalize(q))); } float area (vec3 p, vec3 q) { return length(cross(p, q)) / 2.; } float discreteGaussianCurvature (vec2 uv) { /* Evaluation of the discrete gaussian curvature at point r(u, v). References: 1) S. L. Devadoss, J. O'Rourke in "Discrete and Computational Geometry" 2) Meyer et al. 2003 in "Discrete Differential-Geometry Operators for Triangulated 2-Manifolds" The code below depends on function "vec3 r (vec2 uv)" and variable "vec2 step" defined outside function discreteGaussianCurvature. Function "r" describes the 2D surface immersed in the 3D space, parametrized with coordinates uv. Variable "step" describes the spacing in coordinates uv in the 2D mesh, if for example u domain [a, b] is divided in n subintervals then step.x = (b - a) / n. */ // Vertices positions with respect vertex p0. // Labels are assigned clockwise around vertex p0. vec3 p0 = r(uv); vec3 v1 = r(uv + step * vec2(0., -1.)) - p0; vec3 v2 = r(uv + step * vec2(1., -1.)) - p0; vec3 v3 = r(uv + step * vec2(1., 0.)) - p0; vec3 v4 = r(uv + step * vec2(1., 1.)) - p0; vec3 v5 = r(uv + step * vec2(0., 1.)) - p0; vec3 v6 = r(uv + step * vec2(-1., 1.)) - p0; vec3 v7 = r(uv + step * vec2(-1., 0.)) - p0; vec3 v8 = r(uv + step * vec2(-1., -1.)) - p0; // Variable aij is the angle between vectors vi and vj. float a12 = angle(v1, v2); float a23 = angle(v2, v3); float a34 = angle(v3, v4); float a45 = angle(v4, v5); float a56 = angle(v5, v6); float a67 = angle(v6, v7); float a78 = angle(v7, v8); float a81 = angle(v8, v1); // Variable Aij is the area of the triangle of vertices (0, i, j). float A12 = area(v1, v2); float A23 = area(v2, v3); float A34 = area(v3, v4); float A45 = area(v4, v5); float A56 = area(v5, v6); float A67 = area(v6, v7); float A78 = area(v7, v8); float A81 = area(v8, v1); float twoPI = 6.283185307179586; float angleDeficit = twoPI - (a12 + a23 + a34 + a45 + a56 + a67 + a78 + a81); // Area of the region contained within the 1-ring neighborhood of each vertex. // Note that this is an approximation of what suggested by Meyer et al. float surfaceArea = A12 + A23 + A34 + A45 + A56 + A67 + A78 + A81; return angleDeficit / surfaceArea; } void main() { //curvature = gaussianCurvature(uv); // This is valid only for surface (u, u**2 - v**2, v). curvature = discreteGaussianCurvature(uv); v = vec4(r(uv), 1.); gl_Position = projection * view * model * v; }`, uniforms: { step: mesh.step, }, attributes: { uv: mesh.uv, }, elements: mesh.elements, count: mesh.elements.length * 3, primitive: 'triangles', })35:T47d,function createMesh (opts = {n: 100, uDomain: [0, 1], vDomain: [0, 1]}) { const {n, uDomain, vDomain} = opts; assert(n > 0, 'Variable n is not greater than zero.'); assert(uDomain[1] > uDomain[0], 'Interval uDomain is not valid.'); assert(vDomain[1] > vDomain[0], 'Interval vDomain is not valid.'); const u = i => uDomain[0] + i / n * (uDomain[1] - uDomain[0]); const v = j => vDomain[0] + j / n * (vDomain[1] - vDomain[0]); const elements = []; for (let i = 0; i < n; i++){ for (let j = 0; j < n; j++) { const vertexIndex = i + j * (n + 1); // Add top triangle (clockwise) elements.push([ vertexIndex, vertexIndex + 1, vertexIndex + n + 2 ]); // Add bottom triangle (clockwise) elements.push([ vertexIndex, vertexIndex + n + 2, vertexIndex + n + 1 ]); } } const uv = []; for (let j = 0; j <= n; j++){ for (let i = 0; i <= n; i++) { uv.push([u(i),v(j)]); } } const step = [ (uDomain[1] - uDomain[0]) / n, (vDomain[1] - vDomain[0]) / n, ] return {elements, uv, step}; }36:[] 2f:{"title":"Gaussian curvature in triangulated 2-Manifolds.","cells":[{"id":0,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":1731,"value":"viewof s = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":1844,"value":"md`Use the slider below to visualize narrow/large curvature intervals with red to blue colours. Red and blue colours describe negative and positive curvature values respectively. Intermediate colours, such as violet, describes curvature near zero. The coordinate lines represent values at constant ${tex`u`} and ${tex`v`}.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1816,"value":"viewof e = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":345,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1071,"value":"download = {\n const button = html`