/**------------------------------------------------------------------------
 * @file
 *
 * Copyright 2001 by the University of Toronto (UT).
 *
 *                       All Rights Reserved
 * 
 * Permission to copy and modify this software and its documentation 
 * only for internal use in your organization is hereby granted, 
 * provided that this notice is retained thereon and on all copies.  
 * UT makes no representations as to the suitability and operability 
 * of this software for any purpose.  It is provided "as is" without 
 * express or implied warranty.
 *
 * UT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING
 * ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT
 * SHALL UT BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
 * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA
 * OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
 * TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
 * PERFORMANCE OF THIS SOFTWARE.
 *
 * No other rights, including for example, the right to redistribute 
 * this software and its documentation or the right to prepare 
 * derivative works, are granted unless specifically provided in a 
 * separate license agreement.
 *
 * Author: Diego Macrini
 * 
 * $Log: ImgPreprocess.cpp,v $
 *
 *-----------------------------------------------------------------------*/
#include <stack>
#include <functional>
#include <DirWalker.h>
#include <tools/connected.h>

#include "ImgPreprocess.h"

#include <vector>
#include <field.h>
#include <CImg/CImg.h>
#include <ImageProcessing.h>

#ifdef WIN32
#define drand48() (double(rand()) / RAND_MAX)
#define srand48(A) (srand((unsigned int) (A)))
#else
#define atan2 std::atan2
#endif

#define MAX_NUM_B_N 1
#define BLACK 0
#define WHITE 255

#ifndef PI
#define PI 3.141592653589793238462643
#endif

int g_nInitShockCount = 0;

typedef std::vector<int> IntVector;
typedef std::vector< std::vector<float> > EdgelVector;

using namespace std;
using namespace dml;

