图像融合之拉普拉斯融合（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实现代码如下：

  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 }

View Code

　融合结果如下图：

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

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

  1 #include <opencv2opencv.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 }

View Code

相关阅读:
What Kind of Friends Are You? ZOJ 3960
博弈随笔(未完待续)
Mergeable Stack ZOJ
LIS ZOJ
差分约束 HDU
How far away ？ HDU
wya费用流
 不知道说些什么
 ext大法好啊
 bzoj2348
原文地址：https://www.cnblogs.com/riddick/p/8922381.html