void RGB2Lab(float R, float G, float B, float *L, float *a, float *b)
{
  R = R < 0 ? 0 : R > 255 ? 255 : R;
  G = G < 0 ? 0 : G > 255 ? 255 : G;
  B = B < 0 ? 0 : B > 255 ? 255 : B;

  double X, Y, Z, fX, fY, fZ;

  X = 0.412453*R + 0.357580*G + 0.180423*B;
  Y = 0.212671*R + 0.715160*G + 0.072169*B;
  Z = 0.019334*R + 0.119193*G + 0.950227*B;

  X /= (255 * 0.950456);
  Y /=  255;
  Z /= (255 * 1.088754);

  if (Y > 0.008856)
        {
          fY = pow(Y, 1.0/3.0);
          *L = static_cast<float>(116.0*fY - 16.0 + 0.5);
        }
  else
        {
          fY = 7.787*Y + 16.0/116.0;
          *L = static_cast<float>(903.3*Y + 0.5);
        }

  if (X > 0.008856)
        fX = pow(X, 1.0/3.0);
  else
        fX = 7.787*X + 16.0/116.0;

  if (Z > 0.008856)
        fZ = pow(Z, 1.0/3.0);
  else
        fZ = 7.787*Z + 16.0/116.0;

  *a = static_cast<float>(500.0*(fX - fY) + 0.5);
  *b = static_cast<float>(200.0*(fY - fZ) + 0.5);
}

void Lab2RGB(float L, float a, float b, float *R, float *G, float *B)
{
  double X, Y, Z, fX, fY, fZ;
  float RR, GG, BB;

  fY = pow((L + 16.0) / 116.0, 3.0);
  if (fY < 0.008856)
        fY = L / 903.3;
  Y = fY;

  if (fY > 0.008856)
        fY = pow(fY, 1.0/3.0);
  else
        fY = 7.787 * fY + 16.0/116.0;

  fX = a / 500.0 + fY;
  if (fX > 0.206893)
        X = pow(fX, 3.0);
  else
        X = (fX - 16.0/116.0) / 7.787;

  fZ = fY - b /200.0;
  if (fZ > 0.206893)
        Z = pow(fZ, 3.0);
  else
        Z = (fZ - 16.0/116.0) / 7.787;

  X *= (0.950456 * 255);
  Y *= 255;
  Z *= (1.088754 * 255);

  RR = static_cast<float>(3.240479*X - 1.537150*Y - 0.498535*Z + 0.5);
  GG = static_cast<float>(-0.969256*X + 1.875992*Y + 0.041556*Z + 0.5);
  BB = static_cast<float>(0.055648*X - 0.204043*Y + 1.057311*Z + 0.5);

  *R = RR < 0 ? 0 : RR > 255 ? 255 : RR;
  *G = GG < 0 ? 0 : GG > 255 ? 255 : GG;
  *B = BB < 0 ? 0 : BB > 255 ? 255 : BB;
}

float gaussian(size_t x1,size_t y1,size_t x2,size_t y2,float zigma){
      float distance = (x1-x2) * (x1-x2) + (y1 - y2) * (y1-y2);
      distance = sqrt(distance);
      float x = - (distance * distance) / (2 * zigma * zigma );
      return exp(x);
}
float gaussian(float distance,float zigma){
      float x = - (distance * distance) / (2 * zigma * zigma );
      return exp(x);
}