医学图像处理案例（二十三）——基于cuda的小波变换的3d图像融合
医学图像处理案例（二十三）——基于cuda的小波变换的3d图像融合
今天将介绍使用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、融合规则
规则一：系数绝对值较大法
该融合规则适合高频成分比较丰富，亮度、对比度比较高的源图像，否则在融合图像中只保留一幅源图像的特征，其他的特征被覆盖。小波变换的实际作用是对信号解相关，并将信号的全部信息集中到一部分具有大幅值的小波系数中。这些大的小波系数含有的能量远比小系数含有的能量大，从而在信号的重构中，大的系数比小的系数更重要。
规则二：加权平均法
权重系数可调，适用范围广，可消除部分噪声，源图像信息损失较少，但会造成图像对比度的下降，需要增强图像灰度。
4、基于cuda小波变换的多模态医学图像融合代码实现
将分享python版本代码来融合多模态MR图像，融合策略是低频图像采用平均值法，高频图像采用最大值法。python版本中需要用到ptwt库，可以使用下面命令来安装，具体可以见原文链接。
pip install ptwtpython版本（cuda加速）代码
import pywt
import ptwt
import torch
import numpy as np
import SimpleITK as sitk
import time
import os
os.environ['KMP_DUPLICATE_LIB_OK'] = 'True'
# Use CUDA
os.environ['CUDA_VISIBLE_DEVICES'] = '-1'
os.environ['CUDA_LAUNCH_BLOCKING'] = '0'
use_cuda = torch.cuda.is_available()
device = torch.device('cuda' if use_cuda else 'cpu')
# This function does the coefficient fusing according to the fusion method
def fuseCoeffcuda(cooef1, cooef2, method):
if method == 'mean':
cooef = (cooef1 + cooef2) / 2
elif method == 'min':
cooef = torch.min(cooef1, cooef2)
elif method == 'max':
cooef = torch.max(cooef1, cooef2)
return cooef
def channelTransform3dgpu(ch1, ch2, FUSION_METHOD):
# wavelet transformation
cooef1 = ptwt.wavedec3(torch.as_tensor(ch1).contiguous().to(device=device, dtype=torch.float32),
pywt.Wavelet("haar"), level=1, mode="symmetric")
cooef2 = ptwt.wavedec3(torch.as_tensor(ch2).contiguous().to(device=device, dtype=torch.float32),
pywt.Wavelet("haar"), level=1, mode="symmetric")
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(fuseCoeffcuda(cooef1[0], cooef2[0], 'mean'))
else:
c1 = fuseCoeffcuda(cooef1[i]['aad'], cooef2[i]['aad'], FUSION_METHOD)
c2 = fuseCoeffcuda(cooef1[i]['ada'], cooef2[i]['ada'], FUSION_METHOD)
c3 = fuseCoeffcuda(cooef1[i]['add'], cooef2[i]['add'], FUSION_METHOD)
c4 = fuseCoeffcuda(cooef1[i]['daa'], cooef2[i]['daa'], FUSION_METHOD)
c5 = fuseCoeffcuda(cooef1[i]['dad'], cooef2[i]['dad'], FUSION_METHOD)
c6 = fuseCoeffcuda(cooef1[i]['dda'], cooef2[i]['dda'], FUSION_METHOD)
c7 = fuseCoeffcuda(cooef1[i]['ddd'], cooef2[i]['ddd'], FUSION_METHOD)
dictobj = {'aad': c1, 'ada': c2, 'add': c3, 'daa': c4, 'dad': c5, 'dda': c6, 'ddd': c7}
fusedCooef.append(dictobj)
# wavelet iverse transformation
reconstruction = ptwt.waverec3(fusedCooef, pywt.Wavelet("haar"))
reconstruction = reconstruction.detach().cpu().squeeze().numpy()
return reconstructionpython版本（cpu）代码
import pywt
import ptwt
import torch
import numpy as np
import SimpleITK as sitk
import time
import os
# 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
def channelTransform3d(ch1, ch2, FUSION_METHOD):
# wavelet transformation
cooef1 = pywt.wavedecn(ch1, 'haar', level=1)
cooef2 = pywt.wavedecn(ch2, 'haar', level=1)
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], 'mean'))
else:
c1 = fuseCoeff(cooef1[i]['aad'], cooef2[i]['aad'], FUSION_METHOD)
c2 = fuseCoeff(cooef1[i]['ada'], cooef2[i]['ada'], FUSION_METHOD)
c3 = fuseCoeff(cooef1[i]['add'], cooef2[i]['add'], FUSION_METHOD)
c4 = fuseCoeff(cooef1[i]['daa'], cooef2[i]['daa'], FUSION_METHOD)
c5 = fuseCoeff(cooef1[i]['dad'], cooef2[i]['dad'], FUSION_METHOD)
c6 = fuseCoeff(cooef1[i]['dda'], cooef2[i]['dda'], FUSION_METHOD)
c7 = fuseCoeff(cooef1[i]['ddd'], cooef2[i]['ddd'], FUSION_METHOD)
dictobj = {'aad': c1, 'ada': c2, 'add': c3, 'daa': c4, 'dad': c5, 'dda': c6, 'ddd': c7}
fusedCooef.append(dictobj)
# wavelet iverse transformation
outImageC = pywt.waverecn(fusedCooef, 'haar')
return outImageC5、融合结果
下图是T1图像和T2图像。
cuda和cpu融合结果。
cuda程序运行时间：0.1 秒
cpu程序运行时间：0.29 秒如果碰到任何问题，随时留言，我会尽量去回答的。