openCV中的小波变换

Question

是否有人尝试在 opencv 或 C++ 中实现 DWT？ 我看到了关于这个主题的旧帖子，但我发现它们对我没有用，因为我需要一个近似系数和作为小波变换结果的细节。

我尝试将this 添加到我的项目中，但效果不如预期。

这很简单，因为结果参数我需要近似系数和细节：

void haar1(float *vec, int n, int w)
{
int i=0;
float *vecp = new float[n];
for(i=0;i<n;i++)
    vecp[i] = 0;

    w/=2;
    for(i=0;i<w;i++)
    {
        vecp[i] = (vec[2*i] + vec[2*i+1])/sqrt(2.0);
        vecp[i+w] = (vec[2*i] - vec[2*i+1])/sqrt(2.0);
    }
    
    for(i=0;i<(w*2);i++)
            vec[i] = vecp[i];

        delete [] vecp;
}
void haar2(float **matrix, int rows, int cols)
{
    float *temp_row = new float[cols];
    float *temp_col = new float[rows];

    int i=0,j=0;
    int w = cols, h=rows;
while(w>1 || h>1)
{
    if(w>1)
    {
        for(i=0;i<h;i++)
        {
            for(j=0;j<cols;j++)
                temp_row[j] = matrix[i][j];

            haar1(temp_row,cols,w);
            
            for(j=0;j<cols;j++)
                matrix[i][j] = temp_row[j];
        }
    }

    if(h>1)
    {
        for(i=0;i<w;i++)
        {
            for(j=0;j<rows;j++)
                temp_col[j] = matrix[j][i];
            haar1(temp_col, rows, h);
            for(j=0;j<rows;j++)
                matrix[j][i] = temp_col[j];
        }
    }

    if(w>1)
        w/=2;
    if(h>1)
        h/=2;
}

    delete [] temp_row;
    delete [] temp_col;
}

所以有人可以帮我找到用 C++ 实现的 dwt，或者告诉我如何从上面的代码系数中提取。谢谢

Answer 1

这里是直接和逆 Haar 小波变换（用于滤波）：

#include "opencv2/opencv.hpp"
#include <iostream>
#include <vector>
#include <stdio.h>

using namespace cv;
using namespace std;

