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