Seam Carving – 讓圖片不成比例地縮小

Seam Carving 是我在 Computational Thinking 時學到的演算法，可以針對圖片的內容做縮放。一般將圖片渲染到網頁上時，如果大小不符合容器高、寬，我們會採取一些做法：

1. 將圖片等比例縮放
2. 剪裁掉一部分的圖片內容

Seam Carving 這個演算法可以在改變圖片比例的情況下做到裁剪，其原理是每次縮放時都找出圖片當中最不重要的部分並刪除，然後將圖片縫合起來。今天我們就來談談這個演算法吧！

我們來看看這張圖片

如果我們將他的比例改成這樣

可以發現圖片的內容整個都變形，非常不好看。

Seam Carving 這個演算法可以做到大致保持圖片內容，又可以不按照比例縮放。實際上是怎麼做到的呢？今天就來介紹一下這個演算法。

Seam Carving

Seam Carving 的核心原理是找出圖片中最不重要部分刪除後，把剩下的部分「縫合」起來。那麼重點來了，要怎麼找到圖片中「不重要」的部分？又是怎麼定義不重要的呢？

邊緣檢測

我們先來看看這張圖片

再來看看這張圖片

這張圖片最重要的地方是黑白的交界點，只要能夠維持這個邊界，就算稍微修改比例，也不會感受到太大的差異。用肉眼看是看得懂，那麼在電腦視覺當中要怎麼找到這些點呢？

Sobel 邊緣檢測

Sobel 包含兩個 3x3 的矩陣：

$$K_x = \begin{bmatrix} 1 & 0 & -1 \\ 2 & 0 & -2 \\ 1 & 0 & -1 \end{bmatrix} \qquad K_y = \begin{bmatrix} 1 & 2 & 1 \\ 0 & 0 & 0 \\ -1 & -2 & -1 \end{bmatrix}$$

如果將圖片當作一個矩陣 $A$，把每個像素和周圍 3x3 的鄰居分別與 $K_x$、$K_y$ 相乘後加總，等於對影像做卷積（Convolution），可以得到水平和垂直方向的梯度 $G_x$、$G_y$。Sobel 透過計算離散差分，算出圖像亮度的梯度值：

$$G_x = K_x * A \qquad G_y = K_y * A \qquad G = \sqrt{G_x^2 + G_y^2}$$

從矩陣來看可能不好理解，用程式碼表示會清楚很多。舊版的文章用 JavaScript 跑兩層迴圈，一個像素一個像素算；但每個像素的計算互不相關，非常適合交給 GPU 平行處理。下面用 WebGPU 的 compute shader（WGSL）來寫，每個 thread 負責一個像素：

const KERNEL_X = array<f32, 9>(1.0, 0.0, -1.0, 2.0, 0.0, -2.0, 1.0, 0.0, -1.0);
const KERNEL_Y = array<f32, 9>(1.0, 2.0, 1.0, 0.0, 0.0, 0.0, -1.0, -2.0, -1.0);

fn luminance(x: i32, y: i32) -> f32 {
  let cx = u32(clamp(x, 0, i32(p.width) - 1));
  let cy = u32(clamp(y, 0, i32(p.height) - 1));
  let rgb = unpack4x8unorm(pixels[cy * p.stride + cx]).rgb;
  return dot(rgb, vec3f(0.299, 0.587, 0.114));
}

@compute @workgroup_size(8, 8)
fn computeEnergy(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= p.width || id.y >= p.height) { return; }

  var gx = 0.0;
  var gy = 0.0;
  for (var i = 0; i < 9; i++) {
    let l = luminance(i32(id.x) + i % 3 - 1, i32(id.y) + i / 3 - 1);
    gx += KERNEL_X[i] * l;
    gy += KERNEL_Y[i] * l;
  }
  energy[id.y * p.stride + id.x] = sqrt(gx * gx + gy * gy);
}

luminance 先把 RGB 轉成灰階亮度，再把 3x3 的鄰居分別乘上兩個 kernel，最後算出梯度大小。邊界的像素用 clamp 取最近的像素代替。

這裡的 stride 是原圖的寬度。之後每移除一條 seam 圖片就會少一個像素寬，但 buffer 的大小不變，所以用固定的 stride 計算索引，再用目前的 width 判斷範圍。

接下來看看經過 Sobel 處理的圖片，這是原圖：

經過 Sobel 邊緣檢測後：

經過處理後的圖，可以發現交界處的部分（亮度變化大）的地方，值會越接近白色（也就是越靠近 255）。而亮度變化不明顯的地方則是呈黑色（靠近 0）

計算圖片中的不重要成分

由於邊緣是圖片中比較重要的部分（越白的地方越重要），因此我們可以由上往下從相鄰像素當中找出最小值，然後把它從圖片移除就可以了。

然而，如果我們直接從上往下去找（Greedy），他可能不是全域最小值，例如：

$$\begin{bmatrix} 0.9 & 0.8 & 0.7 \\ 0.1 & 0.1 & 0.05 \\ 0.1 & 1.5 & 1.5 \end{bmatrix}$$

如果從最上排最小的 0.7 開始，每次都往下選最小的鄰居，會得到 0.7+0.05+1.5=2.25。但如果第二步選 0.1，則會得到 0.7+0.1+0.1=0.9。每一步都選最小值，不代表整條路徑最小。我們把經由 Sobel 處理過的矩陣資料叫做 energy，可以先針對此 energy 作預處理。

先從最底部開始，首先將原本的 energy 最底部複製過去，由下往上，每一格都加上它下方三個相鄰元素中的最小值。這樣每一格存的就是「從這格走到底部的最小總 energy」，這就是動態規劃的寫法：

