// TnzCore includes

#include "tmathutil.h"

#include "tcurves.h"

#include "tbezier.h"

#include "tstrokedeformations.h"

#include "tstroke.h"

#include "tcurveutil.h"

#include "tcg_wrap.h"

// tcg includes

#include "tcg/tcg_poly_ops.h"

#define INCLUDE_HPP

#include "tcg/tcg_polylineops.h"

#include "tcg/tcg_cyclic.h"

#undef INCLUDE_HPP

#include "tstrokeutil.h"

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

//    Local namespace  stuff

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

namespace {

typedef std::vector<tthickcubic *=""> TThickCubicArray;</tthickcubic>

typedef std::vector<tthickquadratic *=""> QuadStrokeChunkArray;</tthickquadratic>

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

int getControlPointIndex(const TStroke &stroke, double w) {

  TThickPoint p = stroke.getControlPointAtParameter(w);

  int i                 = 0;

  int controlPointCount = stroke.getControlPointCount();

  for (; i < controlPointCount; ++i)

    if (stroke.getControlPoint(i) == p) return i;

  return controlPointCount - 1;

}

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

double findMinimum(const TStrokeDeformation &def, const TStroke &stroke,

                   double x1, double x2, double xacc, double length = 0,

                   int max_iter = 100)

{

  int j;

  double dx, f, fmid, xmid, rtb;

  f    = def.getDelta(stroke, x1) - length;

  fmid = def.getDelta(stroke, x2) - length;

  if (f == 0) return x1;

  if (fmid == 0) return x2;

  if (f * fmid > 0.0) return -1;

  rtb = f < 0.0 ? (dx = x2 - x1, x1) : (dx = x1 - x2, x2);

  for (j = 1; j <= max_iter; j++) {

    fmid = def.getDelta(stroke, xmid = rtb + (dx *= 0.5)) - length;

    if (fmid <= 0.0) rtb = xmid;

    if (fabs(dx) < xacc || fmid == 0.0) return rtb;

  }

  return -2;

}

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

/**

 * Rationale:

 *  Suppose we want to model a segment (represented by a stroke)  so that it

 * takes the shape of a parabola (the usual case offered by the modifier). We

 * assume that: (o) stroke points lie along the y=-100 axis; (o) the x's that

 * will correspond are x1=-10 and x2=+10 (obvious from the equation).

 *

 *  The parabola may be represented on the left side by a quadratic with control

 * points: P0=(-10,-100), P1=(-5,    0), P2=( 0,    0). If we know the number of

 * linear strokes representing this parabola, we also know how many "samples"

 * are required for its linearization. This parameter can be used to

 * qualitatively determine the value with which to sample the stroke to be

 * tested; there will need to be as many points to move as there are samples in

 * the reference.

 */

double computeIncrement(double strokeLength, double pixelSize) {

  assert(pixelSize > 0 && "Pixel size is negative!!!");

  assert(strokeLength > 0 && "Stroke Length size is negative!!!");

  // height of the parabola (goes downward)

  double height = 100;

  // I suppose I'm doing at least a 100-pixel drag

  assert(height >= 100.0);

  double x = sqrt(height);

  // the point p1 will have to be at the intersection of the tangents to the two

  // extremes. The tangent of the point p2 and the x-axis, the other will have

  // versor given by the gradient at p0,

  //  ie: grad(x,-2 x)

  //  and if y = m x + q

  //  m =

  double m = 2.0 * x;

  double q = m * x - height;

  double p1x = q / m;

  double scale = strokeLength / (2.0 * x);

  TScale scaleAffine(scale, scale);

  TPointD p0 = scaleAffine * TPointD(-x, -height),

          p1 = scaleAffine * TPointD(-p1x, 0.0),

          p2 = scaleAffine * TPointD(0.0, 0.0);

  TQuadratic quadratic(p0, p1, p2);

  double step = computeStep(quadratic, pixelSize);

  //  just to add points even in the worst case.

  if (step >= 1.0) step = 0.1;

  return step;

}

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

void detectEdges(const std::vector<tpointd> &pointArray,</tpointd>

                 std::vector<uint> &edgeIndexArray) {</uint>

  // ASSUMPTION: sharpPointArray does not contain adjacent coincident points

  int size = pointArray.size();

  // I check that there are more than three elements

  if (size < 3) return;

  //  runs pointArray and for each of its points tries to inscribe triangles

  //  (using left and right points) considering potential corners those with

  //  sides l such that dMin <= l <= dMax (actually at the first time that l >

  //  dMax: breack) and with angular aperture alpha <= alphaMax.

  // Then it looks for local maxes among the potential corners in a window of

  // semiamplitude dMax(actually at the first time dMax : breack is exceeded)

  //  default values: dMin = 7; dMax = dMin + 2; alphaMax = 2.6 (150 degrees)

  const double dMin     = 4;

  const double dMax     = dMin + 3;

  const double alphaMax = 2.4;  // ( 137.5 degrees)

  const double dMin2    = dMin * dMin;

  const double dMax2    = dMax * dMax;

  std::vector<double> sharpnessArray;</double>

  sharpnessArray.push_back(M_PI);  //  the first point is a corner

  int nodeCount;

  for (nodeCount = 1; nodeCount < size - 1;

       ++nodeCount) {  //  scrolls the sharpPointArray excluding the extremes

    sharpnessArray.push_back(0);

    TPointD point(pointArray[nodeCount]);

    int leftCount;

    for (leftCount = nodeCount - 1; leftCount >= 0;

         --leftCount) {  //  Calculates the "left" sides of the inscribed

                         //  triangles...

      TPointD left  = pointArray[leftCount];

      double dLeft2 = norm2(left - point);

      if (dLeft2 < dMin2)

        continue;

      else if (dLeft2 > dMax2)

        break;

      int rightCount;

      for (rightCount = nodeCount + 1; rightCount < size;

           ++rightCount) {  // Calculates the "right" sides of the inscribed

                            // triangles...

        TPointD right  = pointArray[rightCount];

        double dRight2 = norm2(right - point);

        if (dRight2 < dMin2)

          continue;

        else if (dMax2 < dRight2)

          break;

        //  Calculates the "center" sides of the inscribed triangles

        double dCenter2 = norm2(left - right);

        assert(dLeft2 != 0.0 && dRight2 != 0.0);

        double cs =

            (dLeft2 + dRight2 - dCenter2) / (2 * sqrt(dLeft2 * dRight2));

        double alpha = acos(cs);

        if (alpha > alphaMax) continue;

        double sharpness = M_PI - alpha;

        if (sharpnessArray[nodeCount] < sharpness)

          sharpnessArray[nodeCount] = sharpness;

      }

    }

  }

  edgeIndexArray.push_back(0);  // the first point is a corner

  // I find local maxima by excluding extremes

  for (nodeCount = 1; nodeCount < size - 1;

       ++nodeCount) {  // scroll through the list excluding the extremes

    bool isCorner = true;

    TPointD point(pointArray[nodeCount]);

    int leftCount;

    for (leftCount = nodeCount - 1; leftCount >= 0;

         --leftCount) {  //  scrolls down the list of sharpPoints to the left of

                         //  node...

      TPointD left  = pointArray[leftCount];

      double dLeft2 = norm2(left - point);

      if (dLeft2 > dMax2) break;

      if (sharpnessArray[leftCount] > sharpnessArray[nodeCount]) {

        isCorner = false;

        break;

      }

    }

    if (isCorner) continue;

    int rightCount;

    for (rightCount = nodeCount + 1; rightCount < size;

         ++rightCount) {  // scrolls the list of sharpPoints to the right of

                          // node..

      TPointD right  = pointArray[rightCount];

      double dRight2 = norm2(right - point);

      if (dRight2 > dMax2) break;

      if (sharpnessArray[rightCount] > sharpnessArray[nodeCount]) {

        isCorner = false;

        break;

      }

    }

    if (isCorner) edgeIndexArray.push_back(nodeCount);

  }

  edgeIndexArray.push_back(size - 1);  //  the last point is a corner

}

}  // namespace

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

