Source code of plot #040 back to plot

///<reference path="../pub/lib/paper.d.ts" />
import {info, init, loadLib, setSketch, dbgRedraw} from "./utils/boilerplate.js";
import {mulberry32, rand, rand_range, rand_select, randn_bm, setRandomGenerator, shuffle} from "./utils/random.js"
import * as THREE from "../pub/lib/three.module.js";
import {Vector3} from "../pub/lib/three.module.js";
import {OrbitControls} from "../pub/lib/OrbitControls.js"
import {WebGLCanvasMasker} from "./utils/webgl-canvas-masker.js";


// https://codepen.io/dissimulate/pen/nYQrNP
// https://steven.codes/blog/cloth-simulation/
// https://istemi-bahceci.github.io/blog/physically%20based%20simulations/2017/03/14/ClothSimulation.html
// https://github.com/matthias-research/pages/blob/master/tenMinutePhysics/14-cloth.html

const pw = 2100;    // Paper width
const ph = 1480;    // Paper height
const w = 1480;     // Drawing width
const h = 1050;     // Drawing height
const margin = 50;

const cpars = {
  // nClothSquares: 120,
  nClothSquares: 120,
  clothSize: 3,
  nSimIters: 8,
  gravity: 1,
  kStructural: 100,
  kShear: 1000,
  kFlex: 200,
  dampening: 0.995,
  dt: 0.016,
};

const nIters1 = 100;
const nIters2 = 100;


let camProps = {
  fov: 30,
  //position: { x: 0, y: -0.5, z: 3}, target: { x: 0, y: 0, z: 0 },
  //"position":{"x":0,"y":2.73,"z":5.29},"target":{"x":0,"y":0.16,"z":0}
  "position":{"x":-3.82,"y":1.64,"z":1.98},"target":{"x":-0.12,"y":-0.58,"z":-0.09}
};
let sceneBgColor = "black";

// Canvas masker & Three JS canvas/machinery
let segLen = 1;
let joinLen = 3;
let rf = 2;                 // Occlusion canvas & three canvas are this many times larger than our area
let elmThreeCanvas;
let renderer, scene, cam, ray;

// If true, we'll create colorizer and colorize faces to sample for vector graphics
// If false, we'll render wireframe-like 3D model
let colorizer = true;

// If true, orbit controls will be created for interactive exploration of 3D model
// If false, we'll render plot
let controls = false;

let seed;     // Random seed
if (window.fxhash) seed = Math.round(fxrand() * 65535);
else seed = Math.round(Math.random() * 65535);
// seed = 57737;

setRandomGenerator(mulberry32(seed));

setSketch(function () {

  info("Seed: " + seed, seed);
  init(w, h, pw, ph);

  // Three JS canvas
  initThree();

  setTimeout(draw, 10);
});

