打开APP
userphoto
未登录

开通VIP,畅享免费电子书等14项超值服

开通VIP
哈尔小波变换的原理及其实现(Haar)

另外参见俄罗斯写的http://www.codeproject.com/Articles/22243/Real-Time-Object-Tracker-in-C

Haar小波在图像处理和数字水印等方面应用较多,这里简单的介绍一下哈尔小波的基本原理以及其实现情况。

一、Haar小波的基本原理

数学理论方面的东西我也不是很熟悉,这边主要用简单的例子来介绍下Haar小波的使用情况。

例如:有a=[8,7,6,9]四个数,并使用b[4]数组来保存结果.

则一级Haar小波变换的结果为:

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

即依次从数组中取两个数字,计算它们的和以及差,并将和一半和差的一半依次保存在数组的前半部分和后半部分。

例如:有a[8],要进行一维Haar小波变换,结果保存在b[8]中

则一级Haar小波变换的结果为:

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

如果需要进行二级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 );
#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;

}

===========================================================================

*****************************************************************************************************************

其他代码实现

*****************************************************************************************************************

// 哈尔小波.cpp : 定义控制台应用程序的入口点。

//

#include "stdafx.h"

#include <iostream>

#include "cv.h"

#include "highgui.h"

#include <ctype.h>

// 二维离散小波变换(单通道浮点图像)

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 _tmain(int argc, _TCHAR* argv[])

