#include "tcenterlinevectP.h"

//==========================================================================

//************************

//*    Polygonization    *

//************************

//--------------------------------------------------------------------------

//===============================

//    Raw Borders Extraction

//===============================

//Raw contour class definition

class RawBorderPoint

{

	TPoint m_position;

	int m_ambiguousTurn; //used to remember cases of multiple turning directions

						 //in a RawBorder extraction.

public:

	RawBorderPoint() : m_ambiguousTurn(0) {}

	RawBorderPoint(int i, int j) : m_position(i, j), m_ambiguousTurn(0) {}

	inline TPoint pos() const { return m_position; }

	inline int x() const { return m_position.x; }

	inline int y() const { return m_position.y; }

	enum { left = 1,

		   right = 2 }; //Direction taken at ambiguous turning point

	inline int getAmbiguous() const { return m_ambiguousTurn; }

	inline void setAmbiguous(int direction) { m_ambiguousTurn = direction; }

};

//--------------------------------------------------------------------------

class RawBorder : public std::vector<rawborderpoint></rawborderpoint>

{

	int m_xExternal; //x coordinate of a specific vertex in the outer

					 //RawBorder which contains this inner one.

	TPointD *m_coordinateSums;

	TPointD *m_coordinateSquareSums;

	double *m_coordinateMixedSums;

public:

	RawBorder() {}

	~RawBorder() {}

	void setXExternalPixel(int a) { m_xExternal = a; }

	int xExternalPixel() { return m_xExternal; }

	TPointD *&sums() { return m_coordinateSums; }

	TPointD *&sums2() { return m_coordinateSquareSums; }

	double *&sumsMix() { return m_coordinateMixedSums; }

};

//--------------------------------------------------------------------------

//Of course we don't want RawBorders to be entirely copied whenever STL

//requires to resize a BorderFamily...

typedef std::vector<rawborder *=""> BorderFamily;</rawborder>

typedef std::vector<borderfamily> BorderList;</borderfamily>

//--------------------------------------------------------------------------

//==========================================================================

//============================

//    Polygonizer Locals

//============================

namespace

{

//Const names

enum { white = 0,

	   black = 1 };

enum { inner = 0,

	   outer = 1,

	   none = 2,

	   invalid = 3 };

}

//=======================================================================================

//-------------------------------

//    Raster Data Functions

//-------------------------------

//NOTA: Il tono di un TPixelCM32 rappresenta la transizione tra colore ink e colore paint.

//      di solito, se il tono e' basso abbiamo un colore ink - che puo' anche essere bianco,

//      in teoria...

//      Sarebbe opportuno che il vettorizzatore riconoscesse - per colormap - il colore

//      delle strokes.

//      Approcci: a)  la Signaturemap diventa *piu'* di una bitmap e si seguono le outline

//                    dei singoli colori.

//                 => Sconnessioni tra i colori adiacenti. Bisogna introdurre la distanza tra

//                    colori per seguire l'outline (ossia, i pixel tenuti a dx della outline

//                    devono essere simili).

//

//                b)  Una volta vettorizzato tutto, si sceglie il colore della stroke.

//                    E' possibile controllare il colore sui vertici delle sequenze semplificate

//                    e fare una media.

//NOTE: Transparency makes colors fade to white. Full transparent black pixels are considered white.

//--------------------------------------------------------------------------

template <typename t=""></typename>

class PixelEvaluator

{

	TRasterPT<t> m_ras;</t>

	int m_threshold;

public:

	PixelEvaluator(const TRasterPT<t> &ras, int threshold)</t>

		: m_ras(ras), m_threshold(threshold) {}

	inline unsigned char getBlackOrWhite(int x, int y);

};

//--------------------------------------------------------------------------

template <>

inline unsigned char PixelEvaluator<tpixel32>::getBlackOrWhite(int x, int y)</tpixel32>

