最近整理了基于Harr小波分解的图像融合技术实现可见光图像和红外图像的简单融合，本文侧重于图像融合，并不注重于Harr小波变换，所以就只做了二层小波分解和重构（分解层数越多，噪声处理起来也越麻烦），下面简要给出理论分析和代码。
一.Harr小波变换
小波变换的基本思想是想用一组小波函数或者基函数表示一个图像信号，其中Harr小波基函数是其中最简单的基函数，数学原理不做推导，有兴趣的同学可以学习小波变换的有关知识。
1.一维Harr小波分解
为了理解什么是一维Harr小波分解，下面用一个具体的例子来说明Harr小波分解的过程。假设有一幅只有8个像素的一维图像，对应像素为[2,4,6,8,10,12,14,16]。
（1）求均值：计算相邻像素对的平均值，得到一幅分辨率变为原来一半的图像[3,7,11,15],该部分储存是图像的整体信息，属于图像的高频分量。
（2）求差值：通过第一步求均值保留了图像的整体信息，但丢失了图像的细节信息，所以还要记录图像的细节信息以便重构时能够恢复图像的全部信息。通过计算原像素[2,4,6,8,10,12,14,16]相邻像素值之间的差值的1/2得到[-1,-1,-1,-1],该部分记录了图像的细节信息，属于低频分量。
（3）将（1）（2）得到的原始图像1/2分辨率的图像组合成新的图像[3,7,11,15,-1,-1,-1,-1],这就是Harr小波第一次分解的结果。如果还想做进一步的分解，可将第一次分解结果作为原始像素值重复第（1）（2）步即可。
2.二维Harr小波分解
有图像处理经验的同学都知道一幅图像可以看作一个二维矩阵，下面就从二维的角度讲解Harr小波分解。假设有一副8x8的灰度图像，表示为：
[64 2 3 61 60 6 7 57
9 55 54 12 13 51 50 16
17 47 46 20 21 43 42 24
40 26 27 37 36 30 31 33
32 34 35 29 28 38 39 25
41 23 22 44 45 19 18 48
49 15 14 52 53 11 10 56
8 58 59 5 4 62 63 1]
（1）经过行分解后得到：
[33 32 33 32 31 -29 27 -25
32 33 32 33 -23 21 -19 17
32 33 32 33 -15 13 -11 9
33 32 33 32 7 -5 3 -1
33 32 33 32 -1 3 -5 7
32 33 32 33 9 -11 13 -15
32 33 32 33 17 -19 21 -23
33 32 33 32 -25 27 -29 31
]
（2）在（1）的基础上经过列分解后得到：
[32.5 32.5 32.5 32.5 4 -4 4 -4
32.5 32.5 32.5 32.5 -4 4 -4 4
32.5 32.5 32.5 32.5 4 -4 4 -4
32.5 32.5 32.5 32.5 -4 4 -4 4
0.5 -0.5 0.5 -0.5 27 -25 23 -21
-0.5 0.5 -0.5 0.5 -11 9 -7 5
0.5 -0.5 0.5 -0.5 -5 7 -9 11
-0.5 0.5 -0.5 0.5 21 -23 25 -27
]。这就是一次二维Harr小波分解的结果，其中高频分量1个（左上角4x4的部分），低频分量3个（剩余部分），若还想多次分解，可重复（1）（2）过程。
二.代码实现
话不多说，直接上代码：
class WaveTransform
{
public:
Mat WDT(const Mat &_src, const string _wname, const int _level);//小波分解
Mat IWDT(const Mat &_src, const string _wname, const int _level);//小波重构
void wavelet_D(const string _wname, Mat &_lowFilter, Mat &_highFilter);//分解包
void wavelet_R(const string _wname, Mat &_lowFilter, Mat &_highFilter);//重构包
Mat waveletDecompose(const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter);
Mat waveletReconstruct(const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter);
};
Mat WaveTransform::WDT(const Mat &_src, const string _wname, const int _level)
{
Mat_<float> src = Mat_<float>(_src);
Mat dst = Mat::zeros(src.rows, src.cols, src.type());
int row = src.rows;
int col = src.cols;
//高通低通滤波器
Mat lowFilter;
Mat highFilter;
wavelet_D(_wname, lowFilter, highFilter);
//小波变换
int t = 1;
while (t <= _level)
{
//先进行 行小波变换
//#pragma omp parallel for
for (int i = 0;i<row;i++)
{
//取出src中要处理的数据的一行
Mat oneRow = Mat::zeros(1, col, src.type());
for (int j = 0;j<col;j++)
{
oneRow.at<float>(0, j) = src.at<float>(i, j);
}
oneRow = waveletDecompose(oneRow, lowFilter, highFilter);
for (int j = 0;j<col;j++)
{
dst.at<float>(i, j) = oneRow.at<float>(0, j);
}
}
//小波列变换
//#pragma omp parallel for
for (int j = 0;j<col;j++)
{
Mat oneCol = Mat::zeros(row, 1, src.type());
for (int i = 0;i<row;i++)
{
oneCol.at<float>(i, 0) = dst.at<float>(i, j);//dst,not src
}
oneCol = (waveletDecompose(oneCol.t(), lowFilter, highFilter)).t();
for (int i = 0;i<row;i++)
{
dst.at<float>(i, j) = oneCol.at<float>(i, 0);
}
}
//更新
row /= 2;
col /= 2;
t++;
src = dst;
}
return dst;
}
Mat WaveTransform::IWDT(const Mat &_src, const string _wname, const int _level)
{
Mat src = Mat_<float>(_src);
Mat dst;
src.copyTo(dst);
int N = src.rows;
int D = src.cols;
//高低通滤波器
Mat lowFilter;
Mat highFilter;
wavelet_R(_wname, lowFilter, highFilter);
//小波变换
int t = 1;
int row = N / std::pow(2., _level - 1);
int col = D / std::pow(2., _level - 1);
while (row <= N && col <= D)
//while(t<=_level)
{
//列逆变换
for (int j = 0;j<col;j++)
{
Mat oneCol = Mat::zeros(row, 1, src.type());
for (int i = 0;i<row;i++)
{
oneCol.at<float>(i, 0) = src.at<float>(i, j);
}
oneCol = (waveletReconstruct(oneCol.t(), lowFilter, highFilter)).t();
for (int i = 0;i<row;i++)
{
dst.at<float>(i, j) = oneCol.at<float>(i, 0);
}
}
//行逆变换
for (int i = 0;i<row;i++)
{
Mat oneRow = Mat::zeros(1, col, src.type());
for (int j = 0;j<col;j++)
{
oneRow.at<float>(0, j) = dst.at<float>(i, j);
}
oneRow = waveletReconstruct(oneRow, lowFilter, highFilter);
for (int j = 0;j<col;j++)
{
dst.at<float>(i, j) = oneRow.at<float>(0, j);
}
}
row *= 2;
col *= 2;
t++;
src = dst;
}
return dst;
}
void WaveTransform::wavelet_D(const string _wname, Mat &_lowFilter, Mat &_highFilter)
{
if (_wname == "haar" || _wname == "db1")
{
int N = 2;
_lowFilter = Mat::zeros(1, N, CV_32F);
_highFilter = Mat::zeros(1, N, CV_32F);
_lowFilter.at<float>(0, 0) = 1 / sqrtf(N);
_lowFilter.at<float>(0, 1) = 1 / sqrtf(N);
_highFilter.at<float>(0, 0) = -1 / sqrtf(N);
_highFilter.at<float>(0, 1) = 1 / sqrtf(N);
}
else if (_wname == "sym2")
{
int N = 4;
float h[] = { -0.4830, 0.8365, -0.2241, -0.1294 };
float l[] = { -0.1294, 0.2241, 0.8365, 0.4830 };
_lowFilter = Mat::zeros(1, N, CV_32F);
_highFilter = Mat::zeros(1, N, CV_32F);
for (int i = 0;i<N;i++)
{
_lowFilter.at<float>(0, i) = l[i];
_highFilter.at<float>(0, i) = h[i];
}
}
}
void WaveTransform::wavelet_R(const string _wname, Mat &_lowFilter, Mat &_highFilter)
{
if (_wname == "haar" || _wname == "db1")
{
int N = 2;
_lowFilter = Mat::zeros(1, N, CV_32F);
_highFilter = Mat::zeros(1, N, CV_32F);
_lowFilter.at<float>(0, 0) = 1 / sqrtf(N);
_lowFilter.at<float>(0, 1) = 1 / sqrtf(N);
_highFilter.at<float>(0, 0) = 1 / sqrtf(N);
_highFilter.at<float>(0, 1) = -1 / sqrtf(N);
}
else if (_wname == "sym2")
{
int N = 4;
float h[] = { -0.1294,-0.2241,0.8365,-0.4830 };
float l[] = { 0.4830, 0.8365, 0.2241, -0.1294 };
_lowFilter = Mat::zeros(1, N, CV_32F);
_highFilter = Mat::zeros(1, N, CV_32F);
for (int i = 0;i<N;i++)
{
_lowFilter.at<float>(0, i) = l[i];
_highFilter.at<float>(0, i) = h[i];
}
}
}
Mat WaveTransform::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);
//下采样
//数据拼接
for (int i = 0, j = 1;i<D / 2;i++, j += 2)
{
src.at<float>(0, i) = dst1.at<float>(0, j);//lowFilter
src.at<float>(0, i + D / 2) = dst2.at<float>(0, j);//highFilter
}
return src;
}
Mat WaveTransform::waveletReconstruct(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);
/// 插值;
Mat Up1 = Mat::zeros(1, D, src.type());
Mat Up2 = Mat::zeros(1, D, src.type());
for (int i = 0, cnt = 0; i < D / 2; i++, cnt += 2)
{
Up1.at<float>(0, cnt) = src.at<float>(0, i); ///< 前一半
Up2.at<float>(0, cnt) = src.at<float>(0, i + D / 2); ///< 后一半
}
/// 前一半低通，后一半高通
Mat dst1 = Mat::zeros(1, D, src.type());
Mat dst2 = Mat::zeros(1, D, src.type());
filter2D(Up1, dst1, -1, lowFilter);
filter2D(Up2, dst2, -1, highFilter);
/// 结果相加
dst1 = dst1 + dst2;
return dst1;
}
void Coef_Fusion(Mat Visible, Mat Infrared, Mat &dst) {
assert(Visible.rows == Infrared.rows&&Visible.cols == Infrared.cols);
Mat Visible_gray,Infrared_gray;
cvtColor(Visible, Visible_gray, CV_RGB2GRAY);
normalize(Visible_gray, Visible_gray, 0, 255, CV_MINMAX);
cvtColor(Infrared, Infrared_gray, CV_RGB2GRAY);
normalize(Infrared_gray, Infrared_gray, 0, 255, CV_MINMAX);
Visible_gray.convertTo(Visible_gray, CV_32F);
Infrared_gray.convertTo(Infrared_gray, CV_32F);
WaveTransform m_waveTransform;
const int level = 2;
Mat visible_dec = m_waveTransform.WDT(Visible_gray, "haar", level);
Mat Infrared_dec = m_waveTransform.WDT(Infrared_gray, "haar", level);
int row = Visible.rows;
int col = Visible.cols;
Mat dec = Mat(Visible.size(), CV_32FC1);
//融合规则：高频部分采用加权平均的方法，低频部分采用模值取大的方法
int halfRow = row / (2 * level);
int halfCol = col / (2 * level);
for (int i = 0;i < row;i++) {
for (int j = 0;j < col;j++) {
if (i < halfRow&&j < halfCol) {
dec.at<float>(i, j) = (visible_dec.at<float>(i, j) + Infrared_dec.at<float>(i, j)) / 2;
}
else {
float p = abs(visible_dec.at<float>(i, j));
float q = abs(Infrared_dec.at<float>(i, j));
if (p > q) {
dec.at<float>(i, j) = visible_dec.at<float>(i, j);
}
else {
dec.at<float>(i, j) = Infrared_dec.at<float>(i, j);
}
}
}
}
dst = m_waveTransform.IWDT(dec, "haar", level);
dst.convertTo(dst, CV_8UC1);
normalize(dst, dst, 0, 255, CV_MINMAX);
}
int main(int argc, char *argv[])
{
Mat Visible= imread("D:/VS2015/image/Visible.bmp"); //右图
Mat Infrared = imread("D:/VS2015/image/Infrared.bmp"); //左图
if (Visible.empty() || Infrared.empty()) {
cout << "could not load the image..." << endl;
return -1;
}
imshow("可见光图像", Visible);
imshow("红外图像", Infrared);
Mat dst;
Coef_Fusion(Visible,Infrared,dst);
if (dst.data) {
imshow("融合结果", dst);
}
waitKey();
destroyAllWindows();
return 0;
}
函数说明：
Mat WDT(const Mat &_src, const string _wname, const int _level);//小波分解（src：输入图像，_wname:小波基，_level:分解层数）
Mat IWDT(const Mat &_src, const string _wname, const int _level);//小波重构（src：输入图像，_wname:小波基，_level:分解层数）
void Coef_Fusion(Mat Visible, Mat Infrared, Mat &dst);//融合函数（Visible：可见光图像，Infrared：红外图像，dst：融合结果）
注意：该方法只能处理灰度图像，也就是融合是基于灰度级的，结果也是灰度图片。
三.实验结果
四：参考
我的小波分解和重构也是参考别人的代码，这里把原作者的博客链接放在这里，大家可以去观看学习。
链接：https://blog.csdn.net/u014426939/article/details/82867874
以上就是所有分析和结果啦，如果有什么问题欢迎大家和我交流，非常乐意和大家交流计算机视觉相关的知识！
扫一扫
906
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
