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