1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198

use std::ops::Range;

#[cfg(feature="ipp")]
mod dwt1d;

#[cfg(feature="ipp")]
mod dwt2d;

#[cfg(feature="ipp")]
pub use dwt1d::*;

#[cfg(feature="ipp")]
pub use dwt2d::*;

// Filter coefficients were extracted from Tan & Jiang (2019) Digital Signal Processing
// Fundamentals and applications (Third Edition); Table 12.2 (p. 621).
pub mod bank {

    // Haar: 1 / sqrt(2) = 0.7...

    pub const HAAR_LOW : [f32; 2] = [1., 1.];

    pub const HAAR_HIGH : [f32; 2] = [-1., 1.];

    pub const HAAR_LOW_INV : [f32; 2] = [1., 1.];

    pub const HAAR_LOW_HIGH : [f32; 2] = [1., -1.];

    pub const DAUB4_LOW : [f32; 4] = [
        0.48296291314453414337487159986,
        0.83651630373780790557529378092,
        0.22414386804201338102597276224,
        -0.12940952255126038117444941881
    ];

    pub const DAUB4_HIGH : [f32; 4] = [
        -0.12940952255126038117444941881,
        -0.22414386804201338102597276224,
        0.83651630373780790557529378092,
        -0.48296291314453414337487159986
    ];

    /*// const fn implementation if const floating point arithmetic were to become stabilized
    // https://github.com/rust-lang/rust/issues/57241
    pub const fn define_coef<const U : usize>(h_k : &[f32; U], t : usize) -> f32 {
        let ix = ((U as i32 - 1 - t as i32) % U as i32) as usize;
        (-1. as f32).powf(t as f32) * h_k[ix]
    }*/

    /*pub const fn advance_level<const U : usize>(h_k : &[f32; 4]) -> [f32; 4] {
        [
            define_coef(h_k, 0),
            define_coef(h_k, 1),
            define_coef(h_k, 2),
            define_coef(h_k, 3)
        ]
    }*/

    // Derive a difference filter from a coarse filter. This is used for testing purposes only,
    // since coefficients are passed as &'static [f32]. If const floating point arithmetic
    // is ever stabilized, we can use this to define the coefficients instead of using
    // hard-coded values.
    pub fn derive_diff<const U : usize>(h_k : &[f32; U]) -> [f32; U] {
        let mut h_kp1 = h_k.clone();
        for t in 0..U {
            let ix = ((U as i32 - 1 - t as i32) % U as i32) as usize;
            h_kp1[t] = (-1. as f32).powf(t as f32) * h_k[ix];
        }
        h_kp1
    }

}

#[test]
fn advance_level() {
    use crate::dwt::bank::*;
    println!("{:?} {:?}", derive_diff(&DAUBECHIES), DAUB4_HIGH);
    /*println!("Original: {:?}", &bank::DAUBECHIES);
    println!("Apply 1: {:?}", bank::derive_diff(&bank::DAUBECHIES));
    println!("Apply 2: {:?}", bank::derive_diff(&bank::derive_diff(&bank::DAUBECHIES)));
    println!("Apply 3: {:?}", bank::derive_diff(&bank::derive_diff(&bank::derive_diff(&bank::DAUBECHIES))));
    println!("Apply 4: {:?}", bank::derive_diff(&bank::derive_diff(&bank::derive_diff(&bank::derive_diff(&bank::DAUBECHIES)))));*/
}

#[derive(Clone, Copy, Debug)]
pub enum Basis {

    // Centered or not
    Haar(bool),

    // k = 4..20 even; centered or not
    Daubechies(usize, bool),

    // k = 103, 105, 202..208 even, 301..309 odd; centered or not
    BSpline(usize, bool)

}

impl Basis {

    fn len(&self) -> usize {
        match self {
            Basis::Daubechies(4, _) => 4,
            _ => unimplemented!()
        }
    }

