|
另外參見俄羅斯寫的http://www./Articles/22243/Real-Time-Object-Tracker-in-C Haar小波在圖像處理和數(shù)字水印等方面應(yīng)用較多,，這里簡單的介紹一下哈爾小波的基本原理以及其實(shí)現(xiàn)情況。 一、Haar小波的基本原理 數(shù)學(xué)理論方面的東西我也不是很熟悉,，這邊主要用簡單的例子來介紹下Haar小波的使用情況。 例如：有a=[8,7,6,9]四個(gè)數(shù),，并使用b[4]數(shù)組來保存結(jié)果. 則一級(jí)Haar小波變換的結(jié)果為: b[0]=(a[0]+a[1])/2, b[2]=(a[0]-a[1])/2 b[1]=(a[2]+a[3])/2, b[3]=(a[2]-a[3])/2 即依次從數(shù)組中取兩個(gè)數(shù)字，計(jì)算它們的和以及差,，并將和一半和差的一半依次保存在數(shù)組的前半部分和后半部分,。 例如：有a[8]，要進(jìn)行一維Haar小波變換,，結(jié)果保存在b[8]中 則一級(jí)Haar小波變換的結(jié)果為: b[0]=(a[0]+a[1])/2, b[4]=(a[0]-a[1])/2 b[1]=(a[2]+a[3])/2, b[5]=(a[2]-a[3])/2 b[2]=(a[4]+a[5])/2, b[6]=(a[4-a[5]])/2 b[3]=(a[6]+a[7])/2, b[7]=(a[6]-a[7])/2 如果需要進(jìn)行二級(jí)Haar小波變換的時(shí)候,，只需要對b[0]-b[3]進(jìn)行Haar小波變換. 對于二維的矩陣來講，每一級(jí)Haar小波變換需要先后進(jìn)行水平方向和豎直方向上的兩次一維小波變換,行和列的先后次序?qū)Y(jié)果不影響,。 二,、Haar小波的實(shí)現(xiàn) 使用opencv來讀取圖片及像素，對圖像的第一個(gè)8*8的矩陣做了一級(jí)小波變換 #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數(shù)據(jù)"<<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級(jí)小波分解數(shù)據(jù)"<<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級(jí)小波逆分解數(shù)據(jù)"<<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 ) {///先進(jìn)行行小波變換 for( int i=0; i<row; i++ ) {/// 取出src中要處理的數(shù)據(jù)的一行 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中的數(shù)據(jù) 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數(shù)據(jù)的一行輸入 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數(shù)據(jù)的一行輸入 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中要處理的數(shù)據(jù)的一行 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中的數(shù)據(jù) 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; } //////////////////////////////////////////////////////////////////////////////////////////// /// 調(diào)用函數(shù) /// 生成不同類型的小波,，現(xiàn)在只有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); /// 頻域?yàn)V波，或時(shí)域卷積,；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() ); /// 數(shù)據(jù)拼接 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 ); /// 結(jié)果相加 dst1 = dst1 + dst2; return dst1; } =========================================================================== ***************************************************************************************************************** 其他代碼實(shí)現(xiàn) ***************************************************************************************************************** // 哈爾小波.cpp : 定義控制臺(tái)應(yīng)用程序的入口點(diǎn)。 // #include "stdafx.h" #include <iostream> #include "cv.h" #include "highgui.h" #include <ctype.h> // 二維離散小波變換（單通道浮點(diǎn)圖像） void DWT(IplImage *pImage, int nLayer) { // 執(zhí)行條件 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; } // 頻帶系數(shù) 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; } // 頻帶系數(shù) for (i = 0; i < nHalfH; i++) { pData[i][x] *= fRadius; pData[nHalfH + i][x] /= fRadius; } } } delete[] pData; delete[] pRow; delete[] pColumn; } } } // 二維離散小波恢復(fù)（單通道浮點(diǎn)圖像） void IDWT(IplImage *pImage, int nLayer) { // 執(zhí)行條件 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); } // 多層小波恢復(fù) for (n = 0; n < nLayer; n++, nWidth *= 2, nHeight *= 2, nHalfW *= 2, nHalfH *= 2) { // 垂直恢復(fù) for (x = 0; x < nWidth; x++) { // 頻帶系數(shù) for (i = 0; i < nHalfH; i++) { pData[i][x] /= fRadius; pData[nHalfH + i][x] *= fRadius; } // 提升小波恢復(fù) 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]; } } // 水平恢復(fù) for (y = 0; y < nHeight; y++) { // 頻帶系數(shù) for (i = 0; i < nHalfW; i++) { pData[y][i] /= fRadius; pData[y][nHalfW + i] *= fRadius; } // 提升小波恢復(fù) 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 _tmain(int argc, _TCHAR* argv[]) { // 小波變換層數(shù) int nLayer = 1; // 輸入彩色圖像 IplImage *pSrc = cvLoadImage("lena.bmp", CV_LOAD_IMAGE_COLOR); // 計(jì)算小波圖象大小 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; } // 創(chuàng)建小波圖象 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); // 二維離散小波恢復(fù) // 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); cvNamedWindow("pWavelet",CV_WINDOW_AUTOSIZE); cvShowImage("pWavelet",pWavelet); cvNamedWindow("pSrc",CV_WINDOW_AUTOSIZE); cvShowImage("pSrc",pSrc); cvWaitKey(0); cvResetImageROI(pWavelet); // 本行代碼有點(diǎn)多余,，但有利用養(yǎng)成良好的編程習(xí)慣 cvReleaseImage(&pWavelet); } // 顯示圖像pSrc // ... cvReleaseImage(&pSrc); return 0; } ************************ -------------------------------------Codeproject Haar ------------------------------------------------------------------------------------------------- void Haar::transrows(char** dest, char** sour, unsigned int w, unsigned int h) const { unsigned int w2 = w / 2; __m64 m00FF; m00FF.m64_u64 = 0x00FF00FF00FF00FF; for (unsigned int y = 0; y < h; y++) { __m64 *mlo = (__m64 *) & dest[y][0]; __m64 *mhi = (__m64 *) & dest[y][w2]; __m64 *msour = (__m64 *) & sour[y][0]; for (unsigned int k = 0; k < w2 / 8; k++) { //k<w2/8 k=8*k __m64 even = _mm_packs_pu16(_mm_and_si64(*msour, m00FF), _mm_and_si64(*(msour + 1), m00FF)); //even coeffs __m64 odd = _mm_packs_pu16(_mm_srli_pi16(*msour, 8), _mm_srli_pi16(*(msour + 1), 8)); //odd coeffs addsub(even, odd, mlo++, mhi++); msour += 2; } if (w2 % 8) { for (unsigned int k = w2 - (w2 % 8); k < w2; k++) { dest[y][k] = char(((int)sour[y][2*k] + (int)sour[y][2*k+1]) / 2); dest[y][k+w2] = char(((int)sour[y][2*k] - (int)sour[y][2*k+1]) / 2); } } } _mm_empty(); } void Haar::transcols(char** dest, char** sour, unsigned int w, unsigned int h) const { unsigned int h2 = h / 2; for (unsigned int k = 0; k < h2; k++) { __m64 *mlo = (__m64 *) & dest[k][0]; __m64 *mhi = (__m64 *) & dest[k+h2][0]; __m64 *even = (__m64 *) & sour[2*k][0]; __m64 *odd = (__m64 *) & sour[2*k+1][0]; for (unsigned int x = 0; x < w / 8; x++) { addsub(*even, *odd, mlo, mhi); even++; odd++; mlo++; mhi++; } } _mm_empty(); //odd remainder for (unsigned int x = w - (w % 8); x < w; x++) { for (unsigned int k = 0; k < h2; k++) { dest[k][x] = char(((int)sour[2*k][x] + (int)sour[2*k+1][x]) / 2); dest[k+h2][x] = char(((int)sour[2*k][x] - (int)sour[2*k+1][x]) / 2); } } } //////// **************************************************************************************************************** 哈爾小波變換示例1.哈爾基函數(shù) 用q_00(x)表示這個(gè)常值函數(shù),，用V0表示由這個(gè)常值函數(shù)生成的矢量空間， 即V0：q_00(x)=1(0<=x<1)或0(其它),它是構(gòu)成矢量空間V0的基,。 如果一幅圖像由2^1=2個(gè)像素組成,，這幅圖像在[0,1) 區(qū)間中有兩個(gè)等間隔的子區(qū)間：[0,1/2) 和[1/2,1) ，每一個(gè)區(qū)間中各有1個(gè)常值函數(shù),，分別用q_10(x) q_11(x)表示,。用V1表示由2個(gè)子區(qū)間中的常值函數(shù)生成的矢量空間，即V1: q_10(x) =1(0<=x<1/2)或0(其它); q_11(x)=1(1/2<=x<1)或0(其它).這2個(gè)常值函數(shù)就是構(gòu)成矢量空間V1的基,。 如果一幅圖像由2^2= 4個(gè)像素組成,，這幅圖像在[0,1)區(qū)間中被分成4個(gè)等間隔的子區(qū)間：[0,1/4),[1/4,1/2),[1/2,3/4)和[3/4,1)，它們的 常值函數(shù)分別用q_20(x) q_21(x) q_22(x) q_23(x)表示,，用V2表示由4個(gè)子區(qū)間中的常值函數(shù)生成的矢量空間,，即V2:q_20(x)=1(0<=x<1/4)或0(其它); q_21(x)=1(1/4<=x<1/2)或0(其它); q_22(x) =1(1/2<=x<3/4)或0(其它); q_23(x)=1(3/4<=x<1)或0(其它).這4個(gè)常值函數(shù)就是構(gòu)成矢量空間V2的基,。 我們可以按照這種方法繼續(xù)定義基函數(shù)和由它生成的矢量空間?？偠灾?，為了表示矢量空間中的矢量，每一個(gè)矢量空間V_i都需要定義一個(gè)基(basis),。
為生成矢量空間V_i而定義的基函數(shù)也叫做尺度函數(shù),這種函數(shù)通常用符號(hào)q_ij(x)表示,。哈爾基函數(shù)定義為q(x)=1(0<=x<1)
或0(其它),哈爾尺度函數(shù)q_ij(x)定義為q_ij(x)=q(2^i·x-j),j=0,1,…,
2^i-1其中,i為尺度因子，改變i使函數(shù)圖形縮小或者放大,j為平移參數(shù),，改變j使函數(shù)沿x軸方向平移,。 小波函數(shù)通常用Ψ_ij(x)表示,。與哈爾基函數(shù)相對應(yīng)的小波稱為基本哈爾小波函數(shù),并由下式定義: Ψ(x)=1(0<=x<1/2)或
-1(1/2<=x<1)或0(其它); 哈爾小波尺度函數(shù)Ψ_ij(x)定義為 Ψ_ij(x)= Ψ(2^i·x-j),j=0,1,…,
2^i-1.用小波函數(shù)構(gòu)成的矢量空間Wi表示,，Wi=span{Ψ_ij(x)}, j=0,1,…, 2^i-1 根據(jù)哈爾小波函數(shù)的定義，可以寫出生成W0,，W1和W2等矢量空間的小波函數(shù),。 生成矢量空間W0的哈爾小波：Ψ_00(x)=1(0<=x<1/2)或-1(1/2<=x<1)或0(其它) 生成矢量空間W1的哈爾小波：Ψ_10(x)=1(0<=x<1/4)或-1(1/4<=x<1/2)或0(其它); Ψ_11(x)=1(1/2<=x<3/4)或-1(3/4<=x<1/2)或0(其它). 生成矢量空間W2的哈爾小波：Ψ_00(x)=1(0<=x<1/2)或-1(1/2<=x<1)或0(其它) Ψ_20(x)=1(0<=x<1/8)或-1(1/8<=x<2/8)或0(其它); Ψ_21(x)=1(2/8<=x<3/8)或-1(3/8<=x<4/8)或0(其它); Ψ_22(x)=1(4/8<=x<5/8)或-1(5/8<=x<6/8)或0(其它); Ψ_00(x)=1(6/8<=x<7/8)或-1(7/8<=x<1)或0(其它). 哈爾基和哈爾小波分別使用下面兩個(gè)式子進(jìn)行規(guī)范化 q_ij(x)=2^(i/2)·q(2^i·x-j), Ψ_ij(x)= 2^(i/2)·Ψ(2^i·x-j), 3.哈爾基的結(jié)構(gòu) 使用哈爾基函數(shù) q_ij(x)和哈爾小波函數(shù)Ψ_ij(x)生成的矢量空間Vi和Wi具有下面的性質(zhì)，V_i+1=V_i⊕W_i,。這就是說,，在矢量空間V_i+1中, 矢量空間W_i中的小波可用來表示一個(gè)函數(shù)在矢量空間V_i 中不能表示的部分。因此,，小波可定義為生成矢量空間W_i的一組線性無關(guān)的函數(shù)Ψ_ij(x)的集合,。 4.哈爾小波變換 小波變換的基本思想是用一組小波函數(shù)或者基函數(shù)表示一個(gè)函數(shù)或者信號(hào),，例如圖像信號(hào),。為了理解什么是小波變換,，下面用一個(gè)具體的例子來說明小波變換的過程,。假設(shè)有一幅分辨率只有4個(gè)像素的一維圖像,，對應(yīng)的像素值分別為[9 7 3 5] 計(jì)算它的哈爾小波變換系數(shù): (1).求均值(averaging)。計(jì)算相鄰像素對的平均值,，得到一幅分辨率比較低的新圖像,，它的像素?cái)?shù)目變成了2個(gè)，即新的圖像的分辨率是原來的1/2,，相應(yīng)的像素值為：[8 4] (2). 求差值(diqqerencing),。很明顯，用2個(gè)像素表示這幅圖像時(shí),，圖像的信息已經(jīng)部分丟失,。為了能夠從由2個(gè)像素組成的圖像重構(gòu)出由4個(gè)像素組成 的原始圖像,，就需要存儲(chǔ)一些圖像的細(xì)節(jié)系數(shù),以便在重構(gòu)時(shí)找回丟失的信息。方法是把像素對的第一個(gè)像素值減去這個(gè)像素對的平均值,，或者使用這個(gè)像素對的差 值除以2,。在這個(gè)例子中，第一個(gè)細(xì)節(jié)系數(shù)是(9-8)=1,，因?yàn)橛?jì)算得到的平均值是8,它比9小1而比7大1,，存儲(chǔ)這個(gè)細(xì)節(jié)系數(shù)就可以恢復(fù)原始圖像的前兩 個(gè)像素值。使用同樣的方法,，第二個(gè)細(xì)節(jié)系數(shù)是(3-4)=-1,，存儲(chǔ)這個(gè)細(xì)節(jié)系數(shù)就可以恢復(fù)后2個(gè)像素值。因此,，原始圖像就可以用下面的兩個(gè)平均值和兩個(gè) 細(xì)節(jié)系數(shù)表示[8 4 1 -1] (3).重復(fù)第1,，2步，把由第一步分解得到的圖像進(jìn)一步分解成分辨率更低的圖像[6 2 1 -1]. 由此可見,，通過上述分解就把由4像素組成的一幅圖像用一個(gè)平均像素值和三個(gè)細(xì)節(jié)系數(shù)表示,，這個(gè)過程就叫做哈爾小波變換,也稱哈爾小波分解(Haar wavelet decomposition)。這個(gè)概念可以推廣到使用其他小波基的變換,。在這個(gè)例子中我們可以看到,，①變換過程中沒有丟失信息，因?yàn)槟軌驈乃涗浀臄?shù)據(jù) 中重構(gòu)出原始圖像,。②對這個(gè)給定的變換,，我們可以從所記錄的數(shù)據(jù)中重構(gòu)出各種分辨率的圖像。例如,，在分辨率為1的圖像基礎(chǔ)上重構(gòu)出分辨率為2的圖像,，在分 辨率為2的圖像基礎(chǔ)上重構(gòu)出分辨率為4的圖像。③通過變換之后產(chǎn)生的細(xì)節(jié)系數(shù)的幅度值比較小,，這就為圖像壓縮提供了一種途徑,，例如去掉一些微不足道的細(xì)節(jié) 系數(shù)而不影響對重構(gòu)圖像的理解。 求均值和差值的過程實(shí)際上就是一維小波變換的過程,，現(xiàn)在用數(shù)學(xué)方法重新描述小波變換的過程,。 (1) 用V2中的哈爾基表示 圖像I(x)=[9 7 3 5]有2^j =2^2=4像素，因此可以用生成矢量空間V2中的基函數(shù)的線性組合表示,， I(x) =9
q_20(x)+
7q_21(x)+3q_22(x)+5q_23(x) (2)用V1和W1中的函數(shù)表示 根據(jù)V2=V1⊕W1,I(x)可表示成 I(x)=8q_10(x)+4q_11(x)+Ψ_10(x)-Ψ_11(x) (3)用V0,W0和W1中的函數(shù)表示 V2=V0⊕W0⊕W1,I(x)可表示成 I(x)
=6 q_00(x)+2Ψ_00(x)+Ψ_10(x)-Ψ_11(x) 5.規(guī)范化算法 規(guī)范化的小波變換與非規(guī)范化的小波變換相比,，唯一的差別是在規(guī)范化的小波變換中需要乘一個(gè)規(guī)范化的系數(shù)。下面用一個(gè)例子說明,。對函數(shù)f(x)= [2,5,8,9,7,4, -1, 1]作哈爾小波變換,。哈爾小波變換實(shí)際上是使用求均值和差值的方法進(jìn)行分解。我們把f(x)看成是矢量空間V3中的一個(gè)向量，尺度因子i= 3,，因此最多可分解為3層.分解過程如下,。 步驟1：c=(2 5,8 9,7,4,-1,1)/√2=(7,17,11,0, -3, -1,3,
-2)/√2 解釋： 平均值（7,17,11,0）/2 差值（-3，-1,3,，-2）/2（與平均值的差值相對于直接差值的1/2 即 A- (A+B)/2 = (A-B)/2 ,我們將除以2提出來） 所以步驟一的結(jié)果是（7,17,11,0）/2 并上（-3,，-1,3，-2）/2 = (7,17,11,0, -3, -1,3, -2)/2 步驟2：c=[(7+17)/√2,(11+0)/√2,(7-17)/√2,(11-0)/√2,-3,-1,3,-2]/√2 =(24/√2,11/√2,-10/√2,11/√2,-3,-1,3,-2)/√2 步驟3:c=[(24+11)/2,(24-11)/2 ,-10/√2,11/√2,-3,-1,3,-2]/√2 =(35/2,13/2,-10/√2,11/√2,-3,-1,3,-2)/√2 這個(gè)結(jié)果實(shí)際上是后面哈爾轉(zhuǎn)換N=8轉(zhuǎn)換矩陣相乘的結(jié)果,。 N= 8,，  。6.二維哈爾小波變換 前面已經(jīng)介紹了一維小波變換的基本原理和變換方法,。這節(jié)將結(jié)合具體的圖像數(shù)據(jù)系統(tǒng)地介紹如何使用小波對圖像進(jìn)行變換,。一幅圖像可看成是由許多像素組成的一 個(gè)矩陣，在進(jìn)行圖像壓縮時(shí),，為降低對存儲(chǔ)器的要求,人們通常把它分成許多小塊,，例如以8×8 個(gè)像素為一塊，并用矩陣表示,，然后分別對每一個(gè)圖像塊進(jìn)行處理,。在小波變換中，由于小波變換中使用的基函數(shù)的長度是可變的,，雖然無須像以離散余弦變換 (DCT)為基礎(chǔ)的JPEG標(biāo)準(zhǔn)算法那樣把輸入圖像進(jìn)行分塊以避免產(chǎn)生JPEG圖像那樣的“塊效應(yīng)”,，但為便于理解小波變換，還是從一個(gè)小的圖像塊入手,， 并且繼續(xù)使用哈爾小波對圖像進(jìn)行變換,。 假設(shè)有一幅灰度圖像，其中的一個(gè)圖像塊用矩陣 A表示為,， [64 2 3 61 60 6 7 57 9 55 54 12 13 51 50 16 17 47 46 20 21 43 42 24 40 26 27 37 36 30 31 33 32 34 35 29 28 38 39 25 41 23 22 44 45 19 18 48 49 15 14 52 53 11 10 56 8 58 59 5 4 62 63 1] 一個(gè)圖像塊是一個(gè)二維的數(shù)據(jù)陣列,，進(jìn)行小波變換時(shí)可以對陣列的每一行進(jìn)行變換，然后對行變換之后的陣列的每一列進(jìn)行變換,，最后對經(jīng)過變換之后的圖像數(shù)據(jù)陣列進(jìn)行編碼,。 從求均值(averaging)與求差值(differencing)開始。在圖像塊矩陣A中,第一行的像素值為[64 2 3 61 60 6 7 57] 步驟1：在第一行上取每一對像素的平均值,，并將結(jié)果放到第一行的前4個(gè)位置,，其余的4個(gè)數(shù)是R0行每一對像素的第一個(gè)數(shù)與其相應(yīng)的平均值之差，并將結(jié)果放到第一行的后4個(gè)位置 步驟2：對第一行前4個(gè)數(shù)使用與第一步相同的方法,，得到兩個(gè)平均值和兩個(gè)差（系數(shù)），并依次放在第一行的前4個(gè)位置,，其余的4個(gè)細(xì)節(jié)系數(shù)位置不動(dòng),。 步驟3：用與第1和2 步相同的方法，對剩余的一對平均值求平均值和差值。 使用求均值和求差值的方法,，對矩陣的每一行進(jìn)行計(jì)算,，得到A’ [32.5 0 0.5 0.5 31 -29 27 -25 32.5 0 -0.5 -0.5 -23 21 -19 17 32.5 0 -0.5 -0.5 -15 13 -11 9 32.5 0 0.5 0.5 7 -5 3 -1 32.5 0 0.5 0.5 -1 3 -5 7 32.5 0 -0.5 -0.5 9 -11 13 -15 32.5 0 -0.5 -0.5 17 -19 21 -23 32.5 0 0.5 0.5 -25 27 -29 31] 其中，每一行的第一個(gè)元素是該行像素值的平均值,，其余的是這行的細(xì)節(jié)系數(shù),。使用同樣的方法，對A’的每一列進(jìn)行計(jì)算,，得到A’’ [32.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 -4 4 -4 0 0 0 0 4 -4 4 -4 0 0 0.5 0.5 27 -25 23 -21 0 0 -0.5 -0.5 -11 9 -7 5 0 0 0.5 0.5 -5 7 -9 11 0 0 -0.5 -0.5 21 -23 25 -27] 其 中,，左上角的元素表示整個(gè)圖像塊的像素值的平均值，其余是該圖像塊的細(xì)節(jié)系數(shù),。根據(jù)這個(gè)事實(shí),，如果從矩陣中去掉表示圖像的某些細(xì)節(jié)系數(shù)，事實(shí)證明重構(gòu)的圖 像質(zhì)量仍然可以接受,。具體做法是設(shè)置一個(gè)閾值d,，例如d>=5的細(xì)節(jié)系數(shù)就把它當(dāng)作“0”看待，這樣經(jīng)過變換之后的矩陣就變成A’’’ [32.5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 27 -25 23 -21 0 0 0 0 -11 9 -7 0 0 0 0 0 0 7 -9 11 0 0 0 0 21 -23 25 -27] “0”的數(shù)目增加了18個(gè),，也就是去掉了18個(gè)細(xì)節(jié)系數(shù),。這樣做的好處是可提高編碼的效率。對A’’’矩陣進(jìn)行逆變換,，得到了重構(gòu)的近似矩陣AA [59.5 5.5 7.5 57.5 55.5 9.5 11.5 53.5 5.5 59.5 57.5 7.5 9.5 55.5 53.5 11.5 21.5 43.5 41.5 23.5 25.5 39.5 32.5 32.5 43.5 21.5 23.5 41.5 39.5 25.5 32.5 32.5 32.5 32.5 39.5 25.5 23.5 41.5 43.5 21.5 32.5 32.5 25.5 39.5 41.5 23.5 21.5 43.5 53.5 11.5 9.5 55.5 57.5 7.5 5.5 59.5 11.5 53.5 55.5 9.5 7.5 57.5 59.5 5.5] 如果矩陣A的數(shù)據(jù)與矩陣AA的數(shù)據(jù)用圖表示,，原圖和重構(gòu)圖差別不是很大 轉(zhuǎn)： http://blog.csdn.net/mmmmmmmmnnnnxx/article/details/5169111 Another example: 小波變換的基本思想是用一組小波函數(shù)或者基函數(shù)表示一個(gè)函數(shù)或者信號(hào)，例如圖像信號(hào),。為了理解什么是小波變換,，下面用一個(gè)具體的例子來說明小波變換的過程。 1. 求有限信號(hào)的均值和差值 [例8. 1]
假設(shè)有一幅分辨率只有4個(gè)像素 的一維圖像,，對應(yīng)的像素值或者叫做圖像位置的系數(shù)分別為： 計(jì)算步驟如下： [8 4] 步 驟2：求差值(differencing),。很明顯，用2個(gè)像素表示這幅圖像時(shí),，圖像的信息已經(jīng)部分丟失,。為了能夠從由2個(gè)像素組成的圖像重構(gòu)出由4個(gè)像 素組成的原始圖像，就需要存儲(chǔ)一些圖像的細(xì)節(jié)系數(shù)(detail coefficient),，以便在重構(gòu)時(shí)找回丟失的信息,。方法是把像素對的第一個(gè)像素值減去這個(gè)像素對的平均值,，或者使用這個(gè)像素對的差值除以2。在這個(gè) 例子中,，第一個(gè)細(xì)節(jié)系數(shù)是(9-8)=1,，因?yàn)橛?jì)算得到的平均值是8,它比9小1而比7大1，存儲(chǔ)這個(gè)細(xì)節(jié)系數(shù)就可以恢復(fù)原始圖像的前兩個(gè)像素值,。使用同 樣的方法,，第二個(gè)細(xì)節(jié)系數(shù)是(3-4)=-1，存儲(chǔ)這個(gè)細(xì)節(jié)系數(shù)就可以恢復(fù)后2個(gè)像素值,。因此,，原始圖像就可以用下面的兩個(gè)平均值和兩個(gè)細(xì)節(jié)系數(shù)表示， [8 4 1 -1] 步驟3：重復(fù)第1,，2步,，把由第一步分解得到的圖像進(jìn)一步分解成分辨率更低的圖像和細(xì)節(jié)系數(shù)。在這個(gè)例子中,，分解到最后,，就用一個(gè)像素的平均值6和三個(gè)細(xì)節(jié)系數(shù)2,1和-1表示整幅圖像。 [6 2 1 -1] 這個(gè)分解過程如表8-1所示,。 表8-1 哈爾變換過程
由此可見,，通過上述分解就把由4像素組成的一幅圖像用一個(gè)平均像素值和三個(gè)細(xì)節(jié)系數(shù)表示，這個(gè)過程就叫做哈爾小波變換(Haar
wavelet transform),，也稱哈爾小波分解(Haar wavelet
decomposition),。這個(gè)概念可以推廣到使用其他小波基的變換。 ======================================================================== =原理 ============================== 哈爾小波變換是于1909年由Alfréd Haar所提出,，是小波變換(Wavelet transform)中最簡單的一種變換，也是最早提出的小波變換,。他是多貝西小波的于N=2的特例,，可稱之為D2 哈爾小波的母小波(mother wavelet)可表示為： 且對應(yīng)的縮放方程式(scaling function)可表示為：
特性哈爾小波具有如下的特性： (1)任一函數(shù)都可以由 (2)任一函數(shù)都可以由常函數(shù)， (3) 正交性(Orthogonal) (4)不同寬度的(不同m)的wavelet/scaling functions之間會(huì)有一個(gè)關(guān)系 φ(t) = φ(2t) + φ(2t? 1) ψ(t) = φ(2t) ? φ(2t? 1) (5)可以用 m+1的 系數(shù)來計(jì)算 m 的系數(shù) 若
哈爾變換Haar Transform最早是由A. Haar在1910年“Zur theorie der orthogonalen funktionensysteme”中所提出,，是一種最簡單又可以反應(yīng)出時(shí)變頻譜（time-variant spectrum）的表示方法,。其觀念與Fourier Transform相近,，F(xiàn)ourier Transform的原理是利用弦波sine與cosine來對信號(hào)進(jìn)行調(diào)制；而Haar Transform則是利用Haar function來對信號(hào)進(jìn)行調(diào)制,。Haar function也含有sine、cosine所擁有的正交性,，也就是說不同的Haar function是互相orthogonal,，其內(nèi)積為零。 以下面N=8的哈爾變換矩陣為例,，我們?nèi)〉谝涣泻偷诙衼碜鰞?nèi)積,，得到的結(jié)果為零；取第二列和第三列來做內(nèi)積,，得到的結(jié)果也是零,。依序下去，我們可以發(fā)現(xiàn)在哈爾變換矩陣任取兩列來進(jìn)行內(nèi)積的運(yùn)算,，所得到的內(nèi)積皆為零,。
Haar Transform有以下幾點(diǎn)特性： 1.不需要乘法（只有相加或加減） 2.輸入與輸出個(gè)數(shù)相同 3.頻率只分為低頻（直流值）與高頻（1和-1）部分 4.可以分析一個(gè)信號(hào)的Localized feature 5.運(yùn)算速度極快,，但不適合用于信號(hào)分析 6.大部分運(yùn)算為0，不用計(jì)算 7.維度小,，使用的memory少 8.因?yàn)榇蟛糠譃楦哳l,，變換較籠統(tǒng) 對一矩陣A做哈爾小波變換的公式為B= HAHT，其中A為一
哈爾小波變換應(yīng)用于圖像壓縮
范例
哈爾小波變換運(yùn)算量比沃爾什變換更少
|
|
|
以及它們的位移函數(shù)所組成
以及它們的位移函數(shù)所組成
,。
的區(qū)塊且H為N點(diǎn)的哈爾小波變換,。而反哈爾小波變換為A= HBHT,。以下為H在2、4及8點(diǎn)時(shí)的值：
