Cut panorama overlaps along seams instead of averaging them

The merge weighted every overlap pixel by its distance from each frame's
edge, a 200 px linear cross-fade. Anything the frames disagreed on —
parallax in the near foreground, grass in the wind, a walker — came out
twice at half strength: a soft double edge at 1:1.

dr_pano::seam picks, per output texel at proxy resolution, which frame a
pixel comes from. Where a new frame overlaps the composite the cost is the
gain-corrected difference plus local detail plus nearness to either
frame's edge, taken as the worst over a small window, and the cut is a
dynamic-programming path across the overlap. merge.wgsl weights each frame
by its tent-filtered share of that map, a 64 px blend that follows the
seam, with the edge feather kept as the fallback. The page's preview uses
the same map, and examples/merge.rs takes --feather-only for comparison.

core/dr-gpu/src/merge.rs

@@ -30,18 +30,27 @@
 //! are the caller's to provide and cache — `source` is asked for frame `k`
 //! as it is needed, and a caller short of memory may demosaic on demand.
 //!
+//! # The blend
+//!
+//! With a seam map (`dr_pano::seam`), a frame's weight at a pixel is its
+//! share of the map about that pixel — whole on its own side of a seam,
+//! nothing on the other, and a ramp across a window `seam_blend` pixels
+//! wide that follows the seam. Without one, or where the map has nothing
+//! to say, the weight is the distance to the frame's edge over `feather`,
+//! which hides exposure steps and does not hide parallax: the average draws
+//! anything the frames disagree on twice.
+//!
 //! # What is not here yet
 //!
-//! A feathered blend, not seams and a Laplacian pyramid: the weight is the
-//! distance to the frame's edge, which hides exposure steps and small
-//! misalignments and does not hide parallax. Gain is a scalar per frame
-//! the caller supplies. Both are panorama.md §10's step 5, after the path
-//! writes a file end to end.
+//! A Laplacian pyramid, which would let the seam's blend be narrow for
+//! detail and wide for exposure at once. Gain is a scalar per frame the
+//! caller supplies.
 
 use std::sync::Arc;
 
 use dr_pano::bundle::Cameras;
 use dr_pano::projection::{Bounds, Projection};
+use dr_pano::seam::SeamMap;
 use wgpu::util::DeviceExt;
 
 use crate::readback::await_mapping;
@@ -58,7 +67,7 @@ pub struct MergeFrame {
 }
 
 /// The output the merge produces.
-#[derive(Debug, Clone, Copy, PartialEq)]
+#[derive(Debug, Clone, PartialEq)]
 pub struct MergeOutput {
     pub projection: Projection,
     /// The projection's scale in output pixels: the cylinder's radius, the
@@ -69,6 +78,11 @@ pub struct MergeOutput {
     pub bounds: Bounds,
     /// Pixels over which a frame's weight ramps up from its edge.
     pub feather: f32,
+    /// Which frame each part of the output is taken from, laid out at the
+    /// proxies' scale; `None` for the feathered average everywhere.
+    pub seams: Option<Arc<SeamMap>>,
+    /// The width, in output pixels, of the blend across a seam.
+    pub seam_blend: f32,
     /// Chunk size: the unit of GPU work and of memory.
     pub chunk: (u32, u32),
     /// Multiplies a normalised sample (1.0 = white) to the sensor's scale.
@@ -115,6 +129,12 @@ struct WarpParams {
     feather: f32,
     clip_onset: f32,
     balance: [f32; 4],
+    seam_origin: [f32; 2],
+    seam_size: [u32; 2],
+    seam_px: f32,
+    seam_radius: f32,
+    frame_index: u32,
+    seam_on: u32,
 }
 
 #[repr(C)]
@@ -182,6 +202,16 @@ impl MergePass {
                         count: None,
                     },
                     storage(2, false),
+                    wgpu::BindGroupLayoutEntry {
+                        binding: 3,
+                        visibility: wgpu::ShaderStages::COMPUTE,
+                        ty: wgpu::BindingType::Texture {
+                            sample_type: wgpu::TextureSampleType::Uint,
+                            view_dimension: wgpu::TextureViewDimension::D2,
+                            multisampled: false,
+                        },
+                        count: None,
+                    },
                 ],
             });
         let resolve_layout =