// Filter type
#define NONE 0  // no filter
#define HARD 1  // hard shrinkage
#define SOFT 2  // soft shrinkage
#define GARROT 3  // garrot filter
//--------------------------------
// signum
//--------------------------------
float sgn(float x)
{
    float res=0;
    if(x==0)
    {
        res=0;
    }
    if(x>0)
    {
        res=1;
    }
    if(x<0)
    {
        res=-1;
    }
    return res;
}
//--------------------------------
// Soft shrinkage
//--------------------------------
float soft_shrink(float d,float T)
{
    float res;
    if(fabs(d)>T)
    {
        res=sgn(d)*(fabs(d)-T);
    }
    else
    {
        res=0;
    }

    return res;
}
//--------------------------------
// Hard shrinkage
//--------------------------------
float hard_shrink(float d,float T)
{
    float res;
    if(fabs(d)>T)
    {
        res=d;
    }
    else
    {
        res=0;
    }

    return res;
}
//--------------------------------
// Garrot shrinkage
//--------------------------------
float Garrot_shrink(float d,float T)
{
    float res;
    if(fabs(d)>T)
    {
        res=d-((T*T)/d);
    }
    else
    {
        res=0;
    }

    return res;
}
//--------------------------------
// Wavelet transform
//--------------------------------
static void cvHaarWavelet(Mat &src,Mat &dst,int NIter)
{
    float c,dh,dv,dd;
    assert( src.type() == CV_32FC1 );
    assert( dst.type() == CV_32FC1 );
    int width = src.cols;
    int height = src.rows;
    for (int k=0;k<NIter;k++) 
    {
        for (int y=0;y<(height>>(k+1));y++)
        {
            for (int x=0; x<(width>>(k+1));x++)
            {
                c=(src.at<float>(2*y,2*x)+src.at<float>(2*y,2*x+1)+src.at<float>(2*y+1,2*x)+src.at<float>(2*y+1,2*x+1))*0.5;
                dst.at<float>(y,x)=c;

                dh=(src.at<float>(2*y,2*x)+src.at<float>(2*y+1,2*x)-src.at<float>(2*y,2*x+1)-src.at<float>(2*y+1,2*x+1))*0.5;
                dst.at<float>(y,x+(width>>(k+1)))=dh;

                dv=(src.at<float>(2*y,2*x)+src.at<float>(2*y,2*x+1)-src.at<float>(2*y+1,2*x)-src.at<float>(2*y+1,2*x+1))*0.5;
                dst.at<float>(y+(height>>(k+1)),x)=dv;

                dd=(src.at<float>(2*y,2*x)-src.at<float>(2*y,2*x+1)-src.at<float>(2*y+1,2*x)+src.at<float>(2*y+1,2*x+1))*0.5;
                dst.at<float>(y+(height>>(k+1)),x+(width>>(k+1)))=dd;
            }
        }
        dst.copyTo(src);
    }   
}
//--------------------------------
//Inverse wavelet transform
//--------------------------------
static void cvInvHaarWavelet(Mat &src,Mat &dst,int NIter, int SHRINKAGE_TYPE=0, float SHRINKAGE_T=50)
{
    float c,dh,dv,dd;
    assert( src.type() == CV_32FC1 );
    assert( dst.type() == CV_32FC1 );
    int width = src.cols;
    int height = src.rows;
    //--------------------------------
    // NIter - number of iterations 
    //--------------------------------
    for (int k=NIter;k>0;k--) 
    {
        for (int y=0;y<(height>>k);y++)
        {
            for (int x=0; x<(width>>k);x++)
            {
                c=src.at<float>(y,x);
                dh=src.at<float>(y,x+(width>>k));
                dv=src.at<float>(y+(height>>k),x);
                dd=src.at<float>(y+(height>>k),x+(width>>k));

               // (shrinkage)
                switch(SHRINKAGE_TYPE)
                {
                case HARD:
                    dh=hard_shrink(dh,SHRINKAGE_T);
                    dv=hard_shrink(dv,SHRINKAGE_T);
                    dd=hard_shrink(dd,SHRINKAGE_T);
                    break;
                case SOFT:
                    dh=soft_shrink(dh,SHRINKAGE_T);
                    dv=soft_shrink(dv,SHRINKAGE_T);
                    dd=soft_shrink(dd,SHRINKAGE_T);
                    break;
                case GARROT:
                    dh=Garrot_shrink(dh,SHRINKAGE_T);
                    dv=Garrot_shrink(dv,SHRINKAGE_T);
                    dd=Garrot_shrink(dd,SHRINKAGE_T);
                    break;
                }

                //-------------------
                dst.at<float>(y*2,x*2)=0.5*(c+dh+dv+dd);
                dst.at<float>(y*2,x*2+1)=0.5*(c-dh+dv-dd);
                dst.at<float>(y*2+1,x*2)=0.5*(c+dh-dv-dd);
                dst.at<float>(y*2+1,x*2+1)=0.5*(c-dh-dv+dd);            
            }
        }
        Mat C=src(Rect(0,0,width>>(k-1),height>>(k-1)));
        Mat D=dst(Rect(0,0,width>>(k-1),height>>(k-1)));
        D.copyTo(C);
    }   
}
//--------------------------------
//
//--------------------------------
int process(VideoCapture& capture)
{
    int n = 0;
    const int NIter=4;
    char filename[200];
    string window_name = "video | q or esc to quit";
    cout << "press space to save a picture. q or esc to quit" << endl;
    namedWindow(window_name, CV_WINDOW_KEEPRATIO); //resizable window;
    Mat frame;
    capture >> frame;

    Mat GrayFrame=Mat(frame.rows, frame.cols, CV_8UC1);
    Mat Src=Mat(frame.rows, frame.cols, CV_32FC1);
    Mat Dst=Mat(frame.rows, frame.cols, CV_32FC1);
    Mat Temp=Mat(frame.rows, frame.cols, CV_32FC1);
    Mat Filtered=Mat(frame.rows, frame.cols, CV_32FC1);
    for (;;) 
    {
        Dst=0;
        capture >> frame;
        if (frame.empty()) continue;
        cvtColor(frame, GrayFrame, CV_BGR2GRAY);
        GrayFrame.convertTo(Src,CV_32FC1);
        cvHaarWavelet(Src,Dst,NIter);

        Dst.copyTo(Temp);

        cvInvHaarWavelet(Temp,Filtered,NIter,GARROT,30);

        imshow(window_name, frame);

        double M=0,m=0;
        //----------------------------------------------------
        // Normalization to 0-1 range (for visualization)
        //----------------------------------------------------
        minMaxLoc(Dst,&m,&M);
        if((M-m)>0) {Dst=Dst*(1.0/(M-m))-m/(M-m);}
        imshow("Coeff", Dst);

        minMaxLoc(Filtered,&m,&M);
        if((M-m)>0) {Filtered=Filtered*(1.0/(M-m))-m/(M-m);}        
        imshow("Filtered", Filtered);

        char key = (char)waitKey(5);
        switch (key) 
        {
        case 'q':
        case 'Q':
        case 27: //escape key
            return 0;
        case ' ': //Save an image
            sprintf(filename,"filename%.3d.jpg",n++);
            imwrite(filename,frame);
            cout << "Saved " << filename << endl;
            break;
        default:
            break;
        }
    }
    return 0;
}

