图像融合之拉普拉斯融合（laplacian blending）

一、拉普拉斯融合基本步骤

 　　1. 两幅图像L，R，以及二值掩模mask，给定金字塔层数level。

　　2. 分别根据L，R构建其对应的拉普拉斯残差金字塔（层数为level），并保留高斯金字塔下采样最顶端的图像（尺寸最小的图像，第level+1层）：

　　　　拉普拉斯残差金字塔构建方法如下，以L图为例：

　　　　(1) 对L进行高斯下采样得到downL，OpenCV中pyrDown()函数可以实现此功能。然后再对downL进行高斯上采样得到upL，OpenCV中pyrUp()函数可以实现此功能。

　　　　(2) 计算原图L与upL之间的残差，得到一幅残差图lapL0。作为残差金字塔最低端的图像。

　　　　(3) 对downL继续进行(1) (2)操作，不断计算残差图lapL1, lap2, lap3.....lapN。这样得到一系列残差图，即为拉普拉斯残差金字塔。

　　　　(4)拉普拉斯 残差金字塔中一共有level幅图像。而我们需要保留第level+1层的高斯下采样图topL，以便后面使用。

　　3. 二值掩模mask下采样构建高斯金字塔，同样利用pyrDown()实现，共有level+1层。

　　4. 利用mask金字塔每一层的mask图，将L图和R图的拉普拉斯残差金字塔对应层的图像合并为一幅图像。这样得到合并后的拉普拉斯残差金字塔。同时利用最顶端的mask将步骤2中保留的topL和topR合并为topLR。

　　5. 以topLR为金字塔最顶端的图像，利用pyrUp()函数对topLR进行高斯上采样，得到upTopLR，并将upTopLR与步骤4中合并后的残差金字塔对应层的图像相加，重建出该层的图像。

　　6. 重复步骤5，直至重建出第0层，也就是金字塔最低端的图像，即blendImg。输出。

 

 二、代码

　　拉普拉斯融合的OpenCV实现代码如下：

#include <opencv2/opencv.hpp>
#include <iostream>
#include <string>

using namespace std;
using namespace cv;

/************************************************************************/
/* 说明：
*金字塔从下到上依次为 [0,1，...，level-1] 层
*blendMask 为图像的掩模
*maskGaussianPyramid为金字塔每一层的掩模
*resultLapPyr 存放每层金字塔中直接用左右两图Laplacian变换拼成的图像
*/
/************************************************************************/


class LaplacianBlending {
private:
    Mat left;
    Mat right;
    Mat blendMask;

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

    int levels;