{

// 小波变换层数

int nLayer = 1;

// 输入彩色图像

IplImage *pSrc = cvLoadImage("lena.bmp", 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);

cvNamedWindow("pWavelet",CV_WINDOW_AUTOSIZE);

cvShowImage("pWavelet",pWavelet);

cvNamedWindow("pSrc",CV_WINDOW_AUTOSIZE);

cvShowImage("pSrc",pSrc);

cvWaitKey(0);

cvResetImageROI(pWavelet); // 本行代码有点多余,但有利用养成良好的编程习惯

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.哈尔基函数
最简单的基函数是哈尔基函数(Haar basis function)。哈尔基函数在1909年提出,它是由一组分段常值函数组成的函数集。这个函数集定义在半开区间[0,1)上,每一个分段常值函数的数 值在一个小范围里是1,其他地方为0,现以图像为例并使用线性代数中的矢量空间来说明哈尔基函数。
如果一幅图像仅由2^0=1个像素组成,这幅图像在整个[0,1) 区间中就是一个常值函数。

用q_00(x)表示这个常值函数,用V0表示由这个常值函数生成的矢量空间,

即V0:q_00(x)=1(0<=x<1)或0(其它),它是构成矢量空间V0的基。

如果一幅图像由2^1=2个像素组成,这幅图像在[0,1) 区间中有两个等间隔的子区间:[0,1/2) 和[1/2,1) ,每一个区间中各有1个常值函数,分别用q_10(x) q_11(x)表示。用V1表示由2个子区间中的常值函数生成的矢量空间,即V1: q_10(x)  =1(0<=x<1/2)或0(其它); q_11(x)=1(1/2<=x<1)或0(其它).这2个常值函数就是构成矢量空间V1的基。

如果一幅图像由2^2= 4个像素组成,这幅图像在[0,1)区间中被分成4个等间隔的子区间:[0,1/4),[1/4,1/2),[1/2,3/4)和[3/4,1),它们的 常值函数分别用q_20(x) q_21(x) q_22(x) q_23(x)表示,用V2表示由4个子区间中的常值函数生成的矢量空间,即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个常值函数就是构成矢量空间V2的基。

我们可以按照这种方法继续定义基函数和由它生成的矢量空间。总而言之,为了表示矢量空间中的矢量,每一个矢量空间V_i都需要定义一个基(basis)。 为生成矢量空间V_i而定义的基函数也叫做尺度函数,这种函数通常用符号q_ij(x)表示。哈尔基函数定义为q(x)=1(0<=x<1) 或0(其它),哈尔尺度函数q_ij(x)定义为q_ij(x)=q(2^i·x-j),j=0,1,…, 2^i-1其中,i为尺度因子,改变i使函数图形缩小或者放大,j为平移参数,改变j使函数沿x轴方向平移。
空间矢量V_i定义为V_i=span{ q_ij(x)},j=0,1,…, 2^i-1.由于定义了基和矢量空间,我们就可以把由2^i个像素组成的一维图像看成为矢量空间V_i中的矢量。由于这些矢量都是在单位区间[0,1)上 定义的函数,所以在V_i矢量空间中的每一个矢量也被包含在V_i+1矢量空间中,即V0∈V1∈...∈V_i∈V_i+1.

2.哈尔小波函数
小波函数通常用Ψ_ij(x)表示。与哈尔基函数相对应的小波称为基本哈尔小波函数,并由下式定义:
Ψ(x)=1(0<=x<1/2)或 -1(1/2<=x<1)或0(其它); 
哈尔小波尺度函数Ψ_ij(x)定义为
Ψ_ij(x)= Ψ(2^i·x-j),j=0,1,…, 2^i-1.用小波函数构成的矢量空间Wi表示,Wi=span{Ψ_ij(x)}, j=0,1,…, 2^i-1
根据哈尔小波函数的定义,可以写出生成W0,W1和W2等矢量空间的小波函数。
生成矢量空间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(其它).
哈尔基和哈尔小波分别使用下面两个式子进行规范化
q_ij(x)=2^(i/2)·q(2^i·x-j),
Ψ_ij(x)= 2^(i/2)·Ψ(2^i·x-j),

3.哈尔基的结构
使用哈尔基函数 q_ij(x)和哈尔小波函数Ψ_ij(x)生成的矢量空间Vi和Wi具有下面的性质,V_i+1=V_i⊕W_i。这就是说,在矢量空间V_i+1中, 矢量空间W_i中的小波可用来表示一个函数在矢量空间V_i 中不能表示的部分。因此,小波可定义为生成矢量空间W_i的一组线性无关的函数Ψ_ij(x)的集合。
4.哈尔小波变换
小波变换的基本思想是用一组小波函数或者基函数表示一个函数或者信号,例如图像信号。为了理解什么是小波变换,下面用一个具体的例子来说明小波变换的过程。假设有一幅分辨率只有4个像素的一维图像,对应的像素值分别为[9 7 3 5]
计算它的哈尔小波变换系数:
(1).求均值(averaging)。计算相邻像素对的平均值,得到一幅分辨率比较低的新图像,它的像素数目变成了2个,即新的图像的分辨率是原来的1/2,相应的像素值为:[8 4]
(2). 求差值(diqqerencing)。很明显,用2个像素表示这幅图像时,图像的信息已经部分丢失。为了能够从由2个像素组成的图像重构出由4个像素组成 的原始图像,就需要存储一些图像的细节系数,以便在重构时找回丢失的信息。方法是把像素对的第一个像素值减去这个像素对的平均值,或者使用这个像素对的差 值除以2。在这个例子中,第一个细节系数是(9-8)=1,因为计算得到的平均值是8,它比9小1而比7大1,存储这个细节系数就可以恢复原始图像的前两 个像素值。使用同样的方法,第二个细节系数是(3-4)=-1,存储这个细节系数就可以恢复后2个像素值。因此,原始图像就可以用下面的两个平均值和两个 细节系数表示[8 4 1 -1]
(3).重复第1,2步,把由第一步分解得到的图像进一步分解成分辨率更低的图像[6 2 1 -1].
由此可见,通过上述分解就把由4像素组成的一幅图像用一个平均像素值和三个细节系数表示,这个过程就叫做哈尔小波变换,也称哈尔小波分解(Haar wavelet decomposition)。这个概念可以推广到使用其他小波基的变换。在这个例子中我们可以看到,①变换过程中没有丢失信息,因为能够从所记录的数据 中重构出原始图像。②对这个给定的变换,我们可以从所记录的数据中重构出各种分辨率的图像。例如,在分辨率为1的图像基础上重构出分辨率为2的图像,在分 辨率为2的图像基础上重构出分辨率为4的图像。③通过变换之后产生的细节系数的幅度值比较小,这就为图像压缩提供了一种途径,例如去掉一些微不足道的细节 系数而不影响对重构图像的理解。
求均值和差值的过程实际上就是一维小波变换的过程,现在用数学方法重新描述小波变换的过程。
(1) 用V2中的哈尔基表示
图像I(x)=[9 7 3 5]有2^j =2^2=4像素,因此可以用生成矢量空间V2中的基函数的线性组合表示,
I(x) =9 q_20(x)+ 7q_21(x)+3q_22(x)+5q_23(x)

(2)用V1和W1中的函数表示
根据V2=V1⊕W1,I(x)可表示成
I(x)=8q_10(x)+4q_11(x)+Ψ_10(x)-Ψ_11(x)

(3)用V0,W0和W1中的函数表示
V2=V0⊕W0⊕W1,I(x)可表示成

I(x) =6 q_00(x)+2Ψ_00(x)+Ψ_10(x)-Ψ_11(x)

5.规范化算法
规范化的小波变换与非规范化的小波变换相比,唯一的差别是在规范化的小波变换中需要乘一个规范化的系数。下面用一个例子说明。对函数f(x)= [2,5,8,9,7,4, -1, 1]作哈尔小波变换。哈尔小波变换实际上是使用求均值和差值的方法进行分解。我们把f(x)看成是矢量空间V3中的一个向量,尺度因子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提出来)

