Poisson Blending

Poisson Blending

from: http://eric-yuan.me/poisson-blending/

By Eric | Published: October 15, 2013

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:

01	for i=1:row number

02	for j=1:col number

03

if(i==j)

04	matrix(i, j)=-4

05	elif(adjacent(pixel(j), pixel(i)))

06	matrix(i, j)=1

07

else

08	matrix(i, j)=0

09

end

10

end

11

end

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

001	#include <iostream>

002	#include <complex>

003	#include "cv.h"

004	#include "cxcore.h"

005	#include "highgui.h"

006	#include "math.h"

007

008	using namespace std;

009	using namespace cv;

010

011	#define pi 3.1416

012

013	void display(const char *name, IplImage *img){

014	cvNamedWindow(name);

015	cvShowImage(name,img);

016

}

017

018	void getGradientx( const IplImage *img, IplImage *gx){

019	int w = img->width;

020	int h = img->height;

021	int channel = img->nChannels;

022

023	cvZero( gx );

024	for(int i=0;i<h;i++)

025	for(int j=0;j<w;j++)

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);

029

}

030

}

031	void getGradienty( const IplImage *img, IplImage *gy){

032	int w = img->width;

033	int h = img->height;

034	int channel = img->nChannels;

035

036	cvZero( gy );

037	for(int i=0;i<h;i++)

038	for(int j=0;j<w;j++)

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);

042

043

}

044

}

045	void lapx( const IplImage *img, IplImage *gxx){

046	int w = img->width;

047	int h = img->height;

048	int channel = img->nChannels;

049

050	cvZero( gxx );

051	for(int i=0;i<h;i++)

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);

056

}

057

}

058	void lapy( const IplImage *img, IplImage *gyy){

059	int w = img->width;

060	int h = img->height;

061	int channel = img->nChannels;

062	cvZero( gyy );

063	for(int i=0;i<h-1;i++)

064	for(int j=0;j<w;j++)

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);

068

069

}

070

}

071

072	void dst(double *gtest, double *gfinal,int h,int w){

073

074	int k,r,z;

075	unsigned long int idx;

076

077	Mat temp = Mat(2*h+2,1,CV_32F);

078	Mat res  = Mat(h,1,CV_32F);

079

080

int p=0;

081	for(int i=0;i<w;i++){

082	temp.at<float>(0,0) = 0.0;

083

084	for(int j=0,r=1;j<h;j++,r++){

085	idx = j*w+i;

086	temp.at<float>(r,0) = gtest[idx];

087

}

088

089	temp.at<float>(h+1,0)=0.0;

090

091	for(int j=h-1, r=h+2;j>=0;j--,r++){

092	idx = j*w+i;

093	temp.at<float>(r,0) = -1*gtest[idx];

094

}

095

096	Mat planes[] = {Mat_<float>(temp), Mat::zeros(temp.size(), CV_32F)};

097

098	Mat complex1;

099	merge(planes, 2, complex1);

100

101	dft(complex1,complex1,0,0);

102

103	Mat planes1[] = {Mat::zeros(complex1.size(), CV_32F), Mat::zeros(complex1.size(), CV_32F)};

104

105	// planes1[0] = Re(DFT(I)), planes1[1] = Im(DFT(I))

106	split(complex1, planes1);

107

108	std::complex<double> two_i = std::sqrt(std::complex<double>(-1));

109

110	double fac = -2*imag(two_i);

111

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;

114

}

115

116	for(int q=0,z=0;q<h;q++,z++){

117	idx = q*w+p;

118	gfinal[idx] =  res.at<float>(z,0);

119

}

120

p++;

121

}

122

123

}

124

125	void idst(double *gtest, double *gfinal,int h,int w){

126	int nn = h+1;

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++){

131	idx = i*w + j;

132	gfinal[idx] = (double) (2*gfinal[idx])/nn;

133

}

134

135

}

136	void transpose(double *mat, double *mat_t,int h,int w){

137

138	Mat tmp = Mat(h,w,CV_32FC1);

139

int p =0;

140	unsigned long int idx;

141	for(int i = 0 ; i < h;i++){

142	for(int j = 0 ; j < w; j++){

143

144	idx = i*(w) + j;

145	tmp.at<float>(i,j) = mat[idx];

146

}

147

}

148	Mat tmp_t = tmp.t();

149

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);

154

}

155

156

}

