小波变换笔记

实现haar小波

参考文章中讲解了haar小波的原理，比较易懂，实现了多级haar小波分解，本文参照它进行了练习，同时添加了小波重建。

# include<iostream>
# include<opencv2/opencv.hpp>

using namespace std;
using namespace cv;

int main()
{
    Mat src = imread("F:\\testdata\\pic\\16.jpg");
    Mat HSV;
    cvtColor(src, HSV, CV_BGR2HSV);
    vector<Mat> channels;
    split(HSV, channels);

    Mat img = channels.at(1);
        //上述是获取HSV中的饱和度图像，按需修改
    namedWindow("img", 0);
    imshow("img", img);
    int height = img.rows;
    int width = img.cols;
    int depth = 2;
    int depthcount = 1;

    Mat tmp = Mat::ones(Size(width, height), CV_32FC1);
    Mat wavelet = Mat::ones(Size(width, height), CV_32FC1);
    Mat imgtmp = img.clone();

    imgtmp.convertTo(imgtmp, CV_32FC1);

    //---------------------------------------小波分解-----------------------------//
    while (depthcount <= depth)
    {
        height = img.rows / pow(2,depthcount-1);
        width = img.cols / pow(2,depthcount-1);

        //水平方向变换
        for (int i = 0; i < height; ++i)
        {
            for (int j = 0; j < width / 2; ++j)
            {
                tmp.at<float>(i, j) = (imgtmp.at<float>(i, 2 * j) + imgtmp.at<float>(i, 2 * j + 1)) / 2;
                tmp.at<float>(i, j+width/2) = (imgtmp.at<float>(i, 2 * j) - imgtmp.at<float>(i, 2 * j + 1)) / 2;
            }
        }
        
        //垂直方向变换
        for (int i = 0; i < height / 2; ++i)
        {
            for (int j = 0; j < width; ++j)
            {
                wavelet.at<float>(i, j) = (tmp.at<float>(2 * i, j) + tmp.at<float>(2 * i + 1, j)) / 2;
                wavelet.at<float>(i+height/2, j) = (tmp.at<float>(2 * i, j) - tmp.at<float>(2 * i + 1, j)) / 2;
            }
        }
        /*
        //低通滤波，选用高斯低通滤波器
        Mat ROI = wavelet(Rect(0, 0, width, height));
        GaussianBlur(ROI, ROI, Size(7, 7), 0.5);
        Mat dst(wavelet, Rect(0, 0, width, height));
        ROI.copyTo(dst);
        */
        imgtmp = wavelet;
        depthcount++;


    }

    
    //------------------------------------------小波重建--------------------------------------//
    while (depth > 0)
    {
        height = img.rows / pow(2, depth - 1);
        width = img.cols / pow(2, depth - 1);
        //列逆变换
        for (int i = 0; i < height/2; ++i)
        {
            for (int j = 0; j < width; ++j)
            {
                float value1 = imgtmp.at<float>(i, j);
                float value2 = imgtmp.at<float>(i+height/2, j);
                tmp.at<float>(2 * i, j) = value1 + value2;
                tmp.at<float>(2 * i+1, j) = value1 - value2;
            }
        }
        
        //行逆变换
        for (int i = 0; i < height; ++i)
        {
            for (int j = 0; j < width / 2; ++j)
            {
                float value1 = tmp.at<float>(i, j);
                float value2 = tmp.at<float>(i , j+width/2);
                wavelet.at<float>(i, 2*j) = value1 + value2;
                wavelet.at<float>(i , 2*j+1) = value1 - value2;
            }
        }

        Mat ROI=wavelet(Rect(0, 0, width, height));
        Mat dst(imgtmp, Rect(0, 0, width, height));
        ROI.copyTo(dst);
        depth--;
    }



    namedWindow("res", 0);
    imgtmp.convertTo(imgtmp, CV_8UC1);
    imshow("res", imgtmp);
    waitKey(0);
    return 0;
}

---

开源代码实现多种小波函数的小波变换

• C++ wavelet library
• 可以选择按照其文档中的配置方法，调用动态链接库，但是我的尝试不成功，出现了内存异常，所以本文直接将其源码放入自己的项目中，和自己的项目一起编译，可以实现同样的功能。需要注意的是，需要配置好FFTW这个依赖，方法见VS2015配置FFTW
• 函数说明
• 利用这个库中的函数可以得到小波分解的近似系数和细节系数

使用示例

重要的只有如下两行函数调用

dwt_2d(data, depth, name, output, flag, length);

idwt_2d(output, flag, name, res, length);

#include <iostream>
#include <fstream>
#include "wavelet2d.h"
#include <vector>
#include <string>
#include <cmath>
#include<opencv2\opencv.hpp>
using namespace std;
using namespace cv;

int main() {
  //读入图片，获取饱和度
  Mat src = imread("F:\\testdata\\poreImgV3\\10.jpg");
  resize(src, src, Size(450, 450));
  Mat HSV;
  cvtColor(src, HSV, CV_BGR2HSV);
  vector<Mat> channels;
  channels.resize(3);
  split(HSV, channels);
  Mat saturation = channels.at(1);
  namedWindow("s", 0);
  imshow("s", saturation);
  waitKey(0);
  //将饱和度图像数据存入vector<vector<double>>中,一行一行放置
  int row = saturation.rows;
  int col = saturation.cols;
  //分配大小
  vector<vector<double>> data;
  data.resize(row);
  for (int i = 0; i < data.size(); ++i)
  {
      data[i].resize(col);
  }
  //放入数据
  for (int i = 0; i < row; ++i)
  {
      for (int j = 0; j < col; ++j)
      {
          data[i][j] = saturation.at<uchar>(i, j);
      }
  }

  int depth = 2;
  string name = "haar";
  vector<double> output;
  //output.resize(row*col);
  vector<double> flag;
  vector<int> length; 
  //小波分解
  dwt_2d(data, depth, name, output, flag, length);

  vector<vector<double>> res;
  res.resize(row);
  for (int i = 0; i < res.size(); ++i)
  {
      res[i].resize(col);
  }
  //----------------------------------------------------------------------------------

  //小波重建
  idwt_2d(output, flag, name, res, length);
  //写回原图像中
  Mat recovery(src.size(), CV_8UC1);
  for (int i = 0; i < row; ++i)
  {
      for (int j = 0; j < col; ++j)
      {
          recovery.at<uchar>(i, j) = (uchar)res[i][j];
      }
  }
  imshow("res", recovery);
  waitKey(0);
  system("pause");
  return 0;
}