/*!
	@brief Find the edge pixels or elements (ie, edgels) in an image
	
	@param img Image to be processed
	@param mag is a return variable. It is given as a reference to CImg<float> image with 
		size equal to img. It returns the gradient magnitudes.
	@param phase is a return variable. It is given as a reference to CImg<float> image with 
		size equal to img. It returns the gradient phase (ie, edge angle at pixel x,y).
	@param x0 is a return variable. It is an int with the x-coordinate of
		the first edge point found. It is zero if no edge point is found.
	@param y0 is a return variable. It is an int with the y-coordinate of
		the first edge point found. It is zero if no edge point is found.
	@return number of edgels found.
*/
int FindEdgels(const cimg_library::CImg<float>& img, cimg_library::CImg<float>& mag, cimg_library::CImg<float>& phase, int& x0, int& y0)
{
	int x,y,numPts;
	long g1, g2, g3, g4;
	float xx, yy, xc, yc, g, gi1, gi2, c, ci1, ci2;

	cimg_library::CImgl<> grad = img.get_gradientXY(3);
	const cimg_library::CImg<float>& dx = grad[0];
	const cimg_library::CImg<float>& dy = grad[1];
	
	// Init first edge point coords and num edgels to zero 
	x0 = y0 = 0;
	numPts = 0;
	
	// Scale floating point magnitudes and angles to 8 bits
	//const float ORI_SCALE = 40.0f;
	//const float MAG_SCALE = 20.0f;

	// Perform non-maximum suppression
	for (y = 1; y < mag.dimy() - 1; y++)
	{
		for (x = 1; x < mag.dimx() - 1; x++)
		{
			mag(x,y) = 0;

			// Treat the x and y derivatives as components of a vector
			xc = dx(x,y);
			yc = dy(x,y);
			
			// Speed-up: ignore pixels with near-zero gradient magnitud
			if (fabs(xc)<0.01 && fabs(yc)<0.01) continue;
	
			// Follow the gradient direction, as indicated by the direction of
			// the vector (xc, yc) to find pair of neighbouring pixels.
			if (fabs(yc) > fabs(xc))
			{
				// The Y component is biggest, so gradient direction is basically UP/DOWN
				xx = (float) (fabs(xc) / fabs(yc));
				yy = 1.0;
	
				g2 = mag.offset(x,y-1);
				g4 = mag.offset(x,y+1);
				if (xc * yc > 0.0) {
					g3 = mag.offset(x+1,y+1);
					g1 = mag.offset(x-1,y-1);
				} 
				else {
					g3 = mag.offset(x-1,y+1);
					g1 = mag.offset(x+1,y-1);
				}
			} 
			else {
				// The X component is biggest, so gradient direction is basically LEFT/RIGHT
				xx = (float) (fabs(yc) / fabs(xc));
				yy = 1.0;
	
				g2 = mag.offset(x+1,y);
				g4 = mag.offset(x-1,y);
				if (xc * yc > 0.0) {
					g3 = mag.offset(x-1,y-1);
					g1 = mag.offset(x+1,y+1);
				}
				else {
					g1 = mag.offset(x+1,y-1);
					g3 = mag.offset(x-1,y+1);
				}
			}

			// Compute gradient magnitud at pixel y,x
			g  = NORM2(xc, yc);
			c = img(x, y);

			// Compute the interpolated value of the gradient magnitude
			gi1 = xx * NORM2(dx[g1],dy[g1]) + (yy - xx) * NORM2(dx[g2],dy[g2]);
			gi2 = xx * NORM2(dx[g3],dy[g3]) + (yy - xx) * NORM2(dx[g4],dy[g4]);

			// Compute the interpolated value of the color intensity
			ci1 = xx * img[g1] + (yy - xx) * img[g2];
			ci2 = xx * img[g3] + (yy - xx) * img[g4];
			
			// Choose local maximum and supress non-max. If there is a tie,
			// choose the pixel with darkest color intensity (black=0, white=255)
			if ((g > gi1 && g > gi2) || 
			    (g == gi1 && g > gi2 && c <= ci1) || 
			    (g > gi1 && g == gi2 && c <= ci2))
			{
				/*if (g * MAG_SCALE <= 255)
					mag(x,y) = (unsigned char)(g * MAG_SCALE);
				else
					mag(x,y) = 255;

				phase(x,y) = (char)(std::atan2(yc, xc) * ORI_SCALE);
				*/
		
				mag(x,y) = g;
				phase(x,y) = atan2(yc, xc);
				
				// Check if it is the first edge point found and update numPts
				if (numPts++ == 0) {
					x0 = x;
					y0 = y;
				}
			} 
			else {
				mag(x,y) = 0;
				phase(x,y) = 0;
			}
		}
	}
	
	/*CImgDisplay main_disp(mag.get_normalize(0,255),"Magnitud");
	
	while (!main_disp.closed)
		main_disp.wait();*/
		
	//mag.save_png("mag.pgm");
	/*dx.save_ascii("loaddx.m");
	dy.save_ascii("loaddy.m");
	
	CImgDisplay main_disp3(dx.get_normalize(0,255),"Matute's window 3");
	CImgDisplay main_disp4(dy.get_normalize(0,255),"Matute's window 4");
	
	while (!main_disp3.closed && !main_disp4.closed) {
		main_disp3.wait();
	}*/
	
	return numPts;
}

void MakeNextEdgel(int& x, int& y, cimg_library::CImg<float>& mag)
{
	const int cx = x, cy = y;
	int i, inc;
	int maxInc = MIN(MIN(x, y), MIN(mag.dimx() - x - 2, mag.dimy() - y - 2));
	
	// Mark current point as visited making its magnitud negative
	mag(cx, cy) = -mag(cx, cy);
	
	for (inc = 1; inc <= 7 && inc <= maxInc; inc++)
	{	
		for (i = 0; i <= inc; i++)
		{	
			if (mag(x = cx - i, y = cy - inc) > 0 || mag(x = cx - i, y = cy + inc) > 0)
					return;
			
			if (mag(x = cx + i, y = cy - inc) > 0 || mag(x = cx + i, y = cy + inc) > 0)
					return;
					
			if (mag(x = cx - inc, y = cy - i) > 0 || mag(x = cx + inc, y = cy - i) > 0)
					return;
					
			if (mag(x = cx - inc, y = cy + i) > 0 || mag(x = cx + inc, y = cy + i) > 0)
					return;
		}
	}
		
	x = 0;
	y = 0;
}

