参考博客

有关小波的几个术语及常见的小波基介绍_jbb0523的博客-CSDN博客_常用小波基

python小波图像融合_koervcor的博客-CSDN博客_python小波图像融合

医学图像处理案例（十四）——基于小波变换的图像融合 - 腾讯云开发者社区-腾讯云 (tencent.com)

图像融合（六）-- 小波融合 - silenceer - 博客园 (cnblogs.com)

常见小波函数：Haar、Daubechies、Biorthogonal、Coiflets、Symlets、Morlet、Mexican Hat、Meyer、Gaus、Dmeyer、ReverseBior、Cgau、Cmor、Fbsp、Shan.

图像融合与小波变换&融合 概述

1、图像融合概述

图像融合（Image Fusion）是指将多源信道所采集到的关于同一目标的图像数据经过图像处理和计算机技术等，最大限度的提取各自信道中的有利信息，最后综合成高质量的图像，以提高图像信息的利用率、改善计算机解译精度和可靠性、提升原始图像的空间分辨率和光谱分辨率，利于监测。

2、小波变换特点介绍

小波变换的固有特性使其在图像处理中有如下优点：

  • 完善的重构能力，保证信号在分解过程中没有信息损失和冗余信息；
  • 把图像分解成低频图像和细节（高频）图像的组合，分别代表了图像的不同结构，因此容易提取原始图像的结构信息和细节信息；
  • 小波分析提供了与人类视觉系统方向相吻合的选择性图像。

一般图像融合的小波分解采用离散小波变换(Discrete Wavelet Transform, DWT)。DWT的函数基由一个称为母小波或分析小波的单一函数通过膨胀和平移获得。因而，DWT同时具有时域和频域分析能力，与一般的金字塔分解相比，DWT图像分解具有以下优势：

  • 1）具有方向性，在提取图像低频信息的同时，还可获得了水平、垂直和对角三个方向的高频信息；
  • 2）通过合理的选择母小波，可使DWT在压缩噪声的同时更有效的提取纹理、边缘等显著信息；
  • 3）金字塔分解各尺度之间具有信息的相关性，而DWT在不同尺度上具有更高的独立性。

3、基于小波变换的图像融合

DWT 融合算法基本思想：首先对源图像进行小波变换，然后按照一定规则对变换系数进行合并；最后对合并后的系数进行小波逆变换得到融合图像。

3.1、小波分解原理简介

（1）小波的简单计算原理

　　[x0，x1，x2，x3]=[90，70，100，70] 为达到压缩 我们可取 (x0+x1)/2  (x0-x1)/2 来代表 x0,x1  这样 [90,70] 可表示为 [80,10] 80即平均数 10是小范围波动数（可想象出一种波的形状） [90,70] --〉[80,10] , [100,70] --〉 [85,15] 可以想象80 和85 都是局部的平均值反映大的总体的状态，是变化相对缓慢的值，可以认为他们是低频部分的值。 而10、15是小范围波动的值局部变换较快，可以认为他们是高频部分的值。

  1. FIRST：把[90,70,100,70] 写成 [80,85,10,15] 即把低频部分写在一起（记频率L） 高频部分写在一起（H) ；
  2. SECOND：而[80,85] 又可经同样的变换--> [82.5, -2.5] 这样 82.5表示更低频的信息(记频率LL) -2.5则表示了频率L上的波动； 
  3. 最后[90,70,100,70] --〉[82.5, -2.5, 10, 15] 这样信息就可被压缩了（数字范围小了）。

　　现在再来扩展一下  [90,70]---> [80,10] 写成矩阵 [90,70] * [ [1/2, 1/2] [1/2 ,-1/2] ] , 矩阵 [1,1；1,-1]/2 为haar转换矩阵。

　　如果是[90,70,100,70]，第一步就可以写成矩阵M1：[0.5,0,0.5,0; 0.5,0,-0.5,0; 0,0.5,0,0.5;0,0.5,0,-0.5]；第二步只对低频L操作，高频不变，可写成M2：[1/2,  1/2, 0, 0; 1/2, -1/2, 0, 0; 0,  0,  1, 0 ;0,  0,  0, 1]。另M= M1*M2，可得到4*4的点阵操作。

　　第一步运算后原图像缩小至左边一半了，右边的是对应波动信息；

　　第二步运算后图像又缩小至左边一半了，对应波动信息。

