医学图像处理案例（二十四）——基于cuda小波变换和cuda脉冲耦合神经网络的图像融合
医学图像处理案例（二十四）——基于cuda小波变换和cuda脉冲耦合神经网络的图像融合
今天将介绍使用cuda小波变换和cuda脉冲耦合神经网络来对多景深图像进行融合。
1、图像融合概述
图像融合（Image Fusion）是指将多源信道所采集到的关于同一目标的图像数据经过图像处理和计算机技术等，最大限度的提取各自信道中的有利信息，最后综合成高质量的图像，以提高图像信息的利用率、改善计算机解译精度和可靠性、提升原始图像的空间分辨率和光谱分辨率，利于监测。
2、小波变换特点介绍
小波变换的固有特性使其在图像处理中有如下优点：完善的重构能力，保证信号在分解过程中没有信息损失和冗余信息；把图像分解成低频图像和细节（高频）图像的组合，分别代表了图像的不同结构，因此容易提取原始图像的结构信息和细节信息；小波分析提供了与人类视觉系统方向相吻合的选择性图像。
一般图像融合的小波分解采用离散小波变换(Discrete Wavelet Transform, DWT)。DWT的函数基由一个称为母小波或分析小波的单一函数通过膨胀和平移获得。因而，DWT同时具有时域和频域分析能力，与一般的金字塔分解相比，DWT图像分解具有以下优势：
1）具有方向性，在提取图像低频信息的同时，还可获得了水平、垂直和对角三个方向的高频信息；
2）通过合理的选择母小波，可使DWT在压缩噪声的同时更有效的提取纹理、边缘等显著信息；
3）金字塔分解各尺度之间具有信息的相关性，而DWT在不同尺度上具有更高的独立性。
3、基于小波变换的图像融合
DWT 融合算法基本思想：首先对源图像进行小波变换，然后按照一定规则对变换系数进行合并；最后对合并后的系数进行小波逆变换得到融合图像。
3.1、小波分解原理简介
LL：水平低频，垂直低频
LH：水平低频，垂直高频
HL：水平高频，垂直低频
HH：水平高频，垂直高频
其中，L表示低频，H表示高频，下标1、2表示一级或二级分解。在每一分解层上，图像均被分解为LL，LH，HH和HL四个频带，下一层的分解仅对低频分量LL进行分解。这四个子图像中的每一个都是由原图与一个小波基函数的内积后，再经过在x和y方向都进行2倍的间隔采样而生成的，这是正变换，也就是图像的分解；逆变换，也就是图像的重建，是通过图像的增频采样和卷积来实现的。
3.2、融合规则
规则一：系数绝对值较大法
该融合规则适合高频成分比较丰富，亮度、对比度比较高的源图像，否则在融合图像中只保留一幅源图像的特征，其他的特征被覆盖。小波变换的实际作用是对信号解相关，并将信号的全部信息集中到一部分具有大幅值的小波系数中。这些大的小波系数含有的能量远比小系数含有的能量大，从而在信号的重构中，大的系数比小的系数更重要。
规则二：加权平均法
权重系数可调，适用范围广，可消除部分噪声，源图像信息损失较少，但会造成图像对比度的下降，需要增强图像灰度。
规则三、脉冲耦合神经网络（PCNN）
PCNN模型是由很多神经元相互连接而形成的单层循环网络，其中单个神经元是由接收区域、耦合调制域和脉冲发生域组成，单个神经元模型如下图所示。
接收域由循环输入Fij线路和线性连通输入Lij线路构成，其中循环输入Fij直接获取来自外部的刺激信号输入Sij，线性连通输入Lij获取来自一定区域内相连接的神经元信号Ykl。耦合调制域是将带有偏置的线性连通输入单元与循环输入单元进行相乘来得到神经元的内部参数Uij。脉冲发生域的构成有脉冲发生器和临界值变化的匹配器，在神经元的内部参数Uij大于该神经元膜电位的动态临界值θij时，神经元会输出一个Yij信息。
当PCNN模型用于处理二维图像时，可以用数学离散形式来描述，如下公式所示。
其中Sij为外部信号，αL和αθ分别是线性连通输入Lij和动态临界值θij的衰减定值，VL和Vθ分别是连通倍数系数和临界值倍数系数，Wijkl是线性连通输入Lij的加权系数，βij为连接强度，决定了线性连通输入Lij对内部参数Uij的贡献。
4、基于cuda小波变换和cuda脉冲耦合神经网络的图像融合代码实现
将分享python版本代码来实现多景深医学图像融合，融合策略是低频图像采用平均值法，高频图像采用PCNN最大值法，PCNN参数设置：链接系数为5，链接参数为0.1，迭代次数为200。python版本中需要用ptwt库，可以使用下面命令来安装，具体可以见原文链接。
pip install ptwtpython版本（cuda加速）代码
import pywt
import ptwt
import torch
import numpy as np
import cv2
import time
import math
from scipy import signal
use_cuda = True
device = torch.device('cuda' if use_cuda else 'cpu')
def pcnn_cuda(img, link=5, beta=0.1, iteration=200):
"""
PCNN generate fire maps image
:param img:source image
:param link:pcnn link parm
:param beta:pcnn link coefficient
:param iteration:pcnn iteration number
:return:fire maps image
"""
b, c = img.shape[0], img.shape[1]
alpha_L = torch.tensor(1).to(device)
alpha_Theta = torch.tensor(0.2).to(device)
vL = torch.tensor(1.0).to(device)
vTheta = torch.tensor(20).to(device)
size = (link, link)
W = torch.zeros(size).to(device)
center = torch.tensor(size) // 2
for i in range(size[0]):
for j in range(size[1]):
W[i][j] = torch.sqrt(torch.sum((torch.tensor([i, j]) - center) ** 2))
W = 1.0 / W
W[center[0], center[1]] = 0
W = W.unsqueeze(0)
W = W.unsqueeze(0)
W = W.repeat(b, c, 1, 1)
imgf = img.to(torch.float)
F = abs(imgf)
L = torch.zeros_like(F).to(device)
Y = torch.zeros_like(F).to(device)
Theta = torch.zeros_like(F).to(device)
img_pcnn = torch.zeros_like(F).to(device)
for i in range(iteration):
K = torch.nn.functional.conv2d(Y, W, padding='same')
L = torch.exp(-alpha_L) * L + vL * K
Theta = torch.exp(-alpha_Theta) * Theta + vTheta * Y
U = torch.multiply(F, 1 + torch.multiply(torch.tensor(beta), L))
Y = U > Theta
Y = Y.to(torch.float)
img_pcnn = img_pcnn + Y
return img_pcnn
def fuseCoefpcnncuda(cooef, link=5, beta=0.1, iteration=200):
pcnn = pcnn_cuda(cooef, link, beta, iteration)
pcnnmaxvalue = torch.max(pcnn, dim=0, keepdim=True)[0]
number = cooef.shape[0]
for i in range(number):
cooef[i][pcnn[i] != pcnnmaxvalue[0]] = 0
sumcooef = torch.sum(cooef, dim=0, keepdim=True)
return sumcooef
def channelTransformgpu(ch_list, FUSION_METHOD):
# wavelet transformation
data = np.array(ch_list)
data = np.transpose(data, (0, 3, 1, 2))
cooef = ptwt.wavedec2(torch.as_tensor(data).contiguous().to(device=device, dtype=torch.float32), "haar", level=1,
mode='symmetric')
cA, (cH, cV, cD) = cooef
# wavelet transformation coor fusion
cA = torch.mean(cA, dim=0, keepdim=True)
if FUSION_METHOD == "max":
cH = torch.max(cH, dim=0, keepdim=True)[0]
cV = torch.max(cV, dim=0, keepdim=True)[0]
cD = torch.max(cD, dim=0, keepdim=True)[0]
elif FUSION_METHOD == "min":
cH = torch.min(cH, dim=0, keepdim=True)[0]
cV = torch.min(cV, dim=0, keepdim=True)[0]
cD = torch.min(cD, dim=0, keepdim=True)[0]
elif FUSION_METHOD == "pcnn":
cH = fuseCoefpcnncuda(cH, 5, 0.1, 200)
cV = fuseCoefpcnncuda(cV, 5, 0.1, 200)
cD = fuseCoefpcnncuda(cD, 5, 0.1, 200)
else:
cH = torch.mean(cH, dim=0, keepdim=True)
cV = torch.mean(cV, dim=0, keepdim=True)
cD = torch.mean(cD, dim=0, keepdim=True)
fincoC = [cA, (cH, cV, cD)]
# wavelet iverse transformation
reconstruction = ptwt.waverec2(fincoC, "haar")
reconstruction = reconstruction.detach().cpu().squeeze().numpy()
reconstruction = np.transpose(reconstruction, (1, 2, 0))
return reconstruction python版本（cpu）代码
import pywt
import ptwt
import torch
import numpy as np
import cv2
import time
import math
from scipy import signal
def pcnn(img, link=5, beta=0.1, iteration=200):
"""
PCNN generate fire maps image
:param img:source image
:param link:pcnn link parm
:param beta:pcnn link coefficient
:param iteration:pcnn iteration number
:return:fire maps image
"""
m, n = np.shape(img)[0], np.shape(img)[1]
alpha_L = 1
alpha_Theta = 0.2
vL = 1.0
vTheta = 20
center_x = round(link / 2)
center_y = round(link / 2)
W = np.zeros((link, link))
for i in range(link):
for j in range(link):
if i == center_x and j == center_y:
W[i, j] = 0
else:
W[i, j] = 1. / math.sqrt(pow(i - center_x, 2) + pow(j - center_y, 2))
imgf = img.astype(np.float)
F = abs(imgf)
L = np.zeros((m, n))
Y = np.zeros((m, n))
Theta = np.zeros((m, n))
img_pcnn = np.zeros((m, n))
for i in range(iteration):
K = signal.convolve2d(Y, W, mode='same')
L = math.exp(-alpha_L) * L + vL * K
Theta = math.exp(-alpha_Theta) * Theta + vTheta * Y
U = np.multiply(F, 1 + np.multiply(beta, L))
Y = U > Theta
Yf = Y.astype(np.float)
img_pcnn = img_pcnn + Yf
return img_pcnn
def fuseCoefpcnn(cooef1, cooef2, cooef3):
pcnn1 = pcnn(cooef1, 5, 0.1, 200)
pcnn2 = pcnn(cooef2, 5, 0.1, 200)
pcnn3 = pcnn(cooef3, 5, 0.1, 200)
pcnn_cooef = np.maximum(pcnn1, pcnn2, pcnn3)
cooef1[pcnn1 != pcnn_cooef] = 0
cooef2[pcnn2 != pcnn_cooef] = 0
cooef3[pcnn3 != pcnn_cooef] = 0
cooef = cooef1 + cooef2 + cooef3
return cooef
def channelTransform(ch1, ch2, ch3, FUSION_METHOD):
# wavelet transformation
cooef1 = pywt.wavedec2(ch1, 'haar', level=1)
cooef2 = pywt.wavedec2(ch2, 'haar', level=1)
cooef3 = pywt.wavedec2(ch3, 'haar', level=1)
cA1, (cH1, cV1, cD1) = cooef1
cA2, (cH2, cV2, cD2) = cooef2
cA3, (cH3, cV3, cD3) = cooef3
# wavelet transformation coor fusion
cA = (cA1 + cA2 + cA3) / 3
if FUSION_METHOD == "max":
cH = np.maximum(cH1, cH2, cH3)
cV = np.maximum(cV1, cV2, cV3)
cD = np.maximum(cD1, cD2, cD3)
elif FUSION_METHOD == "min":
cH = np.minimum(cH1, cH2, cH3)
cV = np.minimum(cV1, cV2, cV3)
cD = np.minimum(cD1, cD2, cD3)
elif FUSION_METHOD == "pcnn":
cH = fuseCoefpcnn(cH1, cH2, cH3)
cV = fuseCoefpcnn(cV1, cV2, cV3)
cD = fuseCoefpcnn(cD1, cD2, cD3)
else:
cH = (cH1 + cH2 + cH3) / 3
cV = (cV1 + cV2 + cV3) / 3
cD = (cD1 + cD2 + cD3) / 3
fincoC = cA, (cH, cV, cD)
# wavelet iverse transformation
outImageC = pywt.waverec2(fincoC, 'haar')
return outImageC5、融合结果
下图是多个景深的宫颈癌细胞病理图像。
下图第一个是cuda小波PCNN融合结果，第二个是cpu小波PCNN融合结果，通过计算图像清晰度数值，可以看到融合后图像比融合前的三张图像质量都好，而计算时间cuda是8.57s，cpu是2135.46s。
