/**
 * This file is adapted by  P. Vincent in 2011
 * from file u_draw_spline.c of xfig (www.xfig.org)
 *
 * FIG : Facility for Interactive Generation of figures
 * Copyright (c) 1985-1988 by Supoj Sutanthavibul
 * Parts Copyright (c) 1989-2002 by Brian V. Smith
 * Parts Copyright (c) 1991 by Paul King
 *
 * Distributed under the GNU General Public Licence
 */


/********************* CURVES FOR SPLINES *****************************

 The following spline drawing routines are from

    "X-splines : A Spline Model Designed for the End User"

    by Carole BLANC and Christophe SCHLICK,
    in Proceedings of SIGGRAPH ' 95

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

#include <cmath.h>

#include "defines.h"
#include "globals.h"
#include "noxprotos.h"
#include "utils.h"

#include "ng_objects.h"


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

#define         MAX_SPLINE_STEP   0.2
/***********************************************************************/

#define Q(s)  (-(s))
#define EQN_NUMERATOR(dim) \
  (A_blend[0]*p0->dim+A_blend[1]*p1->dim+A_blend[2]*p2->dim+A_blend[3]*p3->dim)

static double f_blend (double numerator ,double denominator)
{
  double p = 2 * denominator * denominator;
  double u = numerator / denominator;
  double u2 = u * u;
  return (u * u2 * (10 - p + (2*p - 15)*u + (6 - p)*u2));
}

static double g_blend (double u ,double q)             /* p equals 2 */
{
  return (u*(q + u*(2*q + u*(8 - 12*q + u*(14*q - 11 + u*(4 - 5*q))))));
}

static double h_blend (double u ,double q)
{
  double u2 = u * u;
  return (u * (q + u * (2 * q + u2 * (-2*q - u*q))));
}

static void negative_s1_influence (double t ,double s1 ,double *A0 ,double *A2)
{
  *A0 = h_blend(-t ,Q(s1));
  *A2 = g_blend(t ,Q(s1));
}

static void negative_s2_influence (double t ,double s2 ,double *A1 ,double *A3)
{
  *A1 = g_blend(1-t ,Q(s2));
  *A3 = h_blend(t-1 ,Q(s2));
}

static void positive_s1_influence (int k ,double t ,double s1 ,double *A0 ,double *A2)
{
  double Tk;
  Tk  = k + 1 + s1;
  *A0 = (t+k+1<Tk) ? f_blend(t+k+1-Tk ,k-Tk) : 0.0;
  Tk  = k + 1 - s1;
  *A2 = f_blend(t+k+1-Tk ,k+2-Tk);
}

static void positive_s2_influence (int k ,double t ,double s2 ,double *A1 ,double *A3)
{
  double Tk;
  Tk  = k + 2 + s2;
  *A1 = f_blend(t+k+1-Tk ,k+1-Tk);
  Tk  = k + 2 - s2;
  *A3 = (t+k+1>Tk) ? f_blend(t+k+1-Tk ,k+3-Tk) : 0.0;
}

/******************** for INTERP_XYSPLINE VPoint(s)  ********************/

static int meshmax = 0;
static int j;
#define INCREMENT_ARRAY_SIZE 50

static int point_adding_xy (double *A_blend ,VPoint *p0 , VPoint *p1 , VPoint *p2 , VPoint *p3
			  ,double **pxint ,double **pyint)
 {
   double weights_sum;

   if (j >= meshmax) {
     meshmax += INCREMENT_ARRAY_SIZE;
     *pxint = xrealloc (*pxint ,meshmax*SIZEOF_DOUBLE);
     *pyint = xrealloc (*pyint ,meshmax*SIZEOF_DOUBLE);
     if (*pxint == NULL || *pyint == NULL) return RETURN_FAILURE;
   }

   weights_sum = A_blend[0] + A_blend[1] + A_blend[2] + A_blend[3];
   (*pxint)[j] = EQN_NUMERATOR(x) / (weights_sum);
   (*pyint)[j] = EQN_NUMERATOR(y) / (weights_sum);
   j++;
   return RETURN_SUCCESS;
 }

