random_phase_noise_lib.c

Go to the documentation of this file.

/*
 * Copyright (c) 2010, Bruno Galerne <bruno.galerne@cmla.ens-cachan.fr>
 * All rights reserved.
 *
 * This program is free software: you can use, modify and/or
 * redistribute it under the terms of the GNU General Public
 * License as published by the Free Software Foundation, either
 * version 3 of the License, or (at your option) any later
 * version. You should have received a copy of this license along
 * this program. If not, see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>

#include "mt.h"

/* M_PI is a POSIX definition */
#ifndef M_PI

#define M_PI 3.14159265358979323846
#endif                          /* !M_PI */

/*
* number of threads to use for libfftw
* (uncomment to enable parallel FFT multi-threading)
*/
/* #define FFTW_NTHREADS 4 */

/*
* Preliminary section: some useful standard operations
*/

static float mean(float *data_in, size_t N)
{
    float m = 0.;
    float *ptr, *ptr_end;

    /* Pointers initialization */
    ptr = data_in;
    ptr_end = ptr + N;

    /* Sum loop */
    while (ptr < ptr_end) {
        m += *ptr;
        ptr++;
    }

    /* Normalization */
    m /= ((float) N);

    return (m);
}

static fftwf_complex
    * pointwise_complexfloat_multiplication(fftwf_complex * comp_in,
                                            float *float_in, size_t N)
{
    fftwf_complex *ptr_comp_in, *ptr_comp_end;
    float *ptr_float_in;

    /* check allocaton */
    if (NULL == comp_in || NULL == float_in)
        return (NULL);

    ptr_comp_in = comp_in;
    ptr_float_in = float_in;
    ptr_comp_end = ptr_comp_in + N;
    while (ptr_comp_in < ptr_comp_end) {
        (*ptr_comp_in)[0] *= (*ptr_float_in);
        (*ptr_comp_in)[1] *= (*ptr_float_in);
        ptr_comp_in++;
        ptr_float_in++;
    }
    return (comp_in);
}

static fftwf_complex
    * pointwise_complex_multiplication(fftwf_complex * comp_in1,
                                       fftwf_complex * comp_in2, size_t N)
{
    fftwf_complex *ptr_comp_in1, *ptr_comp_in2, *ptr_comp_end;
    float re1, im1;             /* real and imaginary part of *ptr_comp_in1 */

    /* Check allocaton */
    if (NULL == comp_in1 || NULL == comp_in2)
        return (NULL);

    ptr_comp_in1 = comp_in1;
    ptr_comp_in2 = comp_in2;
    ptr_comp_end = ptr_comp_in1 + N;
    while (ptr_comp_in1 < ptr_comp_end) {
        re1 = (*ptr_comp_in1)[0];
        im1 = (*ptr_comp_in1)[1];
        (*ptr_comp_in1)[0] = re1 * (*ptr_comp_in2)[0]
            - im1 * (*ptr_comp_in2)[1];
        (*ptr_comp_in1)[1] = re1 * (*ptr_comp_in2)[1]
            + im1 * (*ptr_comp_in2)[0];
        ptr_comp_in1++;
        ptr_comp_in2++;
    }
    return (comp_in1);
}

/*
* Periodic component section
*/

static float *discrete_laplacian(float *data_out, float *data_in,
                                 size_t nx, size_t ny)
{
    int i, j;
    float *out_ptr;
    float *in_ptr, *in_ptr_xm1, *in_ptr_xp1, *in_ptr_ym1, *in_ptr_yp1;

    /* check allocation */
    if (NULL == data_in || NULL == data_out)
        return (NULL);

    /* pointers to the data and neighbour values */
    in_ptr = data_in;
    in_ptr_xm1 = data_in - 1;
    in_ptr_xp1 = data_in + 1;
    in_ptr_ym1 = data_in - nx;
    in_ptr_yp1 = data_in + nx;
    out_ptr = data_out;
    /* iterate on j, i, following the array order */
    for (j = 0; j < (int) ny; j++) {
        for (i = 0; i < (int) nx; i++) {
            *out_ptr = 0.;
            /* row differences */
            if (0 < i) {
                *out_ptr += *(in_ptr_xm1) - *in_ptr;
            }
            if ((int) nx - 1 > i) {
                *out_ptr += *(in_ptr_xp1) - *in_ptr;
            }
            /* column differences */
            if (0 < j) {
                *out_ptr += *(in_ptr_ym1) - *in_ptr;
            }
            if ((int) ny - 1 > j) {
                *out_ptr += *(in_ptr_yp1) - *in_ptr;
            }
            in_ptr++;
            in_ptr_xm1++;
            in_ptr_xp1++;
            in_ptr_ym1++;
            in_ptr_yp1++;
            out_ptr++;
        }
    }
    return (data_out);
}