所以步骤一的结果是(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

这个结果实际上是后面哈尔转换N=8转换矩阵相乘的结果。
N= 8,

6.二维哈尔小波变换
前面已经介绍了一维小波变换的基本原理和变换方法。这节将结合具体的图像数据系统地介绍如何使用小波对图像进行变换。一幅图像可看成是由许多像素组成的一 个矩阵,在进行图像压缩时,为降低对存储器的要求,人们通常把它分成许多小块,例如以8×8 个像素为一块,并用矩阵表示,然后分别对每一个图像块进行处理。在小波变换中,由于小波变换中使用的基函数的长度是可变的,虽然无须像以离散余弦变换 (DCT)为基础的JPEG标准算法那样把输入图像进行分块以避免产生JPEG图像那样的“块效应”,但为便于理解小波变换,还是从一个小的图像块入手, 并且继续使用哈尔小波对图像进行变换。
假设有一幅灰度图像,其中的一个图像块用矩阵 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]
一个图像块是一个二维的数据阵列,进行小波变换时可以对阵列的每一行进行变换,然后对行变换之后的阵列的每一列进行变换,最后对经过变换之后的图像数据阵列进行编码。
从求均值(averaging)与求差值(differencing)开始。在图像块矩阵A中,第一行的像素值为[64 2  3  61  60  6  7  57]
步骤1:在第一行上取每一对像素的平均值,并将结果放到第一行的前4个位置,其余的4个数是R0行每一对像素的第一个数与其相应的平均值之差,并将结果放到第一行的后4个位置
步骤2:对第一行前4个数使用与第一步相同的方法,得到两个平均值和两个差(系数),并依次放在第一行的前4个位置,其余的4个细节系数位置不动。
步骤3:用与第1和2 步相同的方法,对剩余的一对平均值求平均值和差值。
使用求均值和求差值的方法,对矩阵的每一行进行计算,得到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]
其中,每一行的第一个元素是该行像素值的平均值,其余的是这行的细节系数。使用同样的方法,对A’的每一列进行计算,得到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]
其 中,左上角的元素表示整个图像块的像素值的平均值,其余是该图像块的细节系数。根据这个事实,如果从矩阵中去掉表示图像的某些细节系数,事实证明重构的图 像质量仍然可以接受。具体做法是设置一个阈值d,例如d>=5的细节系数就把它当作“0”看待,这样经过变换之后的矩阵就变成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”的数目增加了18个,也就是去掉了18个细节系数。这样做的好处是可提高编码的效率。对A’’’矩阵进行逆变换,得到了重构的近似矩阵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的数据与矩阵AA的数据用图表示,原图和重构图差别不是很大

转:

http://blog.csdn.net/mmmmmmmmnnnnxx/article/details/5169111

Another example:

小波变换的基本思想是用一组小波函数或者基函数表示一个函数或者信号,例如图像信号。为了理解什么是小波变换,下面用一个具体的例子来说明小波变换的过程。

1. 求有限信号的均值和差值

