图像算法中的二维傅立叶变换(DCT)及蝶形算法(butterfly algorithm)

今天讲个比较基础的知识，图像算法中的二维傅立叶变换(DCT)及蝶形算法(butterfly algorithm)。

关于离散2维Fourier的基础知识这里就不讲了，需要的可以到这里看一下，

fouri小短裙er.eng.hmc.edu/e101/lectures/Image_Processing/node6.html

我们这里从其展开式开始，解释傅立叶变换在图像压缩（如JPEG）中的算法原理和方法。因为傅立叶变换在图像视频压缩中比重还是比较大的，而快速算法在现代计算机科学中又是重中之重。

注意，图像（也包括视频，因为视频是由一帧帧的图像组成）在处理中，并不是整体进行傅立叶变换，而往往会被分成小块处理，最多的情况是４Ｘ４的像素块，或８Ｘ８的像素块。分成这样小块处理的原因也是因为要达到计算量和质量的平衡。

下面我们以4x4的像素块为例，分析一下二维傅立叶变换(DCT)及蝶形算法(butterfly algorithm)的原理和过程，

考虑

\left\{ \begin{array}{lr} \frac{1}{2}Cos(0) = \frac{1}{2}\\ \\ \frac{1}{2}Cos(\frac{2\pi}{8}) = \frac{1}{2\sqrt{2}} \\ \\ \frac{1}{2}Cos(\frac{3\pi}{8}) = \frac{1}{2}Sin(\frac{\pi}{8})\\ \\ \frac{1}{2}Cos(\frac{5\pi}{8}) = -\frac{1}{2}Sin(\frac{\pi}{8})\\ \\ \frac{1}{2}Cos(\frac{7\pi}{8}) = -\frac{1}{2}Cos(\frac{\pi}{8}) \end{array} \right\}

我们可以得到这样一个4x4矩阵的DCT变换矩阵，

\left( \begin{array}{cccc} \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\ \\ \frac{1}{2} \cos \left(\frac{\pi }{8}\right) & \frac{1}{2} \sin \left(\frac{\pi }{8}\right) & -\frac{1}{2} \sin \left(\frac{\pi }{8}\right) & -\frac{1}{2} \cos \left(\frac{\pi }{8}\right) \\ \\ \frac{1}{2 \sqrt{2}} & -\frac{1}{2 \sqrt{2}} & -\frac{1}{2 \sqrt{2}} & \frac{1}{2 \sqrt{2}} \\ \\ \frac{1}{2} \sin \left(\frac{\pi }{8}\right) & -\frac{1}{2} \cos \left(\frac{\pi }{8}\right) & \frac{1}{2} \cos \left(\frac{\pi }{8}\right) & -\frac{1}{2} \sin \left(\frac{\pi }{8}\right) \\ \end{array} \right)

如果你从傅立叶的标准定义去推导的话，得到的就应该是这个结果。用mathematica的命令FourierDCTMatrix[4]，得到的也是这个样子的表达式。然而大多数的文献不会写在这种形式，第一个系数往往会乘以一个参数 \frac{1}{\sqrt{2}} ，所以展开式也往往不会写成标准形式，而是会写成这样（请注意那个系数C），

F_{uv} = C \sum_{i=0}^{N-1}\sum_{j=0}^{N国际交换生-1} f_{ij} \cos \frac{(2i+1)u\pi}{2N} \cos \frac{(2j+1)v\pi}{2N}

C= \left\{ \begin{array}{lr} \frac{1}{2} \quad (x=0) \\ \\ 1 \quad (x \ne 0) \end{array} \right\}

更多的还会把所有的系数另外再乘以一个（这里我们不考虑） \frac{1}{\sqrt{2}} = \frac{2}{\sqrt{N}} \quad (N=4) 。这样做的目的是使得矩阵可以完全正交，即满足： \textbf{W}^T\textbf{W} = \textbf{W}\textbf{W}^T = \textbf{I} ，从而方便处理和减少计算量。不过这样做就不标准了，为有什么影响呢？我们看一下常数Ｃ的物理意义。

这个常数项直接决定直流系数的大小，比如说在图片的亮度处理中就是亮度的起始点，也就是对应亮度起伏为0的那个幅度。在图大抉择像处理中，这种起始点的变化或比例上的缩放是随处可见的，通常这种变化会在后面的量化中加以考虑（这里涉及的内容比较多，进一步的理解建议参考wiki百科，里面列出了不少相关参考资料：en.wikipedia/wiki/Discrete_cosine_transform ），而且，我们往往需要得到的是图像中的起伏和相对变化的情况，所以对处理的特征结果影响不大。

这样的话，我们就不用顾及形式，可以直接考虑一般情况下使用的DCT变换矩阵，这也是大多数资料里面列出的，对图像变换所采用的形式，

