#include "tcenterlinevectP.h"

//#define _SSDEBUG                                              // Uncomment to enable the debug viewer

//#define _UPDATE                                               // Shows borders updated to current time

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

// Forward declarations

struct VectorizationContext;

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

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

//      Rationale

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

/*

  NOTE: Input is a vector of Contours, representing the borders of a

        polygonal region. First border is the outer one, followed by internal

        counter-borders; each Contour is itself a vector of "ContourNode"s,

        ordered so that the region to be thinned is at the RIGHT of segments

        formed by successive nodes.

        Output is a Graph structure representing the straight skeleton of the

        thinned region. Original contours survive the thinning procedure.

*/

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

//      Straight Skeleton Debugger

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

#ifdef _SSDEBUG

#include <qwidget></qwidget>

#include <qtransform></qtransform>

#include <qeventloop></qeventloop>

#include <qpainter></qpainter>

#include <qmouseevent></qmouseevent>

#include <qwheelevent></qwheelevent>

class SSDebugger : public QWidget

{

public:

	VectorizationContext &m_context;

	QPoint m_pos,

		m_pressPos;

	double m_scale;

	QTransform m_transform;

	QEventLoop m_loop;

public:

	double m_height;

public:

	SSDebugger(VectorizationContext &context);

	~SSDebugger() {}

	void loop() { m_loop.exec(); }

	void paintEvent(QPaintEvent *event);

	void keyPressEvent(QKeyEvent *event);

	void mouseMoveEvent(QMouseEvent *event);

	void mousePressEvent(QMouseEvent *event);

	void mouseReleaseEvent(QMouseEvent *event);

	void wheelEvent(QWheelEvent *event);

	inline QPoint winToWorld(int x, int y);

	inline QPoint worldToWin(double x, double y);

	inline QPointF winToWorldF(int x, int y);

	inline bool isOnScreen(ContourNode *node);

	// Node Updates

	TPointD updated(ContourNode *input);

};

#endif // _SSDEBUG

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

//      Classes

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

class ContourEdge

{

public:

	enum { NOT_OPPOSITE = 0x1 };

public:

	TPointD m_direction;

	unsigned short m_attributes;

public:

	ContourEdge() : m_attributes(0) {}

	ContourEdge(TPointD dir) : m_direction(dir), m_attributes(0) {}

	int hasAttribute(int attr) const { return m_attributes & attr; }

	void setAttribute(int attr) { m_attributes |= attr; }

	void clearAttribute(int attr) { m_attributes &= ~attr; }

};

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

class IndexTable

{

public:

	typedef std::list<contournode *=""> IndexColumn;</contournode>

	std::vector<indexcolumn> m_columns; //!< Countours set by 'column identifier'.</indexcolumn>

	std::vector<int> m_identifiers;		//!< Column identifiers by original contour index.</int>

	// NOTE: Contours are stored in 'comb' structure (vector of lists) since contours may both

	// be SPLIT (new entry in a list) and MERGED (two lists merge).

public:

	IndexTable() {}

	IndexColumn *operator[](int i) { return &m_columns[i]; }

	IndexColumn &columnOfId(int id) { return m_columns[m_identifiers[id]]; }

	// Initialization

	void build(ContourFamily &family);

	void clear();

	// Specific handlers

	IndexColumn::iterator find(ContourNode *index);

	void merge(IndexColumn::iterator index1, IndexColumn::iterator index2);

	void remove(IndexColumn::iterator index);

};

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

class Event

{

public:

	/*! \remark  Values are sorted by preference at simultaneous events.    */

	enum Type //! An event's possible types.

	{

		special,		  //!< A vertex event that is also an edge event (V case).

		edge,			  //!< An edge shrinks to 0 length.

		vertex,			  //!< Two contour nodes clash.

		split_regenerate, //!< Placeholder type for split events that must be regenerated.

		split,			  //!< An edge is split by a clashing contour node.

		failure

	};

public:

	double m_height;

	double m_displacement;

	ContourNode *m_generator;

	ContourNode *m_coGenerator;

	Type m_type;

	unsigned int m_algoritmicTime;

	VectorizationContext *m_context;

public:

	// In-builder event constructor

