#include "toonz/tcenterlinevectorizer.h"
#include "tsystem.h"
#include "tstopwatch.h"
#include "tpalette.h"
//#include "timage_io.h"
//#include "tstrokeutil.h"
//#include "tspecialstyleid.h"
#include "trastercm.h"
#include "ttoonzimage.h"
//#include "dpiscale.h"
#include "tregion.h"
#include "tstroke.h"

#include "trasterimage.h"

#define SEARCH_WINDOW 3
#define JUNCTION_THICKNESS_RATIO 0.7
#define JOIN_LIMIT 0.8
#define THICKNESS_LIMIT 0.2
#define THICKNESS_LIMIT_2 0.04

#undef DEBUG
//---------------------------------------------------------

struct ControlPoint {
	TStroke *m_stroke;
	int m_index;
	ControlPoint(TStroke *stroke, int index)
		: m_stroke(stroke), m_index(index)
	{
	}
	TPointD getPoint() const
	{
		return m_stroke->getControlPoint(m_index);
	}
	void setPoint(const TPointD &p)
	{
		TThickPoint point = m_stroke->getControlPoint(m_index);
		point.x = p.x;
		point.y = p.y;
		m_stroke->setControlPoint(m_index, point);
	}
};

void renormalizeImage(TVectorImage *vi)
{
	int i, j;
	int n = vi->getStrokeCount();
	std::vector<ControlPoint> points;
	points.reserve(n * 2);
	for (i = 0; i < n; i++) {
		TStroke *stroke = vi->getStroke(i);
		int m = stroke->getControlPointCount();
		if (m > 0) {
			if (m == 1)
				points.push_back(ControlPoint(stroke, 0));
			else {
				points.push_back(ControlPoint(stroke, 0));
				points.push_back(ControlPoint(stroke, m - 1));
			}
		}
	}
	int count = points.size();
	for (i = 0; i < count; i++) {
		ControlPoint &pi = points[i];
		TPointD posi = pi.getPoint();
		TPointD center = posi;
		std::vector<int> neighbours;
		neighbours.push_back(i);
		for (j = i + 1; j < count; j++) {
			TPointD posj = points[j].getPoint();
			double d = tdistance(posj, posi);
			if (d < 0.01) {
				/*if(d>0)
             {
              int yy = 123;
             }*/
				neighbours.push_back(j);
				center += posj;
			}
		}
		int m = neighbours.size();
		if (m == 1)
			continue;
		center = center * (1.0 / m);
		for (j = 0; j < m; j++)
			points[neighbours[j]].setPoint(center);
	}
}

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

class Node;

class DataPixel
{
public:
	TPoint m_pos;
	int m_value;
	bool m_ink;
	//int m_index;
	Node *m_node;
	DataPixel() : m_value(0), m_ink(false), /*m_index(0), */ m_node(0) {}
};

//---------------------------------------------------------
#ifdef WIN32
template class DV_EXPORT_API TSmartPointerT<TRasterT<DataPixel>>;
#endif
typedef TRasterPT<DataPixel> DataRasterP;

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

class Junction;

class Node
{
public:
	Node *m_other;
	DataPixel *m_pixel;
	Node *m_prev, *m_next;
	Junction *m_junction;
#ifdef DEBUG
	bool m_flag;
#endif
	bool m_visited;
	Node()
		: m_pixel(0),
		  m_prev(0),
		  m_next(0),
		  m_junction(0),
#ifdef DEBUG
		  m_flag(false),
#endif
		  m_visited(false)
	{
	}
};

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

class ProtoStroke;

class Junction
{
public:
	TThickPoint m_center;
	std::deque<Node *> m_nodes;
	int m_junctionOrder;
	std::vector<ProtoStroke *> m_protoStrokes;
	bool m_locked;
	Junction()
		: m_center(), m_nodes(), m_junctionOrder(0), m_protoStrokes(), m_locked(false) {}
	bool isConvex();
};

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

class ProtoStroke
{
public:
	TPointD m_startDirection, m_endDirection;
	Junction *m_startJunction, *m_endJunction;
	std::deque<TThickPoint> m_points;
	ProtoStroke()
		: m_points(), m_startDirection(), m_endDirection(), m_startJunction(0), m_endJunction(0) {}
	ProtoStroke(std::deque<TThickPoint>::iterator it_b,
				std::deque<TThickPoint>::iterator it_e)
		: m_points(it_b, it_e), m_startDirection(), m_endDirection(), m_startJunction(0), m_endJunction(0) {}
};

class ProtoRegion
{
public:
	TPointD m_center;
	bool m_isConvex;
	std::vector<TThickPoint> m_points;
	ProtoRegion(bool isConvex)
		: m_points(), m_isConvex(isConvex), m_center() {}
};

class JunctionMerge
{
public:
	Junction *m_junction;
	Node *m_startNode, *m_endNode;
	bool m_isNext;
	JunctionMerge(Junction *junction)
		: m_junction(junction), m_startNode(0), m_endNode(0), m_isNext(true) {}
	JunctionMerge(Junction *junction, Node *startNode, Node *endNode, bool
																		  isNext)
		: m_junction(junction), m_startNode(startNode), m_endNode(endNode), m_isNext(isNext) {}
	bool operator==(const JunctionMerge &op2)
	{
		return (m_junction == op2.m_junction);
	}
};

class JunctionLink
{
public:
	Junction *m_first;
	Junction *m_second;
	int m_order;
	JunctionLink(Junction *first, Junction *second)
		: m_first(first), m_second(second), m_order(0) {}
	bool operator==(const JunctionLink &op2) const
	{
		return (m_first == op2.m_first && m_second == op2.m_second) || (m_first ==
																			op2.m_second &&
																		m_second == op2.m_first);
	}
};

class CenterLineVectorizer
{
	TPalette *m_palette;

public:
	//	int m_delta[8];
	//	int m_radius;
	TRasterP m_src;

	VectorizerConfiguration m_configuration;

#ifdef DEBUG
	TRaster32P m_output;
#endif
	DataRasterP m_dataRaster;
	vector<pair<int, DataRasterP>> m_dataRasterArray;
	TVectorImageP m_vimage;
#ifdef DEBUG
	TVectorImageP m_rimage;
#endif

#ifdef DEBUG
	vector<vector<TThickPoint>> m_outlines;
	vector<TPointD> m_unvisitedNodes;
	vector<vector<Node *>> m_junctionPolygons;
#endif

	vector<Node *> m_nodes;
	vector<Junction *> m_junctions;
	vector<JunctionLink> m_links;
	list<ProtoStroke> m_protoStrokes;
	list<ProtoRegion> m_protoRegions;
	//list< std::vector<TThickPoint> > m_protoHoles;

	Node *findOtherSide(Node *node);
	bool testOtherSide(Node *na1, Node *nb1, double &startDist2);

	TPointD computeCenter(Node *na, Node *nb, double &r);

	bool isHole(Node *startNode);
	void traceLine(DataPixel *pix); // , vector<TThickPoint> &points);

	TPoint computeGradient(DataPixel *pix)
	{
		assert(m_dataRaster);
		const int wrap = m_dataRaster->getWrap();

		TPoint g(0, 0);
		int n, s, w, e, nw, sw, ne, se;

		w = pix[-1].m_value;
		nw = pix[-1 + wrap].m_value;
		sw = pix[-1 - wrap].m_value;

		e = pix[+1].m_value;
		ne = pix[+1 + wrap].m_value;
		se = pix[+1 - wrap].m_value;

		n = pix[+wrap].m_value;
		s = pix[-wrap].m_value;

		g.y = -sw + ne - se + nw + 2 * (n - s);
		g.x = -sw + ne + se - nw + 2 * (e - w);
		return g;
	}

	CenterLineVectorizer(const VectorizerConfiguration &configuration, TPalette *palette)
		: m_configuration(configuration), m_palette(palette) {}

	~CenterLineVectorizer();

	void clearNodes();
	void createNodes();

	void clearJunctions();

	void makeOutputRaster();
	void makeDataRaster(const TRasterP &src);

	DataPixel *findUnvisitedInkPixel();

	bool linkNextProtoStroke(
		ProtoStroke *const dstProtoStroke,
		Junction *const currJunction,
		const int k);
	void joinProtoStrokes(ProtoStroke *const dstProtoStroke);

	void resolveUnvisitedNode(Node *node);
	Junction *mergeJunctions(const std::list<JunctionMerge> &junctions);
	void joinJunctions();
	void createJunctionPolygon(Junction *junction);
	void handleLinks();
	void handleJunctions();

	void createStrokes();
	void createRegions();

	void vectorize();
	void init();
	//	void click(const TPoint &p);

	Node *createNode(DataPixel *pix);

	void link(DataPixel *pix, DataPixel *from, DataPixel *to);

	/*#ifdef DEBUG
	inline TPixel32 &getOutPix(DataPixel*pix) {return 
m_output->pixels()[pix->m_index]; }
	void outputNodes(Node *node);
	void outputInks();
#endif*/

private:
	// not implemented
	CenterLineVectorizer(const CenterLineVectorizer &);
	CenterLineVectorizer &operator=(const CenterLineVectorizer &);
};

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

class CompareJunctionNodes
{
	TPointD m_center;

public:
	CompareJunctionNodes(TPointD center)
		: m_center(center) {}

	bool operator()(const Node *n1, const Node *n2)
	{
		TPointD p1 = convert(n1->m_pixel->m_pos) - m_center;
		TPointD p2 = convert(n2->m_pixel->m_pos) - m_center;

		double alpha1, alpha2;

		if (p1.x > 0)
			alpha1 = -p1.y / sqrt(norm2(p1));
		else if (p1.x < 0)
			alpha1 = 2 + p1.y / sqrt(norm2(p1));
		else //(p1.x = 0)
		{
			if (p1.y > 0)
				alpha1 = -1;
			else if (p1.y < 0)
				alpha1 = 1;
			else
				assert(true);
		}

		if (p2.x > 0)
			alpha2 = -p2.y / sqrt(norm2(p2));
		else if (p2.x < 0)
			alpha2 = 2 + p2.y / sqrt(norm2(p2));
		else //(p2.x = 0)
		{
			if (p2.y > 0)
				alpha2 = -1;
			else if (p2.y < 0)
				alpha2 = 1;
			else
				assert(true);
		}

		if (alpha2 - alpha1 > 0)
			return true;
		else if (alpha2 - alpha1 < 0)
			return false;
		else
			return false; // n1->m_pixel->m_pos == n2->m_pixel->m_pos!!
	}
};

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

inline bool collinear(const TPointD &a, const TPointD &b, const TPointD &c)
{
	double area = (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y);
	return area == 0;
}

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

inline bool right(const TPointD &a, const TPointD &b, const TPointD &c)
{
	double area = (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y);
	return area < 0;
}

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

enum RectIntersectionResult {
	parallel,
	coincident,
	intersected
};

RectIntersectionResult rectsIntersect(const TPointD &a1, const TPointD &a2,
									  const TPointD &b1, const TPointD &b2,
									  TPointD &intersection)
{
	double ma, mb, ca, cb;

	if (a2.x == a1.x && b2.x == b1.x) {
		if (collinear(a1, a2, b1) &&
			collinear(a1, a2, b2))
			return coincident;
		else
			return parallel;
	} else if (a1.x == a2.x) {
		mb = (b2.y - b1.y) / (b2.x - b1.x);
		cb = b1.y - mb * b1.x;
		intersection.x = a1.x;
		intersection.y = mb * intersection.x + cb;
		return intersected;
	} else if (b1.x == b2.x) {
		ma = (a2.y - a1.y) / (a2.x - a1.x);
		ca = a1.y - ma * a1.x;
		intersection.x = b1.x;
		intersection.y = ma * intersection.x + ca;
		return intersected;
	} else {
		ma = (a2.y - a1.y) / (a2.x - a1.x);
		mb = (b2.y - b1.y) / (b2.x - b1.x);

		if (ma == mb) {
			if (collinear(a1, a2, b1) &&
				collinear(a1, a2, b2))
				return coincident;
			else
				return parallel;
		}

		ca = a1.y - ma * a1.x;
		cb = b1.y - mb * b1.x;

		intersection.x = (ca - cb) / (mb - ma);
		intersection.y = ma * intersection.x + ca;
		return intersected;
	}
}

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