int main(int ac, char** av) 
{
    VideoCapture capture(0);
    if (!capture.isOpened()) 
    {
        return 1;
    }
    return process(capture);
}

Answer 2

Here 是来自Mahavir 的 OpenCV 中小波变换的另一种实现：

            #include <opencv2\highgui\highgui.hpp>
            #include <opencv2\core\core.hpp>
            #include <opencv2\core\mat.hpp>
            #include <opencv2\imgproc\imgproc.hpp>
            #include<iostream>
            #include<math.h>
            #include<conio.h>
            using namespace std;
            using namespace cv;

            class image
            {
            public:
                Mat im,im1,im2,im3,im4,im5,im6,temp,im11,im12,im13,im14,imi,imd,imr;
                float a,b,c,d;
                int getim();
            };

            int image::getim()
            {
                im=imread("lena.jpg",0); //Load image in Gray Scale
                imi=Mat::zeros(im.rows,im.cols,CV_8U);
                im.copyTo(imi);

                im.convertTo(im,CV_32F,1.0,0.0);
                im1=Mat::zeros(im.rows/2,im.cols,CV_32F);
                im2=Mat::zeros(im.rows/2,im.cols,CV_32F);
                im3=Mat::zeros(im.rows/2,im.cols/2,CV_32F);
                im4=Mat::zeros(im.rows/2,im.cols/2,CV_32F);
                im5=Mat::zeros(im.rows/2,im.cols/2,CV_32F);
                im6=Mat::zeros(im.rows/2,im.cols/2,CV_32F);

                //--------------Decomposition-------------------

                for(int rcnt=0;rcnt<im.rows;rcnt+=2)
                {
                    for(int ccnt=0;ccnt<im.cols;ccnt++)
                    {

                        a=im.at<float>(rcnt,ccnt);
                        b=im.at<float>(rcnt+1,ccnt);
                        c=(a+b)*0.707;
                        d=(a-b)*0.707;
                        int _rcnt=rcnt/2;
                        im1.at<float>(_rcnt,ccnt)=c;
                        im2.at<float>(_rcnt,ccnt)=d;
                    }
                }

                for(int rcnt=0;rcnt<im.rows/2;rcnt++)
                {
                    for(int ccnt=0;ccnt<im.cols;ccnt+=2)
                    {

                        a=im1.at<float>(rcnt,ccnt);
                        b=im1.at<float>(rcnt,ccnt+1);
                        c=(a+b)*0.707;
                        d=(a-b)*0.707;
                        int _ccnt=ccnt/2;
                        im3.at<float>(rcnt,_ccnt)=c;
                        im4.at<float>(rcnt,_ccnt)=d;
                    }
                }

                for(int rcnt=0;rcnt<im.rows/2;rcnt++)
                {
                    for(int ccnt=0;ccnt<im.cols;ccnt+=2)
                    {

                        a=im2.at<float>(rcnt,ccnt);
                        b=im2.at<float>(rcnt,ccnt+1);
                        c=(a+b)*0.707;
                        d=(a-b)*0.707;
                        int _ccnt=ccnt/2;
                        im5.at<float>(rcnt,_ccnt)=c;
                        im6.at<float>(rcnt,_ccnt)=d;
                    }
                }

                imr=Mat::zeros(256,256,CV_32F);
                imd=Mat::zeros(256,256,CV_32F);
                im3.copyTo(imd(Rect(0,0,128,128)));
                im4.copyTo(imd(Rect(0,127,128,128)));
                im5.copyTo(imd(Rect(127,0,128,128)));
                im6.copyTo(imd(Rect(127,127,128,128)));


                //---------------------------------Reconstruction-------------------------------------

                im11=Mat::zeros(im.rows/2,im.cols,CV_32F);
                im12=Mat::zeros(im.rows/2,im.cols,CV_32F);
                im13=Mat::zeros(im.rows/2,im.cols,CV_32F);
                im14=Mat::zeros(im.rows/2,im.cols,CV_32F);

                for(int rcnt=0;rcnt<im.rows/2;rcnt++)
                {
                    for(int ccnt=0;ccnt<im.cols/2;ccnt++)
                    {
                        int _ccnt=ccnt*2;
                        im11.at<float>(rcnt,_ccnt)=im3.at<float>(rcnt,ccnt);     //Upsampling of stage I
                        im12.at<float>(rcnt,_ccnt)=im4.at<float>(rcnt,ccnt);
                        im13.at<float>(rcnt,_ccnt)=im5.at<float>(rcnt,ccnt);
                        im14.at<float>(rcnt,_ccnt)=im6.at<float>(rcnt,ccnt);
                    }
                }


                for(int rcnt=0;rcnt<im.rows/2;rcnt++)
                {
                    for(int ccnt=0;ccnt<im.cols;ccnt+=2)
                    {

                        a=im11.at<float>(rcnt,ccnt);
                        b=im12.at<float>(rcnt,ccnt);
                        c=(a+b)*0.707;
                        im11.at<float>(rcnt,ccnt)=c;
                        d=(a-b)*0.707;                           //Filtering at Stage I
                        im11.at<float>(rcnt,ccnt+1)=d;
                        a=im13.at<float>(rcnt,ccnt);
                        b=im14.at<float>(rcnt,ccnt);
                        c=(a+b)*0.707;
                        im13.at<float>(rcnt,ccnt)=c;
                        d=(a-b)*0.707;
                        im13.at<float>(rcnt,ccnt+1)=d;
                    }
                }

                temp=Mat::zeros(im.rows,im.cols,CV_32F);

                for(int rcnt=0;rcnt<im.rows/2;rcnt++)
                {
                    for(int ccnt=0;ccnt<im.cols;ccnt++)
                    {

                        int _rcnt=rcnt*2;
                        imr.at<float>(_rcnt,ccnt)=im11.at<float>(rcnt,ccnt);     //Upsampling at stage II
                        temp.at<float>(_rcnt,ccnt)=im13.at<float>(rcnt,ccnt); 
                    }
                }

                for(int rcnt=0;rcnt<im.rows;rcnt+=2)
                {
                    for(int ccnt=0;ccnt<im.cols;ccnt++)
                    {

                        a=imr.at<float>(rcnt,ccnt);
                        b=temp.at<float>(rcnt,ccnt);
                        c=(a+b)*0.707;
                        imr.at<float>(rcnt,ccnt)=c;                                      //Filtering at Stage II
                        d=(a-b)*0.707;
                        imr.at<float>(rcnt+1,ccnt)=d;
                    }
                }

                imd.convertTo(imd,CV_8U);
                namedWindow("Input Image",1);
                imshow("Input Image",imi);
                namedWindow("Wavelet Decomposition",1);
                imshow("Wavelet Decomposition",imd);
                imr.convertTo(imr,CV_8U);
                namedWindow("Wavelet Reconstruction",1);
                imshow("Wavelet Reconstruction",imr);
                waitKey(0);
                return 0;
            }

            int main()
            {
                image my;
                my.getim();
                return 0;
            }

