影象融合之拉普拉斯融合（laplacian blending）
阿新 • • 發佈：2019-01-17
一、拉普拉斯融合基本步驟
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 namespaceView Codestd; 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 <