Skip to content

Commit 087668d

Browse files
acharneskiacharneski
committed
neat fractal
1 parent 4f0d5d5 commit 087668d

4 files changed

Lines changed: 368 additions & 131 deletions

File tree

‎experiments/irrational_lattice/fft.js‎

Lines changed: 131 additions & 131 deletions
Original file line numberDiff line numberDiff line change
@@ -1,117 +1,117 @@
11
// Lightweight 2D FFT (power-of-two) and 3D surface rendering of the
2-
// log-magnitude spectrum for the irrational lattice field.
2+
// log-magnitude spectrum for the irrational lattice field.
33

4-
// In-place iterative radix-2 Cooley-Tukey FFT on real/imag arrays.
5-
function fft1d(re, im, inverse) {
4+
// In-place iterative radix-2 Cooley-Tukey FFT on real/imag arrays.
5+
function fft1d(re, im, inverse) {
66
const n = re.length;
77
// Bit reversal permutation.
88
for (let i = 1, j = 0; i < n; i++) {
9-
let bit = n >> 1;
10-
for (; j & bit; bit >>= 1) j ^= bit;
11-
j ^= bit;
12-
if (i < j) {
13-
const tr = re[i]; re[i] = re[j]; re[j] = tr;
14-
const ti = im[i]; im[i] = im[j]; im[j] = ti;
15-
}
9+
let bit = n >> 1;
10+
for (; j & bit; bit >>= 1) j ^= bit;
11+
j ^= bit;
12+
if (i < j) {
13+
const tr = re[i]; re[i] = re[j]; re[j] = tr;
14+
const ti = im[i]; im[i] = im[j]; im[j] = ti;
15+
}
1616
}
1717
for (let len = 2; len <= n; len <<= 1) {
18-
const ang = (inverse ? 2 : -2) * Math.PI / len;
19-
const wpr = Math.cos(ang);
20-
const wpi = Math.sin(ang);
21-
for (let i = 0; i < n; i += len) {
22-
let wr = 1, wi = 0;
23-
for (let k = 0; k < len / 2; k++) {
24-
const a = i + k;
25-
const b = i + k + len / 2;
26-
const tr = wr * re[b] - wi * im[b];
27-
const ti = wr * im[b] + wi * re[b];
28-
re[b] = re[a] - tr;
29-
im[b] = im[a] - ti;
30-
re[a] += tr;
31-
im[a] += ti;
32-
const nwr = wr * wpr - wi * wpi;
33-
wi = wr * wpi + wi * wpr;
34-
wr = nwr;
18+
const ang = (inverse ? 2 : -2) * Math.PI / len;
19+
const wpr = Math.cos(ang);
20+
const wpi = Math.sin(ang);
21+
for (let i = 0; i < n; i += len) {
22+
let wr = 1, wi = 0;
23+
for (let k = 0; k < len / 2; k++) {
24+
const a = i + k;
25+
const b = i + k + len / 2;
26+
const tr = wr * re[b] - wi * im[b];
27+
const ti = wr * im[b] + wi * re[b];
28+
re[b] = re[a] - tr;
29+
im[b] = im[a] - ti;
30+
re[a] += tr;
31+
im[a] += ti;
32+
const nwr = wr * wpr - wi * wpi;
33+
wi = wr * wpi + wi * wpr;
34+
wr = nwr;
35+
}
3536
}
36-
}
3737
}
38-
}
38+
}
3939

40-
function largestPow2LE(n) {
40+
function largestPow2LE(n) {
4141
let p = 1;
4242
while (p * 2 <= n) p *= 2;
4343
return p;
44-
}
44+
}
4545

46-
// Compute the centered log-magnitude 2D FFT of a square field (Float32Array,
47-
// length size*size). Downsamples to a power-of-two grid (max `maxN`).
48-
// Returns { mag: Float32Array(N*N), N, min, max }.
49-
export function computeFFT2D(data, size, maxN = 64) {
46+
// Compute the centered log-magnitude 2D FFT of a square field (Float32Array,
47+
// length size*size). Downsamples to a power-of-two grid (max `maxN`).
48+
// Returns { mag: Float32Array(N*N), N, min, max }.
49+
export function computeFFT2D(data, size, maxN = 64) {
5050
let N = largestPow2LE(Math.min(size, maxN));
5151
const re = new Float64Array(N * N);
5252
const im = new Float64Array(N * N);
5353

5454
// Sample/average the source field down to N x N.
5555
const scale = size / N;
5656
for (let y = 0; y < N; y++) {
57-
for (let x = 0; x < N; x++) {
58-
// Nearest-block average for a smoother spectrum.
59-
const sx = Math.floor(x * scale);
60-
const sy = Math.floor(y * scale);
61-
re[y * N + x] = data[sy * size + sx];
62-
}
57+
for (let x = 0; x < N; x++) {
58+
// Nearest-block average for a smoother spectrum.
59+
const sx = Math.floor(x * scale);
60+
const sy = Math.floor(y * scale);
61+
re[y * N + x] = data[sy * size + sx];
62+
}
6363
}
6464

6565
// Rows.
6666
const rowRe = new Float64Array(N);
6767
const rowIm = new Float64Array(N);
6868
for (let y = 0; y < N; y++) {
69-
for (let x = 0; x < N; x++) {
70-
rowRe[x] = re[y * N + x];
71-
rowIm[x] = im[y * N + x];
72-
}
73-
fft1d(rowRe, rowIm, false);
74-
for (let x = 0; x < N; x++) {
75-
re[y * N + x] = rowRe[x];
76-
im[y * N + x] = rowIm[x];
77-
}
69+
for (let x = 0; x < N; x++) {
70+
rowRe[x] = re[y * N + x];
71+
rowIm[x] = im[y * N + x];
72+
}
73+
fft1d(rowRe, rowIm, false);
74+
for (let x = 0; x < N; x++) {
75+
re[y * N + x] = rowRe[x];
76+
im[y * N + x] = rowIm[x];
77+
}
7878
}
7979
// Columns.
8080
const colRe = new Float64Array(N);
8181
const colIm = new Float64Array(N);
8282
for (let x = 0; x < N; x++) {
83-
for (let y = 0; y < N; y++) {
84-
colRe[y] = re[y * N + x];
85-
colIm[y] = im[y * N + x];
86-
}
87-
fft1d(colRe, colIm, false);
88-
for (let y = 0; y < N; y++) {
89-
re[y * N + x] = colRe[y];
90-
im[y * N + x] = colIm[y];
91-
}
83+
for (let y = 0; y < N; y++) {
84+
colRe[y] = re[y * N + x];
85+
colIm[y] = im[y * N + x];
86+
}
87+
fft1d(colRe, colIm, false);
88+
for (let y = 0; y < N; y++) {
89+
re[y * N + x] = colRe[y];
90+
im[y * N + x] = colIm[y];
91+
}
9292
}
9393

9494
// Log magnitude, fftshift to center the DC component.
9595
const mag = new Float32Array(N * N);
9696
let min = Infinity, max = -Infinity;
9797
const half = N / 2;
9898
for (let y = 0; y < N; y++) {
99-
for (let x = 0; x < N; x++) {
100-
const sx = (x + half) % N;
101-
const sy = (y + half) % N;
102-
const idx = y * N + x;
103-
const m = Math.log1p(Math.hypot(re[idx], im[idx]));
104-
mag[sy * N + sx] = m;
105-
if (m < min) min = m;
106-
if (m > max) max = m;
107-
}
99+
for (let x = 0; x < N; x++) {
100+
const sx = (x + half) % N;
101+
const sy = (y + half) % N;
102+
const idx = y * N + x;
103+
const m = Math.log1p(Math.hypot(re[idx], im[idx]));
104+
mag[sy * N + sx] = m;
105+
if (m < min) min = m;
106+
if (m > max) max = m;
107+
}
108108
}
109109
return { mag, N, min, max };
110-
}
110+
}
111111

112-
// Render the FFT log-magnitude as a 3D surface (isometric/oblique projection).
113-
// opts: { rot (deg), tilt (deg), heightScale }.
114-
export function renderFFT3D(canvas, fft, opts = {}) {
112+
// Render the FFT log-magnitude as a 3D surface (isometric/oblique projection).
113+
// opts: { rot (deg), tilt (deg), heightScale }.
114+
export function renderFFT3D(canvas, fft, opts = {}) {
115115
const ctx = canvas.getContext("2d");
116116
const W = canvas.width, H = canvas.height;
117117
ctx.clearRect(0, 0, W, H);
@@ -129,29 +129,29 @@
129129
// Map grid coords -> screen. We compute then fit to canvas.
130130
const cellH = Math.min(W, H) * 0.32 * heightScale;
131131
const project = (gx, gy, gz) => {
132-
// Normalize grid to [-1,1].
133-
const nx = (gx / (N - 1)) * 2 - 1;
134-
const ny = (gy / (N - 1)) * 2 - 1;
135-
// Rotate around vertical axis.
136-
const rx = nx * cosR - ny * sinR;
137-
const ry = nx * sinR + ny * cosR;
138-
// Oblique tilt projection.
139-
const sx = rx;
140-
const sy = ry * cosT - gz * sinT;
141-
return [sx, sy];
132+
// Normalize grid to [-1,1].
133+
const nx = (gx / (N - 1)) * 2 - 1;
134+
const ny = (gy / (N - 1)) * 2 - 1;
135+
// Rotate around vertical axis.
136+
const rx = nx * cosR - ny * sinR;
137+
const ry = nx * sinR + ny * cosR;
138+
// Oblique tilt projection.
139+
const sx = rx;
140+
const sy = ry * cosT - gz * sinT;
141+
return [sx, sy];
142142
};
143143

144144
// First pass: find bounds for auto-fit.
145145
let minSx = Infinity, maxSx = -Infinity, minSy = Infinity, maxSy = -Infinity;
146146
for (let y = 0; y < N; y++) {
147-
for (let x = 0; x < N; x++) {
148-
const z = ((mag[y * N + x] - min) / range) * (cellH / (Math.min(W, H) * 0.5));
149-
const [sx, sy] = project(x, y, z);
150-
if (sx < minSx) minSx = sx;
151-
if (sx > maxSx) maxSx = sx;
152-
if (sy < minSy) minSy = sy;
153-
if (sy > maxSy) maxSy = sy;
154-
}
147+
for (let x = 0; x < N; x++) {
148+
const z = ((mag[y * N + x] - min) / range) * (cellH / (Math.min(W, H) * 0.5));
149+
const [sx, sy] = project(x, y, z);
150+
if (sx < minSx) minSx = sx;
151+
if (sx > maxSx) maxSx = sx;
152+
if (sy < minSy) minSy = sy;
153+
if (sy > maxSy) maxSy = sy;
154+
}
155155
}
156156
const pad = 14;
157157
const spanX = (maxSx - minSx) || 1;
@@ -161,11 +161,11 @@
161161
const oy = pad - minSy * s;
162162

163163
const toScreen = (gx, gy, z) => {
164-
const [sx, sy] = project(gx, gy, z);
165-
return [sx * s + ox, sy * s + oy];
164+
const [sx, sy] = project(gx, gy, z);
165+
return [sx * s + ox, sy * s + oy];
166166
};
167167
const zAt = (x, y) =>
168-
((mag[y * N + x] - min) / range) * (cellH / (Math.min(W, H) * 0.5));
168+
((mag[y * N + x] - min) / range) * (cellH / (Math.min(W, H) * 0.5));
169169

170170
// Painter's algorithm: draw quads from back to front.
171171
// "Back" = larger projected depth (ry after rotation). We approximate by
@@ -175,46 +175,46 @@
175175
for (let y = 0; y < N - 1; y++) ys.push(y);
176176
// Order rows by their average screen-y so nearer ones paint last.
177177
ys.sort((a, b) => {
178-
const za = zAt(0, a), zb = zAt(0, b);
179-
const [, sya] = project(N / 2, a, za);
180-
const [, syb] = project(N / 2, b, zb);
181-
return sya - syb;
178+
const za = zAt(0, a), zb = zAt(0, b);
179+
const [, sya] = project(N / 2, a, za);
180+
const [, syb] = project(N / 2, b, zb);
181+
return sya - syb;
182182
});
183183

184184
for (const y of ys) {
185-
for (let x = 0; x < N - 1; x++) {
186-
const z00 = zAt(x, y);
187-
const z10 = zAt(x + 1, y);
188-
const z11 = zAt(x + 1, y + 1);
189-
const z01 = zAt(x, y + 1);
190-
const p00 = toScreen(x, y, z00);
191-
const p10 = toScreen(x + 1, y, z10);
192-
const p11 = toScreen(x + 1, y + 1, z11);
193-
const p01 = toScreen(x, y + 1, z01);
194-
195-
const t = (z00 + z10 + z11 + z01) / 4 /
196-
((cellH / (Math.min(W, H) * 0.5)) || 1);
197-
const [r, g, b] = spectrumColor(Math.max(0, Math.min(1, t)));
198-
ctx.fillStyle = `rgb(${r},${g},${b})`;
199-
ctx.strokeStyle = "rgba(0,0,0,0.25)";
200-
ctx.lineWidth = 0.5;
201-
ctx.beginPath();
202-
ctx.moveTo(p00[0], p00[1]);
203-
ctx.lineTo(p10[0], p10[1]);
204-
ctx.lineTo(p11[0], p11[1]);
205-
ctx.lineTo(p01[0], p01[1]);
206-
ctx.closePath();
207-
ctx.fill();
208-
ctx.stroke();
209-
}
185+
for (let x = 0; x < N - 1; x++) {
186+
const z00 = zAt(x, y);
187+
const z10 = zAt(x + 1, y);
188+
const z11 = zAt(x + 1, y + 1);
189+
const z01 = zAt(x, y + 1);
190+
const p00 = toScreen(x, y, z00);
191+
const p10 = toScreen(x + 1, y, z10);
192+
const p11 = toScreen(x + 1, y + 1, z11);
193+
const p01 = toScreen(x, y + 1, z01);
194+
195+
const t = (z00 + z10 + z11 + z01) / 4 /
196+
((cellH / (Math.min(W, H) * 0.5)) || 1);
197+
const [r, g, b] = spectrumColor(Math.max(0, Math.min(1, t)));
198+
ctx.fillStyle = `rgb(${r},${g},${b})`;
199+
ctx.strokeStyle = "rgba(0,0,0,0.25)";
200+
ctx.lineWidth = 0.5;
201+
ctx.beginPath();
202+
ctx.moveTo(p00[0], p00[1]);
203+
ctx.lineTo(p10[0], p10[1]);
204+
ctx.lineTo(p11[0], p11[1]);
205+
ctx.lineTo(p01[0], p01[1]);
206+
ctx.closePath();
207+
ctx.fill();
208+
ctx.stroke();
209+
}
210210
}
211-
}
211+
}
212212

213-
// Small viridis-ish colormap for the 3D surface.
214-
function spectrumColor(t) {
213+
// Small viridis-ish colormap for the 3D surface.
214+
function spectrumColor(t) {
215215
const stops = [
216-
[68, 1, 84], [59, 82, 139], [33, 145, 140],
217-
[94, 201, 98], [253, 231, 37],
216+
[68, 1, 84], [59, 82, 139], [33, 145, 140],
217+
[94, 201, 98], [253, 231, 37],
218218
];
219219
const n = stops.length - 1;
220220
const idx = t * n;
@@ -223,8 +223,8 @@
223223
const f = idx - i0;
224224
const a = stops[i0], b = stops[i1];
225225
return [
226-
Math.round(a[0] + (b[0] - a[0]) * f),
227-
Math.round(a[1] + (b[1] - a[1]) * f),
228-
Math.round(a[2] + (b[2] - a[2]) * f),
226+
Math.round(a[0] + (b[0] - a[0]) * f),
227+
Math.round(a[1] + (b[1] - a[1]) * f),
228+
Math.round(a[2] + (b[2] - a[2]) * f),
229229
];
230-
}
230+
}

‎experiments/irrational_lattice/index.html‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -141,6 +141,20 @@ <h1>Algebraic Colored Lattice Field</h1>
141141
<button id="offsetYPlus">y+</button>
142142
</div>
143143
</div>
144+
<div class="control">
145+
<label>Autocorrelation walk</label>
146+
<div class="offset-buttons">
147+
<button id="acCompute" title="Compute top-2 autocorrelation vectors">analyze</button>
148+
<button id="acStep" title="Random walk one step">step</button>
149+
<button id="acPlay" title="Play / pause random walk">▶ play</button>
150+
<button id="acShow" title="Show / hide autocorrelation vectors on lattice">show vectors</button>
151+
</div>
152+
<output id="acOut">no vectors</output>
153+
</div>
154+
<div class="control">
155+
<label>Walk speed: <output id="acSpeedOut">2.0</output> steps/s</label>
156+
<input type="range" id="acSpeed" min="0.1" max="20" value="2" step="0.1">
157+
</div>
144158

145159

146160
<button id="regen">Regenerate</button>

0 commit comments

Comments
 (0)