希望有人觉得它有用！

Answer 3

我看到java中小波的代码示例很少，特别是如果您使用的是openCV。我不得不在 java 中使用带有 openCV 的小波，我使用来自 @la luvia 的 C 代码并转换为 java。

翻译代码时遇到了很多麻烦，因为它在openCV的方法和使用方式上有很多不同。这本书在这个过程中也帮助了我很多。

我希望这段代码和这本书能提供一些关于如何使用 lib 的观点以及 C 和 Java 之间的一些差异。

代码如下：

import org.opencv.core.Core;
import org.opencv.core.CvType;
import org.opencv.core.Mat;
import org.opencv.imgcodecs.Imgcodecs;
import org.opencv.imgproc.Imgproc;

public class Wavelet {

    //Imperative in java
    static{
        System.loadLibrary(Core.NATIVE_LIBRARY_NAME);
    }

    String pathname = "C:/Users/user/img096";


    public static void main(String[] args) {
        Wavelet wavelet = new Wavelet();
        wavelet.applyHaarFoward();
        wavelet.applyHaarReverse();

        Imgcodecs.imwrite(wavelet.pathname+"imi.jpg", wavelet.imi);
        Imgcodecs.imwrite(wavelet.pathname+"imd.jpg", wavelet.imd);
        Imgcodecs.imwrite(wavelet.pathname+"imr.jpg", wavelet.imr);
    }