　　对一幅图像先进行行变化，在进行列变化，那么就是小波变化了。

 LL：水平低频，垂直低频

LH：水平低频，垂直高频

HL：水平高频，垂直低频

HH：水平高频，垂直高频

其中，L表示低频，H表示高频，下标1、2表示一级或二级分解。在每一分解层上，图像均被分解为LL，LH，HH和HL四个频带，下一层的分解仅对低频分量LL进行分解。这四个子图像中的每一个都是由原图与一个小波基函数的内积后，再经过在x和y方向都进行2倍的间隔采样而生成的，这是正变换，也就是图像的分解；逆变换，也就是图像的重建，是通过图像的增频采样和卷积来实现的。这里有个问题经过处理后，数据或超出255或者出现负数，需要将其归一化到0-255之间，方可显示图像。这里介绍的只是简单的小波计算，小波计算的而不同就在于选取不同的小波系数，一般有haar小波，sym2小波等。

3.2、融合规则

规则一：系数绝对值较大法

该融合规则适合高频成分比较丰富，亮度、对比度比较高的源图像，否则在融合图像中只保留一幅源图像的特征，其他的特征被覆盖。小波变换的实际作用是对信号解相关，并将信号的全部信息集中到一部分具有大幅值的小波系数中。这些大的小波系数含有的能量远比小系数含有的能量大，从而在信号的重构中，大的系数比小的系数更重要。

规则二：加权平均法

权重系数可调，适用范围广，可消除部分噪声，源图像信息损失较少，但会造成图像对比度的下降，需要增强图像灰度。

规则三：局部方差准则

设A(x,y)和B(x,y)分别为高频子图像数据值，F(x,y)为相应高频子图像融合值，将A(x,y)和B(x,y)分成若干个M×N子块图像。对每个子块图像进行数值分布统计，计算其方差。确定A和B图像每个子块图像加权系数K1和K2。如果A图像子块方差大于B图像子块方差,则K1≥K2，否则K1<K2。确定每个子块图像的数据融合数值为：F(i,j)=K1A(i,j)+K2B(i,j)。

4、基于小波变换的图像融合代码实现

matlab和python版本代码来融合红外和可见光图像，融合策略是低频图像采用平均值法，高频图像采用最大值法。python版本中需要用到PyWavelets库 [ pip install PyWavelets ] 。

python版本代码

import pywt
import cv2
import numpy as np
# This function does the coefficient fusing according to the fusion method
def fuseCoeff(cooef1, cooef2, method):
    if (method == 'mean'):
        cooef = (cooef1 + cooef2) / 2
    elif (method == 'min'):
        cooef = np.minimum(cooef1, cooef2)
    elif (method == 'max'):
        cooef = np.maximum(cooef1, cooef2)
    return cooef

# Params
FUSION_METHOD = 'mean'  # Can be 'min' || 'max || anything you choose according theory
FUSION_METHOD1 = 'max'
# Read the two image
I1 = cv2.imread('IR1.png', 0)
I2 = cv2.imread('VIS1.png', 0)
# First: Do wavelet transform on each image
wavelet = 'db2'
cooef1 = pywt.wavedec2(I1[:, :], wavelet, level=1)
cooef2 = pywt.wavedec2(I2[:, :], wavelet, level=1)
# Second: for each level in both image do the fusion according to the desire option
fusedCooef = []
for i in range(len(cooef1)):
    # The first values in each decomposition is the apprximation values of the top level
    if (i == 0):
        fusedCooef.append(fuseCoeff(cooef1[0], cooef2[0], FUSION_METHOD))
    else:
        # For the rest of the levels we have tupels with 3 coeeficents
        c1 = fuseCoeff(cooef1[i][0], cooef2[i][0], FUSION_METHOD1)
        c2 = fuseCoeff(cooef1[i][1], cooef2[i][1], FUSION_METHOD1)
        c3 = fuseCoeff(cooef1[i][2], cooef2[i][2], FUSION_METHOD1)
        fusedCooef.append((c1, c2, c3))
# Third: After we fused the cooefficent we nned to transfor back to get the image
fusedImage = pywt.waverec2(fusedCooef, wavelet)
# Forth: normmalize values to be in uint8
fusedImage1 = np.multiply(np.divide(fusedImage - np.min(fusedImage), (np.max(fusedImage) - np.min(fusedImage))), 255)
fusedImage1 = fusedImage1.astype(np.uint8)
# Fith: Show image
cv2.imwrite("win.bmp", fusedImage1)