[例8. 1] 假设有一幅分辨率只有4个像素 的一维图像,对应的像素值或者叫做图像位置的系数分别为:
[9 7 3 5]
计算它的哈尔小波变换系数。

计算步骤如下:
步骤1:求均值(averaging)。计算相邻像素对的平均值,得到一幅分辨率比较低的新图像,它的像素数目变成了2个,即新的图像的分辨率是原来的1/2,相应的像素值为:

[8 4]

步 骤2:求差值(differencing)。很明显,用2个像素表示这幅图像时,图像的信息已经部分丢失。为了能够从由2个像素组成的图像重构出由4个像 素组成的原始图像,就需要存储一些图像的细节系数(detail coefficient),以便在重构时找回丢失的信息。方法是把像素对的第一个像素值减去这个像素对的平均值,或者使用这个像素对的差值除以2。在这个 例子中,第一个细节系数是(9-8)=1,因为计算得到的平均值是8,它比9小1而比7大1,存储这个细节系数就可以恢复原始图像的前两个像素值。使用同 样的方法,第二个细节系数是(3-4)=-1,存储这个细节系数就可以恢复后2个像素值。因此,原始图像就可以用下面的两个平均值和两个细节系数表示,

[8 4 1 -1]

步骤3:重复第1,2步,把由第一步分解得到的图像进一步分解成分辨率更低的图像和细节系数。在这个例子中,分解到最后,就用一个像素的平均值6和三个细节系数2,1和-1表示整幅图像。

[6 2 1 -1]

这个分解过程如表8-1所示。

表8-1 哈尔变换过程

分辨率

平均值

细节系数

4

[9 7 3 5]


2

[8 4]

[1 -1]

1

[6]

[2]

由此可见,通过上述分解就把由4像素组成的一幅图像用一个平均像素值和三个细节系数表示,这个过程就叫做哈尔小波变换(Haar wavelet transform),也称哈尔小波分解(Haar wavelet decomposition)。这个概念可以推广到使用其他小波基的变换。
从这个例子中我们可以看到:
① 变换过程中没有丢失信息,因为能够从所记录的数据中重构出原始图像。
② 对这个给定的变换,我们可以从所记录的数据中重构出各种分辨率的图像。例如,在分辨率为1的图像基础上重构出分辨率为2的图像,在分辨率为2的图像基础上重构出分辨率为4的图像。
③ 通过变换之后产生的细节系数的幅度值比较小,这就为图像压缩提供了一种途径,例如去掉一些微不足道的细节系数并不影响对重构图像的理解。

========================================================================

=原理

==============================

哈尔小波变换是于1909年由Alfréd Haar所提出,是小波变换(Wavelet transform)中最简单的一种变换,也是最早提出的小波变换。他是多贝西小波的于N=2的特例,可称之为D2

哈尔小波的母小波(mother wavelet)可表示为:

且对应的缩放方程式(scaling function)可表示为:

目录

  • 1 特性

  • 2 哈尔变换

  • 3 哈尔小波变换应用于图像压缩

    • 3.1 说明

    • 3.2 范例

  • 4 哈尔小波变换运算量比沃尔什变换更少

  • 5 参考

特性

哈尔小波具有如下的特性: (1)任一函数都可以由

以及它们的位移函数所组成

(2)任一函数都可以由常函数,

以及它们的位移函数所组成

(3) 正交性(Orthogonal)


(4)不同宽度的(不同m)的wavelet/scaling functions之间会有一个关系

φ(t) = φ(2t) + φ(2t? 1)

ψ(t) = φ(2t) ? φ(2t? 1)

(5)可以用 m+1的 系数来计算 m 的系数 若

哈尔变换

Haar Transform最早是由A. Haar在1910年“Zur theorie der orthogonalen funktionensysteme”中所提出,是一种最简单又可以反应出时变频谱(time-variant spectrum)的表示方法。其观念与Fourier Transform相近,Fourier Transform的原理是利用弦波sine与cosine来对信号进行调制;而Haar Transform则是利用Haar function来对信号进行调制。Haar function也含有sine、cosine所拥有的正交性,也就是说不同的Haar function是互相orthogonal,其内积为零。

