Add pairwise swaps

src/reduction/palette.rs

@@ -168,6 +168,10 @@ pub fn sorted_palette_ezeng(png: &PngImage) -> Option<PngImage> {
     let edges = weighted_edges(&matrix);
     let mut remapping = ezeng_reindex(edges, &matrix);
 
+    // Perform additional optimization with pairwise swaps. 50 is a good value for max_dist to
+    // keep performance reasonable and can actually be better than the full 255 in some cases.
+    pairwise_swap_search(&mut remapping, &matrix, 50);
+
     apply_most_popular_color(png, &mut remapping);
 
     apply_palette_reorder(png, &remapping)
@@ -551,3 +555,50 @@ fn battiato_reindex(num_colors: usize, edges: Vec<(usize, usize)>) -> Vec<usize>
     // Return the completed chain
     chains.swap_remove(0)
 }
+
+// Pairwise swap: for each pair (a, b), swap if it reduces cost.
+// This is an effective means of refining the result of another algorithm. It's currently rather
+// brutish and not very efficient - there are probably ways it could be optimized, or perhaps the
+// ezeng algorithm could be modified to achieve a similar effect without needing this step at all.
+//
+// `max_dist` limits the distance between pairs to consider to keep it performant. A value of 1 is
+// quite fast and still provides a good improvement, but higher values can be a little better.
+fn pairwise_swap_search(remapping: &mut [usize], matrix: &[Vec<u32>], max_dist: usize) {
+    let num_colors = remapping.len();
+
+    // Build a flat weight matrix for cache-friendly inner-loop access
+    let mut weights = vec![0_i64; num_colors * num_colors];
+    for i in 0..num_colors {
+        for j in 0..num_colors {
+            weights[i * num_colors + j] = matrix[i][j] as i64;
+        }
+    }
+
+    // Repeat until no further improvement, but cap the number of passes to num_colors for safety
+    let mut pass = 0;
+    let mut improved = true;
+    while improved && pass < num_colors {
+        pass += 1;
+        improved = false;
+        for a in 0..num_colors - 1 {
+            for b in (a + 1)..(a + 1 + max_dist).min(num_colors) {
+                let va = remapping[a];
+                let vb = remapping[b];
+                let mut delta: i64 = 0;
+                for (i, &vi) in remapping.iter().enumerate() {
+                    if i == a || i == b {
+                        continue;
+                    }
+                    let weight_diff = weights[va * num_colors + vi] - weights[vb * num_colors + vi];
+                    let dist_diff = (b as i64 - i as i64).unsigned_abs() as i64
+                        - (a as i64 - i as i64).unsigned_abs() as i64;
+                    delta += weight_diff * dist_diff;
+                }
+                if delta < 0 {
+                    remapping.swap(a, b);
+                    improved = true;
+                }
+            }
+        }
+    }
+}