//    API  functions

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

bool increaseControlPoints(TStroke &stroke, const TStrokeDeformation &deformer,

                           double pixelSize) {

  if (isAlmostZero(stroke.getLength())) {

    return norm2(deformer.getDisplacement(stroke, 0.0)) > 0;

  }

  // step 1:

  // It's possible to have control point at not null potential

  //  but with delta equal 0 (equipotential control point)

  bool notVoidPotential = false;

  for (int i = 0; i < stroke.getControlPointCount(); ++i) {

    double par = stroke.getParameterAtControlPoint(i);

    if (deformer.getDisplacement(stroke, par) != TThickPoint()) {

      notVoidPotential = true;

      break;

    }

  }

  // step 2:

  //  increase control point checking delta of deformer

  double maxDifference =

      deformer.getMaxDiff();  // above this delta value, points are added

  int strokeControlPoint = stroke.getControlPointCount();

  // pixelSize = sq( pixelSize );

  if (pixelSize < TConsts::epsilon) pixelSize = TConsts::epsilon;

  double length = stroke.getLength(),

         // set the step function of length

      //    step = length > 1.0 ?  pixelSize * 15.0/ length : length,

      // step = 0.01,

      w = 0.0;

  double step = computeIncrement(length, pixelSize);

  double x1, x2, d1, d2, diff, offset, minimum, incr;

  incr = step;

  while (w + incr < 1.0) {

    d1 = deformer.getDelta(stroke, w);

    d2 = deformer.getDelta(stroke, w + incr);

    diff = d2 - d1;

    if (fabs(diff) >= maxDifference)  // if there is a step of potential

    {

      if (tsign(diff) > 0) {

        x1 = w;

        x2 = w + incr;

      } else {

        x1 = w + incr;

        x2 = w;

      }

      offset = (d1 + d2) * 0.5;

      // find the position of step

      minimum = findMinimum(

          deformer, stroke, x1, x2, TConsts::epsilon, offset,

          20);  // A new control point should be put between x1 and x2. where?

      // this function finds the point at which the maxdifference value is

      // exceeded

      // if minimum is not found or is equal to previous value

      //  use an heuristic...

      if (minimum < 0 || w == minimum) {

        minimum = w + incr * 0.5;

        w += step;

      }

      //... else insert a control point in minimum

      w = minimum;  // scanning resumes from the new point, in this way it

                    // thickens ...

      stroke.insertControlPoints(minimum);

      // update of step

      incr = step;

    } else

      incr += step;

  }

  // return true if control point are increased

  return (stroke.getControlPointCount() > strokeControlPoint) ||

         notVoidPotential;

}

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