static float *poisson_complex_filter(float *data_out, size_t nx, size_t ny)
{
    float *cosx, *cosymintwo, *ptr_x, *ptr_y, *ptr_x_end,
        *ptr_y_end, *out_ptr;
    int i, j;
    size_t sx = nx / 2 + 1;
    int N = nx * ny;
    float halfinvN = 0.5 / ((float) N);

    /* check allocation */
    if (NULL == data_out)
        return (NULL);

    /* allocate memory */
    if (NULL == (cosx = (float *) malloc(sx * sizeof(float)))
        || NULL == (cosymintwo = (float *) malloc(ny * sizeof(float))))
        return (NULL);

    /* define cosx and cosymintwo */
    ptr_x = cosx;
    for (i = 0; i < (int) sx; i++) {
        (*ptr_x) = cos(2. * ((float) i) * M_PI / ((float) nx));
        ptr_x++;
    }
    ptr_y = cosymintwo;
    for (j = 0; j < (int) ny; j++) {
        (*ptr_y) = cos(2. * ((float) j) * M_PI / ((float) ny)) - 2.;
        ptr_y++;
    }

    /* define data_out */
    out_ptr = data_out;
    ptr_y = cosymintwo;
    ptr_y_end = cosymintwo + ny;
    ptr_x = cosx;
    ptr_x_end = cosx + sx;
    /* particular case : (0,0) */
    (*out_ptr) = 1.;
    out_ptr++;
    ptr_x++;
    while (ptr_y < ptr_y_end) {
        while (ptr_x < ptr_x_end) {
            (*out_ptr) = halfinvN / ((*ptr_x) + (*ptr_y));
            out_ptr++;
            ptr_x++;
        }
        ptr_x = cosx;
        ptr_y++;
    }

    /* free memory */
    free(cosx);
    free(cosymintwo);

    /* return data_out */
    return (data_out);
}