{

	//return ras->pixels(y)[x].r + 2 * ras->pixels(y)[x].g + ras->pixels(y)[x].b <

	//   threshold * (ras->pixels(y)[x].m / 255.0);

	//NOTE: Green is considered twice brighter than red or blue channel.

	//Using Value of HSV color model

	return tmax(m_ras->pixels(y)[x].r, tmax(m_ras->pixels(y)[x].g, m_ras->pixels(y)[x].b)) <

		   m_threshold * (m_ras->pixels(y)[x].m / 255.0);

	//Using Lightness of HSV color model

	//return (max(ras->pixels(y)[x].r, max(ras->pixels(y)[x].g, ras->pixels(y)[x].b)) +

	//       min(ras->pixels(y)[x].r, min(ras->pixels(y)[x].g, ras->pixels(y)[x].b))) / 2.0 <

	//  threshold * (ras->pixels(y)[x].m / 255.0);

	//Using (relative) Luminance

	//return 0.2126 * ras->pixels(y)[x].r + 0.7152 * ras->pixels(y)[x].g + 0.0722 * ras->pixels(y)[x].b <

	//  threshold * (ras->pixels(y)[x].m / 255.0);

}

template <>

inline unsigned char PixelEvaluator<tpixelgr8>::getBlackOrWhite(int x, int y)</tpixelgr8>

{

	return m_ras->pixels(y)[x].value < m_threshold;

}

template <>

inline unsigned char PixelEvaluator<tpixelcm32>::getBlackOrWhite(int x, int y)</tpixelcm32>

{

	return m_ras->pixels(y)[x].getTone() < m_threshold;

}

//--------------------------------------------------------------------------

//Signaturemap format:

//  stores a map of bytes, whose first bit represents the color (black/white) of corresponding pixel, and

//  the rest its 'signature', used as an int to store information.

//NOTE: given a TRaster32, the corresponding Signaturemap constructed is intended 0(white)-padded

class Signaturemap

{

	unsigned char *m_array;

	int m_rowSize;

	int m_colSize;

public:

	Signaturemap(const TRasterP &ras, int threshold);

	~Signaturemap() { delete[] m_array; }

	template <typename t=""></typename>

	void readRasterData(const TRasterPT<t> &ras, int threshold);</t>

	inline int getRowSize() const { return m_rowSize; }

	inline int getColSize() const { return m_colSize; }

	unsigned char *pixelByte(int x, int y) { return &m_array[(y + 1) * m_rowSize + x + 1]; }

	bool getBitmapColor(int x, int y) const

	{

		return m_array[(y + 1) * m_rowSize + x + 1] & 1;

	}

	inline unsigned char getSignature(int x, int y) const

	{

		return m_array[(y + 1) * m_rowSize + x + 1] >> 1;

	}

	void setSignature(int x, int y, int val)

	{

		unsigned char *pixel = pixelByte(x, y);

		*pixel &= 1;

		*pixel |= (val << 1); //Si puo' fare meglio??

	}

};

//--------------------------------------------------------------------------

Signaturemap::Signaturemap(const TRasterP &ras, int threshold)

{

	//Extrapolate raster type

	TRaster32P rr = (TRaster32P)ras;

	TRasterGR8P rgr = (TRasterGR8P)ras;

	TRasterCM32P rt = (TRasterCM32P)ras;

	assert(rr || rgr || rt);

	//Read raster data

	if (rr) {

		rr->lock();

		readRasterData(rr, threshold);

		rr->unlock();

	} else if (rgr) {

		rgr->lock();

		readRasterData(rgr, threshold);

		rgr->unlock();

	} else {

		rt->lock();

		readRasterData(rt, threshold);

		rt->unlock();

	}

}

//--------------------------------------------------------------------------

template <typename t=""></typename>

void Signaturemap::readRasterData(const TRasterPT<t> &ras, int threshold)</t>

{

	unsigned char *currByte;

	int x, y;

	PixelEvaluator<t> evaluator(ras, threshold);</t>

	m_rowSize = ras->getLx() + 2;

	m_colSize = ras->getLy() + 2;

	m_array = new unsigned char[m_rowSize * m_colSize];

	memset(m_array, none << 1, m_rowSize);

	currByte = m_array + m_rowSize;

	for (y = 0; y < ras->getLy(); ++y) {

		*currByte = none << 1;

		currByte++;

		for (x = 0; x < ras->getLx(); ++x, ++currByte)

			*currByte = evaluator.getBlackOrWhite(x, y) | (none << 1);

		*currByte = none << 1;

		currByte++;

	}

	memset(currByte, none << 1, m_rowSize);

}

//--------------------------------------------------------------------------

//Minority check for amiguous turning directions

inline bool getMinorityCheck(const Signaturemap &ras, int x, int y)

