Skip to content

PhaseNoise.kt

src/jvmMain/kotlin/afield/lang/stacksl/PhaseNoise.kt · 182 lines

package afield.lang.stacksl

StackSL's phase noise family: qnoise, phnoise, phrot — noise whose phase is a control, per frequency.

What it is for

Two fields made from the same random kernels, differing only in phase, are related in a precise way that no amount of offsetting or reseeding gives you. Rotate every kernel's phase by θ turns and the new field correlates with the old one by cos(2πθ): 0 is the same field, 0.25 an uncorrelated partner with the same spectrum, 0.5 its negative. Rotate only the kernels above some frequency and the two fields share their coarse structure exactly while their fine detail decorrelates — stereo pairs, RGB channels or particle species that belong together, and noise that churns at small scales while its large forms hold still.

Value noise (noise, fbm, fbm3) cannot do this. It has no phase to rotate: shifting its output is a translation or an offset, not a phase change. These words are built from kernels that have one — localised sinusoids, which is sparse Gabor noise.

The construction

Each octave is a jittered lattice with one kernel per cell: a Gaussian window exp(-2 r²), truncated to zero at r = 1.5 cells, times a sinusoid with a random direction, a random phase, and a carrier of PhaseNoiseConstants.CARRIER cycles per cell, jittered by up to ±PhaseNoiseConstants.CARRIER_JITTER_OCTAVES/2 octaves per kernel so the spectrum is continuous across a band rather than a comb. The same kernels are summed twice, once through cos and once through sin — a quadrature pair — so the result is (Nc, Ns) and every phase rotation of it is Nc·cos θ − Ns·sin θ, exactly. Each octave is divided by the square root of its own window energy, which is independent of phase, so the identity survives normalisation.

Octave o is the lattice scaled by lacunarity^o and rotated once per octave by the fixed (PhaseNoiseConstants.OCTAVE_ROTATION_COS, PhaseNoiseConstants.OCTAVE_ROTATION_SIN) pair, so a kernel's frequency in cycles per unit of the input coordinate is its carrier times lacunarity^o. That number is what a phase response sees.

Output and range

vec4(Nc, Ns, |N|, 1): the quadrature pair, its envelope (phase-invariant, so still right after a phrot), and an opaque alpha. Nc and Ns are signed and zero-mean, with an RMS near PhaseNoiseConstants.OUTPUT_SCALE·0.7 ≈ 0.35, so they sit mostly inside -1..1; PhaseNoiseTest pins those statistics.

Cost

Sixteen kernels per octave, each a sin/cos pair, two hashes and an exp: five octaves is eighty kernels. That is more than fbm3 and less than curl3. The family is 2D over input.xy; input.w reaches a phase response as time.

Phase is in turns throughout, like phacelle's offset and unlike rotate's radians: 0.25 is a quarter cycle.

object PhaseNoiseConstants

The fixed numbers of the phase noise construction, read by the emitted GLSL.

