图像融合之拉普拉斯融合(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实现代码如下:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
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 //创建mask 151 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写的拉普拉斯融合,仅供参考:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
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 DEBUG 248 cv::pyrDown(currentImg, downImg); 249 cv::pyrUp(downImg, upImg, currentImg.size()); 250 #else 251 myPyrDown(currentImg, downImg, cv::Size()); 252 myPyrUp(downImg, upImg, currentImg.size()); 253 #endif 254 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 DEBUG 282 cv::pyrDown(currentMask, downMask, pyramidImgs[k].size()); 283 #else 284 myPyrDown(currentMask, downMask, pyramidImgs[k].size()); 285 #endif 286 } 287 else 288 { 289 #ifdef DEBUG 290 cv::pyrDown(currentMask, downMask, topLevelImg.size()); 291 #else 292 myPyrDown(currentMask, downMask, topLevelImg.size()); 293 #endif 294 } 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 DEBUG 321 cv::pyrUp(currentImg, upImg, lpPyramidImgs[k].size()); 322 #else 323 myPyrUp(currentImg, upImg, lpPyramidImgs[k].size()); 324 #endif 325 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 }