index.ts

// SPDX-License-Identifier: Apache-2.0
import type { Fn5 } from "@thi.ng/api";
import { range2d } from "@thi.ng/transducers/range2d";
import type { ReadonlyVec, Vec } from "@thi.ng/vectors";

// flattened [to, clear] tuples
// all positive values are given as times 2
// prettier-ignore
const EDGE_INDEX = [
    -1, -1, 4, 0, 2, 0, 2, 0, 0, 0, -1, -1, 0, 0, 0, 0,
    6, 0, 4, 0, -1, -1, 2, 0, 6, 0, 4, 0, 6, 0, -1, -1
];

// flattened coord offsets [x,y] tuples
const NEXT_EDGES = [0, -1, 1, 0, 0, 1, -1, 0];

// flattened [to,clear] tuples (all values times 2)
const S5 = [4, 8, 0, 2, 0, 26, 4, 14];
const S10 = [6, 4, 2, 16, 6, 22, 2, 28];

export const setBorder = (src: Vec, w: number, h: number, val: number) => {
	const w1 = w - 1;
	const h1 = h - 1;
	const idxH1 = h1 * w;
	for (let x = 0; x < w; x++) {
		src[x] = src[idxH1 + x] = val;
	}
	for (let y = 0; y < h; y++) {
		const yy = y * w;
		src[yy] = src[w1 + yy] = val;
	}
	return src;
};

/** @internal */
const __encodeCrossings = (
	src: ReadonlyVec,
	w: number,
	h: number,
	iso: number
) => {
	const out = new Uint8Array(src.length);
	const w1 = w - 1;
	const h1 = h - 1;
	for (let y = 0, i = 0; y < h1; y++) {
		for (let x = 0; x < w1; i++, x++) {
			out[i] =
				(src[i] < iso ? 16 : 0) |
				(src[i + 1] < iso ? 8 : 0) |
				(src[i + 1 + w] < iso ? 4 : 0) |
				(src[i + w] < iso ? 2 : 0);
		}
		i++;
	}
	return out;
};

/** @internal */
const __cellValue = (src: ReadonlyVec, w: number, idx: number) => {
	return (src[idx] + src[idx + 1] + src[idx + w] + src[idx + w + 1]) * 0.25;
};

/** @internal */
const __mix = (
	src: ReadonlyVec,
	w: number,
	x1: number,
	y1: number,
	x2: number,
	y2: number,
	iso: number
) => {
	const a = src[y1 * w + x1];
	const b = src[y2 * w + x2];
	return a === b ? 0 : (a - iso) / (a - b);
};

// prettier-ignore
/** @internal */
const __contourVertex: Fn5<ReadonlyVec, number, number, number, number, Vec>[] = [
    (src, w, x, y, iso) => [x + __mix(src, w, x, y, x + 1, y, iso), y],
    (src, w, x, y, iso) => [x + 1, y + __mix(src, w, x + 1, y, x + 1, y + 1, iso)],
    (src, w, x, y, iso) => [x + __mix(src, w, x, y + 1, x + 1, y + 1, iso), y + 1],
    (src, w, x, y, iso) => [x, y + __mix(src, w, x, y, x, y + 1, iso)]
];

/**
 * Takes an linearized 2d input field of scalars (i.e. `src` array), its
 * `width`, `height` and computes iso lines (contours) for given `iso` threshold
 * value. The computation is implemented as a generator, yielding 1 point array
 * per found contour. The returned points optionally can be scaled using `scale`
 * factor or vector. The default scale of 1.0 will result in coords in these
 * ranges: `[0.5,w-0.5]` (for X) and `[0.5,h-0.5]` (for Y).
 *
 * @param src -
 * @param w -
 * @param h -
 * @param iso -
 * @param scale -
 */
export function* isolines(
	src: ReadonlyVec,
	w: number,
	h: number,
	iso: number,
	scale: ReadonlyVec | number = 1
) {
	const coded = __encodeCrossings(src, w, h, iso);
	let curr: Vec[] = [];
	let from: number;
	let to = -1;
	let clear: number;
	let x!: number;
	let y!: number;
	const w1 = w - 1;
	const h1 = h - 1;
	const [sx, sy] = typeof scale === "number" ? [scale, scale] : scale;
	const cells = range2d(h, w);
	let next: boolean = true;
	let idx: number;
	while (true) {
		from = to;
		if (next) {
			const c = cells.next();
			if (c.done) break;
			[y, x] = c.value;
			next = false;
		}
		if (x >= w1 || y >= h1) {
			next = true;
			continue;
		}
		const i = y * w + x;
		const id = coded[i]; // * 2
		if (id === 10) {
			idx = (__cellValue(src, w, i) > iso ? 0 : 4) + (from === 6 ? 0 : 2);
			to = S5[idx];
			clear = S5[idx + 1];
		} else if (id === 20) {
			idx =
				__cellValue(src, w, i) > iso
					? from === 0
						? 0
						: 2
					: from === 4
					? 4
					: 6;
			to = S10[idx];
			clear = S10[idx + 1];
		} else {
			to = EDGE_INDEX[id];
			clear = EDGE_INDEX[id + 1];
		}
		if (from === -1 && to > -1 && curr.length > 0) {
			yield curr;
			curr = [];
		}
		if (clear !== -1) {
			coded[i] = clear;
		}
		if (to >= 0) {
			const p = __contourVertex[to >> 1](src, w, x, y, iso);
			p[0] = (p[0] + 0.5) * sx;
			p[1] = (p[1] + 0.5) * sy;
			curr.push(p);
			x += NEXT_EDGES[to];
			y += NEXT_EDGES[to + 1];
		} else {
			next = true;
		}
	}
	if (curr.length > 0) {
		yield curr;
	}
}