opencv+小波变换

什么是小波变换

  小波变换，作为图像处理中的一个重要分支，在图像压缩，医疗图像处理中有着比较好的应用，在一定程度上甚至可以取代傅里叶变换。与傅里叶变换不同，由于小波变换是建立在多尺度上的，因而可以有着更高的灵活度。这里为了让各位更好地理解小波，先打个比方。存在一家做产品的公司，该公司等级制度明确，其中有一个负责产品方向的CEO，以及一个负责技术的CTO。此时CEO决定做产品A，跟同为大佬的CTO商量了下能不能做，CTO说能做，于是CEO就把产品的定位，人群跟他的小弟产品经理交代，产品经理就去找CTO的小弟技术大拿商量技术路线，技术大拿确定好了技术路线，于是就给他的技术小弟们开会布置技术任务。这个期间产品经理的小弟们运营就天天去烦这些技术小弟，商讨如何实现一个个功能。最后，经过若干场撕逼，产品A就出炉了。这里面，把一张图比作产品A，那么小波分解就是把产品A分解成整个公司的实现过程。每一层平级的同事就是一个尺度，站在CEO这边的运营人员就是基于尺度的时域部分，他们明确这个产品A是怎么样的，基于客户确定产品A的功能，是能够看得见的部分，而站在CTO这边的技术人员就是基于细节的频域部分，他们实现客户所看不见的部分。

  小波分解

这里先上一张小波分解之后的图

该小波图由两部分组成，一部分为左上角Opencv女神组成的时域部分,其他的就是频域部分（当然这里为了能够让各位看到女神我把浮点型转成整型了）。

进行小波分解的公式有很多种，基于图像的离散小波一般有草帽小波、对称小波、高斯小波、拉普拉斯小波，这里为了后面能够比较顺利地进行还原，采用拉普拉斯小波分解。拉普拉斯小波很容易理解，比方我们有8个数分别为8、6、8、10、12、10、12、8，该数组记为A，那么我们两个1组求平均数求得新的一组为7、9、11、10，记为B1，这B1就是A1低一尺度下的时域部分，有时域自然就得有频域，用B1减去原来数组A的一半，当然得间断地相减，得到-1，-1，-1，-2，记为B2，很明显存在另一组频域部分B3为1,1,1,2，与B2互为相反数，所以可以不保存下来。以此类推，可以再对B1进行分解，得C1为8、10.5，C2为1、-0.5，如下图

8	6	8	10	12	10	12	8
7	9	11	10	-1	-1	-1	-2
8	10.5	1	0.5	-1	-1	-1	-2

那么或许有人会疑问，原来8个数，最后还是8个数，有什么用呢？小波变换当然不可能是脱啥放啥，他能解决的问题可多了。①掌控全局，如果要增大整个数组，在数组中你得分别对8个数进行操作，但在低尺度下，你只需要对两个数进行操作，然后小波合成回去②掌控细节，如果要使整个方差变大，在原数组中是很难操作的，但在低尺度下，你只需要对频域进行操作，比如最下面一行你只需要使大的更大，小的更小，这在图像中就是增强的操作。又能掌控全局，又能掌控细节，你还有什么理由不为小波鼓掌，下面上干货

 

#include <opencv2/opencv.hpp>
#include <math.h>
#include <iostream>
#include <imgproc/imgproc.hpp>  
using namespace std;
using namespace cv;