\left( \begin{array}{cccc} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & \fr小荷才露尖尖角ac{1}{\sqrt{2}} \\ \\ 室内消火栓 \frac{\cos \left(\frac{\pi }{8}\right)}{\sqrt{2}} & \frac{\sin \left(\frac{\pi }{8}\r农村小洋楼ight)}{\sqrt{2}} & -\frac{排卵时间\sin \left(\frac{\pi }{8}\right)}{\sqrt{2}} & -\fr杭州我爱我家ac{\cos \left(\frac{\pi }{8}\right)}{\sqrt{2}} \\ \\ \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} & \frac{1}{2} \\ \\ \frac{\sin \left(\frac{\pi }{8}\right)}{\sqrt{2}} & -\frac{\cos \le环境影响评价工程师报考条件ft(\frac{\pi }{8}\right)}{\sqrt{2}} & \frac{\cos \left(\frac{\pi }{8}\right)}{\sqrt{2}} & -\frac{\sin \left(\frac{\pi }{8}\right)}{\sqrt{2}} \\ \end{array} \right)

令　

\left\{ \begin{array}{lr} a = \frac{1}{\sqrt{2}}Cos(0) = \frac{1}{\sqrt{2}}\\ 陈克礼\\ b = \frac{1}{\sqrt{2}}Cos(\frac{\pi}{8})\\ \\ c = \frac{1}{\sqrt{2}}Cos(\frac{3\pi}{8}) = \frac{1}{\sqrt{2}}Sin(\frac{\pi}{8}) \end{array} \right\}

上面的矩阵就可以简化为下面这种形式，

\textbf{A} = \left( \begin{array}{cccc} a & a & a & a \\ \\ b & c & -c & -b \\ \\ a & -a & -a & a \\ \\ c & -b & b & -c \e汽车椅套nd{array} \right)

不难证明（参考[1]），对一个任一个4x4的矩阵X，我们可以把运算写成下面这样的形式，

\textbf{A}^T \textbf{X} \textbf{A}^T = \textbf{X} \bullet \left( \begin{array}{cccc} a^2 & \frac{ab}{2} & a^2 & \frac{ab}{2} \\ \\ \frac{ab}{2} & \frac{b^2}{4} & \frac{ab}{2} & \frac{b^2}{4} \\ \\ a^2 & \frac{ab}{2} & a^2 & \frac{a异国恋b}{2} \\ \\ \frac{ab}{2} & \frac{b^2}{4} & \frac{ab}{2} & \frac{b^2}{4} \end{array} \right) = \textbf{C}_f \textbf{X} \textbf{C}_f^T \bullet \textbf{E}_f

其中， \bullet 表示点乘，

\textbf{C}_f = \left( \begin{array}{cccc} 1&1&1&1 \\ \\ 2&1&-1&-2\\ \\ 1&-1&-1&1\\ \\ 1&-2&2&-1 \end{array} \right)

\textbf{E}_f = \left( \begin{array}{cccc} a^2 & \frac{ab}{2} & a^2 & \frac{ab}{2} \\ \\ \frac{ab}{2} & \frac{b^2}{4} & \frac{ab}{2} & \frac{b^2}{4} \\ \\ a^2 & \frac{ab}{2} & a^2 & \frac{ab}{2} \\ \\ \frac{ab}{2} & \frac{b^2}{4} & \frac{ab}{2} & \frac{b^2}{4} \end{array} \right)

可以看出，此时在编码的过程中只需要计算 \textbf{C}_f \textbf{X} \textbf{C}_f^T ，这些系数简单， \textbf{E}_f 也完全可以事计算好，这样一来，运算量就大大降低了，速度会很快，我相信现在的计算机直接处理应该是没有问题的。不过，通过适当的变换，我们还能得到更快的计算方式，这就是知名蝶形算法(butterfly algorithm)。

因为矩阵全部扩展开来会比较大，我们先考虑其中一列，

\left( \begin{array}{cccc} 1&1&1&1 \\ \\ 2&1&-1&-2\\ \\ 1&针刺伤-1&-1&1\\ \\ 1&-2&2&-1 \end{array} \right) \left( \begin{array}{cccc} x_{01} \\ \\ x_{11} \\ \\ x_{21} \\ \\ x_{31} \end{array} \right) = \left( \begin{array}{cccc} (x_{01}+x_{31})+ (x_{11}+x_{21}) \\ \\ 2(x_{01}-x_{31}) + (x_{11}-x_{21}) \\ \\ (x_{01}+x_{31}) - (x_{11}+x_{21})\\ \\ (x_{01}-x_{31}) - 2 (x_{11}-x_{21}) \end{array} \right)

