git.maemo.org Git - opencv/blob - cxcore/src/cxmathfuncs.cpp

   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 #include "_cxcore.h"
  43
  44 #ifdef HAVE_CONFIG_H
  45 #include <cvconfig.h>
  46 #endif
  47
  48 #define ICV_MATH_BLOCK_SIZE  256
  49
  50 #define _CV_SQRT_MAGIC     0xbe6f0000
  51
  52 #define _CV_SQRT_MAGIC_DBL CV_BIG_UINT(0xbfcd460000000000)
  53
  54 #define _CV_ATAN_CF0  (-15.8131890796f)
  55 #define _CV_ATAN_CF1  (61.0941945596f)
  56 #define _CV_ATAN_CF2  0.f /*(-0.140500406322f)*/
  57
  58 static const float icvAtanTab[8] = { 0.f + _CV_ATAN_CF2, 90.f - _CV_ATAN_CF2,
  59     180.f - _CV_ATAN_CF2, 90.f + _CV_ATAN_CF2,
  60     360.f - _CV_ATAN_CF2, 270.f + _CV_ATAN_CF2,
  61     180.f + _CV_ATAN_CF2, 270.f - _CV_ATAN_CF2
  62 };
  63
  64 static const int icvAtanSign[8] =
  65     { 0, 0x80000000, 0x80000000, 0, 0x80000000, 0, 0, 0x80000000 };
  66
  67 CV_IMPL float
  68 cvFastArctan( float y, float x )
  69 {
  70     Cv32suf _x, _y;
  71     int ix, iy, ygx, idx;
  72     double z;
  73
  74     _x.f = x; _y.f = y;
  75     ix = _x.i; iy = _y.i;
  76     idx = (ix < 0) * 2 + (iy < 0) * 4;
  77
  78     ix &= 0x7fffffff;
  79     iy &= 0x7fffffff;
  80
  81     ygx = (iy <= ix) - 1;
  82     idx -= ygx;
  83
  84     idx &= ((ix == 0) - 1) | ((iy == 0) - 1);
  85
  86     /* swap ix and iy if ix < iy */
  87     ix ^= iy & ygx;
  88     iy ^= ix & ygx;
  89     ix ^= iy & ygx;
  90
  91     _y.i = iy ^ icvAtanSign[idx];
  92
  93     /* ix = ix != 0 ? ix : 1.f */
  94     _x.i = ((ix ^ CV_1F) & ((ix == 0) - 1)) ^ CV_1F;
  95
  96     z = _y.f / _x.f;
  97     return (float)((_CV_ATAN_CF0*fabs(z) + _CV_ATAN_CF1)*z + icvAtanTab[idx]);
  98 }
  99
 100
 101 IPCVAPI_IMPL( CvStatus, icvFastArctan_32f,
 102     (const float *__y, const float *__x, float *angle, int len ), (__y, __x, angle, len) )
 103 {
 104     int i = 0;
 105     const int *y = (const int*)__y, *x = (const int*)__x;
 106
 107     if( !(y && x && angle && len >= 0) )
 108         return CV_BADFACTOR_ERR;
 109
 110     /* unrolled by 4 loop */
 111     for( ; i <= len - 4; i += 4 )
 112     {
 113         int j, idx[4];
 114         float xf[4], yf[4];
 115         double d = 1.;
 116
 117         /* calc numerators and denominators */
 118         for( j = 0; j < 4; j++ )
 119         {
 120             int ix = x[i + j], iy = y[i + j];
 121             int ygx, k = (ix < 0) * 2 + (iy < 0) * 4;
 122             Cv32suf _x, _y;
 123
 124             ix &= 0x7fffffff;
 125             iy &= 0x7fffffff;
 126
 127             ygx = (iy <= ix) - 1;
 128             k -= ygx;
 129
 130             k &= ((ix == 0) - 1) | ((iy == 0) - 1);
 131
 132             /* swap ix and iy if ix < iy */
 133             ix ^= iy & ygx;
 134             iy ^= ix & ygx;
 135             ix ^= iy & ygx;
 136
 137             _y.i = iy ^ icvAtanSign[k];
 138
 139             /* ix = ix != 0 ? ix : 1.f */
 140             _x.i = ((ix ^ CV_1F) & ((ix == 0) - 1)) ^ CV_1F;
 141             idx[j] = k;
 142             yf[j] = _y.f;
 143             d *= (xf[j] = _x.f);
 144         }
 145
 146         d = 1. / d;
 147
 148         {
 149             double b = xf[2] * xf[3], a = xf[0] * xf[1];
 150
 151             float z0 = (float) (yf[0] * xf[1] * b * d);
 152             float z1 = (float) (yf[1] * xf[0] * b * d);
 153             float z2 = (float) (yf[2] * xf[3] * a * d);
 154             float z3 = (float) (yf[3] * xf[2] * a * d);
 155
 156             z0 = (float)((_CV_ATAN_CF0*fabs(z0) + _CV_ATAN_CF1)*z0 + icvAtanTab[idx[0]]);
 157             z1 = (float)((_CV_ATAN_CF0*fabs(z1) + _CV_ATAN_CF1)*z1 + icvAtanTab[idx[1]]);
 158             z2 = (float)((_CV_ATAN_CF0*fabs(z2) + _CV_ATAN_CF1)*z2 + icvAtanTab[idx[2]]);
 159             z3 = (float)((_CV_ATAN_CF0*fabs(z3) + _CV_ATAN_CF1)*z3 + icvAtanTab[idx[3]]);
 160
 161             angle[i] = z0;
 162             angle[i+1] = z1;
 163             angle[i+2] = z2;
 164             angle[i+3] = z3;
 165         }
 166     }
 167
 168     /* process the rest */
 169     for( ; i < len; i++ )
 170         angle[i] = cvFastArctan( __y[i], __x[i] );
 171
 172     return CV_OK;
 173 }
 174
 175
 176 /* ************************************************************************** *\
 177    Fast cube root by Ken Turkowski
 178    (http://www.worldserver.com/turk/computergraphics/papers.html)
 179 \* ************************************************************************** */
 180 CV_IMPL  float  cvCbrt( float value )
 181 {
 182     float fr;
 183     Cv32suf v, m;
 184     int ix, s;
 185     int ex, shx;
 186
 187     v.f = value;
 188     ix = v.i & 0x7fffffff;
 189     s = v.i & 0x80000000;
 190     ex = (ix >> 23) - 127;
 191     shx = ex % 3;
 192     shx -= shx >= 0 ? 3 : 0;
 193     ex = (ex - shx) / 3; /* exponent of cube root */
 194     v.i = (ix & ((1<<23)-1)) | ((shx + 127)<<23);
 195     fr = v.f;
 196
 197     /* 0.125 <= fr < 1.0 */
 198     /* Use quartic rational polynomial with error < 2^(-24) */
 199     fr = (float)(((((45.2548339756803022511987494 * fr +
 200     192.2798368355061050458134625) * fr +
 201     119.1654824285581628956914143) * fr +
 202     13.43250139086239872172837314) * fr +
 203     0.1636161226585754240958355063)/
 204     ((((14.80884093219134573786480845 * fr +
 205     151.9714051044435648658557668) * fr +
 206     168.5254414101568283957668343) * fr +
 207     33.9905941350215598754191872) * fr +
 208     1.0));
 209
 210     /* fr *= 2^ex * sign */
 211     m.f = value;
 212     v.f = fr;
 213     v.i = (v.i + (ex << 23) + s) & (m.i*2 != 0 ? -1 : 0);
 214     return v.f;
 215 }
 216
 217 //static const double _0_5 = 0.5, _1_5 = 1.5;
 218
 219 IPCVAPI_IMPL( CvStatus, icvInvSqrt_32f, (const float *src, float *dst, int len), (src, dst, len) )
 220 {
 221     int i = 0;
 222
 223     if( !(src && dst && len >= 0) )
 224         return CV_BADFACTOR_ERR;
 225
 226     for( ; i < len; i++ )
 227         dst[i] = (float)(1.f/sqrt(src[i]));
 228
 229     return CV_OK;
 230 }
 231
 232
 233 IPCVAPI_IMPL( CvStatus, icvSqrt_32f, (const float *src, float *dst, int len), (src, dst, len) )
 234 {
 235     int i = 0;
 236
 237     if( !(src && dst && len >= 0) )
 238         return CV_BADFACTOR_ERR;
 239
 240     for( ; i < len; i++ )
 241         dst[i] = (float)sqrt(src[i]);
 242
 243     return CV_OK;
 244 }
 245
 246
 247 IPCVAPI_IMPL( CvStatus, icvSqrt_64f, (const double *src, double *dst, int len), (src, dst, len) )
 248 {
 249     int i = 0;
 250
 251     if( !(src && dst && len >= 0) )
 252         return CV_BADFACTOR_ERR;
 253
 254     for( ; i < len; i++ )
 255         dst[i] = sqrt(src[i]);
 256
 257     return CV_OK;
 258 }
 259
 260
 261 IPCVAPI_IMPL( CvStatus, icvInvSqrt_64f, (const double *src, double *dst, int len), (src, dst, len) )
 262 {
 263     int i = 0;
 264
 265     if( !(src && dst && len >= 0) )
 266         return CV_BADFACTOR_ERR;
 267
 268     for( ; i < len; i++ )
 269         dst[i] = 1./sqrt(src[i]);
 270
 271     return CV_OK;
 272 }
 273
 274
 275 #define ICV_DEF_SQR_MAGNITUDE_FUNC(flavor, arrtype, magtype)\
 276 static CvStatus CV_STDCALL                                  \
 277 icvSqrMagnitude_##flavor(const arrtype* x, const arrtype* y,\
 278                          magtype* mag, int len)             \
 279 {                                                           \
 280     int i;                                                  \
 281                                                             \
 282     for( i = 0; i <= len - 4; i += 4 )                      \
 283     {                                                       \
 284         magtype x0 = (magtype)x[i], y0 = (magtype)y[i];     \
 285         magtype x1 = (magtype)x[i+1], y1 = (magtype)y[i+1]; \
 286                                                             \
 287         x0 = x0*x0 + y0*y0;                                 \
 288         x1 = x1*x1 + y1*y1;                                 \
 289         mag[i] = x0;                                        \
 290         mag[i+1] = x1;                                      \
 291         x0 = (magtype)x[i+2], y0 = (magtype)y[i+2];         \
 292         x1 = (magtype)x[i+3], y1 = (magtype)y[i+3];         \
 293         x0 = x0*x0 + y0*y0;                                 \
 294         x1 = x1*x1 + y1*y1;                                 \
 295         mag[i+2] = x0;                                      \
 296         mag[i+3] = x1;                                      \
 297     }                                                       \
 298                                                             \
 299     for( ; i < len; i++ )                                   \
 300     {                                                       \
 301         magtype x0 = (magtype)x[i], y0 = (magtype)y[i];     \
 302         mag[i] = x0*x0 + y0*y0;                             \
 303     }                                                       \
 304                                                             \
 305     return CV_OK;                                           \
 306 }
 307
 308
 309 ICV_DEF_SQR_MAGNITUDE_FUNC( 32f, float, float )
 310 ICV_DEF_SQR_MAGNITUDE_FUNC( 64f, double, double )
 311
 312 /****************************************************************************************\
 313 *                                  Cartezian -> Polar                                    *
 314 \****************************************************************************************/
 315
 316 CV_IMPL void
 317 cvCartToPolar( const CvArr* xarr, const CvArr* yarr,
 318                CvArr* magarr, CvArr* anglearr,
 319                int angle_in_degrees )
 320 {
 321     CV_FUNCNAME( "cvCartToPolar" );
 322
 323     __BEGIN__;
 324
 325     float* mag_buffer = 0;
 326     float* x_buffer = 0;
 327     float* y_buffer = 0;
 328     int block_size = 0;
 329     CvMat xstub, *xmat = (CvMat*)xarr;
 330     CvMat ystub, *ymat = (CvMat*)yarr;
 331     CvMat magstub, *mag = (CvMat*)magarr;
 332     CvMat anglestub, *angle = (CvMat*)anglearr;
 333     int coi1 = 0, coi2 = 0, coi3 = 0, coi4 = 0;
 334     int depth;
 335     CvSize size;
 336     int x, y;
 337     int cont_flag = CV_MAT_CONT_FLAG;
 338
 339     if( !CV_IS_MAT(xmat))
 340         CV_CALL( xmat = cvGetMat( xmat, &xstub, &coi1 ));
 341
 342     if( !CV_IS_MAT(ymat))
 343         CV_CALL( ymat = cvGetMat( ymat, &ystub, &coi2 ));
 344
 345     if( !CV_ARE_TYPES_EQ( xmat, ymat ) )
 346         CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
 347
 348     if( !CV_ARE_SIZES_EQ( xmat, ymat ) )
 349         CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
 350
 351     depth = CV_MAT_DEPTH( xmat->type );
 352     if( depth < CV_32F )
 353         CV_ERROR( CV_StsUnsupportedFormat, "" );
 354
 355     if( mag )
 356     {
 357         CV_CALL( mag = cvGetMat( mag, &magstub, &coi3 ));
 358
 359         if( !CV_ARE_TYPES_EQ( mag, xmat ) )
 360             CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
 361
 362         if( !CV_ARE_SIZES_EQ( mag, xmat ) )
 363             CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
 364         cont_flag = mag->type;
 365     }
 366
 367     if( angle )
 368     {
 369         CV_CALL( angle = cvGetMat( angle, &anglestub, &coi4 ));
 370
 371         if( !CV_ARE_TYPES_EQ( angle, xmat ) )
 372             CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
 373
 374         if( !CV_ARE_SIZES_EQ( angle, xmat ) )
 375             CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
 376         cont_flag &= angle->type;
 377     }
 378
 379     if( coi1 != 0 || coi2 != 0 || coi3 != 0 || coi4 != 0 )
 380         CV_ERROR( CV_BadCOI, "" );
 381
 382     size = cvGetMatSize(xmat);
 383     size.width *= CV_MAT_CN(xmat->type);
 384
 385     if( CV_IS_MAT_CONT( xmat->type & ymat->type & cont_flag ))
 386     {
 387         size.width *= size.height;
 388         size.height = 1;
 389     }
 390
 391     block_size = MIN( size.width, ICV_MATH_BLOCK_SIZE );
 392     if( depth == CV_64F && angle )
 393     {
 394         x_buffer = (float*)cvStackAlloc( block_size*sizeof(float));
 395         y_buffer = (float*)cvStackAlloc( block_size*sizeof(float));
 396     }
 397     else if( depth == CV_32F && mag )
 398     {
 399         mag_buffer = (float*)cvStackAlloc( block_size*sizeof(float));
 400     }
 401
 402     if( depth == CV_32F )
 403     {
 404         for( y = 0; y < size.height; y++ )
 405         {
 406             float* x_data = (float*)(xmat->data.ptr + xmat->step*y);
 407             float* y_data = (float*)(ymat->data.ptr + ymat->step*y);
 408             float* mag_data = mag ? (float*)(mag->data.ptr + mag->step*y) : 0;
 409             float* angle_data = angle ? (float*)(angle->data.ptr + angle->step*y) : 0;
 410
 411             for( x = 0; x < size.width; x += block_size )
 412             {
 413                 int len = MIN( size.width - x, block_size );
 414
 415                 if( mag )
 416                     icvSqrMagnitude_32f( x_data + x, y_data + x, mag_buffer, len );
 417
 418                 if( angle )
 419                 {
 420                     icvFastArctan_32f( y_data + x, x_data + x, angle_data + x, len );
 421                     if( !angle_in_degrees )
 422                         icvScale_32f( angle_data + x, angle_data + x,
 423                                       len, (float)(CV_PI/180.), 0 );
 424                 }
 425
 426                 if( mag )
 427                     icvSqrt_32f( mag_buffer, mag_data + x, len );
 428             }
 429         }
 430     }
 431     else
 432     {
 433         for( y = 0; y < size.height; y++ )
 434         {
 435             double* x_data = (double*)(xmat->data.ptr + xmat->step*y);
 436             double* y_data = (double*)(ymat->data.ptr + ymat->step*y);
 437             double* mag_data = mag ? (double*)(mag->data.ptr + mag->step*y) : 0;
 438             double* angle_data = angle ? (double*)(angle->data.ptr + angle->step*y) : 0;
 439
 440             for( x = 0; x < size.width; x += block_size )
 441             {
 442                 int len = MIN( size.width - x, block_size );
 443
 444                 if( angle )
 445                 {
 446                     icvCvt_64f32f( x_data + x, x_buffer, len );
 447                     icvCvt_64f32f( y_data + x, y_buffer, len );
 448                 }
 449
 450                 if( mag )
 451                 {
 452                     icvSqrMagnitude_64f( x_data + x, y_data + x, mag_data + x, len );
 453                     icvSqrt_64f( mag_data + x, mag_data + x, len );
 454                 }
 455
 456                 if( angle )
 457                 {
 458                     icvFastArctan_32f( y_buffer, x_buffer, x_buffer, len );
 459                     if( !angle_in_degrees )
 460                         icvScale_32f( x_buffer, x_buffer, len, (float)(CV_PI/180.), 0 );
 461                     icvCvt_32f64f( x_buffer, angle_data + x, len );
 462                 }
 463             }
 464         }
 465     }
 466
 467     __END__;
 468 }
 469
 470
 471 /****************************************************************************************\
 472 *                                  Polar -> Cartezian                                    *
 473 \****************************************************************************************/
 474
 475 static CvStatus CV_STDCALL
 476 icvSinCos_32f( const float *angle,float *sinval, float* cosval,
 477                 int len, int angle_in_degrees )
 478 {
 479     const int N = 64;
 480
 481     static const double sin_table[] =
 482     {
 483      0.00000000000000000000,     0.09801714032956060400,
 484      0.19509032201612825000,     0.29028467725446233000,
 485      0.38268343236508978000,     0.47139673682599764000,
 486      0.55557023301960218000,     0.63439328416364549000,
 487      0.70710678118654746000,     0.77301045336273699000,
 488      0.83146961230254524000,     0.88192126434835494000,
 489      0.92387953251128674000,     0.95694033573220894000,
 490      0.98078528040323043000,     0.99518472667219682000,
 491      1.00000000000000000000,     0.99518472667219693000,
 492      0.98078528040323043000,     0.95694033573220894000,
 493      0.92387953251128674000,     0.88192126434835505000,
 494      0.83146961230254546000,     0.77301045336273710000,
 495      0.70710678118654757000,     0.63439328416364549000,
 496      0.55557023301960218000,     0.47139673682599786000,
 497      0.38268343236508989000,     0.29028467725446239000,
 498      0.19509032201612861000,     0.09801714032956082600,
 499      0.00000000000000012246,    -0.09801714032956059000,
 500     -0.19509032201612836000,    -0.29028467725446211000,
 501     -0.38268343236508967000,    -0.47139673682599764000,
 502     -0.55557023301960196000,    -0.63439328416364527000,
 503     -0.70710678118654746000,    -0.77301045336273666000,
 504     -0.83146961230254524000,    -0.88192126434835494000,
 505     -0.92387953251128652000,    -0.95694033573220882000,
 506     -0.98078528040323032000,    -0.99518472667219693000,
 507     -1.00000000000000000000,    -0.99518472667219693000,
 508     -0.98078528040323043000,    -0.95694033573220894000,
 509     -0.92387953251128663000,    -0.88192126434835505000,
 510     -0.83146961230254546000,    -0.77301045336273688000,
 511     -0.70710678118654768000,    -0.63439328416364593000,
 512     -0.55557023301960218000,    -0.47139673682599792000,
 513     -0.38268343236509039000,    -0.29028467725446250000,
 514     -0.19509032201612872000,    -0.09801714032956050600,
 515     };
 516
 517     static const double k2 = (2*CV_PI)/N;
 518
 519     static const double sin_a0 = -0.166630293345647*k2*k2*k2;
 520     static const double sin_a2 = k2;
 521
 522     static const double cos_a0 = -0.499818138450326*k2*k2;
 523     /*static const double cos_a2 =  1;*/
 524
 525     double k1;
 526     int i;
 527
 528     if( !angle_in_degrees )
 529         k1 = N/(2*CV_PI);
 530     else
 531         k1 = N/360.;
 532
 533     for( i = 0; i < len; i++ )
 534     {
 535         double t = angle[i]*k1;
 536         int it = cvRound(t);
 537         t -= it;
 538         int sin_idx = it & (N - 1);
 539         int cos_idx = (N/4 - sin_idx) & (N - 1);
 540
 541         double sin_b = (sin_a0*t*t + sin_a2)*t;
 542         double cos_b = cos_a0*t*t + 1;
 543
 544         double sin_a = sin_table[sin_idx];
 545         double cos_a = sin_table[cos_idx];
 546
 547         double sin_val = sin_a*cos_b + cos_a*sin_b;
 548         double cos_val = cos_a*cos_b - sin_a*sin_b;
 549
 550         sinval[i] = (float)sin_val;
 551         cosval[i] = (float)cos_val;
 552     }
 553
 554     return CV_OK;
 555 }
 556
 557
 558 CV_IMPL void
 559 cvPolarToCart( const CvArr* magarr, const CvArr* anglearr,
 560                CvArr* xarr, CvArr* yarr, int angle_in_degrees )
 561 {
 562     CV_FUNCNAME( "cvPolarToCart" );
 563
 564     __BEGIN__;
 565
 566     float* x_buffer = 0;
 567     float* y_buffer = 0;
 568     int block_size = 0;
 569     CvMat xstub, *xmat = (CvMat*)xarr;
 570     CvMat ystub, *ymat = (CvMat*)yarr;
 571     CvMat magstub, *mag = (CvMat*)magarr;
 572     CvMat anglestub, *angle = (CvMat*)anglearr;
 573     int coi1 = 0, coi2 = 0, coi3 = 0, coi4 = 0;
 574     int depth;
 575     CvSize size;
 576     int x, y;
 577     int cont_flag;
 578
 579     if( !CV_IS_MAT(angle))
 580         CV_CALL( angle = cvGetMat( angle, &anglestub, &coi4 ));
 581
 582     depth = CV_MAT_DEPTH( angle->type );
 583     if( depth < CV_32F )
 584         CV_ERROR( CV_StsUnsupportedFormat, "" );
 585     cont_flag = angle->type;
 586
 587     if( mag )
 588     {
 589         if( !CV_IS_MAT(mag))
 590             CV_CALL( mag = cvGetMat( mag, &magstub, &coi3 ));
 591
 592         if( !CV_ARE_TYPES_EQ( angle, mag ) )
 593             CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
 594
 595         if( !CV_ARE_SIZES_EQ( angle, mag ) )
 596             CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
 597
 598         cont_flag &= mag->type;
 599     }
 600
 601     if( xmat )
 602     {
 603         if( !CV_IS_MAT(xmat))
 604             CV_CALL( xmat = cvGetMat( xmat, &xstub, &coi1 ));
 605
 606         if( !CV_ARE_TYPES_EQ( angle, xmat ) )
 607             CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
 608
 609         if( !CV_ARE_SIZES_EQ( angle, xmat ) )
 610             CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
 611
 612         cont_flag &= xmat->type;
 613     }
 614
 615     if( ymat )
 616     {
 617         if( !CV_IS_MAT(ymat))
 618             CV_CALL( ymat = cvGetMat( ymat, &ystub, &coi2 ));
 619
 620         if( !CV_ARE_TYPES_EQ( angle, ymat ) )
 621             CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
 622
 623         if( !CV_ARE_SIZES_EQ( angle, ymat ) )
 624             CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
 625
 626         cont_flag &= ymat->type;
 627     }
 628
 629     if( coi1 != 0 || coi2 != 0 || coi3 != 0 || coi4 != 0 )
 630         CV_ERROR( CV_BadCOI, "" );
 631
 632     size = cvGetMatSize(angle);
 633     size.width *= CV_MAT_CN(angle->type);
 634
 635     if( CV_IS_MAT_CONT( cont_flag ))
 636     {
 637         size.width *= size.height;
 638         size.height = 1;
 639     }
 640
 641     block_size = MIN( size.width, ICV_MATH_BLOCK_SIZE );
 642     x_buffer = (float*)cvStackAlloc( block_size*sizeof(float));
 643     y_buffer = (float*)cvStackAlloc( block_size*sizeof(float));
 644
 645     if( depth == CV_32F )
 646     {
 647         for( y = 0; y < size.height; y++ )
 648         {
 649             float* x_data = (float*)(xmat ? xmat->data.ptr + xmat->step*y : 0);
 650             float* y_data = (float*)(ymat ? ymat->data.ptr + ymat->step*y : 0);
 651             float* mag_data = (float*)(mag ? mag->data.ptr + mag->step*y : 0);
 652             float* angle_data = (float*)(angle->data.ptr + angle->step*y);
 653
 654             for( x = 0; x < size.width; x += block_size )
 655             {
 656                 int i, len = MIN( size.width - x, block_size );
 657
 658                 icvSinCos_32f( angle_data+x, y_buffer, x_buffer, len, angle_in_degrees );
 659
 660                 for( i = 0; i < len; i++ )
 661                 {
 662                     float tx = x_buffer[i];
 663                     float ty = y_buffer[i];
 664
 665                     if( mag_data )
 666                     {
 667                         float magval = mag_data[x+i];
 668                         tx *= magval;
 669                         ty *= magval;
 670                     }
 671
 672                     if( xmat )
 673                         x_data[x+i] = tx;
 674                     if( ymat )
 675                         y_data[x+i] = ty;
 676                 }
 677             }
 678         }
 679     }
 680     else
 681     {
 682         for( y = 0; y < size.height; y++ )
 683         {
 684             double* x_data = (double*)(xmat ? xmat->data.ptr + xmat->step*y : 0);
 685             double* y_data = (double*)(ymat ? ymat->data.ptr + ymat->step*y : 0);
 686             double* mag_data = (double*)(mag ? mag->data.ptr + mag->step*y : 0);
 687             double* angle_data = (double*)(angle->data.ptr + angle->step*y);
 688             double C = angle_in_degrees ? CV_PI/180. : 1;
 689
 690             for( x = 0; x < size.width; x++ )
 691             {
 692                 double phi = angle_data[x]*C;
 693                 double magval = mag_data ? mag_data[x] : 1.;
 694                 if( xmat )
 695                     x_data[x] = cos(phi)*magval;
 696                 if( ymat )
 697                     y_data[x] = sin(phi)*magval;
 698             }
 699         }
 700     }
 701
 702     __END__;
 703 }
 704
 705 /****************************************************************************************\
 706 *                                          E X P                                         *
 707 \****************************************************************************************/
 708
 709 typedef union
 710 {
 711     struct {
 712 #if ( defined( WORDS_BIGENDIAN ) && !defined( OPENCV_UNIVERSAL_BUILD ) ) || defined( __BIG_ENDIAN__ )
 713         int hi;
 714         int lo;
 715 #else
 716         int lo;
 717         int hi;
 718 #endif
 719     } i;
 720     double d;
 721 }
 722 DBLINT;
 723
 724 #define EXPTAB_SCALE 6
 725 #define EXPTAB_MASK  (1 << EXPTAB_SCALE) - 1
 726
 727 #define EXPPOLY_32F_A0 .9670371139572337719125840413672004409288e-2
 728
 729 static const double icvExpTab[] = {
 730     1.0 * EXPPOLY_32F_A0,
 731     1.0108892860517004600204097905619 * EXPPOLY_32F_A0,
 732     1.0218971486541166782344801347833 * EXPPOLY_32F_A0,
 733     1.0330248790212284225001082839705 * EXPPOLY_32F_A0,
 734     1.0442737824274138403219664787399 * EXPPOLY_32F_A0,
 735     1.0556451783605571588083413251529 * EXPPOLY_32F_A0,
 736     1.0671404006768236181695211209928 * EXPPOLY_32F_A0,
 737     1.0787607977571197937406800374385 * EXPPOLY_32F_A0,
 738     1.0905077326652576592070106557607 * EXPPOLY_32F_A0,
 739     1.1023825833078409435564142094256 * EXPPOLY_32F_A0,
 740     1.1143867425958925363088129569196 * EXPPOLY_32F_A0,
 741     1.126521618608241899794798643787 * EXPPOLY_32F_A0,
 742     1.1387886347566916537038302838415 * EXPPOLY_32F_A0,
 743     1.151189229952982705817759635202 * EXPPOLY_32F_A0,
 744     1.1637248587775775138135735990922 * EXPPOLY_32F_A0,
 745     1.1763969916502812762846457284838 * EXPPOLY_32F_A0,
 746     1.1892071150027210667174999705605 * EXPPOLY_32F_A0,
 747     1.2021567314527031420963969574978 * EXPPOLY_32F_A0,
 748     1.2152473599804688781165202513388 * EXPPOLY_32F_A0,
 749     1.2284805361068700056940089577928 * EXPPOLY_32F_A0,
 750     1.2418578120734840485936774687266 * EXPPOLY_32F_A0,
 751     1.2553807570246910895793906574423 * EXPPOLY_32F_A0,
 752     1.2690509571917332225544190810323 * EXPPOLY_32F_A0,
 753     1.2828700160787782807266697810215 * EXPPOLY_32F_A0,
 754     1.2968395546510096659337541177925 * EXPPOLY_32F_A0,
 755     1.3109612115247643419229917863308 * EXPPOLY_32F_A0,
 756     1.3252366431597412946295370954987 * EXPPOLY_32F_A0,
 757     1.3396675240533030053600306697244 * EXPPOLY_32F_A0,
 758     1.3542555469368927282980147401407 * EXPPOLY_32F_A0,
 759     1.3690024229745906119296011329822 * EXPPOLY_32F_A0,
 760     1.3839098819638319548726595272652 * EXPPOLY_32F_A0,
 761     1.3989796725383111402095281367152 * EXPPOLY_32F_A0,
 762     1.4142135623730950488016887242097 * EXPPOLY_32F_A0,
 763     1.4296133383919700112350657782751 * EXPPOLY_32F_A0,
 764     1.4451808069770466200370062414717 * EXPPOLY_32F_A0,
 765     1.4609177941806469886513028903106 * EXPPOLY_32F_A0,
 766     1.476826145939499311386907480374 * EXPPOLY_32F_A0,
 767     1.4929077282912648492006435314867 * EXPPOLY_32F_A0,
 768     1.5091644275934227397660195510332 * EXPPOLY_32F_A0,
 769     1.5255981507445383068512536895169 * EXPPOLY_32F_A0,
 770     1.5422108254079408236122918620907 * EXPPOLY_32F_A0,
 771     1.5590044002378369670337280894749 * EXPPOLY_32F_A0,
 772     1.5759808451078864864552701601819 * EXPPOLY_32F_A0,
 773     1.5931421513422668979372486431191 * EXPPOLY_32F_A0,
 774     1.6104903319492543081795206673574 * EXPPOLY_32F_A0,
 775     1.628027421857347766848218522014 * EXPPOLY_32F_A0,
 776     1.6457554781539648445187567247258 * EXPPOLY_32F_A0,
 777     1.6636765803267364350463364569764 * EXPPOLY_32F_A0,
 778     1.6817928305074290860622509524664 * EXPPOLY_32F_A0,
 779     1.7001063537185234695013625734975 * EXPPOLY_32F_A0,
 780     1.7186192981224779156293443764563 * EXPPOLY_32F_A0,
 781     1.7373338352737062489942020818722 * EXPPOLY_32F_A0,
 782     1.7562521603732994831121606193753 * EXPPOLY_32F_A0,
 783     1.7753764925265212525505592001993 * EXPPOLY_32F_A0,
 784     1.7947090750031071864277032421278 * EXPPOLY_32F_A0,
 785     1.8142521755003987562498346003623 * EXPPOLY_32F_A0,
 786     1.8340080864093424634870831895883 * EXPPOLY_32F_A0,
 787     1.8539791250833855683924530703377 * EXPPOLY_32F_A0,
 788     1.8741676341102999013299989499544 * EXPPOLY_32F_A0,
 789     1.8945759815869656413402186534269 * EXPPOLY_32F_A0,
 790     1.9152065613971472938726112702958 * EXPPOLY_32F_A0,
 791     1.9360617934922944505980559045667 * EXPPOLY_32F_A0,
 792     1.9571441241754002690183222516269 * EXPPOLY_32F_A0,
 793     1.9784560263879509682582499181312 * EXPPOLY_32F_A0,
 794 };
 795
 796 static const double exp_prescale = 1.4426950408889634073599246810019 * (1 << EXPTAB_SCALE);
 797 static const double exp_postscale = 1./(1 << EXPTAB_SCALE);
 798 static const double exp_max_val = 3000.*(1 << EXPTAB_SCALE); // log10(DBL_MAX) < 3000
 799
 800 IPCVAPI_IMPL( CvStatus, icvExp_32f, ( const float *_x, float *y, int n ), (_x, y, n) )
 801 {
 802     static const double
 803         EXPPOLY_32F_A4 = 1.000000000000002438532970795181890933776 / EXPPOLY_32F_A0,
 804         EXPPOLY_32F_A3 = .6931471805521448196800669615864773144641 / EXPPOLY_32F_A0,
 805         EXPPOLY_32F_A2 = .2402265109513301490103372422686535526573 / EXPPOLY_32F_A0,
 806         EXPPOLY_32F_A1 = .5550339366753125211915322047004666939128e-1 / EXPPOLY_32F_A0;
 807
 808     #undef EXPPOLY
 809     #define EXPPOLY(x)  \
 810         (((((x) + EXPPOLY_32F_A1)*(x) + EXPPOLY_32F_A2)*(x) + EXPPOLY_32F_A3)*(x) + EXPPOLY_32F_A4)
 811
 812     int i = 0;
 813     DBLINT buf[4];
 814     const Cv32suf* x = (const Cv32suf*)_x;
 815
 816     if( !x || !y )
 817         return CV_NULLPTR_ERR;
 818     if( n <= 0 )
 819         return CV_BADSIZE_ERR;
 820
 821     buf[0].i.lo = buf[1].i.lo = buf[2].i.lo = buf[3].i.lo = 0;
 822
 823     for( ; i <= n - 4; i += 4 )
 824     {
 825         double x0 = x[i].f * exp_prescale;
 826         double x1 = x[i + 1].f * exp_prescale;
 827         double x2 = x[i + 2].f * exp_prescale;
 828         double x3 = x[i + 3].f * exp_prescale;
 829         int val0, val1, val2, val3, t;
 830
 831         if( ((x[i].i >> 23) & 255) > 127 + 10 )
 832             x0 = x[i].i < 0 ? -exp_max_val : exp_max_val;
 833
 834         if( ((x[i+1].i >> 23) & 255) > 127 + 10 )
 835             x1 = x[i+1].i < 0 ? -exp_max_val : exp_max_val;
 836
 837         if( ((x[i+2].i >> 23) & 255) > 127 + 10 )
 838             x2 = x[i+2].i < 0 ? -exp_max_val : exp_max_val;
 839
 840         if( ((x[i+3].i >> 23) & 255) > 127 + 10 )
 841             x3 = x[i+3].i < 0 ? -exp_max_val : exp_max_val;
 842
 843         val0 = cvRound(x0);
 844         val1 = cvRound(x1);
 845         val2 = cvRound(x2);
 846         val3 = cvRound(x3);
 847
 848         x0 = (x0 - val0)*exp_postscale;
 849         x1 = (x1 - val1)*exp_postscale;
 850         x2 = (x2 - val2)*exp_postscale;
 851         x3 = (x3 - val3)*exp_postscale;
 852
 853         t = (val0 >> EXPTAB_SCALE) + 1023;
 854         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 855         buf[0].i.hi = t << 20;
 856
 857         t = (val1 >> EXPTAB_SCALE) + 1023;
 858         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 859         buf[1].i.hi = t << 20;
 860
 861         t = (val2 >> EXPTAB_SCALE) + 1023;
 862         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 863         buf[2].i.hi = t << 20;
 864
 865         t = (val3 >> EXPTAB_SCALE) + 1023;
 866         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 867         buf[3].i.hi = t << 20;
 868
 869         x0 = buf[0].d * icvExpTab[val0 & EXPTAB_MASK] * EXPPOLY( x0 );
 870         x1 = buf[1].d * icvExpTab[val1 & EXPTAB_MASK] * EXPPOLY( x1 );
 871
 872         y[i] = (float)x0;
 873         y[i + 1] = (float)x1;
 874
 875         x2 = buf[2].d * icvExpTab[val2 & EXPTAB_MASK] * EXPPOLY( x2 );
 876         x3 = buf[3].d * icvExpTab[val3 & EXPTAB_MASK] * EXPPOLY( x3 );
 877
 878         y[i + 2] = (float)x2;
 879         y[i + 3] = (float)x3;
 880     }
 881
 882     for( ; i < n; i++ )
 883     {
 884         double x0 = x[i].f * exp_prescale;
 885         int val0, t;
 886
 887         if( ((x[i].i >> 23) & 255) > 127 + 10 )
 888             x0 = x[i].i < 0 ? -exp_max_val : exp_max_val;
 889
 890         val0 = cvRound(x0);
 891         t = (val0 >> EXPTAB_SCALE) + 1023;
 892         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 893
 894         buf[0].i.hi = t << 20;
 895         x0 = (x0 - val0)*exp_postscale;
 896
 897         y[i] = (float)(buf[0].d * icvExpTab[val0 & EXPTAB_MASK] * EXPPOLY(x0));
 898     }
 899
 900     return CV_OK;
 901 }
 902
 903
 904 IPCVAPI_IMPL( CvStatus, icvExp_64f, ( const double *_x, double *y, int n ), (_x, y, n) )
 905 {
 906     static const double
 907         A5 = .99999999999999999998285227504999 / EXPPOLY_32F_A0,
 908         A4 = .69314718055994546743029643825322 / EXPPOLY_32F_A0,
 909         A3 = .24022650695886477918181338054308 / EXPPOLY_32F_A0,
 910         A2 = .55504108793649567998466049042729e-1 / EXPPOLY_32F_A0,
 911         A1 = .96180973140732918010002372686186e-2 / EXPPOLY_32F_A0,
 912         A0 = .13369713757180123244806654839424e-2 / EXPPOLY_32F_A0;
 913
 914     #undef EXPPOLY
 915     #define EXPPOLY(x)  (((((A0*(x) + A1)*(x) + A2)*(x) + A3)*(x) + A4)*(x) + A5)
 916
 917     int i = 0;
 918     DBLINT buf[4];
 919     const Cv64suf* x = (const Cv64suf*)_x;
 920
 921     if( !x || !y )
 922         return CV_NULLPTR_ERR;
 923     if( n <= 0 )
 924         return CV_BADSIZE_ERR;
 925
 926     buf[0].i.lo = buf[1].i.lo = buf[2].i.lo = buf[3].i.lo = 0;
 927
 928     for( ; i <= n - 4; i += 4 )
 929     {
 930         double x0 = x[i].f * exp_prescale;
 931         double x1 = x[i + 1].f * exp_prescale;
 932         double x2 = x[i + 2].f * exp_prescale;
 933         double x3 = x[i + 3].f * exp_prescale;
 934
 935         double y0, y1, y2, y3;
 936         int val0, val1, val2, val3, t;
 937
 938         t = (int)(x[i].i >> 52);
 939         if( (t & 2047) > 1023 + 10 )
 940             x0 = t < 0 ? -exp_max_val : exp_max_val;
 941
 942         t = (int)(x[i+1].i >> 52);
 943         if( (t & 2047) > 1023 + 10 )
 944             x1 = t < 0 ? -exp_max_val : exp_max_val;
 945
 946         t = (int)(x[i+2].i >> 52);
 947         if( (t & 2047) > 1023 + 10 )
 948             x2 = t < 0 ? -exp_max_val : exp_max_val;
 949
 950         t = (int)(x[i+3].i >> 52);
 951         if( (t & 2047) > 1023 + 10 )
 952             x3 = t < 0 ? -exp_max_val : exp_max_val;
 953
 954         val0 = cvRound(x0);
 955         val1 = cvRound(x1);
 956         val2 = cvRound(x2);
 957         val3 = cvRound(x3);
 958
 959         x0 = (x0 - val0)*exp_postscale;
 960         x1 = (x1 - val1)*exp_postscale;
 961         x2 = (x2 - val2)*exp_postscale;
 962         x3 = (x3 - val3)*exp_postscale;
 963
 964         t = (val0 >> EXPTAB_SCALE) + 1023;
 965         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 966         buf[0].i.hi = t << 20;
 967
 968         t = (val1 >> EXPTAB_SCALE) + 1023;
 969         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 970         buf[1].i.hi = t << 20;
 971
 972         t = (val2 >> EXPTAB_SCALE) + 1023;
 973         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 974         buf[2].i.hi = t << 20;
 975
 976         t = (val3 >> EXPTAB_SCALE) + 1023;
 977         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 978         buf[3].i.hi = t << 20;
 979
 980         y0 = buf[0].d * icvExpTab[val0 & EXPTAB_MASK] * EXPPOLY( x0 );
 981         y1 = buf[1].d * icvExpTab[val1 & EXPTAB_MASK] * EXPPOLY( x1 );
 982
 983         y[i] = y0;
 984         y[i + 1] = y1;
 985
 986         y2 = buf[2].d * icvExpTab[val2 & EXPTAB_MASK] * EXPPOLY( x2 );
 987         y3 = buf[3].d * icvExpTab[val3 & EXPTAB_MASK] * EXPPOLY( x3 );
 988
 989         y[i + 2] = y2;
 990         y[i + 3] = y3;
 991     }
 992
 993     for( ; i < n; i++ )
 994     {
 995         double x0 = x[i].f * exp_prescale;
 996         int val0, t;
 997
 998         t = (int)(x[i].i >> 52);
 999         if( (t & 2047) > 1023 + 10 )
1000             x0 = t < 0 ? -exp_max_val : exp_max_val;
1001
1002         val0 = cvRound(x0);
1003         t = (val0 >> EXPTAB_SCALE) + 1023;
1004         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
1005
1006         buf[0].i.hi = t << 20;
1007         x0 = (x0 - val0)*exp_postscale;
1008
1009         y[i] = buf[0].d * icvExpTab[val0 & EXPTAB_MASK] * EXPPOLY( x0 );
1010     }
1011
1012     return CV_OK;
1013 }
1014
1015 #undef EXPTAB_SCALE
1016 #undef EXPTAB_MASK
1017 #undef EXPPOLY_32F_A0
1018
1019 CV_IMPL void cvExp( const CvArr* srcarr, CvArr* dstarr )
1020 {
1021     CV_FUNCNAME( "cvExp" );
1022
1023     __BEGIN__;
1024
1025     CvMat srcstub, *src = (CvMat*)srcarr;
1026     CvMat dststub, *dst = (CvMat*)dstarr;
1027     int coi1 = 0, coi2 = 0, src_depth, dst_depth;
1028     double* buffer = 0;
1029     CvSize size;
1030     int x, y, dx = 0;
1031
1032     if( !CV_IS_MAT(src))
1033         CV_CALL( src = cvGetMat( src, &srcstub, &coi1 ));
1034
1035     if( !CV_IS_MAT(dst))
1036         CV_CALL( dst = cvGetMat( dst, &dststub, &coi2 ));
1037
1038     if( coi1 != 0 || coi2 != 0 )
1039         CV_ERROR( CV_BadCOI, "" );
1040
1041     src_depth = CV_MAT_DEPTH(src->type);
1042     dst_depth = CV_MAT_DEPTH(dst->type);
1043
1044     if( !CV_ARE_CNS_EQ( src, dst ) || src_depth < CV_32F || dst_depth < src_depth )
1045         CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
1046
1047     if( !CV_ARE_SIZES_EQ( src, dst ) )
1048         CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
1049
1050     size = cvGetMatSize(src);
1051     size.width *= CV_MAT_CN(src->type);
1052
1053     if( CV_IS_MAT_CONT( src->type & dst->type ))
1054     {
1055         size.width *= size.height;
1056         size.height = 1;
1057     }
1058
1059     if( !CV_ARE_DEPTHS_EQ( src, dst ))
1060     {
1061         dx = MIN( 1024, size.width );
1062         buffer = (double*)cvStackAlloc( dx*sizeof(buffer[0]) );
1063     }
1064
1065     for( y = 0; y < size.height; y++ )
1066     {
1067         uchar* src_data = src->data.ptr + src->step*y;
1068         uchar* dst_data = dst->data.ptr + dst->step*y;
1069
1070         if( src_depth == CV_64F )
1071         {
1072             icvExp_64f( (double*)src_data, (double*)dst_data, size.width );
1073         }
1074         else if( src_depth == dst_depth )
1075         {
1076             icvExp_32f( (float*)src_data, (float*)dst_data, size.width );
1077         }
1078         else
1079         {
1080             for( x = 0; x < size.width; x += dx )
1081             {
1082                 int len = dx;
1083                 if( x + len > size.width )
1084                     len = size.width - x;
1085                 icvCvt_32f64f( (float*)src_data + x, buffer, len );
1086                 icvExp_64f( buffer, (double*)dst_data + x, len );
1087             }
1088         }
1089     }
1090
1091     __END__;
1092 }
1093
1094
1095 /****************************************************************************************\
1096 *                                          L O G                                         *
1097 \****************************************************************************************/
1098
1099 #define LOGTAB_SCALE    8
1100 #define LOGTAB_MASK         ((1 << LOGTAB_SCALE) - 1)
1101 #define LOGTAB_MASK2        ((1 << (20 - LOGTAB_SCALE)) - 1)
1102 #define LOGTAB_MASK2_32F    ((1 << (23 - LOGTAB_SCALE)) - 1)
1103
1104 static const double icvLogTab[] = {
1105 0.0000000000000000000000000000000000000000,    1.000000000000000000000000000000000000000,
1106 .00389864041565732288852075271279318258166,    .9961089494163424124513618677042801556420,
1107 .00778214044205494809292034119607706088573,    .9922480620155038759689922480620155038760,
1108 .01165061721997527263705585198749759001657,    .9884169884169884169884169884169884169884,
1109 .01550418653596525274396267235488267033361,    .9846153846153846153846153846153846153846,
1110 .01934296284313093139406447562578250654042,    .9808429118773946360153256704980842911877,
1111 .02316705928153437593630670221500622574241,    .9770992366412213740458015267175572519084,
1112 .02697658769820207233514075539915211265906,    .9733840304182509505703422053231939163498,
1113 .03077165866675368732785500469617545604706,    .9696969696969696969696969696969696969697,
1114 .03455238150665972812758397481047722976656,    .9660377358490566037735849056603773584906,
1115 .03831886430213659461285757856785494368522,    .9624060150375939849624060150375939849624,
1116 .04207121392068705056921373852674150839447,    .9588014981273408239700374531835205992509,
1117 .04580953603129420126371940114040626212953,    .9552238805970149253731343283582089552239,
1118 .04953393512227662748292900118940451648088,    .9516728624535315985130111524163568773234,
1119 .05324451451881227759255210685296333394944,    .9481481481481481481481481481481481481481,
1120 .05694137640013842427411105973078520037234,    .9446494464944649446494464944649446494465,
1121 .06062462181643483993820353816772694699466,    .9411764705882352941176470588235294117647,
1122 .06429435070539725460836422143984236754475,    .9377289377289377289377289377289377289377,
1123 .06795066190850773679699159401934593915938,    .9343065693430656934306569343065693430657,
1124 .07159365318700880442825962290953611955044,    .9309090909090909090909090909090909090909,
1125 .07522342123758751775142172846244648098944,    .9275362318840579710144927536231884057971,
1126 .07884006170777602129362549021607264876369,    .9241877256317689530685920577617328519856,
1127 .08244366921107458556772229485432035289706,    .9208633093525179856115107913669064748201,
1128 .08603433734180314373940490213499288074675,    .9175627240143369175627240143369175627240,
1129 .08961215868968712416897659522874164395031,    .9142857142857142857142857142857142857143,
1130 .09317722485418328259854092721070628613231,    .9110320284697508896797153024911032028470,
1131 .09672962645855109897752299730200320482256,    .9078014184397163120567375886524822695035,
1132 .10026945316367513738597949668474029749630,    .9045936395759717314487632508833922261484,
1133 .10379679368164355934833764649738441221420,    .9014084507042253521126760563380281690141,
1134 .10731173578908805021914218968959175981580,    .8982456140350877192982456140350877192982,
1135 .11081436634029011301105782649756292812530,    .8951048951048951048951048951048951048951,
1136 .11430477128005862852422325204315711744130,    .8919860627177700348432055749128919860627,
1137 .11778303565638344185817487641543266363440,    .8888888888888888888888888888888888888889,
1138 .12124924363286967987640707633545389398930,    .8858131487889273356401384083044982698962,
1139 .12470347850095722663787967121606925502420,    .8827586206896551724137931034482758620690,
1140 .12814582269193003360996385708858724683530,    .8797250859106529209621993127147766323024,
1141 .13157635778871926146571524895989568904040,    .8767123287671232876712328767123287671233,
1142 .13499516453750481925766280255629681050780,    .8737201365187713310580204778156996587031,
1143 .13840232285911913123754857224412262439730,    .8707482993197278911564625850340136054422,
1144 .14179791186025733629172407290752744302150,    .8677966101694915254237288135593220338983,
1145 .14518200984449788903951628071808954700830,    .8648648648648648648648648648648648648649,
1146 .14855469432313711530824207329715136438610,    .8619528619528619528619528619528619528620,
1147 .15191604202584196858794030049466527998450,    .8590604026845637583892617449664429530201,
1148 .15526612891112392955683674244937719777230,    .8561872909698996655518394648829431438127,
1149 .15860503017663857283636730244325008243330,    .8533333333333333333333333333333333333333,
1150 .16193282026931324346641360989451641216880,    .8504983388704318936877076411960132890365,
1151 .16524957289530714521497145597095368430010,    .8476821192052980132450331125827814569536,
1152 .16855536102980664403538924034364754334090,    .8448844884488448844884488448844884488449,
1153 .17185025692665920060697715143760433420540,    .8421052631578947368421052631578947368421,
1154 .17513433212784912385018287750426679849630,    .8393442622950819672131147540983606557377,
1155 .17840765747281828179637841458315961062910,    .8366013071895424836601307189542483660131,
1156 .18167030310763465639212199675966985523700,    .8338762214983713355048859934853420195440,
1157 .18492233849401198964024217730184318497780,    .8311688311688311688311688311688311688312,
1158 .18816383241818296356839823602058459073300,    .8284789644012944983818770226537216828479,
1159 .19139485299962943898322009772527962923050,    .8258064516129032258064516129032258064516,
1160 .19461546769967164038916962454095482826240,    .8231511254019292604501607717041800643087,
1161 .19782574332991986754137769821682013571260,    .8205128205128205128205128205128205128205,
1162 .20102574606059073203390141770796617493040,    .8178913738019169329073482428115015974441,
1163 .20421554142869088876999228432396193966280,    .8152866242038216560509554140127388535032,
1164 .20739519434607056602715147164417430758480,    .8126984126984126984126984126984126984127,
1165 .21056476910734961416338251183333341032260,    .8101265822784810126582278481012658227848,
1166 .21372432939771812687723695489694364368910,    .8075709779179810725552050473186119873817,
1167 .21687393830061435506806333251006435602900,    .8050314465408805031446540880503144654088,
1168 .22001365830528207823135744547471404075630,    .8025078369905956112852664576802507836991,
1169 .22314355131420973710199007200571941211830,    .8000000000000000000000000000000000000000,
1170 .22626367865045338145790765338460914790630,    .7975077881619937694704049844236760124611,
1171 .22937410106484582006380890106811420992010,    .7950310559006211180124223602484472049689,
1172 .23247487874309405442296849741978803649550,    .7925696594427244582043343653250773993808,
1173 .23556607131276688371634975283086532726890,    .7901234567901234567901234567901234567901,
1174 .23864773785017498464178231643018079921600,    .7876923076923076923076923076923076923077,
1175 .24171993688714515924331749374687206000090,    .7852760736196319018404907975460122699387,
1176 .24478272641769091566565919038112042471760,    .7828746177370030581039755351681957186544,
1177 .24783616390458124145723672882013488560910,    .7804878048780487804878048780487804878049,
1178 .25088030628580937353433455427875742316250,    .7781155015197568389057750759878419452888,
1179 .25391520998096339667426946107298135757450,    .7757575757575757575757575757575757575758,
1180 .25694093089750041913887912414793390780680,    .7734138972809667673716012084592145015106,
1181 .25995752443692604627401010475296061486000,    .7710843373493975903614457831325301204819,
1182 .26296504550088134477547896494797896593800,    .7687687687687687687687687687687687687688,
1183 .26596354849713793599974565040611196309330,    .7664670658682634730538922155688622754491,
1184 .26895308734550393836570947314612567424780,    .7641791044776119402985074626865671641791,
1185 .27193371548364175804834985683555714786050,    .7619047619047619047619047619047619047619,
1186 .27490548587279922676529508862586226314300,    .7596439169139465875370919881305637982196,
1187 .27786845100345625159121709657483734190480,    .7573964497041420118343195266272189349112,
1188 .28082266290088775395616949026589281857030,    .7551622418879056047197640117994100294985,
1189 .28376817313064456316240580235898960381750,    .7529411764705882352941176470588235294118,
1190 .28670503280395426282112225635501090437180,    .7507331378299120234604105571847507331378,
1191 .28963329258304265634293983566749375313530,    .7485380116959064327485380116959064327485,
1192 .29255300268637740579436012922087684273730,    .7463556851311953352769679300291545189504,
1193 .29546421289383584252163927885703742504130,    .7441860465116279069767441860465116279070,
1194 .29836697255179722709783618483925238251680,    .7420289855072463768115942028985507246377,
1195 .30126133057816173455023545102449133992200,    .7398843930635838150289017341040462427746,
1196 .30414733546729666446850615102448500692850,    .7377521613832853025936599423631123919308,
1197 .30702503529491181888388950937951449304830,    .7356321839080459770114942528735632183908,
1198 .30989447772286465854207904158101882785550,    .7335243553008595988538681948424068767908,
1199 .31275571000389684739317885942000430077330,    .7314285714285714285714285714285714285714,
1200 .31560877898630329552176476681779604405180,    .7293447293447293447293447293447293447293,
1201 .31845373111853458869546784626436419785030,    .7272727272727272727272727272727272727273,
1202 .32129061245373424782201254856772720813750,    .7252124645892351274787535410764872521246,
1203 .32411946865421192853773391107097268104550,    .7231638418079096045197740112994350282486,
1204 .32694034499585328257253991068864706903700,    .7211267605633802816901408450704225352113,
1205 .32975328637246797969240219572384376078850,    .7191011235955056179775280898876404494382,
1206 .33255833730007655635318997155991382896900,    .7170868347338935574229691876750700280112,
1207 .33535554192113781191153520921943709254280,    .7150837988826815642458100558659217877095,
1208 .33814494400871636381467055798566434532400,    .7130919220055710306406685236768802228412,
1209 .34092658697059319283795275623560883104800,    .7111111111111111111111111111111111111111,
1210 .34370051385331840121395430287520866841080,    .7091412742382271468144044321329639889197,
1211 .34646676734620857063262633346312213689100,    .7071823204419889502762430939226519337017,
1212 .34922538978528827602332285096053965389730,    .7052341597796143250688705234159779614325,
1213 .35197642315717814209818925519357435405250,    .7032967032967032967032967032967032967033,
1214 .35471990910292899856770532096561510115850,    .7013698630136986301369863013698630136986,
1215 .35745588892180374385176833129662554711100,    .6994535519125683060109289617486338797814,
1216 .36018440357500774995358483465679455548530,    .6975476839237057220708446866485013623978,
1217 .36290549368936841911903457003063522279280,    .6956521739130434782608695652173913043478,
1218 .36561919956096466943762379742111079394830,    .6937669376693766937669376693766937669377,
1219 .36832556115870762614150635272380895912650,    .6918918918918918918918918918918918918919,
1220 .37102461812787262962487488948681857436900,    .6900269541778975741239892183288409703504,
1221 .37371640979358405898480555151763837784530,    .6881720430107526881720430107526881720430,
1222 .37640097516425302659470730759494472295050,    .6863270777479892761394101876675603217158,
1223 .37907835293496944251145919224654790014030,    .6844919786096256684491978609625668449198,
1224 .38174858149084833769393299007788300514230,    .6826666666666666666666666666666666666667,
1225 .38441169891033200034513583887019194662580,    .6808510638297872340425531914893617021277,
1226 .38706774296844825844488013899535872042180,    .6790450928381962864721485411140583554377,
1227 .38971675114002518602873692543653305619950,    .6772486772486772486772486772486772486772,
1228 .39235876060286384303665840889152605086580,    .6754617414248021108179419525065963060686,
1229 .39499380824086893770896722344332374632350,    .6736842105263157894736842105263157894737,
1230 .39762193064713846624158577469643205404280,    .6719160104986876640419947506561679790026,
1231 .40024316412701266276741307592601515352730,    .6701570680628272251308900523560209424084,
1232 .40285754470108348090917615991202183067800,    .6684073107049608355091383812010443864230,
1233 .40546510810816432934799991016916465014230,    .6666666666666666666666666666666666666667,
1234 .40806588980822172674223224930756259709600,    .6649350649350649350649350649350649350649,
1235 .41065992498526837639616360320360399782650,    .6632124352331606217616580310880829015544,
1236 .41324724855021932601317757871584035456180,    .6614987080103359173126614987080103359173,
1237 .41582789514371093497757669865677598863850,    .6597938144329896907216494845360824742268,
1238 .41840189913888381489925905043492093682300,    .6580976863753213367609254498714652956298,
1239 .42096929464412963239894338585145305842150,    .6564102564102564102564102564102564102564,
1240 .42353011550580327293502591601281892508280,    .6547314578005115089514066496163682864450,
1241 .42608439531090003260516141381231136620050,    .6530612244897959183673469387755102040816,
1242 .42863216738969872610098832410585600882780,    .6513994910941475826972010178117048346056,
1243 .43117346481837132143866142541810404509300,    .6497461928934010152284263959390862944162,
1244 .43370832042155937902094819946796633303180,    .6481012658227848101265822784810126582278,
1245 .43623676677491801667585491486534010618930,    .6464646464646464646464646464646464646465,
1246 .43875883620762790027214350629947148263450,    .6448362720403022670025188916876574307305,
1247 .44127456080487520440058801796112675219780,    .6432160804020100502512562814070351758794,
1248 .44378397241030093089975139264424797147500,    .6416040100250626566416040100250626566416,
1249 .44628710262841947420398014401143882423650,    .6400000000000000000000000000000000000000,
1250 .44878398282700665555822183705458883196130,    .6384039900249376558603491271820448877805,
1251 .45127464413945855836729492693848442286250,    .6368159203980099502487562189054726368159,
1252 .45375911746712049854579618113348260521900,    .6352357320099255583126550868486352357320,
1253 .45623743348158757315857769754074979573500,    .6336633663366336633663366336633663366337,
1254 .45870962262697662081833982483658473938700,    .6320987654320987654320987654320987654321,
1255 .46117571512217014895185229761409573256980,    .6305418719211822660098522167487684729064,
1256 .46363574096303250549055974261136725544930,    .6289926289926289926289926289926289926290,
1257 .46608972992459918316399125615134835243230,    .6274509803921568627450980392156862745098,
1258 .46853771156323925639597405279346276074650,    .6259168704156479217603911980440097799511,
1259 .47097971521879100631480241645476780831830,    .6243902439024390243902439024390243902439,
1260 .47341577001667212165614273544633761048330,    .6228710462287104622871046228710462287105,
1261 .47584590486996386493601107758877333253630,    .6213592233009708737864077669902912621359,
1262 .47827014848147025860569669930555392056700,    .6198547215496368038740920096852300242131,
1263 .48068852934575190261057286988943815231330,    .6183574879227053140096618357487922705314,
1264 .48310107575113581113157579238759353756900,    .6168674698795180722891566265060240963855,
1265 .48550781578170076890899053978500887751580,    .6153846153846153846153846153846153846154,
1266 .48790877731923892879351001283794175833480,    .6139088729016786570743405275779376498801,
1267 .49030398804519381705802061333088204264650,    .6124401913875598086124401913875598086124,
1268 .49269347544257524607047571407747454941280,    .6109785202863961813842482100238663484487,
1269 .49507726679785146739476431321236304938800,    .6095238095238095238095238095238095238095,
1270 .49745538920281889838648226032091770321130,    .6080760095011876484560570071258907363420,
1271 .49982786955644931126130359189119189977650,    .6066350710900473933649289099526066350711,
1272 .50219473456671548383667413872899487614650,    .6052009456264775413711583924349881796690,
1273 .50455601075239520092452494282042607665050,    .6037735849056603773584905660377358490566,
1274 .50691172444485432801997148999362252652650,    .6023529411764705882352941176470588235294,
1275 .50926190178980790257412536448100581765150,    .6009389671361502347417840375586854460094,
1276 .51160656874906207391973111953120678663250,    .5995316159250585480093676814988290398126,
1277 .51394575110223428282552049495279788970950,    .5981308411214953271028037383177570093458,
1278 .51627947444845445623684554448118433356300,    .5967365967365967365967365967365967365967,
1279 .51860776420804555186805373523384332656850,    .5953488372093023255813953488372093023256,
1280 .52093064562418522900344441950437612831600,    .5939675174013921113689095127610208816705,
1281 .52324814376454775732838697877014055848100,    .5925925925925925925925925925925925925926,
1282 .52556028352292727401362526507000438869000,    .5912240184757505773672055427251732101617,
1283 .52786708962084227803046587723656557500350,    .5898617511520737327188940092165898617512,
1284 .53016858660912158374145519701414741575700,    .5885057471264367816091954022988505747126,
1285 .53246479886947173376654518506256863474850,    .5871559633027522935779816513761467889908,
1286 .53475575061602764748158733709715306758900,    .5858123569794050343249427917620137299771,
1287 .53704146589688361856929077475797384977350,    .5844748858447488584474885844748858447489,
1288 .53932196859560876944783558428753167390800,    .5831435079726651480637813211845102505695,
1289 .54159728243274429804188230264117009937750,    .5818181818181818181818181818181818181818,
1290 .54386743096728351609669971367111429572100,    .5804988662131519274376417233560090702948,
1291 .54613243759813556721383065450936555862450,    .5791855203619909502262443438914027149321,
1292 .54839232556557315767520321969641372561450,    .5778781038374717832957110609480812641084,
1293 .55064711795266219063194057525834068655950,    .5765765765765765765765765765765765765766,
1294 .55289683768667763352766542084282264113450,    .5752808988764044943820224719101123595506,
1295 .55514150754050151093110798683483153581600,    .5739910313901345291479820627802690582960,
1296 .55738115013400635344709144192165695130850,    .5727069351230425055928411633109619686801,
1297 .55961578793542265941596269840374588966350,    .5714285714285714285714285714285714285714,
1298 .56184544326269181269140062795486301183700,    .5701559020044543429844097995545657015590,
1299 .56407013828480290218436721261241473257550,    .5688888888888888888888888888888888888889,
1300 .56628989502311577464155334382667206227800,    .5676274944567627494456762749445676274945,
1301 .56850473535266865532378233183408156037350,    .5663716814159292035398230088495575221239,
1302 .57071468100347144680739575051120482385150,    .5651214128035320088300220750551876379691,
1303 .57291975356178548306473885531886480748650,    .5638766519823788546255506607929515418502,
1304 .57511997447138785144460371157038025558000,    .5626373626373626373626373626373626373626,
1305 .57731536503482350219940144597785547375700,    .5614035087719298245614035087719298245614,
1306 .57950594641464214795689713355386629700650,    .5601750547045951859956236323851203501094,
1307 .58169173963462239562716149521293118596100,    .5589519650655021834061135371179039301310,
1308 .58387276558098266665552955601015128195300,    .5577342047930283224400871459694989106754,
1309 .58604904500357812846544902640744112432000,    .5565217391304347826086956521739130434783,
1310 .58822059851708596855957011939608491957200,    .5553145336225596529284164859002169197397,
1311 .59038744660217634674381770309992134571100,    .5541125541125541125541125541125541125541,
1312 .59254960960667157898740242671919986605650,    .5529157667386609071274298056155507559395,
1313 .59470710774669277576265358220553025603300,    .5517241379310344827586206896551724137931,
1314 .59685996110779382384237123915227130055450,    .5505376344086021505376344086021505376344,
1315 .59900818964608337768851242799428291618800,    .5493562231759656652360515021459227467811,
1316 .60115181318933474940990890900138765573500,    .5481798715203426124197002141327623126338,
1317 .60329085143808425240052883964381180703650,    .5470085470085470085470085470085470085470,
1318 .60542532396671688843525771517306566238400,    .5458422174840085287846481876332622601279,
1319 .60755525022454170969155029524699784815300,    .5446808510638297872340425531914893617021,
1320 .60968064953685519036241657886421307921400,    .5435244161358811040339702760084925690021,
1321 .61180154110599282990534675263916142284850,    .5423728813559322033898305084745762711864,
1322 .61391794401237043121710712512140162289150,    .5412262156448202959830866807610993657505,
1323 .61602987721551394351138242200249806046500,    .5400843881856540084388185654008438818565,
1324 .61813735955507864705538167982012964785100,    .5389473684210526315789473684210526315789,
1325 .62024040975185745772080281312810257077200,    .5378151260504201680672268907563025210084,
1326 .62233904640877868441606324267922900617100,    .5366876310272536687631027253668763102725,
1327 .62443328801189346144440150965237990021700,    .5355648535564853556485355648535564853556,
1328 .62652315293135274476554741340805776417250,    .5344467640918580375782881002087682672234,
1329 .62860865942237409420556559780379757285100,    .5333333333333333333333333333333333333333,
1330 .63068982562619868570408243613201193511500,    .5322245322245322245322245322245322245322,
1331 .63276666957103777644277897707070223987100,    .5311203319502074688796680497925311203320,
1332 .63483920917301017716738442686619237065300,    .5300207039337474120082815734989648033126,
1333 .63690746223706917739093569252872839570050,    .5289256198347107438016528925619834710744,
1334 .63897144645792069983514238629140891134750,    .5278350515463917525773195876288659793814,
1335 .64103117942093124081992527862894348800200,    .5267489711934156378600823045267489711934,
1336 .64308667860302726193566513757104985415950,    .5256673511293634496919917864476386036961,
1337 .64513796137358470073053240412264131009600,    .5245901639344262295081967213114754098361,
1338 .64718504499530948859131740391603671014300,    .5235173824130879345603271983640081799591,
1339 .64922794662510974195157587018911726772800,    .5224489795918367346938775510204081632653,
1340 .65126668331495807251485530287027359008800,    .5213849287169042769857433808553971486762,
1341 .65330127201274557080523663898929953575150,    .5203252032520325203252032520325203252033,
1342 .65533172956312757406749369692988693714150,    .5192697768762677484787018255578093306288,
1343 .65735807270835999727154330685152672231200,    .5182186234817813765182186234817813765182,
1344 .65938031808912778153342060249997302889800,    .5171717171717171717171717171717171717172,
1345 .66139848224536490484126716182800009846700,    .5161290322580645161290322580645161290323,
1346 .66341258161706617713093692145776003599150,    .5150905432595573440643863179074446680080,
1347 .66542263254509037562201001492212526500250,    .5140562248995983935742971887550200803213,
1348 .66742865127195616370414654738851822912700,    .5130260521042084168336673346693386773547,
1349 .66943065394262923906154583164607174694550,    .5120000000000000000000000000000000000000,
1350 .67142865660530226534774556057527661323550,    .5109780439121756487025948103792415169661,
1351 .67342267521216669923234121597488410770900,    .5099601593625498007968127490039840637450,
1352 .67541272562017662384192817626171745359900,    .5089463220675944333996023856858846918489,
1353 .67739882359180603188519853574689477682100,    .5079365079365079365079365079365079365079,
1354 .67938098479579733801614338517538271844400,    .5069306930693069306930693069306930693069,
1355 .68135922480790300781450241629499942064300,    .5059288537549407114624505928853754940711,
1356 .68333355911162063645036823800182901322850,    .5049309664694280078895463510848126232742,
1357 .68530400309891936760919861626462079584600,    .5039370078740157480314960629921259842520,
1358 .68727057207096020619019327568821609020250,    .5029469548133595284872298624754420432220,
1359 .68923328123880889251040571252815425395950,    .5019607843137254901960784313725490196078,
1360 .69314718055994530941723212145818, 5.0e-01,
1361 };
1362
1363
1364
1365 #define LOGTAB_TRANSLATE(x,h) (((x) - 1.)*icvLogTab[(h)+1])
1366 static const double ln_2 = 0.69314718055994530941723212145818;
1367
1368 IPCVAPI_IMPL( CvStatus, icvLog_32f, ( const float *_x, float *y, int n ), (_x, y, n) )
1369 {
1370     static const double shift[] = { 0, -1./512 };
1371     static const double
1372         A0 = 0.3333333333333333333333333,
1373         A1 = -0.5,
1374         A2 = 1;
1375
1376     #undef LOGPOLY
1377     #define LOGPOLY(x,k) ((x)+=shift[k],((A0*(x) + A1)*(x) + A2)*(x))
1378
1379     int i = 0;
1380     union
1381     {
1382         int i;
1383         float f;
1384     }
1385     buf[4];
1386
1387     const int* x = (const int*)_x;
1388
1389     if( !x || !y )
1390         return CV_NULLPTR_ERR;
1391     if( n <= 0 )
1392         return CV_BADSIZE_ERR;
1393
1394     for( i = 0; i <= n - 4; i += 4 )
1395     {
1396         double x0, x1, x2, x3;
1397         double y0, y1, y2, y3;
1398         int h0, h1, h2, h3;
1399
1400         h0 = x[i];
1401         h1 = x[i+1];
1402         buf[0].i = (h0 & LOGTAB_MASK2_32F) | (127 << 23);
1403         buf[1].i = (h1 & LOGTAB_MASK2_32F) | (127 << 23);
1404
1405         y0 = (((h0 >> 23) & 0xff) - 127) * ln_2;
1406         y1 = (((h1 >> 23) & 0xff) - 127) * ln_2;
1407
1408         h0 = (h0 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1409         h1 = (h1 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1410
1411         y0 += icvLogTab[h0];
1412         y1 += icvLogTab[h1];
1413
1414         h2 = x[i+2];
1415         h3 = x[i+3];
1416
1417         x0 = LOGTAB_TRANSLATE( buf[0].f, h0 );
1418         x1 = LOGTAB_TRANSLATE( buf[1].f, h1 );
1419
1420         buf[2].i = (h2 & LOGTAB_MASK2_32F) | (127 << 23);
1421         buf[3].i = (h3 & LOGTAB_MASK2_32F) | (127 << 23);
1422
1423         y2 = (((h2 >> 23) & 0xff) - 127) * ln_2;
1424         y3 = (((h3 >> 23) & 0xff) - 127) * ln_2;
1425
1426         h2 = (h2 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1427         h3 = (h3 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1428
1429         y2 += icvLogTab[h2];
1430         y3 += icvLogTab[h3];
1431
1432         x2 = LOGTAB_TRANSLATE( buf[2].f, h2 );
1433         x3 = LOGTAB_TRANSLATE( buf[3].f, h3 );
1434
1435         y0 += LOGPOLY( x0, h0 == 510 );
1436         y1 += LOGPOLY( x1, h1 == 510 );
1437
1438         y[i] = (float) y0;
1439         y[i + 1] = (float) y1;
1440
1441         y2 += LOGPOLY( x2, h2 == 510 );
1442         y3 += LOGPOLY( x3, h3 == 510 );
1443
1444         y[i + 2] = (float) y2;
1445         y[i + 3] = (float) y3;
1446     }
1447
1448     for( ; i < n; i++ )
1449     {
1450         int h0 = x[i];
1451         double x0, y0;
1452
1453         y0 = (((h0 >> 23) & 0xff) - 127) * ln_2;
1454
1455         buf[0].i = (h0 & LOGTAB_MASK2_32F) | (127 << 23);
1456         h0 = (h0 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1457
1458         y0 += icvLogTab[h0];
1459         x0 = LOGTAB_TRANSLATE( buf[0].f, h0 );
1460         y0 += LOGPOLY( x0, h0 == 510 );
1461
1462         y[i] = (float)y0;
1463     }
1464
1465     return CV_OK;
1466 }
1467
1468
1469 IPCVAPI_IMPL( CvStatus, icvLog_64f, ( const double *x, double *y, int n ), (x, y, n) )
1470 {
1471     static const double shift[] = { 0, -1./512 };
1472     static const double
1473         A0 = -.1666666666666666666666666666666666666666,
1474         A1 = +0.2,
1475         A2 = -0.25,
1476         A3 = +0.3333333333333333333333333333333333333333,
1477         A4 = -0.5,
1478         A5 = +1.0;
1479
1480     #undef LOGPOLY
1481     #define LOGPOLY(x,k) ((x)+=shift[k], (xq) = (x)*(x),\
1482         ((A0*(xq) + A2)*(xq) + A4)*(xq) + ((A1*(xq) + A3)*(xq) + A5)*(x))
1483
1484     int i = 0;
1485     DBLINT buf[4];
1486     DBLINT *X = (DBLINT *) x;
1487
1488     if( !x || !y )
1489         return CV_NULLPTR_ERR;
1490     if( n <= 0 )
1491         return CV_BADSIZE_ERR;
1492
1493     for( ; i <= n - 4; i += 4 )
1494     {
1495         double xq;
1496         double x0, x1, x2, x3;
1497         double y0, y1, y2, y3;
1498         int h0, h1, h2, h3;
1499
1500         h0 = X[i].i.lo;
1501         h1 = X[i + 1].i.lo;
1502         buf[0].i.lo = h0;
1503         buf[1].i.lo = h1;
1504
1505         h0 = X[i].i.hi;
1506         h1 = X[i + 1].i.hi;
1507         buf[0].i.hi = (h0 & LOGTAB_MASK2) | (1023 << 20);
1508         buf[1].i.hi = (h1 & LOGTAB_MASK2) | (1023 << 20);
1509
1510         y0 = (((h0 >> 20) & 0x7ff) - 1023) * ln_2;
1511         y1 = (((h1 >> 20) & 0x7ff) - 1023) * ln_2;
1512
1513         h2 = X[i + 2].i.lo;
1514         h3 = X[i + 3].i.lo;
1515         buf[2].i.lo = h2;
1516         buf[3].i.lo = h3;
1517
1518         h0 = (h0 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1519         h1 = (h1 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1520
1521         y0 += icvLogTab[h0];
1522         y1 += icvLogTab[h1];
1523
1524         h2 = X[i + 2].i.hi;
1525         h3 = X[i + 3].i.hi;
1526
1527         x0 = LOGTAB_TRANSLATE( buf[0].d, h0 );
1528         x1 = LOGTAB_TRANSLATE( buf[1].d, h1 );
1529
1530         buf[2].i.hi = (h2 & LOGTAB_MASK2) | (1023 << 20);
1531         buf[3].i.hi = (h3 & LOGTAB_MASK2) | (1023 << 20);
1532
1533         y2 = (((h2 >> 20) & 0x7ff) - 1023) * ln_2;
1534         y3 = (((h3 >> 20) & 0x7ff) - 1023) * ln_2;
1535
1536         h2 = (h2 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1537         h3 = (h3 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1538
1539         y2 += icvLogTab[h2];
1540         y3 += icvLogTab[h3];
1541
1542         x2 = LOGTAB_TRANSLATE( buf[2].d, h2 );
1543         x3 = LOGTAB_TRANSLATE( buf[3].d, h3 );
1544
1545         y0 += LOGPOLY( x0, h0 == 510 );
1546         y1 += LOGPOLY( x1, h1 == 510 );
1547
1548         y[i] = y0;
1549         y[i + 1] = y1;
1550
1551         y2 += LOGPOLY( x2, h2 == 510 );
1552         y3 += LOGPOLY( x3, h3 == 510 );
1553
1554         y[i + 2] = y2;
1555         y[i + 3] = y3;
1556     }
1557
1558     for( ; i < n; i++ )
1559     {
1560         int h0 = X[i].i.hi;
1561         double xq;
1562         double x0, y0 = (((h0 >> 20) & 0x7ff) - 1023) * ln_2;
1563
1564         buf[0].i.hi = (h0 & LOGTAB_MASK2) | (1023 << 20);
1565         buf[0].i.lo = X[i].i.lo;
1566         h0 = (h0 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1567
1568         y0 += icvLogTab[h0];
1569         x0 = LOGTAB_TRANSLATE( buf[0].d, h0 );
1570         y0 += LOGPOLY( x0, h0 == 510 );
1571         y[i] = y0;
1572     }
1573
1574     return CV_OK;
1575 }
1576
1577
1578 CV_IMPL void cvLog( const CvArr* srcarr, CvArr* dstarr )
1579 {
1580     CV_FUNCNAME( "cvLog" );
1581
1582     __BEGIN__;
1583
1584     CvMat srcstub, *src = (CvMat*)srcarr;
1585     CvMat dststub, *dst = (CvMat*)dstarr;
1586     int coi1 = 0, coi2 = 0, src_depth, dst_depth;
1587     double* buffer = 0;
1588     CvSize size;
1589     int x, y, dx = 0;
1590
1591     if( !CV_IS_MAT(src))
1592         CV_CALL( src = cvGetMat( src, &srcstub, &coi1 ));
1593
1594     if( !CV_IS_MAT(dst))
1595         CV_CALL( dst = cvGetMat( dst, &dststub, &coi2 ));
1596
1597     if( coi1 != 0 || coi2 != 0 )
1598         CV_ERROR( CV_BadCOI, "" );
1599
1600     src_depth = CV_MAT_DEPTH(src->type);
1601     dst_depth = CV_MAT_DEPTH(dst->type);
1602
1603     if( !CV_ARE_CNS_EQ( src, dst ) || dst_depth < CV_32F || src_depth < dst_depth )
1604         CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
1605
1606     if( !CV_ARE_SIZES_EQ( src, dst ) )
1607         CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
1608
1609     size = cvGetMatSize(src);
1610     size.width *= CV_MAT_CN(src->type);
1611
1612     if( CV_IS_MAT_CONT( src->type & dst->type ))
1613     {
1614         size.width *= size.height;
1615         size.height = 1;
1616     }
1617
1618     if( !CV_ARE_DEPTHS_EQ( src, dst ))
1619     {
1620         dx = MIN( 1024, size.width );
1621         buffer = (double*)cvStackAlloc( dx*sizeof(buffer[0]) );
1622     }
1623
1624     for( y = 0; y < size.height; y++ )
1625     {
1626         uchar* src_data = src->data.ptr + src->step*y;
1627         uchar* dst_data = dst->data.ptr + dst->step*y;
1628
1629         if( dst_depth == CV_64F )
1630         {
1631             icvLog_64f( (double*)src_data, (double*)dst_data, size.width );
1632         }
1633         else if( src_depth == dst_depth )
1634         {
1635             icvLog_32f( (float*)src_data, (float*)dst_data, size.width );
1636         }
1637         else
1638         {
1639             for( x = 0; x < size.width; x += dx )
1640             {
1641                 int len = dx;
1642                 if( x + len > size.width )
1643                     len = size.width - x;
1644                 icvLog_64f( (double*)src_data + x, buffer, len );
1645                 icvCvt_64f32f( buffer, (float*)dst_data + x, len );
1646             }
1647         }
1648     }
1649
1650     __END__;
1651 }
1652
1653
1654 /****************************************************************************************\
1655 *                                    P O W E R                                           *
1656 \****************************************************************************************/
1657
1658 #define ICV_DEF_IPOW_OP( flavor, arrtype, worktype, cast_macro )                    \
1659 static CvStatus CV_STDCALL                                                          \
1660 icvIPow_##flavor( const arrtype* src, arrtype* dst, int len, int power )            \
1661 {                                                                                   \
1662     int i;                                                                          \
1663                                                                                     \
1664     for( i = 0; i < len; i++ )                                                      \
1665     {                                                                               \
1666         worktype a = 1, b = src[i];                                                 \
1667         int p = power;                                                              \
1668         while( p > 1 )                                                              \
1669         {                                                                           \
1670             if( p & 1 )                                                             \
1671                 a *= b;                                                             \
1672             b *= b;                                                                 \
1673             p >>= 1;                                                                \
1674         }                                                                           \
1675                                                                                     \
1676         a *= b;                                                                     \
1677         dst[i] = cast_macro(a);                                                     \
1678     }                                                                               \
1679                                                                                     \
1680     return CV_OK;                                                                   \
1681 }
1682
1683
1684 ICV_DEF_IPOW_OP( 8u, uchar, int, CV_CAST_8U )
1685 ICV_DEF_IPOW_OP( 16u, ushort, int, CV_CAST_16U )
1686 ICV_DEF_IPOW_OP( 16s, short, int, CV_CAST_16S )
1687 ICV_DEF_IPOW_OP( 32s, int, int, CV_CAST_32S )
1688 ICV_DEF_IPOW_OP( 32f, float, double, CV_CAST_32F )
1689 ICV_DEF_IPOW_OP( 64f, double, double, CV_CAST_64F )
1690
1691 #define icvIPow_8s 0
1692
1693 CV_DEF_INIT_FUNC_TAB_1D( IPow )
1694
1695 typedef CvStatus (CV_STDCALL * CvIPowFunc)( const void* src, void* dst, int len, int power );
1696 typedef CvStatus (CV_STDCALL * CvSqrtFunc)( const void* src, void* dst, int len );
1697
1698 CV_IMPL void cvPow( const CvArr* srcarr, CvArr* dstarr, double power )
1699 {
1700     static CvFuncTable ipow_tab;
1701     static int inittab = 0;
1702
1703     CV_FUNCNAME( "cvPow" );
1704
1705     __BEGIN__;
1706
1707     void* temp_buffer = 0;
1708     int block_size = 0;
1709     CvMat srcstub, *src = (CvMat*)srcarr;
1710     CvMat dststub, *dst = (CvMat*)dstarr;
1711     int coi1 = 0, coi2 = 0;
1712     int depth;
1713     CvSize size;
1714     int x, y;
1715     int ipower = cvRound( power );
1716     int is_ipower = 0;
1717
1718     if( !CV_IS_MAT(src))
1719         CV_CALL( src = cvGetMat( src, &srcstub, &coi1 ));
1720
1721     if( !CV_IS_MAT(dst))
1722         CV_CALL( dst = cvGetMat( dst, &dststub, &coi2 ));
1723
1724     if( coi1 != 0 || coi2 != 0 )
1725         CV_ERROR( CV_BadCOI, "" );
1726
1727     if( !CV_ARE_TYPES_EQ( src, dst ))
1728         CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
1729
1730     if( !CV_ARE_SIZES_EQ( src, dst ) )
1731         CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
1732
1733     depth = CV_MAT_DEPTH( src->type );
1734
1735     if( fabs(ipower - power) < DBL_EPSILON )
1736     {
1737         if( !inittab )
1738         {
1739             icvInitIPowTable( &ipow_tab );
1740             inittab = 1;
1741         }
1742
1743         if( ipower < 0 )
1744         {
1745             CV_CALL( cvDiv( 0, src, dst ));
1746
1747             if( ipower == -1 )
1748                 EXIT;
1749             ipower = -ipower;
1750             src = dst;
1751         }
1752
1753         switch( ipower )
1754         {
1755         case 0:
1756             cvSet( dst, cvScalarAll(1));
1757             EXIT;
1758         case 1:
1759             cvCopy( src, dst );
1760             EXIT;
1761         case 2:
1762             cvMul( src, src, dst );
1763             EXIT;
1764         default:
1765             is_ipower = 1;
1766         }
1767     }
1768     else if( depth < CV_32F )
1769         CV_ERROR( CV_StsUnsupportedFormat,
1770         "Fractional or negative integer power factor can be used "
1771         "with floating-point types only");
1772
1773     size = cvGetMatSize(src);
1774     size.width *= CV_MAT_CN(src->type);
1775
1776     if( CV_IS_MAT_CONT( src->type & dst->type ))
1777     {
1778         size.width *= size.height;
1779         size.height = 1;
1780     }
1781
1782     if( is_ipower )
1783     {
1784         CvIPowFunc pow_func = (CvIPowFunc)ipow_tab.fn_2d[depth];
1785         if( !pow_func )
1786             CV_ERROR( CV_StsUnsupportedFormat, "The data type is not supported" );
1787
1788         for( y = 0; y < size.height; y++ )
1789         {
1790             uchar* src_data = src->data.ptr + src->step*y;
1791             uchar* dst_data = dst->data.ptr + dst->step*y;
1792
1793             pow_func( src_data, dst_data, size.width, ipower );
1794         }
1795     }
1796     else if( fabs(fabs(power) - 0.5) < DBL_EPSILON )
1797     {
1798         CvSqrtFunc sqrt_func = power < 0 ?
1799             (depth == CV_32F ? (CvSqrtFunc)icvInvSqrt_32f : (CvSqrtFunc)icvInvSqrt_64f) :
1800             (depth == CV_32F ? (CvSqrtFunc)icvSqrt_32f : (CvSqrtFunc)icvSqrt_64f);
1801
1802         for( y = 0; y < size.height; y++ )
1803         {
1804             uchar* src_data = src->data.ptr + src->step*y;
1805             uchar* dst_data = dst->data.ptr + dst->step*y;
1806
1807             sqrt_func( src_data, dst_data, size.width );
1808         }
1809     }
1810     else
1811     {
1812         block_size = MIN( size.width, ICV_MATH_BLOCK_SIZE );
1813         temp_buffer = cvStackAlloc( block_size*CV_ELEM_SIZE(depth) );
1814
1815         for( y = 0; y < size.height; y++ )
1816         {
1817             uchar* src_data = src->data.ptr + src->step*y;
1818             uchar* dst_data = dst->data.ptr + dst->step*y;
1819
1820             for( x = 0; x < size.width; x += block_size )
1821             {
1822                 int len = MIN( size.width - x, block_size );
1823                 if( depth == CV_32F )
1824                 {
1825                     icvLog_32f( (float*)src_data + x, (float*)temp_buffer, len );
1826                     icvScale_32f( (float*)temp_buffer, (float*)temp_buffer, len, (float)power, 0 );
1827                     icvExp_32f( (float*)temp_buffer, (float*)dst_data + x, len );
1828                 }
1829                 else
1830                 {
1831                     icvLog_64f( (double*)src_data + x, (double*)temp_buffer, len );
1832                     icvScale_64f( (double*)temp_buffer, (double*)temp_buffer, len, power, 0 );
1833                     icvExp_64f( (double*)temp_buffer, (double*)dst_data + x, len );
1834                 }
1835             }
1836         }
1837     }
1838
1839     __END__;
1840 }
1841
1842
1843 /************************** CheckArray for NaN's, Inf's *********************************/
1844
1845 IPCVAPI_IMPL( CvStatus, icvCheckArray_32f_C1R,
1846     ( const float* src, int srcstep, CvSize size, int flags, double min_val, double max_val ),
1847      (src, srcstep, size, flags, min_val, max_val) )
1848 {
1849     Cv32suf a, b;
1850     int ia, ib;
1851     const int* isrc = (const int*)src;
1852
1853     if( !src )
1854         return CV_NULLPTR_ERR;
1855
1856     if( size.width <= 0 || size.height <= 0 )
1857         return CV_BADSIZE_ERR;
1858
1859     if( flags & CV_CHECK_RANGE )
1860     {
1861         a.f = (float)min_val;
1862         b.f = (float)max_val;
1863     }
1864     else
1865     {
1866         a.f = -FLT_MAX;
1867         b.f = FLT_MAX;
1868     }
1869
1870     ia = CV_TOGGLE_FLT(a.i);
1871     ib = CV_TOGGLE_FLT(b.i);
1872
1873     srcstep /= sizeof(isrc[0]);
1874     for( ; size.height--; isrc += srcstep )
1875     {
1876         int i;
1877         for( i = 0; i < size.width; i++ )
1878         {
1879             int val = isrc[i];
1880             val = CV_TOGGLE_FLT(val);
1881
1882             if( val < ia || val >= ib )
1883                 return CV_BADRANGE_ERR;
1884         }
1885     }
1886
1887     return CV_OK;
1888 }
1889
1890
1891 IPCVAPI_IMPL( CvStatus,  icvCheckArray_64f_C1R,
1892     ( const double* src, int srcstep, CvSize size, int flags, double min_val, double max_val ),
1893     (src, srcstep, size, flags, min_val, max_val) )
1894 {
1895     Cv64suf a, b;
1896     int64 ia, ib;
1897     const int64* isrc = (const int64*)src;
1898
1899     if( !src )
1900         return CV_NULLPTR_ERR;
1901
1902     if( size.width <= 0 || size.height <= 0 )
1903         return CV_BADSIZE_ERR;
1904
1905     if( flags & CV_CHECK_RANGE )
1906     {
1907         a.f = min_val;
1908         b.f = max_val;
1909     }
1910     else
1911     {
1912         a.f = -DBL_MAX;
1913         b.f = DBL_MAX;
1914     }
1915
1916     ia = CV_TOGGLE_DBL(a.i);
1917     ib = CV_TOGGLE_DBL(b.i);
1918
1919     srcstep /= sizeof(isrc[0]);
1920     for( ; size.height--; isrc += srcstep )
1921     {
1922         int i;
1923         for( i = 0; i < size.width; i++ )
1924         {
1925             int64 val = isrc[i];
1926             val = CV_TOGGLE_DBL(val);
1927
1928             if( val < ia || val >= ib )
1929                 return CV_BADRANGE_ERR;
1930         }
1931     }
1932
1933     return CV_OK;
1934 }
1935
1936
1937 CV_IMPL  int  cvCheckArr( const CvArr* arr, int flags,
1938                           double minVal, double maxVal )
1939 {
1940     int result = 0;
1941
1942     CV_FUNCNAME( "cvCheckArr" );
1943
1944     __BEGIN__;
1945
1946     if( arr )
1947     {
1948         CvStatus status = CV_OK;
1949         CvMat stub, *mat = (CvMat*)arr;
1950         int type;
1951         CvSize size;
1952
1953         if( !CV_IS_MAT( mat ))
1954             CV_CALL( mat = cvGetMat( mat, &stub, 0, 1 ));
1955
1956         type = CV_MAT_TYPE( mat->type );
1957         size = cvGetMatSize( mat );
1958
1959         size.width *= CV_MAT_CN( type );
1960
1961         if( CV_IS_MAT_CONT( mat->type ))
1962         {
1963             size.width *= size.height;
1964             size.height = 1;
1965         }
1966
1967         if( CV_MAT_DEPTH(type) == CV_32F )
1968         {
1969             status = icvCheckArray_32f_C1R( mat->data.fl, mat->step, size,
1970                                             flags, minVal, maxVal );
1971         }
1972         else if( CV_MAT_DEPTH(type) == CV_64F )
1973         {
1974             status = icvCheckArray_64f_C1R( mat->data.db, mat->step, size,
1975                                             flags, minVal, maxVal );
1976         }
1977         else
1978         {
1979             CV_ERROR( CV_StsUnsupportedFormat, "" );
1980         }
1981
1982         if( status < 0 )
1983         {
1984             if( status != CV_BADRANGE_ERR || !(flags & CV_CHECK_QUIET))
1985                 CV_ERROR( CV_StsOutOfRange, "CheckArray failed" );
1986             EXIT;
1987         }
1988     }
1989
1990     result = 1;
1991
1992     __END__;
1993
1994     return result;
1995 }
1996
1997
1998 /* End of file. */