diff --git a/filter.blur.inc.c b/filter.blur.inc.c index d8c63bc..4f3449d 100644 --- a/filter.blur.inc.c +++ b/filter.blur.inc.c @@ -1,7 +1,7 @@ int toPowerOf2(int x) { - int xx; + int xx; for(xx = 1; x > 1; x >>= 1) xx <<= 1; return xx; @@ -10,7 +10,7 @@ int toPowerOf2(int x) { void fft(double *real, int realStep, double *img, int imgStep, int count, int invert) { - if (count < 2 || !realStep || !imgStep) return; + if (count < 2 || !realStep || !imgStep) return; count = toPowerOf2(count); // check order @@ -27,8 +27,8 @@ void fft(double *real, int realStep, double *img, int imgStep, int count, int in imgStep = -imgStep; } - double *realEnd = real + count*realStep; - double *imgEnd = img + count*imgStep; + double *realEnd = real + count*realStep; + double *imgEnd = img + count*imgStep; // flip if step was negative if (flipReal) @@ -39,40 +39,40 @@ void fft(double *real, int realStep, double *img, int imgStep, int count, int in { double x = *i; *i = *j; *j = x; } - // bit-reversal permutation - for(int i=0, j=0; i < count; ++i) { - if (j > i) { + // 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, + 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, + // process subsequences + for( double *iR = real + m*realStep, + *iI = img + m*imgStep, + *jR = iR + mmaxRealStep, *jI = iI + mmaxImgStep; iR < realEnd; iR += mmaxRealStep2, @@ -80,18 +80,18 @@ void fft(double *real, int realStep, double *img, int imgStep, int count, int in 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 + // 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; } @@ -99,30 +99,30 @@ void fft(double *real, int realStep, double *img, int imgStep, int count, int in 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; - } + // 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; + 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); + // 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); }