// Function Prototypes; 
//===============================================================================
BOOL laplacian2D(Raw2D* pSrc, int bound=1);
// Find Laplacian of 'src'.
// (resize us as needed) 
// bound%2=0: use Von Neumann boundary cond.
// bound%2=1: use Dirichlet (default).
BOOL poisson2D_naive(Raw2D *pLap, int maxCount, int bound=1, 
					 double maxErr=0.001);
// Slow Poisson Solver for 2D images. Given
// desired Laplacian pLap, iteratively
// find an image that fits it well;
// Convergence rate is O(N^2)
// bound%2=0: use Von Neumann bounds
// bound%2=1: use Dirichlet bounds(DEFAULT)
BOOL poisson2D_GaussSeidel(Raw2D *pLap, int maxCount, int bound=1,
						   double maxErr=0.001);
// Still slow--2X as fast as 'naive' method.
BOOL poisson2D_GaussSeidelRB(Raw2D *pLap, int maxCount, int bound=1,
							 double maxErr=0.001);
// Better: Red-black checkerboard..
BOOL poisson2D_SORfix(Raw2D *pLap, int maxCount, int bound=1,
					  double maxErr=0.001);
// Better: convergence rate approx O(N) with
// a scale factor.
BOOL poisson2D_SOR(Raw2D *pLap, int maxCount, int bound=1,
				   double maxErr=0.001);
// Better: convergence rate approx O(N) 
// with Chebyshev acceleration.
//===============================================================================



// Function Bodies:
//===============================================================================


BOOL Raw2D::laplacian2D(Raw2D* pSrc, int bound)
//------------------------------------------------------------------------------
// Laplacian estimates from N,S,E,W neighbors. Recall that Laplacian for scalar
// function F is another scalar, written del^2 F or del-dot-del F, and is the
// sum of the 2nd derivs in each spatial direction:
// Laplacian(F) = d^2 F/dx + d^2 F/dy.  If we assume our pixel sample positions
// are at integer locations, then we can approximate the Laplacian using forward
// differences.  If we call F(x,y)==C, and define the north,south,east,west,
// neighbors as N,S,E,W, then Laplacian at x,y is (N+S+E+W)-4C.  
//
// bound%2=0: use Von Neumann boundary condition: presume gradient across the
// outer edge of the image is zero (as if image extended to infinity by copying
// the outermost image pixels endlessly)
// bound%2=1: use Dirichlet boundary condition (default): presume all pixels
// outside of the image are zero-valued.  Return FALSE on error.
{

#define JT_BORD		1			// size of border used for bound. conditions

	Raw2D jtmp;
	int i,j,imax,jmax;
	PIXTYPE lap;

	if(pSrc==NULL)
	{
		BEEP;
		DPF("\nRaw2D::laplacian2D() given NULL input!\n");
		return(FALSE);
	}
	if(pSrc->getYsize()<=0 || pSrc->getXsize()<=0)
	{
		BEEP;					// complain; exit.
		DPF("\nRaw2D::laplacian2D() given zero-sized input!\n");
		return(FALSE);
	}
	sizer(pSrc);			// make sure we match the size of the input img.
	// Copy 'src' into 'tmp', but with a JT_BORD-pixel wide border on all sides
	imax = pSrc->getXsize() + JT_BORD;//Find address of 1st pix in far border)
	jmax = pSrc->getYsize() + JT_BORD;
	jtmp.sizer(JT_BORD + imax, JT_BORD + jmax);

	if(bound==1)				// Initialize the buffer:
	{							// for (default) Dirichlet boundary conditions,
		jtmp.putConst((PIXTYPE)0.0);	// set image & surrounds to zero.
	}
	for(j=JT_BORD; j<jmax;j++)			// copy image into center of 'tmp'
	{
		for(i=JT_BORD; i<imax; i++)
		{
			jtmp.put(i,j,pSrc->get(i-JT_BORD,j-JT_BORD));
		}
	}
	if(bound==0)			// for the VonNeumann conditions,
	{						// copy some pixels to fill the borders:
		for(j=0; j<JT_BORD; j++)	// For top & bottom border scanlines,
		{
			for(i=JT_BORD; i<imax; i++)	// and within the horiz. borders,
			{						// copy nearest scanline within border
				jtmp.put(i,j,     jtmp.get(i,JT_BORD));	// top 
				jtmp.put(i,jmax+j,jtmp.get(i,jmax-1));    // bottom
			}
		}
		for(j=JT_BORD; j<jmax; j++)// for left & right side borders,
		{
			for(i=0; i<JT_BORD; i++)
			{						// copy nearest column within border
				jtmp.put(i,j,     jtmp.get(JT_BORD,j));// left side
				jtmp.put(i+imax,j,jtmp.get(imax-1 ,j));// right
			}
		}
		for(j=0; j<JT_BORD; j++)	// Now fill in all the corner values;
		{
			for(i=0; i<JT_BORD; i++)
			{
				jtmp.put(i,     j,     jtmp.get(JT_BORD,JT_BORD));// lower left,
				jtmp.put(i+imax,j,     jtmp.get(imax-1, JT_BORD));// lower right,
				jtmp.put(i,j+jmax,     jtmp.get(JT_BORD,jmax-1));// upper left,
				jtmp.put(i+imax,j+jmax,jtmp.get(imax-1, jmax-1));// upper right.
			}
		}	
	}
	// ---OK, finally:compute the Laplacian for the pixels within the borders:
	for(j=JT_BORD; j<jmax; j++)
	{
		for(i=JT_BORD; i<imax; i++)
		{
			lap = jtmp.get(i+1,j) + 
				jtmp.get(i-1,j) + 
				jtmp.get(i,j+1) + 
				jtmp.get(i,j-1) - (PIXTYPE)4.0 * jtmp.get(i,j);
			put(i-JT_BORD,j-JT_BORD,lap);	// record the laplacian as output.
		}
	}
	jtmp.wipe();							// discard the temporary buffer.
	return(TRUE);
#undef JT_BORD			//locals only!
}

