Pointmap: switch from spheres to surfels (oriented quads aligned to surface normal)

app.py

@@ -236,62 +236,86 @@ def _triangulate_grid(pointmap_hwc: np.ndarray, mask_hw: np.ndarray,
 
 
 def _make_glb(image_pil_texture: Image.Image, pointmap_hwc: np.ndarray,
-              mask_hw: np.ndarray, max_points: int = 25_000) -> str:
-    """Render the pointmap as a cloud of tiny coloured spheres (one per pixel).
+              mask_hw: np.ndarray, max_points: int = 80_000) -> str:
+    """Render the pointmap as **surfels** — each pixel becomes a small flat
+    quad oriented along the surface normal at that point. No triangulation
+    BETWEEN points (so no triangle artifacts), but each surfel still tiles
+    the surface, giving a smooth "skin" appearance instead of granular balls.
 
-    Avoids triangulation artifacts entirely — there's no surface mesh between
-    pixels. Each pixel becomes an icosphere instance colored by its image
-    pixel, radius scaled to the typical inter-pixel distance in 3D so the cloud
-    looks coherent without huge overlaps.
-
-    `image_pil_texture` is sampled (downscaled) to native pixel positions for
-    color picking — high-res isn't useful here since we sample point-wise.
+    4 verts + 2 faces per surfel → much lighter than spheres (~1/10 the verts).
     """
     H, W = pointmap_hwc.shape[:2]
     z = pointmap_hwc[:, :, 2]
     valid = mask_hw & np.isfinite(pointmap_hwc).all(axis=2) & (z > 0.05) & (z < 25.0)
-    yy, xx = np.where(valid)
 
-    rng = np.random.default_rng(0)
+    # Smooth per-pixel normals from the pointmap's spatial gradient.
+    px = np.zeros_like(pointmap_hwc, dtype=np.float32)
+    py = np.zeros_like(pointmap_hwc, dtype=np.float32)
+    px[:, 1:-1] = (pointmap_hwc[:, 2:] - pointmap_hwc[:, :-2]) * 0.5
+    py[1:-1, :] = (pointmap_hwc[2:, :] - pointmap_hwc[:-2, :]) * 0.5
+    n_grid = np.cross(px, py)
+    n_grid /= np.linalg.norm(n_grid, axis=2, keepdims=True).clip(min=1e-8)
+
+    yy, xx = np.where(valid)
     if len(yy) > max_points:
-        idx = rng.choice(len(yy), max_points, replace=False)
+        idx = np.random.default_rng(0).choice(len(yy), max_points, replace=False)
         yy, xx = yy[idx], xx[idx]
 
-    pts = pointmap_hwc[yy, xx].astype(np.float32)                                  # (N, 3)
-    image_native = np.asarray(image_pil_texture.resize((W, H), Image.LANCZOS))      # (H, W, 3)
-    cols = image_native[yy, xx].astype(np.uint8)                                    # (N, 3)
+    pts = pointmap_hwc[yy, xx].astype(np.float32)
+    nrm = n_grid[yy, xx].astype(np.float32)
+    image_native = np.asarray(image_pil_texture.resize((W, H), Image.LANCZOS))
+    cols = image_native[yy, xx].astype(np.uint8)
 
-    # Y-up flip + recenter on centroid for a tight default view.
     flip = np.array([1.0, -1.0, -1.0], dtype=np.float32)
     pts = pts * flip
+    nrm = nrm * flip
     centroid = pts.mean(axis=0).astype(np.float32) if len(pts) else np.zeros(3, np.float32)
     pts = pts - centroid
 
-    # Sphere radius from MEASURED median nearest-neighbour distance.
-    # Why measured (vs. extent / cbrt(N)): our points sit on a ~2D surface
-    # (skin), not in a 3D volume — a kd-tree gives the right spacing directly.
-    # Factor 0.6 → adjacent balls overlap slightly, forming a solid surface
-    # without big gaps and without obvious overlap blobs.
+    # Surfel size from KD-tree median nearest-neighbour distance.
+    # Factor 0.6 → adjacent surfels overlap slightly along the surface so
+    # there are no visible gaps between tiles.
     if len(pts) >= 2:
         from scipy.spatial import cKDTree
-        nn_dist = cKDTree(pts).query(pts, k=2)[0][:, 1]   # k=1 is self
-        radius = max(float(np.median(nn_dist)) * 0.6, 1e-4)
+        nn_dist = cKDTree(pts).query(pts, k=2)[0][:, 1]
+        r = max(float(np.median(nn_dist)) * 0.6, 1e-4)
     else:
-        radius = 0.005  # fallback for empty/tiny clouds
-
-    # Icosphere template (subdivisions=1 → 42 verts, 80 faces — smoother balls).
-    sphere = trimesh.creation.icosphere(subdivisions=1, radius=radius)
-    sv = sphere.vertices.astype(np.float32)         # (V=42, 3)
-    sf = sphere.faces.astype(np.int64)              # (F=80, 3)
-    V, F = len(sv), len(sf)
+        r = 0.005
+
+    # Tangent basis (u, v) per surfel, perpendicular to that surfel's normal.
+    up = np.array([0.0, 1.0, 0.0], dtype=np.float32)
+    u = np.cross(nrm, up)                                      # (N, 3)
+    u_norm = np.linalg.norm(u, axis=1, keepdims=True)
+    # Fallback for any normals nearly parallel to "up".
+    parallel = (u_norm.squeeze(-1) < 1e-4)
+    if parallel.any():
+        ref = np.array([1.0, 0.0, 0.0], dtype=np.float32)
+        u[parallel] = np.cross(nrm[parallel], ref)
+        u_norm[parallel] = np.linalg.norm(u[parallel], axis=1, keepdims=True).clip(min=1e-8)
+    u = u / u_norm.clip(min=1e-8)
+    v = np.cross(nrm, u)                                       # (N, 3) — already unit
+
+    # 4 corners per surfel (CCW when viewed from +normal):
+    #   0 = P − r*u − r*v
+    #   1 = P + r*u − r*v
+    #   2 = P + r*u + r*v
+    #   3 = P − r*u + r*v
+    ru, rv = r * u, r * v
+    c0 = pts - ru - rv
+    c1 = pts + ru - rv
+    c2 = pts + ru + rv
+    c3 = pts - ru + rv
+    verts = np.stack([c0, c1, c2, c3], axis=1).reshape(-1, 3).astype(np.float32)
+
+    # 2 triangles per surfel: (0,1,2) and (0,2,3), offset by 4 per surfel.
     N = len(pts)
+    base = (np.arange(N, dtype=np.int64) * 4)[:, None, None]
+    tri = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int64)[None, :, :]
+    faces = (base + tri).reshape(-1, 3)
 
-    # Vectorized instancing.
-    verts = (sv[None, :, :] + pts[:, None, :]).reshape(-1, 3)              # (N*V, 3)
-    face_off = (np.arange(N, dtype=np.int64) * V)[:, None, None]
-    faces = (sf[None, :, :] + face_off).reshape(-1, 3)                      # (N*F, 3)
-    vc = np.empty((N, V, 4), dtype=np.uint8)
-    vc[..., :3] = cols[:, None, :]                                          # broadcast color to all 12 verts
+    # Uniform colour across each surfel's 4 verts.
+    vc = np.empty((N, 4, 4), dtype=np.uint8)
+    vc[..., :3] = cols[:, None, :]
     vc[..., 3] = 255
     vertex_colors = vc.reshape(-1, 4)
 
@@ -300,6 +324,10 @@ def _make_glb(image_pil_texture: Image.Image, pointmap_hwc: np.ndarray,
         vertex_colors=vertex_colors,
         process=True,
     )
+    # Force double-sided so surfels are visible regardless of normal flip after the Y-flip.
+    if hasattr(mesh.visual, "material") and hasattr(mesh.visual.material, "doubleSided"):
+        mesh.visual.material.doubleSided = True
+
     out_path = tempfile.NamedTemporaryFile(delete=False, suffix=".glb").name
     mesh.export(out_path)
     return out_path