// gFilters.h

#ifndef gFiltersh
#define gFiltersh

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include "../../Global.h"

struct gFilters {

//>>>>>>>>>>>>>> Periodic methods which return a parameter with interval [0,1] <<<<<<<<<<<<<<<<
// x is any value. Period is in the same units as x. return values are in the Interval [0,1]

  static double SawTooth(const double& x, const double& Period) {return Mod(x,Period)/Period;}
  static double Triangle(const double& x, const double& Period) {return fabs(Signed(SawTooth(x,Period)));} // High at zero. To start at 0,0 use 1-Triangle
  static bool   Square  (const double& x, const double& Period) {return Step(SawTooth(x,Period), 0.5);}
  static double Stripes (const double& x, const double& Period) {return HermiteStep(Triangle(x,Period), 0.4, 0.6);} // Smoothed Square
  static double Repeat  (const double& x, const double& Count, bool Mirror=false) {return Mirror ? 1-Triangle(x,1/Count) : SawTooth(x,1/Count);}
// 2 Dimensional:
  static bool   Chess   (const double& x, const double& y, int xCount, int yCount) {return Square(x,1/xCount)==Square(y,1/yCount);} // 0,0 returns true

//>>>>>>>>>>>>>> Step methods which return a parameter with interval [0,1] <<<<<<<<<<<<<<<<
// x is any value, a and b are values that x is likely to pass through, in the same units as x.
// return values are in the Interval [0,1]

  // return is 0 when x<=a and 1 when x>a.
  static bool Step(const double& x, const double& a) {return x>a;}

  // return is 0 when x<=a;  1 when x>b and unchanged in the [a,b] interval:
  static double Threshold(const double& x, const double& a, const double& b) {return x<a ? 0 : (x<b ? x : 1);}

  // returns 0 when x<a, rising linearly to 1 when x>b, holding at 1 thereafter.
  static double Ramp(const double& x, const double& a, const double& b) {
    if(a==b) return 0;
    if(x<a) return 0;
    if(b<x) return 1;
    return Parameterize(x,a,b);
  }
  // returns 0 when x<a, rising to 1 at b with a Hermite curve, holding at 1 thereafter.
  static double HermiteStep(const double& x, const double& a, const double& b) {
    if(a==b) return 0;
    if(x<a) return 0;
    if(b<x) return 1;
    return HermiteCurve(Parameterize(x,a,b)); // 3t^2-2t^3
  }

// Filter an edge transition, Tx, with CatStep(t, Tx-(2.5*Zoom), Tx+(1.5*Zoom)); Large Zoom will blur.
// Results are in the interval [-0.0416667, 1.0416667] which is [-1/24.0, 1/24.0]
  static double CatStep(const double& x, const double& a, const double& b) {
    if(a==b) return 0;
    if(x<a) return 0;
    if(b<x) return 1;
    return CatStepCurve(Parameterize(x,a,b));
  }

//>>>>>>>>>>>>>> Zero-crossing methods which alter a parameter with interval [0,1] <<<<<<<<<<<<<<<<

  static double HermiteCurve(const double& t) {return t*t*(3-2*t);} // 3t²-2t³ Called "Ease" in some art packages

  static double RootCurve(const double& t) {return t==0 ? 0 : sqrt(fabs(t));}

  static double SinusoidalCurve(const double& t) {return sin(M_PI_2*t);}

  static double CircleCurve(const double& t) {return sqrt(fabs(t*(2-t)));} // 1-(1-t)*(1-t)

// Full Gaussian 1D equation is: exp(-(t*t/(2*s*s))/(sqrt(2*M_PI)*s); where s is the standard deviation. This uses s=1/3.0; to make the bell peak at 1.0
  static double GaussianCurve(const double& t) {return exp(-3.125*t*t)*1.0026513098524002009663061139244;} // The number is 1/(sqrt(2*M_PI)/3).
// Full Gaussian 2D equation is: exp(-(x*x+y*y)/(2*s*s))/((2*M_PI)*s*s); where s is the standard deviation.