async function draw() {

  paper.project.currentStyle.strokeColor = "black";
  paper.project.currentStyle.strokeWidth = 2;

  let frame = Path.Rectangle(margin, margin, w - 2 * margin, h - 2 * margin);
  // project.activeLayer.addChild(frame);

  // If colorizer is requested, create it
  if (colorizer) colorizer = new Colorizer(31, 31, 31);

  const cloth = new Cloth(0.5, false);
  // Collision testers
  cloth.objFun = (x, y, z, prevx, prevy, prevz) => {
    // Object 1
    if (x >= -0.3 && x <= 0.2 && z >= -0.4 && z <= -0.2) {
      if (y < 0.45) return true;
    }
    // Object 2
    if (x >= 0.4 && x <= 0.6 && z >= 0.4 && z <= 0.5) {
      if (y < 0.35) return true;
    }
    // Object 3
    if (x >= -0.7 && x <= -0.6 && z >= 0.6 && z <= 0.7) {
      if (y < 0.25) return true;
    }
    // Floor-ish
    if (x >= -1.0 && x <= 1.0 && z >= -1.0 && z <= 1.0) {
      if (y < -0.3) return true;
    }
    return null;
  };

  for (let i = 0; i < nIters1; ++i) {
    cloth.update(cpars.dt);
    if ((i % 100) == 0) console.log(i);
  }
  const mesh = cloth.makeMesh(colorizer);
  scene.add(mesh);

  let iter = 0;
  cloth.updateGeo();
  requestAnimationFrame(animate);

  // If orbit controls are requested, create them, and kick off animation
  if (controls) {
    controls = new OrbitControls(cam, renderer.domElement);
    controls.target.set(camProps.target.x, camProps.target.y, camProps.target.z);
  }

  function renderPlot() {

    let pixels = new Uint8Array(w * rf * h * rf * 4);
    let ctx = elmThreeCanvas.getContext("webgl2", {preserveDrawingBuffer: true});
    ctx.readPixels(0, 0, ctx.drawingBufferWidth, ctx.drawingBufferHeight, ctx.RGBA, ctx.UNSIGNED_BYTE, pixels);

    const lines = cloth.getLines();
    const wcm = new WebGLCanvasMasker(pixels, w, h, rf, true);

    const allLines = lines.horiz.concat(...lines.deep).concat(...lines.diag);
    const visibleLines = wcm.mask(allLines, null, segLen);

    const joinedPaths = [];
    let currPts = [];
    for (const vl of visibleLines) {
      if (currPts.length == 0) {
        currPts.push(vl[0], vl[1]);
      }
      else if (currPts[currPts.length-1].getDistance(vl[0]) < joinLen) {
        currPts.push(vl[1]);
      }
      else {
        joinedPaths.push(currPts);
        currPts = vl;
      }
    }
    if (currPts.length != 0) joinedPaths.push(currPts);
    console.log("Paths: " + joinedPaths.length);

    for (const pts of joinedPaths) {
      let path = new Path({ segments: pts });
      project.activeLayer.addChild(path);
    }
  }

  function animate() {
    let clothChanged = false;
    if (nIters2 <= 0 || iter < nIters2) {
      ++iter;
      cloth.update(cpars.dt);
      cloth.updateGeo();
      clothChanged = true;
      if ((iter % 100) == 0) console.log(iter);
    }
    if (controls) {
      controls.update();
      let newPos = cam.position;
      let newTarget = controls.target;
      let camMoved =
        !nearEq(newPos.x, camProps.position.x) || !nearEq(newPos.y, camProps.position.y) || !nearEq(newPos.z, camProps.position.z) ||
        !nearEq(newTarget.x, camProps.target.x) || !nearEq(newTarget.y, camProps.target.y) || !nearEq(newTarget.z, camProps.target.z);
      if (camMoved) {
        camProps.position.x = twoDecimals(newPos.x);
        camProps.position.y = twoDecimals(newPos.y);
        camProps.position.z = twoDecimals(newPos.z);
        camProps.target.x = twoDecimals(newTarget.x);
        camProps.target.y = twoDecimals(newTarget.y);
        camProps.target.z = twoDecimals(newTarget.z);
        console.log(JSON.stringify(camProps));
      }
    }
    else {
      if (iter == nIters2) {
        renderPlot();
        ++iter;
      }
    }
    if (controls || clothChanged) {
      renderer.render(scene, cam);
      requestAnimationFrame(animate);
    }
  }
}

class FilterableLine {
  constructor(pt1, pt2, clr1, clr2) {
    this.pt1 = pt1;
    this.pt2 = pt2;
    this.clr1 = clr1;
    this.clr2 = clr2;
  }
}

class Cloth {

  constructor(startY, pinCorners) {

    this.colors = null;
    this.mesh = null;
    this.geoPosAttr = null;

    this.nParticles = (cpars.nClothSquares + 1) * (cpars.nClothSquares + 1);
    this.partPosArr = new Float64Array(this.nParticles * 3);
    this.partPrevArr = new Float64Array(this.nParticles * 3);
    this.partForceArr = new Float64Array(this.nParticles * 3);

    this.pins = new Set();

    this.objFun = null;

    this.vtemp = new Vector3();
    this.partDist = cpars.clothSize / cpars.nClothSquares;
    this.diagDist = Math.sqrt(2 * this.partDist * this.partDist);

    let ix = 0;
    for (let z = 0; z <= cpars.nClothSquares; z++) {
      for (let x = 0; x <= cpars.nClothSquares; x++, ix++) {
        let px = -cpars.clothSize / 2 + x * this.partDist;
        let py = startY;
        let pz = -cpars.clothSize / 2 + z * this.partDist;
        this.partPosArr[ix * 3] = px;
        this.partPosArr[ix * 3 + 1] = py;
        this.partPosArr[ix * 3 + 2] = pz;
        this.partPrevArr[ix * 3] = px;
        this.partPrevArr[ix * 3 + 1] = py;
        this.partPrevArr[ix * 3 + 2] = pz;
        this.partForceArr[ix * 3] = 0;
        this.partForceArr[ix * 3 + 1] = 0;
        this.partForceArr[ix * 3 + 2] = 0;
      }
    }
    if (pinCorners) {
      this.pins.add(this.partIx(0, 0));
      this.pins.add(this.partIx(cpars.nClothSquares, 0));
      this.pins.add(this.partIx(0, cpars.nClothSquares));
      this.pins.add(this.partIx(cpars.nClothSquares, cpars.nClothSquares));
    }
  }