EdgelVector FindAndSortEdgels(const cimg_library::CImg<float>& img)
{
	cimg_library::CImg<float> mag(img.dimx(), img.dimy());
	cimg_library::CImg<float> phase(img.dimx(), img.dimy());
	int x = 0, y = 0, numPts, i;
	
	numPts = FindEdgels(img, mag, phase, x, y);
	
	std::vector<float> nilEdgePt(4,0);
	EdgelVector edgels(numPts, nilEdgePt);
	
	for (i = 0; i < numPts && x != 0; i++)
	{
		edgels[i][0] = x;
		edgels[i][1] = y;
		edgels[i][2] = mag(x, y);
		edgels[i][3] = phase(x, y);
		MakeNextEdgel(x, y, mag);
	}
	
	// Update numPts to represent the number of (valid) connected edge pts
	if (i < numPts)
		edgels.resize(i);
	
	// Number the edges in the magnitud matrix for testing
	//numPts = i;
	//for (i = 0; i < numPts; i++)
	//	mag((int)edgels[i][0], (int)edgels[i][1]) = i;
	//mag.save_ascii("dogmag.raw");
	
	return edgels;
}

cimg_library::CImg<float> RGBToGrayscale(const cimg_library::CImg<float>& img)
{
	return img.get_channel(0) * 0.2989 
		+ img.get_channel(1) * 0.5870 
		+ img.get_channel(2) * 0.1140;
}

