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();
}
這是原圖
下面上一張效果圖