#include "tinbetween.h"
#include "tcurves.h"
#include "tstroke.h"
#include "tpalette.h"
//#include "tcolorstyles.h"
//#include "tregion.h"
#include "tmathutil.h"
//#include "tstrokeutil.h"
//#include "tsystem.h"
#include <utility>
#include <limits>
#include <list>

#include "../tvectorimage/tvectorimageP.h"

//#include "tdebugmessage.h"
//--------------------------------------------------------------------------------------

double average(std::vector<double> &values, double range = 2.5) {
  UINT size = values.size();
  if (size == 0) return std::numeric_limits<double>::signaling_NaN();

  if (size == 1) return values[0];

  double sum = 0;
  UINT j     = 0;
  for (; j < size; j++) {
    sum += values[j];
  }

  double average = sum / size;

  double variance = 0;
  for (j = 0; j < size; j++) {
    variance += (average - values[j]) * (average - values[j]);
  }
  variance /= size;

  double err;
  int acceptedNum = 0;
  sum             = 0;
  for (j = 0; j < size; j++) {
    err = values[j] - average;
    err *= err;
    if (err <= range * variance) {
      sum += values[j];
      acceptedNum++;
    }
  }

  assert(acceptedNum > 0);
  return (acceptedNum > 0) ? sum / (double)acceptedNum : average;
}

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

double weightedAverage(std::vector<double> &values,
                       std::vector<double> &weights, double range = 2.5) {
  UINT size = values.size();
  if (size == 0) return std::numeric_limits<double>::signaling_NaN();

  double sum         = 0;
  double totalWeight = 0;
  UINT j             = 0;
  for (; j < size; j++) {
    sum += weights[j] * values[j];
    totalWeight += weights[j];
  }

  assert(totalWeight > 0);
  if (totalWeight == 0) return std::numeric_limits<double>::signaling_NaN();

  double average = sum / totalWeight;

  double variance = 0;
  for (j = 0; j < size; j++) {
    variance += weights[j] * (average - values[j]) * (average - values[j]);
  }
  variance /= totalWeight;

  double err;
  totalWeight = 0;
  sum         = 0;
  for (j = 0; j < size; j++) {
    err = values[j] - average;
    err *= err;
    if (err <= range * variance) {
      sum += weights[j] * values[j];
      totalWeight += weights[j];
    }
  }

  assert(totalWeight > 0);
  return (totalWeight != 0) ? sum / totalWeight : average;
}

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

inline UINT angleNumber(const std::vector<std::pair<int, double>> &corners,
                        double angle) {
  UINT count = 0;
  UINT size  = corners.size();
  for (UINT j = 0; j < size; j++)
    if (corners[j].second >= angle) count++;

  return count;
}

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

inline bool isTooComplex(UINT n1, UINT n2, UINT maxSubsetNumber = 100) {
  UINT n, r;
  if (n1 > n2) {
    n = n1;
    r = n2;
  } else if (n1 < n2) {
    n = n2;
    r = n1;
  } else {
    assert(!"equal angle number");
    return false;
  }

  if (n - r < r) r = n - r;

  /*

n*(n-1)* ...(n-r+1) < n^r that must be <= 2^(num bits of UINT)-1

s:=sizeof(UINT)*8

if
n <= 2^( (s-1)/r)   =>

log n <= (s-1)/r    (the base of log is 2)  =>

r*log n <= s-1      =>

log n^r <= log 2^(s-1)   =>

n^r <= 2^(s-1) <  (2^s)-1

*/

  const UINT one = 1;
  if (n > (one << ((sizeof(UINT) * 8 - 1) / r))) return true;

  UINT product1 = n;  // product1 = n*(n-1)* ...(n-r+1)
  for (UINT i = 1; i < r; i++) product1 *= --n;

  UINT rFact = r;

  while (r > 1) rFact *= --r;

  return (product1 / rFact > maxSubsetNumber);
  // product1/rFact is number of combination
  //  ( n )
  //  ( r )
}

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

void eraseSmallAngles(std::vector<std::pair<int, double>> &corners,
                      double angle) {
  std::vector<std::pair<int, double>>::iterator it = corners.begin();

  while (it != corners.end()) {
    if ((*it).second < angle)
      it = corners.erase(it);
    else
      ++it;
  }
}

//--------------------------------------------------------------------------------------
// output:
// min is the minimum angle greater or equal to minDegree (i.e the minimum angle
// of the corners)
// max is tha maximum angle greater or equal to minDegree

