Source code of plot #060 back to plot

import * as E from "./lib/env.js";
import * as G from "./lib/own/geo2.js"
import {caption} from "./lib/own/caption.js"
import {mulberry32, rand, rand_range, randn_bm, setRandomGenerator, shuffle} from "./lib/own/random.js"
import Clipper2ZFactory from "./lib/thirdparty/clipper2z/clipper2z.js"
import createVoroPP from "./lib/thirdparty/voropp/voropp-module.js"
import {WallPlane, CellData, genVoro} from "./gen-voro.js"

// Declarations below instruct build plugin to copy static files to runtime dir
// STATIC lib/thirdparty/clipper2z/clipper2z.wasm
// STATIC lib/thirdparty/voropp/voropp-module.wasm
// STATIC lib/texture.png

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

const wallA = 1;
const wallB = 0.4;
const wallD = 0.4;
const nPoints = 150;

let seed = Math.round(Math.random() * 65535);
// seed = 48923;

/**
 * @type {MainModule}
 */
let Clipper2Z;

let voroMod;

// Clipper uses integer coordinates.
// Sketch multiplies canvas coordinates by this amount before doing geometry in Clipper.
const clipperMul = 50;

void setup();

async function setup() {
  E.initEnv(w, h, pw, ph);
  console.log(`Seed: ${seed}`);
  E.info("Seed: " + seed);
  setRandomGenerator(mulberry32(seed));
  voroMod = await createVoroPP();
  Clipper2Z = await Clipper2ZFactory();
  await draw();
}

async function draw() {

  caption(w, h, "3d voronoi octahedron", "115");

  const hrange = 2.0;
  const camPos = new G.Vec3(0, 0, -1);
  const proj = new OrthoProjector(hrange, camPos, w, h);
  const yRot = Math.PI * 0.18;
  const xRot = Math.PI * 0.00;
  const maxBlockLevel = 1;
  const equalPtDist = 3;
  const thickLineGap = 2.5;
  const dotLen = 1.5;
  const dotGap = 8;

  const t0 = performance.now();
  // Generate tetrahedron cells via random particles
  const points = genPoints();
  const walls = genTetraWalls();
  const voro = genVoro(voroMod, [-1, 1, -1, 1, -1, 1], walls, points);
  const cells = [];
  for (const cellData of voro) {
    if (cellData.volume < 5e-6) continue;
    const cell = new Cell(cellData);
    cells.push(cell);
    // Rotate
    for (let i = 0; i < cell.vertsAbs.length; ++i) {
      cell.vertsAbs[i].roty(yRot);
      cell.vertsAbs[i].rotx(xRot);
    }
  }
  const t1 = performance.now();
  console.log(`Calculated ${cells.length} Voronoi cells: ${t1-t0} msec`);

  let edgeCount = 0;
  for (const cell of cells) {
    const outlinePts = getProjectedOutline(cell, proj);
    calcCellEdges(proj, cell, outlinePts);
    edgeCount += cell.edges.length;
  }
  const triVerts = [];
  const triIxs = [];
  for (const cell of cells) cell.appendTriangles(triVerts, triIxs);
  console.log(`${cells.length} cells; ${edgeCount} edges; ${triIxs.length} triangles`);

  const triFilter = new TriangleFilter(Math.floor(w/100), Math.floor(h/100));
  for (let i = 0; i < triIxs.length; ++i) {
    const pt1 = proj.project(triVerts[i*3]);
    const pt2 = proj.project(triVerts[i*3+1]);
    const pt3 = proj.project(triVerts[i*3+2]);
    triFilter.addTriangle(pt1, pt2, pt3, i);
  }

  const t2 = performance.now();
  console.log(`Retrieved edges and triangles: ${t2-t1} msec`);

  for (const cell of cells) {
    for (const edge of cell.edges) {
      const visibleSegs = getVisibleSegs(edge, proj, maxBlockLevel, triFilter, triVerts, triIxs);
      for (const vs of visibleSegs) {
        if (vs.isOutline) cell.visibleOutlineSegs.push(vs);
        else cell.visibleInnerSegs.push(vs);
      }
    }
  }
  const t3 = performance.now();
  console.log(`Occlusion testing: ${t3-t2} msec; ${proj.rayTestCount} rays`);

  for (const cell of cells) {
    // Outline as a thick line
    const olSegArrs = joinSegs(cell.visibleOutlineSegs, equalPtDist);
    for (const segs of olSegArrs) {
      const pts = [segs[0].pt1];
      for (const seg of segs) pts.push(seg.pt2);
      drawThickPath(pts, equalPtDist, thickLineGap);
    }
    // Visible and slightly hidden edges
    for (const seg of cell.visibleInnerSegs) {
      if (seg.blockLevel == 0) E.addPath([seg.pt1, seg.pt2]);
      else drawDottedLine(seg.pt1, seg.pt2, dotLen, dotGap);
    }
  }
}