bool Junction::isConvex()
{
	int nonConvexAngles = 0;
	assert(m_nodes.size() > 2);
	std::deque<Node *>::iterator it = m_nodes.begin();
	std::deque<Node *>::iterator it_e = m_nodes.end() - 2;
	for (; it != it_e; ++it) {
		if (right(
				convert((*(it))->m_pixel->m_pos),
				convert((*(it + 1))->m_pixel->m_pos),
				convert((*(it + 2))->m_pixel->m_pos)))
			nonConvexAngles++;
		if (nonConvexAngles > 1)
			return false;
	}
	if (right(
			convert((*(it))->m_pixel->m_pos),
			convert((*(it + 1))->m_pixel->m_pos),
			convert((*(m_nodes.begin()))->m_pixel->m_pos)))
		nonConvexAngles++;
	if (nonConvexAngles > 1)
		return false;
	return true;
}

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

CenterLineVectorizer::~CenterLineVectorizer()
{
	clearNodes();
	clearJunctions();

#ifdef DEBUG
	m_outlines.clear();
	m_unvisitedNodes.clear();
	m_junctionPolygons.clear();
#endif

	m_links.clear();
	m_protoStrokes.clear();
	m_protoRegions.clear();
	//m_protoHoles.clear();
}

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

void CenterLineVectorizer::clearNodes()
{
	for (int i = 0; i < (int)m_nodes.size(); i++)
		delete m_nodes[i];
	m_nodes.clear();
}

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

void CenterLineVectorizer::clearJunctions()
{
	for (int i = 0; i < (int)m_junctions.size(); i++)
		delete m_junctions[i];
	m_junctions.clear();
}

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

void CenterLineVectorizer::link(
	DataPixel *pix,
	DataPixel *srcPix,
	DataPixel *dstPix)
{
	Node *srcNode = 0, *dstNode = 0, *node = 0;
	Node *tmp;
	for (tmp = pix->m_node; tmp; tmp = tmp->m_other) {
		if (tmp->m_pixel == 0)
			continue;
		if (tmp->m_prev && tmp->m_prev->m_pixel == srcPix) {
			assert(srcNode == 0);
			if (node) {
				assert(node->m_next->m_pixel == dstPix);
				assert(node->m_prev == 0);
				node->m_prev = tmp->m_prev;
				tmp->m_prev->m_next = node;
				tmp->m_next = tmp->m_prev = 0;
				tmp->m_pixel = 0;
				return;
			}
			assert(tmp->m_next == 0);
			srcNode = tmp->m_prev;
			node = tmp;
		}
		if (tmp->m_next && tmp->m_next->m_pixel == dstPix) {
			assert(dstNode == 0);
			if (node) {
				assert(node->m_prev->m_pixel == srcPix);
				assert(node->m_next == 0);
				node->m_next = tmp->m_next;
				tmp->m_next->m_prev = node;
				tmp->m_next = tmp->m_prev = 0;
				tmp->m_pixel = 0;
				return;
			}
			assert(tmp->m_prev == 0);
			dstNode = tmp->m_next;
			node = tmp;
		}
	}
	if (!node)
		node = createNode(pix);
	if (!srcNode)
		srcNode = createNode(srcPix);
	if (!dstNode)
		dstNode = createNode(dstPix);

	if (!node->m_next) {
		node->m_next = dstNode;
		assert(dstNode->m_prev == 0);
		dstNode->m_prev = node;
	}
	if (!node->m_prev) {
		node->m_prev = srcNode;
		assert(srcNode->m_next == 0);
		srcNode->m_next = node;
	}

	assert(node->m_next == dstNode);
	assert(node->m_prev == srcNode);
	assert(dstNode->m_prev == node);
	assert(srcNode->m_next == node);
}

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

inline int colorDistance2(const TPixel32 &c0, const TPixel32 &c1)
{
	return ((c0.r - c1.r) * (c0.r - c1.r) + (c0.g - c1.g) * (c0.g - c1.g) + (c0.b - c1.b) * (c0.b - c1.b));
}

//---------------------------------------------------------
#define MAX_TOLERANCE 20

#include "tcolorstyles.h"

void CenterLineVectorizer::makeDataRaster(const TRasterP &src)
{
	m_vimage = new TVectorImage();
	if (!src)
		return;
	m_src = src;

	clearNodes();
	clearJunctions();

	int ii, x, y = 0;
	TRaster32P srcRGBM = (TRaster32P)m_src;
	TRasterCM32P srcCM = (TRasterCM32P)m_src;
	TRasterGR8P srcGR = (TRasterGR8P)m_src;

	DataRasterP dataRaster(m_src->getSize().lx + 2, m_src->getSize().ly + 2);
	int ly = dataRaster->getLy();
	int lx = dataRaster->getLx();
	int wrap = dataRaster->getWrap();
	DataPixel *dataPix0 = dataRaster->pixels(0);
	DataPixel *dataPix1 = dataRaster->pixels(0) + m_src->getLx() + 1;
	for (y = 0; y < ly; y++, dataPix0 += wrap, dataPix1 += wrap) {
		dataPix0->m_pos.x = 0;
		dataPix1->m_pos.x = lx - 1;
		dataPix0->m_pos.y = dataPix1->m_pos.y = y;
		dataPix0->m_value = dataPix1->m_value = 0;
		dataPix0->m_ink = dataPix1->m_ink = false;
		dataPix0->m_node = dataPix1->m_node = 0;
	}
	dataPix0 = dataRaster->pixels(0);
	dataPix1 = dataRaster->pixels(ly - 1);
	for (x = 0; x < lx; x++, dataPix0++, dataPix1++) {
		dataPix0->m_pos.x = dataPix1->m_pos.x = x;
		dataPix0->m_pos.y = 0;
		dataPix1->m_pos.y = ly - 1;
		dataPix0->m_value = dataPix1->m_value = 0;
		dataPix0->m_ink = dataPix1->m_ink = false;
		dataPix0->m_node = dataPix1->m_node = 0;
	}

	if (srcRGBM) {
		assert(m_palette);
		int inkId = m_palette->addStyle(m_configuration.m_inkColor);

		m_dataRasterArray.push_back(pair<int, DataRasterP>(inkId, dataRaster));
		int maxDistance2 = m_configuration.m_threshold * m_configuration.m_threshold;

		for (y = 0; y < m_src->getLy(); y++) {
			TPixel32 *inPix = srcRGBM->pixels(y);
			TPixel32 *inEndPix = inPix + srcRGBM->getLx();
			DataPixel *dataPix = dataRaster->pixels(y + 1) + 1;
			x = 0;
			while (inPix < inEndPix) {
				*dataPix = DataPixel();
				int distance2 = colorDistance2(m_configuration.m_inkColor, *inPix);

				//int value = (inPix->r + 2*inPix->g + inPix->b)>>2;
				if (y == 0 || y == m_src->getLy() - 1 || x == 0 || x == m_src->getLx() - 1) {
					dataPix->m_value = 255;
					dataPix->m_ink = false;
				} else {
					dataPix->m_value = (inPix->r + 2 * inPix->g + inPix->b) >> 2;
					dataPix->m_ink = (distance2 < maxDistance2);
				}
				dataPix->m_pos.x = x++;
				dataPix->m_pos.y = y;
				dataPix->m_node = 0;
				//dataPix->m_index = index++;
				inPix++;
				dataPix++;
			}
		}
	} else if (srcGR) {
		assert(m_palette);
		int inkId = m_palette->addStyle(m_configuration.m_inkColor);

		m_dataRasterArray.push_back(pair<int, DataRasterP>(inkId, dataRaster));
		int threshold = m_configuration.m_threshold;

		for (y = 0; y < m_src->getLy(); y++) {
			TPixelGR8 *inPix = srcGR->pixels(y);
			TPixelGR8 *inEndPix = inPix + srcGR->getLx();
			DataPixel *dataPix = dataRaster->pixels(y + 1) + 1;
			x = 0;
			while (inPix < inEndPix) {
				*dataPix = DataPixel();
				//int distance2 = colorDistance2(m_configuration.m_inkColor, *inPix);

				//int value = (inPix->r + 2*inPix->g + inPix->b)>>2;
				if (y == 0 || y == m_src->getLy() - 1 || x == 0 || x == m_src->getLx() - 1) {
					dataPix->m_value = 255;
					dataPix->m_ink = false;
				} else {
					dataPix->m_value = inPix->value;
					dataPix->m_ink = (inPix->value < threshold);
				}
				dataPix->m_pos.x = x++;
				dataPix->m_pos.y = y;
				dataPix->m_node = 0;
				//dataPix->m_index = index++;
				inPix++;
				dataPix++;
			}
		}
	}

	else if (srcCM) {
		int currInk, nextInk = 0;
		do {
			int threshold = m_configuration.m_threshold; //tolerance: 1->MAX thresh: 1-255
			currInk = nextInk;
			nextInk = 0;
			m_dataRasterArray.push_back(pair<int, DataRasterP>(currInk, dataRaster));

			// inizializza la parte centrale
			for (y = 0; y < srcCM->getLy(); y++) {
				TPixelCM32 *inPix = srcCM->pixels(y);
				TPixelCM32 *inEndPix = inPix + m_src->getLx();
				DataPixel *dataPix = dataRaster->pixels(y + 1) + 1;
				x = 0;
				while (inPix < inEndPix) {
					*dataPix = DataPixel();
					int value = inPix->getTone();
					if (value < 255 && !m_configuration.m_ignoreInkColors) {
						int ink = inPix->getInk();
						if (currInk == 0) {
							currInk = ink;
							m_dataRasterArray.back().first = ink;
						} else if (ink != currInk) {
							value = 255;
							if (nextInk == 0) {
								for (ii = 0; ii < (int)m_dataRasterArray.size() - 1; ii++)
									if (m_dataRasterArray[ii].first == ink)
										break;
								if (ii == (int)m_dataRasterArray.size() - 1)
									nextInk = ink;
							}
						}
					}
					//TPixel col = m_palette->getStyle(inPix->getInk())->getMainColor();
					dataPix->m_pos.x = x++;
					dataPix->m_pos.y = y;
					dataPix->m_value = value;
					dataPix->m_ink = (value < threshold);
					dataPix->m_node = 0;
					//dataPix->m_index = index++;
					inPix++;
					dataPix++;
				}
			}
		} while (nextInk != 0);
		if (m_configuration.m_ignoreInkColors) {
			assert(m_dataRasterArray.size() == 1);
			m_dataRasterArray.back().first = 1;
		}
	} else
		assert(false);
}

