j2k_codec_math/dwt/
linearized53.rs

1// SPDX-License-Identifier: MIT OR Apache-2.0
2
3/// Maximum number of input taps contributing to one linearized 5/3 output.
4pub const DWT53_MAX_LINEAR_TAPS: usize = 5;
5/// Maximum number of input taps contributing to one high-pass 5/3 output.
6pub const DWT53_MAX_HIGH_LINEAR_TAPS: usize = 3;
7
8/// Linearized 5/3 output band.
9#[derive(Clone, Copy, Debug, Eq, PartialEq)]
10pub enum Dwt53Band {
11    /// Low-pass output samples.
12    Low,
13    /// High-pass output samples.
14    High,
15}
16
17/// One non-zero input contribution to a linearized 5/3 output sample.
18#[derive(Clone, Copy, Debug, PartialEq)]
19pub struct Dwt53LinearTap {
20    sample_index: usize,
21    weight: f64,
22}
23
24impl Dwt53LinearTap {
25    const ZERO: Self = Self {
26        sample_index: 0,
27        weight: 0.0,
28    };
29
30    /// Input sample index.
31    #[must_use]
32    pub const fn sample_index(self) -> usize {
33        self.sample_index
34    }
35
36    /// Weight applied to the input sample.
37    #[must_use]
38    pub const fn weight(self) -> f64 {
39        self.weight
40    }
41}
42
43/// Fixed-capacity sparse row for one linearized 5/3 output sample.
44#[derive(Clone, Copy, Debug, PartialEq)]
45pub struct Dwt53LinearRow {
46    taps: [Dwt53LinearTap; DWT53_MAX_LINEAR_TAPS],
47    len: usize,
48}
49
50impl Dwt53LinearRow {
51    const EMPTY: Self = Self {
52        taps: [Dwt53LinearTap::ZERO; DWT53_MAX_LINEAR_TAPS],
53        len: 0,
54    };
55
56    /// Non-zero taps in increasing input-index order.
57    #[must_use]
58    pub fn taps(&self) -> &[Dwt53LinearTap] {
59        &self.taps[..self.len]
60    }
61
62    fn push(&mut self, sample_index: usize, weight: f64) -> Option<()> {
63        if weight == 0.0 {
64            return Some(());
65        }
66        let tap = self.taps.get_mut(self.len)?;
67        *tap = Dwt53LinearTap {
68            sample_index,
69            weight,
70        };
71        self.len += 1;
72        Some(())
73    }
74}
75
76/// Derive one conventional linearized 5/3 analysis row in constant work.
77///
78/// Returns `None` when `output_index` is outside the selected band's output
79/// extent. Symmetric boundary extension is folded into the returned weights.
80#[must_use]
81pub fn linearized_dwt53_row(
82    sample_len: usize,
83    band: Dwt53Band,
84    output_index: usize,
85) -> Option<Dwt53LinearRow> {
86    let output_len = match band {
87        Dwt53Band::Low => low_len(sample_len),
88        Dwt53Band::High => high_len(sample_len),
89    };
90    if output_index >= output_len {
91        return None;
92    }
93
94    let parity = usize::from(matches!(band, Dwt53Band::High));
95    let center = output_index.checked_mul(2)?.checked_add(parity)?;
96    let radius = if matches!(band, Dwt53Band::Low) { 2 } else { 1 };
97    let first = center.saturating_sub(radius);
98    let last = center.saturating_add(radius).min(sample_len - 1);
99    let mut row = Dwt53LinearRow::EMPTY;
100    for sample_index in first..=last {
101        let weight = match band {
102            Dwt53Band::Low => low_weight(sample_len, output_index, sample_index),
103            Dwt53Band::High => high_weight(sample_len, output_index, sample_index),
104        };
105        row.push(sample_index, weight)?;
106    }
107    Some(row)
108}
109
110const fn low_len(sample_len: usize) -> usize {
111    sample_len / 2 + sample_len % 2
112}
113
114const fn high_len(sample_len: usize) -> usize {
115    sample_len / 2
116}
117
118fn low_weight(sample_len: usize, low_index: usize, sample_index: usize) -> f64 {
119    let even_index = low_index * 2;
120    let mut weight = delta(sample_index, even_index);
121    let high_count = high_len(sample_len);
122    let left_high = low_index.checked_sub(1);
123    let right_high = (low_index < high_count).then_some(low_index);
124    match (left_high, right_high) {
125        (Some(left), Some(right)) => {
126            weight += 0.25 * high_weight(sample_len, left, sample_index);
127            weight += 0.25 * high_weight(sample_len, right, sample_index);
128        }
129        (None, Some(right)) => {
130            weight += 0.5 * high_weight(sample_len, right, sample_index);
131        }
132        (Some(left), None) => {
133            weight += 0.5 * high_weight(sample_len, left, sample_index);
134        }
135        (None, None) => {}
136    }
137    weight
138}
139
140fn high_weight(sample_len: usize, high_index: usize, sample_index: usize) -> f64 {
141    let odd_index = high_index * 2 + 1;
142    let left_even = odd_index - 1;
143    let right_even = if odd_index + 1 < sample_len {
144        odd_index + 1
145    } else {
146        left_even
147    };
148    delta(sample_index, odd_index)
149        - 0.5 * (delta(sample_index, left_even) + delta(sample_index, right_even))
150}
151
152fn delta(left: usize, right: usize) -> f64 {
153    f64::from(left == right)
154}
155
156#[cfg(test)]
157mod tests {
158    use super::{linearized_dwt53_row, Dwt53Band, DWT53_MAX_LINEAR_TAPS};
159
160    #[test]
161    fn symbolic_rows_match_the_canonical_linearized_transform() {
162        let samples = [
163            11.0, -7.0, 23.0, 5.0, -19.0, 31.0, 2.0, 17.0, -13.0, 29.0, 3.0, -5.0, 37.0, -11.0,
164            7.0, 41.0, -17.0,
165        ];
166        for sample_len in 1..=samples.len() {
167            for low_index in 0..sample_len.div_ceil(2) {
168                let row = linearized_dwt53_row(sample_len, Dwt53Band::Low, low_index)
169                    .expect("in-range low row");
170                assert_close(
171                    apply(&row, &samples),
172                    direct_low(&samples[..sample_len], low_index),
173                );
174            }
175            for high_index in 0..sample_len / 2 {
176                let row = linearized_dwt53_row(sample_len, Dwt53Band::High, high_index)
177                    .expect("in-range high row");
178                assert_close(
179                    apply(&row, &samples),
180                    direct_high(&samples[..sample_len], high_index),
181                );
182            }
183        }
184    }
185
186    #[test]
187    fn maximum_jpeg_axis_builds_one_fixed_size_row_per_output() {
188        let sample_len = 65_535usize;
189        let mut row_count = 0usize;
190        let mut tap_count = 0usize;
191        for (band, output_len) in [
192            (Dwt53Band::Low, sample_len.div_ceil(2)),
193            (Dwt53Band::High, sample_len / 2),
194        ] {
195            for output_index in 0..output_len {
196                let row = linearized_dwt53_row(sample_len, band, output_index)
197                    .expect("in-range large-axis row");
198                assert!(row.taps().len() <= DWT53_MAX_LINEAR_TAPS);
199                assert!(row.taps().iter().all(|tap| tap.sample_index() < sample_len));
200                row_count += 1;
201                tap_count += row.taps().len();
202            }
203        }
204        assert_eq!(row_count, sample_len);
205        assert!(tap_count <= sample_len * DWT53_MAX_LINEAR_TAPS);
206    }
207
208    fn apply(row: &super::Dwt53LinearRow, samples: &[f64]) -> f64 {
209        row.taps()
210            .iter()
211            .map(|tap| samples[tap.sample_index()] * tap.weight())
212            .sum()
213    }
214
215    fn assert_close(actual: f64, expected: f64) {
216        assert!(
217            (actual - expected).abs() <= f64::EPSILON,
218            "actual={actual} expected={expected}"
219        );
220    }
221
222    fn direct_low(samples: &[f64], low_index: usize) -> f64 {
223        let current = samples[low_index * 2];
224        let left = low_index
225            .checked_sub(1)
226            .map(|index| direct_high(samples, index));
227        let right = (low_index < samples.len() / 2).then(|| direct_high(samples, low_index));
228        current
229            + match (left, right) {
230                (Some(left), Some(right)) => (left + right) * 0.25,
231                (None, Some(right)) => right * 0.5,
232                (Some(left), None) => left * 0.5,
233                (None, None) => 0.0,
234            }
235    }
236
237    fn direct_high(samples: &[f64], high_index: usize) -> f64 {
238        let odd_index = high_index * 2 + 1;
239        let left = samples[odd_index - 1];
240        let right = samples.get(odd_index + 1).copied().unwrap_or(left);
241        samples[odd_index] - (left + right) * 0.5
242    }
243}
j2k_codec_math/dwt/linearized53.rs

j2k_codec_math/dwt/
linearized53.rs
