1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
8 //
9 //
10 // Intel License Agreement
11 // For Open Source Computer Vision Library
12 //
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
15 //
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
18 //
19 // * Redistribution's of source code must retain the above copyright notice,
20 // this list of conditions and the following disclaimer.
21 //
22 // * Redistribution's in binary form must reproduce the above copyright notice,
23 // this list of conditions and the following disclaimer in the documentation
24 // and/or other materials provided with the distribution.
25 //
26 // * The name of Intel Corporation may not be used to endorse or promote products
27 // derived from this software without specific prior written permission.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the Intel Corporation or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
39 //
40 //M*/
41
42 /* ////////////////////////////////////////////////////////////////////
43 //
44 // Filling CvMat/IplImage instances with random numbers
45 //
46 // */
47
48 #include "_cxcore.h"
49
50
51 ///////////////////////////// Functions Declaration //////////////////////////////////////
52
53 /*
54 Multiply-with-carry generator is used here:
55 temp = ( A*X(n) + carry )
56 X(n+1) = temp mod (2^32)
57 carry = temp / (2^32)
58 */
59 #define ICV_RNG_NEXT(x) ((uint64)(unsigned)(x)*1554115554 + ((x) >> 32))
60 #define ICV_CVT_FLT(x) (((unsigned)(x) >> 9)|CV_1F)
61 #define ICV_1D CV_BIG_INT(0x3FF0000000000000)
62 #define ICV_CVT_DBL(x) (((uint64)(unsigned)(x) << 20)|((x) >> 44)|ICV_1D)
63
64 /***************************************************************************************\
65 * Pseudo-Random Number Generators (PRNGs) *
66 \***************************************************************************************/
67
68 #define ICV_IMPL_RAND_BITS( flavor, arrtype, cast_macro ) \
69 static CvStatus CV_STDCALL \
70 icvRandBits_##flavor##_C1R( arrtype* arr, int step, CvSize size, \
71 uint64* state, const int* param ) \
72 { \
73 uint64 temp = *state; \
74 int small_flag = (param[12]|param[13]|param[14]|param[15]) <= 255; \
75 step /= sizeof(arr[0]); \
76 \
77 for( ; size.height--; arr += step ) \
78 { \
79 int i, k = 3; \
80 const int* p = param; \
81 \
82 if( !small_flag ) \
83 { \
84 for( i = 0; i <= size.width - 4; i += 4 ) \
85 { \
86 unsigned t0, t1; \
87 \
88 temp = ICV_RNG_NEXT(temp); \
89 t0 = ((unsigned)temp & p[i + 12]) + p[i]; \
90 temp = ICV_RNG_NEXT(temp); \
91 t1 = ((unsigned)temp & p[i + 13]) + p[i+1]; \
92 arr[i] = cast_macro((int)t0); \
93 arr[i+1] = cast_macro((int)t1); \
94 \
95 temp = ICV_RNG_NEXT(temp); \
96 t0 = ((unsigned)temp & p[i + 14]) + p[i+2]; \
97 temp = ICV_RNG_NEXT(temp); \
98 t1 = ((unsigned)temp & p[i + 15]) + p[i+3]; \
99 arr[i+2] = cast_macro((int)t0); \
100 arr[i+3] = cast_macro((int)t1); \
101 \
102 if( !--k ) \
103 { \
104 k = 3; \
105 p -= 12; \
106 } \
107 } \
108 } \
109 else \
110 { \
111 for( i = 0; i <= size.width - 4; i += 4 ) \
112 { \
113 unsigned t0, t1, t; \
114 \
115 temp = ICV_RNG_NEXT(temp); \
116 t = (unsigned)temp; \
117 t0 = (t & p[i + 12]) + p[i]; \
118 t1 = ((t >> 8) & p[i + 13]) + p[i+1]; \
119 arr[i] = cast_macro((int)t0); \
120 arr[i+1] = cast_macro((int)t1); \
121 \
122 t0 = ((t >> 16) & p[i + 14]) + p[i + 2]; \
123 t1 = ((t >> 24) & p[i + 15]) + p[i + 3]; \
124 arr[i+2] = cast_macro((int)t0); \
125 arr[i+3] = cast_macro((int)t1); \
126 \
127 if( !--k ) \
128 { \
129 k = 3; \
130 p -= 12; \
131 } \
132 } \
133 } \
134 \
135 for( ; i < size.width; i++ ) \
136 { \
137 unsigned t0; \
138 temp = ICV_RNG_NEXT(temp); \
139 \
140 t0 = ((unsigned)temp & p[i + 12]) + p[i]; \
141 arr[i] = cast_macro((int)t0); \
142 } \
143 } \
144 \
145 *state = temp; \
146 return CV_OK; \
147 }
148
149
150 #define ICV_IMPL_RAND( flavor, arrtype, worktype, cast_macro1, cast_macro2 )\
151 static CvStatus CV_STDCALL \
152 icvRand_##flavor##_C1R( arrtype* arr, int step, CvSize size, \
153 uint64* state, const double* param ) \
154 { \
155 uint64 temp = *state; \
156 step /= sizeof(arr[0]); \
157 \
158 for( ; size.height--; arr += step ) \
159 { \
160 int i, k = 3; \
161 const double* p = param; \
162 \
163 for( i = 0; i <= size.width - 4; i += 4 ) \
164 { \
165 worktype f0, f1; \
166 Cv32suf t0, t1; \
167 \
168 temp = ICV_RNG_NEXT(temp); \
169 t0.u = ICV_CVT_FLT(temp); \
170 temp = ICV_RNG_NEXT(temp); \
171 t1.u = ICV_CVT_FLT(temp); \
172 f0 = cast_macro1( t0.f * p[i + 12] + p[i] ); \
173 f1 = cast_macro1( t1.f * p[i + 13] + p[i + 1] ); \
174 arr[i] = cast_macro2(f0); \
175 arr[i+1] = cast_macro2(f1); \
176 \
177 temp = ICV_RNG_NEXT(temp); \
178 t0.u = ICV_CVT_FLT(temp); \
179 temp = ICV_RNG_NEXT(temp); \
180 t1.u = ICV_CVT_FLT(temp); \
181 f0 = cast_macro1( t0.f * p[i + 14] + p[i + 2] ); \
182 f1 = cast_macro1( t1.f * p[i + 15] + p[i + 3] ); \
183 arr[i+2] = cast_macro2(f0); \
184 arr[i+3] = cast_macro2(f1); \
185 \
186 if( !--k ) \
187 { \
188 k = 3; \
189 p -= 12; \
190 } \
191 } \
192 \
193 for( ; i < size.width; i++ ) \
194 { \
195 worktype f0; \
196 Cv32suf t0; \
197 \
198 temp = ICV_RNG_NEXT(temp); \
199 t0.u = ICV_CVT_FLT(temp); \
200 f0 = cast_macro1( t0.f * p[i + 12] + p[i] ); \
201 arr[i] = cast_macro2(f0); \
202 } \
203 } \
204 \
205 *state = temp; \
206 return CV_OK; \
207 }
208
209
210 static CvStatus CV_STDCALL
icvRand_64f_C1R(double * arr,int step,CvSize size,uint64 * state,const double * param)211 icvRand_64f_C1R( double* arr, int step, CvSize size,
212 uint64* state, const double* param )
213 {
214 uint64 temp = *state;
215 step /= sizeof(arr[0]);
216
217 for( ; size.height--; arr += step )
218 {
219 int i, k = 3;
220 const double* p = param;
221
222 for( i = 0; i <= size.width - 4; i += 4 )
223 {
224 double f0, f1;
225 Cv64suf t0, t1;
226
227 temp = ICV_RNG_NEXT(temp);
228 t0.u = ICV_CVT_DBL(temp);
229 temp = ICV_RNG_NEXT(temp);
230 t1.u = ICV_CVT_DBL(temp);
231 f0 = t0.f * p[i + 12] + p[i];
232 f1 = t1.f * p[i + 13] + p[i + 1];
233 arr[i] = f0;
234 arr[i+1] = f1;
235
236 temp = ICV_RNG_NEXT(temp);
237 t0.u = ICV_CVT_DBL(temp);
238 temp = ICV_RNG_NEXT(temp);
239 t1.u = ICV_CVT_DBL(temp);
240 f0 = t0.f * p[i + 14] + p[i + 2];
241 f1 = t1.f * p[i + 15] + p[i + 3];
242 arr[i+2] = f0;
243 arr[i+3] = f1;
244
245 if( !--k )
246 {
247 k = 3;
248 p -= 12;
249 }
250 }
251
252 for( ; i < size.width; i++ )
253 {
254 double f0;
255 Cv64suf t0;
256
257 temp = ICV_RNG_NEXT(temp);
258 t0.u = ICV_CVT_DBL(temp);
259 f0 = t0.f * p[i + 12] + p[i];
260 arr[i] = f0;
261 }
262 }
263
264 *state = temp;
265 return CV_OK;
266 }
267
268
269 /***************************************************************************************\
270 The code below implements algorithm from the paper:
271
272 G. Marsaglia and W.W. Tsang,
273 The Monty Python method for generating random variables,
274 ACM Transactions on Mathematical Software, Vol. 24, No. 3,
275 Pages 341-350, September, 1998.
276 \***************************************************************************************/
277
278 static CvStatus CV_STDCALL
icvRandn_0_1_32f_C1R(float * arr,int len,uint64 * state)279 icvRandn_0_1_32f_C1R( float* arr, int len, uint64* state )
280 {
281 uint64 temp = *state;
282 int i;
283 temp = ICV_RNG_NEXT(temp);
284
285 for( i = 0; i < len; i++ )
286 {
287 double x, y, v, ax, bx;
288
289 for(;;)
290 {
291 x = ((int)temp)*1.167239e-9;
292 temp = ICV_RNG_NEXT(temp);
293 ax = fabs(x);
294 v = 2.8658 - ax*(2.0213 - 0.3605*ax);
295 y = ((unsigned)temp)*2.328306e-10;
296 temp = ICV_RNG_NEXT(temp);
297
298 if( y < v || ax < 1.17741 )
299 break;
300
301 bx = x;
302 x = bx > 0 ? 0.8857913*(2.506628 - ax) : -0.8857913*(2.506628 - ax);
303
304 if( y > v + 0.0506 )
305 break;
306
307 if( log(y) < .6931472 - .5*bx*bx )
308 {
309 x = bx;
310 break;
311 }
312
313 if( log(1.8857913 - y) < .5718733-.5*x*x )
314 break;
315
316 do
317 {
318 v = ((int)temp)*4.656613e-10;
319 x = -log(fabs(v))*.3989423;
320 temp = ICV_RNG_NEXT(temp);
321 y = -log(((unsigned)temp)*2.328306e-10);
322 temp = ICV_RNG_NEXT(temp);
323 }
324 while( y+y < x*x );
325
326 x = v > 0 ? 2.506628 + x : -2.506628 - x;
327 break;
328 }
329
330 arr[i] = (float)x;
331 }
332 *state = temp;
333 return CV_OK;
334 }
335
336
337 #define RAND_BUF_SIZE 96
338
339
340 #define ICV_IMPL_RANDN( flavor, arrtype, worktype, cast_macro1, cast_macro2 ) \
341 static CvStatus CV_STDCALL \
342 icvRandn_##flavor##_C1R( arrtype* arr, int step, CvSize size, \
343 uint64* state, const double* param ) \
344 { \
345 float buffer[RAND_BUF_SIZE]; \
346 step /= sizeof(arr[0]); \
347 \
348 for( ; size.height--; arr += step ) \
349 { \
350 int i, j, len = RAND_BUF_SIZE; \
351 \
352 for( i = 0; i < size.width; i += RAND_BUF_SIZE ) \
353 { \
354 int k = 3; \
355 const double* p = param; \
356 \
357 if( i + len > size.width ) \
358 len = size.width - i; \
359 \
360 icvRandn_0_1_32f_C1R( buffer, len, state ); \
361 \
362 for( j = 0; j <= len - 4; j += 4 ) \
363 { \
364 worktype f0, f1; \
365 \
366 f0 = cast_macro1( buffer[j]*p[j+12] + p[j] ); \
367 f1 = cast_macro1( buffer[j+1]*p[j+13] + p[j+1] ); \
368 arr[i+j] = cast_macro2(f0); \
369 arr[i+j+1] = cast_macro2(f1); \
370 \
371 f0 = cast_macro1( buffer[j+2]*p[j+14] + p[j+2] ); \
372 f1 = cast_macro1( buffer[j+3]*p[j+15] + p[j+3] ); \
373 arr[i+j+2] = cast_macro2(f0); \
374 arr[i+j+3] = cast_macro2(f1); \
375 \
376 if( --k == 0 ) \
377 { \
378 k = 3; \
379 p -= 12; \
380 } \
381 } \
382 \
383 for( ; j < len; j++ ) \
384 { \
385 worktype f0 = cast_macro1( buffer[j]*p[j+12] + p[j] ); \
386 arr[i+j] = cast_macro2(f0); \
387 } \
388 } \
389 } \
390 \
391 return CV_OK; \
392 }
393
394
395 ICV_IMPL_RAND_BITS( 8u, uchar, CV_CAST_8U )
396 ICV_IMPL_RAND_BITS( 16u, ushort, CV_CAST_16U )
397 ICV_IMPL_RAND_BITS( 16s, short, CV_CAST_16S )
398 ICV_IMPL_RAND_BITS( 32s, int, CV_CAST_32S )
399
400 ICV_IMPL_RAND( 8u, uchar, int, cvFloor, CV_CAST_8U )
401 ICV_IMPL_RAND( 16u, ushort, int, cvFloor, CV_CAST_16U )
402 ICV_IMPL_RAND( 16s, short, int, cvFloor, CV_CAST_16S )
403 ICV_IMPL_RAND( 32s, int, int, cvFloor, CV_CAST_32S )
404 ICV_IMPL_RAND( 32f, float, float, CV_CAST_32F, CV_NOP )
405
406 ICV_IMPL_RANDN( 8u, uchar, int, cvRound, CV_CAST_8U )
407 ICV_IMPL_RANDN( 16u, ushort, int, cvRound, CV_CAST_16U )
408 ICV_IMPL_RANDN( 16s, short, int, cvRound, CV_CAST_16S )
409 ICV_IMPL_RANDN( 32s, int, int, cvRound, CV_CAST_32S )
410 ICV_IMPL_RANDN( 32f, float, float, CV_CAST_32F, CV_NOP )
411 ICV_IMPL_RANDN( 64f, double, double, CV_CAST_64F, CV_NOP )
412
icvInitRandTable(CvFuncTable * fastrng_tab,CvFuncTable * rng_tab,CvFuncTable * normal_tab)413 static void icvInitRandTable( CvFuncTable* fastrng_tab,
414 CvFuncTable* rng_tab,
415 CvFuncTable* normal_tab )
416 {
417 fastrng_tab->fn_2d[CV_8U] = (void*)icvRandBits_8u_C1R;
418 fastrng_tab->fn_2d[CV_8S] = 0;
419 fastrng_tab->fn_2d[CV_16U] = (void*)icvRandBits_16u_C1R;
420 fastrng_tab->fn_2d[CV_16S] = (void*)icvRandBits_16s_C1R;
421 fastrng_tab->fn_2d[CV_32S] = (void*)icvRandBits_32s_C1R;
422
423 rng_tab->fn_2d[CV_8U] = (void*)icvRand_8u_C1R;
424 rng_tab->fn_2d[CV_8S] = 0;
425 rng_tab->fn_2d[CV_16U] = (void*)icvRand_16u_C1R;
426 rng_tab->fn_2d[CV_16S] = (void*)icvRand_16s_C1R;
427 rng_tab->fn_2d[CV_32S] = (void*)icvRand_32s_C1R;
428 rng_tab->fn_2d[CV_32F] = (void*)icvRand_32f_C1R;
429 rng_tab->fn_2d[CV_64F] = (void*)icvRand_64f_C1R;
430
431 normal_tab->fn_2d[CV_8U] = (void*)icvRandn_8u_C1R;
432 normal_tab->fn_2d[CV_8S] = 0;
433 normal_tab->fn_2d[CV_16U] = (void*)icvRandn_16u_C1R;
434 normal_tab->fn_2d[CV_16S] = (void*)icvRandn_16s_C1R;
435 normal_tab->fn_2d[CV_32S] = (void*)icvRandn_32s_C1R;
436 normal_tab->fn_2d[CV_32F] = (void*)icvRandn_32f_C1R;
437 normal_tab->fn_2d[CV_64F] = (void*)icvRandn_64f_C1R;
438 }
439
440
441 CV_IMPL void
cvRandArr(CvRNG * rng,CvArr * arr,int disttype,CvScalar param1,CvScalar param2)442 cvRandArr( CvRNG* rng, CvArr* arr, int disttype, CvScalar param1, CvScalar param2 )
443 {
444 static CvFuncTable rng_tab[2], fastrng_tab;
445 static int inittab = 0;
446
447 CV_FUNCNAME( "cvRandArr" );
448
449 __BEGIN__;
450
451 int is_nd = 0;
452 CvMat stub, *mat = (CvMat*)arr;
453 int type, depth, channels;
454 double dparam[2][12];
455 int iparam[2][12];
456 void* param = dparam;
457 int i, fast_int_mode = 0;
458 int mat_step = 0;
459 CvSize size;
460 CvFunc2D_1A2P func = 0;
461 CvMatND stub_nd;
462 CvNArrayIterator iterator_state, *iterator = 0;
463
464 if( !inittab )
465 {
466 icvInitRandTable( &fastrng_tab, &rng_tab[CV_RAND_UNI],
467 &rng_tab[CV_RAND_NORMAL] );
468 inittab = 1;
469 }
470
471 if( !rng )
472 CV_ERROR( CV_StsNullPtr, "Null pointer to RNG state" );
473
474 if( CV_IS_MATND(mat) )
475 {
476 iterator = &iterator_state;
477 CV_CALL( cvInitNArrayIterator( 1, &arr, 0, &stub_nd, iterator ));
478 type = CV_MAT_TYPE(iterator->hdr[0]->type);
479 size = iterator->size;
480 is_nd = 1;
481 }
482 else
483 {
484 if( !CV_IS_MAT(mat))
485 {
486 int coi = 0;
487 CV_CALL( mat = cvGetMat( mat, &stub, &coi ));
488
489 if( coi != 0 )
490 CV_ERROR( CV_BadCOI, "COI is not supported" );
491 }
492
493 type = CV_MAT_TYPE( mat->type );
494 size = cvGetMatSize( mat );
495 mat_step = mat->step;
496
497 if( mat->height > 1 && CV_IS_MAT_CONT( mat->type ))
498 {
499 size.width *= size.height;
500 mat_step = CV_STUB_STEP;
501 size.height = 1;
502 }
503 }
504
505 depth = CV_MAT_DEPTH( type );
506 channels = CV_MAT_CN( type );
507 size.width *= channels;
508
509 if( disttype == CV_RAND_UNI )
510 {
511 if( depth <= CV_32S )
512 {
513 for( i = 0, fast_int_mode = 1; i < channels; i++ )
514 {
515 int t0 = iparam[0][i] = cvCeil( param1.val[i] );
516 int t1 = iparam[1][i] = cvFloor( param2.val[i] ) - t0;
517 double diff = param1.val[i] - param2.val[i];
518
519 fast_int_mode &= INT_MIN <= diff && diff <= INT_MAX && (t1 & (t1 - 1)) == 0;
520 }
521 }
522
523 if( fast_int_mode )
524 {
525 for( i = 0; i < channels; i++ )
526 iparam[1][i]--;
527
528 for( ; i < 12; i++ )
529 {
530 int t0 = iparam[0][i - channels];
531 int t1 = iparam[1][i - channels];
532
533 iparam[0][i] = t0;
534 iparam[1][i] = t1;
535 }
536
537 CV_GET_FUNC_PTR( func, (CvFunc2D_1A2P)(fastrng_tab.fn_2d[depth]));
538 param = iparam;
539 }
540 else
541 {
542 for( i = 0; i < channels; i++ )
543 {
544 double t0 = param1.val[i];
545 double t1 = param2.val[i];
546
547 dparam[0][i] = t0 - (t1 - t0);
548 dparam[1][i] = t1 - t0;
549 }
550
551 CV_GET_FUNC_PTR( func, (CvFunc2D_1A2P)(rng_tab[0].fn_2d[depth]));
552 }
553 }
554 else if( disttype == CV_RAND_NORMAL )
555 {
556 for( i = 0; i < channels; i++ )
557 {
558 double t0 = param1.val[i];
559 double t1 = param2.val[i];
560
561 dparam[0][i] = t0;
562 dparam[1][i] = t1;
563 }
564
565 CV_GET_FUNC_PTR( func, (CvFunc2D_1A2P)(rng_tab[1].fn_2d[depth]));
566 }
567 else
568 {
569 CV_ERROR( CV_StsBadArg, "Unknown distribution type" );
570 }
571
572 if( !fast_int_mode )
573 {
574 for( i = channels; i < 12; i++ )
575 {
576 double t0 = dparam[0][i - channels];
577 double t1 = dparam[1][i - channels];
578
579 dparam[0][i] = t0;
580 dparam[1][i] = t1;
581 }
582 }
583
584 if( !is_nd )
585 {
586 IPPI_CALL( func( mat->data.ptr, mat_step, size, rng, param ));
587 }
588 else
589 {
590 do
591 {
592 IPPI_CALL( func( iterator->ptr[0], CV_STUB_STEP, size, rng, param ));
593 }
594 while( cvNextNArraySlice( iterator ));
595 }
596
597 __END__;
598 }
599
600 /* End of file. */
601