  static double HighContrastCurve(double t) {
              t=Bias(  t, 0.67);
    if(t<0.5) t=Bias(2*t, 0.14)/2;
    return      Bias(  t, 0.35);
  }

  static double CatStepCurve(const double& t) {
    if(t<0.25) return 32*t*t*t*(-1/3.0 + t);
    if(t<0.5)  return t*(t*(t*(-t*96+416/3.0)- 64)+12)- 5/6.0; // Horner's Rule
    if(t<0.75) return t*(t*(t*( t*96-736/3.0)+224)-84)+67/6.0;
    return            t*(t*(t*(-t*32+352/3.0)-160)+96)-61/3.0;
  }

//>>>>>>>>>>>>>> Methods which alter a parameter with interval [0,1] <<<<<<<<<<<<<<<<

  // Gain alters the slope at the mid-point of the interval. Steeper slopes come from higher values:
  static double Gain(const double& t, const double& Gain) {
    return t<0.5 ?   Bias(2*   t , 1-Gain)/2
                 : 1-Bias(2*(1-t), 1-Gain)/2;
  }

  // Bias pulls the interval; <0.5 keeps parameter lower longer
  static double Bias(const double& t, const double& Bias) {
         if(   t<0.001) return 0;
    else if(   t>0.999) return 1;
    else if(Bias<0.001) return 0;
    else if(Bias>0.999) return 1;              // The long number is  1/ln(0.5); [ln in c++ is log; log in c++ is log10]
    else                return pow(t, log(Bias) * -1.4426950408889634073599246810019);
  }

//>>>>>>>>>>>>>> Symmetrical Filters changing a parameter with interval [-1,1] to [0,1] <<<<<<<<<<<<<<<<

// Deals with t in the range [-1,1] (FilterWidth=1)
// Results are in the interval [0,1]
  static double BoxFilter(const double& t) {return ((t>-0.5) && (t<=0.5) ? 1 : 0);}

// Deals with t in the range [-1,1] (FilterWidth=1)
// Results are in the interval [0,1]
  static double TriangleFilter(double t) {
    if(t<0) t=-t;
    if(t<1) return 1-t;
    return 0;
  }

// Deals with t in the range [-1,1] (FilterWidth=1)
// Results are in the interval [0,1]
  static double HermiteFilter(double t) {
    if(t<0) t=-t;
    return (t>=1 ? 0 : 1-HermiteCurve(t));
  }

//>>>>>>> Methods which take (potentially over-sized) signed intervals and return Intervals in the range [0-1] (some return small negatives) <<<<<<<<<<<<

// Deals with t in the range [-1.5,1.5] (FilterWidth=1.5)
// Results are in the interval [0, 0.75]
  static double BellFilter(double t) {
    if(t<0) t=-t;
    if(t<0.5) return 0.75-(t*t);
    if(t<1.5) {
      t-=1.5;
      return 0.5*(t*t);
    }
    return 0;
  }

// Deals with t in the range [-2,2] (FilterWidth=2)
// Results are in the interval [0, 0.6666]
  static double bSplineFilter(double t) {
    if(t<0) t=-t;
    if(t<1) {
      double t2=t*t;
      return 0.5*t2*t-t2+2/3.0;
    }
    if(t<2) {
      t=2-t;
      return 1/6.0*t*t*t;
    }
    return 0;
  }

// Deals with t in the range [-2,2] (FilterWidth=2)
// Results are in the interval [-0.074052, 1]
  static double CatRomFilter(double t) {
    double t2=t*t;
    if(t<0) t=-t;
    if(t<1) return  1.5*t2*t - 2.5*t2       + 1;
    if(t<2) return -0.5*t2*t + 2.5*t2 - 4*t + 2;
    return 0;
  }

// Deals with t in the range [-3,3] (FilterWidth=3)
// Results are in the interval [-0.147233, 1]
  static double LanczosFilter(double t, double FilterWidth=3) {
    if(t<0) t=-t;
    if(t<FilterWidth) return sinc(t) * sinc(t/FilterWidth);
    return 0;
  }

// Deals with t in the range [-2,2] (FilterWidth=2)
// Results are in the interval [-0.0362564, 0.8888]
  static double MitchellFilter(const double& t) {
    double a=1/3.0;
    double b=(  6- 2*a);
    double c=(-18+18*a);
    double d=( 12-15*a);
    double e=(    32*a);
    double f=(   -60*a);
    double g=(    36*a);
    double h=(   - 7*a);
    if (t<-2) return 0;
    if (t<-1) return (-t*(-t*(-t*h+g)+f)+e)*1/6.0; // Horner's Rule
    if (t< 0) return ( t*  t*(-t*d+c)+b)   *1/6.0;
    if (t< 1) return ( t*  t*( t*d+c)+b)   *1/6.0;
    if (t< 2) return ( t*( t*( t*h+g)+f)+e)*1/6.0;
    return 0;
  }

//>>>>>>>>>>>>>>>>>>>>>>>> return any value from an interval [0-1] <<<<<<<<<<<<<<<<<<<<<<<<<<<<

/* Compare the Pixel distance (P) to the Real distance (R) (just give squares: don't use sqrt).
For simple antialias, Blur=3.5*r; where R=r*r;
More generally, if int BlurLevel is between 1 and 5, use double Blur=(1.5*BlurLevel+2)*r;
Returns a Transparency Value (0 is opaque, 1 is transparent) */
  static double Antialias(const double& P, const double& R, const double& Blur) {return (P>R-1) ? 1 : ((P>R-Blur) ? (P-R)/Blur : 0);}