static int point_adding_vp (double *A_blend ,VPoint *p0 , VPoint *p1 , VPoint *p2 , VPoint *p3
			  ,VPoint **svp)
 {
   double weights_sum;

   if (j >= meshmax) {
     meshmax += INCREMENT_ARRAY_SIZE;
     *svp = xrealloc (*svp ,meshmax*sizeof(VPoint));
     if (*svp == NULL) return RETURN_FAILURE;
   }
   weights_sum = A_blend[0] + A_blend[1] + A_blend[2] + A_blend[3];
   (*svp)[j].x = EQN_NUMERATOR(x) / (weights_sum);
   (*svp)[j].y = EQN_NUMERATOR(y) / (weights_sum);
   j++;
   return RETURN_SUCCESS;
 }

static void point_computing (double *A_blend ,VPoint *p0 ,VPoint *p1 ,VPoint *p2 ,VPoint *p3 ,double *x ,double *y)
{
  double weights_sum;
  weights_sum = A_blend[0] + A_blend[1] + A_blend[2] + A_blend[3];
  *x = EQN_NUMERATOR(x) / (weights_sum);
  *y = EQN_NUMERATOR(y) / (weights_sum);
}

/**
 * This function computes the step used to draw the segment (p1 ,p2)
 * (xv1 ,yv1) : coordinates of the vector from middle to origin
 * (xv2 ,yv2) : coordinates of the vector from middle to extremity
 */
static double step_computing (int k ,VPoint *p0 ,VPoint *p1 ,VPoint *p2 ,VPoint *p3 ,double s1 ,double s2 ,double precision)
{
  double A_blend[4];
  double xstart ,ystart ,xend ,yend ,xmid ,ymid ,xlength ,ylength;
  double start_to_end_dist ,number_of_steps;
  double step ,angle_cos ,scal_prod ,xv1 ,xv2 ,yv1 ,yv2 ,sides_length_prod;


  if ((s1 == 0) && (s2 == 0)) return(1.0);  /* only one step in case of linear segment */

  /* compute coordinates of the origin */
  if (s1>0) {
      if (s2<0) {
	  positive_s1_influence (k ,0.0 ,s1 ,&A_blend[0] ,&A_blend[2]);
	  negative_s2_influence (0.0 ,s2 ,&A_blend[1] ,&A_blend[3]);
      } else {
	  positive_s1_influence (k ,0.0 ,s1 ,&A_blend[0] ,&A_blend[2]);
	  positive_s2_influence (k ,0.0 ,s2 ,&A_blend[1] ,&A_blend[3]);
      }
      point_computing (A_blend ,p0 ,p1 ,p2 ,p3 ,&xstart ,&ystart);
  } else {
      xstart = p1->x;
      ystart = p1->y;
  }

  /* compute coordinates  of the extremity */
  if (s2>0) {
      if (s1<0) {
	  negative_s1_influence (1.0 ,s1 ,&A_blend[0] ,&A_blend[2]);
	  positive_s2_influence (k ,1.0 ,s2 ,&A_blend[1] ,&A_blend[3]);
      } else {
	  positive_s1_influence (k ,1.0 ,s1 ,&A_blend[0] ,&A_blend[2]);
	  positive_s2_influence (k ,1.0 ,s2 ,&A_blend[1] ,&A_blend[3]);
      }
      point_computing(A_blend ,p0 ,p1 ,p2 ,p3 ,&xend ,&yend);
  } else {
      xend = p2->x;
      yend = p2->y;
  }

  /* compute coordinates  of the middle */
  if (s2>0) {
      if (s1<0) {
	  negative_s1_influence (0.5 ,s1 ,&A_blend[0] ,&A_blend[2]);
	  positive_s2_influence (k ,0.5 ,s2 ,&A_blend[1] ,&A_blend[3]);
      } else {
	  positive_s1_influence (k ,0.5 ,s1 ,&A_blend[0] ,&A_blend[2]);
	  positive_s2_influence (k ,0.5 ,s2 ,&A_blend[1] ,&A_blend[3]);
	}
  } else if (s1<0) {
      negative_s1_influence (0.5 ,s1 ,&A_blend[0] ,&A_blend[2]);
      negative_s2_influence (0.5 ,s2 ,&A_blend[1] ,&A_blend[3]);
  } else {
      positive_s1_influence (k ,0.5 ,s1 ,&A_blend[0] ,&A_blend[2]);
      negative_s2_influence (0.5 ,s2 ,&A_blend[1] ,&A_blend[3]);
  }
  point_computing (A_blend ,p0 ,p1 ,p2 ,p3 ,&xmid ,&ymid);

  xv1 = xstart - xmid;
  yv1 = ystart - ymid;
  xv2 = xend - xmid;
  yv2 = yend - ymid;

  scal_prod = xv1*xv2 + yv1*yv2;

  sides_length_prod = sqrt ((xv1*xv1 + yv1*yv1)*(xv2*xv2 + yv2*yv2));

  /* compute cosinus of origin-middle-extremity angle ,which approximates the
     curve of the spline segment */
  if (sides_length_prod == 0.0)
    angle_cos = 0.0;
  else
    angle_cos = scal_prod/sides_length_prod;

  xlength = xend - xstart;
  ylength = yend - ystart;

  start_to_end_dist = sqrt(xlength*xlength + ylength*ylength);

  /* more steps if segment's origin and extremity are remote */
  number_of_steps = sqrt (start_to_end_dist)/2.0;

  /* more steps if the curve is high */
  number_of_steps += (int)((1 + angle_cos)*10);

  if (number_of_steps == 0)
    step = 1;
  else
    step = precision / (double)number_of_steps;

  if ((step > MAX_SPLINE_STEP) || (step == 0))
    step = MAX_SPLINE_STEP;
  return (step);
}