	Event(ContourNode *generator, VectorizationContext *context);

	// Event calculators

	inline void calculateEdgeEvent();

	inline void calculateSplitEvent();

	// Auxiliary event calculators

	inline double splitDisplacementWith(ContourNode *plane);

	inline bool tryRayEdgeCollisionWith(ContourNode *edge);

	// Event handlers

	inline bool process();

	inline void processEdgeEvent();

	inline void processMaxEvent();

	inline void processSplitEvent();

	inline void processVertexEvent();

	inline void processSpecialEvent();

private:

	inline bool testRayEdgeCollision(ContourNode *opposite,

									 double &displacement, double &height, double &side1, double &side2);

};

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

struct EventGreater {

	bool operator()(const Event &event1, const Event &event2) const

	{

		return event1.m_height > event2.m_height || (event1.m_height == event2.m_height && event1.m_type > event2.m_type);

	}

};

class Timeline : public std::priority_queue<event, std::vector<event="">, EventGreater></event,>

{

public:

	Timeline() {}

	// NOTE: Timeline construction contains the most complex part of vectorization;

	// progress bar partial notification happens there, so thisVectorizer's signal

	// emission methods must be passed and used.

	void build(ContourFamily &polygons, VectorizationContext &context, VectorizerCore *thisVectorizer);

};

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

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

//    Preliminary methods/functions

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

// IndexTable methods

void IndexTable::build(ContourFamily &family)

{

	unsigned int i;

	m_columns.resize(family.size());

	m_identifiers.resize(family.size());

	//NOTE: At the beginning, m_identifiers= 1, .. , m_columns.size() - 1;

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

		m_identifiers[i] = i;

		m_columns[i].push_back(&family[i][0]);

		//Each node referenced in the Table is signed as 'head' of the cirular list.

		family[i][0].setAttribute(ContourNode::HEAD);

	}

}

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

//Explanation: during the skeletonization process, ContourNodes and calculated

//Events are unaware of global index-changes generated by other events, so

//the position of index stored in one Event has to be retrieved in the

//IndexTable before event processing begins.

//NOTE: Can this be done in a more efficient way?...

inline IndexTable::IndexColumn::iterator

IndexTable::find(ContourNode *sought)

{

	int indexId = m_identifiers[sought->m_ancestorContour];

	IndexColumn::iterator res;

	// Search for the HEAD attribute in index's Contour

	for (; !sought->hasAttribute(ContourNode::HEAD); sought = sought->m_next)

		;

	// Finally, find index through our column

	for (res = m_columns[indexId].begin(); (*res) != sought; ++res)

		;

	return res;

}

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

// Handles active contour merging due to split/vertex events

void IndexTable::merge(IndexColumn::iterator index1, IndexColumn::iterator index2)

{

	IndexColumn::iterator current;

	int identifier1 = m_identifiers[(*index1)->m_ancestorContour],

		identifier2 = m_identifiers[(*index2)->m_ancestorContour];

	remove(index2); // We maintain only one index of the merged contour

	// Now, append columns

	if (!m_columns[identifier2].empty()) {

		append<indextable::indexcolumn, indextable::indexcolumn::reverse_iterator="">(</indextable::indexcolumn,>

			m_columns[identifier1], m_columns[identifier2]);

		m_columns[identifier2].clear();

	}

	// Then, update stored identifiers

	for (unsigned int k = 0; k < m_columns.size(); ++k) {

		if (m_identifiers[k] == identifier2)

			m_identifiers[k] = identifier1;

	}

}

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

// Removes given index in Table

inline void IndexTable::remove(IndexColumn::iterator index)

{

	m_columns[m_identifiers[(*index)->m_ancestorContour]].erase(index);

}

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

inline void IndexTable::clear()

{

	m_columns.clear(), m_identifiers.clear();

}

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

//                    Straight Skeleton Algorithm

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

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

//      Global Variables

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

struct VectorizationContext {

	VectorizerCoreGlobals *m_globals;

	//Globals

	unsigned int m_totalNodes;	 //Number of original contour nodes

	unsigned int m_contoursCount;  //Number of contours in input region

	IndexTable m_activeTable;	  //Index table of active contours

