/****************************************************************************
Copyright (C) 2012 HidraVFX development team
            <https://launchpad.net/~hidravfx-dev-team>

This file is part of the HidraVFX project.

HidraVFX is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

HidraVFX is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with HidraVFX. If not, see <http://www.gnu.org/licenses/>.
****************************************************************************/

#include <math.h>
#include <mathplus.h>

enum tdistancetype {dteuclidian, dtfasteuclidian, dtsqreuclidian, dtmanhattan,
                    dtcamberra, dtminkovsky,dtmaxnorma,dtxmanhattan,dtokta};

double valami(double dx, double dy)
{
   double d,fi;

  d = sqrt(sqr(dx)+sqr(dy));
  fi = atan2(dy,dx);
  return (d*(3.0+cos(fi*d)) );
}

double cubics(double dx, double dy)
{
  double d,fi;

  d = sqrt(sqr(dx) + sqr(dy));
  fi = atan2(dy,dx);
  return ( d*(1.0+0.5*cos(3*fi)*cos(3*fi)*cos(3*fi)) );
}

double heart(double dx, double dy)
{ double d,fi;

  d = sqrt(sqr(dx) + sqr(dy));
  fi = atan2(dy,dx);
  if( d == 0.0 ) return 0.0 ;
  else return (d*(1.0+0.3*cos(1.2*fi/d)));
}

double tentacles(double dx, double dy)
{  double d,fi1,fi2,fi;

   d = sqrt(sqr(dx)+sqr(dy));
   fi1 = atan2(dy,dx);
   fi2 = 4*atan2(fabs(dy),fabs(dx));
   fi = (2*fi1+2*fi2);
   return ( d*(1.0+0.5*cos(4.0*fi+4.0*cos(5.0*d))));
}

double clover(double dx, double dy)
{  double d,fi;

   d = sqrt(sqr(dx)+sqr(dy));
   fi = atan2(dy,dx);
   fi = 4.0*fi;
   return(d*(3.0+cos(d)*sin(fi))/4.0);
}

double zigzag(double dx, double dy)
{  double d,fi;

   d = sqrt(sqr(dx)+sqr(dy));
   fi = atan2(dy,dx);
   return ( d*(1.0+fabs(cos(2.0*fi)))/2.0 );
}

double flower(double dx, double dy)
{ double d,fi;

   d = sqrt(sqr(dx)+sqr(dy));
   fi = 4.0*atan2(fabs(dy),fabs(dx));
   return( d*(1.0+exp(exp(1.0-fabs(cos(fi)))))/7.0);
}

double flower2(double dx,double dy)
{ double d,fi;

   d = sqrt(sqr(dx)+sqr(dy));
   fi = 4.0*atan2(fabs(dy),fabs(dx));
   return( d*(1.0+0.5*cos(3.0*fi)) );
}

double spiral(double dx,double dy)
{ double d,fi;

  d = sqrt(sqr(dx)+sqr(dy));
  fi = atan2(dy,dx);
  return ( d*(1.0+0.5*cos(3.0*fi)));
}

double spiral2(double dx, double dy)
{ double d,fi;

  d = sqrt(sqr(dx)+sqr(dy));
  fi = 4.0*atan2(fabs(dy),fabs(dx));
  return(d*(1.0+0.5*cos(3.0*fi+3.0*cos(5.0*d))));
}

double euclidian(double dx, double dy)
{
   return(sqrt(dx*dx+dy*dy));
}

double baroque(double dx, double dy)
{ double d,fi;

  d = sqrt(sqr(dx)+sqr(dy));
  fi = 4.0*atan2(fabs(dy),fabs(dx));
  return(sin(atan2(d,fi)));
}

double sqreuclidian(double dx, double dy)
{
  return(dx*dx+dy*dy);
}

double manhattan(double dx,double dy)
{
   return(fabs(dx)+fabs(dy));
}

double maxnorma(double dx, double dy)
{
   return(max(fabs(dx),fabs(dy)));
}

double xmanhattan(double dx, double dy)
{ double result;

   result = fabs(dx)*fabs(dy)*500.0+fabs(dx)+fabs(dy);
   if( result > 0.0 ) {  result = log(sqr(result))
   ; }else { result = 0.0; }
   return (result);
}

double okta(double dx, double dy)
{ double result;
   result=max(1.5*(max(fabs(dx),fabs(dy))),(fabs(dx)+fabs(dy)));
   return(result);
}

double circle(double dx, double dy)
{ double result;

   result = sqrt(dx*dx+dy*dy);
   if( result>1.0 ) {  result = result-1.0 ; }else { result = 1.0-result; }
   return(result);
}