static int spline_segment_computing_xy (double step ,int k ,VPoint *p0 ,VPoint *p1 ,VPoint *p2 ,VPoint *p3
				      ,double s1 ,double s2 ,int last
				      ,double **pxint ,double **pyint)
{
  double A_blend[4];
  double t;
  double tfin = last ? 1.0 + 0.5*step : 1.0;
  if (s1 < 0.0) {
    if (s2 < 0.0) {
      for (t=0.0 ; t < tfin ; t+=step) {
	negative_s1_influence (t ,s1 ,&A_blend[0] ,&A_blend[2]);
	negative_s2_influence (t ,s2 ,&A_blend[1] ,&A_blend[3]);
	if (point_adding_xy (A_blend ,p0 ,p1 ,p2 ,p3 ,pxint ,pyint) == RETURN_FAILURE)
	  return RETURN_FAILURE;
      }
    } else {
      for (t=0.0 ; t < tfin ; t+=step) {
	negative_s1_influence (t ,s1 ,&A_blend[0] ,&A_blend[2]);
	positive_s2_influence (k ,t ,s2 ,&A_blend[1] ,&A_blend[3]);
	if (point_adding_xy (A_blend ,p0 ,p1 ,p2 ,p3  ,pxint ,pyint) == RETURN_FAILURE)
	  return RETURN_FAILURE;
      }
    }
  } else if (s2 < 0.0) {
    for (t=0.0 ; t < tfin ; t+=step) {
      positive_s1_influence (k ,t ,s1 ,&A_blend[0] ,&A_blend[2]);
      negative_s2_influence (t ,s2 ,&A_blend[1] ,&A_blend[3]);
      if (point_adding_xy (A_blend ,p0 ,p1 ,p2 ,p3 ,pxint ,pyint) == RETURN_FAILURE)
	  return RETURN_FAILURE;
    }
  } else {
    for (t=0.0 ; t < tfin ; t+=step) {
      positive_s1_influence (k ,t ,s1 ,&A_blend[0] ,&A_blend[2]);
      positive_s2_influence (k ,t ,s2 ,&A_blend[1] ,&A_blend[3]);
      if (point_adding_xy (A_blend ,p0 ,p1 ,p2 ,p3 ,pxint ,pyint) == RETURN_FAILURE)
	return RETURN_FAILURE;
    }
  }
  return RETURN_SUCCESS;
}

