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

 　　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实现代码如下：

1 #include 
      
         2 #include 
       
          3 #include 
        
           4   5 using namespace std;  6 using namespace cv;  7   8 /************************************************************************/  9 /* 说明： 10 *金字塔从下到上依次为 [0,1，...，level-1] 层 11 *blendMask 为图像的掩模 12 *maskGaussianPyramid为金字塔每一层的掩模 13 *resultLapPyr 存放每层金字塔中直接用左右两图Laplacian变换拼成的图像 14 */ 15 /************************************************************************/ 16  17  18 class LaplacianBlending { 19 private: 20     Mat left; 21     Mat right; 22     Mat blendMask; 23  24     //Laplacian Pyramids 25     vector
         
           leftLapPyr, rightLapPyr, resultLapPyr; 26     Mat leftHighestLevel, rightHighestLevel, resultHighestLevel; 27     //mask为三通道方便矩阵相乘 28     vector
          
            maskGaussianPyramid;  29  30     int levels; 31  32     void buildPyramids()  33     { 34         buildLaplacianPyramid(left, leftLapPyr, leftHighestLevel); 35         buildLaplacianPyramid(right, rightLapPyr, rightHighestLevel); 36         buildGaussianPyramid(); 37     } 38  39     void buildGaussianPyramid()  40     { 41         //金字塔内容为每一层的掩模 42         assert(leftLapPyr.size()>0); 43  44         maskGaussianPyramid.clear(); 45         Mat currentImg; 46         cvtColor(blendMask, currentImg, CV_GRAY2BGR); 47         //保存mask金字塔的每一层图像 48         maskGaussianPyramid.push_back(currentImg); //0-level 49  50         currentImg = blendMask; 51         for (int l = 1; l
           
             l) 54 pyrDown(currentImg, _down, leftLapPyr[l].size()); 55 else 56 pyrDown(currentImg, _down, leftHighestLevel.size()); //lowest level 57 58 Mat down; 59 cvtColor(_down, down, CV_GRAY2BGR); 60 //add color blend mask into mask Pyramid 61 maskGaussianPyramid.push_back(down); 62 currentImg = _down; 63 } 64 } 65 66 void buildLaplacianPyramid(const Mat& img, vector
            
             & lapPyr, Mat& HighestLevel) 67 { 68 lapPyr.clear(); 69 Mat currentImg = img; 70 for (int l = 0; l
             
              = 0; l--) 87 { 88 Mat up; 89 pyrUp(currentImg, up, resultLapPyr[l].size()); 90 currentImg = up + resultLapPyr[l]; 91 } 92 return currentImg; 93 } 94 95 void blendLapPyrs() 96 { 97 //获得每层金字塔中直接用左右两图Laplacian变换拼成的图像resultLapPyr 98 resultHighestLevel = leftHighestLevel.mul(maskGaussianPyramid.back()) + 99 rightHighestLevel.mul(Scalar(1.0, 1.0, 1.0) - maskGaussianPyramid.back());100 for (int l = 0; l
             
            
           
          
         
        
       
      

　融合结果如下图：

  

 

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

  

 

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

1 #include 
      
         2 #include 
       
          3 #include 
        
           4   5 using namespace std;  6   7 //#define DEBUG  8   9  10 void borderInterp(cv::Mat &_src, int radius) 11 { 12     int imgH = _src.rows; 13     int imgW = _src.cols; 14     float *pSrc = (float*)_src.data; 15     for (int i = radius; i < imgH-radius*2; i++) 16     { 17         for (int j = 0; j < 2; j++) 18         { 19             int srcIdx = (i*imgW + j + 3) * 3; 20             int dstIdx = (i*imgW + j) * 3; 21             pSrc[dstIdx] = pSrc[srcIdx]; 22             pSrc[dstIdx + 1] = pSrc[srcIdx + 1]; 23             pSrc[dstIdx + 2] = pSrc[srcIdx + 2]; 24         } 25         for (int j = imgW - radius; j < imgW; j++) 26         { 27             int srcIdx = (i*imgW + j - 3) * 3; 28             int dstIdx = (i*imgW + j) * 3; 29             pSrc[dstIdx] = pSrc[srcIdx]; 30             pSrc[dstIdx + 1] = pSrc[srcIdx + 1]; 31             pSrc[dstIdx + 2] = pSrc[srcIdx + 2]; 32  33         } 34     } 35     for (int j = 0; j < imgW; j++) 36     { 37         for (int i = 0; i < 2; i++) 38         { 39             int srcIdx = ((i + 3)*imgW + j) * 3; 40             int dstIdx = (i*imgW + j) * 3; 41             pSrc[dstIdx] = pSrc[srcIdx]; 42             pSrc[dstIdx + 1] = pSrc[srcIdx + 1]; 43             pSrc[dstIdx + 2] = pSrc[srcIdx + 2]; 44         } 45         for (int i = imgH - radius; i < imgH; i++) 46         { 47             int srcIdx = ((i - 3)*imgW + j) * 3; 48             int dstIdx = (i*imgW + j) * 3; 49             pSrc[dstIdx] = pSrc[srcIdx]; 50             pSrc[dstIdx + 1] = pSrc[srcIdx + 1]; 51             pSrc[dstIdx + 2] = pSrc[srcIdx + 2]; 52         } 53     } 54 } 55  56 void myPyrDown(cv::Mat src, cv::Mat &dst, cv::Size dSize) 57 { 58     dSize = dSize.area() == 0 ? cv::Size((src.cols + 1) / 2, (src.rows + 1) / 2) : dSize; 59  60     float scale = 1. / 16; 61  62     int imgH = src.rows; 63     int imgW = src.cols; 64     cv::Mat _src = cv::Mat::zeros(imgH + 4, imgW + 4, CV_32FC3); 65     int _imgH = _src.rows; 66     int _imgW = _src.cols; 67     src.copyTo(_src(cv::Rect(2, 2, imgW, imgH))); 68     borderInterp(_src, 2); 69  70     //高斯卷积 71     cv::Mat gaussImg = cv::Mat::zeros(imgH, imgW, CV_32FC3); 72     cv::Mat tmpRowGaussImg = _src.clone(); 73     float *pSrc = (float*)_src.data; 74     float *pRowGaussImg = (float*)tmpRowGaussImg.data; 75     //行卷积 76     for (int i = 2; i < imgH+2; i++) 77     { 78         for (int j = 2; j < imgW+2; j++) 79         { 80             float val[3] = { 0 }; 81             int idx = i*_imgW + j; 82             for (int chan = 0; chan < 3; chan++) 83             { 84                 val[chan] += pSrc[(idx - 2) * 3 + chan] + pSrc[(idx + 2) * 3 + chan] 85                     + 4 * (pSrc[(idx - 1) * 3 + chan] + pSrc[(idx + 1) * 3 + chan]) 86                     + 6 * pSrc[idx * 3 + chan]; 87             } 88             pRowGaussImg[idx * 3] = val[0] * scale; 89             pRowGaussImg[idx * 3 + 1] = val[1] * scale; 90             pRowGaussImg[idx * 3 + 2] = val[2] * scale; 91         } 92     } 93  94     float *pGaussImg = (float*)gaussImg.data; 95     //列卷积 96     for (int j = 0; j < imgW; j++) 97     { 98         for (int i = 0; i < imgH; i++) 99         {100             int gi = i + 2;101             int gj = j + 2;102             float val[3] = { 0 };103             int idx = gi*_imgW + gj;104             for (int chan = 0; chan < 3; chan++)105             {106                 val[chan] += pRowGaussImg[(idx-2*_imgW) * 3 + chan] + pRowGaussImg[(idx + 2*_imgW) * 3 + chan]107                     + 4 * (pRowGaussImg[(idx - _imgW) * 3 + chan] + pRowGaussImg[(idx + _imgW) * 3 + chan])108                     + 6 * pRowGaussImg[idx * 3 + chan];109             }110             111             int id = (i*imgW + j) * 3;112             pGaussImg[id] = val[0] * scale;113             pGaussImg[id + 1] = val[1] * scale;114             pGaussImg[id + 2] = val[2] * scale;115         }116     }117 118     int downH = dSize.height;119     int downW = dSize.width;120 121     if (abs(downH * 2 - imgH) > 2) downH = imgH*0.5;122     if (abs(downW * 2 - imgW) > 2) downW = imgW*0.5;123     downH = (downH < 1) ? 1 : downH;124     downW = (downW < 1) ? 1 : downW;125 126     dst = cv::Mat::zeros(downH, downW, CV_32FC3);127     float *pDst = (float*)dst.data;128     for (int i = 0; i < imgH; i++)129     {130         for (int j = 0; j < imgW; j++)131         {132             if (i % 2 != 0 || j % 2 != 0) continue;133             int srcIdx = (i*imgW + j) * 3;134             int y = int((i+1) * 0.5);135             int x = int((j+1) * 0.5);136             y = (y >= downH) ? (downH - 1) : y;137             x = (x >= downW) ? (downW - 1) : x;138             int dstIdx = (y*downW + x) * 3;139             pDst[dstIdx] = pGaussImg[srcIdx];140             pDst[dstIdx + 1] = pGaussImg[srcIdx + 1];141             pDst[dstIdx + 2] = pGaussImg[srcIdx + 2];142         }143     }144 }145 146 void myPyrUp(cv::Mat src, cv::Mat &dst, cv::Size dSize)147 {148     dSize = dSize.area() == 0 ? cv::Size(src.cols * 2, src.rows * 2) : dSize;149     cv::Mat _src;150     src.convertTo(_src, CV_32FC3);151 152     float scale = 1. / 8;153 154     int imgH = src.rows;155     int imgW = src.cols;156     int upImgH = dSize.height;157     int upImgW = dSize.width;158 159     if (abs(upImgH - imgH * 2) > upImgH % 2) upImgH = imgH*2;160     if (abs(upImgW - imgW * 2) > upImgW % 2) upImgW = imgW*2;161 162     cv::Mat upImg = cv::Mat::zeros(upImgH, upImgW, CV_32FC3);163     float *pSrc = (float*)_src.data;164     float *pUpImg = (float*)upImg.data;165     for (int i = 0; i < upImgH; i++)166     {167         for (int j = 0; j < upImgW; j++)168         {169             if (i % 2 != 0 || j % 2 != 0) continue;170             int dstIdx = (i*upImgW + j)*3;171             int y = int((i+1)*0.5);172             int x = int((j+1)*0.5);173             y = (y >= imgH) ? (imgH - 1) : y;174             x = (x >= imgW) ? (imgW - 1) : x;175             int srcIdx = (y*imgW + x) * 3;176 177             pUpImg[dstIdx] = pSrc[srcIdx];178             pUpImg[dstIdx + 1] = pSrc[srcIdx + 1];179             pUpImg[dstIdx + 2] = pSrc[srcIdx + 2];180         }181     }182 183     dst = cv::Mat::zeros(dSize, CV_32FC3);184     cv::Mat _upImg = cv::Mat::zeros(upImgH + 4, upImgW + 4, CV_32FC3);185     int _imgH = _upImg.rows;186     int _imgW = _upImg.cols;187     upImg.copyTo(_upImg(cv::Rect(2, 2, upImgW, upImgH)));188     borderInterp(_upImg, 2);189 190     //高斯卷积191     cv::Mat tempRowGaussImg = _upImg.clone();192     float *pUpData = (float*)_upImg.data;193     float *pRowGaussImg = (float*)tempRowGaussImg.data;194     //行卷积195     for (int i = 2; i < upImgH + 2; i++)196     {197         for (int j = 2; j < upImgW + 2; j++)198         {199             float val[3] = { 0 };200             int idx = i*_imgW + j;201             for (int chan = 0; chan < 3; chan++)202             {203                 val[chan] += pUpData[(idx - 2) * 3 + chan] + pUpData[(idx + 2) * 3 + chan]204                     + 4 * (pUpData[(idx - 1) * 3 + chan] + pUpData[(idx + 1) * 3 + chan])205                     + 6 * pUpData[idx * 3 + chan];206             }207 208             pRowGaussImg[idx * 3] = val[0] * scale;209             pRowGaussImg[idx * 3 + 1] = val[1] * scale;210             pRowGaussImg[idx * 3 + 2] = val[2] * scale;211         }212     }213 214 215     //列卷积216     float *pDst = (float*)dst.data;217     for (int j = 0; j < upImgW; j++)218     {219         for (int i = 0; i < upImgH; i++)220         {221             int gi = i + 2;222             int gj = j + 2;223             float val[3] = { 0 };224             int idx = gi*_imgW + gj;225             for (int chan = 0; chan < 3; chan++)226             {227                 val[chan] += pRowGaussImg[(idx - 2 * _imgW) * 3 + chan] + pRowGaussImg[(idx + 2 * _imgW) * 3 + chan]228                     + 4 * (pRowGaussImg[(idx - _imgW) * 3 + chan] + pRowGaussImg[(idx + _imgW) * 3 + chan])229                     + 6 * pRowGaussImg[idx * 3 + chan];230             }231 232             int id = (i*upImgW + j) * 3;233             pDst[id] = val[0] * scale;234             pDst[id + 1] = val[1] * scale;235             pDst[id + 2] = val[2] * scale;236         }237     }238 }239 240 void buildLaplacianPyramid(cv::Mat srcImg, vector
         
           &pyramidImgs, cv::Mat &topLevelImg, int levels)241 {242     cv::Mat currentImg = srcImg;243     for (int k = 0; k < levels; k++)244     {245         cv::Mat downImg, upImg, lpImg;246 247 #ifdef DEBUG248         cv::pyrDown(currentImg, downImg);249         cv::pyrUp(downImg, upImg, currentImg.size());250 #else251         myPyrDown(currentImg, downImg, cv::Size());252         myPyrUp(downImg, upImg, currentImg.size());253 #endif254 255         lpImg = currentImg - upImg;256         pyramidImgs.push_back(lpImg);257         currentImg = downImg;258     }259     currentImg.copyTo(topLevelImg);260 }261 262 void buildGaussPyramid(cv::Mat mask, vector
          
            &maskGaussPyramidImgs, vector
           
             pyramidImgs,cv::Mat topLevelImg, int levels)263 {264 cv::Mat currentMask;265 //mask转3通道266 if (mask.channels() == 1)267 {268 cv::cvtColor(mask, currentMask, CV_GRAY2BGR);269 }270 else if(mask.channels()==3)271 {272 currentMask = mask;273 }274 275 maskGaussPyramidImgs.push_back(currentMask);276 for (int k = 1; k < levels+1; k++)277 {278 cv::Mat downMask;279 if (k < levels)280 {281 #ifdef DEBUG282 cv::pyrDown(currentMask, downMask, pyramidImgs[k].size());283 #else284 myPyrDown(currentMask, downMask, pyramidImgs[k].size());285 #endif286 }287 else288 {289 #ifdef DEBUG290 cv::pyrDown(currentMask, downMask, topLevelImg.size());291 #else292 myPyrDown(currentMask, downMask, topLevelImg.size());293 #endif294 }295 296 maskGaussPyramidImgs.push_back(downMask);297 currentMask = downMask;298 }299 }300 301 void buildResultPyramid(vector
            
              leftPyramidImgs, vector
             
               rightPyramidImgs, vector
              
                maskPyramids, vector
               
                 &resultPyramidImgs, int levels)302 {303 304 for (int k = 0; k < levels; k++)305 {306 cv::Mat left = leftPyramidImgs[k].mul(maskPyramids[k]);307 cv::Mat right = rightPyramidImgs[k].mul(cv::Scalar(1.0,1.0,1.0) - maskPyramids[k]);308 cv::Mat result = left + right;309 resultPyramidImgs.push_back(result);310 }311 312 }313 314 void reconstruct(vector
                
                  lpPyramidImgs, cv::Mat blendTopLevelImg, cv::Mat &blendImg, int levels)315 {316 cv::Mat currentImg = blendTopLevelImg;317 for (int k = levels - 1; k >= 0; k--)318 {319 cv::Mat upImg;320 #ifdef DEBUG321 cv::pyrUp(currentImg, upImg, lpPyramidImgs[k].size());322 #else323 myPyrUp(currentImg, upImg, lpPyramidImgs[k].size());324 #endif325 currentImg = upImg + lpPyramidImgs[k];326 }327 currentImg.copyTo(blendImg);328 }329 330 cv::Mat laplacianBlending(cv::Mat leftImg, cv::Mat rightImg, cv::Mat mask)331 {332 cv::Mat leftImg32f, rightImg32f, mask32f;333 leftImg.convertTo(leftImg32f, CV_32FC1);334 rightImg.convertTo(rightImg32f, CV_32FC1);335 mask.convertTo(mask32f, CV_32FC1);336 337 vector
                 
                   leftLpPyramidImgs, rightLpPyramidImgs, resultLpPyramidImgs, gaussPyramidMaskImgs;338 cv::Mat leftTopLevelImg, rightTopLevelImg;339 int levels =4;340 //拉普拉斯金字塔341 buildLaplacianPyramid(leftImg32f, leftLpPyramidImgs, leftTopLevelImg, levels);342 buildLaplacianPyramid(rightImg32f, rightLpPyramidImgs, rightTopLevelImg, levels);343 //mask创建gauss金字塔344 buildGaussPyramid(mask32f, gaussPyramidMaskImgs, leftLpPyramidImgs, leftTopLevelImg, levels);345 //结合左右两图的laplacian残差图346 buildResultPyramid(leftLpPyramidImgs, rightLpPyramidImgs, gaussPyramidMaskImgs, resultLpPyramidImgs, levels);347 //348 cv::Mat blendImg = cv::Mat::zeros(leftImg.size(), CV_32FC3);349 350 cv::Mat blendTopLevelImg = leftTopLevelImg.mul(gaussPyramidMaskImgs[levels]) + rightTopLevelImg.mul(cv::Scalar(1.0, 1.0, 1.0) - gaussPyramidMaskImgs[levels]);351 reconstruct(resultLpPyramidImgs, blendTopLevelImg, blendImg, levels);352 353 blendImg.convertTo(blendImg, CV_8UC3);354 return blendImg;355 }356 357 void main()358 {359 cv::Mat appleImg = cv::imread("data/apple.jpg");360 cv::Mat pearImg = cv::imread("data/orange.jpg");361 362 int imgH = appleImg.rows;363 int imgW = appleImg.cols;364 cv::Mat mask = cv::Mat::zeros(imgH, imgW, CV_32FC1);365 mask(cv::Range::all(), cv::Range(0, imgW*0.5)) = 1.0;366 cv::Mat blendImg = laplacianBlending(appleImg, pearImg, mask);367 cv::namedWindow("blendImg", 0);368 cv::imshow("blendImg", blendImg);369 cv::imwrite("data/blendImg.png", blendImg);370 cv::waitKey(0);371 }