double correlation(double x1, double x2, double y1, double y2)
{  double k,a1,a2, result;

   a1 = (x1+y1)/2.0;
   a2 = (x2+y2)/2.0;
   k = sqrt((sqr(x1-a1)+sqr(y1-a1))*(sqr(x2-a2)+sqr(y2-a2)));
   if ( k!=0 ) {  result=((x1-a1)*(x2-a2)+(y1-a1)*(y2-a2))/k ; }else { result=0.0; }
   return(result);
}

double kendallsCorrelation(double x1, double x2,double y1,double y2)
{
   return(1.0-sign(y1-x1)*sign(y2-x2));
}

double chisquare(double x1, double x2, double y1, double y2)
{ double result;

   result = sqr(x1/(x1+y1)-x2/(x2+y2))/(x1+x2)
          +sqr(y1/(x1+y1)-y2/(x2+y2))/(y1+y2);
   return(result);
}

double cosin(double x1, double x2, double y1, double y2)
{ double result;

   result = (x1*x2+y1*y2)/( sqrt(x1*x1+y1*y1)*sqrt(x2*x2+y2*y2) );
   return(result);
}

double circleDist(double x1, double x2, double y1, double y2)
{ double result;

  result=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
  if( result>1.0 ) {  result=result-1.0 ; }else { result=1.0-result; }
  return(result);
}

double oktaDist(double x1, double x2, double y1, double y2)
{ double result;

  result = max(1.5*(max(fabs(x1-x2),fabs(y1-y2))),(fabs(x1-x2)+fabs(y1-y2)));
  return(result);
}

double euclidianDist(double x1, double x2, double y1, double y2)
{
  return((sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))));
}

double fasteuclidianDist(double x1, double x2, double y1, double y2)
{
  return( fabs(x1-x2)+0.5*fabs(y1-y2) );
}

double sqreuclidianDist(double x1, double x2, double y1, double y2)
{
  return(((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)));
}

double manhattanDist(double x1, double x2, double y1,double y2)
{
    return(fabs(x1-x2)+fabs(y1-y2));
}

double camberraDist(double x1, double x2, double y1, double y2)
{
    return(fabs(x1-x2)/fabs(x1+x2)+fabs(y1-y2)/fabs(y1+y2));
}

double minkovskyDist(double x1,double x2,double y1,double y2, int n)
{ double result;

  result = sqrt(pow(pow(fabs(x1-x2),n)+pow(fabs(y1-y2),n),1.0/n));
  return(result);
}

double maxnormaDist(double x1,double x2,double y1, double y2)
{
    return( max(fabs(x1-x2),fabs(y1-y2)) );
}

double xmanhattanDist(double x1, double x2, double y1, double y2)
{ double result;

  result = (fabs(x1-x2))*(fabs(y1-y2))*500.0+(fabs(x1-x2))+(fabs(y1-y2));
  if( result>0.0 ) {  result = log(sqr(result)); } else { result=0.0; }
  return(result);
}

double euclidian3Dist(double x1,double x2,double y1,double y2,double z1,double z2)
{
    return(sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)));
}

double manhattan3Dist(double x1,double x2,double y1,double y2,double z1,double z2)
{
    return( fabs(x1-x2)+fabs(y1-y2)+fabs(z1-z2) );
}

double minkovsky3Dist(double x1,double x2,double y1,double y2,double z1,double z2,double n)
{ double result;

    result = sqrt(pow(pow(fabs(x1-x2),n)+pow(fabs(y1-y2)+pow(fabs(z1-z2),n),n),1/n));
    return(result);
}

double camberra3Dist(double x1,double x2,double y1,double y2,double z1,double z2)
{ double result;
  result = fabs(x1-x2)/fabs(x1+x2)+fabs(y1-y2)/fabs(y1+y2)+fabs(z1-z2)/fabs(z1+z2);
  return(result);
}

double fasteuclidian3Dist(double x1,double x2,double y1,double y2,double z1,double z2)
{  double result;

  result = fabs(x1-x2)+0.5*fabs(y1-y2)+0.5*fabs(z1-z2);
  return(result);
}

double sqreuclidian3Dist(double x1,double x2,double y1,double y2,double z1,double z2)
{  double result;
    result = ((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
    return(result);
}

double maxnorma3Dist(double x1,double x2,double y1,double y2,double z1,double z2)
{  double result;

    result = max(fabs(x1-x2),max(fabs(y1-y2),fabs(z1-z2)) );
    return(result);
}

double xmanhattan3Dist(double x1,double x2,double y1,double y2,double z1,double z2)
{  double result;

  result = (fabs(x1-x2))*(fabs(y1-y2))*(fabs(z1-z2))*500.0+(fabs(x1-x2))+(fabs(y1-y2))+(fabs(z1-z2));
  if( result>0.0 ) {  result = log(sqr(result)); }else { result=0; }
  return(result);
}