static float **periodic_component_color(float **data_in_rgb,
                                        size_t nx, size_t ny)
{
    float *laplacian;           /* 2D discrete Laplacian */
    float *pcf;                 /* Poisson complex filter */
    fftwf_plan plan_r2c, plan_c2r;
    fftwf_complex *fft;
    int N = ((int) nx) * ((int) ny);
    float m;                    /* mean of the coefficient of data_in */
    int channel;
    size_t fft_size = (nx / 2 + 1) * ny;        /* physical size of the r2c FFT */

    /* Check allocation */
    if (NULL == data_in_rgb)
        return (NULL);

    /* Start threaded fftw if FFTW_NTHREADS is defined */
#ifdef FFTW_NTHREADS
    if (0 == fftwf_init_threads())
        return (NULL);
    fftwf_plan_with_nthreads(FFTW_NTHREADS);
#endif

    /* Allocate memory */
    if (NULL == (pcf = (float *) malloc(fft_size * sizeof(float)))
        || NULL == (laplacian = (float *) malloc(N * sizeof(float)))
        || NULL == (fft = (fftwf_complex *)
                    fftwf_malloc(fft_size * sizeof(fftwf_complex))))
        return (NULL);

    /* Compute the Poisson complex filter */
    if (NULL == poisson_complex_filter(pcf, nx, ny))
        return (NULL);

    /* For each channel */
    for (channel = 0; channel < 3; channel++) {
        /* Compute m = the mean value of data_in[channel] */
        m = mean(data_in_rgb[channel], N);

        /* Compute the discrete Laplacian of data_in[channel] */
        if (NULL == discrete_laplacian(laplacian,
                                       data_in_rgb[channel], nx, ny))
            return (NULL);

        /* Forward Fourier transform of the Laplacian: create
         * the FFT forward plan and run the FFT */
        plan_r2c = fftwf_plan_dft_r2c_2d((int) (ny), (int) (nx),
                                         laplacian, fft,
                                         FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
        fftwf_execute(plan_r2c);

        /* Inverse the discrete periodic Laplacian in the
         * Fourier domain using the Poisson complex filter */
        if (NULL == pointwise_complexfloat_multiplication(fft, pcf, fft_size))
            return (NULL);
        /* (0,0) frequency : we impose the same mean as the input */
        (*fft)[0] = m;
        (*fft)[1] = 0.;

        /* Backward Fourier transform: the output is stored in
         * data_in_rgb[channel] */
        plan_c2r = fftwf_plan_dft_c2r_2d((int) (ny), (int) (nx),
                                         fft, data_in_rgb[channel],
                                         FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
        fftwf_execute(plan_c2r);
    }

    /* cleanup */
    free(pcf);
    free(laplacian);
    /* destroy the FFT plans and data */
    fftwf_destroy_plan(plan_r2c);
    fftwf_destroy_plan(plan_c2r);
    fftwf_free(fft);
    fftwf_cleanup();
#ifdef FFTW_NTHREADS
    fftwf_cleanup_threads();
#endif

    return (data_in_rgb);
}

/*
* Phase randomization section
*/

static fftwf_complex *random_phase(fftwf_complex * data_out, size_t nx,
                                   size_t ny)
{
    int inx, iny;               /* int version of nx and ny */
    int hinx, hiny;             /* half of nx and ny */
    float invN;                 /* inverse of nx*ny */
    int xev, yev;               /* 0 if x is even, 1 otherwise */
    fftwf_complex *comp_ptr;    /* complex pointer */
    int sx = (nx / 2 + 1);      /* physical x-size of data_out */
    int x, y;                   /* index for for loops */
    int ad;                     /* pixels locations */
    double theta;               /* random phase */
    float costheta, sintheta;   /* cosine and sine of theta */
    float sign;                 /* random sign */

    /* Check allocation */
    if (NULL == data_out)
        return (NULL);

    /* Define sizes and indexes for for loops: depend of the parity */
    inx = (int) nx;
    iny = (int) ny;
    invN = 1. / ((float) (inx * iny));
    xev = inx % 2;
    yev = iny % 2;
    hinx = inx / 2;
    hiny = iny / 2;

    /* Pointer initialization */
    comp_ptr = data_out;

    for (y = 0; y < iny; y++) {
        for (x = 0; x < sx; x++) {
            /* Case 1: (x,y) is its own symmetric: (x,y)
             * corresponds to a real coefficient of the DFT. */
            /* A random sign is drawn (except for (0,0)
             * where we impose sign=1. to conserve the
             * original mean) */
            if (((x == 0) || ((xev == 0) && (x == hinx)))
                && ((y == 0) || ((yev == 0) && (y == hiny)))) {
                sign = (((mt_drand53() < 0.5)
                         || ((x == 0) && (y == 0))) ? 1. : -1.);
                (*comp_ptr)[0] = sign * invN;
                (*comp_ptr)[1] = 0.;
                comp_ptr++;
            }
            /* case 2: Both (x,y) and its symmetric are in
             * the same half-domain, and y > ny/2. */
            /* Then the random phase of this symmetric
             * point has already been drawn. */
            else if (((x == 0) || ((xev == 0) && (x == hinx)))
                     && (y > hiny)) {
                /* copy the symmetric point */
                ad = (iny - y) * sx + x;        /* adress of the
                                                 * symmetric point */
                (*comp_ptr)[0] = data_out[ad][0];
                (*comp_ptr)[1] = data_out[ad][1];
                comp_ptr++;
            }
            else {
                /* draw a random phase */
                theta = (2. * mt_drand53() - 1.) * M_PI;
                costheta = (float) cos(theta);
                sintheta = (float) sin(theta);
                (*comp_ptr)[0] = costheta * invN;
                (*comp_ptr)[1] = sintheta * invN;
                comp_ptr++;
            }
        }
    }
    return (data_out);
}

static float **phase_randomization(float **data_in_rgb, fftwf_complex * rp,
                                   size_t nx, size_t ny)
{
    fftwf_plan plan_r2c, plan_c2r;
    fftwf_complex *fft;
    size_t fft_size = (nx / 2 + 1) * ny;        /* physical size of the r2c FFT */
    int channel;                /* RGB channel index */

    /* Check allocation */
    if (NULL == data_in_rgb || NULL == rp)
        return (NULL);
    /* Start threaded fftw if FFTW_NTHREADS is defined */
#ifdef FFTW_NTHREADS
    if (0 == fftwf_init_threads())
        return (NULL);
    fftwf_plan_with_nthreads(FFTW_NTHREADS);
#endif

    /* Memory allocation */
    if (NULL == (fft = (fftwf_complex *)
                 fftwf_malloc(fft_size * sizeof(fftwf_complex))))
        return (NULL);

    /* For each channel */
    for (channel = 0; channel < 3; channel++) {
        /* Forward Fourier transform : create the FFT forward
         * plan and run the FFT */
        plan_r2c = fftwf_plan_dft_r2c_2d((int) (ny), (int) (nx),
                                         data_in_rgb[channel], fft,
                                         FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
        fftwf_execute(plan_r2c);

        /* Phase randomization */
        if (NULL == pointwise_complex_multiplication(fft, rp, fft_size))
            return (NULL);

        /* Backward Fourier transform */
        plan_c2r = fftwf_plan_dft_c2r_2d((int) (ny), (int) (nx),
                                         fft, data_in_rgb[channel],
                                         FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
        fftwf_execute(plan_c2r);
    }

    /* Free memory : destroy the FFT plans and data */
    fftwf_destroy_plan(plan_r2c);
    fftwf_destroy_plan(plan_c2r);
    fftwf_free(fft);
    fftwf_cleanup();
#ifdef FFTW_NTHREADS
    fftwf_cleanup_threads();
#endif

    return (data_in_rgb);
}

/*
 * Spot extension section
 */

static float *compute_transition(float *transition,
                                 size_t nxin, size_t nyin, float alpha)
{
    float *ptr_trans;           /* pointer towards transition */
    float *bumpx, *bumpy;       /* transition part along
                                 * x-axis and y-axis */
    float *ptr_bumpx, *ptr_bumpy;       /* pointers towards bumpx and bumpy */
    int margx, margy;           /* width of margins */
    int x, y;                   /* indices for for loops */
    float tmp, tmpy, tmpx;      /* temp values */

    /* Check allocation */
    if (NULL == transition)
        return (NULL);

    /* Sizes for memory allocation */
    margx = ceil(alpha * nxin);
    margy = ceil(alpha * nyin);

    /* Memory allocation */
    if (NULL == (bumpx = (float *) malloc(margx * sizeof(float)))
        || NULL == (bumpy = (float *) malloc(margy * sizeof(float))))
        return (NULL);

    /* Computation of the radial transition bumpx: it is the
     * cumulative histogram of the molifier */
    ptr_bumpx = bumpx;
    tmp = 0.;
    for (x = 0; x < margx; x++) {
        tmpx = 2 * x / (margx - 1) - 1; /* affine change of variable
                                         * [0,margx-1] -> [-1,1] */
        tmp += exp(-1. / (1. - tmpx * tmpx));
        *ptr_bumpx = tmp;
        ptr_bumpx++;
    }
    /* Normalization of the cumulative histogram: the total sum
     * must equal 1 */
    ptr_bumpx = bumpx;
    for (x = 0; x < margx; x++) {
        *ptr_bumpx /= tmp;
        ptr_bumpx++;
    }

    /* Same procedure for bumpy */
    ptr_bumpy = bumpy;
    tmp = 0.;
    for (y = 0; y < margy; y++) {
        tmpy = 2 * y / (margy - 1) - 1;
        tmp += exp(-1. / (1. - tmpy * tmpy));
        *ptr_bumpy = tmp;
        ptr_bumpy++;
    }
    ptr_bumpy = bumpy;
    for (y = 0; y < margy; y++) {
        *ptr_bumpy /= tmp;
        ptr_bumpy++;
    }

    /* Computation of the float image transition */
    ptr_trans = transition;
    ptr_bumpy = bumpy;
    ptr_bumpx = bumpx;

    for (y = 0; y < (int) nyin; y++) {
        /* Test to determine the atenuation coefficient tmpy:
         * we go back and forth on the array bumpy */
        if (y < margy) {
            tmpy = *ptr_bumpy;
            ptr_bumpy++;
        }
        else if (y >= (int) nyin - margy) {
            ptr_bumpy--;
            tmpy = *ptr_bumpy;
        }
        else {
            tmpy = 1.;
        }
        for (x = 0; x < (int) nxin; x++) {
            /* Test to determine the atenuation coefficient
             * tmpx: we go back and forth on the array bumpx */
            if (x < margx) {
                tmpx = *ptr_bumpx;
                ptr_bumpx++;
            }
            else if (x >= (int) nxin - margx) {
                ptr_bumpx--;
                tmpx = *ptr_bumpx;
            }
            else {
                tmpx = 1.;
            }
            /* Computation of the coefficient of transition */
            *ptr_trans = tmpy * tmpx;
            ptr_trans++;
        }
    }

    /* Free memory */
    free(bumpx);
    free(bumpy);

    /* Return */
    return (transition);
}

static float **spot_extension_color(float **data_out_rgb, float **data_in_rgb,
                                    size_t nxout, size_t nyout,
                                    size_t nxin, size_t nyin, float alpha)
{
    int channel;
    float *float_ptr_in;        /* pointers towards the
                                 * channels of data_in_rgb */
    float *float_ptr_out;       /* pointers towards the
                                 * channels of data_out_rgb */
    float *float_ptr_out_end, *float_ptr_out_beg;       /* pointers towards the
                                                         * channels of data_out_rgb */
    float *transition;          /* float array for the
                                 * transition function */
    float *ptr_trans, *ptr_trans_end;   /* pointers towards transition */
    float sqsum;                /* sum of the square of the
                                 * coefficient of transition */
    float normf;                /* normalization coefficient */
    float m;                    /* mean of the current color
                                 * channel */
    int lx, ly;                 /* coordinates of the
                                 * top-left corner of the
                                 * inner frame */
    int j;                      /* index for for loops */

    /* Check allocation */
    if (NULL == data_out_rgb || NULL == data_in_rgb)
        return (NULL);

    /* Memory allocation */
    if (NULL == (transition = (float *) malloc(nxin * nyin * sizeof(float))))
        return (NULL);

    /* Compute the transition */
    if (NULL == (transition = compute_transition(transition,
                                                 nxin, nyin, alpha)))
        return (NULL);

    /* Computation of the sum of the square of the coefficient of
     * transition */
    ptr_trans = transition;
    ptr_trans_end = ptr_trans + nxin * nyin;
    sqsum = 0.;
    while (ptr_trans < ptr_trans_end) {
        sqsum += (*ptr_trans) * (*ptr_trans);
        ptr_trans++;
    }
    /* Normalization so that the extended spot has roughly the
     * same variance than the original one */
    normf = sqrt(((float) nxout * nyout) / sqsum);
    ptr_trans = transition;
    while (ptr_trans < ptr_trans_end) {
        *ptr_trans *= normf;
        ptr_trans++;
    }
    /* Copy the normalized and smoothed original spot in the
     * middle of data_out_rgb */
    for (channel = 0; channel < 3; channel++) {
        /* Compute the mean of data_in_rgb[channel] */
        m = mean(data_in_rgb[channel], nxin * nyin);
        /* fill data_out_rgb[channel] with the mean value */
        float_ptr_out = data_out_rgb[channel];
        float_ptr_out_end = float_ptr_out + nxout * nyout;
        while (float_ptr_out < float_ptr_out_end) {
            *float_ptr_out = m;
            float_ptr_out++;
        }
        /* Copy in the middle */
        lx = (int) (0.5 * ((float) (nxout - nxin)));
        ly = (int) (0.5 * ((float) (nyout - nyin)));
        float_ptr_out = data_out_rgb[channel] + ly * nxout + lx;
        float_ptr_in = data_in_rgb[channel];
        ptr_trans = transition;
        for (j = 0; j < (int) nyin; j++) {
            float_ptr_out_beg = float_ptr_out;
            float_ptr_out_end = float_ptr_out + nxin;
            while (float_ptr_out < float_ptr_out_end) {
                *float_ptr_out += (*float_ptr_in - m)
                    * (*ptr_trans);
                float_ptr_out++;
                float_ptr_in++;
                ptr_trans++;
            }
            float_ptr_out = float_ptr_out_beg + nxout;
        }
    }

    /* Free memory */
    free(transition);

    /* Return */
    return (data_out_rgb);
}

/*
* Random Phase Noise section
*/

float **rpn_color_different_size(float **data_out_rgb, float **data_in_rgb,
                                 size_t nxout, size_t nyout,
                                 size_t nxin, size_t nyin, float alpha)
{
    fftwf_complex *rp;          /* random phase */
    size_t rp_size = (nxout / 2 + 1) * nyout;   /* size of the random
                                                 * phase = physical size
                                                 * of the r2c FFT */

    /* Check allocation */
    if (NULL == data_in_rgb || NULL == data_out_rgb)
        return (NULL);

    /* Allocate memory */
    if (NULL == (rp = (fftwf_complex *)
                 fftwf_malloc(rp_size * sizeof(fftwf_complex))))
        return (NULL);

    /* compute the periodic component of data_in_rgb */
    if (NULL == periodic_component_color(data_in_rgb, nxin, nyin))
        return (NULL);

    /* compute the extended spot of the periodic component and
     * store it in data_out_rgb */
    if (NULL == spot_extension_color(data_out_rgb, data_in_rgb,
                                     nxout, nyout, nxin, nyin, alpha))
        return (NULL);

    /* compute a random phase */
    if (NULL == random_phase(rp, nxout, nyout))
        return (NULL);

    /* phase randomization */
    if (NULL == phase_randomization(data_out_rgb, rp, nxout, nyout))
        return (NULL);

    /* free memory */
    fftwf_free(rp);

    return (data_out_rgb);
}

static fftwf_complex
    * random_phase_with_poisson_complex_filter(fftwf_complex * data_out,
                                               size_t nx, size_t ny)
{
    int inx, iny;               /* int version of nx and ny */
    int hinx, hiny;             /* half of nx and ny */
    int xev, yev;               /* 0 if x is even, 1 otherwise */
    fftwf_complex *comp_ptr;    /* complex pointer */
    int sx = (nx / 2 + 1);      /* physical x-size of data_out */
    int x, y;                   /* index for for loops */
    int ad;                     /* pixels locations */
    double theta;               /* random phase */
    float costheta, sintheta;   /* cosine and sine of theta */
    float sign;                 /* random sign */
    float *pcf;                 /* Poisson complex filter */
    float *pcf_ptr;
    size_t fft_size = (nx / 2 + 1) * ny;        /* physical size of the r2c FFT */

    /* Check allocation */
    if (NULL == data_out)
        return (NULL);

    /* Allocate memory */
    if (NULL == (pcf = (float *) malloc(fft_size * sizeof(float))))
        return (NULL);

    /* Compute the Poisson complex filter */
    if (NULL == poisson_complex_filter(pcf, nx, ny))
        return (NULL);

    /* Define sizes and indexes for for loops: depend of the parity */
    inx = (int) nx;
    iny = (int) ny;
    xev = inx % 2;
    yev = iny % 2;
    hinx = inx / 2;
    hiny = iny / 2;

    /* Pointers initialization */
    comp_ptr = data_out;
    pcf_ptr = pcf;

    for (y = 0; y < iny; y++) {
        for (x = 0; x < sx; x++) {
            /* Case 1: (x,y) is its own symmetric: (x,y)
             * corresponds to a real coefficient of the DFT. */
            /* A random sign is drawn (except for (0,0)
             * where we impose sign=1. to conserve the
             * original mean) */
            if (((x == 0) || ((xev == 0) && (x == hinx)))
                && ((y == 0) || ((yev == 0) && (y == hiny)))) {
                sign = (((mt_drand53() < 0.5)
                         || ((x == 0) && (y == 0))) ? 1. : -1.);
                (*comp_ptr)[0] = sign * (*pcf_ptr);
                (*comp_ptr)[1] = 0.;
                comp_ptr++;
                pcf_ptr++;
            }
            /* Case 2: Both (x,y) and its symmetric are in
             * the same half-domain, and y > ny/2. */
            /* Then the random phase of this symmetric
             * point has already been drawn. */
            else if (((x == 0) || ((xev == 0) && (x == hinx)))
                     && (y > hiny)) {
                /* Copy the symmetric point */
                ad = (iny - y) * sx + x;        /* adress of the
                                                 * symmetric point */
                (*comp_ptr)[0] = data_out[ad][0];
                (*comp_ptr)[1] = data_out[ad][1];
                comp_ptr++;
                pcf_ptr++;
            }
            else {
                /* Draw a random phase */
                theta = (2. * mt_drand53() - 1.) * M_PI;
                costheta = (float) cos(theta);
                sintheta = (float) sin(theta);
                (*comp_ptr)[0] = costheta * (*pcf_ptr);
                (*comp_ptr)[1] = sintheta * (*pcf_ptr);
                comp_ptr++;
                pcf_ptr++;
            }
        }
    }
    /* Free memory */
    free(pcf);

    return (data_out);
}

float **rpn_color_same_size(float **data_in_rgb, size_t nx, size_t ny)
{
    float *laplacian;           /* 2D discrete Laplacian */
    fftwf_complex *rp_pcf;      /* random phase with poisson
                                 * complex filter */
    fftwf_plan plan_r2c, plan_c2r;
    fftwf_complex *fft;
    float m;                    /* mean of the coefficient of
                                 * data_in */
    int channel;
    size_t N = nx * ny;
    size_t fft_size = (nx / 2 + 1) * ny;        /* physical size of the r2c FFT */

    /* Check allocation */
    if (NULL == data_in_rgb)
        return (NULL);
    /* Start threaded fftw if FFTW_NTHREADS is defined */
#ifdef FFTW_NTHREADS
    if (0 == fftwf_init_threads())
        return (NULL);
    fftwf_plan_with_nthreads(FFTW_NTHREADS);
#endif
    /* Allocate memory */
    if (NULL == (laplacian = (float *) malloc(N * sizeof(float)))
        || NULL == (rp_pcf = (fftwf_complex *)
                    fftwf_malloc(fft_size * sizeof(fftwf_complex)))
        || NULL == (fft = (fftwf_complex *)
                    fftwf_malloc(fft_size * sizeof(fftwf_complex))))
        return (NULL);

    /* Compute the random phase with poisson complex filter */
    if (NULL == random_phase_with_poisson_complex_filter(rp_pcf, nx, ny))
        return (NULL);

    /* For each channel */
    for (channel = 0; channel < 3; channel++) {
        /* Compute m = the mean value of data_in[channel] */
        m = mean(data_in_rgb[channel], N);

        /* Compute the discrete Laplacian of data_in[channel] */
        if (NULL == discrete_laplacian(laplacian,
                                       data_in_rgb[channel], nx, ny))
            return (NULL);

        /* Forward Fourier transform of the Laplacian: create
         * the FFT forward plan and run the FFT */
        plan_r2c = fftwf_plan_dft_r2c_2d((int) (ny), (int) (nx),
                                         laplacian, fft,
                                         FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
        fftwf_execute(plan_r2c);

        /* Inverse the discrete periodic Laplacian in the
         * Fourier domain using the Poisson complex filter, and
         * randomize the phase */
        if (NULL == pointwise_complex_multiplication(fft, rp_pcf, fft_size))
            return (NULL);

        /* (0,0) frequency : we impose the same mean as
         * data_in[channel] */
        (*fft)[0] = m;
        (*fft)[1] = 0.;

        /* Backward Fourier transform: the output is stored in
         * data_in_rgb[channel] */
        plan_c2r = fftwf_plan_dft_c2r_2d((int) (ny), (int) (nx),
                                         fft, data_in_rgb[channel],
                                         FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
        fftwf_execute(plan_c2r);
    }

    /* Cleanup */
    free(laplacian);

    /* Destroy the FFT plans and data */
    fftwf_destroy_plan(plan_r2c);
    fftwf_destroy_plan(plan_c2r);
    fftwf_free(fft);
    fftwf_free(rp_pcf);
    fftwf_cleanup();
#ifdef FFTW_NTHREADS
    fftwf_cleanup_threads();
#endif

    return (data_in_rgb);
}