	SkeletonGraph *m_output;	   //Output skeleton of input region

	double m_currentHeight;		   //Height of our 'roof-flooding' process

	Timeline m_timeline;		   //Ordered queue of all possible events

	unsigned int m_algoritmicTime; //Number of events precessed up to now

	//Containers

	std::vector<contouredge> m_edgesHeap;</contouredge>

	std::vector<contournode> m_nodesHeap; //of *non-original* nodes only</contournode>

	unsigned int m_nodesHeapCount;		  //number of nodes used in nodesHeap

	//'Linear Axis-added' *pseudo-original* nodes and edges

	std::vector<contournode> m_linearNodesHeap;</contournode>

	std::vector<contouredge> m_linearEdgesHeap;</contouredge>

	unsigned int m_linearNodesHeapCount;

public:

	VectorizationContext(VectorizerCoreGlobals *globals)

		: m_globals(globals) {}

	ContourNode *getNode() { return &m_nodesHeap[m_nodesHeapCount++]; }

	ContourNode *getLinearNode() { return &m_linearNodesHeap[m_linearNodesHeapCount]; }

	ContourEdge *getLinearEdge() { return &m_linearEdgesHeap[m_linearNodesHeapCount++]; }

	inline void addLinearNodeBefore(ContourNode *node);

	inline void repairDegenerations(const std::vector<contournode *=""> °enerates);</contournode>

	inline void prepareGlobals();

	inline void prepareContours(ContourFamily &family);

	inline void newSkeletonLink(unsigned int cur, ContourNode *node);

};

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

//WARNING: To be launched only *after* prepareContours - node countings happen there

inline void VectorizationContext::prepareGlobals()

{

	//NOTE: Let n be the total number of nodes in the family, k the number of split events

	//      effectively happening in the process, m the number of original contours of the family.

	//      Now:

	//        * Each split event eliminates its generating reflex node and introduces

	//          two convex nodes

	//        * Each edge event eliminates its couple of generating nodes and

	//          introduces one new convex node

	//        * Each max event eliminates 3 generating nodes without introducing new ones

	//So, split events introduce 2k non-original nodes, and (k-m+2) is the number of max events

	//necessarily happening, since (m-1) are the *merging* split events.

	//On the n+k-3(k-m+2) nodes remaining for pure edge events, as many non-original nodes are inserted.

	//=> This yields 2k + n-2k+3m-6= n+3m-6 non-original nodes. Contemporaneous events such as

	//vertex and special events can only decrease the number of non-original nodes requested.

	//Initialize non-original nodes container

	m_nodesHeap.resize(m_totalNodes + 3 * m_contoursCount - 6);

	m_nodesHeapCount = 0;

	//Reset time/height variables

	m_currentHeight = 0;

	m_algoritmicTime = 0;

	//Clean IndexTable

	m_activeTable.clear();

}

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

inline void VectorizationContext::newSkeletonLink(unsigned int cur, ContourNode *node)

{

	if (node->hasAttribute(ContourNode::SK_NODE_DROPPED)) {

		SkeletonArc arcCopy(node);

		m_output->newLink(node->m_outputNode, cur, arcCopy);

		arcCopy.turn();

		m_output->newLink(cur, node->m_outputNode, arcCopy);

	}

}

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

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

//    Thinning Functions/Methods

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

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

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

//      Repair Polygon Degenerations

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

//EXPLANATION: After "Polygonizer", there may be simpleness degenerations

//about polygons which are dangerous to deal in the thinning process.

//Typically, these correspond to cases in which node->m_direction.z ~ 0

//(too fast), and are concave.

//We then deal with them *before* the process begins, by splitting one

//such node in two slower ones (known as 'Linear Axis' method).

inline void VectorizationContext::addLinearNodeBefore(ContourNode *node)

