Replace mzeng with new ezeng method

benches/reductions.rs

@@ -262,6 +262,16 @@ fn reductions_palette_sort_mzeng(b: &mut Bencher) {
     b.iter(|| palette::sorted_palette_mzeng(&png.raw));
 }
 
+#[bench]
+fn reductions_palette_sort_ezeng(b: &mut Bencher) {
+    let input = test::black_box(PathBuf::from(
+        "tests/files/palette_8_should_be_palette_8.png",
+    ));
+    let png = PngData::new(&input, &Options::default()).unwrap();
+
+    b.iter(|| palette::sorted_palette_ezeng(&png.raw));
+}
+
 #[bench]
 fn reductions_palette_sort_battiato(b: &mut Bencher) {
     let input = test::black_box(PathBuf::from(

src/reduction/mod.rs

@@ -161,13 +161,13 @@ pub(crate) fn perform_reductions(
             }
         }
 
-        // Attempt to sort the palette using the mzeng method
+        // Attempt to sort the palette using the ezeng method
         if !deadline.passed() {
-            if let Some(reduced) = sorted_palette_mzeng(input) {
+            if let Some(reduced) = sorted_palette_ezeng(input) {
                 if let ColorType::Indexed { palette } = &reduced.ihdr.color_type {
                     if !palettes.contains(palette) {
                         palettes.push(palette.clone());
-                        eval.try_image_with_description(Arc::new(reduced), "Indexed (mzeng sort)");
+                        eval.try_image_with_description(Arc::new(reduced), "Indexed (ezeng sort)");
                         evaluation_added = true;
                     }
                 }

src/reduction/palette.rs

@@ -151,6 +151,28 @@ pub fn sorted_palette_mzeng(png: &PngImage) -> Option<PngImage> {
     apply_palette_reorder(png, &remapping)
 }
 
+/// Sort the colors in the palette using the ezeng technique, returning the sorted image if successful
+#[must_use]
+pub fn sorted_palette_ezeng(png: &PngImage) -> Option<PngImage> {
+    // Interlacing not currently supported
+    if png.ihdr.bit_depth != BitDepth::Eight || png.ihdr.interlaced {
+        return None;
+    }
+    let palette = match &png.ihdr.color_type {
+        // Images with only two colors will remain unchanged from previous luma sort
+        ColorType::Indexed { palette } if palette.len() > 2 => palette,
+        _ => return None,
+    };
+
+    let matrix = co_occurrence_matrix(palette.len(), png);
+    let edges = weighted_edges(&matrix);
+    let mut remapping = ezeng_reindex(edges, &matrix);
+
+    apply_most_popular_color(png, &mut remapping);
+
+    apply_palette_reorder(png, &remapping)
+}
+
 /// Sort the colors in the palette using the battiato technique, returning the sorted image if successful
 #[must_use]
 pub fn sorted_palette_battiato(png: &PngImage) -> Option<PngImage> {
@@ -364,6 +386,91 @@ fn mzeng_reindex(num_colors: usize, edges: Vec<(usize, usize)>, matrix: &[Vec<u3
     remapping
 }
 
+// Apply a different version of Zeng's technique where the best index is inserted at a position
+// minimizing total cost increase, rather than just the ends. This version is significantly more
+// effective, but is more computationally expensive.
+// The "e" here could mean enhanced, extended, exhaustive, or perhaps elephant if you prefer.
+fn ezeng_reindex(edges: Vec<(usize, usize)>, matrix: &[Vec<u32>]) -> Vec<usize> {
+    // Initialize the mapping list with the two best indices.
+    let mut remapping = vec![edges[0].0, edges[0].1];
+
+    // Initialize the sums with the first two remappings and find the best one
+    let mut sums = Vec::new();
+    let mut best_sum_pos = 0;
+    let mut best_sum = (0, 0);
+    for (i, m_row) in matrix.iter().enumerate() {
+        if i == remapping[0] || i == remapping[1] {
+            continue;
+        }
+        let sum = (i, m_row[remapping[0]] + m_row[remapping[1]]);
+        if sum.1 > best_sum.1 {
+            best_sum_pos = sums.len();
+            best_sum = sum;
+        }
+        sums.push(sum);
+    }
+
+    while !sums.is_empty() {
+        let best_index = best_sum.0;
+        // Try all insertion positions and pick the one minimizing total cost increase.
+        // The cost increase of inserting at position p has two components:
+        // 1. New element cost: Σ w(new, placed[k]) * distance(p, k_after_shift)
+        // 2. Cross-pair cost: Σ w(placed[a], placed[b]) for all pairs straddling p
+        //    (these pairs get pushed 1 unit farther apart by the insertion)
+        let m = remapping.len();
+        let mut best_pos = 0;
+        let mut best_cost = i64::MAX;
+        let mut cross_cost: i64 = 0;
+        for p in 0..=m {
+            // New element's weighted distance to all existing elements
+            let new_cost: i64 = (0..m)
+                .map(|k| {
+                    let dist = if k < p { p - k } else { k + 1 - p };
+                    matrix[best_index][remapping[k]] as i64 * dist as i64
+                })
+                .sum();
+
+            let total = new_cost + cross_cost;
+            if total < best_cost {
+                best_cost = total;
+                best_pos = p;
+            }
+
+            // Update cross_cost for position p+1:
+            // Element at position p moves from "right of split" to "left of split"
+            // cross_cost(p+1) - cross_cost(p) = Σ_{b>p} w(p,b) - Σ_{a<p} w(a,p)
+            if p < m {
+                let rp = remapping[p];
+                for b in (p + 1)..m {
+                    cross_cost += matrix[rp][remapping[b]] as i64;
+                }
+                for a in 0..p {
+                    cross_cost -= matrix[remapping[a]][rp] as i64;
+                }
+            }
+        }
+        remapping.insert(best_pos, best_index);
+
+        // Remove best_sum from sums.
+        sums.swap_remove(best_sum_pos);
+        if !sums.is_empty() {
+            // Update all the sums and find the best one.
+            best_sum_pos = 0;
+            best_sum = (0, 0);
+            for (i, sum) in sums.iter_mut().enumerate() {
+                sum.1 += matrix[best_index][sum.0];
+                if sum.1 > best_sum.1 {
+                    best_sum_pos = i;
+                    best_sum = *sum;
+                }
+            }
+        }
+    }
+
+    // Return the completed remapping
+    remapping
+}
+
 // Calculate an approximate solution of the Traveling Salesman Problem using the algorithm
 // from "An efficient Re-indexing algorithm for color-mapped images" by Battiato et al
 // https://ieeexplore.ieee.org/document/1344033