opencv 小波变换与逆变换C++

4 篇文章 0 订阅
3 篇文章 0 订阅

在此记录一下

兴幸在网上看到关于小波变换的代码,但是在逆变换的时候结果跟matlab的有很大差别,因此对照一下matlab的具体代码,对已有的代码进行一点改动。

说明:

1.配置好opencv,就可以运行(附Demo)

2.小波包的选择,代码中只包括了haar,db1,sym2,如需要更多种类的小波包,可以在matlab里查看wfilters函数,对应写进代码中。

 

代码如下:

1.  WaveTransform.h文件

#include<opencv.hpp>
#include <imgproc/imgproc.hpp>
using namespace cv;
class WaveTransform
{
public:
	WaveTransform(void);
	~WaveTransform(void);
	Mat WDT(const Mat &_src,const string _wname,const int _level);//小波分解
	Mat IWDT(const Mat &_src,const string _wname,const int _level);//小波重构
	void wavelet_D( const string _wname, Mat &_lowFilter, Mat &_highFilter );//分解包
	void wavelet_R( const string _wname, Mat &_lowFilter, Mat &_highFilter );//重构包
	Mat waveletDecompose( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter );
	Mat waveletReconstruct(const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter);
};

2. .WaveTransform.cpp文件

#include "WaveTransform.h"


WaveTransform::WaveTransform(void)
{
}


WaveTransform::~WaveTransform(void)
{
}

Mat WaveTransform::WDT(const Mat &_src,const string _wname,const int _level)
{
	//int reValue=THID_ERR_NONE;
    Mat_<float> src=Mat_<float>(_src);
    Mat dst=Mat::zeros(src.rows,src.cols,src.type());
    int row=src.rows;
    int col=src.cols;
    //高通低通滤波器
    Mat lowFilter;
    Mat highFilter;
    wavelet_D(_wname,lowFilter,highFilter);
    //小波变换
    int t=1;
	

    while (t<=_level)
    {
        //先进行 行小波变换
//#pragma omp parallel for
        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);
            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("dst1.jpg",&dstImg1);
#endif
		
        //小波列变换
//#pragma omp parallel for
        for (int j=0;j<col;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);//dst,not src
            }
            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("dst2.jpg",&dstImg2);
#endif
		/*
        char s[10];
        itoa(t,s,10);
        imshow(s,Mat_<uchar>(dst));
        waitKey(1);
		*/
		//归化各子图范围0~255
		/*
		int r_len=row/2,c_len=col/2;
		for(int i=0;i<2;i++)
		{
			for(int j=0;j<2;j++)
			{
				Point p1=Point(i*c_len,j*r_len);
				Point p2=Point((i+1)*c_len,(j+1)*r_len);
				Mat ROI=dst(Rect(p1,p2));
				//ROI=cv::abs(ROI);
				normalize(ROI,ROI,0,255,CV_MINMAX);
				//imshow("a",Mat_<uchar>(ROI));
				//waitKey(0);
			}
		}
		*/

		//
        //更新 
        row/=2;
        col/=2;
        t++;
        src=dst;      
    }
    return dst;
}

Mat WaveTransform::IWDT(const Mat &_src,const string _wname,const int _level)
{
	//int reValue=THID_ERR_NONE;
	Mat src=Mat_<float>(_src);
	Mat dst;//=Mat::zeros(src.rows,src.cols,src.type());
	src.copyTo(dst);
	int N=src.rows;
	int D= src.cols;

	//高低通滤波器
	Mat lowFilter;
	Mat highFilter;
	wavelet_R(_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)
	//while(t<=_level)
	{
		//列逆变换
		for(int j=0;j<col;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
		IplImage dstImg2=IplImage(dst);
		cvSaveImage("dst.jpg",&dstImg2);
#endif

		//行逆变换
		for(int i=0;i<row;i++)
		{
			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);
			for(int j=0;j<col;j++)
			{
				dst.at<float>(i,j)=oneRow.at<float>(0,j);
			}
		}

#if 0
		IplImage dstImg1=IplImage(dst);
		cvSaveImage("dst.jpg",&dstImg1);
#endif
		
        char s[10];
        itoa(t,s,10);
		
		//Rect rrr=Rect(Point(col-1,row-1),Point(src.cols-1,src.rows-1));
		//Rect rrr=Rect(Point(0,0),Point(col-1,row-1));
		/*
		Mat showImg;//=//dst;//(rrr);
		dst.copyTo(showImg);
		Mat showImg1;
		showImg.copyTo(showImg1);
		normalize(showImg1,showImg1,0,255,CV_MINMAX);
        imshow(s,Mat_<uchar>(showImg1));
        waitKey(1);
		*/
		row*=2;
		col*=2;
		t++;
		src=dst;
	}

	return dst;
}