void detectCorners(const TStroke *stroke, double minDegree,
                   std::vector<std::pair<int, double>> &corners, double &min,
                   double &max) {
  const double minSin = fabs(sin(minDegree * M_PI_180));
  double angle, vectorialProduct, metaCornerLen, partialLen;

  UINT quadCount1 = stroke->getChunkCount();
  TPointD speed1, speed2;
  int j;

  TPointD tan1, tan2;
  min = 180;
  max = minDegree;
  for (j = 1; j < (int)quadCount1; j++) {
    speed1 = stroke->getChunk(j - 1)->getSpeed(1);
    speed2 = stroke->getChunk(j)->getSpeed(0);
    if (!(speed1 == TPointD() || speed2 == TPointD())) {
      tan1 = normalize(speed1);
      tan2 = normalize(speed2);

      vectorialProduct = fabs(cross(tan1, tan2));

      if (tan1 * tan2 < 0) {
        angle = 180 - asin(tcrop(vectorialProduct, -1.0, 1.0)) * M_180_PI;
        corners.push_back(std::make_pair(j, angle));

        //------------------------------------------
        //        TDebugMessage::getStream()<<j<<" real angle";
        //        TDebugMessage::flush();
        //------------------------------------------

        if (min > angle) min = angle;
        if (max < angle) max = angle;
      } else if (vectorialProduct >= minSin) {
        angle = asin(tcrop(vectorialProduct, -1.0, 1.0)) * M_180_PI;
        corners.push_back(std::make_pair(j, angle));

        //------------------------------------------
        //        TDebugMessage::getStream()<<j<<" real angle";
        //        TDebugMessage::flush();
        //------------------------------------------

        if (min > angle) min = angle;
        if (max < angle) max = angle;
      }
    }
  }

  const double ratioLen   = 2.5;
  const double ratioAngle = 0.2;
  std::vector<std::pair<int, double>>::iterator it = corners.begin();

  for (j = 1; j < (int)quadCount1;
       j++)  //"meta angoli" ( meta perche' derivabili)
  {
    if (it != corners.end() && j == (*it).first) {
      ++it;
      continue;
    }

    if (j - 2 >= 0 && (corners.empty() || it == corners.begin() ||
                       j - 1 != (*(it - 1)).first) &&
        j + 1 < (int)quadCount1 &&
        (corners.empty() || it == corners.end() || j + 1 != (*it).first)) {
      speed1 = stroke->getChunk(j - 2)->getSpeed(1);
      speed2 = stroke->getChunk(j + 1)->getSpeed(0);
      if (!(speed1 == TPointD() || speed2 == TPointD())) {
        tan1             = normalize(speed1);
        tan2             = normalize(speed2);
        vectorialProduct = fabs(cross(tan1, tan2));

        if (tan1 * tan2 < 0) {
          angle = 180 - asin(tcrop(vectorialProduct, -1.0, 1.0)) * M_180_PI;

          metaCornerLen = ratioLen * (stroke->getChunk(j - 1)->getLength() +
                                      stroke->getChunk(j)->getLength());
          partialLen     = 0;
          bool goodAngle = false;

          for (int i = j - 3;
               i >= 0 && (corners.empty() || it == corners.begin() ||
                          i + 1 != (*(it - 1)).first);
               i--) {
            tan1 = stroke->getChunk(i)->getSpeed(1);
            if (tan1 == TPointD()) continue;
            tan1 = normalize(tan1);

            tan2 = stroke->getChunk(i + 1)->getSpeed(1);
            if (tan2 == TPointD()) continue;
            tan2 = normalize(tan2);

            vectorialProduct = fabs(cross(tan1, tan2));
            double nearAngle =
                asin(tcrop(vectorialProduct, -1.0, 1.0)) * M_180_PI;
            if (tan1 * tan2 < 0) nearAngle = 180 - nearAngle;

            if (nearAngle > ratioAngle * angle) break;

            partialLen += stroke->getChunk(i)->getLength();
            if (partialLen > metaCornerLen) {
              goodAngle = true;
              break;
            }
          }
          if (goodAngle) {
            partialLen = 0;
            for (int i = j + 2; i + 1 < (int)quadCount1 &&
                                (corners.empty() || it == corners.end() ||
                                 i + 1 != (*it).first);
                 i++) {
              tan1 = stroke->getChunk(i)->getSpeed(0);
              if (tan1 == TPointD()) continue;
              tan1 = normalize(tan1);

              tan2 = stroke->getChunk(i + 1)->getSpeed(0);
              if (tan2 == TPointD()) continue;
              tan2 = normalize(tan2);

              vectorialProduct = fabs(cross(tan1, tan2));
              double nearAngle =
                  asin(tcrop(vectorialProduct, -1.0, 1.0)) * M_180_PI;
              if (tan1 * tan2 < 0) nearAngle = 180 - nearAngle;

              if (nearAngle > 0.1 * angle) break;

              partialLen += stroke->getChunk(i)->getLength();
              if (partialLen > metaCornerLen) {
                goodAngle = true;
                break;
              }
            }
          }

          if (goodAngle) {
            // l'angolo viene un po' declassato in quanto meta
            it = corners.insert(it, std::make_pair(j, angle * 0.7)) + 1;
            if (min > angle) min = angle;
            if (max < angle) max = angle;

            //            TDebugMessage::getStream()<<j<<" meta angle";
            //            TDebugMessage::flush();
          }
        }
      }
    }
  }
}

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

double variance(std::vector<double> &values) {
  UINT size = values.size();
  if (size == 0) return std::numeric_limits<double>::signaling_NaN();

  double sum = 0;
  UINT j     = 0;
  for (; j < size; j++) {
    sum += values[j];
  }

  double average = sum / size;

  double variance = 0;
  for (j = 0; j < size; j++) {
    variance += (average - values[j]) * (average - values[j]);
  }

  return variance / size;
}

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

void findBestSolution(const TStroke *stroke1, const TStroke *stroke2,
                      std::pair<int, double> *partialAngles1,
                      UINT partialAngles1Size,
                      const std::vector<std::pair<int, double>> &angles2,
                      UINT r,
                      std::list<std::pair<int, double>> &partialSolution,
                      double &bestValue, std::vector<int> &bestVector) {
  //-------------------------------------------------------------------
  if (r == partialAngles1Size) {
    UINT j;
    std::vector<std::pair<int, double>> angles1;

    std::list<std::pair<int, double>>::iterator it = partialSolution.begin();

    for (; it != partialSolution.end(); ++it) {
      angles1.push_back(*it);
    }
    for (j = 0; j < partialAngles1Size; j++) {
      angles1.push_back(partialAngles1[j]);
    }

    UINT angles1Size = angles1.size();
    std::vector<double> rationAngles(angles1Size), ratioLength(angles1Size + 1);
    std::vector<double> ratioX, ratioY;

    for (j = 0; j < angles1Size; j++) {
      rationAngles[j] = fabs(angles1[j].second - angles2[j].second) /
                        (angles1[j].second + angles2[j].second);
    }

    UINT firstQuad1 = 0;
    UINT firstQuad2 = 0;
    UINT nextQuad1, nextQuad2;

    TRectD bbox1 = stroke1->getBBox();
    TRectD bbox2 = stroke2->getBBox();

    double app, div;
    double invTotalLen1 = stroke1->getLength();
    assert(invTotalLen1 > 0);
    invTotalLen1 = 1.0 / invTotalLen1;

    double invTotalLen2 = stroke2->getLength();
    assert(invTotalLen2 > 0);
    invTotalLen2 = 1.0 / invTotalLen2;

    for (j = 0; j <= angles1Size; j++) {
      if (j < angles1Size) {
        nextQuad1 = angles1[j].first;
        nextQuad2 = angles2[j].first;
      } else {
        nextQuad1 = stroke1->getChunkCount();
        nextQuad2 = stroke2->getChunkCount();
      }

      ratioLength[j] =
          fabs(stroke1->getLengthAtControlPoint(nextQuad1 * 2) * invTotalLen1 -
               stroke2->getLengthAtControlPoint(nextQuad2 * 2) * invTotalLen2);

      TPointD p1(stroke1->getChunk(nextQuad1 - 1)->getP2() -
                 stroke1->getChunk(firstQuad1)->getP0());
      p1.x = fabs(p1.x);
      p1.y = fabs(p1.y);

      TPointD p2(stroke2->getChunk(nextQuad2 - 1)->getP2() -
                 stroke2->getChunk(firstQuad2)->getP0());
      p2.x = fabs(p2.x);
      p2.y = fabs(p2.y);

      app = fabs(bbox1.getLx() * p2.x - bbox2.getLx() * p1.x);
      div = (bbox1.getLx() * p2.x + bbox2.getLx() * p1.x);
      if (div) ratioX.push_back(app / div);

      app = fabs(bbox1.getLy() * p2.y - bbox2.getLy() * p1.y);
      div = (bbox1.getLy() * p2.y + bbox2.getLy() * p1.y);
      if (div) ratioY.push_back(app / div);

      firstQuad1 = nextQuad1;
      firstQuad2 = nextQuad2;
    }

    double varAng, varX, varY, varLen;

    varX   = average(ratioX);
    varY   = average(ratioY);
    varLen = average(ratioLength);
    varAng = average(rationAngles);

    double estimate = varX + varY + varAng + varLen;
    /*
#ifdef _DEBUG
for(UINT dI=0; dI<angles1Size; dI++)
{
TDebugMessage::getStream() << angles1[dI].first<<"  ";
}
TDebugMessage::flush();

TDebugMessage::getStream() <<"estimate "<< estimate<<"=" ;
TDebugMessage::flush();
TDebugMessage::getStream()<<varAng <<"+" ;
TDebugMessage::flush();
TDebugMessage::getStream()<<varX <<"+" ;
TDebugMessage::flush();
TDebugMessage::getStream()<<varY<<"+";
TDebugMessage::flush();
TDebugMessage::getStream()<<varLen;
TDebugMessage::flush();

#endif
*/

    if (estimate < bestValue) {
      bestValue = estimate;
      if (bestVector.size() != angles1Size) {
        assert(!"bad size for bestVector");
        bestVector.resize(angles1Size);
      }

      for (j = 0; j < angles1Size; j++) {
        bestVector[j] = angles1[j].first;
      }
    }

    return;
  }

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

  if (r == 1) {
    for (UINT i = 0; i < partialAngles1Size; i++) {
      findBestSolution(stroke1, stroke2, partialAngles1 + i, 1, angles2, 1,
                       partialSolution, bestValue, bestVector);
    }
    return;
  }

  partialSolution.push_back(partialAngles1[0]);
  findBestSolution(stroke1, stroke2, partialAngles1 + 1, partialAngles1Size - 1,
                   angles2, r - 1, partialSolution, bestValue, bestVector);

  partialSolution.pop_back();
  findBestSolution(stroke1, stroke2, partialAngles1 + 1, partialAngles1Size - 1,
                   angles2, r, partialSolution, bestValue, bestVector);
}

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