/*!
	@brief
				 
*/
//SmartPtr<void> AddBumpsAndNotches(const char* szPPMFileName, const SGRECINFO& info, std::string& newFileName)
/*bool AddBumpsAndNotches(const char* szPPMFileName, const SGRECINFO& info, IMAGEINFO* pImgInfo)
{
//	using namespace cimg_library;
	
	CImg<float>* pImg = new CImg<float>(szPPMFileName);
	
	pImgInfo->ptrImg = pImg;		// for the sg code which deletes it
	CImg<float>& img = *pImg;		// for using it right here
	
	// If color image, we only keep the "red" chanel
	if (img.dimv() > 1)
		img = img.channel(0);
		//img = RGBToGrayscale(img);
	
	// Removes "white holes" and "black satellites" from the image
	RemoveNonMaxComponents(img);
	
	// Find points along the shape boundary
	EdgelVector edgels = FindAndSortEdgels(img);
	
	// If there is no noise to add we just return	
	if (info.m_nBumpSize == 0 && info.m_nNotchSize == 0)
	{
		pImgInfo->strFileName = szPPMFileName;
		pImgInfo->boundaryPts.ReSize(edgels.size());
		
		for (int i = 0; i < edgels.size(); i++)
			pImgInfo->boundaryPts[i].Set(edgels[i][0], edgels[i][1]);
			
		return true;
	}
	
	// Prepare variables, eg, a copy of the original img
	//CImg<float>* pNoisyImg = new CImg<float>(img);
	//pImgInfo->ptrImg = pNoisyImg;		// for the sg code which deletes it
	CImg<float>& noisyImg = *pImg;
	
	int r;
	float colorBump = 0, colorNotch = 255, a;
	double x0, y0, x1, y1, x2, y2;
	int i, x, y, maxRad = MAX(info.m_nBumpSize,info.m_nNotchSize), spacing = 30;
	int delta = 2 * maxRad + spacing;
	int maxNumBN = 10;//edgels.size() / delta * MAX_NUM_B_N;
	int minNumBN = 3;//int(maxNumBN * 0.05);
	
	// Randomly choose some edge pts along the boundary
	IntVector pts(maxNumBN);
	IntVector::iterator start, end, it;
	
	// Make things really random
	srand48((long)time(NULL));

	for (i = 0; i < maxNumBN; i++)
		pts[i] = i * delta + maxRad;
		
	start = pts.begin();
	end = pts.end();
	
	// Shuffle the elements in a random order
   	std::random_shuffle(start, end);
	
	// Choose how many bumps & notches we want to add
	int numBumps = 0;
	int numNotches = 0;
	int addedBumps = 0, addedNotches = 0, numDiffs = 0;
	
	// Set the new beggining of the vector
	it = pts.begin();
	
	// Add bumps
	if (info.m_nBumpSize > 0)
	{
		numBumps = int(drand48() * (maxNumBN - minNumBN)) + minNumBN;
		
		for (; it != end && addedBumps < numBumps; it++, addedBumps++) 
		{
			i = *it;
			x = (int)edgels[i][0];
			y = (int)edgels[i][1];
			a = edgels[i][3];
			r = info.m_nBumpSize; //(int) (1 + drand48() * (info.m_nBumpSize - 1));
			numDiffs += r * 2 + 1;
			
			if (drand48() >= .5)
				noisyImg.draw_circle(x, y, r, &colorBump);
			else {
				x0 =  1.5 * r; y0 = 0;
				x1 = -0.5 * r; y1 = 2.0 * r / 1.5;
				x2 =  x1     ; y2 = -y1;
			
				noisyImg.draw_triangle(
					(int)(x0 * cos(a) + y0 * sin(a)) + x, (int)(x0 * sin(a) - y0 * cos(a)) + y,
					(int)(x1 * cos(a) + y1 * sin(a)) + x, (int)(x1 * sin(a) - y1 * cos(a)) + y,
					(int)(x2 * cos(a) + y2 * sin(a)) + x, (int)(x2 * sin(a) - y2 * cos(a)) + y,
					&colorBump);
			}
		}
	}
	
	// Add notches
	if (info.m_nNotchSize > 0)
	{
		numNotches = int(drand48() * (maxNumBN - minNumBN - numBumps)) + minNumBN;

		for (; it != end && addedNotches < numNotches; it++, addedNotches++) 
		{
			i = *it;
			x = (int)edgels[i][0];
			y = (int)edgels[i][1];
			a = edgels[i][3] + cimg::PI;
			r = info.m_nNotchSize; //(int) (1 + drand48() * (info.m_nNotchSize - 1));
			numDiffs += r * 2 + 1;
			
			if (drand48() >= .5)
				noisyImg.draw_circle(x, y, r, &colorNotch);
			else {
				x0 =  1.5 * r; y0 = 0;
				x1 = -0.5 * r; y1 = 2.0 * r / 1.5;
				x2 =  x1     ; y2 = -y1;
			
				noisyImg.draw_triangle(
					(int)(x0 * cos(a) + y0 * sin(a)) + x, (int)(x0 * sin(a) - y0 * cos(a)) + y,
					(int)(x1 * cos(a) + y1 * sin(a)) + x, (int)(x1 * sin(a) - y1 * cos(a)) + y,
					(int)(x2 * cos(a) + y2 * sin(a)) + x, (int)(x2 * sin(a) - y2 * cos(a)) + y,
					&colorNotch);
			}	
		}
	}
	
	// Remove components that got disconnected after adding noise
	RemoveNonMaxComponents(noisyImg);
	
	// Comp. the level of noise introduced along the boundary
	double dNoiseLevel = (100.0 * numDiffs) / edgels.size();
	
	// Construct a new file name that encodes noise level
	char szSubfix[1024];
	sprintf(szSubfix, "n%.0f.pgm", dNoiseLevel);
	
	DirWalker dw;
	dw.ParseFileName(szPPMFileName);

	std::string newName = dw.path + dw.name + szSubfix; //+ dw.ext;
	cerr << newName << " noise level = " << dNoiseLevel << endl;
	
	// Find the new countour with bumps and notches
	edgels = FindAndSortEdgels(noisyImg);
	
	pImgInfo->boundaryPts.ReSize(edgels.size());
		
	for (int i = 0; i < edgels.size(); i++)
		pImgInfo->boundaryPts[i].Set(edgels[i][0], edgels[i][1]);
			
	
	pImgInfo->strFileName = newName.c_str();
	
	// Save the noisy image
	noisyImg.save_png(newName.c_str());
	
	//CImgDisplay main_disp1(img,"Matute's window 1");
	//CImgDisplay main_disp2(noisyImg,"Matute's window 2");
	//CImgDisplay main_disp3(phase.normalize(0,255),"Matute's window 3");
	
	//while (!main_disp2.closed) {
	//	main_disp2.wait();
	//}
	
	return true;
}*/

