【频域变换】理论与Rust实现



2024年06月01日    Author:Guofei

文章归类: 0x25_CV    文章编号: 1014

版权声明:本文作者是郭飞。转载随意,标明原文链接即可。本人邮箱
原文链接:https://www.guofei.site/2024/06/01/cv6.html


DCT

二维 dct 广泛用于图片频域相关算法,例如 JPEG 压缩

一维离散余弦变换(DCT-II)公式 $F(u) = \alpha(u) \sum\limits_{x=0}^{N-1} f(x) \cos \left( \dfrac{(2x + 1)u\pi}{2N} \right)$

其中,归一化系数 \(\alpha(u) = \begin{cases} \sqrt{\frac{1}{N}} & \text{if } u = 0 \\ \sqrt{\frac{2}{N}} & \text{otherwise} \end{cases}\)

一维离散余弦逆变换(DCT-III)公式 $f(x) = \sum\limits_{u=0}^{N-1} \alpha(u) F(u) \cos \left( \dfrac{(2x + 1)u\pi}{2N} \right)$

二维 DCT 可以由一维 DCT 得出,步骤如下:

  1. 步骤1:对每行应用一维DCT
  2. 步骤2:对每列应用一维DCT

代码

use std::f64::consts::PI;

/// 计算 DCT 的归一化系数 C(u) 或 C(v)
fn c(u: usize, n: usize) -> f64 {
    if u == 0 {
        (1.0 / n as f64).sqrt()
    } else {
        (2.0 / n as f64).sqrt()
    }
}

/// 一维 DCT-II
fn dct_1d(vector: &[f64]) -> Vec<f64> {
    let n = vector.len();
    let mut result = vec![0.0; n];
    for u in 0..n {
        let mut sum = 0.0;
        for x in 0..n {
            sum += vector[x] * ((2 * x + 1) as f64 * u as f64 * PI / (2.0 * n as f64)).cos();
        }
        result[u] = c(u, n) * sum;
    }
    result
}

/// 一维 IDCT-III
fn idct_1d(vector: &[f64]) -> Vec<f64> {
    let n = vector.len();
    let mut result = vec![0.0; n];
    for x in 0..n {
        let mut sum = 0.0;
        for u in 0..n {
            sum += c(u, n) * vector[u] * ((2 * x + 1) as f64 * u as f64 * PI / (2.0 * n as f64)).cos();
        }
        result[x] = sum;
    }
    result
}

/// 二维 DCT
fn dct_2d(matrix: &[Vec<f64>]) -> Vec<Vec<f64>> {
    let n = matrix.len();
    let mut intermediate = vec![vec![0.0; n]; n];
    let mut dct_matrix = vec![vec![0.0; n]; n];

    // 对每一行应用一维 DCT
    for i in 0..n {
        intermediate[i] = dct_1d(&matrix[i]);
    }

    // 对每一列应用一维 DCT
    for j in 0..n {
        let column: Vec<f64> = intermediate.iter().map(|row| row[j]).collect();
        let dct_col = dct_1d(&column);
        for i in 0..n {
            dct_matrix[i][j] = dct_col[i];
        }
    }

    dct_matrix
}

/// 二维 IDCT
fn idct_2d(matrix: &[Vec<f64>]) -> Vec<Vec<f64>> {
    let n = matrix.len();
    let mut intermediate = vec![vec![0.0; n]; n];
    let mut idct_matrix = vec![vec![0.0; n]; n];

    // 对每一列应用一维 IDCT
    for j in 0..n {
        let column: Vec<f64> = matrix.iter().map(|row| row[j]).collect();
        let idct_col = idct_1d(&column);
        for i in 0..n {
            intermediate[i][j] = idct_col[i];
        }
    }

    // 对每一行应用一维 IDCT
    for i in 0..n {
        idct_matrix[i] = idct_1d(&intermediate[i]);
    }

    idct_matrix
}

#[test]
fn tst3() {
    // 示例 4x4 矩阵(例如灰度图像块)
    let matrix = vec![
        vec![52.0, 55.0, 61.0, 66.0],
        vec![70.0, 61.0, 64.0, 73.0],
        vec![63.0, 59.0, 55.0, 90.0],
        vec![67.0, 61.0, 68.0, 104.0],
    ];

    println!("原始矩阵:");
    for row in &matrix {
        println!("{:?}", row);
    }

    // 计算二维 DCT
    let dct = dct_2d(&matrix);
    println!("\nDCT 矩阵:");
    for row in &dct {
        println!("{:?}", row);
    }

    // 计算二维 IDCT
    let idct = idct_2d(&dct);
    println!("\nIDCT 矩阵:");
    for row in &idct {
        println!("{:?}", row);
    }
}

