Lifting-Based Wavelet & Color Reconstruction

Overview

After a JPEG2000 decoder has entropy-decoded the wavelet coefficients from the codestream, it must reconstruct the original pixel values through two stages: inverse discrete wavelet transform (IDWT) and inverse color transform. The JPEG2000 standard specifies two wavelet filters—the reversible 5/3 filter for lossless compression and the irreversible 9/7 filter for lossy compression—both implemented via the lifting scheme rather than convolution.

The lifting scheme decomposes the wavelet transform into a sequence of simple predict-and-update steps that operate in-place, halving the memory requirements compared to convolution-based implementations. Each lifting step modifies either the even-indexed or odd-indexed samples using values from the other set, building up the full filter response incrementally. This structure also makes the inverse trivially derivable: reverse the step order and negate the coefficients.

Color reconstruction follows the IDWT. JPEG2000 defines two color transforms: the Irreversible Color Transform (ICT), which converts YCbCr to RGB using floating-point arithmetic, and the Reversible Color Transform (RCT), which uses integer arithmetic to achieve exact reconstruction. The codec library’s decoder implements both transforms in Zig, handling edge cases like symmetric extension at tile boundaries and proper bit-depth scaling.

Core Concepts

The Lifting Scheme

Traditional wavelet transforms apply a filter bank via convolution, which requires temporary buffers for the output. The lifting scheme, introduced by Sweldens (1996), factors any wavelet filter into a sequence of alternating predict and update steps. The key insight is that each step only needs access to neighboring samples, enabling in-place computation.

For the inverse transform (reconstruction), the lifting steps are applied in reverse order with negated coefficients. If the forward transform applies steps A, B, C, the inverse applies -C, -B, -A. This symmetry means the encoder and decoder share the same code with only the sign and order changed.

Reversible 5/3 Filter

The Le Gall 5/3 filter is used for lossless compression. It uses integer arithmetic exclusively, ensuring perfect reconstruction. The filter has two lifting steps for the inverse transform:

1. Update: odd[n] -= floor((even[n] + even[n+1]) / 2) — remove the prediction from the detail coefficients
2. Predict: even[n] += floor((odd[n-1] + odd[n] + 2) / 4) — reconstruct the approximation coefficients

The floor operations and integer divisors guarantee that the inverse exactly undoes the forward transform, enabling bit-exact lossless reconstruction regardless of platform.

Irreversible 9/7 Filter

The Cohen-Daubechies-Feauveau 9/7 filter is used for lossy compression and provides better energy compaction than the 5/3 filter, yielding higher compression ratios at the same visual quality. It has four lifting steps with irrational coefficients:

1. odd[n] -= 1.586134342 * (even[n] + even[n+1])
2. even[n] -= 0.052980118 * (odd[n-1] + odd[n])
3. odd[n] += 0.882911076 * (even[n] + even[n+1])
4. even[n] += 0.443506852 * (odd[n-1] + odd[n])

After the lifting steps, the even and odd samples are scaled by 1/K and K respectively, where K = 1.230174105.

Symmetric Boundary Extension

At tile and image boundaries, the lifting steps need samples that lie outside the signal. JPEG2000 specifies symmetric (mirror) extension: the signal is reflected at the boundary. For a signal [a, b, c, d], the left extension is [..., c, b, a, b, c, d, ...] and the right extension is [..., a, b, c, d, c, b, ...]. The decoder must handle both even-length and odd-length signals, which have subtly different extension rules.

Color Transforms

JPEG2000 compresses images in a decorrelated color space. The encoder applies a forward color transform before the DWT; the decoder applies the inverse after the IDWT.

Irreversible Color Transform (ICT): Standard YCbCr-to-RGB conversion:

• R = Y + 1.402 * Cr
• G = Y - 0.34413 * Cb - 0.71414 * Cr
• B = Y + 1.772 * Cb

Reversible Color Transform (RCT): Integer-exact conversion:

• G = Y0 - floor((Y2 + Y1) / 4)
• R = Y2 + G
• B = Y1 + G