object PhaseNoiseConstants {

const val MAX_OCTAVES = VolumetricNoiseConstants.MAX_OCTAVES

Loop bound in the emitted GLSL, shared with the volumetric family.

const val CARRIER = 1.0f

Cycles of the carrier per lattice cell, before jitter. With a window of standard deviation half a cell the band is about ±0.32 cycles wide around this, so at lacunarity 2 consecutive octaves' bands meet rather than leaving gaps.

const val CARRIER_JITTER_OCTAVES = 0.5f

Per-kernel carrier jitter, as a total width in octaves (±half of it).

const val WINDOW_FALLOFF = 2.0f

The window exp(-WINDOW_FALLOFF r²).

const val WINDOW_FLOOR = 0.011108996f

exp(-WINDOW_FALLOFF * 1.5²), subtracted from the window so it reaches zero at 1.5 cells. A 4x4 neighbourhood with jitter of ±half a cell is complete for that radius, so no kernel that should contribute is ever skipped and the field has no seams.

const val ENERGY_EPSILON = 1.0e-6f

Added under the energy square root so an octave never divides by zero.

const val OUTPUT_SCALE = 0.5f

Final scale on the quadrature pair; see the class notes for the resulting RMS.

const val OCTAVE_ROTATION_COS = 0.8f

Per-octave lattice rotation (cos, sin) — a 3-4-5 triangle, so it is exactly orthonormal.

const val OCTAVE_ROTATION_SIN = 0.6f

const val OCTAVE_OFFSET_X = 17.17f

Translation after each octave's rotate-and-scale, so the lattices do not share an origin.

const val OCTAVE_OFFSET_Y = 9.53f

const val DEFAULT_OCTAVES = VolumetricNoiseConstants.DEFAULT_OCTAVES

const val DEFAULT_LACUNARITY = VolumetricNoiseConstants.DEFAULT_LACUNARITY

const val DEFAULT_GAIN = VolumetricNoiseConstants.DEFAULT_GAIN

internal fun phaseNoiseHelperSource(prefix: String, helpers: Set<String>): String

Helpers the phase family emits through the shared helper set rather than per call site.

internal fun phaseNoiseHelperSource(prefix: String, helpers: Set<String>): String {
    if ("phrot" !in helpers) return ""
    return buildString {
        appendLine("vec4 ${prefix}_phrot(vec4 q, float turns) {")
        appendLine("    float a = 6.283185307179586 * turns;")
        appendLine("    float c = cos(a);")
        appendLine("    float s = sin(a);")
        appendLine("    return vec4(q.x * c - q.y * s, q.x * s + q.y * c, q.z, q.w);")
        appendLine("}")
        appendLine()
    }
}

internal fun quadratureNoiseFunctionSource

One quadrature noise function, named helperName. shiftExpression is a GLSL float in turns added to every kernel's phase; it may read kernel, a vec4(frequency, octave, random, time) describing the kernel being summed, and any parameter in extraParameters. It is evaluated once per kernel, which is what lets the shift depend on frequency.

Calls ${prefix}_hash33, so the caller must also request the volumetric hash33 helper.

body — 54 lines, from line 120
internal fun quadratureNoiseFunctionSource(
    prefix: String,
    helperName: String,
    shiftExpression: String,
    extraParameters: String,
    wireParamPlaceholder: String
): String {
    val k = PhaseNoiseConstants
    fun g(value: Float) = glslFloatLiteral(value)
    val maxOct = k.MAX_OCTAVES
    return """
        vec4 $helperName(vec4 pValue, float octaves, float lacunarity, float gain$extraParameters$wireParamPlaceholder) {
            const float tau = 6.283185307179586;
            int n = int(clamp(octaves, 1.0, ${g(maxOct.toFloat())}));
            mat2 rotation = mat2(${g(k.OCTAVE_ROTATION_COS)}, ${g(k.OCTAVE_ROTATION_SIN)}, ${g(-k.OCTAVE_ROTATION_SIN)}, ${g(k.OCTAVE_ROTATION_COS)});
            vec2 q = pValue.xy;
            float scale = 1.0;
            float amp = 1.0;
            float norm = 0.0;
            vec2 quad = vec2(0.0);
            for (int o = 0; o < $maxOct; o++) {
                if (o >= n) break;
                vec2 cellBase = floor(q);
                vec2 f = q - cellBase;
                vec2 octaveSum = vec2(0.0);
                float energy = 0.0;
                for (int i = -1; i <= 2; i++) {
                    for (int j = -1; j <= 2; j++) {
                        vec2 cell = cellBase + vec2(float(i), float(j));
                        vec3 h0 = ${prefix}_hash33(vec3(cell, float(o) * 2.0));
                        vec3 h1 = ${prefix}_hash33(vec3(cell, float(o) * 2.0 + 1.0));
                        vec2 r = f - vec2(float(i), float(j)) - (h0.xy - 0.5);
                        float w = max(exp(-${g(k.WINDOW_FALLOFF)} * dot(r, r)) - ${g(k.WINDOW_FLOOR)}, 0.0);
                        float angle = tau * h0.z;
                        vec2 dir = vec2(cos(angle), sin(angle));
                        float cycles = ${g(k.CARRIER)} * exp2((h1.x - 0.5) * ${g(k.CARRIER_JITTER_OCTAVES)});
                        vec4 kernel = vec4(cycles * scale, float(o), h1.z, pValue.w);
                        float shift = $shiftExpression;
                        float phase = tau * (cycles * dot(r, dir) + h1.y + shift);
                        octaveSum += w * vec2(cos(phase), sin(phase));
                        energy += w * w;
                    }
                }
                quad += amp * octaveSum * inversesqrt(energy + ${g(k.ENERGY_EPSILON)});
                norm += amp * amp;
                q = rotation * q * lacunarity + vec2(${g(k.OCTAVE_OFFSET_X)}, ${g(k.OCTAVE_OFFSET_Y)});
                scale *= lacunarity;
                amp *= gain;
            }
            quad *= ${g(k.OUTPUT_SCALE)} * inversesqrt(max(norm, ${g(k.ENERGY_EPSILON)}));
            return vec4(quad, length(quad), 1.0);
        }
    """.trimIndent()
}

internal const val PER_OCTAVE_SHIFT_EXPRESSION =

phnoise's shift: phases holds the turns for octaves 0, 1 and 2, and for octave 3 and everything finer. The weights are a one-hot of the octave index built from abs/clamp, not a variable vector index, so the CPU backend's GLSL subset runs it.

    "dot(phases, clamp(1.0 - abs(vec4(0.0, 1.0, 2.0, 3.0) - min(kernel.y, 3.0)), 0.0, 1.0))"