fretrans.cpp

含《Visual c++数字图像处理典型算法及实现》这本书里所有源代码

 *   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<<k; j++)
		{	
			// 按照蝶形运算图进行运算
			nTemp = 1 << (r-k);
			for(i = 0; i < nTemp / 2; i++)
			{
				m = j * nTemp;
				X2[i + m] = X1[i + m] + X1[i + m + nTemp / 2];
				X2[i + m + nTemp / 2] = X1[i + m] - X1[i + m + nTemp / 2];
			}
		}
		
		// 互换  
		X = X2;
		X2 = X1;
		X1 = X;
	}
	
	// 对系数做调整
	for(j = 0; j < lNum; j++)
	{
		m = 0;
		for(i = 0; i < r; i++)
		{
			if (j & (1<<i))
			{
				m += 1 << (r-i-1);
			}
		}

		dpF[j] = X1[m] / lNum;
	}
	
	// 释放内存
	delete X1;
	delete X2;
}


/*************************************************************************
 *
 * 函数名称：
 *   THREECROSS()
 *
 * 参数:
 *   double  *Matrix     - 用来存放矩阵A
 *   int     Rank        - 矩阵A的阶数
 *   double	 *QMatrix    - 返回householder变换的矩阵Q
 *   double  *MainCross  - 对称三角阵中的主对角元素
 *   double  *HypoCross  - 对称三角阵中的次对角元素
 *
 * 返回值:
 *   BOOL                - 成功返回TRUE，否则返回FALSE。
 *
 * 说明:
 *   该函数用householder变换将n阶实对称矩阵化为对称三角阵
 *
 *
 ***********************************************************************
*/

BOOL WINAPI THREECROSS(double *Matrix, int Rank, double *QMatrix, 
					    double *MainCross, double *HypoCross)
{
	//   循环变量
	int i, j, k, u;
    double h, f, g, h2;
    
	// 对矩阵QMatrix赋值
	for(i = 0; i <= Rank-1; i++)
		for(j = 0; j <= Rank-1; j++)
		{
			u = i*Rank + j; 
			QMatrix[u] = Matrix[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 + 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; n<nHeight ; n++ )
		for(m=0 ; m<nWidth ; m++ )
		{
			*(	pTrRstRpart + n*nWidth + m	) =0;
			*(	pTrRstIpart + n*nWidth + m	) =0;
		}
	double fCosTable;
	double fSinTable;
	int	  nPxValue;

	fCosTable=0 ;
	nPxValue =0;

	double fTmpRstR;
	double fTmpRstI;
	for(n=0; n<nHeight ; n++ )
		for(m=0 ; m<nWidth ; m++ )
		{
			fTmpRstR=0;
			fTmpRstI=0;
			for(y=0; y<nHeight ; y++ )
				for(x=0 ; x<nWidth ; x++ )
				{
					fCosTable= 
						cos(	2*PI*( ((double)m*x)/nWidth + ((double)n*y)/nHeight) ) ;
					fSinTable= 
						sin(	-2*PI*( ((double)m*x)/nWidth + ((double)n*y)/nHeight) ) ;
					nPxValue = *(pImageData+ y*nSaveWidth + x )			;

					fTmpRstR+=fCosTable* nPxValue						;
					fTmpRstI+=fSinTable* nPxValue						;
				}
			*( pTrRstRpart + nWidth * n + m ) = fTmpRstR;
			*( pTrRstIpart + nWidth * n + m ) = fTmpRstI;
		}
}

/*************************************************************************
 *
 * \函数名称：
 *   IDFT_2D()
 *
 * \输入参数:
 *   CDib * pDib				- 指向CDib类的指针，含有图像数据
 *   double * pTrRstRpart		- 指向傅立叶系数实部的指针
 *   double * pTrRstIpart		- 指向傅立叶系数虚部的指针
 *
 * \返回值:
 *   无
 *
 * \说明:
 *   二维傅立叶反变换。
 *
 *************************************************************************
 */
VOID WINAPI IDFT_2D(CDib * pDib,double * pTrRstRpart, double * pTrRstIpart)
{
	double PI = 3.1415926;
	//遍历图象的纵坐标
	int y;

	//遍历图象的横坐标
	int x;

	//频域的横坐标
	int m;

	//频域的纵坐标
	int n;