How It Works

The reconstruction pipeline processes each tile independently. Within a tile, the decoder first dequantizes the wavelet coefficients (multiplying by the quantization step size), then applies the inverse DWT level by level, and finally applies the inverse color transform.

1. Dequantize: Multiply each subband’s coefficients by its quantization step size.
2. IDWT Level N: Combine LL_N, HL_N, LH_N, HH_N to produce LL_{N-1} using the lifting-based inverse filter.
3. Repeat: Apply IDWT for levels N-1 down to 0.
4. Inverse color transform: Convert from YCbCr (or YDbDr) back to RGB.
5. Clamp and bit-shift: Ensure pixel values are within the valid range for the bit depth.

// Inverse 5/3 lifting transform (1D, in-place)
fn inverseLift53(samples: []i32, len: usize) void {
    const half = len / 2;
    // samples[0..half] = even (low-pass), samples[half..len] = odd (high-pass)

    // Step 1: Update odd samples (undo prediction)
    var i: usize = 0;
    while (i < half - 1) : (i += 1) {
        samples[half + i] -= @divFloor(samples[i] + samples[i + 1], 2);
    }
    // Boundary: last odd sample
    samples[half + half - 1] -= samples[half - 1]; // symmetric extension

    // Step 2: Predict even samples (undo update)
    samples[0] += @divFloor(samples[half] + samples[half] + 2, 4); // symmetric ext
    i = 1;
    while (i < half) : (i += 1) {
        samples[i] += @divFloor(samples[half + i - 1] + samples[half + i] + 2, 4);
    }

    // Interleave: rearrange [even, odd] back to [e0, o0, e1, o1, ...]
    deinterleave(samples, len);
}

// Inverse 9/7 lifting transform (1D, in-place)
fn inverseLift97(samples: []f64, len: usize) void {
    const half = len / 2;

    // CDF 9/7 lifting coefficients
    const alpha: f64 = -1.586134342;
    const beta: f64 = -0.052980118;
    const gamma: f64 = 0.882911076;
    const delta: f64 = 0.443506852;
    const K: f64 = 1.230174105;

    // Undo scaling
    for (0..half) |i| samples[i] /= K;        // even
    for (half..len) |i| samples[i] *= K;       // odd

    // Reverse lifting steps (negate coefficients, reverse order)
    // Step 4 (inverse): even -= delta * (odd[n-1] + odd[n])
    liftStep(samples, half, -delta, .even_from_odd);
    // Step 3 (inverse): odd -= gamma * (even[n] + even[n+1])
    liftStep(samples, half, -gamma, .odd_from_even);
    // Step 2 (inverse): even -= beta * (odd[n-1] + odd[n])
    liftStep(samples, half, -beta, .even_from_odd);
    // Step 1 (inverse): odd -= alpha * (even[n] + even[n+1])
    liftStep(samples, half, -alpha, .odd_from_even);

    deinterleave(samples, len);
}

Practical Applications

Progressive Image Rendering

Because the IDWT reconstructs images level by level, a decoder can render a low-resolution preview by stopping after the first few levels. A 4-level decomposition of a 4096×4096 image produces a 256×256 LL approximation at level 4. The decoder can display this thumbnail immediately while continuing to decode higher-resolution levels in the background.

Region-of-Interest Decoding

The tile-based architecture combined with the lifting scheme enables efficient ROI decoding. To reconstruct a specific region, the decoder only needs to process the tiles that intersect the region. Within each tile, only the subbands that contribute to the requested region need to be inverse-transformed. The lifting scheme’s local data access pattern means that reconstructing a small region touches only nearby coefficients.

HDR and Scientific Imaging

JPEG2000 supports bit depths up to 38 bits per component, making it suitable for HDR photography and scientific imaging. The reversible 5/3 filter preserves every bit of precision, and the RCT color transform avoids floating-point rounding. This makes JPEG2000 the format of choice for astronomy (FITS wrapped in JP2), remote sensing, and digital cinema (DCI specification mandates JPEG2000).

