ED.cpp

#include "ED.h"
#include "EDColor.h"
#include <fstream>

using namespace cv;
using namespace std;

ED::ED(Mat _srcImage, GradientOperator _op, int _gradThresh, int _anchorThresh,int _scanInterval, int _minPathLen ,double _sigma, bool _sumFlag)
{	
	// Check parameters for sanity
	if (_gradThresh < 1) _gradThresh = 1;
	if (_anchorThresh < 0) _anchorThresh = 0;
	if (_sigma < 1.0) _sigma = 1.0;

	srcImage = _srcImage;

	height = srcImage.rows;
	width = srcImage.cols;
	
	op = _op;
	gradThresh = _gradThresh;
	anchorThresh = _anchorThresh;
	scanInterval = _scanInterval;
	minPathLen = _minPathLen;
	sigma = _sigma;
	sumFlag = _sumFlag;

	segmentNos = 0;
	segmentPoints.push_back(vector<Point>()); // create empty vector of points for segments

	edgeImage = Mat(height, width, CV_8UC1, Scalar(0)); // initialize edge Image
	smoothImage = Mat(height, width, CV_8UC1);
	gradImage = Mat(height, width, CV_16SC1); // gradImage contains short values 

	srcImg = srcImage.data;

	//// Detect Edges By Edge Drawing Algorithm  ////
	
	/*------------ SMOOTH THE IMAGE BY A GAUSSIAN KERNEL -------------------*/
	if (sigma == 1.0)
		GaussianBlur(srcImage, smoothImage, Size(5, 5), sigma);
	else
		GaussianBlur(srcImage, smoothImage, Size(), sigma); // calculate kernel from sigma

	// Assign Pointers from Mat's data
	smoothImg = smoothImage.data;
	gradImg = (short*)gradImage.data;
	edgeImg = edgeImage.data;

	dirImg = new unsigned char[width*height];

	/*------------ COMPUTE GRADIENT & EDGE DIRECTION MAPS -------------------*/
	ComputeGradient();

	/*------------ COMPUTE ANCHORS -------------------*/
	ComputeAnchorPoints();

	/*------------ JOIN ANCHORS -------------------*/
	JoinAnchorPointsUsingSortedAnchors();
	
	delete[] dirImg;
}

// This constructor for use of EDLines and EDCircle with ED given as constructor argument
// only the necessary attributes are coppied
ED::ED(const ED & cpyObj)
{
	height = cpyObj.height;
	width = cpyObj.width;

	srcImage = cpyObj.srcImage.clone();
	
	op = cpyObj.op;
	gradThresh = cpyObj.gradThresh;
	anchorThresh = cpyObj.anchorThresh;
	scanInterval = cpyObj.scanInterval;
	minPathLen = cpyObj.minPathLen;
	sigma = cpyObj.sigma;
	sumFlag = cpyObj.sumFlag;

	edgeImage = cpyObj.edgeImage.clone();
	smoothImage = cpyObj.smoothImage.clone();
	gradImage = cpyObj.gradImage.clone();

	srcImg = srcImage.data;

	smoothImg = smoothImage.data;
	gradImg = (short*)gradImage.data;
	edgeImg = edgeImage.data;

	segmentPoints = cpyObj.segmentPoints;
	segmentNos = cpyObj.segmentNos;
}