  static double Random(const double& a=0, const double& b=1) {return InterpolateLinear((1402024253UL*rand()+586950981)/double(ULONG_MAX), a,b);}
  static double ParameterizeLinear(const double& x, const double& a, const double& b) {return Parameterize(x,a,b);}
  static double  InterpolateLinear(const double& t, const double& a, const double& b) {return Interpolate (t,a,b);}
  // InterpolateCosine is slow! Use:  InterpolateHermite(t,a,a,b,a); it gives a very similar result!
  static double InterpolateCosine(const double& t, const double& a, const double& b) {return InterpolateLinear((1-cos(t*M_PI))/2, a,b);}

  // Uses four values to find a value between the two MIDDLE ones.
  static double InterpolateBezier(const double& t, const double& a, const double& b, const double& c, const double& d) {return InterpolateBezierStd(t, b, b+(b-a), c-(d-c), c);}
  // Uses four values to find a value between the two END ones (b and c indicate initial direction for a and d).
  static double InterpolateBezierStd(const double& t, const double& a, const double& b, const double& c, const double& d) {
    return t*(t*(t*(-a+3*(b-c)+d)+3*(a-b-b+c))+3*(b-a))+a;
  }
  // Uses four values to find a value between the two MIDDLE ones (passes through b and c):
  // If you don't have four points, a=Interpolate(-1, b,c), d=Interpolate(2, b,c);
  static double InterpolateHermite(const double& t, const double& a, const double& b, const double& c, const double& d) {
    return t*(t*(t*(a+2*(b-c)+d)-(2*a+3*(b-c)+d))+a)+b;
  }//[(2t^3 - 3t^2 +1) (-2t^3 + 3t^2) (t^3 - 2t^2 + t) (t^3 - t^2)]