Implementation Patterns

Separable 2D Transform

The 2D inverse DWT is implemented as a separable filter: apply the 1D inverse transform to all rows, then apply it to all columns (or vice versa). This reduces a 2D operation with O(N²) filter taps to two passes of O(N) taps each. The order (rows-first or columns-first) does not affect the result due to the separability property.

fn inverseDWT2D(
    tile: *TileBuffer,
    level: u32,
    filter: WaveletFilter,
) void {
    const width = tile.widthAtLevel(level);
    const height = tile.heightAtLevel(level);

    // Step 1: Inverse transform all columns (vertical)
    var col: u32 = 0;
    while (col < width) : (col += 1) {
        var column = tile.getColumn(col, level);
        switch (filter) {
            .reversible_5_3 => inverseLift53(column, height),
            .irreversible_9_7 => inverseLift97(column, height),
        }
    }

    // Step 2: Inverse transform all rows (horizontal)
    var row: u32 = 0;
    while (row < height) : (row += 1) {
        var rowdata = tile.getRow(row, level);
        switch (filter) {
            .reversible_5_3 => inverseLift53(rowdata, width),
            .irreversible_9_7 => inverseLift97(rowdata, width),
        }
    }
}

SIMD-Friendly Memory Layout

For the irreversible 9/7 filter, the lifting steps involve multiply-add operations on adjacent samples—a pattern that maps directly to SIMD instructions. Storing coefficients in contiguous arrays aligned to 16 or 32 bytes allows the compiler to auto-vectorize the inner loops. The even/odd interleaving used during reconstruction is naturally amenable to SIMD shuffle operations.

Bit-Depth Clamping

After the inverse color transform, pixel values may exceed the valid range due to floating-point quantization artifacts. The decoder must clamp each component to [0, 2^bitdepth - 1] before writing to the output buffer. For signed components (rare but valid), the range is [-2^(bitdepth-1), 2^(bitdepth-1) - 1].

fn clampAndOutput(
    samples: []f64,
    output: []u8,
    bit_depth: u5,
    is_signed: bool,
) void {
    const max_val = (@as(i64, 1) << bit_depth) - 1;
    const min_val: i64 = if (is_signed) -(@as(i64, 1) << (bit_depth - 1)) else 0;
    const offset: i64 = if (is_signed) @as(i64, 1) << (bit_depth - 1) else 0;

    for (samples, 0..) |s, i| {
        var val = @floatToInt(i64, @round(s));
        val = @max(min_val, @min(val, max_val));
        val += offset; // Shift signed to unsigned for output
        output[i] = @intCast(u8, val);
    }
}

Common Pitfalls

• Mixing up forward and inverse lifting order: The inverse transform reverses the step order and negates the coefficients. Applying forward-order steps during decoding produces plausible-looking but subtly wrong output that passes casual visual inspection but fails bit-exact comparison tests.
• Incorrect symmetric extension for even-length signals: Even-length and odd-length signals use different extension rules at the right boundary. An even-length signal [a,b,c,d] extends as [...,c,d,c,b,...] while an odd-length signal [a,b,c] extends as [...,b,c,b,...]. Getting this wrong introduces boundary artifacts visible as vertical or horizontal lines at tile edges.
• Applying ICT to lossless-decoded data: If the codestream uses the reversible 5/3 filter, the color transform must be RCT, not ICT. The COD marker’s transform flag indicates which was used. Applying the wrong color transform produces severely corrupted colors.
• Integer overflow in 5/3 lifting for high bit depths: The intermediate values during 5/3 lifting can exceed the input bit depth by several bits. For 16-bit input, intermediate values can reach 18 bits. Using i16 instead of i32 causes silent overflow and corrupted coefficients.
• Forgetting the 9/7 scaling step: The 9/7 filter requires a scaling step (multiply by 1/K and K) before the lifting steps. Omitting this produces output with incorrect magnitude—the image appears too bright or too dark depending on which scale factor was missed.