// This constructor for use of EDColor with use of direction and gradient image
// It finds edge image for given gradient and direction image
ED::ED(short *_gradImg, uchar *_dirImg, int _width, int _height, int _gradThresh, int _anchorThresh, int _scanInterval, int _minPathLen, bool selectStableAnchors)
{
	height = _height;
	width = _width;

	gradThresh = _gradThresh;
	anchorThresh = _anchorThresh;
	scanInterval = _scanInterval;
	minPathLen = _minPathLen;

	gradImg = _gradImg;
	dirImg = _dirImg;

	edgeImage = Mat(height, width, CV_8UC1, Scalar(0)); // initialize edge Image

	edgeImg = edgeImage.data;

	if (selectStableAnchors) {

		// Compute anchors with the user supplied parameters
		anchorThresh = 0; // anchorThresh used as zero while computing anchor points if selectStableAnchors set. 
						  // Finding higher number of anchors is OK, because we have following validation steps in selectStableAnchors.
		ComputeAnchorPoints();
		anchorThresh = _anchorThresh; // set it to its initial argument value for further anchor validation.
		anchorPoints.clear(); // considering validation step below, it should constructed again.

		for (int i = 1; i<height - 1; i++) {
			for (int j = 1; j<width - 1; j++) {
				if (edgeImg[i*width + j] != ANCHOR_PIXEL) continue;

				// Take only "stable" anchors
				// 0 degree edge
				if (edgeImg[i*width + j - 1] && edgeImg[i*width + j + 1]) {
					int diff1 = gradImg[i*width + j] - gradImg[(i - 1)*width + j];
					int diff2 = gradImg[i*width + j] - gradImg[(i + 1)*width + j];
					if (diff1 >= anchorThresh && diff2 >= anchorThresh) edgeImg[i*width + j] = 255;

					continue;
				} //end-if

				  // 90 degree edge
				if (edgeImg[(i - 1)*width + j] && edgeImg[(i + 1)*width + j]) {
					int diff1 = gradImg[i*width + j] - gradImg[i*width + j - 1];
					int diff2 = gradImg[i*width + j] - gradImg[i*width + j + 1];
					if (diff1 >= anchorThresh && diff2 >= anchorThresh) edgeImg[i*width + j] = 255;

					continue;
				} //end-if

				  // 135 degree diagonal
				if (edgeImg[(i - 1)*width + j - 1] && edgeImg[(i + 1)*width + j + 1]) {
					int diff1 = gradImg[i*width + j] - gradImg[(i - 1)*width + j + 1];
					int diff2 = gradImg[i*width + j] - gradImg[(i + 1)*width + j - 1];
					if (diff1 >= anchorThresh && diff2 >= anchorThresh) edgeImg[i*width + j] = 255;
					continue;
				} //end-if

				  // 45 degree diagonal
				if (edgeImg[(i - 1)*width + j + 1] && edgeImg[(i + 1)*width + j - 1]) {
					int diff1 = gradImg[i*width + j] - gradImg[(i - 1)*width + j - 1];
					int diff2 = gradImg[i*width + j] - gradImg[(i + 1)*width + j + 1];
					if (diff1 >= anchorThresh && diff2 >= anchorThresh) edgeImg[i*width + j] = 255;
				} //end-if

			} //end-for
		} //end-for

		for (int i = 0; i<width*height; i++) 
			if (edgeImg[i] == ANCHOR_PIXEL) 
				edgeImg[i] = 0; 
			else if (edgeImg[i] == 255) {
				edgeImg[i] = ANCHOR_PIXEL;
				int y = i / width;
				int x = i % width;
				anchorPoints.push_back(Point(x, y)); // push validated anchor point to vector
			}
			
		anchorNos = (int)anchorPoints.size(); // get # of anchor pixels
	}

	else {
		// Compute anchors with the user supplied parameters
		ComputeAnchorPoints(); // anchorThresh used as given as argument. No validation applied. (No stable anchors.)
	} //end-else

	segmentNos = 0;
	segmentPoints.push_back(vector<Point>()); // create empty vector of points for segments

	JoinAnchorPointsUsingSortedAnchors();
}

ED::ED(EDColor &obj) 
{
	width = obj.getWidth();
	height = obj.getHeight();
	segmentPoints = obj.getSegments();
	segmentNos = obj.getSegmentNo();
}

ED::ED()
{
	//
}


Mat ED::getEdgeImage()
{
	return edgeImage;
}

Mat ED::getAnchorImage()
{
	Mat anchorImage = Mat(edgeImage.size(), edgeImage.type(), Scalar(0));

	std::vector<Point>::iterator it;

	for (it = anchorPoints.begin(); it != anchorPoints.end(); it++)
		anchorImage.at<uchar>(*it) = 255;

	return anchorImage;
}

Mat ED::getSmoothImage()
{
	return smoothImage;
}

Mat ED::getGradImage()
{	
	Mat result8UC1;
	convertScaleAbs(gradImage, result8UC1);
	
	return result8UC1;
}

int ED::getSegmentNo()
{
	return segmentNos;
}

int ED::getAnchorNo()
{
	return anchorNos;
}

std::vector<Point> ED::getAnchorPoints()
{
	return anchorPoints;
}

std::vector<std::vector<Point>> ED::getSegments()
{
	return segmentPoints;
}

std::vector<std::vector<Point>> ED::getSortedSegments()
{
		// sort segments from largest to smallest
		std::sort(segmentPoints.begin(), segmentPoints.end(), [](const std::vector<Point> & a, const std::vector<Point> & b) { return a.size() > b.size(); });

		return segmentPoints;
}

Mat ED::drawParticularSegments(std::vector<int> list)
{
	Mat segmentsImage = Mat(edgeImage.size(), edgeImage.type(), Scalar(0));

	std::vector<Point>::iterator it;
	std::vector<int>::iterator itInt;

	for (itInt = list.begin(); itInt != list.end(); itInt++)
		for (it = segmentPoints[*itInt].begin(); it != segmentPoints[*itInt].end(); it++)
			segmentsImage.at<uchar>(*it) = 255;
	
	return segmentsImage;
}


