#include <limits.h>
#include <math.h>
#include <isolines.h>
#include <ps.h>
#include <js.h>

#define MIN(A,B) (A < B ? A : B)
#define MAX(A,B) (A > B ? A : B)

iso_triangle_list* iso_split_cell_in_triangles(dyg_node* node, dyg_dimensions* dim)
{
  iso_coords center;
  iso_triangle_list *ret = 0, *aux = 0;

  center.axis[X] = (dim->min[X] + dim->max[X]) * 0.5;
  center.axis[Y] = (dim->min[Y] + dim->max[Y]) * 0.5;

  //  center.axis[Z] = O valor de f(x, y) quem eh f?

  aux = (iso_triangle_list*)malloc(sizeof(iso_triangle_list));
  aux->triangle =(iso_triangle*)malloc(sizeof(iso_triangle));

  aux->triangle->P.axis[X] = dim->min[X];
  aux->triangle->P.axis[Y] = dim->min[Y];

  aux->triangle->Q.axis[X] = center.axis[X];
  aux->triangle->Q.axis[Y] = dim->min[Y];

  aux->triangle->R.axis[X] = center.axis[X];
  aux->triangle->R.axis[Y] = center.axis[Y];

  aux->next = ret;
  ret = aux;

  aux = (iso_triangle_list*)malloc(sizeof(iso_triangle_list));
  aux->triangle =(iso_triangle*)malloc(sizeof(iso_triangle));

  aux->triangle->P.axis[X] = center.axis[X];
  aux->triangle->P.axis[Y] = dim->min[Y];

  aux->triangle->Q.axis[X] = dim->max[X];
  aux->triangle->Q.axis[Y] = dim->min[Y];

  aux->triangle->R.axis[X] = center.axis[X];
  aux->triangle->R.axis[Y] = center.axis[Y];

  aux->next = ret;
  ret = aux;

  aux = (iso_triangle_list*)malloc(sizeof(iso_triangle_list));
  aux->triangle =(iso_triangle*)malloc(sizeof(iso_triangle));

  aux->triangle->P.axis[X] = dim->max[X];
  aux->triangle->P.axis[Y] = dim->min[Y];

  aux->triangle->Q.axis[X] = dim->max[X];
  aux->triangle->Q.axis[Y] = center.axis[Y];

  aux->triangle->R.axis[X] = center.axis[X];
  aux->triangle->R.axis[Y] = center.axis[Y];

  aux->next = ret;
  ret = aux;

  aux = (iso_triangle_list*)malloc(sizeof(iso_triangle_list));
  aux->triangle =(iso_triangle*)malloc(sizeof(iso_triangle));

  aux->triangle->P.axis[X] = dim->max[X];
  aux->triangle->P.axis[Y] = center.axis[Y];

  aux->triangle->Q.axis[X] = dim->max[X];
  aux->triangle->Q.axis[Y] = dim->max[Y];

  aux->triangle->R.axis[X] = center.axis[X];
  aux->triangle->R.axis[Y] = center.axis[Y];

  aux->next = ret;
  ret = aux;

  aux = (iso_triangle_list*)malloc(sizeof(iso_triangle_list));
  aux->triangle =(iso_triangle*)malloc(sizeof(iso_triangle));

  aux->triangle->P.axis[X] = dim->max[X];
  aux->triangle->P.axis[Y] = dim->max[Y];

  aux->triangle->Q.axis[X] = center.axis[X];
  aux->triangle->Q.axis[Y] = dim->max[Y];

  aux->triangle->R.axis[X] = center.axis[X];
  aux->triangle->R.axis[Y] = center.axis[Y];

  aux->next = ret;
  ret = aux;

  aux = (iso_triangle_list*)malloc(sizeof(iso_triangle_list));
  aux->triangle =(iso_triangle*)malloc(sizeof(iso_triangle));

  aux->triangle->P.axis[X] = center.axis[X];
  aux->triangle->P.axis[Y] = dim->max[Y];

  aux->triangle->Q.axis[X] = dim->min[X];
  aux->triangle->Q.axis[Y] = dim->max[Y];

  aux->triangle->R.axis[X] = center.axis[X];
  aux->triangle->R.axis[Y] = center.axis[Y];

  aux->next = ret;
  ret = aux;

  aux = (iso_triangle_list*)malloc(sizeof(iso_triangle_list));
  aux->triangle =(iso_triangle*)malloc(sizeof(iso_triangle));

  aux->triangle->P.axis[X] = dim->min[X];
  aux->triangle->P.axis[Y] = dim->max[Y];

  aux->triangle->Q.axis[X] = dim->min[X];
  aux->triangle->Q.axis[Y] = center.axis[Y];

  aux->triangle->R.axis[X] = center.axis[X];
  aux->triangle->R.axis[Y] = center.axis[Y];

  aux->next = ret;
  ret = aux;

  aux = (iso_triangle_list*)malloc(sizeof(iso_triangle_list));
  aux->triangle =(iso_triangle*)malloc(sizeof(iso_triangle));

  aux->triangle->P.axis[X] = dim->min[X];
  aux->triangle->P.axis[Y] = center.axis[Y];

  aux->triangle->Q.axis[X] = dim->min[X];
  aux->triangle->Q.axis[Y] = dim->min[Y];

  aux->triangle->R.axis[X] = center.axis[X];
  aux->triangle->R.axis[Y] = center.axis[Y];

  aux->next = ret;
  ret = aux;

  return ret;
}