function genPoints() {

  const res = []
  let i = 0;
  while (i < nPoints) {
    const x = genVal();
    const z = genVal();
    let y = rand();
    y = Math.pow(y, 3);
    if (rand() < 0.5) y *= -1;
    y *= 0.5;
    res.push(new G.Vec3(x, y, z));
    ++i;
  }
  return res;

  function genVal() {
    let val = rand_range(-1, 1);
    val = Math.round(val * 100) / 100;
    return val;
  }
}

function genTetraWalls() {
  return [
    new WallPlane(new G.Vec3(wallA, wallB, 0), wallD),
    new WallPlane(new G.Vec3(wallA, -wallB, 0), wallD),
    new WallPlane(new G.Vec3(-wallA, -wallB, 0), wallD),
    new WallPlane(new G.Vec3(-wallA, wallB, 0, wallD), wallD),
    new WallPlane(new G.Vec3(0, wallB, wallA), wallD),
    new WallPlane(new G.Vec3(0, -wallB, wallA), wallD),
    new WallPlane(new G.Vec3(0, -wallB, -wallA), wallD),
    new WallPlane(new G.Vec3(0, wallB, -wallA), wallD),
  ];
}

class VisibleSeg {
  constructor(pt1, pt2, isOutline, blockLevel) {
    this.pt1 = pt1;
    this.pt2 = pt2;
    this.isOutline = isOutline;
    this.blockLevel = blockLevel;
  }
}

function joinSegs(segs, equalPtDist) {
  if (segs.length == 0) return [];
  const res = [];
  let curr = [segs.pop()];
  while (segs.length > 0) {
    const next = getNext(curr[curr.length-1].pt2);
    if (next != null) curr.push(next);
    const prev = getPrev(curr[0].pt1);
    if (prev != null) curr.unshift(prev);
    if (next == null && prev == null) {
      res.push(curr);
      if (segs.length == 0) break;
      curr = [segs.pop()];
    }
  }
  res.push(curr);
  return res;

  function getNext(pt) {
    for (let i = 0; i < segs.length; ++i) {
      const seg = segs[i];
      if (eq(seg.pt2, pt)) [seg.pt1, seg.pt2] = [seg.pt2, seg.pt1];
      if (eq(seg.pt1, pt)) {
        segs.splice(i, 1);
        return seg;
      }
    }
    return null;
  }

  function getPrev(pt) {
    for (let i = 0; i < segs.length; ++i) {
      const seg = segs[i];
      if (eq(seg.pt1, pt)) [seg.pt1, seg.pt2] = [seg.pt2, seg.pt1];
      if (eq(seg.pt2, pt)) {
        segs.splice(i, 1);
        return seg;
      }
    }
    return null;
  }

  function eq(a, b) {
    return G.dist2(a, b) < equalPtDist;
  }
}