void CenterLineVectorizer::init()
{
	int y;

	/*m_links.clear();
	m_protoStrokes.clear();
	m_protoRegions.clear();
	clearNodes();
	clearJunctions();*/

	DataRasterP dataRaster = m_dataRaster;
	const int wrap = dataRaster->getWrap();
	const int delta[] = {-wrap - 1, -wrap, -wrap + 1, 1, wrap + 1, wrap, wrap - 1, -1};

	for (y = 1; y < dataRaster->getLy() - 1; y++) {
		DataPixel *pix = dataRaster->pixels(y);
		DataPixel *endPix = pix + dataRaster->getLx() - 1;
		pix++;
		for (pix++; pix < endPix; ++pix) {
			if ((pix->m_ink == false) ||
				(pix[-wrap].m_ink && pix[wrap].m_ink &&
				 pix[-1].m_ink && pix[1].m_ink))
				continue;
			int i;
			for (i = 0; i < 8; i++)
				if (pix[delta[i]].m_ink && pix[delta[(i + 1) & 0x7]].m_ink == false)
					break;
			int start = i;
			if (i == 8)
				continue; // punto isolato
			for (;;) {
				int j = (i + 1) & 0x7;
				assert(i < 8 && pix[delta[i]].m_ink);
				assert(j < 8 && pix[delta[j]].m_ink == false);
				do
					j = (j + 1) & 0x7;
				while (pix[delta[j]].m_ink == false);
				assert(j < 8 && pix[delta[j]].m_ink);
				if (((i + 2) & 0x7) != j || (i & 1) == 0) {
					// il bianco comprende anche un fianco
					link(pix, pix + delta[i], pix + delta[j]);
				}
				i = j;
				assert(i < 8);
				while (pix[delta[(i + 1) & 0x7]].m_ink)
					i = (i + 1) & 0x7;
				assert(i < 8 && pix[delta[i]].m_ink);
				assert(pix[delta[(i + 1) & 0x7]].m_ink == false);
				if (i == start)
					break;
			}
		}
	}
}
/*
#ifdef DEBUG
    
    int i;
    for(i = 0;i<(int)m_nodes.size();i++)
	{
		Node *node = m_nodes[i];
		if(node->m_pixel==0) continue;
		assert(node->m_prev);
		assert(node->m_next);
		assert(node->m_prev->m_next == node);
		assert(node->m_next->m_prev == node);
	}

    for(i = 0;i<m_dataRaster->getLx()*m_dataRaster->getLy(); i++)
	{
		DataPixel *pix = m_dataRaster->pixels()+i;
		for(Node *node = pix->m_node; node; node=node->m_other)
		{
			if(node->m_pixel==0) continue;
			assert(node->m_pixel == pix);
			assert(node->m_prev);
			assert(node->m_next);
			assert(node->m_prev->m_next == node);
			assert(node->m_next->m_prev == node);
		}
	}
#endif
*/

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

TPointD CenterLineVectorizer::computeCenter(Node *na, Node *nb, double &r)
{
	TPointD pa = convert(na->m_pixel->m_pos);
	TPointD pb = convert(nb->m_pixel->m_pos);
	r = THICKNESS_LIMIT;
	if (pa == pb)
		return pa;

	TPointD center((pa + pb) * 0.5);
	r = norm(pa - pb) * 0.5;
	return center;

	/*	double vsum = 0.0;
	TPointD d = pa-pb;
	if(fabs(d.x)>fabs(d.y))
    {
		if(pa.x>pb.x) tswap(pa,pb);
		for(int x = pa.x; x<=pb.x; x++)
		{
			int y = pa.y + (pb.y-pa.y)*(x-pa.x)/(pb.x-pa.x);
			int v = 255-m_dataRaster->pixels(y)[x].m_value;
			center += v * TPointD(x,y);
			vsum += v;
        }
    }
	else
    {
		if(pa.y>pb.y) tswap(pa,pb);
		for(int y = pa.y; y<=pb.y; y++)
		{
			int x = pa.x + (pb.x-pa.x)*(y-pa.y)/(pb.y-pa.y);
			int v = 255-m_dataRaster->pixels(y)[x].m_value;
			center += v * TPointD(x,y);
			vsum += v;
		}
    }
	assert(vsum>0);
	r = 0.5*vsum/255.0;
	return center * (1/vsum);
*/
}

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

double computeDistance2(Node *na, Node *nb)
{
	assert(na->m_pixel);
	assert(nb->m_pixel);
	TPointD d = convert(na->m_pixel->m_pos - nb->m_pixel->m_pos);
	return d * d;
}

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

Node *CenterLineVectorizer::findOtherSide(Node *node)
{
	DataPixel *pix = node->m_pixel;

	TPoint dir = -computeGradient(pix);
	if (dir == TPoint(0, 0))
		return 0;
	TPoint d1(tsign(dir.x), 0), d2(0, tsign(dir.y));
	int num = abs(dir.y), den = abs(dir.x);
	if (num > den) {
		tswap(d1, d2);
		tswap(num, den);
	}
	TPoint pos = pix->m_pos;
	int i;
	for (i = 0;; i++) {
		TPoint q(pos.x + d1.x * i + d2.x * num * i / den, pos.y + d1.y * i +
															  d2.y * num * i / den);
		//if (!m_dataRaster->getBounds().contains(q)) break;
		DataPixel *nextPix = m_dataRaster->pixels(q.y) + q.x;
		if (nextPix->m_ink == false)
			break;
		pix = nextPix;
	}
	assert(pix);
	if (!pix->m_node) {
		const int wrap = m_dataRaster->getWrap();
		if (pix[-1].m_node)
			pix--;
		else if (pix[1].m_node)
			pix++;
		else if (pix[wrap].m_node)
			pix += wrap;
		else if (pix[-wrap].m_node)
			pix -= wrap;
		else {
			assert(0);
		}
	}
	if (!pix->m_node)
		return 0;
	Node *q = pix->m_node;
	while (q->m_pixel == 0 && q->m_other)
		q = q->m_other;
	assert(q && q->m_pixel == pix);

	for (i = 0; i < 5; i++) {
		if (!q->m_prev)
			break;
		q = q->m_prev;
	}

	Node *best = q;
	double bestDist2 = computeDistance2(q, node);
	for (i = 0; i < 10; i++) {
		q = q->m_next;
		if (!q)
			break;
		double dist2 = computeDistance2(q, node);
		if (dist2 < bestDist2) {
			bestDist2 = dist2;
			best = q;
		}
	}

	return best;
}

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

TPointD findDirection(Node *na, Node *nb, int skip)
{
	for (int i = 0; i < skip; i++) {
		na = na->m_next;
		nb = nb->m_prev;
	}
	TPointD p0 = 0.5 * (convert(na->m_pixel->m_pos) + convert(nb->m_pixel->m_pos));
	na = na->m_next;
	nb = nb->m_prev;
	TPointD p1 = 0.5 * (convert(na->m_pixel->m_pos) + convert(nb->m_pixel->m_pos));
	na = na->m_next;
	nb = nb->m_prev;
	TPointD p2 = 0.5 * (convert(na->m_pixel->m_pos) + convert(nb->m_pixel->m_pos));
	TPointD dir = 0.25 * (3 * (p1 - p0) + (p2 - p1));
	double dirn = norm(dir);
	if (dirn > 0)
		return dir * (1 / dirn);
	else
		return dir;
}

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

void setEndpointsDirections(ProtoStroke &protoStroke)
{
	assert(!protoStroke.m_points.empty());
	int size = protoStroke.m_points.size();
	double startThickness = 0, endThickness = 0;
	if (size == 1)
		return;
	if (size < 5) {
		protoStroke.m_startDirection = protoStroke.m_points[0] -
									   protoStroke.m_points[1];
		protoStroke.m_endDirection = protoStroke.m_points[size - 1] -
									 protoStroke.m_points[size - 2];
	} else if (size < 10) {
		assert(protoStroke.m_points.front() != protoStroke.m_points.back());
		protoStroke.m_startDirection = protoStroke.m_points[0] -
									   protoStroke.m_points[2];
		protoStroke.m_endDirection = protoStroke.m_points[size - 1] -
									 protoStroke.m_points[size - 3];
		startThickness = 0.5 * (protoStroke.m_points[1].thick +
								protoStroke.m_points[2].thick);
		endThickness = 0.5 * (protoStroke.m_points[size - 2].thick +
							  protoStroke.m_points[size - 3].thick);
	} else {
		protoStroke.m_startDirection = 4 * protoStroke.m_points[0];
		protoStroke.m_endDirection = 4 * protoStroke.m_points[size - 1];

		for (int i = 0; i < 4; ++i) {
			protoStroke.m_startDirection -= protoStroke.m_points[3 + i];
			protoStroke.m_endDirection -= protoStroke.m_points[size - 4 - i];
			startThickness += protoStroke.m_points[i].thick;
			endThickness += protoStroke.m_points[size - i - 1].thick;
		}
		protoStroke.m_startDirection = protoStroke.m_startDirection * (1 / 4.0);
		protoStroke.m_endDirection = protoStroke.m_endDirection * (1 / 4.0);
		startThickness *= (1 / 4.0);
		endThickness *= (1 / 4.0);
	}

	double startNorm = norm(protoStroke.m_startDirection);
	if (startNorm)
		protoStroke.m_startDirection = protoStroke.m_startDirection *
									   (1 / startNorm);

	double endNorm = norm(protoStroke.m_endDirection);
	if (endNorm)
		protoStroke.m_endDirection = protoStroke.m_endDirection * (1 / endNorm);

	if (startThickness != 0)
		protoStroke.m_points[0].thick = protoStroke.m_points[1].thick =
			protoStroke.m_points[2].thick = startThickness;
	if (endThickness != 0)
		protoStroke.m_points[size - 1].thick = protoStroke.m_points[size - 2].thick =
			protoStroke.m_points[size - 3].thick = endThickness;
}

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

bool CenterLineVectorizer::testOtherSide(Node *na1, Node *nb1, double
																   &startDist2)
{
	if (na1->m_pixel->m_pos == nb1->m_pixel->m_pos || na1->m_next == nb1 || na1->m_prev == nb1 || na1->m_next->m_next == nb1 || na1->m_prev->m_prev == nb1)
		return false;

	Node *na2 = na1->m_next;
	for (int i = 1; i < 3; ++i) {
		na2 = na2->m_next;
		Node *nb2 = findOtherSide(na2);
		if (nb2 == 0 && i == 3)
			return true;
		if (nb2 == 0)
			continue;

		double newStartDist2 = computeDistance2(na2, nb2);
		if (newStartDist2 > startDist2)
			startDist2 = newStartDist2;
		for (int j = 0; j < 3 * i + 1; ++j) {
			if (nb2 == nb1)
				return true;
			nb2 = nb2->m_next;
		}
	}
	return false;
}

//---------------------------------------------------------
/*
bool CenterLineVectorizer::isHole(Node *startNode)
{
	Node *node = startNode;
	for(int i=0; i<20; ++i)
	{
		node = node->m_next;
		if(node == startNode)
		{
			std::vector<TThickPoint> points;
			do
			{
				node = node->m_next;
				node->m_visited = true;
				node->m_flag = true;
				points.push_back(convert(node->m_pixel->m_pos));
			}
			while(node != startNode);
			m_protoHoles.push_back(points);
			return true;
		}
	}
	return false;
}
*/

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

bool isJunction(const Node *na, const Node *nb)
{
	TPointD pa = convert(na->m_pixel->m_pos);
	TPointD pb = convert(nb->m_pixel->m_pos);

	Node *na0 = na->m_prev;
	Node *na1 = na->m_next;
	Node *nb0 = nb->m_next;
	Node *nb1 = nb->m_prev;

	TPointD da = convert(na1->m_pixel->m_pos) - pa;
	TPointD db = convert(nb1->m_pixel->m_pos) - pb;

	TPointD da0 = pa - convert(na0->m_pixel->m_pos);
	TPointD db0 = pb - convert(nb0->m_pixel->m_pos);

	if ((da0.x * da.x + da0.y * da.y) / (norm(da0) * norm(da)) <= 0) {
		return true;
	}

	if ((db0.x * db.x + db0.y * db.y) / (norm(db0) * norm(db)) <= 0) {
		return true;
	}

	if (right(TPointD(0, 0), da, db)) {
		if ((da.x * db.x + da.y * db.y) / (norm(da) * norm(db)) <= 0) {
			return true;
		}
	} else {
		if ((da.x * db.x + da.y * db.y) / (norm(da) * norm(db)) <= -0.5) {
			return true;
		}
	}
	return false;
}

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

