图像融合：拉普拉斯金字塔融合算法

继图像拼接的课程设计之后，对这方面依旧十分感兴趣。很巧合的是，数图老师表示刚好手上有这么一个项目，要用到这方面的知识，可以让我去作为毕业设计。虽然距离毕业还远，不过如果能选到一个感兴趣并且有一定深度的题目还是很好的。这些天在看论文的时候，了解到这么一个强大的算法，打算写篇文记录下心得吧。

一些图像融合算法

图像拼接主要可以分为两个步骤：图像配准和图像融合。其中图像配准的目的是将图一场景中不同视角的图像投影到同一平面并进行对准。比如我之前这篇博客中使用SIFT特征检测和单应矩阵的目的，就是进行图像配准。

图像配准
经过图像配准之后，就需要进行图像融合。而图像融合的目的就是使两幅图像的重叠区域过渡自然且平滑。在上图中，可以看到明显的边界，这对拼接来说是无法接受的。这主要是因为外部亮度的变化（天空飘过了一朵萌萌的云彩？）以及曝光时相机参数不一致导致的。要消除或缓和这种现象，就需要进行图像融合。

主流的图像融合算法有：

1）加权平均法。这个很好理解，即简单的使用加权的方式从左边过渡到右边。这种方法效果一般，但算法实现极其简单，速度快。课设时我用的就是这个方法。

2）羽化算法 。这种方法过渡会比加权平均法自然，但会造成不好的模糊效果。

3）拉普拉斯金字塔融合。有的地方也称为多分辨率融合算法。这种方法是将图像建立一个拉普拉斯金字塔，其中金字塔的每一层都包含了图像不同的频段。分开不同频段进行融合效果出奇的好。这也是本文主要介绍的方法。

高斯金字塔、拉普拉斯金字塔

之前在写SIFT相关博客的时候说到过高斯金字塔。图像金字塔的意思无非就是对原图进行下采样，然后塞到一个C++的Vector或者其他什么语言中的数组里。在可视化的时候，最大的图像放在最下面，最小的图像放在最上面，所以称为图像金字塔。

图像金字塔
而高斯金字塔的每一层的构建步骤分为两步：首先对下一层的图像进行高斯模糊。这个步骤相信读者都了解，是图像处理中最基本的概念。然后删除模糊后的图像的偶数行和列，就得到了当前层的图像了。不断进行这个步骤，最终就得到了高斯金字塔。

拉普拉斯金字塔的构造需要用到高斯金字塔。拉普拉斯金字塔第i层的数学定义如下

每一层的定义

意思是拉普拉斯金字塔每一层的图像为同一层高斯金字塔的图像减去上一层的图像进行上采样并高斯模糊的结果。说的有点绕，可以看网上的这幅图进行理解。

高斯金字塔与拉普拉斯金字塔的关系

算法原理

1）首先建立两幅图像高斯金字塔，然后建立一定层数的拉普拉斯金字塔。拉普拉斯金字塔的层数越高，融合效果越好。层数N作为一个参数。

2）传入一个mask掩膜，代表了融合的位置。比如说想在两图的中间进行融合，那么掩膜图像的左半为255，右半为0，反过来是一样的。根据这个mask建立一个高斯金字塔，用于后续融合，层数为N+1。

3）根据mask将两幅图像的拉普拉斯金字塔的图像进行相加，mask为权值。相加的结果即为一个新的金字塔。同时，两幅图像的高斯金字塔的N+1层也进行这个操作，记这个图像为IMG1。

4）根据这个新的金字塔重建出最终的图像，重建的过程跟一般的拉普拉斯金字塔一样。首先对IMG1上采样，然后跟新金字塔的顶层相加，得到IMG2。IMG2进行上采样后跟下一层相加，得到IMG3。重复这个过程，最终得到的结果就是拉普拉斯金字塔融合算法的结果。

因为mask建立金字塔的过程中使用了高斯模糊，所以融合的边缘是比较平滑的。

算法实现

类实现主要参考其他博客

/**************************************************************
 * Created by 杨帮杰 on 1/1/2019
 * Right to use this code in any way you want without
 * warranty, support or any guarantee of it working
 * E-mail: yangbangjie1998@qq.com
 * Association: SCAU 华南农业大学
 **************************************************************/
#include <opencv2/opencv.hpp>
#include <iostream>
#include <string>

using namespace std;
using namespace cv;

#define IMG1_PATH "/home/jacob/图片/blend1.jpg"
#define IMG2_PATH "/home/jacob/图片/blend2.jpg"

/**
 * @brief The LaplacianBlending class
 * @private leftImg & rightImg 用于拼接的左图和右图
 * @private blendMask 用于融合的掩膜，值为加权平均的系数
 * @ref https://blog.csdn.net/abcjennifer/article/details/7628655
 */
class LaplacianBlending {
private:
    Mat leftImg;
    Mat rightImg;
    Mat blendMask;

    //Laplacian Pyramids
    vector<Mat> leftLapPyr, rightLapPyr, resultLapPyr;
    Mat leftHighestLevel, rightHighestLevel, resultHighestLevel;
    //mask为三通道方便矩阵相乘
    vector<Mat> maskGaussianPyramid;