$$\begin{bmatrix} 1.1 & 1.0 & 0.9 \\ 0.2 & 0.2 & 1.55 \\ 0.1 & 1.5 & 1.5 \end{bmatrix}$$

這樣子就可以直接從最上排最小的 0.9 開始找，往下依序是 0.2、0.1，由於已經經過預處理，可以保證這條路線是總 energy 最少的。

接下來就是將對應座標 [x, y] 記錄起來，以上述的例子來說可以得到：

[[2, 0], [1, 1], [0, 2]]

在 GPU 上算動態規劃

Sobel 每個像素都可以獨立計算，但動態規劃不行：每一列都要等下面那一列算完。所以這邊的做法是只開一個 workgroup，讓 256 個 thread 分攤同一列的每一格，算完一列用 storageBarrier() 等所有 thread 寫入完成，再往上算下一列。

找出最小總 energy 的路徑時，每一步都依賴上一步的結果，只能循序執行，所以交給一個 thread 處理：

const THREADS = 256u;

@compute @workgroup_size(THREADS)
fn findSeam(@builtin(local_invocation_index) tid: u32) {
  let last = (p.height - 1u) * p.stride;
  for (var x = tid; x < p.width; x += THREADS) {
    cost[last + x] = energy[last + x];
  }

  for (var y = i32(p.height) - 2; y >= 0; y--) {
    storageBarrier();
    let row = u32(y) * p.stride;
    let below = row + p.stride;
    for (var x = tid; x < p.width; x += THREADS) {
      let left = max(x, 1u) - 1u;
      let right = min(x + 1u, p.width - 1u);
      let minBelow = min(cost[below + left], min(cost[below + x], cost[below + right]));
      cost[row + x] = energy[row + x] + minBelow;
    }
  }
  storageBarrier();

  if (tid != 0u) { return; }

  var x = 0u;
  for (var i = 1u; i < p.width; i++) {
    if (cost[i] < cost[x]) { x = i; }
  }
  seam[0] = x;

  for (var y = 1u; y < p.height; y++) {
    let row = y * p.stride;
    let prev = x;
    if (prev > 0u && cost[row + prev - 1u] < cost[row + x]) { x = prev - 1u; }
    if (prev + 1u < p.width && cost[row + prev + 1u] < cost[row + x]) { x = prev + 1u; }
    seam[y] = x;
  }
}

最後把 seam 右邊的像素往左移一格就完成了。每個像素一樣可以平行處理，寫到另一個 buffer 之後再複製回來，避免讀到其他 thread 已經改過的值：

@compute @workgroup_size(8, 8)
fn removeSeam(@builtin(global_invocation_id) id: vec3u) {
  if (id.x >= p.width - 1u || id.y >= p.height) { return; }

  let source = select(id.x, id.x + 1u, id.x >= seam[id.y]);
  scratch[id.y * p.stride + id.x] = pixels[id.y * p.stride + source];
}

JavaScript 這邊只要依序 dispatch 這三個 pass。每移除一條 seam，寬度就減 1：

const encoder = device.createCommandEncoder()
const pass = encoder.beginComputePass()
pass.setBindGroup(0, bindGroup)
pass.setPipeline(energyPipeline)
pass.dispatchWorkgroups(Math.ceil(width / 8), Math.ceil(height / 8))
pass.setPipeline(seamPipeline)
pass.dispatchWorkgroups(1)
pass.setPipeline(removePipeline)
pass.dispatchWorkgroups(Math.ceil(width / 8), Math.ceil(height / 8))
pass.end()
encoder.copyBufferToBuffer(scratch, 0, pixels, 0, pixels.size)
device.queue.submit([encoder.finish()])

畫面的部分則是用 render pass 直接讀取同一個 storage buffer 畫到 canvas 上，像素從頭到尾都留在 GPU，不需要讀回 CPU。

效果

拉動滑桿就可以縮小圖片的寬度，紅線是下一條要被移除的 seam。切換到 Energy 可以看到 Sobel 處理後的結果，也可以換成自己的圖片（圖片只會在瀏覽器裡處理，不會上傳）。

如果瀏覽器不支援 WebGPU，可以看這部影片：

雖然我們說可以讓圖片不成比例縮放，但根據圖片內容，還是會看起來怪怪的，尤其是當圖片的重要訊息被移除後會越來越奇怪。儘管如此我還是覺得它是個很有趣的演算法。

移除不重要的部分，這個中心思想其實在很多地方都看得到，例如聲音採樣，人耳聽不到的高頻可以移除掉節省空間；矩陣太複雜了，透過 SVD 奇異值分解用重要的特徵向量來表示。雖然看似互不關聯，但觀念相通。

後記

舊版的實作直接用 JavaScript 迴圈處理每個像素，圖片尺寸大一點就有明顯的延遲。改用 WebGPU 之後，Sobel 和移除 seam 都是每個像素平行計算，可以即時拖動滑桿。

比較可惜的是動態規劃那一段：每一列都依賴下一列，只能用一個 workgroup 逐列計算，沒辦法吃滿整張 GPU。另一個可以優化的地方是移除一條 seam 之後，只有 seam 附近的 energy 會改變，其實不需要每次都重算整張圖。

在實際使用上，雖然效果不錯，但其實還是要看圖片的內容為何，如果圖片的細節都很重要，都有非常明顯的亮度變化，那麼就算用 Seam Carving 也還是會失真很多。

相關資源

本篇很多概念是從 Computational Thinking – Seam Carving 這部影片而來，解說是 3Blue1Brown 的作者 Grant Sanderson！裡頭有大量的視覺化解說，非常好理解。