void modifyControlPoints(TStroke &stroke, const TStrokeDeformation &deformer) {

  int cpCount = stroke.getControlPointCount();

  TThickPoint newP;

  for (int i = 0; i < cpCount; ++i) {

    newP = stroke.getControlPoint(i) +

           deformer.getDisplacementForControlPoint(stroke, i);

    if (isAlmostZero(newP.thick, 0.005)) newP.thick = 0;

    stroke.setControlPoint(i, newP);

  }

}

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

void modifyControlPoints(TStroke &stroke, const TStrokeDeformation &deformer,

                         std::vector<double> &controlPointLen) {</double>

  UINT cpCount = stroke.getControlPointCount();

  TThickPoint newP;

#ifdef _DEBUG

  UINT debugVariable = controlPointLen.size();

#endif

  assert(controlPointLen.size() == cpCount);

  for (UINT i = 0; i < cpCount; ++i) {

    newP =

        stroke.getControlPoint(i) +

        deformer.getDisplacementForControlPointLen(stroke, controlPointLen[i]);

    if (isAlmostZero(newP.thick, 0.005)) newP.thick = 0;

    stroke.setControlPoint(i, newP);

  }

}

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

void modifyThickness(TStroke &stroke, const TStrokeDeformation &deformer,

                     std::vector<double> &controlPointLen, bool exponentially) {</double>

  UINT cpCount = stroke.getControlPointCount();

  assert(controlPointLen.size() == cpCount);

  double disp;

  double thick;

  for (UINT i = 0; i < cpCount; ++i) {

    disp =

        (deformer.getDisplacementForControlPointLen(stroke, controlPointLen[i]))

            .thick;

    thick = stroke.getControlPoint(i).thick;

    // The additive version is straightforward.

    // The exponential version is devised to keep derivative 1 at disp == 0;

    // it is typically used when the thickness decreases.

    thick = (exponentially && thick >= 0.005) ? thick * exp(disp / thick)

                                              : thick + disp;

    if (thick < 0.005) thick = 0.0;

    stroke.setControlPoint(i, TThickPoint(stroke.getControlPoint(i), thick));

  }

}

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

