← geometric-flows

Stretch and stack

Stretch-and-stack flow on lattices: the unit cell warps under the flow and restacks, driven by the shared lattice computation.

Source

// Stretch and Stack — driven by the shared JS lattice driver (script.js).
// uLattice (the lattice mat2) and uFlowInv (inverse of the per-frame flow) are
// pushed from script.js; uScale comes from the UI. UnitCell warps the previous
// frame by uFlowInv to animate the stretch, like the original Shadertoy.

#include "../shared/2d_utils.glsl"
#include "../shared/coord_utils.glsl"

bool fadeMoves   = false;   // fade the warped trail instead of hard-clearing it

vec4 backCol = vec4(1., 1., 1., 0.);
vec4 dotCol  = vec4(.8, .05, .04, .1);
vec4 edgeCol = vec4(.651, .669, .95, .1);

// Screen origin for this view: horizontally centred, one third up from the bottom.
vec2 viewOffset(in vec3 res) { return vec2(res.x * 0.5, res.y / 3.); }

// Buffer for whole plane

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec4 col;
    
    vec2 p = coords(fragCoord, iResolution, uScale, viewOffset(iResolution));

    mat2 m = mat2(uLattice);
    
    float a = vec4(m).x;
    float b = -vec4(m).z;
    float c = vec4(m).y;
    float d = vec4(m).w;

    mat2 im = inverse(m);
    vec2 pC = p-m*floor(im*p);

    if(a > b)
    {
        if(pC.x < 0. && pC.y > d) { pC = pC - vec2(-b,d); }
        else if(pC.x > 0.) 
        { 
            pC = pC + floor(pC.x/b)*vec2(-b, d);
            if(pC.y > c ) 
            {
                pC = pC - floor((pC.y-c)/d)*vec2(-b, d);
                pC = pC - vec2(a,c);
            }
        }
    }
    else
    {
        if(pC.x > 0. && pC.y > c) { pC = pC - vec2(a,c); }
        else if(pC.x < 0.) 
        { 
            pC = pC + floor(-pC.x/a)*vec2(a, c);
            if(pC.y > d ) 
            {
                pC = pC - floor((pC.y-d)/c)*vec2(a, c);
                pC = pC - vec2(-b,d);
            }
        }
    }
        
	vec2 uv = iCoords(pC, iResolution, uScale, viewOffset(iResolution)) / iResolution.xy;
    
    if(uv.y > .998) { uv.y = .998;}
    
    col = vec4(texture(UnitCell,uv).xyz,1);

    // Output to screen
    fragColor = col;
}

// Fork of "Modular Surface Flows" by Gelada. https://shadertoy.com/view/wsdXRN
// 2019-10-28 18:45:24

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec4 col;

	vec2 uv = fragCoord.xy / iResolution.xy;
    
    if(uv.y > .998) { uv.y = .998;}
    
    if(uView == 1) { col = vec4(texture(fullLattice,uv).xyz,1); }
    else { col = vec4(texture(UnitCell,uv).xyz,1); }    

    // Output to screen
    fragColor = col;
}

// Unit cell — warps the previous frame back by one flow step (uFlowInv) to
// animate the stretch, then draws the current fundamental domain on top.
// Restored from the original Shadertoy, with uFlowInv replacing gFlow(-1.) and
// uLattice replacing the bufA matrix.

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 p = coords(fragCoord, iResolution, uScale, viewOffset(iResolution));

    vec2 oFrag = iCoords(mat2(uFlowInv) * p, iResolution, uScale, viewOffset(iResolution));
    vec2 uv = oFrag / iResolution.xy;

    vec4 col = backCol;
    vec4 pt = uLattice;       // (b1.x, b1.y, b2.x, b2.y)
    float mixing = 0.;
    float ad = 0.0;



    if (iFrame == 0 || uNewLattice || uv.x > 1. || uv.x < 0. || uv.y > .99 || uv.y < 0.)
    {
        col = drBox(col, p, vec4(-pt.x-ad, 0., 0., pt.w+ad), vec4(edgeCol.xyz, 1.), .0);
        col = drBox(col, p, vec4(0., 0., -pt.z+ad, pt.y+ad), vec4(dotCol.xyz, 1.), .0);
    }
    else
    {
        col = vec4(texture(UnitCell, uv).xyz, 1.);
        vec4 fullCol = vec4(texture(fullLattice, uv).xyz, 1.);

        if (fadeMoves) {
            col = mix(col, vec4(backCol.xyz, 0), .1);
            mixing = .1;
        } else {
        	col = vec4(backCol.xyz,0);
            mixing = 1.;
        }

        if (uFade == 1) {
            mixing = 1.0;
            col = drBox(col, p, vec4(-pt.x-ad, 0., 0., pt.w+ad), vec4(fullCol.xyz, mixing), 0.01);
            col = drBox(col, p, vec4(0., 0., -pt.z+ad, pt.y+ad), vec4(fullCol.xyz, mixing), 0.01);
            mixing = .01;
            col = drBox(col, p, vec4(-pt.x-ad, 0., 0., pt.w+ad), vec4(edgeCol.xyz, mixing), .001);
            col = drBox(col, p, vec4(0., 0., -pt.z+ad, pt.y+ad), vec4(dotCol.xyz, mixing), .001);
        } else {
            col = drBox(col, p, vec4(-pt.x-ad, 0., 0., pt.w+ad), vec4(edgeCol.xyz, mixing), .1);
            col = drBox(col, p, vec4(0., 0., -pt.z+ad, pt.y+ad), vec4(dotCol.xyz, mixing), .1);
        }
    }

    fragColor = col;
}

