From bb9d21dbc4a93009986af5dd59c9d712a098aa18 Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Tue, 30 Jun 2026 22:11:42 +0100 Subject: [PATCH 01/22] =?UTF-8?q?feat:=20streaming=20video=20=E2=86=92=20l?= =?UTF-8?q?arge=20coherent=20point=20clouds=20(+=20DA3=20studio=20demo,=20?= =?UTF-8?q?Vulkan)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Turn time-lapse video into ONE large coherent point cloud via upstream DA3's sliding-window streaming recipe: run the fused multi-view pass on overlapping windows, then stitch each window into a single global frame with a weighted-Umeyama Sim3 on the overlap. No ggml/model changes — pure host orchestration on top of the existing inference + back_project. Streaming (mirrors da3_streaming, minus loop closure): - src/sim3.{hpp,cpp}: weighted-Umeyama Sim3 — Horn quaternion rotation via the existing linalg jacobi_eigen_sym (reflection-free), optimal Umeyama scale (outlier-robust, unlike the RMS-ratio form), Huber IRLS; compose/apply. - src/stream.{hpp,cpp}: overlapping window loop over depth_pose_multi, dense conf-weighted overlap correspondence (matched by pixel, both-valid), cumulative Sim3 by composition, global accumulation with per-window conf filter + budget subsample. counts are per input frame for build-up. - C API da_capi_points_stream (ABI 5→6); removed the unused da_capi_depth_pose_multi. Bound in the Go server as PointsStream. - server: video bake routes through streaming, lifts the 24-frame cap, exposes chunk/overlap/fps; build-up steps strided to ~30 for long clips. - tests (host-only): test_sim3 (exact recovery, reflection trap, Huber, compose) and test_stream_seam (overlap realignment under disjoint masks). Also lands this session's demo + GPU work it builds on: - DA3 studio server (server/): lazy model registry, points/gaussian bake, .splat export + embedded WebGL EWA viewer, model picker, scene gallery. - Vulkan GPU inference: graph-input upload fix for host norm tensors and skip-unreferenced-input fix (src/dino_backbone.cpp, src/backend.cpp); ggml bumped to v0.15.3 (conv_transpose_2d coopmat2 fix, PR #24924). - gs_head 224×224 guard; multi-view unified on preprocess_real (engine). Verified: host unit tests pass; end-to-end flower-bed clip → 32 frames / 4 stitched windows (overlap M=423360, RMS ~0.01) → 996k-point coherent cloud at 60fps (vs ~508k single-window). Co-Authored-By: Claude Opus 4.8 (1M context) --- .gitignore | 5 + CMakeLists.txt | 2 + include/da_capi.h | 53 ++++- scripts/fetch_models.sh | 34 +++ server/README.md | 93 +++++++++ server/convert.go | 162 ++++++++++++++ server/engine.go | 200 ++++++++++++++++++ server/go.mod | 5 + server/go.sum | 2 + server/main.go | 244 ++++++++++++++++++++++ server/models.go | 147 +++++++++++++ server/scenes.go | 397 +++++++++++++++++++++++++++++++++++ server/web/index.html | 417 +++++++++++++++++++++++++++++++++++++ src/backend.cpp | 9 +- src/da_capi.cpp | 196 ++++++++++++++--- src/dino_backbone.cpp | 10 + src/engine.cpp | 5 +- src/gs_head.cpp | 9 + src/sim3.cpp | 190 +++++++++++++++++ src/sim3.hpp | 32 +++ src/stream.cpp | 242 +++++++++++++++++++++ src/stream.hpp | 38 ++++ tests/CMakeLists.txt | 2 + tests/test_capi.cpp | 2 +- tests/test_sim3.cpp | 119 +++++++++++ tests/test_stream_seam.cpp | 102 +++++++++ third_party/ggml | 2 +- 27 files changed, 2680 insertions(+), 39 deletions(-) create mode 100755 scripts/fetch_models.sh create mode 100644 server/README.md create mode 100644 server/convert.go create mode 100644 server/engine.go create mode 100644 server/go.mod create mode 100644 server/go.sum create mode 100644 server/main.go create mode 100644 server/models.go create mode 100644 server/scenes.go create mode 100644 server/web/index.html create mode 100644 src/sim3.cpp create mode 100644 src/sim3.hpp create mode 100644 src/stream.cpp create mode 100644 src/stream.hpp create mode 100644 tests/test_sim3.cpp create mode 100644 tests/test_stream_seam.cpp diff --git a/.gitignore b/.gitignore index 3dfdefc..e926014 100644 --- a/.gitignore +++ b/.gitignore @@ -12,3 +12,8 @@ __pycache__/ # recorder intermediate outputs /out/ + +# demo server: compiled binary + runtime data + model fetch log +/server/da3-server +/.cache/ +/scripts/.fetch.log diff --git a/CMakeLists.txt b/CMakeLists.txt index 732bfd5..e11d2a9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -63,6 +63,8 @@ set(DA_SOURCES src/colmap_export.cpp src/cam_pose.cpp src/linalg.cpp + src/sim3.cpp + src/stream.cpp src/ray_pose.cpp src/nested.cpp src/depth_export.cpp diff --git a/include/da_capi.h b/include/da_capi.h index 8bb9a9c..b085efa 100644 --- a/include/da_capi.h +++ b/include/da_capi.h @@ -6,7 +6,10 @@ extern "C" { #endif typedef struct da_ctx da_ctx; /* ABI version. 3: added da_capi_depth_dense, da_capi_points, da_capi_free_bytes. - 4: added da_capi_load_nested (two-branch metric model). */ + 4: added da_capi_load_nested (two-branch metric model). + 5: added da_capi_points_multi (fused multi-view cloud) + da_capi_gaussians. + 6: added da_capi_points_stream (sliding-window streaming cloud); removed the + unused da_capi_depth_pose_multi. */ int da_capi_abi_version(void); da_ctx* da_capi_load(const char* gguf_path, int n_threads); /* NULL on failure */ /* Load a NESTED metric model from its two branches: the anyview (GIANT) GGUF and @@ -26,12 +29,6 @@ float* da_capi_depth_path(da_ctx* ctx, const char* image_path, int* out_h, int* void da_capi_free_floats(float* p); /* Run pose; fills ext[12] (3x4 row-major) and intr[9] (3x3). Returns 0 ok, -1 error. */ int da_capi_pose_path(da_ctx* ctx, const char* image_path, float out_ext[12], float out_intr[9]); -/* Multi-view depth+pose. n_images paths. Fills, per view i: out_ext[i*12], out_intr[i*9]. - Returns a malloc'd float[n*H*W] depth (view-major), sets *out_h,*out_w,*out_n; NULL on error. - Caller frees the returned buffer via da_capi_free_floats. */ -float* da_capi_depth_pose_multi(da_ctx* ctx, const char** image_paths, int n_images, - int* out_h, int* out_w, int* out_n, - float* out_ext /* n*12 */, float* out_intr /* n*9 */); /* Single-image 3D export. Runs the native depth+pose pipeline, captures the processed-resolution RGB colors, and writes a glTF-2.0 binary point cloud to out_glb. Returns 0 ok, -1 error (see da_capi_last_error). */ @@ -66,6 +63,48 @@ int da_capi_points(da_ctx* ctx, const char* image_path, float conf_thresh, int* out_n, float** out_xyz, unsigned char** out_rgb); /* Free a uint8 buffer returned by da_capi_points (out_rgb). */ void da_capi_free_bytes(unsigned char* p); + +/* Multi-view FUSED colored point cloud (DualDPT/pose-capable DA3 models; returns + -1 for mono/DA2). Runs ONE cross-view depth+pose pass over n_images and back- + projects every view into a single shared world frame, so the cloud is coherent + across frames (no per-frame scale drift). Recovers per-pixel processed RGB and a + perspective-correct per-point radius (~depth/focal * point_size). Keeps pixels + with confidence >= percentile(conf, conf_pct) (conf_pct in [0,100]). Points are + emitted grouped by source view (frames outer); if out_counts != NULL it must + point to int[n_images] and receives the per-view point count (for progressive + build-up rendering). On success sets *out_n and mallocs *out_xyz[3N] (world xyz, + OpenCV frame), *out_rgb[3N] uint8, *out_radius[N]; free xyz/radius via + da_capi_free_floats and rgb via da_capi_free_bytes. Returns 0 ok, -1 error. */ +int da_capi_points_multi(da_ctx* ctx, const char** image_paths, int n_images, + double conf_pct, float point_size, + int* out_n, int* out_counts, + float** out_xyz, unsigned char** out_rgb, float** out_radius); + +/* Sliding-window STREAMING colored point cloud (DualDPT/pose-capable DA3 models; + returns -1 for mono/DA2). For long ordered frame sequences (time-lapse video): + runs the fused cross-view pass on overlapping windows of chunk_size frames + (sharing `overlap` frames) and stitches each window into ONE global frame via a + weighted-Umeyama Sim3 solved on the overlap, so hundreds of frames fuse into a + single coherent cloud at bounded per-pass memory. conf_pct in [0,100]; + global_budget caps total emitted points (<=0 = unlimited, subsampled per window). + out_counts (if non-NULL) is int[n_images] receiving per-INPUT-frame point counts + (frame-major, for progressive build-up). On success sets *out_n and mallocs + *out_xyz[3N] (world xyz, OpenCV frame), *out_rgb[3N] uint8, *out_radius[N]; free + xyz/radius via da_capi_free_floats and rgb via da_capi_free_bytes. 0 ok, -1 err. */ +int da_capi_points_stream(da_ctx* ctx, const char** image_paths, int n_images, + int chunk_size, int overlap, double conf_pct, + float point_size, int global_budget, + int* out_n, int* out_counts, + float** out_xyz, unsigned char** out_rgb, float** out_radius); + +/* Single-image 3D GAUSSIANS (DA3-GIANT / GS-head models only; returns -1 with a + clear last_error otherwise). Returns world-frame (OpenCV) gaussians as parallel + arrays: *out_xyz[3N] means, *out_scale[3N], *out_quat[4N] (wxyz), *out_rgb[3N] + linear colour in [0,1] from the SH DC term, *out_opacity[N]. Sets *out_n. Free + every returned float buffer via da_capi_free_floats. Returns 0 ok, -1 error. */ +int da_capi_gaussians(da_ctx* ctx, const char* image_path, int* out_n, + float** out_xyz, float** out_scale, float** out_quat, + float** out_rgb, float** out_opacity); #ifdef __cplusplus } #endif diff --git a/scripts/fetch_models.sh b/scripts/fetch_models.sh new file mode 100755 index 0000000..6908f6a --- /dev/null +++ b/scripts/fetch_models.sh @@ -0,0 +1,34 @@ +#!/usr/bin/env bash +# Pre-download the curated DA3 model set for the demo app (~13 GB). +# One weight per capability, not the full 31 GB of redundant quants. +# Idempotent: skips files already present at the right size; resumes partials. +set -uo pipefail + +REPO="mudler/depth-anything.cpp-gguf" +BASE="https://huggingface.co/${REPO}/resolve/main" +DIR="$(cd "$(dirname "$0")/../models" && pwd)" + +# capability file +CURATED=( + "relative+pose (fast) depth-anything-small-f32.gguf" + "relative+pose (default) depth-anything-base-q4_k.gguf" + "relative+pose (quality) depth-anything-large-f32.gguf" + "relative+pose+gaussians depth-anything-giant-f32.gguf" + "metric+pose (anyview branch) depth-anything-nested-anyview.gguf" + "metric+pose (metric branch) depth-anything-nested-metric.gguf" +) + +echo "Downloading curated DA3 models into $DIR" +for row in "${CURATED[@]}"; do + cap="${row% *}"; file="${row##* }" + out="$DIR/$file" + if [ -s "$out" ]; then + echo " [skip] $file (present, $(du -h "$out" | cut -f1))" + continue + fi + echo " [get ] $file — $cap" + curl -fL -C - --retry 5 --retry-delay 3 -o "$out" "$BASE/$file" \ + || { echo " [FAIL] $file"; rm -f "$out"; } +done +echo "Done. Models present:" +ls -lh "$DIR"/*.gguf 2>/dev/null diff --git a/server/README.md b/server/README.md new file mode 100644 index 0000000..606c3eb --- /dev/null +++ b/server/README.md @@ -0,0 +1,93 @@ +# DA3 studio — demo server + +A browser app (in the style of free-splatter.cpp's demo) that turns **videos or +photos into a coherent 3D point cloud** — or **3D gaussians** with DA3-Giant — +using the depth-anything.cpp (DA3) engine, rendered live in a WebGL splat viewer. + +* **Browse sample scenes** — a gallery of baked scenes with frame-by-frame + "build-up" playback (`acc_2 … acc_N`). +* **Create** — upload a video (→ one cross-view DA3 pass → fused cloud) or photos + (→ multi-view cloud, or single-image gaussians), pick the **model** and the + **output type** (point cloud / gaussians), and watch it appear. + +The viewer core (EWA splatting, depth-sort worker, camera) is the same proven +renderer as free-splatter.cpp; it consumes antimatter15 `.splat` (32 B/record). + +## How it works + +``` +video ──ffmpeg──▶ frames ──┐ +photos ────────────────────┴─▶ da_capi_points_multi (ONE cross-view pass) + → per-view depth+pose in one shared world frame + → back-project → fused cloud → .splat (coherent) + +photo ─▶ da_capi_gaussians (DA3-Giant, 224²) → anisotropic gaussians → .splat +``` + +Inference is in-process via **purego** FFI to `libdepthanything.so` (no cgo). +Models are loaded lazily and LRU-evicted (`-max-live`, default 1) so the big +~5 GB checkpoints don't all sit in RAM at once. + +## Build & run + +```bash +# 1. shared library (from the repo root) +cmake -B build -DDA_SHARED=ON -DCMAKE_BUILD_TYPE=Release && cmake --build build -j + +# 2. curated model set (~13 GB, one weight per capability) +bash scripts/fetch_models.sh + +# 3. the server (CGO not needed — purego) +cd server +CGO_ENABLED=0 go build -o da3-server . +./da3-server -lib ../build/libdepthanything.so -models-dir ../models + +# open http://localhost:8794 +``` + +### Flags + +| flag | default | meaning | +|------|---------|---------| +| `-addr` | `:8794` | listen address | +| `-lib` | `../build/libdepthanything.so` | shared library path | +| `-models-dir` | `../models` | directory of `.gguf` weights | +| `-scenes-dir` | `../.cache/da3-scenes` | baked/uploaded scenes | +| `-work-dir` | `/tmp/da3-demo` | scratch for uploads/frames | +| `-threads` | `12` | inference threads | +| `-max-live` | `1` | max resident model contexts (LRU) | +| `-max-splats` | `1500000` | cap splats per output | + +## Models + +Catalogued in `models.go`; only those present on disk appear in the picker. + +| model | output | notes | +|-------|--------|-------| +| `da3-small` / `da3-base` / `da3-large` | point cloud | relative depth + pose, increasing quality | +| `da3-giant` | point cloud **or gaussians** | adds the GS head | +| `da3-nested-metric` | point cloud | anyview + metric branches (real-metre scale, single-image) | + +## API + +| method · path | purpose | +|---|---| +| `GET /api/models` | available models + capabilities + default | +| `POST /api/reconstruct` | photos (+`model`,`mode`) → `{id,n,size,seconds}` | +| `GET /api/splat/{id}` | the `.splat` bytes | +| `GET /api/scenes` | list baked/uploaded scenes | +| `POST /api/scene/from-video` | video (+`model`,`mode`,`max_frames`,`conf_pct`) → `{job}` (async) | +| `GET /api/scene/status/{job}` | `{state,total,done,kept,scene}` | +| `GET /scenes-assets/...` | a scene's `manifest.json` + `.splat` + thumbnails | + +## Known limitations + +* **Point-cloud build-up** slices ONE coherent cross-view pass by frame; it is + bounded by cross-view attention (`-max-frames`, default 12 per pass). Long + time-lapses need windowing + chunk alignment (future work). +* **Gaussians** use the engine's GS head, which is fixed to a **224×224** input + (16×16 patches); the server resizes accordingly. Single-image feed-forward + gaussians are a per-pixel depth surface — sharp from the original viewpoint, + thin edge-on. Multi-view gaussian fusion is future work. +* The `da3-nested-metric` model gives metric scale on the single-image path; the + multi-view fuse currently runs its anyview branch (relative, still coherent). diff --git a/server/convert.go b/server/convert.go new file mode 100644 index 0000000..9085531 --- /dev/null +++ b/server/convert.go @@ -0,0 +1,162 @@ +// convert.go — encode engine output into antimatter15 .splat bytes (32 B/record), +// byte-identical to depth-anything.cpp's viewer contract and free-splatter's +// src/splat.h: the OpenCV(y-down,z-fwd)->OpenGL(y-up) flip is applied here. +package main + +import ( + "encoding/binary" + "math" + "sort" +) + +const splatRow = 32 + +func clamp01(v float32) float32 { + if v < 0 { + return 0 + } + if v > 1 { + return 1 + } + return v +} + +func u8(v float32) byte { + if v < 0 { + v = 0 + } + if v > 255 { + v = 255 + } + return byte(v) +} + +// encodeRecord writes one 32-byte splat. pos/scale in world (OpenCV) units; +// quat as (w,x,y,z); rgb and opacity in [0,1]. +func encodeRecord(out []byte, pos [3]float32, scale [3]float32, quat [4]float32, rgb [3]float32, opacity float32) { + p := [3]float32{pos[0], -pos[1], -pos[2]} + q := [4]float32{-quat[1], quat[0], -quat[3], quat[2]} + nrm := float32(math.Sqrt(float64(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]+q[3]*q[3]))) + 1e-12 + binary.LittleEndian.PutUint32(out[0:], math.Float32bits(p[0])) + binary.LittleEndian.PutUint32(out[4:], math.Float32bits(p[1])) + binary.LittleEndian.PutUint32(out[8:], math.Float32bits(p[2])) + binary.LittleEndian.PutUint32(out[12:], math.Float32bits(scale[0])) + binary.LittleEndian.PutUint32(out[16:], math.Float32bits(scale[1])) + binary.LittleEndian.PutUint32(out[20:], math.Float32bits(scale[2])) + out[24] = u8(clamp01(rgb[0]) * 255) + out[25] = u8(clamp01(rgb[1]) * 255) + out[26] = u8(clamp01(rgb[2]) * 255) + out[27] = u8(clamp01(opacity) * 255) + for c := 0; c < 4; c++ { + out[28+c] = u8(q[c]/nrm*128 + 128) + } +} + +// pointsToSplat encodes the first `n` points of a cloud as isotropic dots. +// quat is identity, opacity 1; n<=0 means all points. +func pointsToSplat(c *Cloud, n int, opacity float32) []byte { + if n <= 0 || n > c.N { + n = c.N + } + buf := make([]byte, n*splatRow) + quat := [4]float32{1, 0, 0, 0} + for i := 0; i < n; i++ { + r := c.Rad[i] + encodeRecord(buf[i*splatRow:], + [3]float32{c.XYZ[3*i], c.XYZ[3*i+1], c.XYZ[3*i+2]}, + [3]float32{r, r, r}, quat, + [3]float32{float32(c.RGB[3*i]) / 255, float32(c.RGB[3*i+1]) / 255, float32(c.RGB[3*i+2]) / 255}, + opacity) + } + return buf +} + +func pctile(v []float32, p float64) float32 { + if len(v) == 0 { + return 0 + } + s := append([]float32(nil), v...) + sort.Slice(s, func(a, b int) bool { return s[a] < s[b] }) + i := int(p/100*float64(len(s)-1) + 0.5) + if i < 0 { + i = 0 + } + if i >= len(s) { + i = len(s) - 1 + } + return s[i] +} + +// gaussiansToSplat encodes a GS set, importance-sorted (opacity*volume) and capped +// at maxSplats. Single-image gaussians include sky / low-confidence pixels that +// back-project to huge depths; we drop that long spatial tail (per-axis 1..99 +// percentile box over reasonably-opaque gaussians) and clamp giant blobs so the +// bulk of the surface frames and renders cleanly. +func gaussiansToSplat(g *GS, maxSplats int) []byte { + N := g.N + var xs, ys, zs, ss []float32 + for i := 0; i < N; i++ { + if g.Opacity[i] < 0.03 { + continue + } + xs = append(xs, g.XYZ[3*i]) + ys = append(ys, g.XYZ[3*i+1]) + zs = append(zs, g.XYZ[3*i+2]) + ss = append(ss, (g.Scale[3*i]+g.Scale[3*i+1]+g.Scale[3*i+2])/3) + } + keepAll := len(xs) < 32 + // Robust RADIAL trim: keep gaussians within the p92 radius of the median + // centroid. This catches scattered floaters (huge back-projected depths) that + // per-axis percentiles miss, so the bulk surface frames cleanly. + var cx, cy, cz, rCut float32 + scaleCap := float32(1e30) + if !keepAll { + cx, cy, cz = pctile(xs, 50), pctile(ys, 50), pctile(zs, 50) + rs := make([]float32, len(xs)) + for i := range xs { + dx, dy, dz := xs[i]-cx, ys[i]-cy, zs[i]-cz + rs[i] = dx*dx + dy*dy + dz*dz + } + rCut = pctile(rs, 92) + scaleCap = pctile(ss, 98) * 3 + } + + idx := make([]int, 0, N) + for i := 0; i < N; i++ { + if g.Opacity[i] < 0.02 { + continue + } + if !keepAll { + dx, dy, dz := g.XYZ[3*i]-cx, g.XYZ[3*i+1]-cy, g.XYZ[3*i+2]-cz + if dx*dx+dy*dy+dz*dz > rCut { + continue + } + } + idx = append(idx, i) + } + imp := func(i int) float32 { + return g.Opacity[i] * g.Scale[3*i] * g.Scale[3*i+1] * g.Scale[3*i+2] + } + sort.Slice(idx, func(a, b int) bool { return imp(idx[a]) > imp(idx[b]) }) + n := len(idx) + if maxSplats > 0 && n > maxSplats { + n = maxSplats + } + clampS := func(v float32) float32 { + if v > scaleCap { + return scaleCap + } + return v + } + buf := make([]byte, n*splatRow) + for k := 0; k < n; k++ { + i := idx[k] + encodeRecord(buf[k*splatRow:], + [3]float32{g.XYZ[3*i], g.XYZ[3*i+1], g.XYZ[3*i+2]}, + [3]float32{clampS(g.Scale[3*i]), clampS(g.Scale[3*i+1]), clampS(g.Scale[3*i+2])}, + [4]float32{g.Quat[4*i], g.Quat[4*i+1], g.Quat[4*i+2], g.Quat[4*i+3]}, + [3]float32{g.RGB[3*i], g.RGB[3*i+1], g.RGB[3*i+2]}, + g.Opacity[i]) + } + return buf +} diff --git a/server/engine.go b/server/engine.go new file mode 100644 index 0000000..1bd83c0 --- /dev/null +++ b/server/engine.go @@ -0,0 +1,200 @@ +// engine.go — purego bindings to libdepthanything (the flat C API in +// include/da_capi.h, ABI 6). No cgo: the shared library is dlopen'd once and each +// model is loaded as its own da_ctx. A context is NOT thread-safe (one ggml +// backend + compute graph), so all inference is serialized by the caller. +package main + +import ( + "fmt" + "runtime" + "unsafe" + + "github.com/ebitengine/purego" +) + +// capi is the dlopen'd library with the DA3 C API bound by name. +type capi struct { + handle uintptr + + abiVersion func() int32 + load func(string, int32) uintptr + loadNested func(string, string, int32) uintptr + freeCtx func(uintptr) + lastErrP func(uintptr) uintptr + freeFloats func(uintptr) + freeBytes func(uintptr) + pointsMulti func(ctx uintptr, paths uintptr, n int32, confPct float64, ptSize float32, + outN, outCounts, outXyz, outRgb, outRad unsafe.Pointer) int32 + pointsStream func(ctx uintptr, paths uintptr, n, chunk, overlap int32, confPct float64, + ptSize float32, budget int32, + outN, outCounts, outXyz, outRgb, outRad unsafe.Pointer) int32 + gaussians func(ctx uintptr, path string, + outN, outXyz, outScale, outQuat, outRgb, outOpacity unsafe.Pointer) int32 +} + +func openCAPI(libPath string) (*capi, error) { + h, err := purego.Dlopen(libPath, purego.RTLD_NOW|purego.RTLD_GLOBAL) + if err != nil { + return nil, fmt.Errorf("dlopen %s: %w", libPath, err) + } + a := &capi{handle: h} + purego.RegisterLibFunc(&a.abiVersion, h, "da_capi_abi_version") + purego.RegisterLibFunc(&a.load, h, "da_capi_load") + purego.RegisterLibFunc(&a.loadNested, h, "da_capi_load_nested") + purego.RegisterLibFunc(&a.freeCtx, h, "da_capi_free") + purego.RegisterLibFunc(&a.lastErrP, h, "da_capi_last_error") + purego.RegisterLibFunc(&a.freeFloats, h, "da_capi_free_floats") + purego.RegisterLibFunc(&a.freeBytes, h, "da_capi_free_bytes") + purego.RegisterLibFunc(&a.pointsMulti, h, "da_capi_points_multi") + purego.RegisterLibFunc(&a.pointsStream, h, "da_capi_points_stream") + purego.RegisterLibFunc(&a.gaussians, h, "da_capi_gaussians") + return a, nil +} + +func (a *capi) lastErr(ctx uintptr) string { return cstr(a.lastErrP(ctx)) } + +// cstr reads a NUL-terminated C string at p (0 -> ""). +func cstr(p uintptr) string { + if p == 0 { + return "" + } + var b []byte + for i := uintptr(0); ; i++ { + c := *(*byte)(unsafe.Pointer(p + i)) + if c == 0 { + break + } + b = append(b, c) + } + return string(b) +} + +func cFloats(p uintptr, n int) []float32 { + if p == 0 || n <= 0 { + return nil + } + out := make([]float32, n) + copy(out, unsafe.Slice((*float32)(unsafe.Pointer(p)), n)) + return out +} + +func cBytes(p uintptr, n int) []byte { + if p == 0 || n <= 0 { + return nil + } + out := make([]byte, n) + copy(out, unsafe.Slice((*byte)(unsafe.Pointer(p)), n)) + return out +} + +// Cloud is a fused multi-view colored point cloud (world frame, OpenCV axes). +type Cloud struct { + N int + XYZ []float32 // 3N + RGB []byte // 3N + Rad []float32 // N (per-point world radius) + Counts []int32 // per source view (frames outer; for build-up prefix sums) +} + +// PointsMulti runs ONE cross-view pass over the given image files and returns the +// fused coherent cloud. confPct in [0,100], ptSize multiplies the per-point radius. +func (a *capi) PointsMulti(ctx uintptr, paths []string, confPct float64, ptSize float32) (*Cloud, error) { + if len(paths) == 0 { + return nil, fmt.Errorf("no images") + } + // Build a C char*[] in Go memory. Go's heap is non-moving, and the C side only + // reads the strings during the (synchronous) call, so this is safe with KeepAlive. + cs := make([][]byte, len(paths)) + pp := make([]*byte, len(paths)) + for i, s := range paths { + b := append([]byte(s), 0) + cs[i] = b + pp[i] = &b[0] + } + var n int32 + counts := make([]int32, len(paths)) + var pXyz, pRgb, pRad uintptr + rc := a.pointsMulti(ctx, uintptr(unsafe.Pointer(&pp[0])), int32(len(paths)), confPct, ptSize, + unsafe.Pointer(&n), unsafe.Pointer(&counts[0]), + unsafe.Pointer(&pXyz), unsafe.Pointer(&pRgb), unsafe.Pointer(&pRad)) + runtime.KeepAlive(cs) + runtime.KeepAlive(pp) + if rc != 0 { + return nil, fmt.Errorf("points_multi: %s", a.lastErr(ctx)) + } + np := int(n) + cl := &Cloud{N: np, Counts: counts, + XYZ: cFloats(pXyz, np*3), RGB: cBytes(pRgb, np*3), Rad: cFloats(pRad, np)} + a.freeFloats(pXyz) + a.freeBytes(pRgb) + a.freeFloats(pRad) + return cl, nil +} + +// PointsStream runs the SLIDING-WINDOW streaming pipeline over an ordered frame +// list (time-lapse video): overlapping fused windows of `chunk` frames (sharing +// `overlap`) stitched into one global cloud. budget caps total points (0 = unlimited). +// Counts is per-input-frame (frame-major) for progressive build-up. +func (a *capi) PointsStream(ctx uintptr, paths []string, chunk, overlap int, confPct float64, ptSize float32, budget int) (*Cloud, error) { + if len(paths) == 0 { + return nil, fmt.Errorf("no images") + } + cs := make([][]byte, len(paths)) + pp := make([]*byte, len(paths)) + for i, s := range paths { + b := append([]byte(s), 0) + cs[i] = b + pp[i] = &b[0] + } + var n int32 + counts := make([]int32, len(paths)) + var pXyz, pRgb, pRad uintptr + rc := a.pointsStream(ctx, uintptr(unsafe.Pointer(&pp[0])), int32(len(paths)), + int32(chunk), int32(overlap), confPct, ptSize, int32(budget), + unsafe.Pointer(&n), unsafe.Pointer(&counts[0]), + unsafe.Pointer(&pXyz), unsafe.Pointer(&pRgb), unsafe.Pointer(&pRad)) + runtime.KeepAlive(cs) + runtime.KeepAlive(pp) + if rc != 0 { + return nil, fmt.Errorf("points_stream: %s", a.lastErr(ctx)) + } + np := int(n) + cl := &Cloud{N: np, Counts: counts, + XYZ: cFloats(pXyz, np*3), RGB: cBytes(pRgb, np*3), Rad: cFloats(pRad, np)} + a.freeFloats(pXyz) + a.freeBytes(pRgb) + a.freeFloats(pRad) + return cl, nil +} + +// GS is a single-image 3D-gaussian set (world frame, OpenCV axes). +type GS struct { + N int + XYZ []float32 // 3N means + Scale []float32 // 3N + Quat []float32 // 4N wxyz + RGB []float32 // 3N linear [0,1] + Opacity []float32 // N +} + +// Gaussians runs the GS reconstruction on one image (DA3-GIANT only). +func (a *capi) Gaussians(ctx uintptr, path string) (*GS, error) { + var n int32 + var pXyz, pScale, pQuat, pRgb, pOp uintptr + rc := a.gaussians(ctx, path, unsafe.Pointer(&n), + unsafe.Pointer(&pXyz), unsafe.Pointer(&pScale), unsafe.Pointer(&pQuat), + unsafe.Pointer(&pRgb), unsafe.Pointer(&pOp)) + if rc != 0 { + return nil, fmt.Errorf("gaussians: %s", a.lastErr(ctx)) + } + np := int(n) + g := &GS{N: np, + XYZ: cFloats(pXyz, np*3), Scale: cFloats(pScale, np*3), Quat: cFloats(pQuat, np*4), + RGB: cFloats(pRgb, np*3), Opacity: cFloats(pOp, np)} + a.freeFloats(pXyz) + a.freeFloats(pScale) + a.freeFloats(pQuat) + a.freeFloats(pRgb) + a.freeFloats(pOp) + return g, nil +} diff --git a/server/go.mod b/server/go.mod new file mode 100644 index 0000000..02a74bf --- /dev/null +++ b/server/go.mod @@ -0,0 +1,5 @@ +module da3.local/server + +go 1.26 + +require github.com/ebitengine/purego v0.10.0 diff --git a/server/go.sum b/server/go.sum new file mode 100644 index 0000000..118766f --- /dev/null +++ b/server/go.sum @@ -0,0 +1,2 @@ +github.com/ebitengine/purego v0.10.0 h1:QIw4xfpWT6GWTzaW5XEKy3HXoqrJGx1ijYHzTF0/ISU= +github.com/ebitengine/purego v0.10.0/go.mod h1:iIjxzd6CiRiOG0UyXP+V1+jWqUXVjPKLAI0mRfJZTmQ= diff --git a/server/main.go b/server/main.go new file mode 100644 index 0000000..4481e91 --- /dev/null +++ b/server/main.go @@ -0,0 +1,244 @@ +// da3 demo server — browse sample scenes or upload a video/photos and turn them +// into a coherent point cloud (any pose-capable DA3 model) or gaussians (giant), +// rendered in the embedded WebGL viewer. Mirrors free-splatter's demo UX. +// +// go run . -lib ../build/libdepthanything.so -models-dir ../models +package main + +import ( + "embed" + "flag" + "fmt" + "io" + "io/fs" + "log" + "net/http" + "os" + "path/filepath" + "strconv" + "sync" + "time" +) + +//go:embed web +var webFS embed.FS + +type server struct { + api *capi + reg *Registry + scenesDir string + workDir string + maxSplats int + + inferMu sync.Mutex // serializes ALL inference (single ggml backend) + + // async video bakes: at most one at a time + jobsMu sync.Mutex + jobs map[string]*sceneJob + bakeSem chan struct{} + + // in-memory LRU of reconstruct results (id -> splat bytes) + resMu sync.Mutex + results map[string][]byte + order []string + counter int64 +} + +func (s *server) infer(fn func() error) error { + s.inferMu.Lock() + defer s.inferMu.Unlock() + return fn() +} + +func main() { + addr := flag.String("addr", ":8794", "listen address") + lib := flag.String("lib", "../build/libdepthanything.so", "path to libdepthanything.so") + modelsDir := flag.String("models-dir", "../models", "directory with the .gguf weights") + scenesDir := flag.String("scenes-dir", "../.cache/da3-scenes", "directory of baked/uploaded scenes") + workDir := flag.String("work-dir", "/tmp/da3-demo", "scratch dir for uploads/frames") + threads := flag.Int("threads", 12, "inference threads") + maxLive := flag.Int("max-live", 1, "max resident model contexts (LRU evicted)") + maxSplats := flag.Int("max-splats", 1500000, "cap splats per output") + flag.Parse() + + api, err := openCAPI(*lib) + if err != nil { + log.Fatalf("load library: %v", err) + } + log.Printf("libdepthanything ABI %d", api.abiVersion()) + + _ = os.MkdirAll(*scenesDir, 0o755) + _ = os.MkdirAll(*workDir, 0o755) + + s := &server{ + api: api, + reg: NewRegistry(api, *modelsDir, *threads, *maxLive), + scenesDir: *scenesDir, + workDir: *workDir, + maxSplats: *maxSplats, + jobs: map[string]*sceneJob{}, + bakeSem: make(chan struct{}, 1), + results: map[string][]byte{}, + } + + avail := s.reg.available() + if len(avail) == 0 { + log.Printf("WARNING: no models found in %s — run scripts/fetch_models.sh", *modelsDir) + } else { + names := make([]string, len(avail)) + for i, m := range avail { + names[i] = m.Name + } + log.Printf("models available: %v", names) + } + + mux := http.NewServeMux() + mux.HandleFunc("/api/models", s.handleModels) + mux.HandleFunc("/api/reconstruct", s.handleReconstruct) + mux.HandleFunc("/api/splat/", s.handleSplat) + mux.HandleFunc("/api/scenes", s.handleScenes) + mux.HandleFunc("/api/scene/from-video", s.handleSceneFromVideo) + mux.HandleFunc("/api/scene/status/", s.handleSceneStatus) + mux.Handle("/scenes-assets/", http.StripPrefix("/scenes-assets/", http.FileServer(http.Dir(*scenesDir)))) + + sub, _ := fs.Sub(webFS, "web") + mux.Handle("/", http.FileServer(http.FS(sub))) + + log.Printf("da3 demo on http://localhost%s (scenes=%s)", *addr, *scenesDir) + log.Fatal(http.ListenAndServe(*addr, mux)) +} + +// GET /api/models +func (s *server) handleModels(w http.ResponseWriter, r *http.Request) { + avail := s.reg.available() + def := "" + for _, m := range avail { + if m.Default { + def = m.Name + } + } + if def == "" && len(avail) > 0 { + def = avail[0].Name + } + writeJSON(w, map[string]any{"models": avail, "default": def}) +} + +// GET /api/splat/{id} +func (s *server) handleSplat(w http.ResponseWriter, r *http.Request) { + id := r.URL.Path[len("/api/splat/"):] + s.resMu.Lock() + b, ok := s.results[id] + s.resMu.Unlock() + if !ok { + http.Error(w, "no such result", 404) + return + } + w.Header().Set("Content-Type", "application/octet-stream") + w.Header().Set("Cache-Control", "no-store") + w.Write(b) +} + +func (s *server) storeResult(b []byte) string { + s.resMu.Lock() + defer s.resMu.Unlock() + s.counter++ + id := "r" + strconv.FormatInt(s.counter, 10) + s.results[id] = b + s.order = append(s.order, id) + for len(s.order) > 16 { + old := s.order[0] + s.order = s.order[1:] + delete(s.results, old) + } + return id +} + +// POST /api/reconstruct (multipart: images[], model, mode) +func (s *server) handleReconstruct(w http.ResponseWriter, r *http.Request) { + if err := r.ParseMultipartForm(256 << 20); err != nil { + http.Error(w, "bad form: "+err.Error(), 400) + return + } + model := r.FormValue("model") + if model == "" { + model = "da3-base" + } + mode := r.FormValue("mode") + if mode == "" { + mode = "points" + } + confPct := atofDefault(r.FormValue("conf_pct"), 55) + ptSize := float32(atofDefault(r.FormValue("point_size"), 1.2)) + + files := r.MultipartForm.File["images"] + if len(files) == 0 { + http.Error(w, "no images", 400) + return + } + dir := filepath.Join(s.workDir, "recon", strconv.FormatInt(time.Now().UnixNano(), 10)) + _ = os.MkdirAll(dir, 0o755) + paths := make([]string, 0, len(files)) + for i, fh := range files { + src, err := fh.Open() + if err != nil { + http.Error(w, err.Error(), 400) + return + } + p := filepath.Join(dir, fmt.Sprintf("img_%03d%s", i, filepath.Ext(fh.Filename))) + dst, err := os.Create(p) + if err != nil { + src.Close() + http.Error(w, err.Error(), 500) + return + } + io.Copy(dst, src) + dst.Close() + src.Close() + paths = append(paths, p) + } + + t0 := time.Now() + var splat []byte + var nOut int + err := s.infer(func() error { + return s.reg.WithModel(model, func(ctx uintptr, mi *ModelInfo) error { + if mode == "gaussians" { + if !mi.Gaussians { + return fmt.Errorf("model %q has no gaussian head (use da3-giant)", model) + } + // GS head is fixed at a 224x224 (16x16 patch) input. + g224 := paths[0] + ".g224.jpg" + if e := squareResize(paths[0], g224, 224); e != nil { + return e + } + g, e := s.api.Gaussians(ctx, g224) + if e != nil { + return e + } + splat = gaussiansToSplat(g, s.maxSplats) + nOut = g.N + return nil + } + c, e := s.api.PointsMulti(ctx, paths, confPct, ptSize) + if e != nil { + return e + } + n := c.N + if s.maxSplats > 0 && n > s.maxSplats { + n = s.maxSplats + } + splat = pointsToSplat(c, n, 1.0) + nOut = n + return nil + }) + }) + if err != nil { + http.Error(w, err.Error(), 500) + return + } + id := s.storeResult(splat) + writeJSON(w, map[string]any{ + "id": id, "model": model, "mode": mode, + "n": nOut, "size": len(splat), "seconds": time.Since(t0).Seconds(), + }) +} diff --git a/server/models.go b/server/models.go new file mode 100644 index 0000000..9f57159 --- /dev/null +++ b/server/models.go @@ -0,0 +1,147 @@ +// models.go — the curated DA3 model catalog and a lazy-loading registry. +// Big checkpoints (giant/nested are ~5 GB each) are loaded on first use and +// evicted LRU so we never hold more than -max-live at once. +package main + +import ( + "fmt" + "log" + "os" + "path/filepath" + "sync" +) + +// ModelInfo describes one selectable model and its capabilities. +type ModelInfo struct { + Name string `json:"name"` + Label string `json:"label"` + Backbone string `json:"backbone"` + File string `json:"-"` // primary gguf + File2 string `json:"-"` // metric branch (nested only) + Pose bool `json:"pose"` // has camera pose -> multi-view fusion + Gaussians bool `json:"gaussians"` // has GS head (gaussian output) + Metric bool `json:"metric"` // absolute metre scale + Modes []string `json:"modes"` // subset of {"points","gaussians"} + Default bool `json:"default"` // pre-selected in the UI +} + +// catalog is the curated set fetched by scripts/fetch_models.sh. +var catalog = []ModelInfo{ + {Name: "da3-small", Label: "DA3 Small — fast", Backbone: "ViT-S", + File: "depth-anything-small-f32.gguf", Pose: true, Modes: []string{"points"}}, + {Name: "da3-base", Label: "DA3 Base — default", Backbone: "ViT-B", + File: "depth-anything-base-q4_k.gguf", Pose: true, Modes: []string{"points"}, Default: true}, + {Name: "da3-large", Label: "DA3 Large — quality", Backbone: "ViT-L", + File: "depth-anything-large-f32.gguf", Pose: true, Modes: []string{"points"}}, + {Name: "da3-giant", Label: "DA3 Giant — points + gaussians", Backbone: "ViT-g", + File: "depth-anything-giant-f32.gguf", Pose: true, Gaussians: true, Modes: []string{"points", "gaussians"}}, + {Name: "da3-nested-metric", Label: "DA3 Nested — metric scale", Backbone: "ViT-g + ViT-L", + File: "depth-anything-nested-anyview.gguf", File2: "depth-anything-nested-metric.gguf", + Pose: true, Metric: true, Modes: []string{"points"}}, +} + +func findModel(name string) *ModelInfo { + for i := range catalog { + if catalog[i].Name == name { + return &catalog[i] + } + } + return nil +} + +// loaded is one resident model context plus a serialization lock. +type loaded struct { + ctx uintptr + runMu sync.Mutex +} + +// Registry lazily loads/evicts model contexts. All inference goes through +// WithModel, which holds the registry lock across load+run so eviction can never +// free a context that is in use (max-live defaults to 1 for the big checkpoints). +type Registry struct { + api *capi + dir string + threads int32 + maxLive int + mu sync.Mutex + live map[string]*loaded + order []string // LRU, oldest first +} + +func NewRegistry(api *capi, dir string, threads, maxLive int) *Registry { + if maxLive < 1 { + maxLive = 1 + } + return &Registry{api: api, dir: dir, threads: int32(threads), maxLive: maxLive, + live: map[string]*loaded{}} +} + +// available returns the catalog entries whose files are present on disk. +func (r *Registry) available() []ModelInfo { + var out []ModelInfo + for _, m := range catalog { + if r.present(&m) { + out = append(out, m) + } + } + return out +} + +func (r *Registry) present(m *ModelInfo) bool { + if _, err := os.Stat(filepath.Join(r.dir, m.File)); err != nil { + return false + } + if m.File2 != "" { + if _, err := os.Stat(filepath.Join(r.dir, m.File2)); err != nil { + return false + } + } + return true +} + +// WithModel loads (or reuses) the named model and runs fn with its context held +// under both the registry lock and the model's run lock. +func (r *Registry) WithModel(name string, fn func(ctx uintptr, mi *ModelInfo) error) error { + r.mu.Lock() + defer r.mu.Unlock() + mi := findModel(name) + if mi == nil { + return fmt.Errorf("unknown model %q", name) + } + if !r.present(mi) { + return fmt.Errorf("model %q not downloaded (run scripts/fetch_models.sh)", name) + } + l, ok := r.live[name] + if !ok { + for len(r.order) >= r.maxLive && len(r.order) > 0 { + ev := r.order[0] + r.order = r.order[1:] + r.api.freeCtx(r.live[ev].ctx) + delete(r.live, ev) + log.Printf("model: evicted %s", ev) + } + var ctx uintptr + if mi.File2 != "" { + ctx = r.api.loadNested(filepath.Join(r.dir, mi.File), filepath.Join(r.dir, mi.File2), r.threads) + } else { + ctx = r.api.load(filepath.Join(r.dir, mi.File), r.threads) + } + if ctx == 0 { + return fmt.Errorf("failed to load model %q", name) + } + l = &loaded{ctx: ctx} + r.live[name] = l + r.order = append(r.order, name) + log.Printf("model: loaded %s", name) + } else { + // move to MRU + for i, n := range r.order { + if n == name { + r.order = append(r.order[:i], r.order[i+1:]...) + break + } + } + r.order = append(r.order, name) + } + return fn(l.ctx, mi) +} diff --git a/server/scenes.go b/server/scenes.go new file mode 100644 index 0000000..93e5de2 --- /dev/null +++ b/server/scenes.go @@ -0,0 +1,397 @@ +// scenes.go — sample-scene listing + the async video->scene pipeline. A video is +// turned into ONE coherent cloud via sliding-window streaming (overlapping fused +// windows stitched by a weighted-Umeyama Sim3), then sliced by frame into acc_* +// "build-up" steps for the viewer. +package main + +import ( + "encoding/json" + "fmt" + "io" + "net/http" + "os" + "os/exec" + "path/filepath" + "regexp" + "sort" + "strconv" + "strings" +) + +type manifestStep struct { + Splat string `json:"splat"` + Images []string `json:"images"` + N int `json:"n"` + Label string `json:"label,omitempty"` +} + +type sceneManifest struct { + Model string `json:"model"` + Mode string `json:"mode"` + Steps []manifestStep `json:"steps"` +} + +type sceneInfo struct { + Name string `json:"name"` + Label string `json:"label"` + Steps int `json:"steps"` + Thumb string `json:"thumb"` + Source string `json:"source"` + Model string `json:"model"` + Mode string `json:"mode"` +} + +type sceneJob struct { + State string `json:"state"` // running | done | error + Total int `json:"total"` + Done int `json:"done"` + Kept int `json:"kept"` + Scene string `json:"scene,omitempty"` + Err string `json:"error,omitempty"` +} + +var slugRe = regexp.MustCompile(`[^a-z0-9_-]+`) + +func slug(s string) string { + s = strings.ToLower(strings.TrimSpace(s)) + s = strings.ReplaceAll(s, " ", "-") + s = slugRe.ReplaceAllString(s, "") + if s == "" { + s = "scene" + } + if len(s) > 48 { + s = s[:48] + } + return s +} + +func writeJSON(w http.ResponseWriter, v any) { + w.Header().Set("Content-Type", "application/json") + _ = json.NewEncoder(w).Encode(v) +} + +// GET /api/scenes +func (s *server) handleScenes(w http.ResponseWriter, r *http.Request) { + ents, _ := os.ReadDir(s.scenesDir) + out := []sceneInfo{} + for _, e := range ents { + if !e.IsDir() { + continue + } + mf := filepath.Join(s.scenesDir, e.Name(), "manifest.json") + b, err := os.ReadFile(mf) + if err != nil { + continue + } + var m sceneManifest + if json.Unmarshal(b, &m) != nil || len(m.Steps) == 0 { + continue + } + thumb := "" + if len(m.Steps[0].Images) > 0 { + thumb = "/scenes-assets/" + e.Name() + "/" + m.Steps[0].Images[0] + } + src := "baked" + if _, err := os.Stat(filepath.Join(s.scenesDir, e.Name(), ".uploaded")); err == nil { + src = "uploaded" + } + out = append(out, sceneInfo{ + Name: e.Name(), Label: prettify(e.Name()), Steps: len(m.Steps), + Thumb: thumb, Source: src, Model: m.Model, Mode: m.Mode, + }) + } + sort.Slice(out, func(i, j int) bool { return out[i].Name < out[j].Name }) + writeJSON(w, map[string]any{"scenes": out}) +} + +func prettify(s string) string { + return strings.Title(strings.ReplaceAll(strings.ReplaceAll(s, "-", " "), "_", " ")) +} + +// GET /api/scene/status/{job} +func (s *server) handleSceneStatus(w http.ResponseWriter, r *http.Request) { + id := strings.TrimPrefix(r.URL.Path, "/api/scene/status/") + s.jobsMu.Lock() + j, ok := s.jobs[id] + s.jobsMu.Unlock() + if !ok { + http.Error(w, "no such job", 404) + return + } + writeJSON(w, j) +} + +// POST /api/scene/from-video (multipart: video, name, model, mode, max_frames, conf_pct, point_size) +func (s *server) handleSceneFromVideo(w http.ResponseWriter, r *http.Request) { + if err := r.ParseMultipartForm(512 << 20); err != nil { + http.Error(w, "bad form: "+err.Error(), 400) + return + } + file, hdr, err := r.FormFile("video") + if err != nil { + http.Error(w, "missing video", 400) + return + } + defer file.Close() + name := slug(r.FormValue("name")) + if name == "scene" && hdr != nil { + name = slug(strings.TrimSuffix(hdr.Filename, filepath.Ext(hdr.Filename))) + } + model := r.FormValue("model") + if model == "" { + model = "da3-base" + } + mode := r.FormValue("mode") + if mode == "" { + mode = "points" + } + maxFrames := atoiDefault(r.FormValue("max_frames"), 60) + if maxFrames < 2 { + maxFrames = 2 + } + if maxFrames > 600 { + maxFrames = 600 + } + chunkSize := atoiDefault(r.FormValue("chunk_size"), 12) + if chunkSize < 2 { + chunkSize = 2 + } + if chunkSize > 24 { + chunkSize = 24 + } + overlap := atoiDefault(r.FormValue("overlap"), 3) + if overlap < 0 { + overlap = 0 + } + if overlap > chunkSize-1 { + overlap = chunkSize - 1 + } + fps := atofDefault(r.FormValue("fps"), 6) + if fps < 0.5 { + fps = 0.5 + } + if fps > 30 { + fps = 30 + } + confPct := atofDefault(r.FormValue("conf_pct"), 55) + ptSize := float32(atofDefault(r.FormValue("point_size"), 1.2)) + + // save upload + jobDir := filepath.Join(s.workDir, "uploads", name) + _ = os.MkdirAll(jobDir, 0o755) + vpath := filepath.Join(jobDir, "input"+filepath.Ext(hdr.Filename)) + out, err := os.Create(vpath) + if err != nil { + http.Error(w, err.Error(), 500) + return + } + if _, err := io.Copy(out, file); err != nil { + out.Close() + http.Error(w, err.Error(), 500) + return + } + out.Close() + + s.jobsMu.Lock() + s.jobs[name] = &sceneJob{State: "running"} + s.jobsMu.Unlock() + + go func() { + s.bakeSem <- struct{}{} + defer func() { <-s.bakeSem }() + err := s.bakeVideo(name, vpath, jobDir, model, mode, maxFrames, chunkSize, overlap, fps, confPct, ptSize) + s.jobsMu.Lock() + j := s.jobs[name] + if err != nil { + j.State, j.Err = "error", err.Error() + } else { + j.State, j.Scene = "done", name + } + s.jobsMu.Unlock() + }() + + writeJSON(w, map[string]string{"job": name, "name": name}) +} + +// bakeVideo: ffmpeg -> frames -> DA3 -> acc_*.splat + thumbnails + manifest. +func (s *server) bakeVideo(name, vpath, jobDir, model, mode string, maxFrames, chunkSize, overlap int, fps, confPct float64, ptSize float32) error { + framesDir := filepath.Join(jobDir, "frames") + _ = os.RemoveAll(framesDir) + _ = os.MkdirAll(framesDir, 0o755) + // Bounded extraction (fps cap + downscale), then even stride to maxFrames. + cmd := exec.Command("ffmpeg", "-y", "-loglevel", "error", "-i", vpath, + "-vf", fmt.Sprintf("fps=%g,scale=640:-2", fps), filepath.Join(framesDir, "all%05d.jpg")) + if b, err := cmd.CombinedOutput(); err != nil { + return fmt.Errorf("ffmpeg: %v: %s", err, string(b)) + } + all, _ := filepath.Glob(filepath.Join(framesDir, "all*.jpg")) + sort.Strings(all) + if len(all) < 2 { + return fmt.Errorf("video produced %d frames (need >=2)", len(all)) + } + stride := (len(all) + maxFrames - 1) / maxFrames + if stride < 1 { + stride = 1 + } + var sel []string + for i := 0; i < len(all) && len(sel) < maxFrames; i += stride { + sel = append(sel, all[i]) + } + + sceneDir := filepath.Join(s.scenesDir, name) + _ = os.RemoveAll(sceneDir) + _ = os.MkdirAll(sceneDir, 0o755) + _ = os.WriteFile(filepath.Join(sceneDir, ".uploaded"), []byte("1"), 0o644) + + s.setJob(name, func(j *sceneJob) { j.Total = len(sel) }) + + mi := findModel(model) + if mi == nil { + return fmt.Errorf("unknown model %q", model) + } + + // Gaussian mode: single representative (middle) frame -> GS splat. + if mode == "gaussians" { + if !mi.Gaussians { + return fmt.Errorf("model %q has no gaussian head (use da3-giant)", model) + } + mid := sel[len(sel)/2] + thumb := "view_0.jpg" + _ = copyScaled(mid, filepath.Join(sceneDir, thumb), 360) + // GS head is fixed at a 224x224 (16x16 patch) input. + g224 := filepath.Join(jobDir, "g224.jpg") + if e := squareResize(mid, g224, 224); e != nil { + return e + } + var splat []byte + err := s.infer(func() error { + return s.reg.WithModel(model, func(ctx uintptr, _ *ModelInfo) error { + g, e := s.api.Gaussians(ctx, g224) + if e != nil { + return e + } + splat = gaussiansToSplat(g, s.maxSplats) + return nil + }) + }) + if err != nil { + return err + } + if e := os.WriteFile(filepath.Join(sceneDir, "acc_full.splat"), splat, 0o644); e != nil { + return e + } + s.setJob(name, func(j *sceneJob) { j.Done, j.Kept = len(sel), 1 }) + return writeManifest(sceneDir, sceneManifest{Model: model, Mode: mode, Steps: []manifestStep{ + {Splat: "acc_full.splat", Images: []string{thumb}, N: 1, Label: "gaussians · single frame"}}}) + } + + // Point-cloud mode: sliding-window streaming over all selected frames. + if !mi.Pose { + return fmt.Errorf("model %q has no camera pose; cannot fuse frames", model) + } + // thumbnails first (so progress feels live) + thumbs := make([]string, len(sel)) + for i, f := range sel { + thumbs[i] = fmt.Sprintf("view_%d.jpg", i) + _ = copyScaled(f, filepath.Join(sceneDir, thumbs[i]), 360) + } + + var cloud *Cloud + err := s.infer(func() error { + return s.reg.WithModel(model, func(ctx uintptr, _ *ModelInfo) error { + c, e := s.api.PointsStream(ctx, sel, chunkSize, overlap, confPct, ptSize, s.maxSplats) + if e != nil { + return e + } + cloud = c + return nil + }) + }) + if err != nil { + return err + } + s.setJob(name, func(j *sceneJob) { j.Done = len(sel); j.Kept = len(sel) }) + + // Build-up: acc_k.splat = points from the first k views (frames outer). For long + // clips, stride the steps so the viewer gets ~30 progressive reveals (always + // including the final full cloud), not one file per frame. + steps := []manifestStep{} + prefix := 0 + capN := s.maxSplats + nv := cloud.N0(len(sel)) + bstep := (nv + 29) / 30 + if bstep < 1 { + bstep = 1 + } + for k := 1; k <= nv; k++ { + prefix += int(cloud.Counts[k-1]) + if k < 2 { + continue + } + if k%bstep != 0 && k != nv { + continue + } + n := prefix + if capN > 0 && n > capN { + n = capN + } + fn := fmt.Sprintf("acc_%d.splat", k) + if e := os.WriteFile(filepath.Join(sceneDir, fn), pointsToSplat(cloud, n, 1.0), 0o644); e != nil { + return e + } + label := "" + if k == nv { + label = fmt.Sprintf("full %d-view cloud · %d pts", k, prefix) + } + steps = append(steps, manifestStep{Splat: fn, Images: thumbs[:k], N: k, Label: label}) + } + return writeManifest(sceneDir, sceneManifest{Model: model, Mode: mode, Steps: steps}) +} + +// N0 guards Counts length vs the number of selected views. +func (c *Cloud) N0(nViews int) int { + if len(c.Counts) < nViews { + return len(c.Counts) + } + return nViews +} + +func writeManifest(dir string, m sceneManifest) error { + b, _ := json.MarshalIndent(m, "", " ") + return os.WriteFile(filepath.Join(dir, "manifest.json"), b, 0o644) +} + +func copyScaled(src, dst string, px int) error { + cmd := exec.Command("ffmpeg", "-y", "-loglevel", "error", "-i", src, + "-vf", fmt.Sprintf("scale=%d:-2", px), dst) + return cmd.Run() +} + +// squareResize center-crops to a px*px square (the GS head's fixed input size). +func squareResize(src, dst string, px int) error { + cmd := exec.Command("ffmpeg", "-y", "-loglevel", "error", "-i", src, + "-vf", fmt.Sprintf("scale=%d:%d:force_original_aspect_ratio=increase,crop=%d:%d", px, px, px, px), dst) + return cmd.Run() +} + +func (s *server) setJob(name string, fn func(*sceneJob)) { + s.jobsMu.Lock() + if j := s.jobs[name]; j != nil { + fn(j) + } + s.jobsMu.Unlock() +} + +func atoiDefault(s string, d int) int { + if v, err := strconv.Atoi(strings.TrimSpace(s)); err == nil { + return v + } + return d +} + +func atofDefault(s string, d float64) float64 { + if v, err := strconv.ParseFloat(strings.TrimSpace(s), 64); err == nil { + return v + } + return d +} diff --git a/server/web/index.html b/server/web/index.html new file mode 100644 index 0000000..63dfc82 --- /dev/null +++ b/server/web/index.html @@ -0,0 +1,417 @@ + + + + + + +depth-anything.cpp — DA3 studio + + + + +
+

depth-anything.cpp · DA3 studio

+
video → coherent point cloud / gaussians, in the browser
+ +
+
Scenes
+
Create
+
+ + +
+ +
+ + + — +
+
+ + +
+ + +
+ + +
+ + +
+ + +
+ + +
+
click or drop a video here
+ + +
+
+
+
+
+
+ + +
+ +
+
loaded0%
+
splats0
+
fps—
+
ready
+
+
drag to orbit · scroll to zoom · WASD/arrows to move
+ + + + diff --git a/src/backend.cpp b/src/backend.cpp index 4e98c20..93bc5bf 100644 --- a/src/backend.cpp +++ b/src/backend.cpp @@ -322,8 +322,15 @@ bool Backend::compute(const std::function& build, return false; } - // Inputs are allocated now (->buffer/->data set): push host data in. + // Inputs are allocated now (->buffer/->data set): push host data in. Skip any + // registered input the graph never referenced — the allocator leaves it without + // a buffer (e.g. cam_src_in when forward_mv runs with a single view), and + // uploading to it would assert "tensor buffer not set". Unreachable inputs feed + // nothing, so there is nothing to upload. for (const PendingInput& pi : impl_->pending) { + ggml_backend_buffer_t ib = pi.tensor->view_src ? pi.tensor->view_src->buffer + : pi.tensor->buffer; + if (!ib) continue; ggml_backend_tensor_set(pi.tensor, pi.host, 0, pi.nbytes); } impl_->pending.clear(); diff --git a/src/da_capi.cpp b/src/da_capi.cpp index c5b2081..2ac0db5 100644 --- a/src/da_capi.cpp +++ b/src/da_capi.cpp @@ -5,7 +5,9 @@ #include "glb_export.hpp" #include "colmap_export.hpp" #include "reconstruct.hpp" +#include "stream.hpp" #include +#include #include #include #include @@ -68,7 +70,7 @@ static bool capi_run_nested(da_ctx* c, const char* image_path, } extern "C" { -int da_capi_abi_version(void){ return 4; } +int da_capi_abi_version(void){ return 6; } da_ctx* da_capi_load(const char* path, int n_threads){ if (!path) return nullptr; auto e = da::Engine::load(path, n_threads); @@ -120,33 +122,6 @@ int da_capi_pose_path(da_ctx* c, const char* image_path, float out_ext[12], floa if (out_intr) std::memcpy(out_intr, intr.data(), 9 * sizeof(float)); return 0; } -float* da_capi_depth_pose_multi(da_ctx* c, const char** image_paths, int n_images, - int* out_h, int* out_w, int* out_n, - float* out_ext, float* out_intr){ - if (!c || !c->engine || !image_paths || n_images <= 0){ if (c) c->last_error = "depth_multi: bad args"; return nullptr; } - std::vector imgs(n_images); - for (int i = 0; i < n_images; ++i){ - if (!image_paths[i] || !da::load_image_rgb(image_paths[i], imgs[i])){ - c->last_error = "depth_multi: load image failed"; return nullptr; - } - } - std::vector views; int H = 0, W = 0; - if (!c->engine->depth_pose_multi(imgs, views, H, W)){ c->last_error = "depth_multi: failed"; return nullptr; } - const int n = (int)views.size(); - const size_t per = (size_t)H * W; - float* p = (float*)std::malloc((size_t)n * per * sizeof(float)); - if (!p){ c->last_error = "depth_multi: oom"; return nullptr; } - for (int i = 0; i < n; ++i){ - std::memcpy(p + (size_t)i * per, views[i].depth.data(), per * sizeof(float)); - if (out_ext) std::memcpy(out_ext + (size_t)i * 12, views[i].ext.data(), 12 * sizeof(float)); - if (out_intr) std::memcpy(out_intr + (size_t)i * 9, views[i].intr.data(), 9 * sizeof(float)); - } - if (out_h) *out_h = H; - if (out_w) *out_w = W; - if (out_n) *out_n = n; - return p; -} - // Shared single-image export prep: run native depth+pose, capture processed RGB, // build N=1 exporter inputs. Returns false (with c->last_error set) on failure. static bool capi_export_prep(da_ctx* c, const char* image_path, @@ -300,4 +275,169 @@ int da_capi_points(da_ctx* c, const char* image_path, float conf_thresh, return 0; } void da_capi_free_bytes(unsigned char* p){ std::free(p); } + +int da_capi_points_multi(da_ctx* c, const char** paths, int n_images, + double conf_pct, float point_size, + int* out_n, int* out_counts, + float** out_xyz, unsigned char** out_rgb, float** out_radius){ + if (out_xyz) *out_xyz = nullptr; + if (out_rgb) *out_rgb = nullptr; + if (out_radius) *out_radius = nullptr; + if (!c || !c->engine || !paths || n_images <= 0){ if (c) c->last_error = "points_multi: bad args"; return -1; } + if (c->engine->is_mono() || c->engine->is_da2()){ + c->last_error = "points_multi: model has no camera pose; use a DualDPT DA3 model"; return -1; } + if (!(point_size > 0.f)) point_size = 1.f; + + std::vector imgs(n_images); + for (int i = 0; i < n_images; ++i){ + if (!paths[i] || !da::load_image_rgb(paths[i], imgs[i])){ c->last_error = "points_multi: load image failed"; return -1; } + } + std::vector views; int H = 0, W = 0; + if (!c->engine->depth_pose_multi(imgs, views, H, W)){ c->last_error = "points_multi: depth_pose_multi failed"; return -1; } + const int N = (int)views.size(); + const size_t plane = (size_t)H * (size_t)W; + + std::vector depth_all; depth_all.reserve((size_t)N * plane); + std::vector conf_all; conf_all.reserve((size_t)N * plane); + std::vector> K(N); + std::vector> E(N); + std::vector> rgb_store(N); + std::vector images_u8(N); + bool have_conf = true; + for (int i = 0; i < N; ++i){ + depth_all.insert(depth_all.end(), views[i].depth.begin(), views[i].depth.end()); + if (views[i].conf.size() == plane) conf_all.insert(conf_all.end(), views[i].conf.begin(), views[i].conf.end()); + else have_conf = false; + K[i] = views[i].intr; + std::array e4{}; for (int k = 0; k < 12; ++k) e4[k] = views[i].ext[k]; + e4[12] = 0.f; e4[13] = 0.f; e4[14] = 0.f; e4[15] = 1.f; E[i] = e4; + // Capture the processed RGB (same resize policy depth_pose_multi used) for + // per-point colour; preprocess_real fills rgb_store directly. + da::Preprocessed pp; + if (!da::preprocess_real(imgs[i], c->engine->config(), pp, &rgb_store[i]) || pp.H != H || pp.W != W){ + c->last_error = "points_multi: preprocess color mismatch"; return -1; } + images_u8[i] = rgb_store[i].data(); + } + if (!have_conf) conf_all.clear(); + float conf_thr = -1e30f; + if (have_conf){ + if (conf_pct < 0) conf_pct = 0; if (conf_pct > 100) conf_pct = 100; + conf_thr = (float)da::percentile_linear(conf_all, conf_pct); + } + + da::WorldPoints wp = da::back_project(depth_all, conf_all, K, E, images_u8, H, W, N, conf_thr); + const size_t np = wp.xyz.size() / 3; + if (np == 0){ c->last_error = "points_multi: no points survived (raise conf_pct or check parallax)"; return -1; } + + float* xyz = (float*)std::malloc(np * 3 * sizeof(float)); + unsigned char* rgb = (unsigned char*)std::malloc(np * 3); + float* rad = (float*)std::malloc(np * sizeof(float)); + if (!xyz || !rgb || !rad){ std::free(xyz); std::free(rgb); std::free(rad); c->last_error = "points_multi: oom"; return -1; } + if (out_counts) for (int i = 0; i < n_images; ++i) out_counts[i] = 0; + for (size_t j = 0; j < np; ++j){ + xyz[3*j+0] = wp.xyz[3*j+0]; xyz[3*j+1] = wp.xyz[3*j+1]; xyz[3*j+2] = wp.xyz[3*j+2]; + rgb[3*j+0] = wp.rgb[3*j+0]; rgb[3*j+1] = wp.rgb[3*j+1]; rgb[3*j+2] = wp.rgb[3*j+2]; + int f = wp.frame[j], u = wp.u[j], v = wp.v[j]; + float d = depth_all[(size_t)f * plane + (size_t)v * W + (size_t)u]; + float fx = K[f][0], fy = K[f][4]; + float rr = 0.5f * (d / fx + d / fy) * point_size; + if (!(rr > 0.f) || !std::isfinite(rr)) rr = 1e-4f; + rad[j] = rr; + if (out_counts && f >= 0 && f < n_images) out_counts[f]++; + } + if (out_xyz) *out_xyz = xyz; else std::free(xyz); + if (out_rgb) *out_rgb = rgb; else std::free(rgb); + if (out_radius) *out_radius = rad; else std::free(rad); + if (out_n) *out_n = (int)np; + return 0; +} + +int da_capi_points_stream(da_ctx* c, const char** image_paths, int n_images, + int chunk_size, int overlap, double conf_pct, + float point_size, int global_budget, + int* out_n, int* out_counts, + float** out_xyz, unsigned char** out_rgb, float** out_radius){ + if (out_xyz) *out_xyz = nullptr; + if (out_rgb) *out_rgb = nullptr; + if (out_radius) *out_radius = nullptr; + if (!c || !c->engine || !image_paths || n_images <= 0){ if (c) c->last_error = "points_stream: bad args"; return -1; } + if (c->engine->is_mono() || c->engine->is_da2()){ + c->last_error = "points_stream: model has no camera pose; use a DualDPT DA3 model"; return -1; } + + std::vector paths(n_images); + for (int i = 0; i < n_images; ++i){ + if (!image_paths[i]){ c->last_error = "points_stream: null path"; return -1; } + paths[i] = image_paths[i]; + } + da::StreamParams sp; + if (chunk_size > 0) sp.chunk_size = chunk_size; + sp.overlap = overlap; // clamped inside stream_points + sp.conf_pct = conf_pct; + if (point_size > 0.f) sp.point_size = point_size; + sp.global_budget = global_budget; + + da::StreamCloud sc; + if (!da::stream_points(*c->engine, paths, c->engine->config(), sp, sc, c->last_error)) return -1; + const size_t np = sc.radius.size(); + if (np == 0){ c->last_error = "points_stream: no points"; return -1; } + + float* xyz = (float*)std::malloc(np * 3 * sizeof(float)); + unsigned char* rgb = (unsigned char*)std::malloc(np * 3); + float* rad = (float*)std::malloc(np * sizeof(float)); + if (!xyz || !rgb || !rad){ std::free(xyz); std::free(rgb); std::free(rad); c->last_error = "points_stream: oom"; return -1; } + std::memcpy(xyz, sc.xyz.data(), np * 3 * sizeof(float)); + std::memcpy(rgb, sc.rgb.data(), np * 3); + std::memcpy(rad, sc.radius.data(), np * sizeof(float)); + if (out_counts) for (int i = 0; i < n_images; ++i) out_counts[i] = (i < (int)sc.counts.size()) ? sc.counts[i] : 0; + if (out_xyz) *out_xyz = xyz; else std::free(xyz); + if (out_rgb) *out_rgb = rgb; else std::free(rgb); + if (out_radius) *out_radius = rad; else std::free(rad); + if (out_n) *out_n = (int)np; + return 0; +} + +int da_capi_gaussians(da_ctx* c, const char* path, int* out_n, + float** out_xyz, float** out_scale, float** out_quat, + float** out_rgb, float** out_opacity){ + if (out_xyz) *out_xyz = nullptr; + if (out_scale) *out_scale = nullptr; + if (out_quat) *out_quat = nullptr; + if (out_rgb) *out_rgb = nullptr; + if (out_opacity) *out_opacity = nullptr; + if (!c || !c->engine || !path){ if (c) c->last_error = "gaussians: bad args"; return -1; } + da::Image img; + if (!da::load_image_rgb(path, img)){ c->last_error = "gaussians: load image failed"; return -1; } + da::Gaussians g; int H = 0, W = 0; + if (!c->engine->reconstruct(img, g, H, W)){ + c->last_error = "gaussians: reconstruct failed (needs a GS model, e.g. DA3-GIANT)"; return -1; } + const int N = g.N; + if (N <= 0 || (int)g.means.size() < N*3 || (int)g.scales.size() < N*3 || + (int)g.rotations.size() < N*4 || (int)g.opacities.size() < N || (int)g.harmonics.size() < N*3*9){ + c->last_error = "gaussians: empty/short arrays"; return -1; } + const double SH_C0 = 0.28209479177387814; + float* xyz = (float*)std::malloc((size_t)N * 3 * sizeof(float)); + float* scl = (float*)std::malloc((size_t)N * 3 * sizeof(float)); + float* quat = (float*)std::malloc((size_t)N * 4 * sizeof(float)); + float* rgb = (float*)std::malloc((size_t)N * 3 * sizeof(float)); + float* op = (float*)std::malloc((size_t)N * sizeof(float)); + if (!xyz || !scl || !quat || !rgb || !op){ + std::free(xyz); std::free(scl); std::free(quat); std::free(rgb); std::free(op); + c->last_error = "gaussians: oom"; return -1; } + for (int i = 0; i < N; ++i){ + for (int k = 0; k < 3; ++k){ xyz[3*i+k] = g.means[3*i+k]; scl[3*i+k] = g.scales[3*i+k]; } + for (int k = 0; k < 4; ++k) quat[4*i+k] = g.rotations[4*i+k]; + for (int ch = 0; ch < 3; ++ch){ + double col = 0.5 + SH_C0 * (double)g.harmonics[((size_t)i*3 + ch)*9 + 0]; + rgb[3*i+ch] = (float)(col < 0 ? 0 : (col > 1 ? 1 : col)); + } + op[i] = g.opacities[i]; + } + if (out_xyz) *out_xyz = xyz; else std::free(xyz); + if (out_scale) *out_scale = scl; else std::free(scl); + if (out_quat) *out_quat = quat; else std::free(quat); + if (out_rgb) *out_rgb = rgb; else std::free(rgb); + if (out_opacity) *out_opacity = op; else std::free(op); + if (out_n) *out_n = N; + return 0; +} } diff --git a/src/dino_backbone.cpp b/src/dino_backbone.cpp index e114f21..76376fa 100644 --- a/src/dino_backbone.cpp +++ b/src/dino_backbone.cpp @@ -241,6 +241,16 @@ bool DinoBackbone::build_feats_graph(ggml_context* ctx, const std::vector ggml_tensor* nw = ml_.tensor("vit.norm.weight"); ggml_tensor* nb = ml_.tensor("vit.norm.bias"); + if (be_.is_offloading()) { + // vit.norm.* are kept as host gguf tensors (no backend buffer) for the + // host-side layernorm used by other paths. Here they feed an IN-GRAPH + // layernorm, so on a GPU backend they must be device-resident: a host-only + // operand has no device subbuffer and crashes the op kernel (on Vulkan, + // ggml_vk_mul -> ggml_vk_tensor_subbuffer null-derefs; CUDA happened to + // tolerate it). Upload them as graph inputs so they get a device buffer. + nw = be_.add_graph_input_nd(ctx, pool, (const float*)nw->data, nw->ne, ggml_n_dims(nw)); + nb = be_.add_graph_input_nd(ctx, pool, (const float*)nb->data, nb->ne, ggml_n_dims(nb)); + } // --- prepare tokens (same graph as forward()) --- const int64_t ine[4]={W,H,3,1}; diff --git a/src/engine.cpp b/src/engine.cpp index 010a6e0..cb9e451 100644 --- a/src/engine.cpp +++ b/src/engine.cpp @@ -260,7 +260,10 @@ bool Engine::depth_pose_multi(const std::vector& imgs, std::vector img_resize_target, round to patch), + // same as the single-image native path. Bounds the input resolution, so a + // full-res photo can't blow the graph up to tens of GB. + if (!preprocess_real(imgs[i], ml_.config(), p)) { DA_LOG("depth_pose_multi: preprocess failed"); return false; } if (i == 0) { H = p.H; W = p.W; } else if (p.H != H || p.W != W) { DA_LOG("depth_pose_multi: views differ in H,W"); return false; } views_chw.push_back(std::move(p.chw)); diff --git a/src/gs_head.cpp b/src/gs_head.cpp index bb558c4..785aa9c 100644 --- a/src/gs_head.cpp +++ b/src/gs_head.cpp @@ -2,6 +2,7 @@ #include "dpt_blocks.hpp" #include "uv_posembed.hpp" #include "ggml_extend.hpp" +#include "common.hpp" #include #include @@ -44,6 +45,14 @@ bool GsHead::raw_gaussians(const std::vector>& feats, if ((int)f.size() != N * C) return false; if ((int)image_chw.size() != 3 * H * W) return false; + // The GSDPT refinenet fusion sizes are fixed for a 16x16 patch grid (a 224x224 + // input). Any other resolution mismatches those hardcoded sizes and would abort + // deep in ggml; fail cleanly so callers can resize to 224x224 and retry. + if (pw != 16 || ph != 16) { + DA_LOG("gs_head: GS reconstruction needs a 224x224 input (16x16 patches); got %dx%d (%dx%d patches)", W, H, pw, ph); + return false; + } + auto t = [&](const std::string& n) { return ml_.tensor(n); }; GraphInputPool pool; diff --git a/src/sim3.cpp b/src/sim3.cpp new file mode 100644 index 0000000..e0eb7e1 --- /dev/null +++ b/src/sim3.cpp @@ -0,0 +1,190 @@ +#include "sim3.hpp" +#include "linalg.hpp" + +#include +#include +#include + +namespace da { + +void sim3_apply(const Sim3& T, const double p[3], double out[3]) { + double x = p[0], y = p[1], z = p[2]; + double rx = T.R[0]*x + T.R[1]*y + T.R[2]*z; + double ry = T.R[3]*x + T.R[4]*y + T.R[5]*z; + double rz = T.R[6]*x + T.R[7]*y + T.R[8]*z; + out[0] = T.s*rx + T.t[0]; + out[1] = T.s*ry + T.t[1]; + out[2] = T.s*rz + T.t[2]; +} + +Sim3 sim3_compose(const Sim3& A, const Sim3& B) { + // (A∘B)(x) = sA RA (sB RB x + tB) + tA + // = (sA sB)(RA RB) x + (sA RA tB + tA) + Sim3 C; + C.s = A.s * B.s; + linalg::mat3_mul(A.R, B.R, C.R); + double rb[3] = { + A.R[0]*B.t[0] + A.R[1]*B.t[1] + A.R[2]*B.t[2], + A.R[3]*B.t[0] + A.R[4]*B.t[1] + A.R[5]*B.t[2], + A.R[6]*B.t[0] + A.R[7]*B.t[1] + A.R[8]*B.t[2], + }; + C.t[0] = A.s*rb[0] + A.t[0]; + C.t[1] = A.s*rb[1] + A.t[1]; + C.t[2] = A.s*rb[2] + A.t[2]; + return C; +} + +// Unit quaternion (w,x,y,z) -> row-major active rotation matrix. +static void quat_to_R(const double q[4], double R[9]) { + double w = q[0], x = q[1], y = q[2], z = q[3]; + R[0] = 1 - 2*(y*y + z*z); R[1] = 2*(x*y - w*z); R[2] = 2*(x*z + w*y); + R[3] = 2*(x*y + w*z); R[4] = 1 - 2*(x*x + z*z); R[5] = 2*(y*z - w*x); + R[6] = 2*(x*z - w*y); R[7] = 2*(y*z + w*x); R[8] = 1 - 2*(x*x + y*y); +} + +// One weighted closed-form solve given raw (un-normalized) weights W[M]. +// Returns false if total weight <= 0. On return: out filled; den = Σ w||a||², +// rot_ok=false if the rotation was ill-posed (kept identity). +static bool solve_once(const double* src, const double* tgt, const double* W, int M, + Sim3& out, double& den_out, bool& rot_ok) { + double Wsum = 0.0; + for (int i = 0; i < M; ++i) Wsum += W[i]; + if (!(Wsum > 0.0) || !std::isfinite(Wsum)) return false; + + // Weighted centroids. + double mu_s[3] = {0,0,0}, mu_t[3] = {0,0,0}; + for (int i = 0; i < M; ++i) { + double w = W[i]; + mu_s[0] += w*src[3*i+0]; mu_s[1] += w*src[3*i+1]; mu_s[2] += w*src[3*i+2]; + mu_t[0] += w*tgt[3*i+0]; mu_t[1] += w*tgt[3*i+1]; mu_t[2] += w*tgt[3*i+2]; + } + for (int k = 0; k < 3; ++k) { mu_s[k] /= Wsum; mu_t[k] /= Wsum; } + + // Cross-covariance Cxx = Σ w a b^T (a=src centered, b=tgt centered) and the + // weighted variances of a (den) and b (num) for the RMS-ratio scale. + double C[9] = {0,0,0, 0,0,0, 0,0,0}; + double den = 0.0; + for (int i = 0; i < M; ++i) { + double w = W[i]; + double a0 = src[3*i+0]-mu_s[0], a1 = src[3*i+1]-mu_s[1], a2 = src[3*i+2]-mu_s[2]; + double b0 = tgt[3*i+0]-mu_t[0], b1 = tgt[3*i+1]-mu_t[1], b2 = tgt[3*i+2]-mu_t[2]; + C[0] += w*a0*b0; C[1] += w*a0*b1; C[2] += w*a0*b2; + C[3] += w*a1*b0; C[4] += w*a1*b1; C[5] += w*a1*b2; + C[6] += w*a2*b0; C[7] += w*a2*b1; C[8] += w*a2*b2; + den += w*(a0*a0 + a1*a1 + a2*a2); + } + den_out = den; + + // Rotation via Horn's quaternion method: largest-eigenvalue eigenvector of the + // symmetric 4x4 N built from C. Reflection-free by construction. + double Sxx=C[0], Sxy=C[1], Sxz=C[2]; + double Syx=C[3], Syy=C[4], Syz=C[5]; + double Szx=C[6], Szy=C[7], Szz=C[8]; + double N[16] = { + Sxx+Syy+Szz, Syz-Szy, Szx-Sxz, Sxy-Syx, + Syz-Szy, Sxx-Syy-Szz, Sxy+Syx, Szx+Sxz, + Szx-Sxz, Sxy+Syx, -Sxx+Syy-Szz, Syz+Szy, + Sxy-Syx, Szx+Sxz, Syz+Szy, -Sxx-Syy+Szz, + }; + double evals[4], evecs[16]; + linalg::jacobi_eigen_sym(N, 4, evals, evecs); + int best = 0, second = -1; + for (int i = 1; i < 4; ++i) if (evals[i] > evals[best]) best = i; + for (int i = 0; i < 4; ++i) { + if (i == best) continue; + if (second < 0 || evals[i] > evals[second]) second = i; + } + rot_ok = true; + Sim3 r; // identity by default + double gap = evals[best] - (second >= 0 ? evals[second] : 0.0); + double mag = std::fabs(evals[best]) + 1e-30; + if (gap < 1e-9 * mag) { + rot_ok = false; // ill-posed (collinear / coincident) -> keep R = I + } else { + double q[4] = { evecs[0*4+best], evecs[1*4+best], evecs[2*4+best], evecs[3*4+best] }; + double qn = std::sqrt(q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]); + if (!(qn > 0.0) || !std::isfinite(qn)) { rot_ok = false; } + else { for (int k = 0; k < 4; ++k) q[k] /= qn; quat_to_R(q, r.R); } + } + + // Optimal Umeyama scale s = tr(R·C) / Σ w||a||². Unlike the RMS-ratio form it + // is linear (not quadratic) in the target points, so gross outliers downweighted + // by Huber can't inflate it; identical to the RMS ratio on clean data. Pinned to + // 1 when the source carries no spread. + double trRC = 0.0; + for (int rr = 0; rr < 3; ++rr) + for (int kk = 0; kk < 3; ++kk) + trRC += r.R[rr*3+kk] * C[kk*3+rr]; + r.s = (den > 1e-30) ? trRC / den : 1.0; + if (!std::isfinite(r.s) || r.s <= 0.0) r.s = 1.0; + + // Translation: t = mu_tgt - s*R*mu_src. + double rs[3]; + rs[0] = r.R[0]*mu_s[0] + r.R[1]*mu_s[1] + r.R[2]*mu_s[2]; + rs[1] = r.R[3]*mu_s[0] + r.R[4]*mu_s[1] + r.R[5]*mu_s[2]; + rs[2] = r.R[6]*mu_s[0] + r.R[7]*mu_s[1] + r.R[8]*mu_s[2]; + r.t[0] = mu_t[0] - r.s*rs[0]; + r.t[1] = mu_t[1] - r.s*rs[1]; + r.t[2] = mu_t[2] - r.s*rs[2]; + out = r; + return true; +} + +bool umeyama_sim3_weighted(const double* src, const double* tgt, const double* w, + int M, Sim3& out, double& rms, int min_pts) { + rms = 0.0; + if (M < min_pts) return false; + + std::vector W(w, w + M); // current IRLS weights + std::vector resid(M, 0.0); + Sim3 T; + double den = 0.0; + bool rot_ok = false; + + const int max_iter = 8; + Sim3 prev; + for (int it = 0; it < max_iter; ++it) { + if (!solve_once(src, tgt, W.data(), M, T, den, rot_ok)) return false; + + // residuals r_i = || s*R*src_i + t - tgt_i || + for (int i = 0; i < M; ++i) { + double p[3] = { src[3*i+0], src[3*i+1], src[3*i+2] }; + double q[3]; sim3_apply(T, p, q); + double dx = q[0]-tgt[3*i+0], dy = q[1]-tgt[3*i+1], dz = q[2]-tgt[3*i+2]; + resid[i] = std::sqrt(dx*dx + dy*dy + dz*dz); + } + // Huber band from the median residual; floor relative to scene extent. + std::vector tmp(resid); + size_t mid = tmp.size()/2; + std::nth_element(tmp.begin(), tmp.begin()+mid, tmp.end()); + double med = tmp[mid]; + double scene = (den > 0.0) ? std::sqrt(den / std::max(1.0, (double)M)) : 1.0; + double delta = std::max(1.345 * med, 1e-6 * scene + 1e-30); + // Reweight: conf * huber(r,delta). + for (int i = 0; i < M; ++i) { + double hu = (resid[i] <= delta) ? 1.0 : (delta / resid[i]); + W[i] = w[i] * hu; + } + // Convergence: small change in (s,R,t). + if (it > 0) { + double d = std::fabs(T.s - prev.s); + for (int k = 0; k < 9; ++k) d += std::fabs(T.R[k] - prev.R[k]); + for (int k = 0; k < 3; ++k) d += std::fabs(T.t[k] - prev.t[k]); + if (d < 1e-12) break; + } + prev = T; + } + + // Final weighted RMS with the converged (conf*huber) weights. + double sw = 0.0, se = 0.0; + for (int i = 0; i < M; ++i) { sw += W[i]; se += W[i]*resid[i]*resid[i]; } + rms = (sw > 0.0) ? std::sqrt(se / sw) : 0.0; + + for (int k = 0; k < 9; ++k) if (!std::isfinite(T.R[k])) return false; + for (int k = 0; k < 3; ++k) if (!std::isfinite(T.t[k])) return false; + if (!std::isfinite(T.s)) return false; + out = T; + return true; +} + +} // namespace da diff --git a/src/sim3.hpp b/src/sim3.hpp new file mode 100644 index 0000000..1d071c1 --- /dev/null +++ b/src/sim3.hpp @@ -0,0 +1,32 @@ +#pragma once +// Weighted Umeyama Sim3 alignment (NO ggml, NO engine). Used by the sliding-window +// streaming stitcher to map each window's local point cloud into a single global +// frame. Mirrors da3_streaming/loop_utils/alignment_torch.py: weighted centroids, +// rotation from a cross-covariance, RMS-ratio scale, translation, Huber IRLS. +namespace da { + +// Similarity transform: y = s * R * x + t. R is row-major 3x3. +struct Sim3 { + double s = 1.0; + double R[9] = {1,0,0, 0,1,0, 0,0,1}; + double t[3] = {0,0,0}; +}; + +// out = s*R*p + t. Safe when out aliases p. +void sim3_apply(const Sim3& T, const double p[3], double out[3]); + +// Returns A∘B with (A∘B)(x) = A(B(x)). +Sim3 sim3_compose(const Sim3& A, const Sim3& B); + +// Weighted Umeyama similarity fitting src -> tgt, refined by Huber IRLS. +// src,tgt : length 3*M (x,y,z interleaved). +// w : length M, per-correspondence weight (>=0; e.g. min confidence). +// On success fills `out` (maps src->tgt) and `rms` (weighted RMS residual) and +// returns true. Returns false if degenerate: M < min_pts, all weights ~0, or a +// non-finite result. When the source points carry no scale information (nearly +// coincident) the scale is pinned to 1; when the rotation is ill-posed (collinear/ +// coincident) R is left at identity (translation+scale only). +bool umeyama_sim3_weighted(const double* src, const double* tgt, const double* w, + int M, Sim3& out, double& rms, int min_pts = 50); + +} // namespace da diff --git a/src/stream.cpp b/src/stream.cpp new file mode 100644 index 0000000..069517f --- /dev/null +++ b/src/stream.cpp @@ -0,0 +1,242 @@ +#include "stream.hpp" + +#include "engine.hpp" +#include "preprocess.hpp" +#include "image_io.hpp" +#include "reconstruct.hpp" +#include "sim3.hpp" +#include "common.hpp" + +#include +#include +#include + +namespace da { + +// Dense per-pixel world points (in a window's LOCAL frame) for a handful of +// overlap frames, kept around to align the NEXT window against. +struct OverlapCache { + int H = 0, W = 0, n = 0; // n cached frames (== overlap of a full window) + std::vector valid; // n*H*W : finite,d>0,conf>=thr at snapshot time + std::vector xyz; // n*H*W*3 : LOCAL-frame world points + std::vector conf; // n*H*W + Sim3 G; // LOCAL -> GLOBAL for the cached window +}; + +// Back-project one frame's pixels into LOCAL world space, filling dense valid/xyz/conf. +// Validity is GEOMETRIC only (finite, d>0) — alignment uses conf as a weight, not a +// hard gate, so the overlap correspondence isn't starved at low-confidence window +// edges. conf is stored raw (or 1.0 when the model gives no conf) for that weighting. +static void backproject_frame_dense(const ViewResult& view, + const std::array& K, + const std::array& E4, + int H, int W, bool have_conf, + uint8_t* valid, double* xyz, float* conf) { + std::array Kinv; std::array c2w; + const bool kok = inv3(K, Kinv); + const bool eok = inv4(E4, c2w); + double ki[9]; for (int k=0;k<9;++k) ki[k] = kok ? (double)Kinv[k] : 0.0; + double cw[16]; for (int k=0;k<16;++k) cw[k] = eok ? (double)c2w[k] : 0.0; + const size_t plane = (size_t)H*(size_t)W; + const float* dF = view.depth.data(); + const float* cF = (have_conf && view.conf.size()==plane) ? view.conf.data() : nullptr; + for (int v = 0; v < H; ++v) { + for (int u = 0; u < W; ++u) { + size_t pix = (size_t)v*W + u; + float d = dF[pix]; + float cc = cF ? cF[pix] : 1.0f; + conf[pix] = cc; + bool ok = kok && eok && std::isfinite(d) && d > 0.0f; + if (!ok) { valid[pix] = 0; xyz[3*pix+0]=xyz[3*pix+1]=xyz[3*pix+2]=0.0; continue; } + double rx = ki[0]*u + ki[1]*v + ki[2]; + double ry = ki[3]*u + ki[4]*v + ki[5]; + double rz = ki[6]*u + ki[7]*v + ki[8]; + double xc = rx*d, yc = ry*d, zc = rz*d; + xyz[3*pix+0] = cw[0]*xc + cw[1]*yc + cw[2]*zc + cw[3]; + xyz[3*pix+1] = cw[4]*xc + cw[5]*yc + cw[6]*zc + cw[7]; + xyz[3*pix+2] = cw[8]*xc + cw[9]*yc + cw[10]*zc + cw[11]; + valid[pix] = 1; + } + } +} + +bool stream_points(Engine& eng, const std::vector& frame_paths, + const Config& cfg, const StreamParams& p, + StreamCloud& out, std::string& err) { + const int F = (int)frame_paths.size(); + if (F == 0) { err = "stream: no frames"; return false; } + + int chunk = p.chunk_size < 2 ? 2 : p.chunk_size; + int overlap = p.overlap; + if (overlap < 0) overlap = 0; + if (overlap > chunk - 1) overlap = chunk - 1; + int step = chunk - overlap; + if (step < 1) step = 1; + double conf_pct = p.conf_pct; if (conf_pct < 0) conf_pct = 0; if (conf_pct > 100) conf_pct = 100; + float point_size = (p.point_size > 0.f) ? p.point_size : 1.f; + + out.xyz.clear(); out.rgb.clear(); out.radius.clear(); + out.counts.assign(F, 0); + out.warnings = 0; + + Sim3 G; // current window LOCAL -> GLOBAL + OverlapCache prev; // previous window's tail-overlap geometry + bool have_prev = false; + + for (int w0 = 0; w0 < F; w0 += step) { + const int w1 = std::min(w0 + chunk, F); + const int nv = w1 - w0; + + // --- load + fused inference for this window --- + std::vector imgs(nv); + for (int i = 0; i < nv; ++i) + if (!load_image_rgb(frame_paths[w0+i], imgs[i])) { + err = "stream: load image failed: " + frame_paths[w0+i]; return false; } + + std::vector views; int H = 0, W = 0; + if (!eng.depth_pose_multi(imgs, views, H, W) || (int)views.size() != nv) { + err = "stream: depth_pose_multi failed"; return false; } + const size_t plane = (size_t)H * (size_t)W; + + // K, E4, processed RGB per view. + std::vector> K(nv); + std::vector> E4(nv); + std::vector> rgb_store(nv); + std::vector images_u8(nv); + bool have_conf = true; + std::vector conf_all; conf_all.reserve((size_t)nv*plane); + for (int i = 0; i < nv; ++i) { + K[i] = views[i].intr; + std::array e4{}; for (int k=0;k<12;++k) e4[k]=views[i].ext[k]; + e4[12]=0.f; e4[13]=0.f; e4[14]=0.f; e4[15]=1.f; E4[i]=e4; + Preprocessed pp; + if (!preprocess_real(imgs[i], cfg, pp, &rgb_store[i]) || pp.H!=H || pp.W!=W) { + err = "stream: preprocess color mismatch"; return false; } + images_u8[i] = rgb_store[i].data(); + if (views[i].conf.size() == plane) conf_all.insert(conf_all.end(), views[i].conf.begin(), views[i].conf.end()); + else have_conf = false; + } + if (!have_conf) conf_all.clear(); + float conf_thr = -1e30f; + if (have_conf && !conf_all.empty()) conf_thr = (float)percentile_linear(conf_all, conf_pct); + + // --- stitch: solve G (LOCAL -> GLOBAL) for this window --- + if (!have_prev) { + G = Sim3(); // window 0 defines the global frame + } else { + const int ov = std::min(prev.n, nv); + std::vector srcv, tgtv, wv; // cur-local -> prev-local + srcv.reserve((size_t)ov*plane*3); + for (int o = 0; o < ov; ++o) { + std::array Kinv; std::array c2w; + if (!inv3(K[o], Kinv) || !inv4(E4[o], c2w)) continue; + double ki[9]; for (int k=0;k<9;++k) ki[k]=(double)Kinv[k]; + double cw[16]; for (int k=0;k<16;++k) cw[k]=(double)c2w[k]; + const float* dcur = views[o].depth.data(); + const float* ccur = (have_conf && views[o].conf.size()==plane) ? views[o].conf.data() : nullptr; + const uint8_t* pvalid = prev.valid.data() + (size_t)o*plane; + const double* pxyz = prev.xyz.data() + (size_t)o*plane*3; + const float* pconf = prev.conf.data() + (size_t)o*plane; + for (int v = 0; v < H; ++v) for (int u = 0; u < W; ++u) { + size_t pix = (size_t)v*W + u; + if (!pvalid[pix]) continue; + float d = dcur[pix]; + if (!std::isfinite(d) || d <= 0.0f) continue; + float cc = ccur ? ccur[pix] : 1.0f; // weight only, no hard gate + double rx = ki[0]*u + ki[1]*v + ki[2]; + double ry = ki[3]*u + ki[4]*v + ki[5]; + double rz = ki[6]*u + ki[7]*v + ki[8]; + double xc = rx*d, yc = ry*d, zc = rz*d; + srcv.push_back(cw[0]*xc + cw[1]*yc + cw[2]*zc + cw[3]); + srcv.push_back(cw[4]*xc + cw[5]*yc + cw[6]*zc + cw[7]); + srcv.push_back(cw[8]*xc + cw[9]*yc + cw[10]*zc + cw[11]); + tgtv.push_back(pxyz[3*pix+0]); tgtv.push_back(pxyz[3*pix+1]); tgtv.push_back(pxyz[3*pix+2]); + double wgt = std::min((double)cc, (double)pconf[pix]); + wv.push_back(wgt > 0.0 ? wgt : 1e-6); + } + } + const int M = (int)wv.size(); + Sim3 S_rel; double rms = 0.0; + bool ok = M >= p.min_overlap_pts && + umeyama_sim3_weighted(srcv.data(), tgtv.data(), wv.data(), M, S_rel, rms, p.min_overlap_pts); + if (!ok) { + out.warnings++; + DA_LOG("stream: window @%d degenerate stitch (M=%d) -> identity seam", w0, M); + S_rel = Sim3(); + } else { + DA_LOG("stream: window @%d stitch M=%d s=%.4f rms=%.4g", w0, M, S_rel.s, rms); + } + G = sim3_compose(prev.G, S_rel); // cur-local -> global + } + + // --- emit NEW frames only --- + const int new_lo = have_prev ? std::min(prev.n, nv) : 0; + const int newN = nv - new_lo; + if (newN > 0) { + std::vector depth_slice; depth_slice.reserve((size_t)newN*plane); + std::vector conf_slice; if (have_conf) conf_slice.reserve((size_t)newN*plane); + std::vector> Ks(newN); + std::vector> Es(newN); + std::vector Iu8(newN); + for (int i = new_lo; i < nv; ++i) { + int j = i - new_lo; + depth_slice.insert(depth_slice.end(), views[i].depth.begin(), views[i].depth.end()); + if (have_conf) conf_slice.insert(conf_slice.end(), views[i].conf.begin(), views[i].conf.end()); + Ks[j] = K[i]; Es[j] = E4[i]; Iu8[j] = images_u8[i]; + } + WorldPoints wp = back_project(depth_slice, have_conf ? conf_slice : std::vector{}, + Ks, Es, Iu8, H, W, newN, conf_thr); + const size_t np = wp.frame.size(); + + // Per-window budget: uniform stride over the frame-major point list. + size_t stride = 1; + if (p.global_budget > 0 && F > 0) { + double quota = (double)p.global_budget * (double)newN / (double)F; + if (quota < 1.0) quota = 1.0; + if ((double)np > quota) stride = (size_t)std::ceil((double)np / quota); + if (stride < 1) stride = 1; + } + + for (size_t k = 0; k < np; k += stride) { + int f = wp.frame[k], u = wp.u[k], v = wp.v[k]; + double pL[3] = { wp.xyz[3*k+0], wp.xyz[3*k+1], wp.xyz[3*k+2] }; + double pG[3]; sim3_apply(G, pL, pG); + out.xyz.push_back((float)pG[0]); + out.xyz.push_back((float)pG[1]); + out.xyz.push_back((float)pG[2]); + out.rgb.push_back(wp.rgb[3*k+0]); + out.rgb.push_back(wp.rgb[3*k+1]); + out.rgb.push_back(wp.rgb[3*k+2]); + float d = depth_slice[(size_t)f*plane + (size_t)v*W + u]; + float fx = Ks[f][0], fy = Ks[f][4]; + float rr = 0.5f * (d/fx + d/fy) * point_size * (float)G.s; + if (!(rr > 0.f) || !std::isfinite(rr)) rr = 1e-4f; + out.radius.push_back(rr); + int gframe = w0 + new_lo + f; + if (gframe >= 0 && gframe < F) out.counts[gframe]++; + } + } + + // --- snapshot this window's TAIL overlap (LOCAL frame) for the next stitch --- + const int ovkeep = std::min(overlap, nv); + prev.H = H; prev.W = W; prev.n = ovkeep; prev.G = G; + prev.valid.assign((size_t)ovkeep*plane, 0); + prev.xyz.assign((size_t)ovkeep*plane*3, 0.0); + prev.conf.assign((size_t)ovkeep*plane, 0.f); + for (int o = 0; o < ovkeep; ++o) { + int vi = nv - ovkeep + o; // local index of this tail frame + backproject_frame_dense(views[vi], K[vi], E4[vi], H, W, have_conf, + prev.valid.data() + (size_t)o*plane, + prev.xyz.data() + (size_t)o*plane*3, + prev.conf.data() + (size_t)o*plane); + } + have_prev = true; + + if (w1 >= F) break; // last (possibly truncated) window processed + } + + if (out.radius.empty()) { err = "stream: no points survived (raise conf_pct or check parallax)"; return false; } + return true; +} + +} // namespace da diff --git a/src/stream.hpp b/src/stream.hpp new file mode 100644 index 0000000..6a1e406 --- /dev/null +++ b/src/stream.hpp @@ -0,0 +1,38 @@ +#pragma once +// Sliding-window streaming reconstruction: turn a long ordered frame list into ONE +// coherent global colored point cloud, mirroring upstream da3_streaming. Each +// overlapping window is run through the existing fused multi-view pass +// (Engine::depth_pose_multi); consecutive windows are stitched into a single global +// frame with a weighted-Umeyama Sim3 solved on the overlap frames. No ggml/model +// changes — pure host orchestration on top of inference + back_project. +#include "model_loader.hpp" // da::Config +#include +#include +#include + +namespace da { +class Engine; + +struct StreamParams { + int chunk_size = 12; // frames per fused window (<= ~24) + int overlap = 3; // shared frames between consecutive windows + double conf_pct = 55.0; // per-window confidence percentile gate + float point_size = 1.2f; // per-point radius multiplier + int global_budget = 6'000'000; // <=0 => unlimited; total emitted point cap + int min_overlap_pts = 50; // degenerate guard for the Sim3 solve +}; + +struct StreamCloud { + std::vector xyz; // 3N, GLOBAL frame (OpenCV axes) + std::vector rgb; // 3N + std::vector radius; // N + std::vector counts; // per INPUT frame (length = #frames); build-up prefix sums + int warnings = 0; // # windows that fell back to a degraded seam +}; + +// Stitch frame_paths (in order) into out. Returns false + sets err on failure. +bool stream_points(Engine& eng, const std::vector& frame_paths, + const Config& cfg, const StreamParams& p, + StreamCloud& out, std::string& err); + +} // namespace da diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index fd1b5f5..baae6d4 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -37,6 +37,8 @@ da_add_test(test_gs_head) da_add_test(test_gs_adapter) da_add_test(test_winograd) da_add_test(test_linalg) +da_add_test(test_sim3) +da_add_test(test_stream_seam) da_add_test(test_ray_pose) da_add_test(test_dpt_blocks) da_add_test(test_dpt_head) diff --git a/tests/test_capi.cpp b/tests/test_capi.cpp index b37dbf0..b7fb3df 100644 --- a/tests/test_capi.cpp +++ b/tests/test_capi.cpp @@ -5,7 +5,7 @@ int main(){ const char* gguf = std::getenv("DA_TEST_GGUF"); if (!gguf) return 77; - if (da_capi_abi_version() != 4) return 1; + if (da_capi_abi_version() != 6) return 1; da_ctx* c = da_capi_load(gguf, 1); if (!c) { std::fprintf(stderr, "load failed\n"); return 1; } char* j = da_capi_info_json(c); diff --git a/tests/test_sim3.cpp b/tests/test_sim3.cpp new file mode 100644 index 0000000..e613265 --- /dev/null +++ b/tests/test_sim3.cpp @@ -0,0 +1,119 @@ +// Host-only gate for the weighted-Umeyama Sim3 solver (src/sim3.cpp): exact +// recovery, the reflection trap (proper rotation only), Huber robustness to gross +// outliers, and compose associativity. No GGUF needed. +#include "sim3.hpp" +#include "linalg.hpp" +#include +#include +#include + +using namespace da; + +static unsigned long long g_seed = 88172645463325252ULL; +static double urand() { // xorshift64 -> [0,1) + g_seed ^= g_seed << 13; g_seed ^= g_seed >> 7; g_seed ^= g_seed << 17; + return (double)(g_seed >> 11) / 9007199254740992.0; +} +static double sym() { return urand() * 2.0 - 1.0; } + +static void axis_angle_R(double ax, double ay, double az, double ang, double R[9]) { + double n = std::sqrt(ax*ax + ay*ay + az*az); + if (n < 1e-12) { R[0]=R[4]=R[8]=1; R[1]=R[2]=R[3]=R[5]=R[6]=R[7]=0; return; } + ax/=n; ay/=n; az/=n; + double c = std::cos(ang), s = std::sin(ang), C = 1 - c; + R[0]=c+ax*ax*C; R[1]=ax*ay*C-az*s; R[2]=ax*az*C+ay*s; + R[3]=ay*ax*C+az*s; R[4]=c+ay*ay*C; R[5]=ay*az*C-ax*s; + R[6]=az*ax*C-ay*s; R[7]=az*ay*C+ax*s; R[8]=c+az*az*C; +} +static double rdiff(const double A[9], const double B[9]) { + double s = 0; for (int i = 0; i < 9; ++i) { double d = A[i]-B[i]; s += d*d; } return std::sqrt(s); +} +static double tdiff(const double a[3], const double b[3]) { + double s = 0; for (int i = 0; i < 3; ++i) { double d = a[i]-b[i]; s += d*d; } return std::sqrt(s); +} + +static int fails = 0; +static void check(bool ok, const char* msg) { + if (!ok) { std::fprintf(stderr, "FAIL: %s\n", msg); fails++; } + else std::fprintf(stderr, "ok: %s\n", msg); +} + +// Build M correspondences tgt = GT(src) for random src in [-5,5]^3. +static void make(int M, const Sim3& GT, std::vector& src, std::vector& tgt) { + src.resize(3*M); tgt.resize(3*M); + for (int i = 0; i < M; ++i) { + double p[3] = { sym()*5, sym()*5, sym()*5 }; + double q[3]; sim3_apply(GT, p, q); + for (int k = 0; k < 3; ++k) { src[3*i+k] = p[k]; tgt[3*i+k] = q[k]; } + } +} + +int main() { + // (a) exact recovery. + { + Sim3 GT; GT.s = 1.7; axis_angle_R(0.3,-0.7,0.5, 0.9, GT.R); + GT.t[0]=2.0; GT.t[1]=-1.0; GT.t[2]=0.5; + int M = 300; std::vector src, tgt; make(M, GT, src, tgt); + std::vector w(M, 1.0); + Sim3 out; double rms = 0; + bool ok = umeyama_sim3_weighted(src.data(), tgt.data(), w.data(), M, out, rms); + check(ok, "exact: solver returned true"); + check(rms < 1e-6, "exact: rms ~ 0"); + check(std::fabs(out.s - GT.s) < 1e-6, "exact: scale"); + check(rdiff(out.R, GT.R) < 1e-6, "exact: rotation"); + check(tdiff(out.t, GT.t) < 1e-6, "exact: translation"); + } + // (b) reflection trap: R = diag(-1,-1,1) (180 deg about z) must stay proper. + { + Sim3 GT; GT.s = 1.0; + GT.R[0]=-1; GT.R[4]=-1; GT.R[8]=1; + int M = 250; std::vector src, tgt; make(M, GT, src, tgt); + std::vector w(M, 1.0); + Sim3 out; double rms = 0; + bool ok = umeyama_sim3_weighted(src.data(), tgt.data(), w.data(), M, out, rms); + check(ok, "reflect: solver returned true"); + check(std::fabs(linalg::det3(out.R) - 1.0) < 1e-6, "reflect: det(R) == +1 (proper)"); + check(rdiff(out.R, GT.R) < 1e-5, "reflect: recovered diag(-1,-1,1)"); + } + // (c) Huber robustness: 20% gross outliers, all weights 1. + { + Sim3 GT; GT.s = 1.3; axis_angle_R(-0.2,0.5,0.8, 1.4, GT.R); + GT.t[0]=-3.0; GT.t[1]=0.7; GT.t[2]=4.0; + int M = 400; std::vector src, tgt; make(M, GT, src, tgt); + std::vector w(M, 1.0); + for (int i = 0; i < M; i += 5) { // ~20% corrupted + tgt[3*i+0] += sym()*30; tgt[3*i+1] += sym()*30; tgt[3*i+2] += sym()*30; + } + Sim3 out; double rms = 0; + bool ok = umeyama_sim3_weighted(src.data(), tgt.data(), w.data(), M, out, rms); + check(ok, "huber: solver returned true"); + check(std::fabs(out.s - GT.s) < 5e-2, "huber: scale within 5e-2"); + check(rdiff(out.R, GT.R) < 5e-2, "huber: rotation within 5e-2"); + check(tdiff(out.t, GT.t) < 0.25, "huber: translation within 0.25"); + } + // (d) compose associativity: (A∘B)(x) == A(B(x)). + { + Sim3 A; A.s = 0.8; axis_angle_R(0.1,0.2,0.9, 0.6, A.R); A.t[0]=1; A.t[1]=2; A.t[2]=-1; + Sim3 B; B.s = 1.5; axis_angle_R(0.7,-0.3,0.2, -1.1, B.R); B.t[0]=-2; B.t[1]=0.5; B.t[2]=3; + Sim3 C = sim3_compose(A, B); + double maxd = 0; + for (int i = 0; i < 50; ++i) { + double x[3] = { sym()*4, sym()*4, sym()*4 }; + double bx[3], abx[3], cx[3]; + sim3_apply(B, x, bx); sim3_apply(A, bx, abx); sim3_apply(C, x, cx); + maxd = std::max(maxd, tdiff(abx, cx)); + } + check(maxd < 1e-9, "compose: (A.B)(x) == A(B(x))"); + } + // (e) degenerate: too few points -> false. + { + Sim3 GT; int M = 10; std::vector src, tgt; make(M, GT, src, tgt); + std::vector w(M, 1.0); + Sim3 out; double rms = 0; + bool ok = umeyama_sim3_weighted(src.data(), tgt.data(), w.data(), M, out, rms); + check(!ok, "degenerate: M +#include +#include + +using namespace da; + +static unsigned long long g_seed = 0x9e3779b97f4a7c15ULL; +static double urand() { + g_seed ^= g_seed << 13; g_seed ^= g_seed >> 7; g_seed ^= g_seed << 17; + return (double)(g_seed >> 11) / 9007199254740992.0; +} +static double sym() { return urand() * 2.0 - 1.0; } + +static void axis_angle_R(double ax, double ay, double az, double ang, double R[9]) { + double n = std::sqrt(ax*ax + ay*ay + az*az); ax/=n; ay/=n; az/=n; + double c = std::cos(ang), s = std::sin(ang), C = 1 - c; + R[0]=c+ax*ax*C; R[1]=ax*ay*C-az*s; R[2]=ax*az*C+ay*s; + R[3]=ay*ax*C+az*s; R[4]=c+ay*ay*C; R[5]=ay*az*C-ax*s; + R[6]=az*ax*C-ay*s; R[7]=az*ay*C+ax*s; R[8]=c+az*az*C; +} +static Sim3 sim3_inv(const Sim3& T) { + Sim3 I; I.s = 1.0 / T.s; + I.R[0]=T.R[0]; I.R[1]=T.R[3]; I.R[2]=T.R[6]; + I.R[3]=T.R[1]; I.R[4]=T.R[4]; I.R[5]=T.R[7]; + I.R[6]=T.R[2]; I.R[7]=T.R[5]; I.R[8]=T.R[8]; + double rt[3] = { I.R[0]*T.t[0]+I.R[1]*T.t[1]+I.R[2]*T.t[2], + I.R[3]*T.t[0]+I.R[4]*T.t[1]+I.R[5]*T.t[2], + I.R[6]*T.t[0]+I.R[7]*T.t[1]+I.R[8]*T.t[2] }; + I.t[0]=-I.s*rt[0]; I.t[1]=-I.s*rt[1]; I.t[2]=-I.s*rt[2]; + return I; +} + +static int fails = 0; +static void check(bool ok, const char* msg) { + if (!ok) { std::fprintf(stderr, "FAIL: %s\n", msg); fails++; } + else std::fprintf(stderr, "ok: %s\n", msg); +} + +// Run one seam: prev-window local->global = Gp, cur-window local->global = Sc. +// noise = per-coordinate jitter added to local points. Returns max global residual. +static double run_seam(const Sim3& Gp, const Sim3& Sc, double noise, double& srel_s) { + const Sim3 Gp_inv = sim3_inv(Gp), Sc_inv = sim3_inv(Sc); + const int Ngrid = 3000; + std::vector src, tgt, w; // cur-local -> prev-local + std::vector keep_cur, keep_glob; // both-valid cur-local & true global + for (int i = 0; i < Ngrid; ++i) { + double Pg[3] = { sym()*5, sym()*5, sym()*5 }; + double pl[3], cl[3]; + sim3_apply(Gp_inv, Pg, pl); // prev-local + sim3_apply(Sc_inv, Pg, cl); // cur-local + for (int k = 0; k < 3; ++k) { pl[k] += sym()*noise; cl[k] += sym()*noise; } + bool pv = urand() < 0.7, cv = urand() < 0.7; // independent disjoint masks + if (!(pv && cv)) continue; + double cw = 0.5 + 0.5*urand(), pw = 0.5 + 0.5*urand(); + src.push_back(cl[0]); src.push_back(cl[1]); src.push_back(cl[2]); + tgt.push_back(pl[0]); tgt.push_back(pl[1]); tgt.push_back(pl[2]); + w.push_back(cw < pw ? cw : pw); + keep_cur.push_back(cl[0]); keep_cur.push_back(cl[1]); keep_cur.push_back(cl[2]); + keep_glob.push_back(Pg[0]); keep_glob.push_back(Pg[1]); keep_glob.push_back(Pg[2]); + } + int M = (int)w.size(); + Sim3 S_rel; double rms = 0; + bool ok = umeyama_sim3_weighted(src.data(), tgt.data(), w.data(), M, S_rel, rms); + if (!ok) { srel_s = 0; return 1e30; } + srel_s = S_rel.s; + Sim3 G_cur = sim3_compose(Gp, S_rel); // cur-local -> global + double maxd = 0; + for (size_t i = 0; i < keep_cur.size()/3; ++i) { + double cl[3] = { keep_cur[3*i+0], keep_cur[3*i+1], keep_cur[3*i+2] }; + double g[3]; sim3_apply(G_cur, cl, g); + double dx = g[0]-keep_glob[3*i+0], dy = g[1]-keep_glob[3*i+1], dz = g[2]-keep_glob[3*i+2]; + double d = std::sqrt(dx*dx + dy*dy + dz*dz); + if (d > maxd) maxd = d; + } + std::fprintf(stderr, " seam: M=%d s_rel=%.4f rms=%.3g max_resid=%.3g\n", M, S_rel.s, rms, maxd); + return maxd; +} + +int main() { + Sim3 Gp; Gp.s = 1.2; axis_angle_R(0.2,-0.4,0.6, 0.5, Gp.R); Gp.t[0]=1; Gp.t[1]=-2; Gp.t[2]=0.5; + Sim3 Sc; Sc.s = 0.9; axis_angle_R(-0.5,0.3,0.7, -0.8, Sc.R); Sc.t[0]=-1.5; Sc.t[1]=0.4; Sc.t[2]=2.0; + + // Exact: disjoint masks, no noise -> overlap realigns to machine precision. + double srel_s = 0; + double m0 = run_seam(Gp, Sc, 0.0, srel_s); + check(m0 < 1e-5, "seam exact: overlap realigns (max resid < 1e-5)"); + // S_rel scale should equal Sc.s / Gp.s (cur-local -> prev-local). + check(std::fabs(srel_s - Sc.s/Gp.s) < 1e-4, "seam exact: relative scale = Sc.s/Gp.s"); + + // Noisy: small jitter on local points -> still tightly aligned. + double m1 = run_seam(Gp, Sc, 0.01, srel_s); + check(m1 < 0.2, "seam noisy: overlap realigns within 0.2"); + + std::fprintf(stderr, "%s (%d failures)\n", fails ? "FAILED" : "PASSED", fails); + return fails ? 1 : 0; +} diff --git a/third_party/ggml b/third_party/ggml index 3af5f57..eced84c 160000 --- a/third_party/ggml +++ b/third_party/ggml @@ -1 +1 @@ -Subproject commit 3af5f5760e19a96427f5f7a93b79cbdf3d4b265b +Subproject commit eced84c86f8b012c752c016f7fe789adea168e1e From b4c0f74f8dc536f6fd2f9ba3398fb90f56d9d29d Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Wed, 1 Jul 2026 06:18:16 +0100 Subject: [PATCH 02/22] fix(backend): enlarge graph-node pool for giant multi-view windows MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A full multi-view streaming window on GIANT (40 layers) overflowed the fixed ggml metadata context: 12 views needed just over kGraphSize=16384 nodes and aborted with "not enough space in the context's memory pool" mid-backbone. Base multi-view and single-image giant fit; giant x 12-view did not. Raise kGraphSize to 49152 — sized for the heaviest window the server allows (24-view giant, ~33k nodes) with margin. The pool is only tensor metadata (~408 B/node, ~20 MB), allocated per compute() call, so the cost is trivial. Fixes the server SIGABRT when baking a giant-model video via streaming. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/backend.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/backend.cpp b/src/backend.cpp index 93bc5bf..2f213e4 100644 --- a/src/backend.cpp +++ b/src/backend.cpp @@ -14,10 +14,12 @@ namespace da { namespace { -// Number of graph nodes the metadata context must hold. Ample for the MoonViT -// tower; the per-compute metadata context + graph hash-set scale with this, so -// keep it as small as the largest single graph needs. -constexpr size_t kGraphSize = 16384; +// Number of graph nodes the metadata context must hold. Sized for the heaviest +// supported graph: a full multi-view streaming window on GIANT (40 layers) at the +// max chunk of 24 views (~33k nodes). 16384 fit base multi-view and single-image +// giant but overflowed giant x 12-view by a hair. The per-compute metadata context +// + graph hash-set scale with this (~408 B/node, ~20 MB here) — cheap. +constexpr size_t kGraphSize = 49152; struct PendingInput { ggml_tensor* tensor; From ef4a159ff79e48bd87a16a480cdec6833942604a Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Wed, 1 Jul 2026 06:28:51 +0100 Subject: [PATCH 03/22] fix(ui): source toggle stuck on first file type after a pick MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The drop box held a whose text setSrc() updated, but reflectChosen() overwrites the whole box via drop.textContent when a file is chosen — destroying that span. The next source switch then hit $("droptype").textContent on a null element and threw, aborting setSrc() before it updated the file input's accept filter (and multiple / vidparams). So after picking one file, toggling Video/Photos only moved the highlight and the picker kept filtering to the first type. Drop the span; setSrc() no longer touches it. Set accept/multiple fresh in drop.onclick right before opening the picker (and clear value so re-picking the same file re-fires onchange). A shared dropText() restores the prompt. Co-Authored-By: Claude Opus 4.8 (1M context) --- server/web/index.html | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/server/web/index.html b/server/web/index.html index 63dfc82..40ae5cc 100644 --- a/server/web/index.html +++ b/server/web/index.html @@ -105,7 +105,7 @@

depth-anything.cpp · DA3 studio

-
click or drop a video here
+
click or drop a video here
@@ -324,20 +324,24 @@

depth-anything.cpp · DA3 studio

let mode="points", src="video"; function setMode(m){ mode=m; document.querySelectorAll('#modeseg button').forEach(b=>b.classList.toggle("on",b.dataset.mode===m)); } function setSrc(sv){ src=sv; document.querySelectorAll('#srcseg button').forEach(b=>b.classList.toggle("on",b.dataset.src===sv)); - $("droptype").textContent=sv==="video"?"video":"photo(s)"; $("file").toggleAttribute("multiple",sv==="photos"); - $("file").setAttribute("accept",sv==="video"?"video/*":"image/*"); $("vidparams").classList.toggle("hide",sv!=="video"); chosen=null; refreshGen(); } + $("vidparams").classList.toggle("hide",sv!=="video"); chosen=null; dropText(); refreshGen(); } document.querySelectorAll('#modeseg button').forEach(b=>b.onclick=()=>{ if(!b.disabled) setMode(b.dataset.mode); }); document.querySelectorAll('#srcseg button').forEach(b=>b.onclick=()=>setSrc(b.dataset.src)); // ---- file pick / drop ---- let chosen=null; const drop=$("drop"), fileInput=$("file"); -drop.onclick=()=>fileInput.click(); +// Set the picker filter fresh on every open (a reused otherwise keeps a +// stale accept), and clear value so re-picking the same file still fires onchange. +drop.onclick=()=>{ fileInput.value=""; fileInput.toggleAttribute("multiple",src==="photos"); + fileInput.setAttribute("accept",src==="video"?"video/*":"image/*"); fileInput.click(); }; fileInput.onchange=()=>{ chosen=fileInput.files; reflectChosen(); }; ["dragover","dragenter"].forEach(ev=>drop.addEventListener(ev,e=>{e.preventDefault();drop.classList.add("hl");})); ["dragleave","drop"].forEach(ev=>drop.addEventListener(ev,e=>{e.preventDefault();drop.classList.remove("hl");})); drop.addEventListener("drop",e=>{ chosen=e.dataTransfer.files; reflectChosen(); }); -function reflectChosen(){ if(chosen&&chosen.length){ drop.textContent=chosen.length>1?`${chosen.length} files`:chosen[0].name; } refreshGen(); } +function dropText(){ drop.textContent=(chosen&&chosen.length)?(chosen.length>1?`${chosen.length} files`:chosen[0].name) + :`click or drop a ${src==="video"?"video":"photo(s)"} here`; } +function reflectChosen(){ dropText(); refreshGen(); } function refreshGen(){ $("gen").disabled=!(chosen&&chosen.length); } // ---- generate ---- From f4f4420610d3d490fe4002ae0bab3b75c4b458e7 Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Wed, 1 Jul 2026 07:57:55 +0100 Subject: [PATCH 04/22] feat(ui): present single-image gaussians from the input view MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit DA3's single-image gaussians are view-synthesis primitives: their per-pixel depth is intentionally spread (verified bit-exact against the DA3 reference via the regenerated giant parity dumps — test_gs_head/test_gs_adapter pass at 2e-3, and the reference's own gaussians span the same 100x+ depth range at full res). They only composite into a coherent image from the INPUT camera. The viewer was free-orbiting them and auto-framing to their huge bounding box, which collapsed the scene to a floater dot ("no discernible objects"). For gaussian scenes the viewer now frames from the input camera (the world origin, where the pose head places the single view) looking at the content's median, and clamps orbit (+-0.3 rad) and zoom to stay near that view. WASD panning is disabled there so the fixed input view can't be nudged off. Points/ streaming scenes are untouched: free orbit + bbox framing. Result: gaussian mode now renders the recognizable input view instead of floaters. --- server/web/index.html | 41 +++++++++++++++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 8 deletions(-) diff --git a/server/web/index.html b/server/web/index.html index 40ae5cc..3e97841 100644 --- a/server/web/index.html +++ b/server/web/index.html @@ -233,18 +233,40 @@

