Экранированное решение Пуассона с FFTW3

Question

Я хочу решить экранированное уравнение Пуассона, описанное в этой статье , чтобы вычислить изображение из грубого приближения и его вертикальных и горизонтальных градиентов.Я использую FFTW3 впервые для этого.Однако у меня очень странные результаты.

Мое грубое приближение к изображению выглядит следующим образом:

Но вывод моего решения выглядит следующим образом:

.

Чтобы отладить его, я попробовал несколько вещей -

Сначала я подозревал, что мои вычисления пространственных частот были отключены, поэтому я попытался вычислитьгоризонтальный градиент изображения в области Фурье, в результате чего:

.

Наконец, когда я вычисляю только FT моего изображения и его обратное преобразование, я получаю это:

Но когда я просто копирую FT вновый массив и вычислить его обратный FT, я получаю это:

.

Вот моя функция:

void RayTrace::computeFT(double* &image_in, double* &dx_in, double* &dy_in, double* &result){
    fftw_complex image_FT[width*height];
    fftw_complex dx_FT[width*height];
    fftw_complex dy_FT[width*height];
    fftw_plan image_plan = fftw_plan_dft_r2c_2d(width, height, image_in, image_FT, FFTW_ESTIMATE);
    fftw_plan dx_plan = fftw_plan_dft_r2c_2d(width, height, dx_in, dx_FT, FFTW_ESTIMATE);
    fftw_plan dy_plan = fftw_plan_dft_r2c_2d(width, height, dy_in, dy_FT, FFTW_ESTIMATE);
    //compute the FTs
    fftw_execute(image_plan);
    fftw_execute(dx_plan);
    fftw_execute(dy_plan);
    //delete plans after that
    fftw_destroy_plan(image_plan);
    fftw_destroy_plan(dx_plan);
    fftw_destroy_plan(dy_plan);
    fftw_cleanup();
    //modify the image in the Fourier domain
    fftw_complex result_FT[width*height];
    for(int i=0 ; i<height ; i++){
        for(int j=0 ; j<width ; j++){
            float sx = (i>height/2 ? i-height : i) / float(height);
            float sy = (j>width/2 ? j-width : j) / float(width);
            int idx = i*width+j;
            //These are the things that I tried to do
            //compute the screened Poisson solution
            result_FT[idx][0] = (alpha*image_FT[idx][0] + 2*M_PI*sx*dx_FT[idx][1] + 2*M_PI*sy*dy_FT[idx][1])/(alpha + 4*M_PI*M_PI*(sx*sx+sy*sy));
            result_FT[idx][1] = (alpha*image_FT[idx][1] - 2*M_PI*sx*dx_FT[idx][0] - 2*M_PI*sy*dy_FT[idx][0])/(alpha + 4*M_PI*M_PI*(sx*sx+sy*sy));
            //horizontal gradient
            result_FT[idx][0] = -2*M_PI*sx*image_FT[idx][1];
            result_FT[idx][1] = 2*M_PI*sx*image_FT[idx][0];
            //simple copy
            result_FT[idx][0] = image_FT[idx][0];
            result_FT[idx][1] = image_FT[idx][1];
        }
    }
    //inverse transform
    fftw_plan r_back_plan = fftw_plan_dft_c2r_2d(width, height, image_FT, result, FFTW_ESTIMATE);
    fftw_execute(r_back_plan);
    //delete plan after that
    fftw_destroy_plan(r_back_plan);
    fftw_cleanup();

    //normalize the result
    float norm = height * width;
    for(int i=0 ; i<height ; i++){
        for(int j=0 ; j<width ; j++){
            result[i*width + j] /= norm;
        }
    }
}

Пожалуйста, сообщите мне, если вы видите ошибку в моей реализации

Кроме того, вот метод, из которого я вызываю метод,один раз для каждого цветового компонента в моем изображении:

void RayTrace::computeGradient(glm::vec3* &image, glm::vec3* &xGradient, glm::vec3* &zGradient, glm::vec3* &gradient){
//compute the gradient of the image in the fourier domain
double *r_in = new double[width*height];
double *r_dx_in = new double[width*height];
double *r_dy_in = new double[width*height];
double *r_out = new double[width*height];
double *g_in = new double[width*height];
double *g_dx_in = new double[width*height];
double *g_dy_in = new double[width*height];
double *g_out = new double[width*height];
double *b_in = new double[width*height];
double *b_dx_in = new double[width*height];
double *b_dy_in = new double[width*height];
double *b_out = new double[width*height];
for(int i=0 ; i<height ; i++){
    for(int j=0 ; j<width ; j++){
        r_in[i*width + j] = image[i*width + j].x;
        r_dx_in[i*width + j] = xGradient[i*width + j].x;
        r_dy_in[i*width + j] = zGradient[i*width + j].x;
        g_in[i*width + j] = image[i*width + j].y;
        g_dx_in[i*width + j] = xGradient[i*width + j].y;
        g_dy_in[i*width + j] = zGradient[i*width + j].y;
        b_in[i*width + j] = image[i*width + j].z;
        b_dx_in[i*width + j] = xGradient[i*width + j].z;
        b_dy_in[i*width + j] = zGradient[i*width + j].z;
    }
}
std::cout << width << " " << height << std::endl;
computeFT(r_in, r_dx_in, r_dy_in, r_out);
computeFT(g_in, g_dx_in, g_dy_in, g_out);
computeFT(b_in, b_dx_in, b_dy_in, b_out);
for(int i=0 ; i<height ; i++){
    for(int j=0 ; j<width ; j++){
        gradient[i*width + j] = glm::vec3(r_out[i*width + j], g_out[i*width + j], b_out[i*width + j]);
    }
}
delete r_in;
delete r_out;
delete r_dx_in;
delete r_dy_in;
delete g_in;
delete g_out;
delete g_dx_in;
delete g_dy_in;
delete b_in;
delete b_out;
delete b_dx_in;
delete b_dy_in;

}

РЕДАКТИРОВАТЬ: Оказывается, что мой вывод трассировки пути вообще не был нормализован, поэтому значения в моих массивах варьировались от 1eОт -3 до 1e6, так что, думаю, у FFTW возникли проблемы с масштабированием.Теперь, когда я нормализовал его, я могу вернуть свое изображение после вычисления его FFT, а затем IFFT без изменений, что является некоторым прогрессом.Тем не менее, вычисление градиента дает следующее:

Так что я думаю, что все еще есть проблема с моим вычислением пространственной частоты