function getVisibleSegs(edge, proj, maxBlockLevel, triFilter, triVerts, triIxs) {

  const iv = new G.Vec3();
  const ro = new G.Vec3();
  const ixsToTest = [];
  const stepLen = 1;

  const pt1 = edge.pt1;
  const pt2 = edge.pt2;
  const ab = pt2.clone().sub(pt1);
  const nSteps = Math.round(ab.length() / stepLen);
  if (nSteps < 3) return [];

  const res = [];
  const vw = edge.v2.clone().sub(edge.v1);
  let start = null;
  let nextPt = new G.Vec2(), lastPt = new G.Vec2();
  let blockLevel = -1;
  for (let i = 0; i <= nSteps; ++i) {
    const frac = i / nSteps;
    const e = edge.v1.clone().addMul(vw, frac);
    nextPt.setFrom(pt1).addMul(ab, frac);
    const newBlockLevel = countBlockers(edge, e, nextPt);
    // Current point *very* invisible
    if (newBlockLevel > maxBlockLevel) {
      if (start != null) res.push(new VisibleSeg(start, lastPt.clone(), edge.isOutline, blockLevel));
      start = null;
    }
    // Current point visible
    else {
      // Not building a segment: start now
      if (start == null) start = nextPt.clone();
      // In segment, and blocking level changed: start new segment
      else if (newBlockLevel != blockLevel) {
        res.push(new VisibleSeg(start, lastPt.clone(), edge.isOutline, blockLevel));
        start = nextPt.clone();
      }
    }
    blockLevel = newBlockLevel;
    lastPt.setFrom(nextPt);
  }
  if (start != null && G.dist2(nextPt, start) > 0)
    res.push(new VisibleSeg(start, nextPt.clone(), edge.isOutline, blockLevel));
  return res;

  function countBlockers(edge, v, pt) {

    let count = 0;
    proj.getRayOrigin(pt, ro);
    const vDist = G.dist3(ro, v);
    triFilter.getTriangleIxs(pt, ixsToTest);

    for (let i = 0; i < ixsToTest.length && count <= maxBlockLevel; ++i) {
      const ix = ixsToTest[i];
      const cellId = triIxs[ix].cellId;
      const faceIx = triIxs[ix].faceIx;
      if (edge.cellId == cellId && (edge.faceIx1 == faceIx || edge.faceIx2 == faceIx))
        continue;
      if (!proj.ray(pt, triVerts[ix*3], triVerts[ix*3+1], triVerts[ix*3+2], iv))
        continue;
      const iDist = G.dist3(ro, iv);
      if (iDist < vDist) ++count;
    }
    return count;
  }
}

class ProjectedEdge {
  constructor(v1, v2, pt1, pt2, isOutline, cellId, faceIx1, faceIx2) {
    this.v1 = v1;
    this.v2 = v2;
    this.pt1 = pt1;
    this.pt2 = pt2;
    this.isOutline = isOutline;
    this.cellId = cellId;
    this.faceIx1 = faceIx1;
    this.faceIx2 = faceIx2;
  }
}

function getProjectedOutline(cell, proj) {
  // Project all faces; join resulting 2D polygons
  const projectedFacePolys = [];
  for (const vertIxs of cell.faceVerts) {
    const pts = [];
    vertIxs.forEach(ix => pts.push(proj.project(cell.vertsAbs[ix])));
    projectedFacePolys.push(pts);
  }
  let clUnion = toClipperPolys(projectedFacePolys[0]);
  for (let i = 1; i < projectedFacePolys.length; ++i) {
    const clFace = toClipperPolys(projectedFacePolys[i]);
    clUnion = Clipper2Z.Union64(clUnion, clFace, Clipper2Z.FillRule.NonZero);
  }
  return getClipperPathPoints(clUnion.get(0));
}