void findBestSolution(const TStroke *stroke1, const TStroke *stroke2,
                      std::vector<std::pair<int, double>> &angles1,
                      const std::vector<std::pair<int, double>> &angles2,
                      double &bestValue, std::vector<int> &bestVector) {
  assert(angles1.size() > angles2.size());

  std::list<std::pair<int, double>> partialSolution;

  findBestSolution(stroke1, stroke2, &(angles1[0]), angles1.size(), angles2,
                   angles2.size(), partialSolution, bestValue, bestVector);
}

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

void trivialSolution(const TStroke *stroke1, const TStroke *stroke2,
                     const std::vector<std::pair<int, double>> &angles1,
                     const std::vector<std::pair<int, double>> &angles2,
                     std::vector<int> &solution) {
  assert(angles1.size() > angles2.size());

  UINT j;
  double subStrokeRatio2;

  double invTotalLen2 = stroke2->getLength();
  assert(invTotalLen2);
  invTotalLen2 = 1.0 / invTotalLen2;

  double invTotalLen1 = stroke1->getLength();
  assert(invTotalLen1 > 0);
  invTotalLen1 = 1.0 / invTotalLen1;

  if (solution.size() != angles2.size()) {
    assert(!"bad size for solution");
    solution.resize(angles2.size());
  }

  int toBeErased = angles1.size() - angles2.size();
  UINT count     = 0;

  double diff, ratio, oldRatio = 100;

  subStrokeRatio2 =
      stroke2->getLengthAtControlPoint(angles2[count].first * 2) * invTotalLen2;

  for (j = 0; j < angles1.size() && count < solution.size(); j++) {
    if (toBeErased == 0) {
      solution[count++] = angles1[j].first;
    } else {
      ratio =
          stroke1->getLengthAtControlPoint(angles1[j].first * 2) * invTotalLen1;
      assert(ratio > 0 && ratio <= 1);

      diff = ratio - subStrokeRatio2;
      if (diff >= 0) {
        if (fabs(diff) < fabs(oldRatio - subStrokeRatio2)) {
          solution[count] = angles1[j].first;
          oldRatio        = 100;
        } else {
          assert(j > 0);
          solution[count] = angles1[j - 1].first;
        }
        count++;
        if (angles2.size() < count)
          subStrokeRatio2 =
              stroke2->getLengthAtControlPoint(angles2[count].first * 2) *
              invTotalLen2;
        else
          subStrokeRatio2 = 1;
      } else {
        toBeErased--;
        oldRatio = ratio;
      }
    }
  }
  assert(count == solution.size());
}

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

TStroke *extract(const TStroke *source, UINT firstQuad, UINT lastQuad) {
  UINT quadCount = source->getChunkCount();
  if (firstQuad >= quadCount) {
    assert(!"bad quadric index");
    firstQuad = quadCount - 1;
  }
  if (lastQuad < firstQuad) {
    assert(!"bad quadric index");
    lastQuad = firstQuad;
  }
  if (lastQuad >= quadCount) {
    assert(!"bad quadric index");
    lastQuad = quadCount - 1;
  }

  UINT cpIndex0 = firstQuad * 2;
  UINT cpIndex1 = lastQuad * 2 + 2;

  std::vector<TThickPoint> points(cpIndex1 - cpIndex0 + 1);
  UINT count = 0;
  for (UINT j = cpIndex0; j <= cpIndex1; j++) {
    points[count++] = source->getControlPoint(j);
  }

  return new TStroke(points);
}

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

void sample(const TStroke *stroke, int samplingSize,
            std::vector<TPointD> &sampledPoint) {
  double samplingFrequency = 1.0 / (double)samplingSize;
  sampledPoint.resize(samplingSize);

  double totalLen = stroke->getLength();
  double step     = totalLen * samplingFrequency;
  double len      = 0;

  for (int p = 0; p < samplingSize - 1; p++) {
    sampledPoint[p] = stroke->getPointAtLength(len);
    len += step;
  }
  sampledPoint.back() =
      stroke->getControlPoint(stroke->getControlPointCount() - 1);
}

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

class TInbetween::Imp {
public:
  //----------------------

  struct StrokeTransform {
    typedef enum { EQUAL, POINT, GENERAL } TransformationType;

    TPointD m_translate;
    TPointD m_rotationAndScaleCenter;
    double m_scaleX, m_scaleY, m_rotation;

    TransformationType m_type;

    // saved for optimization
    TAffine m_inverse;

    std::vector<int> m_firstStrokeCornerIndexes;
    std::vector<int> m_secondStrokeCornerIndexes;
  };

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

  TVectorImageP m_firstImage, m_lastImage;

  std::vector<StrokeTransform> m_transformation;

  void computeTransformation();

  void transferColor(const TVectorImageP &destination) const;

  TVectorImageP tween(double t) const;

  Imp(const TVectorImageP firstImage, const TVectorImageP lastImage)
      : m_firstImage(firstImage), m_lastImage(lastImage) {
    computeTransformation();
  }
};
//-------------------------------------------------------------------

TInbetween::TInbetween(const TVectorImageP firstImage,
                       const TVectorImageP lastImage)
    : m_imp(new TInbetween::Imp(firstImage, lastImage)) {}

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

TInbetween::~TInbetween() {}

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