BOOL Raw2D::poisson2D_naive(Raw2D *pLap, int maxCount, int bound, double maxErr)
//------------------------------------------------------------------------------
// Straightforward but slow-to-converge Poisson Solver. Starting with a constant
// valued image, iteratively adjust it until its Laplacian matches pLap, using
// the Jacobi Iteration method. (see Numerical Recipes).
//		Recall that the Laplacian of F is just the sum of the 2nd derivs in x,y
// directions: Laplacian = d^2F/dx^2 + d^2F/dy^2.
//		Given an NxM image of scalar-valued Laplacian of intensity pLap,
// make an output double-buffer of size (N+2*JT_BORD)x(M+2*JT_BORD), big enough 
// to include a JT_BORD-sized border of extra pixels on the image's top, bottom, 
// left and right sides. Initialize to a constant-valued image.  For each non-
// border pixel, find the Laplacian; if it does not match the desired value 
// found in pLap, adjust the buffer pixel's intensity to improve its Laplacian.
//  'Boundary Conditions' let you define the Laplacian for pixels
// along the top,bottom, and sides of images:
//
// If bound%2=1, (DEFAULT):
// Use Dirichlet boundary conditions: Presume the image is surrounded by an
// unbounded field of zero-valued pixels. 
// If bound%2=0 : (Not recommended)
// Uses Von Neumann boundary conditions: assumes gradients across across the
// image boundaries are zero, but unchanged in tangent direction).  
//
//		NOTES: Dirichlet conditions give you a unique solution; Von Neumann
//  conditions have an unknown constant offset. Poisson solutions can 
// have persistent large error near boundaries if the laplacian source image was
//  computed with unknown boundary conditions, or was a section of a larger 
//  Laplacian image. You can improve results by enclosing the Laplacian in a 
// border of zero-valued pixels, and use either VonNeumann or Dirichlet..
// 
//	(fcnGroup::fcnCurv1() finds 2nd derivatives w/ VonNeumann boundaries;)
// **(fcnGroup::fcnCurv2() finds 2nd derivatives w/Dirichlet boundaries.)**
//
// 06/08/2004 Created--J. Tumblin
{
	int i,j,k,imax,jmax;
	Raw2D buf[2];					// double-buffer for computing output;
	int newB,oldB,tmpB;				// indices for buf[]; 'old' is the previously
	// computed output, 'new' is the new output we
	// are computing from values in 'old'.

	PIXTYPE lapNow,lapErr;			// current Laplacian, & how far wrong it is.
	PIXTYPE lapAvg,lapWanted;		// current local average; desired laplacian
	double err2max;					// largest squared error found.

#define JT_BORD   1				// Width of zero-valued border around image.
#define JT_FRAC	  0.5f			// Fraction of Laplacian to use(for non-Jacobi).
#define JT_JACOBI TRUE			// Computing method: If TRUE, use the fast and
	// simple Jacobi iteration, where we just find
	// the avg. of the 4 neighboring pixels and
	// subtract the Laplacian/4 to get a new pixel.
	// If FALSE, use the slower, more-obvious method
	// --find (desired-actual)Laplacian/4,
	// --multiply it by JT_FRAC, subtract from the
	// current pixel value to make new pixel value.

	if(pLap==NULL)				// input does not exist? 
	{
		BEEP;					// complain; exit.
		DPF("Raw2D::Poisson2D_naive() given NULL input!\n");
		return(FALSE);
	}
	if(pLap->getYsize()<=0 || pLap->getXsize()<=0)
	{
		BEEP;					// complain; exit.
		DPF("Raw2D::Poisson2D_naive() given zero-sized input!\n");
		return(FALSE);
	}
	// Make a pair of double-buffered images that
	// can hold pLap plus a border of 'JT_BORD'
	// pixels wide on all 4 sides:
	buf[0].sizer(2*JT_BORD + pLap->getXsize(),2*JT_BORD + pLap->getYsize());
	buf[1].sizer(2*JT_BORD + pLap->getXsize(),2*JT_BORD + pLap->getYsize());
	newB = 0; oldB = 1;			// init. vars that tell us which buffer holds 
	// old results vs. new values we're computing.
	imax = JT_BORD + pLap->getXsize();// Find largest address within the border)
	jmax = JT_BORD + pLap->getYsize();
	if(bound==1)				// Initialize the buffers:
	{							// for (default) Dirichlet boundary conditions,
		buf[newB].putConst((PIXTYPE)0.0);	// set image & surrounds to zero.
		buf[oldB].putConst((PIXTYPE)0.0);
	}
	else
	{							// for (not reccommended) VonNeumann conditions,
		buf[newB].putConst((PIXTYPE)0.5);	// grow darker & lighter from the
		buf[oldB].putConst((PIXTYPE)0.5);	// mid-level gray value of 0.5
	}
	for(k=0; k<maxCount; k++)	// Iteratively modify the output buffer so the
	{							// Laplacian at each pix matches pLap values:
		err2max = 0.0;				// clear the max-error-finder for this pass.
		for(j=JT_BORD; j<jmax; j++)	// For all pixels within the border,
		{
			for(i=JT_BORD; i<imax; i++)
			{
#if JT_JACOBI==FALSE
				// Even slower than Jacobi iteration, but a bit more obvious.
				// find current Laplacian value:
				lapNow = (buf[oldB].get(i+1,j) +//(sum the 4 nearest
					buf[oldB].get(i-1,j) +	// neighbor pixels,
					buf[oldB].get(i,j+1) +
					buf[oldB].get(i,j-1)) -
					(PIXTYPE)4.0 * buf[oldB].get(i,j); // minus 4*center.
				// Find how much it differs from desired value:
				lapErr = pLap->get(i-JT_BORD,j-JT_BORD) - lapNow;
				// correct a fraction of it to make new output:
				buf[newB].put(i,j,buf[oldB].get(i,j) 
					- JT_FRAC*lapErr*(PIXTYPE)0.25);

				//				buf[newB].put(i  ,j,buf[oldB].get(i  ,j)-JT_FRAC*lapErr*(PIXTYPE)0.5);
				//				buf[newB].put(i+1,j,buf[oldB].get(i+1,j)+JT_FRAC*lapErr*(PIXTYPE)0.125);
				//				buf[newB].put(i-1,j,buf[oldB].get(i+1,j)+JT_FRAC*lapErr*(PIXTYPE)0.125);
				//				buf[newB].put(i,j+1,buf[oldB].get(i,j+1)+JT_FRAC*lapErr*(PIXTYPE)0.125);
				//				buf[newB].put(i,j-1,buf[oldB].get(i,j-1)+JT_FRAC*lapErr*(PIXTYPE)0.125);

#else			// Use Jacobi iteration: it's simpler and faster:
				lapWanted = pLap->get(i-JT_BORD, j-JT_BORD);
				lapAvg = (buf[oldB].get(i+1,j) + buf[oldB].get(i-1,j) +
					buf[oldB].get(i,j+1) + buf[oldB].get(i,j-1));
				// before we change things, measure error...
				lapNow = lapAvg - (PIXTYPE)4.0*buf[oldB].get(i,j);
				lapErr = lapWanted - lapNow;
				// Correct exactly 1/4 of the laplacian error;
				// we can ignore current pixel value and get:
				buf[newB].put(i,j, (PIXTYPE)0.25 * (lapAvg - lapWanted));
				// (compute the error).
#endif
				if(err2max < (double)(lapErr*lapErr))
				{				// find image's largest error;
					err2max = (double)(lapErr*lapErr);
				}
			}
			if(j%200==0) DPF(".");// print little dots every 200 scanlines...
		}
		/////Set Boundary conditions-------------------------------------------
		// bound%2==1 means Dirichlet Boundary conditions:-----------
		// ***We already have them!
		//		Border pixels were initialized to zero, and remain unchanged.
		if(bound%2==0)		// VonNeumann Boundary conditions----------
		{						// update the border pixels so that their 
			// values match the adjacent pixel:
			for(j=0; j<JT_BORD; j++)	// For top & bottom border scanlines,
			{
				for(i=JT_BORD; i<imax; i++)	// and within the horiz. borders,
				{						// copy nearest scanline within border
					buf[newB].put(i,j, buf[newB].get(i,JT_BORD));	// top 
					buf[newB].put(i,jmax+j,buf[newB].get(i,jmax-1));// bottom
				}
			}
			for(j=JT_BORD; j<jmax; j++)// for left & right side borders,
			{
				for(i=0; i<JT_BORD; i++)
				{						// copy nearest column within border
					buf[newB].put(i,j,buf[newB].get(JT_BORD,j));// left side
					buf[newB].put(i+imax,j,buf[newB].get(imax-1,j));// right
				}
			}
			for(j=0; j<JT_BORD; j++)	// Now fill in all the corner values;
			{
				for(i=0; i<JT_BORD; i++)
				{
					buf[newB].put(i,j,
						buf[newB].get(JT_BORD,JT_BORD));	// lower left,
					buf[newB].put(i+imax,j,
						buf[newB].get(imax-1,JT_BORD));		// lower right,
					buf[newB].put(i,j+jmax,
						buf[newB].get(JT_BORD,jmax-1));		// upper left,
					buf[newB].put(i+imax,j+jmax,
						buf[newB].get(imax-1,jmax-1));		// upper right.
				}
			}
		}
		//-----END boundary-value setting

		// 'newb' buffer is now complete.
		tmpB = oldB;			// 'newb' buffer is now complete. Swap buffers 
		oldB = newB;			// so that newB becomes the old buffer.
		newB = tmpB;
		if(err2max < maxErr*maxErr)// Are we below the error maximum?
		{
			break;				// Yes. Stop.  Answer is in 'oldB' buffer.
		}						// No; keep going.
	}
	if(bound==0) DPF("VonNeumann--"); else DPF("Dirichlet(DEFAULT)--");
	if(k==maxCount)				// did we stop before error < errMax?
	{							// Yes; report it.
		DPF("Raw2D::Poisson2D_naive() stopped at err=%f, not err=%f!\n",(float)sqrt(err2max),(float)maxErr);
	}
	else
	{
		DPF("Raw2D::Poisson2D_naive() stopped at err=%f after %d counts\n",(float)sqrt(err2max),k);
	}
	///HIDE BORDER/////////////
	sizer(pLap);				// Copy result to output image, without border;
	for(j=JT_BORD; j<jmax; j++)
	{
		for(i=JT_BORD; i<imax; i++)
		{
			put(i-JT_BORD,j-JT_BORD,	buf[oldB].get(i,j));
		}
	}

	//END HIDE BORDER////////////////
	//diagnostic only: Show border://///////////////////
	//	wipecopy(&(buf[oldB]));
	// end diagnostic 'show border'/////////////////////
	return(TRUE);

#undef JT_BORD					// locals only!
#undef JT_FRAC
#undef JT_JACOBI
}