void CenterLineVectorizer::traceLine(DataPixel *pix)
{
	Junction *startJunction = 0;
	Junction *endJunction = 0;
	assert(m_dataRaster);
	const int wrap = m_dataRaster->getWrap();
	//  TRect bounds = m_dataRaster->getBounds();
	if (!pix->m_ink)
		return;
	while (pix->m_node == 0)
		pix -= wrap;

	Node *node = pix->m_node;
	while (node && node->m_pixel == 0 && node->m_other)
		node = node->m_other;
	assert(node && node->m_pixel == pix);

	Node *naa = node;

	bool isSamePixel;
	if (node->m_other && node->m_other->m_pixel && (node->m_pixel->m_pos ==
													node->m_other->m_pixel->m_pos)) {
		node = node->m_other;
		isSamePixel = true;
	} else {
		node = findOtherSide(node);
		isSamePixel = false;
	}
	if (node == 0)
		return;
	Node *nbb = node;

	if (naa->m_visited || nbb->m_visited)
		return;

	double startDist2 = computeDistance2(naa, nbb);
	if (startDist2 > 2000)
		return;
	if (startDist2 < THICKNESS_LIMIT_2)
		startDist2 = THICKNESS_LIMIT_2;

	//try to detect errors in findOtherSide
	if (!isSamePixel && !testOtherSide(naa, nbb, startDist2))
		return;

	double minDist2 = 0;
	double maxDist2 = startDist2 + 3.0;
	if (startDist2 > 36) {
		minDist2 = startDist2 / 2;
		maxDist2 = startDist2 + 0.25 * startDist2;
	}
	if (naa->m_next == nbb || naa->m_prev == nbb)
		return;

	naa->m_visited = nbb->m_visited = true;
	std::deque<TThickPoint> dpoints;

	for (int k = 0; k < 2; k++) {
		Node *na = naa;
		Node *nb = nbb;

		for (;;) {
			TPoint &pa = na->m_pixel->m_pos;
			TPoint &pb = nb->m_pixel->m_pos;
			TPoint ba = pb - pa;

			Node *na1 = na->m_next;
			Node *nb1 = nb->m_prev;
			//if (!na1 || !nb1) break;

			TPoint da = na1->m_pixel->m_pos - pb;
			TPoint db = nb1->m_pixel->m_pos - pa;

			if (na1->m_junction || nb1->m_junction) {
				na1->m_visited = true;
				nb1->m_visited = true;
				TPointD p1;
				if (na1->m_junction && nb1->m_junction) {
					//					assert(na1->m_junction == nb1->m_junction);
				} else {
					if (na1->m_junction) {
						p1 = convert(nb1->m_pixel->m_pos);
						nb1->m_junction = na1->m_junction;
						na1->m_junction->m_nodes.push_front(nb1);
					} else {
						p1 = convert(na1->m_pixel->m_pos);
						na1->m_junction = nb1->m_junction;
						nb1->m_junction->m_nodes.push_back(na1);
					}
				}
				double r;
				TPointD center;
				TThickPoint point;

				center = computeCenter(na, nb, r);
				if (r < THICKNESS_LIMIT)
					r = THICKNESS_LIMIT;
				point = TThickPoint(center, r);
				if (k == 0)
					dpoints.push_front(point);
				else
					dpoints.push_back(point);

				center = computeCenter(na1, nb1, r);
				if (r < THICKNESS_LIMIT)
					r = THICKNESS_LIMIT;
				point = TThickPoint(center, r);
				if (k == 0)
					dpoints.push_front(point);
				else
					dpoints.push_back(point);

				if (k == 0)
					startJunction = na1->m_junction;
				else
					endJunction = na1->m_junction;

				break;
			}

			//			if(da*ba>-db*ba)
			if (norm2(db) >= norm2(da)) {
				if (na1->m_visited)
					break;

				na = na1;
				na->m_visited = true;
#ifdef DEBUG
				getOutPix(na->m_pixel) = TPixel32::Red;
#endif
			} else {
				if (nb1->m_visited)
					break;

				nb = nb1;
				nb->m_visited = true;
#ifdef DEBUG
				getOutPix(nb->m_pixel) = TPixel32::Green;
#endif
			}
			double distance2 = computeDistance2(na, nb);
			if (distance2 < 0.25)
				distance2 = 0.25;

			double r;
			TPointD center = computeCenter(na, nb, r);
			if (r < THICKNESS_LIMIT)
				r = THICKNESS_LIMIT;
			TThickPoint point(center, r);
			TPointI pos = convert(center);
			//TPixel32  *inPix = m_src->pixels(pos.y);
			//inPix += pos.x;
			//int value = (inPix->r + 2*inPix->g + inPix->b)>>2;

			if (distance2 > maxDist2 || distance2 < minDist2 || na->m_next == nb || na->m_prev == nb
				//|| isJunction(na, nb)
				//|| !(value<200)
				) {
				if (na->m_next == nb || na->m_prev == nb || na->m_next->m_next == nb ||
					na->m_prev->m_prev == nb) {
					if (k == 0) {
						if (!dpoints.empty())
							dpoints.pop_front();
						if (!dpoints.empty())
							dpoints.pop_front();
					} else {
						if (!dpoints.empty())
							dpoints.pop_back();
						if (!dpoints.empty())
							dpoints.pop_back();
					}
				}
				//				if(k==0) dpoints.push_front(point);
				//				else dpoints.push_back(point);
				na->m_junction = nb->m_junction = new Junction();
				na->m_junction->m_nodes.push_back(nb);
				nb->m_junction->m_nodes.push_back(na);
				m_junctions.push_back(na->m_junction);

				if (k == 0)
					startJunction = na->m_junction;
				else
					endJunction = na->m_junction;
				break;
			}
			if (k == 0)
				dpoints.push_front(point);
			else
				dpoints.push_back(point);
		}

		tswap(naa, nbb);
	}
	if (dpoints.size() == 0)
		naa->m_visited = nbb->m_visited = false;
	else {
		m_protoStrokes.push_back(ProtoStroke(dpoints.begin(), dpoints.end()));
		setEndpointsDirections(m_protoStrokes.back());
		if (startJunction) {
			m_protoStrokes.back().m_startJunction = startJunction;
			startJunction->m_protoStrokes.push_back(&m_protoStrokes.back());
		}
		if (endJunction) {
			m_protoStrokes.back().m_endJunction = endJunction;
			endJunction->m_protoStrokes.push_back(&m_protoStrokes.back());
		}
	}
}

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

Node *CenterLineVectorizer::createNode(DataPixel *pix)
{
	Node *node = new Node();
	node->m_pixel = pix;
	node->m_other = pix->m_node;
	pix->m_node = node;
	m_nodes.push_back(node);
	return node;
}

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

void CenterLineVectorizer::createNodes()
{
}

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

#ifdef DEBUG
void CenterLineVectorizer::outputNodes(Node *node)
{
	std::vector<TThickPoint> points;
	Node *start = node, *curr = node;
	int i = 0;
	do {
		TThickPoint point(convert(curr->m_pixel->m_pos), 0);
		points.push_back(point);

		curr->m_flag = true;
		assert(curr->m_next);
		curr = curr->m_next;
		i++;
	} while (curr != start);
	TThickPoint point(convert(curr->m_pixel->m_pos), 0);
	points.push_back(point);
	m_outlines.push_back(points);
}
#endif

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

bool CenterLineVectorizer::linkNextProtoStroke(
	ProtoStroke *const dstProtoStroke,
	Junction *const currJunction,
	int k)
{
	assert(dstProtoStroke);
	assert(currJunction);
	//	if(dstProtoStroke->m_startJunction == dstProtoStroke->m_endJunction)
	return false;
	if (k == 0 && dstProtoStroke->m_startJunction->m_locked ||
		k == 1 && dstProtoStroke->m_endJunction->m_locked)
		return false;

	ProtoStroke *srcProtoStroke = 0;
	Junction *nextJunction = 0;

	std::vector<ProtoStroke *>::iterator it;
	std::vector<ProtoStroke *>::iterator it_e =
		currJunction->m_protoStrokes.end();
	std::vector<ProtoStroke *>::iterator it_candidate = it_e;
	double candidateProbability = JOIN_LIMIT;
	int candidate_k;
	for (it = currJunction->m_protoStrokes.begin(); it != it_e; it++) {
		//erase dstProtoStroke in currJunction
		if (*it == dstProtoStroke) {
			it = currJunction->m_protoStrokes.erase(it);
			it_e = currJunction->m_protoStrokes.end();
			if (candidateProbability == JOIN_LIMIT)
				it_candidate = it_e;
			if (it == it_e)
				break;
		}
		double probability;
		int next_k;
		if (((*it))->m_startJunction == currJunction)
			next_k = 0;
		else if ((*it)->m_endJunction == currJunction)
			next_k = 1;
		else
			assert(false);

		if (k == 0) {
			if (next_k == 0)
				probability = 1 - norm2(dstProtoStroke->m_startDirection +
										(*it)->m_startDirection);
			else
				probability = 1 - norm2(dstProtoStroke->m_startDirection +
										(*it)->m_endDirection);
		} else {
			if (next_k == 0)
				probability = 1 - norm2(dstProtoStroke->m_endDirection +
										(*it)->m_startDirection);
			else
				probability = 1 - norm2(dstProtoStroke->m_endDirection +
										(*it)->m_endDirection);
		}

		if (probability > candidateProbability ||
			//			currJunction->m_junctionOrder == 2 ||
			(*it)->m_points.size() < 3) {
			it_candidate = it;
			candidateProbability = probability;
			candidate_k = next_k;
		}
	}

	if (it_candidate == it_e)
		return false;
	srcProtoStroke = *it_candidate;
	assert(srcProtoStroke);

	if ((*it_candidate)->m_startJunction == (*it_candidate)->m_endJunction)
		return false;
	/*	if(k==0)
	{
		if(	candidate_k == 0 && dstProtoStroke->m_startJunction == 
(*it_candidate)->m_endJunction ||
			candidate_k == 1 && dstProtoStroke->m_startJunction == 
(*it_candidate)->m_startJunction	)
			return false;
	}
	else if(k==1)
	{
		if(	candidate_k == 0 && dstProtoStroke->m_endJunction == 
(*it_candidate)->m_endJunction ||
			candidate_k == 1 && dstProtoStroke->m_endJunction == 
(*it_candidate)->m_startJunction	)
			return false;
	}
*/
	if (candidate_k == 0)
		nextJunction = (*it_candidate)->m_endJunction;
	else
		nextJunction = (*it_candidate)->m_startJunction;

	if (k == 0) {
		if (candidate_k == 0) {
			dstProtoStroke->m_points.insert(
				dstProtoStroke->m_points.begin(),
				srcProtoStroke->m_points.rbegin(),
				srcProtoStroke->m_points.rend());
		} else {
			dstProtoStroke->m_points.insert(
				dstProtoStroke->m_points.begin(),
				srcProtoStroke->m_points.begin(),
				srcProtoStroke->m_points.end());
		}
	} else {
		if (candidate_k == 0) {
			dstProtoStroke->m_points.insert(
				dstProtoStroke->m_points.end(),
				srcProtoStroke->m_points.begin(),
				srcProtoStroke->m_points.end());
		} else {
			dstProtoStroke->m_points.insert(
				dstProtoStroke->m_points.end(),
				srcProtoStroke->m_points.rbegin(),
				srcProtoStroke->m_points.rend());
		}
	}

	srcProtoStroke->m_points.clear();
	srcProtoStroke->m_startJunction = 0;
	srcProtoStroke->m_endJunction = 0;

	//refresh currJunction
	currJunction->m_protoStrokes.erase(it_candidate);

	//refresh nextJunction
	if (nextJunction) {
		it = std::find(
			nextJunction->m_protoStrokes.begin(),
			nextJunction->m_protoStrokes.end(),
			srcProtoStroke);
		assert(it != nextJunction->m_protoStrokes.end());
		nextJunction->m_protoStrokes.erase(it);
		nextJunction->m_protoStrokes.push_back(dstProtoStroke);
	}

	if (k == 0) {
		dstProtoStroke->m_startJunction = nextJunction;

		if (candidate_k == 0)
			dstProtoStroke->m_startDirection = srcProtoStroke->m_endDirection;
		else
			dstProtoStroke->m_startDirection = srcProtoStroke->m_startDirection;
	} else {
		dstProtoStroke->m_endJunction = nextJunction;

		if (candidate_k == 0)
			dstProtoStroke->m_endDirection = srcProtoStroke->m_endDirection;
		else
			dstProtoStroke->m_endDirection = srcProtoStroke->m_startDirection;
	}
	if (!nextJunction)
		return false;
	return true;
}

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