    void buildPyramids() 
    {
        buildLaplacianPyramid(left, leftLapPyr, leftHighestLevel);
        buildLaplacianPyramid(right, 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);
            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);
        left(_left), right(_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 img8UL = imread("data/apple.jpg");
    Mat img8UR = imread("data/orange.jpg");

    int imgH = img8UL.rows;
    int imgW = img8UL.cols;

    imshow("left", img8UL);
    imshow("right", img8UR);

    Mat img32fL, img32fR; 
    img8UL.convertTo(img32fL, CV_32F);
    img8UR.convertTo(img32fR, CV_32F);

    //创建mask
    Mat mask = Mat::zeros(imgH, imgW, CV_32FC1);
    mask(Range::all(), Range(0, mask.cols * 0.5)) = 1.0;

    Mat blendImg = LaplacianBlend(img32fL, img32fR, mask);

    blendImg.convertTo(blendImg, CV_8UC3);
    imshow("blended", blendImg);

    waitKey(0);
    return 0;
}

　融合结果如下图：

  

 

金字塔层数level=5                                                金字塔层数level=10

  

 

　　附上自己实现pyrDown和pyrUp写的拉普拉斯融合，仅供参考：

#include <opencv2\opencv.hpp>
#include <iostream>
#include <string>

using namespace std;

//#define DEBUG


void borderInterp(cv::Mat &_src, int radius)
{
    int imgH = _src.rows;
    int imgW = _src.cols;
    float *pSrc = (float*)_src.data;
    for (int i = radius; i < imgH-radius*2; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            int srcIdx = (i*imgW + j + 3) * 3;
            int dstIdx = (i*imgW + j) * 3;
            pSrc[dstIdx] = pSrc[srcIdx];
            pSrc[dstIdx + 1] = pSrc[srcIdx + 1];
            pSrc[dstIdx + 2] = pSrc[srcIdx + 2];
        }
        for (int j = imgW - radius; j < imgW; j++)
        {
            int srcIdx = (i*imgW + j - 3) * 3;
            int dstIdx = (i*imgW + j) * 3;
            pSrc[dstIdx] = pSrc[srcIdx];
            pSrc[dstIdx + 1] = pSrc[srcIdx + 1];
            pSrc[dstIdx + 2] = pSrc[srcIdx + 2];

        }
    }
    for (int j = 0; j < imgW; j++)
    {
        for (int i = 0; i < 2; i++)
        {
            int srcIdx = ((i + 3)*imgW + j) * 3;
            int dstIdx = (i*imgW + j) * 3;
            pSrc[dstIdx] = pSrc[srcIdx];
            pSrc[dstIdx + 1] = pSrc[srcIdx + 1];
            pSrc[dstIdx + 2] = pSrc[srcIdx + 2];
        }
        for (int i = imgH - radius; i < imgH; i++)
        {
            int srcIdx = ((i - 3)*imgW + j) * 3;
            int dstIdx = (i*imgW + j) * 3;
            pSrc[dstIdx] = pSrc[srcIdx];
            pSrc[dstIdx + 1] = pSrc[srcIdx + 1];
            pSrc[dstIdx + 2] = pSrc[srcIdx + 2];
        }
    }
}

void myPyrDown(cv::Mat src, cv::Mat &dst, cv::Size dSize)
{
    dSize = dSize.area() == 0 ? cv::Size((src.cols + 1) / 2, (src.rows + 1) / 2) : dSize;

    float scale = 1. / 16;

    int imgH = src.rows;
    int imgW = src.cols;
    cv::Mat _src = cv::Mat::zeros(imgH + 4, imgW + 4, CV_32FC3);
    int _imgH = _src.rows;
    int _imgW = _src.cols;
    src.copyTo(_src(cv::Rect(2, 2, imgW, imgH)));
    borderInterp(_src, 2);

    //高斯卷积
    cv::Mat gaussImg = cv::Mat::zeros(imgH, imgW, CV_32FC3);
    cv::Mat tmpRowGaussImg = _src.clone();
    float *pSrc = (float*)_src.data;
    float *pRowGaussImg = (float*)tmpRowGaussImg.data;
    //行卷积
    for (int i = 2; i < imgH+2; i++)
    {
        for (int j = 2; j < imgW+2; j++)
        {
            float val[3] = { 0 };
            int idx = i*_imgW + j;
            for (int chan = 0; chan < 3; chan++)
            {
                val[chan] += pSrc[(idx - 2) * 3 + chan] + pSrc[(idx + 2) * 3 + chan]
                    + 4 * (pSrc[(idx - 1) * 3 + chan] + pSrc[(idx + 1) * 3 + chan])
                    + 6 * pSrc[idx * 3 + chan];
            }
            pRowGaussImg[idx * 3] = val[0] * scale;
            pRowGaussImg[idx * 3 + 1] = val[1] * scale;
            pRowGaussImg[idx * 3 + 2] = val[2] * scale;
        }
    }

    float *pGaussImg = (float*)gaussImg.data;
    //列卷积
    for (int j = 0; j < imgW; j++)
    {
        for (int i = 0; i < imgH; i++)
        {
            int gi = i + 2;
            int gj = j + 2;
            float val[3] = { 0 };
            int idx = gi*_imgW + gj;
            for (int chan = 0; chan < 3; chan++)
            {
                val[chan] += pRowGaussImg[(idx-2*_imgW) * 3 + chan] + pRowGaussImg[(idx + 2*_imgW) * 3 + chan]
                    + 4 * (pRowGaussImg[(idx - _imgW) * 3 + chan] + pRowGaussImg[(idx + _imgW) * 3 + chan])
                    + 6 * pRowGaussImg[idx * 3 + chan];
            }
            
            int id = (i*imgW + j) * 3;
            pGaussImg[id] = val[0] * scale;
            pGaussImg[id + 1] = val[1] * scale;
            pGaussImg[id + 2] = val[2] * scale;
        }
    }

    int downH = dSize.height;
    int downW = dSize.width;

    if (abs(downH * 2 - imgH) > 2) downH = imgH*0.5;
    if (abs(downW * 2 - imgW) > 2) downW = imgW*0.5;
    downH = (downH < 1) ? 1 : downH;
    downW = (downW < 1) ? 1 : downW;

    dst = cv::Mat::zeros(downH, downW, CV_32FC3);
    float *pDst = (float*)dst.data;
    for (int i = 0; i < imgH; i++)
    {
        for (int j = 0; j < imgW; j++)
        {
            if (i % 2 != 0 || j % 2 != 0) continue;
            int srcIdx = (i*imgW + j) * 3;
            int y = int((i+1) * 0.5);
            int x = int((j+1) * 0.5);
            y = (y >= downH) ? (downH - 1) : y;
            x = (x >= downW) ? (downW - 1) : x;
            int dstIdx = (y*downW + x) * 3;
            pDst[dstIdx] = pGaussImg[srcIdx];
            pDst[dstIdx + 1] = pGaussImg[srcIdx + 1];
            pDst[dstIdx + 2] = pGaussImg[srcIdx + 2];
        }
    }
}