    int levels;

    void buildPyramids()
    {
        buildLaplacianPyramid(leftImg, leftLapPyr, leftHighestLevel);
        buildLaplacianPyramid(rightImg, rightLapPyr, rightHighestLevel);
        buildGaussianPyramid();
    }

    void buildGaussianPyramid()
    {
        //金字塔内容为每一层的掩模
        assert(leftLapPyr.size()>0);

        maskGaussianPyramid.clear();
        Mat currentImg;
        cvtColor(blendMask, currentImg, CV_GRAY2BGR);
        //保存mask金字塔的每一层图像
        maskGaussianPyramid.push_back(currentImg); //0-level

        currentImg = blendMask;
        for (int l = 1; l<levels + 1; l++) {
            Mat _down;
            if (leftLapPyr.size() > l)
                pyrDown(currentImg, _down, leftLapPyr[l].size());
            else
                pyrDown(currentImg, _down, leftHighestLevel.size()); //lowest level

            Mat down;
            cvtColor(_down, down, CV_GRAY2BGR);
            //add color blend mask into mask Pyramid
            maskGaussianPyramid.push_back(down);
            string winName = to_string(l);
            imshow(winName,down);
            currentImg = _down;
        }
    }

    void buildLaplacianPyramid(const Mat& img, vector<Mat>& lapPyr, Mat& HighestLevel)
    {
        lapPyr.clear();
        Mat currentImg = img;
        for (int l = 0; l<levels; l++) {
            Mat down, up;
            pyrDown(currentImg, down);
            pyrUp(down, up, currentImg.size());
            Mat lap = currentImg - up;
            lapPyr.push_back(lap);
            currentImg = down;
        }
        currentImg.copyTo(HighestLevel);
    }

    Mat reconstructImgFromLapPyramid()
    {
        //将左右laplacian图像拼成的resultLapPyr金字塔中每一层
        //从上到下插值放大并与残差相加，即得blend图像结果
        Mat currentImg = resultHighestLevel;
        for (int l = levels - 1; l >= 0; l--)
        {
            Mat up;
            pyrUp(currentImg, up, resultLapPyr[l].size());
            currentImg = up + resultLapPyr[l];
        }
        return currentImg;
    }

    void blendLapPyrs()
    {
        //获得每层金字塔中直接用左右两图Laplacian变换拼成的图像resultLapPyr
        resultHighestLevel = leftHighestLevel.mul(maskGaussianPyramid.back()) +
            rightHighestLevel.mul(Scalar(1.0, 1.0, 1.0) - maskGaussianPyramid.back());
        for (int l = 0; l<levels; l++)
        {
            Mat A = leftLapPyr[l].mul(maskGaussianPyramid[l]);
            Mat antiMask = Scalar(1.0, 1.0, 1.0) - maskGaussianPyramid[l];
            Mat B = rightLapPyr[l].mul(antiMask);
            Mat blendedLevel = A + B;

            resultLapPyr.push_back(blendedLevel);
        }
    }

public:
    LaplacianBlending(const Mat& _left, const Mat& _right, const Mat& _blendMask, int _levels) ://construct function, used in LaplacianBlending lb(l,r,m,4);
        leftImg(_left), rightImg(_right), blendMask(_blendMask), levels(_levels)
    {
        assert(_left.size() == _right.size());
        assert(_left.size() == _blendMask.size());
        //创建拉普拉斯金字塔和高斯金字塔
        buildPyramids();
        //每层金字塔图像合并为一个
        blendLapPyrs();
    };

    Mat blend()
    {
        //重建拉普拉斯金字塔
        return reconstructImgFromLapPyramid();
    }
};

Mat LaplacianBlend(const Mat &left, const Mat &right, const Mat &mask)
{
    LaplacianBlending laplaceBlend(left, right, mask, 10);
    return laplaceBlend.blend();
}

int main() {
    Mat leftImg = imread(IMG1_PATH);
    Mat rightImg = imread(IMG2_PATH);

    int hight = leftImg.rows;
    int width = leftImg.cols;

    Mat leftImg32f, rightImg32f;
    leftImg.convertTo(leftImg32f, CV_32F);
    rightImg.convertTo(rightImg32f, CV_32F);

    //创建用于混合的掩膜，这里在中间进行混合
    Mat mask = Mat::zeros(hight, width, CV_32FC1);
    mask(Range::all(), Range(0, mask.cols * 0.5)) = 1.0;

    Mat blendImg = LaplacianBlend(leftImg32f, rightImg32f, mask);

    blendImg.convertTo(blendImg, CV_8UC3);

    imshow("left", leftImg);
    imshow("right", rightImg);
    imshow("blended", blendImg);

    waitKey(0);
    return 0;
}

左图 右图 融合结果

后记

可以看到，融合结果简直惊艳。话说写完这篇博客的时候已经是2019的元旦了，如果你能看到这里，就祝你新年快乐吧~~

Reference:
图像融合之拉普拉斯融合
图像拉普拉斯金字塔融合（Laplacian Pyramid Blending）
【OpenCV入门教程之十三】OpenCV图像金字塔：高斯金字塔、拉普拉斯金字塔与图片尺寸缩放