static int spline_segment_computing_vp (double step ,int k ,VPoint *p0 ,VPoint *p1 ,VPoint *p2 ,VPoint *p3
				      ,double s1 ,double s2 ,int last
				      ,VPoint **svp)
{
  double A_blend[4];
  double t;
  double tfin = last ? 1.0 + 0.5*step : 1.0;
  if (s1 < 0.0) {
    if (s2 < 0.0) {
      for (t=0.0 ; t < tfin ; t+=step) {
	negative_s1_influence (t ,s1 ,&A_blend[0] ,&A_blend[2]);
	negative_s2_influence (t ,s2 ,&A_blend[1] ,&A_blend[3]);
	if (point_adding_vp (A_blend ,p0 ,p1 ,p2 ,p3 ,svp) == RETURN_FAILURE)
	  return RETURN_FAILURE;
      }
    } else {
      for (t=0.0 ; t < tfin ; t+=step) {
	negative_s1_influence (t ,s1 ,&A_blend[0] ,&A_blend[2]);
	positive_s2_influence (k ,t ,s2 ,&A_blend[1] ,&A_blend[3]);
	if (point_adding_vp (A_blend ,p0 ,p1 ,p2 ,p3  ,svp) == RETURN_FAILURE)
	  return RETURN_FAILURE;
      }
    }
  } else if (s2 < 0.0) {
    for (t=0.0 ; t < tfin ; t+=step) {
      positive_s1_influence (k ,t ,s1 ,&A_blend[0] ,&A_blend[2]);
      negative_s2_influence (t ,s2 ,&A_blend[1] ,&A_blend[3]);
      if (point_adding_vp (A_blend ,p0 ,p1 ,p2 ,p3 ,svp) == RETURN_FAILURE)
	  return RETURN_FAILURE;
    }
  } else {
    for (t=0.0 ; t < tfin ; t+=step) {
      positive_s1_influence (k ,t ,s1 ,&A_blend[0] ,&A_blend[2]);
      positive_s2_influence (k ,t ,s2 ,&A_blend[1] ,&A_blend[3]);
      if (point_adding_vp (A_blend ,p0 ,p1 ,p2 ,p3 ,svp) == RETURN_FAILURE)
	return RETURN_FAILURE;
    }
  }
  return RETURN_SUCCESS;
}


/****************************** MAIN METHODS ******************************/

/**
 *  X-spline interpolation for sets
 *  input:
 *     xint, yint, meshlen : the allocated set to receive the interpolation
 *               x ,y ,len : the dataset to be interpolated
 *                       s : THe "s" paraemter of X-Splines
 *  output:
 *    xint ,yint ,n_interp : the computed arrays with their actual length
 *
 *  Compare to aspline, seval, interpolate and do_interp
 */
int xspline_eval_xy (double **pxint ,double **pyint ,int meshlen
		     ,double *x ,double *y ,int len ,int *n_interp
		     ,double s ,double pr)
{
  int k ,last = FALSE;
  VPoint p0 ,p1 ,p2 ,p3;
  double step;
  meshmax = meshlen;
  j = 0;
  for (k = 0; k < len - 1; k++) {
    if (k == 0) {
      p0.x = p1.x = x[0]; p2.x = x[1]; p3.x = x[2];
      p0.y = p1.y = y[0]; p2.y = y[1]; p3.y = y[2];
      last = FALSE;
    } else if (k < len - 2) {
      p0.x = x[k-1]; p1.x = x[k]; p2.x = x[k+1]; p3.x = x[k+2];
      p0.y = y[k-1]; p1.y = y[k]; p2.y = y[k+1]; p3.y = y[k+2];
      last = FALSE;
    } else if (k == len - 2) {
      p0.x = x[k-1]; p1.x = x[k]; p2.x = p3.x = x[k+1];
      p0.y = y[k-1]; p1.y = y[k]; p2.y = p3.y = y[k+1];
      last = TRUE;
    }
    step = step_computing (k ,&p0 ,&p1 ,&p2 ,&p3 ,s ,s ,/* precision */ pr);
    if (spline_segment_computing_xy (step ,k ,&p0 ,&p1 ,&p2 ,&p3 ,s ,s ,last ,pxint ,pyint) == RETURN_FAILURE) {
      errmsg ("xspline_eval_xy: error in allocation");
      return RETURN_FAILURE;
    }
  }
  *n_interp = j - 1;
  return RETURN_SUCCESS;
}