BOOL Raw2D::poisson2D_GaussSeidel(Raw2D *pLap, int maxCount, int bound, double maxErr)
//------------------------------------------------------------------------------
// Poisson Solver using GaussSeidel method, raster-scan. It is ALMOST the same 
// solver, but uses in-place calculations without double-buffering. Instead of 
// just 'old' pixel values to compute new ones, it mixes old and new values.
//  'Boundary Conditions' let you define the Laplacian for pixels
// along the top,bottom, and sides of images:
//
// If bound%2=1, (DEFAULT):
// Use Dirichlet boundary conditions: Presume the image is surrounded by an
// unbounded field of zero-valued pixels. 
// If bound%2=0 : (Not recommended)
// Uses VonNeumann boundary conditions: assumes gradients across across the
// image boundaries are zero, but unchanged in tangent direction).  
//
//		NOTES: Dirichlet conditions give you a unique solution;  VonNeumann 
// conditions have an unknown constant offset. Poisson solutions can have 
// persistent large error near boundaries if the laplacian source image was
//  computed with unknown boundary conditions, or was a section of a larger 
//  Laplacian image. You can improve results by enclosing the Laplacian in a 
// border of zero-valued pixels, and use either VonNeumann or Dirichlet..
// 
//	(fcnGroup::fcnCurv1() finds 2nd derivatives w/ VonNeumann boundaries;)
// **(fcnGroup::fcnCurv2() finds 2nd derivatives w/Dirichlet boundaries.)**
//
// 06/08/2004 Created--J. Tumblin
{
	int i,j,k,imax,jmax;
	Raw2D buf;						// single-buffer for computing output;
	PIXTYPE lapNow,lapErr;			// current Laplacian, & how far wrong it is.
	PIXTYPE lapAvg,lapWanted;		// current local average; desired laplacian
	double err2max;					// largest squared error found.

#define JT_BORD   1				// Width of zero-valued border around image.

	if(pLap==NULL)				// input does not exist? 
	{
		BEEP;					// complain; exit.
		DPF("Raw2D::Poisson2D_GaussSeidel() given NULL input!\n");
		return(FALSE);
	}
	if(pLap->getYsize()<=0 || pLap->getXsize()<=0)
	{
		BEEP;					// complain; exit.
		DPF("Raw2D::Poisson2D_GaussSeidel() given zero-sized input!\n");
		return(FALSE);
	}
	// Make a single output image buffer that
	// can hold pLap plus a border of 'JT_BORD'
	// pixels wide on all 4 sides:
	buf.sizer(2*JT_BORD + pLap->getXsize(),2*JT_BORD + pLap->getYsize());
	imax = JT_BORD + pLap->getXsize();// Find largest address within the border)
	jmax = JT_BORD + pLap->getYsize();
	if(bound==1)				// Initialize the buffers:
	{							// for (default) Dirichlet boundary conditions,
		buf.putConst((PIXTYPE)0.0);	// set image & surrounds to zero.
	}
	else
	{							// for (not reccommended) VonNeumann conditions,
		buf.putConst((PIXTYPE)0.5);	// grow darker & lighter from the
	}								// mid-level gray value of 0.5
	for(k=0; k<maxCount; k++)	// Iteratively modify the output buffer so the
	{							// Laplacian at each pix matches pLap values:
		err2max = 0.0;				// clear the max-error-finder for this pass.
		for(j=JT_BORD; j<jmax; j++)	// For all pixels within the border,
		{
			for(i=JT_BORD; i<imax; i++)
			{
				// Using Jacobi iteration in-place: 
				lapWanted = pLap->get(i-JT_BORD, j-JT_BORD);
				lapAvg = (buf.get(i+1,j) + buf.get(i-1,j) +
					buf.get(i,j+1) + buf.get(i,j-1));
				lapNow = lapAvg - (PIXTYPE)4.0*buf.get(i,j);//current Laplacian
				lapErr = lapWanted - lapNow;
				// If we correct 1/4 of the laplacian error,
				// then computing the new pixel value is simple:
				buf.put(i,j, (PIXTYPE)0.25 * (lapAvg - lapWanted));
				// (compute the error).
				if(err2max < (double)(lapErr*lapErr))
				{				// find image's largest error;
					err2max = (double)(lapErr*lapErr);
				}
			}// end of one scanline,
			if(j%200==0) DPF(".");// print little dots every 200 scanlines...
		}// end of one pass,
		/////Set Boundary conditions-------------------------------------------
		// bound%2==1 means Dirichlet Boundary conditions:------------
		// We already have them! Border pixels initialized to zero & unchanged.
		if(bound%2==0)		// Von Neumann Boundary conditions----------
		{						// update the border pixels so that their 
			// values match the adjacent pixel:
			for(j=0; j<JT_BORD; j++)	// For top & bottom border scanlines,
			{
				for(i=JT_BORD; i<imax; i++)	// and within the horiz. borders,
				{						// copy nearest scanline within border
					buf.put(i,j, buf.get(i,JT_BORD));	// top 
					buf.put(i,jmax+j,buf.get(i,jmax-1));// bottom
				}
			}
			for(j=JT_BORD; j<jmax; j++)// for left & right side borders,
			{
				for(i=0; i<JT_BORD; i++)
				{						// copy nearest column within border
					buf.put(i,j,buf.get(JT_BORD,j));// left side
					buf.put(i+imax,j,buf.get(imax-1,j));// right
				}
			}
			for(j=0; j<JT_BORD; j++)	// Now fill in all the corner values;
			{
				for(i=0; i<JT_BORD; i++)
				{
					buf.put(i,j,buf.get(JT_BORD,JT_BORD));		// lower left,
					buf.put(i+imax,j,buf.get(imax-1,JT_BORD));	// lower right,
					buf.put(i,j+jmax,buf.get(JT_BORD,jmax-1));	// upper left,
					buf.put(i+imax,j+jmax,buf.get(imax-1,jmax-1));// upper right.
				}
			}
		}
		//-----END boundary-value setting-----------
		if(err2max < maxErr*maxErr)// Are we below the error maximum?
		{
			break;				// Yes. Stop.  Answer is in 'oldB' buffer.
		}						// No; keep going.
	}	// end of all iterations.

	if(bound==0) DPF("VonNeumann--"); else DPF("Dirichlet(DEFAULT)--");
	if(k==maxCount)				// did we stop before error < errMax?
	{							// Yes; report it.
		DPF("Raw2D::Poisson2D_GaussSeidel() stopped at err=%f, not err=%f!\n",(float)sqrt(err2max),(float)maxErr);
	}
	else
	{
		DPF("Raw2D::Poisson2D_GaussSiedel() stopped at err=%f after %d counts\n",(float)sqrt(err2max),k);
	}
	///HIDE BORDER/////////////
	sizer(pLap);				// Copy result to output image, without border;
	for(j=JT_BORD; j<jmax; j++)
	{
		for(i=JT_BORD; i<imax; i++)
		{
			put(i-JT_BORD,j-JT_BORD,buf.get(i,j));
		}
	}

	//END HIDE BORDER////////////////
	//diagnostic only: Show border://///////////////////
	//	wipecopy(&(buf[oldB]));
	// end diagnostic 'show border'/////////////////////
	return(TRUE);

#undef JT_BORD					// locals only!
#undef JT_FRAC
}
BOOL Raw2D::poisson2D_GaussSeidelRB(Raw2D *pLap, int maxCount, int bound, double maxErr)
//------------------------------------------------------------------------------
// Poisson Solver using Gauss-Seidel method, red-black scan (e.g. each iteration
// split into 2 halves, each updating a checkerboard-like pattern of pixels). 
// It is ALMOST the same solver as naive (Jacobi method), but uses in-place 
// calculations without double-buffering. Instead of 
// just 'old' pixel values to compute new ones, it mixes old and new values.
//  'Boundary Conditions' let you define the Laplacian for pixels
// along the top,bottom, and sides of images:
//
// If bound%2=1, (DEFAULT):
// Use Dirichlet boundary conditions: Presume the image is surrounded by an
// unbounded field of zero-valued pixels. 
// If bound%2=0 : (Not recommended)
// Uses VonNeumann boundary conditions: assumes gradients across across the
// image boundaries are zero, but unchanged in tangent direction).  
//
//		NOTES: Dirichlet conditions give you a unique solution;  VonNeumann 
// conditions have an unknown constant offset. Poisson solutions can 
// have persistent large error near boundaries if the laplacian source image was
//  computed with unknown boundary conditions, or was a section of a larger 
//  Laplacian image. You can improve results by enclosing the Laplacian in a 
// border of zero-valued pixels, and use either VonNeumann or Dirichlet..
// 
//	(fcnGroup::fcnCurv1() finds 2nd derivatives w/ VonNeumann boundaries;)
// **(fcnGroup::fcnCurv2() finds 2nd derivatives w/ Dirichlet boundaries.)**
//
// 06/08/2004 Created--J. Tumblin
{
	int i,j,k,imax,jmax;
	int isOdd,chk;
	Raw2D buf;						// single-buffer for computing output;
	PIXTYPE lapNow,lapErr;			// current Laplacian, & how far wrong it is.
	PIXTYPE lapAvg,lapWanted;		// current local average; desired laplacian
	double err2max;					// largest squared error found.

#define JT_BORD   1				// Width of zero-valued border around image.

	if(pLap==NULL)				// input does not exist? 
	{
		BEEP;					// complain; exit.
		DPF("Raw2D::Poisson2D_GaussSeidelRB() given NULL input!\n");
		return(FALSE);
	}
	if(pLap->getYsize()<=0 || pLap->getXsize()<=0)
	{
		BEEP;					// complain; exit.
		DPF("Raw2D::Poisson2D_GaussSeidelRB() given zero-sized input!\n");
		return(FALSE);
	}
	// Make a single output image buffer that
	// can hold pLap plus a border of 'JT_BORD'
	// pixels wide on all 4 sides:
	buf.sizer(2*JT_BORD + pLap->getXsize(),2*JT_BORD + pLap->getYsize());
	imax = JT_BORD + pLap->getXsize();// Find largest address within the border)
	jmax = JT_BORD + pLap->getYsize();
	if(bound==1)				// Initialize the buffers:
	{							// for (default) Dirichlet boundary conditions,
		buf.putConst((PIXTYPE)0.0);	// set image & surrounds to zero.
	}
	else
	{							// for (not reccommended) VonNeumann conditions,
		buf.putConst((PIXTYPE)0.5);	// grow darker & lighter from the
	}								// mid-level gray value of 0.5
	for(k=0; k<maxCount; k++)	// Iteratively modify the output buffer so the
	{							// Laplacian at each pix matches pLap values:
		err2max = 0.0;				// clear the max-error-finder for this pass.
		for(isOdd=0; isOdd<2; isOdd++)	// for each half-sweep of the image:
		{
			chk = isOdd;			// start on pixel 0;
			for(j=JT_BORD; j<jmax; j++)// update every scanline,
			{
				for(i=JT_BORD+chk; i<imax; i+=2)	// but only half the pixels,
				{	// Using Jacobi iteration in-place: 
					lapWanted = pLap->get(i-JT_BORD, j-JT_BORD);
					lapAvg = (buf.get(i+1,j) + buf.get(i-1,j) +
						buf.get(i,j+1) + buf.get(i,j-1));
					lapNow = lapAvg - (PIXTYPE)4.0*buf.get(i,j);//current Laplacian
					lapErr = lapWanted - lapNow;
					// If we correct 1/4 of the laplacian error,
					// then computing the new pixel value is simple:
					buf.put(i,j, (PIXTYPE)0.25 * (lapAvg - lapWanted));
					// (compute the error).
					if(err2max < (double)(lapErr*lapErr))
					{				// find image's largest error;
						err2max = (double)(lapErr*lapErr);
					}// end of one pixel,
				}// end of one scanline,
				chk = (chk+1)%2;	// toggle 'chk' on each scanline
			}// end of one half-sweep,
			if(j%200==0) DPF(".");// print little dots every 200 scanlines...
			/////Set Boundary conditions---------------------------------------
			// bound%2==1 means Dirichlet Boundary conditions:-----------
			// We already have them! Border pixels initialized to zero & unchanged.
			if(bound%2==0)		// VonNeumann Boundary conditions----------
			{						// update the border pixels so that their 
				// values match the adjacent pixel:
				for(j=0; j<JT_BORD; j++)	// For top & bottom border scanlines,
				{
					for(i=JT_BORD; i<imax; i+=2)
						// and within the horiz. borders,
					{					// copy nearest scanline within border
						buf.put(i,j, buf.get(i,JT_BORD));	// top 
						buf.put(i,jmax+j,buf.get(i,jmax-1));// bottom
					}
				}
				for(j=JT_BORD; j<jmax; j++)
					// for left & right side borders,
				{
					for(i=0; i<JT_BORD; i++)
					{						// copy nearest column within border
						buf.put(i,j,buf.get(JT_BORD,j));// left side
						buf.put(i+imax,j,buf.get(imax-1,j));// right
					}
				}
				for(j=0; j<JT_BORD; j++)	// Now fill in all the corner values;
				{
					for(i=0; i<JT_BORD; i++)
					{
						buf.put(i,j,buf.get(JT_BORD,JT_BORD));		// lower left,
						buf.put(i+imax,j,buf.get(imax-1,JT_BORD));	// lower right,
						buf.put(i,j+jmax,buf.get(JT_BORD,jmax-1));	// upper left,
						buf.put(i+imax,j+jmax,buf.get(imax-1,jmax-1));// upper right.
					}
				}
			}
			//-----END boundary-value setting-----------
		}							// go on to the next iteration.

		if(err2max < maxErr*maxErr)// Are we below the error maximum?
		{
			break;				// Yes. Stop.  Answer is in 'oldB' buffer.
		}						// No; keep going.
	}
	if(bound==0) DPF("VonNeumann--"); else DPF("Dirichlet(DEFAULT)--");
	if(k==maxCount)				// did we stop before error < errMax?
	{							// Yes; report it.
		DPF("Raw2D::Poisson2D_GaussSeidelRB() stopped at err=%f, not err=%f!\n",(float)sqrt(err2max),(float)maxErr);
	}
	else
	{
		DPF("Raw2D::Poisson2D_GaussSiedelRB() stopped at err=%f after %d counts\n",(float)sqrt(err2max),k);
	}
	///HIDE BORDER/////////////
	sizer(pLap);				// Copy result to output image, without border;
	for(j=JT_BORD; j<jmax; j++)
	{
		for(i=JT_BORD; i<imax; i++)
		{
			put(i-JT_BORD,j-JT_BORD,buf.get(i,j));
		}
	}

	//END HIDE BORDER////////////////
	//diagnostic only: Show border://///////////////////
	//	wipecopy(&(buf[oldB]));
	// end diagnostic 'show border'/////////////////////
	return(TRUE);

#undef JT_BORD					// locals only!
#undef JT_FRAC
}

