哈尔小波变换的原理及其实现（Haar）另外参见俄罗斯写的/Articles/22243/Real-Time-Object-Tracker-in-CHaar小波在图像处理和数字水印等方面应用较多，这里简单的介绍一下哈尔小波的基本原理以及其实现情况。
一、Haar小波的基本原理数学理论方面的东西我也不是很熟悉，这边主要用简单的例子来介绍下Haar小波的使用情况。
例如：有a=[8,7,6,9]四个数，并使用b[4]数组来保存结果.则一级Haar小波变换的结果为:b[0]=(a[0]+a[1])/2, b[2]=(a[0]-a[1])/2b[1]=(a[2]+a[3])/2, b[3]=(a[2]-a[3])/2即依次从数组中取两个数字，计算它们的和以及差，并将和一半和差的一半依次保存在数组的前半部分和后半部分。
例如：有a[8]，要进行一维Haar小波变换，结果保存在b[8]中则一级Haar小波变换的结果为:b[0]=(a[0]+a[1])/2, b[4]=(a[0]-a[1])/2b[1]=(a[2]+a[3])/2, b[5]=(a[2]-a[3])/2b[2]=(a[4]+a[5])/2, b[6]=(a[4-a[5]])/2b[3]=(a[6]+a[7])/2, b[7]=(a[6]-a[7])/2如果需要进行二级Haar小波变换的时候，只需要对b[0]-b[3]进行Haar小波变换.对于二维的矩阵来讲，每一级Haar小波变换需要先后进行水平方向和竖直方向上的两次一维小波变换,行和列的先后次序对结果不影响。
二、Haar小波的实现使用opencv来读取图片及像素，对图像的第一个8*8的矩阵做了一级小波变换#include <cv.h>#include <highgui.h>#include <iostream>using namespace std;int main(){IplImage* srcImg;double imgData[8][8];int i,j;srcImg=cvLoadImage("lena.bmp",0);cout<<"原8*8数据"<<endl;for( i=0;i<8;i++){for( j=0;j<8;j++){imgData[i][j]=cvGetReal2D(srcImg,i+256,j+16); cout<<imgData[i][j]<<" ";}cout<<endl;}double tempData[8];//行小波分解for( i=0;i<8;i++){for( j=0;j<4;j++){double temp1=imgData[i][2*j];double temp2=imgData[i][2*j+1]; tempData[j]=(temp1+temp2)/2;tempData[j+4]=(temp1-temp2)/2;}for( j=0;j<8;j++)imgData[i][j]=tempData[j];}//列小波分解for( i=0;i<8;i++){for( j=0;j<4;j++){double temp1=imgData[2*j][i]; double temp2=imgData[2*j+1][i]; tempData[j]=(temp1+temp2)/2; tempData[j+4]=(temp1-temp2)/2; }for( j=0;j<8;j++)imgData[j][i]=tempData[j];}cout<<"1级小波分解数据"<<endl; for( i=0;i<8;i++){for( j=0;j<8;j++){cout<<imgData[i][j]<<" ";}cout<<endl;}//列小波逆分解for( i=0;i<8;i++){for( j=0;j<4;j++){double temp1=imgData[j][i]; double temp2=imgData[j+4][i]; tempData[2*j]=temp1+temp2; tempData[2*j+1]=temp1-temp2; }for( j=0;j<8;j++){imgData[j][i]=tempData[j];}}//行小波逆分解for( i=0;i<8;i++){for( j=0;j<4;j++){double temp1=imgData[i][j]; double temp2=imgData[i][j+4]; tempData[2*j]=temp1+temp2; tempData[2*j+1]=temp1-temp2; }for( j=0;j<2*4;j++){imgData[i][j]=tempData[j];}}cout<<"1级小波逆分解数据"<<endl; for( i=0;i<8;i++){for( j=0;j<8;j++){cout<<imgData[i][j]<<" ";}cout<<endl;}return 0;}====================================== ===================================== /// 小波变换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 );#endifrow *= 2;col *= 2;src = dst;}return dst;}/////////////////////////////////////////////////////////////////// //////////////////////////// 调用函数/// 生成不同类型的小波，现在只有haar，sym2void 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;}====================================== =====================================*************************************************************** **************************************************其他代码实现*************************************************************** **************************************************// 哈尔小波.cpp : 定义控制台应用程序的入口点。