  // Tension: 0 is normal; -1 is low; 1 is high.
  // Bias: 0 is even, positive is towards first segment, negative towards the other.
  static double InterpolateHermite(const double& t, const double& a,const double& b,const double& c,const double& d, const double& Tension, const double& Bias) {
    double t2=t*t;
    double t3=t2*t;
    double m0=(b-a)*(1+Bias)*(1-Tension)
             +(c-b)*(1-Bias)*(1-Tension);
    double m1=(c-b)*(1+Bias)*(1-Tension)
             +(d-c)*(1-Bias)*(1-Tension);
    double a0= 2*t3 - 3*t2 + 1;
    double a1=   t3 - 2*t2 + t;
    double a2=   t3 -   t2;
    double a3=-2*t3 + 3*t2;
    return a0*b+a3*c+(a1*m0+a2*m1)/2;
  }
  // Uses four values to find a value between the two MIDDLE ones (passes through b and c):
  // If you don't have four points, a=Interpolate(-1, b,c), d=Interpolate(2, b,c);
  static double InterpolateCatmullRom(const double& t, const double& a, const double& b, const double& c, const double& d) {
    if(t==0) return b;
//  if((t==0)||(b==c)) return b; Clips it into the Hermite Curve
    if(t==1) return c;
    return 0.5*(t*(t*(t*((d-a)+3*(b-c))+(2*a - 5*b + 4*c -d))-a+c))+b;
  }
  // Uses four values to find a value between the two MIDDLE ones (passes through b and c):
  // If you don't have four points, a=Interpolate(-1, b,c), d=Interpolate(2, b,c);
  static double InterpolateCubic(const double& t, const double& a, const double& b, const double& c, const double& d) {
    return t*(t*(t*(-a+b-c+d)+(2*(a-b)+c-d))+(c-a))+b;
  }
  // Uses four values to find a value between the two MIDDLE ones (passes through b and c):
  // If you don't have four points, a=Interpolate(-1, b,c), d=Interpolate(2, b,c);
  static double InterpolateBSpline(const double& t, const double& a, const double& b, const double& c, const double& d) {
    return 1/6.0*(t*(t*(t*(-a+3*b-3*c+d)+(3*a-6*b+3*c))+3*(c-a))+(a+4*b+c));
  }
};

class CatmullRom1D { // Cubic Spline. Uses four values to find a value between the two middle ones:
  double C[4];
public:
  CatmullRom1D(double v[4]) {
    C[0]=(              v[1]                    );
    C[1]=(-v[0]/2            +     v[2]/2       );
    C[2]=( v[0]   - 2.5*v[1] +   2*v[2] - v[3]/2);
    C[3]=(-v[0]/2 + 1.5*v[1] - 1.5*v[2] + v[3]/2);
  }
  CatmullRom1D(const double& a, const double& b, const double& c,const double& d) {
    C[0]=(           b              );
    C[1]=(-a/2         +     c/2    );
    C[2]=( a   - 2.5*b +   2*c - d/2);
    C[3]=(-a/2 + 1.5*b - 1.5*c + d/2);
  }
  double Interpolate(const double& t) const {return t*(t*(t*C[3]+C[2])+C[1])+C[0];}
};

class CatmullRom2D { // Bicubic Spline. Catmull-Rom interpolation of 4x4 vertex mesh's middle square:
  double C[4][4]; // Coefficients matrix
public:
  CatmullRom2D(double Grid[4][4]) {
    static const double M[4][4]={ // Catmull-Rom basis matrix
      {-0.5,  1.5, -1.5, 0.5},
      { 1  , -2.5,  2  ,-0.5},
      {-0.5,  0  ,  0.5, 0  },
      { 0  ,  1  ,  0  , 0  }
    };
    double T[4][4];
    for(int i=0; i<4; ++i) // T = G * MT
    for(int j=0; j<4; ++j)
    for(int k=0; k<4; ++k) T[i][j] += Grid[i][k] * M[j][k];
    for(int i=0; i<4; ++i) // C = M * T
    for(int j=0; j<4; ++j)
    for(int k=0; k<4; ++k) C[i][j] += M[i][k] * T[k][j];
  }
  double Interpolate(const double& s, const double& t) const {
    const double* X3=C[0];
    const double* X2=C[1];
    const double* X1=C[2];
    const double* X0=C[3];
    return s*(s*(s*(t*(t*(t*X3[0]+X3[1])+X3[2])+X3[3]) // Horner's Rule
                 + (t*(t*(t*X2[0]+X2[1])+X2[2])+X2[3]))
              +    (t*(t*(t*X1[0]+X1[1])+X1[2])+X1[3]))
           +       (t*(t*(t*X0[0]+X0[1])+X0[2])+X0[3]);
  }
};

#endif // ndef gFiltersh