  partIx(x, z) {
    return x + z * (cpars.nClothSquares  + 1);
  }

  makeMesh(colorizer) {

    const nTriangles = 2 * cpars.nClothSquares * cpars.nClothSquares;
    this.colors = [];
    for (let i = 0; i < nTriangles; ++i)
      this.colors.push(colorizer.next());

    const geo = new THREE.BufferGeometry();
    const posArr = new Float32Array(nTriangles * 3 * 3);
    this.geoPosAttr = new THREE.BufferAttribute(posArr, 3, false);
    geo.setAttribute('position', this.geoPosAttr);
    this.updateGeo();
    const clrArr = new Float32Array(nTriangles * 3 * 3);
    for (let i = 0; i < nTriangles; ++i) {
      const clr = this.colors[i];
      for (let j = 0; j < 3; ++j) {
        clrArr[i * 9 + j * 3] = clr.r;
        clrArr[i * 9 + j * 3 + 1] = clr.g;
        clrArr[i * 9 + j * 3 + 2] = clr.b;
      }
    }
    geo.setAttribute('a_color', new THREE.BufferAttribute(clrArr, 3, false));
    const mat = new THREE.ShaderMaterial({
      vertexShader: clothVertShader,
      fragmentShader: clothFragShader,
      transparent: false,
      blending: THREE.NoBlending,
      depthTest: true,
      vertexColors: true,
      side: THREE.DoubleSide,
    });
    // const mat = new THREE.MeshBasicMaterial({vertexColors: THREE.VertexColors, side: THREE.DoubleSide});
    this.mesh = new THREE.Mesh(geo, mat);
    return this.mesh;
  }

  updateGeo() {

    // pXZ
    const p00 = new Vector3();
    const p10 = new Vector3();
    const p01 = new Vector3();
    const p11 = new Vector3();

    for (let z = 0; z < cpars.nClothSquares; ++z) {
      for (let x = 0; x < cpars.nClothSquares; ++x) {

        // Coordinates of four particles around this square
        let partIx = this.partIx(x, z) * 3;
        p00.set(this.partPosArr[partIx], this.partPosArr[partIx + 1], this.partPosArr[partIx + 2]);
        partIx = this.partIx(x + 1, z) * 3;
        p10.set(this.partPosArr[partIx], this.partPosArr[partIx + 1], this.partPosArr[partIx + 2]);
        partIx = this.partIx(x, z + 1) * 3;
        p01.set(this.partPosArr[partIx], this.partPosArr[partIx + 1], this.partPosArr[partIx + 2]);
        partIx = this.partIx(x + 1, z + 1) * 3;
        p11.set(this.partPosArr[partIx], this.partPosArr[partIx + 1], this.partPosArr[partIx + 2]);

        // Start of two triangles. Front facing is clockwise, so:
        // 00 - 10 - 01
        // 10 - 11 - 01
        const triIx = (z * cpars.nClothSquares + x) * 2 * 3;
        this.geoPosAttr.setXYZ(triIx, p00.x, p00.y, p00.z);
        this.geoPosAttr.setXYZ(triIx + 1, p10.x, p10.y, p10.z);
        this.geoPosAttr.setXYZ(triIx + 2, p01.x, p01.y, p01.z);
        this.geoPosAttr.setXYZ(triIx + 3, p10.x, p10.y, p10.z);
        this.geoPosAttr.setXYZ(triIx + 4, p11.x, p11.y, p11.z);
        this.geoPosAttr.setXYZ(triIx + 5, p01.x, p01.y, p01.z);
      }
    }
    this.geoPosAttr.needsUpdate = true;
  }

