summaryrefslogtreecommitdiff
path: root/plugins/FreeImage/Source/FreeImageToolkit/MultigridPoissonSolver.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'plugins/FreeImage/Source/FreeImageToolkit/MultigridPoissonSolver.cpp')
-rw-r--r--plugins/FreeImage/Source/FreeImageToolkit/MultigridPoissonSolver.cpp1010
1 files changed, 505 insertions, 505 deletions
diff --git a/plugins/FreeImage/Source/FreeImageToolkit/MultigridPoissonSolver.cpp b/plugins/FreeImage/Source/FreeImageToolkit/MultigridPoissonSolver.cpp
index 03281a16b5..3b577cbd0a 100644
--- a/plugins/FreeImage/Source/FreeImageToolkit/MultigridPoissonSolver.cpp
+++ b/plugins/FreeImage/Source/FreeImageToolkit/MultigridPoissonSolver.cpp
@@ -1,505 +1,505 @@
-// ==========================================================
-// Poisson solver based on a full multigrid algorithm
-//
-// Design and implementation by
-// - Hervé Drolon (drolon@infonie.fr)
-// Reference:
-// PRESS, W. H., TEUKOLSKY, S. A., VETTERLING, W. T., AND FLANNERY, B. P.
-// 1992. Numerical Recipes in C: The Art of Scientific Computing, 2nd ed. Cambridge University Press.
-//
-// This file is part of FreeImage 3
-//
-// COVERED CODE IS PROVIDED UNDER THIS LICENSE ON AN "AS IS" BASIS, WITHOUT WARRANTY
-// OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES
-// THAT THE COVERED CODE IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE
-// OR NON-INFRINGING. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE COVERED
-// CODE IS WITH YOU. SHOULD ANY COVERED CODE PROVE DEFECTIVE IN ANY RESPECT, YOU (NOT
-// THE INITIAL DEVELOPER OR ANY OTHER CONTRIBUTOR) ASSUME THE COST OF ANY NECESSARY
-// SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL
-// PART OF THIS LICENSE. NO USE OF ANY COVERED CODE IS AUTHORIZED HEREUNDER EXCEPT UNDER
-// THIS DISCLAIMER.
-//
-// Use at your own risk!
-// ==========================================================
-
-#include "FreeImage.h"
-#include "Utilities.h"
-#include "ToneMapping.h"
-
-static const int NPRE = 1; // Number of relaxation sweeps before ...
-static const int NPOST = 1; // ... and after the coarse-grid correction is computed
-static const int NGMAX = 15; // Maximum number of grids
-
-/**
-Copy src into dst
-*/
-static inline void fmg_copyArray(FIBITMAP *dst, FIBITMAP *src) {
- memcpy(FreeImage_GetBits(dst), FreeImage_GetBits(src), FreeImage_GetHeight(dst) * FreeImage_GetPitch(dst));
-}
-
-/**
-Fills src with zeros
-*/
-static inline void fmg_fillArrayWithZeros(FIBITMAP *src) {
- memset(FreeImage_GetBits(src), 0, FreeImage_GetHeight(src) * FreeImage_GetPitch(src));
-}
-
-/**
-Half-weighting restriction. nc is the coarse-grid dimension. The fine-grid solution is input in
-uf[0..2*nc-2][0..2*nc-2], the coarse-grid solution is returned in uc[0..nc-1][0..nc-1].
-*/
-static void fmg_restrict(FIBITMAP *UC, FIBITMAP *UF, int nc) {
- int row_uc, row_uf, col_uc, col_uf;
-
- const int uc_pitch = FreeImage_GetPitch(UC) / sizeof(float);
- const int uf_pitch = FreeImage_GetPitch(UF) / sizeof(float);
-
- float *uc_bits = (float*)FreeImage_GetBits(UC);
- const float *uf_bits = (float*)FreeImage_GetBits(UF);
-
- // interior points
- {
- float *uc_scan = uc_bits + uc_pitch;
- for (row_uc = 1, row_uf = 2; row_uc < nc-1; row_uc++, row_uf += 2) {
- const float *uf_scan = uf_bits + row_uf * uf_pitch;
- for (col_uc = 1, col_uf = 2; col_uc < nc-1; col_uc++, col_uf += 2) {
- // calculate
- // UC(row_uc, col_uc) =
- // 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) ]
- float *uc_pixel = uc_scan + col_uc;
- const float *uf_center = uf_scan + col_uf;
- *uc_pixel = 0.5F * *uf_center + 0.125F * ( *(uf_center + uf_pitch) + *(uf_center - uf_pitch) + *(uf_center + 1) + *(uf_center - 1) );
- }
- uc_scan += uc_pitch;
- }
- }
- // boundary points
- const int ncc = 2*nc-1;
- {
- /*
- calculate the following:
- for (row_uc = 0, row_uf = 0; row_uc < nc; row_uc++, row_uf += 2) {
- UC(row_uc, 0) = UF(row_uf, 0);
- UC(row_uc, nc-1) = UF(row_uf, ncc-1);
- }
- */
- float *uc_scan = uc_bits;
- for (row_uc = 0, row_uf = 0; row_uc < nc; row_uc++, row_uf += 2) {
- const float *uf_scan = uf_bits + row_uf * uf_pitch;
- uc_scan[0] = uf_scan[0];
- uc_scan[nc-1] = uf_scan[ncc-1];
- uc_scan += uc_pitch;
- }
- }
- {
- /*
- calculate the following:
- for (col_uc = 0, col_uf = 0; col_uc < nc; col_uc++, col_uf += 2) {
- UC(0, col_uc) = UF(0, col_uf);
- UC(nc-1, col_uc) = UF(ncc-1, col_uf);
- }
- */
- float *uc_scan_top = uc_bits;
- float *uc_scan_bottom = uc_bits + (nc-1)*uc_pitch;
- const float *uf_scan_top = uf_bits + (ncc-1)*uf_pitch;
- const float *uf_scan_bottom = uf_bits;
- for (col_uc = 0, col_uf = 0; col_uc < nc; col_uc++, col_uf += 2) {
- uc_scan_top[col_uc] = uf_scan_top[col_uf];
- uc_scan_bottom[col_uc] = uf_scan_bottom[col_uf];
- }
- }
-}
-
-/**
-Solution of the model problem on the coarsest grid, where h = 1/2 .
-The right-hand side is input
-in rhs[0..2][0..2] and the solution is returned in u[0..2][0..2].
-*/
-static void fmg_solve(FIBITMAP *U, FIBITMAP *RHS) {
- // fill U with zeros
- fmg_fillArrayWithZeros(U);
- // calculate U(1, 1) = -h*h*RHS(1, 1)/4.0 where h = 1/2
- float *u_scan = (float*)FreeImage_GetScanLine(U, 1);
- const float *rhs_scan = (float*)FreeImage_GetScanLine(RHS, 1);
- u_scan[1] = -rhs_scan[1] / 16;
-}
-
-/**
-Coarse-to-fine prolongation by bilinear interpolation. nf is the fine-grid dimension. The coarsegrid
-solution is input as uc[0..nc-1][0..nc-1], where nc = nf/2 + 1. The fine-grid solution is
-returned in uf[0..nf-1][0..nf-1].
-*/
-static void fmg_prolongate(FIBITMAP *UF, FIBITMAP *UC, int nf) {
- int row_uc, row_uf, col_uc, col_uf;
-
- const int uf_pitch = FreeImage_GetPitch(UF) / sizeof(float);
- const int uc_pitch = FreeImage_GetPitch(UC) / sizeof(float);
-
- float *uf_bits = (float*)FreeImage_GetBits(UF);
- const float *uc_bits = (float*)FreeImage_GetBits(UC);
-
- // do elements that are copies
- {
- const int nc = nf/2 + 1;
-
- float *uf_scan = uf_bits;
- const float *uc_scan = uc_bits;
- for (row_uc = 0; row_uc < nc; row_uc++) {
- for (col_uc = 0, col_uf = 0; col_uc < nc; col_uc++, col_uf += 2) {
- // calculate UF(2*row_uc, col_uf) = UC(row_uc, col_uc);
- uf_scan[col_uf] = uc_scan[col_uc];
- }
- uc_scan += uc_pitch;
- uf_scan += 2 * uf_pitch;
- }
- }
- // do odd-numbered columns, interpolating vertically
- {
- for(row_uf = 1; row_uf < nf-1; row_uf += 2) {
- float *uf_scan = uf_bits + row_uf * uf_pitch;
- for (col_uf = 0; col_uf < nf; col_uf += 2) {
- // calculate UF(row_uf, col_uf) = 0.5 * ( UF(row_uf+1, col_uf) + UF(row_uf-1, col_uf) )
- uf_scan[col_uf] = 0.5F * ( *(uf_scan + uf_pitch + col_uf) + *(uf_scan - uf_pitch + col_uf) );
- }
- }
- }
- // do even-numbered columns, interpolating horizontally
- {
- float *uf_scan = uf_bits;
- for(row_uf = 0; row_uf < nf; row_uf++) {
- for (col_uf = 1; col_uf < nf-1; col_uf += 2) {
- // calculate UF(row_uf, col_uf) = 0.5 * ( UF(row_uf, col_uf+1) + UF(row_uf, col_uf-1) )
- uf_scan[col_uf] = 0.5F * ( uf_scan[col_uf + 1] + uf_scan[col_uf - 1] );
- }
- uf_scan += uf_pitch;
- }
- }
-}
-
-/**
-Red-black Gauss-Seidel relaxation for model problem. Updates the current value of the solution
-u[0..n-1][0..n-1], using the right-hand side function rhs[0..n-1][0..n-1].
-*/
-static void fmg_relaxation(FIBITMAP *U, FIBITMAP *RHS, int n) {
- int row, col, ipass, isw, jsw;
- const float h = 1.0F / (n - 1);
- const float h2 = h*h;
-
- const int u_pitch = FreeImage_GetPitch(U) / sizeof(float);
- const int rhs_pitch = FreeImage_GetPitch(RHS) / sizeof(float);
-
- float *u_bits = (float*)FreeImage_GetBits(U);
- const float *rhs_bits = (float*)FreeImage_GetBits(RHS);
-
- for (ipass = 0, jsw = 1; ipass < 2; ipass++, jsw = 3-jsw) { // Red and black sweeps
- float *u_scan = u_bits + u_pitch;
- const float *rhs_scan = rhs_bits + rhs_pitch;
- for (row = 1, isw = jsw; row < n-1; row++, isw = 3-isw) {
- for (col = isw; col < n-1; col += 2) {
- // Gauss-Seidel formula
- // calculate U(row, col) =
- // 0.25 * [ U(row+1, col) + U(row-1, col) + U(row, col+1) + U(row, col-1) - h2 * RHS(row, col) ]
- float *u_center = u_scan + col;
- const float *rhs_center = rhs_scan + col;
- *u_center = *(u_center + u_pitch) + *(u_center - u_pitch) + *(u_center + 1) + *(u_center - 1);
- *u_center -= h2 * *rhs_center;
- *u_center *= 0.25F;
- }
- u_scan += u_pitch;
- rhs_scan += rhs_pitch;
- }
- }
-}
-
-/**
-Returns minus the residual for the model problem. Input quantities are u[0..n-1][0..n-1] and
-rhs[0..n-1][0..n-1], while res[0..n-1][0..n-1] is returned.
-*/
-static void fmg_residual(FIBITMAP *RES, FIBITMAP *U, FIBITMAP *RHS, int n) {
- int row, col;
-
- const float h = 1.0F / (n-1);
- const float h2i = 1.0F / (h*h);
-
- const int res_pitch = FreeImage_GetPitch(RES) / sizeof(float);
- const int u_pitch = FreeImage_GetPitch(U) / sizeof(float);
- const int rhs_pitch = FreeImage_GetPitch(RHS) / sizeof(float);
-
- float *res_bits = (float*)FreeImage_GetBits(RES);
- const float *u_bits = (float*)FreeImage_GetBits(U);
- const float *rhs_bits = (float*)FreeImage_GetBits(RHS);
-
- // interior points
- {
- float *res_scan = res_bits + res_pitch;
- const float *u_scan = u_bits + u_pitch;
- const float *rhs_scan = rhs_bits + rhs_pitch;
- for (row = 1; row < n-1; row++) {
- for (col = 1; col < n-1; col++) {
- // calculate RES(row, col) =
- // -h2i * [ U(row+1, col) + U(row-1, col) + U(row, col+1) + U(row, col-1) - 4 * U(row, col) ] + RHS(row, col);
- float *res_center = res_scan + col;
- const float *u_center = u_scan + col;
- const float *rhs_center = rhs_scan + col;
- *res_center = *(u_center + u_pitch) + *(u_center - u_pitch) + *(u_center + 1) + *(u_center - 1) - 4 * *u_center;
- *res_center *= -h2i;
- *res_center += *rhs_center;
- }
- res_scan += res_pitch;
- u_scan += u_pitch;
- rhs_scan += rhs_pitch;
- }
- }
-
- // boundary points
- {
- memset(FreeImage_GetScanLine(RES, 0), 0, FreeImage_GetPitch(RES));
- memset(FreeImage_GetScanLine(RES, n-1), 0, FreeImage_GetPitch(RES));
- float *left = res_bits;
- float *right = res_bits + (n-1);
- for(int k = 0; k < n; k++) {
- *left = 0;
- *right = 0;
- left += res_pitch;
- right += res_pitch;
- }
- }
-}
-
-/**
-Does coarse-to-fine interpolation and adds result to uf. nf is the fine-grid dimension. The
-coarse-grid solution is input as uc[0..nc-1][0..nc-1], where nc = nf/2+1. The fine-grid solution
-is returned in uf[0..nf-1][0..nf-1]. res[0..nf-1][0..nf-1] is used for temporary storage.
-*/
-static void fmg_addint(FIBITMAP *UF, FIBITMAP *UC, FIBITMAP *RES, int nf) {
- fmg_prolongate(RES, UC, nf);
-
- const int uf_pitch = FreeImage_GetPitch(UF) / sizeof(float);
- const int res_pitch = FreeImage_GetPitch(RES) / sizeof(float);
-
- float *uf_bits = (float*)FreeImage_GetBits(UF);
- const float *res_bits = (float*)FreeImage_GetBits(RES);
-
- for(int row = 0; row < nf; row++) {
- for(int col = 0; col < nf; col++) {
- // calculate UF(row, col) = UF(row, col) + RES(row, col);
- uf_bits[col] += res_bits[col];
- }
- uf_bits += uf_pitch;
- res_bits += res_pitch;
- }
-}
-
-/**
-Full Multigrid Algorithm for solution of linear elliptic equation, here the model problem (19.0.6).
-On input u[0..n-1][0..n-1] contains the right-hand side ñ, while on output it returns the solution.
-The dimension n must be of the form 2^j + 1 for some integer j. (j is actually the number of
-grid levels used in the solution, called ng below.) ncycle is the number of V-cycles to be
-used at each level.
-*/
-static BOOL fmg_mglin(FIBITMAP *U, int n, int ncycle) {
- int j, jcycle, jj, jpost, jpre, nf, ngrid;
-
- FIBITMAP **IRHO = NULL;
- FIBITMAP **IU = NULL;
- FIBITMAP **IRHS = NULL;
- FIBITMAP **IRES = NULL;
-
- int ng = 0; // number of allocated grids
-
-// --------------------------------------------------------------------------
-
-#define _CREATE_ARRAY_GRID_(array, array_size) \
- array = (FIBITMAP**)malloc(array_size * sizeof(FIBITMAP*));\
- if(!array) throw(1);\
- memset(array, 0, array_size * sizeof(FIBITMAP*))
-
-#define _FREE_ARRAY_GRID_(array, array_size) \
- if(NULL != array) {\
- for(int k = 0; k < array_size; k++) {\
- if(NULL != array[k]) {\
- FreeImage_Unload(array[k]); array[k] = NULL;\
- }\
- }\
- free(array);\
- }
-
-// --------------------------------------------------------------------------
-
- try {
- int nn = n;
- // check grid size and grid levels
- while (nn >>= 1) ng++;
- if (n != 1 + (1L << ng)) {
- FreeImage_OutputMessageProc(FIF_UNKNOWN, "Multigrid algorithm: n = %d, while n-1 must be a power of 2.", n);
- throw(1);
- }
- if (ng > NGMAX) {
- FreeImage_OutputMessageProc(FIF_UNKNOWN, "Multigrid algorithm: ng = %d while NGMAX = %d, increase NGMAX.", ng, NGMAX);
- throw(1);
- }
- // allocate grid arrays
- {
- _CREATE_ARRAY_GRID_(IRHO, ng);
- _CREATE_ARRAY_GRID_(IU, ng);
- _CREATE_ARRAY_GRID_(IRHS, ng);
- _CREATE_ARRAY_GRID_(IRES, ng);
- }
-
- nn = n/2 + 1;
- ngrid = ng - 2;
-
- // allocate storage for r.h.s. on grid (ng - 2) ...
- IRHO[ngrid] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
- if(!IRHO[ngrid]) throw(1);
-
- // ... and fill it by restricting from the fine grid
- fmg_restrict(IRHO[ngrid], U, nn);
-
- // similarly allocate storage and fill r.h.s. on all coarse grids.
- while (nn > 3) {
- nn = nn/2 + 1;
- ngrid--;
- IRHO[ngrid] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
- if(!IRHO[ngrid]) throw(1);
- fmg_restrict(IRHO[ngrid], IRHO[ngrid+1], nn);
- }
-
- nn = 3;
-
- IU[0] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
- if(!IU[0]) throw(1);
- IRHS[0] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
- if(!IRHS[0]) throw(1);
-
- // initial solution on coarsest grid
- fmg_solve(IU[0], IRHO[0]);
- // irho[0] no longer needed ...
- FreeImage_Unload(IRHO[0]); IRHO[0] = NULL;
-
- ngrid = ng;
-
- // nested iteration loop
- for (j = 1; j < ngrid; j++) {
- nn = 2*nn - 1;
-
- IU[j] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
- if(!IU[j]) throw(1);
- IRHS[j] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
- if(!IRHS[j]) throw(1);
- IRES[j] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
- if(!IRES[j]) throw(1);
-
- fmg_prolongate(IU[j], IU[j-1], nn);
-
- // interpolate from coarse grid to next finer grid
-
- // set up r.h.s.
- fmg_copyArray(IRHS[j], j != (ngrid - 1) ? IRHO[j] : U);
-
- // V-cycle loop
- for (jcycle = 0; jcycle < ncycle; jcycle++) {
- nf = nn;
- // downward stoke of the V
- for (jj = j; jj >= 1; jj--) {
- // pre-smoothing
- for (jpre = 0; jpre < NPRE; jpre++) {
- fmg_relaxation(IU[jj], IRHS[jj], nf);
- }
- fmg_residual(IRES[jj], IU[jj], IRHS[jj], nf);
- nf = nf/2 + 1;
- // restriction of the residual is the next r.h.s.
- fmg_restrict(IRHS[jj-1], IRES[jj], nf);
- // zero for initial guess in next relaxation
- fmg_fillArrayWithZeros(IU[jj-1]);
- }
- // bottom of V: solve on coarsest grid
- fmg_solve(IU[0], IRHS[0]);
- nf = 3;
- // upward stroke of V.
- for (jj = 1; jj <= j; jj++) {
- nf = 2*nf - 1;
- // use res for temporary storage inside addint
- fmg_addint(IU[jj], IU[jj-1], IRES[jj], nf);
- // post-smoothing
- for (jpost = 0; jpost < NPOST; jpost++) {
- fmg_relaxation(IU[jj], IRHS[jj], nf);
- }
- }
- }
- }
-
- // return solution in U
- fmg_copyArray(U, IU[ngrid-1]);
-
- // delete allocated arrays
- _FREE_ARRAY_GRID_(IRES, ng);
- _FREE_ARRAY_GRID_(IRHS, ng);
- _FREE_ARRAY_GRID_(IU, ng);
- _FREE_ARRAY_GRID_(IRHO, ng);
-
- return TRUE;
-
- } catch(int) {
- // delete allocated arrays
- _FREE_ARRAY_GRID_(IRES, ng);
- _FREE_ARRAY_GRID_(IRHS, ng);
- _FREE_ARRAY_GRID_(IU, ng);
- _FREE_ARRAY_GRID_(IRHO, ng);
-
- return FALSE;
- }
-}
-
-// --------------------------------------------------------------------------
-
-/**
-Poisson solver based on a multigrid algorithm.
-This routine solves a Poisson equation, remap result pixels to [0..1] and returns the solution.
-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,
-where j is such that 2^j is the nearest larger dimension corresponding to MAX(image width, image height).
-@param Laplacian Laplacian image
-@param ncycle Number of cycles in the multigrid algorithm (usually 2 or 3)
-@return Returns the solved PDE equations if successful, returns NULL otherwise
-*/
-FIBITMAP* DLL_CALLCONV
-FreeImage_MultigridPoissonSolver(FIBITMAP *Laplacian, int ncycle) {
- if(!FreeImage_HasPixels(Laplacian)) return NULL;
-
- int width = FreeImage_GetWidth(Laplacian);
- int height = FreeImage_GetHeight(Laplacian);
-
- // get nearest larger dimension length that is acceptable by the algorithm
- int n = MAX(width, height);
- int size = 0;
- while((n >>= 1) > 0) size++;
- if((1 << size) < MAX(width, height)) {
- size++;
- }
- // size must be of the form 2^j + 1 for some integer j
- size = 1 + (1 << size);
-
- // allocate a temporary square image I
- FIBITMAP *I = FreeImage_AllocateT(FIT_FLOAT, size, size);
- if(!I) return NULL;
-
- // copy Laplacian into I and shift pixels to create a boundary
- FreeImage_Paste(I, Laplacian, 1, 1, 255);
-
- // solve the PDE equation
- fmg_mglin(I, size, ncycle);
-
- // shift pixels back
- FIBITMAP *U = FreeImage_Copy(I, 1, 1, width + 1, height + 1);
- FreeImage_Unload(I);
-
- // remap pixels to [0..1]
- NormalizeY(U, 0, 1);
-
- // copy metadata from src to dst
- FreeImage_CloneMetadata(U, Laplacian);
-
- // return the integrated image
- return U;
-}
-
+// ==========================================================
+// Poisson solver based on a full multigrid algorithm
+//
+// Design and implementation by
+// - Hervé Drolon (drolon@infonie.fr)
+// Reference:
+// PRESS, W. H., TEUKOLSKY, S. A., VETTERLING, W. T., AND FLANNERY, B. P.
+// 1992. Numerical Recipes in C: The Art of Scientific Computing, 2nd ed. Cambridge University Press.
+//
+// This file is part of FreeImage 3
+//
+// COVERED CODE IS PROVIDED UNDER THIS LICENSE ON AN "AS IS" BASIS, WITHOUT WARRANTY
+// OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES
+// THAT THE COVERED CODE IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE
+// OR NON-INFRINGING. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE COVERED
+// CODE IS WITH YOU. SHOULD ANY COVERED CODE PROVE DEFECTIVE IN ANY RESPECT, YOU (NOT
+// THE INITIAL DEVELOPER OR ANY OTHER CONTRIBUTOR) ASSUME THE COST OF ANY NECESSARY
+// SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL
+// PART OF THIS LICENSE. NO USE OF ANY COVERED CODE IS AUTHORIZED HEREUNDER EXCEPT UNDER
+// THIS DISCLAIMER.
+//
+// Use at your own risk!
+// ==========================================================
+
+#include "FreeImage.h"
+#include "Utilities.h"
+#include "ToneMapping.h"
+
+static const int NPRE = 1; // Number of relaxation sweeps before ...
+static const int NPOST = 1; // ... and after the coarse-grid correction is computed
+static const int NGMAX = 15; // Maximum number of grids
+
+/**
+Copy src into dst
+*/
+static inline void fmg_copyArray(FIBITMAP *dst, FIBITMAP *src) {
+ memcpy(FreeImage_GetBits(dst), FreeImage_GetBits(src), FreeImage_GetHeight(dst) * FreeImage_GetPitch(dst));
+}
+
+/**
+Fills src with zeros
+*/
+static inline void fmg_fillArrayWithZeros(FIBITMAP *src) {
+ memset(FreeImage_GetBits(src), 0, FreeImage_GetHeight(src) * FreeImage_GetPitch(src));
+}
+
+/**
+Half-weighting restriction. nc is the coarse-grid dimension. The fine-grid solution is input in
+uf[0..2*nc-2][0..2*nc-2], the coarse-grid solution is returned in uc[0..nc-1][0..nc-1].
+*/
+static void fmg_restrict(FIBITMAP *UC, FIBITMAP *UF, int nc) {
+ int row_uc, row_uf, col_uc, col_uf;
+
+ const int uc_pitch = FreeImage_GetPitch(UC) / sizeof(float);
+ const int uf_pitch = FreeImage_GetPitch(UF) / sizeof(float);
+
+ float *uc_bits = (float*)FreeImage_GetBits(UC);
+ const float *uf_bits = (float*)FreeImage_GetBits(UF);
+
+ // interior points
+ {
+ float *uc_scan = uc_bits + uc_pitch;
+ for (row_uc = 1, row_uf = 2; row_uc < nc-1; row_uc++, row_uf += 2) {
+ const float *uf_scan = uf_bits + row_uf * uf_pitch;
+ for (col_uc = 1, col_uf = 2; col_uc < nc-1; col_uc++, col_uf += 2) {
+ // calculate
+ // UC(row_uc, col_uc) =
+ // 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) ]
+ float *uc_pixel = uc_scan + col_uc;
+ const float *uf_center = uf_scan + col_uf;
+ *uc_pixel = 0.5F * *uf_center + 0.125F * ( *(uf_center + uf_pitch) + *(uf_center - uf_pitch) + *(uf_center + 1) + *(uf_center - 1) );
+ }
+ uc_scan += uc_pitch;
+ }
+ }
+ // boundary points
+ const int ncc = 2*nc-1;
+ {
+ /*
+ calculate the following:
+ for (row_uc = 0, row_uf = 0; row_uc < nc; row_uc++, row_uf += 2) {
+ UC(row_uc, 0) = UF(row_uf, 0);
+ UC(row_uc, nc-1) = UF(row_uf, ncc-1);
+ }
+ */
+ float *uc_scan = uc_bits;
+ for (row_uc = 0, row_uf = 0; row_uc < nc; row_uc++, row_uf += 2) {
+ const float *uf_scan = uf_bits + row_uf * uf_pitch;
+ uc_scan[0] = uf_scan[0];
+ uc_scan[nc-1] = uf_scan[ncc-1];
+ uc_scan += uc_pitch;
+ }
+ }
+ {
+ /*
+ calculate the following:
+ for (col_uc = 0, col_uf = 0; col_uc < nc; col_uc++, col_uf += 2) {
+ UC(0, col_uc) = UF(0, col_uf);
+ UC(nc-1, col_uc) = UF(ncc-1, col_uf);
+ }
+ */
+ float *uc_scan_top = uc_bits;
+ float *uc_scan_bottom = uc_bits + (nc-1)*uc_pitch;
+ const float *uf_scan_top = uf_bits + (ncc-1)*uf_pitch;
+ const float *uf_scan_bottom = uf_bits;
+ for (col_uc = 0, col_uf = 0; col_uc < nc; col_uc++, col_uf += 2) {
+ uc_scan_top[col_uc] = uf_scan_top[col_uf];
+ uc_scan_bottom[col_uc] = uf_scan_bottom[col_uf];
+ }
+ }
+}
+
+/**
+Solution of the model problem on the coarsest grid, where h = 1/2 .
+The right-hand side is input
+in rhs[0..2][0..2] and the solution is returned in u[0..2][0..2].
+*/
+static void fmg_solve(FIBITMAP *U, FIBITMAP *RHS) {
+ // fill U with zeros
+ fmg_fillArrayWithZeros(U);
+ // calculate U(1, 1) = -h*h*RHS(1, 1)/4.0 where h = 1/2
+ float *u_scan = (float*)FreeImage_GetScanLine(U, 1);
+ const float *rhs_scan = (float*)FreeImage_GetScanLine(RHS, 1);
+ u_scan[1] = -rhs_scan[1] / 16;
+}
+
+/**
+Coarse-to-fine prolongation by bilinear interpolation. nf is the fine-grid dimension. The coarsegrid
+solution is input as uc[0..nc-1][0..nc-1], where nc = nf/2 + 1. The fine-grid solution is
+returned in uf[0..nf-1][0..nf-1].
+*/
+static void fmg_prolongate(FIBITMAP *UF, FIBITMAP *UC, int nf) {
+ int row_uc, row_uf, col_uc, col_uf;
+
+ const int uf_pitch = FreeImage_GetPitch(UF) / sizeof(float);
+ const int uc_pitch = FreeImage_GetPitch(UC) / sizeof(float);
+
+ float *uf_bits = (float*)FreeImage_GetBits(UF);
+ const float *uc_bits = (float*)FreeImage_GetBits(UC);
+
+ // do elements that are copies
+ {
+ const int nc = nf/2 + 1;
+
+ float *uf_scan = uf_bits;
+ const float *uc_scan = uc_bits;
+ for (row_uc = 0; row_uc < nc; row_uc++) {
+ for (col_uc = 0, col_uf = 0; col_uc < nc; col_uc++, col_uf += 2) {
+ // calculate UF(2*row_uc, col_uf) = UC(row_uc, col_uc);
+ uf_scan[col_uf] = uc_scan[col_uc];
+ }
+ uc_scan += uc_pitch;
+ uf_scan += 2 * uf_pitch;
+ }
+ }
+ // do odd-numbered columns, interpolating vertically
+ {
+ for(row_uf = 1; row_uf < nf-1; row_uf += 2) {
+ float *uf_scan = uf_bits + row_uf * uf_pitch;
+ for (col_uf = 0; col_uf < nf; col_uf += 2) {
+ // calculate UF(row_uf, col_uf) = 0.5 * ( UF(row_uf+1, col_uf) + UF(row_uf-1, col_uf) )
+ uf_scan[col_uf] = 0.5F * ( *(uf_scan + uf_pitch + col_uf) + *(uf_scan - uf_pitch + col_uf) );
+ }
+ }
+ }
+ // do even-numbered columns, interpolating horizontally
+ {
+ float *uf_scan = uf_bits;
+ for(row_uf = 0; row_uf < nf; row_uf++) {
+ for (col_uf = 1; col_uf < nf-1; col_uf += 2) {
+ // calculate UF(row_uf, col_uf) = 0.5 * ( UF(row_uf, col_uf+1) + UF(row_uf, col_uf-1) )
+ uf_scan[col_uf] = 0.5F * ( uf_scan[col_uf + 1] + uf_scan[col_uf - 1] );
+ }
+ uf_scan += uf_pitch;
+ }
+ }
+}
+
+/**
+Red-black Gauss-Seidel relaxation for model problem. Updates the current value of the solution
+u[0..n-1][0..n-1], using the right-hand side function rhs[0..n-1][0..n-1].
+*/
+static void fmg_relaxation(FIBITMAP *U, FIBITMAP *RHS, int n) {
+ int row, col, ipass, isw, jsw;
+ const float h = 1.0F / (n - 1);
+ const float h2 = h*h;
+
+ const int u_pitch = FreeImage_GetPitch(U) / sizeof(float);
+ const int rhs_pitch = FreeImage_GetPitch(RHS) / sizeof(float);
+
+ float *u_bits = (float*)FreeImage_GetBits(U);
+ const float *rhs_bits = (float*)FreeImage_GetBits(RHS);
+
+ for (ipass = 0, jsw = 1; ipass < 2; ipass++, jsw = 3-jsw) { // Red and black sweeps
+ float *u_scan = u_bits + u_pitch;
+ const float *rhs_scan = rhs_bits + rhs_pitch;
+ for (row = 1, isw = jsw; row < n-1; row++, isw = 3-isw) {
+ for (col = isw; col < n-1; col += 2) {
+ // Gauss-Seidel formula
+ // calculate U(row, col) =
+ // 0.25 * [ U(row+1, col) + U(row-1, col) + U(row, col+1) + U(row, col-1) - h2 * RHS(row, col) ]
+ float *u_center = u_scan + col;
+ const float *rhs_center = rhs_scan + col;
+ *u_center = *(u_center + u_pitch) + *(u_center - u_pitch) + *(u_center + 1) + *(u_center - 1);
+ *u_center -= h2 * *rhs_center;
+ *u_center *= 0.25F;
+ }
+ u_scan += u_pitch;
+ rhs_scan += rhs_pitch;
+ }
+ }
+}
+
+/**
+Returns minus the residual for the model problem. Input quantities are u[0..n-1][0..n-1] and
+rhs[0..n-1][0..n-1], while res[0..n-1][0..n-1] is returned.
+*/
+static void fmg_residual(FIBITMAP *RES, FIBITMAP *U, FIBITMAP *RHS, int n) {
+ int row, col;
+
+ const float h = 1.0F / (n-1);
+ const float h2i = 1.0F / (h*h);
+
+ const int res_pitch = FreeImage_GetPitch(RES) / sizeof(float);
+ const int u_pitch = FreeImage_GetPitch(U) / sizeof(float);
+ const int rhs_pitch = FreeImage_GetPitch(RHS) / sizeof(float);
+
+ float *res_bits = (float*)FreeImage_GetBits(RES);
+ const float *u_bits = (float*)FreeImage_GetBits(U);
+ const float *rhs_bits = (float*)FreeImage_GetBits(RHS);
+
+ // interior points
+ {
+ float *res_scan = res_bits + res_pitch;
+ const float *u_scan = u_bits + u_pitch;
+ const float *rhs_scan = rhs_bits + rhs_pitch;
+ for (row = 1; row < n-1; row++) {
+ for (col = 1; col < n-1; col++) {
+ // calculate RES(row, col) =
+ // -h2i * [ U(row+1, col) + U(row-1, col) + U(row, col+1) + U(row, col-1) - 4 * U(row, col) ] + RHS(row, col);
+ float *res_center = res_scan + col;
+ const float *u_center = u_scan + col;
+ const float *rhs_center = rhs_scan + col;
+ *res_center = *(u_center + u_pitch) + *(u_center - u_pitch) + *(u_center + 1) + *(u_center - 1) - 4 * *u_center;
+ *res_center *= -h2i;
+ *res_center += *rhs_center;
+ }
+ res_scan += res_pitch;
+ u_scan += u_pitch;
+ rhs_scan += rhs_pitch;
+ }
+ }
+
+ // boundary points
+ {
+ memset(FreeImage_GetScanLine(RES, 0), 0, FreeImage_GetPitch(RES));
+ memset(FreeImage_GetScanLine(RES, n-1), 0, FreeImage_GetPitch(RES));
+ float *left = res_bits;
+ float *right = res_bits + (n-1);
+ for(int k = 0; k < n; k++) {
+ *left = 0;
+ *right = 0;
+ left += res_pitch;
+ right += res_pitch;
+ }
+ }
+}
+
+/**
+Does coarse-to-fine interpolation and adds result to uf. nf is the fine-grid dimension. The
+coarse-grid solution is input as uc[0..nc-1][0..nc-1], where nc = nf/2+1. The fine-grid solution
+is returned in uf[0..nf-1][0..nf-1]. res[0..nf-1][0..nf-1] is used for temporary storage.
+*/
+static void fmg_addint(FIBITMAP *UF, FIBITMAP *UC, FIBITMAP *RES, int nf) {
+ fmg_prolongate(RES, UC, nf);
+
+ const int uf_pitch = FreeImage_GetPitch(UF) / sizeof(float);
+ const int res_pitch = FreeImage_GetPitch(RES) / sizeof(float);
+
+ float *uf_bits = (float*)FreeImage_GetBits(UF);
+ const float *res_bits = (float*)FreeImage_GetBits(RES);
+
+ for(int row = 0; row < nf; row++) {
+ for(int col = 0; col < nf; col++) {
+ // calculate UF(row, col) = UF(row, col) + RES(row, col);
+ uf_bits[col] += res_bits[col];
+ }
+ uf_bits += uf_pitch;
+ res_bits += res_pitch;
+ }
+}
+
+/**
+Full Multigrid Algorithm for solution of linear elliptic equation, here the model problem (19.0.6).
+On input u[0..n-1][0..n-1] contains the right-hand side ñ, while on output it returns the solution.
+The dimension n must be of the form 2^j + 1 for some integer j. (j is actually the number of
+grid levels used in the solution, called ng below.) ncycle is the number of V-cycles to be
+used at each level.
+*/
+static BOOL fmg_mglin(FIBITMAP *U, int n, int ncycle) {
+ int j, jcycle, jj, jpost, jpre, nf, ngrid;
+
+ FIBITMAP **IRHO = NULL;
+ FIBITMAP **IU = NULL;
+ FIBITMAP **IRHS = NULL;
+ FIBITMAP **IRES = NULL;
+
+ int ng = 0; // number of allocated grids
+
+// --------------------------------------------------------------------------
+
+#define _CREATE_ARRAY_GRID_(array, array_size) \
+ array = (FIBITMAP**)malloc(array_size * sizeof(FIBITMAP*));\
+ if(!array) throw(1);\
+ memset(array, 0, array_size * sizeof(FIBITMAP*))
+
+#define _FREE_ARRAY_GRID_(array, array_size) \
+ if(NULL != array) {\
+ for(int k = 0; k < array_size; k++) {\
+ if(NULL != array[k]) {\
+ FreeImage_Unload(array[k]); array[k] = NULL;\
+ }\
+ }\
+ free(array);\
+ }
+
+// --------------------------------------------------------------------------
+
+ try {
+ int nn = n;
+ // check grid size and grid levels
+ while (nn >>= 1) ng++;
+ if (n != 1 + (1L << ng)) {
+ FreeImage_OutputMessageProc(FIF_UNKNOWN, "Multigrid algorithm: n = %d, while n-1 must be a power of 2.", n);
+ throw(1);
+ }
+ if (ng > NGMAX) {
+ FreeImage_OutputMessageProc(FIF_UNKNOWN, "Multigrid algorithm: ng = %d while NGMAX = %d, increase NGMAX.", ng, NGMAX);
+ throw(1);
+ }
+ // allocate grid arrays
+ {
+ _CREATE_ARRAY_GRID_(IRHO, ng);
+ _CREATE_ARRAY_GRID_(IU, ng);
+ _CREATE_ARRAY_GRID_(IRHS, ng);
+ _CREATE_ARRAY_GRID_(IRES, ng);
+ }
+
+ nn = n/2 + 1;
+ ngrid = ng - 2;
+
+ // allocate storage for r.h.s. on grid (ng - 2) ...
+ IRHO[ngrid] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
+ if(!IRHO[ngrid]) throw(1);
+
+ // ... and fill it by restricting from the fine grid
+ fmg_restrict(IRHO[ngrid], U, nn);
+
+ // similarly allocate storage and fill r.h.s. on all coarse grids.
+ while (nn > 3) {
+ nn = nn/2 + 1;
+ ngrid--;
+ IRHO[ngrid] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
+ if(!IRHO[ngrid]) throw(1);
+ fmg_restrict(IRHO[ngrid], IRHO[ngrid+1], nn);
+ }
+
+ nn = 3;
+
+ IU[0] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
+ if(!IU[0]) throw(1);
+ IRHS[0] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
+ if(!IRHS[0]) throw(1);
+
+ // initial solution on coarsest grid
+ fmg_solve(IU[0], IRHO[0]);
+ // irho[0] no longer needed ...
+ FreeImage_Unload(IRHO[0]); IRHO[0] = NULL;
+
+ ngrid = ng;
+
+ // nested iteration loop
+ for (j = 1; j < ngrid; j++) {
+ nn = 2*nn - 1;
+
+ IU[j] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
+ if(!IU[j]) throw(1);
+ IRHS[j] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
+ if(!IRHS[j]) throw(1);
+ IRES[j] = FreeImage_AllocateT(FIT_FLOAT, nn, nn);
+ if(!IRES[j]) throw(1);
+
+ fmg_prolongate(IU[j], IU[j-1], nn);
+
+ // interpolate from coarse grid to next finer grid
+
+ // set up r.h.s.
+ fmg_copyArray(IRHS[j], j != (ngrid - 1) ? IRHO[j] : U);
+
+ // V-cycle loop
+ for (jcycle = 0; jcycle < ncycle; jcycle++) {
+ nf = nn;
+ // downward stoke of the V
+ for (jj = j; jj >= 1; jj--) {
+ // pre-smoothing
+ for (jpre = 0; jpre < NPRE; jpre++) {
+ fmg_relaxation(IU[jj], IRHS[jj], nf);
+ }
+ fmg_residual(IRES[jj], IU[jj], IRHS[jj], nf);
+ nf = nf/2 + 1;
+ // restriction of the residual is the next r.h.s.
+ fmg_restrict(IRHS[jj-1], IRES[jj], nf);
+ // zero for initial guess in next relaxation
+ fmg_fillArrayWithZeros(IU[jj-1]);
+ }
+ // bottom of V: solve on coarsest grid
+ fmg_solve(IU[0], IRHS[0]);
+ nf = 3;
+ // upward stroke of V.
+ for (jj = 1; jj <= j; jj++) {
+ nf = 2*nf - 1;
+ // use res for temporary storage inside addint
+ fmg_addint(IU[jj], IU[jj-1], IRES[jj], nf);
+ // post-smoothing
+ for (jpost = 0; jpost < NPOST; jpost++) {
+ fmg_relaxation(IU[jj], IRHS[jj], nf);
+ }
+ }
+ }
+ }
+
+ // return solution in U
+ fmg_copyArray(U, IU[ngrid-1]);
+
+ // delete allocated arrays
+ _FREE_ARRAY_GRID_(IRES, ng);
+ _FREE_ARRAY_GRID_(IRHS, ng);
+ _FREE_ARRAY_GRID_(IU, ng);
+ _FREE_ARRAY_GRID_(IRHO, ng);
+
+ return TRUE;
+
+ } catch(int) {
+ // delete allocated arrays
+ _FREE_ARRAY_GRID_(IRES, ng);
+ _FREE_ARRAY_GRID_(IRHS, ng);
+ _FREE_ARRAY_GRID_(IU, ng);
+ _FREE_ARRAY_GRID_(IRHO, ng);
+
+ return FALSE;
+ }
+}
+
+// --------------------------------------------------------------------------
+
+/**
+Poisson solver based on a multigrid algorithm.
+This routine solves a Poisson equation, remap result pixels to [0..1] and returns the solution.
+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,
+where j is such that 2^j is the nearest larger dimension corresponding to MAX(image width, image height).
+@param Laplacian Laplacian image
+@param ncycle Number of cycles in the multigrid algorithm (usually 2 or 3)
+@return Returns the solved PDE equations if successful, returns NULL otherwise
+*/
+FIBITMAP* DLL_CALLCONV
+FreeImage_MultigridPoissonSolver(FIBITMAP *Laplacian, int ncycle) {
+ if(!FreeImage_HasPixels(Laplacian)) return NULL;
+
+ int width = FreeImage_GetWidth(Laplacian);
+ int height = FreeImage_GetHeight(Laplacian);
+
+ // get nearest larger dimension length that is acceptable by the algorithm
+ int n = MAX(width, height);
+ int size = 0;
+ while((n >>= 1) > 0) size++;
+ if((1 << size) < MAX(width, height)) {
+ size++;
+ }
+ // size must be of the form 2^j + 1 for some integer j
+ size = 1 + (1 << size);
+
+ // allocate a temporary square image I
+ FIBITMAP *I = FreeImage_AllocateT(FIT_FLOAT, size, size);
+ if(!I) return NULL;
+
+ // copy Laplacian into I and shift pixels to create a boundary
+ FreeImage_Paste(I, Laplacian, 1, 1, 255);
+
+ // solve the PDE equation
+ fmg_mglin(I, size, ncycle);
+
+ // shift pixels back
+ FIBITMAP *U = FreeImage_Copy(I, 1, 1, width + 1, height + 1);
+ FreeImage_Unload(I);
+
+ // remap pixels to [0..1]
+ NormalizeY(U, 0, 1);
+
+ // copy metadata from src to dst
+ FreeImage_CloneMetadata(U, Laplacian);
+
+ // return the integrated image
+ return U;
+}
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 19:49:45 +0000