#include <math.h>	/* sqrt() for shock filter equation */
#include <iostream>
#include "ppm.h"
using namespace std;

namespace math144_ppm
{

/* Derivative approximations for h = 1 */
#define X_FWD (old.data[i+1][j] - old.data[i][j])
#define Y_FWD (1) // PROGRAM HERE, replace (1) with the forward difference
#define X_BCK (1) // PROGRAM HERE, replace (1) with the backwards difference
#define Y_BCK (1) // PROGRAM HERE, replace (1) with the backwards difference
#define X_CTR (1) // PROGRAM HERE, replace (1) with the centered difference
#define Y_CTR (1) // PROGRAM HERE, replace (1) with the centered difference
#define X_2ND (1) // PROGRAM HERE, replace (1) with the second derivative
#define Y_2ND (1) // PROGRAM HERE, replace (1) with the second derivative

double m_squared(double a, double b)
{
	/*
	 * PROGRAM HERE
	 * This function should return the square of the "m" function, as
	 * defined for use in the shock filter equation.  See your notes for
	 * details.
	 */
	return 0;
}

void heat(image &img, double dt, unsigned int num_iter)
{
	/*
	 * The heat equation has been programmed as an example for the student.
	 * Note it uses X_2ND and Y_2ND, which must be defined above by the
	 * student.
	 */

	/* Make a copy of the image */
	image old(img);

	/* Run for each iteration */
	for (int count = 0; count < num_iter; count++)
	{
		/* Run the heat equation on the image */
		for (int i = 1; i < img.width - 1; i++) {
			for (int j = 1; j < img.height - 1; j++) {
				img.data[i][j] = old.data[i][j] + dt * (X_2ND + Y_2ND);
			}
		}

		/* Update for the next iteration */
		old = img;
	}
}

void lvlset(image &img, double dt, unsigned int num_iter)
{
	/* PROGRAM THE LEVEL SET EQUATION HERE */
}

void stdylvlset(image &img, double dt, double alpha, double tolerance)
{
	/*
	 * PROGRAM HERE
	 * The steady state level set equation uses the original image (in
	 * your notes, it's called "f(x,y)").  So, in addition to making a copy
	 * of the image for the previous time step, you'll need a copy of the
	 * original image to use as f(x,y).  You should do this now, right
	 * after this comment.
	 */


	/*
	 * The steady state level set equation is different from the other
	 * equations.  Instead of running for a given time, it runs as many
	 * iterations as needed until the image reaches a "steady state,"
	 * meaning that the pixels stop changing.  We program this by seeing
	 * if every pixel is within a tolerance of the previous time step.
	 * If so, then the image is "steady enough" to stop.
	 *
	 * We'll keep track of this by creating a variable called "above_toler."
	 * If above_toler = 1, then we're still above the tolerance and need to
	 *                          keep running the PDE.
	 * If above_toler = 0, then we're below the tolerance and can stop.
	 */
	bool above_toler;
	int iter = 0;
	do {
		/*
		 * Out of curiosity, keep track of how many iterations it takes
		 * to reach steady state.
		 */
		iter++;

		/*
		 * Start off by assuming all the pixels are below the tolerance.
		 * In your program below, you'll change the value of this if
		 * a pixel is above the tolerance.
		 */
		above_toler = 0;

		/*
		 * PROGRAM HERE
		 * Go through each pixel and run the steady state level set
		 * equation.  After you run the PDE for a particular pixel,
		 * check whether it's above the tolerance.  If it is, then
		 * set above_toler = 1 so that we know we haven't reached
		 * steady state yet.  You should do this now, right after
		 * this comment.
		 */


		/*
		 * If we got through all the pixels and each was below the
		 * tolerance, then we've reached steady state.  So, we can
		 * break out of the loop and stop.
		 */
		if (above_toler == 0) {
			break;
		}
	} while (above_toler == 1);

	/* Print the number of iterations it took to reach steady state. */
	cout << "It took " << iter << " iterations to reach steady state."
	     << endl;
}

void shock(image &img, double dt, unsigned int num_iter)
{
	/* PROGRAM THE SHOCK FILTER EQUATION HERE */
}

} // math144_ppm