//import * as math from 'mathjs';
import { matrix, dot, add, subtract, multiply, norm, hypot, subset, index, inv} from 'mathjs';
import * as cfUtils from '../shared/cf_utils.js';

let b1 = [1, 1.618033988749895];
let b2 = [-1.618033988749895, 1];

let p1 = null, p2 = null, hasReturnedFar = false;

let latticeSet = false;

function cDiv(a, b) {
  const d = b[0]*b[0] + b[1]*b[1];
  return [(a[0]*b[0] + a[1]*b[1]) / d, (a[1]*b[0] - a[0]*b[1]) / d];
}

function hpPt(b1, b2) {
  let pt = hypot(...b1) > hypot(...b2) ? cDiv(b1, b2) : cDiv(b2, b1);
  if (pt[1] < 0) pt = [-pt[0], -pt[1]];
  return pt;
}

function checkCF(flow) {
  console.log('checkPeriodSnap called, flow:', flow, 'p1:', p1);
  if (p1 === null) { return; }
  if (flow !== 0) { return; }

  

  const dist1 = norm(subtract(b1,p1));
  const dist2 = norm(subtract(b2,p2));

  const dist = Math.sqrt(dist1*dist1+dist2*dist2);
  

 if (dist < .01) {
    b1 = [...p1];
    b2 = [...p2];
    hasReturnedFar = false;
  }
}

// Input: basis as [[a, b], [c, d]] (rows are basis vectors)
// Output: reduced basis with b1 = shortest vector
function lagrangeGauss(b1, b2) {
  b1 = matrix(b1);
  b2 = matrix(b2);
  const normSq = v => dot(v, v);

  if (normSq(b1) > normSq(b2)) { [b1, b2] = [b2, b1]; }

  while (true) {
    const m = Math.round(dot(b1, b2) / normSq(b1));
    b2 = subtract(b2, multiply(m, b1));
    if (normSq(b2) >= normSq(b1)) break;
    [b1, b2] = [b2, b1];
  }

  if(subset(b1,index(1))<0) { b1 = multiply(-1,b1);}
  if(subset(b2,index(1))<0) { b2 = multiply(-1,b2);}

  return [b1.toArray(), b2.toArray()];
}

