小波变换 opencv

<pre name="code" class="cpp">///  小波变换
Mat WDT( const Mat &_src, const string _wname, const int _level )const
{    
	int reValue = THID_ERR_NONE;  
	Mat src = Mat_<float>(_src);   
	Mat dst = Mat::zeros( src.rows, src.cols, src.type() );   
	int N = src.rows;  
	int D = src.cols;   
	/// 高通低通滤波器  
	Mat lowFilter;    
	Mat highFilter;   
	wavelet( _wname, lowFilter, highFilter );   
	/// 小波变换   
	int t=1;   
	int row = N;   
	int col = D;     
	while( t<=_level ) 
	{    
		///先进行行小波变换      
		for( int i=0; i<row; i++ )      
		{          
			/// 取出src中要处理的数据的一行   
			Mat oneRow = Mat::zeros( 1,col, src.type() );      
			for ( int j=0; j<col; j++ )      
			{            
				oneRow.at<float>(0,j) = src.at<float>(i,j);        
			}         
			oneRow = waveletDecompose( oneRow, lowFilter, highFilter );        
			/// 将src这一行置为oneRow中的数据       
			for ( int j=0; j<col; j++ )     
			{             
				dst.at<float>(i,j) = oneRow.at<float>(0,j);    
			}      
		}
#if 0     
		//normalize( dst, dst, 0, 255, NORM_MINMAX );     
		IplImage dstImg1 = IplImage(dst);      
		cvSaveImage( "dst.jpg", &dstImg1 );
#endif    
			/// 小波列变换      
			for ( int j=0; j<col; j++ )   
			{          
				/// 取出src数据的一行输入     
				Mat oneCol = Mat::zeros( row, 1, src.type() );   
				for ( int i=0; i<row; i++ )        
				{             
					oneCol.at<float>(i,0) = dst.at<float>(i,j);   
				}          
				oneCol = ( waveletDecompose( oneCol.t(), lowFilter, highFilter ) ).t();         
				for ( int i=0; i<row; i++ )        
				{               
					dst.at<float>(i,j) = oneCol.at<float>(i,0);     
				}      
			} 
#if 0       
			//normalize( dst, dst, 0, 255, NORM_MINMAX );     
			IplImage dstImg2 = IplImage(dst);      
			cvSaveImage( "dst.jpg", &dstImg2 );
#endif       
			/// 更新    
			row /= 2;     
			col /=2;      
			t++;       
			src = dst;  
	}   
	return dst;
} 
///  小波逆变换
Mat IWDT( const Mat &_src, const string _wname, const int _level )const
{    
	int reValue = THID_ERR_NONE; 
	Mat src = Mat_<float>(_src);  
	Mat dst = Mat::zeros( src.rows, src.cols, src.type() );  
	int N = src.rows;    int D = src.cols;   
	/// 高通低通滤波器 
	Mat lowFilter;    
	Mat highFilter;   
	wavelet( _wname, lowFilter, highFilter );  
	/// 小波变换
	int t=1;   
	int row = N/std::pow( 2., _level-1);
	int col = D/std::pow(2., _level-1); 
	while ( row<=N && col<=D )   
	{       
		/// 小波列逆变换    
		for ( int j=0; j<col; j++ )    
		{         
			/// 取出src数据的一行输入       
			Mat oneCol = Mat::zeros( row, 1, src.type() );    
			for ( int i=0; i<row; i++ )      
			{           
				oneCol.at<float>(i,0) = src.at<float>(i,j);    
			}         
			oneCol = ( waveletReconstruct( oneCol.t(), lowFilter, highFilter ) ).t();      
			for ( int i=0; i<row; i++ )        
			{            
				dst.at<float>(i,j) = oneCol.at<float>(i,0);          
			}      
		} 
#if 0       
		//normalize( dst, dst, 0, 255, NORM_MINMAX );   
		IplImage dstImg2 = IplImage(dst);    
		cvSaveImage( "dst.jpg", &dstImg2 );
#endif      
		///行小波逆变换      
		for( int i=0; i<row; i++ )    
		{          
			/// 取出src中要处理的数据的一行  
			Mat oneRow = Mat::zeros( 1,col, src.type() );            
			for ( int j=0; j<col; j++ )         
			{              
				oneRow.at<float>(0,j) = dst.at<float>(i,j);     
			}         
			oneRow = waveletReconstruct( oneRow, lowFilter, highFilter );        
			/// 将src这一行置为oneRow中的数据          
			for ( int j=0; j<col; j++ )         
			{              
				dst.at<float>(i,j) = oneRow.at<float>(0,j);     
			}      
		}
#if 0      
		//normalize( dst, dst, 0, 255, NORM_MINMAX );  
		IplImage dstImg1 = IplImage(dst);    
		cvSaveImage( "dst.jpg", &dstImg1 );
#endif      
			row *= 2;    
			col *= 2;      
			src = dst;   
	}    
	return dst;
}  
 
