另外参见俄罗斯写的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.
转:
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的图像。
③
通过变换之后产生的细节系数的幅度值比较小,这就为图像压缩提供了一种途径,例如去掉一些微不足道的细节系数并不影响对重构图像的理解。
========================================================================
=原理
==============================
哈尔小波的母小波(mother wavelet)可表示为:
且对应的缩放方程式(scaling function)可表示为:
目录
|
哈尔小波具有如下的特性: (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做哈尔小波变换的公式为B= HAHT,其中A为一
的区块且H为N点的哈尔小波变换。而反哈尔小波变换为A= HBHT。以下为H在2、4及8点时的值:N= 2,
。N= 4,
。N= 8,
。此外,当N= 2k时,
。其中H除了第0个row为φ(φ=[1 1 1 ... 1]/,共N个1),第2p+ q个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 Time | terms required for NRMSE < 10 ? 5 | |
---|---|---|
离散傅里叶变换 | 9.5秒 | 43 |
沃尔什变换 | 2.2秒 | 65 |
哈尔小波变换 | 0.3秒 | 128 |
联系客服