{

	//Assumes (x,y) is ambiguous case: 2 immediate surrounding pixels are white and 2 black

	return (ras.getBitmapColor(x + 1, y) +

			ras.getBitmapColor(x + 1, y - 1) +

			ras.getBitmapColor(x - 2, y) +

			ras.getBitmapColor(x - 2, y - 1) +

			ras.getBitmapColor(x - 1, y + 1) +

			ras.getBitmapColor(x - 1, y - 2) +

			ras.getBitmapColor(x, y + 1) +

			ras.getBitmapColor(x, y - 2)) > 4;

}

//--------------------------------------------------------------------------

//Sets signature of a given border

inline void setSignature(Signaturemap &ras, const RawBorder &border, int val)

{

	unsigned int j;

	int yOld;

	//Set border's alpha channel

	yOld = border.back().y();

	for (j = 0; j < border.size(); ++j) {

		if (border[j].y() < yOld) {

			ras.setSignature(border[j].x(), border[j].y(), val);

		} else if (border[j].y() > yOld) {

			ras.setSignature(border[j].x(), yOld, val);

		}

		yOld = border[j].y();

	}

}

//==========================================================================

//-------------------------------

//    Raw Borders Extraction

//-------------------------------

//RawBorderPoints correspond to lower-left pixel corners.

//EXAMPLE: (0,0) is the lower-left *corner* of the image, whereas (0,0) also

// represents coordinates of the lower-left *pixel*.

//NOTE: 'Ambiguous turning' vertices are those of kind:

//

//    B|W           W|B

//    -x-    -or-   -x-

//    W|B           B|W

//

//Keeping B on the right of our path-seeking direction, we may either turn

//left or right at these points.

RawBorder *extractPath(Signaturemap &ras, int x0, int y0,

					   int pathType, int xOuterPixel, int despeckling)

{

	RawBorder *path = new RawBorder;

	int x, y;

	short dirX, dirY;

	long int area = 0;

	bool nextLeftPixel, nextRightPixel;

	if (pathType == outer) {

		dirX = 0;

		dirY = 1;

	} else {

		dirX = 1;

		dirY = 0;

		area += y0;

		path->setXExternalPixel(xOuterPixel);

	}

	path->push_back(RawBorderPoint(x0, y0));

	//Check here if (x0, y0) is an ambiguous-direction point

	nextLeftPixel = ras.getBitmapColor(x0 + (dirY - dirX - 1) / 2, y0 + (-dirY - dirX - 1) / 2);

	nextRightPixel = ras.getBitmapColor(x0 + (-dirX - dirY - 1) / 2, y0 + (dirX - dirY - 1) / 2);

	if ((nextRightPixel == black) && (nextLeftPixel == white))

		path->back().setAmbiguous(dirX ? RawBorderPoint::left : RawBorderPoint::right);

	//Begin path extraction

	for (x = x0 + dirX, y = y0 + dirY; !(x == x0 && y == y0); x += dirX, y += dirY) {

		path->push_back(RawBorderPoint(x, y));

		//Calculate next direction

		nextLeftPixel = ras.getBitmapColor(x + (dirX - dirY - 1) / 2, y + (dirY + dirX - 1) / 2);

		nextRightPixel = ras.getBitmapColor(x + (dirX + dirY - 1) / 2, y + (dirY - dirX - 1) / 2);

		if ((nextRightPixel == black) && (nextLeftPixel == black)) {

			//Left Turn

			std::swap(dirY, dirX);

			dirX = -dirX;

		} else if ((nextRightPixel == white) && (nextLeftPixel == white)) {

			//Right Turn

			std::swap(dirY, dirX);

			dirY = -dirY;

		} else if ((nextRightPixel == white) && (nextLeftPixel == black)) {

			//path->back().setAmbiguous();

			//Do a surrounding check and connect minority color

			if (getMinorityCheck(ras, x, y) == black) {

				std::swap(dirY, dirX);

				dirY = -dirY;

				path->back().setAmbiguous(RawBorderPoint::right);

			} //right turn

			else {

				std::swap(dirY, dirX);

				dirX = -dirX;

				path->back().setAmbiguous(RawBorderPoint::left);

			} //left turn

		}

		//Also calculate border area

		area += y * dirX;

		//And sign treated pixel

		if (dirY != 0)

			ras.setSignature(x, y + (dirY - 1) / 2, pathType);

	}

	//If the inner region's overall area is under a given threshold,

	//then erase it (intended as image noise).

	if (abs(area) < despeckling) {

		setSignature(ras, *path, invalid);

		delete path;

		path = 0;

	}

	return path;

}