    fn coefficients(&self) -> (&'static [f32; 4], &'static [f32; 4]) {
        match self {
            Basis::Daubechies(4, _) => (&bank::DAUB4_LOW, &bank::DAUB4_HIGH),
            _ => panic!("Unimplemented basis {:?}", self)
        }
    }

    /*fn taps(&self, n_levels : usize) -> Vec<([f32; 4], [f32; 4])> {
        /*match self {
            Basis::Daubechies(4, _) => {
                let mut coefs = Vec::new();
                let low = self.coefficients().clone();
                let high = bank::derive_diff(&low);
                coefs.push((low, high));
                for i in 1..(n_levels-1) {
                    let low = bank::advance_level(&coefs[i-1].0);
                    let high = bank::advance_level(&coefs[i-1].1);
                    coefs.push((low, high));
                }
                coefs
            },
            _ => unimplemented!()
        }*/
        unimplemented!()
    }*/

}

// Checks if the vector with len given by the argument is a valid DWT argument.
// Returns true if it is an integer power of two. Returns false otherwise. Integer
// powers of two are required for DWT buffers because those are the only
// values that can be divided by two recursively until we arrive at a coarse
// coefficient buffer of size n=filter_len. For any other N even, we might only
// be able to perform one or a few steps of the decomposition, then fall into a value
// with odd number of elements, leaving the transform incomplete. Eventually, this might
// be all that the user needs, but we impose this restriction for now.
pub fn is_valid_dwt_len(len : usize) -> bool {
    len > 1 && (len as f64).log2().fract() == 0.0
}

// Checks if the side of the square image of given buffer len has a valid lenght for dwt decomposition.
pub fn is_valid_dwt_side(len : usize) -> bool {
    is_valid_dwt_len(side_len(len))
}

// Returns maximum number of levels for image of size n.
pub fn dwt_max_levels(n : usize) -> usize {
    (n as f32).log2() as usize - 1
}

// If cascade data is located in a contiguous buffer from coarse to detail, return the range that index
// the transformed buffer, assuming an original signal size n, and that we wish to index the level lvl.
fn dwt_level_index(n : usize, lvl : usize) -> Range<usize> {
    let ix = n / (2usize).pow(lvl as u32 + 1);
    ix..2*ix
}

fn verify_dwt_dimensions(src_len : usize, coarse_len : usize, detail_len : usize) {
    assert!(src_len / 2 == coarse_len);
    assert!(src_len / 2 == detail_len);
    for len in [src_len, coarse_len, detail_len].iter() {
        assert!(is_valid_dwt_len(*len));
    }
}

fn verify_dwt1d_dimensions(src : &[f32], coarse : &[f32], detail : &[f32]) {
    verify_dwt_dimensions(src.len(), coarse.len(), detail.len());
}

fn verify_dwt2d_dimensions(src : &[f32], coarse : &[f32], detail_x : &[f32], detail_y : &[f32], detail_xy : &[f32]) {
    verify_dwt_dimensions(side_len(src.len()), side_len(coarse.len()), side_len(detail_x.len()));
    verify_dwt_dimensions(side_len(src.len()), side_len(coarse.len()), side_len(detail_y.len()));
    verify_dwt_dimensions(side_len(src.len()), side_len(coarse.len()), side_len(detail_xy.len()));
}

/// Tests len is a power of two and returns the side at the same time.
fn side_len(len : usize) -> usize {
    let sq_len = (len as f32).sqrt();
    assert!(sq_len.fract() == 0.0);
    sq_len as usize
}

/// Gets element at index last-1 and element at index last and return them
/// as a tuple of mutable references.
fn index_and_prev_mut<T>(levels : &mut [T], last : usize) -> (&mut T, &mut T) {
    assert!(last >= 1 && last <= levels.len() - 1);
    assert!(levels.len() >= 2);
    let (curr, prev) = levels[last-1..last+1].split_last_mut().unwrap();
    assert!(prev.len() == 1);
    (&mut prev[0], curr)
}