void CenterLineVectorizer::joinProtoStrokes(ProtoStroke *const
												dstProtoStroke)
{
	//std::deque<TThickPoint> &dstPoints = dstProtoStroke->m_points;
	for (int k = 0; k < 2; k++) {
		for (;;) {
			Junction *currJunction;
			if (k == 0 && !dstProtoStroke->m_startJunction)
				break;
			else if (k == 0)
				currJunction = dstProtoStroke->m_startJunction;
			else if (!dstProtoStroke->m_endJunction)
				break;
			else
				currJunction = dstProtoStroke->m_endJunction;

			if (!linkNextProtoStroke(dstProtoStroke, currJunction, k))
				break;
		}
	}
	std::vector<ProtoStroke *>::iterator it;
	if (dstProtoStroke->m_startJunction) {
		it = std::find(
			dstProtoStroke->m_startJunction->m_protoStrokes.begin(),
			dstProtoStroke->m_startJunction->m_protoStrokes.end(),
			dstProtoStroke);
		if (it != dstProtoStroke->m_startJunction->m_protoStrokes.end())
			dstProtoStroke->m_startJunction->m_protoStrokes.erase(it);
		dstProtoStroke->m_startJunction = 0;
	}
	if (dstProtoStroke->m_endJunction) {
		it = std::find(
			dstProtoStroke->m_endJunction->m_protoStrokes.begin(),
			dstProtoStroke->m_endJunction->m_protoStrokes.end(),
			dstProtoStroke);
		if (it != dstProtoStroke->m_endJunction->m_protoStrokes.end())
			dstProtoStroke->m_endJunction->m_protoStrokes.erase(it);
		dstProtoStroke->m_endJunction = 0;
	}
}

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

bool y_intersect(const double y0, const TPointD p1, const TPointD p2,
				 double &x0)
{
	if ((y0 < p1.y && y0 < p2.y) ||
		(y0 > p1.y && y0 > p2.y))
		return false;

	x0 = (p1.x * (p2.y - y0) + p2.x * (y0 - p1.y)) / (p2.y - p1.y);
	return true;
}

bool isContained(TThickPoint point, std::deque<Node *> &polygon)
{
	std::deque<Node *>::iterator it_b = polygon.begin();
	std::deque<Node *>::iterator it_e = polygon.end();
	std::deque<Node *>::iterator it;

	int intersectionsNumber = 0;

	for (it = it_b; it < it_e; ++it) {
		double x;
		TPointD p1, p2;

		p1 = convert((*it)->m_pixel->m_pos);
		if ((it + 1) == it_e)
			p2 = convert((*it_b)->m_pixel->m_pos);
		else
			p2 = convert((*(it + 1))->m_pixel->m_pos);

		if (point.y == p1.y || point.y == p2.y || p1.y == p2.y)
			continue;
		if (y_intersect(point.y, p1, p2, x) && x >= point.x)
			++intersectionsNumber;
	}

	return (intersectionsNumber % 2 == 0) ? false : true;
}

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

Junction *CenterLineVectorizer::mergeJunctions(const std::list<JunctionMerge> &junctions)
{
	Junction *newJunction = new Junction();

	std::list<JunctionMerge>::const_iterator it_junction;
	std::deque<Node *>::iterator it_startNode, it_endNode;

	for (it_junction = junctions.begin(); it_junction != junctions.end();
		 ++it_junction) {
		newJunction->m_center += it_junction->m_junction->m_center;

		std::deque<Node *>::iterator it_node;
		if (it_junction->m_startNode == 0 && it_junction->m_endNode == 0) {
			for (it_node = it_junction->m_junction->m_nodes.begin();
				 it_node != it_junction->m_junction->m_nodes.end(); ++it_node) {
				(*it_node)->m_junction = newJunction;
				newJunction->m_nodes.push_back(*it_node);
			}
		} else {
			//			if( isContained(it_junction->m_junction->m_center, newJunction->m_nodes) ) continue;
			/*			if(	it_junction->m_junction->m_nodes.size() < 3)
			{
				std::vector<ProtoStroke*>::iterator it_protoStroke;
				bool ignoreJunction = false;
				for(it_protoStroke = it_junction->m_junction->m_protoStrokes.begin();
					it_protoStroke!=it_junction->m_junction->m_protoStrokes.end();
					++it_protoStroke)
				{
					if(	(*it_protoStroke)->m_startJunction == newJunction ||
						(*it_protoStroke)->m_endJunction == newJunction )
					{
						ignoreJunction = true;
						break;
					}
				}
				if(ignoreJunction) continue;
			}
*/

			it_startNode = std::find(newJunction->m_nodes.begin(),
									 newJunction->m_nodes.end(),
									 it_junction->m_startNode);
			assert(it_startNode != newJunction->m_nodes.end());
			it_endNode = std::find(it_junction->m_junction->m_nodes.begin(),
								   it_junction->m_junction->m_nodes.end(),
								   it_junction->m_endNode);
			assert(it_endNode != it_junction->m_junction->m_nodes.end());

			for (it_node = it_junction->m_junction->m_nodes.begin();
				 it_node != it_junction->m_junction->m_nodes.end(); ++it_node) {
				(*it_node)->m_junction = newJunction;
			}

			if (it_junction->m_isNext) {
				++it_startNode;
				newJunction->m_nodes.insert(it_startNode,
											it_junction->m_junction->m_nodes.begin(), it_endNode);
				it_startNode = std::find(newJunction->m_nodes.begin(),
										 newJunction->m_nodes.end(),
										 it_junction->m_startNode);
				assert(it_startNode != newJunction->m_nodes.end());
				++it_startNode;
				newJunction->m_nodes.insert(it_startNode, it_endNode,
											it_junction->m_junction->m_nodes.end());
			} else {
				++it_endNode;
				newJunction->m_nodes.insert(it_startNode,
											it_junction->m_junction->m_nodes.begin(), it_endNode);
				it_startNode = std::find(newJunction->m_nodes.begin(),
										 newJunction->m_nodes.end(),
										 it_junction->m_startNode);
				assert(it_startNode != newJunction->m_nodes.end());
				newJunction->m_nodes.insert(it_startNode, it_endNode,
											it_junction->m_junction->m_nodes.end());
			}
		}

		std::vector<ProtoStroke *>::iterator it_protoStroke;
		for (it_protoStroke = it_junction->m_junction->m_protoStrokes.begin();
			 it_protoStroke != it_junction->m_junction->m_protoStrokes.end();
			 ++it_protoStroke) {
			if ((*it_protoStroke)->m_startJunction == it_junction->m_junction) {
				(*it_protoStroke)->m_startJunction = newJunction;
				newJunction->m_protoStrokes.push_back(*it_protoStroke);
			}
			if ((*it_protoStroke)->m_endJunction == it_junction->m_junction) {
				(*it_protoStroke)->m_endJunction = newJunction;
				newJunction->m_protoStrokes.push_back(*it_protoStroke);
			}
		}

		std::vector<JunctionLink>::iterator it_links;
		for (it_links = m_links.begin(); it_links != m_links.end(); ++it_links) {
			if (it_links->m_first == it_junction->m_junction)
				it_links->m_first =
					newJunction;
			if (it_links->m_second == it_junction->m_junction)
				it_links->m_second =
					newJunction;
		}

		//erase old junction
		std::vector<Junction *>::iterator it;
		it = std::find(
			m_junctions.begin(),
			m_junctions.end(),
			it_junction->m_junction);
		assert(it != m_junctions.end());
		delete it_junction->m_junction;
		m_junctions.erase(it);
	}
	double size = junctions.size();
	newJunction->m_center = newJunction->m_center * (1 / size);

	std::vector<ProtoStroke *>::iterator it_protoStroke;
	for (it_protoStroke = newJunction->m_protoStrokes.begin();
		 it_protoStroke != newJunction->m_protoStrokes.end();
		 ++it_protoStroke) {
		if ((*it_protoStroke)->m_points.size() < 6 &&
			(*it_protoStroke)->m_startJunction == newJunction &&
			(*it_protoStroke)->m_endJunction == newJunction) {
			(*it_protoStroke)->m_points.clear();
			(*it_protoStroke)->m_startJunction = (*it_protoStroke)->m_endJunction =
				0;

			ProtoStroke *protoStroke = *it_protoStroke;
			for (;;) {
				std::vector<ProtoStroke *>::iterator it;
				it = std::find(newJunction->m_protoStrokes.begin(),
							   newJunction->m_protoStrokes.end(),
							   protoStroke);
				if (it == newJunction->m_protoStrokes.end())
					break;
				newJunction->m_protoStrokes.erase(it);
			}
			if (newJunction->m_protoStrokes.empty())
				break;
			it_protoStroke = newJunction->m_protoStrokes.begin();
		}
	}

	return newJunction;
}

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

void CenterLineVectorizer::joinJunctions()
{
	for (int i = 0; i < (int)m_junctions.size(); i++) {
		Junction *currJunction = m_junctions[i];

		//		if(currJunction->m_nodes.size() <= 2) continue;

		std::list<JunctionMerge> joinJunctions;
		joinJunctions.push_back(JunctionMerge(currJunction));

		//		std::sort(currJunction->m_nodes.begin(), currJunction->m_nodes.end());
		std::deque<Node *>::iterator it;
		for (it = currJunction->m_nodes.begin(); it != currJunction->m_nodes.end();
			 it++) {
			Node *currNode = *it;
			Node *prevNode = *it;
			Node *nextNode = *it;
			std::vector<JunctionLink>::iterator it_link;
			int j;
			for (j = 0; true; ++j) {
				prevNode = prevNode->m_prev;
				prevNode->m_visited = true;
				if (prevNode->m_junction != 0 && prevNode->m_junction != currJunction) {
					//double junctionMaxDist = 3.5*currJunction->m_center.thick;
					it_link = std::find(m_links.begin(), m_links.end(),
										JunctionLink(currJunction, prevNode->m_junction));
					/*					if(	j >= 1*SEARCH_WINDOW &&
						norm(currJunction->m_center - prevNode->m_junction->m_center) >= 
junctionMaxDist)
					{
						assert(false);
						if(it_link == m_links.end())
							m_links.push_back(JunctionLink(currJunction, prevNode->m_junction));
						else
							++(it_link->m_order);
					}
					else*/ if (it_link == m_links.end() &&
																		std::find(joinJunctions.begin(), joinJunctions.end(),
																				  JunctionMerge(prevNode->m_junction)) == joinJunctions.end()) {
						joinJunctions.push_back(JunctionMerge(prevNode->m_junction, currNode,
															  prevNode, false));
					}
				}
				if (prevNode->m_visited)
					break;
			}
			for (j = 0; true; ++j) {
				nextNode = nextNode->m_next;
				prevNode->m_visited = true;
				if (nextNode->m_junction != 0 && nextNode->m_junction != currJunction) {
					//double junctionMaxDist = 3.5*currJunction->m_center.thick;
					it_link = std::find(m_links.begin(), m_links.end(),
										JunctionLink(currJunction, nextNode->m_junction));
					/*					if(	j >= 1*SEARCH_WINDOW &&
						norm(currJunction->m_center - nextNode->m_junction->m_center) >= 
junctionMaxDist )
					{
						assert(false);
						if(it_link == m_links.end())
							m_links.push_back(JunctionLink(currJunction, nextNode->m_junction));
						else
							++(it_link->m_order);
					}
					else*/ if (it_link == m_links.end() &&
																		std::find(joinJunctions.begin(), joinJunctions.end(),
																				  JunctionMerge(nextNode->m_junction)) == joinJunctions.end()) {
						joinJunctions.push_back(JunctionMerge(nextNode->m_junction, currNode,
															  nextNode, true));
					}
				}
				if (nextNode->m_visited)
					break;
			}
		}
		if (joinJunctions.size() > 1) {
			//erase duplicate junctions
			//			joinJunctions.sort();
			//			std::sort(joinJunctions.begin(), joinJunctions.end());
			//			std::list<JunctionMerge>::iterator it_newEnd;
			//			it_newEnd = std::unique(joinJunctions.begin(), joinJunctions.end());
			//			assert(it_newEnd == joinJunctions.end());
			//			joinJunctions.erase(it_newEnd, joinJunctions.end());

			Junction *newJunction = mergeJunctions(joinJunctions);
			assert(newJunction);
			m_junctions.push_back(newJunction);
		}
		joinJunctions.clear();
	}
}

