第3章 图像正交变换
图像的正交变换是先将数字图像由空域变换到频域，在领域中进行处理，以提高图像的处理速度。在频域法的数字图像处理中，一般采用的变换方式都是线性正交变换，因此图像的正交变换是一种非常关键的预处理方法。本章将以实例的方式系统介绍几种常用的图像正交变换，并以VC代码实现。
3.1 基本正交变换
图像的基本正交变换包括离散傅里叶变换、离散余弦变换、离散沃尔什变换等，下面对这几种正交变换分别加以介绍。
3.1.1 离散傅里叶变换
离散傅里叶变换是傅里叶变换的一种，傅里叶变换的良好性质使离散傅里叶变换在分析处理数字图像等方面具有许多优势，得到普遍的应用。
1. 基本原理
（1）傅里叶变换
傅里叶变换在数学上的定义非常严格，其定义如下，设f (x)为x的函数，如果满足下面的狄里赫莱条件：
▓ 具有有限个间隔点。
▓ 具有有限个极点。
▓ 绝对可积。
则傅里叶变换为
其逆变换公式为
其中，x为时域变量；u为频域变量。通常把以上两个公式称为傅里叶变换对。
一维傅里叶变换可以推广到二维函数。即如果二维函数满足狄里赫莱条件，则下面二维傅里叶变换对存在。
傅里叶变换在一维信号处理中已经得到了广泛应用。
（2）离散傅里叶变换
要在数字图像处理中应用傅里叶变换，还需要解决两个问题：一是在数学中进行傅里叶变换的f (x)为连续（模拟）信号，而计算机处理的是离散（数学）信号（图像数据）；二是数学上采用无穷大概念，而计算机只能进行有限次计算。通常，将受这种限制的傅里叶变换称为离散傅里叶变换（DFT）。
1) 对一维连续函数f (x)的傅里叶变换，如式（3-1）、式（3-2），可以用一个离散序列去逼近这一积分。从另一个角度讲，设{f (x) | f (0), f (1), f (2), ⋯, f (N-1)}为一维信号f (x)的N个抽样，则其一维离散傅里叶变换对为
其中，x，u= 0, 1, 2, ⋯, N－1。可见，离散序列的傅里叶变换仍是一个离散的序列。
2) 在数字图像图理中用到的傅里叶正交变换通常是二维离散函数傅里叶变换，所以这里将重点介绍二维离散函数的傅里叶变换。设一幅二维离散图像f (x, y)的大小为M×N，则二维离散图像傅里叶变换定义为
其中，u= 0, 1, 2, ⋯, M－1；v= 0, 1, 2, ⋯, N－1。
其逆变换为
其中，x= 0, 1, 2,⋯, N－1, y= 0, 1, 2, ⋯, N－1。
在实际处理图像中，通常会选择方形阵列，所以二维离散图像傅里叶变换对在M=N的情况下可表示为
其中，x，u= 0, 1, 2, ⋯, N－1；y，v= 0, 1, 2, ⋯, N－1。
3) 二维离散函数傅里叶变换的性质有以下几种。
▓ 可分离性
由式（3-4）容易得出
可分离性说明，二维离散函数傅里叶变换可以由两次一维离散函数傅里叶变换来实现，即先通过f (x, y)求出F (x, v)，进而求出F (u, v)。对其求逆变换，亦然。
▓ 可平移性
二维离散傅里叶变换和逆变换的位移性质，可表示为
通过二维离散傅里叶变换的可平移性，将f (x, y)的傅里叶变换原点（图像的频域原点）移动至N×N频率方阵的中心(N/2, N/2)，这样可以观察到一个完全周期的傅里叶频谱。
▓ 周期性
二维离散傅里叶变换具有周期性，即
F (u,v) = F (u+ aN, v) = F (u,v+ bN) = F (u+ aN, v+ bN) f(x, y) = f (x+ aN, y+ bN)
由此可知，离散函数f (x, y)经过正交变换后得到的F (u, v)是具有周期性的离散函数。反之，亦然。这也是进行快速傅里叶变换的基础。
▓ 离散卷积定理
二维离散傅里叶卷积定理为
其中，F (u, v)为离散函数f (x, y)的离散傅里叶变换；G (u, v)为离散函数g (x, y)的离散傅里叶变换。
根据此定理，在一些涉及相关卷积的运算中，可将求时域内的卷积转化为求频域内的乘积来简化运算过程。
此外，二维离散傅里叶变换也具有线性、旋转性、比例性和相关定理等性质，这些性质在分析及处理数字图像时有着重要的意义。
（3）快速傅里叶变换
由傅里叶变换的公式可以看到，傅里叶变换的计算量很大，可以证明仅对一个长度为N的一维离散函数进行傅里叶变换来说，其运算次数正比于N2，特别是当N较大时，其运算时间将迅速增长，在某种程度上限制了其使用范围。为此，库里（Cooley）和图基（Tukey）在1965年首先提出了一种逐次加倍的快速傅里叶变换算法（FFT），使用该算法可以将运算的次数降到正比于N。当N很大时，FFT大大减小了计算量，缩短了运算时间，提高了运行效率。当然有关FFT的算法还有其他很多种，FFT并不是一种新的变换，它只是离散傅里叶变换的一种算法。这种算法分析离散傅里叶变换中的多余运算，进而消除这些重复的工作，达到了快速运算的目的。
由于二维离散傅里叶变换具有可分离性，可由两次一维离散傅里叶变换计算得到，因此，只需要研究一维离散傅里叶变换的快速算法即可。
式（3-3）可以写成
其中，
称为旋转因子。
式（3-5）的展开式为
可将此展开式表示为矩阵的形式
其中，由Wux构成的矩阵称为W阵或系数矩阵。
由此系数矩阵和
可知，系数W是以N为周期的，也就是说，W阵中很多系数是相同的，不必进行多次重复计算，且由于W的对称性，可以进一步减小运算量。
例如：在N= 4时，W阵为
由W的周期性得W4＝W0，W6＝W2，W9＝W1。
再由W的对称性可得W3＝-W1，W2＝-W0。
所以，在N= 4时，此W阵变为
可见N= 4的W阵中只需计算W0和W1两个系数即可。说明W阵的系数有许多计算工作是重复的，如果把一个离散序列分解成若干短序列，并充分利用旋转因子W的周期性和对称性来计算离散傅里叶变换，便可以简化运算过程，这就是FFT的基本思想。
设N为2的正整数次幂，即
N = 2n n=1, 2, ⋯
又令M为正整数，且N = 2M，则式（3-5）可写为
由
，化简式（3-6）得
现定义
则式（3-7）可表示为
再由
和
，得
由此，可将一个N点的离散傅里叶变换分解成两个N/2短序列的离散傅里叶变换，即分解为偶数和奇数序列的离散傅里叶变换Fe(u)和Fo(u)。
下面，以计算N= 8的离散傅里叶变换为例进行说明，此时n= 3，M= 4
可以看到，式（3-8）中
如图3-1所示，定义由Fe(0)、Fo(0)、F(0)和F(4)所构成的结构为蝶形运算单元，其左方两个节点为输入节点，代表输入值；右方两个节点为输出节点，代表输入值的叠加，运算由左向右进行，斜线旁的
和-
为加权系数。因此，图3-1所代表的含义即为式（3-9）。
图3-1 蝶形运算单元
下面，继续对式（3-8）进行分解，由于Fe(0)、Fo(0)是N= 4，M= 2离散傅里叶变换，所以再对其按照奇偶进行分组。
由于Fee(u)、Feo(u)、Foe(u)和Foo(u)均为两点的离散傅里叶变换，所以不能再继续分解，可以直接由原始数据f(x)来计算。
由此，可以完成对离散傅里叶变换的快速运算，根据以上公式及各级蝶形运算流程图，对一个8点的FFT，首先将其分成4组，每组为2点的FFT，然后再分成两组，每组为4点的FFT，最后根据已知条件计算出8点的FFT。与DFT相比，FFT计算次数大大减少，运算时间缩短，运算效率得到很大程度提高。
2. 算法描述
实现数字图像在屏幕上的傅里叶变换，可以分为以下5个步骤。
[1] 将源图像保存到缓冲区，并记录下缓冲区的地址。
[2] 获取源图像的高度和宽度，计算图像每行的字节数。
[3] 调用傅里叶变换函数，对图像进行傅里叶变换。
[4] 在傅里叶变换函数过程中，调用快速傅里叶变换函数。
[5] 采用蝶形算法对图像进行快速傅里叶变换，最终显示变换后的图像。
3. 编程实现
下面通过编写VC的函数Fourier ( )实现数字图像的傅里叶变换，其具体代码如下：
/＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊ 函数名称：Fourier(CplexNum ＊ pTd, int lWidth, int lHeight, CplexNum ＊ pFd)
＊ 函数参数：
＊ CplexNum ＊ pTd 指向时域值的指针
＊ int lWidth 图像宽度
＊ int lHeight 图像高度
＊ CplexNum ＊ pFd 指向频域值的指针
＊ 函数功能：二维快速傅里叶变换
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊/
void Fourier(CplexNum ＊ pTd, int lWidth, int lHeight, CplexNum ＊ pFd)
{
int j; // 循环控制变量
int i;
int wid=1; // 进行傅里叶变换的宽度和高度，（2的整数次幂）
int hei=1; // 图像的宽度和高度不一定为2的整数次幂
int widpor=0,heiPor=0; // 2的幂数
while(wid ＊ 2 <= lWidth) // 计算进行傅里叶变换的宽度和高度（2的整数次幂）
{
wid ＊= 2;
widpor++;
}
while(hei ＊ 2 <= lHeight)
{
hei ＊= 2;
heiPor++;
}
for(i = 0; i < hei; i++)
{
FastFourierTran(&pTd[wid ＊ i], &pFd[wid ＊ i], widpor);
// x方向进行快速傅里叶变换
}
// pFd中目前存储了pTd经过行变换的结果，为了直接利用FastFourierTran，需要把
// pFd的二维数据转置，再次利用FastFourierTran进行傅里叶行变换（实际上相当于
// 对列进行傅里叶变换）
for(i = 0; i < hei; i++)
for(j = 0; j < wid; j++)
pTd[hei ＊ j + i] = pFd[wid ＊ i + j];
for(j = 0; j < wid; j++)
{
// 对x方向进行快速傅里叶变换，实际上相当于对原来的图像数据进行列方向的傅
// 里叶变换
FastFourierTran(&pTd[j ＊ hei], &pFd[j ＊ hei], heiPor);
}
// pFd中目前存储了pTd经过二维傅里叶变换的结果，但是为了方便列方向的傅里叶
// 变换，对其进行了转置，现在把结果转置回来。
for(i = 0; i < hei; i++)
for(j = 0; j < wid; j++)
pTd[wid ＊ i + j] = pFd[hei ＊ j + i];
memcpy(pTd, pFd, sizeof(CplexNum) ＊ hei ＊ wid );
}
编写VC的函数FastFourierTran ()来实现图像的快速傅里叶变换，其具体代码如下：
/＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊ 函数名称：FastFourierTran(CplexNum ＊ pTd, CplexNum＊ pFd, int power)
＊ 函数参数：
＊ CplexNum ＊ pTd 指向时域数组的指针
＊ CplexNum ＊ pFd 指向频域数组的指针
＊ int power 2的幂数，即迭代次数
＊lt@span b=1> lt@span b=1> 函数功能：用来实现快速傅里叶变换
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊/
void FastFourierTran(CplexNum ＊ pTd, CplexNum ＊ pFd, int power)
{
long i;
// 行循环变量
long j;
// 列循环变量
long dotCount; // 傅里叶变换点数
int k; // 循环变量
int bfsize,p; // 中间变量
double angle; // 角度
CplexNum ＊pWn,＊temReg1,＊temReg2,＊temReg;
dotCount= 1 <<power; // 计算傅里叶变换点数
pWn= new CplexNum[sizeof(CplexNum)＊dotCount/ 2]; // 分配运算所需存储器
temReg1 = new CplexNum[sizeof(CplexNum)＊dotCount];
temReg2 = new CplexNum[sizeof(CplexNum)＊dotCount];
for(i = 0; i < dotCount/ 2; i++) // 计算加权系数
{
angle = -i ＊ pi＊ 2 / dotCount;
pWn[i].re = cos(angle);
pWn[i].im=sin(angle);
}
memcpy(temReg1, pTd, sizeof(CplexNum)＊dotCount); // 将时域点写入temReg1
for(k = 0; k < power; k++) // 采用蝶形算法进行快速傅里叶变换
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (power-k);
for(i = 0; i < bfsize / 2; i++)
{
p = j ＊ bfsize;
temReg2[i+p]=Add(temReg1[i+p],temReg1[i+p+bfsize/2]);
temReg2[i+p+bfsize/2]=Mul(Sub(temReg1[i+p],temReg1[i+p+bfsize/2]),
pWn[i＊(1<<k)]);
}
}
temReg = temReg1;
temReg1 = temReg2;
temReg2 = temReg;
}
for(j = 0; j <dotCount; j++) // 重新排序
{
p = 0;
for(i = 0; i <power; i++)
{
if (j&(1<<i))
{
p+=1<<(power-i-1);
}
}
pFd[j]=temReg1[p];
}
delete pWn; // 释放内存
delete temReg1;
delete temReg2;
}
可以通过在数字图像处理程序的框架窗口上增加菜单【正交变换】|【离散傅里叶变换】，如图3-2所示。对应视图类中的CDImageProcessView:: OnFourier ( )进行函数的调用。
图3-2 正交变换菜单
void CDImageProcessView::OnFourier()
{
CDImageProcessDoc＊ pDoc = GetDocument();
SetCursor(LoadCursor(NULL, IDC_WAIT)); // 光标显示为处理状态
long lSrcLineBytes; // 图像每行的字节数
long lSrcWidth; // 图像的宽度和高度
long lSrcHeight;
LPSTR lpSrcDib; // 指向源图像的指针
LPSTR lpSrcStartBits; // 指向源像素的指针
lpSrcDib= (LPSTR) ::GlobalLock((HGLOBAL) pDoc->GetHObject()); // 锁定DIB
if (pDoc->m_dib.GetColorNum(lpSrcDib) != 256) // 判断是否是8-bpp位图
{
AfxMessageBox(_T ("对不起，不是256色位图！")); // 警告
::GlobalUnlock((HGLOBAL) pDoc->GetHObject()); // 解除锁定
return; // 返回
}
lpSrcStartBits=pDoc->m_dib.GetBits(lpSrcDib);// 找到DIB图像像素起始位置
lSrcWidth= pDoc->m_dib.GetWidth(lpSrcDib); // 获取图像的宽度
lSrcHeight= pDoc->m_dib.GetHeight(lpSrcDib); // 获取图像的高度
lSrcLineBytes=pDoc->m_dib.GetReqByteWidth(lSrcWidth ＊ 8); // 计算图像每行的字节数
double dTmp; // 临时变量
long wid=1,hei=1; // 进行傅里叶变换的宽度和高度，初始化为1
int widpor=0, heiPor=0; // 2的幂数
unsigned char＊ lpSrcUnChr; // 指向像素的指针
long i; // 行循环变量
long j; // 列循环变量
while(wid ＊ 2 <= lSrcWidth) // 计算进行傅里叶变换的宽度和高度（2的整数次幂）
{
wid ＊= 2;
widpor++;
}
while(hei ＊ 2 <= lSrcHeight)
{
hei ＊= 2;
heiPor++;
}
CplexNum ＊pTd = new CplexNum[sizeof(CplexNum)＊wid ＊ hei]; // 分配内存
CplexNum ＊pFd = new CplexNum[sizeof(CplexNum)＊wid ＊ hei]; // 初始化
// 图像数据的宽和高不一定是2的整数次幂，所以pCTData有一部分数据需要补0
for(i=0; i<hei; i++)
{
for(j=0; j<wid; j++)
{
pTd[i＊wid + j].re=0;
pTd[i＊wid + j].im=0;
}
}
for(i=0; i<lSrcHeight; i++) // 把图像数据传给pCTData
{
for(j=0; j<lSrcWidth; j++)
{
lpSrcUnChr=(unsigned char＊)lpSrcStartBits+
lSrcLineBytes＊(lSrcHeight-1-i)+j;
pTd[i＊wid + j].re=＊lpSrcUnChr;
pTd[i＊wid + j].im=0;
}
}
Fourier(pTd, lSrcWidth, lSrcHeight, pFd) ;
for(i=0; i<lSrcHeight; i++)
{
for(j=0; j<lSrcWidth; j++)
{
dTmp = pFd[i ＊ wid + j].re＊ pFd[i ＊ wid + j].re + pFd[i ＊ wid +
j].im＊ pFd[i ＊ wid + j].im;
dTmp = sqrt(dTmp) ; // 为了显示，需要对幅度的大小进行伸缩
dTmp /= 100 ; // 限制图像数据的大小
dTmp = min(dTmp, 255) ;
lpSrcStartBits[i＊lSrcLineBytes +j] = (unsigned char)(int)dTmp;
}
}
for(i=0; i<lSrcHeight; i++)// 为了在屏幕上显示，我们把幅度值大的部分用黑色显示
{
for(j=0; j<lSrcWidth; j++)
{
lpSrcStartBits[i＊lSrcLineBytes +j] = 255- lpSrcStartBits[i＊lSrcLineBytes +j];
}
}
delete pTd; // 释放内存
delete pFd;
pDoc->SetModifiedFlag(TRUE); // 设置修改标记
pDoc->UpdateAllViews(NULL); // 更新视图
SetCursor(LoadCursor(NULL, IDC_ARROW)); // 光标显示为处理状态
}
4. 效果演示
源图像如图3-3a所示，对其进行离散傅里叶变换，结果如图3-3b所示。
图3-3 离散傅里叶变换效果演示
3.1.2 离散余弦变换
离散余弦变换（DCT）也是一种常用的正交变换，它是从一种特殊形式的傅里叶变换转化而来的，当原始信号满足一定的条件时，就可以将傅里叶变换转化成余弦变换，其变换核为余弦函数。
1. 基本原理
（1）一维离散余弦变换
设{f (x) | f (0), f (1), f (2), ⋯, f (N-1)}为离散序列，由一维离散余弦变换定义如下：
其中，F(u)为第u个余弦变换系数；u为广义频率变量；f(x)为时域中N点序列；x, u = 0, 1, 2, ⋯, N－1。
其逆变换为
其中，x, = 0, 1, 2, ⋯, N－1。
另外，一维离散余弦变换的变换核为
其中，x, u= 0, 1, 2, ⋯, N－1。可见一维离散余弦变换的正变换核与逆变换核是相同的。一维离散余弦变换变换式展开整理后，可以写成矩阵的形式F = Gf，其中
同理，一维离散余弦变换逆变换的矩阵形式为
f = GT F
其中，GT是G的转置矩阵。
（2）二维离散余弦变换
通过一维离散余弦变换的定义，很容易将其推广到二维离散余弦变换。设f (x, y)为M×N的数字图像矩阵，则定义二维离散余弦变换为
其逆变换为
其中，x, u= 0, 1, 2, ⋯, M－1；y, v= 0, 1, 2, ⋯, N－1。
二维离散余弦变换的变换核为
其中，x, u= 0, 1, 2, ⋯, M－1；y, v= 0, 1, 2, ⋯, N－1。
对比一维离散余弦变换的矩阵表示，根据二维余弦变换及其逆变换的定义可以得到二维离散余弦变换的矩阵表达式为
F (u, v) = Gf (x, y) GT
其逆变换的矩阵表达式为
f (x, y) = GTF (u, v) G
（3）离散余弦变换的性质
离散余弦变换与傅里叶变换一样，首先具有变换核可分离性，即二维离散余弦变换及逆变换可以通过一维离散余弦变换进行计算；其次，离散余弦变换也有快速算法，可用FFT进行运算；最后，离散余弦变换来自切比雪夫多项式具有完备的正交性质，可广泛地运用于数字图像处理。另外，离散余弦变换的变换阵的基向量还能很好地描述人类语音信号和图像信号的相关特征。因此，在对语音信号和图像信号的变换中，离散余弦变换被认为是一种准最佳变换。
2. 算法描述
实现数字图像在屏幕上的离散余弦变换，可以分为以下5个步骤：
[1] 将源图像保存到缓冲区，并记录下缓冲区的地址。
[2] 获取源图像的高度和宽度，计算图像每行的字节数。
[3] 调用离散余弦变换函数，对图像进行离散余弦变换。
[4] 在离散余弦变换过程中，调用快速离散余弦变换函数，在调用快速离散余弦变换函数中，又调用快速傅里叶变换函数。
[5] 更新显示变换后的图像。
3. 编程实现
通过编写VC的函数CosTran ( )来实现数字图像的离散余弦变换，其具体代码如下：
/＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊ 函数名称：CosTran(LPSTR lpSrcStartBits, long lWidth, long lHeight,long lLineBytes)
＊ 函数参数：
＊ LPSTR lpSrcStartBits 指向DIB起始像素的指针
＊ long lWidth DIB的宽度
＊ long lHeight DIB的高度
＊ long lLineBytes DIB的行字节数，为4的倍数
＊ 函数功能： 用来对图像进行离散余弦变换
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊/
BOOL CosTran(LPSTR lpSrcStartBits, long lWidth, long lHeight,long lLineBytes)
{
unsigned char＊ lpSrcUnChr; // 指向像素的指针
long i; // 行循环变量
long j; // 列循环变量
long wid=1,hei=1; // 进行傅里叶变换的宽度和高度，初始化为1
double dTemp; // 中间变量
int widpor=0,heiPor=0; // 2的幂数
while(wid ＊ 2 <= lWidth) // 计算进行离散余弦变换的宽度和高度（2的整数次幂）
{
wid ＊= 2;
widpor++;
}
while(hei ＊ 2 <= lHeight)
{
hei ＊= 2;
heiPor++;
}
double ＊pTd= new double[wid ＊ hei]; // 分配内存
double ＊pFd = new double[wid ＊ hei];
for(i = 0; i < hei; i++) // 行
{
for(j = 0; j < wid; j++) // 列
{
// 指向DIB第i行、第j个像素的指针
lpSrcUnChr= (unsigned char＊)lpSrcStartBits + lLineBytes ＊ (lHeight-1- i) + j;
pTd[j + i ＊ wid] = ＊(lpSrcUnChr); // 给时域赋值
}
}
for(i = 0; i < hei; i++)
{
DisFCosTran(&pTd[wid ＊ i], &pFd[wid ＊ i], widpor); // 对y方向进行离散余弦变换
}
for(i = 0; i < hei; i++) // 保存计算结果
for(j = 0; j < wid; j++)
pTd[j ＊ hei + i] = pFd[j + wid ＊ i];
for(j = 0; j < wid; j++)
{
DisFCosTran(&pTd[j ＊ hei], &pFd[j ＊ hei], heiPor); // 对x方向进行离散余弦变换
}
for(i = 0; i < hei; i++) // 行
{
for(j = 0; j < wid; j++) // 列
{
dTemp = fabs(pFd[j＊hei+i]); // 计算频谱
if (dTemp > 255) // 是否超过255，超过255的直接赋值为255
{
dTemp = 255;
}
// 指向DIB第y行、第x个像素的指针
lpSrcUnChr= (unsigned char＊)lpSrcStartBits + lLineBytes ＊ (lHeight-1- i) + j;
＊ (lpSrcUnChr) = (BYTE)(dTemp); // 更新源图像
}
}
delete pTd; // 释放内存
delete pFd;
return TRUE;
}
编写VC的函数DisFCosTran( )来实现图像的快速离散余弦变换，其具体代码如下：
/＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊ 函数名称：DisFCosTran(double ＊pTd, double ＊pFd, int power)
＊ 函数参数：
＊ double ＊ pTd 指向时域值的指针
＊ double ＊ pFd 指向频域值的指针
＊ int power 2的幂数
＊ 函数功能：用来实现快速离散余弦变换
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊/
void DisFCosTran(double ＊pTd, double ＊pFd, int power)
{
long i; // 行循环变量
long dotCount; // 离散余弦变换点数
double dTemp; // 临时变量
CplexNum ＊temReg;
dotCount = 1<<power; // 计算离散余弦变换点数
temReg = new CplexNum[sizeof(CplexNum) ＊dotCount＊2]; // 分配内存
memset(temReg, 0, sizeof(CplexNum) ＊ dotCount ＊ 2); // 赋为0
for(i=0;i<dotCount;i++) // 将时域点写入数组temReg
{
temReg[i].re=pTd[i];
temReg[i].im=0;
}
FastFourierTran(temReg,temReg,power+1); // 调用快速傅里叶变换
dTemp = 1/sqrt((double)dotCount); // 调整系数
pFd[0] = temReg[0].re＊dTemp; // 求pFd[0]
dTemp ＊= sqrt(2.0f);
for(i = 1; i < dotCount; i++) // 求pFd[u]
{
pFd[i]=(temReg[i].re＊ cos(i＊pi/(dotCount＊2)) + temReg[i].im＊ sin(i＊pi/(dotCount＊2)))
＊dTemp;
}
delete temReg; // 释放内存
}
可以通过在数字图像处理程序的框架窗口上增加菜单【正交变换】|【离散余弦变换】。对应视图类中的CDImageProcessView:: OnDct ( )进行函数的调用：
void CDImageProcessView::OnDct()
{
CDImageProcessDoc＊ pDoc = GetDocument();
SetCursor(LoadCursor(NULL, IDC_WAIT)); // 光标显示为处理状态
long lSrcLineBytes; // 图像每行的字节数
long lSrcWidth; // 图像的宽度和高度
long lSrcHeight;
LPSTR lpSrcDib; // 指向源图像的指针
LPSTR lpSrcStartBits; // 指向源像素的指针
lpSrcDib= (LPSTR) ::GlobalLock((HGLOBAL) pDoc->GetHObject()); // 锁定DIB
if (pDoc->m_dib.GetColorNum(lpSrcDib) != 256) // 判断是否是8-bpp位图
{
AfxMessageBox(_T ("对不起，不是256色位图！")); // 警告
::GlobalUnlock((HGLOBAL) pDoc->GetHObject()); // 解除锁定
return; // 返回
}
lpSrcStartBits=pDoc->m_dib.GetBits(lpSrcDib); // 找到DIB图像像素起始位置
lSrcWidth= pDoc->m_dib.GetWidth(lpSrcDib); // 获取图像的宽度
lSrcHeight= pDoc->m_dib.GetHeight(lpSrcDib); // 获取图像的高度
lSrcLineBytes=pDoc->m_dib.GetReqByteWidth(lSrcWidth ＊ 8);// 计算图像每行的字节数
// 调用Dct()函数进行离散余弦变换
if (CosTran(lpSrcStartBits, lSrcWidth, lSrcHeight,lSrcLineBytes))
{
pDoc->SetModifiedFlag(TRUE); // 设置修改标记
pDoc->UpdateAllViews(NULL); // 更新视图
::GlobalUnlock((HGLOBAL) pDoc->GetHObject()); // 解除锁定
}
else
{
AfxMessageBox(_T("分配内存失败！")); // 警告
}
SetCursor(LoadCursor(NULL, IDC_ARROW)); // 光标显示为处理状态
}
4. 效果演示
源图像如图3-4a所示，对其进行离散余弦变换，结果如图3-4b所示。
图3-4 离散余弦变换效果演示
3.1.3 离散沃尔什变换
离散沃尔什变换是由两个数值（+1和-1）为基本函数的级数展开而形成的，它也满足完备的正交特性。由于沃尔什函数是二值正交函数，与数字逻辑中的两个状态相对应，因此更适用于计算机技术、数字信号和数字图像处理。
1. 基本原理
（1）沃尔什函数
沃尔什函数是由美国数学家沃尔什在1923年提出的完整的数字理论。按照排列次序，沃尔什函数可以分为3种定义方法：第1种是按照沃尔什排列来定义（按列率排序）；第2种是按照佩利排列来定义（按自然排序）；第3种是按照哈达玛排列来定义。由于哈达玛排序的沃尔什函数是由2n (n= 0, 1, 2, ⋯)阶哈达玛矩阵得到的，而哈达玛矩阵的最大优点在于其具有简单的递推关系，即高阶矩阵可用两个低阶矩阵的克罗内克积求得，因此在此只介绍哈达玛排列定义的沃尔什变换，即离散沃尔什-哈达玛变换。
哈达玛序的离散沃尔什函数定义为
其中，⊕表示模2加法；tp-1-r为t的二进制形式的第p-1-r位；nr为n的二进制形式的第r位。
例如，在取p=3，计算n = 4时的沃尔什函数的数值上，可以把n写成二进制码的形式，则n=4=n2n1n0=100，t=t2t1t0，从而有
因此，可以得到p = 3时的哈达玛序的离散沃尔什函数的矩阵形式为8×8矩阵(N= 2p= 23= 8)，记为
若设最低阶(N= 2)的哈达玛序的离散沃尔什函数矩阵为
再令HN代表N阶(N= 2)哈达玛序的离散沃尔什函数矩阵，则可以得到以下递推关系式
其中，
为克罗内克积。由此，可以产生任意的2n阶哈达玛矩阵。
（2）一维离散沃尔什哈达玛变换
一维离散沃尔什哈达玛变换定义为
其中，WalH为沃尔什哈达玛函数，u, x= 0, 1, 2, ⋯, N－1。
可用矩阵表示为
其逆变换为
可用矩阵表示为
由哈达玛矩阵的特点可知，离散沃尔什哈达玛变换的本质是将离散序列f (x)的各项值的符号按一定规律改变后，进行加减运算。因此，它比采用复数运算DFT和采用余弦运算的DCT要简单得多。
（3）二维离散沃尔什哈达玛变换
设图像函数f (x, y)的大小为M×N，则哈达玛序的二维离散沃尔什哈达玛变换定义为
其中，u= 0, 1, 2, ⋯, M－1；v= 0, 1, 2, ..., N－1。f (x, y)代表图像像素，x、y分别该像素在空间中的横、纵坐标。
其逆变换为
其中，x= 0, 1, 2, ⋯, M－1；y= 0, 1, 2, ⋯, N－1。
哈达玛变换具有直线性、相乘性、位移定理及平均值定理等性质，与离散傅里叶变换相似，二维离散沃尔什哈达玛变换也是一维情况的推广，也就是二维离散沃尔什哈达玛变换也可以用两次一维离散沃尔什哈达玛变换来计算，且在运算过程中也可以像快速傅里叶变换一样，用蝶形图实现快速离散沃尔什哈达玛变换。即以N = Nx对f (x, y)中Ny个列中的每一列进行一维离散沃尔什哈达玛变换，得出Wx (u, y)，再以N = Ny对Wx (u, y)中Nx个行中的每一行进行一维离散沃尔什哈达玛变换，最终得到Wxy (u, v)。
2. 算法描述
实现数字图像在屏幕上的离散沃尔什哈达玛变换，可以分为以下5个步骤：
[1] 将源图像保存到缓冲区，并记录下缓冲区的地址。
[2] 获取源图像的高度和宽度，计算图像每行的字节数。
[3] 调用沃尔什哈达玛变换函数，对图像进行离散沃尔什哈达玛变换。
[4] 在离散沃尔什哈达玛变换过程中，调用快速沃尔什哈达玛变换函数。
[5] 更新显示变换后的图像。
3. 编程实现
VC函数Walsh_HarTran ( )实现图像的离散沃尔什哈达玛变换，其具体代码如下：
/＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊ 函数名称：Walsh_HarTran(LPSTR lpSrcStartBits, long lWidth, long lHeight,long lLineBytes)
＊ 函数参数：
＊ LPSTR lpSrcStartBits, 指向源DIB图像指针
＊ long lWidth, 源DIB图像宽度
＊ long lHeight, 源DIB图像高度
＊ long lLineBytes, 源DIB图像的行字节数，为4的倍数
＊ 函数功能：用来对图像进行沃尔什哈达玛变换
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊/
BOOL Walsh_HarTran(LPSTR lpSrcStartBits, long lWidth, long lHeight,long lLineBytes)
{
unsigned char＊ lpSrcUnChr; // 指向像素的指针
long i; // 行循环变量
long j; // 列循环变量
long wid=1,hei=1; // 进行傅里叶变换的宽度和高度，初始化为1
double dTemp; // 中间变量
int widpor=0,heiPor=0; // 2的幂数
while(wid ＊ 2 <= lWidth)// 计算进行离散沃尔什哈达玛变换的宽度和高度（2的整数次幂）
{
wid ＊= 2;
widpor++;
}
while(hei ＊ 2 <= lHeight)
{
hei ＊= 2;
heiPor++;
}
double ＊pTd = new double[wid ＊ hei]; // 分配内存
double ＊pFd = new double[wid ＊ hei];
for(i = 0; i < hei; i++) // 行
{
for(j = 0; j < wid; j++) // 列
{
// 指向DIB第i行、第j个像素的指针
lpSrcUnChr= (unsigned char＊)lpSrcStartBits + lLineBytes ＊ (lHeight-1- i) + j;
pTd[j + i ＊ wid] = ＊(lpSrcUnChr); // 给时域赋值
}
}
for(i = 0; i < hei; i++)
{
Walshei_Har(pTd + wid ＊ i, pFd + wid ＊ i, widpor); // 对y方向进行沃尔什哈达玛变换
}
for(i = 0; i < hei; i++) // 保存计算结果
{
for(j = 0; j < wid; j++)
{
pTd[j ＊ hei + i] = pFd[j + wid ＊ i];
}
}
for(j = 0; j < wid; j++)
{
Walshei_Har(pTd + j ＊ hei, pFd+ j ＊ hei, heiPor);// 对x方向进行沃尔什哈达玛变换
}
for(i = 0; i < hei; i++) // 行
{
for(j = 0; j < wid; j++) // 列
{
dTemp = fabs(pFd[j ＊ hei + i] ＊ 1000); // 计算频谱
if (dTemp > 255) // 对于超过255的，直接设置为255
dTemp = 255;
lpSrcUnChr= (unsigned char＊)lpSrcStartBits + lLineBytes ＊ (lHeight-1- i) + j;
＊ (lpSrcUnChr) = (BYTE)(dTemp); // 更新源图像
}
}
delete pTd; // 释放内存
delete pFd;
return TRUE;
}
编写VC的函数Walshei_Har ( )来实现图像的快速离散沃尔什哈达玛变换，代码如下：
/＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊ 函数名称：Walshei_Har(double ＊pTd, double ＊pFd, int power)
＊ 函数参数：
＊ double ＊ pTd 指向时域值的指针
＊ double ＊ pFd 指向频域值的指针
＊ int power 2的幂数
＊ 函数功能： 用来实现快速沃尔什哈达玛变换
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊/
void Walshei_Har(double ＊pTd, double ＊pFd, int power)
{
long dotCount; // 沃尔什哈达玛变换点数
int i,j,k; // 循环变量
int bfsize,p; // 中间变量
double ＊temReg1,＊temReg2,＊temReg;
dotCount = 1 << power; // 计算快速沃尔什变换点数
temReg1 = new double[dotCount]; // 分配运算所需的数组
temReg2 = new double[dotCount];
memcpy(temReg1, pTd, sizeof(double) ＊ dotCount); // 将时域点写入数组temReg1
for(k = 0; k < power; k++) // 蝶形运算
{
for(j = 0; j < 1<<k; j++)
{
bfsize = 1 << (power-k);
for(i = 0; i < bfsize / 2; i++)
{
p = j ＊ bfsize;
temReg2[i + p] = temReg1[i + p] + temReg1[i + p + bfsize / 2];
temReg2[i + p + bfsize / 2] = temReg1[i + p] - temReg1[i + p + bfsize / 2];
}
}
temReg = temReg1; // 互换temReg1和temReg2
temReg1 = temReg2;
temReg2 = temReg;
}
for(j = 0; j < dotCount; j++) // 调整系数
{
p = 0;
for(i = 0; i < power; i++)
{
if (j & (1<<i))
{
p += 1 << (power-i-1);
}
}
pFd[j] = temReg1[p] / dotCount;
}
delete temReg1; // 释放内存
delete temReg2;
}
可以通过在数字图像处理程序的框架窗口上增加菜单【正交变换】|【沃尔什哈达玛变换】。对应视图类中的CDImageProcessView:: OnWalhHar ( )进行函数的调用：
void CDImageProcessView::OnWalhHar()
{
CDImageProcessDoc＊ pDoc = GetDocument();
long lSrcLineBytes; // 图像每行的字节数
long lSrcWidth; // 图像的宽度和高度
long lSrcHeight;
LPSTR lpSrcDib; // 指向源图像的指针
LPSTR lpSrcStartBits; // 指向源像素的指针
lpSrcDib= (LPSTR) ::GlobalLock((HGLOBAL) pDoc->GetHObject()); // 锁定DIB
if (pDoc->m_dib.GetColorNum(lpSrcDib) != 256)
{
AfxMessageBox(_T ("对不起，不是256色位图！")); // 警告
::GlobalUnlock((HGLOBAL) pDoc->GetHObject()); // 解除锁定
return; // 返回
}
lpSrcStartBits=pDoc->m_dib.GetBits(lpSrcDib);// 找到DIB图像像素起始位置
lSrcWidth= pDoc->m_dib.GetWidth(lpSrcDib); // 获取图像的宽度
lSrcHeight= pDoc->m_dib.GetHeight(lpSrcDib); // 获取图像的高度
lSrcLineBytes=pDoc->m_dib.GetReqByteWidth(lSrcWidth ＊ 8);// 计算图像每行的字节数
if (Walsh_HarTran(lpSrcStartBits, lSrcWidth, lSrcHeight,lSrcLineBytes))
{
pDoc->SetModifiedFlag(TRUE); // 设置修改标记
pDoc->UpdateAllViews(NULL); // 更新视图
::GlobalUnlock((HGLOBAL) pDoc->GetHObject()); // 解除锁定
}
else
{
AfxMessageBox(_T("分配内存失败！")); // 警告
}
}
4. 效果演示
源图像如图3-5a所示，其离散沃尔什哈达玛变换结果如图3-5b所示。
图3-5 离散沃尔什哈达玛变换效果演示
3.2 特征变换
数字图像的特征变换包括K-L变换、SVD变换、小波变换等，主要用于图像的特征提取，是图像分析的基础，在图像压缩编码中也会用到。下面对这3种变换进行详细介绍。
3.2.1 K-L变换
当变量之间存在一定的相关关系时，可以通过原始变量的线性组合，构成较少的不相关的新变量代替原始变量，而每个新变量都含有尽量多的原始变量的信息。这种处理问题的方法叫主成分分析（PCA），新变量叫原始变量的主成分。主成分分析的目的是寻找任意统计分布的数据集合主要分量的子集。相应的基向量组满足正交性且由其定义的子空间最优地考虑数据的相关性，从而将原始数据集合变换到主分量空间，使单一数据样本的互相关性降低到最低点。
离散卡胡南-洛夫（Karhunen-Loeve）变换简称为K-L变换。K-L变换基于主成分分析理论，以图像的统计特性为基础的最优正交线性变换。其变换核矩阵由图像阵列的协方差矩阵的特征值和特征向量决定，所以K-L变换又称为特征向量变换。
1. 基本原理
（1）离散K-L变换的概念
假设将某幅N×N的图像f (x, y)在某个传输通道上传输了M次，因会受到各种因素的随机干扰，接收到的图像集合为：
{f1(x, y), f2 (x, y), ⋯, fi(x, y), ⋯, fM (x, y)}
其中，x, y= 0, 1, 2, ⋯, N－1。
将M次传送的图像集合写成M个N2维向量{X1, X2, ⋯, Xi, XM}，生成向量的方法可以采用行堆叠或列堆叠的方法，对第i次获得的图像fi (x, y)，可用N2维向量Xi 表示为
若存在一个合适的正交变换A，使得变换后的图像Y = AX满足如下两个条件：
▓ 图像Y具有M个分量，且M<<N2。
▓ 图像Y经反变换A-1Y恢复出的向量X＇和原始向量X具有最小的均方误差。
此时称正交变换A为K-L变换。也就是说经过一个变换，不仅删除了N2-M个分量，但仍能很好地恢复出源图像（向量）且由变换结果Y重新恢复的图像。
对应式（3-12），X向量的协方差矩阵定义为CX = E{(X-mX) (X-mX)T}
其中E为期望值；T为矩阵的转置；mX为平均值向量。
mX = E{X}
在离散情况下，由于M是有限值，其平均值向量mX和X向量的协方差矩阵CX可近似表示为
设ei和λi是协方差矩阵CX对应的特征向量和特征值，其中，i= 0, 1, 2, ⋯, N2。将特征值按照减序排列，即
λ1>λ2>λ3>⋯>λN2
则K-L变换核矩阵A的行用CX的特征值构成
其中，eij表示第i特征向量的第j个分量；A是N 2× N 2方阵。
由此可得离散K-L变换式为
其中，X－mX是原始图像向量X减去平均值向量mX，称为中心化图像向量。离散K-L变换向量Y是中心化向量X－mX与变换核矩阵相乘所得的结果。
（2）离散K-L变换的性质
▓ Y的平均值向量mY= 0，即为零向量0。
mY= E{Y}= E{A(X-mX)}= A E{X}-AmX = 0
▓Y向量的协方差矩阵为CY= ACX AT。
由CY=E{(Y-mY)(Y-mY)T}=E{Y Y T}，代入离散K-L变换式（3-13）中，可得
▓ 对角性。因为
可知，Y向量的协方差矩阵CX是对角矩阵，对角线上的元素是原始图像向量的协方差矩阵CX对应的特征值λi，也是Y向量的方差。如果非对角线上的元素值（协方差）为0，说明Y向量中各元素之间相关性很小；如果CX矩阵的非对角线上元素不为0，则说明原始图像元素之间相关性强。这就是采用K-L变换进行编码，数据压缩比大的原因。
显然K-L坐标系将矩阵CX对角化了，换句话说，通过K-L变换，消除了原有向量X的各分量之间的相关性，从而可能去掉那些带有较少信息的坐标轴，以达到降低特征空间维数的目的，所以K-L变换的最大优点是去除相关性佳，可用于特征提取、数据压缩和图像旋转等。K-L变换的不足之处在于运用协方差矩阵CX，求特征值λ和特征向量的计算量大，同时K-L变换是非分离的，二维不可分。一般情况下，K-L变换没有快速算法。
2. 算法描述
以图像旋转为例，实现用离散K-L变换完成数字图像在屏幕上的旋转，可以分为以下4个步骤：
[1] 将源图像保存到缓冲区，并记录下缓冲区的地址。
[2] 获取源图像的高度和宽度，计算图像旋转后的最大可能范围。
[3] 将源图像像素的坐标值看作二维随机矢量，调用离散K-L变换函数，对图像进行离散K-L变换。
[4] 更新显示变换后的图像，不显示已经移出源图区域的图像。
3. 编程实现
下面通过编写VC的函数DisK_L ( )来实现图像的离散K-L变换，其具体代码如下：
/＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊ 函数名称：DisK_L(LPSTR lpSrcStartBits, long lWidth, long lHeight,long lLineBytes)
＊函数参数：
＊ LPSTR lpSrcStartBits 指向源DIB图像指针
＊ long lWidth 源DIB图像宽度
＊ long lHeight 源DIB图像高度
＊ long lLineBytes 源DIB图像的行字节数，为4的倍数
＊ 函数功能：用离散K-L变换对图像进行旋转
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊/
BOOL DisK_L(LPSTR lpSrcStartBits, long lWidth, long lHeight,long lLineBytes)
{
unsigned char＊ lpSrcUnChr; // 指向像素的指针
LONG i,j, // 循环变量
lMaxRange, // 经过变换后，图像的最大可能范围
AverEx,AverEy, // 目标坐标均值
ToaCount; // 目标总的像素数
double Matr4C[2][2], // 坐标值的协方差矩阵
QMatrix[2][2], // 存放协方差矩阵的特征向量
MainCross[2],HypoCross[2], // 三对角阵的主对角和次对角线元素
dTemp; // 临时变量
LONG lTempI,lTempJ;
if(lWidth>lHeight)// 估计图像经过旋转后可能最大的宽度和高度
lMaxRange = lWidth;
else
lMaxRange =lHeight;
AverEx=0.0; // 初始化
AverEy=0.0;
ToaCount = 0;
Matr4C[0][0] = Matr4C[0][1] = Matr4C[1][0] = Matr4C[1][1] = 0.0;
double ＊F = new double[lWidth＊lHeight]; // 分配内存
for(i = 0; i < lHeight; i++) // 行
{
for(j = 0; j < lWidth; j++) // 列
{
F[i＊lWidth + j] = 255; // 给旋转后坐标轴的每个点赋零值
// 指向位图第i行、第j列像素的指针
lpSrcUnChr= (unsigned char＊)lpSrcStartBits + lLineBytes＊i + j;
// 值小于255（非背景色白色）的像素认为目标的一部分，并将其坐标值x和y
// 看作二维随机矢量
if((＊lpSrcUnChr) < 255)
{
AverEx=AverEx+i; // 属于目标像素的y坐标和x坐标累计值
AverEy=AverEy+j;
ToaCount++; // 目标总的像素数加1
// 随机矢量协方差矩阵的累计值
Matr4C[0][0] = Matr4C[0][0] + i＊i;
Matr4C[0][1] = Matr4C[0][1] + i＊j;
Matr4C[1][0] = Matr4C[1][0] + j＊i;
Matr4C[1][1] = Matr4C[1][1] + j＊j;
}
}
}
AverEx = AverEx/ToaCount; // 计算随机矢量的均值
AverEy = AverEy/ToaCount;
Matr4C[0][0] = Matr4C[0][0]/ToaCount - AverEx＊AverEx; // 计算随机矢量的协方差矩阵
Matr4C[0][1] = Matr4C[0][1]/ToaCount - AverEx＊AverEy;
Matr4C[1][0] = Matr4C[1][0]/ToaCount - AverEx＊AverEy;
Matr4C[1][1] = Matr4C[1][1]/ToaCount - AverEy＊AverEy;
double Precision = 0.000001; // 规定迭代的计算精度
ThreeCrossMat(＊Matr4C, 2, ＊QMatrix, MainCross, HypoCross);
// 将协方差矩阵化作三对角阵
EigenvalueVector(2, MainCross,HypoCross, ＊Matr4C, Precision, 50);
// 求协方差矩阵的特征值和特征矢向量
dTemp = Matr4C[0][1]; // 将特征列向量转化成特征列向量
Matr4C[0][1] = Matr4C[1][0];
Matr4C[1][0] = dTemp;
for(i=0;i<=1;i++)
{
// 对特征列向量进行归一化
dTemp = pow(Matr4C[i][0],2) + pow(Matr4C[i][1],2);
dTemp = sqrt(dTemp);
Matr4C[i][0] = Matr4C[i][0]/dTemp;
Matr4C[i][1] = Matr4C[i][1]/dTemp;
}
// 查找经离散K-L变换后的坐标点在原坐标系中的坐标
for(i = -lMaxRange+1; i < lMaxRange; i++)
{
for(j = -lMaxRange+1; j < lMaxRange; j++)
{
// 将新坐标值映射到原坐标系
int Cx = (int)(i＊Matr4C[0][0]-j＊Matr4C[0][1])+AverEx;
int Cy = (int)(-i＊Matr4C[1][0]+j＊Matr4C[1][1])+AverEy;
// 映射值是否属于源图像
if( Cx>=0 && Cx<lHeight && Cy>=0 && Cy<lWidth )
{
lpSrcUnChr= (unsigned char＊)lpSrcStartBits + lLineBytes＊Cx + Cy;
// 映射值是否属于原来的目标
if(＊(lpSrcUnChr)<255)
{
// 将新坐标系原点平移到中心，以便显示
lTempI=(LONG)(lHeight/2)+j;
lTempJ=(LONG)(lWidth/2)+i;
// 看如果能够进行显示，赋值给数组，进行存储
if( lTempI>=0 && lTempI<lHeight && lTempJ>=0 && lTempJ<lWidth )
F[lTempJ+ (lTempI) ＊ lWidth]=＊(lpSrcUnChr);
}
}
}
}
for(i = 0; i < lMaxRange; i++) // 行
{
for(j = 0; j < lMaxRange; j++) // 列
{
dTemp = F[i ＊ lMaxRange + j] ; // 离散K-L变换后的像素值
// 指向位图第i行、第j列像素的指针
lpSrcUnChr= (unsigned char＊)lpSrcStartBits + lLineBytes ＊ (lHeight -1- i) + j;
＊ (lpSrcUnChr) = (BYTE)(dTemp); // 更新源图像
}
}
delete F; // 释放内存
return TRUE; // 返回
}
可以通过在数字图像处理程序的框架窗口上增加菜单【正交变换】|【离散K-L变换】。对应视图类中的CDImageProcessView:: OnKL ( )进行函数的调用：
void CDImageProcessView::OnKL()
{
CDImageProcessDoc＊ pDoc = GetDocument(); // 获取文档
long lSrcLineBytes; // 图像每行的字节数
long lSrcWidth; // 图像的宽度和高度
long lSrcHeight;
LPSTR lpSrcDib; // 指向源图像的指针
LPSTR lpSrcStartBits; // 指向源像素的指针
lpSrcDib= (LPSTR) ::GlobalLock((HGLOBAL) pDoc->GetHObject()); // 锁定DIB
if (pDoc->m_dib.GetColorNum(lpSrcDib) != 256)
{
AfxMessageBox(_T ("对不起，不是256色位图！")); // 警告
::GlobalUnlock((HGLOBAL) pDoc->GetHObject()); // 解除锁定
return; // 返回
}
lpSrcStartBits=pDoc->m_dib.GetBits(lpSrcDib);// 找到DIB图像像素起始位置
lSrcWidth= pDoc->m_dib.GetWidth(lpSrcDib); // 获取图像的宽度
lSrcHeight= pDoc->m_dib.GetHeight(lpSrcDib); // 获取图像的高度
lSrcLineBytes=pDoc->m_dib.GetReqByteWidth(lSrcWidth ＊ 8);// 计算图像每行的字节数
if (DisK_L(lpSrcStartBits,lSrcWidth,lSrcHeight,lSrcLineBytes))
{
pDoc->SetModifiedFlag(TRUE); // 设置修改标记
pDoc->UpdateAllViews(NULL); // 更新视图
::GlobalUnlock((HGLOBAL) pDoc->GetHObject()); // 解除锁定
}
else
{
AfxMessageBox(_T("分配内存失败！")); // 警告
}
}
4. 效果演示
源图像如图3-6a所示，对其进行离散K-L变换后的效果如图3-6b所示。
图3-6 离散K-L变换效果演示
3.2.2 SVD变换
SVD变换即奇异值分解变换，是一种对非对称矩阵的正交变换，变换的结果矩阵是有一子块为对角矩阵，而其余阵元均为零。
1. 基本原理
从矩阵论的角度来看，奇异值分解可定义如下：若A是m×n的实矩阵，且A的秩为r(r∈[0, min(m, n)])，则存在m阶正交矩阵U和n阶正交矩阵V，使得
其中，∑是非负对角矩阵，且∑=diag(σ1,σ2,⋯,σr)，σ1≥σ2≥⋯≥σr>0是A的全部非零奇异值。
矩阵ATA和AAT具有相同的特征值λ，且λ的平方根就是A的奇异值，即
。
若A为数字图像，则一个M×N实值图像A可以看作是M×N维空间中的一个向量，也可将此图像表示在r维子空间，为r矩阵A的秩。由奇异值分解的定义可得
式（3-14）和（3-15）所表示的正反变换，叫做奇异值分解（SVD）变换。其中，σi是A的非零奇异值，U和V分别是M×r和N×r矩阵，它们的列向量分别是ATA和AAT的正交特征向量，
可以看作是一个基图像，则数字图像A表示为r个秩为1的基图像的加权和，权系数是奇异值σi。
U和V的核矩阵取决于图像A，进行图像的SVD变换一般需要计算图像A的ATA和AAT特征向量。图像按奇异值分解后，有如下性质：
▓图像的奇异值代表图像的能量信息，因而具有稳定性。
▓图像的奇异值具有比例不变性。
▓图像的奇异值具有旋转不变性。
因此，图像的奇异值变换在图像特征提取、压缩、检索、匹配、数字水印等领域应用广泛。以图像压缩为例，数字图像A经奇异值分解后，其纹理和几何信息都集中在U、V之中，而∑中的奇异值则代表它刻画了图像的能量信息，若N≤M，则对角阵Λ至多有N个非零元素，至少可以获得M倍的无损压缩。同时，奇异值中总有一些小到可以被忽略的值，忽略这些值依然可得到满意的视觉效果，选择合适的奇异值个数来实现图像的有损压缩，可以得到更高的压缩比。奇异值个数选取的越少，图像的压缩比率就越大。调节压缩率可以达到不同的压缩效果，是一种实现过程简单、行之有效的压缩方法，具有一定的应用前景。
2. 算法描述
以图像压缩为例，实现用SVD变换完成数字图像在屏幕上的压缩，达到既大幅度减小图像的容量，又基本满足视觉效果的目的。可以分为以下6个步骤：
[1] 将源图像保存到缓冲区，并记录下缓冲区的地址。
[2] 获取源图像的高度和宽度，创建图像（M×N）的数据矩阵A。
[3] 根据式（3-14）对图像矩阵A进行SVD变换，得到矩阵U、V和Λ。
[4] 使用对话框设置图像的压缩比率，根据压缩比率选择合适的奇异值个数K，舍弃多余的奇异值，得到K×K的矩阵Λ，M×K的矩阵UK和N×K矩阵VK。
[5] 由公式
，得到压缩后的矩阵AK (M×N)。
[6] 将变换后的矩阵数据赋值给图像，更新显示变换后的图像。
3. 编程实现
由于数字图像的SVD变换中使用了大量的矩阵运算,编程和调试工作比较复杂且极易出错，所以在VC中使用Matrix<LIB>C++库，来完成图像变换的操作。
Matrix<LIB>C++数学库是MathTools公司利用MatCom技术开发的一个面向专业从事工程技术和科学计算人员的矩阵运算动态链接库。该库提供了绝大多数的关于矩阵类、矩阵操作及数值计算等函数的定义。
具体实现步骤如下：
[1] 将v4501v.lib和matlib.h复制到所需用到的工程文件路径中。
v4501v.lib和matlib.h这两个文件可以通过安装MATCOM软件（其最新版本为4.5）后，在其安装路径下的lib文件夹中得到。
[2] 单击【项目】|【属性】项，显示【工程属性页】对话框，如图3-7所示。在【工程属性页】对话框中的【配置属性】【链接器】|【命令行】右侧的“附加选项（D）”中，添加v4501v.lib文件。
图3-7 【工程属性页】对话框
[3] 在需要使用Matrix<LIB>C++数学库的头文件（.h）和实现文件（.cpp）中添加头文件，即包含如下代码：
#include "matlib.h"
[4] 初如化Matrix<LIB>C++数学库，在MainFrm.cpp文件中加入如下代码：
CMainFrame::CMainFrame()
{
initM(MATCOM_VERSION); //初始化类库调用
}
CMainFrame::～CMainFrame()
{
exitM(); //结束类库调用
}
[5] 将Matrix<LIB>C++数学库的动态链接库文件ago4501.dll、v4501v.dll、opengl132.dll和glu32.dll复制到Windows的system32文件夹下。
如果系统已经安装了MATCOM 4.5，则这4个动态链接库文件会自动安装在相应的目录下。
[6] 通过在数字图像处理程序的框架窗口上增加菜单【正交变换】|【SVD变换】。
[7] 编写对应视图类中的CDImageProcessView:: OnSvd ( )进行图像变换的代码如下：
void CDImageProcessView::OnSvd()
{
CDImageProcessDoc＊ pDoc = GetDocument(); // 获取文档
LPSTR lpDIB; // 指向源图像的指针
LPSTR lDstLineBytes; // 新图像每行的字节数
LPSTR lpSrcStartBits; // 指向DIB像素的指针
Mm m_matBits;
lpDIB = (LPSTR) ::GlobalLock((HGLOBAL) pDoc->GetHObject()); // 锁定DIB
lpSrcStartBits = pDoc->m_dib.GetBits(lpDIB);// 找到DIB图像像素起始位置
int nWidth; // 图像宽度
int nHeight; // 图像高度
long nWidthBytes;
nWidth=(int)pDoc->m_dib.GetWidth(lpDIB); // 获取图像宽度
nHeight=(int)pDoc->m_dib.GetHeight(lpDIB); // 获取图像高度
nWidthBytes=pDoc->m_dib.GetReqByteWidth(nWidth ＊ 8);// 计算图像每行的字节数
DWORD SizeImage = nWidthBytes ＊ nHeight; // 计算图像数据的总数
m_matBits = zeros(1,SizeImage); //创建一个1×SizeImage的一维0矩阵
m_matBits = muint8(m_matBits);//将矩阵由整型转换为double型，以便将图像数据赋给矩阵
memcpy(m_matBits.addr(),lpSrcStartBits,SizeImage);
m_matBits = rot90(reshape(m_matBits,nWidthBytes,nHeight));
if (nWidthBytes != nWidth)
{
m_matBits = m_matBits(c_p,colon(1,1,nWidth));
}
m_matBits = mdouble(m_matBits); // 将矩阵由unsigned char型转换为double型
::GlobalUnlock((HGLOBAL) pDoc->GetHObject()); // 内存解锁
Mm mSize ;
Mm Usvd; // 使用Mm定义矩阵
Mm Wsvd;
Mm Ssvd;
Mm ff;
Mm uu;
Mm ww;
mSize = size(m_matBits); // 获得矩阵的行数和列数
Usvd = svd_U(m_matBits); // 矩阵SVD变换
Wsvd = svd_W(m_matBits);
Ssvd = svd_S(m_matBits);
int fsvdAngle;
fsvdAngle = columns(m_matBits); // 得到矩阵列长度
CDlgsvd svdPara; // 创建对话框
if (svdPara.DoModal() != IDOK) // 显示对话框，提示用户设定量
return;
int temsvd = svdPara.m_svd;
int k; // 奇异值的个数
k = fsvdAngle ＊ temsvd /100;
Mm ffn;
ff = Ssvd(colon(1.0,1.0, k),c_p);
ff = transpose(ff);
ffn = ff(colon(1.0,1.0, k),c_p);
ffn = transpose(ffn);
Usvd = transpose(Usvd);
uu = Usvd(colon(1.0,1.0, k),c_p);
uu = transpose(uu);
Wsvd = transpose(Wsvd);
ww = Wsvd(colon(1.0,1.0, k),c_p);
m_matBits = uu ＊ ffn ＊ ww; // 变换后的压缩矩阵
m_matBits = muint8(m_matBits); // 将矩阵赋值给图像
int nTmp = (int)rem(nWidth,4);
int nPadding;
if (nTmp > 0)
{
nPadding = 4- nTmp;
m_matBits = cat(2,m_matBits,repmat(muint8(0),(BR(size(m_matBits)),nPadding)));
}
m_matBits = rot90(m_matBits,-1);
memcpy(lpSrcStartBits,m_matBits.addr(),(nWidthBytes＊nHeight)＊sizeof(unsigned char));
pDoc->SetDib(); // 更新DIB大小和调色板
pDoc->SetModifiedFlag(TRUE); // 设置修改标记
pDoc->UpdateAllViews(NULL); // 更新视图
::GlobalUnlock((HGLOBAL) pDoc->GetHObject());
}
4. 效果演示
源图像如图3-8a所示，对其进行SVD变换，分别取压缩比率为50%、30%和20%，其结果分别如图3-8b、图3-8c和图3-8d所示。可以看到，当压缩比率大于30%时，均可得到比较满意的视觉效果。
图3-8 SVD变换效果演示
3.2.3 小波变换
由于傅里叶分析是一种全局的变换，要么完全在时间域，要么完全在频率域，它以两个方向上都无限伸展的正弦曲线波作为正交基函数，因此无法表述信号的时频局部性质。而时频局部性质恰好是非平稳信号最基本和最关键的性质。为了分析和处理非平稳信号，能够对瞬态信号或高度局部化的信号（如边缘）进行分析，在傅里叶分析理论基础上，定义一组能够表现出信号瞬时性的基函数来帮助计算信号的瞬时傅里叶变换。这组有限宽度的基函数不仅在频率上是变化的，而且在位置上也是变化的，我们把有限宽度的波称为小波，基于这种波的变换就称为小波变换。
小波变换是在克服傅里叶变换缺点的基础上发展而来的强有力的时频分析(处理)工具。它的一个重要性质是在时域和频域均具有很好的局部化特征，能够提供目标信号各个频率子段的频率信息。这种信息对于信号分类非常有用。小波变换已成功应用于信号处理、图像处理、模式识别、量子力学、理论物理、计算机分类与识别、音乐与语言的人工合成、医学成像与诊断、地震勘探数据处理以及大型机械的故障诊断等许多方面。
小波变换通过缩放母小波的宽度来获得信号的频率特征，通过平移母小波来获得信号的时间信息。对母小波的缩放和平移操作是为了计算小波系数，这些小波系数反映了小波和局部信号之间的相关程度。
1. 基本原理
（1）连续小波变换（CWT）
小波分析是把一个信号分解为将母小波经过缩放和平移之后的一系列小波，因此小波是小波变换的基函数。小波变换可以理解为用经过缩放和平移的一系列小波函数代替傅里叶变换的正弦波和余弦波进行傅里叶变换的结果。
正弦波和小波不同，正弦波从负无穷一直延续到正无穷，是平滑而且是可预测的，而小波是在有限区间内快速衰减到0的函数，其平均值为0，小波趋于不规则、不对称。从小波和正弦波的形状可以看出，变化剧烈的信号，用不规则的小波进行分析比用平滑的正弦波更好，即用小波更能描述信号的局部特征。
CWT是将任意L2(R)空间下的信号f (t)在小波基下进行展开，定义为
其中，scale是缩放因子；position表示平移参数；ψ( )表示被缩放和平移的小波函数。
基本小波函数ψ( )的缩放操作是指压缩或伸展基本小波，缩放系数越小，则小波越窄。其平移操作是指波的延迟或超前。在数学上，函数f (t)延迟k的表达式为f (t-k)。
CWT的计算过程主要有如下5个步骤：
[1] 将小波ψ(t)和原始信号f (t)的开始部分进行比较。
[2] 计算系数C。该系数表示该部分信号与小波的近似程度。系数C的值越高表示信号与小波越相似，因此系数C可以反映这种波形的相关程度。
[3] 把小波向右移，重复步骤[1]和步骤[2]，直到覆盖整个信号f (t)。若向右移动距离为k，得到的小波函数为ψ(t-k)，再把小波向右移，得到小波ψ(t-2k)。
[4] 扩展小波ψ(t)，重复步骤[1]和步骤[2]。若扩展一倍，则得到的小波函数为ψ(t/2)。
[5] 对所有的缩放，重复步骤[1]～步骤[4]。
小波的缩放因子与信号频率之间的关系是：缩放因子越小，表示小波越窄，度量的是信号的细节变化，表示信号频率越高；缩放因子越大，表示小波越宽，度量的是信号的粗糙程度，表示信号频率越低。
（2）离散小波变换（DWT）
将连续小波变换的缩放因子scale离散化，得到二进小波变换；再将其平移参数position也离散化，就得到离散小波变换。离散小波变换是针对连续小波变换中的缩放因子和平移参数进行离散化，而不是针对时间变量t的离散化。
在计算连续小波变换时，实际上也是用离散的数据进行计算的，只是所用的缩放因子和平移参数比较小而已。不难想象，连续小波变换的计算量是惊人的，而且很多数据是无用的。为了解决计算量的问题，缩放因子和平移参数一般都选择2j (j>0且为整数)的倍数。即选择部分缩放因子和平移参数来进行计算。从而使分析的数据量大大减少。使用这样的缩放因子和平移参数的小波变换叫双尺度小波变换，它是离散小波变换的一种形式。通常离散小波变换就是指双尺度小波变换。
执行离散小波变换的有效方法是使用滤波器。该方法是Mallat在1988年开发的，叫Mallat算法。这种方法实际上是一种信号的分解方法，在数字信号处理中称双通道子带编码。用滤波器执行离散小波变换的概念如图3-9所示。
图3-9 小波分解示意图
图中，S表示原始的输入信号，通过两个互补的滤波器产生A和D信号，A表示信号的近似值（approximation），D表示信号的细节值（detail）。在许多应用中，信号的低频部分是最重要的，而高频部分起一个“添加剂”的作用。比如人的声音，把高频分量去掉之后，听起来声音确实是变了，但还能够听清楚说的是什么内容。相反，如果把低频部分去掉，就什么内容也听不出来。在小波分析中，近似值是大缩放因子产生的系数，表示信号的低频分量。而细节值是小缩放因子产生的系数，表示信号的高频分量。
原始信号经过一对互补的滤波器组进行的分解称为一级分解，信号的分解过程可以迭代，也就是可以进行多级分解。如果对信号的高频分量不再分解，而对低频分量进行连续分解，就可以得到信号不同分辨率下的低频分量，这也称为信号的多分辨率分析。如此进行下去，就会形成的一棵比较大的分解树，如图3-10所示，称为信号的小波分解树。实际中，分解的级数取决于要分析的信号数据特征及用户的具体需要。
图3-10 多级信号分解示意图
图中，Lo_D为低通滤波器；Hi_D为高通滤波器。对于一个信号，如采用多级信号分解的方法，理论上产生的数据量将是原始数据的两倍。所以，根据奈奎斯特（Nyquist）采样定理，可用下采样的方法来减少数据量，即在每个通道内（高通和低通通道）每两个样本数据取一个，便可得到离散小波变换的系数（coefficient），分别用cA和cD表示，如图3-11所示。
图3-11 小波分解示意图
（3）小波重构（IDWT）
小波重构也叫小波合成，就是将信号的小波分解的分量进行处理后，再根据需要把信号恢复出来，也就是利用信号的小波分解的系数还原出原始信号，这一合成过程的数学运算叫做逆离散小波变换（Inverse Discrete Wavelet Transform，IDWT）。
在使用滤波器做小波变换时包含滤波和下采样两个过程，那么在小波重构时要包含上采样和滤波过程。
1) 由小波分解的近似系数和细节系数可以重构出原始信号，方法如图3-12所示。
图3-12 小波重构算法示意图
2) 由近似系数和细节系数分别重构出信号的近似值和细节值，方法如图3-13所示。这时只要近似系数和细节系数置为零即可。
图3-13 重构近似和细节信号示意图
3) 多层重构
根据小波变换多级信号分解的逆运算，在重构出第一层信号的近似值A1与细节值D1之后，则原信号可用A1＋D1＝S重构出来。而近似值A1又可以由第2层的近似值A2与细节值D2相加求出。依次类推，可以把小波变换多层信号重构出来。
在信号重构中，滤波器的选择非常重要，关系到能否重构出满意的原始信号。低通分解滤波器（L）、高通分解滤波器（H）及重构滤波器组（L＇和H＇）构成一个系统，这个系统称为正交镜像滤波器系统。
4）二维离散小波变换
二维离散小波变换是一维离散小波变换的推广，其实质是将二维信号在不同尺度上的分解，得到原始信号的近似值和细节值。由于信号是二维的，因此分解也是二维的。分解的结果为：近似分量cA、水平细节分量cH、垂直细节分量cV和对角细节分量cD。同样也可以利用二维小波分解的结果在不同尺度上重构信号。
对静态二维数字图像的小波变换，一般可先对其进行若干次二维DWT变换，将图像信息分解为高频成分H、V、D和低频成分A。然后针对各部分信息进行相应的处理，最后将各部分处理结果再进行相应的二维变换，得到处理后的图像数据。
2. 算法描述
采用常用的Daubechies系列的小波基，实现离散二维小波变换，完成数字图像的一层小波分解，可以分为以下5个步骤：
[1] 将源图像保存到缓冲区，并记录下缓冲区的地址。
[2] 获取源图像的高度和宽度。
[3] 从Daubechies系列的小波基中，选取一种小波基实现小波变换的单步运算，即每次实现小波的一层变换。
[4] 实现小波变换的单步运算中，调用函数来实现单步二维小波变换，而单步二维小波变换又是通过调单步一维小波变换函数实现的。
[5] 通过小波变换的反变换转换数据，重构原始图像。
3. 编程实现
为实现小波变换，首先采用了Daubechies系列的小波基，共10个系列，其数据列表如下：
const double hCoef[10][20] =
{
{ .707106781187, .707106781187},
{ .482962913145, .836516303738, .224143868042, -.129409522551 },
{ .332670552950, .806891509311, .459877502118, -.135011020010, -85441273882, .035226291882 },
{ .230377813309, .714846570553, .630880767930, -.027983769417,
-.187034811719, .030841381836, .032883011667, -.010597401785 },
{ .160102397974, .603829269797, .724308528438, .138428145901, -.242294887066,
.032244869585, .077571493840, -.006241490213, -.012580751999, .003335725285 },
{ .111540743350, .494623890398, .751133908021, .315250351709, -.226264693965,
-.129766867567, .097501605587, .027522865530, -.031582039318, .000553842201,
.004777257511, -.001077301085 },
{ .077852054085, .396539319482, .729132090846, .469782287405, -.143906003929,
-.224036184994, .071309219267, .080612609151, -.038029936935, -.016574541631,
.012550998556, .000429577973, -.001801640704, .000353713800 },
{ .054415842243, .312871590914, .675630736297, .585354683654, -.015829105256,
-.284015542962, .000472484574, .128747426620, -.017369301002, -.044088253931,
.013981027917, .008746094047, -.004870352993, -.000391740373, .000675449406,
-.000117476784 },
{ .038077947364, .243834674613, .604823123690, .657288078051, .133197385825,
-.293273783279, -.096840783223, .148540749338, .030725681479, -.067632829061,
.000250947115, .022361662124, -.004723204758, -.004281503682, .001847646883,
.000230385764, -.000251963189, .000039347320 },
{ .026670057901, .188176800078, .527201188932, .688459039454, .281172343661,
-.249846424327, -.195946274377, .127369340336, .093057364604, -.071394147166,
-.029457536822, .033212674059, .003606553567, -.010733175483, .001395351747,
.001992405295, -.000685856695, -.000116466855, .000093588670, -.000013264203 }
};
下面通过编写VC的函数DIBDWTStep( )来实现图像的小波变换，其具体代码如下：
/＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊ 函数名称： DIBDWTStep
＊ 函数参数：
＊ LPSTR lpDIBBits 指向源数据的指针
＊ double ＊m_pDbImage 小波变换的存储空间
＊ int nWidth 图像的宽度
＊ int nHeight 图像的高度
＊ int nInv 是否为DWT，1表示为IDWT，0表示为DWT
＊ int m_nDWTCurDepth 当前已经分解的层数
＊ int nSupp 小波基的紧支集的长度
＊ 函数功能：该函数对图像进行单步小波分解，分解结果数据仍然保存在lpDIBBits中
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊/
BOOL DIBDWTStep(LPSTR lpDIBBits, double＊m_pDbImage, int nWidth, int nHeight,
int nInv, int m_nDWTCurDepth, int m_nSupp)
{
int i, j; // 循环变量
// 获取变换的最大层数
int nMaxWLevel = Log2(nWidth);
int nMaxHLevel = Log2(nHeight);
int nMaxLevel;
if (nWidth == 1<<nMaxWLevel && nHeight == 1<<nMaxHLevel)
nMaxLevel = min(nMaxWLevel, nMaxHLevel);
double ＊pDbTemp; // 临时变量
BYTE ＊pBits;
// 如果小波变换的存储内存还没有分配，则分配此内存
if(!m_pDbImage)
{
m_pDbImage = new double[nWidth＊nHeight];
if (!m_pDbImage)return FALSE; // 将图像数据放入m_pDbImage中
for (j=0; j<nHeight; j++)
{
pDbTemp = m_pDbImage + j＊nWidth;
pBits = (unsigned char ＊)lpDIBBits + (nHeight-1-j)＊nWidth;
for (i=0; i<nWidth; i++)
pDbTemp[i] = pBits[i];
}
}
// 进行小波变换（或反变换）
if (!DWTStep_2D(m_pDbImage, nMaxWLevel-m_nDWTCurDepth, nMaxHLevel
-m_nDWTCurDepth,nMaxWLevel, nMaxHLevel, nInv, 1, m_nSupp))
return FALSE;
if (nInv)
m_nDWTCurDepth --; // 如果是反变换，则当前层数减1
else
m_nDWTCurDepth ++; // 否则加1
// 然后，将数据复制回原CDib中，并进行相应的数据转换
int lfw = nWidth>>m_nDWTCurDepth, lfh = nHeight>>m_nDWTCurDepth;
for (j=0; j<nHeight; j++)
{
pDbTemp = m_pDbImage + j＊nWidth;
pBits = (unsigned char ＊)lpDIBBits + (nHeight-1-j)＊nWidth;
for (i=0; i<nWidth; i++)
{
if (j<lfh && i<lfw)
pBits[i] = FloatToByte(pDbTemp[i]);
else
pBits[i] = BYTE(FloatToChar(pDbTemp[i]) ^ 0x80);
}
}
return TRUE;
}
其中,单步二维小波变换是通过编写VC的函数DWTStep_2D ( )实现的，具体代码如下：
/＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊ 函数名称： DWTStep_2D
＊ 输入参数：
＊ double ＊ pDbSrc 指向源数据的指针
＊ int nCurWLevel x方向上当前分解的层数
＊ int nCurHLevel y方向上当前分解的层数
＊ int nMaxWLevel x方向上最大可分解的层数
＊ int nMaxHLevel y方向上最大可分解的层数
＊ int nInv 是否为DWT，1表示为IDWT，0表示DWT
＊ int nStep 当前的计算层数
＊ int nSupp 小波基的紧支集的长度
＊ 说明：
＊ 该函数用于对存放在pDbSrc中的数据进行一层的二维DWT或者IDWT。
＊ 计算后数据仍存放在pDbSrc中
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊/
BOOL DWTStep_2D(double＊ pDbSrc, int nCurWLevel, int nCurHLevel,int nMaxWLevel,
int nMaxHLevel, int nInv, int nStep, int nSupp)
{
int W = 1<<nMaxWLevel, H = 1<<nMaxHLevel; //计算图像的长度和宽度（2次幂对齐）
int CurW = 1<<nCurWLevel, CurH = 1<<nCurHLevel; //计算当前分解的图像的长度和宽度
if (!nInv) // 判断是进行DWT，还是IDWT
{
int i = 0;
for (i=0; i<CurH; i++) // 对行进行一维DWT
if (!DWTStep_1D(pDbSrc+(int)i＊W＊nStep, nCurWLevel, nInv, nStep, nSupp))
return FALSE;
for (i=0; i<CurW; i++) // 对列进行一维DWT
if (!DWTStep_1D(pDbSrc+i＊nStep, nCurHLevel, nInv, W＊nStep, nSupp))
return FALSE;
}
else // 否则进行IDWT
{
// 计算当前变换的图像的长度和宽度
CurW <<= 1;
CurH <<= 1;
int i = 0;
for (i=0; i<CurW; i++) // 对列进行IDWT
if (!DWTStep_1D(pDbSrc+i＊nStep, nCurHLevel, nInv, W＊nStep, nSupp))
return FALSE;
for (i=0; i<CurH; i++) // 对行进行IDWT
if (!DWTStep_1D(pDbSrc+(int)i＊W＊nStep, nCurWLevel, nInv, nStep, nSupp))
return FALSE;
}
return TRUE;
}
其中，单步一维小波变换是通过编写VC的函数DWTStep_1D ( )实现，具体代码如下：
/＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
＊ 函数名称： DWTStep_1D
＊ 输入参数：
＊ double ＊ pDbSrc 指向源数据的指针
＊ int nCurLevel 当前分界的层数
＊ int nInv 是否为DWT，1表示为IDWT，0表示DWT
＊ int nStep 当前的计算层数
＊ int nSupp 小波基的紧支集的长度
＊ 说明：
＊ 该函数用于对存放在pDBSrc中的数据进行一层的一维DWT或者IDWT。其中，nInv为表示进行
＊ DWT或者IDWT的标志；nCurLevel为当前需要进行分界的层数；nStep为已经分界的层数
＊ 计算后数据仍存放在pDbSrc中
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊/
BOOL DWTStep_1D(double＊ pDbSrc, int nCurLevel, int nInv, int nStep,int nSupp)
{
double s = sqrt((double)2); // 获得小波基的指针
double＊ h = (double＊)hCoef[nSupp-1]; // 确认当前层数有效
ASSERT(nCurLevel>=0);
int CurN = 1<<nCurLevel; // 计算当前层数的长度
if (nInv) CurN <<= 1;
if (nSupp<1 || nSupp>10 || CurN<2＊nSupp)
return FALSE; // 确认所选择的小波基和当前层数的长度有效
double ＊ptemp = new double[CurN]; // 分配临时内存用于存放结果
if (!ptemp)
return FALSE;
double s1, s2;
int Index1, Index2;
if (!nInv) // 判断是进行DWT，还是IDWT
{
Index1=0; // DWT
Index2=2＊nSupp-1;
// 进行卷积，其中s1为低频部分的结果，s2为高频部分的结果
for (int i=0; i<CurN/2; i++)
{
s1 = s2 = 0;
double t = -1;
for (int j=0; j<2＊nSupp; j++, t=-t)
{
s1 += h[j]＊pDbSrc[(Index1 & CurN-1) ＊ nStep];
s2 += t＊h[j]＊pDbSrc[(Index2 & CurN-1) ＊ nStep];
Index1++;
Index2--;
}
ptemp[i] = s1/s; // 将结果存放在临时内存中
ptemp[i+CurN/2] = s2/s;
Index1-= 2＊nSupp;
Index2 += 2＊nSupp;
Index1 += 2;
Index2 += 2;
}
}
else // 否则进行IDWT
{
Index1 = CurN/2; // IDWT
Index2 = CurN/2-nSupp+1;
// 进行卷积，其中s1为低频部分的结果，s2为高频部分的结果
for (int i=0; i<CurN/2; i++)
{
s1 = s2 = 0;
int Index3 = 0;
for (int j=0; j<nSupp; j++)
{
s1 += h[Index3]＊pDbSrc[(Index1 & CurN/2-1) ＊ nStep]
+h[Index3+1]＊pDbSrc[((Index2 & CurN/2-1) + CurN/2) ＊ nStep];
s2 += h[Index3+1]＊pDbSrc[(Index1 & CurN/2-1) ＊ nStep]
-h[Index3]＊pDbSrc[((Index2 & CurN/2-1) + CurN/2) ＊ nStep];
Index3+=2;
Index1--, Index2++;
}
ptemp[2＊i] = s1＊s; // 将结果存入临时内存
ptemp[2＊i+1] = s2＊s;
Index1 += nSupp;
Index2-= nSupp;
Index1++;
Index2++;
}
}
for (int i=0; i<CurN; i++) // 将结果存入源图像中
pDbSrc[i＊nStep] = ptemp[i];
delete[] ptemp; // 释放临时内存，并返回
return TRUE;
}
可以通过在数字图像处理程序的框架窗口上增加菜单【正交变换】|【小波变换】。对应视图类中的CDImageProcessView:: OnWl ( )进行函数的调用：
void CDImageProcessView::OnWl()
{
CDImageProcessDoc＊ pDoc = GetDocument(); // 获取文档
LPSTR lpDIB; // 指向源图像的指针
LPSTR lpDIBBits; // 指向DIB像素指针
lpDIB = (LPSTR) ::GlobalLock((HGLOBAL) pDoc->GetHObject());//锁定DIB
lpDIBBits = pDoc->m_dib.GetBits(lpDIB); // 找到DIB图像像素起始位置
if (pDoc->m_dib.GetColorNum(lpDIB) != 256) // 判断是否是256色位图
{
MessageBox("目前只支持256色位图！", "系统提示" , MB_ICONINFORMATION | MB_OK);
::GlobalUnlock((HGLOBAL) pDoc->GetHObject());// 解除锁定
return; // 返回
}
int nWidth; // 图像宽度
int nHeight; // 图像高度
nWidth=(int)pDoc->m_dib.GetWidth(lpDIB); // 获取图像宽度
nHeight=(int)pDoc->m_dib.GetHeight(lpDIB); // 获取图像高度
DIBDWTStep(lpDIBBits,m_pDbImage, nWidth,nHeight, 0, m_nDWTCurDepth,1);
pDoc->UpdateAllViews(NULL);
::GlobalUnlock((HGLOBAL) pDoc->GetHObject());
}
4. 效果演示
源图像如图3-14a所示，对其进行一层小波分解，结果如图3-14b所示。
图3-14 一层小波分解效果演示
3.3 综合实例—特征提取
本节将通过一个具体的VC++工程实例—特征提取，实现图像正交变换中的基本正交变换和特征变换过程，其主要程序是在第2章实例程序的框架基础之上编制的。
【例3-1】图像正交变换综合实例—特征提取
[1] 打开Visual Studio 2005，打开第2章工程项目“DImageProcess”。
[2] 为程序添加【正交变换】菜单。在【正交变换】菜单下分别设置了【离散傅里叶变换】、【离散余弦变换】、【沃尔什哈达玛】、【离散K-L变换】、【SVD变换】和【小波变换】6个子菜单，如图3-15所示。
图3-15 【正交变换】菜单
[3] 更改各个菜单项的ID号。新增菜单对应的ID号设定如表3-1所示。
表3-1 新增菜单对应的ID号
[4] 设计对话框。根据程序设置，设计【SVD变换】菜单所对应的【使用SVD算法实现图像压缩】对话框，如图3-16所示。
图3-16 【使用SVD算法实现图像压缩】对话框
[5] 设定【使用SVD算法实现图像压缩】对话框的类名以及对话框上主要控件所对应ID号，如表3-2所示。
表3-2 各对话框名称及其主要控件对应的ID号
[6] 为对话框添加变量及编程。
为CDlgsvd类添加int型变量m_svd，并且增加如下代码初始化变量：
m_svd = 50;
为CDlgsvd类添加OnInitDialog函数，并且增加如下代码：
CSpinButtonCtrl＊ pSpinsvd=(CSpinButtonCtrl＊)GetDlgItem(IDC_SPIN_svd);
pSpinsvd->SetRange(1,100);
为CDlgsvd类添加DoDataExchange函数，并且增加如下代码：
DDX_Text(pDX, IDC_rot_svd, m_svd);
[7] 在“function.h”头文件中加入图像正交变换各个函数以及一些需要调用的计算函数，其具体代码如下：
// 函数功能：将n阶实对称矩阵化为对称三角阵
BOOL ThreeCrossMat(double ＊pMatrix, int rank, double ＊pQMatrix, double ＊pMainCross,
double ＊pHypoCross)
{
int i, j, k, u; // 变量声明
double h, f, g, h2;
for(i = 0; i <= rank-1; i++) // 将矩阵pQMatrix初始化
{
for(j = 0; j <= rank-1; j++)
{
u = i＊rank + j;
pQMatrix[u] = pMatrix[u];
}
}
for (i = rank-1; i >= 1; i--)
{
h = 0.0;
if (i > 1)
for (k = 0; k <= i-1; k++)
{
u = i＊rank + k;
h = h + pQMatrix[u]＊pQMatrix[u];
}
if (h + 1.0 == 1.0) // 如果一行全部为零
{
pHypoCross[i] = 0.0;
if (i == 1)
{
pHypoCross[i] = pQMatrix[i＊rank+i-1];
}
pMainCross[i] = 0.0;
}
else // 否则求正交矩阵的值
{
pHypoCross[i] = sqrt(h); // 求次对角元素的值
u = i＊rank + i -1;
if (pQMatrix[u] > 0.0) // 判断第i行、第i-1列元素是否大于零
{
pHypoCross[i] = -pHypoCross[i];
}
h = h - pQMatrix[u]＊pHypoCross[i];
pQMatrix[u] = pQMatrix[u] - pHypoCross[i];
f = 0.0;
for (j = 0; j <= i-1; j++) // householder变换
{
pQMatrix[j＊rank+i] = pQMatrix[i＊rank+j] / h;
g = 0.0;
for (k = 0; k <= j; k++)
{
g = g + pQMatrix[j＊rank+k]＊pQMatrix[i＊rank+k];
}
if (j+1 <= i-1)
for (k = j+1; k <= i-1; k++)
{
g = g + pQMatrix[k＊rank+j]＊pQMatrix[i＊rank+k];
}
pHypoCross[j] = g / h;
f = f + g＊pQMatrix[j＊rank+i];
}
h2 = f / (h + h);
for (j = 0; j <= i-1; j++) // 求正交矩阵的值
{
f = pQMatrix[i＊rank + j];
= pHypoCross[j] - h2＊f;
pHypoCross[j] = g;
for (k = 0; k <= j; k++)
{
u = j＊rank + k;
pQMatrix[u] = pQMatrix[u] - f＊pHypoCross[k]-g＊pQMatrix[i＊rank + k];
}
}
pMainCross[i] = h;
}
}
for (i = 0; i <= rank-2; i++) // 赋零值
{
pHypoCross[i] = pHypoCross[i + 1];
}
pHypoCross[rank -1] = 0.0;
pMainCross[0] = 0.0;
for (i = 0; i <= rank-1; i++)
{
if ((pMainCross[i] != 0.0) && (i-1 >= 0)) // 主对角元素的计算
for (j = 0; j <= i-1; j++)
{
g = 0.0;
for (k = 0; k <= i-1; k++)
g = g + pQMatrix[i＊rank + k]＊pQMatrix[k＊rank + j];
for (k = 0; k <= i-1; k++)
{
u = k＊rank + j;
pQMatrix[u] = pQMatrix[u] - g＊pQMatrix[k＊rank + i];
}
}
u = i＊rank + i; // 存储主对角线的元素
pMainCross[i] = pQMatrix[u];
pQMatrix[u] = 1.0;
if (i-1 >= 0) // 将三对角外所有的元素赋零值
for (j = 0; j <= i-1; j++)
{
pQMatrix[i＊rank + j] = 0.0;
pQMatrix[j＊rank+i] = 0.0;
}
}
return(TRUE);
}
// 函数功能：计算实对称三角矩阵的全部特征值及相应的特征向量
BOOL EigenvalueVector(int rank, double ＊pMainCross, double ＊pHypoCross,
double ＊pMatrix, double Precision, int MaxT)
{
int i, j, k, m, it, u, v; // 变量声明
double d, f, h, g, p, r, e, s;
pHypoCross[rank -1] = 0.0; // 初始化
d = 0.0;
f = 0.0;
for(j = 0; j <= rank-1; j++) // 迭代精度的控制
{
it = 0;
h = Precision ＊ (fabs(pMainCross[j]) + fabs(pHypoCross[j]));
if(h > d)
{
d = h;
}
m = j;
while((m <= rank-1) && (fabs(pHypoCross[m]) > d))
{
m = m + 1;
}
if(m != j) // 迭代求矩阵A的特征值和特征向量
{
do
{
if(it == MaxT) // 超过迭代次数，迭代失败
{
return(FALSE);
}
it = it + 1;
g = pMainCross[j];
p = (pMainCross[j + 1] - g) / (2.0 ＊ pHypoCross[j]);
r = sqrt(p＊p + 1.0);
if (p >= 0.0) // 如果p大于等于0
{
pMainCross[j] = pHypoCross[j]/(p + r);
}
else
{
pMainCross[j] = pHypoCross[j]/(p - r);
}
h = g - pMainCross[j]; // 计算主对角线的元素
for (i = j + 1; i <= rank -1; i++)
{
pMainCross[i] = pMainCross[i] - h;
}
f = f + h; // 赋值
p = pMainCross[m];
e = 1.0; s = 0.0;
for(i = m -1; i >= j; i--)
{
g = e ＊ pHypoCross[i];
h = e ＊ p;
// 主对角线元素的绝对值是否大于次对角线元素的绝对值
if(fabs(p) >= fabs(pHypoCross[i]))
{
e = pHypoCross[i] / p;
r = sqrt(e＊e + 1.0);
pHypoCross[i + 1] = s＊p＊r;
s = e / r; e = 1.0 / r;
}
else
{
e = p / pHypoCross[i];
r = sqrt(e＊e + 1.0);
pHypoCross[i+1] = s ＊ pHypoCross[i] ＊ r;
s = 1.0 / r; e = e / r;
}
p = e＊pMainCross[i] - s＊g;
pMainCross[i + 1] = h + s＊(e＊g + s＊pMainCross[i]);
for(k = 0; k <= rank -1; k++) // 重新存储特征向量
{
u = k＊rank + i + 1; v = u -1;
h = pMatrix[u];
pMatrix[u] = s＊pMatrix[v] + e＊h;
pMatrix[v] = e＊pMatrix[v] - s＊h;
}
}
pHypoCross[j] = s ＊ p; // 将主对角线和次对角线元素重新赋值
pMainCross[j] = e ＊ p;
} while (fabs(pHypoCross[j]) > d);
}
pMainCross[j] = pMainCross[j] + f;
}
for (i = 0; i <= rank-1; i++)
{
k = i; p = pMainCross[i]; // 返回A的特征值
if(i+1 <= rank-1) // 将A的特征值赋给p
{
j = i + 1;
while((j <= rank-1) && (pMainCross[j] <= p))
{
k = j;
p = pMainCross[j];
j = j+1;
}
}
if (k != i)
{
pMainCross[k] = pMainCross[i]; // 存储A的特征值和特征向量
pMainCross[i] = p;
for(j = 0; j <= rank-1; j++)
{
u = j＊rank + i;
v = j＊rank + k;
p = pMatrix[u];
pMatrix[u] = pMatrix[v];
pMatrix[v] = p;
}
}
}
return(TRUE);
}
// 该函数求取输入参数的以2为底的对数，并转换为整型输出。
int Log2(int n)
{
int rsl = 0;
while (n >>= 1) rsl++;
return rsl;
}
void FastFourierTran(CplexNum ＊ pTd, CplexNum ＊ pFd, int power)
{⋯⋯}
void DisFCosTran(double ＊pTd, double ＊pFd, int power)
{⋯⋯}
BOOL CosTran(LPSTR lpSrcStartBits, long lWidth, long lHeight,long lLineBytes)
{⋯⋯}
void Walshei_Har(double ＊pTd, double ＊pFd, int power)
{⋯⋯}
BOOL Walsh_HarTran(LPSTR lpSrcStartBits, long lWidth, long lHeight,long lLineBytes)
{⋯⋯}
BOOL DisK_L(LPSTR lpSrcStartBits, long lWidth, long lHeight,long lLineBytes)
{⋯⋯}
BOOL DIBDWTStep(LPSTR lpDIBBits,double＊m_pDbImage, int nWidth,int nHeight, int nInv,
int m_nDWTCurDepth,int m_nSupp)
{⋯⋯}
BOOL DWTStep_1D(double＊ pDbSrc, int nCurLevel, int nInv, int nStep,int nSupp)
{⋯⋯}
BOOL DWTStep_2D(double＊ pDbSrc, int nCurWLevel, int nCurHLevel,
int nMaxWLevel, int nMaxHLevel, int nInv, int nStep, int nSupp)
{⋯⋯}
[8] 完善视类文件CDImageProcessView的代码。
在视类实现文件DImageProcessView.cpp中加入对本章所需要的头文件的引用，其具体代码如下：
#include "Dlgsvd.h"
#include "matlib.h"
[9] 为增加的菜单添加消息处理函数，将这些函数添加到DimageProcessView.cpp视图类中。具体形式如下：
void CDImageProcessView::OnDct() {⋯⋯}
void CDImageProcessView::OnWalhHar() {⋯⋯}
void CDImageProcessView::OnFourier() {⋯⋯}
void CDImageProcessView::OnKL() {⋯⋯}
void CDImageProcessView::OnWl() {⋯⋯}
void CDImageProcessView::Svd() {⋯⋯}
消息处理函数的具体内容已在图像正交变换的相关节中给出。
[10] 按F5键或工具栏中的运行按钮
运行程序。
打开一幅图片，此时程序界面上出现该图片，并且在程序菜单栏上出现【正交变换】菜单，如图3-17所示。在【正交变换】菜单下选择要对该图片进行的操作，例如选择【小波变换】，对图片进行小波变换，效果如图3-18所示。
图3-17 打开图像窗口
图3-18 图像小波变换
3.4 实践拓展
1. 怎么样才能避免重复包含同一头文件
在程序运行过程中，头文件被多次包含是比较常见的现象，这直接导致了程序运行效率的下降，为了防止此类问题的发生，Visual Studio 2005一般会在头文件中加入如下代码：
#pragma once
#pragma once是由编译器提供保证同一个头文件不会被包含多次的一种方式。注意这里所说的“同一个文件”是指物理上的一个文件，而不是指内容相同的两个文件。所以不必考虑重新定义宏名，也就不会出现宏名相同引发的问题。其缺点是如果某个头文件有多份复制，则不能保证不被重复包含。
在一些VC程序中还可能看到，类似以下的代码：
#if !defined(AFX_FUNCTION_H_6E194843_FEB3_491F_8062_765AA3465CBC__INCLUDED_)
#define AFX_FUNCTION_H__6E194843_FEB3_491F_8062_765AA3465CBC__INCLUDED_
#if _MSC_VER > 1000
#pragma once
#endif
这段代码多用于VC以前的版本，作用也是防止该头文件被多次包含，提高程序运行效率。它的优点是：#ifndef的方式依赖于宏名字不能冲突，不仅可以保证同一个文件不会被包含多次，也能保证内容完全相同的两个文件不会被同时包含。它的缺点是：如果不同头文件的宏名相同，可能就会导致头文件明明存在，编译器却显示找不到声明的情况。
2. 怎样提高变换后图像的质量
数字图像变换后所得图像的质量取决于变换算法使用和计算值精确程度。对于离散的数字图像来说，无论图像的几何变换还是正交变换都不可能得到绝对正确的变换后图像。例如，在图像的几何变换中进行图像缩放时，计算缩放后图像的实际宽度和实际高度代码如下：
lDstWidth= (long) (lSrcWidth＊fX + 0.5); // 计算缩放后的图像实际宽度
lDstHeight= (long) (lSrcHeight ＊ fY + 0.5); // 计算缩放后的图像高度
可以看到代码中加上0.5后取整，就是为了在强制类型转换时不使用四舍五入，而是直接截去小数部分，保留整数部分。
在图像的正交变换中，许多变换（如离散傅里叶变换）本身就是一种近似的运算。在SVD变换中，经过多次试验可以得出：一般情况下，对于256≤n≤2048的图像进行压缩，选取的奇异值个数在25≤k≤100时，可以得到比较满意的视觉效果，这样既满足了图像压缩后的质量要求，又节省了大量的矩阵运算时间。
总之，在使用VC程序进行数字图像处理时，要保持程序运行速度与变换后图像质量之间的一个平衡，不能顾此失彼。确切地说，对于一个算法，会有很多编程的方法，程序的数据结构和流程不同，最终得到的结果也不同。比如在需要运用插值算法进行图像变换时，如果运用双线性插值或双三次插值的方法，会得到更好的变换后的图像，但同时也增加了程序运行的时间。本章对于图像变换所给出的程序只是针对以各种算法为基础的众多编程方法中的一种，仅是给读者提供了一个图像变换编程方法的思路，读者可以结合自身的编程水平，对数字图像变换的VC程序进行更有效的扩展和实践，编制出运行时间更短、运行结果更好、运行效率更高的程序。