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


#pragma hdrstop

#include "Vector_F.h"
#include "BezierCurveApproximator.h"
#include "math.h"
//---------------------------------------------------------------------------

#pragma package(smart_init)
BezierCurveApproximator::BezierCurveApproximator():Error_Threshold(1) {
}
bool BezierCurveApproximator::Approximate(const int start_index, const int end_index, BezierCurvePoints& result_curve_points) const
{
  P_ITER start, end;
  start = end = freeDrawingPointList.begin();
	std::advance(start, start_index);
  std::advance(end, end_index);

  int list_size = end_index - start_index;

  if(list_size < 2)
    return true;

	if(list_size == 3)
  {
    result_curve_points.start = *start;
    ++start;
    result_curve_points.control1 = *start;
    ++start;
    result_curve_points.control2 = *start;
    result_curve_points.end = *start;
//    result_curve_points.end = *end;
    result_curve_points.start_cross_product = 0;
    result_curve_points.end_cross_product = 0;
    result_curve_points.is_straight_start = result_curve_points.is_straight_end = true;

    return true;
  }

  else if(list_size == 2)
  {
    result_curve_points.start = *start;
    result_curve_points.control1 = *start;
    ++start;
    result_curve_points.control2 = *start;
    result_curve_points.end = *start;
    result_curve_points.start_cross_product = 0;
    result_curve_points.end_cross_product = 0;
    result_curve_points.is_straight_start = result_curve_points.is_straight_end = true;

    return true;
  }

  double* x = new double[list_size + 1];
  double* y = new double[list_size + 1];

  P_ITER temp_iter = start;
  for(int i = 0 ; i < list_size ; ++i, ++temp_iter)
  {
    x[i] = temp_iter->x;
    y[i] = temp_iter->y;
  }

  if(end != freeDrawingPointList.end())
  {
    x[list_size] = end->x;
    y[list_size] = end->y;
  }

  double* interval = new double[list_size];
  interval[0] = 0;

  for(int i = 1 ; i < list_size ; ++i)
  {
    double diff_x = x[i] - x[i - 1];
    double diff_y = y[i] - y[i - 1];
    interval[i] = sqrt(diff_x * diff_x + diff_y * diff_y) + interval[i - 1];
  }


  double** mat_A_for_x = new double*[list_size];
  double* temp = new double[list_size * 4];
  for(int i = 0 ; i < list_size ; ++i)
    mat_A_for_x[i] = temp + 4 * i;
  double** mat_A_for_y = new double*[list_size];
  temp = new double[list_size * 4];
  for(int i = 0 ; i < list_size ; ++i)
    mat_A_for_y[i] = temp + 4 * i;


  for(int i = 0 ; i < list_size ; ++i)
  {
    double ratio_x = interval[i] / interval[list_size - 1];
    double ratio_y = interval[i] / interval[list_size - 1];

    mat_A_for_x[i][3] = 1;
    mat_A_for_x[i][2] = mat_A_for_x[i][3] * ratio_x;
    mat_A_for_x[i][1] = mat_A_for_x[i][2] * ratio_x;
    mat_A_for_x[i][0] = mat_A_for_x[i][1] * ratio_x;

    mat_A_for_y[i][3] = 1;
    mat_A_for_y[i][2] = mat_A_for_y[i][3] * ratio_y;
    mat_A_for_y[i][1] = mat_A_for_y[i][2] * ratio_y;
    mat_A_for_y[i][0] = mat_A_for_y[i][1] * ratio_y;
  }

  double mat_A_ast_A_for_x[4][4];
  double mat_A_ast_A_for_y[4][4];

  for(int i = 0 ; i < 4 ; ++i)
    for(int j = 0 ; j < 4 ; ++j)
    {
      mat_A_ast_A_for_x[i][j] = 0;
      mat_A_ast_A_for_y[i][j] = 0;
      for(int k = 0 ; k < list_size ; ++k)
      {
        mat_A_ast_A_for_x[i][j] += (mat_A_for_x[k][i] * mat_A_for_x[k][j]);
        mat_A_ast_A_for_y[i][j] += (mat_A_for_y[k][i] * mat_A_for_y[k][j]);
      }
    }

  double mat_A_ast_A_inv_for_x[4][4];
  double mat_A_ast_A_inv_for_y[4][4];

  bool is_able_to_approximate_for_x = MakeInverse(mat_A_ast_A_inv_for_x, mat_A_ast_A_for_x);
  bool is_able_to_approximate_for_y = MakeInverse(mat_A_ast_A_inv_for_y, mat_A_ast_A_for_y);

  double** mat_A_ast_A_inv_A_ast_for_x = new double*[4];
  temp = new double[4 * list_size];
  for(int i = 0 ; i < 4 ; ++i)
    mat_A_ast_A_inv_A_ast_for_x[i] = temp + list_size * i;
  double** mat_A_ast_A_inv_A_ast_for_y = new double*[4];
  temp = new double[4 * list_size];
  for(int i = 0 ; i < 4 ; ++i)
    mat_A_ast_A_inv_A_ast_for_y[i] = temp + list_size * i;

  bool is_qualified = false;
  if(is_able_to_approximate_for_x && is_able_to_approximate_for_y)
  {
    for(int i = 0 ; i < 4 ; ++i)
      for(int j = 0 ; j < list_size ; ++j)
      {
        mat_A_ast_A_inv_A_ast_for_x[i][j] = 0;
        mat_A_ast_A_inv_A_ast_for_y[i][j] = 0;
        for(int k = 0 ; k < 4 ; ++k)
        {
          mat_A_ast_A_inv_A_ast_for_x[i][j] += (mat_A_ast_A_inv_for_x[i][k] * mat_A_for_x[j][k]);
          mat_A_ast_A_inv_A_ast_for_y[i][j] += (mat_A_ast_A_inv_for_y[i][k] * mat_A_for_y[j][k]);
        }
      }

    double result_x[4];
    double result_y[4];

    for(int i = 0 ; i < 4 ; ++i)
		{
      result_x[i] = 0;
      result_y[i] = 0;
      for(int j = 0 ; j < list_size ; ++j)
      {
        result_x[i] += (mat_A_ast_A_inv_A_ast_for_x[i][j] * x[j]);
        result_y[i] += (mat_A_ast_A_inv_A_ast_for_y[i][j] * y[j]);
      }
    }


		is_qualified = (CalculateError((const double**)mat_A_for_x, result_x, x, list_size) < Error_Threshold &&
                    CalculateError((const double**)mat_A_for_y, result_y, y, list_size) < Error_Threshold
										);

		if(is_qualified)
		{
			double start_diff1_x = x[1] - x[0];
			double start_diff1_y = y[1] - y[0];
			double start_diff2_x = x[3] - x[0];
			double start_diff2_y = y[3] - y[0];
			double end_diff1_x = x[list_size - 2] - x[list_size - 3];
			double end_diff1_y = y[list_size - 2] - y[list_size - 3];
			double end_diff2_x = x[list_size - 1] - x[list_size - 2];
			double end_diff2_y = y[list_size - 1] - y[list_size - 2];

			double start_cross_product = abs(start_diff1_x * start_diff2_y - start_diff2_x * start_diff1_y);
			double end_cross_product = abs(end_diff1_x * end_diff2_y - end_diff2_x * end_diff1_y);

			double temp1 = sqrt(start_diff1_x * start_diff1_x + start_diff1_y * start_diff1_y +
												 start_diff2_x * start_diff2_x + start_diff2_y * start_diff2_y);
			double temp2 = sqrt(end_diff1_x * end_diff1_x + end_diff1_y * end_diff1_y +
												 end_diff2_x * end_diff2_x + end_diff2_y * end_diff2_y);

			result_curve_points.is_straight_start = false;
			result_curve_points.is_straight_end = false;
			if(start_cross_product / temp1 < 0.5)
				result_curve_points.is_straight_start = true;
			if(end_cross_product / temp2 < 0.5)
				result_curve_points.is_straight_end = true;

			if(end != freeDrawingPointList.end())
			{
				end_diff1_x = x[list_size] - x[list_size - 1];
				end_diff1_y = y[list_size] - y[list_size - 1];
				end_cross_product = abs(end_diff2_x * end_diff1_y - end_diff1_x * end_diff2_y);
				temp2 = sqrt(end_diff1_x * end_diff1_x + end_diff1_y * end_diff1_y +
										 end_diff2_x * end_diff2_x + end_diff2_y * end_diff2_y);
				result_curve_points.end_sin_degree = end_cross_product / temp2;
			}



			start_diff1_x = x[1] - x[0];
			start_diff1_y = y[1] - y[0];
			start_diff2_x = x[2] - x[1];
			start_diff2_y = y[2] - y[1];
			end_diff1_x = x[list_size - 2] - x[list_size - 3];
			end_diff1_y = y[list_size - 2] - y[list_size - 3];
			end_diff2_x = x[list_size - 1] - x[list_size - 2];
			end_diff2_y = y[list_size - 1] - y[list_size - 2];

			result_curve_points.start_cross_product = abs(start_diff1_x * start_diff2_y - start_diff2_x * start_diff1_y);
			result_curve_points.end_cross_product = abs(end_diff1_x * end_diff2_y - end_diff2_x * end_diff1_y);


			MakeResultToPoint(result_x);
			MakeResultToPoint(result_y);

			result_curve_points.start.DPoint(result_x[3], result_y[3]);
			result_curve_points.control1.DPoint(result_x[2], result_y[2]);
			result_curve_points.control2.DPoint(result_x[1], result_y[1]);


/*      int real_start = start_index - ei;
			result_curve_points.start.DPoint(x[ei], y[ei]);

			double start_ratio = interval[ei] / interval[list_size - 1];
			double temp_1_2_x = result_x[2] + (result_x[1] - result_x[2]) * start_ratio;
			double temp_1_2_y = result_y[2] + (result_y[1] - result_y[2]) * start_ratio;
			double temp_2_end_x = result_x[1] + (result_x[0] - result_x[1]) * start_ratio;
			double temp_2_end_y = result_y[1] + (result_y[0] - result_y[1]) * start_ratio;
			double control1_by_real_start_x = temp_1_2_x + (temp_2_end_x - temp_1_2_x) * start_ratio;
			double control1_by_real_start_y = temp_1_2_y + (temp_2_end_y - temp_1_2_y) * start_ratio;

			result_curve_points.control1.DPoint(control1_by_real_start_x, control1_by_real_start_y);
			result_curve_points.control2.DPoint(temp_2_end_x, temp_2_end_y);*/
			result_curve_points.end.DPoint(result_x[0], result_y[0]);
//      result_curve_points.end.DPoint(x[list_size - 1], y[list_size - 1]);
		}
		else {
			int error_x = CalculateError((const double**)mat_A_for_x, result_x, x, list_size);
			int error_y = CalculateError((const double**)mat_A_for_y, result_y, y, list_size);
		}
	}


	delete[] x;
	delete[] y;
	delete[] interval;
	delete[] mat_A_for_x[0];
	delete[] mat_A_for_x;
	delete[] mat_A_for_y[0];
	delete[] mat_A_for_y;
  delete[] mat_A_ast_A_inv_A_ast_for_x[0];
  delete[] mat_A_ast_A_inv_A_ast_for_x;
  delete[] mat_A_ast_A_inv_A_ast_for_y[0];
  delete[] mat_A_ast_A_inv_A_ast_for_y;

  return is_qualified;
}
//---------------------------------------------------------------------------
double BezierCurveApproximator::CalculateError(const double** t_val, const double* result,
                                               const double* original, const int size) const
{
  double error_value = 0;

  for(int i = 0 ; i < size ; ++i)
  {
    double approximated = 0;
    for(int j = 0 ; j < 4 ; ++j)
      approximated += t_val[i][j] * result[j];
    error_value += (int)((approximated - original[i]) * (approximated - original[i]));
  }
  return error_value;
}
//---------------------------------------------------------------------------
bool BezierCurveApproximator::MakeInverse(double inv[4][4], const double x[4][4]) const
{
  inv[0][0] =
  x[1][1]*x[2][2]*x[3][3] + x[1][2]*x[2][3]*x[3][1] + x[1][3]*x[2][1]*x[3][2] - x[1][1]*x[2][3]*x[3][2] - x[1][2]*x[2][1]*x[3][3] - x[1][3]*x[2][2]*x[3][1];
  inv[0][1] =
  x[0][1]*x[2][3]*x[3][2] + x[0][2]*x[2][1]*x[3][3] + x[0][3]*x[2][2]*x[3][1] - x[0][1]*x[2][2]*x[3][3] - x[0][2]*x[2][3]*x[3][1] - x[0][3]*x[2][1]*x[3][2];
  inv[0][2] =
  x[0][1]*x[1][2]*x[3][3] + x[0][2]*x[1][3]*x[3][1] + x[0][3]*x[1][1]*x[3][2] - x[0][1]*x[1][3]*x[3][2] - x[0][2]*x[1][1]*x[3][3] - x[0][3]*x[1][2]*x[3][1];
  inv[0][3] =
  x[0][1]*x[1][3]*x[2][2] + x[0][2]*x[1][1]*x[2][3] + x[0][3]*x[1][2]*x[2][1] - x[0][1]*x[1][2]*x[2][3] - x[0][2]*x[1][3]*x[2][1] - x[0][3]*x[1][1]*x[2][2];
  inv[1][0] =
  x[1][0]*x[2][3]*x[3][2] + x[1][2]*x[2][0]*x[3][3] + x[1][3]*x[2][2]*x[3][0] - x[1][0]*x[2][2]*x[3][3] - x[1][2]*x[2][3]*x[3][0] - x[1][3]*x[2][0]*x[3][2];
  inv[1][1] =
  x[0][0]*x[2][2]*x[3][3] + x[0][2]*x[2][3]*x[3][0] + x[0][3]*x[2][0]*x[3][2] - x[0][0]*x[2][3]*x[3][2] - x[0][2]*x[2][0]*x[3][3] - x[0][3]*x[2][2]*x[3][0];
  inv[1][2] =
  x[0][0]*x[1][3]*x[3][2] + x[0][2]*x[1][0]*x[3][3] + x[0][3]*x[1][2]*x[3][0] - x[0][0]*x[1][2]*x[3][3] - x[0][2]*x[1][3]*x[3][0] - x[0][3]*x[1][0]*x[3][2];
  inv[1][3] =
  x[0][0]*x[1][2]*x[2][3] + x[0][2]*x[1][3]*x[2][0] + x[0][3]*x[1][0]*x[2][2] - x[0][0]*x[1][3]*x[2][2] - x[0][2]*x[1][0]*x[2][3] - x[0][3]*x[1][2]*x[2][0];
  inv[2][0] =
  x[1][0]*x[2][1]*x[3][3] + x[1][1]*x[2][3]*x[3][0] + x[1][3]*x[2][0]*x[3][1] - x[1][0]*x[2][3]*x[3][1] - x[1][1]*x[2][0]*x[3][3] - x[1][3]*x[2][1]*x[3][0];
  inv[2][1] =
  x[0][0]*x[2][3]*x[3][1] + x[0][1]*x[2][0]*x[3][3] + x[0][3]*x[2][1]*x[3][0] - x[0][0]*x[2][1]*x[3][3] - x[0][1]*x[2][3]*x[3][0] - x[0][3]*x[2][0]*x[3][1];
  inv[2][2] =
  x[0][0]*x[1][1]*x[3][3] + x[0][1]*x[1][3]*x[3][0] + x[0][3]*x[1][0]*x[3][1] - x[0][0]*x[1][3]*x[3][1] - x[0][1]*x[1][0]*x[3][3] - x[0][3]*x[1][1]*x[3][0];
  inv[2][3] =
  x[0][0]*x[1][3]*x[2][1] + x[0][1]*x[1][0]*x[2][3] + x[0][3]*x[1][1]*x[2][0] - x[0][0]*x[1][1]*x[2][3] - x[0][1]*x[1][3]*x[2][0] - x[0][3]*x[1][0]*x[2][1];
  inv[3][0] =
  x[1][0]*x[2][2]*x[3][1] + x[1][1]*x[2][0]*x[3][2] + x[1][2]*x[2][1]*x[3][0] - x[1][0]*x[2][1]*x[3][2] - x[1][1]*x[2][2]*x[3][0] - x[1][2]*x[2][0]*x[3][1];
  inv[3][1] =
  x[0][0]*x[2][1]*x[3][2] + x[0][1]*x[2][2]*x[3][0] + x[0][2]*x[2][0]*x[3][1] - x[0][0]*x[2][2]*x[3][1] - x[0][1]*x[2][0]*x[3][2] - x[0][2]*x[2][1]*x[3][0];
  inv[3][2] =
  x[0][0]*x[1][2]*x[3][1] + x[0][1]*x[1][0]*x[3][2] + x[0][2]*x[1][1]*x[3][0] - x[0][0]*x[1][1]*x[3][2] - x[0][1]*x[1][2]*x[3][0] - x[0][2]*x[1][0]*x[3][1];
  inv[3][3] =
  x[0][0]*x[1][1]*x[2][2] + x[0][1]*x[1][2]*x[2][0] + x[0][2]*x[1][0]*x[2][1] - x[0][0]*x[1][2]*x[2][1] - x[0][1]*x[1][0]*x[2][2] - x[0][2]*x[1][1]*x[2][0];

  double det = MakeDet(x);
  if(det == 0)
    return false;

  for(int i = 0 ; i < 4 ; ++i)
    for(int j = 0 ; j < 4 ; ++j)
      inv[i][j] /= det;

  return true;
}
//---------------------------------------------------------------------------
double BezierCurveApproximator::MakeDet(const double x[4][4]) const
{
  double det =
  x[0][0]*x[1][1]*x[2][2]*x[3][3] + x[0][0]*x[1][2]*x[2][3]*x[3][1] + x[0][0]*x[1][3]*x[2][1]*x[3][2] +
  x[0][1]*x[1][0]*x[2][3]*x[3][2] + x[0][1]*x[1][2]*x[2][0]*x[3][3] + x[0][1]*x[1][3]*x[2][2]*x[3][0] +
  x[0][2]*x[1][0]*x[2][1]*x[3][3] + x[0][2]*x[1][1]*x[2][3]*x[3][0] + x[0][2]*x[1][3]*x[2][0]*x[3][1] +
  x[0][3]*x[1][0]*x[2][2]*x[3][1] + x[0][3]*x[1][1]*x[2][0]*x[3][2] + x[0][3]*x[1][2]*x[2][1]*x[3][0] -
  x[0][0]*x[1][1]*x[2][3]*x[3][2] - x[0][0]*x[1][2]*x[2][1]*x[3][3] - x[0][0]*x[1][3]*x[2][2]*x[3][1] -
  x[0][1]*x[1][0]*x[2][2]*x[3][3] - x[0][1]*x[1][2]*x[2][3]*x[3][0] - x[0][1]*x[1][3]*x[2][0]*x[3][2] -
  x[0][2]*x[1][0]*x[2][3]*x[3][1] - x[0][2]*x[1][1]*x[2][0]*x[3][3] - x[0][2]*x[1][3]*x[2][1]*x[3][0] -
  x[0][3]*x[1][0]*x[2][1]*x[3][2] - x[0][3]*x[1][1]*x[2][2]*x[3][0] - x[0][3]*x[1][2]*x[2][0]*x[3][1] ;

  return det;
}
//---------------------------------------------------------------------------
void BezierCurveApproximator::MakeResultToPoint(double result[4]) const
{
  result[2] /= 3;
  result[1] /= 3;
  result[0] += result[1];
  result[1] += result[2];
  result[2] += result[3];
  result[0] += result[1];
  result[1] += result[2];
  result[0] += result[1];
}
//---------------------------------------------------------------------------

