j2k_codec_math/
dwt.rs

1//! JPEG 2000 discrete wavelet transform constants.
2
3mod linearized53;
4pub use linearized53::{
5    linearized_dwt53_row, Dwt53Band, Dwt53LinearRow, Dwt53LinearTap, DWT53_MAX_HIGH_LINEAR_TAPS,
6    DWT53_MAX_LINEAR_TAPS,
7};
8
9/// Return the maximum number of DWT decomposition levels supported by an
10/// image geometry.
11///
12/// The shared encoder policy is `floor(log2(min(width, height)))`; a zero or
13/// unit-length axis supports no decomposition levels.
14#[must_use]
15pub const fn max_decomposition_levels(width: u32, height: u32) -> u8 {
16    j2k_types::encode_geometry::maximum_decomposition_levels(width, height)
17}
18
19/// Forward irreversible 9/7 lifting step alpha, rounded for existing CPU/GPU paths.
20pub const DWT97_ALPHA_F32: f32 = -1.586_134_3;
21/// Forward irreversible 9/7 lifting step beta, rounded for existing CPU/GPU paths.
22pub const DWT97_BETA_F32: f32 = -0.052_980_117;
23/// Forward irreversible 9/7 lifting step gamma, rounded for existing CPU/GPU paths.
24pub const DWT97_GAMMA_F32: f32 = 0.882_911_1;
25/// Forward irreversible 9/7 lifting step delta, rounded for existing CPU/GPU paths.
26pub const DWT97_DELTA_F32: f32 = 0.443_506_87;
27/// Irreversible 9/7 scaling factor, rounded for existing CPU/GPU paths.
28pub const DWT97_KAPPA_F32: f32 = 1.230_174_1;
29/// Inverse irreversible 9/7 scaling factor, computed the same way as existing paths.
30pub const DWT97_INV_KAPPA_F32: f32 = 1.0 / DWT97_KAPPA_F32;
31/// Historical high-pass synthesis scale used by `OpenJPEG`.
32///
33/// `OpenJPEG` pairs this value with uncompensated irreversible subband step
34/// sizes.  The product is intentionally not identical to `2 / KAPPA`.
35pub const IDWT97_OPENJPEG_TWO_INV_KAPPA_F32: f32 = 1.625_732_4;
36
37/// Inverse 9/7 alpha step used by synthesis paths.
38pub const IDWT97_NEG_ALPHA_F32: f32 = -DWT97_ALPHA_F32;
39/// Inverse 9/7 beta step used by synthesis paths.
40pub const IDWT97_NEG_BETA_F32: f32 = -DWT97_BETA_F32;
41/// Inverse 9/7 gamma step used by synthesis paths.
42pub const IDWT97_NEG_GAMMA_F32: f32 = -DWT97_GAMMA_F32;
43/// Inverse 9/7 delta step used by synthesis paths.
44pub const IDWT97_NEG_DELTA_F32: f32 = -DWT97_DELTA_F32;
45
46/// Forward irreversible 9/7 lifting step alpha at f64 precision.
47pub const DWT97_ALPHA_F64: f64 = -1.586_134_342_059_924;
48/// Forward irreversible 9/7 lifting step beta at f64 precision.
49pub const DWT97_BETA_F64: f64 = -0.052_980_118_572_961;
50/// Forward irreversible 9/7 lifting step gamma at f64 precision.
51pub const DWT97_GAMMA_F64: f64 = 0.882_911_075_530_934;
52/// Forward irreversible 9/7 lifting step delta at f64 precision.
53pub const DWT97_DELTA_F64: f64 = 0.443_506_852_043_971;
54/// Irreversible 9/7 scaling factor at f64 precision.
55pub const DWT97_KAPPA_F64: f64 = 1.230_174_104_914_001;
56/// Inverse irreversible 9/7 scaling factor at f64 precision.
57pub const DWT97_INV_KAPPA_F64: f64 = 1.0 / DWT97_KAPPA_F64;
58
59/// Output columns one block of the tiled horizontal inverse-DWT kernel writes.
60pub const IDWT_HORIZONTAL_TILE_COLUMNS: u32 = 120;
61/// Columns read on each side of a horizontal tile for the lifting steps.
62pub const IDWT_HORIZONTAL_TILE_HALO: u32 = 4;
63/// Rows one block of the tiled horizontal inverse-DWT kernel processes.
64pub const IDWT_HORIZONTAL_TILE_ROWS: u32 = 16;
65/// Threads per tiled horizontal block: one per tile and halo column.
66pub const IDWT_HORIZONTAL_TILE_THREADS: u32 =
67    IDWT_HORIZONTAL_TILE_COLUMNS + 2 * IDWT_HORIZONTAL_TILE_HALO;
68/// Columns one block of the vertical strip inverse-DWT kernels processes.
69pub const IDWT_VERTICAL_STRIP_COLUMNS: u32 = 8;
70/// Columns one block of the fused final-vertical RGB8 store writes.
71pub const FUSED_VERTICAL_TILE_COLUMNS: u32 = 32;
72/// Rows one block of the fused final-vertical RGB8 store writes.
73pub const FUSED_VERTICAL_TILE_ROWS: u32 = 32;
74/// Threads per fused final-vertical RGB8 store block.
75pub const FUSED_VERTICAL_TILE_THREADS: u32 = 256;
76
77/// Whole-sample symmetric extension used at wavelet boundaries: reflects
78/// `index` about the first and last samples into `0..len`. `len` must be at
79/// least 2; a one-sample signal has no reflection period.
80#[must_use]
81#[expect(
82    clippy::cast_possible_truncation,
83    clippy::cast_sign_loss,
84    reason = "every returned index lies in 0..len, which fits u32"
85)]
86pub const fn reflect_index(index: i64, len: u32) -> u32 {
87    let len = len as i64;
88    let positive = if index < 0 { -index } else { index };
89    if positive < len {
90        return positive as u32;
91    }
92    let period = 2 * (len - 1);
93    if positive <= period {
94        return (period - positive) as u32;
95    }
96    let wrapped = positive % period;
97    if wrapped < len {
98        wrapped as u32
99    } else {
100        (period - wrapped) as u32
101    }
102}
103
104#[cfg(test)]
105mod tests {
106    use super::*;
107
108    /// Reflect one edge at a time until `index` lands inside the signal.
109    fn reflect_by_steps(mut index: i64, len: u32) -> u32 {
110        let len = i64::from(len);
111        while index < 0 || index >= len {
112            index = if index < 0 {
113                -index
114            } else {
115                2 * (len - 1) - index
116            };
117        }
118        u32::try_from(index).expect("reflected index is in range")
119    }
120
121    #[test]
122    fn reflect_index_matches_repeated_edge_reflection() {
123        for len in 2..=40 {
124            for index in -300..=300 {
125                assert_eq!(
126                    reflect_index(index, len),
127                    reflect_by_steps(index, len),
128                    "index {index}, len {len}"
129                );
130            }
131        }
132    }
133
134    const MAX_LEVELS_FOR_U32_GEOMETRY: u8 = max_decomposition_levels(u32::MAX, u32::MAX);
135
136    #[test]
137    fn maximum_decomposition_level_compatibility_export_is_const() {
138        assert_eq!(MAX_LEVELS_FOR_U32_GEOMETRY, 31);
139    }
140
141    #[test]
142    fn f32_constants_match_existing_backend_rounding() {
143        assert_eq!(DWT97_ALPHA_F32.to_bits(), (-1.586_134_3f32).to_bits());
144        assert_eq!(DWT97_BETA_F32.to_bits(), (-0.052_980_117f32).to_bits());
145        assert_eq!(DWT97_GAMMA_F32.to_bits(), 0.882_911_1f32.to_bits());
146        assert_eq!(DWT97_DELTA_F32.to_bits(), 0.443_506_87f32.to_bits());
147        assert_eq!(DWT97_KAPPA_F32.to_bits(), 1.230_174_1f32.to_bits());
148        assert_eq!(
149            DWT97_INV_KAPPA_F32.to_bits(),
150            (1.0f32 / 1.230_174_1f32).to_bits()
151        );
152    }
153
154    #[test]
155    fn inverse_constants_are_exact_negations_of_forward_steps() {
156        assert_eq!(IDWT97_NEG_ALPHA_F32.to_bits(), (-DWT97_ALPHA_F32).to_bits());
157        assert_eq!(IDWT97_NEG_BETA_F32.to_bits(), (-DWT97_BETA_F32).to_bits());
158        assert_eq!(IDWT97_NEG_GAMMA_F32.to_bits(), (-DWT97_GAMMA_F32).to_bits());
159        assert_eq!(IDWT97_NEG_DELTA_F32.to_bits(), (-DWT97_DELTA_F32).to_bits());
160    }
161}
j2k_codec_math/dwt.rs

j2k_codec_math/
dwt.rs