矩阵化

发现事实:

  1. 用到的 cos 值是可以预先计算的,可以做一个 cosine_table
  2. 归一化系数 alpha 也是可以预先计算的
  3. 进而发现循环可以转化为矩阵积

看 dct 公式:$F(u) = \alpha(u) \sum\limits_{x=0}^{N-1} f(x) \cos \left( \dfrac{(2x + 1)u\pi}{2N} \right)$,后半部分的 $\sum$ 可以写成矩阵积,而前半部分的乘法可以写成 矩阵元素积(Hadamard),如下:

\[[F(0), F(1),..., F(N)] \\ = ([f(0), f(1), ..., f(N)] \left[ \begin{array}{l} \cos \dfrac{(2\cdot 0 + 1)\cdot 0\pi}{2N} & \cos \dfrac{(2\cdot 0 + 1)\cdot 1\pi}{2N} & ... & \cos \dfrac{(2\cdot 0 + 1)\cdot u\pi}{2N} & ... & \cos \dfrac{(2\cdot 0 + 1)\cdot N\pi}{2N} \\ \cos \dfrac{(2\cdot 1 + 1)\cdot 0\pi}{2N} & \cos \dfrac{(2\cdot 1 + 1)\cdot 1\pi}{2N} & ... & \cos \dfrac{(2\cdot 1 + 1)\cdot u\pi}{2N} & ... & \cos \dfrac{(2\cdot 1 + 1)\cdot N\pi}{2N} \\ ...&...&...&...&...&...\\ \cos \dfrac{(2x + 1)\cdot 0\pi}{2N} & \cos \dfrac{(2x + 1)\cdot 1\pi}{2N} & ... & \cos \dfrac{(2x + 1)\cdot u\pi}{2N} & ... & \cos \dfrac{(2x + 1)\cdot N\pi}{2N} \\ ...&...&...&...&...&...\\ \cos \dfrac{(2N + 1)\cdot 0\pi}{2N} & \cos \dfrac{(2N + 1)\cdot 1\pi}{2N} & ... & \cos \dfrac{(2N + 1)\cdot u\pi}{2N} & ... & \cos \dfrac{(2N + 1)\cdot N\pi}{2N} \\ \end{array}\right] ) \\ \circ [\sqrt{\frac{1}{N}},\sqrt{\frac{2}{N}},...,\sqrt{\frac{2}{N}}]\]

说明

  • $\circ$ 表示对应元素积(Hadamard 积)
  • 因为 $\alpha$ 的特殊性,这里面的 $\circ$ 运算实际上是初等列变换,是每列乘以常数。
    • 这对于 2d-dct 也成立

上面的式子可以简写为:

  • dct: $F_{1\times n} = (f_{1\times n}\times \cos_{n\times n}) \circ \alpha_{1\times n}$ (矩阵运算)
  • idct $f= (F_{1\times n} \circ \alpha_{1\times n})\times \cos_{n\times n}^T$

而把 Hadamard 积 写成矩阵积,如下:

  • dct: $F_{1\times n} = f_{1\times n}\times \cos_{n\times n} \times \Lambda_{n\times n}$
  • idct $f_{1\times n}= F_{1\times n} \times \Lambda_{n\times n} \times \cos_{n\times n}^T$

说明:

  • 其中 $\Lambda = \mathrm{diag}(\alpha(0),\alpha(1),…,\alpha(N))$,是个对角矩阵
  • 然后我们发现,可以证明 $\cos_{n\times n} \times \Lambda_{n\times n}$ 是一个正交矩阵。
  • 若记为 $E = \cos_{n\times n} \times \Lambda_{n\times n}$,则有 $EE^T=I$
  • 进而转化为
    • dct: $F=f\times E$
    • idct: $f=F\times E^T$

对于 2d-dct,是先对每一行做 dct,然后对每一列做 dct:

  • 对每一行做dct:$T_1=f_{k\times n} D_{n\times n}$
  • 对每一列做dct,我们借用对每一行做 dct 的方式,先转置,对每一行做dct,然后再转置回来:$T_2=(T_1^T D_{k\times k})^T$
  • $T_2 = D_{k\times k}^T f_{k\times n} D_{n\times n}$
  • 因此:
    • dct-2d:$D_{k\times k}^T f_{k\times n} D_{n\times n}$
    • idct-2d: $D_{k\times k} f_{k\times n} D_{n\times n}^T$
  • 说明
    • 上面的 $D_{k\times k}$ 和 $D_{n\times n}$ 是不同的矩阵。
    • 如果 f 是方阵(k=n),那么 $D_{k\times k} = D_{n\times n}$

您的支持将鼓励我继续创作!