void transform_thickness(TStroke &stroke, const double poly[], int deg) {

  int cp, cpCount = stroke.getControlPointCount();

  for (cp = 0; cp != cpCount; ++cp) {

    TThickPoint cpPoint = stroke.getControlPoint(cp);

    cpPoint.thick =

        std::max(tcg::poly_ops::evaluate(poly, deg, cpPoint.thick), 0.0);

    stroke.setControlPoint(cp, cpPoint);

  }

}

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

TStroke *Toonz::merge(const std::vector<tstroke *=""> &strokes) {</tstroke>

  if (strokes.empty()) return 0;

  std::vector<tthickpoint> new_stroke_cp;</tthickpoint>

  int size_stroke_array = strokes.size();

  int size_cp;

  const TStroke *ref;

  TThickPoint last = TConsts::natp;

  if (!strokes[0]) return 0;

  new_stroke_cp.push_back(strokes[0]->getControlPoint(0));

  int i, j;

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

    ref = strokes[i];

    if (!ref) return 0;

    size_cp = ref->getControlPointCount();

    for (j = 0; j < size_cp - 1; j++) {

      const TThickPoint &pnt = ref->getControlPoint(j);

      if (last != TConsts::natp && j == 0) {

        // new_stroke_cp.push_back( (last+pnt)*0.5 );

        new_stroke_cp.push_back(last);

      }

      if (j > 0) new_stroke_cp.push_back(pnt);

    }

    // last point needs to be merged

    last = ref->getControlPoint(size_cp - 1);

  }

  new_stroke_cp.push_back(ref->getControlPoint(size_cp - 1));

  TStroke *out = new TStroke(new_stroke_cp);

  return out;

}

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

namespace {

class CpsReader {

  std::vector<tthickpoint> &m_cps;</tthickpoint>

public:

  typedef TPointD value_type;

public:

  CpsReader(std::vector<tthickpoint> &cps) : m_cps(cps) {}</tthickpoint>

  void openContainer(const TPointD &point) { addElement(point); }

  void addElement(const TPointD &point) {

    m_cps.push_back(TThickPoint(point, 0.0));

  }

  void closeContainer() {}

};

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

//    Triplet to Quadratics

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

template <typename iter_type=""></typename>

double buildLength(const iter_type &begin, const iter_type &end, double tol) {

  // Build direction

  iter_type it = begin, jt;

  ++it;

  const TPointD &a = *begin, &b = *it;

  TPointD dir(normalize(b - a)), segDir;

  double dist;

  for (jt = it, ++it; it != end; jt = it, ++it) {

    segDir = *it - *jt;

    if (dir * segDir < 0) break;

    dist = tcg::point_ops::lineSignedDist(*it, a, dir);

    if (fabs(dist) > tol) {

      double s, t;

      if (dist > 0) {

        tcg::point_ops::intersectionCoords(

            *jt, segDir, a + tol * tcg::point_ops::ortLeft(dir), dir, s, t);

      } else {

        tcg::point_ops::intersectionCoords(

            *jt, segDir, a + tol * tcg::point_ops::ortRight(dir), dir, s, t);

      }

      s = tcrop(s, 0.0, 1.0);

      return (*jt + s * segDir - a) * dir;

    }

  }

  return (*jt - a) * dir;

}

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

/*

  Converts the specified points triplet into a sequence of quadratics' CPs

  (point

  a is not included, whereas c is).

  Conversion takes 4 parameters:

   - Adherence:     How much quadratics bend toward corners

   - Angle:         Inner product of corner's edges - full corners threshold

   - Relative:      Curvature radius/edge length    - full corners threshold

   - RelativeDist:  Tolerance about edge length build-ups

  See below for extended explanation.

*/

class TripletsConverter {

  typedef std::vector<tpointd>::const_iterator iter_type;</tpointd>

  typedef std::reverse_iterator<iter_type> riter_type;</iter_type>

  typedef tcg::cyclic_iterator<iter_type> cyclic_iter_type;</iter_type>

  typedef std::reverse_iterator<cyclic_iter_type> rcyclic_iter_type;</cyclic_iter_type>

  bool m_circular;

  iter_type m_first, m_end, m_last;

  double m_adherenceTol, m_angleTol, m_relativeTol, m_relativeDistTol;

public:

  TripletsConverter(const iter_type &begin, const iter_type &end,

                    double adherenceTol, double angleTol, double relativeTol,

                    double relativeDistTol)

      : m_circular(*begin == *(end - 1))

