一、拉普拉斯融合基本步骤
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 //创建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写的拉普拉斯融合,仅供参考:
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 }