不难发现，其右边小括号内的每一组合项都出现了２次，我把这4个组合项列在下面

\begin{array}{cccc} (x_{01}+x_{31})\\ (x_{婚戒价格01}-x_{31})\\ (x_{11}+x_{21})\\ (x_{11}-x_{21}) \end{array}

根据上面这个原理得来的结果，就是蝶形算法。其中心思想就是，把那些能组合起来的结果尽量先组合起来，得到一个组合值，当后面出现这样的同样组合时，就不必再去计算这个组合值了；从而避免了同样组合的重复计算。

下面我们解释一下常见到的4x4的蝶形算法图，令

\textbf{X} = \left( \begin{array}{cccc} \text{x0} & \text{x1} & \text{x2} & \text{x3} \\ \text{x4} & \text{x5} & \text{x6} & \text{x7} \\ \text{x8} & \text{x9} & \text{x10} & \text{x11} \什么叫公关\ \text{x12} & \text{x13} & \text{x14} & \text{x15} \\ \end{array} \right)

这次把 \textbf{C}_f \textbf{X} \textbf{C}_f^去皱纹的最好方法T 全部展开，得到下面的形式，注意一行就是一个４x4矩阵中的元素，

第1行：

{{x0+x1+x10+x11+x12+x13+x14+x15+x2+x3+x4+x5+x6+x7+x8+x9,

x1-x10+x13-x14-x2+x5-x6-2 (x11+x15+x3+x7)妙心+2 (x0+x12+x4+x8)+x9,

x0-x1-x10+x11+x12-x13-x14+x15-x2+x3+x4-x5-x6+x7+x8-x9,

x0-x11+x12-x15-x3+x4+2 (x10+x14+x2+x6)-x7+x8-2 (x1+x13+x5+x9)},

第2行：

{2 x0+2 x1-x10-x11-2 x12-2 x13-2 x14-2 x15+2 x2+2 x3+x4+x5+x6+x7-x8-x9,

2 x1+x10-2 x13+2 x14-2 x2+x5-x6-2 (-x11-2 x15+2 x3+x7)+2 (2 x0-2 x12+x4-x8)-x9,

2 x0-2 x1+x10-x11-2 x12+2 x13+2 x14-2 x15-2 x2+2 x3+x4-x5-x余额宝有风险吗6+x7-x8+x9,

2 x0+x11-2 x12+2 x15-2 x3+x4+2 (-x10-2 x14+2 x2+x6)-x7-x8-2 (2 x1-2 x13+x5-x9)},

第3行：

{x0+x1-x10-x11+x12+范伟金马奖x13+x14+x15+x2+x3-x4-x5-x6-x7-x8-x9,

x1+x10+x13-x14-x2-x5+x6-2 (-x11+x15+x3-x7)+2 (x0+x12-x4-x8)-x9,

x0-x1+x10-x11+x12-x13-x14+x15-x2+x3-x4+x5+x6-x7-x8+x9,

x0+x11+x12-x15-x3-x4+2 (-x10+x14+x2-x6)+x7-x8-2 (x1+x13-x5-x9)},

第4行：

{x0+x1+2 x10+2 x11-x12-x13-x14-x15+x2+x3-2 x4-2 x5-2 x6-2 x7+2 x8+2 x9,

x1-2 x10-x13+x14-x2-2 x5+2 x6-2 (2 x11-x15+x3-2 x7)+2 (x0-x12-2 x4+2 x8)+2 x9,

x0-x1-2 x10+2 x11-x12+x13+x14-x15-x2+x3-2 x4+2 x5+2 x6-2 x7+2 x8-2 x9,

x0-2 x11-x12+x15-x3-2 x4+2 (2 x10-x14+x2-2 x6)+2 x7+2 x8-2 (x1-x13-2 x5+2 x9)}}

对应的蝶形算法的图形是这样的：

怎么看图？

比如说对X(0)，有4条路径指向X(0),　也就是图中标注红色的(1),(2),(3),(4),

路径(1) +路径 (2) => [(x0+x8)+(x4+x12)]

同理，路径(3) => [(x2+x10) + (x6+x14)]

算法同上，可得路径(4) => {[(x1+x9)+(x5+x13)]+[(x3+x11)+(x7+x15)]}

那么这4条路径的和，就是最终的X(0)＝路径(1) +路径 (2) +路径(3)+路径(４)。

从图中不难看出，后面的计算可以直接采用前面计评价诸葛亮算的结果，一级一级地形成组合值的依赖关系，从而大大减少运算量。

参考资料

[１]　Optimization of 4x4 Integer DCT in H.264/AVC Encoder ，Charles S. Lubobya1, Mqele M. Dlodlo1, Gerhard De. Jager1 and Keith L.Ferguson2