int xspline_eval_vp (VPoint **svp ,VPoint *vp ,int len ,int *n_interp ,double s ,double pr ,int closed)
{
  int k ,last = FALSE;
  VPoint p0 ,p1 ,p2 ,p3;
  double step = 0.1;
  meshmax = len;
  j = 0;
  for (k = 0; k < len - 1; k++) {
    if (k == 0 && closed) {
      p0 = vp[len-1]; p1 = vp[0]; p2 = vp[1];  p3 = vp[2];
    } else if (k == 0) {
      p0 = p1 = vp[0]; 		  p2 = vp[1];  p3 = vp[2];
      last = FALSE;
    } else if (k < len - 2) {
      p0 = vp[k-1]; p1 = vp[k]; p2 = vp[k+1]; p3 = vp[k+2];
      last = FALSE;
    } else if (k == len - 2 && closed) {
      p0 = vp[k-1]; p1 = vp[k]; p2 =vp[k+1];  p3 = vp[0];
    } else if (k == len - 2) {
      p0 = vp[k-1]; p1 = vp[k]; p2 = p3 = vp[k+1];
      last = TRUE;
    }
    step = step_computing (k ,&p0 ,&p1 ,&p2 ,&p3 ,s ,s ,/* precision */ pr);
    if (spline_segment_computing_vp (step ,k ,&p0 ,&p1 ,&p2 ,&p3 ,s ,s ,last ,svp) == RETURN_FAILURE) {
      errmsg ("xspline_eval_xy: error in allocation");
      return RETURN_FAILURE;
    }
  }
  if (closed) {
    p0 = vp[len-2]; p1 = vp[len-1]; p2 = vp[0]; p3 = vp[1];
    if (spline_segment_computing_vp (step ,len-1 ,&p0 ,&p1 ,&p2 ,&p3 ,s ,s ,last ,svp) == RETURN_FAILURE) {
      errmsg ("xspline_eval_xy: error in allocation");
      return RETURN_FAILURE;
    }
  }
  *n_interp = j - 1;
  return RETURN_SUCCESS;
}

/****************************** METHODS for  INTERP_XSPLINE : y = f(x) ******************************/

static double point_adding_y (double *A_blend ,double y0 , double y1 , double y2 , double y3)
 {
   double weights_sum ,numerator ,yint;
   weights_sum =  A_blend[0]    + A_blend[1]    + A_blend[2]    + A_blend[3];
   numerator   = (A_blend[0]*y0 + A_blend[1]*y1 + A_blend[2]*y2 + A_blend[3]*y3);
   yint =  numerator / weights_sum;
   return yint;
 }


/**   ulen
 *  X-spline interpolation for sets
 *  input:  mesh ,meshlen ,x ,y ,len
 *  output: yint
 *  Compare to aspline, seval, interpolate and do_interp
 */
void xspline_eval_y (double *mesh ,double *yint ,int meshlen ,double *x ,double *y ,int len ,double xspline_s)
{
  int j ,k ,kk ,m;
  double y0 ,y1 ,y2 ,y3;
  double t ,A_blend[4];
  m = monotonicity (x ,len ,FALSE);
  y0 = y1 = y2 = y3 = y[0];
  kk = -1;
  for (j = 0; j < meshlen; j++) {
    k = find_span_index (x ,len ,m ,mesh[j]);
    if (k < 0) {
      k = 0;
      t = (mesh[j] - x[0]) / (x[1] - x[0]);
    } else if (k > len-2) {
      k = len - 1;
      t = (mesh[j] - x[len-2]) / (x[len-1] - x[len-2]);
    } else {
      t = (mesh[j] - x[k]) / (x[k+1] - x[k]);
    }
    if (k != kk) {
      /* Control points have chenged */
      if (k == 0) {
	y0 = y1 = y[0];         y2 = y[1];   y3 = y[2];
      } else if (k < len - 2) {
	y0 = y[k-1]; y1 = y[k]; y2 = y[k+1]; y3 = y[k+2];
      } else if (k == len - 2) {
	y0 = y[k-1]; y1 = y[k]; y2 = y[k+1]; y3 = y[k+1];
      }
    }
    if (xspline_s < 0.0) {
      negative_s1_influence (t ,xspline_s ,&A_blend[0] ,&A_blend[2]);
      negative_s2_influence (t ,xspline_s ,&A_blend[1] ,&A_blend[3]);
    } else {
      positive_s1_influence (k ,t ,xspline_s ,&A_blend[0] ,&A_blend[2]);
      positive_s2_influence (k ,t ,xspline_s ,&A_blend[1] ,&A_blend[3]);
    }
    yint[j] = point_adding_y (A_blend ,y0 ,y1 ,y2 ,y3);
  }

}
