vcs.maemo.org Git - opencv/blob - src/cxcore/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 //                           License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
  14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
  15 // Third party copyrights are property of their respective owners.
  16 //
  17 // Redistribution and use in source and binary forms, with or without modification,
  18 // are permitted provided that the following conditions are met:
  19 //
  20 //   * Redistribution's of source code must retain the above copyright notice,
  21 //     this list of conditions and the following disclaimer.
  22 //
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,
  24 //     this list of conditions and the following disclaimer in the documentation
  25 //     and/or other materials provided with the distribution.
  26 //
  27 //   * The name of the copyright holders may not be used to endorse or promote products
  28 //     derived from this software without specific prior written permission.
  29 //
  30 // This software is provided by the copyright holders and contributors "as is" and
  31 // any express or implied warranties, including, but not limited to, the implied
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,
  34 // indirect, incidental, special, exemplary, or consequential damages
  35 // (including, but not limited to, procurement of substitute goods or services;
  36 // loss of use, data, or profits; or business interruption) however caused
  37 // and on any theory of liability, whether in contract, strict liability,
  38 // or tort (including negligence or otherwise) arising in any way out of
  39 // the use of this software, even if advised of the possibility of such damage.
  40 //
  41 //M*/
  42
  43 #include "_cxcore.h"
  44
  45 namespace cv
  46 {
  47
  48 static const int MAX_BLOCK_SIZE = 1024;
  49
  50 typedef CvStatus (CV_STDCALL * MathFunc)(const void* src, void* dst, int len);
  51
  52 #define ICV_MATH_BLOCK_SIZE  256
  53
  54 #define _CV_SQRT_MAGIC     0xbe6f0000
  55
  56 #define _CV_SQRT_MAGIC_DBL CV_BIG_UINT(0xbfcd460000000000)
  57
  58 #define _CV_ATAN_CF0  (-15.8131890796f)
  59 #define _CV_ATAN_CF1  (61.0941945596f)
  60 #define _CV_ATAN_CF2  0.f /*(-0.140500406322f)*/
  61
  62 static const float icvAtanTab[8] = { 0.f + _CV_ATAN_CF2, 90.f - _CV_ATAN_CF2,
  63     180.f - _CV_ATAN_CF2, 90.f + _CV_ATAN_CF2,
  64     360.f - _CV_ATAN_CF2, 270.f + _CV_ATAN_CF2,
  65     180.f + _CV_ATAN_CF2, 270.f - _CV_ATAN_CF2
  66 };
  67
  68 static const int icvAtanSign[8] =
  69     { 0, 0x80000000, 0x80000000, 0, 0x80000000, 0, 0, 0x80000000 };
  70
  71 float fastAtan2( float y, float x )
  72 {
  73         double a, x2 = (double)x*x, y2 = (double)y*y;
  74         if( y2 <= x2 )
  75         {
  76                 a = (180./CV_PI)*x*y/(x2 + 0.28*y2 + DBL_EPSILON);
  77                 return (float)(x < 0 ? a + 180 : y >= 0 ? a : 360+a);
  78         }
  79         a = (180./CV_PI)*x*y/(y2 + 0.28*x2 + DBL_EPSILON);
  80         return (float)(y >= 0 ? 90 - a : 270 - a);
  81 }
  82
  83 static CvStatus CV_STDCALL FastAtan2_32f(const float *Y, const float *X, float *angle, int len, bool angleInDegrees=true )
  84 {
  85     if( !Y || !X || !angle || len < 0 )
  86         return CV_BADFACTOR_ERR;
  87
  88         int i = 0;
  89         float scale = angleInDegrees ? (float)(180/CV_PI) : 1.f;
  90
  91 #if CV_SSE2
  92         static const int CV_DECL_ALIGNED(16) iabsmask[] = {0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff};
  93         __m128 eps = _mm_set1_ps((float)DBL_EPSILON), absmask = _mm_load_ps((const float*)iabsmask);
  94         __m128 _90 = _mm_set1_ps((float)(CV_PI*0.5)), _180 = _mm_set1_ps((float)CV_PI), _360 = _mm_set1_ps((float)(CV_PI*2));
  95         __m128 zero = _mm_setzero_ps(), _0_28 = _mm_set1_ps(0.28f), scale4 = _mm_set1_ps(scale);
  96
  97         for( ; i <= len - 4; i += 4 )
  98         {
  99                 __m128 x4 = _mm_loadu_ps(X + i), y4 = _mm_loadu_ps(Y + i);
 100                 __m128 xq4 = _mm_mul_ps(x4, x4), yq4 = _mm_mul_ps(y4, y4);
 101                 __m128 xly = _mm_cmplt_ps(xq4, yq4);
 102                 __m128 z4 = _mm_div_ps(_mm_mul_ps(x4, y4), _mm_add_ps(_mm_add_ps(_mm_max_ps(xq4, yq4),
 103                         _mm_mul_ps(_mm_min_ps(xq4, yq4), _0_28)), eps));
 104
 105                 // a4 <- x < y ? 90 : 0;
 106                 __m128 a4 = _mm_and_ps(xly, _90);
 107                 // a4 <- (y < 0 ? 360 - a4 : a4) == ((x < y ? y < 0 ? 270 : 90) : (y < 0 ? 360 : 0))
 108                 __m128 mask = _mm_cmplt_ps(y4, zero);
 109                 a4 = _mm_or_ps(_mm_and_ps(_mm_sub_ps(_360, a4), mask), _mm_andnot_ps(mask, a4));
 110                 // a4 <- (x < 0 && !(x < y) ? 180 : a4)
 111                 mask = _mm_andnot_ps(xly, _mm_cmplt_ps(x4, zero));
 112                 a4 = _mm_or_ps(_mm_and_ps(_180, mask), _mm_andnot_ps(mask, a4));
 113
 114                 // a4 <- (x < y ? a4 - z4 : a4 + z4)
 115                 a4 = _mm_mul_ps(_mm_add_ps(_mm_xor_ps(z4, _mm_andnot_ps(absmask, xly)), a4), scale4);
 116                 _mm_storeu_ps(angle + i, a4);
 117         }
 118 #endif
 119
 120         for( ; i < len; i++ )
 121         {
 122                 float x = X[i], y = Y[i];
 123                 float a, x2 = x*x, y2 = y*y;
 124                 if( y2 <= x2 )
 125                         a = x*y/(x2 + 0.28f*y2 + (float)DBL_EPSILON) + (float)(x < 0 ? CV_PI : y >= 0 ? 0 : CV_PI*2);
 126                 else
 127                         a = (float)(y >= 0 ? CV_PI*0.5 : CV_PI*1.5) - x*y/(y2 + 0.28f*x2 + (float)DBL_EPSILON);
 128                 angle[i] = a*scale;
 129         }
 130
 131     return CV_OK;
 132 }
 133
 134
 135 /* ************************************************************************** *\
 136    Fast cube root by Ken Turkowski
 137    (http://www.worldserver.com/turk/computergraphics/papers.html)
 138 \* ************************************************************************** */
 139 float  cubeRoot( float value )
 140 {
 141     float fr;
 142     Cv32suf v, m;
 143     int ix, s;
 144     int ex, shx;
 145
 146     v.f = value;
 147     ix = v.i & 0x7fffffff;
 148     s = v.i & 0x80000000;
 149     ex = (ix >> 23) - 127;
 150     shx = ex % 3;
 151     shx -= shx >= 0 ? 3 : 0;
 152     ex = (ex - shx) / 3; /* exponent of cube root */
 153     v.i = (ix & ((1<<23)-1)) | ((shx + 127)<<23);
 154     fr = v.f;
 155
 156     /* 0.125 <= fr < 1.0 */
 157     /* Use quartic rational polynomial with error < 2^(-24) */
 158     fr = (float)(((((45.2548339756803022511987494 * fr +
 159     192.2798368355061050458134625) * fr +
 160     119.1654824285581628956914143) * fr +
 161     13.43250139086239872172837314) * fr +
 162     0.1636161226585754240958355063)/
 163     ((((14.80884093219134573786480845 * fr +
 164     151.9714051044435648658557668) * fr +
 165     168.5254414101568283957668343) * fr +
 166     33.9905941350215598754191872) * fr +
 167     1.0));
 168
 169     /* fr *= 2^ex * sign */
 170     m.f = value;
 171     v.f = fr;
 172     v.i = (v.i + (ex << 23) + s) & (m.i*2 != 0 ? -1 : 0);
 173     return v.f;
 174 }
 175
 176 template<typename T> static CvStatus CV_STDCALL InvSqrt(const T* src, T* dst, int len)
 177 {
 178     for( int i = 0; i < len; i++ )
 179         dst[i] = 1/std::sqrt(src[i]);
 180     return CV_OK;
 181
 182 }
 183
 184 template<typename T> static CvStatus CV_STDCALL Sqrt(const T* src, T* dst, int len)
 185 {
 186     for( int i = 0; i < len; i++ )
 187         dst[i] = std::sqrt(src[i]);
 188     return CV_OK;
 189 }
 190
 191 template<typename T> static CvStatus CV_STDCALL
 192 Magnitude(const T* x, const T* y, T* mag, int len)
 193 {
 194     int i;
 195     for( i = 0; i <= len - 4; i += 4 )
 196     {
 197         T x0 = x[i], y0 = y[i], x1 = x[i+1], y1 = y[i+1];
 198
 199         x0 = x0*x0 + y0*y0;
 200         x1 = x1*x1 + y1*y1;
 201         mag[i] = x0;
 202         mag[i+1] = x1;
 203         x0 = x[i+2], y0 = y[i+2];
 204         x1 = x[i+3], y1 = y[i+3];
 205         x0 = x0*x0 + y0*y0;
 206         x1 = x1*x1 + y1*y1;
 207         mag[i+2] = x0;
 208         mag[i+3] = x1;
 209     }
 210
 211     for( ; i < len; i++ )
 212     {
 213         T x0 = x[i], y0 = y[i];
 214         mag[i] = x0*x0 + y0*y0;
 215     }
 216     Sqrt( mag, mag, len );
 217
 218     return CV_OK;
 219 }
 220
 221
 222 #if CV_SSE2
 223 template<> CvStatus CV_STDCALL InvSqrt(const float* src, float* dst, int len)
 224 {
 225     int i = 0;
 226     __m128 _0_5 = _mm_set1_ps(0.5f), _1_5 = _mm_set1_ps(1.5f);
 227     if( (((size_t)src|(size_t)dst) & 15) == 0 )
 228         for( ; i <= len - 8; i += 8 )
 229         {
 230             __m128 t0 = _mm_load_ps(src + i), t1 = _mm_load_ps(src + i + 4);
 231             __m128 h0 = _mm_mul_ps(t0, _0_5), h1 = _mm_mul_ps(t1, _0_5);
 232             t0 = _mm_rsqrt_ps(t0); t1 = _mm_rsqrt_ps(t1);
 233             t0 = _mm_mul_ps(t0, _mm_sub_ps(_1_5, _mm_mul_ps(_mm_mul_ps(t0,t0),h0)));
 234             t1 = _mm_mul_ps(t1, _mm_sub_ps(_1_5, _mm_mul_ps(_mm_mul_ps(t1,t1),h1)));
 235             _mm_store_ps(dst + i, t0); _mm_store_ps(dst + i + 4, t1);
 236         }
 237     else
 238         for( ; i <= len - 8; i += 8 )
 239         {
 240             __m128 t0 = _mm_loadu_ps(src + i), t1 = _mm_loadu_ps(src + i + 4);
 241             __m128 h0 = _mm_mul_ps(t0, _0_5), h1 = _mm_mul_ps(t1, _0_5);
 242             t0 = _mm_rsqrt_ps(t0); t1 = _mm_rsqrt_ps(t1);
 243             t0 = _mm_mul_ps(t0, _mm_sub_ps(_1_5, _mm_mul_ps(_mm_mul_ps(t0,t0),h0)));
 244             t1 = _mm_mul_ps(t1, _mm_sub_ps(_1_5, _mm_mul_ps(_mm_mul_ps(t1,t1),h1)));
 245             _mm_storeu_ps(dst + i, t0); _mm_storeu_ps(dst + i + 4, t1);
 246         }
 247     for( ; i < len; i++ )
 248         dst[i] = 1/std::sqrt(src[i]);
 249     return CV_OK;
 250 }
 251
 252 template<> CvStatus CV_STDCALL Sqrt(const float* src, float* dst, int len)
 253 {
 254     int i = 0;
 255     if( (((size_t)src|(size_t)dst) & 15) == 0 )
 256         for( ; i <= len - 8; i += 8 )
 257         {
 258             __m128 t0 = _mm_load_ps(src + i), t1 = _mm_load_ps(src + i + 4);
 259             t0 = _mm_sqrt_ps(t0); t1 = _mm_sqrt_ps(t1);
 260             _mm_store_ps(dst + i, t0); _mm_store_ps(dst + i + 4, t1);
 261         }
 262     else
 263         for( ; i <= len - 8; i += 8 )
 264         {
 265             __m128 t0 = _mm_loadu_ps(src + i), t1 = _mm_loadu_ps(src + i + 4);
 266             t0 = _mm_sqrt_ps(t0); t1 = _mm_sqrt_ps(t1);
 267             _mm_storeu_ps(dst + i, t0); _mm_storeu_ps(dst + i + 4, t1);
 268         }
 269     for( ; i < len; i++ )
 270         dst[i] = std::sqrt(src[i]);
 271     return CV_OK;
 272 }
 273
 274 template<> CvStatus CV_STDCALL Sqrt(const double* src, double* dst, int len)
 275 {
 276     int i = 0;
 277     if( (((size_t)src|(size_t)dst) & 15) == 0 )
 278         for( ; i <= len - 4; i += 4 )
 279         {
 280             __m128d t0 = _mm_load_pd(src + i), t1 = _mm_load_pd(src + i + 2);
 281             t0 = _mm_sqrt_pd(t0); t1 = _mm_sqrt_pd(t1);
 282             _mm_store_pd(dst + i, t0); _mm_store_pd(dst + i + 2, t1);
 283         }
 284     else
 285         for( ; i <= len - 4; i += 4 )
 286         {
 287             __m128d t0 = _mm_loadu_pd(src + i), t1 = _mm_loadu_pd(src + i + 2);
 288             t0 = _mm_sqrt_pd(t0); t1 = _mm_sqrt_pd(t1);
 289             _mm_storeu_pd(dst + i, t0); _mm_storeu_pd(dst + i + 2, t1);
 290         }
 291     for( ; i < len; i++ )
 292         dst[i] = std::sqrt(src[i]);
 293     return CV_OK;
 294 }
 295
 296 template<>  CvStatus CV_STDCALL
 297 Magnitude(const float* x, const float* y, float* mag, int len)
 298 {
 299     int i = 0;
 300     for( ; i <= len - 8; i += 8 )
 301     {
 302         __m128 x0 = _mm_loadu_ps(x + i), x1 = _mm_loadu_ps(x + i + 4);
 303         __m128 y0 = _mm_loadu_ps(y + i), y1 = _mm_loadu_ps(y + i + 4);
 304         x0 = _mm_add_ps(_mm_mul_ps(x0, x0), _mm_mul_ps(y0, y0));
 305         x1 = _mm_add_ps(_mm_mul_ps(x1, x1), _mm_mul_ps(y1, y1));
 306         x0 = _mm_sqrt_ps(x0); x1 = _mm_sqrt_ps(x1);
 307         _mm_storeu_ps(mag + i, x0); _mm_storeu_ps(mag + i + 4, x1);
 308     }
 309     for( ; i < len; i++ )
 310     {
 311         float x0 = x[i], y0 = y[i];
 312         mag[i] = std::sqrt(x0*x0 + y0*y0);
 313     }
 314     return CV_OK;
 315 }
 316 #endif
 317
 318 static CvStatus CV_STDCALL Sqrt_32f(const float* src, float* dst, int len)
 319 {
 320     return Sqrt( src, dst, len );
 321 }
 322
 323 static CvStatus CV_STDCALL InvSqrt_32f(const float* src, float* dst, int len)
 324 {
 325     return InvSqrt( src, dst, len );
 326 }
 327
 328 static CvStatus CV_STDCALL Sqrt_64f(const double* src, double* dst, int len)
 329 {
 330     return Sqrt( src, dst, len );
 331 }
 332
 333 static CvStatus CV_STDCALL InvSqrt_64f(const double* src, double* dst, int len)
 334 {
 335     return InvSqrt( src, dst, len );
 336 }
 337
 338
 339 /****************************************************************************************\
 340 *                                  Cartezian -> Polar                                    *
 341 \****************************************************************************************/
 342
 343 void magnitude( const Mat& X, const Mat& Y, Mat& Mag )
 344 {
 345         int type = X.type(), depth = X.depth(), cn = X.channels();
 346         CV_Assert( X.size() == Y.size() && type == Y.type() && (depth == CV_32F || depth == CV_64F));
 347     Mag.create( X.size(), type );
 348
 349     Size size = getContinuousSize( X, Y, Mag, cn );
 350
 351     if( depth == CV_32F )
 352     {
 353         const float *x = (const float*)X.data, *y = (const float*)Y.data;
 354         float *mag = (float*)Mag.data;
 355         size_t xstep = X.step/sizeof(x[0]), ystep = Y.step/sizeof(y[0]);
 356         size_t mstep = Mag.step/sizeof(mag[0]);
 357
 358         for( ; size.height--; x += xstep, y += ystep, mag += mstep )
 359             Magnitude( x, y, mag, size.width );
 360     }
 361     else
 362     {
 363         const double *x = (const double*)X.data, *y = (const double*)Y.data;
 364         double *mag = (double*)Mag.data;
 365         size_t xstep = X.step/sizeof(x[0]), ystep = Y.step/sizeof(y[0]);
 366         size_t mstep = Mag.step/sizeof(mag[0]);
 367
 368         for( ; size.height--; x += xstep, y += ystep, mag += mstep )
 369             Magnitude( x, y, mag, size.width );
 370     }
 371 }
 372
 373 void phase( const Mat& X, const Mat& Y, Mat& Angle, bool angleInDegrees )
 374 {
 375     float buf[2][MAX_BLOCK_SIZE];
 376     int i, j, type = X.type(), depth = X.depth(), cn = X.channels();
 377
 378     CV_Assert( X.size() == Y.size() && type == Y.type() && (depth == CV_32F || depth == CV_64F));
 379     Angle.create( X.size(), type );
 380
 381     Size size = getContinuousSize( X, Y, Angle, cn );
 382
 383     if( depth == CV_32F )
 384     {
 385         const float *x = (const float*)X.data, *y = (const float*)Y.data;
 386         float *angle = (float*)Angle.data;
 387         size_t xstep = X.step/sizeof(x[0]), ystep = Y.step/sizeof(y[0]);
 388         size_t astep = Angle.step/sizeof(angle[0]);
 389
 390         for( ; size.height--; x += xstep, y += ystep, angle += astep )
 391             FastAtan2_32f( y, x, angle, size.width, angleInDegrees );
 392     }
 393     else
 394     {
 395         const double *x = (const double*)X.data, *y = (const double*)Y.data;
 396         double *angle = (double*)Angle.data;
 397         size_t xstep = X.step/sizeof(x[0]), ystep = Y.step/sizeof(y[0]);
 398         size_t astep = Angle.step/sizeof(angle[0]);
 399
 400         for( ; size.height--; x += xstep, y += ystep, angle += astep )
 401         {
 402             for( i = 0; i < size.width; i += MAX_BLOCK_SIZE )
 403             {
 404                 int block_size = std::min(MAX_BLOCK_SIZE, size.width - i);
 405                 for( j = 0; j < block_size; j++ )
 406                 {
 407                     buf[0][j] = (float)x[i + j];
 408                     buf[1][j] = (float)y[i + j];
 409                 }
 410                 FastAtan2_32f( buf[1], buf[0], buf[0], block_size, angleInDegrees );
 411                 for( j = 0; j < block_size; j++ )
 412                                         angle[i + j] = buf[0][j];
 413             }
 414         }
 415     }
 416 }
 417
 418 void cartToPolar( const Mat& X, const Mat& Y, Mat& Mag, Mat& Angle, bool angleInDegrees )
 419 {
 420     float buf[2][MAX_BLOCK_SIZE];
 421     int i, j, type = X.type(), depth = X.depth(), cn = X.channels();
 422
 423     CV_Assert( X.size() == Y.size() && type == Y.type() && (depth == CV_32F || depth == CV_64F));
 424     Mag.create( X.size(), type );
 425     Angle.create( X.size(), type );
 426
 427     Size size = getContinuousSize( X, Y, Mag, Angle, cn );
 428     bool inplace = Mag.data == X.data || Mag.data == Y.data;
 429
 430     if( depth == CV_32F )
 431     {
 432         const float *x = (const float*)X.data, *y = (const float*)Y.data;
 433         float *mag = (float*)Mag.data, *angle = (float*)Angle.data;
 434         size_t xstep = X.step/sizeof(x[0]), ystep = Y.step/sizeof(y[0]);
 435         size_t mstep = Mag.step/sizeof(mag[0]), astep = Angle.step/sizeof(angle[0]);
 436
 437         for( ; size.height--; x += xstep, y += ystep, mag += mstep, angle += astep )
 438         {
 439             for( i = 0; i < size.width; i += MAX_BLOCK_SIZE )
 440             {
 441                 int block_size = std::min(MAX_BLOCK_SIZE, size.width - i);
 442                 Magnitude( x + i, y + i, inplace ? buf[0] : mag + i, block_size );
 443                 FastAtan2_32f( y + i, x + i, angle + i, block_size, angleInDegrees );
 444                 if( inplace )
 445                     for( j = 0; j < block_size; j++ )
 446                         mag[i + j] = buf[0][j];
 447             }
 448         }
 449     }
 450     else
 451     {
 452         const double *x = (const double*)X.data, *y = (const double*)Y.data;
 453         double *mag = (double*)Mag.data, *angle = (double*)Angle.data;
 454         size_t xstep = X.step/sizeof(x[0]), ystep = Y.step/sizeof(y[0]);
 455         size_t mstep = Mag.step/sizeof(mag[0]), astep = Angle.step/sizeof(angle[0]);
 456
 457         for( ; size.height--; x += xstep, y += ystep, mag += mstep, angle += astep )
 458         {
 459             for( i = 0; i < size.width; i += MAX_BLOCK_SIZE )
 460             {
 461                 int block_size = std::min(MAX_BLOCK_SIZE, size.width - i);
 462                 for( j = 0; j < block_size; j++ )
 463                 {
 464                     buf[0][j] = (float)x[i + j];
 465                     buf[1][j] = (float)y[i + j];
 466                 }
 467                 FastAtan2_32f( buf[1], buf[0], buf[0], block_size, angleInDegrees );
 468                 Magnitude( x + i, y + i, mag + i, block_size );
 469                                 for( j = 0; j < block_size; j++ )
 470                                         angle[i + j] = buf[0][j];
 471             }
 472         }
 473     }
 474 }
 475
 476
 477 /****************************************************************************************\
 478 *                                  Polar -> Cartezian                                    *
 479 \****************************************************************************************/
 480
 481 static CvStatus CV_STDCALL
 482 SinCos_32f( const float *angle,float *sinval, float* cosval,
 483             int len, int angle_in_degrees )
 484 {
 485     const int N = 64;
 486
 487     static const double sin_table[] =
 488     {
 489      0.00000000000000000000,     0.09801714032956060400,
 490      0.19509032201612825000,     0.29028467725446233000,
 491      0.38268343236508978000,     0.47139673682599764000,
 492      0.55557023301960218000,     0.63439328416364549000,
 493      0.70710678118654746000,     0.77301045336273699000,
 494      0.83146961230254524000,     0.88192126434835494000,
 495      0.92387953251128674000,     0.95694033573220894000,
 496      0.98078528040323043000,     0.99518472667219682000,
 497      1.00000000000000000000,     0.99518472667219693000,
 498      0.98078528040323043000,     0.95694033573220894000,
 499      0.92387953251128674000,     0.88192126434835505000,
 500      0.83146961230254546000,     0.77301045336273710000,
 501      0.70710678118654757000,     0.63439328416364549000,
 502      0.55557023301960218000,     0.47139673682599786000,
 503      0.38268343236508989000,     0.29028467725446239000,
 504      0.19509032201612861000,     0.09801714032956082600,
 505      0.00000000000000012246,    -0.09801714032956059000,
 506     -0.19509032201612836000,    -0.29028467725446211000,
 507     -0.38268343236508967000,    -0.47139673682599764000,
 508     -0.55557023301960196000,    -0.63439328416364527000,
 509     -0.70710678118654746000,    -0.77301045336273666000,
 510     -0.83146961230254524000,    -0.88192126434835494000,
 511     -0.92387953251128652000,    -0.95694033573220882000,
 512     -0.98078528040323032000,    -0.99518472667219693000,
 513     -1.00000000000000000000,    -0.99518472667219693000,
 514     -0.98078528040323043000,    -0.95694033573220894000,
 515     -0.92387953251128663000,    -0.88192126434835505000,
 516     -0.83146961230254546000,    -0.77301045336273688000,
 517     -0.70710678118654768000,    -0.63439328416364593000,
 518     -0.55557023301960218000,    -0.47139673682599792000,
 519     -0.38268343236509039000,    -0.29028467725446250000,
 520     -0.19509032201612872000,    -0.09801714032956050600,
 521     };
 522
 523     static const double k2 = (2*CV_PI)/N;
 524
 525     static const double sin_a0 = -0.166630293345647*k2*k2*k2;
 526     static const double sin_a2 = k2;
 527
 528     static const double cos_a0 = -0.499818138450326*k2*k2;
 529     /*static const double cos_a2 =  1;*/
 530
 531     double k1;
 532     int i;
 533
 534     if( !angle_in_degrees )
 535         k1 = N/(2*CV_PI);
 536     else
 537         k1 = N/360.;
 538
 539     for( i = 0; i < len; i++ )
 540     {
 541         double t = angle[i]*k1;
 542         int it = cvRound(t);
 543         t -= it;
 544         int sin_idx = it & (N - 1);
 545         int cos_idx = (N/4 - sin_idx) & (N - 1);
 546
 547         double sin_b = (sin_a0*t*t + sin_a2)*t;
 548         double cos_b = cos_a0*t*t + 1;
 549
 550         double sin_a = sin_table[sin_idx];
 551         double cos_a = sin_table[cos_idx];
 552
 553         double sin_val = sin_a*cos_b + cos_a*sin_b;
 554         double cos_val = cos_a*cos_b - sin_a*sin_b;
 555
 556         sinval[i] = (float)sin_val;
 557         cosval[i] = (float)cos_val;
 558     }
 559
 560     return CV_OK;
 561 }
 562
 563
 564 void polarToCart( const Mat& Mag, const Mat& Angle, Mat& X, Mat& Y, bool angleInDegrees )
 565 {
 566     int i, j, type = Angle.type(), depth = Angle.depth();
 567     Size size;
 568
 569     CV_Assert( depth == CV_32F || depth == CV_64F );
 570     X.create( Angle.size(), type );
 571     Y.create( Angle.size(), type );
 572
 573     if( Mag.data )
 574     {
 575         CV_Assert( Mag.size() == Angle.size() && Mag.type() == Angle.type() );
 576         size = getContinuousSize( Mag, Angle, X, Y, Angle.channels() );
 577     }
 578     else
 579         size = getContinuousSize( Angle, X, Y, Angle.channels() );
 580
 581     if( depth == CV_32F )
 582     {
 583         float *x = (float*)X.data, *y = (float*)Y.data;
 584         const float *mag = (const float*)Mag.data, *angle = (const float*)Angle.data;
 585         size_t xstep = X.step/sizeof(x[0]), ystep = Y.step/sizeof(y[0]);
 586         size_t mstep = Mag.step/sizeof(mag[0]), astep = Angle.step/sizeof(angle[0]);
 587
 588         for( ; size.height--; x += xstep, y += ystep, mag += mstep, angle += astep )
 589         {
 590             for( i = 0; i < size.width; i += MAX_BLOCK_SIZE )
 591             {
 592                 int block_size = std::min(MAX_BLOCK_SIZE, size.width - i);
 593                 SinCos_32f( angle + i, y + i, x + i, block_size, angleInDegrees );
 594                 for( j = 0; j < block_size; j++ )
 595                 {
 596                     float m = mag ? mag[i + j] : 1.f;
 597                     float t0 = x[i + j]*m, t1 = y[i + j]*m;
 598                     x[i + j] = t0; y[i + j] = t1;
 599                 }
 600             }
 601         }
 602     }
 603     else
 604     {
 605         double *x = (double*)X.data, *y = (double*)Y.data;
 606         const double *mag = (const double*)Mag.data, *angle = (const double*)Angle.data;
 607         size_t xstep = X.step/sizeof(x[0]), ystep = Y.step/sizeof(y[0]);
 608         size_t mstep = Mag.step/sizeof(mag[0]), astep = Angle.step/sizeof(angle[0]);
 609         double ascale = angleInDegrees ? CV_PI/180. : 1;
 610
 611         for( ; size.height--; x += xstep, y += ystep, mag += mstep, angle += astep )
 612         {
 613             for( j = 0; j < size.width; j++ )
 614             {
 615                 double alpha = angle[j]*ascale, m = mag ? mag[j] : 1.;
 616                 double a = cos(alpha), b = sin(alpha);
 617                 x[j] = m*a; y[j] = m*b;
 618             }
 619         }
 620     }
 621 }
 622
 623 /****************************************************************************************\
 624 *                                          E X P                                         *
 625 \****************************************************************************************/
 626
 627 typedef union
 628 {
 629     struct {
 630 #if ( defined( WORDS_BIGENDIAN ) && !defined( OPENCV_UNIVERSAL_BUILD ) ) || defined( __BIG_ENDIAN__ )
 631         int hi;
 632         int lo;
 633 #else
 634         int lo;
 635         int hi;
 636 #endif
 637     } i;
 638     double d;
 639 }
 640 DBLINT;
 641
 642 #ifndef HAVE_IPP
 643
 644 #define EXPTAB_SCALE 6
 645 #define EXPTAB_MASK  ((1 << EXPTAB_SCALE) - 1)
 646
 647 #define EXPPOLY_32F_A0 .9670371139572337719125840413672004409288e-2
 648
 649 static const double expTab[] = {
 650     1.0 * EXPPOLY_32F_A0,
 651     1.0108892860517004600204097905619 * EXPPOLY_32F_A0,
 652     1.0218971486541166782344801347833 * EXPPOLY_32F_A0,
 653     1.0330248790212284225001082839705 * EXPPOLY_32F_A0,
 654     1.0442737824274138403219664787399 * EXPPOLY_32F_A0,
 655     1.0556451783605571588083413251529 * EXPPOLY_32F_A0,
 656     1.0671404006768236181695211209928 * EXPPOLY_32F_A0,
 657     1.0787607977571197937406800374385 * EXPPOLY_32F_A0,
 658     1.0905077326652576592070106557607 * EXPPOLY_32F_A0,
 659     1.1023825833078409435564142094256 * EXPPOLY_32F_A0,
 660     1.1143867425958925363088129569196 * EXPPOLY_32F_A0,
 661     1.126521618608241899794798643787 * EXPPOLY_32F_A0,
 662     1.1387886347566916537038302838415 * EXPPOLY_32F_A0,
 663     1.151189229952982705817759635202 * EXPPOLY_32F_A0,
 664     1.1637248587775775138135735990922 * EXPPOLY_32F_A0,
 665     1.1763969916502812762846457284838 * EXPPOLY_32F_A0,
 666     1.1892071150027210667174999705605 * EXPPOLY_32F_A0,
 667     1.2021567314527031420963969574978 * EXPPOLY_32F_A0,
 668     1.2152473599804688781165202513388 * EXPPOLY_32F_A0,
 669     1.2284805361068700056940089577928 * EXPPOLY_32F_A0,
 670     1.2418578120734840485936774687266 * EXPPOLY_32F_A0,
 671     1.2553807570246910895793906574423 * EXPPOLY_32F_A0,
 672     1.2690509571917332225544190810323 * EXPPOLY_32F_A0,
 673     1.2828700160787782807266697810215 * EXPPOLY_32F_A0,
 674     1.2968395546510096659337541177925 * EXPPOLY_32F_A0,
 675     1.3109612115247643419229917863308 * EXPPOLY_32F_A0,
 676     1.3252366431597412946295370954987 * EXPPOLY_32F_A0,
 677     1.3396675240533030053600306697244 * EXPPOLY_32F_A0,
 678     1.3542555469368927282980147401407 * EXPPOLY_32F_A0,
 679     1.3690024229745906119296011329822 * EXPPOLY_32F_A0,
 680     1.3839098819638319548726595272652 * EXPPOLY_32F_A0,
 681     1.3989796725383111402095281367152 * EXPPOLY_32F_A0,
 682     1.4142135623730950488016887242097 * EXPPOLY_32F_A0,
 683     1.4296133383919700112350657782751 * EXPPOLY_32F_A0,
 684     1.4451808069770466200370062414717 * EXPPOLY_32F_A0,
 685     1.4609177941806469886513028903106 * EXPPOLY_32F_A0,
 686     1.476826145939499311386907480374 * EXPPOLY_32F_A0,
 687     1.4929077282912648492006435314867 * EXPPOLY_32F_A0,
 688     1.5091644275934227397660195510332 * EXPPOLY_32F_A0,
 689     1.5255981507445383068512536895169 * EXPPOLY_32F_A0,
 690     1.5422108254079408236122918620907 * EXPPOLY_32F_A0,
 691     1.5590044002378369670337280894749 * EXPPOLY_32F_A0,
 692     1.5759808451078864864552701601819 * EXPPOLY_32F_A0,
 693     1.5931421513422668979372486431191 * EXPPOLY_32F_A0,
 694     1.6104903319492543081795206673574 * EXPPOLY_32F_A0,
 695     1.628027421857347766848218522014 * EXPPOLY_32F_A0,
 696     1.6457554781539648445187567247258 * EXPPOLY_32F_A0,
 697     1.6636765803267364350463364569764 * EXPPOLY_32F_A0,
 698     1.6817928305074290860622509524664 * EXPPOLY_32F_A0,
 699     1.7001063537185234695013625734975 * EXPPOLY_32F_A0,
 700     1.7186192981224779156293443764563 * EXPPOLY_32F_A0,
 701     1.7373338352737062489942020818722 * EXPPOLY_32F_A0,
 702     1.7562521603732994831121606193753 * EXPPOLY_32F_A0,
 703     1.7753764925265212525505592001993 * EXPPOLY_32F_A0,
 704     1.7947090750031071864277032421278 * EXPPOLY_32F_A0,
 705     1.8142521755003987562498346003623 * EXPPOLY_32F_A0,
 706     1.8340080864093424634870831895883 * EXPPOLY_32F_A0,
 707     1.8539791250833855683924530703377 * EXPPOLY_32F_A0,
 708     1.8741676341102999013299989499544 * EXPPOLY_32F_A0,
 709     1.8945759815869656413402186534269 * EXPPOLY_32F_A0,
 710     1.9152065613971472938726112702958 * EXPPOLY_32F_A0,
 711     1.9360617934922944505980559045667 * EXPPOLY_32F_A0,
 712     1.9571441241754002690183222516269 * EXPPOLY_32F_A0,
 713     1.9784560263879509682582499181312 * EXPPOLY_32F_A0,
 714 };
 715
 716 static const double exp_prescale = 1.4426950408889634073599246810019 * (1 << EXPTAB_SCALE);
 717 static const double exp_postscale = 1./(1 << EXPTAB_SCALE);
 718 static const double exp_max_val = 3000.*(1 << EXPTAB_SCALE); // log10(DBL_MAX) < 3000
 719
 720 static CvStatus CV_STDCALL Exp_32f( const float *_x, float *y, int n )
 721 {
 722     static const double
 723         EXPPOLY_32F_A4 = 1.000000000000002438532970795181890933776 / EXPPOLY_32F_A0,
 724         EXPPOLY_32F_A3 = .6931471805521448196800669615864773144641 / EXPPOLY_32F_A0,
 725         EXPPOLY_32F_A2 = .2402265109513301490103372422686535526573 / EXPPOLY_32F_A0,
 726         EXPPOLY_32F_A1 = .5550339366753125211915322047004666939128e-1 / EXPPOLY_32F_A0;
 727
 728     #undef EXPPOLY
 729     #define EXPPOLY(x)  \
 730         (((((x) + EXPPOLY_32F_A1)*(x) + EXPPOLY_32F_A2)*(x) + EXPPOLY_32F_A3)*(x) + EXPPOLY_32F_A4)
 731
 732     int i = 0;
 733     DBLINT buf[4];
 734     const Cv32suf* x = (const Cv32suf*)_x;
 735
 736     if( !x || !y )
 737         return CV_NULLPTR_ERR;
 738     if( n <= 0 )
 739         return CV_BADSIZE_ERR;
 740
 741     buf[0].i.lo = buf[1].i.lo = buf[2].i.lo = buf[3].i.lo = 0;
 742
 743     for( ; i <= n - 4; i += 4 )
 744     {
 745         double x0 = x[i].f * exp_prescale;
 746         double x1 = x[i + 1].f * exp_prescale;
 747         double x2 = x[i + 2].f * exp_prescale;
 748         double x3 = x[i + 3].f * exp_prescale;
 749         int val0, val1, val2, val3, t;
 750
 751         if( ((x[i].i >> 23) & 255) > 127 + 10 )
 752             x0 = x[i].i < 0 ? -exp_max_val : exp_max_val;
 753
 754         if( ((x[i+1].i >> 23) & 255) > 127 + 10 )
 755             x1 = x[i+1].i < 0 ? -exp_max_val : exp_max_val;
 756
 757         if( ((x[i+2].i >> 23) & 255) > 127 + 10 )
 758             x2 = x[i+2].i < 0 ? -exp_max_val : exp_max_val;
 759
 760         if( ((x[i+3].i >> 23) & 255) > 127 + 10 )
 761             x3 = x[i+3].i < 0 ? -exp_max_val : exp_max_val;
 762
 763         val0 = cvRound(x0);
 764         val1 = cvRound(x1);
 765         val2 = cvRound(x2);
 766         val3 = cvRound(x3);
 767
 768         x0 = (x0 - val0)*exp_postscale;
 769         x1 = (x1 - val1)*exp_postscale;
 770         x2 = (x2 - val2)*exp_postscale;
 771         x3 = (x3 - val3)*exp_postscale;
 772
 773         t = (val0 >> EXPTAB_SCALE) + 1023;
 774         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 775         buf[0].i.hi = t << 20;
 776
 777         t = (val1 >> EXPTAB_SCALE) + 1023;
 778         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 779         buf[1].i.hi = t << 20;
 780
 781         t = (val2 >> EXPTAB_SCALE) + 1023;
 782         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 783         buf[2].i.hi = t << 20;
 784
 785         t = (val3 >> EXPTAB_SCALE) + 1023;
 786         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 787         buf[3].i.hi = t << 20;
 788
 789         x0 = buf[0].d * expTab[val0 & EXPTAB_MASK] * EXPPOLY( x0 );
 790         x1 = buf[1].d * expTab[val1 & EXPTAB_MASK] * EXPPOLY( x1 );
 791
 792         y[i] = (float)x0;
 793         y[i + 1] = (float)x1;
 794
 795         x2 = buf[2].d * expTab[val2 & EXPTAB_MASK] * EXPPOLY( x2 );
 796         x3 = buf[3].d * expTab[val3 & EXPTAB_MASK] * EXPPOLY( x3 );
 797
 798         y[i + 2] = (float)x2;
 799         y[i + 3] = (float)x3;
 800     }
 801
 802     for( ; i < n; i++ )
 803     {
 804         double x0 = x[i].f * exp_prescale;
 805         int val0, t;
 806
 807         if( ((x[i].i >> 23) & 255) > 127 + 10 )
 808             x0 = x[i].i < 0 ? -exp_max_val : exp_max_val;
 809
 810         val0 = cvRound(x0);
 811         t = (val0 >> EXPTAB_SCALE) + 1023;
 812         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 813
 814         buf[0].i.hi = t << 20;
 815         x0 = (x0 - val0)*exp_postscale;
 816
 817         y[i] = (float)(buf[0].d * expTab[val0 & EXPTAB_MASK] * EXPPOLY(x0));
 818     }
 819
 820     return CV_OK;
 821 }
 822
 823
 824 static CvStatus CV_STDCALL Exp_64f( const double *_x, double *y, int n )
 825 {
 826     static const double
 827         A5 = .99999999999999999998285227504999 / EXPPOLY_32F_A0,
 828         A4 = .69314718055994546743029643825322 / EXPPOLY_32F_A0,
 829         A3 = .24022650695886477918181338054308 / EXPPOLY_32F_A0,
 830         A2 = .55504108793649567998466049042729e-1 / EXPPOLY_32F_A0,
 831         A1 = .96180973140732918010002372686186e-2 / EXPPOLY_32F_A0,
 832         A0 = .13369713757180123244806654839424e-2 / EXPPOLY_32F_A0;
 833
 834     #undef EXPPOLY
 835     #define EXPPOLY(x)  (((((A0*(x) + A1)*(x) + A2)*(x) + A3)*(x) + A4)*(x) + A5)
 836
 837     int i = 0;
 838     DBLINT buf[4];
 839     const Cv64suf* x = (const Cv64suf*)_x;
 840
 841     if( !x || !y )
 842         return CV_NULLPTR_ERR;
 843     if( n <= 0 )
 844         return CV_BADSIZE_ERR;
 845
 846     buf[0].i.lo = buf[1].i.lo = buf[2].i.lo = buf[3].i.lo = 0;
 847
 848     for( ; i <= n - 4; i += 4 )
 849     {
 850         double x0 = x[i].f * exp_prescale;
 851         double x1 = x[i + 1].f * exp_prescale;
 852         double x2 = x[i + 2].f * exp_prescale;
 853         double x3 = x[i + 3].f * exp_prescale;
 854
 855         double y0, y1, y2, y3;
 856         int val0, val1, val2, val3, t;
 857
 858         t = (int)(x[i].i >> 52);
 859         if( (t & 2047) > 1023 + 10 )
 860             x0 = t < 0 ? -exp_max_val : exp_max_val;
 861
 862         t = (int)(x[i+1].i >> 52);
 863         if( (t & 2047) > 1023 + 10 )
 864             x1 = t < 0 ? -exp_max_val : exp_max_val;
 865
 866         t = (int)(x[i+2].i >> 52);
 867         if( (t & 2047) > 1023 + 10 )
 868             x2 = t < 0 ? -exp_max_val : exp_max_val;
 869
 870         t = (int)(x[i+3].i >> 52);
 871         if( (t & 2047) > 1023 + 10 )
 872             x3 = t < 0 ? -exp_max_val : exp_max_val;
 873
 874         val0 = cvRound(x0);
 875         val1 = cvRound(x1);
 876         val2 = cvRound(x2);
 877         val3 = cvRound(x3);
 878
 879         x0 = (x0 - val0)*exp_postscale;
 880         x1 = (x1 - val1)*exp_postscale;
 881         x2 = (x2 - val2)*exp_postscale;
 882         x3 = (x3 - val3)*exp_postscale;
 883
 884         t = (val0 >> EXPTAB_SCALE) + 1023;
 885         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 886         buf[0].i.hi = t << 20;
 887
 888         t = (val1 >> EXPTAB_SCALE) + 1023;
 889         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 890         buf[1].i.hi = t << 20;
 891
 892         t = (val2 >> EXPTAB_SCALE) + 1023;
 893         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 894         buf[2].i.hi = t << 20;
 895
 896         t = (val3 >> EXPTAB_SCALE) + 1023;
 897         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 898         buf[3].i.hi = t << 20;
 899
 900         y0 = buf[0].d * expTab[val0 & EXPTAB_MASK] * EXPPOLY( x0 );
 901         y1 = buf[1].d * expTab[val1 & EXPTAB_MASK] * EXPPOLY( x1 );
 902
 903         y[i] = y0;
 904         y[i + 1] = y1;
 905
 906         y2 = buf[2].d * expTab[val2 & EXPTAB_MASK] * EXPPOLY( x2 );
 907         y3 = buf[3].d * expTab[val3 & EXPTAB_MASK] * EXPPOLY( x3 );
 908
 909         y[i + 2] = y2;
 910         y[i + 3] = y3;
 911     }
 912
 913     for( ; i < n; i++ )
 914     {
 915         double x0 = x[i].f * exp_prescale;
 916         int val0, t;
 917
 918         t = (int)(x[i].i >> 52);
 919         if( (t & 2047) > 1023 + 10 )
 920             x0 = t < 0 ? -exp_max_val : exp_max_val;
 921
 922         val0 = cvRound(x0);
 923         t = (val0 >> EXPTAB_SCALE) + 1023;
 924         t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
 925
 926         buf[0].i.hi = t << 20;
 927         x0 = (x0 - val0)*exp_postscale;
 928
 929         y[i] = buf[0].d * expTab[val0 & EXPTAB_MASK] * EXPPOLY( x0 );
 930     }
 931
 932     return CV_OK;
 933 }
 934
 935 #undef EXPTAB_SCALE
 936 #undef EXPTAB_MASK
 937 #undef EXPPOLY_32F_A0
 938
 939 #else
 940
 941 #define Exp_32f ippsExp_32f_A21
 942 #define Exp_64f ippsExp_64f_A50
 943
 944 #endif
 945
 946 void exp( const Mat& src, Mat& dst )
 947 {
 948     int depth = src.depth();
 949     dst.create( src.size(), src.type() );
 950     Size size = getContinuousSize( src, dst, src.channels() );
 951
 952     MathFunc func = depth == CV_32F ? (MathFunc)Exp_32f : depth == CV_64F ? (MathFunc)Exp_64f : 0;
 953     CV_Assert(func != 0);
 954
 955     for( int y = 0; y < size.height; y++ )
 956         func( src.data + src.step*y, dst.data + dst.step*y, size.width );
 957 }
 958
 959
 960 /****************************************************************************************\
 961 *                                          L O G                                         *
 962 \****************************************************************************************/
 963
 964 #ifndef HAVE_IPP
 965
 966 #define LOGTAB_SCALE    8
 967 #define LOGTAB_MASK         ((1 << LOGTAB_SCALE) - 1)
 968 #define LOGTAB_MASK2        ((1 << (20 - LOGTAB_SCALE)) - 1)
 969 #define LOGTAB_MASK2_32F    ((1 << (23 - LOGTAB_SCALE)) - 1)
 970
 971 static const double icvLogTab[] = {
 972 0.0000000000000000000000000000000000000000,    1.000000000000000000000000000000000000000,
 973 .00389864041565732288852075271279318258166,    .9961089494163424124513618677042801556420,
 974 .00778214044205494809292034119607706088573,    .9922480620155038759689922480620155038760,
 975 .01165061721997527263705585198749759001657,    .9884169884169884169884169884169884169884,
 976 .01550418653596525274396267235488267033361,    .9846153846153846153846153846153846153846,
 977 .01934296284313093139406447562578250654042,    .9808429118773946360153256704980842911877,
 978 .02316705928153437593630670221500622574241,    .9770992366412213740458015267175572519084,
 979 .02697658769820207233514075539915211265906,    .9733840304182509505703422053231939163498,
 980 .03077165866675368732785500469617545604706,    .9696969696969696969696969696969696969697,
 981 .03455238150665972812758397481047722976656,    .9660377358490566037735849056603773584906,
 982 .03831886430213659461285757856785494368522,    .9624060150375939849624060150375939849624,
 983 .04207121392068705056921373852674150839447,    .9588014981273408239700374531835205992509,
 984 .04580953603129420126371940114040626212953,    .9552238805970149253731343283582089552239,
 985 .04953393512227662748292900118940451648088,    .9516728624535315985130111524163568773234,
 986 .05324451451881227759255210685296333394944,    .9481481481481481481481481481481481481481,
 987 .05694137640013842427411105973078520037234,    .9446494464944649446494464944649446494465,
 988 .06062462181643483993820353816772694699466,    .9411764705882352941176470588235294117647,
 989 .06429435070539725460836422143984236754475,    .9377289377289377289377289377289377289377,
 990 .06795066190850773679699159401934593915938,    .9343065693430656934306569343065693430657,
 991 .07159365318700880442825962290953611955044,    .9309090909090909090909090909090909090909,
 992 .07522342123758751775142172846244648098944,    .9275362318840579710144927536231884057971,
 993 .07884006170777602129362549021607264876369,    .9241877256317689530685920577617328519856,
 994 .08244366921107458556772229485432035289706,    .9208633093525179856115107913669064748201,
 995 .08603433734180314373940490213499288074675,    .9175627240143369175627240143369175627240,
 996 .08961215868968712416897659522874164395031,    .9142857142857142857142857142857142857143,
 997 .09317722485418328259854092721070628613231,    .9110320284697508896797153024911032028470,
 998 .09672962645855109897752299730200320482256,    .9078014184397163120567375886524822695035,
 999 .10026945316367513738597949668474029749630,    .9045936395759717314487632508833922261484,
1000 .10379679368164355934833764649738441221420,    .9014084507042253521126760563380281690141,
1001 .10731173578908805021914218968959175981580,    .8982456140350877192982456140350877192982,
1002 .11081436634029011301105782649756292812530,    .8951048951048951048951048951048951048951,
1003 .11430477128005862852422325204315711744130,    .8919860627177700348432055749128919860627,
1004 .11778303565638344185817487641543266363440,    .8888888888888888888888888888888888888889,
1005 .12124924363286967987640707633545389398930,    .8858131487889273356401384083044982698962,
1006 .12470347850095722663787967121606925502420,    .8827586206896551724137931034482758620690,
1007 .12814582269193003360996385708858724683530,    .8797250859106529209621993127147766323024,
1008 .13157635778871926146571524895989568904040,    .8767123287671232876712328767123287671233,
1009 .13499516453750481925766280255629681050780,    .8737201365187713310580204778156996587031,
1010 .13840232285911913123754857224412262439730,    .8707482993197278911564625850340136054422,
1011 .14179791186025733629172407290752744302150,    .8677966101694915254237288135593220338983,
1012 .14518200984449788903951628071808954700830,    .8648648648648648648648648648648648648649,
1013 .14855469432313711530824207329715136438610,    .8619528619528619528619528619528619528620,
1014 .15191604202584196858794030049466527998450,    .8590604026845637583892617449664429530201,
1015 .15526612891112392955683674244937719777230,    .8561872909698996655518394648829431438127,
1016 .15860503017663857283636730244325008243330,    .8533333333333333333333333333333333333333,
1017 .16193282026931324346641360989451641216880,    .8504983388704318936877076411960132890365,
1018 .16524957289530714521497145597095368430010,    .8476821192052980132450331125827814569536,
1019 .16855536102980664403538924034364754334090,    .8448844884488448844884488448844884488449,
1020 .17185025692665920060697715143760433420540,    .8421052631578947368421052631578947368421,
1021 .17513433212784912385018287750426679849630,    .8393442622950819672131147540983606557377,
1022 .17840765747281828179637841458315961062910,    .8366013071895424836601307189542483660131,
1023 .18167030310763465639212199675966985523700,    .8338762214983713355048859934853420195440,
1024 .18492233849401198964024217730184318497780,    .8311688311688311688311688311688311688312,
1025 .18816383241818296356839823602058459073300,    .8284789644012944983818770226537216828479,
1026 .19139485299962943898322009772527962923050,    .8258064516129032258064516129032258064516,
1027 .19461546769967164038916962454095482826240,    .8231511254019292604501607717041800643087,
1028 .19782574332991986754137769821682013571260,    .8205128205128205128205128205128205128205,
1029 .20102574606059073203390141770796617493040,    .8178913738019169329073482428115015974441,
1030 .20421554142869088876999228432396193966280,    .8152866242038216560509554140127388535032,
1031 .20739519434607056602715147164417430758480,    .8126984126984126984126984126984126984127,
1032 .21056476910734961416338251183333341032260,    .8101265822784810126582278481012658227848,
1033 .21372432939771812687723695489694364368910,    .8075709779179810725552050473186119873817,
1034 .21687393830061435506806333251006435602900,    .8050314465408805031446540880503144654088,
1035 .22001365830528207823135744547471404075630,    .8025078369905956112852664576802507836991,
1036 .22314355131420973710199007200571941211830,    .8000000000000000000000000000000000000000,
1037 .22626367865045338145790765338460914790630,    .7975077881619937694704049844236760124611,
1038 .22937410106484582006380890106811420992010,    .7950310559006211180124223602484472049689,
1039 .23247487874309405442296849741978803649550,    .7925696594427244582043343653250773993808,
1040 .23556607131276688371634975283086532726890,    .7901234567901234567901234567901234567901,
1041 .23864773785017498464178231643018079921600,    .7876923076923076923076923076923076923077,
1042 .24171993688714515924331749374687206000090,    .7852760736196319018404907975460122699387,
1043 .24478272641769091566565919038112042471760,    .7828746177370030581039755351681957186544,
1044 .24783616390458124145723672882013488560910,    .7804878048780487804878048780487804878049,
1045 .25088030628580937353433455427875742316250,    .7781155015197568389057750759878419452888,
1046 .25391520998096339667426946107298135757450,    .7757575757575757575757575757575757575758,
1047 .25694093089750041913887912414793390780680,    .7734138972809667673716012084592145015106,
1048 .25995752443692604627401010475296061486000,    .7710843373493975903614457831325301204819,
1049 .26296504550088134477547896494797896593800,    .7687687687687687687687687687687687687688,
1050 .26596354849713793599974565040611196309330,    .7664670658682634730538922155688622754491,
1051 .26895308734550393836570947314612567424780,    .7641791044776119402985074626865671641791,
1052 .27193371548364175804834985683555714786050,    .7619047619047619047619047619047619047619,
1053 .27490548587279922676529508862586226314300,    .7596439169139465875370919881305637982196,
1054 .27786845100345625159121709657483734190480,    .7573964497041420118343195266272189349112,
1055 .28082266290088775395616949026589281857030,    .7551622418879056047197640117994100294985,
1056 .28376817313064456316240580235898960381750,    .7529411764705882352941176470588235294118,
1057 .28670503280395426282112225635501090437180,    .7507331378299120234604105571847507331378,
1058 .28963329258304265634293983566749375313530,    .7485380116959064327485380116959064327485,
1059 .29255300268637740579436012922087684273730,    .7463556851311953352769679300291545189504,
1060 .29546421289383584252163927885703742504130,    .7441860465116279069767441860465116279070,
1061 .29836697255179722709783618483925238251680,    .7420289855072463768115942028985507246377,
1062 .30126133057816173455023545102449133992200,    .7398843930635838150289017341040462427746,
1063 .30414733546729666446850615102448500692850,    .7377521613832853025936599423631123919308,
1064 .30702503529491181888388950937951449304830,    .7356321839080459770114942528735632183908,
1065 .30989447772286465854207904158101882785550,    .7335243553008595988538681948424068767908,
1066 .31275571000389684739317885942000430077330,    .7314285714285714285714285714285714285714,
1067 .31560877898630329552176476681779604405180,    .7293447293447293447293447293447293447293,
1068 .31845373111853458869546784626436419785030,    .7272727272727272727272727272727272727273,
1069 .32129061245373424782201254856772720813750,    .7252124645892351274787535410764872521246,
1070 .32411946865421192853773391107097268104550,    .7231638418079096045197740112994350282486,
1071 .32694034499585328257253991068864706903700,    .7211267605633802816901408450704225352113,
1072 .32975328637246797969240219572384376078850,    .7191011235955056179775280898876404494382,
1073 .33255833730007655635318997155991382896900,    .7170868347338935574229691876750700280112,
1074 .33535554192113781191153520921943709254280,    .7150837988826815642458100558659217877095,
1075 .33814494400871636381467055798566434532400,    .7130919220055710306406685236768802228412,
1076 .34092658697059319283795275623560883104800,    .7111111111111111111111111111111111111111,
1077 .34370051385331840121395430287520866841080,    .7091412742382271468144044321329639889197,
1078 .34646676734620857063262633346312213689100,    .7071823204419889502762430939226519337017,
1079 .34922538978528827602332285096053965389730,    .7052341597796143250688705234159779614325,
1080 .35197642315717814209818925519357435405250,    .7032967032967032967032967032967032967033,
1081 .35471990910292899856770532096561510115850,    .7013698630136986301369863013698630136986,
1082 .35745588892180374385176833129662554711100,    .6994535519125683060109289617486338797814,
1083 .36018440357500774995358483465679455548530,    .6975476839237057220708446866485013623978,
1084 .36290549368936841911903457003063522279280,    .6956521739130434782608695652173913043478,
1085 .36561919956096466943762379742111079394830,    .6937669376693766937669376693766937669377,
1086 .36832556115870762614150635272380895912650,    .6918918918918918918918918918918918918919,
1087 .37102461812787262962487488948681857436900,    .6900269541778975741239892183288409703504,
1088 .37371640979358405898480555151763837784530,    .6881720430107526881720430107526881720430,
1089 .37640097516425302659470730759494472295050,    .6863270777479892761394101876675603217158,
1090 .37907835293496944251145919224654790014030,    .6844919786096256684491978609625668449198,
1091 .38174858149084833769393299007788300514230,    .6826666666666666666666666666666666666667,
1092 .38441169891033200034513583887019194662580,    .6808510638297872340425531914893617021277,
1093 .38706774296844825844488013899535872042180,    .6790450928381962864721485411140583554377,
1094 .38971675114002518602873692543653305619950,    .6772486772486772486772486772486772486772,
1095 .39235876060286384303665840889152605086580,    .6754617414248021108179419525065963060686,
1096 .39499380824086893770896722344332374632350,    .6736842105263157894736842105263157894737,
1097 .39762193064713846624158577469643205404280,    .6719160104986876640419947506561679790026,
1098 .40024316412701266276741307592601515352730,    .6701570680628272251308900523560209424084,
1099 .40285754470108348090917615991202183067800,    .6684073107049608355091383812010443864230,
1100 .40546510810816432934799991016916465014230,    .6666666666666666666666666666666666666667,
1101 .40806588980822172674223224930756259709600,    .6649350649350649350649350649350649350649,
1102 .41065992498526837639616360320360399782650,    .6632124352331606217616580310880829015544,
1103 .41324724855021932601317757871584035456180,    .6614987080103359173126614987080103359173,
1104 .41582789514371093497757669865677598863850,    .6597938144329896907216494845360824742268,
1105 .41840189913888381489925905043492093682300,    .6580976863753213367609254498714652956298,
1106 .42096929464412963239894338585145305842150,    .6564102564102564102564102564102564102564,
1107 .42353011550580327293502591601281892508280,    .6547314578005115089514066496163682864450,
1108 .42608439531090003260516141381231136620050,    .6530612244897959183673469387755102040816,
1109 .42863216738969872610098832410585600882780,    .6513994910941475826972010178117048346056,
1110 .43117346481837132143866142541810404509300,    .6497461928934010152284263959390862944162,
1111 .43370832042155937902094819946796633303180,    .6481012658227848101265822784810126582278,
1112 .43623676677491801667585491486534010618930,    .6464646464646464646464646464646464646465,
1113 .43875883620762790027214350629947148263450,    .6448362720403022670025188916876574307305,
1114 .44127456080487520440058801796112675219780,    .6432160804020100502512562814070351758794,
1115 .44378397241030093089975139264424797147500,    .6416040100250626566416040100250626566416,
1116 .44628710262841947420398014401143882423650,    .6400000000000000000000000000000000000000,
1117 .44878398282700665555822183705458883196130,    .6384039900249376558603491271820448877805,
1118 .45127464413945855836729492693848442286250,    .6368159203980099502487562189054726368159,
1119 .45375911746712049854579618113348260521900,    .6352357320099255583126550868486352357320,
1120 .45623743348158757315857769754074979573500,    .6336633663366336633663366336633663366337,
1121 .45870962262697662081833982483658473938700,    .6320987654320987654320987654320987654321,
1122 .46117571512217014895185229761409573256980,    .6305418719211822660098522167487684729064,
1123 .46363574096303250549055974261136725544930,    .6289926289926289926289926289926289926290,
1124 .46608972992459918316399125615134835243230,    .6274509803921568627450980392156862745098,
1125 .46853771156323925639597405279346276074650,    .6259168704156479217603911980440097799511,
1126 .47097971521879100631480241645476780831830,    .6243902439024390243902439024390243902439,
1127 .47341577001667212165614273544633761048330,    .6228710462287104622871046228710462287105,
1128 .47584590486996386493601107758877333253630,    .6213592233009708737864077669902912621359,
1129 .47827014848147025860569669930555392056700,    .6198547215496368038740920096852300242131,
1130 .48068852934575190261057286988943815231330,    .6183574879227053140096618357487922705314,
1131 .48310107575113581113157579238759353756900,    .6168674698795180722891566265060240963855,
1132 .48550781578170076890899053978500887751580,    .6153846153846153846153846153846153846154,
1133 .48790877731923892879351001283794175833480,    .6139088729016786570743405275779376498801,
1134 .49030398804519381705802061333088204264650,    .6124401913875598086124401913875598086124,
1135 .49269347544257524607047571407747454941280,    .6109785202863961813842482100238663484487,
1136 .49507726679785146739476431321236304938800,    .6095238095238095238095238095238095238095,
1137 .49745538920281889838648226032091770321130,    .6080760095011876484560570071258907363420,
1138 .49982786955644931126130359189119189977650,    .6066350710900473933649289099526066350711,
1139 .50219473456671548383667413872899487614650,    .6052009456264775413711583924349881796690,
1140 .50455601075239520092452494282042607665050,    .6037735849056603773584905660377358490566,
1141 .50691172444485432801997148999362252652650,    .6023529411764705882352941176470588235294,
1142 .50926190178980790257412536448100581765150,    .6009389671361502347417840375586854460094,
1143 .51160656874906207391973111953120678663250,    .5995316159250585480093676814988290398126,
1144 .51394575110223428282552049495279788970950,    .5981308411214953271028037383177570093458,
1145 .51627947444845445623684554448118433356300,    .5967365967365967365967365967365967365967,
1146 .51860776420804555186805373523384332656850,    .5953488372093023255813953488372093023256,
1147 .52093064562418522900344441950437612831600,    .5939675174013921113689095127610208816705,
1148 .52324814376454775732838697877014055848100,    .5925925925925925925925925925925925925926,
1149 .52556028352292727401362526507000438869000,    .5912240184757505773672055427251732101617,
1150 .52786708962084227803046587723656557500350,    .5898617511520737327188940092165898617512,
1151 .53016858660912158374145519701414741575700,    .5885057471264367816091954022988505747126,
1152 .53246479886947173376654518506256863474850,    .5871559633027522935779816513761467889908,
1153 .53475575061602764748158733709715306758900,    .5858123569794050343249427917620137299771,
1154 .53704146589688361856929077475797384977350,    .5844748858447488584474885844748858447489,
1155 .53932196859560876944783558428753167390800,    .5831435079726651480637813211845102505695,
1156 .54159728243274429804188230264117009937750,    .5818181818181818181818181818181818181818,
1157 .54386743096728351609669971367111429572100,    .5804988662131519274376417233560090702948,
1158 .54613243759813556721383065450936555862450,    .5791855203619909502262443438914027149321,
1159 .54839232556557315767520321969641372561450,    .5778781038374717832957110609480812641084,
1160 .55064711795266219063194057525834068655950,    .5765765765765765765765765765765765765766,
1161 .55289683768667763352766542084282264113450,    .5752808988764044943820224719101123595506,
1162 .55514150754050151093110798683483153581600,    .5739910313901345291479820627802690582960,
1163 .55738115013400635344709144192165695130850,    .5727069351230425055928411633109619686801,
1164 .55961578793542265941596269840374588966350,    .5714285714285714285714285714285714285714,
1165 .56184544326269181269140062795486301183700,    .5701559020044543429844097995545657015590,
1166 .56407013828480290218436721261241473257550,    .5688888888888888888888888888888888888889,
1167 .56628989502311577464155334382667206227800,    .5676274944567627494456762749445676274945,
1168 .56850473535266865532378233183408156037350,    .5663716814159292035398230088495575221239,
1169 .57071468100347144680739575051120482385150,    .5651214128035320088300220750551876379691,
1170 .57291975356178548306473885531886480748650,    .5638766519823788546255506607929515418502,
1171 .57511997447138785144460371157038025558000,    .5626373626373626373626373626373626373626,
1172 .57731536503482350219940144597785547375700,    .5614035087719298245614035087719298245614,
1173 .57950594641464214795689713355386629700650,    .5601750547045951859956236323851203501094,
1174 .58169173963462239562716149521293118596100,    .5589519650655021834061135371179039301310,
1175 .58387276558098266665552955601015128195300,    .5577342047930283224400871459694989106754,
1176 .58604904500357812846544902640744112432000,    .5565217391304347826086956521739130434783,
1177 .58822059851708596855957011939608491957200,    .5553145336225596529284164859002169197397,
1178 .59038744660217634674381770309992134571100,    .5541125541125541125541125541125541125541,
1179 .59254960960667157898740242671919986605650,    .5529157667386609071274298056155507559395,
1180 .59470710774669277576265358220553025603300,    .5517241379310344827586206896551724137931,
1181 .59685996110779382384237123915227130055450,    .5505376344086021505376344086021505376344,
1182 .59900818964608337768851242799428291618800,    .5493562231759656652360515021459227467811,
1183 .60115181318933474940990890900138765573500,    .5481798715203426124197002141327623126338,
1184 .60329085143808425240052883964381180703650,    .5470085470085470085470085470085470085470,
1185 .60542532396671688843525771517306566238400,    .5458422174840085287846481876332622601279,
1186 .60755525022454170969155029524699784815300,    .5446808510638297872340425531914893617021,
1187 .60968064953685519036241657886421307921400,    .5435244161358811040339702760084925690021,
1188 .61180154110599282990534675263916142284850,    .5423728813559322033898305084745762711864,
1189 .61391794401237043121710712512140162289150,    .5412262156448202959830866807610993657505,
1190 .61602987721551394351138242200249806046500,    .5400843881856540084388185654008438818565,
1191 .61813735955507864705538167982012964785100,    .5389473684210526315789473684210526315789,
1192 .62024040975185745772080281312810257077200,    .5378151260504201680672268907563025210084,
1193 .62233904640877868441606324267922900617100,    .5366876310272536687631027253668763102725,
1194 .62443328801189346144440150965237990021700,    .5355648535564853556485355648535564853556,
1195 .62652315293135274476554741340805776417250,    .5344467640918580375782881002087682672234,
1196 .62860865942237409420556559780379757285100,    .5333333333333333333333333333333333333333,
1197 .63068982562619868570408243613201193511500,    .5322245322245322245322245322245322245322,
1198 .63276666957103777644277897707070223987100,    .5311203319502074688796680497925311203320,
1199 .63483920917301017716738442686619237065300,    .5300207039337474120082815734989648033126,
1200 .63690746223706917739093569252872839570050,    .5289256198347107438016528925619834710744,
1201 .63897144645792069983514238629140891134750,    .5278350515463917525773195876288659793814,
1202 .64103117942093124081992527862894348800200,    .5267489711934156378600823045267489711934,
1203 .64308667860302726193566513757104985415950,    .5256673511293634496919917864476386036961,
1204 .64513796137358470073053240412264131009600,    .5245901639344262295081967213114754098361,
1205 .64718504499530948859131740391603671014300,    .5235173824130879345603271983640081799591,
1206 .64922794662510974195157587018911726772800,    .5224489795918367346938775510204081632653,
1207 .65126668331495807251485530287027359008800,    .5213849287169042769857433808553971486762,
1208 .65330127201274557080523663898929953575150,    .5203252032520325203252032520325203252033,
1209 .65533172956312757406749369692988693714150,    .5192697768762677484787018255578093306288,
1210 .65735807270835999727154330685152672231200,    .5182186234817813765182186234817813765182,
1211 .65938031808912778153342060249997302889800,    .5171717171717171717171717171717171717172,
1212 .66139848224536490484126716182800009846700,    .5161290322580645161290322580645161290323,
1213 .66341258161706617713093692145776003599150,    .5150905432595573440643863179074446680080,
1214 .66542263254509037562201001492212526500250,    .5140562248995983935742971887550200803213,
1215 .66742865127195616370414654738851822912700,    .5130260521042084168336673346693386773547,
1216 .66943065394262923906154583164607174694550,    .5120000000000000000000000000000000000000,
1217 .67142865660530226534774556057527661323550,    .5109780439121756487025948103792415169661,
1218 .67342267521216669923234121597488410770900,    .5099601593625498007968127490039840637450,
1219 .67541272562017662384192817626171745359900,    .5089463220675944333996023856858846918489,
1220 .67739882359180603188519853574689477682100,    .5079365079365079365079365079365079365079,
1221 .67938098479579733801614338517538271844400,    .5069306930693069306930693069306930693069,
1222 .68135922480790300781450241629499942064300,    .5059288537549407114624505928853754940711,
1223 .68333355911162063645036823800182901322850,    .5049309664694280078895463510848126232742,
1224 .68530400309891936760919861626462079584600,    .5039370078740157480314960629921259842520,
1225 .68727057207096020619019327568821609020250,    .5029469548133595284872298624754420432220,
1226 .68923328123880889251040571252815425395950,    .5019607843137254901960784313725490196078,
1227 .69314718055994530941723212145818, 5.0e-01,
1228 };
1229
1230
1231
1232 #define LOGTAB_TRANSLATE(x,h) (((x) - 1.)*icvLogTab[(h)+1])
1233 static const double ln_2 = 0.69314718055994530941723212145818;
1234
1235 static CvStatus CV_STDCALL Log_32f( const float *_x, float *y, int n )
1236 {
1237     static const double shift[] = { 0, -1./512 };
1238     static const double
1239         A0 = 0.3333333333333333333333333,
1240         A1 = -0.5,
1241         A2 = 1;
1242
1243     #undef LOGPOLY
1244     #define LOGPOLY(x,k) ((x)+=shift[k],((A0*(x) + A1)*(x) + A2)*(x))
1245
1246     int i = 0;
1247     union
1248     {
1249         int i;
1250         float f;
1251     }
1252     buf[4];
1253
1254     const int* x = (const int*)_x;
1255
1256     if( !x || !y )
1257         return CV_NULLPTR_ERR;
1258     if( n <= 0 )
1259         return CV_BADSIZE_ERR;
1260
1261     for( i = 0; i <= n - 4; i += 4 )
1262     {
1263         double x0, x1, x2, x3;
1264         double y0, y1, y2, y3;
1265         int h0, h1, h2, h3;
1266
1267         h0 = x[i];
1268         h1 = x[i+1];
1269         buf[0].i = (h0 & LOGTAB_MASK2_32F) | (127 << 23);
1270         buf[1].i = (h1 & LOGTAB_MASK2_32F) | (127 << 23);
1271
1272         y0 = (((h0 >> 23) & 0xff) - 127) * ln_2;
1273         y1 = (((h1 >> 23) & 0xff) - 127) * ln_2;
1274
1275         h0 = (h0 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1276         h1 = (h1 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1277
1278         y0 += icvLogTab[h0];
1279         y1 += icvLogTab[h1];
1280
1281         h2 = x[i+2];
1282         h3 = x[i+3];
1283
1284         x0 = LOGTAB_TRANSLATE( buf[0].f, h0 );
1285         x1 = LOGTAB_TRANSLATE( buf[1].f, h1 );
1286
1287         buf[2].i = (h2 & LOGTAB_MASK2_32F) | (127 << 23);
1288         buf[3].i = (h3 & LOGTAB_MASK2_32F) | (127 << 23);
1289
1290         y2 = (((h2 >> 23) & 0xff) - 127) * ln_2;
1291         y3 = (((h3 >> 23) & 0xff) - 127) * ln_2;
1292
1293         h2 = (h2 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1294         h3 = (h3 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1295
1296         y2 += icvLogTab[h2];
1297         y3 += icvLogTab[h3];
1298
1299         x2 = LOGTAB_TRANSLATE( buf[2].f, h2 );
1300         x3 = LOGTAB_TRANSLATE( buf[3].f, h3 );
1301
1302         y0 += LOGPOLY( x0, h0 == 510 );
1303         y1 += LOGPOLY( x1, h1 == 510 );
1304
1305         y[i] = (float) y0;
1306         y[i + 1] = (float) y1;
1307
1308         y2 += LOGPOLY( x2, h2 == 510 );
1309         y3 += LOGPOLY( x3, h3 == 510 );
1310
1311         y[i + 2] = (float) y2;
1312         y[i + 3] = (float) y3;
1313     }
1314
1315     for( ; i < n; i++ )
1316     {
1317         int h0 = x[i];
1318         double x0, y0;
1319
1320         y0 = (((h0 >> 23) & 0xff) - 127) * ln_2;
1321
1322         buf[0].i = (h0 & LOGTAB_MASK2_32F) | (127 << 23);
1323         h0 = (h0 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1324
1325         y0 += icvLogTab[h0];
1326         x0 = LOGTAB_TRANSLATE( buf[0].f, h0 );
1327         y0 += LOGPOLY( x0, h0 == 510 );
1328
1329         y[i] = (float)y0;
1330     }
1331
1332     return CV_OK;
1333 }
1334
1335
1336 static CvStatus CV_STDCALL Log_64f( const double *x, double *y, int n )
1337 {
1338     static const double shift[] = { 0, -1./512 };
1339     static const double
1340         A7 = 1.0,
1341         A6 = -0.5,
1342         A5 = 0.333333333333333314829616256247390992939472198486328125,
1343         A4 = -0.25,
1344         A3 = 0.2,
1345         A2 = -0.1666666666666666574148081281236954964697360992431640625,
1346         A1 = 0.1428571428571428769682682968777953647077083587646484375,
1347         A0 = -0.125;
1348
1349     #undef LOGPOLY
1350     #define LOGPOLY(x,k) ((x)+=shift[k], xq = (x)*(x),\
1351         (((A0*xq + A2)*xq + A4)*xq + A6)*xq + \
1352         (((A1*xq + A3)*xq + A5)*xq + A7)*(x))
1353
1354     int i = 0;
1355     DBLINT buf[4];
1356     DBLINT *X = (DBLINT *) x;
1357
1358     if( !x || !y )
1359         return CV_NULLPTR_ERR;
1360     if( n <= 0 )
1361         return CV_BADSIZE_ERR;
1362
1363     for( ; i <= n - 4; i += 4 )
1364     {
1365         double xq;
1366         double x0, x1, x2, x3;
1367         double y0, y1, y2, y3;
1368         int h0, h1, h2, h3;
1369
1370         h0 = X[i].i.lo;
1371         h1 = X[i + 1].i.lo;
1372         buf[0].i.lo = h0;
1373         buf[1].i.lo = h1;
1374
1375         h0 = X[i].i.hi;
1376         h1 = X[i + 1].i.hi;
1377         buf[0].i.hi = (h0 & LOGTAB_MASK2) | (1023 << 20);
1378         buf[1].i.hi = (h1 & LOGTAB_MASK2) | (1023 << 20);
1379
1380         y0 = (((h0 >> 20) & 0x7ff) - 1023) * ln_2;
1381         y1 = (((h1 >> 20) & 0x7ff) - 1023) * ln_2;
1382
1383         h2 = X[i + 2].i.lo;
1384         h3 = X[i + 3].i.lo;
1385         buf[2].i.lo = h2;
1386         buf[3].i.lo = h3;
1387
1388         h0 = (h0 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1389         h1 = (h1 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1390
1391         y0 += icvLogTab[h0];
1392         y1 += icvLogTab[h1];
1393
1394         h2 = X[i + 2].i.hi;
1395         h3 = X[i + 3].i.hi;
1396
1397         x0 = LOGTAB_TRANSLATE( buf[0].d, h0 );
1398         x1 = LOGTAB_TRANSLATE( buf[1].d, h1 );
1399
1400         buf[2].i.hi = (h2 & LOGTAB_MASK2) | (1023 << 20);
1401         buf[3].i.hi = (h3 & LOGTAB_MASK2) | (1023 << 20);
1402
1403         y2 = (((h2 >> 20) & 0x7ff) - 1023) * ln_2;
1404         y3 = (((h3 >> 20) & 0x7ff) - 1023) * ln_2;
1405
1406         h2 = (h2 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1407         h3 = (h3 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1408
1409         y2 += icvLogTab[h2];
1410         y3 += icvLogTab[h3];
1411
1412         x2 = LOGTAB_TRANSLATE( buf[2].d, h2 );
1413         x3 = LOGTAB_TRANSLATE( buf[3].d, h3 );
1414
1415         y0 += LOGPOLY( x0, h0 == 510 );
1416         y1 += LOGPOLY( x1, h1 == 510 );
1417
1418         y[i] = y0;
1419         y[i + 1] = y1;
1420
1421         y2 += LOGPOLY( x2, h2 == 510 );
1422         y3 += LOGPOLY( x3, h3 == 510 );
1423
1424         y[i + 2] = y2;
1425         y[i + 3] = y3;
1426     }
1427
1428     for( ; i < n; i++ )
1429     {
1430         int h0 = X[i].i.hi;
1431         double xq;
1432         double x0, y0 = (((h0 >> 20) & 0x7ff) - 1023) * ln_2;
1433
1434         buf[0].i.hi = (h0 & LOGTAB_MASK2) | (1023 << 20);
1435         buf[0].i.lo = X[i].i.lo;
1436         h0 = (h0 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1437
1438         y0 += icvLogTab[h0];
1439         x0 = LOGTAB_TRANSLATE( buf[0].d, h0 );
1440         y0 += LOGPOLY( x0, h0 == 510 );
1441         y[i] = y0;
1442     }
1443
1444     return CV_OK;
1445 }
1446
1447 #else
1448
1449 #define Log_32f ippsLn_32f_A21
1450 #define Log_64f ippsLn_64f_A50
1451
1452 #endif
1453
1454 void log( const Mat& src, Mat& dst )
1455 {
1456     int depth = src.depth();
1457     dst.create( src.size(), src.type() );
1458     Size size = getContinuousSize( src, dst, src.channels() );
1459
1460     MathFunc func = depth == CV_32F ? (MathFunc)Log_32f : depth == CV_64F ? (MathFunc)Log_64f : 0;
1461     CV_Assert(func != 0);
1462
1463     for( int y = 0; y < size.height; y++ )
1464         func( src.data + src.step*y, dst.data + dst.step*y, size.width );
1465 }
1466
1467
1468 /****************************************************************************************\
1469 *                                    P O W E R                                           *
1470 \****************************************************************************************/
1471
1472 template<typename T, typename WT>
1473 static CvStatus CV_STDCALL
1474 IPow( const void* _src, void* _dst, int len, int power )
1475 {
1476     int i;
1477     const T* src = (const T*)_src;
1478     T* dst = (T*)_dst;
1479     for( i = 0; i < len; i++ )
1480     {
1481         WT a = 1, b = src[i];
1482         int p = power;
1483         while( p > 1 )
1484         {
1485             if( p & 1 )
1486                 a *= b;
1487             b *= b;
1488             p >>= 1;
1489         }
1490
1491         a *= b;
1492         dst[i] = saturate_cast<T>(a);
1493     }
1494     return CV_OK;
1495 }
1496
1497 typedef CvStatus (CV_STDCALL * IPowFunc)( const void* src, void* dst, int len, int power );
1498
1499 void pow( const Mat& _src, double power, Mat& dst )
1500 {
1501     int ipower = cvRound( power ), i, j;
1502     bool is_ipower = 0;
1503     int depth = _src.depth();
1504     const Mat* src = &_src;
1505
1506     dst.create( _src.size(), _src.type() );
1507
1508     if( fabs(ipower - power) < DBL_EPSILON )
1509     {
1510         if( ipower < 0 )
1511         {
1512             divide( 1., _src, dst );
1513             if( ipower == -1 )
1514                 return;
1515             ipower = -ipower;
1516             src = &dst;
1517         }
1518
1519         switch( ipower )
1520         {
1521         case 0:
1522             dst = Scalar::all(1);
1523             return;
1524         case 1:
1525             src->copyTo(dst);
1526             return;
1527         case 2:
1528             multiply(*src, *src, dst);
1529             return;
1530         default:
1531             is_ipower = true;
1532         }
1533     }
1534     else
1535         CV_Assert( depth == CV_32F || depth == CV_64F );
1536
1537     Size size = getContinuousSize( *src, dst, src->channels() );
1538
1539     if( is_ipower )
1540     {
1541         static IPowFunc tab[] =
1542         {
1543             IPow<uchar, int>, 0, IPow<ushort, int>, IPow<short, int>, IPow<int, int>,
1544             IPow<float, float>, IPow<double, double>, 0
1545         };
1546
1547         IPowFunc func = tab[depth];
1548         CV_Assert( func != 0 );
1549         for( i = 0; i < size.height; i++ )
1550             func( src->data + src->step*i, dst.data + dst.step*i, size.width, ipower );
1551     }
1552     else if( fabs(fabs(power) - 0.5) < DBL_EPSILON )
1553     {
1554         MathFunc func = power < 0 ?
1555             (depth == CV_32F ? (MathFunc)InvSqrt_32f : (MathFunc)InvSqrt_64f) :
1556             (depth == CV_32F ? (MathFunc)Sqrt_32f : (MathFunc)Sqrt_64f);
1557
1558         for( i = 0; i < size.height; i++ )
1559             func( src->data + src->step*i, dst.data + dst.step*i, size.width );
1560     }
1561     else if( depth == CV_32F )
1562     {
1563         const float *x = (const float*)src->data;
1564         float *y = (float*)dst.data;
1565         size_t xstep = src->step/sizeof(x[0]), ystep = dst.step/sizeof(y[0]);
1566         float p = (float)power;
1567
1568         for( ; size.height--; x += xstep, y += ystep )
1569         {
1570             for( i = 0; i < size.width; i += MAX_BLOCK_SIZE )
1571             {
1572                 int block_size = std::min(MAX_BLOCK_SIZE, size.width - i);
1573                 Log_32f(x + i, y + i, block_size);
1574                 for( j = 0; j < block_size; j++ )
1575                     y[i + j] *= p;
1576                 Exp_32f(y + i, y + i, block_size);
1577             }
1578         }
1579     }
1580     else
1581     {
1582         const double *x = (const double*)src->data;
1583         double *y = (double*)dst.data;
1584         size_t xstep = src->step/sizeof(x[0]), ystep = dst.step/sizeof(y[0]);
1585
1586         for( ; size.height--; x += xstep, y += ystep )
1587         {
1588             for( i = 0; i < size.width; i += MAX_BLOCK_SIZE )
1589             {
1590                 int block_size = std::min(MAX_BLOCK_SIZE, size.width - i);
1591                 Log_64f(x + i, y + i, block_size);
1592                 for( j = 0; j < block_size; j++ )
1593                     y[i + j] *= power;
1594                 Exp_64f(y + i, y + i, block_size);
1595             }
1596         }
1597     }
1598 }
1599
1600 void sqrt(const Mat& a, Mat& b)
1601 {
1602     pow(a, 0.5, b);
1603 }
1604
1605 /************************** CheckArray for NaN's, Inf's *********************************/
1606
1607 bool checkRange(const Mat& src, bool quiet, Point* pt,
1608                 double minVal, double maxVal)
1609 {
1610     int depth = src.depth();
1611     Point badPt(-1, -1);
1612     double badValue = 0;
1613
1614     if( depth < CV_32F )
1615     {
1616         double m = 0, M = 0, badValue = 0;
1617         Point mp, MP, badPt(-1,-1);
1618         minMaxLoc(src.reshape(1,0), &m, &M, &mp, &MP);
1619         if( M >= maxVal )
1620         {
1621             badPt = MP;
1622             badValue = M;
1623         }
1624         else if( m < minVal )
1625         {
1626             badPt = mp;
1627             badValue = m;
1628         }
1629     }
1630     else
1631     {
1632         int i, loc = 0;
1633         Size size = getContinuousSize( src, src.channels() );
1634
1635         if( depth == CV_32F )
1636         {
1637             Cv32suf a, b;
1638             int ia, ib;
1639             const int* isrc = (const int*)src.data;
1640             size_t step = src.step/sizeof(isrc[0]);
1641
1642             a.f = (float)std::max(minVal, (double)-FLT_MAX);
1643             b.f = (float)std::min(maxVal, (double)FLT_MAX);
1644
1645             ia = CV_TOGGLE_FLT(a.i);
1646             ib = CV_TOGGLE_FLT(b.i);
1647
1648             for( ; badPt.x < 0 && size.height--; loc += size.width, isrc += step )
1649             {
1650                 for( i = 0; i < size.width; i++ )
1651                 {
1652                     int val = isrc[i];
1653                     val = CV_TOGGLE_FLT(val);
1654
1655                     if( val < ia || val >= ib )
1656                     {
1657                         badPt = Point((loc + i) % src.cols, (loc + i) / src.cols);
1658                         badValue = ((const float*)isrc)[i];
1659                         break;
1660                     }
1661                 }
1662             }
1663         }
1664         else
1665         {
1666             Cv64suf a, b;
1667             int64 ia, ib;
1668             const int64* isrc = (const int64*)src.data;
1669             size_t step = src.step/sizeof(isrc[0]);
1670
1671             a.f = minVal;
1672             b.f = maxVal;
1673
1674             ia = CV_TOGGLE_DBL(a.i);
1675             ib = CV_TOGGLE_DBL(b.i);
1676
1677             for( ; badPt.x < 0 && size.height--; loc += size.width, isrc += step )
1678             {
1679                 for( i = 0; i < size.width; i++ )
1680                 {
1681                     int64 val = isrc[i];
1682                     val = CV_TOGGLE_DBL(val);
1683
1684                     if( val < ia || val >= ib )
1685                     {
1686                         badPt = Point((loc + i) % src.cols, (loc + i) / src.cols);
1687                         badValue = ((const double*)isrc)[i];
1688                         break;
1689                     }
1690                 }
1691             }
1692         }
1693     }
1694
1695     if( badPt.x >= 0 )
1696     {
1697         if( pt )
1698             *pt = badPt;
1699         if( !quiet )
1700             CV_Error_( CV_StsOutOfRange,
1701             ("the value at (%d, %d)=%g is out of range", badPt.x, badPt.y, badValue));
1702     }
1703     return badPt.x < 0;
1704 }
1705
1706 }
1707
1708 CV_IMPL float cvCbrt(float value) { return cv::cubeRoot(value); }
1709 CV_IMPL float cvFastArctan(float y, float x) { return cv::fastAtan2(y, x); }
1710
1711 CV_IMPL void
1712 cvCartToPolar( const CvArr* xarr, const CvArr* yarr,
1713                CvArr* magarr, CvArr* anglearr,
1714                int angle_in_degrees )
1715 {
1716     cv::Mat X = cv::cvarrToMat(xarr), Y = cv::cvarrToMat(yarr), Mag, Angle;
1717     if( magarr )
1718     {
1719         Mag = cv::cvarrToMat(magarr);
1720         CV_Assert( Mag.size() == X.size() && Mag.type() == X.type() );
1721     }
1722     if( anglearr )
1723     {
1724         Angle = cv::cvarrToMat(anglearr);
1725         CV_Assert( Angle.size() == X.size() && Angle.type() == X.type() );
1726     }
1727         if( magarr )
1728         {
1729                 if( anglearr )
1730                         cv::cartToPolar( X, Y, Mag, Angle, angle_in_degrees != 0 );
1731                 else
1732                         cv::magnitude( X, Y, Mag );
1733         }
1734         else
1735                 cv::phase( X, Y, Angle, angle_in_degrees != 0 );
1736 }
1737
1738 CV_IMPL void
1739 cvPolarToCart( const CvArr* magarr, const CvArr* anglearr,
1740                CvArr* xarr, CvArr* yarr, int angle_in_degrees )
1741 {
1742     cv::Mat X, Y, Angle = cv::cvarrToMat(anglearr), Mag;
1743     if( magarr )
1744     {
1745         Mag = cv::cvarrToMat(magarr);
1746         CV_Assert( Mag.size() == Angle.size() && Mag.type() == Angle.type() );
1747     }
1748     if( xarr )
1749     {
1750         X = cv::cvarrToMat(xarr);
1751         CV_Assert( X.size() == Angle.size() && X.type() == Angle.type() );
1752     }
1753     if( yarr )
1754     {
1755         Y = cv::cvarrToMat(yarr);
1756         CV_Assert( Y.size() == Angle.size() && Y.type() == Angle.type() );
1757     }
1758
1759     cv::polarToCart( Mag, Angle, X, Y, angle_in_degrees != 0 );
1760 }
1761
1762 CV_IMPL void cvExp( const CvArr* srcarr, CvArr* dstarr )
1763 {
1764     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
1765     CV_Assert( src.type() == dst.type() && src.size() == dst.size() );
1766     cv::exp( src, dst );
1767 }
1768
1769 CV_IMPL void cvLog( const CvArr* srcarr, CvArr* dstarr )
1770 {
1771     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
1772     CV_Assert( src.type() == dst.type() && src.size() == dst.size() );
1773     cv::log( src, dst );
1774 }
1775
1776 CV_IMPL void cvPow( const CvArr* srcarr, CvArr* dstarr, double power )
1777 {
1778     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
1779     CV_Assert( src.type() == dst.type() && src.size() == dst.size() );
1780     cv::pow( src, power, dst );
1781 }
1782
1783 CV_IMPL int cvCheckArr( const CvArr* arr, int flags,
1784                         double minVal, double maxVal )
1785 {
1786     if( (flags & CV_CHECK_RANGE) == 0 )
1787         minVal = -DBL_MAX, maxVal = DBL_MAX;
1788     return cv::checkRange(cv::cvarrToMat(arr), (flags & CV_CHECK_QUIET) != 0, 0, minVal, maxVal );
1789 }
1790
1791
1792 /*
1793   Finds real roots of cubic, quadratic or linear equation.
1794   The original code has been taken from Ken Turkowski web page
1795   (http://www.worldserver.com/turk/opensource/) and adopted for OpenCV.
1796   Here is the copyright notice.
1797
1798   -----------------------------------------------------------------------
1799   Copyright (C) 1978-1999 Ken Turkowski. <turk@computer.org>
1800
1801     All rights reserved.
1802
1803     Warranty Information
1804       Even though I have reviewed this software, I make no warranty
1805       or representation, either express or implied, with respect to this
1806       software, its quality, accuracy, merchantability, or fitness for a
1807       particular purpose.  As a result, this software is provided "as is,"
1808       and you, its user, are assuming the entire risk as to its quality
1809       and accuracy.
1810
1811     This code may be used and freely distributed as long as it includes
1812     this copyright notice and the above warranty information.
1813   -----------------------------------------------------------------------
1814 */
1815 CV_IMPL int
1816 cvSolveCubic( const CvMat* coeffs, CvMat* roots )
1817 {
1818     int n = 0;
1819
1820     double a0 = 1., a1, a2, a3;
1821     double x0 = 0., x1 = 0., x2 = 0.;
1822     size_t step = 1;
1823     int coeff_count;
1824
1825     if( !CV_IS_MAT(coeffs) )
1826         CV_Error( !coeffs ? CV_StsNullPtr : CV_StsBadArg, "Input parameter is not a valid matrix" );
1827
1828     if( !CV_IS_MAT(roots) )
1829         CV_Error( !roots ? CV_StsNullPtr : CV_StsBadArg, "Output parameter is not a valid matrix" );
1830
1831     if( (CV_MAT_TYPE(coeffs->type) != CV_32FC1 && CV_MAT_TYPE(coeffs->type) != CV_64FC1) ||
1832         (CV_MAT_TYPE(roots->type) != CV_32FC1 && CV_MAT_TYPE(roots->type) != CV_64FC1) )
1833         CV_Error( CV_StsUnsupportedFormat,
1834         "Both matrices should be floating-point (single or double precision)" );
1835
1836     coeff_count = coeffs->rows + coeffs->cols - 1;
1837
1838     if( (coeffs->rows != 1 && coeffs->cols != 1) || (coeff_count != 3 && coeff_count != 4) )
1839         CV_Error( CV_StsBadSize,
1840         "The matrix of coefficients must be 1-dimensional vector of 3 or 4 elements" );
1841
1842     if( (roots->rows != 1 && roots->cols != 1) ||
1843         roots->rows + roots->cols - 1 != 3 )
1844         CV_Error( CV_StsBadSize,
1845         "The matrix of roots must be 1-dimensional vector of 3 elements" );
1846
1847     if( CV_MAT_TYPE(coeffs->type) == CV_32FC1 )
1848     {
1849         const float* c = coeffs->data.fl;
1850         if( coeffs->rows > 1 )
1851             step = coeffs->step/sizeof(c[0]);
1852         if( coeff_count == 4 )
1853             a0 = c[0], c += step;
1854         a1 = c[0];
1855         a2 = c[step];
1856         a3 = c[step*2];
1857     }
1858     else
1859     {
1860         const double* c = coeffs->data.db;
1861         if( coeffs->rows > 1 )
1862             step = coeffs->step/sizeof(c[0]);
1863         if( coeff_count == 4 )
1864             a0 = c[0], c += step;
1865         a1 = c[0];
1866         a2 = c[step];
1867         a3 = c[step*2];
1868     }
1869
1870     if( a0 == 0 )
1871     {
1872         if( a1 == 0 )
1873         {
1874             if( a2 == 0 )
1875                 n = a3 == 0 ? -1 : 0;
1876             else
1877             {
1878                 // linear equation
1879                 x0 = a3/a2;
1880                 n = 1;
1881             }
1882         }
1883         else
1884         {
1885             // quadratic equation
1886             double d = a2*a2 - 4*a1*a3;
1887             if( d >= 0 )
1888             {
1889                 d = sqrt(d);
1890                 double q = (-a2 + (a2 < 0 ? -d : d)) * 0.5;
1891                 x0 = q / a1;
1892                 x1 = a3 / q;
1893                 n = d > 0 ? 2 : 1;
1894             }
1895         }
1896     }
1897     else
1898     {
1899         a0 = 1./a0;
1900         a1 *= a0;
1901         a2 *= a0;
1902         a3 *= a0;
1903
1904         double Q = (a1 * a1 - 3 * a2) * (1./9);
1905         double R = (2 * a1 * a1 * a1 - 9 * a1 * a2 + 27 * a3) * (1./54);
1906         double Qcubed = Q * Q * Q;
1907         double d = Qcubed - R * R;
1908
1909         if( d >= 0 )
1910         {
1911             double theta = acos(R / sqrt(Qcubed));
1912             double sqrtQ = sqrt(Q);
1913             double t0 = -2 * sqrtQ;
1914             double t1 = theta * (1./3);
1915             double t2 = a1 * (1./3);
1916             x0 = t0 * cos(t1) - t2;
1917             x1 = t0 * cos(t1 + (2.*CV_PI/3)) - t2;
1918             x2 = t0 * cos(t1 + (4.*CV_PI/3)) - t2;
1919             n = 3;
1920         }
1921         else
1922         {
1923             double e;
1924             d = sqrt(-d);
1925             e = pow(d + fabs(R), 0.333333333333);
1926             if( R > 0 )
1927                 e = -e;
1928             x0 = (e + Q / e) - a1 * (1./3);
1929             n = 1;
1930         }
1931     }
1932
1933     step = 1;
1934
1935     if( CV_MAT_TYPE(roots->type) == CV_32FC1 )
1936     {
1937         float* r = roots->data.fl;
1938         if( roots->rows > 1 )
1939             step = roots->step/sizeof(r[0]);
1940         r[0] = (float)x0;
1941         r[step] = (float)x1;
1942         r[step*2] = (float)x2;
1943     }
1944     else
1945     {
1946         double* r = roots->data.db;
1947         if( roots->rows > 1 )
1948             step = roots->step/sizeof(r[0]);
1949         r[0] = x0;
1950         r[step] = x1;
1951         r[step*2] = x2;
1952     }
1953
1954     return n;
1955 }
1956
1957
1958 /*
1959   Finds real and complex roots of polynomials of any degree with real
1960   coefficients. The original code has been taken from Ken Turkowski's web
1961   page (http://www.worldserver.com/turk/opensource/) and adopted for OpenCV.
1962   Here is the copyright notice.
1963
1964   -----------------------------------------------------------------------
1965   Copyright (C) 1981-1999 Ken Turkowski. <turk@computer.org>
1966
1967   All rights reserved.
1968
1969   Warranty Information
1970   Even though I have reviewed this software, I make no warranty
1971   or representation, either express or implied, with respect to this
1972   software, its quality, accuracy, merchantability, or fitness for a
1973   particular purpose.  As a result, this software is provided "as is,"
1974   and you, its user, are assuming the entire risk as to its quality
1975   and accuracy.
1976
1977   This code may be used and freely distributed as long as it includes
1978   this copyright notice and the above warranty information.
1979 */
1980
1981
1982 /*******************************************************************************
1983  * FindPolynomialRoots
1984  *
1985  * The Bairstow and Newton correction formulae are used for a simultaneous
1986  * linear and quadratic iterated synthetic division.  The coefficients of
1987  * a polynomial of degree n are given as a[i] (i=0,i,..., n) where a[0] is
1988  * the constant term.  The coefficients are scaled by dividing them by
1989  * their geometric mean.  The Bairstow or Newton iteration method will
1990  * nearly always converge to the number of figures carried, fig, either to
1991  * root values or to their reciprocals.  If the simultaneous Newton and
1992  * Bairstow iteration fails to converge on root values or their
1993  * reciprocals in maxiter iterations, the convergence requirement will be
1994  * successively reduced by one decimal figure.  This program anticipates
1995  * and protects against loss of significance in the quadratic synthetic
1996  * division.  (Refer to "On Programming the Numerical Solution of
1997  * Polynomial Equations," by K. W. Ellenberger, Commun. ACM 3 (Dec. 1960),
1998  * 644-647.)  The real and imaginary part of each root is stated as u[i]
1999  * and v[i], respectively.
2000  *
2001  * ACM algorithm #30 - Numerical Solution of the Polynomial Equation
2002  * K. W. Ellenberger
2003  * Missle Division, North American Aviation, Downey, California
2004  * Converted to C, modified, optimized, and structured by
2005  * Ken Turkowski
2006  * CADLINC, Inc., Palo Alto, California
2007  *******************************************************************************/
2008
2009 #define MAXN 16
2010
2011 static void icvFindPolynomialRoots(const double *a, double *u, int n, int maxiter, int fig)
2012 {
2013   int i;
2014   int j;
2015   double h[MAXN + 3], b[MAXN + 3], c[MAXN + 3], d[MAXN + 3], e[MAXN + 3];
2016   // [-2 : n]
2017   double K, ps, qs, pt, qt, s, rev, r = 0;
2018   int t;
2019   double p = 0, q = 0, qq;
2020
2021   // Zero elements with negative indices
2022   b[2 + -1] = b[2 + -2] =
2023     c[2 + -1] = c[2 + -2] =
2024     d[2 + -1] = d[2 + -2] =
2025     e[2 + -1] = e[2 + -2] =
2026     h[2 + -1] = h[2 + -2] = 0.0;
2027
2028   // Copy polynomial coefficients to working storage
2029   for (j = n; j >= 0; j--)
2030     h[2 + j] = *a++; // Note reversal of coefficients
2031
2032   t = 1;
2033   K = pow(10.0, (double)(fig)); // Relative accuracy
2034
2035   for (; h[2 + n] == 0.0; n--) { // Look for zero high-order coeff.
2036     *u++ = 0.0;
2037     *u++ = 0.0;
2038   }
2039
2040  INIT:
2041   if (n == 0)
2042     return;
2043
2044   ps = qs = pt = qt = s = 0.0;
2045   rev = 1.0;
2046   K = pow(10.0, (double)(fig));
2047
2048   if (n == 1) {
2049     r = -h[2 + 1] / h[2 + 0];
2050     goto LINEAR;
2051   }
2052
2053   for (j = n; j >= 0; j--) // Find geometric mean of coeff's
2054     if (h[2 + j] != 0.0)
2055       s += log(fabs(h[2 + j]));
2056   s = exp(s / (n + 1));
2057
2058   for (j = n; j >= 0; j--) // Normalize coeff's by mean
2059     h[2 + j] /= s;
2060
2061   if (fabs(h[2 + 1] / h[2 + 0]) < fabs(h[2 + n - 1] / h[2 + n])) {
2062   REVERSE:
2063     t = -t;
2064     for (j = (n - 1) / 2; j >= 0; j--) {
2065       s = h[2 + j];
2066       h[2 + j] = h[2 + n - j];
2067       h[2 + n - j] = s;
2068     }
2069   }
2070   if (qs != 0.0) {
2071     p = ps;
2072     q = qs;
2073   } else {
2074     if (h[2 + n - 2] == 0.0) {
2075       q = 1.0;
2076       p = -2.0;
2077     } else {
2078       q = h[2 + n] / h[2 + n - 2];
2079       p = (h[2 + n - 1] - q * h[2 + n - 3]) / h[2 + n - 2];
2080     }
2081     if (n == 2)
2082       goto QADRTIC;
2083     r = 0.0;
2084   }
2085  ITERATE:
2086   for (i = maxiter; i > 0; i--) {
2087
2088     for (j = 0; j <= n; j++) { // BAIRSTOW
2089       b[2 + j] = h[2 + j] - p * b[2 + j - 1] - q * b[2 + j - 2];
2090       c[2 + j] = b[2 + j] - p * c[2 + j - 1] - q * c[2 + j - 2];
2091     }
2092     if ((h[2 + n - 1] != 0.0) && (b[2 + n - 1] != 0.0)) {
2093       if (fabs(h[2 + n - 1] / b[2 + n - 1]) >= K) {
2094         b[2 + n] = h[2 + n] - q * b[2 + n - 2];
2095       }
2096       if (b[2 + n] == 0.0)
2097         goto QADRTIC;
2098       if (K < fabs(h[2 + n] / b[2 + n]))
2099         goto QADRTIC;
2100     }
2101
2102     for (j = 0; j <= n; j++) { // NEWTON
2103       d[2 + j] = h[2 + j] + r * d[2 + j - 1]; // Calculate polynomial at r
2104       e[2 + j] = d[2 + j] + r * e[2 + j - 1]; // Calculate derivative at r
2105     }
2106     if (d[2 + n] == 0.0)
2107       goto LINEAR;
2108     if (K < fabs(h[2 + n] / d[2 + n]))
2109       goto LINEAR;
2110
2111     c[2 + n - 1] = -p * c[2 + n - 2] - q * c[2 + n - 3];
2112     s = c[2 + n - 2] * c[2 + n - 2] - c[2 + n - 1] * c[2 + n - 3];
2113     if (s == 0.0) {
2114       p -= 2.0;
2115       q *= (q + 1.0);
2116     } else {
2117       p += (b[2 + n - 1] * c[2 + n - 2] - b[2 + n] * c[2 + n - 3]) / s;
2118       q += (-b[2 + n - 1] * c[2 + n - 1] + b[2 + n] * c[2 + n - 2]) / s;
2119     }
2120     if (e[2 + n - 1] == 0.0)
2121       r -= 1.0; // Minimum step
2122     else
2123       r -= d[2 + n] / e[2 + n - 1]; // Newton's iteration
2124   }
2125   ps = pt;
2126   qs = qt;
2127   pt = p;
2128   qt = q;
2129   if (rev < 0.0)
2130     K /= 10.0;
2131   rev = -rev;
2132   goto REVERSE;
2133
2134  LINEAR:
2135   if (t < 0)
2136     r = 1.0 / r;
2137   n--;
2138   *u++ = r;
2139   *u++ = 0.0;
2140
2141   for (j = n; j >= 0; j--) { // Polynomial deflation by lin-nomial
2142     if ((d[2 + j] != 0.0) && (fabs(h[2 + j] / d[2 + j]) < K))
2143       h[2 + j] = d[2 + j];
2144     else
2145       h[2 + j] = 0.0;
2146   }
2147
2148   if (n == 0)
2149     return;
2150   goto ITERATE;
2151
2152  QADRTIC:
2153   if (t < 0) {
2154     p /= q;
2155     q = 1.0 / q;
2156   }
2157   n -= 2;
2158
2159   if (0.0 < (q - (p * p / 4.0))) { // Two complex roots
2160     s = sqrt(q - (p * p / 4.0));
2161     *u++ = -p / 2.0;
2162     *u++ = s;
2163     *u++ = -p / 2.0;
2164     *u++ = -s;
2165   } else { // Two real roots
2166     s = sqrt(((p * p / 4.0)) - q);
2167     if (p < 0.0)
2168       *u++ = qq = -p / 2.0 + s;
2169     else
2170       *u++ = qq = -p / 2.0 - s;
2171     *u++ = 0.0;
2172     *u++ = q / qq;
2173     *u++ = 0.0;
2174   }
2175
2176   for (j = n; j >= 0; j--) { // Polynomial deflation by quadratic
2177     if ((b[2 + j] != 0.0) && (fabs(h[2 + j] / b[2 + j]) < K))
2178       h[2 + j] = b[2 + j];
2179     else
2180       h[2 + j] = 0.0;
2181   }
2182   goto INIT;
2183 }
2184
2185 #undef MAXN
2186
2187 void cvSolvePoly(const CvMat* a, CvMat *r, int maxiter, int fig)
2188 {
2189     __BEGIN__;
2190
2191     int m, n;
2192     double *ad = 0, *rd = 0;
2193
2194     CV_FUNCNAME("cvSolvePoly");
2195
2196     CV_ASSERT(maxiter > 0);
2197     if (CV_MAT_TYPE(a->type) != CV_32FC1 &&
2198         CV_MAT_TYPE(a->type) != CV_64FC1)
2199         CV_Error(CV_StsUnsupportedFormat, "coeffs must be either CV_32FC1 or CV_64FC1");
2200     if (CV_MAT_TYPE(r->type) != CV_32FC2 &&
2201         CV_MAT_TYPE(r->type) != CV_64FC2)
2202         CV_Error(CV_StsUnsupportedFormat, "roots must be either CV_32FC2 or CV_64FC2");
2203     m = a->rows * a->cols;
2204     n = r->rows * r->cols;
2205
2206     if (m - 1 != n)
2207         CV_Error(CV_StsUnmatchedFormats, "must have n + 1 coefficients");
2208
2209     if( CV_MAT_DEPTH(a->type) == CV_32F || !CV_IS_MAT_CONT(a->type))
2210     {
2211         ad = (double*)cvStackAlloc(m*sizeof(ad[0]));
2212         CvMat _a = cvMat( a->rows, a->cols, CV_64F, ad );
2213         cvConvert( a, &_a );
2214     }
2215     else
2216         ad = a->data.db;
2217
2218     if( CV_MAT_DEPTH(r->type) == CV_32F || !CV_IS_MAT_CONT(r->type))
2219         rd = (double*)cvStackAlloc(n*sizeof(ad[0]));
2220     else
2221         rd = r->data.db;
2222
2223     icvFindPolynomialRoots( ad, rd, n, maxiter, fig);
2224     if( rd != r->data.db )
2225     {
2226         CvMat _r = cvMat( r->rows, r->cols, CV_64FC2, rd );
2227         cvConvert( &_r, r );
2228     }
2229
2230     __END__;
2231 }
2232
2233
2234 /* End of file. */