  getLines() {

    const res = {
      // Line segments along the X axis, each line left to right
      horiz: [],
      // Line segments along the Z axiz, each line back to front
      deep: [],
      // Diagonal line segments, each line front left to back right
      diag: [],
    };

    const vert1 = new Vector3();
    const vert2 = new Vector3();

    const getFL = (x1, z1, x2, z2) => {
      const ix1 = this.partIx(x1, z1) * 3;
      const ix2 = this.partIx(x2, z2) * 3;
      vert1.set(this.partPosArr[ix1], this.partPosArr[ix1 + 1], this.partPosArr[ix1 + 2]);
      vert2.set(this.partPosArr[ix2], this.partPosArr[ix2 + 1], this.partPosArr[ix2 + 2]);
      this.mesh.localToWorld(vert1);
      this.mesh.localToWorld(vert2);
      const fl = new FilterableLine;
      fl.pt1 = new Point();
      fl.pt2 = new Point();
      projInPlace(vert1, fl.pt1);
      projInPlace(vert2, fl.pt2);
      return fl;
    }

    const addColors = (fl, clr1Ix, clr2Ix) => {
      if (clr1Ix == -1) clr1Ix = clr2Ix;
      else if (clr2Ix == -1) clr2Ix = clr1Ix;
      fl.clr1 = clrTo8Bit(this.colors[clr1Ix]);
      fl.clr2 = clrTo8Bit(this.colors[clr2Ix]);

      function clrTo8Bit(clr) {
        return {
          r: Math.floor(clr.r >= 1 ? 255 : clr.r * 256),
          g: Math.floor(clr.g >= 1 ? 255 : clr.g * 256),
          b: Math.floor(clr.b >= 1 ? 255 : clr.b * 256),
        };
      }

    }

    // X axis lines
    for (let z = 0; z <= cpars.nClothSquares; ++z) {
      for (let x = 0; x < cpars.nClothSquares; ++x) {
        const fl = getFL(x, z, x + 1, z);
        let clr1Ix = -1, clr2Ix = -1;
        if (z > 0) clr1Ix = (z - 1) * cpars.nClothSquares * 2 + x * 2 + 1;
        if (z < cpars.nClothSquares) clr2Ix = z * cpars.nClothSquares * 2 + x * 2;
        addColors(fl, clr1Ix, clr2Ix);
        res.horiz.push(fl);
      }
    }
    // Z axis lines
    for (let x = 0; x <= cpars.nClothSquares; ++x) {
      for (let z = 0; z < cpars.nClothSquares; ++z) {
        const fl = getFL(x, z, x, z + 1);
        let clr1Ix = -1, clr2Ix = -1;
        if (x > 0) clr1Ix = z * cpars.nClothSquares * 2 + x * 2 - 1;
        if (x < cpars.nClothSquares) clr2Ix = z * cpars.nClothSquares * 2 + x * 2;
        addColors(fl, clr1Ix, clr2Ix);
        res.deep.push(fl);
      }
    }

    // Diagonals
    for (let n = 1; n < 2 * cpars.nClothSquares - 1; ++n) {
      // TODO
    }

    return res;
  }

  updateParticle(ix, dtSq) {
    const ixX = ix * 3;
    const ixY = ixX + 1;
    const ixZ = ixX + 2;
    const px = this.partPosArr[ixX];
    const py = this.partPosArr[ixY];
    const pz = this.partPosArr[ixZ];
    const prevx = this.partPrevArr[ixX];
    const prevy = this.partPrevArr[ixY];
    const prevz = this.partPrevArr[ixZ];
    const forcex = this.partForceArr[ixX];
    const forcey = this.partForceArr[ixY];
    const forcez = this.partForceArr[ixZ];
    const newX = px + (px - prevx) * cpars.dampening + forcex * dtSq;
    const newY = py + (py - prevy) * cpars.dampening + forcey * dtSq;
    const newZ = pz + (pz - prevz) * cpars.dampening + forcez * dtSq;
    this.partPrevArr[ixX] = px;
    this.partPrevArr[ixY] = py;
    this.partPrevArr[ixZ] = pz;
    this.partPosArr[ixX] = newX;
    this.partPosArr[ixY] = newY;
    this.partPosArr[ixZ] = newZ;
    this.partForceArr[ixX] = 0;
    this.partForceArr[ixY] = 0;
    this.partForceArr[ixZ] = 0;
  }

