Rudin-Osher-Fatemi Total Variation Denoising using Split Bregman
randmt.c
Go to the documentation of this file.
1 
65 #include <stdio.h>
66 #include <stdlib.h>
67 #include <limits.h>
68 #include <math.h>
69 #include <time.h>
70 #include "randmt.h"
71 
72 /* Get the high-precision timer specified by POSIX.1-2001 if available */
73 #if defined(USE_GETTIMEOFDAY) || defined(unix) || defined(__unix__) || defined(__unix)
74 #include <unistd.h>
75 #if defined(USE_GETTIMEOFDAY) || (defined(_POSIX_VERSION) && (_POSIX_VERSION >= 200112L))
76 #include <sys/time.h>
77 #define USE_GETTIMEOFDAY
78 #endif
79 #endif
80 
81 #ifndef M_LOGSQRT2PI
82 
83 #define M_LOGSQRT2PI 0.9189385332046727417803297
84 #endif
85 
86 /* Period parameters */
87 #define MT_N 624
88 #define MT_M 397
89 #define MT_MATRIX_A 0x9908B0DFUL
90 #define MT_UPPER_MASK 0x80000000UL
91 #define MT_LOWER_MASK 0x7FFFFFFFUL
95 typedef struct randmtstruct_t
96 {
97  unsigned long mt[MT_N];
98  int mti;
99 } randmttype_t;
100 
103 
104 
105 /* Create a new randmt_t */
107 {
108  randmt_t *generator;
109 
110  if((generator = (randmt_t *)calloc(1, sizeof(randmt_t))))
111  {
112  /* mti == MT_N+1 means mt[MT_N] is not initialized */
113  generator->mti = MT_N + 1;
114  }
115 
116  return generator;
117 }
118 
119 
120 /* Free a randmt_t */
121 void free_randmt(randmt_t *generator)
122 {
123  if(generator)
124  free(generator);
125  return;
126 }
127 
128 
129 /* Initialize randmt_t with a seed */
130 void init_randmt_r(randmt_t *generator, unsigned long seed)
131 {
132  int mti;
133  unsigned long *mt;
134 
135  if(!generator)
136  return;
137 
138 #if ULONG_MAX > 0xFFFFFFFFUL
139  seed &= 0xFFFFFFFFUL; /* for WORDSIZE > 32bit */
140 #endif
141  (mt = generator->mt)[0] = seed;
142 
143  for(mti = 1; mti < MT_N; mti++)
144  {
145  mt[mti] = (1812433253UL * (mt[mti - 1] ^ (mt[mti - 1] >> 30)) + mti);
146  /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
147  * In the previous versions, MSBs of the seed affect
148  * only MSBs of the array mt[].
149  * 2002/01/09 modified by Makoto Matsumoto
150  */
151 #if ULONG_MAX > 0xFFFFFFFFUL
152  mt[mti] &= 0xFFFFFFFFUL; /* for WORDSIZE > 32bit machines */
153 #endif
154  }
155 
156  generator->mti = MT_N;
157  return;
158 }
159 
160 
161 /* Initialize generator with the current time and memory address */
162 void init_randmt_auto_r(randmt_t *generator)
163 {
164  unsigned long int seed;
165 
166  /* The basic (weak) initialization uses the current time, in seconds */
167  seed = (unsigned long int) time(NULL);
168  /* Add to the generator's memory address so that different generators
169  use different seeds. */
170  seed += ((unsigned long int) generator);
171 
172 #if defined(USE_GETTIMEOFDAY)
173  /* gettimeofday() provides microsecond resolution */
174  {
175  struct timeval tp;
176  (void) gettimeofday(&tp, NULL);
177 
178  seed *= 1000000;
179  seed += tp.tv_usec;
180  }
181 #endif
182 
183  init_randmt(seed);
184  return;
185 }
186 
187 
188 /* Generate a random 32-bit unsigned integer (reentrant) */
189 unsigned long rand_uint32_r(randmt_t *generator)
190 {
191  /* mag01[x] = x * MT_MATRIX_A for x=0,1 */
192  const unsigned long mag01[2] = { 0x0UL, MT_MATRIX_A };
193  unsigned long y, *mt;
194  int mti;
195 
196  if(!generator)
197  return 0;
198 
199  mt = generator->mt;
200  mti = generator->mti;
201 
202  if (mti >= MT_N) /* generate MT_N words at one time */
203  {
204  int kk;
205 
206  /* If init_randmt_r() has not been called, use a default seed. */
207  if (mti == MT_N + 1)
208  init_randmt_r(generator, 5489UL);
209 
210  for(kk = 0; kk < MT_N - MT_M; kk++)
211  {
212  y = (mt[kk] & MT_UPPER_MASK) | (mt[kk + 1] & MT_LOWER_MASK);
213  mt[kk] = mt[kk + MT_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
214  }
215  for(; kk < MT_N - 1; kk++)
216  {
217  y = (mt[kk] & MT_UPPER_MASK) | (mt[kk + 1] & MT_LOWER_MASK);
218  mt[kk] = mt[kk + (MT_M - MT_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
219  }
220 
221  y = (mt[MT_N - 1] & MT_UPPER_MASK) | (mt[0] & MT_LOWER_MASK);
222  mt[MT_N - 1] = mt[MT_M - 1] ^ (y >> 1) ^ mag01[y & 0x1UL];
223  mti = 0;
224  }
225 
226  y = mt[mti++];
227 
228  /* Tempering */
229  y ^= (y >> 11);
230  y ^= (y << 7) & 0x9D2C5680UL;
231  y ^= (y << 15) & 0xEFC60000UL;
232  y ^= (y >> 18);
233 
234  generator->mti = mti;
235 
236  return y;
237 }
238 
239 
240 /* Generate a Gamma-distributed number (reentrant) */
241 double rand_gamma_r(randmt_t *generator, double a, double b)
242 {
243  const double d = a - 1.0/3.0;
244  const double c = 1.0/sqrt(9*d);
245  double x, v, u;
246 
247  if(a >= 1)
248  {
249  do
250  {
251  do
252  {
253  x = rand_normal_r(generator);
254  v = 1 + c*x;
255  }while(v <= 0);
256 
257  v = v*v*v;
258  u = rand_unif_r(generator);
259  x = x*x;
260  }while(u >= 1 - 0.0331*x*x
261  && log(u) >= x/2 + d*(1 - v + log(v)));
262 
263  return b*d*v;
264  }
265  else if(a > 0)
266  return rand_gamma_r(generator, 1 + a, b)
267  * pow(rand_unif_r(generator), 1/a);
268  else
269  return 0;
270 }
271 
272 
273 /* Compute the natural logarithm of the factorial (used by rand_poisson_r) */
274 static double logfactorial(double n)
275 {
276  /* Look-up table for log(n!) for n = 2, ..., 11 */
277  static const double Table[10] = {
278  0.69314718055994530942, 1.7917594692280550008,
279  3.1780538303479456197, 4.7874917427820459943,
280  6.5792512120101009951, 8.5251613610654143002,
281  10.604602902745250228, 12.801827480081469611,
282  15.104412573075515295, 17.502307845873885839};
283  /* Asymptotic approximation coefficients */
284  static const double C[7] = { -1.910444077728e-03,
285  8.4171387781295e-04, -5.952379913043012e-04,
286  7.93650793500350248e-04, -2.777777777777681622553e-03,
287  8.333333333333333331554247e-02, 5.7083835261e-03};
288  double x, xsq, res, corr;
289  int i;
290 
291  if(n > 11) /* Asymptotic approximation, based on the public domain */
292  { /* NETLIB (Fortran) code by W. J. Cody and L. Stoltz */
293  x = 1 + floor(n + 0.5);
294 
295  if(x <= 2.25e76)
296  {
297  res = C[6];
298  xsq = x*x;
299 
300  for(i = 0; i < 6; i++)
301  res = res/xsq + C[i];
302  }
303  else
304  res = 0;
305 
306  corr = log(x);
307  return res/x + M_LOGSQRT2PI - corr/2 + x*(corr - 1);
308  }
309  else /* Use look-up table */
310  {
311  i = (int)floor(n + 0.5);
312  return (i >= 2) ? Table[i - 2] : 0;
313  }
314 }
315 
316 
317 /* Generate a Poisson-distributed number (reentrant) */
318 double rand_poisson_r(randmt_t *generator, double mu)
319 {
320  double k, b, a, U, V, us, vr, invalpha, logmu;
321 
322  if(mu < 10)
323  { /* Use simple direct algorthm for small mu */
324  vr = exp(-mu);
325  k = -1;
326  V = 1;
327 
328  do
329  {
330  k += 1;
331  V *= rand_unif_r(generator);
332  }while(V > vr);
333  }
334  else
335  { /* Use the "PTRS" algorithm of Hormann */
336  b = 0.931 + 2.53*sqrt(mu);
337  a = -0.059 + 0.02483*b;
338  vr = 0.9277 - 3.6224/(b - 2);
339  invalpha = 1.1239 + 1.1328/(b - 3.4);
340  logmu = log(mu);
341 
342  do
343  {
344  U = rand_unif_r(generator) - 0.5;
345  V = rand_unif_r(generator);
346  us = 0.5 - fabs(U);
347  k = floor((2*a/us + b)*U + mu + 0.43);
348  }while((us < 0.07 || V > vr)
349  && (k < 0 || (us < 0.013 && V > us)
350  || log(V*invalpha/(a/(us*us) + b))
351  > -mu + k*logmu - logfactorial(k)));
352  }
353 
354  return k;
355 }