void DrawEllipse(cimg_library::CImg<float>& img, int x, int y, float a, int r, float color)
{
	img.draw_circle(x, y, r, &color);
}

void DrawTriangle(cimg_library::CImg<float>& img, int x, int y, float a, int r, float color)
{
	double x0, y0, x1, y1, x2, y2;
	
	x0 =  1.5 * r; y0 = 0;
	x1 = -0.5 * r; y1 = 2.0 * r / 1.5;
	x2 =  x1     ; y2 = -y1;

	img.draw_triangle(
		(int)(x0 * cos(a) + y0 * sin(a)) + x, (int)(x0 * sin(a) - y0 * cos(a)) + y,
		(int)(x1 * cos(a) + y1 * sin(a)) + x, (int)(x1 * sin(a) - y1 * cos(a)) + y,
		(int)(x2 * cos(a) + y2 * sin(a)) + x, (int)(x2 * sin(a) - y2 * cos(a)) + y,
		&color);
}

int CountShockNodes(cimg_library::CImg<float>& img, const SGRECINFO& info)
{
	IMAGEINFO imgInfo;
	SGCP sgparams;
	ShockGraph sg;
	
	//imgInfo.ptrImg = new cimg_library::CImg<float>(img);
	imgInfo.pField = CopyImageToField(img);
	imgInfo.strFileName = "./Obj/Test10.ppm";
	
	sgparams = info.sgparams;
	sgparams.tsTau1 = 0; // no simplification
	sgparams.tsTau2 = 10;
	
	sg.ComputeFromPPMFile3(imgInfo, sgparams);
	
	
	//instead of returning the node count we are
	//now returning the edge count, because what
	//we really want is to have more branches
	//return sg.GetNumberOfBranches();
	return sg.GetNodeCount();
}

double AddNoiseSpot(cimg_library::CImg<float>& img, const SGRECINFO& info, int x, int y, float a, int size, int type, bool bTrySmaller)
{
	if (size < 2)
		return 0;
	
	cimg_library::CImg<float> noisyImg(img);
		
	if (type == 0) 
	{
		if (drand48() >= .5)
			DrawEllipse(noisyImg, x, y, a, size, BLACK);
		else
			DrawTriangle(noisyImg, x, y, a, size, BLACK);
	}
	else
	{
		if (drand48() >= .5)
			DrawEllipse(noisyImg, x, y, a + PI, size, WHITE);
		else
			DrawTriangle(noisyImg, x, y, a + PI, size, WHITE);
	}
	
	// Remove components that got disconnected after adding noise
	RemoveNonMaxComponents(noisyImg);
		
	double totalDiff = 0;
	cimg_mapXY(img,x0,y0) {
			if (img(x0,y0) != noisyImg(x0,y0))
				totalDiff++;
	}
	
	if (totalDiff > PI * size * size)
	{
		// The noise spot is too big for its context
		if (bTrySmaller)
			return AddNoiseSpot(img, info, x, y, a, size - 1, type, true);
		else
			return 0;
	}
	
	int nShockCount = CountShockNodes(noisyImg, info);
// 	DBG_MSG1("NEW Node Count", nShockCount);
	
	if (nShockCount <= g_nInitShockCount)
		return 0;
	else {
// 		cerr << "\n\n\n\n\nYES!!!!!!\n\n\n\n\n\n" << endl;
		g_nInitShockCount = nShockCount;
	}
	
	img = noisyImg;
	
	return size * size * PI / 2;
}