157	void poisson_solver(const IplImage *img, IplImage *gxx , IplImage *gyy, Mat &result){

158

159	int w = img->width;

160	int h = img->height;

161	int channel = img->nChannels;

162

163	unsigned long int idx,idx1;

164

165	IplImage *lap  = cvCreateImage(cvGetSize(img), 32, 1);

166

167	for(int i =0;i<h;i++)

168	for(int j=0;j<w;j++)

169	CV_IMAGE_ELEM(lap,float,i,j)=CV_IMAGE_ELEM(gyy,float,i,j)+CV_IMAGE_ELEM(gxx,float,i,j);

170

171	Mat bound(img);

172

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;

176

}

177

178	double *f_bp = new double[h*w];

179

180

//div(G)

181	for(int i =1;i<h-1;i++)

182	for(int j=1;j<w-1;j++){

183	idx=i*w + 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);

189

}

190

191	Mat diff = Mat(h,w,CV_32FC1);

192	for(int i =0;i<h;i++){

193	for(int j=0;j<w;j++){

194	idx = i*w+j;

195	diff.at<float>(i,j) = (CV_IMAGE_ELEM(lap,float,i,j) - f_bp[idx]);

196

}

197

}

198

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++){

202	idx = i*(w-2) + j;

203	gtest[idx] = diff.at<float>(i+1,j+1);

204

205

}

206

}

207	/ Find DST  /

208

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)];

215

216	dst(gtest,gfinal,h-2,w-2);

217

218	transpose(gfinal,gfinal_t,h-2,w-2);

219

220	dst(gfinal_t,gfinal,w-2,h-2);

221

222	transpose(gfinal,gfinal_t,w-2,h-2);

223

224

int cx=1;

225

int cy=1;

226

227	for(int i = 0 ; i < w-2;i++,cy++){

228	for(int j = 0,cx = 1; j < h-2; j++,cx++){

229	idx = j*(w-2) + i;

230	denom[idx] = (float) 2*cos(pi*cy/( (double) (w-1))) - 2 + 2*cos(pi*cx/((double) (h-1))) - 2;

231

232

}

233

}

234

235	for(idx = 0 ; idx < (w-2)*(h-2) ;idx++){

236	gfinal_t[idx] = gfinal_t[idx]/denom[idx];

237

}

238

239	idst(gfinal_t,f3,h-2,w-2);

240

241	transpose(f3,f3_t,h-2,w-2);

242

243	idst(f3_t,f3,w-2,h-2);

244

245	transpose(f3,f3_t,w-2,h-2);

246

247	for(int i = 0 ; i < h;i++){

248	for(int j = 0 ; j < w; j++){

249	idx = i*w + j;

250	img_d[idx] = (double)CV_IMAGE_ELEM(img,uchar,i,j);

251

}

252

}

253	for(int i = 1 ; i < h-1;i++)

254

{

255	for(int j = 1 ; j < w-1; j++)

256

{

257	idx = i*w + j;

258	img_d[idx] = 0.0;

259

}

260

}

261	int id1,id2;

262	for(int i = 1,id1=0 ; i < h-1;i++,id1++){

263	for(int j = 1,id2=0 ; j < w-1; j++,id2++){

264	idx = i*w + j;

265	idx1= id1*(w-2) + id2;

266	img_d[idx] = f3_t[idx1];

267

}

268

}

269	for(int i = 0 ; i < h;i++){

270	for(int j = 0 ; j < w; j++){

271	idx = i*w + j;

272	if(img_d[idx] < 0.0)

273	result.at<uchar>(i,j) = 0;

274	else if(img_d[idx] > 255.0)

275	result.at<uchar>(i,j) = 255.0;

276

else

277	result.at<uchar>(i,j) = img_d[idx];

278

}

279

}

280

281

}

282

