1/*M///////////////////////////////////////////////////////////////////////////////////////
2//
3//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4//
5//  By downloading, copying, installing or using the software you agree to this license.
6//  If you do not agree to this license, do not download, install,
7//  copy or use the software.
8//
9//
10//                        Intel License Agreement
11//                For Open Source Computer Vision Library
12//
13// Copyright (C) 2000, Intel Corporation, all rights reserved.
14// Third party copyrights are property of their respective owners.
15//
16// Redistribution and use in source and binary forms, with or without modification,
17// are permitted provided that the following conditions are met:
18//
19//   * Redistribution's of source code must retain the above copyright notice,
20//     this list of conditions and the following disclaimer.
21//
22//   * Redistribution's in binary form must reproduce the above copyright notice,
23//     this list of conditions and the following disclaimer in the documentation
24//     and/or other materials provided with the distribution.
25//
26//   * The name of Intel Corporation may not be used to endorse or promote products
27//     derived from this software without specific prior written permission.
28//
29// This software is provided by the copyright holders and contributors "as is" and
30// any express or implied warranties, including, but not limited to, the implied
31// warranties of merchantability and fitness for a particular purpose are disclaimed.
32// In no event shall the Intel Corporation or contributors be liable for any direct,
33// indirect, incidental, special, exemplary, or consequential damages
34// (including, but not limited to, procurement of substitute goods or services;
35// loss of use, data, or profits; or business interruption) however caused
36// and on any theory of liability, whether in contract, strict liability,
37// or tort (including negligence or otherwise) arising in any way out of
38// the use of this software, even if advised of the possibility of such damage.
39//
40//M*/
41
42#include "_cxcore.h"
43
44#ifdef HAVE_CONFIG_H
45#include <cvconfig.h>
46#endif
47
48#define ICV_MATH_BLOCK_SIZE  256
49
50#define _CV_SQRT_MAGIC     0xbe6f0000
51
52#define _CV_SQRT_MAGIC_DBL CV_BIG_UINT(0xbfcd460000000000)
53
54#define _CV_ATAN_CF0  (-15.8131890796f)
55#define _CV_ATAN_CF1  (61.0941945596f)
56#define _CV_ATAN_CF2  0.f /*(-0.140500406322f)*/
57
58static const float icvAtanTab[8] = { 0.f + _CV_ATAN_CF2, 90.f - _CV_ATAN_CF2,
59    180.f - _CV_ATAN_CF2, 90.f + _CV_ATAN_CF2,
60    360.f - _CV_ATAN_CF2, 270.f + _CV_ATAN_CF2,
61    180.f + _CV_ATAN_CF2, 270.f - _CV_ATAN_CF2
62};
63
64static const int icvAtanSign[8] =
65    { 0, 0x80000000, 0x80000000, 0, 0x80000000, 0, 0, 0x80000000 };
66
67CV_IMPL float
68cvFastArctan( float y, float x )
69{
70    Cv32suf _x, _y;
71    int ix, iy, ygx, idx;
72    double z;
73
74    _x.f = x; _y.f = y;
75    ix = _x.i; iy = _y.i;
76    idx = (ix < 0) * 2 + (iy < 0) * 4;
77
78    ix &= 0x7fffffff;
79    iy &= 0x7fffffff;
80
81    ygx = (iy <= ix) - 1;
82    idx -= ygx;
83
84    idx &= ((ix == 0) - 1) | ((iy == 0) - 1);
85
86    /* swap ix and iy if ix < iy */
87    ix ^= iy & ygx;
88    iy ^= ix & ygx;
89    ix ^= iy & ygx;
90
91    _y.i = iy ^ icvAtanSign[idx];
92
93    /* ix = ix != 0 ? ix : 1.f */
94    _x.i = ((ix ^ CV_1F) & ((ix == 0) - 1)) ^ CV_1F;
95
96    z = _y.f / _x.f;
97    return (float)((_CV_ATAN_CF0*fabs(z) + _CV_ATAN_CF1)*z + icvAtanTab[idx]);
98}
99
100
101IPCVAPI_IMPL( CvStatus, icvFastArctan_32f,
102    (const float *__y, const float *__x, float *angle, int len ), (__y, __x, angle, len) )
103{
104    int i = 0;
105    const int *y = (const int*)__y, *x = (const int*)__x;
106
107    if( !(y && x && angle && len >= 0) )
108        return CV_BADFACTOR_ERR;
109
110    /* unrolled by 4 loop */
111    for( ; i <= len - 4; i += 4 )
112    {
113        int j, idx[4];
114        float xf[4], yf[4];
115        double d = 1.;
116
117        /* calc numerators and denominators */
118        for( j = 0; j < 4; j++ )
119        {
120            int ix = x[i + j], iy = y[i + j];
121            int ygx, k = (ix < 0) * 2 + (iy < 0) * 4;
122            Cv32suf _x, _y;
123
124            ix &= 0x7fffffff;
125            iy &= 0x7fffffff;
126
127            ygx = (iy <= ix) - 1;
128            k -= ygx;
129
130            k &= ((ix == 0) - 1) | ((iy == 0) - 1);
131
132            /* swap ix and iy if ix < iy */
133            ix ^= iy & ygx;
134            iy ^= ix & ygx;
135            ix ^= iy & ygx;
136
137            _y.i = iy ^ icvAtanSign[k];
138
139            /* ix = ix != 0 ? ix : 1.f */
140            _x.i = ((ix ^ CV_1F) & ((ix == 0) - 1)) ^ CV_1F;
141            idx[j] = k;
142            yf[j] = _y.f;
143            d *= (xf[j] = _x.f);
144        }
145
146        d = 1. / d;
147
148        {
149            double b = xf[2] * xf[3], a = xf[0] * xf[1];
150
151            float z0 = (float) (yf[0] * xf[1] * b * d);
152            float z1 = (float) (yf[1] * xf[0] * b * d);
153            float z2 = (float) (yf[2] * xf[3] * a * d);
154            float z3 = (float) (yf[3] * xf[2] * a * d);
155
156            z0 = (float)((_CV_ATAN_CF0*fabs(z0) + _CV_ATAN_CF1)*z0 + icvAtanTab[idx[0]]);
157            z1 = (float)((_CV_ATAN_CF0*fabs(z1) + _CV_ATAN_CF1)*z1 + icvAtanTab[idx[1]]);
158            z2 = (float)((_CV_ATAN_CF0*fabs(z2) + _CV_ATAN_CF1)*z2 + icvAtanTab[idx[2]]);
159            z3 = (float)((_CV_ATAN_CF0*fabs(z3) + _CV_ATAN_CF1)*z3 + icvAtanTab[idx[3]]);
160
161            angle[i] = z0;
162            angle[i+1] = z1;
163            angle[i+2] = z2;
164            angle[i+3] = z3;
165        }
166    }
167
168    /* process the rest */
169    for( ; i < len; i++ )
170        angle[i] = cvFastArctan( __y[i], __x[i] );
171
172    return CV_OK;
173}
174
175
176/* ************************************************************************** *\
177   Fast cube root by Ken Turkowski
178   (http://www.worldserver.com/turk/computergraphics/papers.html)
179\* ************************************************************************** */
180CV_IMPL  float  cvCbrt( float value )
181{
182    float fr;
183    Cv32suf v, m;
184    int ix, s;
185    int ex, shx;
186
187    v.f = value;
188    ix = v.i & 0x7fffffff;
189    s = v.i & 0x80000000;
190    ex = (ix >> 23) - 127;
191    shx = ex % 3;
192    shx -= shx >= 0 ? 3 : 0;
193    ex = (ex - shx) / 3; /* exponent of cube root */
194    v.i = (ix & ((1<<23)-1)) | ((shx + 127)<<23);
195    fr = v.f;
196
197    /* 0.125 <= fr < 1.0 */
198    /* Use quartic rational polynomial with error < 2^(-24) */
199    fr = (float)(((((45.2548339756803022511987494 * fr +
200    192.2798368355061050458134625) * fr +
201    119.1654824285581628956914143) * fr +
202    13.43250139086239872172837314) * fr +
203    0.1636161226585754240958355063)/
204    ((((14.80884093219134573786480845 * fr +
205    151.9714051044435648658557668) * fr +
206    168.5254414101568283957668343) * fr +
207    33.9905941350215598754191872) * fr +
208    1.0));
209
210    /* fr *= 2^ex * sign */
211    m.f = value;
212    v.f = fr;
213    v.i = (v.i + (ex << 23) + s) & (m.i*2 != 0 ? -1 : 0);
214    return v.f;
215}
216
217//static const double _0_5 = 0.5, _1_5 = 1.5;
218
219IPCVAPI_IMPL( CvStatus, icvInvSqrt_32f, (const float *src, float *dst, int len), (src, dst, len) )
220{
221    int i = 0;
222
223    if( !(src && dst && len >= 0) )
224        return CV_BADFACTOR_ERR;
225
226    for( ; i < len; i++ )
227        dst[i] = (float)(1.f/sqrt(src[i]));
228
229    return CV_OK;
230}
231
232
233IPCVAPI_IMPL( CvStatus, icvSqrt_32f, (const float *src, float *dst, int len), (src, dst, len) )
234{
235    int i = 0;
236
237    if( !(src && dst && len >= 0) )
238        return CV_BADFACTOR_ERR;
239
240    for( ; i < len; i++ )
241        dst[i] = (float)sqrt(src[i]);
242
243    return CV_OK;
244}
245
246
247IPCVAPI_IMPL( CvStatus, icvSqrt_64f, (const double *src, double *dst, int len), (src, dst, len) )
248{
249    int i = 0;
250
251    if( !(src && dst && len >= 0) )
252        return CV_BADFACTOR_ERR;
253
254    for( ; i < len; i++ )
255        dst[i] = sqrt(src[i]);
256
257    return CV_OK;
258}
259
260
261IPCVAPI_IMPL( CvStatus, icvInvSqrt_64f, (const double *src, double *dst, int len), (src, dst, len) )
262{
263    int i = 0;
264
265    if( !(src && dst && len >= 0) )
266        return CV_BADFACTOR_ERR;
267
268    for( ; i < len; i++ )
269        dst[i] = 1./sqrt(src[i]);
270
271    return CV_OK;
272}
273
274
275#define ICV_DEF_SQR_MAGNITUDE_FUNC(flavor, arrtype, magtype)\
276static CvStatus CV_STDCALL                                  \
277icvSqrMagnitude_##flavor(const arrtype* x, const arrtype* y,\
278                         magtype* mag, int len)             \
279{                                                           \
280    int i;                                                  \
281                                                            \
282    for( i = 0; i <= len - 4; i += 4 )                      \
283    {                                                       \
284        magtype x0 = (magtype)x[i], y0 = (magtype)y[i];     \
285        magtype x1 = (magtype)x[i+1], y1 = (magtype)y[i+1]; \
286                                                            \
287        x0 = x0*x0 + y0*y0;                                 \
288        x1 = x1*x1 + y1*y1;                                 \
289        mag[i] = x0;                                        \
290        mag[i+1] = x1;                                      \
291        x0 = (magtype)x[i+2], y0 = (magtype)y[i+2];         \
292        x1 = (magtype)x[i+3], y1 = (magtype)y[i+3];         \
293        x0 = x0*x0 + y0*y0;                                 \
294        x1 = x1*x1 + y1*y1;                                 \
295        mag[i+2] = x0;                                      \
296        mag[i+3] = x1;                                      \
297    }                                                       \
298                                                            \
299    for( ; i < len; i++ )                                   \
300    {                                                       \
301        magtype x0 = (magtype)x[i], y0 = (magtype)y[i];     \
302        mag[i] = x0*x0 + y0*y0;                             \
303    }                                                       \
304                                                            \
305    return CV_OK;                                           \
306}
307
308
309ICV_DEF_SQR_MAGNITUDE_FUNC( 32f, float, float )
310ICV_DEF_SQR_MAGNITUDE_FUNC( 64f, double, double )
311
312/****************************************************************************************\
313*                                  Cartezian -> Polar                                    *
314\****************************************************************************************/
315
316CV_IMPL void
317cvCartToPolar( const CvArr* xarr, const CvArr* yarr,
318               CvArr* magarr, CvArr* anglearr,
319               int angle_in_degrees )
320{
321    CV_FUNCNAME( "cvCartToPolar" );
322
323    __BEGIN__;
324
325    float* mag_buffer = 0;
326    float* x_buffer = 0;
327    float* y_buffer = 0;
328    int block_size = 0;
329    CvMat xstub, *xmat = (CvMat*)xarr;
330    CvMat ystub, *ymat = (CvMat*)yarr;
331    CvMat magstub, *mag = (CvMat*)magarr;
332    CvMat anglestub, *angle = (CvMat*)anglearr;
333    int coi1 = 0, coi2 = 0, coi3 = 0, coi4 = 0;
334    int depth;
335    CvSize size;
336    int x, y;
337    int cont_flag = CV_MAT_CONT_FLAG;
338
339    if( !CV_IS_MAT(xmat))
340        CV_CALL( xmat = cvGetMat( xmat, &xstub, &coi1 ));
341
342    if( !CV_IS_MAT(ymat))
343        CV_CALL( ymat = cvGetMat( ymat, &ystub, &coi2 ));
344
345    if( !CV_ARE_TYPES_EQ( xmat, ymat ) )
346        CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
347
348    if( !CV_ARE_SIZES_EQ( xmat, ymat ) )
349        CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
350
351    depth = CV_MAT_DEPTH( xmat->type );
352    if( depth < CV_32F )
353        CV_ERROR( CV_StsUnsupportedFormat, "" );
354
355    if( mag )
356    {
357        CV_CALL( mag = cvGetMat( mag, &magstub, &coi3 ));
358
359        if( !CV_ARE_TYPES_EQ( mag, xmat ) )
360            CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
361
362        if( !CV_ARE_SIZES_EQ( mag, xmat ) )
363            CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
364        cont_flag = mag->type;
365    }
366
367    if( angle )
368    {
369        CV_CALL( angle = cvGetMat( angle, &anglestub, &coi4 ));
370
371        if( !CV_ARE_TYPES_EQ( angle, xmat ) )
372            CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
373
374        if( !CV_ARE_SIZES_EQ( angle, xmat ) )
375            CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
376        cont_flag &= angle->type;
377    }
378
379    if( coi1 != 0 || coi2 != 0 || coi3 != 0 || coi4 != 0 )
380        CV_ERROR( CV_BadCOI, "" );
381
382    size = cvGetMatSize(xmat);
383    size.width *= CV_MAT_CN(xmat->type);
384
385    if( CV_IS_MAT_CONT( xmat->type & ymat->type & cont_flag ))
386    {
387        size.width *= size.height;
388        size.height = 1;
389    }
390
391    block_size = MIN( size.width, ICV_MATH_BLOCK_SIZE );
392    if( depth == CV_64F && angle )
393    {
394        x_buffer = (float*)cvStackAlloc( block_size*sizeof(float));
395        y_buffer = (float*)cvStackAlloc( block_size*sizeof(float));
396    }
397    else if( depth == CV_32F && mag )
398    {
399        mag_buffer = (float*)cvStackAlloc( block_size*sizeof(float));
400    }
401
402    if( depth == CV_32F )
403    {
404        for( y = 0; y < size.height; y++ )
405        {
406            float* x_data = (float*)(xmat->data.ptr + xmat->step*y);
407            float* y_data = (float*)(ymat->data.ptr + ymat->step*y);
408            float* mag_data = mag ? (float*)(mag->data.ptr + mag->step*y) : 0;
409            float* angle_data = angle ? (float*)(angle->data.ptr + angle->step*y) : 0;
410
411            for( x = 0; x < size.width; x += block_size )
412            {
413                int len = MIN( size.width - x, block_size );
414
415                if( mag )
416                    icvSqrMagnitude_32f( x_data + x, y_data + x, mag_buffer, len );
417
418                if( angle )
419                {
420                    icvFastArctan_32f( y_data + x, x_data + x, angle_data + x, len );
421                    if( !angle_in_degrees )
422                        icvScale_32f( angle_data + x, angle_data + x,
423                                      len, (float)(CV_PI/180.), 0 );
424                }
425
426                if( mag )
427                    icvSqrt_32f( mag_buffer, mag_data + x, len );
428            }
429        }
430    }
431    else
432    {
433        for( y = 0; y < size.height; y++ )
434        {
435            double* x_data = (double*)(xmat->data.ptr + xmat->step*y);
436            double* y_data = (double*)(ymat->data.ptr + ymat->step*y);
437            double* mag_data = mag ? (double*)(mag->data.ptr + mag->step*y) : 0;
438            double* angle_data = angle ? (double*)(angle->data.ptr + angle->step*y) : 0;
439
440            for( x = 0; x < size.width; x += block_size )
441            {
442                int len = MIN( size.width - x, block_size );
443
444                if( angle )
445                {
446                    icvCvt_64f32f( x_data + x, x_buffer, len );
447                    icvCvt_64f32f( y_data + x, y_buffer, len );
448                }
449
450                if( mag )
451                {
452                    icvSqrMagnitude_64f( x_data + x, y_data + x, mag_data + x, len );
453                    icvSqrt_64f( mag_data + x, mag_data + x, len );
454                }
455
456                if( angle )
457                {
458                    icvFastArctan_32f( y_buffer, x_buffer, x_buffer, len );
459                    if( !angle_in_degrees )
460                        icvScale_32f( x_buffer, x_buffer, len, (float)(CV_PI/180.), 0 );
461                    icvCvt_32f64f( x_buffer, angle_data + x, len );
462                }
463            }
464        }
465    }
466
467    __END__;
468}
469
470
471/****************************************************************************************\
472*                                  Polar -> Cartezian                                    *
473\****************************************************************************************/
474
475static CvStatus CV_STDCALL
476icvSinCos_32f( const float *angle,float *sinval, float* cosval,
477                int len, int angle_in_degrees )
478{
479    const int N = 64;
480
481    static const double sin_table[] =
482    {
483     0.00000000000000000000,     0.09801714032956060400,
484     0.19509032201612825000,     0.29028467725446233000,
485     0.38268343236508978000,     0.47139673682599764000,
486     0.55557023301960218000,     0.63439328416364549000,
487     0.70710678118654746000,     0.77301045336273699000,
488     0.83146961230254524000,     0.88192126434835494000,
489     0.92387953251128674000,     0.95694033573220894000,
490     0.98078528040323043000,     0.99518472667219682000,
491     1.00000000000000000000,     0.99518472667219693000,
492     0.98078528040323043000,     0.95694033573220894000,
493     0.92387953251128674000,     0.88192126434835505000,
494     0.83146961230254546000,     0.77301045336273710000,
495     0.70710678118654757000,     0.63439328416364549000,
496     0.55557023301960218000,     0.47139673682599786000,
497     0.38268343236508989000,     0.29028467725446239000,
498     0.19509032201612861000,     0.09801714032956082600,
499     0.00000000000000012246,    -0.09801714032956059000,
500    -0.19509032201612836000,    -0.29028467725446211000,
501    -0.38268343236508967000,    -0.47139673682599764000,
502    -0.55557023301960196000,    -0.63439328416364527000,
503    -0.70710678118654746000,    -0.77301045336273666000,
504    -0.83146961230254524000,    -0.88192126434835494000,
505    -0.92387953251128652000,    -0.95694033573220882000,
506    -0.98078528040323032000,    -0.99518472667219693000,
507    -1.00000000000000000000,    -0.99518472667219693000,
508    -0.98078528040323043000,    -0.95694033573220894000,
509    -0.92387953251128663000,    -0.88192126434835505000,
510    -0.83146961230254546000,    -0.77301045336273688000,
511    -0.70710678118654768000,    -0.63439328416364593000,
512    -0.55557023301960218000,    -0.47139673682599792000,
513    -0.38268343236509039000,    -0.29028467725446250000,
514    -0.19509032201612872000,    -0.09801714032956050600,
515    };
516
517    static const double k2 = (2*CV_PI)/N;
518
519    static const double sin_a0 = -0.166630293345647*k2*k2*k2;
520    static const double sin_a2 = k2;
521
522    static const double cos_a0 = -0.499818138450326*k2*k2;
523    /*static const double cos_a2 =  1;*/
524
525    double k1;
526    int i;
527
528    if( !angle_in_degrees )
529        k1 = N/(2*CV_PI);
530    else
531        k1 = N/360.;
532
533    for( i = 0; i < len; i++ )
534    {
535        double t = angle[i]*k1;
536        int it = cvRound(t);
537        t -= it;
538        int sin_idx = it & (N - 1);
539        int cos_idx = (N/4 - sin_idx) & (N - 1);
540
541        double sin_b = (sin_a0*t*t + sin_a2)*t;
542        double cos_b = cos_a0*t*t + 1;
543
544        double sin_a = sin_table[sin_idx];
545        double cos_a = sin_table[cos_idx];
546
547        double sin_val = sin_a*cos_b + cos_a*sin_b;
548        double cos_val = cos_a*cos_b - sin_a*sin_b;
549
550        sinval[i] = (float)sin_val;
551        cosval[i] = (float)cos_val;
552    }
553
554    return CV_OK;
555}
556
557
558CV_IMPL void
559cvPolarToCart( const CvArr* magarr, const CvArr* anglearr,
560               CvArr* xarr, CvArr* yarr, int angle_in_degrees )
561{
562    CV_FUNCNAME( "cvPolarToCart" );
563
564    __BEGIN__;
565
566    float* x_buffer = 0;
567    float* y_buffer = 0;
568    int block_size = 0;
569    CvMat xstub, *xmat = (CvMat*)xarr;
570    CvMat ystub, *ymat = (CvMat*)yarr;
571    CvMat magstub, *mag = (CvMat*)magarr;
572    CvMat anglestub, *angle = (CvMat*)anglearr;
573    int coi1 = 0, coi2 = 0, coi3 = 0, coi4 = 0;
574    int depth;
575    CvSize size;
576    int x, y;
577    int cont_flag;
578
579    if( !CV_IS_MAT(angle))
580        CV_CALL( angle = cvGetMat( angle, &anglestub, &coi4 ));
581
582    depth = CV_MAT_DEPTH( angle->type );
583    if( depth < CV_32F )
584        CV_ERROR( CV_StsUnsupportedFormat, "" );
585    cont_flag = angle->type;
586
587    if( mag )
588    {
589        if( !CV_IS_MAT(mag))
590            CV_CALL( mag = cvGetMat( mag, &magstub, &coi3 ));
591
592        if( !CV_ARE_TYPES_EQ( angle, mag ) )
593            CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
594
595        if( !CV_ARE_SIZES_EQ( angle, mag ) )
596            CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
597
598        cont_flag &= mag->type;
599    }
600
601    if( xmat )
602    {
603        if( !CV_IS_MAT(xmat))
604            CV_CALL( xmat = cvGetMat( xmat, &xstub, &coi1 ));
605
606        if( !CV_ARE_TYPES_EQ( angle, xmat ) )
607            CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
608
609        if( !CV_ARE_SIZES_EQ( angle, xmat ) )
610            CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
611
612        cont_flag &= xmat->type;
613    }
614
615    if( ymat )
616    {
617        if( !CV_IS_MAT(ymat))
618            CV_CALL( ymat = cvGetMat( ymat, &ystub, &coi2 ));
619
620        if( !CV_ARE_TYPES_EQ( angle, ymat ) )
621            CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
622
623        if( !CV_ARE_SIZES_EQ( angle, ymat ) )
624            CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
625
626        cont_flag &= ymat->type;
627    }
628
629    if( coi1 != 0 || coi2 != 0 || coi3 != 0 || coi4 != 0 )
630        CV_ERROR( CV_BadCOI, "" );
631
632    size = cvGetMatSize(angle);
633    size.width *= CV_MAT_CN(angle->type);
634
635    if( CV_IS_MAT_CONT( cont_flag ))
636    {
637        size.width *= size.height;
638        size.height = 1;
639    }
640
641    block_size = MIN( size.width, ICV_MATH_BLOCK_SIZE );
642    x_buffer = (float*)cvStackAlloc( block_size*sizeof(float));
643    y_buffer = (float*)cvStackAlloc( block_size*sizeof(float));
644
645    if( depth == CV_32F )
646    {
647        for( y = 0; y < size.height; y++ )
648        {
649            float* x_data = (float*)(xmat ? xmat->data.ptr + xmat->step*y : 0);
650            float* y_data = (float*)(ymat ? ymat->data.ptr + ymat->step*y : 0);
651            float* mag_data = (float*)(mag ? mag->data.ptr + mag->step*y : 0);
652            float* angle_data = (float*)(angle->data.ptr + angle->step*y);
653
654            for( x = 0; x < size.width; x += block_size )
655            {
656                int i, len = MIN( size.width - x, block_size );
657
658                icvSinCos_32f( angle_data+x, y_buffer, x_buffer, len, angle_in_degrees );
659
660                for( i = 0; i < len; i++ )
661                {
662                    float tx = x_buffer[i];
663                    float ty = y_buffer[i];
664
665                    if( mag_data )
666                    {
667                        float magval = mag_data[x+i];
668                        tx *= magval;
669                        ty *= magval;
670                    }
671
672                    if( xmat )
673                        x_data[x+i] = tx;
674                    if( ymat )
675                        y_data[x+i] = ty;
676                }
677            }
678        }
679    }
680    else
681    {
682        for( y = 0; y < size.height; y++ )
683        {
684            double* x_data = (double*)(xmat ? xmat->data.ptr + xmat->step*y : 0);
685            double* y_data = (double*)(ymat ? ymat->data.ptr + ymat->step*y : 0);
686            double* mag_data = (double*)(mag ? mag->data.ptr + mag->step*y : 0);
687            double* angle_data = (double*)(angle->data.ptr + angle->step*y);
688            double C = angle_in_degrees ? CV_PI/180. : 1;
689
690            for( x = 0; x < size.width; x++ )
691            {
692                double phi = angle_data[x]*C;
693                double magval = mag_data ? mag_data[x] : 1.;
694                if( xmat )
695                    x_data[x] = cos(phi)*magval;
696                if( ymat )
697                    y_data[x] = sin(phi)*magval;
698            }
699        }
700    }
701
702    __END__;
703}
704
705/****************************************************************************************\
706*                                          E X P                                         *
707\****************************************************************************************/
708
709typedef union
710{
711    struct {
712#if ( defined( WORDS_BIGENDIAN ) && !defined( OPENCV_UNIVERSAL_BUILD ) ) || defined( __BIG_ENDIAN__ )
713        int hi;
714        int lo;
715#else
716        int lo;
717        int hi;
718#endif
719    } i;
720    double d;
721}
722DBLINT;
723
724#define EXPTAB_SCALE 6
725#define EXPTAB_MASK  ((1 << EXPTAB_SCALE) - 1)
726
727#define EXPPOLY_32F_A0 .9670371139572337719125840413672004409288e-2
728
729static const double icvExpTab[] = {
730    1.0 * EXPPOLY_32F_A0,
731    1.0108892860517004600204097905619 * EXPPOLY_32F_A0,
732    1.0218971486541166782344801347833 * EXPPOLY_32F_A0,
733    1.0330248790212284225001082839705 * EXPPOLY_32F_A0,
734    1.0442737824274138403219664787399 * EXPPOLY_32F_A0,
735    1.0556451783605571588083413251529 * EXPPOLY_32F_A0,
736    1.0671404006768236181695211209928 * EXPPOLY_32F_A0,
737    1.0787607977571197937406800374385 * EXPPOLY_32F_A0,
738    1.0905077326652576592070106557607 * EXPPOLY_32F_A0,
739    1.1023825833078409435564142094256 * EXPPOLY_32F_A0,
740    1.1143867425958925363088129569196 * EXPPOLY_32F_A0,
741    1.126521618608241899794798643787 * EXPPOLY_32F_A0,
742    1.1387886347566916537038302838415 * EXPPOLY_32F_A0,
743    1.151189229952982705817759635202 * EXPPOLY_32F_A0,
744    1.1637248587775775138135735990922 * EXPPOLY_32F_A0,
745    1.1763969916502812762846457284838 * EXPPOLY_32F_A0,
746    1.1892071150027210667174999705605 * EXPPOLY_32F_A0,
747    1.2021567314527031420963969574978 * EXPPOLY_32F_A0,
748    1.2152473599804688781165202513388 * EXPPOLY_32F_A0,
749    1.2284805361068700056940089577928 * EXPPOLY_32F_A0,
750    1.2418578120734840485936774687266 * EXPPOLY_32F_A0,
751    1.2553807570246910895793906574423 * EXPPOLY_32F_A0,
752    1.2690509571917332225544190810323 * EXPPOLY_32F_A0,
753    1.2828700160787782807266697810215 * EXPPOLY_32F_A0,
754    1.2968395546510096659337541177925 * EXPPOLY_32F_A0,
755    1.3109612115247643419229917863308 * EXPPOLY_32F_A0,
756    1.3252366431597412946295370954987 * EXPPOLY_32F_A0,
757    1.3396675240533030053600306697244 * EXPPOLY_32F_A0,
758    1.3542555469368927282980147401407 * EXPPOLY_32F_A0,
759    1.3690024229745906119296011329822 * EXPPOLY_32F_A0,
760    1.3839098819638319548726595272652 * EXPPOLY_32F_A0,
761    1.3989796725383111402095281367152 * EXPPOLY_32F_A0,
762    1.4142135623730950488016887242097 * EXPPOLY_32F_A0,
763    1.4296133383919700112350657782751 * EXPPOLY_32F_A0,
764    1.4451808069770466200370062414717 * EXPPOLY_32F_A0,
765    1.4609177941806469886513028903106 * EXPPOLY_32F_A0,
766    1.476826145939499311386907480374 * EXPPOLY_32F_A0,
767    1.4929077282912648492006435314867 * EXPPOLY_32F_A0,
768    1.5091644275934227397660195510332 * EXPPOLY_32F_A0,
769    1.5255981507445383068512536895169 * EXPPOLY_32F_A0,
770    1.5422108254079408236122918620907 * EXPPOLY_32F_A0,
771    1.5590044002378369670337280894749 * EXPPOLY_32F_A0,
772    1.5759808451078864864552701601819 * EXPPOLY_32F_A0,
773    1.5931421513422668979372486431191 * EXPPOLY_32F_A0,
774    1.6104903319492543081795206673574 * EXPPOLY_32F_A0,
775    1.628027421857347766848218522014 * EXPPOLY_32F_A0,
776    1.6457554781539648445187567247258 * EXPPOLY_32F_A0,
777    1.6636765803267364350463364569764 * EXPPOLY_32F_A0,
778    1.6817928305074290860622509524664 * EXPPOLY_32F_A0,
779    1.7001063537185234695013625734975 * EXPPOLY_32F_A0,
780    1.7186192981224779156293443764563 * EXPPOLY_32F_A0,
781    1.7373338352737062489942020818722 * EXPPOLY_32F_A0,
782    1.7562521603732994831121606193753 * EXPPOLY_32F_A0,
783    1.7753764925265212525505592001993 * EXPPOLY_32F_A0,
784    1.7947090750031071864277032421278 * EXPPOLY_32F_A0,
785    1.8142521755003987562498346003623 * EXPPOLY_32F_A0,
786    1.8340080864093424634870831895883 * EXPPOLY_32F_A0,
787    1.8539791250833855683924530703377 * EXPPOLY_32F_A0,
788    1.8741676341102999013299989499544 * EXPPOLY_32F_A0,
789    1.8945759815869656413402186534269 * EXPPOLY_32F_A0,
790    1.9152065613971472938726112702958 * EXPPOLY_32F_A0,
791    1.9360617934922944505980559045667 * EXPPOLY_32F_A0,
792    1.9571441241754002690183222516269 * EXPPOLY_32F_A0,
793    1.9784560263879509682582499181312 * EXPPOLY_32F_A0,
794};
795
796static const double exp_prescale = 1.4426950408889634073599246810019 * (1 << EXPTAB_SCALE);
797static const double exp_postscale = 1./(1 << EXPTAB_SCALE);
798static const double exp_max_val = 3000.*(1 << EXPTAB_SCALE); // log10(DBL_MAX) < 3000
799
800IPCVAPI_IMPL( CvStatus, icvExp_32f, ( const float *_x, float *y, int n ), (_x, y, n) )
801{
802    static const double
803        EXPPOLY_32F_A4 = 1.000000000000002438532970795181890933776 / EXPPOLY_32F_A0,
804        EXPPOLY_32F_A3 = .6931471805521448196800669615864773144641 / EXPPOLY_32F_A0,
805        EXPPOLY_32F_A2 = .2402265109513301490103372422686535526573 / EXPPOLY_32F_A0,
806        EXPPOLY_32F_A1 = .5550339366753125211915322047004666939128e-1 / EXPPOLY_32F_A0;
807
808    #undef EXPPOLY
809    #define EXPPOLY(x)  \
810        (((((x) + EXPPOLY_32F_A1)*(x) + EXPPOLY_32F_A2)*(x) + EXPPOLY_32F_A3)*(x) + EXPPOLY_32F_A4)
811
812    int i = 0;
813    DBLINT buf[4];
814    const Cv32suf* x = (const Cv32suf*)_x;
815
816    if( !x || !y )
817        return CV_NULLPTR_ERR;
818    if( n <= 0 )
819        return CV_BADSIZE_ERR;
820
821    buf[0].i.lo = buf[1].i.lo = buf[2].i.lo = buf[3].i.lo = 0;
822
823    for( ; i <= n - 4; i += 4 )
824    {
825        double x0 = x[i].f * exp_prescale;
826        double x1 = x[i + 1].f * exp_prescale;
827        double x2 = x[i + 2].f * exp_prescale;
828        double x3 = x[i + 3].f * exp_prescale;
829        int val0, val1, val2, val3, t;
830
831        if( ((x[i].i >> 23) & 255) > 127 + 10 )
832            x0 = x[i].i < 0 ? -exp_max_val : exp_max_val;
833
834        if( ((x[i+1].i >> 23) & 255) > 127 + 10 )
835            x1 = x[i+1].i < 0 ? -exp_max_val : exp_max_val;
836
837        if( ((x[i+2].i >> 23) & 255) > 127 + 10 )
838            x2 = x[i+2].i < 0 ? -exp_max_val : exp_max_val;
839
840        if( ((x[i+3].i >> 23) & 255) > 127 + 10 )
841            x3 = x[i+3].i < 0 ? -exp_max_val : exp_max_val;
842
843        val0 = cvRound(x0);
844        val1 = cvRound(x1);
845        val2 = cvRound(x2);
846        val3 = cvRound(x3);
847
848        x0 = (x0 - val0)*exp_postscale;
849        x1 = (x1 - val1)*exp_postscale;
850        x2 = (x2 - val2)*exp_postscale;
851        x3 = (x3 - val3)*exp_postscale;
852
853        t = (val0 >> EXPTAB_SCALE) + 1023;
854        t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
855        buf[0].i.hi = t << 20;
856
857        t = (val1 >> EXPTAB_SCALE) + 1023;
858        t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
859        buf[1].i.hi = t << 20;
860
861        t = (val2 >> EXPTAB_SCALE) + 1023;
862        t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
863        buf[2].i.hi = t << 20;
864
865        t = (val3 >> EXPTAB_SCALE) + 1023;
866        t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
867        buf[3].i.hi = t << 20;
868
869        x0 = buf[0].d * icvExpTab[val0 & EXPTAB_MASK] * EXPPOLY( x0 );
870        x1 = buf[1].d * icvExpTab[val1 & EXPTAB_MASK] * EXPPOLY( x1 );
871
872        y[i] = (float)x0;
873        y[i + 1] = (float)x1;
874
875        x2 = buf[2].d * icvExpTab[val2 & EXPTAB_MASK] * EXPPOLY( x2 );
876        x3 = buf[3].d * icvExpTab[val3 & EXPTAB_MASK] * EXPPOLY( x3 );
877
878        y[i + 2] = (float)x2;
879        y[i + 3] = (float)x3;
880    }
881
882    for( ; i < n; i++ )
883    {
884        double x0 = x[i].f * exp_prescale;
885        int val0, t;
886
887        if( ((x[i].i >> 23) & 255) > 127 + 10 )
888            x0 = x[i].i < 0 ? -exp_max_val : exp_max_val;
889
890        val0 = cvRound(x0);
891        t = (val0 >> EXPTAB_SCALE) + 1023;
892        t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
893
894        buf[0].i.hi = t << 20;
895        x0 = (x0 - val0)*exp_postscale;
896
897        y[i] = (float)(buf[0].d * icvExpTab[val0 & EXPTAB_MASK] * EXPPOLY(x0));
898    }
899
900    return CV_OK;
901}
902
903
904IPCVAPI_IMPL( CvStatus, icvExp_64f, ( const double *_x, double *y, int n ), (_x, y, n) )
905{
906    static const double
907        A5 = .99999999999999999998285227504999 / EXPPOLY_32F_A0,
908        A4 = .69314718055994546743029643825322 / EXPPOLY_32F_A0,
909        A3 = .24022650695886477918181338054308 / EXPPOLY_32F_A0,
910        A2 = .55504108793649567998466049042729e-1 / EXPPOLY_32F_A0,
911        A1 = .96180973140732918010002372686186e-2 / EXPPOLY_32F_A0,
912        A0 = .13369713757180123244806654839424e-2 / EXPPOLY_32F_A0;
913
914    #undef EXPPOLY
915    #define EXPPOLY(x)  (((((A0*(x) + A1)*(x) + A2)*(x) + A3)*(x) + A4)*(x) + A5)
916
917    int i = 0;
918    DBLINT buf[4];
919    const Cv64suf* x = (const Cv64suf*)_x;
920
921    if( !x || !y )
922        return CV_NULLPTR_ERR;
923    if( n <= 0 )
924        return CV_BADSIZE_ERR;
925
926    buf[0].i.lo = buf[1].i.lo = buf[2].i.lo = buf[3].i.lo = 0;
927
928    for( ; i <= n - 4; i += 4 )
929    {
930        double x0 = x[i].f * exp_prescale;
931        double x1 = x[i + 1].f * exp_prescale;
932        double x2 = x[i + 2].f * exp_prescale;
933        double x3 = x[i + 3].f * exp_prescale;
934
935        double y0, y1, y2, y3;
936        int val0, val1, val2, val3, t;
937
938        t = (int)(x[i].i >> 52);
939        if( (t & 2047) > 1023 + 10 )
940            x0 = t < 0 ? -exp_max_val : exp_max_val;
941
942        t = (int)(x[i+1].i >> 52);
943        if( (t & 2047) > 1023 + 10 )
944            x1 = t < 0 ? -exp_max_val : exp_max_val;
945
946        t = (int)(x[i+2].i >> 52);
947        if( (t & 2047) > 1023 + 10 )
948            x2 = t < 0 ? -exp_max_val : exp_max_val;
949
950        t = (int)(x[i+3].i >> 52);
951        if( (t & 2047) > 1023 + 10 )
952            x3 = t < 0 ? -exp_max_val : exp_max_val;
953
954        val0 = cvRound(x0);
955        val1 = cvRound(x1);
956        val2 = cvRound(x2);
957        val3 = cvRound(x3);
958
959        x0 = (x0 - val0)*exp_postscale;
960        x1 = (x1 - val1)*exp_postscale;
961        x2 = (x2 - val2)*exp_postscale;
962        x3 = (x3 - val3)*exp_postscale;
963
964        t = (val0 >> EXPTAB_SCALE) + 1023;
965        t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
966        buf[0].i.hi = t << 20;
967
968        t = (val1 >> EXPTAB_SCALE) + 1023;
969        t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
970        buf[1].i.hi = t << 20;
971
972        t = (val2 >> EXPTAB_SCALE) + 1023;
973        t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
974        buf[2].i.hi = t << 20;
975
976        t = (val3 >> EXPTAB_SCALE) + 1023;
977        t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
978        buf[3].i.hi = t << 20;
979
980        y0 = buf[0].d * icvExpTab[val0 & EXPTAB_MASK] * EXPPOLY( x0 );
981        y1 = buf[1].d * icvExpTab[val1 & EXPTAB_MASK] * EXPPOLY( x1 );
982
983        y[i] = y0;
984        y[i + 1] = y1;
985
986        y2 = buf[2].d * icvExpTab[val2 & EXPTAB_MASK] * EXPPOLY( x2 );
987        y3 = buf[3].d * icvExpTab[val3 & EXPTAB_MASK] * EXPPOLY( x3 );
988
989        y[i + 2] = y2;
990        y[i + 3] = y3;
991    }
992
993    for( ; i < n; i++ )
994    {
995        double x0 = x[i].f * exp_prescale;
996        int val0, t;
997
998        t = (int)(x[i].i >> 52);
999        if( (t & 2047) > 1023 + 10 )
1000            x0 = t < 0 ? -exp_max_val : exp_max_val;
1001
1002        val0 = cvRound(x0);
1003        t = (val0 >> EXPTAB_SCALE) + 1023;
1004        t = (t | ((t < 2047) - 1)) & (((t < 0) - 1) & 2047);
1005
1006        buf[0].i.hi = t << 20;
1007        x0 = (x0 - val0)*exp_postscale;
1008
1009        y[i] = buf[0].d * icvExpTab[val0 & EXPTAB_MASK] * EXPPOLY( x0 );
1010    }
1011
1012    return CV_OK;
1013}
1014
1015#undef EXPTAB_SCALE
1016#undef EXPTAB_MASK
1017#undef EXPPOLY_32F_A0
1018
1019CV_IMPL void cvExp( const CvArr* srcarr, CvArr* dstarr )
1020{
1021    CV_FUNCNAME( "cvExp" );
1022
1023    __BEGIN__;
1024
1025    CvMat srcstub, *src = (CvMat*)srcarr;
1026    CvMat dststub, *dst = (CvMat*)dstarr;
1027    int coi1 = 0, coi2 = 0, src_depth, dst_depth;
1028    double* buffer = 0;
1029    CvSize size;
1030    int x, y, dx = 0;
1031
1032    if( !CV_IS_MAT(src))
1033        CV_CALL( src = cvGetMat( src, &srcstub, &coi1 ));
1034
1035    if( !CV_IS_MAT(dst))
1036        CV_CALL( dst = cvGetMat( dst, &dststub, &coi2 ));
1037
1038    if( coi1 != 0 || coi2 != 0 )
1039        CV_ERROR( CV_BadCOI, "" );
1040
1041    src_depth = CV_MAT_DEPTH(src->type);
1042    dst_depth = CV_MAT_DEPTH(dst->type);
1043
1044    if( !CV_ARE_CNS_EQ( src, dst ) || src_depth < CV_32F || dst_depth < src_depth )
1045        CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
1046
1047    if( !CV_ARE_SIZES_EQ( src, dst ) )
1048        CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
1049
1050    size = cvGetMatSize(src);
1051    size.width *= CV_MAT_CN(src->type);
1052
1053    if( CV_IS_MAT_CONT( src->type & dst->type ))
1054    {
1055        size.width *= size.height;
1056        size.height = 1;
1057    }
1058
1059    if( !CV_ARE_DEPTHS_EQ( src, dst ))
1060    {
1061        dx = MIN( 1024, size.width );
1062        buffer = (double*)cvStackAlloc( dx*sizeof(buffer[0]) );
1063    }
1064
1065    for( y = 0; y < size.height; y++ )
1066    {
1067        uchar* src_data = src->data.ptr + src->step*y;
1068        uchar* dst_data = dst->data.ptr + dst->step*y;
1069
1070        if( src_depth == CV_64F )
1071        {
1072            icvExp_64f( (double*)src_data, (double*)dst_data, size.width );
1073        }
1074        else if( src_depth == dst_depth )
1075        {
1076            icvExp_32f( (float*)src_data, (float*)dst_data, size.width );
1077        }
1078        else
1079        {
1080            for( x = 0; x < size.width; x += dx )
1081            {
1082                int len = dx;
1083                if( x + len > size.width )
1084                    len = size.width - x;
1085                icvCvt_32f64f( (float*)src_data + x, buffer, len );
1086                icvExp_64f( buffer, (double*)dst_data + x, len );
1087            }
1088        }
1089    }
1090
1091    __END__;
1092}
1093
1094
1095/****************************************************************************************\
1096*                                          L O G                                         *
1097\****************************************************************************************/
1098
1099#define LOGTAB_SCALE    8
1100#define LOGTAB_MASK         ((1 << LOGTAB_SCALE) - 1)
1101#define LOGTAB_MASK2        ((1 << (20 - LOGTAB_SCALE)) - 1)
1102#define LOGTAB_MASK2_32F    ((1 << (23 - LOGTAB_SCALE)) - 1)
1103
1104static const double icvLogTab[] = {
11050.0000000000000000000000000000000000000000,    1.000000000000000000000000000000000000000,
1106.00389864041565732288852075271279318258166,    .9961089494163424124513618677042801556420,
1107.00778214044205494809292034119607706088573,    .9922480620155038759689922480620155038760,
1108.01165061721997527263705585198749759001657,    .9884169884169884169884169884169884169884,
1109.01550418653596525274396267235488267033361,    .9846153846153846153846153846153846153846,
1110.01934296284313093139406447562578250654042,    .9808429118773946360153256704980842911877,
1111.02316705928153437593630670221500622574241,    .9770992366412213740458015267175572519084,
1112.02697658769820207233514075539915211265906,    .9733840304182509505703422053231939163498,
1113.03077165866675368732785500469617545604706,    .9696969696969696969696969696969696969697,
1114.03455238150665972812758397481047722976656,    .9660377358490566037735849056603773584906,
1115.03831886430213659461285757856785494368522,    .9624060150375939849624060150375939849624,
1116.04207121392068705056921373852674150839447,    .9588014981273408239700374531835205992509,
1117.04580953603129420126371940114040626212953,    .9552238805970149253731343283582089552239,
1118.04953393512227662748292900118940451648088,    .9516728624535315985130111524163568773234,
1119.05324451451881227759255210685296333394944,    .9481481481481481481481481481481481481481,
1120.05694137640013842427411105973078520037234,    .9446494464944649446494464944649446494465,
1121.06062462181643483993820353816772694699466,    .9411764705882352941176470588235294117647,
1122.06429435070539725460836422143984236754475,    .9377289377289377289377289377289377289377,
1123.06795066190850773679699159401934593915938,    .9343065693430656934306569343065693430657,
1124.07159365318700880442825962290953611955044,    .9309090909090909090909090909090909090909,
1125.07522342123758751775142172846244648098944,    .9275362318840579710144927536231884057971,
1126.07884006170777602129362549021607264876369,    .9241877256317689530685920577617328519856,
1127.08244366921107458556772229485432035289706,    .9208633093525179856115107913669064748201,
1128.08603433734180314373940490213499288074675,    .9175627240143369175627240143369175627240,
1129.08961215868968712416897659522874164395031,    .9142857142857142857142857142857142857143,
1130.09317722485418328259854092721070628613231,    .9110320284697508896797153024911032028470,
1131.09672962645855109897752299730200320482256,    .9078014184397163120567375886524822695035,
1132.10026945316367513738597949668474029749630,    .9045936395759717314487632508833922261484,
1133.10379679368164355934833764649738441221420,    .9014084507042253521126760563380281690141,
1134.10731173578908805021914218968959175981580,    .8982456140350877192982456140350877192982,
1135.11081436634029011301105782649756292812530,    .8951048951048951048951048951048951048951,
1136.11430477128005862852422325204315711744130,    .8919860627177700348432055749128919860627,
1137.11778303565638344185817487641543266363440,    .8888888888888888888888888888888888888889,
1138.12124924363286967987640707633545389398930,    .8858131487889273356401384083044982698962,
1139.12470347850095722663787967121606925502420,    .8827586206896551724137931034482758620690,
1140.12814582269193003360996385708858724683530,    .8797250859106529209621993127147766323024,
1141.13157635778871926146571524895989568904040,    .8767123287671232876712328767123287671233,
1142.13499516453750481925766280255629681050780,    .8737201365187713310580204778156996587031,
1143.13840232285911913123754857224412262439730,    .8707482993197278911564625850340136054422,
1144.14179791186025733629172407290752744302150,    .8677966101694915254237288135593220338983,
1145.14518200984449788903951628071808954700830,    .8648648648648648648648648648648648648649,
1146.14855469432313711530824207329715136438610,    .8619528619528619528619528619528619528620,
1147.15191604202584196858794030049466527998450,    .8590604026845637583892617449664429530201,
1148.15526612891112392955683674244937719777230,    .8561872909698996655518394648829431438127,
1149.15860503017663857283636730244325008243330,    .8533333333333333333333333333333333333333,
1150.16193282026931324346641360989451641216880,    .8504983388704318936877076411960132890365,
1151.16524957289530714521497145597095368430010,    .8476821192052980132450331125827814569536,
1152.16855536102980664403538924034364754334090,    .8448844884488448844884488448844884488449,
1153.17185025692665920060697715143760433420540,    .8421052631578947368421052631578947368421,
1154.17513433212784912385018287750426679849630,    .8393442622950819672131147540983606557377,
1155.17840765747281828179637841458315961062910,    .8366013071895424836601307189542483660131,
1156.18167030310763465639212199675966985523700,    .8338762214983713355048859934853420195440,
1157.18492233849401198964024217730184318497780,    .8311688311688311688311688311688311688312,
1158.18816383241818296356839823602058459073300,    .8284789644012944983818770226537216828479,
1159.19139485299962943898322009772527962923050,    .8258064516129032258064516129032258064516,
1160.19461546769967164038916962454095482826240,    .8231511254019292604501607717041800643087,
1161.19782574332991986754137769821682013571260,    .8205128205128205128205128205128205128205,
1162.20102574606059073203390141770796617493040,    .8178913738019169329073482428115015974441,
1163.20421554142869088876999228432396193966280,    .8152866242038216560509554140127388535032,
1164.20739519434607056602715147164417430758480,    .8126984126984126984126984126984126984127,
1165.21056476910734961416338251183333341032260,    .8101265822784810126582278481012658227848,
1166.21372432939771812687723695489694364368910,    .8075709779179810725552050473186119873817,
1167.21687393830061435506806333251006435602900,    .8050314465408805031446540880503144654088,
1168.22001365830528207823135744547471404075630,    .8025078369905956112852664576802507836991,
1169.22314355131420973710199007200571941211830,    .8000000000000000000000000000000000000000,
1170.22626367865045338145790765338460914790630,    .7975077881619937694704049844236760124611,
1171.22937410106484582006380890106811420992010,    .7950310559006211180124223602484472049689,
1172.23247487874309405442296849741978803649550,    .7925696594427244582043343653250773993808,
1173.23556607131276688371634975283086532726890,    .7901234567901234567901234567901234567901,
1174.23864773785017498464178231643018079921600,    .7876923076923076923076923076923076923077,
1175.24171993688714515924331749374687206000090,    .7852760736196319018404907975460122699387,
1176.24478272641769091566565919038112042471760,    .7828746177370030581039755351681957186544,
1177.24783616390458124145723672882013488560910,    .7804878048780487804878048780487804878049,
1178.25088030628580937353433455427875742316250,    .7781155015197568389057750759878419452888,
1179.25391520998096339667426946107298135757450,    .7757575757575757575757575757575757575758,
1180.25694093089750041913887912414793390780680,    .7734138972809667673716012084592145015106,
1181.25995752443692604627401010475296061486000,    .7710843373493975903614457831325301204819,
1182.26296504550088134477547896494797896593800,    .7687687687687687687687687687687687687688,
1183.26596354849713793599974565040611196309330,    .7664670658682634730538922155688622754491,
1184.26895308734550393836570947314612567424780,    .7641791044776119402985074626865671641791,
1185.27193371548364175804834985683555714786050,    .7619047619047619047619047619047619047619,
1186.27490548587279922676529508862586226314300,    .7596439169139465875370919881305637982196,
1187.27786845100345625159121709657483734190480,    .7573964497041420118343195266272189349112,
1188.28082266290088775395616949026589281857030,    .7551622418879056047197640117994100294985,
1189.28376817313064456316240580235898960381750,    .7529411764705882352941176470588235294118,
1190.28670503280395426282112225635501090437180,    .7507331378299120234604105571847507331378,
1191.28963329258304265634293983566749375313530,    .7485380116959064327485380116959064327485,
1192.29255300268637740579436012922087684273730,    .7463556851311953352769679300291545189504,
1193.29546421289383584252163927885703742504130,    .7441860465116279069767441860465116279070,
1194.29836697255179722709783618483925238251680,    .7420289855072463768115942028985507246377,
1195.30126133057816173455023545102449133992200,    .7398843930635838150289017341040462427746,
1196.30414733546729666446850615102448500692850,    .7377521613832853025936599423631123919308,
1197.30702503529491181888388950937951449304830,    .7356321839080459770114942528735632183908,
1198.30989447772286465854207904158101882785550,    .7335243553008595988538681948424068767908,
1199.31275571000389684739317885942000430077330,    .7314285714285714285714285714285714285714,
1200.31560877898630329552176476681779604405180,    .7293447293447293447293447293447293447293,
1201.31845373111853458869546784626436419785030,    .7272727272727272727272727272727272727273,
1202.32129061245373424782201254856772720813750,    .7252124645892351274787535410764872521246,
1203.32411946865421192853773391107097268104550,    .7231638418079096045197740112994350282486,
1204.32694034499585328257253991068864706903700,    .7211267605633802816901408450704225352113,
1205.32975328637246797969240219572384376078850,    .7191011235955056179775280898876404494382,
1206.33255833730007655635318997155991382896900,    .7170868347338935574229691876750700280112,
1207.33535554192113781191153520921943709254280,    .7150837988826815642458100558659217877095,
1208.33814494400871636381467055798566434532400,    .7130919220055710306406685236768802228412,
1209.34092658697059319283795275623560883104800,    .7111111111111111111111111111111111111111,
1210.34370051385331840121395430287520866841080,    .7091412742382271468144044321329639889197,
1211.34646676734620857063262633346312213689100,    .7071823204419889502762430939226519337017,
1212.34922538978528827602332285096053965389730,    .7052341597796143250688705234159779614325,
1213.35197642315717814209818925519357435405250,    .7032967032967032967032967032967032967033,
1214.35471990910292899856770532096561510115850,    .7013698630136986301369863013698630136986,
1215.35745588892180374385176833129662554711100,    .6994535519125683060109289617486338797814,
1216.36018440357500774995358483465679455548530,    .6975476839237057220708446866485013623978,
1217.36290549368936841911903457003063522279280,    .6956521739130434782608695652173913043478,
1218.36561919956096466943762379742111079394830,    .6937669376693766937669376693766937669377,
1219.36832556115870762614150635272380895912650,    .6918918918918918918918918918918918918919,
1220.37102461812787262962487488948681857436900,    .6900269541778975741239892183288409703504,
1221.37371640979358405898480555151763837784530,    .6881720430107526881720430107526881720430,
1222.37640097516425302659470730759494472295050,    .6863270777479892761394101876675603217158,
1223.37907835293496944251145919224654790014030,    .6844919786096256684491978609625668449198,
1224.38174858149084833769393299007788300514230,    .6826666666666666666666666666666666666667,
1225.38441169891033200034513583887019194662580,    .6808510638297872340425531914893617021277,
1226.38706774296844825844488013899535872042180,    .6790450928381962864721485411140583554377,
1227.38971675114002518602873692543653305619950,    .6772486772486772486772486772486772486772,
1228.39235876060286384303665840889152605086580,    .6754617414248021108179419525065963060686,
1229.39499380824086893770896722344332374632350,    .6736842105263157894736842105263157894737,
1230.39762193064713846624158577469643205404280,    .6719160104986876640419947506561679790026,
1231.40024316412701266276741307592601515352730,    .6701570680628272251308900523560209424084,
1232.40285754470108348090917615991202183067800,    .6684073107049608355091383812010443864230,
1233.40546510810816432934799991016916465014230,    .6666666666666666666666666666666666666667,
1234.40806588980822172674223224930756259709600,    .6649350649350649350649350649350649350649,
1235.41065992498526837639616360320360399782650,    .6632124352331606217616580310880829015544,
1236.41324724855021932601317757871584035456180,    .6614987080103359173126614987080103359173,
1237.41582789514371093497757669865677598863850,    .6597938144329896907216494845360824742268,
1238.41840189913888381489925905043492093682300,    .6580976863753213367609254498714652956298,
1239.42096929464412963239894338585145305842150,    .6564102564102564102564102564102564102564,
1240.42353011550580327293502591601281892508280,    .6547314578005115089514066496163682864450,
1241.42608439531090003260516141381231136620050,    .6530612244897959183673469387755102040816,
1242.42863216738969872610098832410585600882780,    .6513994910941475826972010178117048346056,
1243.43117346481837132143866142541810404509300,    .6497461928934010152284263959390862944162,
1244.43370832042155937902094819946796633303180,    .6481012658227848101265822784810126582278,
1245.43623676677491801667585491486534010618930,    .6464646464646464646464646464646464646465,
1246.43875883620762790027214350629947148263450,    .6448362720403022670025188916876574307305,
1247.44127456080487520440058801796112675219780,    .6432160804020100502512562814070351758794,
1248.44378397241030093089975139264424797147500,    .6416040100250626566416040100250626566416,
1249.44628710262841947420398014401143882423650,    .6400000000000000000000000000000000000000,
1250.44878398282700665555822183705458883196130,    .6384039900249376558603491271820448877805,
1251.45127464413945855836729492693848442286250,    .6368159203980099502487562189054726368159,
1252.45375911746712049854579618113348260521900,    .6352357320099255583126550868486352357320,
1253.45623743348158757315857769754074979573500,    .6336633663366336633663366336633663366337,
1254.45870962262697662081833982483658473938700,    .6320987654320987654320987654320987654321,
1255.46117571512217014895185229761409573256980,    .6305418719211822660098522167487684729064,
1256.46363574096303250549055974261136725544930,    .6289926289926289926289926289926289926290,
1257.46608972992459918316399125615134835243230,    .6274509803921568627450980392156862745098,
1258.46853771156323925639597405279346276074650,    .6259168704156479217603911980440097799511,
1259.47097971521879100631480241645476780831830,    .6243902439024390243902439024390243902439,
1260.47341577001667212165614273544633761048330,    .6228710462287104622871046228710462287105,
1261.47584590486996386493601107758877333253630,    .6213592233009708737864077669902912621359,
1262.47827014848147025860569669930555392056700,    .6198547215496368038740920096852300242131,
1263.48068852934575190261057286988943815231330,    .6183574879227053140096618357487922705314,
1264.48310107575113581113157579238759353756900,    .6168674698795180722891566265060240963855,
1265.48550781578170076890899053978500887751580,    .6153846153846153846153846153846153846154,
1266.48790877731923892879351001283794175833480,    .6139088729016786570743405275779376498801,
1267.49030398804519381705802061333088204264650,    .6124401913875598086124401913875598086124,
1268.49269347544257524607047571407747454941280,    .6109785202863961813842482100238663484487,
1269.49507726679785146739476431321236304938800,    .6095238095238095238095238095238095238095,
1270.49745538920281889838648226032091770321130,    .6080760095011876484560570071258907363420,
1271.49982786955644931126130359189119189977650,    .6066350710900473933649289099526066350711,
1272.50219473456671548383667413872899487614650,    .6052009456264775413711583924349881796690,
1273.50455601075239520092452494282042607665050,    .6037735849056603773584905660377358490566,
1274.50691172444485432801997148999362252652650,    .6023529411764705882352941176470588235294,
1275.50926190178980790257412536448100581765150,    .6009389671361502347417840375586854460094,
1276.51160656874906207391973111953120678663250,    .5995316159250585480093676814988290398126,
1277.51394575110223428282552049495279788970950,    .5981308411214953271028037383177570093458,
1278.51627947444845445623684554448118433356300,    .5967365967365967365967365967365967365967,
1279.51860776420804555186805373523384332656850,    .5953488372093023255813953488372093023256,
1280.52093064562418522900344441950437612831600,    .5939675174013921113689095127610208816705,
1281.52324814376454775732838697877014055848100,    .5925925925925925925925925925925925925926,
1282.52556028352292727401362526507000438869000,    .5912240184757505773672055427251732101617,
1283.52786708962084227803046587723656557500350,    .5898617511520737327188940092165898617512,
1284.53016858660912158374145519701414741575700,    .5885057471264367816091954022988505747126,
1285.53246479886947173376654518506256863474850,    .5871559633027522935779816513761467889908,
1286.53475575061602764748158733709715306758900,    .5858123569794050343249427917620137299771,
1287.53704146589688361856929077475797384977350,    .5844748858447488584474885844748858447489,
1288.53932196859560876944783558428753167390800,    .5831435079726651480637813211845102505695,
1289.54159728243274429804188230264117009937750,    .5818181818181818181818181818181818181818,
1290.54386743096728351609669971367111429572100,    .5804988662131519274376417233560090702948,
1291.54613243759813556721383065450936555862450,    .5791855203619909502262443438914027149321,
1292.54839232556557315767520321969641372561450,    .5778781038374717832957110609480812641084,
1293.55064711795266219063194057525834068655950,    .5765765765765765765765765765765765765766,
1294.55289683768667763352766542084282264113450,    .5752808988764044943820224719101123595506,
1295.55514150754050151093110798683483153581600,    .5739910313901345291479820627802690582960,
1296.55738115013400635344709144192165695130850,    .5727069351230425055928411633109619686801,
1297.55961578793542265941596269840374588966350,    .5714285714285714285714285714285714285714,
1298.56184544326269181269140062795486301183700,    .5701559020044543429844097995545657015590,
1299.56407013828480290218436721261241473257550,    .5688888888888888888888888888888888888889,
1300.56628989502311577464155334382667206227800,    .5676274944567627494456762749445676274945,
1301.56850473535266865532378233183408156037350,    .5663716814159292035398230088495575221239,
1302.57071468100347144680739575051120482385150,    .5651214128035320088300220750551876379691,
1303.57291975356178548306473885531886480748650,    .5638766519823788546255506607929515418502,
1304.57511997447138785144460371157038025558000,    .5626373626373626373626373626373626373626,
1305.57731536503482350219940144597785547375700,    .5614035087719298245614035087719298245614,
1306.57950594641464214795689713355386629700650,    .5601750547045951859956236323851203501094,
1307.58169173963462239562716149521293118596100,    .5589519650655021834061135371179039301310,
1308.58387276558098266665552955601015128195300,    .5577342047930283224400871459694989106754,
1309.58604904500357812846544902640744112432000,    .5565217391304347826086956521739130434783,
1310.58822059851708596855957011939608491957200,    .5553145336225596529284164859002169197397,
1311.59038744660217634674381770309992134571100,    .5541125541125541125541125541125541125541,
1312.59254960960667157898740242671919986605650,    .5529157667386609071274298056155507559395,
1313.59470710774669277576265358220553025603300,    .5517241379310344827586206896551724137931,
1314.59685996110779382384237123915227130055450,    .5505376344086021505376344086021505376344,
1315.59900818964608337768851242799428291618800,    .5493562231759656652360515021459227467811,
1316.60115181318933474940990890900138765573500,    .5481798715203426124197002141327623126338,
1317.60329085143808425240052883964381180703650,    .5470085470085470085470085470085470085470,
1318.60542532396671688843525771517306566238400,    .5458422174840085287846481876332622601279,
1319.60755525022454170969155029524699784815300,    .5446808510638297872340425531914893617021,
1320.60968064953685519036241657886421307921400,    .5435244161358811040339702760084925690021,
1321.61180154110599282990534675263916142284850,    .5423728813559322033898305084745762711864,
1322.61391794401237043121710712512140162289150,    .5412262156448202959830866807610993657505,
1323.61602987721551394351138242200249806046500,    .5400843881856540084388185654008438818565,
1324.61813735955507864705538167982012964785100,    .5389473684210526315789473684210526315789,
1325.62024040975185745772080281312810257077200,    .5378151260504201680672268907563025210084,
1326.62233904640877868441606324267922900617100,    .5366876310272536687631027253668763102725,
1327.62443328801189346144440150965237990021700,    .5355648535564853556485355648535564853556,
1328.62652315293135274476554741340805776417250,    .5344467640918580375782881002087682672234,
1329.62860865942237409420556559780379757285100,    .5333333333333333333333333333333333333333,
1330.63068982562619868570408243613201193511500,    .5322245322245322245322245322245322245322,
1331.63276666957103777644277897707070223987100,    .5311203319502074688796680497925311203320,
1332.63483920917301017716738442686619237065300,    .5300207039337474120082815734989648033126,
1333.63690746223706917739093569252872839570050,    .5289256198347107438016528925619834710744,
1334.63897144645792069983514238629140891134750,    .5278350515463917525773195876288659793814,
1335.64103117942093124081992527862894348800200,    .5267489711934156378600823045267489711934,
1336.64308667860302726193566513757104985415950,    .5256673511293634496919917864476386036961,
1337.64513796137358470073053240412264131009600,    .5245901639344262295081967213114754098361,
1338.64718504499530948859131740391603671014300,    .5235173824130879345603271983640081799591,
1339.64922794662510974195157587018911726772800,    .5224489795918367346938775510204081632653,
1340.65126668331495807251485530287027359008800,    .5213849287169042769857433808553971486762,
1341.65330127201274557080523663898929953575150,    .5203252032520325203252032520325203252033,
1342.65533172956312757406749369692988693714150,    .5192697768762677484787018255578093306288,
1343.65735807270835999727154330685152672231200,    .5182186234817813765182186234817813765182,
1344.65938031808912778153342060249997302889800,    .5171717171717171717171717171717171717172,
1345.66139848224536490484126716182800009846700,    .5161290322580645161290322580645161290323,
1346.66341258161706617713093692145776003599150,    .5150905432595573440643863179074446680080,
1347.66542263254509037562201001492212526500250,    .5140562248995983935742971887550200803213,
1348.66742865127195616370414654738851822912700,    .5130260521042084168336673346693386773547,
1349.66943065394262923906154583164607174694550,    .5120000000000000000000000000000000000000,
1350.67142865660530226534774556057527661323550,    .5109780439121756487025948103792415169661,
1351.67342267521216669923234121597488410770900,    .5099601593625498007968127490039840637450,
1352.67541272562017662384192817626171745359900,    .5089463220675944333996023856858846918489,
1353.67739882359180603188519853574689477682100,    .5079365079365079365079365079365079365079,
1354.67938098479579733801614338517538271844400,    .5069306930693069306930693069306930693069,
1355.68135922480790300781450241629499942064300,    .5059288537549407114624505928853754940711,
1356.68333355911162063645036823800182901322850,    .5049309664694280078895463510848126232742,
1357.68530400309891936760919861626462079584600,    .5039370078740157480314960629921259842520,
1358.68727057207096020619019327568821609020250,    .5029469548133595284872298624754420432220,
1359.68923328123880889251040571252815425395950,    .5019607843137254901960784313725490196078,
1360.69314718055994530941723212145818, 5.0e-01,
1361};
1362
1363
1364
1365#define LOGTAB_TRANSLATE(x,h) (((x) - 1.)*icvLogTab[(h)+1])
1366static const double ln_2 = 0.69314718055994530941723212145818;
1367
1368IPCVAPI_IMPL( CvStatus, icvLog_32f, ( const float *_x, float *y, int n ), (_x, y, n) )
1369{
1370    static const double shift[] = { 0, -1./512 };
1371    static const double
1372        A0 = 0.3333333333333333333333333,
1373        A1 = -0.5,
1374        A2 = 1;
1375
1376    #undef LOGPOLY
1377    #define LOGPOLY(x,k) ((x)+=shift[k],((A0*(x) + A1)*(x) + A2)*(x))
1378
1379    int i = 0;
1380    union
1381    {
1382        int i;
1383        float f;
1384    }
1385    buf[4];
1386
1387    const int* x = (const int*)_x;
1388
1389    if( !x || !y )
1390        return CV_NULLPTR_ERR;
1391    if( n <= 0 )
1392        return CV_BADSIZE_ERR;
1393
1394    for( i = 0; i <= n - 4; i += 4 )
1395    {
1396        double x0, x1, x2, x3;
1397        double y0, y1, y2, y3;
1398        int h0, h1, h2, h3;
1399
1400        h0 = x[i];
1401        h1 = x[i+1];
1402        buf[0].i = (h0 & LOGTAB_MASK2_32F) | (127 << 23);
1403        buf[1].i = (h1 & LOGTAB_MASK2_32F) | (127 << 23);
1404
1405        y0 = (((h0 >> 23) & 0xff) - 127) * ln_2;
1406        y1 = (((h1 >> 23) & 0xff) - 127) * ln_2;
1407
1408        h0 = (h0 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1409        h1 = (h1 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1410
1411        y0 += icvLogTab[h0];
1412        y1 += icvLogTab[h1];
1413
1414        h2 = x[i+2];
1415        h3 = x[i+3];
1416
1417        x0 = LOGTAB_TRANSLATE( buf[0].f, h0 );
1418        x1 = LOGTAB_TRANSLATE( buf[1].f, h1 );
1419
1420        buf[2].i = (h2 & LOGTAB_MASK2_32F) | (127 << 23);
1421        buf[3].i = (h3 & LOGTAB_MASK2_32F) | (127 << 23);
1422
1423        y2 = (((h2 >> 23) & 0xff) - 127) * ln_2;
1424        y3 = (((h3 >> 23) & 0xff) - 127) * ln_2;
1425
1426        h2 = (h2 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1427        h3 = (h3 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1428
1429        y2 += icvLogTab[h2];
1430        y3 += icvLogTab[h3];
1431
1432        x2 = LOGTAB_TRANSLATE( buf[2].f, h2 );
1433        x3 = LOGTAB_TRANSLATE( buf[3].f, h3 );
1434
1435        y0 += LOGPOLY( x0, h0 == 510 );
1436        y1 += LOGPOLY( x1, h1 == 510 );
1437
1438        y[i] = (float) y0;
1439        y[i + 1] = (float) y1;
1440
1441        y2 += LOGPOLY( x2, h2 == 510 );
1442        y3 += LOGPOLY( x3, h3 == 510 );
1443
1444        y[i + 2] = (float) y2;
1445        y[i + 3] = (float) y3;
1446    }
1447
1448    for( ; i < n; i++ )
1449    {
1450        int h0 = x[i];
1451        double x0, y0;
1452
1453        y0 = (((h0 >> 23) & 0xff) - 127) * ln_2;
1454
1455        buf[0].i = (h0 & LOGTAB_MASK2_32F) | (127 << 23);
1456        h0 = (h0 >> (23 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1457
1458        y0 += icvLogTab[h0];
1459        x0 = LOGTAB_TRANSLATE( buf[0].f, h0 );
1460        y0 += LOGPOLY( x0, h0 == 510 );
1461
1462        y[i] = (float)y0;
1463    }
1464
1465    return CV_OK;
1466}
1467
1468
1469IPCVAPI_IMPL( CvStatus, icvLog_64f, ( const double *x, double *y, int n ), (x, y, n) )
1470{
1471    static const double shift[] = { 0, -1./512 };
1472    static const double
1473        A0 = -.1666666666666666666666666666666666666666,
1474        A1 = +0.2,
1475        A2 = -0.25,
1476        A3 = +0.3333333333333333333333333333333333333333,
1477        A4 = -0.5,
1478        A5 = +1.0;
1479
1480    #undef LOGPOLY
1481    #define LOGPOLY(x,k) ((x)+=shift[k], (xq) = (x)*(x),\
1482        ((A0*(xq) + A2)*(xq) + A4)*(xq) + ((A1*(xq) + A3)*(xq) + A5)*(x))
1483
1484    int i = 0;
1485    DBLINT buf[4];
1486    DBLINT *X = (DBLINT *) x;
1487
1488    if( !x || !y )
1489        return CV_NULLPTR_ERR;
1490    if( n <= 0 )
1491        return CV_BADSIZE_ERR;
1492
1493    for( ; i <= n - 4; i += 4 )
1494    {
1495        double xq;
1496        double x0, x1, x2, x3;
1497        double y0, y1, y2, y3;
1498        int h0, h1, h2, h3;
1499
1500        h0 = X[i].i.lo;
1501        h1 = X[i + 1].i.lo;
1502        buf[0].i.lo = h0;
1503        buf[1].i.lo = h1;
1504
1505        h0 = X[i].i.hi;
1506        h1 = X[i + 1].i.hi;
1507        buf[0].i.hi = (h0 & LOGTAB_MASK2) | (1023 << 20);
1508        buf[1].i.hi = (h1 & LOGTAB_MASK2) | (1023 << 20);
1509
1510        y0 = (((h0 >> 20) & 0x7ff) - 1023) * ln_2;
1511        y1 = (((h1 >> 20) & 0x7ff) - 1023) * ln_2;
1512
1513        h2 = X[i + 2].i.lo;
1514        h3 = X[i + 3].i.lo;
1515        buf[2].i.lo = h2;
1516        buf[3].i.lo = h3;
1517
1518        h0 = (h0 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1519        h1 = (h1 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1520
1521        y0 += icvLogTab[h0];
1522        y1 += icvLogTab[h1];
1523
1524        h2 = X[i + 2].i.hi;
1525        h3 = X[i + 3].i.hi;
1526
1527        x0 = LOGTAB_TRANSLATE( buf[0].d, h0 );
1528        x1 = LOGTAB_TRANSLATE( buf[1].d, h1 );
1529
1530        buf[2].i.hi = (h2 & LOGTAB_MASK2) | (1023 << 20);
1531        buf[3].i.hi = (h3 & LOGTAB_MASK2) | (1023 << 20);
1532
1533        y2 = (((h2 >> 20) & 0x7ff) - 1023) * ln_2;
1534        y3 = (((h3 >> 20) & 0x7ff) - 1023) * ln_2;
1535
1536        h2 = (h2 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1537        h3 = (h3 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1538
1539        y2 += icvLogTab[h2];
1540        y3 += icvLogTab[h3];
1541
1542        x2 = LOGTAB_TRANSLATE( buf[2].d, h2 );
1543        x3 = LOGTAB_TRANSLATE( buf[3].d, h3 );
1544
1545        y0 += LOGPOLY( x0, h0 == 510 );
1546        y1 += LOGPOLY( x1, h1 == 510 );
1547
1548        y[i] = y0;
1549        y[i + 1] = y1;
1550
1551        y2 += LOGPOLY( x2, h2 == 510 );
1552        y3 += LOGPOLY( x3, h3 == 510 );
1553
1554        y[i + 2] = y2;
1555        y[i + 3] = y3;
1556    }
1557
1558    for( ; i < n; i++ )
1559    {
1560        int h0 = X[i].i.hi;
1561        double xq;
1562        double x0, y0 = (((h0 >> 20) & 0x7ff) - 1023) * ln_2;
1563
1564        buf[0].i.hi = (h0 & LOGTAB_MASK2) | (1023 << 20);
1565        buf[0].i.lo = X[i].i.lo;
1566        h0 = (h0 >> (20 - LOGTAB_SCALE - 1)) & LOGTAB_MASK * 2;
1567
1568        y0 += icvLogTab[h0];
1569        x0 = LOGTAB_TRANSLATE( buf[0].d, h0 );
1570        y0 += LOGPOLY( x0, h0 == 510 );
1571        y[i] = y0;
1572    }
1573
1574    return CV_OK;
1575}
1576
1577
1578CV_IMPL void cvLog( const CvArr* srcarr, CvArr* dstarr )
1579{
1580    CV_FUNCNAME( "cvLog" );
1581
1582    __BEGIN__;
1583
1584    CvMat srcstub, *src = (CvMat*)srcarr;
1585    CvMat dststub, *dst = (CvMat*)dstarr;
1586    int coi1 = 0, coi2 = 0, src_depth, dst_depth;
1587    double* buffer = 0;
1588    CvSize size;
1589    int x, y, dx = 0;
1590
1591    if( !CV_IS_MAT(src))
1592        CV_CALL( src = cvGetMat( src, &srcstub, &coi1 ));
1593
1594    if( !CV_IS_MAT(dst))
1595        CV_CALL( dst = cvGetMat( dst, &dststub, &coi2 ));
1596
1597    if( coi1 != 0 || coi2 != 0 )
1598        CV_ERROR( CV_BadCOI, "" );
1599
1600    src_depth = CV_MAT_DEPTH(src->type);
1601    dst_depth = CV_MAT_DEPTH(dst->type);
1602
1603    if( !CV_ARE_CNS_EQ( src, dst ) || dst_depth < CV_32F || src_depth < dst_depth )
1604        CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
1605
1606    if( !CV_ARE_SIZES_EQ( src, dst ) )
1607        CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
1608
1609    size = cvGetMatSize(src);
1610    size.width *= CV_MAT_CN(src->type);
1611
1612    if( CV_IS_MAT_CONT( src->type & dst->type ))
1613    {
1614        size.width *= size.height;
1615        size.height = 1;
1616    }
1617
1618    if( !CV_ARE_DEPTHS_EQ( src, dst ))
1619    {
1620        dx = MIN( 1024, size.width );
1621        buffer = (double*)cvStackAlloc( dx*sizeof(buffer[0]) );
1622    }
1623
1624    for( y = 0; y < size.height; y++ )
1625    {
1626        uchar* src_data = src->data.ptr + src->step*y;
1627        uchar* dst_data = dst->data.ptr + dst->step*y;
1628
1629        if( dst_depth == CV_64F )
1630        {
1631            icvLog_64f( (double*)src_data, (double*)dst_data, size.width );
1632        }
1633        else if( src_depth == dst_depth )
1634        {
1635            icvLog_32f( (float*)src_data, (float*)dst_data, size.width );
1636        }
1637        else
1638        {
1639            for( x = 0; x < size.width; x += dx )
1640            {
1641                int len = dx;
1642                if( x + len > size.width )
1643                    len = size.width - x;
1644                icvLog_64f( (double*)src_data + x, buffer, len );
1645                icvCvt_64f32f( buffer, (float*)dst_data + x, len );
1646            }
1647        }
1648    }
1649
1650    __END__;
1651}
1652
1653
1654/****************************************************************************************\
1655*                                    P O W E R                                           *
1656\****************************************************************************************/
1657
1658#define ICV_DEF_IPOW_OP( flavor, arrtype, worktype, cast_macro )                    \
1659static CvStatus CV_STDCALL                                                          \
1660icvIPow_##flavor( const arrtype* src, arrtype* dst, int len, int power )            \
1661{                                                                                   \
1662    int i;                                                                          \
1663                                                                                    \
1664    for( i = 0; i < len; i++ )                                                      \
1665    {                                                                               \
1666        worktype a = 1, b = src[i];                                                 \
1667        int p = power;                                                              \
1668        while( p > 1 )                                                              \
1669        {                                                                           \
1670            if( p & 1 )                                                             \
1671                a *= b;                                                             \
1672            b *= b;                                                                 \
1673            p >>= 1;                                                                \
1674        }                                                                           \
1675                                                                                    \
1676        a *= b;                                                                     \
1677        dst[i] = cast_macro(a);                                                     \
1678    }                                                                               \
1679                                                                                    \
1680    return CV_OK;                                                                   \
1681}
1682
1683
1684ICV_DEF_IPOW_OP( 8u, uchar, int, CV_CAST_8U )
1685ICV_DEF_IPOW_OP( 16u, ushort, int, CV_CAST_16U )
1686ICV_DEF_IPOW_OP( 16s, short, int, CV_CAST_16S )
1687ICV_DEF_IPOW_OP( 32s, int, int, CV_CAST_32S )
1688ICV_DEF_IPOW_OP( 32f, float, double, CV_CAST_32F )
1689ICV_DEF_IPOW_OP( 64f, double, double, CV_CAST_64F )
1690
1691#define icvIPow_8s 0
1692
1693CV_DEF_INIT_FUNC_TAB_1D( IPow )
1694
1695typedef CvStatus (CV_STDCALL * CvIPowFunc)( const void* src, void* dst, int len, int power );
1696typedef CvStatus (CV_STDCALL * CvSqrtFunc)( const void* src, void* dst, int len );
1697
1698CV_IMPL void cvPow( const CvArr* srcarr, CvArr* dstarr, double power )
1699{
1700    static CvFuncTable ipow_tab;
1701    static int inittab = 0;
1702
1703    CV_FUNCNAME( "cvPow" );
1704
1705    __BEGIN__;
1706
1707    void* temp_buffer = 0;
1708    int block_size = 0;
1709    CvMat srcstub, *src = (CvMat*)srcarr;
1710    CvMat dststub, *dst = (CvMat*)dstarr;
1711    int coi1 = 0, coi2 = 0;
1712    int depth;
1713    CvSize size;
1714    int x, y;
1715    int ipower = cvRound( power );
1716    int is_ipower = 0;
1717
1718    if( !CV_IS_MAT(src))
1719        CV_CALL( src = cvGetMat( src, &srcstub, &coi1 ));
1720
1721    if( !CV_IS_MAT(dst))
1722        CV_CALL( dst = cvGetMat( dst, &dststub, &coi2 ));
1723
1724    if( coi1 != 0 || coi2 != 0 )
1725        CV_ERROR( CV_BadCOI, "" );
1726
1727    if( !CV_ARE_TYPES_EQ( src, dst ))
1728        CV_ERROR_FROM_CODE( CV_StsUnmatchedFormats );
1729
1730    if( !CV_ARE_SIZES_EQ( src, dst ) )
1731        CV_ERROR_FROM_CODE( CV_StsUnmatchedSizes );
1732
1733    depth = CV_MAT_DEPTH( src->type );
1734
1735    if( fabs(ipower - power) < DBL_EPSILON )
1736    {
1737        if( !inittab )
1738        {
1739            icvInitIPowTable( &ipow_tab );
1740            inittab = 1;
1741        }
1742
1743        if( ipower < 0 )
1744        {
1745            CV_CALL( cvDiv( 0, src, dst ));
1746
1747            if( ipower == -1 )
1748                EXIT;
1749            ipower = -ipower;
1750            src = dst;
1751        }
1752
1753        switch( ipower )
1754        {
1755        case 0:
1756            cvSet( dst, cvScalarAll(1));
1757            EXIT;
1758        case 1:
1759            cvCopy( src, dst );
1760            EXIT;
1761        case 2:
1762            cvMul( src, src, dst );
1763            EXIT;
1764        default:
1765            is_ipower = 1;
1766        }
1767    }
1768    else if( depth < CV_32F )
1769        CV_ERROR( CV_StsUnsupportedFormat,
1770        "Fractional or negative integer power factor can be used "
1771        "with floating-point types only");
1772
1773    size = cvGetMatSize(src);
1774    size.width *= CV_MAT_CN(src->type);
1775
1776    if( CV_IS_MAT_CONT( src->type & dst->type ))
1777    {
1778        size.width *= size.height;
1779        size.height = 1;
1780    }
1781
1782    if( is_ipower )
1783    {
1784        CvIPowFunc pow_func = (CvIPowFunc)ipow_tab.fn_2d[depth];
1785        if( !pow_func )
1786            CV_ERROR( CV_StsUnsupportedFormat, "The data type is not supported" );
1787
1788        for( y = 0; y < size.height; y++ )
1789        {
1790            uchar* src_data = src->data.ptr + src->step*y;
1791            uchar* dst_data = dst->data.ptr + dst->step*y;
1792
1793            pow_func( src_data, dst_data, size.width, ipower );
1794        }
1795    }
1796    else if( fabs(fabs(power) - 0.5) < DBL_EPSILON )
1797    {
1798        CvSqrtFunc sqrt_func = power < 0 ?
1799            (depth == CV_32F ? (CvSqrtFunc)icvInvSqrt_32f : (CvSqrtFunc)icvInvSqrt_64f) :
1800            (depth == CV_32F ? (CvSqrtFunc)icvSqrt_32f : (CvSqrtFunc)icvSqrt_64f);
1801
1802        for( y = 0; y < size.height; y++ )
1803        {
1804            uchar* src_data = src->data.ptr + src->step*y;
1805            uchar* dst_data = dst->data.ptr + dst->step*y;
1806
1807            sqrt_func( src_data, dst_data, size.width );
1808        }
1809    }
1810    else
1811    {
1812        block_size = MIN( size.width, ICV_MATH_BLOCK_SIZE );
1813        temp_buffer = cvStackAlloc( block_size*CV_ELEM_SIZE(depth) );
1814
1815        for( y = 0; y < size.height; y++ )
1816        {
1817            uchar* src_data = src->data.ptr + src->step*y;
1818            uchar* dst_data = dst->data.ptr + dst->step*y;
1819
1820            for( x = 0; x < size.width; x += block_size )
1821            {
1822                int len = MIN( size.width - x, block_size );
1823                if( depth == CV_32F )
1824                {
1825                    icvLog_32f( (float*)src_data + x, (float*)temp_buffer, len );
1826                    icvScale_32f( (float*)temp_buffer, (float*)temp_buffer, len, (float)power, 0 );
1827                    icvExp_32f( (float*)temp_buffer, (float*)dst_data + x, len );
1828                }
1829                else
1830                {
1831                    icvLog_64f( (double*)src_data + x, (double*)temp_buffer, len );
1832                    icvScale_64f( (double*)temp_buffer, (double*)temp_buffer, len, power, 0 );
1833                    icvExp_64f( (double*)temp_buffer, (double*)dst_data + x, len );
1834                }
1835            }
1836        }
1837    }
1838
1839    __END__;
1840}
1841
1842
1843/************************** CheckArray for NaN's, Inf's *********************************/
1844
1845IPCVAPI_IMPL( CvStatus, icvCheckArray_32f_C1R,
1846    ( const float* src, int srcstep, CvSize size, int flags, double min_val, double max_val ),
1847     (src, srcstep, size, flags, min_val, max_val) )
1848{
1849    Cv32suf a, b;
1850    int ia, ib;
1851    const int* isrc = (const int*)src;
1852
1853    if( !src )
1854        return CV_NULLPTR_ERR;
1855
1856    if( size.width <= 0 || size.height <= 0 )
1857        return CV_BADSIZE_ERR;
1858
1859    if( flags & CV_CHECK_RANGE )
1860    {
1861        a.f = (float)min_val;
1862        b.f = (float)max_val;
1863    }
1864    else
1865    {
1866        a.f = -FLT_MAX;
1867        b.f = FLT_MAX;
1868    }
1869
1870    ia = CV_TOGGLE_FLT(a.i);
1871    ib = CV_TOGGLE_FLT(b.i);
1872
1873    srcstep /= sizeof(isrc[0]);
1874    for( ; size.height--; isrc += srcstep )
1875    {
1876        int i;
1877        for( i = 0; i < size.width; i++ )
1878        {
1879            int val = isrc[i];
1880            val = CV_TOGGLE_FLT(val);
1881
1882            if( val < ia || val >= ib )
1883                return CV_BADRANGE_ERR;
1884        }
1885    }
1886
1887    return CV_OK;
1888}
1889
1890
1891IPCVAPI_IMPL( CvStatus,  icvCheckArray_64f_C1R,
1892    ( const double* src, int srcstep, CvSize size, int flags, double min_val, double max_val ),
1893    (src, srcstep, size, flags, min_val, max_val) )
1894{
1895    Cv64suf a, b;
1896    int64 ia, ib;
1897    const int64* isrc = (const int64*)src;
1898
1899    if( !src )
1900        return CV_NULLPTR_ERR;
1901
1902    if( size.width <= 0 || size.height <= 0 )
1903        return CV_BADSIZE_ERR;
1904
1905    if( flags & CV_CHECK_RANGE )
1906    {
1907        a.f = min_val;
1908        b.f = max_val;
1909    }
1910    else
1911    {
1912        a.f = -DBL_MAX;
1913        b.f = DBL_MAX;
1914    }
1915
1916    ia = CV_TOGGLE_DBL(a.i);
1917    ib = CV_TOGGLE_DBL(b.i);
1918
1919    srcstep /= sizeof(isrc[0]);
1920    for( ; size.height--; isrc += srcstep )
1921    {
1922        int i;
1923        for( i = 0; i < size.width; i++ )
1924        {
1925            int64 val = isrc[i];
1926            val = CV_TOGGLE_DBL(val);
1927
1928            if( val < ia || val >= ib )
1929                return CV_BADRANGE_ERR;
1930        }
1931    }
1932
1933    return CV_OK;
1934}
1935
1936
1937CV_IMPL  int  cvCheckArr( const CvArr* arr, int flags,
1938                          double minVal, double maxVal )
1939{
1940    int result = 0;
1941
1942    CV_FUNCNAME( "cvCheckArr" );
1943
1944    __BEGIN__;
1945
1946    if( arr )
1947    {
1948        CvStatus status = CV_OK;
1949        CvMat stub, *mat = (CvMat*)arr;
1950        int type;
1951        CvSize size;
1952
1953        if( !CV_IS_MAT( mat ))
1954            CV_CALL( mat = cvGetMat( mat, &stub, 0, 1 ));
1955
1956        type = CV_MAT_TYPE( mat->type );
1957        size = cvGetMatSize( mat );
1958
1959        size.width *= CV_MAT_CN( type );
1960
1961        if( CV_IS_MAT_CONT( mat->type ))
1962        {
1963            size.width *= size.height;
1964            size.height = 1;
1965        }
1966
1967        if( CV_MAT_DEPTH(type) == CV_32F )
1968        {
1969            status = icvCheckArray_32f_C1R( mat->data.fl, mat->step, size,
1970                                            flags, minVal, maxVal );
1971        }
1972        else if( CV_MAT_DEPTH(type) == CV_64F )
1973        {
1974            status = icvCheckArray_64f_C1R( mat->data.db, mat->step, size,
1975                                            flags, minVal, maxVal );
1976        }
1977        else
1978        {
1979            CV_ERROR( CV_StsUnsupportedFormat, "" );
1980        }
1981
1982        if( status < 0 )
1983        {
1984            if( status != CV_BADRANGE_ERR || !(flags & CV_CHECK_QUIET))
1985                CV_ERROR( CV_StsOutOfRange, "CheckArray failed" );
1986            EXIT;
1987        }
1988    }
1989
1990    result = 1;
1991
1992    __END__;
1993
1994    return result;
1995}
1996
1997
1998/* End of file. */
1999