void WaveTransform::wavelet_D( const string _wname, Mat &_lowFilter, Mat &_highFilter )
{
	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);
    }
    else if (_wname=="sym2")
    {
        int N=4;
        float h[]={-0.4830, 0.8365, -0.2241, -0.1294};
        float l[]={-0.1294, 0.2241,  0.8365, 0.4830};
        
        _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];
        }
    }
}
void WaveTransform::wavelet_R( const string _wname, Mat &_lowFilter, Mat &_highFilter )
{
	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);
    }
    else if (_wname=="sym2")
    {
        int N=4;
        float h[]={-0.1294,-0.2241,0.8365,-0.4830};
        float l[]={0.4830, 0.8365, 0.2241, -0.1294};
        
        _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 WaveTransform::waveletDecompose( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter )
{
    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);

    //下采样
    //数据拼接
    for (int i=0,j=1;i<D/2;i++,j+=2)
    {
        src.at<float>(0,i)=dst1.at<float>(0,j);//lowFilter
        src.at<float>(0,i+D/2)=dst2.at<float>(0,j);//highFilter
    }
    return src;
}


Mat WaveTransform::waveletReconstruct(const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter)
{
	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());

	//Mat roi1(src, Rect(0, 0, D / 2, 1));
	//Mat roi2(src, Rect(D / 2, 0, D / 2, 1));
	/// 插值为0
	
	for ( int i=0, cnt=0; 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 ); ///< 后一半
	}
	//std::cout<<Up1<<std::endl;

	/// 线性插值
	
	//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;

}

3. Demo

#include "WaveTransform.h"
#include<iostream>
#include<opencv.hpp>
using namespace cv;
using namespace std;
int main()
{
    char*filename="1.jpg";
	Mat src=imread(filename);
	
	int level=3;//分解阶次
	double dishu=2;
	int result=(int)pow(dishu,level);

	WaveTransform m_waveTransform;
	
	//double a=clock();
	
	//resize(src,src,Size((512/result)*result,(512/result)*result));
	
	Mat img;
	cvtColor(src,img,CV_RGB2GRAY);	
	normalize(img,img,0,255,CV_MINMAX);
	imshow("img",img);
	Mat float_src;
	img.convertTo(float_src,CV_32F);
	
	Mat imgWave=m_waveTransform.WDT(float_src,"sym2",level);	//haar,sym2
	imgWave.convertTo(float_src,CV_32F);
	Mat imgIWave=m_waveTransform.IWDT(float_src,"sym2",level);
	imshow("imgWave",Mat_<uchar>(imgWave));
	normalize(imgIWave,imgIWave,0,255,CV_MINMAX);
	imshow("IWDT",Mat_<uchar>(imgIWave));
	
	waitKey(0)
    return 0;
}

4. 效果图

原图:

小波分解图:(为了显示效果,对每个子图灰度值进行了normalization,WDT程序中注释掉部分)

小波重构图:

 

完!

  • 20
    点赞
  • 84
    收藏
    觉得还不错? 一键收藏
  • 28
    评论
小波变换是一种时频分析方法,可以将信号分解成不同频率的子带,从而更好地理解信号的特征。下面是使用C++OpenCV进行一维连续小波变换的示例代码: ```cpp #include <opencv2/opencv.hpp> #include <opencv2/core.hpp> #include <opencv2/highgui.hpp> #include <opencv2/imgproc.hpp> #include <opencv2/imgcodecs.hpp> #include <opencv2/core/utility.hpp> #include <iostream> #include <vector> #include <cmath> using namespace cv; using namespace std; // 小波变换函数 void wavelet(Mat& src, Mat& dst, int level) { Mat temp = src.clone(); for (int i = 0; i < level; i++) { int rows = temp.rows; int cols = temp.cols; Mat LL, LH, HL, HH; pyrDown(temp(Rect(0, 0, cols, rows)), LL); pyrDown(temp(Rect(0, rows / 2, cols, rows / 2)), LH); pyrDown(temp(Rect(cols / 2, 0, cols / 2, rows)), HL); pyrDown(temp(Rect(cols / 2, rows / 2, cols / 2, rows / 2)), HH); Mat LHHL; pyrUp(LH, LHHL, LL.size()); Mat HHLH; pyrUp(HL, HHLH, LL.size()); Mat HHLL; pyrUp(HH, HHLL, LL.size()); Mat LL_new; LL.copyTo(LL_new); for (int i = 0; i < LL.rows; i++) { for (int j = 0; j < LL.cols; j++) { LL_new.at<float>(i, j) = (LL.at<float>(i, j) + LHHL.at<float>(i, j) + HHLH.at<float>(i, j) + HHLL.at<float>(i, j)) / 4; } } temp = LL_new.clone(); } dst = temp.clone(); } int main() { // 生成测试信号 int N = 512; Mat signal(N, 1, CV_32FC1); for (int i = 0; i < N; i++) { signal.at<float>(i, 0) = sin(2 * CV_PI * i / 32) + sin(2 * CV_PI * i / 16); } // 进行小波变换 Mat dst; wavelet(signal, dst, 3); // 显示结果 Mat dst_show; normalize(dst, dst_show, 0, 255, NORM_MINMAX); dst_show.convertTo(dst_show, CV_8UC1); imshow("Wavelet Transform", dst_show); waitKey(0); return 0; } ``` 上述代码中,我们首先生成了一个长度为512的测试信号,然后调用wavelet函数进行小波变换,最后将结果进行归一化并显示出来。在wavelet函数中,我们使用了OpenCV自带的pyrDown和pyrUp函数进行了一维小波变换

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 28
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值