raycast.ts

examples/spark/js/raycast.ts · raw · open the example · all examples · site home

Exports: SplatPicker

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
// Ray picking against splats on the CPU, from the loaded ArrayData (Radix
// has no splat raycast; Spark's SplatMesh.raycast does this). Each splat
// near the ray adds its Gaussian's opacity where the ray passes closest to
// it; the hit is where the accumulated opacity first reaches `coverage`, so
// many faint splats (grass, hair) stop the ray like one opaque one.
import { Splats, type ArrayData, type ExtSplat } from "@sparkjsdev/radix";

const SIGMAS = 3;

interface Vec3Like {
  x: number;
  y: number;
  z: number;
}

interface Approach {
  t: number;
  alpha: number;
  index: number;
}

export class SplatPicker {
  decode: (i: number, out: ExtSplat) => ExtSplat;
  coverage: number;
  count: number;
  // Per splat: the center and a bounding radius (-1: skipped)
  centers: Float32Array;
  radii: Float32Array;
  splat: ExtSplat;

  constructor(data: ArrayData, { minOpacity = 0.02, coverage = 0.5 } = {}) {
    const { bundle } = Splats.get(data);
    this.decode = Splats.decoder(bundle);
    this.coverage = coverage;
    this.count = data.count;
    this.centers = new Float32Array(this.count * 3);
    this.radii = new Float32Array(this.count);
    const s = Splats.makeExtSplat();
    for (let i = 0; i < this.count; i += 1) {
      this.decode(i, s);
      this.centers[i * 3] = s.center.x;
      this.centers[i * 3 + 1] = s.center.y;
      this.centers[i * 3 + 2] = s.center.z;
      this.radii[i] = s.opacity < minOpacity ? -1 : SIGMAS * Math.max(s.scales.x, s.scales.y, s.scales.z);
    }
    this.splat = s;
  }

  // The hit of the ray origin + t dir (object space; dir need not be unit)
  // as { t, index }, or undefined
  intersect(origin: Vec3Like, dir: Vec3Like, near = 0, far = Infinity): { t: number; index: number } | undefined {
    const { centers, radii, splat: s } = this;
    const dd = dir.x * dir.x + dir.y * dir.y + dir.z * dir.z;
    const length = Math.sqrt(dd);
    const hits: Approach[] = [];
    for (let i = 0; i < this.count; i += 1) {
      const r = radii[i];
      if (r < 0) {
        continue;
      }
      // Bounding sphere first
      const ox = centers[i * 3] - origin.x;
      const oy = centers[i * 3 + 1] - origin.y;
      const oz = centers[i * 3 + 2] - origin.z;
      const along = (ox * dir.x + oy * dir.y + oz * dir.z) / dd;
      const px = ox - along * dir.x;
      const py = oy - along * dir.y;
      const pz = oz - along * dir.z;
      if (px * px + py * py + pz * pz > r * r || along + r / length < near || along - r / length > far) {
        continue;
      }
      const hit = this.closest(i, s, origin, dir);
      if (hit && hit.t >= near && hit.t <= far) {
        hits.push(hit);
      }
    }
    if (hits.length === 0) {
      return undefined;
    }
    hits.sort((a, b) => a.t - b.t);
    let transmittance = 1;
    let best = hits[0];
    for (const hit of hits) {
      transmittance *= 1 - hit.alpha;
      if (hit.alpha > best.alpha) {
        best = hit;
      }
      if (1 - transmittance >= this.coverage) {
        return { t: hit.t, index: hit.index };
      }
    }
    // A thin cloud: its most opaque point
    return { t: best.t, index: best.index };
  }

  // Where the ray passes closest to splat i, in its Gaussian's metric, and
  // the opacity it meets there: { t, alpha, index }, or undefined. The
  // center is the picker's (edited centers move the splat)
  closest(i: number, s: ExtSplat, origin: Vec3Like, dir: Vec3Like): Approach | undefined {
    this.decode(i, s);
    const { centers } = this;
    const { x: qx, y: qy, z: qz, w: qw } = s.quaternion;
    const sx = Math.max(s.scales.x, 1e-9);
    const sy = Math.max(s.scales.y, 1e-9);
    const sz = Math.max(s.scales.z, 1e-9);
    // Into the splat's frame: rotate by the conjugate, divide by the scales
    const p = rotateConjugate(origin.x - centers[i * 3], origin.y - centers[i * 3 + 1], origin.z - centers[i * 3 + 2], qx, qy, qz, qw);
    const v = rotateConjugate(dir.x, dir.y, dir.z, qx, qy, qz, qw);
    const [px, py, pz] = [p[0] / sx, p[1] / sy, p[2] / sz];
    const [vx, vy, vz] = [v[0] / sx, v[1] / sy, v[2] / sz];
    const vv = vx * vx + vy * vy + vz * vz;
    const t = -(px * vx + py * vy + pz * vz) / vv;
    const [cx, cy, cz] = [px + t * vx, py + t * vy, pz + t * vz];
    const m2 = cx * cx + cy * cy + cz * cz;
    if (m2 > SIGMAS * SIGMAS) {
      return undefined;
    }
    return { t, alpha: Math.min(0.99, s.opacity * Math.exp(-0.5 * m2)), index: i };
  }
}

// v rotated by the conjugate of the unit quaternion q
function rotateConjugate(vx: number, vy: number, vz: number, qx: number, qy: number, qz: number, qw: number) {
  const [x, y, z] = [-qx, -qy, -qz];
  // t = 2 cross(q.xyz, v); v + w t + cross(q.xyz, t)
  const tx = 2 * (y * vz - z * vy);
  const ty = 2 * (z * vx - x * vz);
  const tz = 2 * (x * vy - y * vx);
  return [vx + qw * tx + (y * tz - z * ty), vy + qw * ty + (z * tx - x * tz), vz + qw * tz + (x * ty - y * tx)];
}