#include <math.h>

#include "vector2d.h"

#include "Common.h"
#include "Kernel.h"


Kernel::Kernel()
{
  h = h2 = h2inv = 1.0;
}


Kernel* Kernel::newKernel(DEKernelShape ks)
{
  switch (ks)
    {
    case DE_CYLINDER_KERNEL:
      return new CylinderKernel();

    case DE_CONE_KERNEL:
      return new ConeKernel();

    case DE_K2_KERNEL:
      return new K2Kernel();

    case DE_EPANECHNIKOV_KERNEL:
      return new EpanechnikovKernel();

    default:
      break;
    }

  return NULL;
}


void Kernel::setSize(const float newh)
{
  // Check precondition when debugging.
  DebugAssert(h>0.0f,"Kernel size is a stricly positive number");

  h     = newh;
  h2    = h*h;
  h2inv = 1.0f / h2;
}


int Kernel::isInside(const VEC2D &point,const VEC2D &center)
{
  VEC2D  d;
  float  f;

  // Use square norm, faster...
  VEC2DSUBTRACT(point,center,d);
  f = VEC2DNORM2(d);
  
  return (f<h2);
}


int Kernel::intersectSegment(const VEC2D &lorg,const VEC2D &ldir,const float ls,const VEC2D &cc)
{
  VEC2D  eo;
  float  v,disc,l,d,i0,i1;
  
  VEC2DSUBTRACT(cc,lorg,eo);
  v = VEC2DDOTPRODUCT(eo,ldir);
  l = VEC2DDOTPRODUCT(eo,eo);
  disc = h2 - (l - (v*v));
  
  // How much intersections ?
  if (disc > EPSILON)
    {
      // The line intersects the kernel.
      d = sqrt(disc);
      i0 = v-d;
      i1 = v+d;

      // Check if intersection is on line segment.
      if (((i0>0.0f)&&(i0<ls))||((i1>0.0f)&&(i1<ls)))
	return 1;
    }
     
  return 0;
}



float CylinderKernel::evaluate(const VEC2D &point,const VEC2D &center)
{
  VEC2D aux;
  float tp;
  
  // Find distance
  VEC2DSUBTRACT(point,center,aux);
  tp = VEC2DNORM2(aux);
  if (tp<h2)
    {
      // Point inside kernel
      return h2inv*M_1_PI;
    }
  else
    return 0.0f;
}


float ConeKernel::evaluate(const VEC2D &point,const VEC2D &center)
{
  VEC2D aux;
  float tp;
  
  // Find distance
  VEC2DSUBTRACT(point,center,aux);
  tp = VEC2DNORM2(aux);
  if (tp<h2)
    {
      // Point inside kernel
      return 3.0f*M_1_PI*h2inv*(1.0f-sqrt(tp*h2inv));
    }
  else
    return 0.0f;
}


float K2Kernel::evaluate(const VEC2D &point,const VEC2D &center)
{
  VEC2D aux;
  float tp;
  
  // Find distance
  VEC2DSUBTRACT(point,center,aux);
  tp = VEC2DNORM2(aux);
  if (tp<h2)
    {
      // Point inside kernel
      tp = 1.0f-(tp*h2inv);
      return 3.0f*M_1_PI*tp*tp*h2inv;
    }
  else
    return 0.0f;
}


float EpanechnikovKernel::evaluate(const VEC2D &point,const VEC2D &center)
{
  VEC2D aux;
  float tp;
  
  // Find distance
  VEC2DSUBTRACT(point,center,aux);
  tp = VEC2DNORM2(aux);
  if (tp<h2)
    {
      // Point inside kernel
      tp = (1.0f-(tp*h2inv));
      tp = M_2_PI*tp*h2inv;
      return tp;
    }
  else
    return 0.0f;
}