{

	ContourNode *newNode = getLinearNode();

	ContourEdge *newEdge = getLinearEdge();

	newNode->m_position = node->m_position;

	//Build new edge

	if (node->m_direction.z < 0.1)

		newEdge->m_direction = rotate270(node->m_edge->m_direction);

	else

		newEdge->m_direction = normalize(

			node->m_edge->m_direction + node->m_prev->m_edge->m_direction);

	newNode->m_edge = newEdge;

	//Link newNode

	newNode->m_prev = node->m_prev;

	newNode->m_next = node;

	node->m_prev->m_next = newNode;

	node->m_prev = newNode;

	//Build remaining infos

	node->buildNodeInfos(); //Rebuild

	newNode->buildNodeInfos();

	newNode->m_updateTime = 0;

	newNode->m_ancestor = node->m_ancestor;

	newNode->m_ancestorContour = node->m_ancestorContour;

	//Set node and newNode's edges not to be recognized as possible

	//opposites by the other (could happen in *future* instants)

	//                *DO NOT REMOVE!*

	node->m_notOpposites.push_back(newNode->m_edge);

	node->m_notOpposites.push_back(newNode->m_prev->m_edge);

	newNode->m_notOpposites.push_back(node->m_edge);

	//Further sign newly added node

	newNode->setAttribute(ContourNode::LINEAR_ADDED);

}

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

inline void VectorizationContext::repairDegenerations(const std::vector<contournode *=""> °enerates)</contournode>

{

	unsigned int i;

	m_linearNodesHeap.resize(degenerates.size());

	m_linearEdgesHeap.resize(degenerates.size());

	m_linearNodesHeapCount = 0;

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

		if (!degenerates[i]->hasAttribute(ContourNode::AMBIGUOUS_LEFT)) {

			addLinearNodeBefore(degenerates[i]);

			m_totalNodes++;

		}

	}

}

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

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

//    Node Infos Construction

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

inline void VectorizationContext::prepareContours(ContourFamily &family)

{

	std::vector<contournode *=""> degenerateNodes;</contournode>

	//Build circular links

	unsigned int i, j, k;

	unsigned int current;

	m_contoursCount = family.size();

	m_totalNodes = 0;

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

		for (j = 0, k = family[i].size() - 1; j < family[i].size(); k = j, ++j) {

			family[i][k].m_next = &family[i][j];

			family[i][j].m_prev = &family[i][k];

		}

		m_totalNodes += family[i].size();

	}

	//Build node edges

	m_edgesHeap.resize(m_totalNodes);

	current = 0;

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

		for (j = 0, k = family[i].size() - 1; j < family[i].size(); k = j, ++j) {

			m_edgesHeap[current].m_direction =

				normalize(planeProjection(family[i][j].m_position - family[i][k].m_position));

			family[i][k].m_edge = &m_edgesHeap[current];

			current++;

		}

	}

	bool maxThicknessNotZero = m_globals->currConfig->m_maxThickness > 0.0;

	//Now build remaining infos

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

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

			family[i][j].buildNodeInfos();

			family[i][j].m_updateTime = 0;

			family[i][j].m_ancestor = j;

			family[i][j].m_ancestorContour = i;

			//Check the following degeneration

			if (family[i][j].m_concave && family[i][j].m_direction.z < 0.3) {

				//Push this node among degenerate ones

				degenerateNodes.push_back(&family[i][j]);

			}

			//Insert output node in sharp angles

			if (!family[i][j].m_concave && family[i][j].m_direction.z < 0.6 && maxThicknessNotZero) {

				family[i][j].setAttribute(ContourNode::SK_NODE_DROPPED);

				family[i][j].m_outputNode = m_output->newNode(family[i][j].m_position);

			}

			//Push on nodes having AMBIGUOUS_RIGHT attribute

			if (family[i][j].hasAttribute(ContourNode::AMBIGUOUS_RIGHT))

				family[i][j].m_position += 0.02 * family[i][j].m_direction;

		}

	}

	//Finally, ensure polygon degenerations found are solved

	if (maxThicknessNotZero)

		repairDegenerations(degenerateNodes);

}

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

//WARNING: m_edge field of *this* and *previous* node must already be defined.

inline void ContourNode::buildNodeInfos(bool forceConvex)