// Constrained Gauss–Lagrange reduction.
//   Keeps b1 in upper right quadrant
//   Keeps b2 in upper left quadrant
function coneReduce(b1, b2) {
  // Accept plain arrays or mathjs matrices (onFrame passes multiply results).
  b1 = (b1.toArray ? b1.toArray() : b1).slice();
  b2 = (b2.toArray ? b2.toArray() : b2).slice();
  const clamp = (k, lo, hi) => Math.max(lo, Math.min(hi, k));

  let changed = true;
  while (changed) {
    changed = false;

    // shorten b2 along b1, keeping b2 in QII (x <= 0, y >= 0)
    {
      const q = dot(b1, b1);
      if (q > 0) {
        let k = Math.round(-dot(b1, b2) / q);          // unconstrained minimiser
        const hi = b1[0] > 0 ? Math.floor(-b2[0] / b1[0]) :  Infinity; // x <= 0
        const lo = b1[1] > 0 ? Math.ceil (-b2[1] / b1[1]) : -Infinity; // y >= 0
        k = clamp(k, lo, hi);
        if (k !== 0) { b2 = add(b2, multiply(k, b1)); changed = true; }
      }
    }

    // shorten b1 along b2, keeping b1 in QI (x >= 0, y >= 0)
    {
      const q = dot(b2, b2);
      if (q > 0) {
        let k = Math.round(-dot(b1, b2) / q);
        const hi = b2[0] < 0 ? Math.floor(-b1[0] / b2[0]) :  Infinity; // x >= 0
        const lo = b2[1] > 0 ? Math.ceil (-b1[1] / b2[1]) : -Infinity; // y >= 0
        k = clamp(k, lo, hi);
        if (k !== 0) { b1 = add(b1, multiply(k, b2)); changed = true; }
      }
    }
  }
  return [b1, b2]; // b1/b2 are plain arrays (add/multiply preserve arrays)
}

const fps = 60;

// Geodesic (hyperbolic diagonal) flow matrix — stretches one axis, shrinks the other
function gFlow(dTime) {
  return matrix([
    [Math.exp(dTime / (3 * fps)), 0],
    [0, Math.exp(-dTime / (3 * fps))]
  ]);
}

// Horocyclic (shear) flow matrix (vertical)
function hFlow(dTime) {
  return matrix([
    [1, 0],
    [dTime / fps, 1]
  ]);
}

// Horocyclic (shear) flow matrix (horizontal)
function hFlowH(dTime) {
  return matrix([
    [1, dTime / fps],
    [0, 1]
  ]);
}

export function onSetLattice(inputs, handle) {
  const alpha = cfUtils.quadCFStrings(inputs.cfPrefix ?? '', inputs.cfPeriod ?? '0');
  const beta = cfUtils.quadCFStrings('', cfUtils.reverseCFString(inputs.cfPeriod) ?? '0');
  b1 = [1, beta];
  b2 = [-alpha, 1];

  [b1,b2] = coneReduce(b1,b2);
  
  

  //handle.setUniform('uLattice', [1, beta, -alpha, 1]);
  latticeSet = true;
}

export function setup(engine) {
  b1 = [1, 1.618033988749895];
  b2 = [-1.618033988749895, 1];
}

export function onFrame(engine, time, deltaTime, frame) {
  const flow = engine.getUniformValue('uFlow');
  const speed = engine.getUniformValue('uSpeed');

  let m;
  switch (flow) {
    case 1:  m = hFlow(speed); break;
    case 2:  m = hFlowH(speed); break;
    default: m = gFlow(speed * 2.);
  }

  const mb1 = multiply(m, b1);
  const mb2 = multiply(m, b2);

  [b1, b2] = coneReduce(mb1, mb2);

  checkCF(flow);

  // Live CF of the flow's forward endpoint — the future itinerary. This endpoint
  // ratio is invariant under the flow (the per-frame scale factors cancel), so it
  // stays stable and only advances one CF digit per coneReduce reduction (the
  // Gauss map). The invariant ratio depends on flow type/direction:
  //   geodesic (0): x-ratio forward when speed>=0, else y-ratio (reversed flow)
  //   horocyclic vertical (1):   x-ratio (parabolic fixed endpoint)
  //   horocyclic horizontal (2): y-ratio
  let endNum, endDen;
  if (flow === 2 || (flow === 0 && speed < 0)) {
    endNum = b2[1]; endDen = b1[1];   // y-ratio
  } else {
    endNum = b2[0]; endDen = b1[0];   // x-ratio
  }
  if (Math.abs(endDen) > 1e-9) {
    engine.emit('cf', cfUtils.cfExpansion(Math.abs(endNum / endDen), 12).join(', '));
  }

  // Inverse of this frame's smooth flow (pre-reduction), flattened column-major
  // so GLSL `mat2(uFlowInv) * p` walks a point back one flow step. Used by the
  // stretch-warp trails; shaders that don't read it simply ignore it.
  const mi = inv(m);
  engine.setUniformValue('uFlowInv', [mi.get([0, 0]), mi.get([1, 0]), mi.get([0, 1]), mi.get([1, 1])]);

  engine.setUniformValue('uLattice', [b1[0], b1[1], b2[0], b2[1]]);
  engine.setUniformValue('uNewLattice', latticeSet);
  latticeSet = false;
}