diff --git a/src/core/shaders.ts b/src/core/shaders.ts index 09ef681..127a103 100644 --- a/src/core/shaders.ts +++ b/src/core/shaders.ts @@ -168,11 +168,17 @@ export const displayShader = /* glsl */ ` uniform int uAlgorithm; uniform int uEnableAlpha; + float getDensity(vec2 uv) { + float raw = max(texture2D(uTexture, uv).r, 0.0); + float e2 = exp(2.0 * min(raw, 10.0)); + return (e2 - 1.0) / (e2 + 1.0); + } + void main () { float obs = texture2D(uObstacle, vUv).r; // Smooth fade: density falls off over the blurred obstacle boundary zone // rather than cutting off at a hard step, eliminating the jagged fringe. - float density = max(texture2D(uTexture, vUv).r, 0.0) * (1.0 - obs); + float density = getDensity(vUv) * (1.0 - obs); float coverage = texture2D(uCoverage, vUv).r; // 8-tap Sobel normal — computes gradient via 3×3 kernel (no centre sample needed). @@ -180,14 +186,14 @@ export const displayShader = /* glsl */ ` // ≈ 3 sim texels per axis. The Sobel kernel properly averages the gradient in // every direction — no 4-tap cross bias that creates 45° circuit-board artefacts. float sx = texelSize.x * 6.0, sy = texelSize.y * 6.0; - float d00 = max(texture2D(uTexture, vUv + vec2(-sx, -sy)).r, 0.0); - float d10 = max(texture2D(uTexture, vUv + vec2(0.0, -sy)).r, 0.0); - float d20 = max(texture2D(uTexture, vUv + vec2( sx, -sy)).r, 0.0); - float d01 = max(texture2D(uTexture, vUv + vec2(-sx, 0.0)).r, 0.0); - float d21 = max(texture2D(uTexture, vUv + vec2( sx, 0.0)).r, 0.0); - float d02 = max(texture2D(uTexture, vUv + vec2(-sx, sy)).r, 0.0); - float d12 = max(texture2D(uTexture, vUv + vec2(0.0, sy)).r, 0.0); - float d22 = max(texture2D(uTexture, vUv + vec2( sx, sy)).r, 0.0); + float d00 = getDensity(vUv + vec2(-sx, -sy)); + float d10 = getDensity(vUv + vec2(0.0, -sy)); + float d20 = getDensity(vUv + vec2( sx, -sy)); + float d01 = getDensity(vUv + vec2(-sx, 0.0)); + float d21 = getDensity(vUv + vec2( sx, 0.0)); + float d02 = getDensity(vUv + vec2(-sx, sy)); + float d12 = getDensity(vUv + vec2(0.0, sy)); + float d22 = getDensity(vUv + vec2( sx, sy)); float gx = (d20 + 2.0*d21 + d22) - (d00 + 2.0*d01 + d02); float gy = (d02 + 2.0*d12 + d22) - (d00 + 2.0*d10 + d20); diff --git a/src/core/wgsl-shaders.ts b/src/core/wgsl-shaders.ts index d36a27b..3c71d7f 100644 --- a/src/core/wgsl-shaders.ts +++ b/src/core/wgsl-shaders.ts @@ -12,7 +12,7 @@ // render-to-texture round-trips consistent (no Y flip in simulation FBOs). // Source textures are uploaded WITHOUT flipY so their row 0 (visual top) // also lands at uv.y=0. -const SHARED_VS_STRUCT = /* wgsl */` +const SHARED_VS_STRUCT = /* wgsl */ ` struct VSOut { @builtin(position) pos : vec4f, @location(0) uv : vec2f, @@ -24,7 +24,7 @@ struct VSOut { // ─── Advection ─────────────────────────────────────────────────────────────── -export const advectionWGSL = /* wgsl */` +export const advectionWGSL = /* wgsl */ ` ${SHARED_VS_STRUCT} struct U { @@ -60,7 +60,7 @@ struct U { // ─── Divergence ────────────────────────────────────────────────────────────── -export const divergenceWGSL = /* wgsl */` +export const divergenceWGSL = /* wgsl */ ` ${SHARED_VS_STRUCT} struct U { texelSize: vec2f, _pad: vec2f } @@ -91,7 +91,7 @@ struct U { texelSize: vec2f, _pad: vec2f } // ─── Pressure ──────────────────────────────────────────────────────────────── -export const pressureWGSL = /* wgsl */` +export const pressureWGSL = /* wgsl */ ` ${SHARED_VS_STRUCT} struct U { texelSize: vec2f, _pad: vec2f } @@ -125,7 +125,7 @@ struct U { texelSize: vec2f, _pad: vec2f } // ─── Gradient subtract ─────────────────────────────────────────────────────── -export const gradientSubtractWGSL = /* wgsl */` +export const gradientSubtractWGSL = /* wgsl */ ` ${SHARED_VS_STRUCT} struct U { texelSize: vec2f, _pad: vec2f } @@ -160,7 +160,7 @@ struct U { texelSize: vec2f, _pad: vec2f } // ─── Splat ─────────────────────────────────────────────────────────────────── -export const splatWGSL = /* wgsl */` +export const splatWGSL = /* wgsl */ ` ${SHARED_VS_STRUCT} // texelSize occupies bytes 0-7 for the shared vertex stage; aspectRatio/radius fill the rest. @@ -197,7 +197,7 @@ struct U { // ─── Curl ──────────────────────────────────────────────────────────────────── -export const curlWGSL = /* wgsl */` +export const curlWGSL = /* wgsl */ ` ${SHARED_VS_STRUCT} struct U { texelSize: vec2f, _pad: vec2f } @@ -227,7 +227,7 @@ struct U { texelSize: vec2f, _pad: vec2f } // ─── Vorticity ─────────────────────────────────────────────────────────────── -export const vorticityWGSL = /* wgsl */` +export const vorticityWGSL = /* wgsl */ ` ${SHARED_VS_STRUCT} struct U { @@ -278,7 +278,7 @@ struct U { // 56 algorithm i32 // 60 enableAlpha i32 (1 = premultiplied alpha output, 0 = opaque) -export const displayWGSL = /* wgsl */` +export const displayWGSL = /* wgsl */ ` ${SHARED_VS_STRUCT} struct U { @@ -300,6 +300,11 @@ struct U { @group(0) @binding(5) var uCov : texture_2d; @group(0) @binding(6) var uVel : texture_2d; +fn getDensity(uv: vec2f) -> f32 { + let raw = max(textureSample(uTex, samp, uv).r, 0.0); + return tanh(raw); +} + @vertex fn vs(@location(0) a: vec2f) -> VSOut { var o: VSOut; o.uv = vec2f(a.x * 0.5 + 0.5, 0.5 - a.y * 0.5); @@ -313,19 +318,19 @@ struct U { @fragment fn fs(i: VSOut) -> @location(0) vec4f { let obs = textureSample(uObs, samp, i.uv).r; - let density = max(textureSample(uTex, samp, i.uv).r, 0.0) * (1.0 - obs); + let density = getDensity(i.uv) * (1.0 - obs); let cov = textureSample(uCov, samp, i.uv).r; let sx = u.texelSize.x * 6.0; let sy = u.texelSize.y * 6.0; - let d00 = max(textureSample(uTex, samp, i.uv + vec2f(-sx, -sy)).r, 0.0); - let d10 = max(textureSample(uTex, samp, i.uv + vec2f(0.0, -sy)).r, 0.0); - let d20 = max(textureSample(uTex, samp, i.uv + vec2f( sx, -sy)).r, 0.0); - let d01 = max(textureSample(uTex, samp, i.uv + vec2f(-sx, 0.0)).r, 0.0); - let d21 = max(textureSample(uTex, samp, i.uv + vec2f( sx, 0.0)).r, 0.0); - let d02 = max(textureSample(uTex, samp, i.uv + vec2f(-sx, sy)).r, 0.0); - let d12 = max(textureSample(uTex, samp, i.uv + vec2f(0.0, sy)).r, 0.0); - let d22 = max(textureSample(uTex, samp, i.uv + vec2f( sx, sy)).r, 0.0); + let d00 = getDensity(i.uv + vec2f(-sx, -sy)); + let d10 = getDensity(i.uv + vec2f(0.0, -sy)); + let d20 = getDensity(i.uv + vec2f( sx, -sy)); + let d01 = getDensity(i.uv + vec2f(-sx, 0.0)); + let d21 = getDensity(i.uv + vec2f( sx, 0.0)); + let d02 = getDensity(i.uv + vec2f(-sx, sy)); + let d12 = getDensity(i.uv + vec2f(0.0, sy)); + let d22 = getDensity(i.uv + vec2f( sx, sy)); let gx = (d20 + 2.0*d21 + d22) - (d00 + 2.0*d01 + d02); let gy = (d02 + 2.0*d12 + d22) - (d00 + 2.0*d10 + d20); diff --git a/tests/core/shaders.test.js b/tests/core/shaders.test.js new file mode 100644 index 0000000..089df6c --- /dev/null +++ b/tests/core/shaders.test.js @@ -0,0 +1,65 @@ +import { describe, expect, it } from 'vitest'; + +import { displayShader } from '../../src/core/shaders.ts'; +import { displayWGSL } from '../../src/core/wgsl-shaders.ts'; + +describe('displayShader density saturation', () => { + it('GLSL displayShader contains getDensity helper', () => { + expect(displayShader).toContain('getDensity'); + }); + + it('GLSL getDensity uses tanh-equivalent soft saturation', () => { + // The implementation uses (e2-1)/(e2+1) which equals tanh(raw). + // Verify the shader does NOT use the old unbounded max(..., 0.0) directly + // for the main density or Sobel samples. + const lines = displayShader.split('\n'); + const mainDensityLine = lines.find((l) => l.includes('float density') && l.includes('1.0 - obs')); + expect(mainDensityLine).toContain('getDensity'); + }); + + it('WGSL displayWGSL contains getDensity helper', () => { + expect(displayWGSL).toContain('getDensity'); + }); + + it('WGSL getDensity uses tanh', () => { + expect(displayWGSL).toContain('tanh('); + }); + + it('WGSL main density line uses getDensity', () => { + const lines = displayWGSL.split('\n'); + const mainDensityLine = lines.find((l) => l.includes('let density') && l.includes('1.0 - obs')); + expect(mainDensityLine).toContain('getDensity'); + }); +}); + +describe('getDensity saturation behaviour (JS simulation)', () => { + // Verify the tanh-like formula used in GLSL: (e^2x - 1) / (e^2x + 1) + const glslGetDensity = (raw) => { + const r = Math.max(raw, 0.0); + const e2 = Math.exp(2.0 * Math.min(r, 10.0)); + return (e2 - 1.0) / (e2 + 1.0); + }; + + it('maps 0 to 0', () => { + expect(glslGetDensity(0)).toBeCloseTo(0, 5); + }); + + it('maps 1 to < 1', () => { + expect(glslGetDensity(1)).toBeLessThan(1.0); + }); + + it('maps values > 1 to < 1 (prevents glow saturation)', () => { + expect(glslGetDensity(2.0)).toBeLessThan(1.0); + expect(glslGetDensity(5.0)).toBeLessThan(1.0); + expect(glslGetDensity(100.0)).toBeLessThan(1.0); + }); + + it('is monotonically increasing', () => { + expect(glslGetDensity(0.5)).toBeLessThan(glslGetDensity(1.0)); + expect(glslGetDensity(1.0)).toBeLessThan(glslGetDensity(2.0)); + }); + + it('clamps negative raw values to 0', () => { + expect(glslGetDensity(-1.0)).toBeCloseTo(0, 5); + }); +});