void iso_lines_rec
      (
      FILE* f, 
      dyg_node* node, dyg_node_list* list, dyg_dimensions* dim, int inicial_axis, int param, 
      double fStep, double fStart, bool isLinear
     )
{
  iso_triangle_list* triangles, *tritemp;
  dyg_dimensions ret_dim;
  int log2;

  dyg_get_cell_dimension( node->index, dim, inicial_axis, &ret_dim);
  if( !iso_split_criterion(&ret_dim) && node->c0 == NULL && node->c1 == NULL)
    {
      triangles = iso_split_cell_in_triangles(node, &ret_dim);
      tritemp = 0;
      for( ; triangles != NULL; triangles = triangles->next)
	{
	  if( tritemp )
	    {
	      free(tritemp->triangle);
	      free(tritemp);
	    }
	  triangles->triangle->P.axis[Z] = iso_avail(node, list, &triangles->triangle->P, param, dim, inicial_axis, false);
	  triangles->triangle->Q.axis[Z] = iso_avail(node, list, &triangles->triangle->Q, param, dim, inicial_axis, false);
	  triangles->triangle->R.axis[Z] = iso_avail(node, list, &triangles->triangle->R, param, dim, inicial_axis, false);

	  iso_process_triangle(f, fStart, fStep, isLinear, triangles->triangle);

	  tritemp = triangles;
     	}
      if( tritemp )
	{
	  free(tritemp->triangle);
	  free(tritemp);
	}
  }
  else
    {
      if( node->c0 != NULL || node->c1 != NULL )
	{
	  iso_lines_rec(f, node->c0, list, dim, inicial_axis, param, fStep, fStart, isLinear);
	  iso_lines_rec(f, node->c1, list, dim, inicial_axis, param, fStep, fStart, isLinear);
	}
      else
	{
	  log2 = (int)(log(node->index) / log(2));
	  if( log2 )
	    log2 += (int)floor(node->index / pow(2, log2 + 1));

	  dyg_shatter_node( node, MAX_INT, iso_split_criterion, iso_virtual_data_adjust, &param, dim, 
			    (log2 % 2) ? next_axis(inicial_axis) : inicial_axis, true);

	  iso_lines_rec(f, node->c0, list, dim, inicial_axis, param, fStep, fStart, isLinear );
	  iso_lines_rec(f, node->c1, list, dim, inicial_axis, param, fStep, fStart, isLinear );
	  
	  dyg_free_tree_rec(node->c0);
	  dyg_free_tree_rec(node->c1);
	  
	  node->c0 = NULL;
	  node->c1 = NULL;
	}
    }
}