  applyToSprings(x, z, fun) {
    let ix = this.partIx(x, z);
    if (x > 0) {
      fun(ix, this.partIx(x - 1, z), this.partDist, cpars.kStructural);
      if (z < cpars.nClothSquares - 1)
        fun(ix, this.partIx(x - 1, z + 1), this.diagDist, cpars.kShear);
    }
    if (x > 1) {
      fun(ix, this.partIx(x - 2, z), 2 * this.partDist, cpars.kFlex);
      // if (z < cpars.nClothSquares - 2)
      //   fun(ix, this.partIx(x - 2, z + 2), 2 * this.diagDist, kShear);
    }
    if (z > 0) {
      fun(ix, this.partIx(x, z - 1), this.partDist, cpars.kStructural);
      if (x > 0)
        fun(ix, this.partIx(x - 1, z - 1), this.diagDist, cpars.kShear);
    }
    if (z > 1) {
      fun(ix, this.partIx(x, z - 2), 2 * this.partDist, cpars.kFlex);
      // if (x > 1)
      //   fun(ix, this.partIx(x - 2, z - 2), 2 * this.diagDist, kShear);
    }
  }

  resolveSpring(ix1, ix2, length, k) {
    const ix1x = ix1 * 3;
    const ix1y = ix1x + 1;
    const ix1z = ix1x + 2;
    const ix2x = ix2 * 3;
    const ix2y = ix2x + 1;
    const ix2z = ix2x + 2;
    const x1 = this.partPosArr[ix1x];
    const y1 = this.partPosArr[ix1y];
    const z1 = this.partPosArr[ix1z];
    const x2 = this.partPosArr[ix2x];
    const y2 = this.partPosArr[ix2y];
    const z2 = this.partPosArr[ix2z];
    this.vtemp.x = x1 - x2;
    this.vtemp.y = y1 - y2;
    this.vtemp.z = z1 - z2;
    const dist = this.vtemp.length();
    const forceScalar = length / dist - 1; // Negative to bring closer, positive to repel
    this.vtemp.multiplyScalar(k * forceScalar / dist); // 1 / dist is unit length; multiplied by scalar force * spring constant
    this.partForceArr[ix1x] += this.vtemp.x;
    this.partForceArr[ix1y] += this.vtemp.y;
    this.partForceArr[ix1z] += this.vtemp.z;
    this.partForceArr[ix2x] -= this.vtemp.x;
    this.partForceArr[ix2y] -= this.vtemp.y;
    this.partForceArr[ix2z] -= this.vtemp.z;
  }

  update(dt) {

    const dtSq = (dt / cpars.nSimIters) * (dt / cpars.nSimIters);

    for (let i = 0; i < cpars.nSimIters; ++i) {

      for (let ix = 0; ix < this.nParticles; ++ix)
        this.partForceArr[ix * 3 + 1] += -cpars.gravity;

      for (let z = 0; z <= cpars.nClothSquares; ++z) {
        for (let x = 0; x <= cpars.nClothSquares; ++x) {
          this.applyToSprings(x, z, (ix1, ix2, length, k) => {
            this.resolveSpring(ix1, ix2, length, k);
          });
        }
      }

      for (let ix = 0; ix < this.nParticles; ++ix) {

        if (this.pins.has(ix)) continue;

        this.updateParticle(ix, dtSq);
        if (!this.objFun) continue;

        const ixX = ix * 3;
        const ixY = ixX + 1;
        const ixZ = ixX + 2;
        let upd = this.objFun(this.partPosArr[ixX], this.partPosArr[ixY], this.partPosArr[ixZ], this.vtemp);
        if (upd) {
          this.partPosArr[ixX] = this.partPrevArr[ixX];
          this.partPosArr[ixY] = this.partPrevArr[ixY];
          this.partPosArr[ixZ] = this.partPrevArr[ixZ];
        }
      }
    }
  }
}

const clothVertShader = `
precision mediump float;
precision mediump int;
attribute vec4 a_color;
varying vec4 vColor;
void main() {
  gl_Position = projectionMatrix * modelViewMatrix * vec4(position, 1.0);
  vColor = a_color;
}
`;