void myPyrUp(cv::Mat src, cv::Mat &dst, cv::Size dSize)
{
    dSize = dSize.area() == 0 ? cv::Size(src.cols * 2, src.rows * 2) : dSize;
    cv::Mat _src;
    src.convertTo(_src, CV_32FC3);

    float scale = 1. / 8;

    int imgH = src.rows;
    int imgW = src.cols;
    int upImgH = dSize.height;
    int upImgW = dSize.width;

    if (abs(upImgH - imgH * 2) > upImgH % 2) upImgH = imgH*2;
    if (abs(upImgW - imgW * 2) > upImgW % 2) upImgW = imgW*2;

    cv::Mat upImg = cv::Mat::zeros(upImgH, upImgW, CV_32FC3);
    float *pSrc = (float*)_src.data;
    float *pUpImg = (float*)upImg.data;
    for (int i = 0; i < upImgH; i++)
    {
        for (int j = 0; j < upImgW; j++)
        {
            if (i % 2 != 0 || j % 2 != 0) continue;
            int dstIdx = (i*upImgW + j)*3;
            int y = int((i+1)*0.5);
            int x = int((j+1)*0.5);
            y = (y >= imgH) ? (imgH - 1) : y;
            x = (x >= imgW) ? (imgW - 1) : x;
            int srcIdx = (y*imgW + x) * 3;

            pUpImg[dstIdx] = pSrc[srcIdx];
            pUpImg[dstIdx + 1] = pSrc[srcIdx + 1];
            pUpImg[dstIdx + 2] = pSrc[srcIdx + 2];
        }
    }

    dst = cv::Mat::zeros(dSize, CV_32FC3);
    cv::Mat _upImg = cv::Mat::zeros(upImgH + 4, upImgW + 4, CV_32FC3);
    int _imgH = _upImg.rows;
    int _imgW = _upImg.cols;
    upImg.copyTo(_upImg(cv::Rect(2, 2, upImgW, upImgH)));
    borderInterp(_upImg, 2);

    //高斯卷积
    cv::Mat tempRowGaussImg = _upImg.clone();
    float *pUpData = (float*)_upImg.data;
    float *pRowGaussImg = (float*)tempRowGaussImg.data;
    //行卷积
    for (int i = 2; i < upImgH + 2; i++)
    {
        for (int j = 2; j < upImgW + 2; j++)
        {
            float val[3] = { 0 };
            int idx = i*_imgW + j;
            for (int chan = 0; chan < 3; chan++)
            {
                val[chan] += pUpData[(idx - 2) * 3 + chan] + pUpData[(idx + 2) * 3 + chan]
                    + 4 * (pUpData[(idx - 1) * 3 + chan] + pUpData[(idx + 1) * 3 + chan])
                    + 6 * pUpData[idx * 3 + chan];
            }

            pRowGaussImg[idx * 3] = val[0] * scale;
            pRowGaussImg[idx * 3 + 1] = val[1] * scale;
            pRowGaussImg[idx * 3 + 2] = val[2] * scale;
        }
    }


    //列卷积
    float *pDst = (float*)dst.data;
    for (int j = 0; j < upImgW; j++)
    {
        for (int i = 0; i < upImgH; i++)
        {
            int gi = i + 2;
            int gj = j + 2;
            float val[3] = { 0 };
            int idx = gi*_imgW + gj;
            for (int chan = 0; chan < 3; chan++)
            {
                val[chan] += pRowGaussImg[(idx - 2 * _imgW) * 3 + chan] + pRowGaussImg[(idx + 2 * _imgW) * 3 + chan]
                    + 4 * (pRowGaussImg[(idx - _imgW) * 3 + chan] + pRowGaussImg[(idx + _imgW) * 3 + chan])
                    + 6 * pRowGaussImg[idx * 3 + chan];
            }

            int id = (i*upImgW + j) * 3;
            pDst[id] = val[0] * scale;
            pDst[id + 1] = val[1] * scale;
            pDst[id + 2] = val[2] * scale;
        }
    }
}

