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:对每行应用一维DCT
- 步骤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);
}
}
矩阵化
发现事实:
- 用到的 cos 值是可以预先计算的,可以做一个
cosine_table
- 归一化系数
alpha
也是可以预先计算的 - 进而发现循环可以转化为矩阵积
看 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}$