void GetNextEdgel(int& x, int& y, int& dir, 
	cimg_library::CImg<float>& img, int nDimx, int nDimy)
{
	for( int k = 0; k < 8; ++k)
	{
		switch( (dir + k) % 8)
		{
			case 0:	//try at (x+1, y)
				if( x < nDimx && img(x+1,y) == 0)
				{
					++x;
					dir = 5;
					return;
				}
			case 1:	//try at (x+1, y+1)
				if( x < nDimx && y < nDimy && img(x+1,y+1) == 0)
				{
					++x;
					++y;
					dir = 6;
					return;
				}
			case 2:	//try at (x, y+1)
				if( y < nDimy && img(x,y+1) == 0)
				{
					++y;
					dir = 7;
					return;
				}
			case 3:	//try at (x-1, y+1)
				if( x > 0 && y < nDimy && img(x-1,y+1) == 0)
				{
					--x;
					++y;
					dir = 0;
					return;
				}
			case 4:	//try at (x-1, y)
				if( x > 0 && img(x-1,y) == 0)
				{
					--x;
					dir = 1;
					return;
				}
			case 5:	//try at (x-1, y-1)
				if( x > 0 && y > 0 && img(x-1,y-1) == 0)
				{
					--x;
					--y;
					dir = 2;
					return;
				}
			case 6:	//try at (x, y-1)
				if( y > 0 && img(x,y-1) == 0)
				{
					--y;
					dir = 3;
					return;
				}
			case 7:	//try at (x+1, y-1)
				if( x < nDimx && y > 0 && img(x+1,y-1) == 0)
				{
					++x;
					--y;
					dir = 4;
					return;
				}
		}
	}
}



/*
	Finds the direction we should walk in given
	the first boundary point we use to create the boundary
*/
int FirstDirectionBoundary(int& x, int& y, cimg_library::CImg<float>& img)
{
	//check whether this is a boundary pixel
	if(img(x,y) != 0) return -1;
	
	//Find first neighbor outside
	//Check for border pixels
	if(x+1 >= img.dimx()) return 0;
	if(y+1 >= img.dimy()) return 2;
	if(x == 0) return 4;
	if(y == 0) return 6;
	
	if(img(x+1, y) != 0)	return 0;
	if(img(x+1, y+1) != 0)	return 1;
	if(img(x, y+1) != 0)	return 2;
	if(img(x-1, y+1) != 0)	return 3;
	if(img(x-1, y) != 0)	return 4;
	if(img(x-1, y-1) != 0)	return 5;
	if(img(x, y-1) != 0)	return 6;
	if(img(x+1, y-1) != 0)	return 7;
	
	return -1; //no neighbors are outside this point
}