void BezierCurveApproximator::Reset()
{
  freeDrawingPointList.clear();
  bezierCurvePointsList.clear();
}
//---------------------------------------------------------------------------
void BezierCurveApproximator::AddPoint(const DPOINT& point)
{
  DPOINT last = *(freeDrawingPointList.rbegin());
  if(last.x != point.x || last.y != point.y)
    freeDrawingPointList.push_back(point);
}
//---------------------------------------------------------------------------
void BezierCurveApproximator::MakeBezierCurvePointsList(const int start_index, const int end_index)
{
//  const int size = freeDrawingPointList.size();
  const int size = end_index;
  int low_search_bound, high_search_bound, start, end;
//  low_search_bound = start = 0;
//  high_search_bound = end = size;

  low_search_bound = start = start_index;
  high_search_bound = end = end_index;

	int i = 0;
	while(start < end)
	{
		BezierCurvePoints result_curve_points;
//    if(high_search_bound - low_search_bound > 100)
//      end = high_search_bound  = low_search_bound + 101;
		while(1)
		{
			if (i++ % 100 == 1)
				Application->ProcessMessages();
			int next_bounce = (high_search_bound - low_search_bound) / 2;
      if(start + 1 == end && end == size) break;
			if(Approximate(start, end, result_curve_points))
      {
        break;                                              // ¹Ù·Î break Çصµ »ó°ü ¾øÀ»µí by tmddn85 2009.11.24
/*        if(next_bounce == 0 || end == size)
          break;
        low_search_bound = high_search_bound;
        high_search_bound += next_bounce;
        if(high_search_bound > size)
          high_search_bound = size;
        end = high_search_bound;
*/
      }
      else
      {
				if(next_bounce == 0)
          break;

        high_search_bound -= next_bounce;
        end = high_search_bound;
      }
    }
//    if(start + 1 == end && end == size)
//      break;
    bezierCurvePointsList.push_back(result_curve_points);
    if(end == size) break;
    low_search_bound = start = end - 1;
    high_search_bound = end = size;
  }
}