//--------------------------------------------------------------------------

BorderList *extractBorders(const TRasterP &ras, int threshold, int despeckling)

{

	Signaturemap byteImage(ras, threshold);

	BorderList *borderHierarchy = new BorderList;

	vector<rawborder *=""> outerBorders;</rawborder>

	list<rawborder *=""> innerBorders;</rawborder>

	RawBorder *foundPath;

	int x, y;

	bool Color, oldColor;

	int xOuterPixel = 0;

	bool enteredRegionType;

	unsigned char signature;

	//Traverse image to extract raw borders

	for (y = 0; y < ras->getLy(); ++y) {

		oldColor = white;

		enteredRegionType = outer;

		for (x = 0; x < ras->getLx(); ++x) {

			if (oldColor ^ (Color = byteImage.getBitmapColor(x, y))) {

				//Region type changes

				enteredRegionType = !enteredRegionType;

				if ((signature = byteImage.getSignature(x, y)) == none) {

					//We've found a border

					if (foundPath = extractPath(byteImage, x, y, !enteredRegionType, xOuterPixel, despeckling))

						if (enteredRegionType == outer)

							innerBorders.push_back(foundPath);

						else

							outerBorders.push_back(foundPath);

				}

				//If leaving a white region, remember it - in order to establish

				//border hierarchy in the future

				if (enteredRegionType == inner && signature != invalid)

					xOuterPixel = x;

				//Invalid pixels got signed by a cut-out path, due to insufficient area

				if (signature == invalid)

					byteImage.setSignature(x, y, none); //Restore them now

				oldColor = Color;

			}

		}

	}

	//Now, we have all borders found, but no hierarchy between them.

	unsigned int i;

	list<rawborder *="">::iterator l;</rawborder>

	//Build hierarchy

	innerBorders.push_front(0); //Just to keep a fixed list head

	for (i = 0; i < outerBorders.size(); ++i) {

		//Initialize a border family

		borderHierarchy->push_back(BorderFamily());

		borderHierarchy->back().push_back(outerBorders[i]);

		//Reset outerBorders[i]'s signature

		setSignature(byteImage, *outerBorders[i], none);

		//Now check inner borders for insideness - check if the outerPixel

		//remembered in path extraction has been cleared

		for (l = innerBorders.begin(), ++l; l != innerBorders.end(); ++l) {

			if (byteImage.getSignature((*l)->xExternalPixel(), (**l)[0].y()) == none) {

				borderHierarchy->back().push_back(*l);

				setSignature(byteImage, **l, none);

				l = innerBorders.erase(l);

				--l;

			}

		}

	}

	return borderHierarchy;

}

//--------------------------------------------------------------------------

//==================================

//    Calculate optimal polygons

//==================================

//The optimal polygon for a given original border is found like:

// 1) Find couples (i,k(i)), so that k(i) be the largest k:

//        d(j,ik) <= 1;  for *all* i

//    It can be shown that such a condition is equivalent to:

//      exists line l : d(l,j)<=1/2, for all i<=j<=k(i).

// 2) Clean the above couples - find couples (i,l(i)):

//        l(i)=min{k(j)},  j=i..n.

// 3) Calculate clipped couples (i',l'); where i'=i+1, l'=l(i)-1.

// 4) Calculate sums for path penalties.

// 5) Apply optimality algorithm.

//NOTE: Weak simpleness reads like: a set of polygons is weak-simple if no edge

//      *crosses* another edge. Superposition and collision of edges with vertices

//      are still admitted.

//        =>  It can be shown that due to 1) and special conditions on ambiguous

//            turnings applied in both 1) and 3), weak simpleness is insured in

//            our polygonization.

//--------------------------------------------------------------------------

//Helper functions/classes: circular-indexed vectors

//returns 1 whenever the triple (a,b,c) is 'circular' mod n.

//NOTE: We'll find useful taking (i,i,j) as 1 and (i,j,j) as 0.

inline bool isCircular(int a, int b, int c)

{

	return a <= c ? a <= b && b < c : c > b || b >= a;

}

//--------------------------------------------------------------------------

//Extracts a 'next corner' array - helps improving overall speed

inline int *findNextCorners(RawBorder &path)

