#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></utility>

#include <limits></limits>

#include <list></list>

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

//#include "tdebugmessage.h"

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

double average(std::vector<double> &values, double range = 2.5)</double>

{

	UINT size = values.size();

	if (size == 0)

		return std::numeric_limits<double>::signaling_NaN();</double>

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

{

	UINT size = values.size();

	if (size == 0)

		return std::numeric_limits<double>::signaling_NaN();</double>

	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>

	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)</std::pair<int,>

{

	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;

	register UINT product1 = n; //product1 = n*(n-1)* ...(n-r+1)

	for (register UINT i = 1; i < r; i++)

		product1 *= --n;

	register 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::pair<int,>

{

	std::vector<std::pair<int, double="">>::iterator it = corners.begin();</std::pair<int,>

	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)</std::pair<int,>

{

	const double minSin = fabs(sin(minDegree * TConsts::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)) * TConsts::invOf_pi_180;

				corners.push_back(std::make_pair(j, angle));

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

				//        TDebugMessage::getStream()<

				//        TDebugMessage::flush();

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

				if (min > angle)

					min = angle;

				if (max < angle)

					max = angle;

			} else if (vectorialProduct >= minSin) {

				angle = asin(tcrop(vectorialProduct, -1.0, 1.0)) * TConsts::invOf_pi_180;

				corners.push_back(std::make_pair(j, angle));

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

				//        TDebugMessage::getStream()<

				//        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();</std::pair<int,>

	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)) * TConsts::invOf_pi_180;

					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)) * TConsts::invOf_pi_180;

						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)) * TConsts::invOf_pi_180;

							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()<

						//            TDebugMessage::flush();

					}

				}

			}

		}

	}

}

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

double variance(std::vector<double> &values)</double>

{

	UINT size = values.size();

	if (size == 0)

		return std::numeric_limits<double>::signaling_NaN();</double>

	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,</int,>

					  const std::vector<std::pair<int, double="">> &angles2,</std::pair<int,>

					  UINT r, std::list<std::pair<int, double="">> &partialSolution,</std::pair<int,>

					  double &bestValue,

					  std::vector<int> &bestVector)</int>

{

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

	if (r == partialAngles1Size) {

		UINT j;

		vector<std::pair<int, double="">> angles1;</std::pair<int,>

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

		for (; it != partialSolution.end(); ++it) {

			angles1.push_back(*it);

		}

		for (j = 0; j < partialAngles1Size; j++) {

			angles1.push_back(partialAngles1[j]);

		}

		UINT angles1Size = angles1.size();

		vector<double> rationAngles(angles1Size), ratioLength(angles1Size + 1);</double>

		vector<double> ratioX, ratioY;</double>

		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

    {

      TDebugMessage::getStream() << angles1[dI].first<<"  ";      

    }

    TDebugMessage::flush();

    TDebugMessage::getStream() <<"estimate "<< estimate<<"=" ;

    TDebugMessage::flush();

    TDebugMessage::getStream()<

    TDebugMessage::flush();

    TDebugMessage::getStream()<

    TDebugMessage::flush();

    TDebugMessage::getStream()<

    TDebugMessage::flush();

    TDebugMessage::getStream()<

    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,</std::pair<int,>

					  const std::vector<std::pair<int, double="">> &angles2,</std::pair<int,>

					  double &bestValue,

					  std::vector<int> &bestVector)</int>

{

	assert(angles1.size() > angles2.size());

	list<std::pair<int, double="">> partialSolution;</std::pair<int,>

	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,</std::pair<int,>

					 const std::vector<std::pair<int, double="">> &angles2,</std::pair<int,>

					 std::vector<int> &solution)</int>

{

	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;

	vector<tthickpoint> points(cpIndex1 - cpIndex0 + 1);</tthickpoint>

	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, vector<tpointd> &sampledPoint)</tpointd>

{

	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;

		vector<int> m_firstStrokeCornerIndexes;</int>

		vector<int> m_secondStrokeCornerIndexes;</int>

	};

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

	TVectorImageP m_firstImage, m_lastImage;

	std::vector<stroketransform> m_transformation;</stroketransform>

	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()

{

	delete m_imp;

}

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

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;

	vector<tpointd> samplingPoint1(samplingPointNumber), samplingPoint2(samplingPointNumber);</tpointd>

	TStroke *stroke1, *stroke2;

	vector<double> ratioSampling, weigths, subStrokeXScaling, subStrokeYScaling;</double>

	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;

				vector<pair<int, double="">> angles1, angles2;</pair<int,>

				transform.m_type = StrokeTransform::GENERAL;

				double minAngle, maxAngle;

				int minAngle1, maxAngle1, minAngle2, maxAngle2;

				angles1.clear();

				angles2.clear();

				//        TDebugMessage::getStream()<

				//        TDebugMessage::flush();

				detectCorners(stroke1, startMinAngle, angles1, minAngle, maxAngle);

				minAngle1 = (int)minAngle;

				maxAngle1 = (int)maxAngle;

				//        TDebugMessage::getStream()<

				//        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)();</double>

				if (angles1.size() != angles2.size()) {

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

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

					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;

							vector<double> sortedAngles1, sortedAngles2;</double>

							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);

							}

							vector<double> sortedAngles(sortedAngles1.size() + sortedAngles2.size());</double>

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