bool BezierCurveApproximator::Approximate2(const int start_index, const int end_index)
{
  P_ITER start, end;
  start = end = freeDrawingPointList.begin();
  std::advance(start, start_index);
  std::advance(end, end_index);

  int list_size = end_index - start_index;

  if(list_size < 2)
    return true;
  if(list_size == 3)
  {
    BezierCurvePoints result_curve_points;

    result_curve_points.start = *start;
    ++start;
    result_curve_points.control1 = *start;
    ++start;
    result_curve_points.control2 = *start;
    result_curve_points.end = *start;

    bezierCurvePointsList.push_back(result_curve_points);

    return true;
  }
  else if(list_size == 2)
  {
    BezierCurvePoints result_curve_points;

    result_curve_points.start = *start;
    result_curve_points.control1 = *start;
    ++start;
    result_curve_points.control2 = *start;
    result_curve_points.end = *start;

    bezierCurvePointsList.push_back(result_curve_points);

    return true;
  }

  double* x = new double[list_size];
  double* y = new double[list_size];

  P_ITER temp_iter = start;
  for(int i = 0 ; i < list_size ; ++i, ++temp_iter)
  {
    x[i] = temp_iter->x;
    y[i] = temp_iter->y;
  }

  double* interval = new double[list_size];
  interval[0] = 0;

  for(int i = 1 ; i < list_size ; ++i)
  {
    double diff_x = x[i] - x[i - 1];
    double diff_y = y[i] - y[i - 1];
    interval[i] = sqrt(diff_x * diff_x + diff_y * diff_y) + interval[i - 1];
  }


  double** mat_A_for_x = new double*[list_size];
  double* temp = new double[list_size * 4];
  for(int i = 0 ; i < list_size ; ++i)
    mat_A_for_x[i] = temp + 4 * i;
  double** mat_A_for_y = new double*[list_size];
  temp = new double[list_size * 4];
  for(int i = 0 ; i < list_size ; ++i)
    mat_A_for_y[i] = temp + 4 * i;


  for(int i = 0 ; i < list_size ; ++i)
  {
    double ratio_x = interval[i] / interval[list_size - 1];
    double ratio_y = interval[i] / interval[list_size - 1];

    mat_A_for_x[i][3] = 1;
    mat_A_for_x[i][2] = mat_A_for_x[i][3] * ratio_x;
    mat_A_for_x[i][1] = mat_A_for_x[i][2] * ratio_x;
    mat_A_for_x[i][0] = mat_A_for_x[i][1] * ratio_x;

    mat_A_for_y[i][3] = 1;
    mat_A_for_y[i][2] = mat_A_for_y[i][3] * ratio_y;
    mat_A_for_y[i][1] = mat_A_for_y[i][2] * ratio_y;
    mat_A_for_y[i][0] = mat_A_for_y[i][1] * ratio_y;
  }

  double mat_A_ast_A_for_x[4][4];
  double mat_A_ast_A_for_y[4][4];

  for(int i = 0 ; i < 4 ; ++i)
    for(int j = 0 ; j < 4 ; ++j)
    {
      mat_A_ast_A_for_x[i][j] = 0;
      mat_A_ast_A_for_y[i][j] = 0;
      for(int k = 0 ; k < list_size ; ++k)
      {
        mat_A_ast_A_for_x[i][j] += (mat_A_for_x[k][i] * mat_A_for_x[k][j]);
        mat_A_ast_A_for_y[i][j] += (mat_A_for_y[k][i] * mat_A_for_y[k][j]);
      }
    }

  double mat_A_ast_A_inv_for_x[4][4];
  double mat_A_ast_A_inv_for_y[4][4];

  bool is_able_to_approximate_for_x = MakeInverse(mat_A_ast_A_inv_for_x, mat_A_ast_A_for_x);
  bool is_able_to_approximate_for_y = MakeInverse(mat_A_ast_A_inv_for_y, mat_A_ast_A_for_y);

  double** mat_A_ast_A_inv_A_ast_for_x = new double*[4];
  temp = new double[4 * list_size];
  for(int i = 0 ; i < 4 ; ++i)
    mat_A_ast_A_inv_A_ast_for_x[i] = temp + list_size * i;
  double** mat_A_ast_A_inv_A_ast_for_y = new double*[4];
  temp = new double[4 * list_size];
  for(int i = 0 ; i < 4 ; ++i)
    mat_A_ast_A_inv_A_ast_for_y[i] = temp + list_size * i;

  bool is_qualified = false;
  if(is_able_to_approximate_for_x && is_able_to_approximate_for_y)
  {
    for(int i = 0 ; i < 4 ; ++i)
      for(int j = 0 ; j < list_size ; ++j)
      {
        mat_A_ast_A_inv_A_ast_for_x[i][j] = 0;
        mat_A_ast_A_inv_A_ast_for_y[i][j] = 0;
        for(int k = 0 ; k < 4 ; ++k)
        {
          mat_A_ast_A_inv_A_ast_for_x[i][j] += (mat_A_ast_A_inv_for_x[i][k] * mat_A_for_x[j][k]);
          mat_A_ast_A_inv_A_ast_for_y[i][j] += (mat_A_ast_A_inv_for_y[i][k] * mat_A_for_y[j][k]);
        }
      }

    double result_x[4];
    double result_y[4];

    for(int i = 0 ; i < 4 ; ++i)
    {
      result_x[i] = 0;
      result_y[i] = 0;
      for(int j = 0 ; j < list_size ; ++j)
      {
        result_x[i] += (mat_A_ast_A_inv_A_ast_for_x[i][j] * x[j]);
        result_y[i] += (mat_A_ast_A_inv_A_ast_for_y[i][j] * y[j]);
      }
    }


		bool is_qualified = (CalculateError((const double**)mat_A_for_x, result_x, x, list_size) < Error_Threshold + 1 &&
												 CalculateError((const double**)mat_A_for_y, result_y, y, list_size) < Error_Threshold + 1
                         );

    if(is_qualified)
		{
			MakeResultToPoint(result_x);
      MakeResultToPoint(result_y);

      BezierCurvePoints result_curve_points;

      result_curve_points.start.DPoint(result_x[3], result_y[3]);
      result_curve_points.control1.DPoint(result_x[2], result_y[2]);
      result_curve_points.control2.DPoint(result_x[1], result_y[1]);
      result_curve_points.end.DPoint(result_x[0], result_y[0]);

      bezierCurvePointsList.push_back(result_curve_points);
    }
    else
    {
			int half_index = (end_index + start_index) / 2 + 1;
      Approximate2(start_index, half_index);
      Approximate2(half_index, end_index);
    }
  }


  delete[] x;
  delete[] y;
  delete[] interval;
  delete[] mat_A_for_x[0];
  delete[] mat_A_for_x;
  delete[] mat_A_for_y[0];
  delete[] mat_A_for_y;
  delete[] mat_A_ast_A_inv_A_ast_for_x[0];
  delete[] mat_A_ast_A_inv_A_ast_for_x;
  delete[] mat_A_ast_A_inv_A_ast_for_y[0];
  delete[] mat_A_ast_A_inv_A_ast_for_y;

  return is_qualified;
}