{

	int i, currentCorner;

	int *corners = new int[path.size()];

	//NOTE: 0 is a corner, due to the path extraction procedure.

	currentCorner = 0;

	for (i = path.size() - 1; i >= 0; --i) {

		if (path[currentCorner].x() != path[i].x() &&

			path[currentCorner].y() != path[i].y())

			currentCorner = i + 1;

		corners[i] = currentCorner;

	}

	return corners;

}

//--------------------------------------------------------------------------

//Calculate furthest k satisfying 1) for all fixed i.

inline int *furthestKs(RawBorder &path, int *&nextCorners)

{

	int n = path.size();

	int *kVector = new int[n];

	enum { left,

		   up,

		   right,

		   down };

	int directionsOccurred[4];

	nextCorners = findNextCorners(path);

	int i, j, k;

	TPoint shift;

	TPoint leftConstraint, rightConstraint, violatedConstraint;

	TPoint newLeftConstraint, newRightConstraint;

	TPoint jPoint, jNextPoint, iPoint, direction;

	int directionSignature;

	for (i = 0; i < n; ++i) {

		//Initialize search

		leftConstraint = rightConstraint = TPoint();

		directionsOccurred[0] = directionsOccurred[1] = directionsOccurred[2] =

			directionsOccurred[3] = 0;

		j = i;

		jNextPoint = iPoint = path[i].pos();

		//Search for k(i)

		while (1) {

			//NOTE: Here using TPoint::operator= is less effective than setting

			//its x and y components directly...

			jPoint = jNextPoint;

			jNextPoint = path[nextCorners[j]].pos();

			//Update directions count

			directionSignature = jNextPoint.x > jPoint.x ? right : jNextPoint.x < jPoint.x ? left

																						   : jNextPoint.y > jPoint.y ? up : down;

			directionsOccurred[directionSignature] = 1;

			//If all 4 axis directions occurred, quit

			if (directionsOccurred[left] && directionsOccurred[right] &&

				directionsOccurred[up] && directionsOccurred[down]) {

				k = j;

				goto foundK;

			}

			//Update displacement from i

			shift = jNextPoint - iPoint;

			//Test j against constraints

			//if(cross(shift, leftConstraint)<0 || cross(shift, rightConstraint)>0)

			if (cross(shift, leftConstraint) < 0) {

				violatedConstraint = leftConstraint;

				break;

			}

			if (cross(shift, rightConstraint) > 0) {

				violatedConstraint = rightConstraint;

				break;

			}

			//Update constraints

			if (abs(shift.x) > 1 || abs(shift.y) > 1) {

				newLeftConstraint.x = shift.x +

									  (shift.y < 0 || (shift.y == 0 && shift.x < 0) ? 1 : -1);

				newLeftConstraint.y = shift.y +

									  (shift.x > 0 || (shift.x == 0 && shift.y < 0) ? 1 : -1);

				if (cross(newLeftConstraint, leftConstraint) >= 0)

					leftConstraint = newLeftConstraint;

				newRightConstraint.x = shift.x +

									   (shift.y > 0 || (shift.y == 0 && shift.x < 0) ? 1 : -1);

				newRightConstraint.y = shift.y +

									   (shift.x < 0 || (shift.x == 0 && shift.y < 0) ? 1 : -1);

				if (cross(newRightConstraint, rightConstraint) <= 0)

					rightConstraint = newRightConstraint;

			}

			//Imposing strict constraint for ambiguous turnings, to ensure polygons' weak simpleness.

			//Has to be defined *outside* abs checks.

			if (path[nextCorners[j]].getAmbiguous()) {

				if (path[nextCorners[j]].getAmbiguous() == RawBorderPoint::left)

					rightConstraint = shift;

				else

					leftConstraint = shift;

			}

			j = nextCorners[j];

		}

		//At this point, constraints are violated by the next corner.

		//Then, search for the last k between j and corners[j] not violating them.

		direction = convert(normalize(convert(jNextPoint - jPoint)));

		k = (j + cross(jPoint - iPoint, violatedConstraint) / cross(violatedConstraint, direction)) % n;

	foundK:

		kVector[i] = k;

	}

	return kVector;

}

//--------------------------------------------------------------------------

// Now find the effective intervals inside which we can define possible

// arcs approximating the given raw border:

//    for every a in [i,res[i]], the arc connecting border[i] and