/// 调用函数
/// 生成不同类型的小波,现在只有haar,sym2
void wavelet( const string _wname, Mat &_lowFilter, Mat &_highFilter )const
{   
	if ( _wname=="haar" || _wname=="db1" )  
	{       
		int N = 2;    
		_lowFilter = Mat::zeros( 1, N, CV_32F );    

		_highFilter = Mat::zeros( 1, N, CV_32F );             
		_lowFilter.at<float>(0, 0) = 1/sqrtf(N);    
		_lowFilter.at<float>(0, 1) = 1/sqrtf(N);        
		_highFilter.at<float>(0, 0) = -1/sqrtf(N);        
		_highFilter.at<float>(0, 1) = 1/sqrtf(N);   
	}   
	if ( _wname =="sym2" )  
	{       
		int N = 4;     
		float h[] = {-0.483, 0.836, -0.224, -0.129 };  
		float l[] = {-0.129, 0.224,    0.837, 0.483 };  
		_lowFilter = Mat::zeros( 1, N, CV_32F );     
		_highFilter = Mat::zeros( 1, N, CV_32F );    
		for ( int i=0; i<N; i++ )      
		{         
			_lowFilter.at<float>(0, i) = l[i];     
			_highFilter.at<float>(0, i) = h[i];   
		}    
	}
} 
/// 小波分解
Mat waveletDecompose( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter )const
{   
	assert( _src.rows==1 && _lowFilter.rows==1 && _highFilter.rows==1 );    
	assert( _src.cols>=_lowFilter.cols && _src.cols>=_highFilter.cols );  
	Mat &src = Mat_<float>(_src);   
	int D = src.cols;   
	Mat &lowFilter = Mat_<float>(_lowFilter); 
	Mat &highFilter = Mat_<float>(_highFilter);   
	/// 频域滤波,或时域卷积;ifft( fft(x) * fft(filter)) = cov(x,filter)    
	Mat dst1 = Mat::zeros( 1, D, src.type() );  
	Mat dst2 = Mat::zeros( 1, D, src.type()  );  
	filter2D( src, dst1, -1, lowFilter ); 
	filter2D( src, dst2, -1, highFilter );     
	/// 下采样    
	Mat downDst1 = Mat::zeros( 1, D/2, src.type() ); 
	Mat downDst2 = Mat::zeros( 1, D/2, src.type() );   
	resize( dst1, downDst1, downDst1.size() );  
	resize( dst2, downDst2, downDst2.size() );     
	/// 数据拼接   
	for ( int i=0; i<D/2; i++ )  
	{     
		src.at<float>(0, i) = downDst1.at<float>( 0, i );    
		src.at<float>(0, i+D/2) = downDst2.at<float>( 0, i );   
	}    
	return src;
}
/// 小波重建
Mat waveletReconstruct( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter )const
{    
	assert( _src.rows==1 && _lowFilter.rows==1 && _highFilter.rows==1 );  
	assert( _src.cols>=_lowFilter.cols && _src.cols>=_highFilter.cols ); 
	Mat &src = Mat_<float>(_src);   
	int D = src.cols;  
	Mat &lowFilter = Mat_<float>(_lowFilter); 
	Mat &highFilter = Mat_<float>(_highFilter);    
	/// 插值;  
	Mat Up1 = Mat::zeros( 1, D, src.type() );  
	Mat Up2 = Mat::zeros( 1, D, src.type() );  
	/// 插值为0 
	//for ( int i=0, cnt=1; i<D/2; i++,cnt+=2 ) 
	//{ 
	//    Up1.at<float>( 0, cnt ) = src.at<float>( 0, i );   ///< 前一半   
	//    Up2.at<float>( 0, cnt ) = src.at<float>( 0, i+D/2 ); ///< 后一半  
	//}   
	/// 线性插值 
	Mat roi1( src, Rect(0, 0, D/2, 1) );   
	Mat roi2( src, Rect(D/2, 0, D/2, 1) );   
	resize( roi1, Up1, Up1.size(), 0, 0, INTER_CUBIC ); 
	resize( roi2, Up2, Up2.size(), 0, 0, INTER_CUBIC ); 
	/// 前一半低通,后一半高通  
	Mat dst1 = Mat::zeros( 1, D, src.type() );   
	Mat dst2= Mat::zeros( 1, D, src.type() );   
	filter2D( Up1, dst1, -1, lowFilter );
	filter2D( Up2, dst2, -1, highFilter );    
	/// 结果相加   
	dst1 = dst1 + dst2;   
	return dst1; 
} 


 









