Blame simple/imgfilters/filter.blur.inc.c

608c44
608c44
608c44
int toPowerOf2(int x) {
608c44
  int xx;
608c44
  for(xx = 1; x > 1; x >>= 1)
608c44
    xx <<= 1;
608c44
  return xx;
608c44
}
608c44
608c44
608c44
608c44
void fft(double *real, int realStep, double *img, int imgStep, int count, int invert) {
608c44
  if (count < 2 || !realStep || !imgStep) return;
608c44
  count = toPowerOf2(count);
608c44
608c44
  // check order
608c44
  int flipReal = 0;
608c44
  if (realStep < 0) {
608c44
    flipReal = 1;
608c44
    real += realStep*(count - 1);
608c44
    realStep = -realStep;
608c44
  }
608c44
  int flipImg = 0;
608c44
  if (imgStep < 0) {
608c44
    flipImg = 1;
608c44
    img += imgStep*(count - 1);
608c44
    imgStep = -imgStep;
608c44
  }
608c44
608c44
  double *realEnd = real + count*realStep;
608c44
  double *imgEnd  = img  + count*imgStep;
608c44
608c44
  // flip if step was negative
608c44
  if (flipReal)
608c44
    for(double *i = real, *j = realEnd - realStep; i < j; i += realStep, j -= realStep)
608c44
      { double x = *i; *i = *j; *j = x; }
608c44
  if (flipImg)
608c44
    for(double *i = img, *j = imgEnd - imgStep; i < j; i += imgStep, j -= imgStep)
608c44
      { double x = *i; *i = *j; *j = x; }
608c44
608c44
608c44
  // bit-reversal permutation
608c44
  for(int i=0, j=0; i < count; ++i) {
608c44
    if (j > i) {
608c44
      double x, *a, *b;
608c44
      a = real + i*realStep; b = real + j*realStep; x = *a; *a = *b; *b = x;
608c44
      a = img  + i*imgStep;  b = img  + j*imgStep;  x = *a; *a = *b; *b = x;
608c44
    }
608c44
    int m = count/2;
608c44
    while(m >= 1 && j >= m)
608c44
      { j -= m; m /= 2; }
608c44
    j += m;
608c44
  }
608c44
608c44
  double k = invert ? PI : -PI;
608c44
  for(int mmax = 1; mmax < count; mmax *= 2) {
608c44
    // rotation coefficients
608c44
    double angle = k/mmax;
608c44
    double wpR = sin(0.5*angle), wpI = sin(angle);
608c44
    wpR *= wpR*2;
608c44
608c44
    double wR = 1, wI = 0;
608c44
    int mmaxRealStep  = mmax*realStep;
608c44
    int mmaxRealStep2 = 2*mmaxRealStep;
608c44
    int mmaxImgStep   = mmax*imgStep;
608c44
    int mmaxImgStep2  = 2*mmaxImgStep;
608c44
    for(int m = 0; m < mmax; ++m) {
608c44
      // rotate w
608c44
      double wwR = wR;
608c44
      wR = wR - wwR*wpR - wI*wpI,
608c44
      wI = wI + wwR*wpI - wI*wpR;
608c44
      // process subsequences
608c44
      for( double *iR = real + m*realStep,
608c44
                  *iI = img  + m*imgStep,
608c44
                  *jR = iR + mmaxRealStep,
608c44
                  *jI = iI + mmaxImgStep;
608c44
           iR < realEnd;
608c44
           iR += mmaxRealStep2,
608c44
           iI += mmaxImgStep2,
608c44
           jR += mmaxRealStep2,
608c44
           jI += mmaxImgStep2 )
608c44
      {
608c44
        // radix
608c44
        double tR = *jR*wR - *jI*wI;
608c44
        double tI = *jR*wI + *jI*wR;
608c44
        *jR = *iR - tR;
608c44
        *jI = *iI - tI;
608c44
        *iR += tR;
608c44
        *iI += tI;
608c44
      }
608c44
    }
608c44
  }
608c44
608c44
  // reverse order
608c44
  if (!flipReal)
608c44
    for(double *i = real, *j = realEnd - realStep; i < j; i += realStep, j -= realStep)
608c44
      { double x = *i; *i = *j; *j = x; }
608c44
  if (!flipImg)
608c44
    for(double *i = img, *j = imgEnd - imgStep; i < j; i += imgStep, j -= imgStep)
608c44
      { double x = *i; *i = *j; *j = x; }
608c44
608c44
  // divide by count to complete back-FFT
608c44
  if (invert) {
608c44
    double k = 1.0/count;
608c44
    for(double *i = real; i < realEnd; i += realStep)
608c44
      *i *= k;
608c44
    for(double *i = img; i < imgEnd; i += imgStep)
608c44
      *i *= k;
608c44
  }
608c44
}
608c44
608c44
608c44
void fft2d(double *real, int realColStep, int realRowStep, double *img, int imgColStep, int imgRowStep, int rows, int cols, int invert) {
608c44
  if (rows < 1 || cols < 1 || !realColStep || !imgColStep || !realRowStep || imgRowStep) return;
608c44
  rows = toPowerOf2(rows);
608c44
  cols = toPowerOf2(cols);
608c44
608c44
  // fft rows
608c44
  if (cols > 1)
608c44
    for(double *r = real, *i = img, *end = real + rows*realRowStep; r < end; r += realRowStep, i += imgRowStep)
608c44
      fft(r, realColStep, i, imgColStep, cols, invert);
608c44
  // fft cols
608c44
  if (rows > 1)
608c44
    for(double *r = real, *i = img, *end = real + cols*realColStep; r < end; r += realColStep, i += imgColStep)
608c44
      fft(r, realRowStep, i, imgRowStep, rows, invert);
608c44
}
608c44
608c44
608c44
double preGauss(double x, double r) {
608c44
  if (fabs(r) < IMG_PRECISION)
608c44
    return fabs(x) < IMG_PRECISION ? 1 : 0;
608c44
  return exp(-0.5*x*x/(r*r));
608c44
}
608c44
608c44
608c44
double gauss(double x, double r) {
608c44
  static const double k = 1/sqrt(2*PI);
608c44
  return k*preGauss(x, r)/fabs(r);
608c44
}
608c44
608c44
608c44
void gaussBlurRows(double *row, int colStep, int rowStep, int cols, int rows, double size) {
608c44
  if (!(size > IMG_PRECISION)) return;
608c44
  int cnt = cols + 2*(int)ceil(size*4 - IMG_PRECISION);
608c44
  int p2cnt = toPowerOf2(cnt);
608c44
  if (p2cnt < cnt) p2cnt <<= 1;
608c44
  cnt = p2cnt;
608c44
608c44
  double *buf = (double*)calloc(1, sizeof(double)*4*cnt);
608c44
  double *pat = buf + 2*cnt;
608c44
  pat[0]   = preGauss(0,     size);
608c44
  pat[cnt] = preGauss(cnt/2, size);
608c44
  for(int i = 2; i < cnt; i += 2)
608c44
    pat[i] = pat[cnt*2 - i] = preGauss(i, size);
608c44
  fft(pat, 2, pat+1, 2, cnt, FALSE);
608c44
608c44
  double *buf0 = buf + (cnt - cols);
608c44
  double *buf1 = buf0 + (cols - 1)*2;
608c44
  for(int i = 0; i < rows; ++i) {
608c44
    for(int ch = 0; ch < 4; ++ch) {
608c44
      memset(buf, 0, sizeof(double)*2*cnt);
608c44
      double *chrow = row + i*rowStep + ch;
608c44
      for(int j = 0; j < cols; ++j)
608c44
        buf0[j*2] = chrow[j*colStep];
608c44
      for(double v = *buf0, *b = buf;      b < buf0; b += 2) *b = v;
608c44
      for(double v = *buf1, *b = buf1 + 2; b < pat;  b += 2) *b = v;
608c44
      fft(buf, 2, buf+1, 2, cnt, FALSE);
608c44
      for(double *b = buf, *p = pat; b < pat; b += 2, p += 2) {
608c44
        double bR = b[0], bI = b[1];
608c44
        b[0] = bR*p[0] - bI*p[1];
608c44
        b[1] = bR*p[1] + bI*p[0];
608c44
      }
608c44
      fft(buf, 2, buf+1, 2, cnt, TRUE);
608c44
      for(int j = 0; j < cols; ++j) {
608c44
        double *b = buf0 + j*2;
608c44
        chrow[j*colStep] = sqrt(b[0]*b[0] + b[1]*b[1]);
608c44
      }
608c44
    }
608c44
  }
608c44
  free(buf);
608c44
}
608c44
608c44
608c44
void filterBlurGauss(Img *img, double sx, double sy) {
608c44
  if (!img->data) return;
608c44
  sx = fabs(sx);
608c44
  sy = fabs(sy);
608c44
  if (!(sx > IMG_PRECISION) && !(sy > IMG_PRECISION)) return;
608c44
608c44
  filterMulAlpha(img);
608c44
  gaussBlurRows(img->data, 4, 4*img->w, img->w, img->h, sx);
608c44
  gaussBlurRows(img->data, 4*img->w, 4, img->h, img->w, sy);
608c44
  filterDivAlpha(img);
608c44
}
608c44
608c44