matlab版本代码

close all
clear
clc
%% Read Images
% the size of images must be equal
a=imread('IR1.png');
b=imread('VIS1.png');
%%   Wavelet Transform 
[a1,b1,c1,d1]=dwt2(a,'db2');
[a2,b2,c2,d2]=dwt2(b,'db2');
[k1,k2]=size(a1);
%% Fusion Rules
%% Average Rule
for i=1:k1
    for j=1:k2
        a3(i,j)=(a1(i,j)+a2(i,j))/2;
   end
end
%% Max Rule
for i=1:k1
    for j=1:k2
        b3(i,j)=max(b1(i,j),b2(i,j));
        c3(i,j)=max(c1(i,j),c2(i,j));
        d3(i,j)=max(d1(i,j),d2(i,j));
    end
end
%% Inverse Wavelet Transform 
c=idwt2(a3,b3,c3,d3,'db2');
imshow(a)
title('First Image')
figure,imshow(b)
title('Second Image')
figure,imshow(c,[])
title('Fused Image')
fusedimage=(c-min(c(:)))/(max(c(:))-min(c(:)))*255;
fusedimage1=uint8(fusedimage);
imwrite(fusedimage1,'fusedimage.bmp')

效果展示

左上：红外图像；左下：可见光图像；右上matlab融合，右下python融合结果。

小波变换融合实例

python小波图像融合_koervcor的博客-CSDN博客_python小波图像融合

对两张原图进行小波融合，得到一个清晰图。 

两张原图，左图左边模糊，右图右边模糊。

实施流程

1. 对两张图像分别对水平方向和垂直方向进行小波分解，效果如下（L：代表低频，H：代表高频）

2. 对两张图像高频进行绝对值最大法，低频进行局部方差准则的图像进行融合，（低频：加权平均，高频：绝对值取大。）效果如下：

3. 图像重构（小波逆变化），得到融合后的图像，效果如下：

python实现代码

# -*- coding:utf-8 -*-
# @auth hdy
# 哈尔小波变换
import cv2
import numpy as np 
 
imgA = cv2.imread("a.tif") #载入图片A
imgB = cv2.imread("b.tif") #载入图片B
heigh, wide, channel = imgA.shape #获取图像的高、宽、通道数
###############################
#临时变量、存储哈尔小波处理后的数据#
tempA1 = []                   #
tempA2 = []                   #
tempB1 = []                   #
tempB2 = []                   #
###############################
 
waveImgA = np.zeros((heigh, wide, channel), np.float32) #存储A图片小波处理后数据的变量
waveImgB = np.zeros((heigh, wide, channel), np.float32) #存储B图片小波处理后数据的变量
 
#水平方向的哈尔小波处理，对图片的B、G、R三个通道分别遍历进行
for c in range(channel):
	for x in range(heigh):
		for y in range(0,wide,2):
			tempA1.append((float(imgA[x,y,c]) + float(imgA[x,y+1,c]))/2) #将图片A小波处理后的低频存储在tempA1中
			tempA2.append((float(imgA[x,y,c]) + float(imgA[x,y+1,c]))/2 - float(imgA[x,y,c])) #将图片A小波处理后的高频存储在tempA2中
			tempB1.append((float(imgB[x,y,c]) + float(imgB[x,y+1,c]))/2) #将图片B小波处理后的低频存储在tempB1中
			tempB2.append((float(imgB[x,y,c]) + float(imgB[x,y+1,c]))/2 - float(imgB[x,y,c])) #将图片B小波处理后的高频存储在tempB2中
		tempA1 = tempA1 + tempA2 #小波处理完图片A每一个水平方向数据统一保存在tempA1中
		tempB1 = tempB1 + tempB2 #小波处理完图片B每一个水平方向数据统一保存在tempB1中
		for i in range(len(tempA1)):
			waveImgA[x,i,c] = tempA1[i] #图片A水平方向前半段存储低频，后半段存储高频
			waveImgB[x,i,c] = tempB1[i] #图片B水平方向前半段存储低频，后半段存储高频
		tempA1 = [] #当前水平方向数据处理完之后，临时变量重置
		tempA2 = [] 
		tempB1 = []
		tempB2 = []
 
