|
|
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 |
|