@@ -279,6 +309,32 @@ impl MergePass {
         let mut band_cov = vec![false; (out_w * ch) as usize];
         let mut chunk_px: Vec<u32> = Vec::new();
 
+        // The seam map, once for the whole output, and where it sits in
+        // this output's coordinates. A one-texel stand-in when there is
+        // none, because the binding is not optional.
+        let (seam_tex, seam_origin, seam_px, seam_radius, seam_size) = match &output.seams {
+            Some(m) => {
+                let ((ou, ov), px) = m.at_scale(output.scale);
+                let radius = m.blend_radius(output.scale, f64::from(output.seam_blend));
+                (
+                    self.label_texture(m.width as u32, m.height as u32, &m.labels),
+                    [ou as f32, ov as f32],
+                    px as f32,
+                    radius as f32,
+                    [m.width as u32, m.height as u32],
+                )
+            }
+            None => (
+                self.label_texture(1, 1, &[dr_pano::seam::NONE]),
+                [0.0; 2],
+                1.0,
+                1.0,
+                [1, 1],
+            ),
+        };
+        let seam_view = seam_tex.create_view(&Default::default());
+        let seam_on = u32::from(output.seams.is_some());
+
         let mut y = 0u32;
         while y < out_h {
             let rows = ch.min(out_h - y);
@@ -346,8 +402,14 @@ impl MergePass {
                             output.balance[2].max(1e-3),
                             0.0,
                         ],
+                        seam_origin,
+                        seam_size,
+                        seam_px,
+                        seam_radius,
+                        frame_index: k as u32,
+                        seam_on,
                     };
-                    self.accumulate(&params, tile);
+                    self.accumulate(&params, tile, &seam_view);
                 }
 
                 self.resolve_chunk((cols, rows), output.sample_scale, &mut chunk_px)?;
@@ -385,7 +447,42 @@ impl MergePass {
         self.ctx.queue.submit(Some(enc.finish()));
     }
 
-    fn accumulate(&mut self, params: &WarpParams, tile: &wgpu::Texture) {
+    /// The seam map's labels as an `r8uint` texture.
+    fn label_texture(&self, width: u32, height: u32, labels: &[u8]) -> wgpu::Texture {
+        let size = wgpu::Extent3d {
+            width,
+            height,
+            depth_or_array_layers: 1,
+        };
+        let tex = self.ctx.device.create_texture(&wgpu::TextureDescriptor {
+            label: Some("merge-seams"),
+            size,
+            mip_level_count: 1,
+            sample_count: 1,
+            dimension: wgpu::TextureDimension::D2,
+            format: wgpu::TextureFormat::R8Uint,
+            usage: wgpu::TextureUsages::TEXTURE_BINDING | wgpu::TextureUsages::COPY_DST,
+            view_formats: &[],
+        });
+        self.ctx.queue.write_texture(
+            wgpu::TexelCopyTextureInfo {
+                texture: &tex,
+                mip_level: 0,
+                origin: wgpu::Origin3d::ZERO,
+                aspect: wgpu::TextureAspect::All,
+            },
+            labels,
+            wgpu::TexelCopyBufferLayout {
+                offset: 0,
+                bytes_per_row: Some(width),
+                rows_per_image: Some(height),
+            },
+            size,
+        );
+        tex
+    }
+
+    fn accumulate(&mut self, params: &WarpParams, tile: &wgpu::Texture, seams: &wgpu::TextureView) {
         let chunk = (params.chunk_size[0], params.chunk_size[1]);
         let uniforms = self
             .ctx
@@ -417,6 +514,10 @@ impl MergePass {
                         binding: 2,
                         resource: acc.as_entire_binding(),
                     },
+                    wgpu::BindGroupEntry {
+                        binding: 3,
+                        resource: wgpu::BindingResource::TextureView(seams),
+                    },
                 ],
             });
         let mut enc = self.ctx.device.create_command_encoder(&Default::default());

core/dr-gpu/src/shaders/merge.wgsl

@@ -5,8 +5,9 @@
 // pixel it asks which direction that pixel looks along, turns the
 // direction into the frame's camera, projects it to a source pixel, and
 // if that pixel is inside the tile that was rendered for this chunk,
-// samples it and adds it — weighted by its distance from the frame's edge
-// — into the accumulator. `resolve` runs once per chunk after every frame
+// samples it and adds it — weighted by the frame's share of the seam map
+// there, or by its distance from the frame's edge where there is no map —
+// into the accumulator. `resolve` runs once per chunk after every frame
 // has been added: divides the sums by the weights and packs the result as
 // sixteen-bit samples at the sensor's scale (FR-MRG-3).
 //