以下面N=8的哈尔变换矩阵为例,我们取第一列和第二列来做内积,得到的结果为零;取第二列和第三列来做内积,得到的结果也是零。依序下去,我们可以发现在哈尔变换矩阵任取两列来进行内积的运算,所得到的内积皆为零。

  • N= 8,


在此前提下,利用Fourier Transform的观念,假设所要分析的信号可以使用多个频率与位移不同的Haar function来组合而成,进行Haar Transform时,因为Haar function的正交性,便可求出信号在不同Haar function(不同频率)的情况下所占有的比例。

Haar Transform有以下几点特性:

1.不需要乘法(只有相加或加减)

2.输入与输出个数相同

3.频率只分为低频(直流值)与高频(1和-1)部分

4.可以分析一个信号的Localized feature

5.运算速度极快,但不适合用于信号分析

6.大部分运算为0,不用计算

7.维度小,使用的memory少

8.因为大部分为高频,变换较笼统

对一矩阵A做哈尔小波变换的公式为BHAHT,其中A为一

的区块且HN点的哈尔小波变换。而反哈尔小波变换为AHBHT。以下为H在2、4及8点时的值:

  • N= 2,

  • N= 4,

  • N= 8,

  • 此外,当N= 2k时,

    。其中H除了第0个row为φ(φ=[1 1 1 ... 1]/
    ,共N个1),第2pq个row为hp,q

哈尔小波变换应用于图像压缩


  • 由于数字图片档案过大,因此我们往往会对图片做图像压缩,压缩过后的档案大小不仅存放于电脑中不会占到过大容量,也方便我们于网络上传送。哈尔小波变换其中一种应用便是用来压缩图像。压缩图像的基本概念为将图像存成到一矩阵,矩阵中的每一元素则代表是每一图像的某画素值,介于0到255间。例如256x256大小的图片会存成256x256大小的矩阵。JPEG影像压缩的概念为先将图像切成8x8大小的区块,每一区块为一8x8的矩阵。示意图可见右图。

  • 在处理8x8二维矩阵前,先试着对一维矩阵

    作哈尔小波变换,

  • 公式为

范例

  • 对8x8的二维矩阵A作哈尔小波变换,由于AH是对A的每一行作哈尔小波变换,作完后还要对A的每一列作哈尔小波变换,因此公式为HTAH。以下为一简单的例子:

  • 列哈尔小波变换(row Haar wavelet transform)

  • 行哈尔小波变换(column Haar wavelet transform)

  • 由以上例子可以看出哈尔小波变换的效果,原本矩阵中变化量不大的元素经过变换后会趋近零,再配合适当量化便可以达到压缩的效果了。此外若一矩阵作完哈尔小波变换后所含的零元素非常多的话,称此矩阵叫稀疏,若一矩阵越稀疏压缩效果越好。因此可对定一临界值ε若矩阵中元素的绝对值小于此临界值ε,可将该元素令成零,可得到更大的压缩率。然而ε取过大的话会造成图像严重失真,因此如何取适当的ε也是值得讨论的议题。

哈尔小波变换运算量比沃尔什变换更少

  • 若应用于区域的频谱分析及侦测边缘的话,离散傅立叶变换、Walsh-Hadamard变换及哈尔小波变换的计算量见下表


Running Timeterms required for NRMSE < 10 ? 5
离散傅里叶变换9.5秒43
沃尔什变换2.2秒65
哈尔小波变换0.3秒128
本站仅提供存储服务,所有内容均由用户发布,如发现有害或侵权内容,请点击举报
打开APP,阅读全文并永久保存 查看更多类似文章
猜你喜欢
类似文章
【热】打开小程序,算一算2024你的财运
小波变换-我的理解(图像Haar小波变换)
OpenCV-高斯低通&高通滤波器(C++)
图像处理基础:图像的灰度变换
图像的灰度直方图、直方图均衡化、直方图规定化(匹配)
OpenCV学习笔记(五十六)
使用openCV与C++求图片特征值及特征向量并进行图片处理
更多类似文章 >>
生活服务
热点新闻
分享 收藏 导长图 关注 下载文章
绑定账号成功
后续可登录账号畅享VIP特权!
如果VIP功能使用有故障,
可点击这里联系客服!

联系客服