#include <opencv2\highgui\highgui.hpp>
#include <math.h>

#include <iostream>
using namespace std;
//二维离散变换小波源代码
void DWT(IplImage *pImage, int nLayer)
{
   // 执行条件
   if (pImage)
   {
      if (pImage->nChannels == 1 &&
         pImage->depth == IPL_DEPTH_32F &&
         ((pImage->width >> nLayer) << nLayer) == pImage->width &&
         ((pImage->height >> nLayer) << nLayer) == pImage->height)
      {
         int     i, x, y, n;
         float   fValue   = 0;
         float   fRadius  = sqrt(2.0f);
         int     nWidth   = pImage->width;
         int     nHeight  = pImage->height;
         int     nHalfW   = nWidth / 2;
         int     nHalfH   = nHeight / 2;
         float **pData    = new float*[pImage->height];
         float  *pRow     = new float[pImage->width];
         float  *pColumn  = new float[pImage->height];
         for (i = 0; i < pImage->height; i++)   pData[i] = (float*) (pImage->imageData + pImage->widthStep * i);
         // 多层小波变换
         for (n = 0; n < nLayer; n++, nWidth /= 2, nHeight /= 2, nHalfW /= 2, nHalfH /= 2)
         {
            // 水平变换
            for (y = 0; y < nHeight; y++)
            {
               // 奇偶分离
               memcpy(pRow, pData[y], sizeof(float) * nWidth);
               for (i = 0; i < nHalfW; i++)
               {
                  x = i * 2;
                  pData[y][i] = pRow[x];
                  pData[y][nHalfW + i] = pRow[x + 1];
               }
               // 提升小波变换
               for (i = 0; i < nHalfW - 1; i++)
               {
                  fValue = (pData[y][i] + pData[y][i + 1]) / 2;
                  pData[y][nHalfW + i] -= fValue;
               }
               fValue = (pData[y][nHalfW - 1] + pData[y][nHalfW - 2]) / 2;
               pData[y][nWidth - 1] -= fValue;
               fValue = (pData[y][nHalfW] + pData[y][nHalfW + 1]) / 4;
               pData[y][0] += fValue;
               for (i = 1; i < nHalfW; i++)
               {
                  fValue = (pData[y][nHalfW + i] + pData[y][nHalfW + i - 1]) / 4;
                  pData[y][i] += fValue;
               }
               // 频带系数
               for (i = 0; i < nHalfW; i++)
               {
                  pData[y][i] *= fRadius;
                  pData[y][nHalfW + i] /= fRadius;
               }
            }
            // 垂直变换
            for (x = 0; x < nWidth; x++)
            {
               // 奇偶分离
               for (i = 0; i < nHalfH; i++)
               {
                  y = i * 2;
                  pColumn[i] = pData[y][x];
                  pColumn[nHalfH + i] = pData[y + 1][x];
               }
               for (i = 0; i < nHeight; i++)
               {
                  pData[i][x] = pColumn[i];
               }
               // 提升小波变换
               for (i = 0; i < nHalfH - 1; i++)
               {
                  fValue = (pData[i][x] + pData[i + 1][x]) / 2;
                  pData[nHalfH + i][x] -= fValue;
               }
               fValue = (pData[nHalfH - 1][x] + pData[nHalfH - 2][x]) / 2;
               pData[nHeight - 1][x] -= fValue;
               fValue = (pData[nHalfH][x] + pData[nHalfH + 1][x]) / 4;
               pData[0][x] += fValue;
               for (i = 1; i < nHalfH; i++)
               {
                  fValue = (pData[nHalfH + i][x] + pData[nHalfH + i - 1][x]) / 4;
                  pData[i][x] += fValue;
               }
               // 频带系数
               for (i = 0; i < nHalfH; i++)
               {
                  pData[i][x] *= fRadius;
                  pData[nHalfH + i][x] /= fRadius;
               }
            }
         }
         delete[] pData;
         delete[] pRow;
         delete[] pColumn;
      }
   }
}