@@ -50,12 +51,83 @@ struct Params {
     // white balance the composite will be developed with.
     clip_onset: f32,
     balance: vec4<f32>,
+    // The seam map (`dr_pano::seam`): where its texel (0, 0)'s corner sits
+    // in this output's centred coordinates, its size, output pixels per
+    // texel, the blend's radius in texels, which frame this dispatch is,
+    // and whether there is a map at all.
+    seam_origin: vec2<f32>,
+    seam_size: vec2<u32>,
+    seam_px: f32,
+    seam_radius: f32,
+    frame_index: u32,
+    seam_on: u32,
 };
 
 @group(0) @binding(0) var<uniform> p: Params;
 @group(0) @binding(1) var tile: texture_2d<f32>;
 // rgb·w summed, then w: four floats per chunk pixel.
 @group(0) @binding(2) var<storage, read_write> acc: array<vec4<f32>>;
+// One frame index per texel, 255 for none.
+@group(0) @binding(3) var seams: texture_2d<u32>;
+
+const NO_FRAME: u32 = 255u;
+
+fn label(i: i32, j: i32) -> u32 {
+    if (i < 0 || j < 0 || i >= i32(p.seam_size.x) || j >= i32(p.seam_size.y)) {
+        return NO_FRAME;
+    }
+    return textureLoad(seams, vec2<i32>(i, j), 0).r;
+}
+
+// This frame's share of the seam map about output point (u, v): the
+// tent-weighted fraction of the texels within the radius that it owns, and
+// the weight of the texels owned by anyone (zero where the map has nothing
+// to say). `SeamMap::share` verbatim.
+fn seam_share(u: f32, v: f32) -> vec2<f32> {
+    let x = (u - p.seam_origin.x) / p.seam_px - 0.5;
+    let y = (v - p.seam_origin.y) / p.seam_px - 0.5;
+    let r = max(p.seam_radius, 1.0);
+    let x0 = i32(ceil(x - r));
+    let x1 = i32(floor(x + r));
+    let y0 = i32(ceil(y - r));
+    let y1 = i32(floor(y + r));
+    // Most pixels are nowhere near a seam: if the window's corners, edge
+    // midpoints and centre agree, so does the window. A seam crossing it
+    // has to cross its border, between two of those.
+    let xm = i32(round(x));
+    let ym = i32(round(y));
+    let c = label(xm, ym);
+    if (label(x0, y0) == c && label(x1, y0) == c && label(x0, y1) == c && label(x1, y1) == c
+        && label(xm, y0) == c && label(xm, y1) == c && label(x0, ym) == c && label(x1, ym) == c) {
+        if (c == NO_FRAME) {
+            return vec2<f32>(0.0, 0.0);
+        }
+        return vec2<f32>(select(0.0, 1.0, c == p.frame_index), 1.0);
+    }
+    var mine = 0.0;
+    var owned = 0.0;
+    for (var j = y0; j <= y1; j = j + 1) {
+        let wy = 1.0 - abs(y - f32(j)) / r;
+        if (wy <= 0.0) {
+            continue;
+        }
+        for (var i = x0; i <= x1; i = i + 1) {
+            let wx = 1.0 - abs(x - f32(i)) / r;
+            let l = label(i, j);
+            if (wx <= 0.0 || l == NO_FRAME) {
+                continue;
+            }
+            owned = owned + wx * wy;
+            if (l == p.frame_index) {
+                mine = mine + wx * wy;
+            }
+        }
+    }
+    if (owned <= 0.0) {
+        return vec2<f32>(0.0, 0.0);
+    }
+    return vec2<f32>(mine / owned, 1.0);
+}
 
 fn to_direction(u: f32, v: f32) -> vec3<f32> {
     let s = p.proj_scale;
@@ -98,7 +170,17 @@ fn warp(@builtin(global_invocation_id) gid: vec3<u32>) {
     if (edge <= 0.0) {
         return;
     }
-    let w = clamp(edge / max(p.feather, 1.0), 0.0, 1.0);
+    var w = clamp(edge / max(p.feather, 1.0), 0.0, 1.0);
+    // With seams, the share of the map scales it. The small floor keeps
+    // the feather underneath as the answer wherever no frame that reaches
+    // this pixel owns it — the map is coarser than the output, so at the
+    // frames' outer edges it can name a frame that falls just short.
+    if (p.seam_on != 0u) {
+        let s = seam_share(u, v);
+        if (s.y > 0.0) {
+            w = w * (s.x + 1e-4);
+        }
+    }
     // Into the tile.
     let tx = sx - p.tile_origin.x;
     let ty = sy - p.tile_origin.y;