背景与现状：
在现实生活中，人们使用摄像机拍摄时，希望可以获得同一场景中所有景物都清晰的图像。但是摄像机镜头受景深的限制，无法同时聚焦所有目标，因此拍摄的照片中部分区域清晰，部分区域模糊。多聚焦图像融合技术可以将多幅同一场景下聚焦区域不同的图像融合成一幅全清晰的图像，从而有效地解决这个问题，提高图像的信息利用率。准确地识别和提取源图像中的聚焦区域是多聚焦图像融合算法中的难题，如果聚焦区域提取不完整，会导致融合结果中出现伪影及边缘轮廓丢失的现象。图像融合技术包含三个不同的级别，即像素级、特征级和决策级。像素级的图像融合方法直接将源图像中的原始信息进行组合，其保留边缘、纹理等细节信息的能力要优于其他两个级别。现有的多聚焦图像融合技术都是在像素级上展开的。
当前，多聚焦图像融合要解决的核心问题是：如何将源图像中的聚焦区域识别和提取出来组成一幅全清晰的图像。解决这个问题的方法是准确找到聚焦和非聚焦区域的边界。许多学者通过构造不同的像素清晰度测量方法，从而正确区分源图像中哪些是聚焦像素，哪些是模糊像素，提高融合图像的性能。近年来，越来越多的研究人员致力于多聚焦图像融合方法的研究，并提出了许多有效的融合算法。现有的像素级多聚焦图像融合方法主要分为三类：基于空间域的方法、基于变换域的方法和基于深度学习的方法。
基于空间域的图像融合方法一般是用固定大小的窗口或者图像分割技术将源图像分成一定大小的图像块，然后通过某种清晰度测量方法分辨出源图像中的清晰区域和模糊区域，最后将这些清晰的图像块组合成融合图像。
基于变换域的图像融合方法通常是利用某种变换将源图像分解成不同频带的子图像，然后选取不同的融合规则得到融合后的子图像，最后通过逆变换生成融合图像。常用的变换域融合方法有多尺度几何变换法、稀疏表示法、替代法等。
基于深度学习的图像融合方法包括利用脉冲耦合神经网络（Pulse Coupled Neural Network，PCNN）、卷积神经网络（Convolutional NeuralNetwork，CNN）等方法提取图像特征并构建融合权重图，提高了聚焦区域和非聚焦区域分类的准确性，依靠后处理可以有效提高图像融合的质量。
算法思路：
目前比较常用的图像融合算法有线性融合、羽化融合、多频段融合，其中多频段融合通过在变换域分别对图像的不同频率成分做融合，既能保留高频轮廓，又能消除伪影。本文算法使用多频段融合做为基础。
此外，对于多张不同景深的图像，融合时需要有效地选择合适的掩码。关于掩码选择的方法，我大致描述如下：
基于空间域将每张图像分成固定大小的网格，然后计算每个网格的平均梯度，形成这个一个平均梯度矩阵数组。然后遍历每一个网格，将所有图像相同位置的平均梯度对对比，选取值最大的那张设置该网格掩码为1，其余的设为0。当所有网格遍历完成，所有图像的掩码矩阵也基本形成。然后，对所有图像及其掩码构建高斯金字塔，对图像构建拉普拉斯金字塔，并逐层融合所有图像的拉普拉斯图像，形成最后融合效果。
高斯金字塔：
常用的建立高斯金字塔的方法是自底向上一层一层地迭代计算生成。该方法通常分为两个步骤：第一，使用高斯低通滤波器对图像进行平滑（高斯模糊），第二，对经高斯平滑处理后的图像进行下采样，得到一系列尺寸的图像。减少一幅图像的尺寸，仅仅靠采样会丢失许多信息。根据采样定理，需要让所有小于最短波长的1/4采样而得到的精细结构能通过平滑滤波器来消除掉，这样才能获得一幅正确的采样图像，故采样前先进行高斯平滑。如图1所示为高斯金字塔模型的建立过程。从图中可以看出高斯金字塔的第K层图像通过高斯模糊和降采样就可得到第K+1层的高斯图像。这边的降采样为隔行隔列降采样，即可将所有的偶数行列去除，得到上层图像的像素为下层图像的1/4。通过这两个步骤反复迭代就能得到一个完整的高斯金字塔。
通过不断迭代计算就可以生成一系列的图像G0，G1…，GN,构成一个完整的高斯金字塔，G0即为原图像即为高斯金字塔的底层，GN为高斯金字塔的顶层。除了上述这种构建方法外，还可以直接将原始图像即金字塔底层与一组等效加权函数进行卷积生成高斯金字塔的各个层。如下图所示即为图像通过高斯金字塔处理后生成的一系列不同尺度的图像。
拉普拉斯金字塔：
拉普拉斯金字塔可以认为是残差金字塔，用来存储下采样后图片与原始图片的差异。我们知道，如果高斯金字塔中任意一张图Gi（比如G0为最初的高分辨率图像）先进行下采样得到图Down(Gi)，再进行上采样得到图Up(Down(Gi))，得到的Up(Down(Gi))与Gi是存在差异的，因为下采样过程丢失的信息不能通过上采样来完全恢复，也就是说下采样是不可逆的。原始图片下采样后得到的小尺寸图片虽然保留了视觉效果，但是将该小尺寸图像再次上采样也不能完整的恢复出原始图像。为了能够从下采样图像Down(Gi)中还原原始图像Gi，我们需要记录再次上采样得到Up(Down(Gi))与原始图片Gi之间的差异，这就是拉普拉斯金字塔的核心思想。下图是一张拉普拉斯金字塔过程的描述：
算法实现：
opencv自带了多频段融合算法（MultiBandBlender）的实现，经初步使用，它对传入的图像格式做了较多的要求，如要求必须为三通道等等，而且性能也是个问题。下面的实现同时支持CV_8UC3和CV_8UC1两种格式，因而也可以扩展到使用BAYER格式的场景。
主要代码如下：
// 图像金字塔信息对象
class CImagePyramidInfo
{
private:
cv::Mat m_matMeanFocus; // 平均梯度网格矩阵
std::vector<cv::Mat> m_vecGaussian; // 高斯金字塔
std::vector<cv::Mat> m_vecLaplac; // 拉普拉斯金字塔
std::vector<cv::Mat> m_vecMask; // 掩码金字塔
public:
CImagePyramidInfo(const cv::Mat& matImage, int nGradSize = 100)
{
__makeMeanFocus(matImage, m_matMeanFocus, nGradSize);
cv::Mat matImageFloat;
matImage.convertTo(matImageFloat, CV_32F, 1.0 / 255);
__makeGaussian(matImageFloat, m_vecGaussian, 5);
__makeLaplac(m_vecGaussian, m_vecLaplac);
}
virtual ~CImagePyramidInfo()
{
}
// 构建掩码金字塔
void makeMask(const cv::Mat& matMask)
{
__makeMask(matMask, m_vecMask, m_vecGaussian.size());
}
// 获取平均梯度网格矩阵
const cv::Mat& getMeanFocus() const
{
return m_matMeanFocus;
}
// 获取高斯金字塔
const std::vector<cv::Mat>& getGaussian() const
{
return m_vecGaussian;
}
// 获取拉普拉斯金字塔
const std::vector<cv::Mat>& getLaplac() const
{
return m_vecLaplac;
}
// 获取掩码金字塔
const std::vector<cv::Mat>& getMask() const
{
return m_vecMask;
}
private:
// 构建平均梯度网格矩阵
static void __makeMeanFocus(const cv::Mat& matImage, cv::Mat& matMeanFocus, int nGradSize = 100)
{
// 只接受CV_8UC3和CV_8UC1格式的图像
if (matImage.type() != CV_8UC3 && matImage.type() != CV_8UC1)
return;
// 计算网格数
int nGradRows = matImage.rows / nGradSize;
if (matImage.rows % nGradSize)
++nGradRows;
int nGradCols = matImage.cols / nGradSize;
if (matImage.cols % nGradSize)
++nGradCols;
matMeanFocus = cv::Mat(nGradRows, nGradCols, CV_32FC1, cv::Scalar::all(0.0));
// 计算每个网格的平均梯度
for (int nRow = 0; nRow < nGradRows; ++nRow)
{
for (int nCol = 0; nCol < nGradCols; ++nCol)
{
int nOffsetX = nCol * nGradSize;
int nOffsetY = nRow * nGradSize;
int nWidth = nOffsetX + nGradSize > matImage.cols ? matImage.cols - nOffsetX : nGradSize;
int nHeight = nOffsetY + nGradSize > matImage.rows ? matImage.rows - nOffsetY : nGradSize;
cv::Mat matImageROISrc = matImage(cv::Rect(nOffsetX, nOffsetY, nWidth, nHeight));
cv::Mat matImageROIDst;
// 将BGR转换为灰度图像
if (matImage.type() == CV_8UC3)
{
matImageROIDst = cv::Mat(matImageROISrc.size(), CV_8UC1);
cv::cvtColor(matImageROISrc, matImageROIDst, cv::COLOR_BGR2GRAY);
}
else
{
matImageROIDst = matImageROISrc.clone();
}
cv::GaussianBlur(matImageROIDst, matImageROIDst, cv::Size(3, 3), 0, 0);
cv::Mat matImageROISobel;
cv::Sobel(matImageROIDst, matImageROISobel, CV_8UC1, 1, 1, 3);
matMeanFocus.at<float>(nRow, nCol) = cv::mean(matImageROISobel)[0];
}
}
}
// 构建高斯金字塔
static void __makeGaussian(const cv::Mat& matImage, std::vector<cv::Mat>& vecGaussian, int nLayerCount = 5)
{
vecGaussian.resize(nLayerCount);
for (int nLayerIndex = 0; nLayerIndex < nLayerCount; ++nLayerIndex)
{
if (nLayerIndex == 0)
{
vecGaussian[nLayerIndex] = matImage.clone();
}
else
{
cv::Mat matImageDown;
cv::pyrDown(vecGaussian[nLayerIndex - 1], matImageDown);
vecGaussian[nLayerIndex] = matImageDown;
}
}
}
// 构建拉普拉斯金字塔
static void __makeLaplac(const std::vector<cv::Mat>& vecGaussian, std::vector<cv::Mat>& vecLaplac)
{
vecLaplac.resize(vecGaussian.size());
for (int nLayerIndex = vecLaplac.size() - 1; nLayerIndex >= 0; --nLayerIndex)
{
if (nLayerIndex == vecLaplac.size() - 1)
{
vecLaplac[nLayerIndex] = vecGaussian[nLayerIndex].clone();
}
else
{
cv::Mat matImageUp;
cv::pyrUp(vecGaussian[nLayerIndex + 1], matImageUp, vecGaussian[nLayerIndex].size());
vecLaplac[nLayerIndex] = vecGaussian[nLayerIndex] - matImageUp;
}
}
}
// 构建掩码金字塔
static void __makeMask(const cv::Mat& matMask, std::vector<cv::Mat>& vecMask, int nLayerCount = 5)
{
vecMask.resize(nLayerCount);
for (int nLayerIndex = 0; nLayerIndex < vecMask.size(); ++nLayerIndex)
{
if (nLayerIndex == 0)
{
vecMask[nLayerIndex] = matMask.clone();
}
else
{
cv::Mat matMaskDown;
cv::pyrDown(vecMask[nLayerIndex - 1], matMaskDown);
vecMask[nLayerIndex] = matMaskDown;
}
}
}
};
// 多聚焦多频段图像融合算法
static bool multiFocusBandBlender(const std::vector<cv::Mat>& vecImage, cv::Mat& matResult, int nGradSize = 100)
{
if (vecImage.size() < 2)
return false;
// 检查每张图片的大小和格式是否一致
for (int i = 1; i < vecImage.size(); ++i)
{
if (vecImage[i].size() != vecImage[i - 1].size())
return false;
if (vecImage[i].type() != vecImage[i - 1].type())
return false;
}
// 检查每张图片的格式是否为CV_8UC3和CV_8UC1格式
for (int i = 0; i < vecImage.size(); ++i)
{
if (vecImage[i].type() != CV_8UC3 && vecImage[i].type() != CV_8UC1)
return false;
}
// 为每张图像建立高斯和拉普拉斯金字塔
std::vector<CImagePyramidInfo*> vecImagePyramidInfo(vecImage.size());
for (int i = 0; i < vecImage.size(); ++i)
{
vecImagePyramidInfo[i] = new CImagePyramidInfo(vecImage[i], nGradSize);
}
// 合成聚焦网格矩阵
std::vector<cv::Mat> vecMeanFocus(vecImagePyramidInfo.size());
std::vector<cv::Mat> vecBestIndex(vecImagePyramidInfo.size());
for (int i = 0; i < vecImagePyramidInfo.size(); ++i)
{
vecMeanFocus[i] = vecImagePyramidInfo[i]->getMeanFocus();
vecBestIndex[i] = cv::Mat(vecMeanFocus[i].size(), CV_8UC1, cv::Scalar::all(0));
}
for (int nRow = 0; nRow < vecMeanFocus[0].rows; ++nRow)
{
for (int nCol = 0; nCol < vecMeanFocus[0].cols; ++nCol)
{
// 找梯度最大图片
float fMax = -1.0;
int nIndex = -1;
for (int i = 0; i < vecMeanFocus.size(); ++i)
{
float fMeanFocus = vecMeanFocus[i].at<float>(nRow, nCol);
if (fMeanFocus > fMax)
{
fMax = fMeanFocus;
nIndex = i;
}
}
// 标记该图片当前网格为选中
vecBestIndex[nIndex].at<uchar>(nRow, nCol) = 1;
}
}
// 构建掩码金字塔
std::vector<cv::Mat> vecMaskPyramid(vecBestIndex.size());
for (int i = 0; i < vecBestIndex.size(); ++i)
{
vecMaskPyramid[i] = cv::Mat(vecImage[i].size(), (vecImage[i].channels() == 1 ? CV_32FC1 : CV_32FC3), cv::Scalar::all(0.0));
for (int nRow = 0; nRow < vecBestIndex[i].rows; ++nRow)
{
for (int nCol = 0; nCol < vecBestIndex[i].cols; ++nCol)
{
if (vecBestIndex[i].at<uchar>(nRow, nCol) == 1)
{
int nOffsetX = nCol * nGradSize;
int nOffsetY = nRow * nGradSize;
int nWidth = nOffsetX + nGradSize > vecImage[0].cols ? vecImage[0].cols - nOffsetX : nGradSize;
int nHeight = nOffsetY + nGradSize > vecImage[0].rows ? vecImage[0].rows - nOffsetY : nGradSize;
cv::Mat matROI = vecMaskPyramid[i](cv::Rect(nOffsetX, nOffsetY, nWidth, nHeight));
matROI.setTo(cv::Scalar::all(1.0));
}
}
}
}
for (int i = 0; i < vecImagePyramidInfo.size(); ++i)
{
vecImagePyramidInfo[i]->makeMask(vecMaskPyramid[i]);
}
// 逐层融合拉普拉斯金字塔
// 层次遍历
for (int nLayerIndex = vecImagePyramidInfo[0]->getLaplac().size() - 1; nLayerIndex >= 0; --nLayerIndex)
{
cv::Mat matBlendLayer;
// 图片遍历
for (int i = 0; i < vecImagePyramidInfo.size(); ++i)
{
cv::Mat matLaplac = vecImagePyramidInfo[i]->getLaplac()[nLayerIndex];
cv::Mat matMask = vecImagePyramidInfo[i]->getMask()[nLayerIndex];
if (matBlendLayer.empty())
{
matBlendLayer = matLaplac.mul(matMask);
}
else
{
matBlendLayer += matLaplac.mul(matMask);
}
}
if (matResult.empty())
{
matResult = matBlendLayer;
}
else
{
cv::pyrUp(matResult, matResult, matBlendLayer.size());
matResult += matBlendLayer;
}
}
cv::convertScaleAbs(matResult, matResult, 255);
for (auto iter = vecImagePyramidInfo.begin(); iter != vecImagePyramidInfo.end(); ++iter)
{
delete (*iter);
(*iter) = NULL;
}
return true;
}
调用的测试代码：
// 多聚焦多频段图像融合算法
{
cv::Mat matImg0 = cv::imread("full2.jpg");
cv::Mat matImg1 = matImg0.clone();
cv::Mat matImg2 = matImg0.clone();
cv::Mat matImg3 = matImg0.clone();
//cv::cvtColor(matImg1, matImg1, cv::COLOR_BGR2GRAY);
//cv::cvtColor(matImg2, matImg2, cv::COLOR_BGR2GRAY);
//cv::cvtColor(matImg3, matImg3, cv::COLOR_BGR2GRAY);
const int nGradSize = 50;
// 随机模糊一些区域
int nGradRows = matImg0.rows / nGradSize;
int nGradCols = matImg0.cols / nGradSize;
cv::RNG rng1(1);
for (int i = 0; i < 100; ++i)
{
int nRow = rng1.uniform(0, nGradRows);
int nCol = rng1.uniform(0, nGradCols);
cv::Mat matROI = matImg1(cv::Rect(nCol * nGradSize, nRow * nGradSize, nGradSize, nGradSize));
cv::GaussianBlur(matROI, matROI, cv::Size(5, 5), 2.0);
}
cv::RNG rng2(2);
for (int i = 0; i < 100; ++i)
{
int nRow = rng2.uniform(0, nGradRows);
int nCol = rng2.uniform(0, nGradCols);
cv::Mat matROI = matImg2(cv::Rect(nCol * nGradSize, nRow * nGradSize, nGradSize, nGradSize));
cv::GaussianBlur(matROI, matROI, cv::Size(5, 5), 2.0);
}
cv::RNG rng3(3);
for (int i = 0; i < 100; ++i)
{
int nRow = rng3.uniform(0, nGradRows);
int nCol = rng3.uniform(0, nGradCols);
cv::Mat matROI = matImg3(cv::Rect(nCol * nGradSize, nRow * nGradSize, nGradSize, nGradSize));
cv::GaussianBlur(matROI, matROI, cv::Size(5, 5), 2.0);
}
int64 llTickStart = cv::getTickCount();
// 这里用于性能测试
for (int n = 0; n < 1000; ++n)
{
std::vector<cv::Mat> vecImage;
vecImage.push_back(matImg1);
vecImage.push_back(matImg2);
vecImage.push_back(matImg3);
cv::Mat matResult;
multiFocusBandBlender(vecImage, matResult, nGradSize);
cv::imwrite("full2_result.jpg", matResult);
}
int64 llTickEnd = cv::getTickCount();
printf("use %fs \n", (llTickEnd - llTickStart) / cv::getTickFrequency());
return 0;
}
运行效果：
原图：
生成的三张随机区域模糊图：
图1：
图2：
图3：
合成的效果图：