void iso_isolines
     (
      FILE* f, 
      dyg_header* header, int param, 
      double fStep, double fStart, bool isLinear
     )
{
  dyg_node_list* list;
  affirm(header, "Header is empty");
  affirm(header->root, "Tree is empty");

 list = iso_find_leafs(header->root, header->inicial_dimension, header->inicial_axis, 1, MAX_INT, false);  
 iso_lines_rec(f, header->root, list, header->inicial_dimension, header->inicial_axis, param,
		fStep, fStart, isLinear);
 dyg_flush_list(list);
}

void iso_process_triangle
     (
      FILE* f, 
      double fStart, double fStep, bool isLinear,
      iso_triangle* T
     )
{

  int kMax, kMin, k;
  double fMin, fMax, fLevel;
  /*
  ps_set_pen(f, 0.0, 0.0, 1.0, 0.1, 0.0, 0.0);
 	ps_draw_segment(f, T->P.axis[X], T->P.axis[Y], T->Q.axis[X], T->Q.axis[Y]);
	ps_draw_segment(f, T->R.axis[X], T->R.axis[Y], T->Q.axis[X], T->Q.axis[Y]);
	ps_draw_segment(f, T->P.axis[X], T->P.axis[Y], T->R.axis[X], T->R.axis[Y]);
  ps_set_pen(f, 0.0, 0.0, 0.0, 0.1, 0.0, 0.0);
  */
  fMin = MIN(MIN(T->P.axis[Z], T->Q.axis[Z]), T->R.axis[Z]);
  fMax = MAX(MAX(T->P.axis[Z], T->Q.axis[Z]), T->R.axis[Z]);

  kMin = (int)(floor((MAX(0.0, fMin) - fStart) / fStep)) + 1;
  kMax = (int)(ceil((fMax - fStart) / fStep));

  for( k = kMin; k < kMax; k++ )
    {
      if( isLinear )
	      fLevel = fStart + k * fStep;
      else
	      fLevel = fStart + pow(fStep, k);

      iso_plot_zeros_in_triangle
	    (
	       f, 
	       T->P.axis[X], T->P.axis[Y], T->P.axis[Z] - fLevel,
	       T->Q.axis[X], T->Q.axis[Y], T->Q.axis[Z] - fLevel,
	       T->R.axis[X], T->R.axis[Y], T->R.axis[Z] - fLevel
      );
    }
}