function calcCellEdges(proj, cell, outlinePts) {
  const ixMul = 100_000;
  const ixMap = new Map();
  for (let faceIx = 0; faceIx < cell.faceVerts.length; ++faceIx) {
    const indexes = cell.faceVerts[faceIx];
    for (let i = 0; i < indexes.length; ++i) {
      const j = (i+1) % indexes.length;
      let vertIx1 = indexes[i];
      let vertIx2 = indexes[j];
      if (vertIx1 > vertIx2) [vertIx1, vertIx2] = [vertIx2, vertIx1];
      const ixVal = vertIx1 * ixMul + vertIx2;
      if (ixMap.has(ixVal)) {
        const edge = ixMap.get(ixVal);
        edge.faceIx2 = faceIx;
      }
      else ixMap.set(ixVal, { faceIx1: faceIx });
    }
  }
  for (const [ixVal, edge] of ixMap) {
    const ix1 = Math.floor(ixVal / ixMul);
    const ix2 = ixVal - ix1 * ixMul;
    const v1 = cell.vertsAbs[ix1];
    const v2 = cell.vertsAbs[ix2];
    const pt1 = proj.project(v1);
    const pt2 = proj.project(v2);
    const edgeIsOutline = isOutline(pt1, pt2);
    const pe = new ProjectedEdge(v1, v2, pt1, pt2, edgeIsOutline, cell.id, edge.faceIx1, edge.faceIx2);
    cell.edges.push(pe);
  }
  function isOutline(pt1, pt2) {
    const ep = 1e-1;
    for (let i = 0; i < outlinePts.length; ++i) {
      const o1 = outlinePts[i];
      const o2 = outlinePts[(i+1)%outlinePts.length];
      let [d1, d2] = [G.dist2(pt1, o1), G.dist2(pt2, o2)];
      if (d1 < ep && d2 < ep) return true;
      [d1, d2] = [G.dist2(pt1, o2), G.dist2(pt2, o1)];
      if (d1 < ep && d2 < ep) return true;
    }
    return false;
  }
}

// Creates Clipper Paths64 from closed polygon as Vec2 points
function toClipperPolys(pts) {
  const flatPtArr = [];
  for (const pt of pts) flatPtArr.push(Math.round(pt.x * clipperMul), Math.round(pt.y * clipperMul));
  const polys = new Clipper2Z.Paths64();
  polys.push_back(Clipper2Z.MakePath64(flatPtArr));
  return polys;
}

// Converts Clipper Path64 to Vec2 points
function getClipperPathPoints(clPath) {
  const pts = [];
  for (let i = 0; i < clPath.size(); i++) {
    const point = clPath.get(i);
    pts.push(new G.Vec2(Number(point.x) / clipperMul, Number(point.y) / clipperMul));
  }
  return pts;
}

class Cell {

  constructor(data) {

    this.id = data.id;
    this.particle = data.particle;
    this.vertsRel = data.vertices;
    this.vertsAbs = [];
    this.faceVerts = data.faceVertIxs;

    this.edges = [];
    this.visibleOutlineSegs = [];
    this.visibleInnerSegs = [];

    // Offset shards
    if (true) {
      const pnorm = this.particle.clone();
      pnorm.y = 0;
      pnorm.mul(0.2);
      this.particle.add(pnorm);
      this.particle.y *= 1.5;
    }

    // Calculate absolute vertex positiions
    this.vertsRel.forEach(v => this.vertsAbs.push(v.clone().add(this.particle)));
  }

  getFacePts(faceIx) {
    const vertIxs = this.faceVerts[faceIx];
    const res = [];
    for (const vertIx of vertIxs) {
      res.push(this.vertsAbs[vertIx].clone());
    }
    return res;
  }

  appendTriangles(triVerts, triIxs) {
    for (let faceIx = 0; faceIx < this.faceVerts.length; ++faceIx) {
      const indexes = this.faceVerts[faceIx];
      if (indexes.length == 3) {
        triVerts.push(this.vertsAbs[indexes[0]], this.vertsAbs[indexes[1]], this.vertsAbs[indexes[2]]);
        triIxs.push({ cellId: this.id, faceIx: faceIx, triIx: 0});
        continue;
      }
      const center = calcCenter(this.vertsAbs, indexes);
      for (let i = 0; i < indexes.length; ++i) {
        const j = (i+1)%indexes.length;
        triVerts.push(center, this.vertsAbs[indexes[i]], this.vertsAbs[indexes[j]]);
        triIxs.push({ cellId: this.id, faceIx: faceIx, triIx: i});
      }
    }
  }
}

