Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Wavelet transform in openCV

did someone tried to implement DWT in opencv or in C++? I saw older posts on this subject and i didn't find them useful for me, because I need a approximation coefficient and details as a result of wavelet transformation.

I tried to add this to my project but it's not working as well as planned.

And this is to simple, because as a result parameters i need approximation coefficient and details:

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;
}

So can someone help me find dwt implemented in C++ or point me how to extract from above code coefficients. Thanks

like image 938
la lluvia Avatar asked Nov 19 '13 12:11

la lluvia


People also ask

What is a wavelet transform in image processing?

Wavelet transform is one of the important methods of compressing image data so that it takes up less memory. Wavelet based compression techniques have advantages such as multi-resolution, scalability and tolerable degradation over other techniques. Read more. Science.

What is wavelet transform used for?

The wavelet transform (WT) can be used to analyze signals in time–frequency space and reduce noise, while retaining the important components in the original signals. In the past 20 years, WT has become a very effective tool in signal processing.

What is meant by wavelet transform?

Wavelet transform offers a generalization of STFT. From a signal theory point of view, similar to DFT and STFT, wavelet transform can be viewed as the projection of a signal into a set of basis functions named wavelets. Such basis functions offer localization in the frequency domain.


3 Answers

Here is direct and inverse Haar Wavelet transform (used for filtering):

#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);
}
like image 180
Andrey Smorodov Avatar answered Oct 13 '22 16:10

Andrey Smorodov


Here is another implementation of Wavelet transform in OpenCV from Mahavir:

            #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;
            }

Hope someone finds it useful!

like image 35
la lluvia Avatar answered Oct 13 '22 14:10

la lluvia


I see that there's very few code examples for wavelet in java, especially if you're using openCV. I had to use wavelet in java with openCV and I used the C code from @la luvia and converted to java.

There was a lot of trouble while translating the code, because it had a lot of diferences in the openCV methods and ways of using it. This book helped me a lot in the process too.

I hope this code and the book give some perspective of how to use the lib and few diferences between the C and Java.

Here's the code:

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);
                    }
                }

    }
}

Hope it's useful.

like image 1
RochaRF Avatar answered Oct 13 '22 14:10

RochaRF