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
29 static const int NPRE = 1; // Number of relaxation sweeps before ...
30 static const int NPOST = 1; // ... and after the coarse-grid correction is computed
31 static const int NGMAX = 15; // Maximum number of grids
32
33 /**
34 Copy src into dst
35 */
fmg_copyArray(FIBITMAP * dst,FIBITMAP * src)36 static 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 /**
41 Fills src with zeros
42 */
fmg_fillArrayWithZeros(FIBITMAP * src)43 static inline void fmg_fillArrayWithZeros(FIBITMAP *src) {
44 memset(FreeImage_GetBits(src), 0, FreeImage_GetHeight(src) * FreeImage_GetPitch(src));
45 }
46
47 /**
48 Half-weighting restriction. nc is the coarse-grid dimension. The fine-grid solution is input in
49 uf[0..2*nc-2][0..2*nc-2], the coarse-grid solution is returned in uc[0..nc-1][0..nc-1].
50 */
fmg_restrict(FIBITMAP * UC,FIBITMAP * UF,int nc)51 static 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 /**
114 Solution of the model problem on the coarsest grid, where h = 1/2 .
115 The right-hand side is input
116 in rhs[0..2][0..2] and the solution is returned in u[0..2][0..2].
117 */
fmg_solve(FIBITMAP * U,FIBITMAP * RHS)118 static 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 /**
128 Coarse-to-fine prolongation by bilinear interpolation. nf is the fine-grid dimension. The coarsegrid
129 solution is input as uc[0..nc-1][0..nc-1], where nc = nf/2 + 1. The fine-grid solution is
130 returned in uf[0..nf-1][0..nf-1].
131 */
fmg_prolongate(FIBITMAP * UF,FIBITMAP * UC,int nf)132 static 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 /**
180 Red-black Gauss-Seidel relaxation for model problem. Updates the current value of the solution
181 u[0..n-1][0..n-1], using the right-hand side function rhs[0..n-1][0..n-1].
182 */
fmg_relaxation(FIBITMAP * U,FIBITMAP * RHS,int n)183 static 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 /**
215 Returns minus the residual for the model problem. Input quantities are u[0..n-1][0..n-1] and
216 rhs[0..n-1][0..n-1], while res[0..n-1][0..n-1] is returned.
217 */
fmg_residual(FIBITMAP * RES,FIBITMAP * U,FIBITMAP * RHS,int n)218 static 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 /**
270 Does coarse-to-fine interpolation and adds result to uf. nf is the fine-grid dimension. The
271 coarse-grid solution is input as uc[0..nc-1][0..nc-1], where nc = nf/2+1. The fine-grid solution
272 is returned in uf[0..nf-1][0..nf-1]. res[0..nf-1][0..nf-1] is used for temporary storage.
273 */
fmg_addint(FIBITMAP * UF,FIBITMAP * UC,FIBITMAP * RES,int nf)274 static 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 /**
294 Full Multigrid Algorithm for solution of linear elliptic equation, here the model problem (19.0.6).
295 On input u[0..n-1][0..n-1] contains the right-hand side �, while on output it returns the solution.
296 The dimension n must be of the form 2^j + 1 for some integer j. (j is actually the number of
297 grid levels used in the solution, called ng below.) ncycle is the number of V-cycles to be
298 used at each level.
299 */
fmg_mglin(FIBITMAP * U,int n,int ncycle)300 static 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 /**
457 Poisson solver based on a multigrid algorithm.
458 This routine solves a Poisson equation, remap result pixels to [0..1] and returns the solution.
459 NB: The input image is first stored inside a square image whose size is (2^j + 1)x(2^j + 1) for some integer j,
460 where 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 */
465 FIBITMAP* DLL_CALLCONV
FreeImage_MultigridPoissonSolver(FIBITMAP * Laplacian,int ncycle)466 FreeImage_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
506