数字图像处理-离散余弦变换的蝶形运算

转自http://blog.csdn.net/sno_guo/article/details/8039659,并自己做了一些内容上的删减。

###为什么要进行变换

空间图像数据通常是很难压缩的：相邻的采样点具有很强的相关性（相互关联的），而且能量一般平均分布在一幅图像中，从而要想丢掉某些数据和降低数据精度而不明显影响图像质量，就要选择合适的变换，方法，使图像易于被压缩。适合压缩的变换方法要有这样几个性质：

1.可以聚集图像的能量（将能量集中到少数有意义的数值上）；如下图：

三维DCT示意图

2.可以除去数据之间的相关性（以使丢掉不重要的数据对图像的质量影响很少）。
3.变换方法应该适合用软硬件实现。

###DCT算法的实例 H264视频压缩

###4x4整数DCT推导过程

$y_k=C_k\sum_{n=0}^{N-1}x_ncos\frac{(2n+1)k\pi}{2N}$

将上面的DCT公式转变为矩阵公式，为了说明标准矩阵中整数的变换和反变换，先假设 $$d=\frac{c}{b}$ $。

$A= \begin{bmatrix} a & a & a & a\\\ b & c & -c &-b \\\ a & -a & -a & a \\\ c & -b & b & -c \end{bmatrix}= \begin{bmatrix} a & & & \\\ & b & & \\\ & & a & \\\ & & & b \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 & 1\\\ 1 & d & -d & -1 \\\ 1 & -1 & -1 & 1 \\\ d & -1 & 1 & -d \end{bmatrix}$ $A^T= \begin{bmatrix} a & b & a & c \\\ a & c & -a & -b \\\ a & -c & -a & b \\\ a & -b & a & -c \end{bmatrix}= \begin{bmatrix} 1 & 1 & 1 & d\\\ 1 & d & -1 & -1 \\\ 1 & -d & -1 & 1 \\\ 1 & -1 & 1 & -d \end{bmatrix} \begin{bmatrix} a & & & \\\ & b & & \\\ & & a & \\\ & & & b \end{bmatrix}$

这里两个[]之间是矩阵相乘的关系

我们可以把DCT变换写成下式：

$Y=AXA^T= \begin{bmatrix} a & a & a & a\\\ b & c & -c &-b \\\ a & -a & -a & a \\\ c & -b & b & -c \end{bmatrix} [X] \begin{bmatrix} a & b & a & c \\\ a & c & -a & -b \\\ a & -c & -a & b \\\ a & -b & a & -c \end{bmatrix}$

三个矩阵相乘，就是DCT转换的公式

$$Y=C_f * X * C_f^t$ $

为了保证正交，a,b,c,d取值如下

$a=\frac{1}{2},b=\sqrt{\frac{1}{2}}cos \left( \frac{\pi}{8} \right),c=\sqrt{\frac{1}{2}}cos \left( frac{3\pi}{8} \right),d=\frac{c}{b}$

DCT变换公式可改为:

$Y= \begin{bmatrix} a & & & \\\ & b & & \\\ & & a & \\\ & & & b \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 & 1\\\ 1 & d & -d & -1 \\\ 1 & -1 & -1 & 1 \\\ d & -1 & 1 & -d \end{bmatrix} [X] \begin{bmatrix} 1 & 1 & 1 & d\\\ 1 & d & -1 & -1 \\\ 1 & -d & -1 & 1 \\\ 1 & -1 & 1 & -d \end{bmatrix} \begin{bmatrix} a & & & \\\ & b & & \\\ & & a & \\\ & & & b \end{bmatrix}$

把两边的对称矩阵移到左边可得：(对角阵移项)

$\begin{bmatrix} a^{-1} & & & \\\ & b^{-1} & & \\\ & & a^{-1} & \\\ & & & b^{-1} \end{bmatrix} [Y] \begin{bmatrix} a^{-1} & & & \\\ & b^{-1} & & \\\ & & a^{-1} & \\\ & & & b^{-1} \end{bmatrix}= \\\ Y\otimes \begin{bmatrix} a^{-2} & (ab)^{-1} & a^{-2} & (ab)^{-1} \\\ (ab)^{-1} & b^{-2} & (ab)^{-1} & b^{-2} \\\ a^{-2} & (ab)^{-1} & a^{-2} & (ab)^{-1} \\\ (ab)^{-1} & b^{-2} & (ab)^{-1} & b^{-2} \end{bmatrix}= \begin{bmatrix} 1 & 1 & 1 & 1\\\ 1 & d & -d & -1 \\\ 1 & -1 & -1 & 1 \\\ d & -1 & 1 & -d \end{bmatrix} [X] \begin{bmatrix} 1 & 1 & 1 & d\\\ 1 & d & -1 & -1 \\\ 1 & -d & -1 & 1 \\\ 1 & -1 & 1 & -d \end{bmatrix}$