bool AddBumpsAndNotches(const char* szPPMFileName, const SGRECINFO& info, IMAGEINFO* pImgInfo)
{
	/*cimg_library::CImg<float>* pImg = new cimg_library::CImg<float>(szPPMFileName);

	pImgInfo->ptrImg = pImg;		// for the sg code which deletes it
	cimg_library::CImg<float>& img = *pImg;		// for using right here*/

	cimg_library::CImg<float> img(szPPMFileName);

	// If color image, we only keep the "red" chanel
	if (img.dimv() > 1)
		img = img.channel(0);
		//img = RGBToGrayscale(img);

	// Removes "white holes" and "black satellites" from the image
	RemoveNonMaxComponents(img);

	// Find points along the shape boundary
	EdgelVector edgels = FindAndSortEdgels(img);

// //BEGIN: for now...
// 
// // 	int xInit = 0, yInit = 0;
// // 	
// // 	//Find initial starting point
// // 	cimg_mapXY(img,x,y)
// // 	{
// // 		if(img(x,y) == 0) // image foreground point
// // 		{
// // 			xInit = x;
// // 			yInit = y;
// // 			break;
// // 		}
// // 	}
// // 	
// // 	//Find the direction we should walk in first
// // 	int dir = FirstDirectionBoundary(xInit, yInit, img);
// // 	
// // 	if(dir < 0)
// // 		std::cerr << "We did not start with a boundary pixel\n";
// 		
// 	
// 	
// // 	CImg<float> mag(img.dimx(), img.dimy());
// // 	CImg<float> bndry(img.dimx(), img.dimy());
// // 	
// // 	cimg_mapXY(bndry, x, y)
// // 	{
// // 		bndry(x,y) = 0;
// // 	}
// // 	
// // 	CImg<float> phase(img.dimx(), img.dimy());
// // 	int x = 0, y = 0, numPts, i;
// // 	
// // 	numPts = FindEdgels(img, mag, phase, x, y);
// // 	//get some extra space...
// // 	numPts *= 3;
// // 	
// // 	x = xInit;
// // 	y = yInit;
// // 	i = 0;
// // 	
// // 	std::vector<float> nilEdgePt(4,0);
// // 	EdgelVector edgels(numPts, nilEdgePt);
// // 	
// // 	do
// // 	{
// // 		edgels[i][0] = x;
// // 		edgels[i][1] = y;
// // 		edgels[i][2] = mag(x,y);
// // 		edgels[i][3] = phase(x,y);
// // 		bndry(x,y) = 255;
// // 		GetNextEdgel(x, y, dir, img, img.dimx(), img.dimy());
// // 		i++;
// // 	}
// // 	while(x != xInit || y != yInit);
// // 	
// // 	// Update numPts to represent the number of (valid) connected edge pts
// // 	if (i < numPts)
// // 		edgels.resize(i);
// 
// 
// /*
// 	Just for now, resize and store resized images:
// */
// 
// // 	int resizeX = int(img.dimx() * 1.1);
// // 	int resizeY = int(img.dimy() * 1.1);
// // 	CImg<float>size(resizeX, resizeY);
// // 	CImg<float> resized = img.get_resize(size, 1);
// 
// // 	CImg<float> resized(img);
// 	CImg<float> copyImg(img);
// 	
// 	
// 	cimg_mapXY(copyImg,x,y)
// 	{
// 		if(copyImg(x,y) == 0)
// 			copyImg(x,y) = 150;
// 	}
// 	
// 	CImg<float>resized = copyImg.rotate(10,2);
// 
// //	CImg<float>resized = copyImg.get_blur_median();
// 
// // 	CImg<float> resized(img);
// 	
// // 	unsigned int nBlur = 3;
// 
// 	cimg_mapXY(resized,x,y)
// 	{
// 		if(resized(x,y) == 150)
// 			resized(x,y) = 0;
// 		else
// 			resized(x,y) = 255;
// 	}
// 		
// 	//resized.blur(0.2,0.2,0.2,1);
// // 	resized.blur(3,1);
// 	
// // 	cimg_mapXY(resized,x,y)
// // 	{
// // 		if(resized(x,y) == 0)
// // 			resized(x,y) = 255;
// // 	}
// // 	
// // 	cimg_mapXY(resized,x,y)
// // 	{
// // 		if(resized(x,y) < 130)
// // 			resized(x,y) = 0;
// // 		else
// // 			resized(x,y) = 255;
// // 	}
// 	
// 	// Construct a new file name that encodes noise level
// 	char szSubfix2[1024];
// 	sprintf(szSubfix2, "n10.ppm");
// 	DirWalker dw2;
// 	dw2.ParseFileName(szPPMFileName);
// 	std::string newPath1 = info.outDir;
// 	std::string newName1 = newPath1 + dw2.name + szSubfix2; //+ dw.ext;
// 	cerr << newName1 << endl;	
// 	//pImgInfo->strFileName = newName.c_str();
// 	// Save the noisy image
// 	resized.save_pnm(newName1.c_str());
// 	
// //END for now...

	// If there is no noise to add we just return	
	if (info.m_nBumpSize == 0 && info.m_nNotchSize == 0)
	{
		pImgInfo->pField = CopyImageToField(img);
		pImgInfo->strFileName = szPPMFileName;
		pImgInfo->boundaryPts.ReSize(edgels.size());

		int numEdgels = edgels.size();

		for (int i = 0; i < numEdgels; i++)
			pImgInfo->boundaryPts[i].Set(edgels[i][0], edgels[i][1]);
	
		return true;
	}

	g_nInitShockCount = CountShockNodes(img, info);
// 	DBG_MSG1("INIT Node Count", g_nInitShockCount);

	int maxRad = MAX(info.m_nBumpSize,info.m_nNotchSize);
	int spacing = 5;
	int delta = 2 * maxRad + spacing;
	int maxNumBN = edgels.size() / delta * MAX_NUM_B_N;
	int minNumBN = 1;

	// Randomly choose some edge pts along the boundary
	IntVector pts(maxNumBN);
	IntVector::iterator start, end, it;

	// Make things really random
	srand48((long)time(NULL));
	int i;

	for (i = 0; i < maxNumBN; i++)
		pts[i] = i * delta + maxRad;

	start = pts.begin();
	end = pts.end();

	// Shuffle the elements in a random order
   	std::random_shuffle(start, end);

	// Choose how many bumps & notches we want to add
	int numBumps = int(drand48() * (maxNumBN - minNumBN)) + minNumBN;
	int numNotches = int(drand48() * (maxNumBN - minNumBN - numBumps)) + minNumBN;
	int addedSpots = 0;
	double totalSize = 0, progresSize = 0;

	if (info.m_nBumpSize == 0)
		numBumps = 0;

	if (info.m_nNotchSize == 0)
		numNotches = 0;

	ASSERT(numBumps + numNotches <= maxNumBN);
	
// 	DBG_MSG1("Num bumps to add", numBumps)
// 	DBG_MSG1("Num notches to add", numNotches)	
	

	it = pts.begin();
	
	for (int nSwitch = 0; it != pts.end() && addedSpots < numBumps + numNotches; it++, nSwitch++) 
	{
		i = *it;
		ASSERT(i < edgels.size());
		
		totalSize += AddNoiseSpot(img, info, (int)edgels[i][0], (int)edgels[i][1], edgels[i][3], 
			(/*addedSpots < numBumps*/ nSwitch%2 == 0) ? info.m_nBumpSize:info.m_nNotchSize, 
			(/*addedSpots < numBumps*/ nSwitch%2 == 0) ? 0:1, true);
			
		if(progresSize < totalSize)
		{
			addedSpots++;
			progresSize = totalSize;
		}
	}
	
	
	//DBG_VAL(addedSpots)
	
	// Comp. the level of noise introduced along the boundary
	double dNoiseLevel = (100.0 * totalSize) / edgels.size();
	
	// Construct a new file name that encodes noise level
	char szSubfix[1024];
	sprintf(szSubfix, "n%.0f.pgm", dNoiseLevel);
	
	DirWalker dw;
	dw.ParseFileName(szPPMFileName);

	std::string newPath = info.outDir;
	std::string newName = newPath + dw.name + szSubfix; //+ dw.ext;
	cerr << newName << " noise level = " << dNoiseLevel << endl;
	
	// Find the new countour with bumps and notches
	edgels = FindAndSortEdgels(img);
	int numEdgels = edgels.size();
	pImgInfo->boundaryPts.ReSize(numEdgels);
		
	for (i = 0; i < numEdgels; i++)
		pImgInfo->boundaryPts[i].Set(edgels[i][0], edgels[i][1]);
	
	pImgInfo->strFileName = newName.c_str();
	pImgInfo->pField = CopyImageToField(img);

	// Save the noisy image
	img.save_png(newName.c_str());
	
	return true;
}