depth-anything.cpp · DA3 studio

const _p=new URLSearchParams(location.search); let target=[0,0,0], yaw=+(_p.get("yaw")??0), pitch=+(_p.get("pitch")??0.05), dist=2.0, distLock=false; -function frameToData(){ if(framed||count===0) return; let lo=[1e30,1e30,1e30],hi=[-1e30,-1e30,-1e30]; +let curMode="points", orbitLimit=null; // orbitLimit set for gaussian (input-view) scenes +function frameToData(){ if(framed||count===0) return; + if(curMode==="gaussians"){ frameInputView(); framed=true; sortDirty=true; return; } + let lo=[1e30,1e30,1e30],hi=[-1e30,-1e30,-1e30]; for(let i=0;ihi[k])hi[k]=v; } const diag=Math.hypot(hi[0]-lo[0],hi[1]-lo[1],hi[2]-lo[2]); target=[(lo[0]+hi[0])/2,(lo[1]+hi[1])/2,(lo[2]+hi[2])/2]; if(!distLock) dist=Math.max(0.5,0.7*diag); - framed=true; sortDirty=true; } + orbitLimit=null; framed=true; sortDirty=true; } +// DA3 single-image gaussians are view-synthesis primitives: their depth is spread, +// but they composite into the sharp INPUT view (the camera sits at the world origin). +// So frame from the origin looking at the content's median, and clamp the orbit to +// stay near that view instead of orbiting into the (expected) depth-spread floaters. +function frameInputView(){ + const step=Math.max(1,Math.floor(count/4000)); const xs=[],ys=[],zs=[]; + for(let i=0;i{ a.sort((p,q)=>p-q); return a[a.length>>1]; }; + const c=[med(xs),med(ys),med(zs)]; target=c; + const D=Math.hypot(c[0],c[1],c[2])||1; dist=D; // eye lands exactly at the origin + const dir=[-c[0]/D,-c[1]/D,-c[2]/D]; // target -> eye(origin) + pitch=Math.asin(Math.max(-1,Math.min(1,dir[1]))); yaw=Math.atan2(dir[0],dir[2]); + orbitLimit={yaw0:yaw,pitch0:pitch,d:0.3,dmin:D*0.5,dmax:D*1.25}; +} const keys={}; addEventListener("keydown",e=>keys[e.key.toLowerCase()]=true); addEventListener("keyup",e=>keys[e.key.toLowerCase()]=false); let drag=false,lx=0,ly=0; canvas.addEventListener("pointerdown",e=>{drag=true;lx=e.clientX;ly=e.clientY;canvas.setPointerCapture(e.pointerId);}); canvas.addEventListener("pointerup",()=>drag=false); canvas.addEventListener("pointermove",e=>{ if(!drag)return; yaw-=(e.clientX-lx)*0.005; pitch-=(e.clientY-ly)*0.005; - pitch=Math.max(-1.5,Math.min(1.5,pitch)); lx=e.clientX; ly=e.clientY; sortDirty=true; }); -canvas.addEventListener("wheel",e=>{ dist*=Math.exp(e.deltaY*0.001); dist=Math.max(0.05,Math.min(50,dist)); distLock=true; sortDirty=true; e.preventDefault(); },{passive:false}); + pitch=Math.max(-1.5,Math.min(1.5,pitch)); + if(orbitLimit){ const L=orbitLimit; + yaw=Math.max(L.yaw0-L.d,Math.min(L.yaw0+L.d,yaw)); pitch=Math.max(L.pitch0-L.d,Math.min(L.pitch0+L.d,pitch)); } + lx=e.clientX; ly=e.clientY; sortDirty=true; }); +canvas.addEventListener("wheel",e=>{ dist*=Math.exp(e.deltaY*0.001); dist=Math.max(0.05,Math.min(50,dist)); + if(orbitLimit) dist=Math.max(orbitLimit.dmin,Math.min(orbitLimit.dmax,dist)); + distLock=true; sortDirty=true; e.preventDefault(); },{passive:false}); function perspective(f,a,n,fr){ const t=1/Math.tan(f/2); return new Float32Array([t/a,0,0,0,0,t,0,0,0,0,(fr+n)/(n-fr),-1,0,0,(2*fr*n)/(n-fr),0]); } const sub=(a,b)=>[a[0]-b[0],a[1]-b[1],a[2]-b[2]],dot=(a,b)=>a[0]*b[0]+a[1]*b[1]+a[2]*b[2]; const cross=(a,b)=>[a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]]; @@ -257,10 +279,12 @@

depth-anything.cpp · DA3 studio

function render(){ const speed=dist*0.02; const fwd=[Math.sin(yaw)*Math.cos(pitch),Math.sin(pitch),Math.cos(yaw)*Math.cos(pitch)]; const right=norm(cross([0,1,0],fwd)); - if(keys["w"]||keys["arrowup"]){for(let i=0;i<3;i++)target[i]-=fwd[i]*speed;sortDirty=true;} - if(keys["s"]||keys["arrowdown"]){for(let i=0;i<3;i++)target[i]+=fwd[i]*speed;sortDirty=true;} - if(keys["a"]||keys["arrowleft"]){for(let i=0;i<3;i++)target[i]-=right[i]*speed;sortDirty=true;} - if(keys["d"]||keys["arrowright"]){for(let i=0;i<3;i++)target[i]+=right[i]*speed;sortDirty=true;} + if(!orbitLimit){ // WASD panning would break the fixed input-view for gaussians + if(keys["w"]||keys["arrowup"]){for(let i=0;i<3;i++)target[i]-=fwd[i]*speed;sortDirty=true;} + if(keys["s"]||keys["arrowdown"]){for(let i=0;i<3;i++)target[i]+=fwd[i]*speed;sortDirty=true;} + if(keys["a"]||keys["arrowleft"]){for(let i=0;i<3;i++)target[i]-=right[i]*speed;sortDirty=true;} + if(keys["d"]||keys["arrowright"]){for(let i=0;i<3;i++)target[i]+=right[i]*speed;sortDirty=true;} + } const eye=[target[0]+dist*Math.sin(yaw)*Math.cos(pitch),target[1]+dist*Math.sin(pitch),target[2]+dist*Math.cos(yaw)*Math.cos(pitch)]; const view=lookAt(eye,target); const aspect=canvas.width/canvas.height; const fovy=1.0; const proj=perspective(fovy,aspect,0.05,100); const fy=canvas.height/(2*Math.tan(fovy/2)),fx=fy; @@ -394,6 +418,7 @@

depth-anything.cpp · DA3 studio

async function openScene(name){ curScene=name; document.querySelectorAll(".card").forEach(c=>c.classList.toggle("active",c.dataset.name===name)); const m=await (await fetch(`/scenes-assets/${name}/manifest.json`)).json(); + curMode=m.mode||"points"; steps=m.steps||[]; if(!steps.length) return; $("playrow").classList.remove("hide"); $("steprange").max=steps.length-1; stepI=steps.length-1; $("steprange").value=stepI; distLock=false; From cfcd641da0ff32c765973071fb3762d3208bb504 Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Wed, 1 Jul 2026 13:41:39 +0100 Subject: [PATCH 05/22] feat(gs): render single-image gaussians from the input camera + photo colour DA3's single-image gaussians only resolve into a sharp picture from the exact input viewpoint, and their learned SH-DC colour is a deliberately flat, near-grey base (higher SH bands masked x0.025). The demo baked SH-DC colour and viewed the depth-spread gaussians through a wrong, too-wide, non-anamorphic projection -- a scattered grey blob cloud ("extremely washed out"). - Colour each gaussian by its own source pixel (de-normalized p.chw, the exact model input), same convention as the point cloud -> faithful, saturated colour. - Expose the input camera intrinsics + processed size from da_capi_gaussians (ABI 6 -> 7); thread through the Go server into the scene manifest and the /api/reconstruct response as `cam`. - Viewer renders gaussian scenes from the origin (pose is canonical identity) through the input K (anamorphic fx/fy, pillarboxed to height), so every gaussian lands on its source pixel and the input view resolves. Small clamped orbit gives parallax. Points scenes unchanged (free orbit + bbox framing). Result: gaussian scenes composite into the recognizable, correctly-coloured input image (soft, as inherent to DA3's 224^2 single-image gaussians) instead of a washed-out blob cloud. Verified by headless GPU render on two scenes; points path unregressed. Host tests pass (sim3/stream_seam/linalg/reconstruct/gs_adapter, and capi on CPU with ABI 7; the ctest Vulkan aborts here are the known env issue). Co-Authored-By: Claude Opus 4.8 (1M context) --- include/da_capi.h | 15 +++++++++++---- server/engine.go | 17 ++++++++++++----- server/main.go | 4 +++- server/scenes.go | 17 ++++++++++++++++- server/web/index.html | 38 +++++++++++++++++++++++--------------- src/da_capi.cpp | 14 +++++++++++--- src/engine.cpp | 18 ++++++++++++++++++ src/gs_adapter.hpp | 11 +++++++++++ tests/test_capi.cpp | 2 +- 9 files changed, 106 insertions(+), 30 deletions(-) diff --git a/include/da_capi.h b/include/da_capi.h index b085efa..cfe3093 100644 --- a/include/da_capi.h +++ b/include/da_capi.h @@ -9,7 +9,9 @@ typedef struct da_ctx da_ctx; 4: added da_capi_load_nested (two-branch metric model). 5: added da_capi_points_multi (fused multi-view cloud) + da_capi_gaussians. 6: added da_capi_points_stream (sliding-window streaming cloud); removed the - unused da_capi_depth_pose_multi. */ + unused da_capi_depth_pose_multi. + 7: da_capi_gaussians also returns the input camera intrinsics + processed size + (out_intr[9], out_w, out_h) so the viewer can render from the input view. */ int da_capi_abi_version(void); da_ctx* da_capi_load(const char* gguf_path, int n_threads); /* NULL on failure */ /* Load a NESTED metric model from its two branches: the anyview (GIANT) GGUF and @@ -100,11 +102,16 @@ int da_capi_points_stream(da_ctx* ctx, const char** image_paths, int n_images, /* Single-image 3D GAUSSIANS (DA3-GIANT / GS-head models only; returns -1 with a clear last_error otherwise). Returns world-frame (OpenCV) gaussians as parallel arrays: *out_xyz[3N] means, *out_scale[3N], *out_quat[4N] (wxyz), *out_rgb[3N] - linear colour in [0,1] from the SH DC term, *out_opacity[N]. Sets *out_n. Free - every returned float buffer via da_capi_free_floats. Returns 0 ok, -1 error. */ + linear colour in [0,1] (input-photo colour; see note), *out_opacity[N]. Sets + *out_n. Also, when non-NULL, writes the input camera intrinsics to out_intr[9] + (K 3x3 row-major, pixel units) and the processed resolution to *out_w,*out_h -- + the pose is canonical (identity) so a viewer placed at the origin looking down + +z with this K reproduces the input view. Free every returned float buffer via + da_capi_free_floats. Returns 0 ok, -1 error. */ int da_capi_gaussians(da_ctx* ctx, const char* image_path, int* out_n, float** out_xyz, float** out_scale, float** out_quat, - float** out_rgb, float** out_opacity); + float** out_rgb, float** out_opacity, + float* out_intr, int* out_w, int* out_h); #ifdef __cplusplus } #endif diff --git a/server/engine.go b/server/engine.go index 1bd83c0..3744baa 100644 --- a/server/engine.go +++ b/server/engine.go @@ -1,5 +1,5 @@ // engine.go — purego bindings to libdepthanything (the flat C API in -// include/da_capi.h, ABI 6). No cgo: the shared library is dlopen'd once and each +// include/da_capi.h, ABI 7). No cgo: the shared library is dlopen'd once and each // model is loaded as its own da_ctx. A context is NOT thread-safe (one ggml // backend + compute graph), so all inference is serialized by the caller. package main @@ -29,7 +29,8 @@ type capi struct { ptSize float32, budget int32, outN, outCounts, outXyz, outRgb, outRad unsafe.Pointer) int32 gaussians func(ctx uintptr, path string, - outN, outXyz, outScale, outQuat, outRgb, outOpacity unsafe.Pointer) int32 + outN, outXyz, outScale, outQuat, outRgb, outOpacity, + outIntr, outW, outH unsafe.Pointer) int32 } func openCAPI(libPath string) (*capi, error) { @@ -175,22 +176,28 @@ type GS struct { Quat []float32 // 4N wxyz RGB []float32 // 3N linear [0,1] Opacity []float32 // N + // Input camera (pose is canonical identity): K at the processed resolution. + Fx, Fy, Cx, Cy float32 + W, H int } // Gaussians runs the GS reconstruction on one image (DA3-GIANT only). func (a *capi) Gaussians(ctx uintptr, path string) (*GS, error) { - var n int32 + var n, w, h int32 var pXyz, pScale, pQuat, pRgb, pOp uintptr + var intr [9]float32 rc := a.gaussians(ctx, path, unsafe.Pointer(&n), unsafe.Pointer(&pXyz), unsafe.Pointer(&pScale), unsafe.Pointer(&pQuat), - unsafe.Pointer(&pRgb), unsafe.Pointer(&pOp)) + unsafe.Pointer(&pRgb), unsafe.Pointer(&pOp), + unsafe.Pointer(&intr[0]), unsafe.Pointer(&w), unsafe.Pointer(&h)) if rc != 0 { return nil, fmt.Errorf("gaussians: %s", a.lastErr(ctx)) } np := int(n) g := &GS{N: np, XYZ: cFloats(pXyz, np*3), Scale: cFloats(pScale, np*3), Quat: cFloats(pQuat, np*4), - RGB: cFloats(pRgb, np*3), Opacity: cFloats(pOp, np)} + RGB: cFloats(pRgb, np*3), Opacity: cFloats(pOp, np), + Fx: intr[0], Fy: intr[4], Cx: intr[2], Cy: intr[5], W: int(w), H: int(h)} a.freeFloats(pXyz) a.freeFloats(pScale) a.freeFloats(pQuat) diff --git a/server/main.go b/server/main.go index 4481e91..2c375f0 100644 --- a/server/main.go +++ b/server/main.go @@ -200,6 +200,7 @@ func (s *server) handleReconstruct(w http.ResponseWriter, r *http.Request) { t0 := time.Now() var splat []byte var nOut int + var cam *Cam err := s.infer(func() error { return s.reg.WithModel(model, func(ctx uintptr, mi *ModelInfo) error { if mode == "gaussians" { @@ -217,6 +218,7 @@ func (s *server) handleReconstruct(w http.ResponseWriter, r *http.Request) { } splat = gaussiansToSplat(g, s.maxSplats) nOut = g.N + cam = &Cam{Fx: g.Fx, Fy: g.Fy, Cx: g.Cx, Cy: g.Cy, W: g.W, H: g.H} return nil } c, e := s.api.PointsMulti(ctx, paths, confPct, ptSize) @@ -238,7 +240,7 @@ func (s *server) handleReconstruct(w http.ResponseWriter, r *http.Request) { } id := s.storeResult(splat) writeJSON(w, map[string]any{ - "id": id, "model": model, "mode": mode, + "id": id, "model": model, "mode": mode, "cam": cam, "n": nOut, "size": len(splat), "seconds": time.Since(t0).Seconds(), }) } diff --git a/server/scenes.go b/server/scenes.go index 93e5de2..aa66f91 100644 --- a/server/scenes.go +++ b/server/scenes.go @@ -25,10 +25,23 @@ type manifestStep struct { Label string `json:"label,omitempty"` } +// Cam is the input camera for gaussian scenes: the pose is canonical (identity), +// so the viewer places its eye at the origin looking down +z and uses this K +// (pixel units at W×H) to reproduce the exact input view the gaussians composite to. +type Cam struct { + Fx float32 `json:"fx"` + Fy float32 `json:"fy"` + Cx float32 `json:"cx"` + Cy float32 `json:"cy"` + W int `json:"w"` + H int `json:"h"` +} + type sceneManifest struct { Model string `json:"model"` Mode string `json:"mode"` Steps []manifestStep `json:"steps"` + Cam *Cam `json:"cam,omitempty"` } type sceneInfo struct { @@ -264,6 +277,7 @@ func (s *server) bakeVideo(name, vpath, jobDir, model, mode string, maxFrames, c return e } var splat []byte + var cam *Cam err := s.infer(func() error { return s.reg.WithModel(model, func(ctx uintptr, _ *ModelInfo) error { g, e := s.api.Gaussians(ctx, g224) @@ -271,6 +285,7 @@ func (s *server) bakeVideo(name, vpath, jobDir, model, mode string, maxFrames, c return e } splat = gaussiansToSplat(g, s.maxSplats) + cam = &Cam{Fx: g.Fx, Fy: g.Fy, Cx: g.Cx, Cy: g.Cy, W: g.W, H: g.H} return nil }) }) @@ -281,7 +296,7 @@ func (s *server) bakeVideo(name, vpath, jobDir, model, mode string, maxFrames, c return e } s.setJob(name, func(j *sceneJob) { j.Done, j.Kept = len(sel), 1 }) - return writeManifest(sceneDir, sceneManifest{Model: model, Mode: mode, Steps: []manifestStep{ + return writeManifest(sceneDir, sceneManifest{Model: model, Mode: mode, Cam: cam, Steps: []manifestStep{ {Splat: "acc_full.splat", Images: []string{thumb}, N: 1, Label: "gaussians · single frame"}}}) } diff --git a/server/web/index.html b/server/web/index.html index 3e97841..6b66716 100644 --- a/server/web/index.html +++ b/server/web/index.html @@ -234,6 +234,7 @@

depth-anything.cpp · DA3 studio

const _p=new URLSearchParams(location.search); let target=[0,0,0], yaw=+(_p.get("yaw")??0), pitch=+(_p.get("pitch")??0.05), dist=2.0, distLock=false; let curMode="points", orbitLimit=null; // orbitLimit set for gaussian (input-view) scenes +let gsCam=null; // {fx,fy,cx,cy,w,h} input camera for gaussian scenes (pose = identity) function frameToData(){ if(framed||count===0) return; if(curMode==="gaussians"){ frameInputView(); framed=true; sortDirty=true; return; } let lo=[1e30,1e30,1e30],hi=[-1e30,-1e30,-1e30]; @@ -241,19 +242,15 @@

depth-anything.cpp · DA3 studio

const diag=Math.hypot(hi[0]-lo[0],hi[1]-lo[1],hi[2]-lo[2]); target=[(lo[0]+hi[0])/2,(lo[1]+hi[1])/2,(lo[2]+hi[2])/2]; if(!distLock) dist=Math.max(0.5,0.7*diag); orbitLimit=null; framed=true; sortDirty=true; } -// DA3 single-image gaussians are view-synthesis primitives: their depth is spread, -// but they composite into the sharp INPUT view (the camera sits at the world origin). -// So frame from the origin looking at the content's median, and clamp the orbit to -// stay near that view instead of orbiting into the (expected) depth-spread floaters. +// DA3 single-image gaussians are view-synthesis primitives: their per-pixel depth +// is deliberately spread, but they composite into the sharp INPUT view. The pose is +// canonical (identity), so the input camera sits at the world origin looking down +z +// (=-z after the OpenCV->GL flip). Put the eye exactly there and render through the +// input intrinsics (gsCam.fx/fy, anamorphic) so every gaussian lands on its own pixel +// and the picture resolves. A small orbit around a point 1 unit ahead gives parallax. function frameInputView(){ - const step=Math.max(1,Math.floor(count/4000)); const xs=[],ys=[],zs=[]; - for(let i=0;i{ a.sort((p,q)=>p-q); return a[a.length>>1]; }; - const c=[med(xs),med(ys),med(zs)]; target=c; - const D=Math.hypot(c[0],c[1],c[2])||1; dist=D; // eye lands exactly at the origin - const dir=[-c[0]/D,-c[1]/D,-c[2]/D]; // target -> eye(origin) - pitch=Math.asin(Math.max(-1,Math.min(1,dir[1]))); yaw=Math.atan2(dir[0],dir[2]); - orbitLimit={yaw0:yaw,pitch0:pitch,d:0.3,dmin:D*0.5,dmax:D*1.25}; + target=[0,0,-1]; yaw=0; pitch=0; dist=1; // eye at the origin, looking down -z + orbitLimit={yaw0:0,pitch0:0,d:0.35,dmin:0.35,dmax:1.6}; } const keys={}; addEventListener("keydown",e=>keys[e.key.toLowerCase()]=true); addEventListener("keyup",e=>keys[e.key.toLowerCase()]=false); let drag=false,lx=0,ly=0; @@ -286,8 +283,18 @@

depth-anything.cpp · DA3 studio

if(keys["d"]||keys["arrowright"]){for(let i=0;i<3;i++)target[i]+=right[i]*speed;sortDirty=true;} } const eye=[target[0]+dist*Math.sin(yaw)*Math.cos(pitch),target[1]+dist*Math.sin(pitch),target[2]+dist*Math.cos(yaw)*Math.cos(pitch)]; - const view=lookAt(eye,target); const aspect=canvas.width/canvas.height; const fovy=1.0; - const proj=perspective(fovy,aspect,0.05,100); const fy=canvas.height/(2*Math.tan(fovy/2)),fx=fy; + const view=lookAt(eye,target); const cw=canvas.width,ch=canvas.height; const aspect=cw/ch; + let proj,fx,fy; + if(gsCam && curMode==="gaussians"){ + // Render through the input camera's (anamorphic) intrinsics, fitting the WxH + // image to canvas height (pillarboxed). p00/p11 map camera x/y to NDC exactly + // as the model's K would to pixels, so gaussians land on their source pixels. + const p11=2*gsCam.fy/gsCam.h, p00=2*gsCam.fx*ch/(gsCam.h*cw); + const n=0.05,f=500; proj=new Float32Array([p00,0,0,0, 0,p11,0,0, 0,0,(f+n)/(n-f),-1, 0,0,(2*f*n)/(n-f),0]); + fx=p00*cw/2; fy=p11*ch/2; + }else{ + const fovy=1.0; proj=perspective(fovy,aspect,0.05,100); fy=ch/(2*Math.tan(fovy/2)); fx=fy; + } if(sortDirty){ requestSort(view); sortDirty=false; } gl.clearColor(0.043,0.051,0.063,1); gl.clear(gl.COLOR_BUFFER_BIT); gl.useProgram(prog); gl.uniformMatrix4fv(U.view,false,view); gl.uniformMatrix4fv(U.proj,false,proj); @@ -390,6 +397,7 @@

depth-anything.cpp · DA3 studio

if(!resp.ok) throw new Error(await resp.text()); const r=await resp.json(); setStatus(`done · ${r.n.toLocaleString()} splats · ${r.seconds.toFixed(1)}s`); + curMode=mode; gsCam=r.cam||null; await loadSplat("/api/splat/"+r.id); } }catch(err){ setStatus(err.message,true); } @@ -418,7 +426,7 @@