//---------------------------------------------------------
/*
void CenterLineVectorizer::resolveUnvisitedNode(Node *node)
{
	assert(node->m_visited == false);
	assert(node->m_junction == 0);

	JunctionLink link(0,0);

	Node* currNode = node;
	for(int i=0; i<4*SEARCH_WINDOW; i++)
	{
		currNode = currNode->m_next;
		assert(currNode);
		if(currNode->m_junction)
		{
			link.m_first = currNode->m_junction;
			break;
		}
		if(currNode->m_visited) break;
	}

	if(!link.m_first) return;

	int steps = 0;
	for(i=0; i<4*SEARCH_WINDOW; i++)
	{
		currNode = currNode->m_prev;
		assert(currNode);
		if(currNode->m_junction)
		{
			link.m_second = currNode->m_junction;
			break;
		}
		if(currNode->m_visited) break;
		currNode->m_visited = true;
		++steps;
	}

	if(	(link.m_first) && (link.m_second) && (link.m_first != link.m_second)
		&& steps >= 2*SEARCH_WINDOW && norm(link.m_first->m_center - 
link.m_second->m_center) > 2.5 &&
		std::find(m_links.begin(), m_links.end(), link) == m_links.end() &&
		std::find(m_links.begin(), m_links.end(), JunctionLink(link.m_second, 
link.m_first) ) == m_links.end() )
	{
		m_links.push_back(link);
	}
	else
	{
		for(int i=0; i<4*SEARCH_WINDOW; i++)
		{
			currNode = currNode->m_next;
			assert(currNode);
			if(currNode->m_junction) break;
			currNode->m_visited = false;
		}
	}
}
*/
//---------------------------------------------------------

void CenterLineVectorizer::createJunctionPolygon(Junction *junction)
{
	//int size = junction->m_nodes.size();
	//	std::sort(junction->m_nodes.begin(), junction->m_nodes.end(), CompareJunctionNodes(junction->m_center));
	//	std::reverse(junction->m_nodes.begin(), junction->m_nodes.end());	std::deque<Node*>::iterator it, it_temp;
	std::deque<Node *>::iterator it, it_temp;
	for (it = junction->m_nodes.begin(); it != junction->m_nodes.end(); it++) {
		int step = 0;

		Node *currNode;
		Node *otherNode;

		if (it + 1 == junction->m_nodes.end())
			otherNode = *junction->m_nodes.begin();
		else
			otherNode = *(it + 1);

		/*		currNode = (*it)->m_prev;
		if(!currNode->m_visited)
		{
			step=0;
			while(!currNode->m_visited && currNode->m_junction != junction)
			{
				it = junction->m_nodes.insert(it+1, currNode);
				currNode->m_visited = true;
				currNode = currNode->m_prev;
				++step;
			}
			if(currNode->m_junction != junction)
				for(int i=0; i<step; ++i)
				{
					currNode = currNode->m_next;
					currNode->m_visited = false;
					it = junction->m_nodes.erase(it);
					--it;
				}
			else
			{
				if(currNode != otherNode)
				{
					it_temp = std::find(junction->m_nodes.begin(), junction->m_nodes.end(), 
currNode);
					if(it+1 == junction->m_nodes.end())
						std::swap(*junction->m_nodes.begin(), *it_temp);
					else
						std::swap(*(it+1), *it_temp);
				}
				continue;
			}
		}
*/
		currNode = (*it)->m_next;
		if (!currNode->m_visited) {
			step = 0;
			while (!currNode->m_visited && currNode->m_junction != junction) {
				it = junction->m_nodes.insert(it + 1, currNode);
				currNode->m_visited = true;
				currNode = currNode->m_next;
				++step;
			}
			if (currNode->m_junction != junction) //|| currNode != otherNode)
				for (int i = 0; i < step; ++i) {
					currNode = currNode->m_prev;
					currNode->m_visited = false;
					it = junction->m_nodes.erase(it);
					--it;
				}
			else {
				//				assert(*(it+1) == otherNode || (it+1 == junction->m_nodes.end() && *junction->m_nodes.begin() == otherNode) );
				if (currNode != otherNode) {
					Node tempNode = *otherNode;
					*otherNode = *currNode;
					*currNode = tempNode;
					//					it_temp = std::find(junction->m_nodes.begin(), junction->m_nodes.end(), currNode);
					//					assert(it_temp != junction->m_nodes.end());
					//					if(it+1 == junction->m_nodes.end())
					//						std::swap(*junction->m_nodes.begin(), *it_temp);
					//					else
					//						std::swap(*(it+1), *it_temp);
				}
				continue;
			}
		}
	}
}

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

void CenterLineVectorizer::handleLinks()
{
	std::vector<JunctionLink>::iterator it_links;
	for (it_links = m_links.begin(); it_links != m_links.end(); ++it_links) {
		ProtoStroke linkingProtoStroke;

		//initialize points
		TThickPoint p1 = it_links->m_first->m_center;
		TThickPoint p2 = it_links->m_second->m_center - p1;
		linkingProtoStroke.m_points.push_back(TThickPoint(p1,
														  JUNCTION_THICKNESS_RATIO * p1.thick));
		linkingProtoStroke.m_points.push_back(TThickPoint(p1 + 0.5 * p2,
														  JUNCTION_THICKNESS_RATIO * (p1.thick + 0.5 * p2.thick)));
		linkingProtoStroke.m_points.push_back(TThickPoint(p1 + p2,
														  JUNCTION_THICKNESS_RATIO * (p1.thick + p2.thick)));

		//initialize directions
		setEndpointsDirections(linkingProtoStroke);

		//initialize junctions
		/*		linkingProtoStroke.m_startJunction = it_links->first;
		linkingProtoStroke.m_endJunction = it_links->second;*/
		m_protoStrokes.push_back(linkingProtoStroke);
		/*		it_links->first->m_protoStrokes.push_back(&m_protoStrokes.back());
		it_links->second->m_protoStrokes.push_back(&m_protoStrokes.back());*/
	}
}

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

bool isFalseBranch(int k, ProtoStroke *protoStroke, int recursionOrder,
				   TThickPoint &edge)
{
	if (recursionOrder > 3)
		return false;
	++recursionOrder;
	if (k == 0) {
		if (protoStroke->m_points.front().thick > 3.0)
			return false;
		if (protoStroke->m_points.size() <= 20) {
			if (protoStroke->m_endJunction) {
				if (protoStroke->m_endJunction->m_protoStrokes.size() == 1) {
					edge = protoStroke->m_points.back();
					protoStroke->m_endJunction->m_protoStrokes.clear();
					protoStroke->m_points.clear();
					protoStroke->m_startJunction = 0;
					protoStroke->m_endJunction = 0;
					return true;
				} else if (protoStroke->m_endJunction->m_protoStrokes.size() == 2) {
					ProtoStroke *nextProtoStroke;
					int next_k;

					if (protoStroke->m_endJunction->m_protoStrokes[0] == protoStroke)
						nextProtoStroke = protoStroke->m_endJunction->m_protoStrokes[1];
					else if (protoStroke->m_endJunction->m_protoStrokes[1] == protoStroke)
						nextProtoStroke = protoStroke->m_endJunction->m_protoStrokes[0];
					if (nextProtoStroke->m_startJunction == protoStroke->m_endJunction)
						next_k = 0;
					else if (nextProtoStroke->m_endJunction == protoStroke->m_endJunction)
						next_k = 1;

					if (isFalseBranch(next_k, nextProtoStroke, recursionOrder, edge)) {
						protoStroke->m_endJunction->m_protoStrokes.clear();
						protoStroke->m_points.clear();
						protoStroke->m_startJunction = 0;
						protoStroke->m_endJunction = 0;
						return true;
					} else
						return false;
				} else
					return false;
			}
			edge = protoStroke->m_points.back();
			protoStroke->m_points.clear();
			protoStroke->m_startJunction = 0;
			protoStroke->m_endJunction = 0;
			return true;
		}
	} else if (k == 1) {
		if (protoStroke->m_points.back().thick > 3.0)
			return false;
		if (protoStroke->m_points.size() <= 20) {
			if (protoStroke->m_startJunction) {
				if (protoStroke->m_startJunction->m_protoStrokes.size() == 1) {
					edge = protoStroke->m_points.front();
					protoStroke->m_startJunction->m_protoStrokes.clear();
					protoStroke->m_points.clear();
					protoStroke->m_startJunction = 0;
					protoStroke->m_endJunction = 0;
					return true;
				} else if (protoStroke->m_startJunction->m_protoStrokes.size() == 2) {
					ProtoStroke *nextProtoStroke;
					int next_k;

					if (protoStroke->m_startJunction->m_protoStrokes[0] == protoStroke)
						nextProtoStroke = protoStroke->m_startJunction->m_protoStrokes[1];
					else if (protoStroke->m_startJunction->m_protoStrokes[1] == protoStroke)
						nextProtoStroke = protoStroke->m_startJunction->m_protoStrokes[0];
					if (nextProtoStroke->m_startJunction == protoStroke->m_startJunction)
						next_k = 0;
					else if (nextProtoStroke->m_endJunction == protoStroke->m_startJunction)
						next_k = 1;

					if (isFalseBranch(next_k, nextProtoStroke, recursionOrder, edge)) {
						protoStroke->m_startJunction->m_protoStrokes.clear();
						protoStroke->m_points.clear();
						protoStroke->m_startJunction = 0;
						protoStroke->m_endJunction = 0;
						return true;
					} else
						return false;
				} else
					return false;
			}
			edge = protoStroke->m_points.front();
			protoStroke->m_points.clear();
			protoStroke->m_startJunction = 0;
			protoStroke->m_endJunction = 0;
			return true;
		}
	}
	return false;
}

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