void buildLaplacianPyramid(cv::Mat srcImg, vector<cv::Mat> &pyramidImgs, cv::Mat &topLevelImg, int levels)
{
    cv::Mat currentImg = srcImg;
    for (int k = 0; k < levels; k++)
    {
        cv::Mat downImg, upImg, lpImg;

#ifdef DEBUG
        cv::pyrDown(currentImg, downImg);
        cv::pyrUp(downImg, upImg, currentImg.size());
#else
        myPyrDown(currentImg, downImg, cv::Size());
        myPyrUp(downImg, upImg, currentImg.size());
#endif

        lpImg = currentImg - upImg;
        pyramidImgs.push_back(lpImg);
        currentImg = downImg;
    }
    currentImg.copyTo(topLevelImg);
}

void buildGaussPyramid(cv::Mat mask, vector<cv::Mat> &maskGaussPyramidImgs, vector<cv::Mat> pyramidImgs,cv::Mat topLevelImg, int levels)
{
    cv::Mat currentMask;
    //mask转3通道
    if (mask.channels() == 1)
    {
        cv::cvtColor(mask, currentMask, CV_GRAY2BGR);
    }
    else if(mask.channels()==3)
    {
        currentMask = mask;
    }
    
    maskGaussPyramidImgs.push_back(currentMask);
    for (int k = 1; k < levels+1; k++)
    {
        cv::Mat downMask;
        if (k < levels)
        {
#ifdef DEBUG
            cv::pyrDown(currentMask, downMask, pyramidImgs[k].size());
#else
            myPyrDown(currentMask, downMask, pyramidImgs[k].size());
#endif
        }
        else
        {
#ifdef DEBUG
            cv::pyrDown(currentMask, downMask, topLevelImg.size());
#else
            myPyrDown(currentMask, downMask, topLevelImg.size());
#endif
        }

        maskGaussPyramidImgs.push_back(downMask);
        currentMask = downMask;
    }
}

