1. 程式人生 > >opencv+小波變換

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

這是原圖

下面上一張效果圖