      , m_first(m_circular ? begin + 1 : begin)

      , m_end(end)

      , m_adherenceTol(adherenceTol)

      , m_angleTol(angleTol)

      , m_relativeTol(relativeTol)

      , m_relativeDistTol(relativeDistTol) {}

  // Using bisector to convert a triplet

  void operator()(const TPointD &a, const iter_type &bt, const TPointD &c,

                  tcg::sequential_reader<std::vector<tpointd>> &output) {</std::vector<tpointd>

    const TPointD &b = *bt;

    double prod =

        tcg::point_ops::direction(b, a) * tcg::point_ops::direction(b, c);

    if (prod > m_angleTol) {

      // Full corner

      output.addElement(0.5 * (a + b));

      output.addElement(b);

      output.addElement(0.5 * (b + c));

    } else {

      // Build the angle bisector

      TPointD a_b(a - b);

      TPointD c_b(c - b);

      double norm_a_b = norm(a_b);

      double norm_c_b = norm(c_b);

      a_b = a_b * (1.0 / norm_a_b);

      c_b = c_b * (1.0 / norm_c_b);

      TPointD v(tcg::point_ops::normalized(a_b + c_b));

      double cos_v_dir = fabs(a_b * v);

      double t1 = tcrop(m_adherenceTol / (cos_v_dir * norm_a_b), 0.0, 0.5);

      double t2 = tcrop(m_adherenceTol / (cos_v_dir * norm_c_b), 0.0, 0.5);

      if (t1 == 0.5 && t2 == 0.5) {

        // Direct conversion

        output.addElement(b);

      } else {

        // Build the quadratic split

        TPointD d(b + t1 * (a - b)), f(b + t2 * (c - b)), e(0.5 * (d + f));

        // Build curvature radiuses at the corner

        // NOTE: Both speed and acceleration would hold 2.0 as multiplier, which

        // is calculated implicitly.

        TPointD speed(f - d);

        double num = norm(speed);

        if (num <= TConsts::epsilon) {

          // Curvature radius is 0 - full corner

          output.addElement(0.5 * (a + b));

          output.addElement(b);

          output.addElement(0.5 * (b + c));

        } else {

          num = 2.0 * num * num *

                num;  // would be * 8 = 2^3, divided by the 4 below

          double den1 =

              fabs(cross(speed, a - d));  // * 4, from both args of the cross

          double den2 = fabs(cross(speed, c - f));

          double radius1 = (den1 == 0.0) ? 0.0 : num / den1;

          double radius2 = (den1 == 0.0) ? 0.0 : num / den2;

          // Build edges length

          double length1, length2;

          if (m_circular) {

            cyclic_iter_type it(bt, m_first, m_end, 0);

            cyclic_iter_type it1(bt, m_first, m_end, 1);

            cyclic_iter_type it_1(bt, m_first, m_end, -1);

            rcyclic_iter_type rit(it + 1), rit1(it_1 + 1);

            length1 = buildLength(rit, rit1, 0.25);

            length2 = buildLength(it, it1, 0.25);

          } else {

            riter_type rit(bt + 1), rend(m_first);

            length1 = buildLength(rit, rend, m_relativeDistTol);

            length2 = buildLength(bt, m_end, m_relativeDistTol);

          }

          // Test curvature radiuses against edge length

          if (radius1 / length1 < m_relativeTol &&  // both must hold

              radius2 / length2 < m_relativeTol) {

            // Full corner

            output.addElement(0.5 * (a + b));

            output.addElement(b);

            output.addElement(0.5 * (b + c));

          } else {

            // Quadratic split

            output.addElement(d);

            output.addElement(e);

            output.addElement(f);

          }

        }

      }

    }

    output.addElement(c);

  }

};

}  // namespace

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

void polylineToQuadratics(const std::vector<tpointd> &polyline,</tpointd>

                          std::vector<tthickpoint> &cps, double adherenceTol,</tthickpoint>

                          double angleTol, double relativeTol,

                          double relativeDistTol, double mergeTol) {

  CpsReader cpsReader(cps);

  TripletsConverter op(polyline.begin(), polyline.end(), adherenceTol, angleTol,

                       relativeTol, relativeDistTol);

  tcg::polyline_ops::toQuadratics(polyline.begin(), polyline.end(), cpsReader,

                                  op, mergeTol);

}