//    border[a] will be a possible one.

inline int *calculateForwardArcs(RawBorder &border, bool ambiguitiesCheck)

{

	int i, j, n = (int)border.size();

	int *nextCorners;

	int *k = furthestKs(border, nextCorners);

	int *K = new int[n];

	int *res = new int[n];

	//find K[i]= min {k[j]}, j=i..n-1.

	for (i = 0; i < n; ++i) {

		for (j = i, K[i] = k[i]; isCircular(i, j, K[i]); j = (j + 1) % n)

			if (isCircular(j, k[j], K[i]))

				K[i] = k[j];

	}

	//Finally, we perform the following clean-up operations:

	//  first, extremities of [i,K[i]] are clipped away, to obtain a

	//    smoother optimal polygon (and deal with cases like the unitary

	//    square);

	//  second, arcs of the kind [i,j] with j

	for (i = n - 1, j = 0; j < n; i = j, ++j) {

		res[j] = K[i] < j ? (K[i] == 0 ? n - 1 : n) : K[i] - 1;

	}

	//Amibiguities check for vertex and edge superpositions. Prevent problems in the forecoming

	//straight-skeleton thinning process.

	if (ambiguitiesCheck) {

		for (i = 1; nextCorners[i] > 0; i = nextCorners[i]) {

			if (border[i].getAmbiguous() == RawBorderPoint::right) {

				//Check vertices from i (excluded) to res[res[i]]; if in it there exists vertex k so that pos(k)==pos(i)...

				//This prevents the existence of 0 degree angles in the optimal polygon.

				int rrPlus1 = (res[res[i] % n] + 1) % n;

				for (j = nextCorners[i];

					 isCircular(i, j, rrPlus1) && j != i; //remember that isCircular(a,a,b) == 1 ...

					 j = nextCorners[j]) {

					if (border[j].getAmbiguous() && (border[j].pos() == border[i].pos())) {

						res[res[i] % n] = j - 1;

						assert((res[i] % n) != j - 1);

						//Further, ensure res is increasing

						for (int k = res[i] % n; res[k] >= j - 1 && k >= 0; --k) {

							res[k] = j - 1;

							assert(k != j - 1);

						}

						break;

					}

				}

			}

		}

	}

	delete[] k;

	delete[] K;

	delete[] nextCorners;

	return res;

}

//--------------------------------------------------------------------------

//Let sum[i] and sum2[i] be respectively the sums of vertex coordinates

//from 0 to i, and the sums of their squares; sumsMix contain sums of

//xy terms.

inline void calculateSums(RawBorder &path)

{

	unsigned int i, n = path.size();

	TPointD currentRelativePos;

	path.sums() = new TPointD[n + 1];

	path.sums2() = new TPointD[n + 1];

	path.sumsMix() = new double[n + 1];

	path.sums()[0].x = path.sums()[0].y = path.sums2()[0].x = path.sums2()[0].y = 0;

	for (i = 1; i < path.size(); ++i) {

		currentRelativePos = convert(path[i].pos() - path[0].pos());

		path.sums()[i] = path.sums()[i - 1] + currentRelativePos;

		path.sums2()[i].x = path.sums2()[i - 1].x + currentRelativePos.x * currentRelativePos.x;

		path.sums2()[i].y = path.sums2()[i - 1].y + currentRelativePos.y * currentRelativePos.y;

		path.sumsMix()[i] = path.sumsMix()[i - 1] + currentRelativePos.x * currentRelativePos.y;

	}

	// path[n] is virtually intended as path[0], but we prefer to introduce

	// it in the optimality algorithm's count

	path.sums()[n].x = path.sums()[n].y = path.sums2()[n].x = path.sums2()[n].y = 0;

}

//--------------------------------------------------------------------------

//Let a,b the index-extremities of an arc of this path.

//Then return its penalty.

inline double penalty(RawBorder &path, int a, int b)

{

	int n = b - a + 1;

	TPointD v = convert(rotate90(path[b == path.size() ? 0 : b].pos() - path[a].pos()));

	TPointD sum = path.sums()[b] - path.sums()[a];

	TPointD sum2 = path.sums2()[b] - path.sums2()[a];

	double sumMix = path.sumsMix()[b] - path.sumsMix()[a];

	double F1 = sum2.x - 2 * sum.x * path[a].x() + n * path[a].x() * path[a].x();

	double F2 = sum2.y - 2 * sum.y * path[a].y() + n * path[a].y() * path[a].y();

	double F3 = sumMix - sum.x * path[a].y() - sum.y * path[a].x() + n * path[a].x() * path[a].y();

	return sqrt((v.y * v.y * F1 + v.x * v.x * F2 - 2 * v.x * v.y * F3) / n);

}

