/*
 * smooth.c
 *
 * smooth(double *invals, double *outvals, int size, int window, double max_gradient, double min_gradient);
 *
 * Smooth the values in the array invals, writing out the smoothed values to array outvals, of size size,
 * by using a sliding window average with a window size of window. Further limit the average to have maximal 
 * slope of max_gradient +ve and min_gradient -ve
 *
 */


/*
 * gradient filter 
 *
 */


static int g_filter ;
static double l_t ;
static double l_v ;
static double l_min ;
static double l_max ;
static int gcount ;

static void
init_gradient_filter(double min, double max)
{
  l_t = 0 ;
  l_v = 0.0 ;
  g_filter =  !((min == 0.0) && (max == 0.0)) ;
  l_min = min ;
  l_max = max ;
  gcount = 0 ;
}


static double
gradient_filter(unsigned long t, double v, double *ming, double *maxg)
{
  double slope ;


  if (!l_t) {
    l_v = v ;
    l_t = t ;
    }
  else if (l_t < t) {
    slope = (v - l_v) / (double) (t - l_t) ;
    if (!gcount++) {
      *ming = slope ;
      *maxg = slope ;
      }
    if (slope < *ming) *ming = slope ;
    if (slope > *maxg) *maxg = slope ;
    //    printf("slope = %5.3f min=%5.3f Max = %5.3f\n",slope,*ming,*maxg) ;

    if ((l_max != 0.0) && (slope > l_max)) slope = l_max;
    if ((l_min != 0.0) && (slope < l_min)) slope = l_min;
    l_v +=  (slope * (t - l_t)); 
    l_t = t ;
    }
  if (!g_filter) return(v) ;  
  return(l_v) ;
}


/*
 * entry point: smooth
 *
 */

void
smooth(double *invals, double *outvals, int size, int window, 
       double min_gradient, double max_gradient, double *ming, double *maxg)
{
  int index = 0;
  int lhs, rhs ;
  int w ;
  double sum ;
  int count ;
  int startindex = -1 ;
  int swin ;

  init_gradient_filter(min_gradient,max_gradient) ;

  for (index = 0 ; index < size ; ++index) {
    swin = window ;
    if ((startindex < 0) && (invals[index] == 0.0))
      outvals[index] = 0.0 ;
    else {
      if (startindex < 0) startindex = index ;
      count = 0 ;
      sum = 0.0 ;
      // now scale the window back
      // if lhs < 0 then window += (2 *lhs)
      // if rhs >= size then window += (2 * (size - rhs - 1))
      if ((index - (swin / 2)) < 0) swin += (2 * (index - (swin / 2))) ;
      if ((index + (swin / 2)) >= size) swin -= (2 * (index + (swin / 2) - size  + 1)) ;
      lhs = index - (swin / 2) ;
      rhs = index + (swin / 2) ;

      for (w = lhs ; w <= rhs ; ++w) {
        if ((w >= startindex) && (w < size)) {
          sum += invals[w] ;
          ++count ;
          }
        }
      outvals[index]  = sum / (double) count ;
      outvals[index] = gradient_filter(index +1, outvals[index],ming,maxg) ;
      }
    }
}