void TInbetween::Imp::computeTransformation() {
  const UINT samplingPointNumber = 10;
  const UINT bboxSamplingWeight  = samplingPointNumber / 2;

  StrokeTransform transform;
  double cs, sn, totalLen1, totalLen2, constK, constQ, constB, constD, constA;
  UINT cpCount1, cpCount2;
  TPointD stroke1Centroid, stroke2Centroid, stroke1Begin, stroke2Begin,
      stroke1End, stroke2End, versor1, versor2;
  std::vector<TPointD> samplingPoint1(samplingPointNumber),
      samplingPoint2(samplingPointNumber);
  TStroke *stroke1, *stroke2;
  std::vector<double> ratioSampling, weigths, subStrokeXScaling,
      subStrokeYScaling;

  UINT strokeCount1 = m_firstImage->getStrokeCount();
  UINT strokeCount2 = m_lastImage->getStrokeCount();
  if (strokeCount1 > strokeCount2) strokeCount1 = strokeCount2;

  m_transformation.clear();
  m_transformation.reserve(strokeCount1);

  const int maxSubSetNum = (strokeCount1) ? 1000 / strokeCount1 : 1;

  UINT j, p;

  for (UINT i = 0; i < strokeCount1; i++) {
    stroke1 = m_firstImage->getStroke(i);
    stroke2 = m_lastImage->getStroke(i);

    // check if the strokes are equal
    cpCount1 = stroke1->getControlPointCount();
    cpCount2 = stroke2->getControlPointCount();

    transform.m_firstStrokeCornerIndexes.clear();
    transform.m_secondStrokeCornerIndexes.clear();
    transform.m_translate              = TPointD();
    transform.m_rotationAndScaleCenter = TPointD();
    transform.m_scaleX                 = 0;
    transform.m_scaleY                 = 0;
    transform.m_rotation               = 0;

    bool isEqual = true;

    if (cpCount1 == cpCount2) {
      for (j = 0; j < cpCount1 && isEqual; j++) {
        isEqual = (stroke1->getControlPoint(j) == stroke2->getControlPoint(j));
      }
    } else
      isEqual = false;

    if (isEqual) {
      transform.m_type = StrokeTransform::EQUAL;
    } else {
      totalLen1 = stroke1->getLength();
      totalLen2 = stroke2->getLength();

      if (totalLen1 == 0 || totalLen2 == 0) {
        if (totalLen1 == 0 && totalLen2 == 0) {
          transform.m_type = StrokeTransform::POINT;
        } else {
          transform.m_inverse = TAffine();
          transform.m_firstStrokeCornerIndexes.resize(2);
          transform.m_firstStrokeCornerIndexes[0] = 0;
          transform.m_firstStrokeCornerIndexes[1] = stroke1->getChunkCount();
          transform.m_secondStrokeCornerIndexes.resize(2);
          transform.m_secondStrokeCornerIndexes[0] = 0;
          transform.m_secondStrokeCornerIndexes[1] = stroke2->getChunkCount();
        }
      } else {
        const double startMinAngle = 30.0;
        std::vector<std::pair<int, double>> angles1, angles2;

        transform.m_type = StrokeTransform::GENERAL;

        double minAngle, maxAngle;
        int minAngle1, maxAngle1, minAngle2, maxAngle2;

        angles1.clear();
        angles2.clear();

        //        TDebugMessage::getStream()<<j<<" stroke1 corner detection";
        //        TDebugMessage::flush();

        detectCorners(stroke1, startMinAngle, angles1, minAngle, maxAngle);
        minAngle1 = (int)minAngle;
        maxAngle1 = (int)maxAngle;

        //        TDebugMessage::getStream()<<j<<" stroke2 corner detection";
        //        TDebugMessage::flush();

        detectCorners(stroke2, startMinAngle, angles2, minAngle, maxAngle);
        minAngle2 = (int)minAngle;
        maxAngle2 = (int)maxAngle;

        if (angles1.empty()) angles2.clear();
        if (angles2.empty()) angles1.clear();

        /*
debugStream.open("c:\\temp\\inbetween.txt", ios_base::out);
debugStream <<"num angoli 1: "<< angles1.size() << endl;
debugStream <<"num angoli 2: "<< angles2.size() << endl;
*/

        double bestValue = (std::numeric_limits<double>::max)();

        if (angles1.size() != angles2.size()) {
          bestValue = (std::numeric_limits<double>::max)();

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

          if (isTooComplex(angles1.size(), angles2.size(), maxSubSetNum)) {
            // debugStream <<"is too complex" << endl;
            int firstAngle =
                (int)((angles1.size() < angles2.size()) ? minAngle2
                                                        : minAngle1);
            int lastAngle =
                (int)((angles1.size() < angles2.size()) ? maxAngle1
                                                        : maxAngle2);
            int bestAngle = (int)startMinAngle;

            if ((int)(angles1.size() + angles2.size()) <
                lastAngle - firstAngle + 1) {
              double tempAngle;
              std::vector<double> sortedAngles1, sortedAngles2;
              sortedAngles1.reserve(angles1.size());
              sortedAngles2.reserve(angles2.size());
              for (j = 0; j < angles1.size(); j++) {
                tempAngle = angles1[j].second;
                if (tempAngle >= firstAngle && tempAngle <= lastAngle)
                  sortedAngles1.push_back(tempAngle);
              }
              for (j = 0; j < angles2.size(); j++) {
                tempAngle = angles2[j].second;
                if (tempAngle >= firstAngle && tempAngle <= lastAngle)
                  sortedAngles2.push_back(tempAngle);
              }
              std::vector<double> sortedAngles(sortedAngles1.size() +
                                               sortedAngles2.size());

              std::sort(sortedAngles1.begin(), sortedAngles1.end());
              std::sort(sortedAngles2.begin(), sortedAngles2.end());
              std::merge(sortedAngles1.begin(), sortedAngles1.end(),
                         sortedAngles2.begin(), sortedAngles2.end(),
                         sortedAngles.begin());

              for (j = 0; j < sortedAngles.size(); j++) {
                int numAng1 = angleNumber(angles1, sortedAngles[j]);
                int numAng2 = angleNumber(angles2, sortedAngles[j]);
                double val  = (numAng1 == numAng2)
                                 ? 0
                                 : fabs((float)(numAng1 - numAng2)) /
                                       (numAng1 + numAng2);
                if (val < bestValue) {
                  bestValue = val;
                  bestAngle = (int)(sortedAngles[j]);
                  if (bestValue == 0 ||
                      !isTooComplex(numAng1, numAng2, maxSubSetNum))
                    break;
                }
              }

            } else  //-----------------------------------------------------
            {
              for (int angle = firstAngle; angle <= lastAngle; angle++) {
                int numAng1 = angleNumber(angles1, angle);
                int numAng2 = angleNumber(angles2, angle);
                double val  = (numAng1 == numAng2)
                                 ? 0
                                 : fabs((float)(numAng1 - numAng2)) /
                                       (numAng1 + numAng2);
                if (val < bestValue) {
                  bestValue = val;
                  bestAngle = angle;
                  if (bestValue == 0 ||
                      !isTooComplex(numAng1, numAng2, maxSubSetNum))
                    break;
                }
              }
            }

            eraseSmallAngles(angles1, bestAngle);
            eraseSmallAngles(angles2, bestAngle);

            /*
debugStream <<"bestAngle: "<< bestAngle << endl;
debugStream <<"num angoli 1: "<< angles1.size() << endl;
debugStream <<"num angoli 2: "<< angles2.size() << endl;
*/
          }
          //--------------------------------------------------------------------------

          bestValue = (std::numeric_limits<double>::max)();

          if (angles1.size() == angles2.size()) {
            transform.m_firstStrokeCornerIndexes.push_back(0);
            for (j = 0; j < angles1.size(); j++)
              transform.m_firstStrokeCornerIndexes.push_back(angles1[j].first);
            transform.m_firstStrokeCornerIndexes.push_back(
                stroke1->getChunkCount());

            transform.m_secondStrokeCornerIndexes.push_back(0);
            for (j = 0; j < angles2.size(); j++)
              transform.m_secondStrokeCornerIndexes.push_back(angles2[j].first);
            transform.m_secondStrokeCornerIndexes.push_back(
                stroke2->getChunkCount());
          } else {
            if (isTooComplex(angles1.size(), angles2.size(), maxSubSetNum)) {
              if (angles1.size() > angles2.size()) {
                transform.m_firstStrokeCornerIndexes.resize(angles2.size());
                trivialSolution(stroke1, stroke2, angles1, angles2,
                                transform.m_firstStrokeCornerIndexes);
                transform.m_firstStrokeCornerIndexes.insert(
                    transform.m_firstStrokeCornerIndexes.begin(), 0);
                transform.m_firstStrokeCornerIndexes.push_back(
                    stroke1->getChunkCount());

                transform.m_secondStrokeCornerIndexes.push_back(0);
                for (j = 0; j < angles2.size(); j++)
                  transform.m_secondStrokeCornerIndexes.push_back(
                      angles2[j].first);
                transform.m_secondStrokeCornerIndexes.push_back(
                    stroke2->getChunkCount());
              } else {
                transform.m_firstStrokeCornerIndexes.push_back(0);
                for (j = 0; j < angles1.size(); j++)
                  transform.m_firstStrokeCornerIndexes.push_back(
                      angles1[j].first);
                transform.m_firstStrokeCornerIndexes.push_back(
                    stroke1->getChunkCount());

                transform.m_secondStrokeCornerIndexes.resize(angles1.size());
                trivialSolution(stroke2, stroke1, angles2, angles1,
                                transform.m_secondStrokeCornerIndexes);
                transform.m_secondStrokeCornerIndexes.insert(
                    transform.m_secondStrokeCornerIndexes.begin(), 0);
                transform.m_secondStrokeCornerIndexes.push_back(
                    stroke2->getChunkCount());
              }
            } else {
              if (angles1.size() > angles2.size()) {
                transform.m_firstStrokeCornerIndexes.resize(angles2.size());
                findBestSolution(stroke1, stroke2, angles1, angles2, bestValue,
                                 transform.m_firstStrokeCornerIndexes);
                transform.m_firstStrokeCornerIndexes.insert(
                    transform.m_firstStrokeCornerIndexes.begin(), 0);
                transform.m_firstStrokeCornerIndexes.push_back(
                    stroke1->getChunkCount());

                transform.m_secondStrokeCornerIndexes.push_back(0);
                for (j = 0; j < angles2.size(); j++)
                  transform.m_secondStrokeCornerIndexes.push_back(
                      angles2[j].first);
                transform.m_secondStrokeCornerIndexes.push_back(
                    stroke2->getChunkCount());
              } else {
                transform.m_firstStrokeCornerIndexes.push_back(0);
                for (j = 0; j < angles1.size(); j++)
                  transform.m_firstStrokeCornerIndexes.push_back(
                      angles1[j].first);
                transform.m_firstStrokeCornerIndexes.push_back(
                    stroke1->getChunkCount());

                transform.m_secondStrokeCornerIndexes.resize(angles1.size());
                findBestSolution(stroke2, stroke1, angles2, angles1, bestValue,
                                 transform.m_secondStrokeCornerIndexes);
                transform.m_secondStrokeCornerIndexes.insert(
                    transform.m_secondStrokeCornerIndexes.begin(), 0);
                transform.m_secondStrokeCornerIndexes.push_back(
                    stroke2->getChunkCount());
              }
            }
          }
        } else {
          transform.m_firstStrokeCornerIndexes.push_back(0);
          for (j = 0; j < angles1.size(); j++)
            transform.m_firstStrokeCornerIndexes.push_back(angles1[j].first);
          transform.m_firstStrokeCornerIndexes.push_back(
              stroke1->getChunkCount());

          transform.m_secondStrokeCornerIndexes.push_back(0);
          for (j = 0; j < angles2.size(); j++)
            transform.m_secondStrokeCornerIndexes.push_back(angles2[j].first);
          transform.m_secondStrokeCornerIndexes.push_back(
              stroke2->getChunkCount());
        }

        UINT cornerSize = transform.m_firstStrokeCornerIndexes.size();
        assert(cornerSize == transform.m_secondStrokeCornerIndexes.size());
        assert(cornerSize >= 2);

        double totalRadRotation = 0;

        TStroke *subStroke1 = 0;
        TStroke *subStroke2 = 0;

        stroke1Centroid = stroke1->getCentroid();
        stroke2Centroid = stroke2->getCentroid();

#ifdef _DEBUG
        assert(transform.m_firstStrokeCornerIndexes[0] == 0);
        assert(transform.m_secondStrokeCornerIndexes[0] == 0);
        assert(transform.m_firstStrokeCornerIndexes[cornerSize - 1] ==
               stroke1->getChunkCount());
        assert(transform.m_secondStrokeCornerIndexes[cornerSize - 1] ==
               stroke2->getChunkCount());
        for (j = 0; j < cornerSize - 1; j++) {
          assert(transform.m_firstStrokeCornerIndexes[j] <
                 transform.m_firstStrokeCornerIndexes[j + 1]);
          assert(transform.m_secondStrokeCornerIndexes[j] <
                 transform.m_secondStrokeCornerIndexes[j + 1]);

          assert(transform.m_firstStrokeCornerIndexes[j] <
                 stroke1->getChunkCount());
          assert(transform.m_secondStrokeCornerIndexes[j] <
                 stroke2->getChunkCount());
        }
#endif

        for (j = 0; j < cornerSize - 1; j++) {
          ///////////////////////////////////////// sampling

          subStroke1 = extract(stroke1, transform.m_firstStrokeCornerIndexes[j],
                               transform.m_firstStrokeCornerIndexes[j + 1] - 1);
          sample(subStroke1, samplingPointNumber, samplingPoint1);
          subStroke2 =
              extract(stroke2, transform.m_secondStrokeCornerIndexes[j],
                      transform.m_secondStrokeCornerIndexes[j + 1] - 1);
          sample(subStroke2, samplingPointNumber, samplingPoint2);

          ///////////////////////////////////////// compute Rotation

          ratioSampling.clear();
          ratioSampling.reserve(samplingPointNumber);
          weigths.clear();
          weigths.reserve(samplingPointNumber);

          TPointD pOld, pNew;
          // double totalW=0;

          for (p = 0; p < samplingPointNumber; p++) {
            pOld = samplingPoint1[p];
            pNew = samplingPoint2[p];
            if (pOld == stroke1Centroid) continue;
            if (pNew == stroke2Centroid) continue;
            versor1 = normalize(pOld - stroke1Centroid);
            versor2 = normalize(pNew - stroke2Centroid);
            weigths.push_back(tdistance(pOld, stroke1Centroid) +
                              tdistance(pNew, stroke2Centroid));
            cs       = versor1 * versor2;
            sn       = cross(versor1, versor2);
            double v = atan2(sn, cs);
            ratioSampling.push_back(v);
          }
          delete subStroke1;
          delete subStroke2;
          subStroke1 = 0;
          subStroke2 = 0;

          double radRotation = weightedAverage(ratioSampling, weigths);

          totalRadRotation += radRotation;
        }
        totalRadRotation /= (cornerSize - 1);
        transform.m_rotation = totalRadRotation * M_180_PI;

        if (isAlmostZero(transform.m_rotation, 2)) {
          transform.m_rotation = 0.0;
          totalRadRotation     = 0.0;
        }

#ifdef _DEBUG
//        TDebugMessage::getStream()<<"rotation "<< transform.m_rotation;
//        TDebugMessage::flush();
#endif

        ///////////////////////////////////////// compute Scale

        if (transform.m_rotation == 0.0) {
          subStrokeXScaling.clear();
          subStrokeXScaling.reserve(cornerSize - 1);
          subStrokeYScaling.clear();
          subStrokeYScaling.reserve(cornerSize - 1);

          for (j = 0; j < cornerSize - 1; j++) {
            ///////////////////////////////////////// sampling

            subStroke1 =
                extract(stroke1, transform.m_firstStrokeCornerIndexes[j],
                        transform.m_firstStrokeCornerIndexes[j + 1] - 1);
            sample(subStroke1, samplingPointNumber, samplingPoint1);
            subStroke2 =
                extract(stroke2, transform.m_secondStrokeCornerIndexes[j],
                        transform.m_secondStrokeCornerIndexes[j + 1] - 1);
            sample(subStroke2, samplingPointNumber, samplingPoint2);

            ///////////////////////////////////////// compute X Scale

            ratioSampling.clear();
            ratioSampling.reserve(samplingPointNumber + bboxSamplingWeight);
            double appX, appY;

            TPointD appPoint;
            double bboxXMin, bboxYMin, bboxXMax, bboxYMax;
            bboxXMin = bboxYMin = (std::numeric_limits<double>::max)();
            bboxXMax = bboxYMax = -(std::numeric_limits<double>::max)();
            int h;

            for (h = 0; h < subStroke1->getControlPointCount(); ++h) {
              appPoint = subStroke1->getControlPoint(h);
              if (appPoint.x < bboxXMin) bboxXMin = appPoint.x;
              if (appPoint.x > bboxXMax) bboxXMax = appPoint.x;
              if (appPoint.y < bboxYMin) bboxYMin = appPoint.y;
              if (appPoint.y > bboxYMax) bboxYMax = appPoint.y;
            }

            appX = bboxXMax - bboxXMin;
            appY = bboxYMax - bboxYMin;

            if (appX) {
              bboxXMin = (std::numeric_limits<double>::max)();
              bboxXMax = -(std::numeric_limits<double>::max)();
              for (h = 0; h < subStroke2->getControlPointCount(); ++h) {
                appPoint = subStroke2->getControlPoint(h);
                if (appPoint.x < bboxXMin) bboxXMin = appPoint.x;
                if (appPoint.x > bboxXMax) bboxXMax = appPoint.x;
              }
              appX = (isAlmostZero(appX, 1e-01)) ? -1
                                                 : (bboxXMax - bboxXMin) / appX;
              for (UINT tms = 0; tms < bboxSamplingWeight && appX >= 0; tms++)
                ratioSampling.push_back(appX);

              for (p = 0; p < samplingPointNumber; p++) {
                appX = fabs(samplingPoint1[p].x - stroke1Centroid.x);
                if (appX)
                  ratioSampling.push_back(
                      fabs(samplingPoint2[p].x - stroke2Centroid.x) / appX);
              }

              if (!ratioSampling.empty()) {
                subStrokeXScaling.push_back(average(ratioSampling));
              }
            }

            ///////////////////////////////////////// compute Y Scale

            ratioSampling.clear();
            ratioSampling.reserve(samplingPointNumber + bboxSamplingWeight);

            if (appY) {
              bboxYMin = (std::numeric_limits<double>::max)();
              bboxYMax = -(std::numeric_limits<double>::max)();
              for (h = 0; h < subStroke2->getControlPointCount(); ++h) {
                appPoint = subStroke2->getControlPoint(h);
                if (appPoint.y < bboxYMin) bboxYMin = appPoint.y;
                if (appPoint.y > bboxYMax) bboxYMax = appPoint.y;
              }

              appY = (isAlmostZero(appY, 1e-01)) ? -1
                                                 : (bboxYMax - bboxYMin) / appY;
              for (UINT tms = 0; tms < bboxSamplingWeight && appY >= 0; tms++)
                ratioSampling.push_back(appY);

              for (p = 0; p < samplingPointNumber; p++) {
                appY = fabs(samplingPoint1[p].y - stroke1Centroid.y);
                if (appY)
                  ratioSampling.push_back(
                      fabs(samplingPoint2[p].y - stroke2Centroid.y) / appY);
              }

              if (!ratioSampling.empty()) {
                subStrokeYScaling.push_back(average(ratioSampling));
              }
            }

            delete subStroke1;
            delete subStroke2;
            subStroke1 = 0;
            subStroke2 = 0;
          }

          if (subStrokeXScaling.empty()) {
            transform.m_scaleX = 1.0;
          } else {
            transform.m_scaleX = average(subStrokeXScaling);
            if (isAlmostZero(transform.m_scaleX - 1.0))
              transform.m_scaleX = 1.0;
          }

          if (subStrokeYScaling.empty()) {
            transform.m_scaleY = 1.0;
          } else {
            transform.m_scaleY = average(subStrokeYScaling);
            if (isAlmostZero(transform.m_scaleY - 1.0))
              transform.m_scaleY = 1.0;
          }
          /*
#ifdef _DEBUG

TDebugMessage::getStream()<<"x scale "<< transform.m_scaleX ;
TDebugMessage::flush();
TDebugMessage::getStream()<<"y scale "<< transform.m_scaleY ;
TDebugMessage::flush();
#endif
*/
        } else {
          subStrokeXScaling.clear();
          subStrokeXScaling.reserve(cornerSize - 1);

          for (j = 0; j < cornerSize - 1; j++) {
            ///////////////////////////////////////// sampling

            subStroke1 =
                extract(stroke1, transform.m_firstStrokeCornerIndexes[j],
                        transform.m_firstStrokeCornerIndexes[j + 1] - 1);
            sample(subStroke1, samplingPointNumber, samplingPoint1);
            subStroke2 =
                extract(stroke2, transform.m_secondStrokeCornerIndexes[j],
                        transform.m_secondStrokeCornerIndexes[j + 1] - 1);
            sample(subStroke2, samplingPointNumber, samplingPoint2);

            ///////////////////////////////////////// compute Scale

            ratioSampling.clear();
            ratioSampling.reserve(samplingPointNumber + bboxSamplingWeight);

            TRectD bbox1 = subStroke1->getBBox();
            double app   = tdistance2(bbox1.getP00(), bbox1.getP11());
            if (app) {
              TRectD bbox2 =
                  TRotation(transform.m_rotation).inv() * subStroke2->getBBox();
              app = sqrt(tdistance2(bbox2.getP00(), bbox2.getP11()) / app);

              for (UINT tms = 0; tms < bboxSamplingWeight; tms++)
                ratioSampling.push_back(app);

              double app;
              for (p = 0; p < samplingPointNumber; p++) {
                app = tdistance(samplingPoint1[p], stroke1Centroid);
                if (app) {
                  ratioSampling.push_back(
                      tdistance(samplingPoint2[p], stroke2Centroid) / app);
                }
              }
            }
            if (!ratioSampling.empty()) {
              subStrokeXScaling.push_back(average(ratioSampling));
            }

            delete subStroke1;
            delete subStroke2;
            subStroke1 = 0;
            subStroke2 = 0;
          }

          if (subStrokeXScaling.empty()) {
            transform.m_scaleX = transform.m_scaleY = 1.0;
          } else {
            transform.m_scaleX = transform.m_scaleY =
                average(subStrokeXScaling);
            if (isAlmostZero(transform.m_scaleX - 1.0, 0.00001))
              transform.m_scaleX = transform.m_scaleY = 1.0;
          }
          /*
#ifdef _DEBUG

TDebugMessage::getStream()<<"scale "<< transform.m_scaleX ;
TDebugMessage::flush();
#endif
*/
        }

        ///////////////////////////////////////// compute centre of Rotation and
        ///Scaling

        std::vector<TPointD> vpOld(cornerSize), vpNew(cornerSize);
        TPointD pOld, pNew;

        for (j = 0; j < cornerSize - 1; j++) {
          vpOld[j] = stroke1->getChunk(transform.m_firstStrokeCornerIndexes[j])
                         ->getP0();
          vpNew[j] = stroke2->getChunk(transform.m_secondStrokeCornerIndexes[j])
                         ->getP0();
        }
        vpOld[j] = stroke1->getControlPoint(stroke1->getControlPointCount());
        vpNew[j] = stroke2->getControlPoint(stroke2->getControlPointCount());

        if (transform.m_rotation == 0.0) {
          if (transform.m_scaleX == 1.0 && transform.m_scaleY == 1.0) {
            transform.m_translate = stroke2Centroid - stroke1Centroid;
            transform.m_rotationAndScaleCenter = TPointD();
          } else {
            if (transform.m_scaleX == 1.0) {
              transform.m_rotationAndScaleCenter.x = 0;
              transform.m_translate.x              = 0;
              for (j = 0; j < cornerSize; j++) {
                transform.m_translate.x += vpNew[j].x - vpOld[j].x;
              }
              transform.m_translate.x = transform.m_translate.x / cornerSize;
            } else {
              transform.m_rotationAndScaleCenter.x = 0;

              for (j = 0; j < cornerSize; j++) {
                pOld = vpOld[j];
                pNew = vpNew[j];
                transform.m_rotationAndScaleCenter.x +=
                    (transform.m_scaleX * pOld.x - pNew.x) /
                    (transform.m_scaleX - 1.0);
              }
              transform.m_rotationAndScaleCenter.x =
                  transform.m_rotationAndScaleCenter.x / cornerSize;
              transform.m_translate.x = 0;
            }

            if (transform.m_scaleY == 1.0) {
              transform.m_rotationAndScaleCenter.y = 0;
              transform.m_translate.y              = 0;
              for (j = 0; j < cornerSize; j++) {
                transform.m_translate.y += vpNew[j].y - vpOld[j].y;
              }
              transform.m_translate.y = transform.m_translate.y / cornerSize;
            } else {
              transform.m_rotationAndScaleCenter.y = 0;

              for (j = 0; j < cornerSize; j++) {
                pOld = vpOld[j];
                pNew = vpNew[j];
                transform.m_rotationAndScaleCenter.y +=
                    (transform.m_scaleY * pOld.y - pNew.y) /
                    (transform.m_scaleY - 1.0);
              }
              transform.m_rotationAndScaleCenter.y =
                  transform.m_rotationAndScaleCenter.y / cornerSize;
              transform.m_translate.y = 0;
            }
          }

        } else {
          assert(transform.m_scaleX == transform.m_scaleY);

          cs = transform.m_scaleX * cos(totalRadRotation);
          sn = transform.m_scaleX * sin(totalRadRotation);

          // scelgo punti da usare come vincolo, per ottenere la translazione,
          // dato un centro di rotazione

          // dato il punto pOld e pNew si calcola analiticamnete il punto di
          // rotazione e scala
          // che minimizza la traslazione aggiuntiva e la traslazione stessa

          for (j = 0; j < cornerSize; j++) {
            pOld   = vpOld[j];
            pNew   = vpNew[j];
            constK = pNew.x - cs * pOld.x + sn * pOld.y;
            constQ = pNew.y - sn * pOld.x - cs * pOld.y;
            constB = 2 * (constK * (cs - 1.) + sn * constQ);
            constD = 2 * (constQ * (cs - 1.) - sn * constK);
            constA = transform.m_scaleX * transform.m_scaleX + 1 - 2 * cs;
            assert(constA > 0);
            constA = 1.0 / (2 * constA);
            transform.m_rotationAndScaleCenter.x += -constB * constA;
            transform.m_rotationAndScaleCenter.y += -constD * constA;
          }

          transform.m_rotationAndScaleCenter =
              transform.m_rotationAndScaleCenter * (1.0 / (double)cornerSize);

          transform.m_translate.x =
              (cs - 1.0) * transform.m_rotationAndScaleCenter.x -
              sn * transform.m_rotationAndScaleCenter.y + constK;

          transform.m_translate.y =
              sn * transform.m_rotationAndScaleCenter.x +
              (cs - 1.0) * transform.m_rotationAndScaleCenter.y + constQ;
        }

        /////////////////////////////////////////

        transform.m_inverse = (TTranslation(transform.m_translate) *
                               TScale(transform.m_rotationAndScaleCenter,
                                      transform.m_scaleX, transform.m_scaleY) *
                               TRotation(transform.m_rotationAndScaleCenter,
                                         transform.m_rotation))
                                  .inv();

        // debugStream.close();
      }  // end if !isPoint

    }  // end if !isEqual

    m_transformation.push_back(transform);

  }  // end for each stroke
}

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