//--------------------------------------------------------------------------

//NOTA: Il seguente algoritmo di riduzione assicura la semplicita' (debole) dei poligoni prodotti.

//

inline void reduceBorder(RawBorder &border, Contour &res, bool ambiguitiesCheck)

{

	int i, j, k, a, *b, n = border.size(), m;

	double newPenalty;

	int minPenaltyNext;

	int *minPenaltyNextArray = new int[n];

	//Calculate preliminary infos

	int *longestArcFrom = calculateForwardArcs(border, ambiguitiesCheck);

	calculateSums(border);

	double *penaltyToEnd = new double[n + 1];

	//EXPLANATION:

	//The fastest way to extract the optimal reduced border is based on the

	//weakly monotonic property of longestArc[].

	//The minimal number of its vertices 'm' is easily found by

	//traversing the path with the longest step allowed. Let b[] be that

	//succession; then, given res[i], it has to be reached by a vertex in

	//the interval: {a[i-1], .. , b[i-1]}, where longestArc[a[i-1]]=a[i],

	//longestArc[a[i-1]-1]

	//Calculate m

	for (i = 0, m = 0; i < n; i = longestArcFrom[i])

		++m;

	//Calculate b[]

	b = new int[m + 1];

	b[m] = n;

	for (i = 0, j = 0; j < m; i = longestArcFrom[i], ++j)

		b[j] = i;

	//NOTE: a[] need not be completely found - we just remember the

	//a=a[j+1] currently needed.

	//Now, build the optimal polygon

	for (j = m - 1, a = n; j >= 0; --j) {

		for (k = b[j]; k >= 0 && longestArcFrom[k] >= a; --k) {

			penaltyToEnd[k] = infinity;

			for (i = a; i <= longestArcFrom[k]; ++i) {

				newPenalty = penaltyToEnd[i] + penalty(border, k, i);

				if (newPenalty < penaltyToEnd[k])

					penaltyToEnd[k] = newPenalty;

				minPenaltyNext = i;

			}

			minPenaltyNextArray[k] = minPenaltyNext;

		}

		a = k + 1;

	}

	//Finally, read off the optimal polygon

	res.resize(m);

	for (i = j = 0; i < n; i = minPenaltyNextArray[i], ++j) {

		res[j] = ContourNode(border[i].x(), border[i].y());

		//Ambiguities are still remembered in the output polygon.

		if (border[i].getAmbiguous() == RawBorderPoint::left)

			res[j].setAttribute(ContourNode::AMBIGUOUS_LEFT);

		if (border[i].getAmbiguous() == RawBorderPoint::right)

			res[j].setAttribute(ContourNode::AMBIGUOUS_RIGHT);

	}

	delete[] longestArcFrom;

	delete[] minPenaltyNextArray;

	delete[] border.sums();

	delete[] border.sums2();

	delete[] border.sumsMix();

	delete[] penaltyToEnd;

	delete[] b;

}

//--------------------------------------------------------------------------

//Reduction caller and list copier.

inline void reduceBorders(BorderList &borders, Contours &result, bool ambiguitiesCheck)

{

	unsigned int i, j;

	//Initialize output container

	result.resize(borders.size());

	//Copy results

	for (i = 0; i < borders.size(); ++i) {

		result[i].resize(borders[i].size());

		for (j = 0; j < borders[i].size(); ++j) {

			reduceBorder(*borders[i][j], result[i][j], ambiguitiesCheck);

			delete borders[i][j];

		}

	}

}

//--------------------------------------------------------------------------

//===========================

//    Polygonization Main

//===========================

//Extracts a polygonal, minimal yet faithful representation of image contours

//Contours* polygonize(const TRasterP &ras){

void polygonize(const TRasterP &ras, Contours &polygons, VectorizerCoreGlobals &g)

{

	BorderList *borders;

	borders = extractBorders(ras, g.currConfig->m_threshold, g.currConfig->m_despeckling);

	reduceBorders(*borders, polygons, g.currConfig->m_maxThickness > 0.0);