depth-anything.cpp · DA3 studio

async function openScene(name){ curScene=name; document.querySelectorAll(".card").forEach(c=>c.classList.toggle("active",c.dataset.name===name)); const m=await (await fetch(`/scenes-assets/${name}/manifest.json`)).json(); - curMode=m.mode||"points"; + curMode=m.mode||"points"; gsCam=m.cam||null; steps=m.steps||[]; if(!steps.length) return; $("playrow").classList.remove("hide"); $("steprange").max=steps.length-1; stepI=steps.length-1; $("steprange").value=stepI; distLock=false; diff --git a/src/da_capi.cpp b/src/da_capi.cpp index 2ac0db5..62f7690 100644 --- a/src/da_capi.cpp +++ b/src/da_capi.cpp @@ -70,7 +70,7 @@ static bool capi_run_nested(da_ctx* c, const char* image_path, } extern "C" { -int da_capi_abi_version(void){ return 6; } +int da_capi_abi_version(void){ return 7; } da_ctx* da_capi_load(const char* path, int n_threads){ if (!path) return nullptr; auto e = da::Engine::load(path, n_threads); @@ -398,7 +398,8 @@ int da_capi_points_stream(da_ctx* c, const char** image_paths, int n_images, int da_capi_gaussians(da_ctx* c, const char* path, int* out_n, float** out_xyz, float** out_scale, float** out_quat, - float** out_rgb, float** out_opacity){ + float** out_rgb, float** out_opacity, + float* out_intr, int* out_w, int* out_h){ if (out_xyz) *out_xyz = nullptr; if (out_scale) *out_scale = nullptr; if (out_quat) *out_quat = nullptr; @@ -410,6 +411,9 @@ int da_capi_gaussians(da_ctx* c, const char* path, int* out_n, da::Gaussians g; int H = 0, W = 0; if (!c->engine->reconstruct(img, g, H, W)){ c->last_error = "gaussians: reconstruct failed (needs a GS model, e.g. DA3-GIANT)"; return -1; } + if (out_intr) for (int i = 0; i < 9; ++i) out_intr[i] = g.intr[i]; + if (out_w) *out_w = W; + if (out_h) *out_h = H; const int N = g.N; if (N <= 0 || (int)g.means.size() < N*3 || (int)g.scales.size() < N*3 || (int)g.rotations.size() < N*4 || (int)g.opacities.size() < N || (int)g.harmonics.size() < N*3*9){ @@ -423,11 +427,15 @@ int da_capi_gaussians(da_ctx* c, const char* path, int* out_n, if (!xyz || !scl || !quat || !rgb || !op){ std::free(xyz); std::free(scl); std::free(quat); std::free(rgb); std::free(op); c->last_error = "gaussians: oom"; return -1; } + // Prefer the input-photo colour baked by Engine::reconstruct; SH-DC is a + // near-grey flat base (see gs_adapter.hpp). Fall back to SH-DC if absent. + const bool have_col = (int)g.colors.size() >= N*3; for (int i = 0; i < N; ++i){ for (int k = 0; k < 3; ++k){ xyz[3*i+k] = g.means[3*i+k]; scl[3*i+k] = g.scales[3*i+k]; } for (int k = 0; k < 4; ++k) quat[4*i+k] = g.rotations[4*i+k]; for (int ch = 0; ch < 3; ++ch){ - double col = 0.5 + SH_C0 * (double)g.harmonics[((size_t)i*3 + ch)*9 + 0]; + double col = have_col ? (double)g.colors[3*i+ch] + : 0.5 + SH_C0 * (double)g.harmonics[((size_t)i*3 + ch)*9 + 0]; rgb[3*i+ch] = (float)(col < 0 ? 0 : (col > 1 ? 1 : col)); } op[i] = g.opacities[i]; diff --git a/src/engine.cpp b/src/engine.cpp index cb9e451..c499a1b 100644 --- a/src/engine.cpp +++ b/src/engine.cpp @@ -307,6 +307,24 @@ bool Engine::reconstruct(const Image& img, Gaussians& g, int& H, int& W){ if (!gs.raw_gaussians(feats, p.chw, H, W, raw_gs, gs_conf)) { DA_LOG("reconstruct: gs_head failed"); return false; } GsAdapter ad; if (!ad.build(raw_gs, depth, gs_conf, ext, intr, H, W, g)) { DA_LOG("reconstruct: gs_adapter failed"); return false; } + g.ext = ext; g.intr = intr; g.H = H; g.W = W; // input camera, for input-view rendering + // Colour each gaussian by its own source pixel (de-normalize p.chw, the exact + // model input). DA3's learned SH-DC colour is a near-grey flat base; the true + // photo colour gives a faithful input-view presentation (same convention as the + // point cloud). g.colors empty => callers fall back to SH-DC. + { + const auto& mean = ml_.config().img_mean; + const auto& std = ml_.config().img_std; + if (mean.size() >= 3 && std.size() >= 3 && (int)p.chw.size() >= 3*H*W && g.N == H*W) { + const size_t HW = (size_t)H * W; + g.colors.resize((size_t)g.N * 3); + for (size_t pix = 0; pix < HW; ++pix) + for (int c = 0; c < 3; ++c) { + float v = p.chw[(size_t)c*HW + pix] * std[c] + mean[c]; + g.colors[pix*3 + c] = v < 0.f ? 0.f : (v > 1.f ? 1.f : v); + } + } + } return true; } bool Engine::reconstruct_path(const std::string& image_path, Gaussians& g, int& H, int& W){ diff --git a/src/gs_adapter.hpp b/src/gs_adapter.hpp index 476264a..d556da1 100644 --- a/src/gs_adapter.hpp +++ b/src/gs_adapter.hpp @@ -13,6 +13,17 @@ struct Gaussians { std::vector rotations; // [N*4] wxyz (pixel*4 + c) std::vector harmonics; // [N*3*9] (pixel*3 + color)*9 + coeff std::vector opacities; // [N] + // Optional per-gaussian input-photo RGB in [0,1] (pixel*3 + color), filled by + // Engine::reconstruct from the processed image. DA3's learned SH-DC colour is a + // deliberately flat, near-grey base (higher SH bands masked x0.025), so for the + // single-image input-view presentation we colour each gaussian by its own source + // pixel -- same convention as the point cloud. Empty => fall back to SH-DC. + std::vector colors; // [N*3] + // Input camera (from the model's pose head), so the demo can render the + // gaussians from the exact viewpoint they composite into. + std::array ext{}; // world->camera 3x4 (w2c), row-major + std::array intr{}; // K 3x3 (pixel units at processed W,H), row-major + int H = 0, W = 0; // processed resolution the gaussians live at }; // GaussianAdapter: pure host math (NO learned weights). Converts the GSDPT diff --git a/tests/test_capi.cpp b/tests/test_capi.cpp index b7fb3df..813a810 100644 --- a/tests/test_capi.cpp +++ b/tests/test_capi.cpp @@ -5,7 +5,7 @@ int main(){ const char* gguf = std::getenv("DA_TEST_GGUF"); if (!gguf) return 77; - if (da_capi_abi_version() != 6) return 1; + if (da_capi_abi_version() != 7) return 1; da_ctx* c = da_capi_load(gguf, 1); if (!c) { std::fprintf(stderr, "load failed\n"); return 1; } char* j = da_capi_info_json(c); From a899ad45a7412c7ded1205bd498381bb3d4a0e9b Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Sun, 5 Jul 2026 10:21:13 +0100 Subject: [PATCH 06/22] =?UTF-8?q?feat(geom):=20streaming=20de-ghosting=20?= =?UTF-8?q?=E2=80=94=20surface=20fusion,=20per-seam=20ICP,=20loop-closure?= =?UTF-8?q?=20pose-graph?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Collapse the doubled sheets and drift that appear when stitching many sliding-window reconstructions into one coherent cloud: - A: voxel surface fusion (fuse.*) — de-densify + merge coincident sheets. - B: per-seam point-to-plane ICP (icp.*, spatial_hash.hpp) — tighten each window-to-window seam before it accumulates. - C: loop-closure detection + Sim3 pose-graph optimisation (posegraph.*, sim3.*) — remove global drift when the trajectory revisits a place. All params are scene-relative (monocular DA3 lives at arbitrary scale). CPU geometry loops are OpenMP-parallelised. Includes analytic self-tests, synthetic scene/trajectory generators, open3d/gtsam reference-oracle scripts, and standalone dump harnesses. Co-Authored-By: Claude Opus 4.8 --- CMakeLists.txt | 15 +++ scripts/oracle.sh | 68 ++++++++++ scripts/oracle_fuse.py | 89 +++++++++++++ scripts/oracle_icp.py | 114 ++++++++++++++++ scripts/oracle_posegraph.py | 134 +++++++++++++++++++ scripts/stream_regress.sh | 27 ++++ src/fuse.cpp | 133 +++++++++++++++++++ src/fuse.hpp | 39 ++++++ src/icp.cpp | 186 ++++++++++++++++++++++++++ src/icp.hpp | 61 +++++++++ src/posegraph.cpp | 207 +++++++++++++++++++++++++++++ src/posegraph.hpp | 62 +++++++++ src/sim3.cpp | 17 +++ src/sim3.hpp | 3 + src/spatial_hash.hpp | 64 +++++++++ tests/fuse_dump.cpp | 37 ++++++ tests/geom_metrics.hpp | 158 ++++++++++++++++++++++ tests/icp_dump.cpp | 54 ++++++++ tests/posegraph_dump.cpp | 56 ++++++++ tests/stream_regress.cpp | 144 ++++++++++++++++++++ tests/synth_scene.hpp | 254 ++++++++++++++++++++++++++++++++++++ tests/test_fuse.cpp | 114 ++++++++++++++++ tests/test_geom_metrics.cpp | 155 ++++++++++++++++++++++ tests/test_icp.cpp | 193 +++++++++++++++++++++++++++ tests/test_posegraph.cpp | 149 +++++++++++++++++++++ tests/test_stream_loop.cpp | 143 ++++++++++++++++++++ tests/test_stream_synth.cpp | 217 ++++++++++++++++++++++++++++++ tests/test_synth_scene.cpp | 111 ++++++++++++++++ 28 files changed, 3004 insertions(+) create mode 100755 scripts/oracle.sh create mode 100644 scripts/oracle_fuse.py create mode 100644 scripts/oracle_icp.py create mode 100644 scripts/oracle_posegraph.py create mode 100755 scripts/stream_regress.sh create mode 100644 src/fuse.cpp create mode 100644 src/fuse.hpp create mode 100644 src/icp.cpp create mode 100644 src/icp.hpp create mode 100644 src/posegraph.cpp create mode 100644 src/posegraph.hpp create mode 100644 src/spatial_hash.hpp create mode 100644 tests/fuse_dump.cpp create mode 100644 tests/geom_metrics.hpp create mode 100644 tests/icp_dump.cpp create mode 100644 tests/posegraph_dump.cpp create mode 100644 tests/stream_regress.cpp create mode 100644 tests/synth_scene.hpp create mode 100644 tests/test_fuse.cpp create mode 100644 tests/test_geom_metrics.cpp create mode 100644 tests/test_icp.cpp create mode 100644 tests/test_posegraph.cpp create mode 100644 tests/test_stream_loop.cpp create mode 100644 tests/test_stream_synth.cpp create mode 100644 tests/test_synth_scene.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index e11d2a9..31dd5de 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -65,6 +65,9 @@ set(DA_SOURCES src/linalg.cpp src/sim3.cpp src/stream.cpp + src/fuse.cpp + src/icp.cpp + src/posegraph.cpp src/ray_pose.cpp src/nested.cpp src/depth_export.cpp @@ -94,6 +97,18 @@ if(DA_HAS_AVX512F) endif() target_link_libraries(depthanything PUBLIC ggml) +# OpenMP for the single-threaded CPU geometry post-processing (surface fusion, +# per-seam ICP, normal estimation). ggml already pulls OpenMP into the toolchain; +# link it into OUR target too so the #pragma omp loops in fuse.cpp/icp.cpp use all +# cores instead of one. Optional: without it the pragmas compile to no-ops. +find_package(OpenMP QUIET) +if(OpenMP_CXX_FOUND) + target_link_libraries(depthanything PUBLIC OpenMP::OpenMP_CXX) + message(STATUS "depth-anything.cpp: OpenMP enabled for CPU geometry") +else() + message(STATUS "depth-anything.cpp: OpenMP NOT found — CPU geometry stays single-threaded") +endif() + if(DA_BUILD_CLI AND EXISTS ${DA_ROOT}/examples/cli/CMakeLists.txt) add_subdirectory(examples/cli) endif() diff --git a/scripts/oracle.sh b/scripts/oracle.sh new file mode 100755 index 0000000..dd387d7 --- /dev/null +++ b/scripts/oracle.sh @@ -0,0 +1,68 @@ +#!/usr/bin/env bash +# Reference-oracle Python runner for the streaming-alignment validation suite +# (tasks A/B/C). Provides open3d (ICP + voxel fusion + pose-graph), gtsam +# (Sim3 pose-graph), and scipy/numpy — the independent implementations we +# differentially test our C++ against. +# +# NixOS note: the manylinux open3d wheel dlopen's X11/GL at import, which are +# not on the default loader path here. We assemble LD_LIBRARY_PATH from the +# store (cached next to the venv) so the wheel loads. gtsam pins numpy<2, so +# the whole stack is pinned to a numpy 1.26 / scipy 1.13 line. +# +# Usage: +# scripts/oracle.sh setup # (re)create venv + install pinned deps +# scripts/oracle.sh python script.py # run a script with the oracle env +# scripts/oracle.sh -c "import open3d" # inline +set -euo pipefail + +ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +VENV="$ROOT/.oracle-venv" +PY="$VENV/bin/python" +LDP_CACHE="$VENV/.ldpath" + +# Libraries the open3d wheel needs at import (X11/GL cluster + toolchain). +NEED_LIBS=( + libX11.so.6 libXext.so.6 libXrandr.so.2 libXinerama.so.1 libXcursor.so.1 + libXi.so.6 libXfixes.so.3 libXrender.so.1 libGL.so.1 libEGL.so.1 + libGLX.so.0 libGLdispatch.so.0 libgomp.so.1 libstdc++.so.6 libgcc_s.so.1 + libxcb.so.1 libXau.so.6 libXdmcp.so.6 +) + +# True (0) iff $1 is a 64-bit ELF (byte 4 of the header == 2). The store also +# holds 32-bit multilib copies (lib32/) that must not leak onto the path. +is_elf64() { + [ "$(od -An -t u1 -j4 -N1 "$1" 2>/dev/null | tr -d ' ')" = "2" ] +} + +compute_ldpath() { + local dirs="" cand + for lib in "${NEED_LIBS[@]}"; do + while IFS= read -r cand; do + if is_elf64 "$cand"; then dirs="$dirs:$(dirname "$cand")"; break; fi + done < <(find /nix/store -maxdepth 3 -name "$lib" 2>/dev/null) + done + echo "$dirs" | tr ':' '\n' | sort -u | grep -v '^$' | paste -sd: +} + +ensure_ldpath() { + if [ ! -s "$LDP_CACHE" ]; then compute_ldpath > "$LDP_CACHE"; fi + cat "$LDP_CACHE" +} + +case "${1:-python}" in + setup) + uv venv --python 3.12 "$VENV" + uv pip install --python "$PY" \ + "numpy==1.26.4" "scipy==1.13.1" "open3d==0.19.0" "gtsam==4.2.1" + compute_ldpath > "$LDP_CACHE" + echo "oracle env ready: $VENV" + ;; + python) + shift + LD_LIBRARY_PATH="$(ensure_ldpath):${LD_LIBRARY_PATH:-}" exec "$PY" "$@" + ;; + *) + # treat all args as python args (e.g. -c "...") + LD_LIBRARY_PATH="$(ensure_ldpath):${LD_LIBRARY_PATH:-}" exec "$PY" "$@" + ;; +esac diff --git a/scripts/oracle_fuse.py b/scripts/oracle_fuse.py new file mode 100644 index 0000000..f25d549 --- /dev/null +++ b/scripts/oracle_fuse.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 +# Differential oracle for task A (src/fuse.cpp). Diffs our C++ fusion against +# reference implementations on a SURFACE-concentrated cloud (planes + sphere): +# * voxel stage vs open3d.voxel_down_sample +# * radius stage vs an independent scipy cKDTree neighbourhood-centroid +# Run via: scripts/oracle.sh python scripts/oracle_fuse.py +import os, struct, subprocess, sys, tempfile +import numpy as np +import open3d as o3d +from scipy.spatial import cKDTree + +DUMP = sys.argv[1] if len(sys.argv) > 1 else "./build/tests/fuse_dump" + +def surface_cloud(seed=0): + rng = np.random.default_rng(seed) + parts = [] + # floor z=0, +x wall x=2, -y wall y=-2 (each a thin noisy sheet) + xy = rng.uniform(-2, 2, (20000, 2)); parts.append(np.c_[xy, rng.normal(0, 0.004, 20000)]) + yz = rng.uniform(-2, 2, (20000, 2)); parts.append(np.c_[2 + rng.normal(0, 0.004, 20000), yz]) + xz = rng.uniform(-2, 2, (16000, 2)); parts.append(np.c_[xz[:,0], -2 + rng.normal(0, 0.004, 16000), xz[:,1]]) + d = rng.normal(size=(14000, 3)); d /= np.linalg.norm(d, axis=1, keepdims=True) # sphere r=1 + parts.append(d) + return np.vstack(parts).astype(np.float32) + +def run_dump(P, voxel, radius=0.0): + with tempfile.TemporaryDirectory() as td: + ip, op = os.path.join(td, "in.bin"), os.path.join(td, "out.bin") + with open(ip, "wb") as f: + f.write(struct.pack(" order preserved, N out +tree = cKDTree(Ps) +ref_s = np.empty_like(Ps, dtype=np.float64) +nbrs = tree.query_ball_point(Ps, r) +for i, nb in enumerate(nbrs): + ref_s[i] = Ps[nb].mean(axis=0) +maxerr = np.max(np.linalg.norm(ours_s - ref_s, axis=1)) +print(f"[radius] N={len(Ps)} max |ours - scipy_centroid| = {maxerr*1e6:.2f} um") +check(len(ours_s) == len(Ps), "radius stage preserves point count/order") +check(maxerr < 1e-4, "radius smoothing matches scipy neighbourhood centroid") + +print(("ORACLE PASS" if fails == 0 else f"ORACLE FAIL ({fails})")) +sys.exit(1 if fails else 0) diff --git a/scripts/oracle_icp.py b/scripts/oracle_icp.py new file mode 100644 index 0000000..3f6efef --- /dev/null +++ b/scripts/oracle_icp.py @@ -0,0 +1,114 @@ +#!/usr/bin/env python3 +# Differential oracle for task B (src/icp.cpp). On SURFACE-concentrated clouds it +# diffs our point-to-plane ICP against open3d.registration_icp (also point-to- +# plane) AND against the known ground-truth transform: +# * curvature-rich corner+sphere: both recover a KNOWN rigid transform, and our +# result agrees with open3d's to <0.1 deg / sub-mm; +# * single flat wall: both drive the point-to-plane residual to ~0 (the +# observable normal error), confirming the fit is correct where it is +# constrained (our rank-truncation of the unobservable tangential slide is +# asserted separately in tests/test_icp.cpp). +# Run via: scripts/oracle.sh python scripts/oracle_icp.py +import os, struct, subprocess, sys, tempfile +import numpy as np +import open3d as o3d +from scipy.spatial.transform import Rotation + +DUMP = sys.argv[1] if len(sys.argv) > 1 else "./build/tests/icp_dump" +MAXC, NRAD = 0.10, 0.06 + +fails = 0 +def check(ok, msg): + global fails + print(("ok: " if ok else "FAIL: ") + msg) + if not ok: fails += 1 + +def corner_scene(seed=0): + rng = np.random.default_rng(seed) + xy = rng.uniform(-0.5, 0.5, (3000, 2)); floor = np.c_[xy[:,0]+0.5, xy[:,1], np.zeros(3000)] + yz = rng.uniform(-0.5, 0.5, (3000, 2)); wall = np.c_[np.zeros(3000), yz[:,0], yz[:,1]+0.5] + d = rng.normal(size=(2500,3)); d /= np.linalg.norm(d, axis=1, keepdims=True) + sph = np.c_[d[:,0]*0.22+0.5, d[:,1]*0.22, d[:,2]*0.22+0.7] + return np.vstack([floor, wall, sph]).astype(np.float32) + +def wall_scene(seed=1): + rng = np.random.default_rng(seed) + xy = rng.uniform(-0.6, 0.6, (8000, 2)) + return np.c_[xy[:,0], xy[:,1], np.zeros(8000)].astype(np.float32) + +def write_cloud(path, P): + with open(path, "wb") as f: + f.write(struct.pack(" T ~ G^-1 +Ginv = np.linalg.inv(G) + +with tempfile.TemporaryDirectory() as td: + T_ours, s, rms, meta = run_ours(src, tgt, td) +T_o3d, reg = run_o3d(src, tgt) + +ang_ours_gt = rot_angle_deg(T_ours[:3,:3], Ginv[:3,:3]); t_ours_gt = np.linalg.norm(T_ours[:3,3]-Ginv[:3,3]) +ang_o3d_gt = rot_angle_deg(T_o3d[:3,:3], Ginv[:3,:3]); t_o3d_gt = np.linalg.norm(T_o3d[:3,3]-Ginv[:3,3]) +ang_pair = rot_angle_deg(T_ours[:3,:3], T_o3d[:3,:3]); t_pair = np.linalg.norm(T_ours[:3,3]-T_o3d[:3,3]) +print(f"[corner] ours vs GT : {ang_ours_gt:.4f} deg {t_ours_gt*1e3:.3f} mm (rank={meta[0]} iters={meta[1]} s={s:.4f})") +print(f"[corner] o3d vs GT : {ang_o3d_gt:.4f} deg {t_o3d_gt*1e3:.3f} mm (fitness={reg.fitness:.3f} rmse={reg.inlier_rmse*1e3:.3f}mm)") +print(f"[corner] ours vs o3d: {ang_pair:.4f} deg {t_pair*1e3:.3f} mm") +check(ang_ours_gt < 0.1 and t_ours_gt < 1e-3, "ours recovers known transform (<0.1deg, <1mm)") +check(ang_o3d_gt < 0.1 and t_o3d_gt < 1e-3, "open3d recovers known transform (<0.1deg, <1mm)") +check(ang_pair < 0.1 and t_pair < 1e-3, "ours agrees with open3d (<0.1deg, <1mm)") +check(rms[1] < 1e-3, "ours point-to-plane residual converged to ~0") + +# ================= 2. flat wall: both drive plane residual to ~0 ============ +wt = wall_scene() +Rw = Rotation.from_rotvec(np.radians(1.0) * np.array([1.0,0,0])) +Gw = np.eye(4); Gw[:3,:3] = Rw.as_matrix(); Gw[:3,3] = [0.05, 0.0, 0.03] # tangential x + normal z + tilt +ws = apply4(Gw, wt.astype(np.float64)).astype(np.float32) + +with tempfile.TemporaryDirectory() as td: + Tw_ours, _, rms_w, meta_w = run_ours(ws, wt, td) +Tw_o3d, reg_w = run_o3d(ws, wt) + +# residual to the true plane z=0 after each transform (the OBSERVABLE error). +def plane_resid_z(T, P): + return float(np.sqrt(np.mean(apply4(T, P.astype(np.float64))[:,2]**2))) +rz_ours, rz_o3d = plane_resid_z(Tw_ours, ws), plane_resid_z(Tw_o3d, ws) +print(f"[wall] plane resid ours={rz_ours*1e3:.3f}mm o3d={rz_o3d*1e3:.3f}mm (ours rank={meta_w[0]}, o3d rmse={reg_w.inlier_rmse*1e3:.3f}mm)") +check(rz_ours < 2e-3, "ours corrects the observable normal error (plane residual ~0)") +check(rz_o3d < 2e-3, "open3d corrects the observable normal error (plane residual ~0)") +check(meta_w[0] == 3, "ours flags the wall as rank-3 (aperture problem)") + +print("ORACLE PASS" if fails == 0 else f"ORACLE FAIL ({fails})") +sys.exit(1 if fails else 0) diff --git a/scripts/oracle_posegraph.py b/scripts/oracle_posegraph.py new file mode 100644 index 0000000..550bb76 --- /dev/null +++ b/scripts/oracle_posegraph.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python3 +# Differential oracle for task C (src/posegraph.cpp). Cross-checks our Sim3 pose- +# graph optimizer against gtsam's Pose3 pose-graph (BetweenFactorPose3 + +# PriorFactorPose3 + LevenbergMarquardt) on the SE3 (scale=1) case, which is what +# gtsam 4.2's Python wrapper can express (it has no BetweenFactorSimilarity3). +# 1. consistent ring + loop + perturbed init -> both recover the KNOWN poses, +# and our result agrees with gtsam node-for-node; +# 2. systematic drift + strong loop closure -> both cut ATE the same way. +# The Sim3 SCALE dimension (which gtsam Python can't graph-optimize) is validated +# analytically in tests/test_posegraph.cpp. +# scripts/oracle.sh python scripts/oracle_posegraph.py +import os, struct, subprocess, sys, tempfile +import numpy as np +import gtsam +from scipy.spatial.transform import Rotation + +DUMP = sys.argv[1] if len(sys.argv) > 1 else "./build/tests/posegraph_dump" +np.random.seed(7) + +fails = 0 +def check(ok, msg): + global fails + print(("ok: " if ok else "FAIL: ") + msg) + if not ok: fails += 1 + +# ---- SE3 helpers (4x4 homogeneous) --------------------------------------- +def se3(R, t): + M = np.eye(4); M[:3,:3] = R; M[:3,3] = t; return M +def inv(M): + R = M[:3,:3]; t = M[:3,3]; o = np.eye(4); o[:3,:3] = R.T; o[:3,3] = -R.T @ t; return o +def rand_se3(tsc, rsc): + return se3(Rotation.from_rotvec(rsc*np.random.randn(3)).as_matrix(), tsc*np.random.randn(3)) +def pose_diff(A, B): + dt = np.linalg.norm(A[:3,3] - B[:3,3]) + Rr = A[:3,:3].T @ B[:3,:3] + ang = np.degrees(np.arccos(np.clip((np.trace(Rr)-1)/2, -1, 1))) + return dt, ang +def gt_loop(N, radius): + P = [] + for i in range(N): + a = 2*np.pi*i/N + R = Rotation.from_rotvec(a*np.array([0,0,1.0])).as_matrix() + P.append(se3(R, [radius*np.cos(a), radius*np.sin(a), 0.2*np.sin(2*a)])) + return P + +# ---- pose-graph (de)serialization for our dump --------------------------- +def pack_T(M): # SE3 as Sim3 with s=1: 13 f64 [s, R0..8, t0..2] + R = M[:3,:3].reshape(-1); t = M[:3,3] + return struct.pack("<13d", 1.0, *R, *t) +def unpack_T(buf): + b = struct.unpack("<13d", buf); R = np.array(b[1:10]).reshape(3,3) + M = np.eye(4); M[:3,:3] = R; M[:3,3] = b[10:13]; return M +def run_ours(nodes, edges, fixed, td): + ip, op = os.path.join(td,"in.bin"), os.path.join(td,"out.bin") + with open(ip,"wb") as f: + f.write(struct.pack("{c:.2e}, it={it})") +print(f"[exact] gtsam vs GT: {d_gt[0]*1e3:.4f} mm {d_gt[1]:.5f} deg") +print(f"[exact] ours vs gts: {d_pair[0]*1e3:.4f} mm {d_pair[1]:.5f} deg") +check(d_ours[0] < 1e-4 and d_ours[1] < 1e-3, "ours recovers known poses (<0.1mm, <1e-3deg)") +check(d_gt[0] < 1e-4 and d_gt[1] < 1e-3, "gtsam recovers known poses (<0.1mm, <1e-3deg)") +check(d_pair[0] < 1e-4 and d_pair[1] < 1e-3, "ours agrees with gtsam node-for-node") + +# ================= 2. systematic drift + strong loop closure ================ +N = 16; GT = gt_loop(N, 3.0) +bias = se3(Rotation.from_rotvec(0.010*np.array([0,0,1.0])).as_matrix(), [0.010, 0, 0]) +edges = []; nodes = [GT[0].copy()] +for i in range(N-1): + Z = rel(GT[i], GT[i+1]) @ bias @ rand_se3(0.004, 0.004) # systematic + small iid + edges.append((i, i+1, 1.0, Z)) + nodes.append(nodes[i] @ Z) # chain => accumulates drift +edges.append((N-1, 0, 25.0, rel(GT[N-1], GT[0]))) # accurate loop closure +fixed = [0] +with tempfile.TemporaryDirectory() as td: + ours2, _, _, _ = run_ours(nodes, edges, fixed, td) +gt2 = run_gtsam(nodes, edges, fixed) +def ate(P): return np.sqrt(np.mean([np.linalg.norm(P[i][:3,3]-GT[i][:3,3])**2 for i in range(N)])) +a_open, a_ours, a_gt = ate(nodes), ate(ours2), ate(gt2) +d_pair2 = max(pose_diff(ours2[i], gt2[i]) for i in range(N)) +# On INCONSISTENT data the two optima differ by the parametrization of the +# translation error (we use E.t directly; gtsam's Pose3 Logmap folds in the SE3 +# left-Jacobian), so we cross-check EQUIVALENT BEHAVIOUR — both cut the drift to +# the same degree and agree to a small fraction of the corrected drift — not a +# bit-identical trajectory. +print(f"[drift] ATE open={a_open*1e3:.2f}mm ours={a_ours*1e3:.2f}mm gtsam={a_gt*1e3:.2f}mm") +print(f"[drift] ours vs gtsam: {d_pair2[0]*1e3:.3f} mm ({100*d_pair2[0]/a_open:.1f}% of drift) {d_pair2[1]:.4f} deg") +check(a_ours < 0.3*a_open, "ours: loop closure cuts systematic-drift ATE by >70%") +check(a_gt < 0.3*a_open, "gtsam: loop closure cuts systematic-drift ATE by >70%") +check(abs(a_ours-a_gt) < 0.2*max(a_ours,a_gt), "ours & gtsam fix the drift to the same degree (ATE within 20%)") +check(d_pair2[0] < 0.05*a_open and d_pair2[1] < 0.1, "ours agrees with gtsam to <5% of the corrected drift, <0.1deg") + +print("ORACLE PASS" if fails == 0 else f"ORACLE FAIL ({fails})") +sys.exit(1 if fails else 0) diff --git a/scripts/stream_regress.sh b/scripts/stream_regress.sh new file mode 100755 index 0000000..27d699d --- /dev/null +++ b/scripts/stream_regress.sh @@ -0,0 +1,27 @@ +#!/usr/bin/env bash +# Real-clip regression for the streaming de-ghosting toggles (tasks A/B/C). +# Extracts frames from a video (ffmpeg, mirroring the demo server's fps/scale) and +# runs build/tests/stream_regress over the A/B/C variants, printing per-variant +# point count, GT-free surface thickness (ghosting proxy), and wall-clock time. +# +# Needs a pose-capable DA3 model (e.g. models/depth-anything-giant-f32.gguf). +# Build the harness first: +# cmake --build build --target stream_regress +# Usage: +# scripts/stream_regress.sh
@@ -504,6 +513,13 @@

depth-anything.cpp · DA3 studio

document.querySelectorAll('#modeseg button').forEach(b=>b.onclick=()=>{ if(!b.disabled) setMode(b.dataset.mode); }); document.querySelectorAll('#srcseg button').forEach(b=>b.onclick=()=>setSrc(b.dataset.src)); +// Surface-fusion (TSDF) sliders: reveal only when fusion is on; live value readouts. +const fuseCb=$("fuse"); +function syncFuse(){ $("fuseparams").classList.toggle("hide",!fuseCb.checked); } +fuseCb.onchange=syncFuse; syncFuse(); +$("fusedetail").oninput=()=>$("fusedetailval").textContent=(+$("fusedetail").value*100).toFixed(2)+"%"; +$("fusetrunc").oninput=()=>$("fusetruncval").textContent=$("fusetrunc").value+"× voxel"; + // ---- file pick / drop ---- let chosen=null; const drop=$("drop"), fileInput=$("file"); @@ -533,6 +549,8 @@

depth-anything.cpp · DA3 studio

fd.append("icp_refine",$("icprefine").checked?"1":"0"); fd.append("loop_close",$("loopclose").checked?"1":"0"); fd.append("fuse",$("fuse").checked?"1":"0"); + fd.append("fuse_voxel_frac",$("fusedetail").value); + fd.append("fuse_trunc_mult",$("fusetrunc").value); fd.append("metric",$("metric").checked?"1":"0"); fd.append("near_bias",$("nearbias").value); fd.append("surfels",$("surfels").checked?"1":"0"); diff --git a/src/da_capi.cpp b/src/da_capi.cpp index 77ac232..0258ff6 100644 --- a/src/da_capi.cpp +++ b/src/da_capi.cpp @@ -27,6 +27,10 @@ struct da_ctx { // axes, 3F each). da_capi_points_stream is at the purego arg-count ceiling, so these // are retrieved via da_capi_stream_last_poses instead of as out-params. std::vector frame_pos, frame_fwd; + // Scene-relative TSDF fusion knobs for the next da_capi_points_stream (set via + // da_capi_set_fuse_params, which dodges the stream's arg-count ceiling). 0 => the + // fuse_tsdf defaults (voxel = 0.4% of bbox diag, truncation = 4 voxels). + double fuse_voxel_frac = 0.0, fuse_trunc_mult = 0.0; }; static char* dup_cstr(const std::string& s){ @@ -81,7 +85,18 @@ static bool capi_run_nested(da_ctx* c, const char* image_path, } extern "C" { -int da_capi_abi_version(void){ return 9; } +int da_capi_abi_version(void){ return 10; } + +// Scene-relative TSDF fusion knobs applied by the NEXT da_capi_points_stream when +// its fuse flag is set. voxel_frac = voxel edge as a fraction of the bbox diagonal +// ("fusion detail"); trunc_mult = truncation as a multiple of the voxel ("merge +// range", max merged gap = 2*trunc_mult voxels). Either <=0 keeps the built-in +// default. Split out from da_capi_points_stream, which is at the purego arg ceiling. +void da_capi_set_fuse_params(da_ctx* c, double voxel_frac, double trunc_mult){ + if (!c) return; + c->fuse_voxel_frac = voxel_frac; + c->fuse_trunc_mult = trunc_mult; +} da_ctx* da_capi_load(const char* path, int n_threads){ if (!path) return nullptr; auto e = da::Engine::load(path, n_threads); @@ -406,7 +421,9 @@ int da_capi_points_stream(da_ctx* c, const char** image_paths, int n_images, // the truncation band (max merge gap = 2*trunc) defaults to 4 voxels inside. if (fuse != 0){ da::TsdfParams tp; - tp.voxel = (float)fuse_voxel_m; // <=0 => fuse_tsdf derives ~0.4% of bbox diag + tp.voxel = (float)fuse_voxel_m; // absolute override (harness); <=0 => use frac + tp.voxel_frac = (float)c->fuse_voxel_frac; // scene-relative "detail" knob + tp.trunc_mult = (float)c->fuse_trunc_mult; // scene-relative "merge range" knob long pre = 0; for (int cc : sc.counts) pre += cc; size_t pre_pts = sc.radius.size(); auto t_fuse0 = std::chrono::steady_clock::now(); diff --git a/src/tsdf.cpp b/src/tsdf.cpp index 08454b9..d912ba2 100644 --- a/src/tsdf.cpp +++ b/src/tsdf.cpp @@ -53,9 +53,11 @@ int fuse_tsdf(std::vector& xyz, std::vector& rgb, const double diag = std::sqrt((double)(hi[0]-lo[0])*(hi[0]-lo[0]) + (double)(hi[1]-lo[1])*(hi[1]-lo[1]) + (double)(hi[2]-lo[2])*(hi[2]-lo[2])); - double vox = (p.voxel > 0.f) ? p.voxel : (diag > 0 ? 0.004*diag : 0.03); + const double vfrac = (p.voxel_frac > 0.f) ? p.voxel_frac : 0.004; + double vox = (p.voxel > 0.f) ? p.voxel : (diag > 0 ? vfrac*diag : 0.03); if (!(vox > 0)) vox = 0.03; - double trunc = (p.trunc > 0.f) ? p.trunc : 4.0*vox; + const double tmult = (p.trunc_mult > 0.f) ? p.trunc_mult : 4.0; + double trunc = (p.trunc > 0.f) ? p.trunc : tmult*vox; double nr = (p.normal_radius > 0.f) ? p.normal_radius : 2.5*vox; const int min_hits = std::max(1, p.min_hits); diff --git a/src/tsdf.hpp b/src/tsdf.hpp index e39d1b9..d58302b 100644 --- a/src/tsdf.hpp +++ b/src/tsdf.hpp @@ -21,12 +21,20 @@ namespace da { struct TsdfParams { - // Voxel edge (world units). <=0 => derived as 0.004 * bbox-diagonal, because - // monocular reconstructions live at an arbitrary metric scale. + // Voxel edge (world units). <=0 => voxel_frac * bbox-diagonal, because monocular + // reconstructions live at an arbitrary metric scale (an absolute cm-voxel would + // either no-op or collapse the whole cloud). float voxel = 0.f; + // Voxel edge as a FRACTION of the bbox diagonal (used only when voxel<=0). This + // is the scene-relative "fusion detail" knob. <=0 => 0.004 (0.4% of extent). + float voxel_frac = 0.f; // Truncation half-width of the SDF band (world units). Two sheets merge only - // where their bands overlap, i.e. gaps up to 2*trunc collapse. <=0 => 4*voxel. + // where their bands overlap, i.e. gaps up to 2*trunc collapse. <=0 => trunc_mult*voxel. float trunc = 0.f; + // Truncation as a MULTIPLE of the voxel edge (used only when trunc<=0). This is + // the scene-relative "merge range" knob (max merged gap = 2*trunc_mult voxels). + // <=0 => 4. + float trunc_mult = 0.f; // PCA normal-estimation radius (world units). <=0 => 2.5*voxel. float normal_radius = 0.f; // Drop zero-crossing voxels backed by fewer than this many point deposits From 7b1f5e599a997ad38deaa26563d237690fb18d48 Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Sun, 5 Jul 2026 21:08:46 +0100 Subject: [PATCH 16/22] fix(fuse): keep the flythrough reveal frame-major after TSDF fusion MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit TSDF (like the old voxel fusion) re-sorts points spatially and the fuse block only rescaled the per-frame counts proportionally, so the progressive/ flythrough reveal — which reveals a frame-major prefix — showed a spatial voxel slab unrelated to the camera position. During a flythrough large chunks of already-passed geometry read as missing until the slab swept to them. fuse_tsdf now tags each output voxel with the FIRST input frame that observed it (via a per-input-point frame array built from the stream counts), emits frame-major by that frame, and returns the per-point frames; da_capi rebuilds real per-frame counts from them. The reveal now reveals each surface as soon as any camera has seen it and never hides it again, tracking the capture path (clustered by sliding window but monotonic, cumulating to the full cloud). (Nearest-camera assignment was tried and rejected: it back-loads the reveal by hiding already-seen surfaces until the camera's closest approach.) Co-Authored-By: Claude Opus 4.8 --- src/da_capi.cpp | 30 ++++++++++++++++++++++++------ src/tsdf.cpp | 20 +++++++++++++++++--- src/tsdf.hpp | 11 ++++++++++- 3 files changed, 51 insertions(+), 10 deletions(-) diff --git a/src/da_capi.cpp b/src/da_capi.cpp index 0258ff6..d8e1f5c 100644 --- a/src/da_capi.cpp +++ b/src/da_capi.cpp @@ -424,18 +424,36 @@ int da_capi_points_stream(da_ctx* c, const char** image_paths, int n_images, tp.voxel = (float)fuse_voxel_m; // absolute override (harness); <=0 => use frac tp.voxel_frac = (float)c->fuse_voxel_frac; // scene-relative "detail" knob tp.trunc_mult = (float)c->fuse_trunc_mult; // scene-relative "merge range" knob - long pre = 0; for (int cc : sc.counts) pre += cc; - size_t pre_pts = sc.radius.size(); + const int F = (int)sc.counts.size(); + const size_t pre_pts = sc.radius.size(); + // Per-input-point capture frame. The emitted cloud is frame-major, so the + // per-frame counts define the frame of each point index. Handing this to + // fuse_tsdf lets it tag each output voxel with its first-observing frame and + // emit FRAME-MAJOR, so the progressive/flythrough reveal (which reveals a + // frame prefix) still lines up after fusion. + std::vector in_frame; in_frame.reserve(pre_pts); + for (int f = 0; f < F; ++f) for (int k = 0; k < sc.counts[f]; ++k) in_frame.push_back(f); + if (in_frame.size() != pre_pts) in_frame.clear(); // only use if it lines up + std::vector out_frame; auto t_fuse0 = std::chrono::steady_clock::now(); // frame_pos orients normals toward the observing camera (correct SDF sign); // weights default to 1/radius inside (near/high-confidence points dominate). - int nf = da::fuse_tsdf(sc.xyz, sc.rgb, sc.radius, tp, nullptr, &sc.frame_pos); + int nf = da::fuse_tsdf(sc.xyz, sc.rgb, sc.radius, tp, nullptr, &sc.frame_pos, + in_frame.empty() ? nullptr : &in_frame, &out_frame); DA_LOG("stream timing: tsdf(cpu)=%.2fs %zu->%d pts", std::chrono::duration(std::chrono::steady_clock::now()-t_fuse0).count(), pre_pts, nf); - // Fusion reorders points spatially, so per-frame attribution is lost; rescale - // the build-up counts to sum to the fused N (keeps the progressive reveal valid). - if (pre > 0) for (int& cc : sc.counts) cc = (int)((long long)cc * nf / pre); + // Rebuild real per-frame counts from the fused points' first-observing frame + // (the cloud is now frame-major), so the reveal shows each surface as its + // capture frame is reached. Fall back to a proportional rescale if unavailable. + if ((int)out_frame.size() == nf && F > 0){ + std::vector cc(F, 0); + for (int f : out_frame) if (f >= 0 && f < F) cc[f]++; + sc.counts.swap(cc); + } else { + long pre = 0; for (int cc : sc.counts) pre += cc; + if (pre > 0) for (int& cc : sc.counts) cc = (int)((long long)cc * nf / pre); + } } const size_t np = sc.radius.size(); diff --git a/src/tsdf.cpp b/src/tsdf.cpp index d912ba2..abfa279 100644 --- a/src/tsdf.cpp +++ b/src/tsdf.cpp @@ -34,17 +34,21 @@ struct Cell { double cr = 0, cg = 0, cb = 0; // Σ wi*rgb double rad = 0; // Σ wi*radius int hits = 0; + int fmin = INT32_MAX; // first (min) input frame that fed this voxel }; int fuse_tsdf(std::vector& xyz, std::vector& rgb, std::vector& radius, const TsdfParams& p, const std::vector* weights, - const std::vector* cameras) { + const std::vector* cameras, + const std::vector* in_frame, + std::vector* out_frame) { const int N = (int)(xyz.size()/3); if (N < 8) return N; const bool have_rgb = ((int)rgb.size() == 3*N); const bool have_rad = ((int)radius.size() == N); const bool have_w = (weights && (int)weights->size() == N); + const bool have_fr = (in_frame && (int)in_frame->size() == N); // --- scene-relative defaults (arbitrary monocular scale) --- float lo[3] = {xyz[0],xyz[1],xyz[2]}, hi[3] = {xyz[0],xyz[1],xyz[2]}; @@ -121,12 +125,13 @@ int fuse_tsdf(std::vector& xyz, std::vector& rgb, c.nx += w*nx; c.ny += w*ny; c.nz += w*nz; c.cr += w*cr; c.cg += w*cg; c.cb += w*cb; c.rad += w*ri; c.hits++; + if (have_fr) { int f = (*in_frame)[i]; if (f < c.fmin) c.fmin = f; } } } // --- extract the zero-crossing shell: voxels whose centre is within half a cell // of the surface. Move each onto the level set along its averaged normal. --- - struct Out { VKey k; float x,y,z; uint8_t r,g,b; float rad; }; + struct Out { VKey k; int frame; float x,y,z; uint8_t r,g,b; float rad; }; std::vector outs; outs.reserve(field.size()/2 + 16); const double half = 0.5*vox; @@ -143,14 +148,21 @@ int fuse_tsdf(std::vector& xyz, std::vector& rgb, Out o; o.k = kv.first; o.x = (float)(cx - D*nx); o.y = (float)(cy - D*ny); o.z = (float)(cz - D*nz); + // Reveal frame = the FIRST input frame that observed this voxel, so the + // additive flythrough reveals each surface as soon as any camera has seen it + // and never hides it again (nearest-camera instead back-loads the reveal). + o.frame = (c.fmin == INT32_MAX) ? 0 : c.fmin; double iw = 1.0/c.w; auto c8 = [](double v){ double x=std::max(0.0,std::min(255.0,v)); return (uint8_t)std::lround(x); }; o.r = c8(c.cr*iw); o.g = c8(c.cg*iw); o.b = c8(c.cb*iw); o.rad = (float)(c.rad*iw); outs.push_back(o); } - // Deterministic order (hash iteration is unspecified): sort by voxel key. + // Order FRAME-MAJOR (so the flythrough/build-up reveal, which reveals a frame + // prefix, shows points as their first-observing camera is reached), with the + // voxel key as a deterministic tiebreak (hash iteration is unspecified). std::sort(outs.begin(), outs.end(), [](const Out& a, const Out& b){ + if (a.frame != b.frame) return a.frame < b.frame; if (a.k.x != b.k.x) return a.k.x < b.k.x; if (a.k.y != b.k.y) return a.k.y < b.k.y; return a.k.z < b.k.z; @@ -159,10 +171,12 @@ int fuse_tsdf(std::vector& xyz, std::vector& rgb, const int M = (int)outs.size(); if (M == 0) return N; // nothing extracted (degenerate) -> leave input untouched std::vector ox(3*M), orad(M); std::vector orgb(3*M); + if (out_frame) out_frame->resize(M); for (int i = 0; i < M; ++i) { ox[3*i]=outs[i].x; ox[3*i+1]=outs[i].y; ox[3*i+2]=outs[i].z; orgb[3*i]=outs[i].r; orgb[3*i+1]=outs[i].g; orgb[3*i+2]=outs[i].b; orad[i]=outs[i].rad; + if (out_frame) (*out_frame)[i] = outs[i].frame; } xyz.swap(ox); if (have_rgb) rgb.swap(orgb); else rgb.clear(); diff --git a/src/tsdf.hpp b/src/tsdf.hpp index d58302b..2b4aaad 100644 --- a/src/tsdf.hpp +++ b/src/tsdf.hpp @@ -50,9 +50,18 @@ struct TsdfParams { // centres, world frame) orients each normal toward its nearest camera so the SDF // sign is globally consistent — REQUIRED for correct fusion; if empty, normals are // left with PCA's arbitrary sign and fusion degrades to a denoise. +// +// `in_frame` (length N, optional): the input capture frame of each point. When +// given, each output surface voxel inherits the MINIMUM (first-observing) frame of +// the points that fed it, the output is sorted frame-major by that frame, and (if +// `out_frame` is non-null) the per-output-point frame is written there. This keeps +// the progressive/flythrough reveal — which reveals a frame-major prefix — valid +// after fusion; without it the output is voxel-sorted and reveal order is spatial. int fuse_tsdf(std::vector& xyz, std::vector& rgb, std::vector& radius, const TsdfParams& p, const std::vector* weights = nullptr, - const std::vector* cameras = nullptr); + const std::vector* cameras = nullptr, + const std::vector* in_frame = nullptr, + std::vector* out_frame = nullptr); } // namespace da From 4deec354e202bd0529cfc2fef720384c8da9fcad Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Sun, 5 Jul 2026 21:28:02 +0100 Subject: [PATCH 17/22] fix(fuse): size fused splats to the voxel (dark/sparse surface) + linear colour MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The TSDF emits one point per surface voxel but was keeping the averaged input radius (~0.0014 on a scene where voxels sit ~0.011 apart), so splat scale/spacing collapsed from 0.39 (dense cloud) to 0.07 — the splats covered 7% of the gap and the dark background showed through, reading as a dark, sparse surface. Size each fused splat to 0.6*voxel so the one-point-per-cell surface tiles solidly again (scale/spacing back to ~0.34). Also average colour in LINEAR light (sRGB<->linear) rather than gamma-encoded sRGB. Gamma-space averaging biases merged colour toward black at any local contrast; the global mean barely moved here (so it was not the cause of the perceived darkening) but it is the correct way to blend and helps texture edges. Co-Authored-By: Claude Opus 4.8 --- src/tsdf.cpp | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/src/tsdf.cpp b/src/tsdf.cpp index abfa279..e768ca4 100644 --- a/src/tsdf.cpp +++ b/src/tsdf.cpp @@ -26,6 +26,19 @@ static inline VKey key_at(double x, double y, double z, double cell) { return { (int32_t)std::floor(x/cell), (int32_t)std::floor(y/cell), (int32_t)std::floor(z/cell) }; } +// Colours must be averaged in LINEAR light, not in gamma-encoded sRGB: averaging +// sRGB values pulls the mean toward black (the transfer curve is concave), so a +// voxel merging any colour contrast comes out visibly darker. Convert sRGB u8 -> +// linear before accumulating and back on output so merged colour keeps its brightness. +static inline double srgb_to_lin(double v) { + double c = v / 255.0; + return c <= 0.04045 ? c / 12.92 : std::pow((c + 0.055) / 1.055, 2.4); +} +static inline double lin_to_srgb(double c) { + double s = c <= 0.0031308 ? 12.92 * c : 1.055 * std::pow(c, 1.0 / 2.4) - 0.055; + return s * 255.0; +} + // Per-voxel signed-distance accumulator. sdf/nrm/col are confidence-weighted sums; // divide by w at the end. hits counts contributing deposits (speckle cull). struct Cell { @@ -105,9 +118,9 @@ int fuse_tsdf(std::vector& xyz, std::vector& rgb, double w = have_w ? (double)(*weights)[i] : (have_rad ? 1.0/((double)radius[i] + eps) : 1.0); if (!(w > 0)) w = eps; - const double cr = have_rgb ? rgb[3*i] : 200.0; - const double cg = have_rgb ? rgb[3*i+1] : 200.0; - const double cb = have_rgb ? rgb[3*i+2] : 200.0; + const double cr = srgb_to_lin(have_rgb ? rgb[3*i] : 200.0); // accumulate in + const double cg = srgb_to_lin(have_rgb ? rgb[3*i+1] : 200.0); // linear light so + const double cb = srgb_to_lin(have_rgb ? rgb[3*i+2] : 200.0); // merging keeps brightness const double ri = have_rad ? radius[i] : (float)vox; // March along ±normal one voxel at a time; deposit into each distinct voxel // the SIGNED distance from the surface (point) to that voxel's centre. @@ -154,8 +167,11 @@ int fuse_tsdf(std::vector& xyz, std::vector& rgb, o.frame = (c.fmin == INT32_MAX) ? 0 : c.fmin; double iw = 1.0/c.w; auto c8 = [](double v){ double x=std::max(0.0,std::min(255.0,v)); return (uint8_t)std::lround(x); }; - o.r = c8(c.cr*iw); o.g = c8(c.cg*iw); o.b = c8(c.cb*iw); - o.rad = (float)(c.rad*iw); + o.r = c8(lin_to_srgb(c.cr*iw)); o.g = c8(lin_to_srgb(c.cg*iw)); o.b = c8(lin_to_srgb(c.cb*iw)); + // Size each splat to its voxel so the one-point-per-cell surface tiles solidly. + // Keeping the (much smaller) averaged input radius leaves dark gaps between the + // now voxel-spaced points, which reads as a dark, sparse surface. + o.rad = (float)(0.6 * vox); outs.push_back(o); } // Order FRAME-MAJOR (so the flythrough/build-up reveal, which reveals a frame From 0036040ec66fdb3c133b4f46346e1a810fb7abbc Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Sun, 5 Jul 2026 21:42:32 +0100 Subject: [PATCH 18/22] feat(demo): voxel cubes follow the fusion grid; hide unused inputs per mode MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Voxel mode: when Surface fusion (A) is on, render one cube per FUSION voxel (cube cell = the C-side fuse_voxel_frac of the bbox diagonal) instead of re-blocking the fused cloud on the separate Voxel-res grid. The single "Fuse voxel" slider now drives both the fusion detail and the block size, so the two never disagree; Voxel-res is used only when fusion is off. UI: a single updateControls() shows only the inputs a given mode/model/source uses — streaming controls are video-only; window/overlap/conf/near-bias and the de-ghost toggles are hidden in single-image gaussians mode; surfels hide in voxel + gaussians modes; metric shows only for a metric-capable model; Voxel res hides when fusion drives the grid. Wired to mode/src/model/fuse changes. Co-Authored-By: Claude Opus 4.8 --- server/scenes.go | 21 +++++++++++++++++---- server/web/index.html | 44 ++++++++++++++++++++++++++++++------------- 2 files changed, 48 insertions(+), 17 deletions(-) diff --git a/server/scenes.go b/server/scenes.go index e394fb2..b7f2fb5 100644 --- a/server/scenes.go +++ b/server/scenes.go @@ -383,11 +383,24 @@ func (s *server) bakeVideo(name, vpath, jobDir, model, mode string, maxFrames, c asVoxels := mode == "voxels" var cell float64 if asVoxels { - if voxelRes < 1 { - voxelRes = 100 + if fuse { + // Follow the fusion grid: the cloud is already one point per fusion voxel, + // so render one cube per fusion voxel instead of re-blocking on a second + // grid. Use the same scene-relative fraction the C-side TSDF used (see + // da_capi/tsdf fuse_voxel_frac; 0 => its 0.4%-of-bbox-diagonal default). + frac := fuseVoxelFrac + if frac <= 0 { + frac = 0.004 + } + cell = frac * cloudBBoxDiag(cloud, cloud.N) + log.Printf("scene %s: voxel mode follows fusion cell=%.5g (frac=%.4g)", name, cell, frac) + } else { + if voxelRes < 1 { + voxelRes = 100 + } + cell = cloudBBoxDiag(cloud, cloud.N) / voxelRes + log.Printf("scene %s: voxel mode res=%.0f cell=%.5g", name, voxelRes, cell) } - cell = cloudBBoxDiag(cloud, cloud.N) / voxelRes - log.Printf("scene %s: voxel mode res=%.0f cell=%.5g", name, voxelRes, cell) } else { // Additive reveal: overlapping frames re-observe the same surface, so revealing a // new frame otherwise re-blends already-shown regions. Collapse coincident points diff --git a/server/web/index.html b/server/web/index.html index a9bb65f..89029f9 100644 --- a/server/web/index.html +++ b/server/web/index.html @@ -122,13 +122,13 @@

depth-anything.cpp · DA3 studio

-
+
-
-
-
+
+
+
-
+
@@ -139,10 +139,10 @@

depth-anything.cpp · DA3 studio

Loop closure (C) - -
@@ -504,21 +504,39 @@

depth-anything.cpp · DA3 studio

vBtn.disabled=!(curModel&&curModel.pose); // voxels need pose (same as points) const cur=document.querySelector('#modeseg button.on'); if(cur&&cur.disabled) setMode("points"); + updateControls(); } // ---- mode / source segmented ---- let mode="points", src="video"; -function setMode(m){ mode=m; document.querySelectorAll('#modeseg button').forEach(b=>b.classList.toggle("on",b.dataset.mode===m)); } +function setMode(m){ mode=m; document.querySelectorAll('#modeseg button').forEach(b=>b.classList.toggle("on",b.dataset.mode===m)); updateControls(); } function setSrc(sv){ src=sv; document.querySelectorAll('#srcseg button').forEach(b=>b.classList.toggle("on",b.dataset.src===sv)); - $("vidparams").classList.toggle("hide",sv!=="video"); chosen=null; dropText(); refreshGen(); } + chosen=null; dropText(); refreshGen(); updateControls(); } document.querySelectorAll('#modeseg button').forEach(b=>b.onclick=()=>{ if(!b.disabled) setMode(b.dataset.mode); }); document.querySelectorAll('#srcseg button').forEach(b=>b.onclick=()=>setSrc(b.dataset.src)); -// Surface-fusion (TSDF) sliders: reveal only when fusion is on; live value readouts. -const fuseCb=$("fuse"); -function syncFuse(){ $("fuseparams").classList.toggle("hide",!fuseCb.checked); } -fuseCb.onchange=syncFuse; syncFuse(); +// Show only the inputs a given mode/model/source actually uses. Streaming controls +// are video-only; window/overlap/conf/near-bias/de-ghost toggles don't apply to the +// single-image gaussians mode; surfels shape the smooth point surface (cubes ignore +// them); metric needs a metric-capable model; and Voxel res is unused when fusion +// drives the cube grid. +function show(id,on){ const e=$(id); if(e) e.classList.toggle("hide",!on); } +function updateControls(){ + const vid=(src==='video'), gs=(mode==='gaussians'), vox=(mode==='voxels'); + const fuseOn=$('fuse').checked, metricModel=!!(curModel&&curModel.metric); + $('vidparams').classList.toggle('hide',!vid); + $('vidclean').classList.toggle('hide',!vid||gs); + $('fuseparams').classList.toggle('hide',!vid||gs||!fuseOn); + show('fld-window',!gs); show('fld-overlap',!gs); show('fld-conf',!gs); show('fld-nearbias',!gs); + show('fld-voxelres', vox && !fuseOn); // fusion drives the cube grid when on + show('fld-surfels', !gs && !vox); + show('fld-metric', !gs && metricModel); +} + +// Surface-fusion (TSDF) sliders: live value readouts (visibility via updateControls). +$("fuse").onchange=updateControls; $("fusedetail").oninput=()=>$("fusedetailval").textContent=(+$("fusedetail").value*100).toFixed(2)+"%"; $("fusetrunc").oninput=()=>$("fusetruncval").textContent=$("fusetrunc").value+"× voxel"; +updateControls(); // ---- file pick / drop ---- let chosen=null; From 2d48816f604ffa0bf25bb25d86354fd129bf2b7c Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Mon, 6 Jul 2026 08:33:52 +0100 Subject: [PATCH 19/22] feat(ui): runtime brightness slider (overbright to offset cube face-shading) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Voxel cubes are face-lit (0.40 ambient .. 1.0 lit), so away from the light they read much darker than the unshaded points. Add a ☀ brightness slider to the viewer that multiplies the fragment output (u_bright uniform in both the cube and splat programs, default 1.0, range 0.5..3). Cranking it overbrightens the shaded faces to match the points; it doubles as a general exposure control for point/gaussian scenes. Runtime only — the continuous render loop picks up the value each frame, no rebake. Co-Authored-By: Claude Opus 4.8 --- server/web/index.html | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/server/web/index.html b/server/web/index.html index 89029f9..811bdff 100644 --- a/server/web/index.html +++ b/server/web/index.html @@ -96,6 +96,9 @@

depth-anything.cpp · DA3 studio

— +
@@ -216,8 +219,8 @@

depth-anything.cpp · DA3 studio

v_pos = a_corner; gl_Position = vec4(ndc+offset, clip.z/clip.w, 1.0); }`; const FRAG = `#version 300 es -precision highp float; in vec4 v_color; in vec2 v_pos; out vec4 o; -void main(){ float p=-dot(v_pos,v_pos); if(p<-4.0) discard; float a=exp(p)*v_color.a; o=vec4(v_color.rgb*a,a); }`; +precision highp float; uniform float u_bright; in vec4 v_color; in vec2 v_pos; out vec4 o; +void main(){ float p=-dot(v_pos,v_pos); if(p<-4.0) discard; float a=exp(p)*v_color.a; o=vec4(v_color.rgb*a*u_bright,a); }`; function shader(t,src){ const s=gl.createShader(t); gl.shaderSource(s,src); gl.compileShader(s); if(!gl.getShaderParameter(s,gl.COMPILE_STATUS)) throw new Error(gl.getShaderInfoLog(s)); return s; } @@ -232,7 +235,8 @@

depth-anything.cpp · DA3 studio

const idxBuf = gl.createBuffer(); const aIndex = gl.getAttribLocation(prog,"a_index"); const tex = gl.createTexture(); const TEX_W = 2048; const U = { proj:gl.getUniformLocation(prog,"u_proj"), view:gl.getUniformLocation(prog,"u_view"), - focal:gl.getUniformLocation(prog,"u_focal"), viewport:gl.getUniformLocation(prog,"u_viewport") }; + focal:gl.getUniformLocation(prog,"u_focal"), viewport:gl.getUniformLocation(prog,"u_viewport"), + bright:gl.getUniformLocation(prog,"u_bright") }; gl.uniform1i(gl.getUniformLocation(prog,"u_tex"),0); gl.disable(gl.DEPTH_TEST); gl.enable(gl.BLEND); gl.blendFuncSeparate(gl.ONE,gl.ONE_MINUS_SRC_ALPHA,gl.ONE,gl.ONE_MINUS_SRC_ALPHA); @@ -259,14 +263,18 @@

depth-anything.cpp · DA3 studio

v_shade = 0.40 + 0.60*max(dot(a_cnrm,L),0.0); }`; const CFRAG=`#version 300 es -precision highp float; in vec3 v_rgb; in float v_shade; out vec4 o; -void main(){ o = vec4(v_rgb*v_shade, 1.0); }`; +precision highp float; uniform float u_bright; in vec3 v_rgb; in float v_shade; out vec4 o; +void main(){ o = vec4(v_rgb*v_shade*u_bright, 1.0); }`; const cubeProg = gl.createProgram(); gl.attachShader(cubeProg, shader(gl.VERTEX_SHADER,CVERT)); gl.attachShader(cubeProg, shader(gl.FRAGMENT_SHADER,CFRAG)); gl.linkProgram(cubeProg); if(!gl.getProgramParameter(cubeProg,gl.LINK_STATUS)){ setStatus(gl.getProgramInfoLog(cubeProg),true); throw new Error("cube link"); } gl.useProgram(cubeProg); gl.uniform1i(gl.getUniformLocation(cubeProg,"u_tex"),0); -const cubeU = { proj:gl.getUniformLocation(cubeProg,"u_proj"), view:gl.getUniformLocation(cubeProg,"u_view") }; +const cubeU = { proj:gl.getUniformLocation(cubeProg,"u_proj"), view:gl.getUniformLocation(cubeProg,"u_view"), + bright:gl.getUniformLocation(cubeProg,"u_bright") }; +// Viewer brightness/exposure (runtime): compensates for cube face-shading in voxel +// mode (and general exposure for splats). 1 = as-authored; >1 overbrightens. +let viewBright = 1.0; const aCpos = gl.getAttribLocation(cubeProg,"a_cpos"), aCnrm = gl.getAttribLocation(cubeProg,"a_cnrm"); // 36 verts (6 faces × 2 tris), each = pos(±0.5) + outward normal. No culling, so winding is irrelevant. const cubeGeom = (()=>{ const V=[[-.5,-.5,-.5],[.5,-.5,-.5],[.5,.5,-.5],[-.5,.5,-.5],[-.5,-.5,.5],[.5,-.5,.5],[.5,.5,.5],[-.5,.5,.5]]; @@ -435,6 +443,7 @@

depth-anything.cpp · DA3 studio

// Opaque, depth-tested cubes — no back-to-front sort needed. gl.enable(gl.DEPTH_TEST); gl.disable(gl.BLEND); gl.useProgram(cubeProg); gl.uniformMatrix4fv(cubeU.view,false,view); gl.uniformMatrix4fv(cubeU.proj,false,proj); + gl.uniform1f(cubeU.bright, viewBright); if(aCorner>=0) gl.disableVertexAttribArray(aCorner); if(aIndex>=0) gl.disableVertexAttribArray(aIndex); gl.bindBuffer(gl.ARRAY_BUFFER,cubeBuf); @@ -446,6 +455,7 @@

depth-anything.cpp · DA3 studio

gl.disable(gl.DEPTH_TEST); gl.enable(gl.BLEND); gl.useProgram(prog); gl.uniformMatrix4fv(U.view,false,view); gl.uniformMatrix4fv(U.proj,false,proj); gl.uniform2f(U.focal,fx,fy); gl.uniform2f(U.viewport,canvas.width,canvas.height); + gl.uniform1f(U.bright, viewBright); if(aCpos>=0) gl.disableVertexAttribArray(aCpos); if(aCnrm>=0) gl.disableVertexAttribArray(aCnrm); gl.bindBuffer(gl.ARRAY_BUFFER,quadBuf); gl.enableVertexAttribArray(aCorner); gl.vertexAttribPointer(aCorner,2,gl.FLOAT,false,0,0); gl.vertexAttribDivisor(aCorner,0); @@ -644,6 +654,7 @@

depth-anything.cpp · DA3 studio

revealTarget = (st.npts!=null&&st.npts>0) ? Math.min(st.npts,count) : count; // same file: ease the reveal } } +$("brightness").oninput=()=>{ viewBright=+$("brightness").value; }; // runtime; continuous render picks it up $("steprange").oninput=()=>{ stopPlay(); showStep(+$("steprange").value,false); }; $("playbtn").onclick=()=>{ if(playTimer) stopPlay(); else playFrom(stepI>=steps.length-1?0:stepI); }; $("flybtn").onclick=()=>{ if($("flybtn").disabled) return; From 4c36138898b111958aad6b532a629f387c11922e Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Mon, 6 Jul 2026 08:36:24 +0100 Subject: [PATCH 20/22] fix(ui): stop the play/brightness row overflowing horizontally The brightness slider pushed the control row past its box: #steprange is flex:1 but range inputs default to min-width:auto, so it wouldn't shrink to make room. Give it min-width:60px and let .play wrap, so the brightness control drops to a second line on narrow panels instead of overflowing. Co-Authored-By: Claude Opus 4.8 --- server/web/index.html | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/server/web/index.html b/server/web/index.html index 811bdff..5d21fcf 100644 --- a/server/web/index.html +++ b/server/web/index.html @@ -66,11 +66,11 @@ .stat span:first-child { color: #8b93a1; } #status { margin-top: 8px; color: #9fd3a0; font-size: 11.5px; min-height: 1.2em; } #status.err { color: #ff8a8a; } - .play { display: flex; align-items: center; gap: 8px; margin-top: 8px; } + .play { display: flex; flex-wrap: wrap; align-items: center; gap: 8px; margin-top: 8px; } .play button { padding: 5px 10px; border-radius: 7px; cursor: pointer; background: rgba(255,255,255,.07); border: 1px solid rgba(255,255,255,.12); color: #cdd3dc; } .play button.on { background: rgba(90,140,255,.25); color: #cfe0ff; border-color: rgba(90,140,255,.45); } .play button:disabled { opacity: .35; cursor: not-allowed; } - #steprange { flex: 1; } + #steprange { flex: 1; min-width: 60px; } .hide { display: none !important; } #hint { position: fixed; bottom: 14px; left: 50%; transform: translateX(-50%); z-index: 10; color: #6b7280; font-size: 11.5px; background: rgba(16,19,24,.7); padding: 6px 12px; border-radius: 8px; } From 906aea6c057c3349bfe7bb727cce78cb0cec9151 Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Mon, 6 Jul 2026 08:43:01 +0100 Subject: [PATCH 21/22] fix(ui): brightness on its own line; shorten frame readout to the number MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The step readout printed the full label ("60 · full 60-view cloud · 81147 pts") on the last step, which crowded the sliders and reflowed the brightness control onto the frame row as the label appeared/disappeared. Show just the frame number, and give the brightness control flex-basis:100% so it always sits on its own line with a flex:1 slider that fills the width. Co-Authored-By: Claude Opus 4.8 --- server/web/index.html | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/server/web/index.html b/server/web/index.html index 5d21fcf..4ffa547 100644 --- a/server/web/index.html +++ b/server/web/index.html @@ -97,8 +97,8 @@

depth-anything.cpp · DA3 studio

— + style="flex-basis:100%;display:flex;align-items:center;gap:6px;color:#8b93a1;font-size:11px"> + ☀ @@ -642,7 +642,7 @@

depth-anything.cpp · DA3 studio

await showStep(stepI,true); } async function showStep(i,reframe){ stepI=i; $("steprange").value=i; const st=steps[i]; - $("stepnum").textContent=`${st.n}${st.label?" · "+st.label:""}`; + $("stepnum").textContent=`${st.n}`; // just the frame number; keeps the row short + stable const url=`/scenes-assets/${curScene}/${st.splat}?v=${curBust}`; if(url!==curSplatURL){ // New scenes share ONE cloud file across all steps -> load once, then only the reveal From f56e9be43a22c12ef575584d2fa57a6a5d5be7ae Mon Sep 17 00:00:00 2001 From: Richard Palethorpe Date: Tue, 7 Jul 2026 16:09:36 +0100 Subject: [PATCH 22/22] feat(diag): log flash-attn placement + backbone shape to explain VRAM Two one-time stderr diagnostics so a VRAM/latency complaint can be attributed without guesswork: - backend.cpp: reports whether the fused flash-attention nodes actually run on the GPU or get silently offloaded to CPU (Vulkan gates flash on coopmat2 / subgroup shuffle+vote + head_dim%8; when unmet, supports_op returns false and the scheduler runs attention on CPU, materializing the O(N^2) scores there). Uses the same supports_op the scheduler uses, so it reports the real placement. Guarded to log once per backend. - dino_backbone.cpp: reports resolved embed/heads/head_dim/depth (=model size) and the img/patch -> Ntok/view, S views -> Ntok*S global-attn token count, which is what per-frame activation memory scales with. Re-fires only when a larger window arrives (the VRAM high-water mark), so streaming's many windows don't spam. Validated on Vulkan (RTX 5070 Ti, NV_coopmat2): flash runs on-device, da3-base resolves to embed=768/heads=12/depth=12, 504x280 -> Ntok=721, S=6 -> 4326 tokens. No behavior change; pure DA_LOG output. Co-Authored-By: Claude Opus 4.8 --- src/backend.cpp | 30 ++++++++++++++++++++++++++++++ src/dino_backbone.cpp | 11 +++++++++++ src/dino_backbone.hpp | 3 +++ 3 files changed, 44 insertions(+) diff --git a/src/backend.cpp b/src/backend.cpp index 69eb5bb..801def1 100644 --- a/src/backend.cpp +++ b/src/backend.cpp @@ -40,6 +40,7 @@ struct Backend::Impl { ggml_backend_sched_t sched = nullptr; // GPU path: schedules over {backend, cpu_backend} size_t sched_nodes = 0; // graph_size the sched was created with (recreated when a bigger graph arrives) bool use_sched = false; // true only when `backend` is a GPU device + bool flash_diag_logged = false; // one-time flash-attn placement diagnostic (GPU path) // Inputs registered by the build lambda for the in-flight compute. Copied // into the allocated tensors after the graph is allocated, then cleared. // Never overlaps across calls (compute is not re-entrant). @@ -280,6 +281,35 @@ bool Backend::compute(const std::function& build, } } + // One-time diagnostic (GPU path): does the fused flash-attention node actually + // run on THIS device, or does the scheduler silently offload it to CPU? On + // Vulkan the flash kernel is gated on coopmat2 / subgroup shuffle+vote and the + // head-dim being a multiple of 8; when the device lacks those, supports_op() + // returns false and attention falls back to CPU (correct output, but the O(N^2) + // score matrix materializes CPU-side and the VRAM/latency win is lost). This + // fallback is otherwise invisible, so surface it once per backend instance. + if (impl_->use_sched && !impl_->flash_diag_logged) { + const int n_nodes = ggml_graph_n_nodes(gf); + int n_flash = 0, n_flash_off = 0; + for (int i = 0; i < n_nodes; ++i) { + ggml_tensor* node = ggml_graph_node(gf, i); + if (node->op != GGML_OP_FLASH_ATTN_EXT) continue; + ++n_flash; + if (!ggml_backend_supports_op(impl_->backend, node)) ++n_flash_off; + } + if (n_flash > 0) { + impl_->flash_diag_logged = true; + if (n_flash_off == 0) + DA_LOG("flash-attn: %d node(s) run on %s (fused, no N^2 score matrix)", + n_flash, device_name_.c_str()); + else + DA_LOG("flash-attn: %d/%d node(s) UNSUPPORTED on %s -> offloaded to CPU " + "(scores materialize CPU-side; check coopmat2 / subgroup " + "shuffle+vote support, or run DA_ATTN=manual)", + n_flash_off, n_flash, device_name_.c_str()); + } + } + bool alloc_ok = false; if (need_sched) { // GPU path: schedule across {GPU, CPU}. Unsupported ops fall back to CPU. diff --git a/src/dino_backbone.cpp b/src/dino_backbone.cpp index 2511b4f..0cad35e 100644 --- a/src/dino_backbone.cpp +++ b/src/dino_backbone.cpp @@ -465,6 +465,17 @@ bool DinoBackbone::forward_mv_ordered(const std::vector>& vie const int S=(int)views_chw.size(); const float eps=c.ln_eps; + // One-time shape diagnostic (re-fires only for a larger window, which sets the + // VRAM high-water mark): the cross-view/global attention runs over Ntok*S tokens + // in one pass, so per-frame activation memory scales with this product. Pairs + // with backend.cpp's flash-attn placement line to explain the frame ceiling. + if (S > diag_max_s_) { + diag_max_s_ = S; + DA_LOG("backbone: embed=%d heads=%d head_dim=%d depth=%d | %dx%d img, patch=%d " + "-> Ntok=%d/view, S=%d views -> global attn over %d tokens", + embed, heads, hd, (int)c.depth, W, H, patch, Ntok, S, Ntok*S); + } + // Per-view RoPE position sets (identical across views: same patch grid). std::vector pos_local(2*Ntok, 0.f), pos_nodiff(2*Ntok, 0.f); for (int t=1;t>& cls, int embed) const; ModelLoader& ml_; Backend& be_; + // Largest view-count logged by the one-time shape diagnostic so far; the log + // re-fires only when a bigger window arrives (that window sets the VRAM peak). + int diag_max_s_ = 0; }; }