TVectorImageP TInbetween::Imp::tween(double t) const {
  const double step             = 5.0;
  const double interpolateError = 1.0;

  TVectorImageP vi = new TVectorImage;

  UINT strokeCount1 = m_firstImage->getStrokeCount();
  UINT strokeCount2 = m_lastImage->getStrokeCount();
  vi->setPalette(m_firstImage->getPalette());
  if (strokeCount1 > strokeCount2) strokeCount1 = strokeCount2;

  assert(m_transformation.size() == strokeCount1);

  double totalLen1, totalLen2, len1, len2, step1, step2;
  std::vector<TThickPoint> points;
  TStroke *stroke1, *stroke2, *subStroke1, *subStroke2, *stroke;

  TAffine mt, invMatrix;
  TThickPoint point2, finalPoint;
  UINT i, j, cp, cpSize;

  for (i = 0; i < strokeCount1; i++) {
    stroke1 = m_firstImage->getStroke(i);
    stroke2 = m_lastImage->getStroke(i);

    if (m_transformation[i].m_type == StrokeTransform::EQUAL) {
      stroke = new TStroke(*stroke1);
    } else {
      points.clear();
      totalLen1 = stroke1->getLength();
      totalLen2 = stroke2->getLength();
      ;

      if (stroke1->getControlPointCount() == stroke2->getControlPointCount() &&
          stroke1->isSelfLoop() && stroke2->isSelfLoop()) {
        for (int i = 0; i < stroke1->getControlPointCount(); i++) {
          TThickPoint p0 = stroke1->getControlPoint(i);
          TThickPoint p1 = stroke2->getControlPoint(i);
          points.push_back(p0 * (1 - t) + p1 * t);
        }
        stroke = new TStroke(points);
      } else if (m_transformation[i].m_type == StrokeTransform::POINT) {
        TThickPoint pOld = stroke1->getThickPointAtLength(0.5 * totalLen1);
        TThickPoint pNew = stroke2->getThickPointAtLength(0.5 * totalLen2);
        points.push_back(pOld * (1 - t) + pNew * t);
        points.push_back(points[0]);
        points.push_back(points[0]);
        stroke = new TStroke(points);
      } else {
        mt = (m_transformation[i].m_inverse.isIdentity())
                 ? TAffine()
                 : TTranslation(m_transformation[i].m_translate * t) *
                       TScale(m_transformation[i].m_rotationAndScaleCenter,
                              (1 - t) + t * m_transformation[i].m_scaleX,
                              (1 - t) + t * m_transformation[i].m_scaleY) *
                       TRotation(m_transformation[i].m_rotationAndScaleCenter,
                                 m_transformation[i].m_rotation * t);

        UINT cornerSize = m_transformation[i].m_firstStrokeCornerIndexes.size();
        assert(cornerSize ==
               m_transformation[i].m_secondStrokeCornerIndexes.size());
        if (cornerSize > m_transformation[i].m_secondStrokeCornerIndexes.size())
          cornerSize = m_transformation[i].m_secondStrokeCornerIndexes.size();

        assert(cornerSize >= 2);

        std::vector<TThickPoint> controlPoints;

        // if not m_transformation[i].m_findCorners => detect corner return
        // different size =>cornerSize==2
        // assert(!m_transformation[i].m_findCorners || cornerSize==2);

        for (j = 0; j < cornerSize - 1; j++) {
          points.clear();
          subStroke1 = extract(
              stroke1, m_transformation[i].m_firstStrokeCornerIndexes[j],
              m_transformation[i].m_firstStrokeCornerIndexes[j + 1] - 1);
          subStroke2 = extract(
              stroke2, m_transformation[i].m_secondStrokeCornerIndexes[j],
              m_transformation[i].m_secondStrokeCornerIndexes[j + 1] - 1);

          totalLen1 = subStroke1->getLength();
          totalLen2 = subStroke2->getLength();

          if (totalLen1 > totalLen2) {
            step1 = step;
            step2 = (totalLen2 / totalLen1) * step;
          } else {
            step1 = (totalLen1 / totalLen2) * step;
            step2 = step;
          }

          len1 = 0;
          len2 = 0;

          while (len1 <= totalLen1 && len2 <= totalLen2) {
            point2 = subStroke2->getThickPointAtLength(len2);
            point2 = TThickPoint(m_transformation[i].m_inverse *
                                     subStroke2->getThickPointAtLength(len2),
                                 point2.thick);
            finalPoint =
                subStroke1->getThickPointAtLength(len1) * (1 - t) + t * point2;

            points.push_back(TThickPoint(mt * (finalPoint), finalPoint.thick));
            len1 += step1;
            len2 += step2;
          }
          point2 = subStroke2->getThickPointAtLength(totalLen2);
          point2 = TThickPoint(m_transformation[i].m_inverse *
                                   subStroke2->getThickPointAtLength(totalLen2),
                               point2.thick);
          finalPoint = subStroke1->getThickPointAtLength(totalLen1) * (1 - t) +
                       t * point2;

          points.push_back(TThickPoint(mt * (finalPoint), finalPoint.thick));

          stroke = TStroke::interpolate(points, interpolateError, false
                                        /*m_transformation[i].m_findCorners*/);

          if (j == 0)
            controlPoints.push_back(stroke->getControlPoint(0));
          else
            controlPoints.back() =
                (controlPoints.back() + stroke->getControlPoint(0)) * 0.5;

          cpSize = stroke->getControlPointCount();
          for (cp = 1; cp < cpSize; cp++) {
            controlPoints.push_back(stroke->getControlPoint(cp));
          }

          delete subStroke1;
          delete subStroke2;
          delete stroke;
          subStroke1 = 0;
          subStroke2 = 0;
          stroke     = 0;
        }

        stroke = new TStroke(controlPoints);
      }
    }

    if (stroke1->isSelfLoop() && stroke2->isSelfLoop()) stroke->setSelfLoop();

    stroke->setStyle(stroke1->getStyle());
    stroke->outlineOptions() = stroke1->outlineOptions();
    VIStroke *vs =
        new VIStroke(stroke, m_firstImage->getVIStroke(i)->m_groupId);
    vi->m_imp->m_strokes.push_back(vs);

  }  // end for each stroke

  if (m_firstImage->isComputedRegionAlmostOnce()) transferColor(vi);

  return vi;
}

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

void TInbetween::Imp::transferColor(const TVectorImageP &destination) const {
  const TVectorImageP &original = m_firstImage;
  destination->setPalette(original->getPalette());

  destination->findRegions();

  if (destination->getRegionCount()) {
    UINT strokeCount1 = original->getStrokeCount();
    UINT strokeCount2 = destination->getStrokeCount();
    if (strokeCount1 > strokeCount2) strokeCount1 = strokeCount2;

    for (UINT i = 0; i < strokeCount1; i++) {
      TVectorImage::transferStrokeColors(original, i, destination, i);
    }
  }
}

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

TVectorImageP TInbetween::tween(double t) const { return m_imp->tween(t); }

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