{

	TPointD direction;

	double parameter;

	//Calculate node convexity

	if (forceConvex)

		m_concave = 0;

	else if (cross(m_edge->m_direction, m_prev->m_edge->m_direction) < 0) {

		m_concave = 1;

	} else

		m_concave = 0;

	//Build node direction

	direction = m_edge->m_direction - m_prev->m_edge->m_direction;

	parameter = norm(direction);

	if (parameter > 0.01) {

		direction = direction * (1 / parameter);

		if (m_concave)

			direction = -direction;

	} else

		direction = rotate270(m_edge->m_direction);

	m_direction.x = direction.x;

	m_direction.y = direction.y;

	//Calculate node speed

	m_direction.z = cross(planeProjection(m_direction), m_edge->m_direction);

	if (m_direction.z < 0)

		m_direction.z = 0;

	//Calculate angular momentum

	m_AngularMomentum = cross(m_position, m_direction);

	if (m_concave) {

		m_AuxiliaryMomentum1 = m_AuxiliaryMomentum2 = m_AngularMomentum;

	} else {

		m_AuxiliaryMomentum1 =

			cross(m_position, T3DPointD(m_edge->m_direction.y, -m_edge->m_direction.x, 1));

		m_AuxiliaryMomentum2 =

			cross(m_position, T3DPointD(m_prev->m_edge->m_direction.y, -m_prev->m_edge->m_direction.x, 1));

	}

}

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

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

//      Timeline Construction

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

//NOTE: In the following, we achieve these results:

//        * Build the timeline - events priority queue

//        * Process those split events which *necessarily* happen

//Pre-processing of split events is useful in order to lower execution times.

//Each node is first associated to a random integer; then a referencing

//vector is ordered according to those integers - events are calculated

//following this order. Split events are therefore calculated sparsely

//along the polygons, allowing a significant time reduction effect.

class RandomizedNode

{

public:

	ContourNode *m_node;

	int m_number;

	RandomizedNode() {}

	RandomizedNode(ContourNode *node) : m_node(node), m_number(rand()) {}

	inline ContourNode *operator->(void)

	{

		return m_node;

	}

};

class RandomizedNodeLess

{

public:

	RandomizedNodeLess() {}

	inline bool operator()(RandomizedNode a, RandomizedNode b)

	{

		return (a.m_number < b.m_number);

	}

};

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

void Timeline::build(ContourFamily &polygons, VectorizationContext &context,

					 VectorizerCore *thisVectorizer)

{

	unsigned int i, j, current;

	std::vector<randomizednode> nodesToBeTreated(context.m_totalNodes);</randomizednode>

	T3DPointD momentum, ray;

	//Build casual ordered node-array

	for (i = 0, current = 0; i < polygons.size(); ++i)

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

			nodesToBeTreated[current++] = RandomizedNode(&polygons[i][j]);

	//Same for linear-added nodes

	for (i = 0; i < context.m_linearNodesHeapCount; ++i)

		nodesToBeTreated[current++] = RandomizedNode(&context.m_linearNodesHeap[i]);

	double maxThickness = context.m_globals->currConfig->m_maxThickness;

	//Compute events generated by nodes

	//NOTE: are edge events to be computed BEFORE split ones?

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

		//Break calculation at user cancel press

		if (thisVectorizer->isCanceled())

			break;

		Event currentEvent(nodesToBeTreated[i].m_node, &context);

		//Notify event calculation

		if (!nodesToBeTreated[i].m_node->hasAttribute(ContourNode::LINEAR_ADDED))

			thisVectorizer->emitPartialDone();

		if (currentEvent.m_type != Event::failure && currentEvent.m_height < maxThickness)

#ifdef _PREPROCESS

			if (currentEvent.m_type == Event::split) {

				if (currentEvent.m_coGenerator->m_concave) {

					ray = T3DPointD(currentEvent.m_coGenerator->m_edge->m_direction.y,

									-currentEvent.m_coGenerator->m_edge->m_direction.x, 1);

					momentum =

						cross(currentEvent.m_coGenerator->m_position, ray);

					if (currentEvent.m_generator->m_direction * momentum +

							ray * currentEvent.m_generator->m_AngularMomentum <

						0) {

						timeline.push(currentEvent);

						continue;

					}

				}

				if (currentEvent.m_coGenerator->m_next->m_concave) {

					ray = T3DPointD(currentEvent.m_coGenerator->m_edge->m_direction.y,

									-currentEvent.m_coGenerator->m_edge->m_direction.x, 1);

					momentum =

						cross(currentEvent.m_coGenerator->m_next->m_position, ray);

					if (currentEvent.m_generator->m_direction * momentum +

							ray * currentEvent.m_generator->m_AngularMomentum >

						0) {

						timeline.push(currentEvent);

						continue;

					}

				}

				if (cross(currentEvent.m_generator->m_edge->m_direction,

						  currentEvent.m_coGenerator->m_edge->m_direction) > 0.02 &&

					cross(currentEvent.m_coGenerator->m_edge->m_direction,

						  currentEvent.m_generator->m_prev->m_edge->m_direction) > 0.02) // 0.02 in comparison with 'parameter' in buildNodeInfos

				{

					//Pre-processing succeeded

					currentEvent.process();

					continue;

				}

			}

#endif

		push(currentEvent);

	}

}

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

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