void ED::ComputeGradient()
{	
	// Initialize gradient image for row = 0, row = height-1, column=0, column=width-1 
	for (int j = 0; j<width; j++) { gradImg[j] = gradImg[(height - 1)*width + j] = gradThresh - 1; }
	for (int i = 1; i<height - 1; i++) { gradImg[i*width] = gradImg[(i + 1)*width - 1] = gradThresh - 1; }

	for (int i = 1; i<height - 1; i++) {
		for (int j = 1; j<width - 1; j++) {
			// Prewitt Operator in horizontal and vertical direction
			// A B C
			// D x E
			// F G H
			// gx = (C-A) + (E-D) + (H-F)
			// gy = (F-A) + (G-B) + (H-C)
			//
			// To make this faster: 
			// com1 = (H-A)
			// com2 = (C-F)
			// 
			// For Prewitt
			// Then: gx = com1 + com2 + (E-D) = (H-A) + (C-F) + (E-D) = (C-A) + (E-D) + (H-F)
			//       gy = com1 - com2 + (G-B) = (H-A) - (C-F) + (G-B) = (F-A) + (G-B) + (H-C)
			// 
			// For Sobel
			// Then: gx = com1 + com2 + 2*(E-D) = (H-A) + (C-F) + 2*(E-D) = (C-A) + 2*(E-D) + (H-F)
			//       gy = com1 - com2 + 2*(G-B) = (H-A) - (C-F) + 2*(G-B) = (F-A) + 2*(G-B) + (H-C)
			//
			// For Scharr
			// Then: gx = 3*(com1 + com2) + 10*(E-D) = 3*(H-A) + 3*(C-F) + 10*(E-D) = 3*(C-A) + 10*(E-D) + 3*(H-F)
			//       gy = 3*(com1 - com2) + 10*(G-B) = 3*(H-A) - 3*(C-F) + 10*(G-B) = 3*(F-A) + 10*(G-B) + 3*(H-C)
			//
			// For LSD
			// A B
			// C D
			// gx = (B-A) + (D-C)
			// gy = (C-A) + (D-B)
			//
			// To make this faster: 
			// com1 = (D-A)
			// com2 = (B-C)
			// Then: gx = com1 + com2 = (D-A) + (B-C) = (B-A) + (D-C)
			//       gy = com1 - com2 = (D-A) - (B-C) = (C-A) + (D-B)

			int com1 = smoothImg[(i + 1)*width + j + 1] - smoothImg[(i - 1)*width + j - 1];
			int com2 = smoothImg[(i - 1)*width + j + 1] - smoothImg[(i + 1)*width + j - 1];

			int gx;
			int gy;
			
			switch (op)
			{
			case PREWITT_OPERATOR:
				gx = abs(com1 + com2 + (smoothImg[i*width + j + 1] - smoothImg[i*width + j - 1]));
				gy = abs(com1 - com2 + (smoothImg[(i + 1)*width + j] - smoothImg[(i - 1)*width + j]));
				break;
			case SOBEL_OPERATOR:
				gx = abs(com1 + com2 + 2 * (smoothImg[i*width + j + 1] - smoothImg[i*width + j - 1]));
				gy = abs(com1 - com2 + 2 * (smoothImg[(i + 1)*width + j] - smoothImg[(i - 1)*width + j]));
				break;
			case SCHARR_OPERATOR:
				gx = abs(3 * (com1 + com2) + 10 * (smoothImg[i*width + j + 1] - smoothImg[i*width + j - 1]));
				gy = abs(3 * (com1 - com2) + 10 * (smoothImg[(i + 1)*width + j] - smoothImg[(i - 1)*width + j]));
				break;
			case LSD_OPERATOR:
				// com1 and com2 differs from previous operators, because LSD has 2x2 kernel
				int com1 = smoothImg[(i + 1)*width + j + 1] - smoothImg[i*width + j];
				int com2 = smoothImg[i*width + j + 1] - smoothImg[(i + 1)*width + j];

				gx = abs(com1 + com2);
				gy = abs(com1 - com2);
			}
			
			int sum;

			if(sumFlag)
				sum = gx + gy;
			else
				sum = (int)sqrt((double)gx*gx + gy*gy);

			int index = i*width + j;
			gradImg[index] = sum;

			if (sum >= gradThresh) {
				if (gx >= gy) dirImg[index] = EDGE_VERTICAL;
				else          dirImg[index] = EDGE_HORIZONTAL;
			} //end-if
		} // end-for
	} // end-for
}