//小波分解
void laplace_decompose(Mat src,int s,Mat &wave)
{
    Mat full_src(src.rows, src.cols, CV_32FC1);
    Mat dst = src.clone();
    dst.convertTo(dst, CV_32FC1);
    
    for (int m = 0; m < s; m++)
    {    
        dst.convertTo(dst, CV_32FC1);
        Mat wave_src(dst.rows, dst.cols, CV_32FC1);
        //列变换，detail=mean-original
        for (int i = 0; i < wave_src.rows; i++)
        {
            for (int j = 0; j < wave_src.cols/2; j++)
            {
                wave_src.at<float>(i, j) = (dst.at<float>(i, 2 * j) + dst.at<float>(i, 2 * j + 1)) / 2;
                wave_src.at<float>(i, j + wave_src.cols/2) = wave_src.at<float>(i, j) - dst.at<float>(i, 2 * j);
            }
        }
        Mat temp = wave_src.clone();
        for (int i = 0; i < wave_src.rows/2; i++)
        {
            for (int j = 0; j < wave_src.cols / 2; j++)
            {
                wave_src.at<float>(i, j) = (temp.at<float>(2 * i, j) + temp.at<float>(2 * i + 1, j)) / 2;
                wave_src.at<float>(i + wave_src.rows / 2, j) = wave_src.at<float>(i, j) - temp.at<float>(2 * i, j);
            }
        }    
        dst.release();
        dst = wave_src(Rect(0, 0, wave_src.cols / 2, wave_src.rows /2));
        wave_src.copyTo(full_src(Rect(0, 0, wave_src.cols, wave_src.rows)));
    }
    wave = full_src.clone();
}
//小波复原
void wave_recover(Mat full_scale, Mat &original,int level)
{
    //每一个循环中把一个级数的小波还原
    for (int m = 0; m < level; m++)
    {    
        Mat temp = full_scale(Rect(0, 0, full_scale.cols / pow(2, level - m-1), full_scale.rows / pow(2, level - m-1)));
        
        //先恢复左边
        Mat recover_src(temp.rows, temp.cols, CV_32FC1);
        for (int i = 0; i < recover_src.rows; i++)
        {
            for (int j = 0; j < recover_src.cols/2; j++)
            {
                if (i % 2 == 0)
                    recover_src.at<float>(i, j) = temp.at <float>(i / 2, j) - temp.at<float>(i / 2 + recover_src.rows / 2, j);
                else
                    recover_src.at<float>(i, j) = temp.at <float>(i / 2, j)+ temp.at<float>(i / 2 + recover_src.rows / 2, j);
            }
        }
        Mat temp2 = recover_src.clone();
        //再恢复整个
        for (int i = 0; i < recover_src.rows; i++)
        {
            for (int j = 0; j < recover_src.cols; j++)
            {
                if (j % 2 == 0)
                    recover_src.at<float>(i, j) = temp2.at<float>(i, j / 2) - temp.at<float>(i, j / 2 + temp.cols / 2);
                else
                    recover_src.at<float>(i, j) = temp2.at<float>(i, j / 2) + temp.at<float>(i, j / 2 + temp.cols / 2);
            }
        }
        recover_src.copyTo(temp);
    }
    original = full_scale.clone();
    original.convertTo(original, CV_8UC1);

    
}

//小波操作
void ware_operate(Mat &full_scale, int level)
{
    //取出最低尺度的那一层，对其进行操作，仅最低尺度那层可以对时域进行操作，其他层只能对频域进行操作
    Mat temp = full_scale(Rect(0, 0, full_scale.cols / 4, full_scale.rows /4));
    temp = temp(Rect(0, 0, temp.cols / 2, temp.rows / 2));
    Mat temp2 = temp.clone();
    //这里对时域操作，降低灰度
    for (int i = 0; i < temp2.rows;i++)
    for (int j = 0; j < temp2.cols; j++)
        temp2.at<float>(i, j) -= 20;
    temp2.copyTo(temp);
    
    //这里对频域操作，拉伸细节
    //先处理左下角
    for (int i = temp.rows / 2; i < temp.rows; i++)
    {
        for (int j = 0; j < temp.cols / 2; j++)
        {
            if (temp.at<float>(i, j)>0)
                temp.at<float>(i, j) += 5;
            if (temp.at<float>(i, j) < 0)
                temp.at<float>(i, j) -= 5;
        }
    }
    //再处理右半边
    for (int i = 0; i < temp.rows; i++)
    {
        for (int j = 0; j < temp.cols; j++)
        {
            if (temp.at<float>(i, j)>0)
                temp.at<float>(i, j) += 5;
            if (temp.at<float>(i, j)<0)
                temp.at<float>(i, j) -= 5;
        }
    }
}

//小波分解
Mat 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);

    //下采样
    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;
}
void main()
{
    
    Mat src = imread("timg.jpg",0);
    imshow("original", src);
    Mat full_src;
    laplace_decompose(src, 3, full_src);
    ware_operate(full_src, 3);
    Mat src_recover;
    wave_recover(full_src, src_recover, 3);
    imshow("recover", src_recover);
    
    waitKey();
}

这是原图

下面上一张效果图