再移项变形可得：

$Y= \begin{bmatrix} 1 & 1 & 1 & 1\\\ 1 & d & -d & -1 \\\ 1 & -1 & -1 & 1 \\\ d & -1 & 1 & -d \end{bmatrix} [X] \begin{bmatrix} 1 & 1 & 1 & d\\\ 1 & d & -1 & -1 \\\ 1 & -d & -1 & 1 \\\ 1 & -1 & 1 & -d \end{bmatrix} \otimes \begin{bmatrix} a^2 & ab & a^2 & ab \\\ ab & b^2 & ab & b^2 \\\ a^2 & ab & a^2 & ab \\\ ab & b^2 & ab & b^2 \end{bmatrix}$

此时，d的值为0.4142。这样的话，还是实数运算。如我们令 $$d=\frac{1}{2}$ $,则 $$b=\sqrt{\frac{2}{5}}$ $, $$c=\sqrt{\frac{1}{10}}$ $,同样，可以保证矩阵的正交，同时，可以把运算变为整数运算。

$Y=\left( \begin{bmatrix} 1 & 1 & 1 & 1\\\ 1 & \frac{1}{2} & -\frac{1}{2} & -1 \\\ 1 & -1 & -1 & 1 \\\ \frac{1}{2} & -1 & 1 & -\frac{1}{2} \end{bmatrix} [X] \begin{bmatrix} 1 & 1 & 1 & \frac{1}{2}\\\ 1 & \frac{1}{2} & -1 & -1 \\\ 1 & -\frac{1}{2} & -1 & 1 \\\ 1 & -1 & 1 & -\frac{1}{2} \end{bmatrix} \right) \otimes \begin{bmatrix} a^2 & ab & a^2 & ab \\\ ab & b^2 & ab & b^2 \\\ a^2 & ab & a^2 & ab \\\ ab & b^2 & ab & b^2 \end{bmatrix}$

把 $$\frac{1}{2}$ $提到矩阵之外，并与右边的矩阵点乘合并，得:

$Y=(C_fXC_f^T)\otimes E_f \\\= \left( \begin{bmatrix} 1 & 1 & 1 & 1\\\ 2 & 1 & -1 & -2 \\\ 1 & -1 & -1 & 1 \\\ 1 & -2 & 2 & -1 \end{bmatrix} [X] \begin{bmatrix} 1 & 2 & 1 & 1\\\ 1 & 1 & -1 & -2 \\\ 1 & -1 & -1 & 2 \\\ 1 & -2 & 1 & -1 \end{bmatrix} \right) \otimes \begin{bmatrix} 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{bmatrix}$

其中，

$C_f= \begin{bmatrix} 1 & 1 & 1 & 1\\\ 2 & 1 & -1 & -2 \\\ 1 & -1 & -1 & 1 \\\ 1 & -2 & 2 & -1 \end{bmatrix}, E_f= \begin{bmatrix} 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{bmatrix}$

在H264编码的JM编码器中，变换只包括：

$\begin{bmatrix} 1 & 1 & 1 & 1\\\ 2 & 1 & -1 & -2 \\\ 1 & -1 & -1 & 1 \\\ 1 & -2 & 2 & -1 \end{bmatrix} [X] \begin{bmatrix} 1 & 2 & 1 & 1\\\ 1 & 1 & -1 & -2 \\\ 1 & -1 & -1 & 2 \\\ 1 & -2 & 1 & -1 \end{bmatrix}$

后面的点乘实际上是在量化过程中进行，因为后面的点乘还有实数运算，实数运算将不可避免地产生精度误差,而且运算量巨大。而量化本身就会丢失一些信号，因些，这些实数运算放在量化过程中将大大的降低变换的运算率同时又不明显影响精度？[博主sno的解释: 之所有把后面的Ef和这里分割开，是因为在转换和量化这一步，不可避免的会遇到乘法操作，这里说的只用加减和移位操作，是在这三个矩阵乘法大数据量运算之间只使用加减和移位，可以把64次的乘法运算转为只有两次的乘法运算，极大的减少消耗，在量化的时候，再对每个W[i][j]分别进行乘法操作。]

