www.pudn.com > 图像识别.rar > FreTrans.cpp
#include "stdafx.h"
#include "GlobalApi.h"
#include "Cdib.h"
#include
#include
#include
using namespace std;
// FOURBYTES就是用来计算离4最近的整倍数
#define FOURBYTES(bits) (((bits) + 31) / 32 * 4)
/**************************************************************************
* 文件名：FreTrans.cpp
*
* 正交变换API函数库：
*
* THREECROSS() - 将实对称矩阵化作三对角矩阵
* BSTQ() - 求三对角对称矩阵的特征值和特征向量
* FFT_1D() - 快速一维傅立叶变换
* IFFT_1D() - 快速一维傅立叶反变换
* IFFT_2D() - 快速二维傅立叶反变换
* DCT() - 离散余弦变换
* IDCT() - 离散余弦反变换
* WALSH() - 沃尔什－哈达玛变换
* IWALSH() - 沃尔什－哈达玛反变换
*
*
* DIBFFT_2D() - 图象的二维离散傅立叶快速变换
* DIBDFT_2D() - 图象的二维离散傅立叶变换
* DIBHOTELLING() - 图象的霍特林变换
* DIBDct() - 图像的离散余弦变换
* DIBWalsh() - 图像的沃尔什－哈达玛变换
*
************************************************************************
*/
/*************************************************************************
*
* \函数名称：
* FFT_1D()
*
* \输入参数:
* complex * pCTData - 指向时域数据的指针，输入的需要变换的数据
* complex * pCFData - 指向频域数据的指针，输出的经过变换的数据
* nLevel －傅立叶变换蝶形算法的级数，2的幂数，
*
* \返回值:
* 无
*
* \说明:
* 一维快速傅立叶变换。
*
*************************************************************************
*/
VOID WINAPI FFT_1D(complex * pCTData, complex * pCFData, int nLevel)
{
// 循环控制变量
int i;
int j;
int k;
double PI = 3.1415926;
// 傅立叶变换点数
int nCount =0 ;
// 计算傅立叶变换点数
nCount =(int)pow(2,nLevel) ;
// 某一级的长度
int nBtFlyLen;
nBtFlyLen = 0 ;
// 变换系数的角度 ＝2 * PI * i / nCount
double dAngle;
complex *pCW ;
// 分配内存，存储傅立叶变化需要的系数表
pCW = new complex[nCount / 2];
// 计算傅立叶变换的系数
for(i = 0; i < nCount / 2; i++)
{
dAngle = -2 * PI * i / nCount;
pCW[i] = complex ( cos(dAngle), sin(dAngle) );
}
// 变换需要的工作空间
complex *pCWork1,*pCWork2;
// 分配工作空间
pCWork1 = new complex[nCount];
pCWork2 = new complex[nCount];
// 临时变量
complex *pCTmp;
// 初始化，写入数据
memcpy(pCWork1, pCTData, sizeof(complex) * nCount);
// 临时变量
int nInter;
nInter = 0;
// 蝶形算法进行快速傅立叶变换
for(k = 0; k < nLevel; k++)
{
for(j = 0; j < (int)pow(2,k); j++)
{
//计算长度
nBtFlyLen = (int)pow( 2,(nLevel-k) );
//倒序重排，加权计算
for(i = 0; i < nBtFlyLen/2; i++)
{
nInter = j * nBtFlyLen;
pCWork2[i + nInter] =
pCWork1[i + nInter] + pCWork1[i + nInter + nBtFlyLen / 2];
pCWork2[i + nInter + nBtFlyLen / 2] =
(pCWork1[i + nInter] - pCWork1[i + nInter + nBtFlyLen / 2])
* pCW[(int)(i * pow(2,k))];
}
}
// 交换 pCWork1和pCWork2的数据
pCTmp = pCWork1 ;
pCWork1 = pCWork2 ;
pCWork2 = pCTmp ;
}
// 重新排序
for(j = 0; j < nCount; j++)
{
nInter = 0;
for(i = 0; i < nLevel; i++)
{
if ( j&(1< * pCTData - 指向时域数据的指针，输入的需要反变换的数据
* complex * pCFData - 指向频域数据的指针，输出的经过反变换的数据
* nLevel －傅立叶变换蝶形算法的级数，2的幂数，
*
* \返回值:
* 无
*
* \说明:
* 一维快速傅立叶反变换。
*
************************************************************************
*/
VOID WINAPI IFFT_1D(complex * pCFData, complex * pCTData, int nLevel)
{
// 循环控制变量
int i;
// 傅立叶反变换点数
int nCount;
// 计算傅立叶变换点数
nCount = (int)pow(2,nLevel) ;
// 变换需要的工作空间
complex *pCWork;
// 分配工作空间
pCWork = new complex[nCount];
// 将需要反变换的数据写入工作空间pCWork
memcpy(pCWork, pCFData, sizeof(complex) * nCount);
// 为了利用傅立叶正变换,可以把傅立叶频域的数据取共轭
// 然后直接利用正变换，输出结果就是傅立叶反变换结果的共轭
for(i = 0; i < nCount; i++)
{
pCWork[i] = complex (pCWork[i].real(), -pCWork[i].imag());
}
// 调用快速傅立叶变换实现反变换，结果存储在pCTData中
FFT_1D(pCWork, pCTData, nLevel);
// 求时域点的共轭，求得最终结果
// 根据傅立叶变换原理，利用这样的方法求得的结果和实际的时域数据
// 相差一个系数nCount
for(i = 0; i < nCount; i++)
{
pCTData[i] =
complex (pCTData[i].real() / nCount, -pCTData[i].imag() / nCount);
}
// 释放内存
delete pCWork;
pCWork = NULL;
}
/*************************************************************************
*
* \函数名称：
* FFT_2D()
*
* \输入参数:
* complex * pCTData - 图像数据
* int nWidth - 数据宽度
* int nHeight - 数据高度
* complex * pCFData - 傅立叶变换后的结果
*
* \返回值:
* 无
*
* \说明:
* 二维傅立叶变换。
*
************************************************************************
*/
VOID WINAPI DIBFFT_2D(complex * pCTData, int nWidth, int nHeight, complex * pCFData)
{
// 循环控制变量
int x;
int y;
// 临时变量
double dTmpOne;
double dTmpTwo;
// 进行傅立叶变换的宽度和高度，（2的整数次幂）
// 图像的宽度和高度不一定为2的整数次幂
int nTransWidth;
int nTransHeight;
// 计算进行傅立叶变换的宽度 （2的整数次幂）
dTmpOne = log(nWidth)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransWidth = (int) dTmpTwo ;
// 计算进行傅立叶变换的高度 （2的整数次幂）
dTmpOne = log(nHeight)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransHeight = (int) dTmpTwo ;
// x，y（行列）方向上的迭代次数
int nXLev;
int nYLev;
// 计算x，y（行列）方向上的迭代次数
nXLev = (int) ( log(nTransWidth)/log(2) + 0.5 );
nYLev = (int) ( log(nTransHeight)/log(2) + 0.5 );
for(y = 0; y < nTransHeight; y++)
{
// x方向进行快速傅立叶变换
FFT_1D(&pCTData[nTransWidth * y], &pCFData[nTransWidth * y], nXLev);
}
// pCFData中目前存储了pCTData经过行变换的结果
// 为了直接利用FFT_1D，需要把pCFData的二维数据转置，再一次利用FFT_1D进行
// 傅立叶行变换（实际上相当于对列进行傅立叶变换）
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTData[nTransHeight * x + y] = pCFData[nTransWidth * y + x];
}
}
for(x = 0; x < nTransWidth; x++)
{
// 对x方向进行快速傅立叶变换，实际上相当于对原来的图象数据进行列方向的
// 傅立叶变换
FFT_1D(&pCTData[x * nTransHeight], &pCFData[x * nTransHeight], nYLev);
}
// pCFData中目前存储了pCTData经过二维傅立叶变换的结果，但是为了方便列方向
// 的傅立叶变换，对其进行了转置，现在把结果转置回来
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTData[nTransWidth * y + x] = pCFData[nTransHeight * x + y];
}
}
memcpy(pCTData, pCFData, sizeof(complex) * nTransHeight * nTransWidth );
}
/*************************************************************************
*
* \函数名称：
* IFFT_2D()
*
* \输入参数:
* complex * pCFData - 频域数据
* complex * pCTData - 时域数据
* int nWidth - 图像数据宽度
* int nHeight - 图像数据高度
*
* \返回值:
* 无
*
* \说明:
* 二维傅立叶反变换。
*
************************************************************************
*/
VOID WINAPI IFFT_2D(complex * pCFData, complex * pCTData, int nWidth, int nHeight)
{
// 循环控制变量
int x;
int y;
// 临时变量
double dTmpOne;
double dTmpTwo;
// 进行傅立叶变换的宽度和高度，（2的整数次幂）
// 图像的宽度和高度不一定为2的整数次幂
int nTransWidth;
int nTransHeight;
// 计算进行傅立叶变换的宽度 （2的整数次幂）
dTmpOne = log(nWidth)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransWidth = (int) dTmpTwo ;
// 计算进行傅立叶变换的高度 （2的整数次幂）
dTmpOne = log(nHeight)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransHeight = (int) dTmpTwo ;
// 分配工作需要的内存空间
complex *pCWork= new complex[nTransWidth * nTransHeight];
//临时变量
complex *pCTmp ;
// 为了利用傅立叶正变换,可以把傅立叶频域的数据取共轭
// 然后直接利用正变换，输出结果就是傅立叶反变换结果的共轭
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTmp = &pCFData[nTransWidth * y + x] ;
pCWork[nTransWidth * y + x] = complex( pCTmp->real() , -pCTmp->imag() );
}
}
// 调用傅立叶正变换
::DIBFFT_2D(pCWork, nWidth, nHeight, pCTData) ;
// 求时域点的共轭，求得最终结果
// 根据傅立叶变换原理，利用这样的方法求得的结果和实际的时域数据
// 相差一个系数
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTmp = &pCTData[nTransWidth * y + x] ;
pCTData[nTransWidth * y + x] =
complex( pCTmp->real()/(nTransWidth*nTransHeight),
-pCTmp->imag()/(nTransWidth*nTransHeight) );
}
}
delete pCWork ;
pCWork = NULL ;
}
/*************************************************************************
*
* 函数名称：
* DCT()
*
* 参数:
* double * dpf - 指向时域值的指针
* double * dpF - 指向频域值的指针
* r －2的幂数
*
* 返回值:
* 无。
*
* 说明:
* 实现一维快速离散余弦变换。
*
***********************************************************************
*/
VOID WINAPI DCT(double *dpf, double *dpF, int r)
{
double PI = 3.1415926;
// 离散余弦变换点数
LONG lNum;
// 循环变量
int i;
// 中间变量
double dTemp;
complex *comX;
// 离散余弦变换点数
lNum = 1<[lNum*2];
// 赋初值0
memset(comX, 0, sizeof(complex) * lNum * 2);
// 将时域点转化为复数形式
for(i=0;i (dpf[i], 0);
}
// 调用快速付立叶变换
FFT_1D(comX,comX,r+1);
// 调整系数
dTemp = 1/sqrt(lNum);
// 求dpF[0]
dpF[0] = comX[0].real() * dTemp;
dTemp *= sqrt(2);
// 求dpF[u]
for(i = 1; i < lNum; i++)
{
dpF[i]=(comX[i].real() * cos(i*PI/(lNum*2))
+ comX[i].imag() * sin(i*PI/(lNum*2))) * dTemp;
}
// 释放内存
delete comX;
}
/*************************************************************************
*
* 函数名称：
* IDCT()
*
* 参数:
* double * dpF - 指向频域值的指针
* double * dpf - 指向时域值的指针
* r －2的幂数
*
* 返回值:
* 无。
*
* 说明:
* 实现一维快速离散余弦反变换。
*
************************************************************************/
VOID WINAPI IDCT(double *dpF, double *dpf, int r)
{
double PI = 3.1415926;
// 离散余弦反变换点数
LONG lNum;
// 计算离散余弦变换点数
lNum = 1< *comX;
// 给复数变量分配内存
comX = new complex[lNum*2];
// 赋初值0
memset(comX, 0, sizeof(complex) * lNum * 2);
// 将频域变换后点写入数组comX
for(i=0;i (cos(i*PI/(lNum*2)) * dpF[i], sin(i*PI/(lNum*2) * dpF[i]));
}
// 调用快速付立叶反变换
IFFT_1D(comX,comX,r+1);
// 调整系数
dTemp = sqrt(2.0/lNum);
d0 = (sqrt(1.0/lNum) - dTemp) * dpF[0];
// 计算dpf(x)
for(i = 0; i < lNum; i++)
{
dpf[i] = d0 + 2 * lNum * comX[i].real()* dTemp ;
}
// 释放内存
delete comX;
}
/*************************************************************************
*
* 函数名称：
* WALSH()
*
* 参数:
* double * dpf - 指向时域值的指针
* double * dpF - 指向频域值的指针
* r －2的幂数
*
* 返回值:
* 无。
*
* 说明:
* 该函数用来实现一维快速沃尔什-哈达玛变换。
*
***********************************************************************
*/
VOID WINAPI WALSH(double *dpf, double *dpF, int r)
{
// 沃尔什-哈达玛变换点数
LONG lNum;
// 快速沃尔什变换点数
lNum = 1 << r;
// 循环变量
int i,j,k;
// 中间变量
int nTemp,m;
double *X1,*X2,*X;
// 分配运算所需的数组
X1 = new double[lNum];
X2 = new double[lNum];
// 将时域点写入数组X1
memcpy(X1, dpf, sizeof(double) * lNum);
for(k = 0; k < r; k++)
{
for(j = 0; j < 1<= 1; i--)
{
h = 0.0;
if (i > 1)
for (k = 0; k <= i-1; k++)
{
u = i*Rank + k;
h = h + QMatrix[u]*QMatrix[u];
}
// 如果一行全部为零
if (h + 1.0 == 1.0)
{
HypoCross[i] = 0.0;
if (i == 1)
HypoCross[i] = QMatrix[i*Rank+i-1];
MainCross[i] = 0.0;
}
// 否则进行householder变换，求正交矩阵的值
else
{
// 求次对角元素的值
HypoCross[i] = sqrt(h);
// 判断i行i-1列元素是不是大于零
u = i*Rank + i - 1;
if (QMatrix[u] > 0.0)
HypoCross[i] = -HypoCross[i];
h = h - QMatrix[u]*HypoCross[i];
QMatrix[u] = QMatrix[u] - HypoCross[i];
f = 0.0;
// householder变换
for (j = 0; j <= i-1; j++)
{
QMatrix[j*Rank+i] = QMatrix[i*Rank+j] / h;
g = 0.0;
for (k = 0; k <= j; k++)
g = g + QMatrix[j*Rank+k]*QMatrix[i*Rank+k];
if (j+1 <= i-1)
for (k = j+1; k <= i-1; k++)
g = g + QMatrix[k*Rank+j]*QMatrix[i*Rank+k];
HypoCross[j] = g / h;
f = f + g*QMatrix[j*Rank+i];
}
h2 = f / (h + h);
// 求正交矩阵的值
for (j = 0; j <= i-1; j++)
{
f = QMatrix[i*Rank + j];
g = HypoCross[j] - h2*f;
HypoCross[j] = g;
for (k = 0; k <= j; k++)
{
u = j*Rank + k;
QMatrix[u] = QMatrix[u] - f*HypoCross[k] - g*QMatrix[i*Rank + k];
}
}
MainCross[i] = h;
}
}
// 赋零值
for (i = 0; i <= Rank-2; i++)
HypoCross[i] = HypoCross[i + 1];
HypoCross[Rank - 1] = 0.0;
MainCross[0] = 0.0;
for (i = 0; i <= Rank-1; i++)
{
// 主对角元素的计算
if ((MainCross[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 + QMatrix[i*Rank + k]*QMatrix[k*Rank + j];
for (k = 0; k <= i-1; k++)
{
u = k*Rank + j;
QMatrix[u] = QMatrix[u] - g*QMatrix[k*Rank + i];
}
}
// 将主对角线的元素进行存储，以便返回
u = i*Rank + i;
MainCross[i] = QMatrix[u];
QMatrix[u] = 1.0;
// 将三对角外所有的元素赋零值
if (i-1 >= 0)
for (j = 0; j <= i-1; j++)
{
QMatrix[i*Rank + j] = 0.0;
QMatrix[j*Rank+i] = 0.0;
}
}
// 返回值
return(TRUE);
}
/*************************************************************************
*
* 函数名称：
* BSTQ()
*
* 参数:
* int Rank - 矩阵A的阶数
* double *MainCross - 对称三角阵中的主对角元素，返回时存放A的特征值
* double *HypoCross - 对称三角阵中的次对角元素
* double *Matrix - 返回对称矩阵A的特征向量
* double Eps - 控制精度
* int MaxT - 最大迭代次数
*
* 返回值:
* BOOL - 成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数用变形QR方法计算实对称三角矩阵的全部特征值以及相应的特征向量
*
*
***********************************************************************
*/
BOOL WINAPI BSTQ(int Rank, double *MainCross, double *HypoCross,
double *Matrix, double Eps, int MaxT)
{
// 变量的定义
int i, j, k, m, it, u, v;
double d, f, h, g, p, r, e, s;
// 赋零值
HypoCross[Rank - 1] = 0.0;
d = 0.0;
f = 0.0;
for(j = 0; j <= Rank-1; j++)
{
// 迭代精度的控制
it = 0;
h = Eps * (fabs(MainCross[j]) + fabs(HypoCross[j]));
if(h > d)
d = h;
m = j;
while((m <= Rank-1) && (fabs(HypoCross[m]) > d))
m = m + 1;
if(m != j)
{
// 进行迭代，求得矩阵A的特征值和特征向量
do
{
// 超过迭代次数，返回迭代失败
if(it == MaxT)
return(FALSE);
it = it + 1;
g = MainCross[j];
p = (MainCross[j + 1] - g) / (2.0 * HypoCross[j]);
r = sqrt(p*p + 1.0);
// 如果p大于0
if (p >= 0.0)
MainCross[j] = HypoCross[j]/(p + r);
else
MainCross[j] = HypoCross[j]/(p - r);
h = g - MainCross[j];
// 计算主对角线的元素
for (i = j + 1; i <= Rank - 1; i++)
MainCross[i] = MainCross[i] - h;
// 赋值
f = f + h;
p = MainCross[m];
e = 1.0; s = 0.0;
for(i = m - 1; i >= j; i--)
{
g = e * HypoCross[i];
h = e * p;
// 主对角线元素的绝对值是否大于次对角线元素的
if(fabs(p) >= fabs(HypoCross[i]))
{
e = HypoCross[i] / p;
r = sqrt(e*e + 1.0);
HypoCross[i + 1] = s*p*r;
s = e / r; e = 1.0 / r;
}
else
{
e = p / HypoCross[i];
r = sqrt(e*e + 1.0);
HypoCross[i+1] = s * HypoCross[i] * r;
s = 1.0 / r; e = e / r;
}
p = e*MainCross[i] - s*g;
MainCross[i + 1] = h + s*(e*g + s*MainCross[i]);
// 重新存储特征向量
for(k = 0; k <= Rank - 1; k++)
{
u = k*Rank + i + 1; v = u - 1;
h = Matrix[u];
Matrix[u] = s*Matrix[v] + e*h;
Matrix[v] = e*Matrix[v] - s*h;
}
}
// 将主对角线和次对角线元素重新赋值
HypoCross[j] = s * p;
MainCross[j] = e * p;
}
while (fabs(HypoCross[j]) > d);
}
MainCross[j] = MainCross[j] + f;
}
// 返回A的特征值
for (i = 0; i <= Rank-1; i++)
{
k = i; p = MainCross[i];
// 将A特征值赋给p
if(i+1 <= Rank-1)
{
j = i + 1;
while((j <= Rank-1) && (MainCross[j] <= p))
{ k = j;
p = MainCross[j];
j = j+1;
}
}
// 存储A的特征值和特征向量
if (k != i)
{
MainCross[k] = MainCross[i];
MainCross[i] = p;
for(j = 0; j <= Rank-1; j++)
{
u = j*Rank + i;
v = j*Rank + k;
p = Matrix[u];
Matrix[u] = Matrix[v];
Matrix[v] = p;
}
}
}
// 返回值
return(TRUE);
}
/*************************************************************************
*
* 函数名称：
* IWALSH()
*
* 参数:
* double * dpF - 指向频域值的指针
* double * dpf - 指向时域值的指针
* n －2的幂数
*
* 返回值:
* 无。
*
* 说明:
* 该函数用来实现一维快速沃尔什-哈达玛反变换。
*
***********************************************************************
*/
VOID WINAPI IWALSH(double *dpF, double *dpf, int n)
{
// 变换点数
LONG lNum;
// 循环变量
int i;
// 计算变换点数
lNum = 1 << n;
// 用快速沃尔什－哈达玛变换进行反变换
WALSH(dpF, dpf, n);
// 对系数进行调整
for(i = 0; i < lNum; i++)
{
dpf[i] *= lNum;
}
}
/*************************************************************************
*
* \函数名称：
* DFT_2D()
*
* \输入参数:
* CDib * pDib - 指向CDib类的指针，含有图像数据
* double * pTrRstRpart - 指向傅立叶系数实部的指针
* double * pTrRstIpart - 指向傅立叶系数虚部的指针
*
* \返回值:
* 无
*
* \说明:
* 二维傅立叶变换。
*
*************************************************************************
*/
VOID WINAPI DIBDFT_2D(CDib * pDib,double * pTrRstRpart, double * pTrRstIpart)
{
double PI = 3.14159;
//遍历图象的纵坐标
int y;
//遍历图象的横坐标
int x;
//频域的横坐标
int m;
//频域的纵坐标
int n;
//图象的长宽大小
CSize sizeImage = pDib->GetDimensions();
int nWidth = sizeImage.cx ;
int nHeight = sizeImage.cy ;
//图像在计算机在存储中的实际大小
CSize sizeImageSave = pDib->GetDibSaveDim();
int nSaveWidth = sizeImageSave.cx;
//图像数据的指针
LPBYTE pImageData = pDib->m_lpImage;
//初始化结果数据
for(n=0; nGetDimensions();
int nWidth = sizeImage.cx ;
int nHeight = sizeImage.cy ;
//图像在计算机在存储中的实际大小
CSize sizeImageSave = pDib->GetDibSaveDim();
int nSaveWidth = sizeImageSave.cx;
//图像数据的指针
LPBYTE pImageData = pDib->m_lpImage;
// 正弦和余弦值
double fCosTable;
double fSinTable;
fCosTable=0 ;
fSinTable=0 ;
// 临时变量
double fTmpPxValue;
double fRpartValue;
double fIpartValue;
fTmpPxValue=0;
fRpartValue=0;
fIpartValue=0;
for(y=0; yGetDimensions();
lWidth = SizeDim.cx;
lHeight = SizeDim.cy;
//得到实际的图象存储大小
CSize SizeRealDim;
SizeRealDim = pDib->GetDibSaveDim();
// 图像每行的字节数
LONG lLineBytes;
// 计算图像每行的字节数
lLineBytes = SizeRealDim.cx;
//图像数据的指针
LPBYTE lpDIBBits = pDib->m_lpImage;
// 保证离散余弦变换的宽度和高度为2的整数次方
while(lW * 2 <= lWidth)
{
lW = lW * 2;
wp++;
}
while(lH * 2 <= lHeight)
{
lH = lH * 2;
hp++;
}
// 分配内存
double *dpf = new double[lW * lH];
double *dpF = new double[lW * lH];
// 时域赋值
for(i = 0; i < lH; i++)
{
// 列
for(j = 0; j < lW; j++)
{
// 指向DIBi行j列象素的指针
lpSrc = lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
// 将象素值赋值给时域数组
dpf[j + i * lW] = *(lpSrc);
}
}
for(i = 0; i < lH; i++)
// 对y方向进行沃尔什-哈达玛变换
WALSH(dpf + lW * i, dpF + lW * i, wp);
// 保存计算结果
for(i = 0; i < lH; i++)
{
for(j = 0; j < lW; j++)
{
dpf[j * lH + i] = dpF[j + lW * i];
}
}
for(j = 0; j < lW; j++)
// 对x方向进行沃尔什-哈达玛变换
WALSH(dpf + j * lH, dpF + j * lH, hp);
// 行
for(i = 0; i < lH; i++)
{
// 列
for(j = 0; j < lW; j++)
{
// 计算频谱
dTemp = fabs(dpF[j * lH + i] * 1000);
if (dTemp > 255)
{
// 超过255直接设置为255
dTemp = 255;
}
// 指向DIBi行j列象素的指针
lpSrc = lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
}
//释放内存
delete dpf;
delete dpF;
// 返回
return TRUE;
}
/*************************************************************************
*
* 函数名称：
* DIBDct()
*
* 参数:
* CDib *pDib - 指向CDib类的指针
*
* 返回值:
* BOOL - 成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数用来对图像进行二维离散余弦变换。
*
***********************************************************************
*/
BOOL WINAPI DIBDct(CDib *pDib)
{
// 指向源图像的指针
unsigned char* lpSrc;
//图象的宽度和高度
LONG lWidth;
LONG lHeight;
// 循环变量
LONG i;
LONG j;
// 离散余弦变换的宽度和高度，必须为2的整数次方
LONG lW = 1;
LONG lH = 1;
int wp = 0;
int hp = 0;
// 中间变量
double dTemp;
//得到图象的宽度和高度
CSize SizeDim;
SizeDim = pDib->GetDimensions();
lWidth = SizeDim.cx;
lHeight = SizeDim.cy;
// 图像每行的字节数
LONG lLineBytes;
//得到实际的Dib图象存储大小
CSize SizeRealDim;
SizeRealDim = pDib->GetDibSaveDim();
// 计算图像每行的字节数
lLineBytes = SizeRealDim.cx;
//图像数据的指针
LPBYTE lpDIBBits = pDib->m_lpImage;
// 保证离散余弦变换的宽度和高度为2的整数次方
while(lW * 2 <= lWidth)
{
lW = lW * 2;
wp++;
}
while(lH * 2 <= lHeight)
{
lH = lH * 2;
hp++;
}
// 分配内存
double *dpf = new double[lW * lH];
double *dpF = new double[lW * lH];
// 时域赋值
for(i = 0; i < lH; i++)
{
for(j = 0; j < lW; j++)
{
// 指针指向位图i行j列的象素
lpSrc = lpDIBBits + lLineBytes * ( i) + j;
// lpSrc = lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
// 将象素值赋给时域数组
dpf[j + i * lW] = *(lpSrc);
}
}
for(i = 0; i < lH; i++)
// y方向进行离散余弦变换
DCT(&dpf[lW * i], &dpF[lW * i], wp);
// 保存计算结果
for(i = 0; i < lH; i++)
for(j = 0; j < lW; j++)
dpf[j * lH + i] = dpF[j + lW * i];
for(j = 0; j < lW; j++)
// x方向进行离散余弦变换
DCT(&dpf[j * lH], &dpF[j * lH], hp);
// 频谱的计算
for(i = 0; i < lH; i++)
{
for(j = 0; j < lW; j++)
{
dTemp = fabs(dpF[j*lH+i]);
// 是否超过255
if (dTemp > 255)
// 如果超过，设置为255
dTemp = 255;
// 指针指向位图i行j列的象素
lpSrc = lpDIBBits + lLineBytes * ( i) + j;
// lpSrc = lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
}
// 释放内存
delete dpf;
delete dpF;
// 返回
return TRUE;
}
/*************************************************************************
*
* 函数名称：
* DIBOHOTELLING()
*
* 参数:
* CDib *pDib - 指向CDib类的指针
*
* 返回值:
* BOOL - 成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数用霍特林变换来对图像进行旋转。
*
***********************************************************************
*/
BOOL WINAPI DIBHOTELLING(CDib *pDib)
{
// 指向源图像的指针
unsigned char* lpSrc;
// 循环变量
LONG i;
LONG j;
//图象的宽度和高度
LONG lWidth;
LONG lHeight;
// 图像每行的字节数
LONG lLineBytes;
// 经过变换后图象最大可能范围
LONG lMaxRange;
// 物体坐标的均值
LONG AverEx;
LONG AverEy;
// 物体总的象素数
LONG ToaCount;
// 坐标值的协方差矩阵
double Matr4C[2][2];
// 存放协方差矩阵的特征向量
double QMatrix[2][2];
// 三对角阵的主对角和次对角线元素
double MainCross[2];
double HypoCross[2];
// 中间变量
double dTemp;
LONG lTempI;
LONG lTempJ;
//得到图象的宽度和高度
CSize SizeDim;
SizeDim = pDib->GetDimensions();
lWidth = SizeDim.cx;
lHeight = SizeDim.cy;
//得到实际的Dib图象存储大小
CSize SizeRealDim;
SizeRealDim = pDib->GetDibSaveDim();
// 计算图像每行的字节数
lLineBytes = SizeRealDim.cx;
//图像数据的指针
LPBYTE lpDIBBits = pDib->m_lpImage;
// 估计图象经过旋转后可能最大的宽度和高度
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++)
{
// 给旋转后坐标轴的每个点赋零值（灰度值255对应显示白色）
F[i*lWidth + j] = 255;
// 指向位图i行j列象素的指针
lpSrc = lpDIBBits + lLineBytes*i + j;
// 值小于255（非背景色白色）的象素认为物体的一部分
// 并将其坐标值x和y看作二维随机矢量
if((*lpSrc) < 255)
{
// 属于物体象素的Y坐标和X坐标累计值
AverEx=AverEx+i;
AverEy=AverEy+j;
// 物体总的象素数加一
ToaCount++;
// 随机矢量协方差矩阵的累计值
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 Eps = 0.000001;
// 将协方差矩阵化作三对角对称阵
THREECROSS(*Matr4C, 2, *QMatrix, MainCross, HypoCross);
// 求协方差矩阵的特征值和特征矢向量
BSTQ(2, MainCross,HypoCross, *Matr4C, Eps, 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;
}
// 查找经霍特林变换后的坐标点在原坐标系中的坐标
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=0 && Cy=0 && lTempI=0 && lTempJ