function.c 1.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100
  1. #include <math.h>
  2. #include <grass/gis.h>
  3. #include <grass/vector.h>
  4. #include "global.h"
  5. /* probability for gaussian distribution */
  6. double gaussian2dBySigma(double d, double sigma)
  7. {
  8. double res;
  9. res =
  10. 1. / (2. * M_PI * sigma * sigma) * exp(-d * d / (2. * sigma * sigma));
  11. return (res);
  12. }
  13. double gaussianFunction(double x, double sigma, double dimension)
  14. {
  15. return ((1. / (pow(2. * M_PI, dimension / 2.) * pow(sigma, dimension))) *
  16. exp(-0.5 * pow(x / sigma, 2.)));
  17. }
  18. /* probability for gaussian distribution */
  19. double gaussianKernel(double x, double term)
  20. {
  21. return (term * exp(-(x * x) / 2.));
  22. }
  23. /*
  24. term1 = 1./(2.*M_PI*sigma*sigma)
  25. term2 = 2.*sigma*sigma;
  26. */
  27. double gaussian2dByTerms(double d, double term1, double term2)
  28. {
  29. double res;
  30. res = term1 * exp(-d * d / term2);
  31. return (res);
  32. }
  33. double segno(double x)
  34. {
  35. double y;
  36. y = (x > 0 ? 1. : 0.) + (x < 0 ? -1. : 0.);
  37. return y;
  38. }
  39. double kernel1(double d, double rs, double lambda)
  40. {
  41. double res;
  42. double a = lambda - 1.;
  43. if (lambda == 1.) {
  44. res = 1. / (M_PI * (d * d + rs * rs));
  45. }
  46. else {
  47. res =
  48. segno(a) * (a / M_PI) * (pow(rs, 2. * a)) * (1 /
  49. pow(d * d + rs * rs,
  50. lambda));
  51. }
  52. /* res=1./(M_PI*(d*d+rs*rs)); */
  53. return (res);
  54. }
  55. double invGaussian2d(double sigma, double prob)
  56. {
  57. double d;
  58. d = sqrt(-2 * sigma * sigma * log(prob * M_PI * 2 * sigma * sigma));
  59. return (d);
  60. }
  61. /* euclidean distance between vectors x and y of length n */
  62. double euclidean_distance(double *x, double *y, int n)
  63. {
  64. int j;
  65. double out = 0.0;
  66. double tmp;
  67. for (j = 0; j < n; j++) {
  68. tmp = x[j] - y[j];
  69. out += tmp * tmp;
  70. }
  71. return sqrt(out);
  72. }