void ED::ComputeAnchorPoints()
{
	//memset(edgeImg, 0, width*height);
	for (int i = 2; i<height - 2; i++) {
		int start = 2;
		int inc = 1;
		if (i%scanInterval != 0) { start = scanInterval; inc = scanInterval; }

		for (int j = start; j<width - 2; j += inc) {
			if (gradImg[i*width + j] < gradThresh) continue;

			if (dirImg[i*width + j] == EDGE_VERTICAL) {
				// vertical edge
				int diff1 = gradImg[i*width + j] - gradImg[i*width + j - 1];
				int diff2 = gradImg[i*width + j] - gradImg[i*width + j + 1];
				if (diff1 >= anchorThresh && diff2 >= anchorThresh) {
					edgeImg[i*width + j] = ANCHOR_PIXEL;
					anchorPoints.push_back(Point(j, i)); 
				}

			}
			else {
				// horizontal edge
				int diff1 = gradImg[i*width + j] - gradImg[(i - 1)*width + j];
				int diff2 = gradImg[i*width + j] - gradImg[(i + 1)*width + j];
				if (diff1 >= anchorThresh && diff2 >= anchorThresh) {
					edgeImg[i*width + j] = ANCHOR_PIXEL;
					anchorPoints.push_back(Point(j, i)); 
				}
			} // end-else
		} //end-for-inner
	} //end-for-outer

	anchorNos = (int)anchorPoints.size(); // get the total number of anchor points
}