#垂直方向哈尔小波处理，与水平方向同理
for c in range(channel):
	for y in range(wide):
		for x in range(0,heigh-1,2):
			tempA1.append((float(waveImgA[x,y,c]) + float(waveImgA[x+1,y,c]))/2)
			tempA2.append((float(waveImgA[x,y,c]) + float(waveImgA[x+1,y,c]))/2 - float(waveImgA[x,y,c]))
			tempB1.append((float(waveImgB[x,y,c]) + float(waveImgB[x+1,y,c]))/2)
			tempB2.append((float(waveImgB[x,y,c]) + float(waveImgB[x+1,y,c]))/2 - float(waveImgB[x,y,c]))
		tempA1 = tempA1 + tempA2
		tempB1 = tempB1 + tempB2
		for i in range(len(tempA1)):
			waveImgA[i,y,c] = tempA1[i]
			waveImgB[i,y,c] = tempB1[i]
		tempA1 = []
		tempA2 = []
		tempB1 = []
		tempB2 = []
 
#求以x,y为中心的5x5矩阵的方差，  “//”在python3中表示整除，没有小数，“/”在python3中会有小数，  python2中“/”即可，“//”也行都表示整除
varImgA = np.zeros((heigh//2, wide//2, channel), np.float32) #将图像A中低频数据求方差之后存储的变量
varImgB = np.zeros((heigh//2, wide//2, channel), np.float32) #将图像B中低频数据求方差之后存储的变量
for c in range(channel):
	for x in range(heigh//2):
		for y in range(wide//2):
			#############################
			#对图片边界(或临近)的像素点进行处理
			if x - 3    <   0:
				up      =   0
			else:
				up      =   x - 3
			if x + 3    >   heigh//2:
				down    =   heigh//2
			else:
				down    =   x + 3
			if y - 3    <   0:
				left    =   0
			else:
				left    =   y - 3
			if y + 3    >   wide//2:
				right   =   wide//2
			else:
				right   =   y + 3
			#############################
			meanA, varA = cv2.meanStdDev(waveImgA[up:down,left:right,c]) #求图片A以x,y为中心的5x5矩阵的方差，mean表示平均值，var表示方差
			meanB, varB = cv2.meanStdDev(waveImgB[up:down,left:right,c]) #求图片B以x,y为中心的5x5矩阵的方差，
 
			varImgA[x,y,c] = varA #将图片A对应位置像素的方差存储在变量中
			varImgB[x,y,c] = varB #将图片B对应位置像素的方差存储在变量中
 
#求两图的权重
weightImgA = np.zeros((heigh//2, wide//2, channel), np.float32) #图像A存储权重的变量
weightImgB = np.zeros((heigh//2, wide//2, channel), np.float32) #图像B存储权重的变量
for c in range(channel):
	for x in range(heigh//2):
		for y in range(wide//2):
			weightImgA[x,y,c] = varImgA[x,y,c] / (varImgA[x,y,c]+varImgB[x,y,c]+0.00000001) #分别求得图片A与图片B的权重
			weightImgB[x,y,c] = varImgB[x,y,c] / (varImgA[x,y,c]+varImgB[x,y,c]+0.00000001) #“0.00000001”防止零除
 
#进行融合，高频————系数绝对值最大化，低频————局部方差准则
reImgA = np.zeros((heigh, wide, channel), np.float32) #图像融合后的存储数据的变量
reImgB = np.zeros((heigh, wide, channel), np.float32) #临时变量
for c in range(channel):
	for x in range(heigh):
		for y in range(wide):
			if x < heigh//2 and y < wide//2:
				reImgA[x,y,c] = weightImgA[x,y,c]*waveImgA[x,y,c] + weightImgB[x,y,c]*waveImgB[x,y,c] #对两张图片低频的地方进行权值融合数据
			else:
				reImgA[x,y,c] = waveImgA[x,y,c] if abs(waveImgA[x,y,c]) >= abs(waveImgB[x,y,c]) else waveImgB[x,y,c] #对两张图片高频的进行绝对值系数最大规则融合
 
#由于是先进行水平方向小波处理，因此重构是先进行垂直方向
#垂直方向进行重构
for c in range(channel):
	for y in range(wide):
		for x in range(heigh):
			if x%2 == 0:
				reImgB[x,y,c] = reImgA[x//2,y,c] - reImgA[x//2 + heigh//2,y,c] #根据哈尔小波原理，将重构后的数据存储在临时变量中
			else:
				reImgB[x,y,c] = reImgA[x//2,y,c] + reImgA[x//2 + heigh//2,y,c] #图片的前半段是低频后半段是高频，除以2余数为0相减，不为0相加
 
#水平方向进行重构，与垂直方向同理
for c in range(channel):
	for x in range(heigh):
		for y in range(wide):
			if y%2 ==0:
				reImgA[x,y,c] = reImgB[x,y//2,c] - reImgB[x,y//2 + wide//2,c]
			else:
				reImgA[x,y,c] = reImgB[x,y//2,c] + reImgB[x,y//2 + wide//2,c]
#限制图像的范围(0-255),若不限制，根据np.astype(np.uint8)的规则，会对图片产生噪声
reImgA[reImgA[:, :, :] < 0] = 0
reImgA[reImgA[:, :, :] > 255] = 255
 
cv2.imshow("reImg", reImgA.astype(np.uint8))
cv2.waitKey(0)
cv2.destroyAllWindows()

附录

有关小波的几个术语及常见的小波基介绍_jbb0523的博客-CSDN博客_常用小波基

一、小波基选择标准

小波变换不同于傅里叶变换，根据小波母函数的不同，小波变换的结果也不尽相同。现实中到底选择使用哪一种小波的标准一般有以下几点：

1、支撑长度

小波函数Ψ(t)、Ψ(ω)、尺度函数φ(t)和φ(ω)的支撑区间，是当时间或频率趋向于无穷大时，Ψ(t)、Ψ(ω)、φ(t)和φ(ω)从一个有限值收敛到0的长度。支撑长度越长，一般需要耗费更多的计算时间，且产生更多高幅值的小波系数。大部分应用选择支撑长度为5~9之间的小波，因为支撑长度太长会产生边界问题，支撑长度太短消失矩太低，不利于信号能量的集中。

这里常常见到“紧支撑”的概念，通俗来讲，对于函数f(x)，如果自变量x在0附近的取值范围内，f(x)能取到值；而在此之外，f(x)取值为0，那么这个函数f(x)就是紧支撑函数，而这个0附近的取值范围就叫做紧支撑集。总结为一句话就是“除在一个很小的区域外，函数为零，即函数有速降性”。

2、对称性

具有对称性的小波，在图像处理中可以很有效地避免相位畸变，因为该小波对应的滤波器具有线性相位的特点。

3、消失矩

在实际中，对基本小波往往不仅要求满足容许条件，对还要施加所谓的消失矩（Vanishing Moments）条件，使尽量多的小波系数为零或者产生尽量少的非零小波系数，这样有利于数据压缩和消除噪声。消失矩越大，就使更多的小波系数为零。但在一般情况下，消失矩越高，支撑长度也越长。所以在支撑长度和消失矩上，我们必须要折衷处理。

        小波的消失矩的定义为，若 

 其中，Ψ(t)为基本小波，0<=p<N。则称小波函数具有N阶消失矩。从上式还可以得出，同任意n-1阶多项式正交。在频域内表示就是Ψ(ω)在ω=0处有高阶零点（一阶零点就是容许条件）。

4、正则性

在量化或者舍入小波系数时，为了减小重构误差对人眼的影响，我们必须尽量增大小波的光滑性或者连续可微性。因为人眼对“不规则”(irregular)误差比“平滑”误差更加敏感。换句话说，我们需要强加“正则性”(regularity)条件。也就是说正则性好的小波，能在信号或图像的重构中获得较好的平滑效果，减小量化或舍入误差的视觉影响。但在一般情况下，正则性好，支撑长度就长，计算时间也就越大。因此正则性和支撑长度上，我们也要有所权衡。

消失矩和正则性之间有很大关系，对很多重要的小波（比如，样条小波，Daubechies小波等）来说，随着消失矩的增加，小波的正则性变大，但是，并不能说随着小波消失矩的增加，小波的正则性一定增加，有的反而变小。

5、相似性

选择和信号波形相似的小波，这对于压缩和消噪是有参考价值的。
 

二、常见的小波基

以下列出的15种小波基是Matlab中支持的15种。

# python # 图像处理 # 算法
Logo
GitCode 开源社区

旨在为数千万中国开发者提供一个无缝且高效的云端环境，以支持学习、使用和贡献开源项目。

更多推荐

cover

GitTalk | DevUI Admin 前端项目构建

avatar GitCode 开源社区
cover

GitTalk | 使用面向业务的狮偶编程语言提升开发效率

avatar GitCode 开源社区
cover

GitTalk | DevUI Suits 场景解决方案

avatar GitCode 开源社区