void CenterLineVectorizer::handleJunctions()
{
	std::vector<Junction *>::iterator it_junction = m_junctions.begin();
	for (it_junction = m_junctions.begin(); it_junction < m_junctions.end();
		 it_junction++) {
		Junction *currJunction = *it_junction;

		if (currJunction->m_protoStrokes.empty()) {
			std::deque<Node *>::iterator it_node;
			for (it_node = currJunction->m_nodes.begin();
				 it_node != currJunction->m_nodes.end(); ++it_node) {
				(*it_node)->m_junction = 0;
				(*it_node)->m_visited = false;
			}
			delete currJunction;
			it_junction = m_junctions.erase(it_junction);
			--it_junction;
			continue;
		}

		std::vector<ProtoStroke *>::iterator it;

		//set polygon center
		/*		double size = currJunction->m_nodes.size();
		std::deque<Node*>::iterator it_node;
		for(it_node = currJunction->m_nodes.begin(); 
it_node!=currJunction->m_nodes.end(); it_node++)
		{
			currJunction->m_center += convert( (*it_node)->m_pixel->m_pos );
		}
		currJunction->m_center = currJunction->m_center * (1/size);

		size = currJunction->m_protoStrokes.size();
		if(size == 0)
		{
			currJunction->m_center.thick = THICKNESS_LIMIT;
			continue;
		}
*/
		double size = currJunction->m_protoStrokes.size();
		for (it = currJunction->m_protoStrokes.begin();
			 it != currJunction->m_protoStrokes.end(); it++) {
			if (((*it))->m_startJunction == currJunction && (*it)->m_endJunction == currJunction) {
				currJunction->m_center += 0.5 * (*it)->m_points.front();
				currJunction->m_center += 0.5 * (*it)->m_points.back();
			} else if (((*it))->m_startJunction == currJunction) {
				currJunction->m_center += (*it)->m_points.front();
			} else if ((*it)->m_endJunction == currJunction) {
				currJunction->m_center += (*it)->m_points.back();
			}
		}
		currJunction->m_center = currJunction->m_center * (1 / size);
	}

#ifdef DEBUG
	int j = 0;
	for (; j < (int)m_nodes.size(); j++) {
		Node *node = m_nodes[j];
		if (!node->m_visited && node->m_pixel) {
			m_unvisitedNodes.push_back(convert(node->m_pixel->m_pos));
			//			resolveUnvisitedNode(node);
		}
	}

	for (j = 0; j < (int)m_nodes.size(); j++) {
		Node *node = m_nodes[j];
		//		assert(node->m_visited || !node->m_pixel);
		if (node->m_pixel == 0 || node->m_flag)
			continue;
		outputNodes(node);
	}
#endif

	joinJunctions();
	//	return;
	for (it_junction = m_junctions.begin(); it_junction < m_junctions.end();
		 it_junction++) {
		Junction *currJunction = *it_junction;
		if (currJunction->m_center.thick < 2 * THICKNESS_LIMIT)
			currJunction->m_center.thick = 2 * THICKNESS_LIMIT;

		std::vector<ProtoStroke *>::iterator it, it_next;

		createJunctionPolygon(currJunction);

		currJunction->m_junctionOrder = currJunction->m_protoStrokes.size();

		if (currJunction->m_protoStrokes.size() == 3) {
			std::vector<ProtoStroke *>::iterator it_e =
				currJunction->m_protoStrokes.end();
			std::vector<ProtoStroke *>::iterator it_candidate1 = it_e, it_candidate2 =
																		   it_e;
			int candidate1_k, candidate2_k;
			double candidateProbability = JOIN_LIMIT;
			int k, next_k;
			for (it = currJunction->m_protoStrokes.begin(); it != it_e; it++) {
				if (((*it))->m_startJunction == currJunction)
					k = 0;
				else if ((*it)->m_endJunction == currJunction)
					k = 1;
				else
					assert(false);

				if (isFalseBranch(k, *it, 0, currJunction->m_center)) {
					currJunction->m_protoStrokes.erase(it);
					currJunction->m_locked = true;
					break;
				}

				it_next = it;
				it_next++;
				for (; it_next != it_e; it_next++) {
					double probability;
					if (((*it_next))->m_startJunction == currJunction)
						next_k = 0;
					else if ((*it_next)->m_endJunction == currJunction)
						next_k = 1;
					else
						assert(false);

					if (k == 0) {
						if (next_k == 0)
							probability = 1 - norm2((*it)->m_startDirection +
													(*it_next)->m_startDirection);
						else
							probability = 1 - norm2((*it)->m_startDirection +
													(*it_next)->m_endDirection);
					} else {
						if (next_k == 0)
							probability = 1 - norm2((*it)->m_endDirection +
													(*it_next)->m_startDirection);
						else
							probability = 1 - norm2((*it)->m_endDirection +
													(*it_next)->m_endDirection);
					}
					if (probability > candidateProbability ||
						(*it)->m_points.size() < 3) {
						candidateProbability = probability;
						it_candidate1 = it;
						it_candidate2 = it_next;
						candidate1_k = k;
						candidate2_k = next_k;
					}
				}
			}
			if (currJunction->m_protoStrokes.size() == 3 && it_candidate1 != it_e &&
				it_candidate2 != it_e) {
				TThickPoint p1 = (candidate1_k == 0) ? (*it_candidate1)->m_points.front()
													 : (*it_candidate1)->m_points.back();
				TThickPoint p2 = (candidate2_k == 0) ? (*it_candidate2)->m_points.front()
													 : (*it_candidate2)->m_points.back();
				currJunction->m_center = 0.5 * (p1 + p2);
			}
		}
		if (false || currJunction->m_protoStrokes.size() == 2) {
			TPointD a1, a2, b1, b2;

			if (currJunction->m_protoStrokes[0]->m_startJunction == currJunction) {
				a1 = currJunction->m_protoStrokes[0]->m_points.front();
				a2 = currJunction->m_protoStrokes[0]->m_startDirection - a1;
			} else {
				assert(currJunction->m_protoStrokes[0]->m_endJunction == currJunction);
				a1 = currJunction->m_protoStrokes[0]->m_points.back();
				a2 = currJunction->m_protoStrokes[0]->m_endDirection - a1;
			}

			if (currJunction->m_protoStrokes[1]->m_startJunction == currJunction) {
				b1 = currJunction->m_protoStrokes[1]->m_points.front();
				b2 = currJunction->m_protoStrokes[1]->m_startDirection - b1;
			} else {
				assert(currJunction->m_protoStrokes[1]->m_endJunction == currJunction);
				b1 = currJunction->m_protoStrokes[1]->m_points.back();
				b2 = currJunction->m_protoStrokes[1]->m_endDirection - b1;
			}

/*			TPointD intersection;
			RectIntersectionResult res;
			res = rectsIntersect(a1, a2, b1, b2, intersection);
//			assert(res == intersected);
			if( norm(intersection - currJunction->m_center) < 5 && res == 
intersected)
			{
				currJunction->m_center.x = intersection.x;
				currJunction->m_center.y = intersection.y;
			}
*/		}

if (true || currJunction->m_nodes.size() <= 12 ||
	currJunction->isConvex()) {
	TThickPoint junctionPoint(currJunction->m_center);
	for (it = currJunction->m_protoStrokes.begin();
		 it != currJunction->m_protoStrokes.end(); it++) {
		assert((*it)->m_startJunction == currJunction || (*it)->m_endJunction == currJunction);
		if ((*it)->m_startJunction == currJunction) {
			TPointD deltaPoint = (*it)->m_points.front() - junctionPoint;
			double deltaThick = (*it)->m_points.front().thick -
								junctionPoint.thick;
			if (norm2(deltaPoint) > 0.2) {
				if (norm2(deltaPoint) > 2) {
					(*it)->m_points.push_front(
						TThickPoint(junctionPoint + 0.7 * deltaPoint, junctionPoint.thick +
																		  0.7 * deltaThick));
					(*it)->m_points.push_front(
						TThickPoint(junctionPoint + 0.4 * deltaPoint, junctionPoint.thick +
																		  0.4 * deltaThick));
				}
				(*it)->m_points.push_front(
					TThickPoint(junctionPoint + 0.0 * deltaPoint, junctionPoint.thick +
																	  0.0 * deltaThick));
			}
		}
		if ((*it)->m_endJunction == currJunction) {
			TPointD deltaPoint = (*it)->m_points.back() - junctionPoint;
			double deltaThick = (*it)->m_points.back().thick - junctionPoint.thick;
			if (norm2(deltaPoint) > 0.2) {
				if (norm2(deltaPoint) > 2) {
					(*it)->m_points.push_back(
						TThickPoint(junctionPoint + 0.7 * deltaPoint, junctionPoint.thick +
																		  0.7 * deltaThick));
					(*it)->m_points.push_back(
						TThickPoint(junctionPoint + 0.4 * deltaPoint, junctionPoint.thick +
																		  0.4 * deltaThick));
				}
				(*it)->m_points.push_back(
					TThickPoint(junctionPoint + 0.0 * deltaPoint, junctionPoint.thick +
																	  0.0 * deltaThick));
			}
		}
	}
} else {
	ProtoRegion region(currJunction->isConvex());
	region.m_center = currJunction->m_center;

	std::deque<Node *>::iterator it_nodes;
	for (it_nodes = currJunction->m_nodes.begin();
		 it_nodes != currJunction->m_nodes.end(); it_nodes++)
		region.m_points.push_back(
			TThickPoint(convert((*it_nodes)->m_pixel->m_pos), 0.0));

	region.m_points.push_back(
		TThickPoint(convert((*currJunction->m_nodes.begin())->m_pixel->m_pos), 0.0));

	m_protoRegions.push_back(region);

	for (it = currJunction->m_protoStrokes.begin();
		 it != currJunction->m_protoStrokes.end(); it++) {
		if ((*it)->m_startJunction == currJunction && (*it)->m_endJunction == currJunction) {
			if ((*it)->m_points.size() < 10)
				(*it)->m_points.clear();
		}

		if ((*it)->m_startJunction == currJunction) {
			(*it)->m_startJunction = 0;
		}
		if ((*it)->m_endJunction == currJunction) {
			(*it)->m_endJunction = 0;
		}
	}
	currJunction->m_protoStrokes.clear();
}
	}

	//handleLinks();
}

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

void CenterLineVectorizer::createStrokes()
{
	std::list<ProtoStroke>::iterator it_centerlines = m_protoStrokes.begin();
	for (; it_centerlines != m_protoStrokes.end(); it_centerlines++) {
		if (!it_centerlines->m_points.empty()) {
			joinProtoStrokes(&(*it_centerlines));
			if (it_centerlines->m_points.size() <= 1)
				continue;

			std::vector<TThickPoint> points;
			std::deque<TThickPoint>::iterator it;

			if (m_configuration.m_smoothness > 0) // gmt: prima era true
			{
				// gmt, 30-3-2006:
				/*
				it = it_centerlines->m_points.begin()+1;
				for(;;)
				{
					if( it >= it_centerlines->m_points.end()-(m_configuration.m_smoothness+1)) break;
					for(int j=0; j<m_configuration.m_smoothness;j++)
						it = it_centerlines->m_points.erase(it);
					++it;
				}
        */
				// gmt. credo che l'idea sia di lasciare un punto e
				// toglierne m_configuration.m_smoothness e cosi' via fino alla fine
				// n.b. ho modificato il comportamento sugli eventuali ultimi punti
				// (<m_configuration.m_smoothness) che prima venivano lasciati e adesso
				// vengono tolti
				std::deque<TThickPoint> &points = it_centerlines->m_points;
				it = points.begin();
				for (;;) {
					if (it == points.end())
						break;
					++it;
					int maxD = std::distance(it, points.end());
					int d = tmin(maxD, m_configuration.m_smoothness);
					if (d == 0)
						break;
					it = points.erase(it, it + d);
				}
			}
			it = it_centerlines->m_points.begin();
			points.push_back(*it);
			assert(it->thick >= THICKNESS_LIMIT);
			TThickPoint old = *it;
			for (it++; it != it_centerlines->m_points.end(); ++it) {
				assert(it->thick >= THICKNESS_LIMIT);
				TThickPoint point(0.5 * (*it + old), 0.5 * (it->thick + old.thick));
				points.push_back(point);
				old = *it;
			}
			points.push_back(it_centerlines->m_points.back());

			int n = points.size();
			if (n > 5) {
				double *rr = new double[n];
				int i;
				for (i = 0; i < n; i++)
					rr[i] = points[i].thick;
				for (i = 2; i < n - 2; i++)
					points[i].thick = 0.8 * (1 / 5.0) *
									  (/*rr[i-4]+rr[i-3]*/ +rr[i - 2] + rr[i - 1] + rr[i] + rr[i + 1] + rr[i + 2] /*+rr[i+3]+rr[i+4]*/);
				points[0].thick = points[1].thick = points[2].thick;			 //=points[3].thick;//=points[4].thick;
				points[n - 1].thick = points[n - 2].thick = points[n - 3].thick; //=points[n-4].thick;//=points[n-5].thick;
				delete[] rr;
			}
			TStroke *stroke = TStroke::interpolate(points, m_configuration.m_interpolationError);
			stroke->setStyle(m_configuration.m_strokeStyleId); // FirstUserStyle+2);
			m_vimage->addStroke(stroke);
		}
	}
}

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