function calcCenter(verts, indexes) {
  const v = verts[indexes[0]].clone();
  for (let i = 1; i < indexes.length; ++i)
    v.add(verts[indexes[i]]);
  return v.mul(1 / indexes.length);
}

function parseInts(str) {
  str = str.replace("(", "").replace(")", "");
  const parts = str.split(",");
  const res = [];
  for (const part of parts)
    res.push(Number.parseInt(part));
  return res;
}

function parseVec3(str) {
  str = str.replace("(", "").replace(")", "");
  const parts = str.split(",");
  return new G.Vec3(
    Number.parseFloat(parts[0]),
    Number.parseFloat(parts[1]),
    Number.parseFloat(parts[2]),
  );
}

function drawDottedLine(a, b, dotLen, dotGap) {
  const len = G.dist2(a, b) - dotLen;
  const nSegs = Math.round(len / (dotLen + dotGap));
  if (nSegs < 2) return null;
  const ab = b.clone().sub(a);
  const dot = ab.clone().mul(dotLen / ab.length());
  const lines = [];
  for (let i = 0; i <= nSegs; ++i) {
    const pt1 = a.clone();
    pt1.x += ab.x * i / nSegs;
    pt1.y += ab.y * i / nSegs;
    const pt2 = pt1.clone().add(dot);
    lines.push([pt1, pt2]);
  }
  return E.addLinesAsPath(lines);
}

function drawThickPath(pts, equalPtDist, lineGap) {
  const closed = pts.length > 2 && G.dist2(pts[0], pts[pts.length-1]) < equalPtDist;
  const clPoly = toClipperPolys(pts);
  if (closed) {
    const clExt = Clipper2Z.InflatePaths64(clPoly, 0.5 * lineGap * clipperMul,
      Clipper2Z.JoinType.Miter, Clipper2Z.EndType.Polygon, 2, 0);
    const ptsExt = getClipperPathPoints(clExt.get(0));
    E.addPath(ptsExt, true);
    const clInt = Clipper2Z.InflatePaths64(clPoly, -0.5 * lineGap * clipperMul,
      Clipper2Z.JoinType.Miter, Clipper2Z.EndType.Polygon, 2, 0);
    const ptsInt = getClipperPathPoints(clInt.get(0));
    E.addPath(ptsInt, true);
  }
  else {
    const clAround = Clipper2Z.InflatePaths64(clPoly, 0.5 * lineGap * clipperMul,
      Clipper2Z.JoinType.Miter, Clipper2Z.EndType.Butt, 2, 0);
    const ptsAround = getClipperPathPoints(clAround.get(0));
    E.addPath(ptsAround, true);
  }
}

class OrthoProjector {
  /**
   * Initializes a new projector.
   * @param range The camera's horizontal range.
   * @param camPos The camera's position. Camera always points at origin.
   * @param canvasWidth Projection canvas width.
   * @param canvasHeight Projection canvas height.
   */
  constructor(range, camPos, canvasWidth, canvasHeight) {
    this.range = range;
    this.camPos = camPos;
    this.canvasWidth = canvasWidth;
    this.canvasHeight = canvasHeight;
    this.aspect = canvasWidth / canvasHeight;
    this.rangeY = range / this.aspect;
    this.camDir = new G.Vec3(0, 0, 0).sub(this.camPos).normalize();
    this.camRight = G.cross3(new G.Vec3(0, 1, 0), this.camDir).normalize();
    this.camUp = G.cross3(this.camDir, this.camRight).normalize();

    // Interim variables for ray function
    this.rayTestCount = 0;
    this._ro = new G.Vec3();
    this._edge1 = new G.Vec3();
    this._edge2 = new G.Vec3();
    this._s = new G.Vec3();
    this._h = new G.Vec3();
    this._q = new G.Vec3();
  }

