int toPowerOf2(int x) {
  int xx;
  for(xx = 1; x > 1; x >>= 1)
    xx <<= 1;
  return xx;
}



void fft(double *real, int realStep, double *img, int imgStep, int count, int invert) {
  if (count < 2 || !realStep || !imgStep) return;
  count = toPowerOf2(count);

  // check order
  int flipReal = 0;
  if (realStep < 0) {
    flipReal = 1;
    real += realStep*(count - 1);
    realStep = -realStep;
  }
  int flipImg = 0;
  if (imgStep < 0) {
    flipImg = 1;
    img += imgStep*(count - 1);
    imgStep = -imgStep;
  }

  double *realEnd = real + count*realStep;
  double *imgEnd  = img  + count*imgStep;

  // flip if step was negative
  if (flipReal)
    for(double *i = real, *j = realEnd - realStep; i < j; i += realStep, j -= realStep)
      { double x = *i; *i = *j; *j = x; }
  if (flipImg)
    for(double *i = img, *j = imgEnd - imgStep; i < j; i += imgStep, j -= imgStep)
      { double x = *i; *i = *j; *j = x; }


  // bit-reversal permutation
  for(int i=0, j=0; i < count; ++i) {
    if (j > i) {
      double x, *a, *b;
      a = real + i*realStep; b = real + j*realStep; x = *a; *a = *b; *b = x;
      a = img  + i*imgStep;  b = img  + j*imgStep;  x = *a; *a = *b; *b = x;
    }
    int m = count/2;
    while(m >= 1 && j >= m)
      { j -= m; m /= 2; }
    j += m;
  }

  double k = invert ? PI : -PI;
  for(int mmax = 1; mmax < count; mmax *= 2) {
    // rotation coefficients
    double angle = k/mmax;
    double wpR = sin(0.5*angle), wpI = sin(angle);
    wpR *= wpR*2;

    double wR = 1, wI = 0;
    int mmaxRealStep  = mmax*realStep;
    int mmaxRealStep2 = 2*mmaxRealStep;
    int mmaxImgStep   = mmax*imgStep;
    int mmaxImgStep2  = 2*mmaxImgStep;
    for(int m = 0; m < mmax; ++m) {
      // rotate w
      double wwR = wR;
      wR = wR - wwR*wpR - wI*wpI,
      wI = wI + wwR*wpI - wI*wpR;
      // process subsequences
      for( double *iR = real + m*realStep,
                  *iI = img  + m*imgStep,
                  *jR = iR + mmaxRealStep,
                  *jI = iI + mmaxImgStep;
           iR < realEnd;
           iR += mmaxRealStep2,
           iI += mmaxImgStep2,
           jR += mmaxRealStep2,
           jI += mmaxImgStep2 )
      {
        // radix
        double tR = *jR*wR - *jI*wI;
        double tI = *jR*wI + *jI*wR;
        *jR = *iR - tR;
        *jI = *iI - tI;
        *iR += tR;
        *iI += tI;
      }
    }
  }

  // reverse order
  if (!flipReal)
    for(double *i = real, *j = realEnd - realStep; i < j; i += realStep, j -= realStep)
      { double x = *i; *i = *j; *j = x; }
  if (!flipImg)
    for(double *i = img, *j = imgEnd - imgStep; i < j; i += imgStep, j -= imgStep)
      { double x = *i; *i = *j; *j = x; }

  // divide by count to complete back-FFT
  if (invert) {
    double k = 1.0/count;
    for(double *i = real; i < realEnd; i += realStep)
      *i *= k;
    for(double *i = img; i < imgEnd; i += imgStep)
      *i *= k;
  }
}


void fft2d(double *real, int realColStep, int realRowStep, double *img, int imgColStep, int imgRowStep, int rows, int cols, int invert) {
  if (rows < 1 || cols < 1 || !realColStep || !imgColStep || !realRowStep || imgRowStep) return;
  rows = toPowerOf2(rows);
  cols = toPowerOf2(cols);

  // fft rows
  if (cols > 1)
    for(double *r = real, *i = img, *end = real + rows*realRowStep; r < end; r += realRowStep, i += imgRowStep)
      fft(r, realColStep, i, imgColStep, cols, invert);
  // fft cols
  if (rows > 1)
    for(double *r = real, *i = img, *end = real + cols*realColStep; r < end; r += realColStep, i += imgColStep)
      fft(r, realRowStep, i, imgRowStep, rows, invert);
}


double preGauss(double x, double r) {
  if (fabs(r) < IMG_PRECISION)
    return fabs(x) < IMG_PRECISION ? 1 : 0;
  return exp(-0.5*x*x/(r*r));
}


double gauss(double x, double r) {
  static const double k = 1/sqrt(2*PI);
  return k*preGauss(x, r)/fabs(r);
}


void gaussBlurRows(double *row, int colStep, int rowStep, int cols, int rows, double size) {
  if (!(size > IMG_PRECISION)) return;
  int cnt = cols + 2*(int)ceil(size*4 - IMG_PRECISION);
  int p2cnt = toPowerOf2(cnt);
  if (p2cnt < cnt) p2cnt <<= 1;
  cnt = p2cnt;

  double *buf = (double*)calloc(1, sizeof(double)*4*cnt);
  double *pat = buf + 2*cnt;
  pat[0]   = preGauss(0,     size);
  pat[cnt] = preGauss(cnt/2, size);
  for(int i = 2; i < cnt; i += 2)
    pat[i] = pat[cnt*2 - i] = preGauss(i, size);
  fft(pat, 2, pat+1, 2, cnt, FALSE);

  double *buf0 = buf + (cnt - cols);
  double *buf1 = buf0 + (cols - 1)*2;
  for(int i = 0; i < rows; ++i) {
    for(int ch = 0; ch < 4; ++ch) {
      memset(buf, 0, sizeof(double)*2*cnt);
      double *chrow = row + i*rowStep + ch;
      for(int j = 0; j < cols; ++j)
        buf0[j*2] = chrow[j*colStep];
      for(double v = *buf0, *b = buf;      b < buf0; b += 2) *b = v;
      for(double v = *buf1, *b = buf1 + 2; b < pat;  b += 2) *b = v;
      fft(buf, 2, buf+1, 2, cnt, FALSE);
      for(double *b = buf, *p = pat; b < pat; b += 2, p += 2) {
        double bR = b[0], bI = b[1];
        b[0] = bR*p[0] - bI*p[1];
        b[1] = bR*p[1] + bI*p[0];
      }
      fft(buf, 2, buf+1, 2, cnt, TRUE);
      for(int j = 0; j < cols; ++j) {
        double *b = buf0 + j*2;
        chrow[j*colStep] = sqrt(b[0]*b[0] + b[1]*b[1]);
      }
    }
  }
  free(buf);
}


void filterBlurGauss(Img *img, double sx, double sy) {
  if (!img->data) return;
  sx = fabs(sx);
  sy = fabs(sy);
  if (!(sx > IMG_PRECISION) && !(sy > IMG_PRECISION)) return;

  filterMulAlpha(img);
  gaussBlurRows(img->data, 4, 4*img->w, img->w, img->h, sx);
  gaussBlurRows(img->data, 4*img->w, 4, img->h, img->w, sy);
  filterDivAlpha(img);
}