void buildResultPyramid(vector<cv::Mat> leftPyramidImgs, vector<cv::Mat> rightPyramidImgs, vector<cv::Mat> maskPyramids, vector<cv::Mat> &resultPyramidImgs, int levels)
{

    for (int k = 0; k < levels; k++)
    {
        cv::Mat left = leftPyramidImgs[k].mul(maskPyramids[k]);
        cv::Mat right = rightPyramidImgs[k].mul(cv::Scalar(1.0,1.0,1.0) - maskPyramids[k]);
        cv::Mat result = left + right;
        resultPyramidImgs.push_back(result);
    }

}

void reconstruct(vector<cv::Mat> lpPyramidImgs, cv::Mat blendTopLevelImg, cv::Mat &blendImg, int levels)
{
    cv::Mat currentImg = blendTopLevelImg;
    for (int k = levels - 1; k >= 0; k--)
    {
        cv::Mat upImg;
#ifdef DEBUG
        cv::pyrUp(currentImg, upImg, lpPyramidImgs[k].size());
#else
        myPyrUp(currentImg, upImg, lpPyramidImgs[k].size());
#endif
        currentImg = upImg + lpPyramidImgs[k];
    }
    currentImg.copyTo(blendImg);
}

cv::Mat laplacianBlending(cv::Mat leftImg, cv::Mat rightImg, cv::Mat mask)
{
    cv::Mat leftImg32f, rightImg32f, mask32f;
    leftImg.convertTo(leftImg32f, CV_32FC1);
    rightImg.convertTo(rightImg32f, CV_32FC1);
    mask.convertTo(mask32f, CV_32FC1);

    vector<cv::Mat> leftLpPyramidImgs, rightLpPyramidImgs, resultLpPyramidImgs, gaussPyramidMaskImgs;
    cv::Mat leftTopLevelImg, rightTopLevelImg;
    int levels =4;
    //拉普拉斯金字塔
    buildLaplacianPyramid(leftImg32f, leftLpPyramidImgs, leftTopLevelImg, levels);
    buildLaplacianPyramid(rightImg32f, rightLpPyramidImgs, rightTopLevelImg, levels);
    //mask创建gauss金字塔
    buildGaussPyramid(mask32f, gaussPyramidMaskImgs, leftLpPyramidImgs, leftTopLevelImg, levels);
    //结合左右两图的laplacian残差图
    buildResultPyramid(leftLpPyramidImgs, rightLpPyramidImgs, gaussPyramidMaskImgs, resultLpPyramidImgs, levels);
    //
    cv::Mat blendImg = cv::Mat::zeros(leftImg.size(), CV_32FC3);

    cv::Mat blendTopLevelImg = leftTopLevelImg.mul(gaussPyramidMaskImgs[levels]) + rightTopLevelImg.mul(cv::Scalar(1.0, 1.0, 1.0) - gaussPyramidMaskImgs[levels]);
    reconstruct(resultLpPyramidImgs, blendTopLevelImg, blendImg, levels);

    blendImg.convertTo(blendImg, CV_8UC3);
    return blendImg;
}

void main()
{
    cv::Mat appleImg = cv::imread("data/apple.jpg");
    cv::Mat pearImg = cv::imread("data/orange.jpg");

    int imgH = appleImg.rows;
    int imgW = appleImg.cols;
    cv::Mat mask = cv::Mat::zeros(imgH, imgW, CV_32FC1);
    mask(cv::Range::all(), cv::Range(0, imgW*0.5)) = 1.0;
    cv::Mat blendImg = laplacianBlending(appleImg, pearImg, mask);
    cv::namedWindow("blendImg", 0);
    cv::imshow("blendImg", blendImg);
    cv::imwrite("data/blendImg.png", blendImg);
    cv::waitKey(0);
}