  project(v) {
    const relativePt = new G.sub3(v, this.camPos);
    const projectedX = relativePt.dot(this.camRight);
    const projectedY = relativePt.dot(this.camUp);
    const canvasX = ((projectedX + this.range / 2) / this.range) * this.canvasWidth;
    const rangeY = this.range / this.aspect;
    const canvasY = ((rangeY / 2 - projectedY) / rangeY) * this.canvasHeight;
    return new G.Vec2(canvasX, canvasY);
  }

  getRayOrigin(pt, ro) {
    // Calculate the 3D coordinates of the canvas point in the orthographic projection
    const prx = this.range * (pt.x / this.canvasWidth - 0.5);
    const pry = this.rangeY * (0.5 - pt.y / this.canvasHeight);
    ro.setFrom(this.camPos);
    ro.x += this.camRight.x * prx + this.camUp.x * pry;
    ro.y += this.camRight.y * prx + this.camUp.y * pry;
    ro.z += this.camRight.z * prx + this.camUp.z * pry;
  }

  ray(pt, v1, v2, v3, iv) {

    ++this.rayTestCount;
    this.getRayOrigin(pt, this._ro);

    // Möller–Trumbore intersection algorithm
    this._edge1.setFrom(v2).sub(v1);
    this._edge2.setFrom(v3).sub(v1);
    this._h.setFrom(this.camDir).cross(this._edge2);
    const a = this._edge1.dot(this._h);

    // Ray is parallel to triangle
    if (a > -1e-6 && a < 1e-6) return false;

    const f = 1 / a;
    this._s.setFrom(this._ro).sub(v1);
    const u = f * this._s.dot(this._h);

    // Intersection point is outside of the triangle
    if (u < 0 || u > 1) return false;

    this._q.setFrom(this._s).cross(this._edge1);
    const v = f * this.camDir.dot(this._q);

    // Intersection point is outside of the triangle
    if (v < 0 || u + v > 1) return false;

    const t = f * this._edge2.dot(this._q);

    // Intersection point is along the ray and in front of the camera
    if (t > 1e-6) {
      iv.setFrom(this._ro);
      iv.x += this.camDir.x * t;
      iv.y += this.camDir.y * t;
      iv.z += this.camDir.z * t;
      return true;
    }
    // Intersection point is behind the camera
    else return false;
  }
}

class TriangleFilter {

  constructor(nHoriz, nVert) {
    this.bucketW = w / nHoriz;
    this.bucketH = h / nVert;
    this.cols = [];
    for (let i = 0; i < nHoriz; ++i) {
      const row = [];
      this.cols.push(row);
      for (let j = 0; j < nVert; ++j) {
        row.push([])
      }
    }
  }

  addTriangle(pt1, pt2, pt3, ix) {

    const getRanges = bounds => {
      let firstColIx = Math.floor(bounds.left / this.bucketW);
      let lastColIx = Math.floor(bounds.right / this.bucketW);
      let firstRowIx = Math.floor(bounds.top / this.bucketH);
      let lastRowIx = Math.floor(bounds.bottom / this.bucketH);
      firstColIx = Math.min(Math.max(firstColIx, 0), this.cols.length - 1);
      lastColIx = Math.min(Math.max(lastColIx, 0), this.cols.length - 1);
      firstRowIx = Math.min(Math.max(firstRowIx, 0), this.cols[0].length - 1);
      lastRowIx = Math.min(Math.max(lastRowIx, 0), this.cols[0].length - 1);
      return [firstColIx, lastColIx, firstRowIx, lastRowIx];
    }

    const bounds = G.bounds2([pt1, pt2, pt3]);
    const [firstColIx, lastColIx, firstRowIx, lastRowIx] = getRanges(bounds);
    for (let colIx = firstColIx; colIx <= lastColIx; ++colIx) {
      const col = this.cols[colIx];
      for (let rowIx = firstRowIx; rowIx <= lastRowIx; ++rowIx) {
        col[rowIx].push(ix);
      }
    }
  }

  getTriangleIxs(pt, ixs) {
    ixs.length = 0;
    const colIx = Math.floor(pt.x / this.bucketW);
    const rowIx = Math.floor(pt.y / this.bucketH);
    ixs.push(...this.cols[colIx][rowIx]);
  }
}