Acceleration

CHANGELOG.md

@@ -4,6 +4,13 @@ All notable WIMF changes are recorded here. The project follows semantic version
 
 ## 2.2.0 - Unreleased
 
+- First slice of the context-modeled entropy stage (known-flaws A3): lossy
+  wavelet tiles now use marker-free (run, zigzag) coefficient tokens flagged by
+  the tile's reversible byte value 2, removing one mandatory byte per token
+  (about a quarter of packed lossy streams). Pre-2.2 decoders reject the new
+  flag cleanly; the native and Python decoders accept all layouts, and legacy
+  files decode bit-identically. The Python encoder keeps emitting legacy
+  packing, which remains fully valid.
 - Filled the lossy quality dead zone: tile scoring now evaluates a second
   wavelet candidate at 0.9x quantizer scale. Each payload stores its own
   quantizer, so no format change is needed and the decoder is untouched; the

docs/known-flaws.md

@@ -27,8 +27,11 @@ determinism, corruption rejection, and memory guards all audited clean
   per channel). A reversible RGB→YCoCg transform is the standard first win on
   photographic content. Landed: native encoder/decoder plus the Python decoder mirror (the mirror was exactly what run #180s failures exposed). YCoCg-with-offsets refinement pending.
 - **[P0] A3 Generic entropy stage.** Tile payloads are Zstandard bytes of raw
-  prediction residuals or zigzag varint coefficients. No context modeling of
-  residuals/subbands - the structural advantage modern image codecs exploit.
+  prediction residuals or zigzag varint coefficients. Partially improved:
+  lossy wavelet tiles now pack marker-free (run, zigzag) tokens (reversible
+  byte value 2), cutting about a byte per token from packed coefficient
+  streams. Still missing: context modeling of residuals/subbands - the
+  structural advantage modern image codecs exploit.
 - **[P1] A4 Coarse quantizer dead zone.** Largely fixed: monotonic scoring
   (A4), the divisor retune 8.0→16.0, and the 0.9x quantizer sub-step in tile
   scoring removed the non-monotonicity and the one-ladder-notch gap that showed

src/v2_core.cpp

@@ -375,6 +375,44 @@ std::vector<int64_t> unpack_coefficients(const uint8_t* data, size_t size, size_
     return output;
 }
 
+// Marker-free token packing (reversible flag value 2). Legacy packing spends a
+// mandatory 0x00 marker byte on every run/value token, roughly a quarter of
+// the packed lossy stream. V2 tokens are strict (run, zigzag) varint pairs;
+// the reversible byte distinguishes the layouts and pre-v2.2 decoders reject
+// the value cleanly.
+std::vector<uint8_t> pack_coefficients_v2(const std::vector<int64_t>& coefficients) {
+    std::vector<uint8_t> output;
+    size_t run = 0;
+    for (const int64_t value : coefficients) {
+        if (value == 0) {
+            ++run;
+            continue;
+        }
+        append_varint(output, run);
+        run = 0;
+        const uint64_t zigzag = (static_cast<uint64_t>(value) << 1) ^ static_cast<uint64_t>(value >> 63);
+        append_varint(output, zigzag);
+    }
+    // Trailing zero-run is omitted: the decoder zero-fills the remainder.
+    return output;
+}
+
+std::vector<int64_t> unpack_coefficients_v2(const uint8_t* data, size_t size, size_t count) {
+    std::vector<int64_t> output(count);
+    size_t position = 0, index = 0;
+    while (index < count) {
+        if (position >= size) break; // stream ended, remaining coefficients are zero
+        const uint64_t run = read_varint(data, size, position);
+        if (run > count - index) throw std::runtime_error("coefficient zero run exceeds tile");
+        index += static_cast<size_t>(run);
+        if (index == count) break;
+        const uint64_t zigzag = read_varint(data, size, position);
+        output[index++] = static_cast<int64_t>((zigzag >> 1) ^ (0 - (zigzag & 1)));
+    }
+    if (position != size) throw std::runtime_error("trailing coefficient data");
+    return output;
+}
+
 uint32_t next_power_of_two(uint32_t value) {
     uint32_t output = 1;
     while (output < std::max(2u, value)) output <<= 1;
@@ -451,7 +489,9 @@ std::vector<uint8_t> encode_wavelet_tile(const ImageView& tile, uint8_t quality,
     put16(output, static_cast<uint16_t>(padded_height));
     put16(output, static_cast<uint16_t>(padded_width));
     output.push_back(static_cast<uint8_t>(levels | 0x80));
-    output.push_back(lossless ? 1 : 0);
+    // reversible byte doubles as the coefficient-packing selector: 1 lossless
+    // (legacy packing), 2 lossy with marker-free v2 packing, 0 lossy legacy.
+    output.push_back(lossless ? 1 : 2);
     append_float(output, quantizer);
     if (reconstructed) reconstructed->assign(static_cast<size_t>(tile.width) * tile.height * tile.channels * tile.bytes_per_sample, 0);
 
@@ -465,7 +505,8 @@ std::vector<uint8_t> encode_wavelet_tile(const ImageView& tile, uint8_t quality,
         }
         const auto coefficients = wavelet_forward(plane.data(), padded_width, padded_height,
                                                   tile.bytes_per_sample, lossless, levels, quantizer);
-        const auto packed = pack_coefficients(reorder_subbands(coefficients, padded_width, padded_height, levels));
+        const auto ordered = reorder_subbands(coefficients, padded_width, padded_height, levels);
+        const auto packed = lossless ? pack_coefficients(ordered) : pack_coefficients_v2(ordered);
         put32(output, static_cast<uint32_t>(packed.size()));
         output.insert(output.end(), packed.begin(), packed.end());
         if (reconstructed) {
@@ -489,7 +530,7 @@ std::vector<uint8_t> decode_wavelet_tile(const uint8_t* data, size_t size, uint3
     const uint8_t levels = data[4] & 0x7F, reversible = data[5];
     const float quantizer = read_float(data + 6);
     if (padded_width > 256 || padded_height > 256 || padded_width < width || padded_height < height ||
-        levels > 8 || reversible > 1 || !std::isfinite(quantizer) || quantizer <= 0)
+        levels > 8 || reversible > 2 || !std::isfinite(quantizer) || quantizer <= 0)
         throw std::runtime_error("invalid wavelet dimensions");
     size_t position = 10;
     std::vector<uint8_t> output(static_cast<size_t>(width) * height * channels * bytes_per_sample);
@@ -498,12 +539,14 @@ std::vector<uint8_t> decode_wavelet_tile(const uint8_t* data, size_t size, uint3
         const uint32_t packed_size = read32(data + position);
         position += 4;
         if (packed_size > size - position) throw std::runtime_error("truncated wavelet coefficients");
-        auto coefficients = unpack_coefficients(data + position, packed_size,
-                                                static_cast<size_t>(padded_width) * padded_height);
+        auto coefficients = reversible == 2 ? unpack_coefficients_v2(data + position, packed_size,
+                                                                      static_cast<size_t>(padded_width) * padded_height)
+                                            : unpack_coefficients(data + position, packed_size,
+                                                                  static_cast<size_t>(padded_width) * padded_height);
         if (subband) coefficients = restore_raster_order(std::move(coefficients), padded_width, padded_height, levels);
         position += packed_size;
         const auto plane = wavelet_inverse(coefficients.data(), coefficients.size(), padded_width,
-                                           padded_height, bytes_per_sample, reversible != 0, levels, quantizer);
+                                           padded_height, bytes_per_sample, reversible == 1, levels, quantizer);
         for (uint32_t y = 0; y < height; ++y) for (uint32_t x = 0; x < width; ++x) {
             const size_t source = (static_cast<size_t>(y) * padded_width + x) * bytes_per_sample;
             const size_t target = (static_cast<size_t>(y) * width * channels + x * channels + channel) * bytes_per_sample;

wiki/Known-Flaws.md

@@ -27,8 +27,11 @@ determinism, corruption rejection, and memory guards all audited clean
   per channel). A reversible RGB→YCoCg transform is the standard first win on
   photographic content. Landed: native encoder/decoder plus the Python decoder mirror (the mirror was exactly what run #180s failures exposed). YCoCg-with-offsets refinement pending.
 - **[P0] A3 Generic entropy stage.** Tile payloads are Zstandard bytes of raw
-  prediction residuals or zigzag varint coefficients. No context modeling of
-  residuals/subbands - the structural advantage modern image codecs exploit.
+  prediction residuals or zigzag varint coefficients. Partially improved:
+  lossy wavelet tiles now pack marker-free (run, zigzag) tokens (reversible
+  byte value 2), cutting about a byte per token from packed coefficient
+  streams. Still missing: context modeling of residuals/subbands - the
+  structural advantage modern image codecs exploit.
 - **[P1] A4 Coarse quantizer dead zone.** Largely fixed: monotonic scoring
   (A4), the divisor retune 8.0→16.0, and the 0.9x quantizer sub-step in tile
   scoring removed the non-monotonicity and the one-ladder-notch gap that showed

wimf/hybrid.py

@@ -289,6 +289,30 @@ def _varints_decode(data, count):
     return values
 
 
+def _varints_decode_v2(data, count):
+    """Marker-free (run, zigzag) token pairs - mirrors native unpack_coefficients_v2.
+
+    Used for lossy wavelet tiles flagged with reversible == 2; the legacy
+    decoder handles reversible == 0/1 streams unchanged."""
+    values = np.zeros(count, dtype=np.int64)
+    pos = index = 0
+    while index < count:
+        if pos >= len(data):
+            break
+        run, pos = _read_varint(data, pos)
+        if run > count - index:
+            raise ValueError("coefficient zero run exceeds tile")
+        index += run
+        if index == count:
+            break
+        zz, pos = _read_varint(data, pos)
+        values[index] = (zz >> 1) ^ -(zz & 1)
+        index += 1
+    if pos != len(data):
+        raise ValueError("trailing coefficient data")
+    return values
+
+
 def _reorder_subbands(coeff, pw, ph, levels):
     """Raster-order DWT coefficients into dyadic subband sequence (mirrors native)."""
     plane = coeff.reshape(ph, pw)
@@ -322,6 +346,9 @@ def _restore_subbands(flat, pw, ph, levels):
 
 
 def _wavelet_encode(tile, quality, lossless):
+    # The Python encoder intentionally keeps the legacy reversible == 0/1
+    # packing; native encoders emit reversible == 2 (marker-free pairs) for
+    # lossy tiles and both decoders accept every layout.
     h, w, channels = tile.shape
     ph = 1 << int(np.ceil(np.log2(max(2, h))))
     pw = 1 << int(np.ceil(np.log2(max(2, w))))
@@ -354,7 +381,7 @@ def _wavelet_decode(data, h, w, channels, dtype):
         or ph < h
         or pw < w
         or not 0 <= levels <= 8
-        or reversible not in (0, 1)
+        or reversible not in (0, 1, 2)
         or not np.isfinite(q)
         or q <= 0
     ):
@@ -369,18 +396,21 @@ def _wavelet_decode(data, h, w, channels, dtype):
         pos += 4
         if pos + size > len(data):
             raise ValueError("truncated wavelet coefficients")
-        coeff = _varints_decode(data[pos : pos + size], ph * pw)
+        if reversible == 2:
+            coeff = _varints_decode_v2(data[pos : pos + size], ph * pw)
+        else:
+            coeff = _varints_decode(data[pos : pos + size], ph * pw)
         if subband:
             coeff = _restore_subbands(coeff, pw, ph, levels)
         coeff = coeff.reshape(ph, pw)
         pos += size
         if native is not None:
             decoded = native.wavelet_inverse(
-                np.ascontiguousarray(coeff), pw, ph, np.dtype(dtype).itemsize, bool(reversible), levels, q
+                np.ascontiguousarray(coeff), pw, ph, np.dtype(dtype).itemsize, reversible == 1, levels, q
             )
             plane = np.frombuffer(decoded, dtype=dtype).reshape(ph, pw)[:h, :w]
         else:
-            plane = _wavelet_inverse_2d(coeff * q, levels, bool(reversible))[:h, :w]
+            plane = _wavelet_inverse_2d(coeff * q, levels, reversible == 1)[:h, :w]
         planes.append(np.clip(np.rint(plane), 0, max_value).astype(dtype))
     if pos != len(data):
         raise ValueError("trailing wavelet tile data")

wimf/test_v2.py

@@ -22,6 +22,20 @@ def test_all_lossy_quality_and_preset_combinations(quality, preset):
     assert info["width"] == 40 and info["height"] == 24
 
 
+def test_varints_decode_v2_marker_free_pairs():
+    """reversible == 2 wavelet tiles use strict (run, zigzag) varint pairs.
+
+    Vector: [0, 0, 5, -3] encodes as run=2, zigzag(5)=10, run=0, zigzag(-3)=5;
+    the decoder zero-fills the trailing remainder."""
+    stream = bytes([2, 10, 0, 5])
+    assert list(hybrid._varints_decode_v2(stream, 5)) == [0, 0, 5, -3, 0]
+    assert list(hybrid._varints_decode_v2(b"", 3)) == [0, 0, 0]
+    with pytest.raises(ValueError):
+        hybrid._varints_decode_v2(bytes([9, 1]), 3)
+    with pytest.raises(ValueError):
+        hybrid._varints_decode_v2(bytes([1]), 3)
+
+
 def test_v2_lossless_rgb_odd_dimensions(tmp_path):
     arr = np.random.default_rng(1).integers(0, 256, (133, 259, 3), dtype=np.uint8)
     path = tmp_path / "odd.wimf"