libs/filters/GP_Rand.c

   1 /*****************************************************************************
   2  * This file is part of gfxprim library.                                     *
   3  *                                                                           *
   4  * Gfxprim is free software; you can redistribute it and/or                  *
   5  * modify it under the terms of the GNU Lesser General Public                *
   6  * License as published by the Free Software Foundation; either              *
   7  * version 2.1 of the License, or (at your option) any later version.        *
   8  *                                                                           *
   9  * Gfxprim is distributed in the hope that it will be useful,                *
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU         *
  12  * Lesser General Public License for more details.                           *
  13  *                                                                           *
  14  * You should have received a copy of the GNU Lesser General Public          *
  15  * License along with gfxprim; if not, write to the Free Software            *
  16  * Foundation, Inc., 51 Franklin Street, Fifth Floor,                        *
  17  * Boston, MA  02110-1301  USA                                               *
  18  *                                                                           *
  19  * Copyright (C) 2009-2012 Cyril Hrubis <metan@ucw.cz>                       *
  20  *                                                                           *
  21  *****************************************************************************/
  22
  23 #include <math.h>
  24
  25 #include "core/GP_Common.h"
  26
  27 #include "GP_Rand.h"
  28
  29 void GP_NormInt(int *arr, unsigned int size, int sigma, int mu)
  30 {
  31         unsigned int i = 0;
  32         float a, b, rsq;
  33
  34         GP_ASSERT(size%2 == 0);
  35
  36         while (i < size) {
  37                 /* Sample two point inside of the unit circle */
  38                 do {
  39                         a = (float)random() / (RAND_MAX/2) - 1;
  40                         b = (float)random() / (RAND_MAX/2) - 1;
  41
  42                         rsq = a * a + b * b;
  43
  44                 } while (rsq >= 1 || rsq == 0);
  45
  46                 /* Polar form of the Box-Muller transformation */
  47                 float mul = sqrtf(-2.0 * logf(rsq)/rsq);
  48
  49                 arr[i++] = mu + sigma * a * mul;
  50                 arr[i++] = mu + sigma * b * mul;
  51         }
  52 }