- Main Page
- Related Pages
- Files

- File List
- Globals

# random_phase_noise_lib.c

Go to the documentation of this file.

    00001 /*
    00002  * Copyright (c) 2010, Bruno Galerne <bruno.galerne@cmla.ens-cachan.fr>
    00003  * All rights reserved.
    00004  *
    00005  * This program is free software: you can use, modify and/or
    00006  * redistribute it under the terms of the GNU General Public
    00007  * License as published by the Free Software Foundation, either
    00008  * version 3 of the License, or (at your option) any later
    00009  * version. You should have received a copy of this license along
    00010  * this program. If not, see <http://www.gnu.org/licenses/>.
    00011  */
    00012
    00026 #include <stdio.h>
    00027 #include <stdlib.h>
    00028 #include <math.h>
    00029 #include <fftw3.h>
    00030
    00031 #include "mt.h"
    00032
    00033 /* M_PI is a POSIX definition */
    00034 #ifndef M_PI
    00035
    00036 #define M_PI 3.14159265358979323846
    00037 #endif                          /* !M_PI */
    00038
    00039 /*
    00040 * number of threads to use for libfftw
    00041 * (uncomment to enable parallel FFT multi-threading)
    00042 */
    00043 /* #define FFTW_NTHREADS 4 */
    00044
    00045 /*
    00046 * Preliminary section: some useful standard operations
    00047 */
    00048
    00056 static float mean(float *data_in, size_t N)
    00057 {
    00058     float m = 0.;
    00059     float *ptr, *ptr_end;
    00060
    00061     /* Pointers initialization */
    00062     ptr = data_in;
    00063     ptr_end = ptr + N;
    00064
    00065     /* Sum loop */
    00066     while (ptr < ptr_end) {
    00067         m += *ptr;
    00068         ptr++;
    00069     }
    00070
    00071     /* Normalization */
    00072     m /= ((float) N);
    00073
    00074     return (m);
    00075 }
    00076
    00087 static fftwf_complex
    00088     * pointwise_complexfloat_multiplication(fftwf_complex * comp_in,
    00089                                             float *float_in, size_t N)
    00090 {
    00091     fftwf_complex *ptr_comp_in, *ptr_comp_end;
    00092     float *ptr_float_in;
    00093
    00094     /* check allocaton */
    00095     if (NULL == comp_in || NULL == float_in)
    00096         return (NULL);
    00097
    00098     ptr_comp_in = comp_in;
    00099     ptr_float_in = float_in;
    00100     ptr_comp_end = ptr_comp_in + N;
    00101     while (ptr_comp_in < ptr_comp_end) {
    00102         (*ptr_comp_in)[0] *= (*ptr_float_in);
    00103         (*ptr_comp_in)[1] *= (*ptr_float_in);
    00104         ptr_comp_in++;
    00105         ptr_float_in++;
    00106     }
    00107     return (comp_in);
    00108 }
    00109
    00120 static fftwf_complex
    00121     * pointwise_complex_multiplication(fftwf_complex * comp_in1,
    00122                                        fftwf_complex * comp_in2, size_t N)
    00123 {
    00124     fftwf_complex *ptr_comp_in1, *ptr_comp_in2, *ptr_comp_end;
    00125     float re1, im1;             /* real and imaginary part of *ptr_comp_in1 */
    00126
    00127     /* Check allocaton */
    00128     if (NULL == comp_in1 || NULL == comp_in2)
    00129         return (NULL);
    00130
    00131     ptr_comp_in1 = comp_in1;
    00132     ptr_comp_in2 = comp_in2;
    00133     ptr_comp_end = ptr_comp_in1 + N;
    00134     while (ptr_comp_in1 < ptr_comp_end) {
    00135         re1 = (*ptr_comp_in1)[0];
    00136         im1 = (*ptr_comp_in1)[1];
    00137         (*ptr_comp_in1)[0] = re1 * (*ptr_comp_in2)[0]
    00138             - im1 * (*ptr_comp_in2)[1];
    00139         (*ptr_comp_in1)[1] = re1 * (*ptr_comp_in2)[1]
    00140             + im1 * (*ptr_comp_in2)[0];
    00141         ptr_comp_in1++;
    00142         ptr_comp_in2++;
    00143     }
    00144     return (comp_in1);
    00145 }
    00146
    00147 /*
    00148 * Periodic component section
    00149 */
    00150
    00170 static float *discrete_laplacian(float *data_out, float *data_in,
    00171                                  size_t nx, size_t ny)
    00172 {
    00173     int i, j;
    00174     float *out_ptr;
    00175     float *in_ptr, *in_ptr_xm1, *in_ptr_xp1, *in_ptr_ym1, *in_ptr_yp1;
    00176
    00177     /* check allocation */
    00178     if (NULL == data_in || NULL == data_out)
    00179         return (NULL);
    00180
    00181     /* pointers to the data and neighbour values */
    00182     in_ptr = data_in;
    00183     in_ptr_xm1 = data_in - 1;
    00184     in_ptr_xp1 = data_in + 1;
    00185     in_ptr_ym1 = data_in - nx;
    00186     in_ptr_yp1 = data_in + nx;
    00187     out_ptr = data_out;
    00188     /* iterate on j, i, following the array order */
    00189     for (j = 0; j < (int) ny; j++) {
    00190         for (i = 0; i < (int) nx; i++) {
    00191             *out_ptr = 0.;
    00192             /* row differences */
    00193             if (0 < i) {
    00194                 *out_ptr += *(in_ptr_xm1) - *in_ptr;
    00195             }
    00196             if ((int) nx - 1 > i) {
    00197                 *out_ptr += *(in_ptr_xp1) - *in_ptr;
    00198             }
    00199             /* column differences */
    00200             if (0 < j) {
    00201                 *out_ptr += *(in_ptr_ym1) - *in_ptr;
    00202             }
    00203             if ((int) ny - 1 > j) {
    00204                 *out_ptr += *(in_ptr_yp1) - *in_ptr;
    00205             }
    00206             in_ptr++;
    00207             in_ptr_xm1++;
    00208             in_ptr_xp1++;
    00209             in_ptr_ym1++;
    00210             in_ptr_yp1++;
    00211             out_ptr++;
    00212         }
    00213     }
    00214     return (data_out);
    00215 }
    00216
    00242 static float *poisson_complex_filter(float *data_out, size_t nx, size_t ny)
    00243 {
    00244     float *cosx, *cosymintwo, *ptr_x, *ptr_y, *ptr_x_end,
    00245         *ptr_y_end, *out_ptr;
    00246     int i, j;
    00247     size_t sx = nx / 2 + 1;
    00248     int N = nx * ny;
    00249     float halfinvN = 0.5 / ((float) N);
    00250
    00251     /* check allocation */
    00252     if (NULL == data_out)
    00253         return (NULL);
    00254
    00255     /* allocate memory */
    00256     if (NULL == (cosx = (float *) malloc(sx * sizeof(float)))
    00257         || NULL == (cosymintwo = (float *) malloc(ny * sizeof(float))))
    00258         return (NULL);
    00259
    00260     /* define cosx and cosymintwo */
    00261     ptr_x = cosx;
    00262     for (i = 0; i < (int) sx; i++) {
    00263         (*ptr_x) = cos(2. * ((float) i) * M_PI / ((float) nx));
    00264         ptr_x++;
    00265     }
    00266     ptr_y = cosymintwo;
    00267     for (j = 0; j < (int) ny; j++) {
    00268         (*ptr_y) = cos(2. * ((float) j) * M_PI / ((float) ny)) - 2.;
    00269         ptr_y++;
    00270     }
    00271
    00272     /* define data_out */
    00273     out_ptr = data_out;
    00274     ptr_y = cosymintwo;
    00275     ptr_y_end = cosymintwo + ny;
    00276     ptr_x = cosx;
    00277     ptr_x_end = cosx + sx;
    00278     /* particular case : (0,0) */
    00279     (*out_ptr) = 1.;
    00280     out_ptr++;
    00281     ptr_x++;
    00282     while (ptr_y < ptr_y_end) {
    00283         while (ptr_x < ptr_x_end) {
    00284             (*out_ptr) = halfinvN / ((*ptr_x) + (*ptr_y));
    00285             out_ptr++;
    00286             ptr_x++;
    00287         }
    00288         ptr_x = cosx;
    00289         ptr_y++;
    00290     }
    00291
    00292     /* free memory */
    00293     free(cosx);
    00294     free(cosymintwo);
    00295
    00296     /* return data_out */
    00297     return (data_out);
    00298 }
    00299
    00332 static float **periodic_component_color(float **data_in_rgb,
    00333                                         size_t nx, size_t ny)
    00334 {
    00335     float *laplacian;           /* 2D discrete Laplacian */
    00336     float *pcf;                 /* Poisson complex filter */
    00337     fftwf_plan plan_r2c, plan_c2r;
    00338     fftwf_complex *fft;
    00339     int N = ((int) nx) * ((int) ny);
    00340     float m;                    /* mean of the coefficient of data_in */
    00341     int channel;
    00342     size_t fft_size = (nx / 2 + 1) * ny;        /* physical size of the r2c FFT */
    00343
    00344     /* Check allocation */
    00345     if (NULL == data_in_rgb)
    00346         return (NULL);
    00347
    00348     /* Start threaded fftw if FFTW_NTHREADS is defined */
    00349 #ifdef FFTW_NTHREADS
    00350     if (0 == fftwf_init_threads())
    00351         return (NULL);
    00352     fftwf_plan_with_nthreads(FFTW_NTHREADS);
    00353 #endif
    00354
    00355     /* Allocate memory */
    00356     if (NULL == (pcf = (float *) malloc(fft_size * sizeof(float)))
    00357         || NULL == (laplacian = (float *) malloc(N * sizeof(float)))
    00358         || NULL == (fft = (fftwf_complex *)
    00359                     fftwf_malloc(fft_size * sizeof(fftwf_complex))))
    00360         return (NULL);
    00361
    00362     /* Compute the Poisson complex filter */
    00363     if (NULL == poisson_complex_filter(pcf, nx, ny))
    00364         return (NULL);
    00365
    00366     /* For each channel */
    00367     for (channel = 0; channel < 3; channel++) {
    00368         /* Compute m = the mean value of data_in[channel] */
    00369         m = mean(data_in_rgb[channel], N);
    00370
    00371         /* Compute the discrete Laplacian of data_in[channel] */
    00372         if (NULL == discrete_laplacian(laplacian,
    00373                                        data_in_rgb[channel], nx, ny))
    00374             return (NULL);
    00375
    00376         /* Forward Fourier transform of the Laplacian: create
    00377          * the FFT forward plan and run the FFT */
    00378         plan_r2c = fftwf_plan_dft_r2c_2d((int) (ny), (int) (nx),
    00379                                          laplacian, fft,
    00380                                          FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
    00381         fftwf_execute(plan_r2c);
    00382
    00383         /* Inverse the discrete periodic Laplacian in the
    00384          * Fourier domain using the Poisson complex filter */
    00385         if (NULL == pointwise_complexfloat_multiplication(fft, pcf, fft_size))
    00386             return (NULL);
    00387         /* (0,0) frequency : we impose the same mean as the input */
    00388         (*fft)[0] = m;
    00389         (*fft)[1] = 0.;
    00390
    00391         /* Backward Fourier transform: the output is stored in
    00392          * data_in_rgb[channel] */
    00393         plan_c2r = fftwf_plan_dft_c2r_2d((int) (ny), (int) (nx),
    00394                                          fft, data_in_rgb[channel],
    00395                                          FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
    00396         fftwf_execute(plan_c2r);
    00397     }
    00398
    00399     /* cleanup */
    00400     free(pcf);
    00401     free(laplacian);
    00402     /* destroy the FFT plans and data */
    00403     fftwf_destroy_plan(plan_r2c);
    00404     fftwf_destroy_plan(plan_c2r);
    00405     fftwf_free(fft);
    00406     fftwf_cleanup();
    00407 #ifdef FFTW_NTHREADS
    00408     fftwf_cleanup_threads();
    00409 #endif
    00410
    00411     return (data_in_rgb);
    00412 }
    00413
    00414 /*
    00415 * Phase randomization section
    00416 */
    00417
    00433 static fftwf_complex *random_phase(fftwf_complex * data_out, size_t nx,
    00434                                    size_t ny)
    00435 {
    00436     int inx, iny;               /* int version of nx and ny */
    00437     int hinx, hiny;             /* half of nx and ny */
    00438     float invN;                 /* inverse of nx*ny */
    00439     int xev, yev;               /* 0 if x is even, 1 otherwise */
    00440     fftwf_complex *comp_ptr;    /* complex pointer */
    00441     int sx = (nx / 2 + 1);      /* physical x-size of data_out */
    00442     int x, y;                   /* index for for loops */
    00443     int ad;                     /* pixels locations */
    00444     double theta;               /* random phase */
    00445     float costheta, sintheta;   /* cosine and sine of theta */
    00446     float sign;                 /* random sign */
    00447
    00448     /* Check allocation */
    00449     if (NULL == data_out)
    00450         return (NULL);
    00451
    00452     /* Define sizes and indexes for for loops: depend of the parity */
    00453     inx = (int) nx;
    00454     iny = (int) ny;
    00455     invN = 1. / ((float) (inx * iny));
    00456     xev = inx % 2;
    00457     yev = iny % 2;
    00458     hinx = inx / 2;
    00459     hiny = iny / 2;
    00460
    00461     /* Pointer initialization */
    00462     comp_ptr = data_out;
    00463
    00464     for (y = 0; y < iny; y++) {
    00465         for (x = 0; x < sx; x++) {
    00466             /* Case 1: (x,y) is its own symmetric: (x,y)
    00467              * corresponds to a real coefficient of the DFT. */
    00468             /* A random sign is drawn (except for (0,0)
    00469              * where we impose sign=1. to conserve the
    00470              * original mean) */
    00471             if (((x == 0) || ((xev == 0) && (x == hinx)))
    00472                 && ((y == 0) || ((yev == 0) && (y == hiny)))) {
    00473                 sign = (((mt_drand53() < 0.5)
    00474                          || ((x == 0) && (y == 0))) ? 1. : -1.);
    00475                 (*comp_ptr)[0] = sign * invN;
    00476                 (*comp_ptr)[1] = 0.;
    00477                 comp_ptr++;
    00478             }
    00479             /* case 2: Both (x,y) and its symmetric are in
    00480              * the same half-domain, and y > ny/2. */
    00481             /* Then the random phase of this symmetric
    00482              * point has already been drawn. */
    00483             else if (((x == 0) || ((xev == 0) && (x == hinx)))
    00484                      && (y > hiny)) {
    00485                 /* copy the symmetric point */
    00486                 ad = (iny - y) * sx + x;        /* adress of the
    00487                                                  * symmetric point */
    00488                 (*comp_ptr)[0] = data_out[ad][0];
    00489                 (*comp_ptr)[1] = data_out[ad][1];
    00490                 comp_ptr++;
    00491             }
    00492             else {
    00493                 /* draw a random phase */
    00494                 theta = (2. * mt_drand53() - 1.) * M_PI;
    00495                 costheta = (float) cos(theta);
    00496                 sintheta = (float) sin(theta);
    00497                 (*comp_ptr)[0] = costheta * invN;
    00498                 (*comp_ptr)[1] = sintheta * invN;
    00499                 comp_ptr++;
    00500             }
    00501         }
    00502     }
    00503     return (data_out);
    00504 }
    00505
    00518 static float **phase_randomization(float **data_in_rgb, fftwf_complex * rp,
    00519                                    size_t nx, size_t ny)
    00520 {
    00521     fftwf_plan plan_r2c, plan_c2r;
    00522     fftwf_complex *fft;
    00523     size_t fft_size = (nx / 2 + 1) * ny;        /* physical size of the r2c FFT */
    00524     int channel;                /* RGB channel index */
    00525
    00526     /* Check allocation */
    00527     if (NULL == data_in_rgb || NULL == rp)
    00528         return (NULL);
    00529     /* Start threaded fftw if FFTW_NTHREADS is defined */
    00530 #ifdef FFTW_NTHREADS
    00531     if (0 == fftwf_init_threads())
    00532         return (NULL);
    00533     fftwf_plan_with_nthreads(FFTW_NTHREADS);
    00534 #endif
    00535
    00536     /* Memory allocation */
    00537     if (NULL == (fft = (fftwf_complex *)
    00538                  fftwf_malloc(fft_size * sizeof(fftwf_complex))))
    00539         return (NULL);
    00540
    00541     /* For each channel */
    00542     for (channel = 0; channel < 3; channel++) {
    00543         /* Forward Fourier transform : create the FFT forward
    00544          * plan and run the FFT */
    00545         plan_r2c = fftwf_plan_dft_r2c_2d((int) (ny), (int) (nx),
    00546                                          data_in_rgb[channel], fft,
    00547                                          FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
    00548         fftwf_execute(plan_r2c);
    00549
    00550         /* Phase randomization */
    00551         if (NULL == pointwise_complex_multiplication(fft, rp, fft_size))
    00552             return (NULL);
    00553
    00554         /* Backward Fourier transform */
    00555         plan_c2r = fftwf_plan_dft_c2r_2d((int) (ny), (int) (nx),
    00556                                          fft, data_in_rgb[channel],
    00557                                          FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
    00558         fftwf_execute(plan_c2r);
    00559     }
    00560
    00561     /* Free memory : destroy the FFT plans and data */
    00562     fftwf_destroy_plan(plan_r2c);
    00563     fftwf_destroy_plan(plan_c2r);
    00564     fftwf_free(fft);
    00565     fftwf_cleanup();
    00566 #ifdef FFTW_NTHREADS
    00567     fftwf_cleanup_threads();
    00568 #endif
    00569
    00570     return (data_in_rgb);
    00571 }
    00572
    00573 /*
    00574  * Spot extension section
    00575  */
    00576
    00603 static float *compute_transition(float *transition,
    00604                                  size_t nxin, size_t nyin, float alpha)
    00605 {
    00606     float *ptr_trans;           /* pointer towards transition */
    00607     float *bumpx, *bumpy;       /* transition part along
    00608                                  * x-axis and y-axis */
    00609     float *ptr_bumpx, *ptr_bumpy;       /* pointers towards bumpx and bumpy */
    00610     int margx, margy;           /* width of margins */
    00611     int x, y;                   /* indices for for loops */
    00612     float tmp, tmpy, tmpx;      /* temp values */
    00613
    00614     /* Check allocation */
    00615     if (NULL == transition)
    00616         return (NULL);
    00617
    00618     /* Sizes for memory allocation */
    00619     margx = ceil(alpha * nxin);
    00620     margy = ceil(alpha * nyin);
    00621
    00622     /* Memory allocation */
    00623     if (NULL == (bumpx = (float *) malloc(margx * sizeof(float)))
    00624         || NULL == (bumpy = (float *) malloc(margy * sizeof(float))))
    00625         return (NULL);
    00626
    00627     /* Computation of the radial transition bumpx: it is the
    00628      * cumulative histogram of the molifier */
    00629     ptr_bumpx = bumpx;
    00630     tmp = 0.;
    00631     for (x = 0; x < margx; x++) {
    00632         tmpx = 2 * x / (margx - 1) - 1; /* affine change of variable
    00633                                          * [0,margx-1] -> [-1,1] */
    00634         tmp += exp(-1. / (1. - tmpx * tmpx));
    00635         *ptr_bumpx = tmp;
    00636         ptr_bumpx++;
    00637     }
    00638     /* Normalization of the cumulative histogram: the total sum
    00639      * must equal 1 */
    00640     ptr_bumpx = bumpx;
    00641     for (x = 0; x < margx; x++) {
    00642         *ptr_bumpx /= tmp;
    00643         ptr_bumpx++;
    00644     }
    00645
    00646     /* Same procedure for bumpy */
    00647     ptr_bumpy = bumpy;
    00648     tmp = 0.;
    00649     for (y = 0; y < margy; y++) {
    00650         tmpy = 2 * y / (margy - 1) - 1;
    00651         tmp += exp(-1. / (1. - tmpy * tmpy));
    00652         *ptr_bumpy = tmp;
    00653         ptr_bumpy++;
    00654     }
    00655     ptr_bumpy = bumpy;
    00656     for (y = 0; y < margy; y++) {
    00657         *ptr_bumpy /= tmp;
    00658         ptr_bumpy++;
    00659     }
    00660
    00661     /* Computation of the float image transition */
    00662     ptr_trans = transition;
    00663     ptr_bumpy = bumpy;
    00664     ptr_bumpx = bumpx;
    00665
    00666     for (y = 0; y < (int) nyin; y++) {
    00667         /* Test to determine the atenuation coefficient tmpy:
    00668          * we go back and forth on the array bumpy */
    00669         if (y < margy) {
    00670             tmpy = *ptr_bumpy;
    00671             ptr_bumpy++;
    00672         }
    00673         else if (y >= (int) nyin - margy) {
    00674             ptr_bumpy--;
    00675             tmpy = *ptr_bumpy;
    00676         }
    00677         else {
    00678             tmpy = 1.;
    00679         }
    00680         for (x = 0; x < (int) nxin; x++) {
    00681             /* Test to determine the atenuation coefficient
    00682              * tmpx: we go back and forth on the array bumpx */
    00683             if (x < margx) {
    00684                 tmpx = *ptr_bumpx;
    00685                 ptr_bumpx++;
    00686             }
    00687             else if (x >= (int) nxin - margx) {
    00688                 ptr_bumpx--;
    00689                 tmpx = *ptr_bumpx;
    00690             }
    00691             else {
    00692                 tmpx = 1.;
    00693             }
    00694             /* Computation of the coefficient of transition */
    00695             *ptr_trans = tmpy * tmpx;
    00696             ptr_trans++;
    00697         }
    00698     }
    00699
    00700     /* Free memory */
    00701     free(bumpx);
    00702     free(bumpy);
    00703
    00704     /* Return */
    00705     return (transition);
    00706 }
    00707
    00729 static float **spot_extension_color(float **data_out_rgb, float **data_in_rgb,
    00730                                     size_t nxout, size_t nyout,
    00731                                     size_t nxin, size_t nyin, float alpha)
    00732 {
    00733     int channel;
    00734     float *float_ptr_in;        /* pointers towards the
    00735                                  * channels of data_in_rgb */
    00736     float *float_ptr_out;       /* pointers towards the
    00737                                  * channels of data_out_rgb */
    00738     float *float_ptr_out_end, *float_ptr_out_beg;       /* pointers towards the
    00739                                                          * channels of data_out_rgb */
    00740     float *transition;          /* float array for the
    00741                                  * transition function */
    00742     float *ptr_trans, *ptr_trans_end;   /* pointers towards transition */
    00743     float sqsum;                /* sum of the square of the
    00744                                  * coefficient of transition */
    00745     float normf;                /* normalization coefficient */
    00746     float m;                    /* mean of the current color
    00747                                  * channel */
    00748     int lx, ly;                 /* coordinates of the
    00749                                  * top-left corner of the
    00750                                  * inner frame */
    00751     int j;                      /* index for for loops */
    00752
    00753     /* Check allocation */
    00754     if (NULL == data_out_rgb || NULL == data_in_rgb)
    00755         return (NULL);
    00756
    00757     /* Memory allocation */
    00758     if (NULL == (transition = (float *) malloc(nxin * nyin * sizeof(float))))
    00759         return (NULL);
    00760
    00761     /* Compute the transition */
    00762     if (NULL == (transition = compute_transition(transition,
    00763                                                  nxin, nyin, alpha)))
    00764         return (NULL);
    00765
    00766     /* Computation of the sum of the square of the coefficient of
    00767      * transition */
    00768     ptr_trans = transition;
    00769     ptr_trans_end = ptr_trans + nxin * nyin;
    00770     sqsum = 0.;
    00771     while (ptr_trans < ptr_trans_end) {
    00772         sqsum += (*ptr_trans) * (*ptr_trans);
    00773         ptr_trans++;
    00774     }
    00775     /* Normalization so that the extended spot has roughly the
    00776      * same variance than the original one */
    00777     normf = sqrt(((float) nxout * nyout) / sqsum);
    00778     ptr_trans = transition;
    00779     while (ptr_trans < ptr_trans_end) {
    00780         *ptr_trans *= normf;
    00781         ptr_trans++;
    00782     }
    00783     /* Copy the normalized and smoothed original spot in the
    00784      * middle of data_out_rgb */
    00785     for (channel = 0; channel < 3; channel++) {
    00786         /* Compute the mean of data_in_rgb[channel] */
    00787         m = mean(data_in_rgb[channel], nxin * nyin);
    00788         /* fill data_out_rgb[channel] with the mean value */
    00789         float_ptr_out = data_out_rgb[channel];
    00790         float_ptr_out_end = float_ptr_out + nxout * nyout;
    00791         while (float_ptr_out < float_ptr_out_end) {
    00792             *float_ptr_out = m;
    00793             float_ptr_out++;
    00794         }
    00795         /* Copy in the middle */
    00796         lx = (int) (0.5 * ((float) (nxout - nxin)));
    00797         ly = (int) (0.5 * ((float) (nyout - nyin)));
    00798         float_ptr_out = data_out_rgb[channel] + ly * nxout + lx;
    00799         float_ptr_in = data_in_rgb[channel];
    00800         ptr_trans = transition;
    00801         for (j = 0; j < (int) nyin; j++) {
    00802             float_ptr_out_beg = float_ptr_out;
    00803             float_ptr_out_end = float_ptr_out + nxin;
    00804             while (float_ptr_out < float_ptr_out_end) {
    00805                 *float_ptr_out += (*float_ptr_in - m)
    00806                     * (*ptr_trans);
    00807                 float_ptr_out++;
    00808                 float_ptr_in++;
    00809                 ptr_trans++;
    00810             }
    00811             float_ptr_out = float_ptr_out_beg + nxout;
    00812         }
    00813     }
    00814
    00815     /* Free memory */
    00816     free(transition);
    00817
    00818     /* Return */
    00819     return (data_out_rgb);
    00820 }
    00821
    00822 /*
    00823 * Random Phase Noise section
    00824 */
    00825
    00850 float **rpn_color_different_size(float **data_out_rgb, float **data_in_rgb,
    00851                                  size_t nxout, size_t nyout,
    00852                                  size_t nxin, size_t nyin, float alpha)
    00853 {
    00854     fftwf_complex *rp;          /* random phase */
    00855     size_t rp_size = (nxout / 2 + 1) * nyout;   /* size of the random
    00856                                                  * phase = physical size
    00857                                                  * of the r2c FFT */
    00858
    00859     /* Check allocation */
    00860     if (NULL == data_in_rgb || NULL == data_out_rgb)
    00861         return (NULL);
    00862
    00863     /* Allocate memory */
    00864     if (NULL == (rp = (fftwf_complex *)
    00865                  fftwf_malloc(rp_size * sizeof(fftwf_complex))))
    00866         return (NULL);
    00867
    00868     /* compute the periodic component of data_in_rgb */
    00869     if (NULL == periodic_component_color(data_in_rgb, nxin, nyin))
    00870         return (NULL);
    00871
    00872     /* compute the extended spot of the periodic component and
    00873      * store it in data_out_rgb */
    00874     if (NULL == spot_extension_color(data_out_rgb, data_in_rgb,
    00875                                      nxout, nyout, nxin, nyin, alpha))
    00876         return (NULL);
    00877
    00878     /* compute a random phase */
    00879     if (NULL == random_phase(rp, nxout, nyout))
    00880         return (NULL);
    00881
    00882     /* phase randomization */
    00883     if (NULL == phase_randomization(data_out_rgb, rp, nxout, nyout))
    00884         return (NULL);
    00885
    00886     /* free memory */
    00887     fftwf_free(rp);
    00888
    00889     return (data_out_rgb);
    00890 }
    00891
    00903 static fftwf_complex
    00904     * random_phase_with_poisson_complex_filter(fftwf_complex * data_out,
    00905                                                size_t nx, size_t ny)
    00906 {
    00907     int inx, iny;               /* int version of nx and ny */
    00908     int hinx, hiny;             /* half of nx and ny */
    00909     int xev, yev;               /* 0 if x is even, 1 otherwise */
    00910     fftwf_complex *comp_ptr;    /* complex pointer */
    00911     int sx = (nx / 2 + 1);      /* physical x-size of data_out */
    00912     int x, y;                   /* index for for loops */
    00913     int ad;                     /* pixels locations */
    00914     double theta;               /* random phase */
    00915     float costheta, sintheta;   /* cosine and sine of theta */
    00916     float sign;                 /* random sign */
    00917     float *pcf;                 /* Poisson complex filter */
    00918     float *pcf_ptr;
    00919     size_t fft_size = (nx / 2 + 1) * ny;        /* physical size of the r2c FFT */
    00920
    00921     /* Check allocation */
    00922     if (NULL == data_out)
    00923         return (NULL);
    00924
    00925     /* Allocate memory */
    00926     if (NULL == (pcf = (float *) malloc(fft_size * sizeof(float))))
    00927         return (NULL);
    00928
    00929     /* Compute the Poisson complex filter */
    00930     if (NULL == poisson_complex_filter(pcf, nx, ny))
    00931         return (NULL);
    00932
    00933     /* Define sizes and indexes for for loops: depend of the parity */
    00934     inx = (int) nx;
    00935     iny = (int) ny;
    00936     xev = inx % 2;
    00937     yev = iny % 2;
    00938     hinx = inx / 2;
    00939     hiny = iny / 2;
    00940
    00941     /* Pointers initialization */
    00942     comp_ptr = data_out;
    00943     pcf_ptr = pcf;
    00944
    00945     for (y = 0; y < iny; y++) {
    00946         for (x = 0; x < sx; x++) {
    00947             /* Case 1: (x,y) is its own symmetric: (x,y)
    00948              * corresponds to a real coefficient of the DFT. */
    00949             /* A random sign is drawn (except for (0,0)
    00950              * where we impose sign=1. to conserve the
    00951              * original mean) */
    00952             if (((x == 0) || ((xev == 0) && (x == hinx)))
    00953                 && ((y == 0) || ((yev == 0) && (y == hiny)))) {
    00954                 sign = (((mt_drand53() < 0.5)
    00955                          || ((x == 0) && (y == 0))) ? 1. : -1.);
    00956                 (*comp_ptr)[0] = sign * (*pcf_ptr);
    00957                 (*comp_ptr)[1] = 0.;
    00958                 comp_ptr++;
    00959                 pcf_ptr++;
    00960             }
    00961             /* Case 2: Both (x,y) and its symmetric are in
    00962              * the same half-domain, and y > ny/2. */
    00963             /* Then the random phase of this symmetric
    00964              * point has already been drawn. */
    00965             else if (((x == 0) || ((xev == 0) && (x == hinx)))
    00966                      && (y > hiny)) {
    00967                 /* Copy the symmetric point */
    00968                 ad = (iny - y) * sx + x;        /* adress of the
    00969                                                  * symmetric point */
    00970                 (*comp_ptr)[0] = data_out[ad][0];
    00971                 (*comp_ptr)[1] = data_out[ad][1];
    00972                 comp_ptr++;
    00973                 pcf_ptr++;
    00974             }
    00975             else {
    00976                 /* Draw a random phase */
    00977                 theta = (2. * mt_drand53() - 1.) * M_PI;
    00978                 costheta = (float) cos(theta);
    00979                 sintheta = (float) sin(theta);
    00980                 (*comp_ptr)[0] = costheta * (*pcf_ptr);
    00981                 (*comp_ptr)[1] = sintheta * (*pcf_ptr);
    00982                 comp_ptr++;
    00983                 pcf_ptr++;
    00984             }
    00985         }
    00986     }
    00987     /* Free memory */
    00988     free(pcf);
    00989
    00990     return (data_out);
    00991 }
    00992
    01013 float **rpn_color_same_size(float **data_in_rgb, size_t nx, size_t ny)
    01014 {
    01015     float *laplacian;           /* 2D discrete Laplacian */
    01016     fftwf_complex *rp_pcf;      /* random phase with poisson
    01017                                  * complex filter */
    01018     fftwf_plan plan_r2c, plan_c2r;
    01019     fftwf_complex *fft;
    01020     float m;                    /* mean of the coefficient of
    01021                                  * data_in */
    01022     int channel;
    01023     size_t N = nx * ny;
    01024     size_t fft_size = (nx / 2 + 1) * ny;        /* physical size of the r2c FFT */
    01025
    01026     /* Check allocation */
    01027     if (NULL == data_in_rgb)
    01028         return (NULL);
    01029     /* Start threaded fftw if FFTW_NTHREADS is defined */
    01030 #ifdef FFTW_NTHREADS
    01031     if (0 == fftwf_init_threads())
    01032         return (NULL);
    01033     fftwf_plan_with_nthreads(FFTW_NTHREADS);
    01034 #endif
    01035     /* Allocate memory */
    01036     if (NULL == (laplacian = (float *) malloc(N * sizeof(float)))
    01037         || NULL == (rp_pcf = (fftwf_complex *)
    01038                     fftwf_malloc(fft_size * sizeof(fftwf_complex)))
    01039         || NULL == (fft = (fftwf_complex *)
    01040                     fftwf_malloc(fft_size * sizeof(fftwf_complex))))
    01041         return (NULL);
    01042
    01043     /* Compute the random phase with poisson complex filter */
    01044     if (NULL == random_phase_with_poisson_complex_filter(rp_pcf, nx, ny))
    01045         return (NULL);
    01046
    01047     /* For each channel */
    01048     for (channel = 0; channel < 3; channel++) {
    01049         /* Compute m = the mean value of data_in[channel] */
    01050         m = mean(data_in_rgb[channel], N);
    01051
    01052         /* Compute the discrete Laplacian of data_in[channel] */
    01053         if (NULL == discrete_laplacian(laplacian,
    01054                                        data_in_rgb[channel], nx, ny))
    01055             return (NULL);
    01056
    01057         /* Forward Fourier transform of the Laplacian: create
    01058          * the FFT forward plan and run the FFT */
    01059         plan_r2c = fftwf_plan_dft_r2c_2d((int) (ny), (int) (nx),
    01060                                          laplacian, fft,
    01061                                          FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
    01062         fftwf_execute(plan_r2c);
    01063
    01064         /* Inverse the discrete periodic Laplacian in the
    01065          * Fourier domain using the Poisson complex filter, and
    01066          * randomize the phase */
    01067         if (NULL == pointwise_complex_multiplication(fft, rp_pcf, fft_size))
    01068             return (NULL);
    01069
    01070         /* (0,0) frequency : we impose the same mean as
    01071          * data_in[channel] */
    01072         (*fft)[0] = m;
    01073         (*fft)[1] = 0.;
    01074
    01075         /* Backward Fourier transform: the output is stored in
    01076          * data_in_rgb[channel] */
    01077         plan_c2r = fftwf_plan_dft_c2r_2d((int) (ny), (int) (nx),
    01078                                          fft, data_in_rgb[channel],
    01079                                          FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
    01080         fftwf_execute(plan_c2r);
    01081     }
    01082
    01083     /* Cleanup */
    01084     free(laplacian);
    01085
    01086     /* Destroy the FFT plans and data */
    01087     fftwf_destroy_plan(plan_r2c);
    01088     fftwf_destroy_plan(plan_c2r);
    01089     fftwf_free(fft);
    01090     fftwf_free(rp_pcf);
    01091     fftwf_cleanup();
    01092 #ifdef FFTW_NTHREADS
    01093     fftwf_cleanup_threads();
    01094 #endif
    01095
    01096     return (data_in_rgb);
    01097 }