Poisson blending is one of the topics that spent me days trying to understand recently (not fully understand yet), it is a wonderful method, and using wonderful maths, as well. Now I’ll try to explain this method, with as less Math formulae as I can.
ABOUT IMAGE BLENDING
As an image blending method, our goal is to seamlessly blend an object or texture from a source image into a target image, here is an example:
we are going to put the dog (or bear??) in 1st image, the two children in 2nd image into the regions in the 3rd image. The simplest way to do that is copy and paste the pixels from 1st and 2nd images into the 3rd image, the result will like this:
Weird, we can see very noticeable difference in these regions, because the color of water in the three sources are different, even if the backgrounds are matched, we can also see seams between these regions. Why?
Our people sight is much more sensitive to gradient than the overall intensity of an image, that means, we don’t know the exact value of each pixel, but if value changes somewhere, we can easily know that it changes, and this in fact leads to our goal: clone pixels to another image, but let color not change abruptly, i.e. maximally preserve the gradient of the source region. This result is much better than the above one:
POISSON BLENDING OVERVIEW
Poisson blending is one of gradient domain image processing methods, maybe this figure can explain it well:
- v – Gradient of a region in an image (as vector)
- g – Selected region of source (dog, children…)
- f* – Known functions that exist in domain S (water in 3rd source above)
- f – Unknown functions that exist in domain Ω
- Ω – Region g that is now placed on domain S (target background)
- ∂Ω – Boundaries between the source and target regions
So Actually what we are going to do is nothing but given vector field v, find the value of f in unknown region that optimize.
whose solution is the solution of Poisson equation:
where div v = ∂v/∂x + ∂v/∂y, is the Divergence Property of gradient field, and ∆ is the Laplacian operator:
Applying the Laplacian operator to the Poisson equation, we can get:
div G = -4 f(x,y) + f(x-1,y) + f(x,y-1) + f(x+1,y) + f(x,y+1)
1-D EXAMPLE
Say we are going to copy the red part of left image (4, 3, 5, 4) into the right image.
Fist, we got the gradient of left image, as showed above the red bars, now we are going to move these pixels to the right hand side image, without abruptly change of value, so we want:
With f1 = 6, f6 = 1.
that is,
Let’s denote it Q, so we can get:
we can solve it by using matrix:
f2 = 6, f3 = 4, f4 = 5, f5 = 3
We can find that, the matrix above is a convolution of kernel (-2 4 -2), because every pixel value is influenced by the two adjacent pixels, this is very similar to the Laplacian operator in 2D situation, every pixel value is influenced by the four adjacent pixels (up, down, left, right).
2-D EXAMPLE
Say we are going to copy some pixels into the white areas of the above image (red pixels are boundary pixels), we doing that by using the equation as same as 1D situation:
Ax = b
in which A is a N*N matrix, N is the number of pixels we are going to copy, in this case, N=10.
We can create the 10*10 matrix as this:
05 | elif(adjacent(pixel(j), pixel(i))) |
So the matrix in this case is like:
In the equation, b is a N-elements vector, and is created by:
b[i] = div ( G( Source(x,y) ) ) + Neighbor(target i) ; i=1..N
div v can be calculated by the formula in front of this article, and Neighbor means pixel i’s 4 neighbors which belong to the boundary, for example:
Neighbor(pixel 1) are pixels left, up, right;
Neighbor(pixel 8) are pixels left, up;
Neighbor(pixel 5) is pixels up.
and if we know matrix A and vector b, we can calculate vector x, which are the value of pixels in target image inside.
x = A \ b (‘\’ means matrix left division)
or, x = A-1 * b
MORE ABOUT POISSON SOLVER
Prof. Gilbert Strang in MIT tells more about FAST POISSON SOLVER, but this only works when the region we want to copy is a square, i.e. with same height and width.
In the above case, the matrix A (as Gilbert’s interpretation, matrix K2D) is in a very beautiful form that we can use the same method of FFT to deal with it, and it can save a lot of time. More details is here:
Lecture 25: Fast Poisson Solver (part 1)
Lecture 26: Fast Poisson Solver (part 2)
TEST RESULTS
Since Chris Martin has charming eyes, so I copied his eyes to everyone in his band…
C++ CODE
013 | void display( const char *name, IplImage *img){ |
015 | cvShowImage(name,img); |
018 | void getGradientx( const IplImage *img, IplImage *gx){ |
021 | int channel = img->nChannels; |
026 | for ( int c=0;c<channel;++c){ |
027 | CV_IMAGE_ELEM(gx, float ,i,j*channel+c) = |
028 | ( float )CV_IMAGE_ELEM(img,uchar,i,(j+1)*channel+c) - ( float )CV_IMAGE_ELEM(img,uchar,i,j*channel+c); |
031 | void getGradienty( const IplImage *img, IplImage *gy){ |
034 | int channel = img->nChannels; |
039 | for ( int c=0;c<channel;++c){ |
040 | CV_IMAGE_ELEM(gy, float ,i,j*channel+c) = |
041 | ( float )CV_IMAGE_ELEM(img,uchar,(i+1),j*channel+c) - ( float )CV_IMAGE_ELEM(img,uchar,i,j*channel+c); |
045 | void lapx( const IplImage *img, IplImage *gxx){ |
048 | int channel = img->nChannels; |
052 | for ( int j=0;j<w-1;j++) |
053 | for ( int c=0;c<channel;++c){ |
054 | CV_IMAGE_ELEM(gxx, float ,i,(j+1)*channel+c) = |
055 | ( float )CV_IMAGE_ELEM(img, float ,i,(j+1)*channel+c) - ( float )CV_IMAGE_ELEM(img, float ,i,j*channel+c); |
058 | void lapy( const IplImage *img, IplImage *gyy){ |
061 | int channel = img->nChannels; |
063 | for ( int i=0;i<h-1;i++) |
065 | for ( int c=0;c<channel;++c){ |
066 | CV_IMAGE_ELEM(gyy, float ,i+1,j*channel+c) = |
067 | ( float )CV_IMAGE_ELEM(img, float ,(i+1),j*channel+c) - ( float )CV_IMAGE_ELEM(img, float ,i,j*channel+c); |
072 | void dst( double *gtest, double *gfinal, int h, int w){ |
075 | unsigned long int idx; |
077 | Mat temp = Mat(2*h+2,1,CV_32F); |
078 | Mat res = Mat(h,1,CV_32F); |
081 | for ( int i=0;i<w;i++){ |
082 | temp.at< float >(0,0) = 0.0; |
084 | for ( int j=0,r=1;j<h;j++,r++){ |
086 | temp.at< float >(r,0) = gtest[idx]; |
089 | temp.at< float >(h+1,0)=0.0; |
091 | for ( int j=h-1, r=h+2;j>=0;j--,r++){ |
093 | temp.at< float >(r,0) = -1*gtest[idx]; |
096 | Mat planes[] = {Mat_< float >(temp), Mat::zeros(temp.size(), CV_32F)}; |
099 | merge(planes, 2, complex1); |
101 | dft(complex1,complex1,0,0); |
103 | Mat planes1[] = {Mat::zeros(complex1.size(), CV_32F), Mat::zeros(complex1.size(), CV_32F)}; |
106 | split(complex1, planes1); |
108 | std::complex< double > two_i = std:: sqrt (std::complex< double >(-1)); |
110 | double fac = -2*imag(two_i); |
112 | for ( int c=1,z=0;c<h+1;c++,z++){ |
113 | res.at< float >(z,0) = planes1[1].at< float >(c,0)/fac; |
116 | for ( int q=0,z=0;q<h;q++,z++){ |
118 | gfinal[idx] = res.at< float >(z,0); |
125 | void idst( double *gtest, double *gfinal, int h, int w){ |
127 | unsigned long int idx; |
128 | dst(gtest,gfinal,h,w); |
129 | for ( int i= 0;i<h;i++) |
130 | for ( int j=0;j<w;j++){ |
132 | gfinal[idx] = ( double ) (2*gfinal[idx])/nn; |
136 | void transpose( double *mat, double *mat_t, int h, int w){ |
138 | Mat tmp = Mat(h,w,CV_32FC1); |
140 | unsigned long int idx; |
141 | for ( int i = 0 ; i < h;i++){ |
142 | for ( int j = 0 ; j < w; j++){ |
145 | tmp.at< float >(i,j) = mat[idx]; |
150 | for ( int i = 0;i < tmp_t.size().height; i++) |
151 | for ( int j=0;j<tmp_t.size().width;j++){ |
152 | idx = i*tmp_t.size().width + j; |
153 | mat_t[idx] = tmp_t.at< float >(i,j); |
157 | void poisson_solver( const IplImage *img, IplImage *gxx , IplImage *gyy, Mat &result){ |
161 | int channel = img->nChannels; |
163 | unsigned long int idx,idx1; |
165 | IplImage *lap = cvCreateImage(cvGetSize(img), 32, 1); |
167 | for ( int i =0;i<h;i++) |
169 | CV_IMAGE_ELEM(lap, float ,i,j)=CV_IMAGE_ELEM(gyy, float ,i,j)+CV_IMAGE_ELEM(gxx, float ,i,j); |
173 | for ( int i =1;i<h-1;i++) |
174 | for ( int j=1;j<w-1;j++){ |
175 | bound.at<uchar>(i,j) = 0.0; |
178 | double *f_bp = new double [h*w]; |
181 | for ( int i =1;i<h-1;i++) |
182 | for ( int j=1;j<w-1;j++){ |
184 | f_bp[idx] = -4*( int )bound.at<uchar>(i,j) |
185 | +( int )bound.at<uchar>(i,(j+1)) |
186 | +( int )bound.at<uchar>(i,(j-1)) |
187 | +( int )bound.at<uchar>((i-1),j) |
188 | +( int )bound.at<uchar>((i+1),j); |
191 | Mat diff = Mat(h,w,CV_32FC1); |
192 | for ( int i =0;i<h;i++){ |
193 | for ( int j=0;j<w;j++){ |
195 | diff.at< float >(i,j) = (CV_IMAGE_ELEM(lap, float ,i,j) - f_bp[idx]); |
199 | double *gtest = new double [(h-2)*(w-2)]; |
200 | for ( int i = 0 ; i < h-2;i++){ |
201 | for ( int j = 0 ; j < w-2; j++){ |
203 | gtest[idx] = diff.at< float >(i+1,j+1); |
209 | double *gfinal = new double [(h-2)*(w-2)]; |
210 | double *gfinal_t = new double [(h-2)*(w-2)]; |
211 | double *denom = new double [(h-2)*(w-2)]; |
212 | double *f3 = new double [(h-2)*(w-2)]; |
213 | double *f3_t = new double [(h-2)*(w-2)]; |
214 | double *img_d = new double [(h)*(w)]; |
216 | dst(gtest,gfinal,h-2,w-2); |
218 | transpose(gfinal,gfinal_t,h-2,w-2); |
220 | dst(gfinal_t,gfinal,w-2,h-2); |
222 | transpose(gfinal,gfinal_t,w-2,h-2); |
227 | for ( int i = 0 ; i < w-2;i++,cy++){ |
228 | for ( int j = 0,cx = 1; j < h-2; j++,cx++){ |
230 | denom[idx] = ( float ) 2* cos (pi*cy/( ( double ) (w-1))) - 2 + 2* cos (pi*cx/(( double ) (h-1))) - 2; |
235 | for (idx = 0 ; idx < (w-2)*(h-2) ;idx++){ |
236 | gfinal_t[idx] = gfinal_t[idx]/denom[idx]; |
239 | idst(gfinal_t,f3,h-2,w-2); |
241 | transpose(f3,f3_t,h-2,w-2); |
243 | idst(f3_t,f3,w-2,h-2); |
245 | transpose(f3,f3_t,w-2,h-2); |
247 | for ( int i = 0 ; i < h;i++){ |
248 | for ( int j = 0 ; j < w; j++){ |
250 | img_d[idx] = ( double )CV_IMAGE_ELEM(img,uchar,i,j); |
253 | for ( int i = 1 ; i < h-1;i++) |
255 | for ( int j = 1 ; j < w-1; j++) |
262 | for ( int i = 1,id1=0 ; i < h-1;i++,id1++){ |
263 | for ( int j = 1,id2=0 ; j < w-1; j++,id2++){ |
265 | idx1= id1*(w-2) + id2; |
266 | img_d[idx] = f3_t[idx1]; |
269 | for ( int i = 0 ; i < h;i++){ |
270 | for ( int j = 0 ; j < w; j++){ |
273 | result.at<uchar>(i,j) = 0; |
274 | else if (img_d[idx] > 255.0) |
275 | result.at<uchar>(i,j) = 255.0; |
277 | result.at<uchar>(i,j) = img_d[idx]; |
283 | IplImage *poisson_blend(IplImage *I, IplImage *mask, int posx, int posy){ |
285 | unsigned long int idx; |
286 | if (I->nChannels < 3){ |
287 | printf ( "Enter RGB image\n" ); |
290 | IplImage *grx = cvCreateImage(cvGetSize(I), 32, 3); |
291 | IplImage *gry = cvCreateImage(cvGetSize(I), 32, 3); |
293 | IplImage *sgx = cvCreateImage(cvGetSize(mask), 32, 3); |
294 | IplImage *sgy = cvCreateImage(cvGetSize(mask), 32, 3); |
296 | IplImage *S = cvCreateImage(cvGetSize(I), 8, 3); |
297 | IplImage *ero = cvCreateImage(cvGetSize(I), 8, 3); |
298 | IplImage *res = cvCreateImage(cvGetSize(I), 8, 3); |
303 | IplImage *O = cvCreateImage(cvGetSize(I), 8, 3); |
304 | IplImage *error= cvCreateImage(cvGetSize(I), 8, 3); |
308 | int channel = I->nChannels; |
310 | int wMask = mask->width; |
311 | int hMask = mask->height; |
312 | int channel1 = mask->nChannels; |
318 | getGradientx(mask,sgx); |
319 | getGradienty(mask,sgy); |
322 | for ( int i=posx, ii =0;i<posx+wMask;i++,ii++) |
323 | for ( int j=0,jj=posy;j<hMask;j++,jj++) |
324 | for ( int c=0;c<channel;++c){ |
325 | CV_IMAGE_ELEM(S,uchar,jj,i*channel+c) = 255; |
328 | IplImage* bmaskx = cvCreateImage(cvGetSize(ero),32,3); |
329 | cvConvertScale(S,bmaskx,1.0/255.0,0.0); |
331 | IplImage* bmasky = cvCreateImage(cvGetSize(ero),32,3); |
332 | cvConvertScale(S,bmasky,1.0/255.0,0.0); |
335 | for ( int i=posx, ii =0;i<posx+wMask;i++,ii++) |
336 | for ( int j=0,jj=posy;j<hMask;j++,jj++) |
337 | for ( int c=0;c<channel;++c){ |
338 | CV_IMAGE_ELEM(bmaskx, float ,jj,i*channel+c) = CV_IMAGE_ELEM(sgx, float ,j,ii*channel+c); |
339 | CV_IMAGE_ELEM(bmasky, float ,jj,i*channel+c) = CV_IMAGE_ELEM(sgy, float ,j,ii*channel+c); |
342 | cvErode(S,ero,NULL,1); |
344 | IplImage* smask = cvCreateImage(cvGetSize(ero),32,3); |
345 | cvConvertScale(ero,smask,1.0/255.0,0.0); |
347 | IplImage* srx32 = cvCreateImage(cvGetSize(res),32,3); |
348 | cvConvertScale(res,srx32,1.0/255.0,0.0); |
350 | IplImage* sry32 = cvCreateImage(cvGetSize(res),32,3); |
351 | cvConvertScale(res,sry32,1.0/255.0,0.0); |
353 | for ( int i=0;i < h; i++) |
354 | for ( int j=0; j < w; j++) |
355 | for ( int c=0;c<channel;++c){ |
356 | CV_IMAGE_ELEM(srx32, float ,i,j*channel+c) = |
357 | (CV_IMAGE_ELEM(bmaskx, float ,i,j*channel+c)*CV_IMAGE_ELEM(smask, float ,i,j*channel+c)); |
358 | CV_IMAGE_ELEM(sry32, float ,i,j*channel+c) = |
359 | (CV_IMAGE_ELEM(bmasky, float ,i,j*channel+c)*CV_IMAGE_ELEM(smask, float ,i,j*channel+c)); |
364 | IplImage* smask1 = cvCreateImage(cvGetSize(ero),32,3); |
365 | cvConvertScale(ero,smask1,1.0/255.0,0.0); |
367 | IplImage* grx32 = cvCreateImage(cvGetSize(res),32,3); |
368 | cvConvertScale(res,grx32,1.0/255.0,0.0); |
370 | IplImage* gry32 = cvCreateImage(cvGetSize(res),32,3); |
371 | cvConvertScale(res,gry32,1.0/255.0,0.0); |
373 | for ( int i=0;i < h; i++) |
374 | for ( int j=0; j < w; j++) |
375 | for ( int c=0;c<channel;++c){ |
376 | CV_IMAGE_ELEM(grx32, float ,i,j*channel+c) = |
377 | (CV_IMAGE_ELEM(grx, float ,i,j*channel+c)*CV_IMAGE_ELEM(smask1, float ,i,j*channel+c)); |
378 | CV_IMAGE_ELEM(gry32, float ,i,j*channel+c) = |
379 | (CV_IMAGE_ELEM(gry, float ,i,j*channel+c)*CV_IMAGE_ELEM(smask1, float ,i,j*channel+c)); |
382 | IplImage* fx = cvCreateImage(cvGetSize(res),32,3); |
383 | IplImage* fy = cvCreateImage(cvGetSize(res),32,3); |
385 | for ( int i=0;i < h; i++) |
386 | for ( int j=0; j < w; j++) |
387 | for ( int c=0;c<channel;++c){ |
388 | CV_IMAGE_ELEM(fx, float ,i,j*channel+c) = |
389 | (CV_IMAGE_ELEM(grx32, float ,i,j*channel+c)+CV_IMAGE_ELEM(srx32, float ,i,j*channel+c)); |
390 | CV_IMAGE_ELEM(fy, float ,i,j*channel+c) = |
391 | (CV_IMAGE_ELEM(gry32, float ,i,j*channel+c)+CV_IMAGE_ELEM(sry32, float ,i,j*channel+c)); |
394 | IplImage *gxx = cvCreateImage(cvGetSize(I), 32, 3); |
395 | IplImage *gyy = cvCreateImage(cvGetSize(I), 32, 3); |
400 | IplImage *rx_channel = cvCreateImage(cvGetSize(I), IPL_DEPTH_32F, 1 ); |
401 | IplImage *gx_channel = cvCreateImage(cvGetSize(I), IPL_DEPTH_32F, 1 ); |
402 | IplImage *bx_channel = cvCreateImage(cvGetSize(I), IPL_DEPTH_32F, 1 ); |
404 | cvCvtPixToPlane(gxx, rx_channel, gx_channel, bx_channel,0); |
406 | IplImage *ry_channel = cvCreateImage(cvGetSize(I), IPL_DEPTH_32F, 1 ); |
407 | IplImage *gy_channel = cvCreateImage(cvGetSize(I), IPL_DEPTH_32F, 1 ); |
408 | IplImage *by_channel = cvCreateImage(cvGetSize(I), IPL_DEPTH_32F, 1 ); |
410 | cvCvtPixToPlane(gyy, ry_channel, gy_channel, by_channel,0); |
412 | IplImage *r_channel = cvCreateImage(cvGetSize(I), 8, 1 ); |
413 | IplImage *g_channel = cvCreateImage(cvGetSize(I), 8, 1 ); |
414 | IplImage *b_channel = cvCreateImage(cvGetSize(I), 8, 1 ); |
416 | cvCvtPixToPlane(I, r_channel, g_channel, b_channel,0); |
418 | Mat resultr = Mat(h,w,CV_8UC1); |
419 | Mat resultg = Mat(h,w,CV_8UC1); |
420 | Mat resultb = Mat(h,w,CV_8UC1); |
421 | clock_t tic = clock (); |
423 | poisson_solver(r_channel,rx_channel, ry_channel,resultr); |
424 | poisson_solver(g_channel,gx_channel, gy_channel,resultg); |
425 | poisson_solver(b_channel,bx_channel, by_channel,resultb); |
427 | clock_t toc = clock (); |
428 | printf ( "Execution time: %f seconds\n" , ( double )(toc - tic) / CLOCKS_PER_SEC); |
430 | IplImage *final = cvCreateImage(cvGetSize(I), 8, 3 ); |
433 | for ( int j=0;j<w;j++){ |
434 | CV_IMAGE_ELEM(final,uchar,i,j*3+0) = resultr.at<uchar>(i,j); |
435 | CV_IMAGE_ELEM(final,uchar,i,j*3+1) = resultg.at<uchar>(i,j); |
436 | CV_IMAGE_ELEM(final,uchar,i,j*3+2) = resultb.at<uchar>(i,j); |
441 | Mat pbmethod(Mat img1, Mat img2, Rect ROI, int posX, int posY){ |
443 | int width1, width2, height1, height2; |
453 | IplImage abc = IplImage(img2); |
454 | IplImage subimg = IplImage(roiImg); |
456 | myresult = poisson_blend(&abc,&subimg,posX,posY); |
457 | Mat result(myresult); |
This code is slightly altered from Siddharth Kherada’s version
Additionally, a Matlab code (Hu Rui in University of Delaware) can be found HERE
NEW: I re-wrote the algorithm introduced above (not the FFT version), and you can find it HERE.