const clothFragShader = `
precision mediump float;
precision mediump int;
varying vec4 vColor;
void main() {
  vec4 color = vec4(vColor);
  gl_FragColor = vColor;
}
`;


class Colorizer {

  constructor(nHues, nSats, nLights) {
    this.currIx = 0;
    const minSat = 30;
    const maxSat = 80;
    const minLight = 30;
    const maxLight = 80;
    this.colors = [];
    for (let iHue = 0; iHue < nHues; ++iHue) {
      for (let iSat = 0; iSat < nSats; ++iSat) {
        for (let iLight = 0; iLight < nLights; ++iLight) {
          let hue = 360 * iHue / nHues;
          let sat = minSat + (maxSat - minSat) * iSat / nSats;
          let light = minLight + (maxLight - minLight) * iLight / nLights;
          hue = Math.round(hue);
          sat = Math.round(sat);
          light = Math.round(light);
          let str = "hsl(" + hue + ", " + sat + "%, " + light + "%)";
          this.colors.push(new THREE.Color(str));
        }
      }
    }
    shuffle(this.colors);
  }

  next() {
    let res = this.colors[this.currIx];
    this.currIx = (this.currIx + 1) % this.colors.length;
    return res;
  }

  colorFaces(geo) {
    const positionAttribute = geo.getAttribute('position');
    const colors = [];
    let color;
    for (let i = 0; i < positionAttribute.count; ++i) {
      if ((i % 6) == 0) color = this.next();
      colors.push(color.r, color.g, color.b);
    }
    // define the new attribute
    geo.setAttribute('color', new THREE.Float32BufferAttribute(colors, 3));
  }
}

// ===========================================================================
// Three JS machinery
// ===========================================================================


function initThree() {

  const elmPaperCanvas = document.getElementById("paper-canvas");
  const elmCanvasHost = document.getElementById("canvasHost");
  const canvasWidth = elmPaperCanvas.clientWidth;
  const canvasHeight = canvasWidth * h / w;
  const asprat = w / h;

  renderer = new THREE.WebGLRenderer({preserveDrawingBuffer: true});
  elmCanvasHost.appendChild(renderer.domElement);
  elmThreeCanvas = renderer.domElement;
  elmThreeCanvas.id = "three-canvas";
  renderer.setSize(w * rf, h * rf);

  elmCanvasHost.style.width = (canvasWidth * 2) + "px";
  elmPaperCanvas.style.width = canvasWidth + "px";
  elmPaperCanvas.style.position = "relative";
  elmThreeCanvas.style.position = "relative";
  elmThreeCanvas.style.float = "right";
  elmThreeCanvas.style.width = canvasWidth + "px";
  elmThreeCanvas.style.height = canvasHeight + "px";

  let D = w;
  // cam = new THREE.OrthographicCamera(-D, D, D / asprat, -D / asprat, 1, 10000);
  cam = new THREE.PerspectiveCamera(camProps.fov, asprat, 1, 2000);

  cam.position.set(camProps.position.x, camProps.position.y, camProps.position.z);
  cam.lookAt(camProps.target.x, camProps.target.y, camProps.target.z);
  cam.updateProjectionMatrix();

  scene = new THREE.Scene();
  scene.background = new THREE.Color(sceneBgColor);

  ray = new THREE.Raycaster();
}


function proj(vec) {
  let projected = vec.clone().project(cam);
  return [new Point((projected.x + 1) * w / rf, (1 - projected.y) * h / rf), projected.z];
}

function projInPlace(vec, pt) {
  vec.project(cam);
  pt.x = (vec.x + 1) * w / rf;
  pt.y = (1 - vec.y) * h / rf;
}

function addWF(mesh, geo) {
  let wfg = new THREE.WireframeGeometry(geo);
  let wmat = new THREE.LineBasicMaterial({color: 0xeffffff});
  let wf = new THREE.LineSegments(wfg, wmat);
  mesh.add(wf);
}


function nearEq(f, g) {
  if (f == g) return true;
  f = twoDecimals(f);
  g = twoDecimals(g);
  let ratio = f / g;
  return ratio > 0.99 && ratio < 1.01;
}

function twoDecimals(x) {
  return Math.round(x * 100) / 100;
}