    Mat im,im1,im2,im3,im4,im5,im6,temp,im11,im12,im13,im14,imi,imd,imr;
    float a,b,c,d;

    private void applyHaarFoward(){
        try{
            im = Imgcodecs.imread(pathname+".jpg", 0);
            imi = new Mat(im.rows(), im.cols(), CvType.CV_8U);
            im.copyTo(imi);

            //in CvType. If the number of channels is omitted, it evaluates to 1. 
            im.convertTo(im, CvType.CV_32F, 1.0, 0.0);
            im1 = new Mat(im.rows()/2, im.cols(), CvType.CV_32F);
            im2 = new Mat(im.rows()/2, im.cols(), CvType.CV_32F);

            im3 = new Mat(im.rows()/2, im.cols()/2, CvType.CV_32F);
            im4 = new Mat(im.rows()/2, im.cols()/2, CvType.CV_32F);

            im5 = new Mat(im.rows()/2, im.cols()/2, CvType.CV_32F);
            im6 = new Mat(im.rows()/2, im.cols()/2, CvType.CV_32F);


            // ------------------- Decomposition ------------------- 

            for (int rcnt = 0; rcnt < im.rows(); rcnt+=2) {
                for (int ccnt = 0; ccnt < im.cols(); ccnt++) {
                    //even though the CvType is float with only one channel
                    //the method Mat.get() return a double array 
                    //with only one position, [0].
                    a = (float) im.get(rcnt, ccnt)[0];
                    b = (float) im.get(rcnt+1, ccnt)[0];
                    c = (float) ((a+b)*0.707);
                    d = (float) ((a-b)*0.707);
                    int _rcnt= rcnt/2;
                    im1.put(_rcnt, ccnt, c);
                    im2.put(_rcnt, ccnt, d);

                }
            }

            for (int rcnt = 0; rcnt < im.rows()/2; rcnt++) {
                for (int ccnt = 0; ccnt < im.cols() - 2; ccnt+=2) {
                    a = (float) im1.get(rcnt, ccnt)[0];
                    b = (float) im1.get(rcnt, ccnt+1)[0];
                    c = (float) ((a+b)*0.707);
                    d = (float) ((a-b)*0.707);
                    int _ccnt = ccnt/2;
                    im3.put(rcnt, _ccnt, c);
                    im4.put(rcnt, _ccnt, d);
                }
            }

            for (int rcnt = 0; rcnt < im.rows()/2; rcnt++) {
                for (int ccnt = 0; ccnt < im.cols() - 2; ccnt+=2) {
                    a = (float) im2.get(rcnt, ccnt)[0];
                    b = (float) im2.get(rcnt, ccnt+1)[0];
                    c = (float) ((a+b)*0.707);
                    d = (float) ((a-b)*0.707);
                    int _ccnt = ccnt/2;
                    im5.put(rcnt, _ccnt, c);
                    im6.put(rcnt, _ccnt, d);
                }
            }

            imr = Mat.zeros(im.rows(), im.cols(), CvType.CV_32F);//imr = Mat.zeros(512, 512, CvType.CV_32F);
            imd = Mat.zeros(512, 512, CvType.CV_32F);
            im3.copyTo(imd.adjustROI(0, 0, 256, 256));
            im4.copyTo(imd.adjustROI(0, 255, 256, 256));
            im5.copyTo(imd.adjustROI(255, 0, 256, 256));
            im6.copyTo(imd.adjustROI(255, 255, 256, 256));





        }catch(Exception ex){
            System.err.println(ex.getLocalizedMessage());
            ex.printStackTrace();
        }
    }