然而，4X4的矩阵运算如果按常规算法的话仍要进行64次乘法运算和48次加法，这将耗费较多的时间，于是在H.264中，有一种改进的算法（蝶形算法）可以减少运算的次数。这种矩阵运算算法构造非常巧妙，利用构造的矩阵的整数性质和对称性，可完全将乘法运算转化为加法运算。

蝶形运算示意图

变换过程在JM编码器中的实现

//  Horizontal transform水平变换，其实就是左乘Cf,

for (j=0; j < BLOCK_SIZE && !lossless_qpprime; j++) {

        for (i=0; i < 2; i++)

{

  i1=3-i;

  m5[i]=img->m7[i][j]+img->m7[i1][j];

  m5[i1]=img->m7[i][j]-img->m7[i1][j];

}

img->m7[0][j]=(m5[0]+m5[1]);

img->m7[2][j]=(m5[0]-m5[1]);

img->m7[1][j]=m5[3]*2+m5[2];

img->m7[3][j]=m5[3]-m5[2]*2;

} // Vertical transform垂直变换，其实就是右乘CfT for (i=0; i < BLOCK_SIZE && !lossless_qpprime; i++) {

        for (j=0; j < 2; j++)

{

  j1=3-j;

  m5[j]=img->m7[i][j]+img->m7[i][j1];

  m5[j1]=img->m7[i][j]-img->m7[i][j1];

}

img->m7[i][0]=(m5[0]+m5[1]);

img->m7[i][2]=(m5[0]-m5[1]);

img->m7[i][1]=m5[3]*2+m5[2];

img->m7[i][3]=m5[3]-m5[2]*2;

}

###分析该蝶形算法

上面的JM代码就是计算下面三个4x4矩阵的过程。

分析一下前两个矩阵的乘法，只分析他们结果矩阵的第一行。有什么办法可以减少运算量呢？首先采用传统方法计算，得到结果：

X[0] = x[00]+x[10]+x[20]+x[30]

X[1] = 2 * x[00]+x[10]-x[20]-2 * x[30]

X[2]= x[00]-x[10]-x[20]+x[30]

X[3] = x[00]-2x[10]+2x[20]-x[30]

计算代价是16次乘法12次加法，考虑到矩阵的1的乘法可以省略，去除8个乘1，还需要8次乘法和12次加法。那么我们再仔细思考他们的相关性，从一般算法意义上来说，可以用空间代价换时间代价，比如设置中间变量来减少计算次数。用不同的颜色把需要重复运算的部分标上，作为中间变量。

X[0] = x[00]+x[10]+x[20]+x[30]

X[1] = 2 * x[00]+x[10]-x[20]-2 * 2 * x[30]

X[2]= x[00]-</font>+x[10]-x[20]+x[30]

X[3] = x[00]-2x[10]+2x[20]-x[30]

那么提取出来的中间变量将是：

x[00]+x[30] x[00]-x[30] x[10]+x[20]

x[10]-x[20]

存储了这四个中间变量，我们对比看看蝶形图，和图中第一层的算式相符合。用这些中间变量来组合，就可以把最终的X[0]..X[3], 计算出来。这样，就把运算量降低到2个乘法和8个加法，剩余的运算就是叠代这个算法。

所以，可以得出以下结论：

这个蝶形图和一般意义的FFT或FDCT蝶形图不同，是对H.264在整数DCT基础上的具体算法优化，只对于以上Cf矩阵。
计算过程是把上面的三个4x4矩阵乘法分成两两矩阵相乘。再把残差矩阵和后来的中间结果Cf x X一行行分别输入蝶形图进行一维整数DCT计算。
蝶形图优化思想就是提取矩阵的相关部分，定义中间变量，减少运算次数。［sno博主添加：这里很重要，之所以使用蝶形算法，就是统计一下多次运行的规律，不重复计算相关的量。］
sno: 矩阵相乘 A*B 就是把A矩阵当前元素所在的这一行乘以B矩阵当前元素所在的这一列。代码实现如下


	int  matrix_multiply(int a[4][4],int b[4][4],int c[4][4])
	{
		int
    	int i,j,k;
    	for(i=0;i<4;i++)
            for(j=0;j<4;j++)
            {
            	c[i][j]=0;
                for(k=0;k<4;k++)
                	c[i][j]+=a[i][k]*b[k][j];  ///c[i][j]=a的i这一行*b的j这一列的总和 ,这里是累计和                          
            }
        return 0;        
	}

菊长的菊花田

Qinghuai,Jiecao,Neihan,Xiaxian,we got them all!

数字图像处理-离散余弦变换的蝶形运算

关于我

站内搜索

版权声明

近期文章

分类目录

标签

良师益友

存档