// 二维离散小波恢复(单通道浮点图像)
void IDWT(IplImage *pImage, int nLayer)
{
   // 执行条件
   if (pImage)
   {
      if (pImage->nChannels == 1 &&
         pImage->depth == IPL_DEPTH_32F &&
         ((pImage->width >> nLayer) << nLayer) == pImage->width &&
         ((pImage->height >> nLayer) << nLayer) == pImage->height)
      {
         int     i, x, y, n;
         float   fValue   = 0;
         float   fRadius  = sqrt(2.0f);
         int     nWidth   = pImage->width >> (nLayer - 1);
         int     nHeight  = pImage->height >> (nLayer - 1);
         int     nHalfW   = nWidth / 2;
         int     nHalfH   = nHeight / 2;
         float **pData    = new float*[pImage->height];
         float  *pRow     = new float[pImage->width];
         float  *pColumn  = new float[pImage->height];
         for (i = 0; i < pImage->height; i++)
         {
            pData[i] = (float*) (pImage->imageData + pImage->widthStep * i);
         }
         // 多层小波恢复
         for (n = 0; n < nLayer; n++, nWidth *= 2, nHeight *= 2, nHalfW *= 2, nHalfH *= 2)
         {
            // 垂直恢复
            for (x = 0; x < nWidth; x++)
            {
               // 频带系数
               for (i = 0; i < nHalfH; i++)
               {
                  pData[i][x] /= fRadius;
                  pData[nHalfH + i][x] *= fRadius;
               }
               // 提升小波恢复
               fValue = (pData[nHalfH][x] + pData[nHalfH + 1][x]) / 4;
               pData[0][x] -= fValue;
               for (i = 1; i < nHalfH; i++)
               {
                  fValue = (pData[nHalfH + i][x] + pData[nHalfH + i - 1][x]) / 4;
                  pData[i][x] -= fValue;
               }
               for (i = 0; i < nHalfH - 1; i++)
               {
                  fValue = (pData[i][x] + pData[i + 1][x]) / 2;
                  pData[nHalfH + i][x] += fValue;
               }
               fValue = (pData[nHalfH - 1][x] + pData[nHalfH - 2][x]) / 2;
               pData[nHeight - 1][x] += fValue;
               // 奇偶合并
               for (i = 0; i < nHalfH; i++)
               {
                  y = i * 2;
                  pColumn[y] = pData[i][x];
                  pColumn[y + 1] = pData[nHalfH + i][x];
               }
               for (i = 0; i < nHeight; i++)   pData[i][x] = pColumn[i];
            }
            // 水平恢复
            for (y = 0; y < nHeight; y++)
            {
               // 频带系数
               for (i = 0; i < nHalfW; i++)
               {
                  pData[y][i] /= fRadius;
                  pData[y][nHalfW + i] *= fRadius;
               }
               // 提升小波恢复
               fValue = (pData[y][nHalfW] + pData[y][nHalfW + 1]) / 4;
               pData[y][0] -= fValue;
               for (i = 1; i < nHalfW; i++)
               {
                  fValue = (pData[y][nHalfW + i] + pData[y][nHalfW + i - 1]) / 4;
                  pData[y][i] -= fValue;
               }
               for (i = 0; i < nHalfW - 1; i++)
               {
                  fValue = (pData[y][i] + pData[y][i + 1]) / 2;
                  pData[y][nHalfW + i] += fValue;
               }
               fValue = (pData[y][nHalfW - 1] + pData[y][nHalfW - 2]) / 2;
               pData[y][nWidth - 1] += fValue;
               // 奇偶合并
               for (i = 0; i < nHalfW; i++)
               {
                  x = i * 2;
                  pRow[x] = pData[y][i];
                  pRow[x + 1] = pData[y][nHalfW + i];
               }
               memcpy(pData[y], pRow, sizeof(float) * nWidth);
            }
         }
         delete[] pData;
         delete[] pRow;
         delete[] pColumn;
      }
   }
} 