void iso_plot_zeros_in_triangle
     (
      FILE* f, 				
      double Px, double Py, double Pf,
      double Qx, double Qy, double Qf,
      double Rx, double Ry, double Rf
      )
{
  /*
    Given a triangle "P,Q,R" on the plane,
    and function values at its vertices, interpolates a 
    linear function through that data, and plots the set of 
    points where that function is zero. 
  */

  double temp;
  double SPR, SRP, PRx, PRy, SQR, SRQ, QRx, QRy;

  if( Pf == 0.0 && Qf == 0.0 && Rf == 0.0 )
    return;
     /*
        f.segment(Px, Py, Qx, Qy);
        f.segment(Qx, Qy, Rx, Ry);
        f.segment(Rx, Ry, Px, Py);
        f.setFillColor(PSPlot.Black);
        f.triangle(Px, Py, Qx, Qy, Rx, Ry);
     */

  if( Pf == 0.0 && Qf == 0.0 )
    ps_draw_segment(f,Px, Py, Qx, Qy);
  if( Qf == 0.0 && Rf == 0.0 )
    ps_draw_segment(f,Qx, Qy, Rx, Ry);
  if( Pf == 0.0 && Rf == 0.0 )
    ps_draw_segment(f,Px, Py, Rx, Ry);
  if( ((Pf == 0.0) + (Qf == 0.0) + (Rf == 0.0)) >= 2 )
      return;
  if( Pf >= 0.0 && Qf >= 0.0 && Rf >= 0.0 )
      return;
  if( Pf <= 0.0 && Qf <= 0.0 && Rf <= 0.0 )
      return;

  /* Now there is at least one positive corner and one negative corner */
  /* Permute corners so that Rf <= Pf <= Qf: */
  
  if( Rf > Pf )
    {
        temp = Px; Px = Rx; Rx = temp;      
        temp = Py; Py = Ry; Ry = temp;
        temp = Pf; Pf = Rf; Rf = temp;
    }
  if( Pf > Qf )
    {
        temp= Px; Px = Qx; Qx = temp;
        temp= Py; Py = Qy; Qy = temp;
        temp= Pf; Pf = Qf; Qf = temp; 
    }
  if( Rf > Pf )
    {
      temp = Px; Px = Rx; Rx = temp;      
      temp = Py; Py = Ry; Ry = temp;
        temp = Pf; Pf = Rf; Rf = temp;
    }

  affirm( Rf <= Pf, "Rf > Pf" );
  affirm( Pf <= Qf, "Pf > Qf" );
  affirm( Rf < 0.0, "Rf >= 0.0" );
  affirm( Qf > 0.0, "Qf <= 0.0" );

  /* Now swap q,r so that sign(Pf) # sign(Rf): */

  if( Pf < 0.0 )
    {
      temp = Qx; Qx = Rx; Rx = temp;
      temp = Qy; Qy = Ry; Ry = temp;
      temp = Qf; Qf = Rf; Rf = temp; 
    }
  
  /*now Qf # 0.0d0, Rf # 0.0d0 */
  /* also Pf = 0.0d0, OR (sign(Pf) # sig(Qf) AND sig(Pf) # sig(Rf)) */

  affirm( Qf != 0.0, "Qf == 0.0" );
  affirm( Rf != 0.0, "Rf == 0.0" );
  affirm( (Qf > 0.0) != (Rf > 0.0), "Qf sign == Rf sign" );
  affirm( (Pf == 0.0) || ((Pf > 0.0) == (Qf > 0.0)) && ((Pf > 0.0) != (Rf > 0.0)), "Pf sign == Rf sign" );

  if( Pf != 0.0 )
    {
      SPR = Rf / (Rf - Pf);
      SRP = Pf / (Pf - Rf);
      PRx = Px * SPR + Rx * SRP;
      PRy = Py * SPR + Ry * SRP;
      SQR = Rf / (Rf - Qf);
      SRQ = Qf / (Qf - Rf);
      QRx = Qx * SQR + Rx * SRQ;
      QRy = Qy * SQR + Ry * SRQ;
      ps_draw_segment(f,PRx, PRy, QRx, QRy);	  
    }
  else
    {
      SQR = Rf / (Rf - Qf);
      SRQ = Qf / (Qf - Rf);
      QRx = Qx * SQR + Rx * SRQ;
      QRy = Qy * SQR + Ry * SRQ;
      affirm( (Qf > 0.0) != (Rf > 0.0), "Qf sign == Rf sign");
      ps_draw_segment(f,Px, Py, QRx, QRy);	  
    }
}

double iso_avail(dyg_node* node,dyg_node_list* list, iso_coords* q, int param,dyg_dimensions* inicial_dim, int inicial_axis, bool virtual){
  dyg_node* ptr, *root;
  dyg_node_list *aux;
  double ret = 0.0, f;
  int i;
  
  if(node->virtual)
    ptr = *((dyg_node**)(node->data));
  else
    ptr = node;
      
  for( aux = list; aux != NULL; aux = aux->next )
    {
      bool flag = true;
      for( i = 0; i < DIM - 1 && flag; i++ )
	{
	  double di = (aux->dim.max[i] - aux->dim.min[i]) * 0.5;
	  if( aux->dim.max[i] + di < q->axis[i] || aux->dim.min[i] - di > q->axis[i] )
	    flag = false;
	}

      if(!flag)
	continue;

      f = 1.0;
      for(i = 0; i < DIM - 1; i++)
	f *= 1 - (fabs(((aux->dim.max[i] + aux->dim.min[i]) * 0.5) - q->axis[i]) / ((aux->dim.max[i] - aux->dim.min[i])));
      
      ret += ((double*)(aux->node->data))[param] * f;
    }
  return ret;
}