void CenterLineVectorizer::createRegions()
{
	std::list<ProtoRegion>::iterator it_regions;
	for (it_regions = m_protoRegions.begin(); it_regions != m_protoRegions.end(); it_regions++) {
		assert(it_regions->m_points.size() > 4);
		TStroke *stroke = TStroke::interpolate(it_regions->m_points, 0.1);
		stroke->setSelfLoop();
		stroke->setStyle(m_configuration.m_strokeStyleId);
		m_vimage->addStroke(stroke);

		/*
		if(it_regions->m_isConvex)
			m_vimage->fill(it_regions->m_center, TPalette::DefaultStrokeStyle);
		else
		{
			TPointD internalPoint;
			internalPoint = rotate90(it_regions->m_points[1] - 
it_regions->m_points[0]);
			assert(norm2(internalPoint) != 0);
			normalize(internalPoint);
			internalPoint = internalPoint*0.001 + 0.5*(it_regions->m_points[1] + 
it_regions->m_points[0]);
			m_vimage->fill(internalPoint, TPalette::DefaultStrokeStyle);
		}*/
	}
}

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

void CenterLineVectorizer::vectorize()
{
	int j = 0;

	for (j = 0; j < (int)m_nodes.size(); j++) {
		Node *node = m_nodes[j];
		if (node->m_pixel == 0 || node->m_visited)
			continue;
		traceLine(node->m_pixel);
	}

	handleJunctions();
	createStrokes();
	createRegions();
}

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

#ifdef DEBUG
void printTime(TStopWatch &sw, string name)
{
	ostrstream ss;
	ss << name << " : ";
	sw.print(ss);
	ss << '\n' << '\0';
	string s = ss.str();
	ss.freeze(false);
	TSystem::outputDebug(s);
}
#endif

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

#ifdef DEBUG
void CenterLineVectorizer::outputInks()
{
	for (int i = 0; i < m_dataRaster->getLx() * m_dataRaster->getLy(); i++) {
		DataPixel *pix = m_dataRaster->pixels() + i;
		if (pix->m_ink)
			getOutPix(pix) = TPixel32::Black;
		else
			getOutPix(pix) = TPixel32::White;
	}
}
#endif

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

TVectorImageP centerlineVectorize(const TImageP &image, const VectorizerConfiguration &c, TPalette *palette)
{
	CenterLineVectorizer vectorizer(c, palette);

	TRasterImageP ri = image;
	TToonzImageP vi = image;
	if (ri)
		vectorizer.makeDataRaster(ri->getRaster());
	else
		vectorizer.makeDataRaster(vi->getRaster());
	for (int i = 0; i < (int)vectorizer.m_dataRasterArray.size(); i++) {
		vectorizer.m_dataRaster = vectorizer.m_dataRasterArray[i].second;
		vectorizer.m_configuration.m_strokeStyleId = vectorizer.m_dataRasterArray[i].first;
		vectorizer.init();
		vectorizer.vectorize();
	}

	renormalizeImage(vectorizer.m_vimage.getPointer());

	return vectorizer.m_vimage;
}

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

class OutlineVectorizer : public CenterLineVectorizer
{

public:
	list<vector<TThickPoint>> m_protoOutlines;

	void traceOutline(Node *initialNode);

	OutlineVectorizer(const VectorizerConfiguration &configuration, TPalette *palette) : CenterLineVectorizer(configuration, palette) {}
	~OutlineVectorizer();

	void createOutlineStrokes();

private:
	// not implemented
	OutlineVectorizer(const OutlineVectorizer &);
	OutlineVectorizer &operator=(const OutlineVectorizer &);
};

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

OutlineVectorizer::~OutlineVectorizer()
{
	m_protoOutlines.clear();
}

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

void OutlineVectorizer::traceOutline(Node *initialNode)
{
	Node *startNode = initialNode;
	Node *node;
	do {
		if (!startNode)
			break;
		node = findOtherSide(startNode);
		if (!node)
			break;

		double startDist2 = computeDistance2(startNode, node);
		if (startDist2 > 0.1)
			break;

		startNode = startNode->m_next;
		//		assert(startNode != initialNode);
	} while (startNode != initialNode);

	if (!startNode)
		return;
	node = startNode;
	std::vector<TThickPoint> points;
	do {
		node = node->m_next;
		if (!node)
			break;
		node->m_visited = true;
		points.push_back(TThickPoint(convert(node->m_pixel->m_pos), 0));
	} while (node != startNode);
	m_protoOutlines.push_back(points);
}

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

void OutlineVectorizer::createOutlineStrokes()
{
	m_vimage->enableRegionComputing(true, false);
	int j;

	for (j = 0; j < (int)m_nodes.size(); j++) {
		Node *node = m_nodes[j];
		if (node->m_pixel == 0 || node->m_visited)
			continue;
		traceOutline(node);
	}

#ifdef DEBUG
	for (j = 0; j < (int)m_nodes.size(); j++) {
		Node *node = m_nodes[j];
		//		assert(node->m_visited || !node->m_pixel);
		if (node->m_pixel == 0 || node->m_flag)
			continue;
		outputNodes(node);
	}
#endif

	std::list<std::vector<TThickPoint>>::iterator it_outlines = m_protoOutlines.begin();
	for (; it_outlines != m_protoOutlines.end(); it_outlines++) {
		if (it_outlines->size() > 3) {
			std::vector<TThickPoint> points;
			std::vector<TThickPoint>::iterator it;

			if (it_outlines->size() > 10) {
				it = it_outlines->begin() + 1;
				for (;;) {
					if (it >= it_outlines->end() - (m_configuration.m_smoothness + 1))
						break;
					for (j = 0; j < m_configuration.m_smoothness; j++)
						it = it_outlines->erase(it);
					++it;
				}
			}
			points.push_back(it_outlines->front());
			it = it_outlines->begin();
			TThickPoint old = *it;
			++it;
			for (; it != it_outlines->end(); ++it) {
				TThickPoint point((1 / 2.0) * (*it + old));
				points.push_back(point);
				old = *it;
			}
			points.push_back(it_outlines->back());
			points.push_back(it_outlines->front());

			/*TPointD internalPoint;
			internalPoint = rotate90(points[1] - points[0]);
			assert(norm2(internalPoint) != 0);
			normalize(internalPoint);
			internalPoint = internalPoint*0.001 + points[0];

			TPointD externalPoint;
			externalPoint = rotate270(points[1] - points[0]);
			assert(norm2(externalPoint) != 0);
			normalize(externalPoint);
			externalPoint = externalPoint*0.001 + points[0];*/
			TStroke *stroke = TStroke::interpolate(points, m_configuration.m_interpolationError);
			stroke->setStyle(m_configuration.m_strokeStyleId);
			stroke->setSelfLoop();
			m_vimage->addStroke(stroke);
			//			TStroke *stroke = TStroke::interpolate(points, m_configuration.interpolationError);
			//     		stroke->setStyle(strokeStyleId);
			//			m_vimage->addStroke(stroke);//, TPalette::DefaultStrokeStyle, TPalette::DefaultFillStyle);
			//internalPoints.push_back(internalPoint);
			//m_vimage->fill(externalPoint, TPalette::DefaultFillStyle);
		}
	}

	//m_vimage->autoFill(m_configuration.m_strokeStyleId);
	//for (UINT i=0; i<internalPoints.size(); i++)
	//  m_vimage->fill(internalPoints[i],fillStyleId);
}

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

TVectorImageP outlineVectorize(const TImageP &image, const VectorizerConfiguration &configuration, TPalette *palette)
{
	TVectorImageP out;

	OutlineVectorizer vectorizer(configuration, palette);

	TRasterImageP ri = image;
	TToonzImageP vi = image;
	if (ri)
		vectorizer.makeDataRaster(ri->getRaster());
	else
		vectorizer.makeDataRaster(vi->getRaster());
	int layersCount = vectorizer.m_dataRasterArray.size();
	if (layersCount > 1) {
		out = new TVectorImage();
		out->setPalette(palette);
	}
	for (int i = 0; i < (int)layersCount; i++) {
		vectorizer.m_dataRaster = vectorizer.m_dataRasterArray[i].second;
		vectorizer.m_configuration.m_strokeStyleId = vectorizer.m_dataRasterArray[i].first;
		vectorizer.m_protoOutlines.clear();
		vectorizer.init();
		vectorizer.createOutlineStrokes();
		renormalizeImage(vectorizer.m_vimage.getPointer());
		vectorizer.m_vimage->setPalette(palette);
		if (layersCount > 1)
			out->mergeImage(vectorizer.m_vimage, TAffine());

		if (i != (int)layersCount - 1)
			vectorizer.m_vimage = new TVectorImage();
	}

	return (layersCount == 1) ? vectorizer.m_vimage : out;
}

//---------------------------------------------------------
void applyFillColors(TRegion *r, const TRasterP &ras, TPalette *palette, const TPoint &offset, bool leaveUnpainted)
{
	TPointD pd;
	TRasterCM32P rt = (TRasterCM32P)ras;
	TRaster32P rr = (TRaster32P)ras;
	TRasterGR8P rgr = (TRasterGR8P)ras;
	assert(rt || rr || rgr);

	if (/*r->getStyle()==0 && */ r->getInternalPoint(pd)) {
		TPoint p = convert(pd) + offset;
		if (ras->getBounds().contains(p)) {
			int styleId = 0;
			if (rt) {
				TPixelCM32 col = rt->pixels(p.y)[p.x];
				int tone = col.getTone();
				if (tone < 100) {
					styleId = col.getInk();
					if (styleId == 0 && !leaveUnpainted)
						styleId = col.getPaint();
				} else if (!leaveUnpainted || col.getPaint() == 0) {
					styleId = col.getPaint();
					if (styleId == 0)
						styleId = col.getInk();
				}
			} else {
				TPixel32 color;
				if (rr)
					color = rr->pixels(p.y)[p.x];
				else {
					int val = rgr->pixels(p.y)[p.x].value;
					color = (val < 80) ? TPixel32::Black : TPixel32::White;
				}
				styleId = palette->getClosestStyle(color);
				if (palette->getStyle(styleId)->getMainColor() != color)
					styleId = palette->addStyle(color);
			}
			r->fill(pd, styleId);
		}
	}

	for (int i = 0; i < (int)r->getSubregionCount(); i++)
		applyFillColors(r->getSubregion(i), ras, palette, offset, leaveUnpainted);
}

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

void applyFillColors(TVectorImageP vi, const TImageP &img, TPalette *palette, bool leaveUnpainted)
{
	TToonzImageP ti = (TToonzImageP)img;
	TRasterImageP ri = (TRasterImageP)img;
	TRectD vb = vi->getBBox();
	TRectD tb = img->getBBox();

	TRasterP ras;

	if (ti)
		ras = ti->getRaster();
	else if (ri)
		ras = ri->getRaster();
	else
		assert(false);

	vi->findRegions();
	for (int i = 0; i < (int)vi->getRegionCount(); i++)
		applyFillColors(vi->getRegion(i), ras, palette, convert(tb.getP00() - vb.getP00()), leaveUnpainted);
}

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

TVectorImageP vectorize(const TImageP &img, const VectorizerConfiguration &c, TPalette *plt)
{
	TVectorImageP vi = c.m_outline
						   ? outlineVectorize(img, c, plt)
						   : centerlineVectorize(img, c, plt);

	applyFillColors(vi, img, plt, c.m_leaveUnpainted);

	//else if (c.m_outline)
	//  vi->autoFill(c.m_strokeStyleId);

	TPointD center;
	if (TToonzImageP ti = img)
		center = ti->getRaster()->getCenterD();
	else if (TRasterImageP ri = img)
		center = ri->getRaster()->getCenterD();

	TAffine aff = TTranslation(-center);
	vi->transform(aff);

	return vi;
}

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