int main()
{
	// 小波变换层数
	int nLayer = 2;
	// 输入彩色图像
	IplImage *pSrc = cvLoadImage("lena.jpg", CV_LOAD_IMAGE_COLOR);
	// 计算小波图象大小
	CvSize size = cvGetSize(pSrc);
	if ((pSrc->width >> nLayer) << nLayer != pSrc->width)	size.width = ((pSrc->width >> nLayer) + 1) << nLayer;
	if ((pSrc->height >> nLayer) << nLayer != pSrc->height)  size.height = ((pSrc->height >> nLayer) + 1) << nLayer;
	// 创建小波图象
	IplImage *pWavelet = cvCreateImage(size, IPL_DEPTH_32F, pSrc->nChannels);
	if (pWavelet) {
		// 小波图象赋值
		cvSetImageROI(pWavelet, cvRect(0, 0, pSrc->width, pSrc->height));
		cvConvertScale(pSrc, pWavelet, 1, -128);
		cvResetImageROI(pWavelet);
		// 彩色图像小波变换
		IplImage *pImage = cvCreateImage(cvGetSize(pWavelet), IPL_DEPTH_32F, 1);
		if (pImage) {
			for (int i = 1; i <= pWavelet->nChannels; i++)   {
				cvSetImageCOI(pWavelet, i);
				cvCopy(pWavelet, pImage, NULL);
				// 二维离散小波变换
				DWT(pImage, nLayer);
				// 二维离散小波恢复
				// IDWT(pImage, nLayer);
				cvCopy(pImage, pWavelet, NULL);
			}
			cvSetImageCOI(pWavelet, 0);
			cvReleaseImage(&pImage);
		}
		// 小波变换图象
		cvSetImageROI(pWavelet, cvRect(0, 0, pSrc->width, pSrc->height));
		cvConvertScale(pWavelet, pSrc, 1, 128);
		cvResetImageROI(pWavelet);
		cvReleaseImage(&pWavelet);
	}
	// 显示图像pSrc
	cv::namedWindow("result");
	cvShowImage("result",pSrc);
	cv::waitKey(0);
	cvReleaseImage(&pSrc); 
}


  • 0
    点赞
  • 9
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值