int iso_dummy_data_adjust(dyg_node* node)
{
  return 0;
}

int iso_virtual_data_adjust(dyg_node* node)
{
  if( node->virtual && node->p->virtual )
    {
      node->data = malloc(sizeof(dyg_node*));
      *((dyg_node**)(node->data)) = *((dyg_node**)(node->p->data));
    }
  else if( node->virtual )
    {
      node->data = malloc(sizeof(dyg_node**));
      *((dyg_node**)(node->data)) = node->p;
    }

  return 0;
}

double iso_interpole(double* val, dyg_dimensions* dim, iso_coords* p)
{
  double ty1 = ((p->axis[Y] - dim->min[Y]) * val[2] + (dim->max[Y] - p->axis[Y]) * val[0]) / (dim->max[Y] - dim->min[Y]);
  double ty2 = ((p->axis[Y] - dim->min[Y]) * val[3] + (dim->max[Y] - p->axis[Y]) * val[1]) / (dim->max[Y] - dim->min[Y]);
  return       ((p->axis[X] - dim->min[X]) * ty2    + (dim->max[X] - p->axis[X]) * ty1)    / (dim->max[X] - dim->min[X]);
}

bool iso_split_criterion(dyg_dimensions* dim)
/*
  An arbitrary splitting criterion. */
{
  return GT((dim->max[X] - dim->min[X]) * (dim->max[Y] - dim->min[Y]), 1e-2 );
}

dyg_node_list* iso_find_leafs(dyg_node* root, dyg_dimensions* dim, int naxis, int deep, int max_deep, bool virtual_cells)
{
  dyg_dimensions d;
  dyg_node_list* ret = 0;
  int i;

  if( deep > max_deep )
    return ret;

  if( (root->c0 == NULL && root->c1 == NULL) )
    {
      dyg_insert_node(&ret, root);
      for(i = 0; i < DIM - 1; i++)
	ret->dim.min[i] = dim->min[i], ret->dim.max[i] = dim->max[i];
      ret->next = 0;
    }
  else if(!virtual_cells && root->c0->virtual && root->c1->virtual)// É folha
    {
      dyg_insert_node(&ret, root);
      for(i = 0; i < DIM - 1; i++)
	ret->dim.min[i] = dim->min[i], ret->dim.max[i] = dim->max[i];
      ret->next = 0;
    }
  else
    {
      for(i = 0; i < DIM; i++ )
	d.max[i] = dim->max[i];
      for(i = 0; i < DIM; i++ )
	d.min[i] = dim->min[i];
      
      if( root->c1 != NULL && (virtual_cells || !root->c1->virtual) )
	{
	  d.min[naxis] = (dim->min[naxis] + dim->max[naxis]) * 0.5;
	  ret = iso_find_leafs(root->c1, &d, next_axis(naxis), deep + 1, max_deep, virtual_cells);
	  d.min[naxis] = dim->min[naxis];           
	}
      
      d.max[naxis] = (dim->min[naxis] + dim->max[naxis]) * 0.5;
      if( root->c0 != NULL && (virtual_cells || !root->c0->virtual) )
	{
	  if( !ret )
	    ret = iso_find_leafs(root->c0, &d, next_axis(naxis), deep + 1, max_deep, virtual_cells);
	  else
	    {
	      dyg_node_list* next = ret;
	      while (next->next != NULL)
		next = next->next;
	      next->next = iso_find_leafs(root->c0, &d, next_axis(naxis), deep + 1, max_deep, virtual_cells);
	    }
	}
    }
  return ret;
}