    private void applyHaarReverse(){
        // ------------------- Reconstruction ------------------- 
                im11 = Mat.zeros(im.rows()/2, im.cols(), CvType.CV_32F);
                im12 = Mat.zeros(im.rows()/2, im.cols(), CvType.CV_32F);
                im13 = Mat.zeros(im.rows()/2, im.cols(), CvType.CV_32F);
                im14 = Mat.zeros(im.rows()/2, im.cols(), CvType.CV_32F);

                for (int rcnt = 0; rcnt < im.rows()/2; rcnt++) {
                    for (int ccnt = 0; ccnt < im.cols()/2; ccnt++) {
                        int _ccnt  = ccnt*2;
                        im11.put(rcnt, _ccnt, im3.get(rcnt, ccnt));
                        im12.put(rcnt, _ccnt, im4.get(rcnt, ccnt));
                        im13.put(rcnt, _ccnt, im5.get(rcnt, ccnt));
                        im14.put(rcnt, _ccnt, im6.get(rcnt, ccnt));
                    }
                }

                for (int rcnt = 0; rcnt < im.rows()/2; rcnt++) {
                    for (int ccnt = 0; ccnt < im.cols() - 2; ccnt+=2) {
                        a = (float) im11.get(rcnt, ccnt)[0];
                        b = (float) im12.get(rcnt, ccnt)[0];
                        c = (float) ((a+b)*0.707);
                        im11.put(rcnt, ccnt, c);
                        d = (float) ((a-b)*0.707);
                        im11.put(rcnt, ccnt+1, d);

                        a = (float) im13.get(rcnt, ccnt)[0];
                        b = (float) im14.get(rcnt, ccnt)[0];
                        c = (float) ((a+b)*0.707);
                        im13.put(rcnt, ccnt, c);
                        d = (float) ((a-b)*0.707);
                        im13.put(rcnt, ccnt+1, d);
                    }
                }

                temp = Mat.zeros(im.rows(), im.cols(), CvType.CV_32F);

                for (int rcnt = 0; rcnt < im.rows()/2; rcnt++) {
                    for (int ccnt = 0; ccnt < im.cols(); ccnt++) {
                        int _rcnt = rcnt*2;
                        imr.put(_rcnt, ccnt, im11.get(rcnt, ccnt));
                        temp.put(_rcnt, ccnt, im13.get(rcnt, ccnt));
                    }
                }


                for (int rcnt = 0; rcnt < im.rows()-2; rcnt+=2) {
                    for (int ccnt = 0; ccnt < im.cols(); ccnt++) {
                        a = (float) imr.get(rcnt, ccnt)[0];
                        b = (float) temp.get(rcnt, ccnt)[0];
                        c = (float) ((a+b)*0.707);
                        imr.put(rcnt, ccnt, c);
                        d = (float) ((a-b)*0.707);
                        imr.put(rcnt+1, ccnt, d);
                    }
                }

    }
}

希望有用。

Answer 4

建议： 我不建议你从头实现 dwt，很难满足你的要求。 如果你在工作中确实需要它的cpp版本，我推荐你试试PyWavelets，它的dwt的基本功能是用C实现的，所以你可以很容易地把它移植到你的程序中。 否则，你也可以先用python验证你的想法，而不是冒险付出无用的努力。

---

Here is one of my implementation of dwt 支持多种小波滤波器，效果好但效果不好。随着层次的增加，重构越来越模糊。

如果你想改变小波滤波器，你可以使用matlab（wfilters(wname)）或者从PyWavelets's source code中挑选，它也会从这些系数中得到the rule to get totall 4 filter。