-
Notifications
You must be signed in to change notification settings - Fork 320
XUASTC LDR ‐ Weight Grid DCT
Copyright (C) 2025-2026 Binomial LLC. All rights reserved except as granted under the Apache 2.0 LICENSE.
- Intro
- Weight Decoding
- Zigzag Order Generation - All Profiles
- Weight Grid DCT Symbol Decoding - Full ZStd Profile
- Adaptive Quantization - compute_level_scale
- Adaptive Quantization - get_max_span_len
- Baseline JPEG Quantization Table Sampling
- Weight Grid Decoding - decode_block_weights
- 2D IDCT
In addition to simple DPCM coding, XUASTC LDR provides an optional transform-based compression mode for ASTC weight grids. This mode applies a lightweight, JPEG-inspired 2D DCT to each ASTC weight grid (for each interpolation plane), enabling significantly higher compression.
The specific DCT type used by transcoders is the orthonormal inverse DCT-II (or the forward DCT-III), and it is applied in the normalized [0,64] ASTC weight grid domain independent of the block's weight grid BISE level.
Weight grid DCT can be applied to any XUASTC LDR mode configuration: single or dual plane, 1-3 subsets, any endpoint or weight BISE range, any partition pattern, and any supported CEM. A strong encoder can evaluate weight grid DCT error while evaluating trial encodings and partition patterns, and choose the configuration that results in the lowest overall perceptual error.
Each ASTC block explicitly signals whether it uses DPCM or DCT for its weight data. When DCT is selected, the grid is encoded as a small number of symbols rather than a full set of weights: a single DC value, followed by a run-length–encoded sequence of AC coefficients ordered via a weight grid size specific zigzag pattern. In the full Zstd profile, these symbols are stored in independent Zstd streams and decoded per block as needed.
To reconstruct the weights, the decoder performs three main steps:
1. Coefficient decoding
The process is quite similar to baseline JPEG: The weight DCT DC value is decoded first, and handled specially.
Weight AC symbols are read in a precomputed zigzag order specific to the grid size. Runs of zero coefficients (extremely common) are handled using RLE coding. AC coefficient magnitudes are always 8-bits (1-256), and the sign bits are stored uncompressed. AC coefficients larger than 256 (which is extremely rare) can either be clamped during encoding, or the compressor can fall back to DPCM.
2. Dequantization with Adaptive Quantization (AQ)
Unlike JPEG, which uses a fixed quantization matrix, XUASTC LDR weight-grid DCT uses a hybrid of a classic JPEG 8×8 luma quantization table and an adaptive scale factor computed per block. The adaptive scale incorporates two block-dependent factors:
* Weight plane endpoint span length, representing how much color variation that block's weight plane can express, and
* Weight precision, derived from the block’s ASTC weight ISE range.
These factors modulate the effective AC coefficient quantization to compensate for the decoded block's pixel dynamic range defined by the block's endpoints. The quantization table is resampled (via bilinear interpolation) to account for the grid’s actual dimensions, the block's dimensions, and weight grid bilinear upsampling. Each weight plane in dual plane blocks is decompressed separately.
3. Orthonormal inverse DCT (DCT-III)
After deadzone dequantization, the AC and DC coefficients form a complete DCT block of the same dimensions as the ASTC weight grid. The decoder applies a 2D inverse DCT to recover a continuous-valued weight grid. The specification allows this IDCT step to be implemented using either floating-point or integer arithmetic. The reference transcoder uses a separable float IDCT, with automatically generated functions for 2-12 samples, with special fast paths to skip zero coefficients.
4. Weight Re-Quantization
The reconstructed normalized [0,64] weight grid values are then clamped and quantized into the block’s ASTC weight ISE domain. Each weight is mapped through a corresponding ASTC weight value quantization table to yield the exact ISE code used in the final ASTC block. This ensures full compatibility with the ASTC format: once the weight grid has been reconstructed, no further adjustments are needed before constructing the final ASTC block bit pattern.
The DCT mode is entirely local to each block - reconstructed weights do not influence prediction of neighboring blocks and do not feed into any other XUASTC decoding stage. This property makes the DCT transform safe to implement in floating point, while enabling aggressive compression on blocks where DPCM would be suboptimal.
The following pseudo-code describes how XUASTC LDR constructs a DCT coefficient zigzag sequence ordering for an arbitrary width × height ASTC weight grid size. This zigzag pattern is used when encoding and decoding DCT AC coefficients. The reference transcoder precomputes this zigzag order sequence for all valid ASTC weight grid dimensions at initialization. This zigzag pattern is used to map between 1D coefficient indices and 2D DCT positions.
function generate_zigzag_order(width: int, height: int) -> int[] or null
{
assert((width > 0) and (height > 0));
const total: int = width * height;
# Allocate output array of size `total`
pOrder: int[] = allocate_int_array(total);
idx: int = 0;
# Walk over all anti-diagonals with sum index s = x + y
for (s: int = 0; s < (width + height - 1); ++s)
{
# Start x at max(0, s - height + 1), end at min(s, width - 1)
x_start: int = (s < height) ? 0 : (s - height + 1);
x_end: int = (s < width) ? s : (width - 1);
# Diagonal size
diag_size: int = x_end - x_start + 1;
# Temporary array for this diagonal
pDiag: int[] = allocate_int_array(diag_size);
j: int = 0;
for (x: int = x_start; x <= x_end; ++x)
{
y: int = s - x;
assert(j < diag_size);
pDiag[j] = x + y * width;
++j;
}
# Reverse direction on odd diagonals to achieve zigzag ordering
if ((s & 1) == 1)
{
for (k: int = diag_size - 1; k >= 0; --k)
{
assert(idx < total);
pOrder[idx] = pDiag[k];
++idx;
}
}
else
{
for (k: int = 0; k < diag_size; ++k)
{
assert(idx < total);
pOrder[idx] = pDiag[k];
++idx;
}
}
free(pDiag);
}
return pOrder;
}
The following pseudo-code describes how the decoder retrieves the per-block DCT symbols (DC value + run-length encoded AC coefficients) once it determines that weight-grid DCT is active for the current block.
This logic occurs inside the main block-decoding loop of the Full Zstd profile,
after reading mode_byte and verifying the global use_dct flag.
uint32_t get_num_weight_dc_levels(uint32_t weight_ise_range)
{
float scaled_weight_coding_scale = SCALED_WEIGHT_BASE_CODING_SCALE;
if (weight_ise_range <= astc_helpers::BISE_8_LEVELS)
scaled_weight_coding_scale = 1.0f / 8.0f;
return (uint32_t)(64.0f * scaled_weight_coding_scale) + 1;
}
const DCT_RUN_LEN_EOB_SYM_INDEX = 64;
const DCT_MEAN_LEVELS0 = 9, DCT_MEAN_LEVELS1 = 33;
# Determine if the current block actually uses DCT
block_used_dct = false
# global use_dct flag, decoded from the Zstd profile System Header
if use_dct:
block_used_dct = ((mode_byte & XUASTC_LDR_MODE_BYTE_USE_DCT) != 0)
if block_used_dct:
# Determine how many DC-level symbols exist for this block's weight ISE range
num_dc_levels = get_num_weight_dc_levels(log_blk.weight_ise_range)
syms.num_dc_levels = num_dc_levels
# Each ASTC block may have 1 or 2 interpolation planes
for plane_iter in range(total_planes):
# Reset AC coefficient list
syms.coeffs.clear()
# Decode the DC coefficient symbol.
# Some blocks use a 1-byte DC symbol, others a 4-bit DC symbol.
if num_dc_levels == DCT_MEAN_LEVELS1:
syms.dc_sym = mean1_bytes.get_bits8()
else:
syms.dc_sym = mean0_bits.get_bits4()
# AC decoding uses zigzag order, starting at index 1 (0 is DC)
cur_zig_ofs = 1
# Repeatedly read run-lengths and AC coefficient magnitudes
while cur_zig_ofs < total_weights:
run_len = run_bytes.get_bits8()
# If run_len equals the EOB marker, AC decoding stops
if run_len == DCT_RUN_LEN_EOB_SYM_INDEX:
break
# Advance by the run length (skip zeros)
cur_zig_ofs += run_len
# If we run past the coefficient count, this is an error
if cur_zig_ofs >= total_weights:
ERROR("DCT decode error: illegal zigzag offset")
return false
# Read sign and magnitude of the next nonzero AC coefficient
sign_bit = sign_bits.get_bits1()
magnitude = coeff_bytes.get_bits8() + 1
if sign_bit != 0:
magnitude = -magnitude
# Store the (run_len, coefficient) pair
syms.coeffs.push_back( coeff(run_len, magnitude) )
# Advance to the next zigzag position
cur_zig_ofs += 1
# At this point syms contains:
# - syms.dc_sym (DC coefficient symbol)
# - syms.coeffs[] (RLE-coded AC coefficient list)
# Reconstruct the weight grid via inverse quantization and IDCT
if not grid_dct.decode_block_weights(
dct_q, # global DCT q factor, read from the Zstd profile System Header
plane_iter, # weight plane index (0 or 1)
log_blk, # the destination ASTC logical block
dct_work, # IDCT temporary work buffer
syms ): # decoded symbols
ERROR("DCT decode failed")
return false
# If block_used_dct == false, fallback to weight-grid DPCM decoding (handled elsewhere)
The following pseudo-code computes the adaptive per-block quantization scale (“level_scale”) used for DCT coefficient dequantization in the weight-grid DCT path.
# Precomputed table used for adjusting quantization strength based on
# the ASTC weight ISE range (i.e., number of weight levels).
# These factors compensate for the quantization noise introduced at low weight BISE levels,
# which becomes spread across multiple AC coefficients in the frequency domain.
#
# g_scale_quant_steps[BISE_2_LEVELS .. BISE_32_LEVELS]
float g_scale_quant_steps[12] = {
1.51333141f, 1.41198814f, 1.35588217f, 1.31743157f,
1.28835952f, 1.24573100f, 1.21481407f, 1.19067919f,
1.15431654f, 1.12734985f, 1.10601568f, 1.07348967f
};
function compute_level_scale(
q: float, # global DCT quality (1..100)
span_len: float, # endpoint span length of the weight plane
weight_ise_range: int # ASTC weight ISE range (BISE_2_LEVELS..BISE_32_LEVELS)
) -> float
{
# Ensure weight_ise_range is within the legal ASTC domain.
assert(weight_ise_range >= BISE_2_LEVELS and weight_ise_range <= BISE_32_LEVELS)
#--------------------------------------------------------------------------------------
# 1. Standard JPEG-style global quality factor scaling
#--------------------------------------------------------------------------------------
# Clamp q into [1,100]
q = clamp(q, 1.0, 100.0)
float level_scale
if q < 50.0:
# JPEG convention for q<50
level_scale = 5000.0 / q
else:
# JPEG convention for q>=50
level_scale = 200.0 - 2.0 * q
# Convert JPEG’s quant scaling (0..100) to a normalized scale
level_scale = level_scale * (1.0 / 100.0)
#--------------------------------------------------------------------------------------
# 2. Adaptive quantization based on block endpoint span
#--------------------------------------------------------------------------------------
# Span floor prevents overly aggressive scaling for tiny spans.
const float span_floor = 14.0
# span_len is the decoded HDR/LD range the weights operate over.
# Smaller spans → higher adaptive_factor → more aggressive quant.
float adaptive_factor = 64.0 / max(span_len, span_floor)
#--------------------------------------------------------------------------------------
# 3. Adaptive scaling based on number of ASTC weight levels
#--------------------------------------------------------------------------------------
# The fewer ISE levels available, the more quant noise will spread
# into multiple AC coefficients. g_scale_quant_steps compensates.
float weight_quant_adaptive_factor = g_scale_quant_steps[weight_ise_range]
adaptive_factor = adaptive_factor * weight_quant_adaptive_factor
#--------------------------------------------------------------------------------------
# 4. Combine JPEG global scaling with adaptive scaling
#--------------------------------------------------------------------------------------
level_scale = level_scale * adaptive_factor
# Larger level_scale → heavier quantization of AC coefficients.
return level_scale
}
This function computes the endpoint span length used by Adaptive Quantization (AQ)
during weight-grid DCT decoding.
The span measures the magnitude of color variation represented by a block’s
interpolation plane(s). Dual-plane and single-plane blocks use different rules.
function get_max_span_len(log_blk: astc_logical_block, plane_index: uint32_t) -> float
{
span_len: float = 0.0
#----------------------------------------------------------------------
# Case 1: Dual-plane block
# A dual-plane block encodes two independent interpolation planes.
#
# plane_index == 0 → use all color channels EXCEPT the CCS selector
# plane_index == 1 → use ONLY the CCS selector channel
#
# Spatial span = sqrt( Σ (ΔC)^2 ) over the selected channels.
#----------------------------------------------------------------------
if log_blk.dual_plane == true:
{
# Decode low/high endpoints for the first subset (subset 0) by following the ASTC
# specification for endpoint dequantization and decoding to RGBA
color32 l, h
decode_endpoints(
log_blk.color_endpoint_modes[0],
log_blk.endpoints,
log_blk.endpoint_ise_range,
l, h)
# Accumulate squared differences depending on plane_index
for c in 0 .. 3: # c = R,G,B,A channels
if plane_index == 1:
# Only include the selected component
if c == log_blk.color_component_selector:
span_len += square( float(h[c]) - float(l[c]) )
else:
# Include all components EXCEPT the selected one
if c != log_blk.color_component_selector:
span_len += square( float(h[c]) - float(l[c]) )
# Final span length is Euclidean magnitude
span_len = sqrt(span_len)
}
#----------------------------------------------------------------------
# Case 2: Single-plane block
# Compute the span for each partition separately. The block's
# effective span is the maximum span over all partitions.
#----------------------------------------------------------------------
else:
{
for i in 0 .. (log_blk.num_partitions - 1):
# Decode endpoints for partition i
color32 l, h
offset = i * get_num_cem_values(log_blk.color_endpoint_modes[0])
decode_endpoints(
log_blk.color_endpoint_modes[0],
log_blk.endpoints + offset,
log_blk.endpoint_ise_range,
l, h)
# Compute 4-channel Euclidean span for this partition
part_span_len = sqrt(
square( float(h.r) - float(l.r) )
+ square( float(h.g) - float(l.g) )
+ square( float(h.b) - float(l.b) )
+ square( float(h.a) - float(l.a) )
)
# Take maximum over partitions
span_len = max(span_len, part_span_len)
}
return span_len
}
The following pseudo-code describes how XUASTC LDR samples a quantization value from the modified JPEG baseline luma quantization matrix, scaled by block-dependent Adaptive Quantization (AQ). This value is used to dequantize each AC DCT coefficient in the ASTC weight-grid DCT path (but not DC coefficients which are handled specially).
struct sample_quant_table_state
{
float m_q # JPEG-style quality factor (1..100)
float m_sx # horizontal scale for mapping grid freq → JPEG freq
float m_sy # vertical scale
float m_level_scale # adaptive per-block quantization scale (AQ)
function init(q: float,
block_width: uint32_t,
block_height: uint32_t,
level_scale: float)
{
m_q = q
m_level_scale = level_scale
Bx = block_width
By = block_height
# JPEG matrix is 8×8, so scale by 8/Bx and 8/By
m_sx = 8.0 / float(Bx)
m_sy = 8.0 / float(By)
}
}
# The first entry (0,0) has been modified from the original JPEG table to slightly
# lower the DC value when sampling via bilinear interpolation, improving low-frequency
# retention for ASTC grids. Note the DC coefficient is handled specially, and is not
# dequantized using this matrix.
int g_baseline_jpeg_y[8][8] =
{
{ 4, 11, 10, 16, 24, 40, 51, 61 },
{ 12, 12, 14, 19, 26, 58, 60, 55 },
{ 14, 13, 16, 24, 40, 57, 69, 56 },
{ 14, 17, 22, 29, 51, 87, 80, 62 },
{ 18, 22, 37, 56, 68,109,103, 77 },
{ 24, 35, 55, 64, 81,104,113, 92 },
{ 49, 64, 78, 87,103,121,120,101 },
{ 72, 92, 95, 98,112,100,103, 99 }
};
# This function bilinearly samples the JPEG luma quant table at a 2D grid-scaled
# location (rx, ry) and multiplies by the block's AQ level_scale.
# The result determines the quant step for a specific DCT coefficient (x,y).
function sample_quant_table(state: sample_quant_table_state,
x: uint32_t,
y: uint32_t) -> int
{
# Reference rule: (x,y) must not both be zero (DC is handled separately)
assert(x != 0 or y != 0)
# If quality >= 100 → maximal quality → quant scale = 1
if state.m_q >= 100.0:
return 1
# ---------------------------------------------------------------------
# Compute scaled sampling coordinates in JPEG frequency space
# ---------------------------------------------------------------------
ny = float(y)
ry = ny * state.m_sy
nx = float(x)
rx = nx * state.m_sx
assert(x != 0 or y != 0)
# Clamp sampling coordinates to JPEG’s 8×8 domain
i = clamp(rx, 0.0, 7.0)
j = clamp(ry, 0.0, 7.0)
# Integer bounding box for bilinear interpolation
i0 = int(i)
j0 = int(j)
i1 = min(i0 + 1, 7)
j1 = min(j0 + 1, 7)
# Fractional parts for bilinear interpolation
ti = i - float(i0)
tj = j - float(j0)
# ---------------------------------------------------------------------
# Bilinear interpolation of JPEG quant table
# ---------------------------------------------------------------------
a = (1 - ti) * g_baseline_jpeg_y[j0][i0] +
ti * g_baseline_jpeg_y[j0][i1]
b = (1 - ti) * g_baseline_jpeg_y[j1][i0] +
ti * g_baseline_jpeg_y[j1][i1]
base = (1 - tj) * a + tj * b
# ---------------------------------------------------------------------
# Apply adaptive per-block quant scaling, round to int
# ---------------------------------------------------------------------
quant_scale = int(base * state.m_level_scale + 0.5)
# Clamp to minimum of 1
quant_scale = max(1, quant_scale)
return quant_scale
}
This function reconstructs a grid of ASTC block weights, for a single plane, from their compressed representation. It begins by determining the block’s grid dimensions, adaptive quantization parameters, and the appropriate DCT transform tables. The function then decodes and dequantizes the DC and AC DCT coefficients, applying deadzone quantization for most AC frequencies and standard linear dequantization for the lowest-frequency AC terms.
These coefficients are arranged in zigzag order, expanded through run-length decoding, and assembled into a full 2D DCT coefficient grid. A 2D inverse DCT (orthonormal DCT-III) is then applied to convert the frequency-domain information back into spatial-domain weight values, which always range from [0,64] independent of the block's BISE weight range.
Finally, the function adds back the DC mean, clamps the results to the valid dequantized [0,64] ASTC weight range, maps them through the BISE quantization table, and writes the completed ASTC BISE encoded weights into the block’s output structure. This procedure fully reverses the weight-grid compression process, restoring the quantized and BISE encoded texel interpolation weights needed for ASTC texture decoding.
const DEADZONE_ALPHA = .5f;
const SCALED_WEIGHT_BASE_CODING_SCALE = .5f;
# q: global DCT q factor [1,100] - read from system header
# plane_index: [0,1]
# log_blk: Current block's ASTC logical block struct (where decoded weights are placed)
# dct_work: Temporary buffers for IDCT
# syms: Decoded DC/AC symbols from compressed stream
function decode_block_weights(q, plane_index, log_blk, dct_work, syms):
# ------------------------------------------------------------
# 1. Retrieve weight grid dimensions and transform parameters
# ------------------------------------------------------------
grid_w = log_blk.grid_width
grid_h = log_blk.grid_height
total = grid_w * grid_h
num_planes = (log_blk.dual_plane) ? 2 : 1
zigzag = precomputed zigzag order for (grid_w × grid_h)
DCT = precomputed orthonormal inverse DCT-II / forward DCT-III object for this grid size
# ------------------------------------------------------------
# 2. Compute adaptive quantization scaling
# ------------------------------------------------------------
span_len = get_max_span_len(log_blk, plane_index)
level_scale = compute_level_scale(q, span_len, log_blk.weight_ise_range)
scaled_dc_scale =
if (log_blk.weight_ise_range <= astc_helpers::BISE_8_LEVELS) then (1/8)
else SCALED_WEIGHT_BASE_CODING_SCALE
mean_weight = syms.dc_sym / scaled_dc_scale
# ------------------------------------------------------------
# 3. Dequantize all AC coefficients using zigzag traversal
# ------------------------------------------------------------
dct_coeffs = array(total) filled with 0
zig = 1 # DC is index 0
for each (run_len, value) in syms.coeffs:
if ((zig + run_len) > total)
return false; # decode error
zig += run_len
if zig >= total:
break
idx = zigzag[zig]
x = idx % grid_w
y = idx / grid_w
quant_step = sample_quant_table(q, x, y, level_scale)
# normal or deadzone dequantization of AC coefficient
if (x,y) == (1,0) or (x,y) == (0,1):
dct_coeffs[idx] = value * quant_step
else if value == 0 or quant_step <= 0:
dct_coeffs[idx] = 0
else:
tau = DEADZONE_ALPHA * quant_step
mag = tau + abs(value) * quant_step
dct_coeffs[idx] = (value < 0) ? -mag : mag
zig += 1
# ------------------------------------------------------------
# 4. Apply 2D IDCT (orthonormal DCT-III)
# ------------------------------------------------------------
idct_coeffs = DCT.inverse(dct_coeffs, dct_work)
# ------------------------------------------------------------
# 5. Reconstruct final grid weights using BISE mapping
# ------------------------------------------------------------
# quant_tab[] converts [0,64] values to the nearest quantized ASTC BISE index
quant_tab = astc_helpers::lookup_quant_table(log_blk.weight_ise_range)
for y in 0 .. grid_h-1:
for x in 0 .. grid_w-1:
# Add mean_weight (this could also be done via placing it into the DC coefficient before IDCT), round to int and clamp to [0,64]
w = mean_weight + idct_coeffs[x,y]
w = round(w)
w = clamp(w, 0, 64)
# Write plane weight
out_index = (x + y * grid_w) * num_planes + plane_index
log_blk.weights[out_index] = quant_tab[w]
return true
This function implements an orthonormal 2D inverse DCT (DCT-III) used to reconstruct ASTC decoder weight grids from their frequency-domain representation. It operates as a separable transform: a vertical 1D IDCT is applied first, followed by a horizontal 1D IDCT, using cosine and α-scaling tables generated during initialization.
Because 2D ASTC weight grids may have several different dimensions (2-12 on each dimension, potentially non-square), the cosine tables and scaling factors for every supported grid size should be precomputed once before any blocks are decoded, avoiding repeated trigonometric evaluation during runtime. The implementation shown here is the straightforward O(N²) formulation of the IDCT, chosen for clarity and simplicity. Faster, factorized IDCT algorithms - such as Loeffler, AAN, Chen, or other well-known DCT decompositions - could be substituted in performance-critical builds, but this version provides a simple and reliable baseline transform with predictable numerical behavior. All-integer implementations are also possible.
class dct2f:
constant MAX_SIZE = 12 # max grid dimension for ASTC is 12
members:
rows # number of DCT rows (height)
cols # number of DCT columns (width)
c_col[] # cosine table for vertical dimension: c_col[u * rows + x]
c_row[] # cosine table for horizontal dimension: c_row[v * cols + y]
a_col[] # alpha(u) scaling for rows
a_row[] # alpha(v) scaling for cols
constructor:
rows = 0
cols = 0
-------------------------------------------------------------------------
function init(rows, cols)
-------------------------------------------------------------------------
if rows < 2 or rows > MAX_SIZE or cols < 2 or cols > MAX_SIZE:
assert false
return false
self.rows = rows
self.cols = cols
allocate c_col of size rows * rows, zero-filled
allocate c_row of size cols * cols, zero-filled
allocate a_col of size rows, zero-filled
allocate a_row of size cols, zero-filled
pi = 3.14159265358979323846
# alpha scaling for orthonormal DCT
inv_m = 1.0 / rows
a_col[0] = sqrt(inv_m)
for u in 1 .. rows-1:
a_col[u] = sqrt(2.0 * inv_m)
inv_n = 1.0 / cols
a_row[0] = sqrt(inv_n)
for v in 1 .. cols-1:
a_row[v] = sqrt(2.0 * inv_n)
# precompute cosine terms for vertical dimension
for u in 0 .. rows-1:
for x in 0 .. rows-1:
angle = pi * (2*x + 1) * u / (2 * rows)
c_col[u * rows + x] = cos(angle)
# precompute cosine terms for horizontal dimension
for v in 0 .. cols-1:
for y in 0 .. cols-1:
angle = pi * (2*y + 1) * v / (2 * cols)
c_row[v * cols + y] = cos(angle)
return true
-------------------------------------------------------------------------
# convenience wrapper: stride is cols
-------------------------------------------------------------------------
function inverse(src, dst, work)
call inverse(src, cols, dst, cols, work)
-------------------------------------------------------------------------
# Full 2D inverse transform (orthonormal DCT-III), separable:
# First pass – vertical (rows)
# Second pass – horizontal (cols)
#
# src_stride and dst_stride must equal cols.
-------------------------------------------------------------------------
function inverse(src, src_stride, dst, dst_stride, work)
assert rows > 0 and cols > 0
resize work to (rows * cols)
m = rows
n = cols
work_matrix = work[]
---------------------------------------------------------------------
# Pass 1: Vertical IDCT
#
# For each column v:
# sums[x] = sum over u: ( src[u,v] * alpha(u) * cos term )
#
# Result is stored into work_matrix[x,v].
---------------------------------------------------------------------
for v in 0 .. n-1:
declare sums[0 .. MAX_SIZE-1] = all zeros
for u in 0 .. m-1:
# Skip coefficient if zero
yU = src[u*src_stride + v]
if yU == 0:
continue
yU = yU * a_col[u]
for x in 0 .. m-1:
cU = c_col[u * m + x]
sums[x] += yU * cU
for x in 0 .. m-1:
work_matrix[x * n + v] = sums[x]
---------------------------------------------------------------------
# Pass 2: Horizontal IDCT
#
# For each row x:
# For each output y:
# sum over v: ( work[x,v] * alpha(v) * cos term )
#
# Result is written into dst[x,y].
---------------------------------------------------------------------
for x in 0 .. m-1:
row_temp = pointer to work_matrix[x * n]
row_out = pointer to dst[x * dst_stride]
for y in 0 .. n-1:
s = 0.0
for v in 0 .. n-1:
cV = c_row[v * n + y]
s += (row_temp[v] * a_row[v]) * cV
row_out[y] = s