Вперед FFT изображение и назад FFT изображение, чтобы получить тот же результат - PullRequest
Новая
       19

Вперед FFT изображение и назад FFT изображение, чтобы получить тот же результат

8 голосов
/ 17 октября 2011

Я пытаюсь БПФ изображение, используя библиотеку от http://www.fftw.org/, чтобы я мог сделать свертку в частотной области. Но я не могу понять, как заставить это работать. Чтобы понять, как это сделать, я пытаюсь перенаправить БПФ изображение в виде массива пикселей, а затем обратно БПФ, чтобы получить тот же массив пикселей. Вот что я делаю: 

fftw_plan planR, planG, planB;
fftw_complex *inR, *inG, *inB, *outR, *outG, *outB, *resultR, *resultG, *resultB;

//Allocate arrays.
inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);

outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);

resultR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
resultG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
resultB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);

//Fill in arrays with the pixelcolors.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        inR[y * width + x][0] = pixelColors[currentIndex];
        inG[y * width + x][0] = pixelColors[currentIndex + 1];
        inB[y * width + x][0] = pixelColors[currentIndex + 2];
    }
}

//Forward plans.
planR = fftw_plan_dft_2d(width, width, inR, outR, FFTW_FORWARD, FFTW_MEASURE);
planG = fftw_plan_dft_2d(width, width, inG, outG, FFTW_FORWARD, FFTW_MEASURE);
planB = fftw_plan_dft_2d(width, width, inB, outB, FFTW_FORWARD, FFTW_MEASURE);

//Forward FFT.
fftw_execute(planR);
fftw_execute(planG);
fftw_execute(planB);

//Backward plans.
planR = fftw_plan_dft_2d(width, width, outR, resultR, FFTW_BACKWARD, FFTW_MEASURE);
planG = fftw_plan_dft_2d(width, width, outG, resultG, FFTW_BACKWARD, FFTW_MEASURE);
planB = fftw_plan_dft_2d(width, width, outB, resultB, FFTW_BACKWARD, FFTW_MEASURE);

//Backward fft
fftw_execute(planR);
fftw_execute(planG);
fftw_execute(planB);

//Overwrite the pixelcolors with the result.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        pixelColors[currentIndex] = resultR[y * width + x][0];
        pixelColors[currentIndex + 1] = resultG[y * width + x][0];
        pixelColors[currentIndex + 2] = resultB[y * width + x][0];
    }
}

Может кто-нибудь показать мне пример того, как переслать FFT изображение, а затем назад FFT изображение, используя FFTW, чтобы получить тот же результат? Я смотрел на множество примеров, показывающих, как использовать FFTW для FFT, но я не могу понять, как это применимо к моей ситуации, когда у меня есть массив пиксельных цветов, представляющих изображение.

Ответы [ 2 ]

16 голосов
/ 17 октября 2011

Одна важная вещь, которую следует отметить, когда вы выполняете прямое БПФ с последующим обратным БПФ, это то, что это обычно приводит к тому, что к конечному результату применяется масштабный коэффициент N, то есть результирующие значения пикселей изображения должны быть разделены на N по порядку.чтобы соответствовать исходным значениям пикселей.(N - размер БПФ.) Таким образом, ваш цикл вывода должен выглядеть примерно так: 

//Overwrite the pixelcolors with the result.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        pixelColors[currentIndex] = resultR[y * width + x][0] / (width * height);
        pixelColors[currentIndex + 1] = resultG[y * width + x][0] / (width * height);
        pixelColors[currentIndex + 2] = resultB[y * width + x][0] / (width * height);
    }
}

Также обратите внимание, что вы, вероятно, хотите сделать БПФ от реального к сложному, за которым следует комплексноеРеальный IFFT (несколько более эффективный с точки зрения памяти и производительности).Пока что, похоже, вы делаете сложное в обоих направлениях, что нормально, но вы не правильно заполняете свои входные массивы.Если вы собираетесь придерживаться комплекса «комплекс», то вы, вероятно, захотите изменить свой цикл ввода на что-то вроде этого: 

//Fill in arrays with the pixelcolors.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        inR[y * width + x][0] = (double)pixelColors[currentIndex];
        inR[y * width + x][1] = 0.0;
        inG[y * width + x][0] = (double)pixelColors[currentIndex + 1];
        inG[y * width + x][1] = 0.0;
        inB[y * width + x][0] = (double)pixelColors[currentIndex + 2];
        inB[y * width + x][1] = 0.0;
    }
}

т.е. значения пикселов входят в реальные части комплексных значений вводаи воображаемые части должны быть обнулены.

Еще одна вещь, на которую стоит обратить внимание: когда вы в конечном итоге получите эту работу, вы обнаружите, что производительность ужасна - создание плана относительно времени, затраченного на это, занимает много временидля фактического БПФ.Идея состоит в том, что вы создаете план только один раз, но используете его для выполнения многих БПФ.Таким образом, вы захотите отделить создание плана от фактического кода FFT и поместить его в подпрограмму инициализации, конструктор или что-то еще.

2 голосов
/ 18 октября 2011

Но если вы используете realToComplex или ComplexToRealFunction, обратите внимание на тот факт, что изображение будет сохранено в матрице измерений [высота x (ширина / 2 +1)], и если вы хотите выполнить некоторые промежуточные вычисления вВ частотной области они станут немного сложнее ...

...