void ED::JoinAnchorPointsUsingSortedAnchors()
{
	int *chainNos = new int[(width + height) * 8];

	Point *pixels = new Point[width*height];
	StackNode *stack = new StackNode[width*height];
	Chain *chains = new Chain[width*height];

	// sort the anchor points by their gradient value in decreasing order
	int *A = sortAnchorsByGradValue1();

	// Now join the anchors starting with the anchor having the greatest gradient value
	int totalPixels = 0;

	for (int k = anchorNos - 1; k >= 0; k--) {
		int pixelOffset = A[k];

		int i = pixelOffset / width;
		int j = pixelOffset % width;

		//int i = anchorPoints[k].y;
		//int j = anchorPoints[k].x;

		if (edgeImg[i*width + j] != ANCHOR_PIXEL) continue;

		chains[0].len = 0;
		chains[0].parent = -1;
		chains[0].dir = 0;
		chains[0].children[0] = chains[0].children[1] = -1;
		chains[0].pixels = NULL;


		int noChains = 1;
		int len = 0;
		int duplicatePixelCount = 0;
		int top = -1;  // top of the stack 

		if (dirImg[i*width + j] == EDGE_VERTICAL) {
			stack[++top].r = i;
			stack[top].c = j;
			stack[top].dir = ED_DOWN;
			stack[top].parent = 0;

			stack[++top].r = i;
			stack[top].c = j;
			stack[top].dir = ED_UP;
			stack[top].parent = 0;

		}
		else {
			stack[++top].r = i;
			stack[top].c = j;
			stack[top].dir = ED_RIGHT;
			stack[top].parent = 0;

			stack[++top].r = i;
			stack[top].c = j;
			stack[top].dir = ED_LEFT;
			stack[top].parent = 0;
		} //end-else

		  // While the stack is not empty
	StartOfWhile:
		while (top >= 0) {
			int r = stack[top].r;
			int c = stack[top].c;
			int dir = stack[top].dir;
			int parent = stack[top].parent;
			top--;

			if (edgeImg[r*width + c] != EDGE_PIXEL) duplicatePixelCount++;

			chains[noChains].dir = dir;   // traversal direction
			chains[noChains].parent = parent;
			chains[noChains].children[0] = chains[noChains].children[1] = -1;


			int chainLen = 0;

			chains[noChains].pixels = &pixels[len];

			pixels[len].y = r;
			pixels[len].x = c;
			len++;
			chainLen++;

			if (dir == ED_LEFT) {
				while (dirImg[r*width + c] == EDGE_HORIZONTAL) {
					edgeImg[r*width + c] = EDGE_PIXEL;

					// The edge is horizontal. Look LEFT
					//
					//   A
					//   B x 
					//   C 
					//
					// cleanup up & down pixels
					if (edgeImg[(r - 1)*width + c] == ANCHOR_PIXEL) edgeImg[(r - 1)*width + c] = 0;
					if (edgeImg[(r + 1)*width + c] == ANCHOR_PIXEL) edgeImg[(r + 1)*width + c] = 0;

					// Look if there is an edge pixel in the neighbors
					if (edgeImg[r*width + c - 1] >= ANCHOR_PIXEL) { c--; }
					else if (edgeImg[(r - 1)*width + c - 1] >= ANCHOR_PIXEL) { r--; c--; }
					else if (edgeImg[(r + 1)*width + c - 1] >= ANCHOR_PIXEL) { r++; c--; }
					else {
						// else -- follow max. pixel to the LEFT
						int A = gradImg[(r - 1)*width + c - 1];
						int B = gradImg[r*width + c - 1];
						int C = gradImg[(r + 1)*width + c - 1];

						if (A > B) {
							if (A > C) r--;
							else       r++;
						}
						else  if (C > B) r++;
						c--;
					} //end-else

					if (edgeImg[r*width + c] == EDGE_PIXEL || gradImg[r*width + c] < gradThresh) {
						if (chainLen > 0) {
							chains[noChains].len = chainLen;
							chains[parent].children[0] = noChains;
							noChains++;
						} // end-if
						goto StartOfWhile;
					} //end-else


					pixels[len].y = r;
					pixels[len].x = c;
					len++;
					chainLen++;
				} //end-while

				stack[++top].r = r;
				stack[top].c = c;
				stack[top].dir = ED_DOWN;
				stack[top].parent = noChains;

				stack[++top].r = r;
				stack[top].c = c;
				stack[top].dir = ED_UP;
				stack[top].parent = noChains;

				len--;
				chainLen--;

				chains[noChains].len = chainLen;
				chains[parent].children[0] = noChains;
				noChains++;

			}
			else if (dir == ED_RIGHT) {
				while (dirImg[r*width + c] == EDGE_HORIZONTAL) {
					edgeImg[r*width + c] = EDGE_PIXEL;

					// The edge is horizontal. Look RIGHT
					//
					//     A
					//   x B
					//     C
					//
					// cleanup up&down pixels
					if (edgeImg[(r + 1)*width + c] == ANCHOR_PIXEL) edgeImg[(r + 1)*width + c] = 0;
					if (edgeImg[(r - 1)*width + c] == ANCHOR_PIXEL) edgeImg[(r - 1)*width + c] = 0;

					// Look if there is an edge pixel in the neighbors
					if (edgeImg[r*width + c + 1] >= ANCHOR_PIXEL) { c++; }
					else if (edgeImg[(r + 1)*width + c + 1] >= ANCHOR_PIXEL) { r++; c++; }
					else if (edgeImg[(r - 1)*width + c + 1] >= ANCHOR_PIXEL) { r--; c++; }
					else {
						// else -- follow max. pixel to the RIGHT
						int A = gradImg[(r - 1)*width + c + 1];
						int B = gradImg[r*width + c + 1];
						int C = gradImg[(r + 1)*width + c + 1];

						if (A > B) {
							if (A > C) r--;       // A
							else       r++;       // C
						}
						else if (C > B) r++;  // C
						c++;
					} //end-else

					if (edgeImg[r*width + c] == EDGE_PIXEL || gradImg[r*width + c] < gradThresh) {
						if (chainLen > 0) {
							chains[noChains].len = chainLen;
							chains[parent].children[1] = noChains;
							noChains++;
						} // end-if
						goto StartOfWhile;
					} //end-else


					pixels[len].y = r;
					pixels[len].x = c;
					len++;
					chainLen++;
				} //end-while

				stack[++top].r = r;
				stack[top].c = c;
				stack[top].dir = ED_DOWN;  // Go down
				stack[top].parent = noChains;

				stack[++top].r = r;
				stack[top].c = c;
				stack[top].dir = ED_UP;   // Go up
				stack[top].parent = noChains;

				len--;
				chainLen--;

				chains[noChains].len = chainLen;
				chains[parent].children[1] = noChains;
				noChains++;

			}
			else if (dir == ED_UP) {
				while (dirImg[r*width + c] == EDGE_VERTICAL) {
					edgeImg[r*width + c] = EDGE_PIXEL;

					// The edge is vertical. Look UP
					//
					//   A B C
					//     x
					//
					// Cleanup left & right pixels
					if (edgeImg[r*width + c - 1] == ANCHOR_PIXEL) edgeImg[r*width + c - 1] = 0;
					if (edgeImg[r*width + c + 1] == ANCHOR_PIXEL) edgeImg[r*width + c + 1] = 0;

					// Look if there is an edge pixel in the neighbors
					if (edgeImg[(r - 1)*width + c] >= ANCHOR_PIXEL) { r--; }
					else if (edgeImg[(r - 1)*width + c - 1] >= ANCHOR_PIXEL) { r--; c--; }
					else if (edgeImg[(r - 1)*width + c + 1] >= ANCHOR_PIXEL) { r--; c++; }
					else {
						// else -- follow the max. pixel UP
						int A = gradImg[(r - 1)*width + c - 1];
						int B = gradImg[(r - 1)*width + c];
						int C = gradImg[(r - 1)*width + c + 1];

						if (A > B) {
							if (A > C) c--;
							else       c++;
						}
						else if (C > B) c++;
						r--;
					} //end-else

					if (edgeImg[r*width + c] == EDGE_PIXEL || gradImg[r*width + c] < gradThresh) {
						if (chainLen > 0) {
							chains[noChains].len = chainLen;
							chains[parent].children[0] = noChains;
							noChains++;
						} // end-if
						goto StartOfWhile;
					} //end-else


					pixels[len].y = r;
					pixels[len].x = c;

					len++;
					chainLen++;
				} //end-while

				stack[++top].r = r;
				stack[top].c = c;
				stack[top].dir = ED_RIGHT;
				stack[top].parent = noChains;

				stack[++top].r = r;
				stack[top].c = c;
				stack[top].dir = ED_LEFT;
				stack[top].parent = noChains;

				len--;
				chainLen--;

				chains[noChains].len = chainLen;
				chains[parent].children[0] = noChains;
				noChains++;

			}
			else { // dir == DOWN
				while (dirImg[r*width + c] == EDGE_VERTICAL) {
					edgeImg[r*width + c] = EDGE_PIXEL;

					// The edge is vertical
					//
					//     x
					//   A B C
					//
					// cleanup side pixels
					if (edgeImg[r*width + c + 1] == ANCHOR_PIXEL) edgeImg[r*width + c + 1] = 0;
					if (edgeImg[r*width + c - 1] == ANCHOR_PIXEL) edgeImg[r*width + c - 1] = 0;

					// Look if there is an edge pixel in the neighbors
					if (edgeImg[(r + 1)*width + c] >= ANCHOR_PIXEL) { r++; }
					else if (edgeImg[(r + 1)*width + c + 1] >= ANCHOR_PIXEL) { r++; c++; }
					else if (edgeImg[(r + 1)*width + c - 1] >= ANCHOR_PIXEL) { r++; c--; }
					else {
						// else -- follow the max. pixel DOWN
						int A = gradImg[(r + 1)*width + c - 1];
						int B = gradImg[(r + 1)*width + c];
						int C = gradImg[(r + 1)*width + c + 1];

						if (A > B) {
							if (A > C) c--;       // A
							else       c++;       // C
						}
						else if (C > B) c++;  // C
						r++;
					} //end-else

					if (edgeImg[r*width + c] == EDGE_PIXEL || gradImg[r*width + c] < gradThresh) {
						if (chainLen > 0) {
							chains[noChains].len = chainLen;
							chains[parent].children[1] = noChains;
							noChains++;
						} // end-if
						goto StartOfWhile;
					} //end-else

					pixels[len].y = r;
					pixels[len].x = c;

					len++;
					chainLen++;
				} //end-while

				stack[++top].r = r;
				stack[top].c = c;
				stack[top].dir = ED_RIGHT;
				stack[top].parent = noChains;

				stack[++top].r = r;
				stack[top].c = c;
				stack[top].dir = ED_LEFT;
				stack[top].parent = noChains;

				len--;
				chainLen--;

				chains[noChains].len = chainLen;
				chains[parent].children[1] = noChains;
				noChains++;
			} //end-else

		} //end-while


		if (len - duplicatePixelCount < minPathLen) {
			for (int k = 0; k<len; k++) {

				edgeImg[pixels[k].y*width + pixels[k].x] = 0;
				edgeImg[pixels[k].y*width + pixels[k].x] = 0;

			} //end-for

		}
		else {

			int noSegmentPixels = 0;

			int totalLen = LongestChain(chains, chains[0].children[1]);

			if (totalLen > 0) {
				// Retrieve the chainNos
				int count = RetrieveChainNos(chains, chains[0].children[1], chainNos);

				// Copy these pixels in the reverse order
				for (int k = count - 1; k >= 0; k--) {
					int chainNo = chainNos[k];

#if 1
					/* See if we can erase some pixels from the last chain. This is for cleanup */

					int fr = chains[chainNo].pixels[chains[chainNo].len - 1].y;
					int fc = chains[chainNo].pixels[chains[chainNo].len - 1].x;

					int index = noSegmentPixels - 2;
					while (index >= 0) {
						int dr = abs(fr - segmentPoints[segmentNos][index].y);
						int dc = abs(fc - segmentPoints[segmentNos][index].x);

						if (dr <= 1 && dc <= 1) {
							// neighbors. Erase last pixel
							segmentPoints[segmentNos].pop_back();
							noSegmentPixels--;
							index--;
						}
						else break;
					} //end-while

					if (chains[chainNo].len > 1 && noSegmentPixels > 0) {
						fr = chains[chainNo].pixels[chains[chainNo].len - 2].y;
						fc = chains[chainNo].pixels[chains[chainNo].len - 2].x;

						int dr = abs(fr - segmentPoints[segmentNos][noSegmentPixels - 1].y);
						int dc = abs(fc - segmentPoints[segmentNos][noSegmentPixels - 1].x);

						if (dr <= 1 && dc <= 1) chains[chainNo].len--;
					} //end-if
#endif

					for (int l = chains[chainNo].len - 1; l >= 0; l--) {
						segmentPoints[segmentNos].push_back(chains[chainNo].pixels[l]);
						noSegmentPixels++;
					} //end-for

					chains[chainNo].len = 0;  // Mark as copied
				} //end-for
			} //end-if

			totalLen = LongestChain(chains, chains[0].children[0]);
			if (totalLen > 1) {
				// Retrieve the chainNos
				int count = RetrieveChainNos(chains, chains[0].children[0], chainNos);

				// Copy these chains in the forward direction. Skip the first pixel of the first chain
				// due to repetition with the last pixel of the previous chain
				int lastChainNo = chainNos[0];
				chains[lastChainNo].pixels++;
				chains[lastChainNo].len--;

				for (int k = 0; k<count; k++) {
					int chainNo = chainNos[k];

#if 1
					/* See if we can erase some pixels from the last chain. This is for cleanup */
					int fr = chains[chainNo].pixels[0].y;
					int fc = chains[chainNo].pixels[0].x;

					int index = noSegmentPixels - 2;
					while (index >= 0) {
						int dr = abs(fr - segmentPoints[segmentNos][index].y);
						int dc = abs(fc - segmentPoints[segmentNos][index].x);

						if (dr <= 1 && dc <= 1) {
							// neighbors. Erase last pixel
							segmentPoints[segmentNos].pop_back();
							noSegmentPixels--;
							index--;
						}
						else break;
					} //end-while

					int startIndex = 0;
					int chainLen = chains[chainNo].len;
					if (chainLen > 1 && noSegmentPixels > 0) {
						int fr = chains[chainNo].pixels[1].y;
						int fc = chains[chainNo].pixels[1].x;

						int dr = abs(fr - segmentPoints[segmentNos][noSegmentPixels - 1].y);
						int dc = abs(fc - segmentPoints[segmentNos][noSegmentPixels - 1].x);

						if (dr <= 1 && dc <= 1) { startIndex = 1; }
					} //end-if
#endif

					  /* Start a new chain & copy pixels from the new chain */
					for (int l = startIndex; l<chains[chainNo].len; l++) {
						segmentPoints[segmentNos].push_back(chains[chainNo].pixels[l]);
						noSegmentPixels++;
					} //end-for

					chains[chainNo].len = 0;  // Mark as copied
				} //end-for
			} //end-if


			  // See if the first pixel can be cleaned up
			int fr = segmentPoints[segmentNos][1].y;
			int fc = segmentPoints[segmentNos][1].x;


			int dr = abs(fr - segmentPoints[segmentNos][noSegmentPixels - 1].y);
			int dc = abs(fc - segmentPoints[segmentNos][noSegmentPixels - 1].x);


			if (dr <= 1 && dc <= 1) {
				segmentPoints[segmentNos].erase(segmentPoints[segmentNos].begin());
				noSegmentPixels--;
			} //end-if

			segmentNos++;
			segmentPoints.push_back(vector<Point>()); // create empty vector of points for segments

													  // Copy the rest of the long chains here
			for (int k = 2; k<noChains; k++) {
				if (chains[k].len < 2) continue;

				totalLen = LongestChain(chains, k);

				if (totalLen >= 10) {

					// Retrieve the chainNos
					int count = RetrieveChainNos(chains, k, chainNos);

					// Copy the pixels
					noSegmentPixels = 0;
					for (int k = 0; k<count; k++) {
						int chainNo = chainNos[k];

#if 1					
						/* See if we can erase some pixels from the last chain. This is for cleanup */
						int fr = chains[chainNo].pixels[0].y;
						int fc = chains[chainNo].pixels[0].x;

						int index = noSegmentPixels - 2;
						while (index >= 0) {
							int dr = abs(fr - segmentPoints[segmentNos][index].y);
							int dc = abs(fc - segmentPoints[segmentNos][index].x);

							if (dr <= 1 && dc <= 1) {
								// neighbors. Erase last pixel
								segmentPoints[segmentNos].pop_back();
								noSegmentPixels--;
								index--;
							}
							else break;
						} //end-while

						int startIndex = 0;
						int chainLen = chains[chainNo].len;
						if (chainLen > 1 && noSegmentPixels > 0) {
							int fr = chains[chainNo].pixels[1].y;
							int fc = chains[chainNo].pixels[1].x;

							int dr = abs(fr - segmentPoints[segmentNos][noSegmentPixels - 1].y);
							int dc = abs(fc - segmentPoints[segmentNos][noSegmentPixels - 1].x);

							if (dr <= 1 && dc <= 1) { startIndex = 1; }
						} //end-if
#endif
						  /* Start a new chain & copy pixels from the new chain */
						for (int l = startIndex; l<chains[chainNo].len; l++) {
							segmentPoints[segmentNos].push_back(chains[chainNo].pixels[l]);
							noSegmentPixels++;
						} //end-for

						chains[chainNo].len = 0;  // Mark as copied
					} //end-for
					segmentPoints.push_back(vector<Point>()); // create empty vector of points for segments
					segmentNos++;
				} //end-if          
			} //end-for

		} //end-else

	} //end-for-outer

	// pop back last segment from vector
	// because of one preallocation in the beginning, it will always empty
	segmentPoints.pop_back();

	// Clean up
	delete[] A;
	delete[] chains;
	delete[] stack;
	delete[] chainNos;
	delete[] pixels;
}

