- PVSM.RU - https://www.pvsm.ru -
Многие слышали о великом и ужасном быстром преобразовании Фурье (БПФ / FFT — fast fourier transform) — но как его можно применять для решения практических задач за исключением JPEG/MPEG сжатия — зачастую остается неясным вопросом.
Недавно я наткнулся на интересную реализацию игры «Жизнь» Конвея, использующую быстрое преобразование Фурье(!!!) — и надеюсь, оно поможет вам понять всю силу и универсальность этого алгоритма.
Вспомним правила классической «Жизни» — на поле с квадратными клетками, живая клетка погибает если у неё больше 3 или меньше 2 соседей, и если у пустой клетки ровно 3 соседей — она рождается. Соответственно, для эффективной реализации алгоритма нужно быстро считать количество соседей вокруг клетки.
Алгоритмов для этого существует целая куча (в том числе и моя JS реализация [1]). Но есть у задачи и математическое решение, которое может давать хорошую скорость для плотно заполненных полей, и быстро уходит в отрыв с ростом сложности правил и площади/объема суммирования (например в Smoothlife [2]-подобных, и 3D вариантах)
Идея алгоритма следующая:
Весь этот процесс называется сверткой / сonvolution.
//Author: Mark VandeWettering http://brainwagon.org/2012/10/11/crazy-programming-experiment-of-the-evening/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <unistd.h>
#include <fftw3.h>
#define SIZE (512)
#define SHIFT (18)
fftw_complex *filter ;
fftw_complex *state ;
fftw_complex *tmp ;
fftw_complex *sum ;
int
main(int argc, char *argv[])
{
fftw_plan fwd, rev, flt ;
fftw_complex *ip, *jp ;
int x, y, g ;
srand48(getpid()) ;
filter = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
state = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
tmp = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
sum = (fftw_complex *) fftw_malloc(SIZE * SIZE * sizeof(fftw_complex)) ;
flt = fftw_plan_dft_2d(SIZE, SIZE,
filter, filter, FFTW_FORWARD, FFTW_ESTIMATE) ;
fwd = fftw_plan_dft_2d(SIZE, SIZE,
state, tmp, FFTW_FORWARD, FFTW_ESTIMATE) ;
rev = fftw_plan_dft_2d(SIZE, SIZE,
tmp, sum, FFTW_BACKWARD, FFTW_ESTIMATE) ;
/* initialize the state */
for (y=0, ip=state; y<SIZE; y++) {
for (x=0; x<SIZE; x++, ip++) {
*ip = (fftw_complex) (lrand48() % 2) ;
}
}
/* initialize and compute the filter */
for (y=0, ip=filter; y<SIZE; y++, ip++) {
for (x=0; x<SIZE; x++) {
*ip = 0. ;
}
}
#define IDX(x, y) (((x + SIZE) % SIZE) + ((y+SIZE) % SIZE) * SIZE)
filter[IDX(-1, -1)] = 1. ;
filter[IDX( 0, -1)] = 1. ;
filter[IDX( 1, -1)] = 1. ;
filter[IDX(-1, 0)] = 1. ;
filter[IDX( 1, 0)] = 1. ;
filter[IDX(-1, 1)] = 1. ;
filter[IDX( 0, 1)] = 1. ;
filter[IDX( 1, 1)] = 1. ;
fftw_execute(flt) ;
for (g = 0; g < 1000; g++) {
fprintf(stderr, "generation %03dr", g) ;
fftw_execute(fwd) ;
/* convolve */
for (y=0, ip=tmp, jp=filter; y<SIZE; y++) {
for (x=0; x<SIZE; x++, ip++, jp++) {
*ip *= *jp ;
}
}
/* go back to the sums */
fftw_execute(rev) ;
for (y=0, ip=state, jp=sum; y<SIZE; y++) {
for (x=0; x<SIZE; x++, ip++, jp++) {
int s = (int) round(creal(*ip)) ;
int t = ((int) round(creal(*jp))) >> SHIFT ;
if (s)
*ip = (t == 2 || t == 3) ;
else
*ip = (t == 3) ;
}
}
/* that's it! dump the frame! */
char fname[80] ;
sprintf(fname, "frame.%04d.pgm", g) ;
FILE *fp = fopen(fname, "wb") ;
fprintf(fp, "P5n%d %dn%dn", SIZE, SIZE, 255) ;
for (y=0, ip=state; y<SIZE; y++) {
for (x=0; x<SIZE; x++, ip++) {
int s = ((int) creal(*ip)) ;
fputc(255*s, fp) ;
}
}
fclose(fp) ;
}
fprintf(stderr, "n") ;
return 0 ;
}
Для сборки нужна библиотека FFTW [3]. Ключи для сборки в gcc:
gcc life.cpp -lfftw3 -lm -lstdc++
в Visual Studio нужны изменения в работе с комплексными числами.
Результат вполне ожидаемый:
«Закольцовывание» поля получается автоматически из-за БПФ.
Надеюсь, мистики в преобразовании Фурье теперь стало немного меньше.
Автор: BarsMonster
Источник [4]
Сайт-источник PVSM.RU: https://www.pvsm.ru
Путь до страницы источника: https://www.pvsm.ru/programmirovanie/34602
Ссылки в тексте:
[1] JS реализация: http://habrahabr.ru/post/144237/
[2] Smoothlife: http://habrahabr.ru/post/154509/
[3] FFTW: http://www.fftw.org/
[4] Источник: http://habrahabr.ru/post/180135/
Нажмите здесь для печати.