//      Event Calculation

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

//Calculates event generated by input node

Event::Event(ContourNode *generator, VectorizationContext *context)

	: m_height(infinity), m_displacement(infinity), m_generator(generator), m_type(failure), m_algoritmicTime(context->m_algoritmicTime), m_context(context)

{

	if (generator->m_concave)

		calculateSplitEvent();

	else

		calculateEdgeEvent();

}

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

// The edge event *generated by a node* is defined as the earliest edge event

// generated by its adjacent edges. Remember that 'edge events' correspond to

// those in which one edge gets 0 length.

inline void Event::calculateEdgeEvent()

{

	struct locals {

		static inline void buildDisplacements(ContourNode *edgeFirst, double &d1, double &d2)

		{

			ContourNode *edgeSecond = edgeFirst->m_next;

			// If bisectors are almost opposite, avoid: there must be another bisector

			// colliding with m_generator *before* coGenerator - allowing a positive

			// result here may interfere with it.

			if ((edgeFirst->m_concave && edgeSecond->m_concave) || edgeFirst->m_direction * edgeSecond->m_direction < -0.9) {

				d1 = d2 = -1.0;

				return;

			}

			double det = edgeFirst->m_direction.y * edgeSecond->m_direction.x - edgeFirst->m_direction.x * edgeSecond->m_direction.y;

			double cx = edgeSecond->m_position.x - edgeFirst->m_position.x,

				   cy = edgeSecond->m_position.y - edgeFirst->m_position.y;

			d1 = (edgeSecond->m_direction.x * cy - edgeSecond->m_direction.y * cx) / det;

			d2 = (edgeFirst->m_direction.x * cy - edgeFirst->m_direction.y * cx) / det;

		}

		static inline double height(ContourNode *node, double displacement)

		{

			return node->m_position.z + displacement * node->m_direction.z;

		}

	}; // locals

	double minHeight, minDisplacement;

	bool positiveEdgeDispl;

	m_type = edge;

	// Calculate the two possible displacement parameters

	double firstDisplacement, prevDisplacement, nextDisplacement, lastDisplacement;

	// right == prev

	locals::buildDisplacements(m_generator, nextDisplacement, lastDisplacement);

	locals::buildDisplacements(m_generator->m_prev, firstDisplacement, prevDisplacement);

	// Take the smallest positive between them and assign the co-generator

	// NOTE: In a missed vertex event, the threshold value is to compare with the possible pushes at the end of processSplit

	// However, admitting slightly negative displacements should be ok: due to the weak linear axis imposed for concave

	// vertices, it is impossible to have little negative displacements apart from the above mentioned pushed case.

	//   ..currently almost true..

	static const double minusTol = -0.03;

	bool prevDispPositive = (prevDisplacement > minusTol);

	bool nextDispPositive = (nextDisplacement > minusTol);

	if (nextDispPositive) {

		if (!prevDispPositive || nextDisplacement < prevDisplacement) {

			m_coGenerator = m_generator;

			minDisplacement = nextDisplacement;

			minHeight = locals::height(m_coGenerator, nextDisplacement);

			positiveEdgeDispl = (nextDispPositive && lastDisplacement > minusTol);

		} else {

			m_coGenerator = m_generator->m_prev;

			minDisplacement = prevDisplacement;

			minHeight = locals::height(m_coGenerator, firstDisplacement);			// Height is built on the edge's first

			positiveEdgeDispl = (prevDispPositive && firstDisplacement > minusTol); // endpoint to have the same values on adjacent

		}																			// generators. It's important for SPECIAL events.

	} else if (prevDispPositive) {

		m_coGenerator = m_generator->m_prev;

		minDisplacement = prevDisplacement;

		minHeight = locals::height(m_coGenerator, firstDisplacement); // Same here

		positiveEdgeDispl = (prevDispPositive && firstDisplacement > minusTol);

	} else {

		m_type = failure;

		return;

	}

	// NOTA: Le const di tolleranza sono da confrontare tra:

	// a) i push alla fine di processSplit per evitare rogne coi vertex multipli

	// b) Le condizioni di esclusione su side1 e side2 in trySplit che evitano a)

	// c) Le condizioni di riconoscimento di vertex e special events - perche' se mancano...

	// Special cases: (forse da raffinare le condizioni - comunque ora sono efficaci)

	if (nextDispPositive && !m_generator->m_concave) {

		if (m_generator->m_prev->m_concave && m_generator->m_next->m_concave &&

			fabs(nextDisplacement - prevDisplacement) < 0.1) //condizione debole per escludere subito i casi evidentemente innocenti

		{

			// Check 'V' (special) event - can generate a new concave vertex

			ContourNode *prevRay = m_generator->m_prev,

						*nextRay = m_generator->m_next;

			double side = prevRay->m_direction * nextRay->m_AngularMomentum +

						  nextRay->m_direction * prevRay->m_AngularMomentum;

			// NOTE: fabs(side) / || prevRay->dir x nextRay->dir ||  is the distance between the two rays.

			if (fabs(side) < 0.03 * norm(cross(prevRay->m_direction, nextRay->m_direction)))

				m_type = special, m_coGenerator = m_generator;

		} else if (fabs(nextDisplacement - prevDisplacement) < 0.01) {

			// Then choose to make the event with a concave vertex (to resemble the 'LL' case)

			m_coGenerator = m_generator->m_next->m_concave ? m_generator : m_generator->m_prev;

		}

	}

	// Now, if calculated height is coherent, this Event is valid.

	if (positiveEdgeDispl										// Edges shrinking to a point after a FORWARD

		|| minHeight > m_context->m_currentHeight - 0.01)		// displacement are processable - this dominates

		m_height = minHeight, m_displacement = minDisplacement; // height considerations which may be affected by

	else														// numerical errors

		m_type = failure;

}

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

