一、拉普拉斯融合基本步骤
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 <opencv2/opencv.hpp> 2 #include <iostream> 3 #include <string> 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<Mat> leftLapPyr, rightLapPyr, resultLapPyr; 26 Mat leftHighestLevel, rightHighestLevel, resultHighestLevel; 27 //mask为三通道方便矩阵相乘 28 vector<Mat> 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<levels + 1; l++) { 52 Mat _down; 53 if (leftLapPyr.size() > 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<Mat>& lapPyr, Mat& HighestLevel) 67 { 68 lapPyr.clear(); 69 Mat currentImg = img; 70 for (int l = 0; l<levels; l++) { 71 Mat down, up; 72 pyrDown(currentImg, down); 73 pyrUp(down, up, currentImg.size()); 74 Mat lap = currentImg - up; 75 lapPyr.push_back(lap); 76 currentImg = down; 77 } 78 currentImg.copyTo(HighestLevel); 79 } 80 81 Mat reconstructImgFromLapPyramid() 82 { 83 //将左右laplacian图像拼成的resultLapPyr金字塔中每一层 84 //从上到下插值放大并与残差相加，即得blend图像结果 85 Mat currentImg = resultHighestLevel; 86 for (int l = levels - 1; 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<levels; l++) 101 {102 Mat A = leftLapPyr[l].mul(maskGaussianPyramid[l]);103 Mat antiMask = Scalar(1.0, 1.0, 1.0) - maskGaussianPyramid[l];104 Mat B = rightLapPyr[l].mul(antiMask);105 Mat blendedLevel = A + B;106 107 resultLapPyr.push_back(blendedLevel);108 }109 }110 111 public:112 LaplacianBlending(const Mat& _left, const Mat& _right, const Mat& _blendMask, int _levels) ://construct function, used in LaplacianBlending lb(l,r,m,4);113 left(_left), right(_right), blendMask(_blendMask), levels(_levels)114 {115 assert(_left.size() == _right.size());116 assert(_left.size() == _blendMask.size());117 //创建拉普拉斯金字塔和高斯金字塔118 buildPyramids();119 //每层金字塔图像合并为一个120 blendLapPyrs();121 };122 123 Mat blend() 124 {125 //重建拉普拉斯金字塔126 return reconstructImgFromLapPyramid();127 }128 };129 130 Mat LaplacianBlend(const Mat &left, const Mat &right, const Mat &mask) 131 {132 LaplacianBlending laplaceBlend(left, right, mask, 10);133 return laplaceBlend.blend();134 }135 136 int main() {137 Mat img8UL = imread("data/apple.jpg");138 Mat img8UR = imread("data/orange.jpg");139 140 int imgH = img8UL.rows;141 int imgW = img8UL.cols;142 143 imshow("left", img8UL);144 imshow("right", img8UR);145 146 Mat img32fL, img32fR; 147 img8UL.convertTo(img32fL, CV_32F);148 img8UR.convertTo(img32fR, CV_32F);149 150 //创建mask151 Mat mask = Mat::zeros(imgH, imgW, CV_32FC1);152 mask(Range::all(), Range(0, mask.cols * 0.5)) = 1.0;153 154 Mat blendImg = LaplacianBlend(img32fL, img32fR, mask);155 156 blendImg.convertTo(blendImg, CV_8UC3);157 imshow("blended", blendImg);158 159 waitKey(0);160 return 0;161 }融合结果如下图：
金字塔层数level=5 金字塔层数level=10
附上自己实现pyrDown和pyrUp写的拉普拉斯融合，仅供参考：
1 #include <opencv2\opencv.hpp> 2 #include <iostream> 3 #include <string> 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<cv::Mat> &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<cv::Mat> &maskGaussPyramidImgs, vector<cv::Mat> 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<cv::Mat> leftPyramidImgs, vector<cv::Mat> rightPyramidImgs, vector<cv::Mat> maskPyramids, vector<cv::Mat> &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<cv::Mat> 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<cv::Mat> 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 }
联系客服
