/* === S Y N F I G ========================================================= */

/*!	\file centerlineskeletonizer.cpp

**

**	$Id$

**

**	\legal

**	Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley

**

**	This package is free software; you can redistribute it and/or

**	modify it under the terms of the GNU General Public License as

**	published by the Free Software Foundation; either version 2 of

**	the License, or (at your option) any later version.

**

**	This package is distributed in the hope that it will be useful,

**	but WITHOUT ANY WARRANTY; without even the implied warranty of

**	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

**	General Public License for more details.

**	\endlegal

*/

/* ========================================================================= */

/* === H E A D E R S ======================================================= */

#include "polygonizerclasses.h"

#include <queue></queue>

#include <synfig vector.h=""></synfig>

/* === U S I N G =========================================================== */

using namespace std;

using namespace etl;

using namespace studio;

using namespace synfig;

//<---------------------------Some Useful functions----------------------------->

inline double cross(const synfig::Point &a, const synfig::Point &b) 

{

  return a[0] * b[1] - a[1] * b[0];

}

inline synfig::Point3 cross(const synfig::Point3 &a, const synfig::Point3 &b) {

  return synfig::Point3(a[1] * b[2] - b[1] * a[2], a[2] * b[0] - b[2] * a[0],

                      a[0] * b[1] - b[0] * a[1]);

}

inline bool angleLess(const synfig::Point &a, const synfig::Point &b) {

  return a[1] >= 0 ? b[1] >= 0 ? a[0] > b[0] : 1 : b[1] < 0 ? a[0] < b[0] : 0;

}

// a, b, ref assumed normalized

inline bool angleLess(const synfig::Point &a, const synfig::Point &b, const synfig::Point &ref) {

  return angleLess(a, ref) ? angleLess(b, ref) ? angleLess(a, b) : 0

                           : angleLess(b, ref) ? 1 : angleLess(a, b);

}

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

//      Classes

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

struct VectorizationContext;

class studio::ContourEdge 

{

public:

  enum { NOT_OPPOSITE = 0x1 };

public:

  synfig::Point m_direction;

  unsigned short m_attributes;

public:

  ContourEdge() : m_attributes(0) {}

  ContourEdge(synfig::Point 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></indexcolumn>

      m_columns;  //!< Contours set by 'column identifier'.

  std::vector<int></int>

      m_identifiers;  //!< Column identifiers by original contour index.

  // 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_algorithmicTime;

  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(studio::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 final

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

public:

  Timeline() {}

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

  // vectorization;

  void build(ContourFamily &polygons, VectorizationContext &context);

};

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

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

//    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_algorithmicTime; // Number of events processed 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];

  }

  studio::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_algorithmicTime = 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[2] ~ 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[2] < 0.1)

    newEdge->m_direction = (node->m_edge->m_direction).perp();

  else

    newEdge->m_direction = (node->m_edge->m_direction + node->m_prev->m_edge->m_direction).norm();

  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 = ((family[i][j].m_position - family[i][k].m_position).to_2d()).norm();

      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[2] < 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[2] < 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 += (family[i][j].m_direction) * 0.02;

    }

  }

  // 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) {

  synfig::Point 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 = direction.mag();

  if (parameter > 0.01) 

  {

    direction                = direction * (1 / parameter);

    if (m_concave) direction = -direction;

  }

  else

    direction = (m_edge->m_direction).perp();

  m_direction[0] = direction[0];

  m_direction[1] = direction[1];

  // Calculate node speed

  m_direction[2] = cross(m_direction.to_2d(), m_edge->m_direction);

  if (m_direction[2] < 0) m_direction[2] = 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,

              synfig::Point3(m_edge->m_direction[1], -(m_edge->m_direction[0]), 1));

    m_AuxiliaryMomentum2 =

        cross(m_position, synfig::Point3(m_prev->m_edge->m_direction[1],

                                    -(m_prev->m_edge->m_direction[0]), 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) {

  unsigned int i, j, current;

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

  synfig::Point3 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) 

  {

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

    // Notify event calculation

    if (currentEvent.m_type != Event::failure &&

        currentEvent.m_height < maxThickness)

    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_algorithmicTime(context->m_algorithmicTime)

    , 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[1] * edgeSecond->m_direction[0] -

                   edgeFirst->m_direction[0] * edgeSecond->m_direction[1];

      double cx = edgeSecond->m_position[0] - edgeFirst->m_position[0],

             cy = edgeSecond->m_position[1] - edgeFirst->m_position[1];

      d1 = (edgeSecond->m_direction[0] * cy - edgeSecond->m_direction[1] * cx) /

           det;

      d2 =

          (edgeFirst->m_direction[0] * cy - edgeFirst->m_direction[1] * cx) / det;

    }

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

      return node->m_position[2] + displacement * node->m_direction[2];

    }

  };  // 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;

  }

  if (nextDispPositive && !m_generator->m_concave) 

  {

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

        fabs(nextDisplacement - prevDisplacement) < 0.1)  

    {

      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;

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

        m_type = special, m_coGenerator = m_generator;

    } 

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

    {

      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; !opposite->hasAttribute(ContourNode::HEAD)  // opposite!=first

             || (forceFirst ? forceFirst = 0, 1 : 0); opposite = opposite->m_next) 

        {

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

            tryRayEdgeCollisionWith(opposite);

        }

      }

    }

  }

  // Restore edge attributes

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

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

}

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

inline bool Event::testRayEdgeCollision(ContourNode *opposite,

                                        double &displacement, double &height,

                                        double &side1, double &side2) {

  // Initialize test vectors

  // NOTE: In the convex case, slab guards MUST be orthogonal to the edge, due

  // to this case:

  //

  //    ______/|                the ray would not hit the edge - AND THUS FOREGO

  //    INTERACTION

  //           |                WITH IT COMPLETELY

  //     ->    |

  synfig::Point3 firstSlabGuard =

      opposite->m_concave ? opposite->m_direction

                          : synfig::Point3(opposite->m_edge->m_direction[1],

                                      -(opposite->m_edge->m_direction[0]), 1);

  synfig::Point3 lastSlabGuard =

      opposite->m_next->m_concave

          ? opposite->m_next->m_direction

          : synfig::Point3(opposite->m_edge->m_direction[1],

                      -(opposite->m_edge->m_direction[0]), 1);

  synfig::Point3 roofSlabOrthogonal(-(opposite->m_edge->m_direction[1]),

                               opposite->m_edge->m_direction[0], 1);