inline void Event::calculateSplitEvent()

{

	unsigned int i;

	bool forceFirst;

	ContourNode *opposite, *first, *last;

	std::list<contournode *="">::iterator currentContour;</contournode>

	// Sign *edges* not to be taken as possible opposites

	for (i = 0; i < m_generator->m_notOpposites.size(); ++i)

		m_generator->m_notOpposites[i]->setAttribute(ContourEdge::NOT_OPPOSITE);

	// Check adjacent edge events

	calculateEdgeEvent(); // DO NOT REMOVE - adjacent convexes may have

						  // been calculated too earlier

	// First check opposites in the m_generator active contour

	first = m_generator->m_next->m_next; // Adjacent edges were already considered

	last = m_generator->m_prev->m_prev;  // by calculateEdgeEvents()

	for (opposite = first; opposite != last; opposite = opposite->m_next) {

		if (!opposite->m_edge->hasAttribute(ContourEdge::NOT_OPPOSITE))

			tryRayEdgeCollisionWith(opposite);

	}

	IndexTable &activeTable = m_context->m_activeTable;

	//Then, try in the remaining active contours whose identifier is != our

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

		for (currentContour = activeTable[i]->begin();

			 currentContour != activeTable[i]->end(); currentContour++) {

			//Da spostare sopra il 2o for

			if (activeTable.m_identifiers[(*currentContour)->m_ancestorContour] !=

				activeTable.m_identifiers[m_generator->m_ancestorContour]) {

				first = *currentContour;

				for (opposite = first, forceFirst = 1;

					 //Better the first next cond. - in case of thinning errors, at least it does not get loop'd.

					 !opposite->hasAttribute(ContourNode::HEAD) //opposite!=first