view release on metacpan or  search on metacpan

src/Source/FreeImageToolkit/MultigridPoissonSolver.cpp  view on Meta::CPAN

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