PageRenderTime 56ms CodeModel.GetById 16ms app.highlight 34ms RepoModel.GetById 1ms app.codeStats 0ms

/src/FreeImage/Source/FreeImageToolkit/MultigridPoissonSolver.cpp

https://bitbucket.org/cabalistic/ogredeps/
C++ | 505 lines | 297 code | 67 blank | 141 comment | 46 complexity | fb22ce7e555c960bcdc5780976223ec4 MD5 | raw file
  1// ==========================================================
  2// Poisson solver based on a full multigrid algorithm
  3//
  4// Design and implementation by
  5// - Hervé Drolon (drolon@infonie.fr)
  6// Reference:
  7// PRESS, W. H., TEUKOLSKY, S. A., VETTERLING, W. T., AND FLANNERY, B. P.
  8// 1992. Numerical Recipes in C: The Art of Scientific Computing, 2nd ed. Cambridge University Press.
  9//
 10// This file is part of FreeImage 3
 11//
 12// COVERED CODE IS PROVIDED UNDER THIS LICENSE ON AN "AS IS" BASIS, WITHOUT WARRANTY
 13// OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES
 14// THAT THE COVERED CODE IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE
 15// OR NON-INFRINGING. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE COVERED
 16// CODE IS WITH YOU. SHOULD ANY COVERED CODE PROVE DEFECTIVE IN ANY RESPECT, YOU (NOT
 17// THE INITIAL DEVELOPER OR ANY OTHER CONTRIBUTOR) ASSUME THE COST OF ANY NECESSARY
 18// SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL
 19// PART OF THIS LICENSE. NO USE OF ANY COVERED CODE IS AUTHORIZED HEREUNDER EXCEPT UNDER
 20// THIS DISCLAIMER.
 21//
 22// Use at your own risk!
 23// ==========================================================
 24
 25#include "FreeImage.h"
 26#include "Utilities.h"
 27#include "ToneMapping.h"
 28
 29static const int NPRE	= 1;		// Number of relaxation sweeps before ...
 30static const int NPOST	= 1;		// ... and after the coarse-grid correction is computed
 31static const int NGMAX	= 15;		// Maximum number of grids
 32
 33/**
 34Copy src into dst
 35*/
 36static inline void fmg_copyArray(FIBITMAP *dst, FIBITMAP *src) {
 37	memcpy(FreeImage_GetBits(dst), FreeImage_GetBits(src), FreeImage_GetHeight(dst) * FreeImage_GetPitch(dst));
 38}
 39
 40/**
 41Fills src with zeros
 42*/
 43static inline void fmg_fillArrayWithZeros(FIBITMAP *src) {
 44	memset(FreeImage_GetBits(src), 0, FreeImage_GetHeight(src) * FreeImage_GetPitch(src));
 45}
 46
 47/**
 48Half-weighting restriction. nc is the coarse-grid dimension. The fine-grid solution is input in
 49uf[0..2*nc-2][0..2*nc-2], the coarse-grid solution is returned in uc[0..nc-1][0..nc-1].
 50*/
 51static void fmg_restrict(FIBITMAP *UC, FIBITMAP *UF, int nc) {
 52	int row_uc, row_uf, col_uc, col_uf;
 53
 54	const int uc_pitch  = FreeImage_GetPitch(UC) / sizeof(float);
 55	const int uf_pitch  = FreeImage_GetPitch(UF) / sizeof(float);
 56	
 57	float *uc_bits = (float*)FreeImage_GetBits(UC);
 58	const float *uf_bits = (float*)FreeImage_GetBits(UF);
 59
 60	// interior points
 61	{
 62		float *uc_scan = uc_bits + uc_pitch;
 63		for (row_uc = 1, row_uf = 2; row_uc < nc-1; row_uc++, row_uf += 2) {
 64			const float *uf_scan = uf_bits + row_uf * uf_pitch;
 65			for (col_uc = 1, col_uf = 2; col_uc < nc-1; col_uc++, col_uf += 2) { 
 66				// calculate 
 67				// UC(row_uc, col_uc) = 
 68				// 0.5 * UF(row_uf, col_uf) + 0.125 * [ UF(row_uf+1, col_uf) + UF(row_uf-1, col_uf) + UF(row_uf, col_uf+1) + UF(row_uf, col_uf-1) ]
 69				float *uc_pixel = uc_scan + col_uc;
 70				const float *uf_center = uf_scan + col_uf;
 71				*uc_pixel = 0.5F * *uf_center + 0.125F * ( *(uf_center + uf_pitch) + *(uf_center - uf_pitch) + *(uf_center + 1) + *(uf_center - 1) );
 72			}
 73			uc_scan += uc_pitch;
 74		}
 75	}
 76	// boundary points
 77	const int ncc = 2*nc-1;
 78	{
 79		/*
 80		calculate the following: 
 81		for (row_uc = 0, row_uf = 0; row_uc < nc; row_uc++, row_uf += 2) { 
 82			UC(row_uc, 0) = UF(row_uf, 0);		
 83			UC(row_uc, nc-1) = UF(row_uf, ncc-1);
 84		}
 85		*/
 86		float *uc_scan = uc_bits;
 87		for (row_uc = 0, row_uf = 0; row_uc < nc; row_uc++, row_uf += 2) { 
 88			const float *uf_scan = uf_bits + row_uf * uf_pitch;
 89			uc_scan[0] = uf_scan[0];
 90			uc_scan[nc-1] = uf_scan[ncc-1];
 91			uc_scan += uc_pitch;
 92		}
 93	}
 94	{
 95		/*
 96		calculate the following: 
 97		for (col_uc = 0, col_uf = 0; col_uc < nc; col_uc++, col_uf += 2) {
 98			UC(0, col_uc) = UF(0, col_uf);
 99			UC(nc-1, col_uc) = UF(ncc-1, col_uf);
100		}
101		*/
102		float *uc_scan_top = uc_bits;
103		float *uc_scan_bottom = uc_bits + (nc-1)*uc_pitch;
104		const float *uf_scan_top = uf_bits + (ncc-1)*uf_pitch;
105		const float *uf_scan_bottom = uf_bits;
106		for (col_uc = 0, col_uf = 0; col_uc < nc; col_uc++, col_uf += 2) {
107			uc_scan_top[col_uc] = uf_scan_top[col_uf];
108			uc_scan_bottom[col_uc] = uf_scan_bottom[col_uf];
109		}
110	}
111}
112
113/**
114Solution of the model problem on the coarsest grid, where h = 1/2 . 
115The right-hand side is input
116in rhs[0..2][0..2] and the solution is returned in u[0..2][0..2].
117*/
118static void fmg_solve(FIBITMAP *U, FIBITMAP *RHS) {
119	// fill U with zeros
120	fmg_fillArrayWithZeros(U);
121	// calculate U(1, 1) = -h*h*RHS(1, 1)/4.0 where h = 1/2
122	float *u_scan = (float*)FreeImage_GetScanLine(U, 1);
123	const float *rhs_scan = (float*)FreeImage_GetScanLine(RHS, 1);
124	u_scan[1] = -rhs_scan[1] / 16;
125}
126
127/**
128Coarse-to-fine prolongation by bilinear interpolation. nf is the fine-grid dimension. The coarsegrid
129solution is input as uc[0..nc-1][0..nc-1], where nc = nf/2 + 1. The fine-grid solution is
130returned in uf[0..nf-1][0..nf-1].
131*/
132static void fmg_prolongate(FIBITMAP *UF, FIBITMAP *UC, int nf) {
133	int row_uc, row_uf, col_uc, col_uf;
134
135	const int uf_pitch  = FreeImage_GetPitch(UF) / sizeof(float);
136	const int uc_pitch  = FreeImage_GetPitch(UC) / sizeof(float);
137	
138	float *uf_bits = (float*)FreeImage_GetBits(UF);
139	const float *uc_bits = (float*)FreeImage_GetBits(UC);
140	
141	// do elements that are copies
142	{
143		const int nc = nf/2 + 1;
144
145		float *uf_scan = uf_bits;
146		const float *uc_scan = uc_bits;		
147		for (row_uc = 0; row_uc < nc; row_uc++) {
148			for (col_uc = 0, col_uf = 0; col_uc < nc; col_uc++, col_uf += 2) {
149				// calculate UF(2*row_uc, col_uf) = UC(row_uc, col_uc);
150				uf_scan[col_uf] = uc_scan[col_uc];
151			}
152			uc_scan += uc_pitch;
153			uf_scan += 2 * uf_pitch;
154		}
155	}
156	// do odd-numbered columns, interpolating vertically
157	{		
158		for(row_uf = 1; row_uf < nf-1; row_uf += 2) {
159			float *uf_scan = uf_bits + row_uf * uf_pitch;
160			for (col_uf = 0; col_uf < nf; col_uf += 2) {
161				// calculate UF(row_uf, col_uf) = 0.5 * ( UF(row_uf+1, col_uf) + UF(row_uf-1, col_uf) )
162				uf_scan[col_uf] = 0.5F * ( *(uf_scan + uf_pitch + col_uf) + *(uf_scan - uf_pitch + col_uf) );
163			}
164		}
165	}
166	// do even-numbered columns, interpolating horizontally
167	{
168		float *uf_scan = uf_bits;
169		for(row_uf = 0; row_uf < nf; row_uf++) {
170			for (col_uf = 1; col_uf < nf-1; col_uf += 2) {
171				// calculate UF(row_uf, col_uf) = 0.5 * ( UF(row_uf, col_uf+1) + UF(row_uf, col_uf-1) )
172				uf_scan[col_uf] = 0.5F * ( uf_scan[col_uf + 1] + uf_scan[col_uf - 1] );
173			}
174			uf_scan += uf_pitch;
175		}
176	}
177}
178
179/**
180Red-black Gauss-Seidel relaxation for model problem. Updates the current value of the solution
181u[0..n-1][0..n-1], using the right-hand side function rhs[0..n-1][0..n-1].
182*/
183static void fmg_relaxation(FIBITMAP *U, FIBITMAP *RHS, int n) {
184	int row, col, ipass, isw, jsw;
185	const float h = 1.0F / (n - 1);
186	const float h2 = h*h;
187
188	const int u_pitch  = FreeImage_GetPitch(U) / sizeof(float);
189	const int rhs_pitch  = FreeImage_GetPitch(RHS) / sizeof(float);
190	
191	float *u_bits = (float*)FreeImage_GetBits(U);
192	const float *rhs_bits = (float*)FreeImage_GetBits(RHS);
193
194	for (ipass = 0, jsw = 1; ipass < 2; ipass++, jsw = 3-jsw) { // Red and black sweeps
195		float *u_scan = u_bits + u_pitch;
196		const float *rhs_scan = rhs_bits + rhs_pitch;
197		for (row = 1, isw = jsw; row < n-1; row++, isw = 3-isw) {
198			for (col = isw; col < n-1; col += 2) { 
199				// Gauss-Seidel formula
200				// calculate U(row, col) = 
201				// 0.25 * [ U(row+1, col) + U(row-1, col) + U(row, col+1) + U(row, col-1) - h2 * RHS(row, col) ]		 
202				float *u_center = u_scan + col;
203				const float *rhs_center = rhs_scan + col;
204				*u_center = *(u_center + u_pitch) + *(u_center - u_pitch) + *(u_center + 1) + *(u_center - 1);
205				*u_center -= h2 * *rhs_center;
206				*u_center *= 0.25F;
207			}
208			u_scan += u_pitch;
209			rhs_scan += rhs_pitch;
210		}
211	}
212}
213
214/**
215Returns minus the residual for the model problem. Input quantities are u[0..n-1][0..n-1] and
216rhs[0..n-1][0..n-1], while res[0..n-1][0..n-1] is returned.
217*/
218static void fmg_residual(FIBITMAP *RES, FIBITMAP *U, FIBITMAP *RHS, int n) {
219	int row, col;
220
221	const float h = 1.0F / (n-1);	
222	const float h2i = 1.0F / (h*h);
223
224	const int res_pitch  = FreeImage_GetPitch(RES) / sizeof(float);
225	const int u_pitch  = FreeImage_GetPitch(U) / sizeof(float);
226	const int rhs_pitch  = FreeImage_GetPitch(RHS) / sizeof(float);
227	
228	float *res_bits = (float*)FreeImage_GetBits(RES);
229	const float *u_bits = (float*)FreeImage_GetBits(U);
230	const float *rhs_bits = (float*)FreeImage_GetBits(RHS);
231
232	// interior points
233	{
234		float *res_scan = res_bits + res_pitch;
235		const float *u_scan = u_bits + u_pitch;
236		const float *rhs_scan = rhs_bits + rhs_pitch;
237		for (row = 1; row < n-1; row++) {
238			for (col = 1; col < n-1; col++) {
239				// calculate RES(row, col) = 
240				// -h2i * [ U(row+1, col) + U(row-1, col) + U(row, col+1) + U(row, col-1) - 4 * U(row, col) ] + RHS(row, col);
241				float *res_center = res_scan + col;
242				const float *u_center = u_scan + col;
243				const float *rhs_center = rhs_scan + col;
244				*res_center = *(u_center + u_pitch) + *(u_center - u_pitch) + *(u_center + 1) + *(u_center - 1) - 4 * *u_center;
245				*res_center *= -h2i;
246				*res_center += *rhs_center;
247			}
248			res_scan += res_pitch;
249			u_scan += u_pitch;
250			rhs_scan += rhs_pitch;
251		}
252	}
253
254	// boundary points
255	{
256		memset(FreeImage_GetScanLine(RES, 0), 0, FreeImage_GetPitch(RES));
257		memset(FreeImage_GetScanLine(RES, n-1), 0, FreeImage_GetPitch(RES));
258		float *left = res_bits;
259		float *right = res_bits + (n-1);
260		for(int k = 0; k < n; k++) {
261			*left = 0;
262			*right = 0;
263			left += res_pitch;
264			right += res_pitch;
265		}
266	}
267}
268
269/**
270Does coarse-to-fine interpolation and adds result to uf. nf is the fine-grid dimension. The
271coarse-grid solution is input as uc[0..nc-1][0..nc-1], where nc = nf/2+1. The fine-grid solution
272is returned in uf[0..nf-1][0..nf-1]. res[0..nf-1][0..nf-1] is used for temporary storage.
273*/
274static void fmg_addint(FIBITMAP *UF, FIBITMAP *UC, FIBITMAP *RES, int nf) {
275	fmg_prolongate(RES, UC, nf);
276
277	const int uf_pitch  = FreeImage_GetPitch(UF) / sizeof(float);
278	const int res_pitch  = FreeImage_GetPitch(RES) / sizeof(float);	
279
280	float *uf_bits = (float*)FreeImage_GetBits(UF);
281	const float *res_bits = (float*)FreeImage_GetBits(RES);
282
283	for(int row = 0; row < nf; row++) {
284		for(int col = 0; col < nf; col++) {
285			// calculate UF(row, col) = UF(row, col) + RES(row, col);
286			uf_bits[col] += res_bits[col];
287		}
288		uf_bits += uf_pitch;
289		res_bits += res_pitch;
290	}
291}
292
293/**
294Full Multigrid Algorithm for solution of linear elliptic equation, here the model problem (19.0.6).
295On input u[0..n-1][0..n-1] contains the right-hand side ñ, while on output it returns the solution.
296The dimension n must be of the form 2^j + 1 for some integer j. (j is actually the number of
297grid levels used in the solution, called ng below.) ncycle is the number of V-cycles to be
298used at each level.
299*/
300static BOOL fmg_mglin(FIBITMAP *U, int n, int ncycle) {
301	int j, jcycle, jj, jpost, jpre, nf, ngrid;
302
303	FIBITMAP **IRHO = NULL;
304	FIBITMAP **IU   = NULL;
305	FIBITMAP **IRHS = NULL;
306	FIBITMAP **IRES = NULL;
307	
308	int ng = 0;		// number of allocated grids
309
310// --------------------------------------------------------------------------
311
312#define _CREATE_ARRAY_GRID_(array, array_size) \
313	array = (FIBITMAP**)malloc(array_size * sizeof(FIBITMAP*));\
314	if(!array) throw(1);\
315	memset(array, 0, array_size * sizeof(FIBITMAP*))
316
317#define _FREE_ARRAY_GRID_(array, array_size) \
318	if(NULL != array) {\
319		for(int k = 0; k < array_size; k++) {\
320			if(NULL != array[k]) {\
321				FreeImage_Unload(array[k]); array[k] = NULL;\
322			}\
323		}\
324		free(array);\
325	}
326
327// --------------------------------------------------------------------------
328
329	try {
330		int nn = n;
331		// check grid size and grid levels
332		while (nn >>= 1) ng++;
333		if (n != 1 + (1L << ng)) {
334			FreeImage_OutputMessageProc(FIF_UNKNOWN, "Multigrid algorithm: n = %d, while n-1 must be a power of 2.", n);
335			throw(1);
336		}
337		if (ng > NGMAX) {
338			FreeImage_OutputMessageProc(FIF_UNKNOWN, "Multigrid algorithm: ng = %d while NGMAX = %d, increase NGMAX.", ng, NGMAX);
339			throw(1);
340		}
341		// allocate grid arrays
342		{
343			_CREATE_ARRAY_GRID_(IRHO, ng);
344			_CREATE_ARRAY_GRID_(IU, ng);
345			_CREATE_ARRAY_GRID_(IRHS, ng);
346			_CREATE_ARRAY_GRID_(IRES, ng);
347		}
348
349		nn = n/2 + 1;
350		ngrid = ng - 2;
351
352		// allocate storage for r.h.s. on grid (ng - 2) ...
353		IRHO[ngrid] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
354		if(!IRHO[ngrid]) throw(1);
355
356		// ... and fill it by restricting from the fine grid
357		fmg_restrict(IRHO[ngrid], U, nn);	
358
359		// similarly allocate storage and fill r.h.s. on all coarse grids.
360		while (nn > 3) {
361			nn = nn/2 + 1; 
362			ngrid--;
363			IRHO[ngrid] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
364			if(!IRHO[ngrid]) throw(1);
365			fmg_restrict(IRHO[ngrid], IRHO[ngrid+1], nn);
366		}
367
368		nn = 3;
369
370		IU[0] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
371		if(!IU[0]) throw(1);
372		IRHS[0] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
373		if(!IRHS[0]) throw(1);
374
375		// initial solution on coarsest grid
376		fmg_solve(IU[0], IRHO[0]);
377		// irho[0] no longer needed ...
378		FreeImage_Unload(IRHO[0]); IRHO[0] = NULL;
379
380		ngrid = ng;
381
382		// nested iteration loop
383		for (j = 1; j < ngrid; j++) {
384			nn = 2*nn - 1;
385
386			IU[j] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
387			if(!IU[j]) throw(1);
388			IRHS[j] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
389			if(!IRHS[j]) throw(1);
390			IRES[j] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
391			if(!IRES[j]) throw(1);
392
393			fmg_prolongate(IU[j], IU[j-1], nn);
394			
395			// interpolate from coarse grid to next finer grid
396
397			// set up r.h.s.
398			fmg_copyArray(IRHS[j], j != (ngrid - 1) ? IRHO[j] : U);
399			
400			// V-cycle loop
401			for (jcycle = 0; jcycle < ncycle; jcycle++) {
402				nf = nn;
403				// downward stoke of the V
404				for (jj = j; jj >= 1; jj--) {
405					// pre-smoothing
406					for (jpre = 0; jpre < NPRE; jpre++) {
407						fmg_relaxation(IU[jj], IRHS[jj], nf);
408					}
409					fmg_residual(IRES[jj], IU[jj], IRHS[jj], nf);
410					nf = nf/2 + 1;
411					// restriction of the residual is the next r.h.s.
412					fmg_restrict(IRHS[jj-1], IRES[jj], nf);				
413					// zero for initial guess in next relaxation
414					fmg_fillArrayWithZeros(IU[jj-1]);
415				}
416				// bottom of V: solve on coarsest grid
417				fmg_solve(IU[0], IRHS[0]); 
418				nf = 3; 
419				// upward stroke of V.
420				for (jj = 1; jj <= j; jj++) { 
421					nf = 2*nf - 1;
422					// use res for temporary storage inside addint
423					fmg_addint(IU[jj], IU[jj-1], IRES[jj], nf);				
424					// post-smoothing
425					for (jpost = 0; jpost < NPOST; jpost++) {
426						fmg_relaxation(IU[jj], IRHS[jj], nf);
427					}
428				}
429			}
430		}
431
432		// return solution in U
433		fmg_copyArray(U, IU[ngrid-1]);
434
435		// delete allocated arrays
436		_FREE_ARRAY_GRID_(IRES, ng);
437		_FREE_ARRAY_GRID_(IRHS, ng);
438		_FREE_ARRAY_GRID_(IU, ng);
439		_FREE_ARRAY_GRID_(IRHO, ng);
440
441		return TRUE;
442
443	} catch(int) {
444		// delete allocated arrays
445		_FREE_ARRAY_GRID_(IRES, ng);
446		_FREE_ARRAY_GRID_(IRHS, ng);
447		_FREE_ARRAY_GRID_(IU, ng);
448		_FREE_ARRAY_GRID_(IRHO, ng);
449
450		return FALSE;
451	}
452}
453
454// --------------------------------------------------------------------------
455
456/**
457Poisson solver based on a multigrid algorithm. 
458This routine solves a Poisson equation, remap result pixels to [0..1] and returns the solution. 
459NB: The input image is first stored inside a square image whose size is (2^j + 1)x(2^j + 1) for some integer j, 
460where j is such that 2^j is the nearest larger dimension corresponding to MAX(image width, image height). 
461@param Laplacian Laplacian image
462@param ncycle Number of cycles in the multigrid algorithm (usually 2 or 3)
463@return Returns the solved PDE equations if successful, returns NULL otherwise
464*/
465FIBITMAP* DLL_CALLCONV 
466FreeImage_MultigridPoissonSolver(FIBITMAP *Laplacian, int ncycle) {
467	if(!FreeImage_HasPixels(Laplacian)) return NULL;
468
469	int width = FreeImage_GetWidth(Laplacian);
470	int height = FreeImage_GetHeight(Laplacian);
471
472	// get nearest larger dimension length that is acceptable by the algorithm
473	int n = MAX(width, height);
474	int size = 0;
475	while((n >>= 1) > 0) size++;
476	if((1 << size) < MAX(width, height)) {
477		size++;
478	}
479	// size must be of the form 2^j + 1 for some integer j
480	size = 1 + (1 << size);
481
482	// allocate a temporary square image I
483	FIBITMAP *I = FreeImage_AllocateT(FIT_FLOAT, size, size);
484	if(!I) return NULL;
485
486	// copy Laplacian into I and shift pixels to create a boundary
487	FreeImage_Paste(I, Laplacian, 1, 1, 255);
488
489	// solve the PDE equation
490	fmg_mglin(I, size, ncycle);
491
492	// shift pixels back
493	FIBITMAP *U = FreeImage_Copy(I, 1, 1, width + 1, height + 1);
494	FreeImage_Unload(I);
495
496	// remap pixels to [0..1]
497	NormalizeY(U, 0, 1);
498
499	// copy metadata from src to dst
500	FreeImage_CloneMetadata(U, Laplacian);
501
502	// return the integrated image
503	return U;
504}
505