void ED::sortAnchorsByGradValue()
{
	auto sortFunc = [&](const Point &a, const Point &b)
	{
		return gradImg[a.y * width + a.x] > gradImg[b.y * width + b.x];
	};
	
	std::sort(anchorPoints.begin(), anchorPoints.end(), sortFunc);

	/*
	ofstream myFile;
	myFile.open("anchorsNew.txt");
	for (int i = 0; i < anchorPoints.size(); i++) {
		int x = anchorPoints[i].x;
		int y = anchorPoints[i].y;

		myFile << i << ". value: " << gradImg[y*width + x] << "  Cord: (" << x << "," << y << ")" << endl;
	}
	myFile.close(); 
	
	
	vector<Point> temp(anchorNos);

	int x, y, i = 0;
	char c;
	std::ifstream infile("cords.txt");
	while (infile >> x >> c >> y && c == ',') {
		temp[i] = Point(x, y);
		i++;
	} 

	anchorPoints = temp; 
	*/
}

int * ED::sortAnchorsByGradValue1()
{
	int SIZE = 128 * 256;
	int *C = new int[SIZE];
	memset(C, 0, sizeof(int)*SIZE);

	// Count the number of grad values
	for (int i = 1; i<height - 1; i++) {
		for (int j = 1; j<width - 1; j++) {
			if (edgeImg[i*width + j] != ANCHOR_PIXEL) continue;

			int grad = gradImg[i*width + j];
			C[grad]++;
		} //end-for
	} //end-for 

	// Compute indices
	for (int i = 1; i<SIZE; i++) C[i] += C[i - 1];

	int noAnchors = C[SIZE - 1];
	int *A = new int[noAnchors];
	memset(A, 0, sizeof(int)*noAnchors);


	for (int i = 1; i<height - 1; i++) {
		for (int j = 1; j<width - 1; j++) {
			if (edgeImg[i*width + j] != ANCHOR_PIXEL) continue;

			int grad = gradImg[i*width + j];
			int index = --C[grad];
			A[index] = i*width + j;    // anchor's offset 
		} //end-for
	} //end-for  

	delete[] C;

	/*
	ofstream myFile;
	myFile.open("aNew.txt");
	for (int i = 0; i < noAnchors; i++)
		myFile << A[i] << endl;

	myFile.close(); */

	return A;

}


int ED::LongestChain(Chain *chains, int root) {
	if (root == -1 || chains[root].len == 0) return 0;

	int len0 = 0;
	if (chains[root].children[0] != -1) len0 = LongestChain(chains, chains[root].children[0]);

	int len1 = 0;
	if (chains[root].children[1] != -1) len1 = LongestChain(chains, chains[root].children[1]);

	int max = 0;

	if (len0 >= len1) {
		max = len0;
		chains[root].children[1] = -1;

	}
	else {
		max = len1;
		chains[root].children[0] = -1;
	} //end-else

	return chains[root].len + max;
} //end-LongestChain

int ED::RetrieveChainNos(Chain * chains, int root, int chainNos[])
{
	int count = 0;

	while (root != -1) {
		chainNos[count] = root;
		count++;

		if (chains[root].children[0] != -1) root = chains[root].children[0];
		else                                root = chains[root].children[1];
	} //end-while

	return count;
}