283	IplImage *poisson_blend(IplImage *I, IplImage *mask, int posx, int posy){

284

285	unsigned long int idx;

286	if(I->nChannels < 3){

287	printf("Enter RGB image\n");

288

exit(0);

289

}

290	IplImage *grx  = cvCreateImage(cvGetSize(I), 32, 3);

291	IplImage *gry  = cvCreateImage(cvGetSize(I), 32, 3);

292

293	IplImage *sgx  = cvCreateImage(cvGetSize(mask), 32, 3);

294	IplImage *sgy  = cvCreateImage(cvGetSize(mask), 32, 3);

295

296	IplImage *S    = cvCreateImage(cvGetSize(I), 8, 3);

297	IplImage *ero  = cvCreateImage(cvGetSize(I), 8, 3);

298	IplImage *res  = cvCreateImage(cvGetSize(I), 8, 3);

299

300	cvZero(S);

301	cvZero(res);

302

303	IplImage *O    = cvCreateImage(cvGetSize(I), 8, 3);

304	IplImage *error= cvCreateImage(cvGetSize(I), 8, 3);

305

306	int w = I->width;

307	int h = I->height;

308	int channel = I->nChannels;

309

310	int wMask = mask->width;

311	int hMask = mask->height;

312	int channel1 = mask->nChannels;

313

314	//get gradiant map of both I and mask

315	getGradientx(I,grx);

316	getGradienty(I,gry);

317

318	getGradientx(mask,sgx);

319	getGradienty(mask,sgy);

320

321	//create a mask map which white pixels is where mask in I

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;

326

}

327

328	IplImage* bmaskx = cvCreateImage(cvGetSize(ero),32,3);

329	cvConvertScale(S,bmaskx,1.0/255.0,0.0);

330

331	IplImage* bmasky = cvCreateImage(cvGetSize(ero),32,3);

332	cvConvertScale(S,bmasky,1.0/255.0,0.0);

333

334	//copy gradiant map of mask into mask map S

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);

340

}

341

342	cvErode(S,ero,NULL,1);

343

344	IplImage* smask = cvCreateImage(cvGetSize(ero),32,3);

345	cvConvertScale(ero,smask,1.0/255.0,0.0);

346

347	IplImage* srx32 = cvCreateImage(cvGetSize(res),32,3);

348	cvConvertScale(res,srx32,1.0/255.0,0.0);

349

350	IplImage* sry32 = cvCreateImage(cvGetSize(res),32,3);

351	cvConvertScale(res,sry32,1.0/255.0,0.0);

352

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));

360

}

361

362	cvNot(ero,ero);

363

364	IplImage* smask1 = cvCreateImage(cvGetSize(ero),32,3);

365	cvConvertScale(ero,smask1,1.0/255.0,0.0);

366

367	IplImage* grx32 = cvCreateImage(cvGetSize(res),32,3);

368	cvConvertScale(res,grx32,1.0/255.0,0.0);

369

370	IplImage* gry32 = cvCreateImage(cvGetSize(res),32,3);

371	cvConvertScale(res,gry32,1.0/255.0,0.0);

372

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));

380

}

381

382	IplImage* fx = cvCreateImage(cvGetSize(res),32,3);

383	IplImage* fy = cvCreateImage(cvGetSize(res),32,3);

384	//create a new gradient map, with I's gradient and mask's gradient in ROI

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));

392

}

393

394	IplImage *gxx  = cvCreateImage(cvGetSize(I), 32, 3);

395	IplImage *gyy  = cvCreateImage(cvGetSize(I), 32, 3);

396

397	lapx(fx,gxx);

398	lapy(fy,gyy);

399

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 );

403

404	cvCvtPixToPlane(gxx, rx_channel, gx_channel, bx_channel,0);

405

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 );

409

410	cvCvtPixToPlane(gyy, ry_channel, gy_channel, by_channel,0);

411

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 );

415

416	cvCvtPixToPlane(I, r_channel, g_channel, b_channel,0);

417

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();

422

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);

426

427	clock_t toc = clock();

428	printf("Execution time: %f seconds\n", (double)(toc - tic) / CLOCKS_PER_SEC);

429

430	IplImage *final = cvCreateImage(cvGetSize(I), 8, 3 );

431

432	for(int i=0;i<h;i++)

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);

437

}

438	return final;

439

}

440

441	Mat pbmethod(Mat img1, Mat img2, Rect ROI, int posX, int posY){

442

443	int width1, width2, height1, height2;

444	width1 = img1.cols;

445	width2 = img2.cols;

446	height1 = img1.rows;

447	height2 = img2.rows;

448

449	Mat roiImg;

450	roiImg = img1(ROI);

451

452	IplImage* myresult;

453	IplImage abc = IplImage(img2);

454	IplImage subimg = IplImage(roiImg);

455

456	myresult = poisson_blend(&abc,&subimg,posX,posY);

457	Mat result(myresult);

458	return result;

459

}

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.