BOOL Raw2D::poisson2D_SORfix(Raw2D *pLap, int maxCount, int bound, double maxErr)
//------------------------------------------------------------------------------
// Poisson Solver using Successive Over-Relaxation method--converges O(N), while
// 'naive' and GaussSeidel converges O(N^2). It is very similar to the Gauss-
// SeidelRB() solver, but the correction for Laplacian error is exaggerated by
// an a fixed factor.
//  
//If bound%2=1, (DEFAULT):
// Use Dirichlet boundary conditions: Presume the image is surrounded by an
// unbounded field of zero-valued pixels. 
// If bound%2=0 : (Not recommended)
// Uses VonNeumann boundary conditions: assumes gradients across across the
// image boundaries are zero, but unchanged in tangent direction).  
//
//		NOTES: Dirichlet conditions give you a unique solution; VonNeumann 
// conditions have an unknown constant offset. Poisson solutions can 
// have persistent large error near boundaries if the Laplacian source image was
//  computed with unknown boundary conditions, or was a section of a larger 
//  Laplacian image. You can improve results by enclosing the Laplacian in a 
// border of zero-valued pixels, and use either VonNeumann or Dirichlet..
// 
//	(fcnGroup::fcnCurv1() finds 2nd derivatives w/ VonNeumann boundaries;)
// **(fcnGroup::fcnCurv2() finds 2nd derivatives w/Dirichlet boundaries.)**
//
// 06/08/2004 Created--J. Tumblin
{
	int i,j,k,imax,jmax;
	int isOdd,chk;					// Checkerboard generator.
	Raw2D buf;						// single-buffer for computing output;
	PIXTYPE lapNow,lapErr;			// current Laplacian, & how far wrong it is.
	PIXTYPE lapAvg,lapWanted;		// current local average; desired laplacian
	double err2max;					// largest squared error found.
	PIXTYPE chebW;					// Chebyshev polynom. weight (see Num. Recipes)
#define JT_BORD   1				// Width of zero-valued border around image.

	if(pLap==NULL)				// input does not exist? 
	{
		BEEP;					// complain; exit.
		DPF("Raw2D::Poisson2D_SOR() given NULL input!\n");
		return(FALSE);
	}
	if(pLap->getYsize()<=0 || pLap->getXsize()<=0)
	{
		BEEP;					// complain; exit.
		DPF("Raw2D::Poisson2D_SOR() given zero-sized input!\n");
		return(FALSE);
	}
	// Make a single output image buffer that
	// can hold pLap plus a border of 'JT_BORD'
	// pixels wide on all 4 sides:
	buf.sizer(2*JT_BORD + pLap->getXsize(),2*JT_BORD + pLap->getYsize());
	imax = JT_BORD + pLap->getXsize();// Find largest address within the border)
	jmax = JT_BORD + pLap->getYsize();
	if(bound==1)				// Initialize the buffers:
	{							// for (default) Dirichlet boundary conditions,
		buf.putConst((PIXTYPE)0.0);	// set image & surrounds to zero.
	}
	else
	{							// for (not reccommended) VonNeumann conditions,
		buf.putConst((PIXTYPE)0.5);	// grow darker & lighter from the
	}								// mid-level gray value of 0.5
	chebW = (PIXTYPE)1.99;
	for(k=0; k<maxCount; k++)	// Iteratively modify the output buffer so the
	{							// Laplacian at each pix matches pLap values:
		err2max = 0.0;				// clear the max-error-finder for this pass.
		for(isOdd=0; isOdd<2; isOdd++)	// for each half-sweep of the image:
		{
			chk = isOdd;
			for(j=JT_BORD; j<jmax; j++)// update every scanline,
			{
				for(i=JT_BORD+chk; i<imax; i+=2)	// but only half the pixels,
				{	// Using Jacobi iteration in-place: 
					lapWanted = pLap->get(i-JT_BORD, j-JT_BORD);
					lapAvg = (buf.get(i+1,j) + buf.get(i-1,j) +
						buf.get(i,j+1) + buf.get(i,j-1));
					lapNow = lapAvg - (PIXTYPE)4.0*buf.get(i,j);//current Laplacian
					lapErr = lapWanted - lapNow;
					// If we correct 1/4 of the laplacian error,
					// then computing the new pixel value is simple:
					buf.put(i,j, buf.get(i,j) - (PIXTYPE)0.25 *chebW * lapErr);
					// (compute the error).
					//					if(err2max < (double)(lapErr*lapErr))
					//					{				// find image's largest error;
					//						err2max = (double)(lapErr*lapErr);
					//					}// end of one pixel,
					err2max += (double)(lapErr*lapErr);		// sum-of-square err
				}// end of one scanline,
				chk = (chk+1)%2;	// toggle 'chk' on each scanline
			}// end of one half-sweep,
			//if(j%100==0) DPF(".");// print little dots every 100 scanlines...
			/////Set Boundary conditions---------------------------------------
			// bound%2==1 means Dirichlet Boundary conditions:
			// We already have them! Border pixels initialized to zero & unchanged.
			if(bound%2==0)		// VonNeumann Boundary conditions----------
			{						// update the border pixels so that their 
				// values match the adjacent pixel:
				for(j=0; j<JT_BORD; j++)	// For top & bottom border scanlines,
				{
					for(i=JT_BORD; i<imax; i+=2)
						// and within the horiz. borders,
					{					// copy nearest scanline within border
						buf.put(i,j, buf.get(i,JT_BORD));	// top 
						buf.put(i,jmax+j,buf.get(i,jmax-1));// bottom
					}
				}
				for(j=JT_BORD; j<jmax; j++)
					// for left & right side borders,
				{
					for(i=0; i<JT_BORD; i++)
					{						// copy nearest column within border
						buf.put(i,j,buf.get(JT_BORD,j));// left side
						buf.put(i+imax,j,buf.get(imax-1,j));// right
					}
				}
				for(j=0; j<JT_BORD; j++)	// Now fill in all the corner values;
				{
					for(i=0; i<JT_BORD; i++)
					{
						buf.put(i,j,buf.get(JT_BORD,JT_BORD));		// lower left,
						buf.put(i+imax,j,buf.get(imax-1,JT_BORD));	// lower right,
						buf.put(i,j+jmax,buf.get(JT_BORD,jmax-1));	// upper left,
						buf.put(i+imax,j+jmax,buf.get(imax-1,jmax-1));// upper right.
					}
				}
			}
			//-----END boundary-value setting-----------
			DPF("%f\n", (float)err2max);
		}							// go on to the next iteration.

		if(err2max < maxErr*maxErr)// Are we below the error maximum?
		{
			break;				// Yes. Stop.  Answer is in 'oldB' buffer.
		}						// No; keep going.
	}
	if(bound==0) DPF("VonNeumann--"); else DPF("Dirichlet(DEFAULT)--");
	if(k==maxCount)				// did we stop before error < errMax?
	{							// Yes; report it.
		DPF("Raw2D::Poisson2D_SOR() stopped at err=%f, not err=%f!\n",(float)sqrt(err2max),(float)maxErr);
	}
	else
	{
		DPF("Raw2D::Poisson2D_SOR() stopped at err=%f after %d counts\n",(float)sqrt(err2max),k);
	}
	///HIDE BORDER/////////////
	sizer(pLap);				// Copy result to output image, without border;
	for(j=JT_BORD; j<jmax; j++)
	{
		for(i=JT_BORD; i<imax; i++)
		{
			put(i-JT_BORD,j-JT_BORD,buf.get(i,j));
		}
	}

	//END HIDE BORDER////////////////
	//diagnostic only: Show border://///////////////////
	//	wipecopy(&(buf[oldB]));
	// end diagnostic 'show border'/////////////////////
	return(TRUE);

#undef JT_BORD					// locals only!
#undef JT_FRAC
}

BOOL Raw2D::poisson2D_SOR(Raw2D *pLap, int maxCount, int bound, double maxErr)
//------------------------------------------------------------------------------
// Poisson Solver using Successive Over-Relaxation method--converges O(N), while
// 'naive' and GaussSeidel converges O(N^2). It is very similar to the solver
// above (SORfix), but the correction for Laplacian error is exaggerated by
// an amount given by a Chebyshev polynomial. Rather than use a fixed value <2,
// (e.g. 1.99) where error can grow for the first N iterations before converging
// using the Chebyshev acceleration guarantees monotonically shrinking error.
//
//If bound%2=1, (DEFAULT):
// Use Dirichlet boundary conditions: Presume the image is surrounded by an
// unbounded field of zero-valued pixels. 
// If bound%2=0 : (Not recommended)
// Uses VonNeumann boundary conditions: assumes gradients across across the
// image boundaries are zero, but unchanged in tangent direction).  
//
//		NOTES: Dirichlet conditions give you a unique solution; VonNeumann
// boundary conditions have an unknown constant offset. Poisson solutions can 
// have persistent large error near boundaries if the laplacian source image was
//  computed with unknown boundary conditions, or was a section of a larger 
//  Laplacian image. You can improve results by enclosing the Laplacian in a 
// border of zero-valued pixels, and use either VonNeumann or Dirichlet..
// 
//	(fcnGroup::fcnCurv1() finds 2nd derivatives w/ VonNeumann boundaries;)
// **(fcnGroup::fcnCurv2() finds 2nd derivatives w/Dirichlet boundaries.)**
//
// 06/10/2004 Created--J. Tumblin
{
	int i,j,k,imax,jmax;
	int isOdd,chk;					// Checkerboard generator.
	Raw2D buf;						// single-buffer for computing output;
	PIXTYPE lapNow,lapErr;			// current Laplacian, & how far wrong it is.
	PIXTYPE lapAvg,lapWanted;		// current local average; desired laplacian
	double err2sum;					// squared error sum in output's Laplacian
	double rjac,rjac2;				// Radius of Convergence for the Jacobi method
	// which is known analytically for the Poisson
	// solver (pg 869, eqn 19.5.24, Num. Recipes)
	PIXTYPE chebW;					// Chebyshev acceleration weight.
#define JT_BORD   1				// Width of zero-valued border around image.

	if(pLap==NULL)				// input does not exist? 
	{
		BEEP;					// complain; exit.
		DPF("Raw2D::Poisson2D_SORcheby() given NULL input!\n");
		return(FALSE);
	}
	if(pLap->getYsize()<=0 || pLap->getXsize()<=0)
	{
		BEEP;					// complain; exit.
		DPF("Raw2D::Poisson2D_SORcheby() given zero-sized input!\n");
		return(FALSE);
	}
	// Make a single output image buffer that
	// can hold pLap plus a border of 'JT_BORD'
	// pixels wide on all 4 sides:
	buf.sizer(2*JT_BORD + pLap->getXsize(),2*JT_BORD + pLap->getYsize());
	imax = JT_BORD + pLap->getXsize();// Find largest address within the border)
	jmax = JT_BORD + pLap->getYsize();
	if(bound==1)				// Initialize the buffers:
	{							// for (default) Dirichlet boundary conditions,
		buf.putConst((PIXTYPE)0.0);	// set image & surrounds to zero.
	}
	else
	{							// for (not reccommended)VonNeumann conditions,
		buf.putConst((PIXTYPE)0.5);	// grow darker & lighter from the
	}								// mid-level gray value of 0.5
	// To compute the Chebyshev acceleration weight chebW we need to know the
	// radius of convergence for Jacobi iterations for our solver. Rectangular
	// boundaries (like ours) with either Dirichlet or VonNeumann conditions
	// used in solution to a Poisson equation have a radius of convergence given
	// by an analytical expression (eqn. 19.5.24, pg. 869, Numerical Recipes)
	rjac = (cos(M_PI/imax) + cos(M_PI/jmax)) / 2.0;
	rjac2 = rjac*rjac;
	// (Note rjac is asymptotically 1.0 as imax and jmax grows.)
	chebW = (PIXTYPE)1.0;			// for 1st pass, no acceleration.
	for(k=0; k<maxCount; k++)	// Iteratively modify the output buffer so the
	{							// Laplacian at each pix matches pLap values:
		err2sum = 0.0;				// clear the max-error-finder for this pass.
		for(isOdd=0; isOdd<2; isOdd++)	// for each checkerboard (red-black)
		{								// half-sweep of the image:		
			chk = isOdd;				// (starting pixel for current scanline)
			for(j=JT_BORD; j<jmax; j++)// update every scanline,
			{
				for(i=JT_BORD+chk; i<imax; i+=2)	// but only half the pixels,
				{	
					lapWanted = pLap->get(i-JT_BORD, j-JT_BORD); 
					lapAvg = (buf.get(i+1,j) + buf.get(i-1,j) + // sum NSEW neighbors
						buf.get(i,j+1) + buf.get(i,j-1));
					lapNow = lapAvg - (PIXTYPE)4.0*buf.get(i,j);//current Laplacian
					lapErr = lapWanted - lapNow;
					// in-place correction of 1/4 of the Laplacian error
					// but exaggerate the correction by chebW factor,
					// which grows as large as 2.0.
					buf.put(i,j, buf.get(i,j) - (PIXTYPE)0.25 *chebW * lapErr);
					// (compute the error).
					err2sum += (double)(lapErr*lapErr);	// sum-of-square err
				}// end of one scanline,
				chk = (chk+1)%2;	// toggle 'chk' on each scanline
				if(k==0 && isOdd==0)// Initial pass on chekerboard's 2nd half?
				{					// calc the Chebyshev polynomial:
					chebW = (PIXTYPE)( 1.0 / (1.0 - 0.5*rjac2));
				}					// (chebW should be ~2.0)
				else				// No? asymptotically go to optimal value:
				{
					chebW = (PIXTYPE)(1.0 / (1.0 - 0.25*rjac2*chebW));
				}					// (chebW should be ~2.0)
			}// end of one half-sweep,
			//if(j%100==0) DPF(".");// print little dots every 100 scanlines...
			/////Set Boundary conditions---------------------------------------
			// bound%2==1 means Dirichlet Boundary conditions:
			// We already have them! Border pixels initialized to zero & unchanged.
			if(bound%2==0)		// VonNeumann Boundary conditions----------
			{						// update the border pixels so that their 
				// values match the adjacent pixel:
				for(j=0; j<JT_BORD; j++)	// For top & bottom border scanlines,
				{
					for(i=JT_BORD; i<imax; i+=2)
						// and within the horiz. borders,
					{					// copy nearest scanline within border
						buf.put(i,j, buf.get(i,JT_BORD));	// top 
						buf.put(i,jmax+j,buf.get(i,jmax-1));// bottom
					}
				}
				for(j=JT_BORD; j<jmax; j++)
					// for left & right side borders,
				{
					for(i=0; i<JT_BORD; i++)
					{						// copy nearest column within border
						buf.put(i,j,buf.get(JT_BORD,j));// left side
						buf.put(i+imax,j,buf.get(imax-1,j));// right
					}
				}
				for(j=0; j<JT_BORD; j++)	// Now fill in all the corner values;
				{
					for(i=0; i<JT_BORD; i++)
					{
						buf.put(i,j,buf.get(JT_BORD,JT_BORD));		// lower left,
						buf.put(i+imax,j,buf.get(imax-1,JT_BORD));	// lower right,
						buf.put(i,j+jmax,buf.get(JT_BORD,jmax-1));	// upper left,
						buf.put(i+imax,j+jmax,buf.get(imax-1,jmax-1));// upper right.
					}
				}
			}
			//-----END boundary-value setting-----------
			//DPF("%f\n", (float)err2sum);
		}							// go on to the next iteration.

		if(err2sum < maxErr*maxErr)// Are we below the error maximum?
		{
			break;				// Yes. Stop.  Answer is in 'oldB' buffer.
		}						// No; keep going.
	}
	if(bound==0) DPF("VonNeumann--"); else DPF("Dirichlet(DEFAULT)--");
	if(k==maxCount)				// did we stop before error < errMax?
	{							// Yes; report it.
		DPF("Raw2D::Poisson2D_SOR() stopped at err=%f, not err=%f!\n",(float)sqrt(err2sum),(float)maxErr);
	}
	else
	{
		DPF("Raw2D::Poisson2D_SOR() stopped at err=%f after %d counts\n",(float)sqrt(err2sum),k);
	}
	///HIDE BORDER/////////////
	sizer(pLap);				// Copy result to output image, without border;
	for(j=JT_BORD; j<jmax; j++)
	{
		for(i=JT_BORD; i<imax; i++)
		{
			put(i-JT_BORD,j-JT_BORD,buf.get(i,j));
		}
	}

	//END HIDE BORDER////////////////
	//diagnostic only: Show border://///////////////////
	//	wipecopy(&(buf[oldB]));
	// end diagnostic 'show border'/////////////////////
	return(TRUE);

#undef JT_BORD					// locals only!
#undef JT_FRAC
}