CoCalc -- mathfuncs.cpp

GitHub Repository: Tetragramm/opencv
Path: blob/master/modules/core/src/mathfuncs.cpp
¹⁶³³⁷ views
1
/*M///////////////////////////////////////////////////////////////////////////////////////
2
//
3
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4
//
5
//  By downloading, copying, installing or using the software you agree to this license.
6
//  If you do not agree to this license, do not download, install,
7
//  copy or use the software.
8
//
9
//
10
//                           License Agreement
11
//                For Open Source Computer Vision Library
12
//
13
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
14
// Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
15
// Copyright (C) 2014-2015, Itseez Inc., all rights reserved.
16
// Third party copyrights are property of their respective owners.
17
//
18
// Redistribution and use in source and binary forms, with or without modification,
19
// are permitted provided that the following conditions are met:
20
//
21
//   * Redistribution's of source code must retain the above copyright notice,
22
//     this list of conditions and the following disclaimer.
23
//
24
//   * Redistribution's in binary form must reproduce the above copyright notice,
25
//     this list of conditions and the following disclaimer in the documentation
26
//     and/or other materials provided with the distribution.
27
//
28
//   * The name of the copyright holders may not be used to endorse or promote products
29
//     derived from this software without specific prior written permission.
30
//
31
// This software is provided by the copyright holders and contributors "as is" and
32
// any express or implied warranties, including, but not limited to, the implied
33
// warranties of merchantability and fitness for a particular purpose are disclaimed.
34
// In no event shall the Intel Corporation or contributors be liable for any direct,
35
// indirect, incidental, special, exemplary, or consequential damages
36
// (including, but not limited to, procurement of substitute goods or services;
37
// loss of use, data, or profits; or business interruption) however caused
38
// and on any theory of liability, whether in contract, strict liability,
39
// or tort (including negligence or otherwise) arising in any way out of
40
// the use of this software, even if advised of the possibility of such damage.
41
//
42
//M*/
43

44
#include "precomp.hpp"
45
#include "opencl_kernels_core.hpp"
46
#include <limits>
47
#include <iostream>
48
#include "mathfuncs.hpp"
49

50
namespace cv
51
{
52

53
typedef void (*MathFunc)(const void* src, void* dst, int len);
54

55
#ifdef HAVE_OPENCL
56

57
enum { OCL_OP_LOG=0, OCL_OP_EXP=1, OCL_OP_MAG=2, OCL_OP_PHASE_DEGREES=3, OCL_OP_PHASE_RADIANS=4 };
58

59
static const char* oclop2str[] = { "OP_LOG", "OP_EXP", "OP_MAG", "OP_PHASE_DEGREES", "OP_PHASE_RADIANS", 0 };
60

61
static bool ocl_math_op(InputArray _src1, InputArray _src2, OutputArray _dst, int oclop)
62
{
63
    int type = _src1.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
64
    int kercn = oclop == OCL_OP_PHASE_DEGREES ||
65
            oclop == OCL_OP_PHASE_RADIANS ? 1 : ocl::predictOptimalVectorWidth(_src1, _src2, _dst);
66

67
    const ocl::Device d = ocl::Device::getDefault();
68
    bool double_support = d.doubleFPConfig() > 0;
69
    if (!double_support && depth == CV_64F)
70
        return false;
71
    int rowsPerWI = d.isIntel() ? 4 : 1;
72

73
    ocl::Kernel k("KF", ocl::core::arithm_oclsrc,
74
                  format("-D %s -D %s -D dstT=%s -D DEPTH_dst=%d -D rowsPerWI=%d%s", _src2.empty() ? "UNARY_OP" : "BINARY_OP",
75
                         oclop2str[oclop], ocl::typeToStr(CV_MAKE_TYPE(depth, kercn)), depth, rowsPerWI,
76
                         double_support ? " -D DOUBLE_SUPPORT" : ""));
77
    if (k.empty())
78
        return false;
79

80
    UMat src1 = _src1.getUMat(), src2 = _src2.getUMat();
81
    _dst.create(src1.size(), type);
82
    UMat dst = _dst.getUMat();
83

84
    ocl::KernelArg src1arg = ocl::KernelArg::ReadOnlyNoSize(src1),
85
            src2arg = ocl::KernelArg::ReadOnlyNoSize(src2),
86
            dstarg = ocl::KernelArg::WriteOnly(dst, cn, kercn);
87

88
    if (src2.empty())
89
        k.args(src1arg, dstarg);
90
    else
91
        k.args(src1arg, src2arg, dstarg);
92

93
    size_t globalsize[] = { (size_t)src1.cols * cn / kercn, ((size_t)src1.rows + rowsPerWI - 1) / rowsPerWI };
94
    return k.run(2, globalsize, 0, false);
95
}
96

97
#endif
98

99
/* ************************************************************************** *\
100
   Fast cube root by Ken Turkowski
101
   (http://www.worldserver.com/turk/computergraphics/papers.html)
102
\* ************************************************************************** */
103
float  cubeRoot( float value )
104
{
105
    CV_INSTRUMENT_REGION();
106

107
    float fr;
108
    Cv32suf v, m;
109
    int ix, s;
110
    int ex, shx;
111

112
    v.f = value;
113
    ix = v.i & 0x7fffffff;
114
    s = v.i & 0x80000000;
115
    ex = (ix >> 23) - 127;
116
    shx = ex % 3;
117
    shx -= shx >= 0 ? 3 : 0;
118
    ex = (ex - shx) / 3; /* exponent of cube root */
119
    v.i = (ix & ((1<<23)-1)) | ((shx + 127)<<23);
120
    fr = v.f;
121

122
    /* 0.125 <= fr < 1.0 */
123
    /* Use quartic rational polynomial with error < 2^(-24) */
124
    fr = (float)(((((45.2548339756803022511987494 * fr +
125
    192.2798368355061050458134625) * fr +
126
    119.1654824285581628956914143) * fr +
127
    13.43250139086239872172837314) * fr +
128
    0.1636161226585754240958355063)/
129
    ((((14.80884093219134573786480845 * fr +
130
    151.9714051044435648658557668) * fr +
131
    168.5254414101568283957668343) * fr +
132
    33.9905941350215598754191872) * fr +
133
    1.0));
134

135
    /* fr *= 2^ex * sign */
136
    m.f = value;
137
    v.f = fr;
138
    v.i = (v.i + (ex << 23) + s) & (m.i*2 != 0 ? -1 : 0);
139
    return v.f;
140
}
141

142
/****************************************************************************************\
143
*                                  Cartezian -> Polar                                    *
144
\****************************************************************************************/
145

146
void magnitude( InputArray src1, InputArray src2, OutputArray dst )
147
{
148
    CV_INSTRUMENT_REGION();
149

150
    int type = src1.type(), depth = src1.depth(), cn = src1.channels();
151
    CV_Assert( src1.size() == src2.size() && type == src2.type() && (depth == CV_32F || depth == CV_64F));
152

153
    CV_OCL_RUN(dst.isUMat() && src1.dims() <= 2 && src2.dims() <= 2,
154
               ocl_math_op(src1, src2, dst, OCL_OP_MAG))
155

156
    Mat X = src1.getMat(), Y = src2.getMat();
157
    dst.create(X.dims, X.size, X.type());
158
    Mat Mag = dst.getMat();
159

160
    const Mat* arrays[] = {&X, &Y, &Mag, 0};
161
    uchar* ptrs[3] = {};
162
    NAryMatIterator it(arrays, ptrs);
163
    int len = (int)it.size*cn;
164

165
    for( size_t i = 0; i < it.nplanes; i++, ++it )
166
    {
167
        if( depth == CV_32F )
168
        {
169
            const float *x = (const float*)ptrs[0], *y = (const float*)ptrs[1];
170
            float *mag = (float*)ptrs[2];
171
            hal::magnitude32f( x, y, mag, len );
172
        }
173
        else
174
        {
175
            const double *x = (const double*)ptrs[0], *y = (const double*)ptrs[1];
176
            double *mag = (double*)ptrs[2];
177
            hal::magnitude64f( x, y, mag, len );
178
        }
179
    }
180
}
181

182
void phase( InputArray src1, InputArray src2, OutputArray dst, bool angleInDegrees )
183
{
184
    CV_INSTRUMENT_REGION();
185

186
    int type = src1.type(), depth = src1.depth(), cn = src1.channels();
187
    CV_Assert( src1.size() == src2.size() && type == src2.type() && (depth == CV_32F || depth == CV_64F));
188

189
    CV_OCL_RUN(dst.isUMat() && src1.dims() <= 2 && src2.dims() <= 2,
190
               ocl_math_op(src1, src2, dst, angleInDegrees ? OCL_OP_PHASE_DEGREES : OCL_OP_PHASE_RADIANS))
191

192
    Mat X = src1.getMat(), Y = src2.getMat();
193
    dst.create( X.dims, X.size, type );
194
    Mat Angle = dst.getMat();
195

196
    const Mat* arrays[] = {&X, &Y, &Angle, 0};
197
    uchar* ptrs[3] = {};
198
    NAryMatIterator it(arrays, ptrs);
199
    int j, total = (int)(it.size*cn), blockSize = total;
200
    size_t esz1 = X.elemSize1();
201
    for( size_t i = 0; i < it.nplanes; i++, ++it )
202
    {
203
        for( j = 0; j < total; j += blockSize )
204
        {
205
            int len = std::min(total - j, blockSize);
206
            if( depth == CV_32F )
207
            {
208
                const float *x = (const float*)ptrs[0], *y = (const float*)ptrs[1];
209
                float *angle = (float*)ptrs[2];
210
                hal::fastAtan32f( y, x, angle, len, angleInDegrees );
211
            }
212
            else
213
            {
214
                const double *x = (const double*)ptrs[0], *y = (const double*)ptrs[1];
215
                double *angle = (double*)ptrs[2];
216
                hal::fastAtan64f(y, x, angle, len, angleInDegrees);
217
            }
218
            ptrs[0] += len*esz1;
219
            ptrs[1] += len*esz1;
220
            ptrs[2] += len*esz1;
221
        }
222
    }
223
}
224

225
#ifdef HAVE_OPENCL
226

227
static bool ocl_cartToPolar( InputArray _src1, InputArray _src2,
228
                             OutputArray _dst1, OutputArray _dst2, bool angleInDegrees )
229
{
230
    const ocl::Device & d = ocl::Device::getDefault();
231
    int type = _src1.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type),
232
            rowsPerWI = d.isIntel() ? 4 : 1;
233
    bool doubleSupport = d.doubleFPConfig() > 0;
234

235
    if ( !(_src1.dims() <= 2 && _src2.dims() <= 2 &&
236
           (depth == CV_32F || depth == CV_64F) && type == _src2.type()) ||
237
         (depth == CV_64F && !doubleSupport) )
238
        return false;
239

240
    ocl::Kernel k("KF", ocl::core::arithm_oclsrc,
241
                  format("-D BINARY_OP -D dstT=%s -D DEPTH_dst=%d -D rowsPerWI=%d -D OP_CTP_%s%s",
242
                         ocl::typeToStr(CV_MAKE_TYPE(depth, 1)), depth,
243
                         rowsPerWI, angleInDegrees ? "AD" : "AR",
244
                         doubleSupport ? " -D DOUBLE_SUPPORT" : ""));
245
    if (k.empty())
246
        return false;
247

248
    UMat src1 = _src1.getUMat(), src2 = _src2.getUMat();
249
    Size size = src1.size();
250
    CV_Assert( size == src2.size() );
251

252
    _dst1.create(size, type);
253
    _dst2.create(size, type);
254
    UMat dst1 = _dst1.getUMat(), dst2 = _dst2.getUMat();
255

256
    k.args(ocl::KernelArg::ReadOnlyNoSize(src1),
257
           ocl::KernelArg::ReadOnlyNoSize(src2),
258
           ocl::KernelArg::WriteOnly(dst1, cn),
259
           ocl::KernelArg::WriteOnlyNoSize(dst2));
260

261
    size_t globalsize[2] = { (size_t)dst1.cols * cn, ((size_t)dst1.rows + rowsPerWI - 1) / rowsPerWI };
262
    return k.run(2, globalsize, NULL, false);
263
}
264

265
#endif
266

267
void cartToPolar( InputArray src1, InputArray src2,
268
                  OutputArray dst1, OutputArray dst2, bool angleInDegrees )
269
{
270
    CV_INSTRUMENT_REGION();
271

272
    CV_OCL_RUN(dst1.isUMat() && dst2.isUMat(),
273
            ocl_cartToPolar(src1, src2, dst1, dst2, angleInDegrees))
274

275
    Mat X = src1.getMat(), Y = src2.getMat();
276
    int type = X.type(), depth = X.depth(), cn = X.channels();
277
    CV_Assert( X.size == Y.size && type == Y.type() && (depth == CV_32F || depth == CV_64F));
278
    dst1.create( X.dims, X.size, type );
279
    dst2.create( X.dims, X.size, type );
280
    Mat Mag = dst1.getMat(), Angle = dst2.getMat();
281

282
    const Mat* arrays[] = {&X, &Y, &Mag, &Angle, 0};
283
    uchar* ptrs[4] = {};
284
    NAryMatIterator it(arrays, ptrs);
285
    int j, total = (int)(it.size*cn), blockSize = std::min(total, ((BLOCK_SIZE+cn-1)/cn)*cn);
286
    size_t esz1 = X.elemSize1();
287

288
    for( size_t i = 0; i < it.nplanes; i++, ++it )
289
    {
290
        for( j = 0; j < total; j += blockSize )
291
        {
292
            int len = std::min(total - j, blockSize);
293
            if( depth == CV_32F )
294
            {
295
                const float *x = (const float*)ptrs[0], *y = (const float*)ptrs[1];
296
                float *mag = (float*)ptrs[2], *angle = (float*)ptrs[3];
297
                hal::magnitude32f( x, y, mag, len );
298
                hal::fastAtan32f( y, x, angle, len, angleInDegrees );
299
            }
300
            else
301
            {
302
                const double *x = (const double*)ptrs[0], *y = (const double*)ptrs[1];
303
                double *angle = (double*)ptrs[3];
304
                hal::magnitude64f(x, y, (double*)ptrs[2], len);
305
                hal::fastAtan64f(y, x, angle, len, angleInDegrees);
306
            }
307
            ptrs[0] += len*esz1;
308
            ptrs[1] += len*esz1;
309
            ptrs[2] += len*esz1;
310
            ptrs[3] += len*esz1;
311
        }
312
    }
313
}
314

315

316
/****************************************************************************************\
317
*                                  Polar -> Cartezian                                    *
318
\****************************************************************************************/
319

320
static void SinCos_32f( const float *angle, float *sinval, float* cosval,
321
                        int len, int angle_in_degrees )
322
{
323
    const int N = 64;
324

325
    static const double sin_table[] =
326
    {
327
     0.00000000000000000000,     0.09801714032956060400,
328
     0.19509032201612825000,     0.29028467725446233000,
329
     0.38268343236508978000,     0.47139673682599764000,
330
     0.55557023301960218000,     0.63439328416364549000,
331
     0.70710678118654746000,     0.77301045336273699000,
332
     0.83146961230254524000,     0.88192126434835494000,
333
     0.92387953251128674000,     0.95694033573220894000,
334
     0.98078528040323043000,     0.99518472667219682000,
335
     1.00000000000000000000,     0.99518472667219693000,
336
     0.98078528040323043000,     0.95694033573220894000,
337
     0.92387953251128674000,     0.88192126434835505000,
338
     0.83146961230254546000,     0.77301045336273710000,
339
     0.70710678118654757000,     0.63439328416364549000,
340
     0.55557023301960218000,     0.47139673682599786000,
341
     0.38268343236508989000,     0.29028467725446239000,
342
     0.19509032201612861000,     0.09801714032956082600,
343
     0.00000000000000012246,    -0.09801714032956059000,
344
    -0.19509032201612836000,    -0.29028467725446211000,
345
    -0.38268343236508967000,    -0.47139673682599764000,
346
    -0.55557023301960196000,    -0.63439328416364527000,
347
    -0.70710678118654746000,    -0.77301045336273666000,
348
    -0.83146961230254524000,    -0.88192126434835494000,
349
    -0.92387953251128652000,    -0.95694033573220882000,
350
    -0.98078528040323032000,    -0.99518472667219693000,
351
    -1.00000000000000000000,    -0.99518472667219693000,
352
    -0.98078528040323043000,    -0.95694033573220894000,
353
    -0.92387953251128663000,    -0.88192126434835505000,
354
    -0.83146961230254546000,    -0.77301045336273688000,
355
    -0.70710678118654768000,    -0.63439328416364593000,
356
    -0.55557023301960218000,    -0.47139673682599792000,
357
    -0.38268343236509039000,    -0.29028467725446250000,
358
    -0.19509032201612872000,    -0.09801714032956050600,
359
    };
360

361
    static const double k2 = (2*CV_PI)/N;
362

363
    static const double sin_a0 = -0.166630293345647*k2*k2*k2;
364
    static const double sin_a2 = k2;
365

366
    static const double cos_a0 = -0.499818138450326*k2*k2;
367
    /*static const double cos_a2 =  1;*/
368

369
    double k1;
370
    int i = 0;
371

372
    if( !angle_in_degrees )
373
        k1 = N/(2*CV_PI);
374
    else
375
        k1 = N/360.;
376

377
#if CV_AVX2
378
    if (USE_AVX2)
379
    {
380
        __m128d v_k1 = _mm_set1_pd(k1);
381
        __m128d v_1 = _mm_set1_pd(1);
382
        __m128i v_N1 = _mm_set1_epi32(N - 1);
383
        __m128i v_N4 = _mm_set1_epi32(N >> 2);
384
        __m128d v_sin_a0 = _mm_set1_pd(sin_a0);
385
        __m128d v_sin_a2 = _mm_set1_pd(sin_a2);
386
        __m128d v_cos_a0 = _mm_set1_pd(cos_a0);
387

388
        for ( ; i <= len - 4; i += 4)
389
        {
390
            __m128 v_angle = _mm_loadu_ps(angle + i);
391

392
            // 0-1
393
            __m128d v_t = _mm_mul_pd(_mm_cvtps_pd(v_angle), v_k1);
394
            __m128i v_it = _mm_cvtpd_epi32(v_t);
395
            v_t = _mm_sub_pd(v_t, _mm_cvtepi32_pd(v_it));
396

397
            __m128i v_sin_idx = _mm_and_si128(v_it, v_N1);
398
            __m128i v_cos_idx = _mm_and_si128(_mm_sub_epi32(v_N4, v_sin_idx), v_N1);
399

400
            __m128d v_t2 = _mm_mul_pd(v_t, v_t);
401
            __m128d v_sin_b = _mm_mul_pd(_mm_add_pd(_mm_mul_pd(v_sin_a0, v_t2), v_sin_a2), v_t);
402
            __m128d v_cos_b = _mm_add_pd(_mm_mul_pd(v_cos_a0, v_t2), v_1);
403

404
            __m128d v_sin_a = _mm_i32gather_pd(sin_table, v_sin_idx, 8);
405
            __m128d v_cos_a = _mm_i32gather_pd(sin_table, v_cos_idx, 8);
406

407
            __m128d v_sin_val_0 = _mm_add_pd(_mm_mul_pd(v_sin_a, v_cos_b),
408
                                             _mm_mul_pd(v_cos_a, v_sin_b));
409
            __m128d v_cos_val_0 = _mm_sub_pd(_mm_mul_pd(v_cos_a, v_cos_b),
410
                                             _mm_mul_pd(v_sin_a, v_sin_b));
411

412
            // 2-3
413
            v_t = _mm_mul_pd(_mm_cvtps_pd(_mm_castsi128_ps(_mm_srli_si128(_mm_castps_si128(v_angle), 8))), v_k1);
414
            v_it = _mm_cvtpd_epi32(v_t);
415
            v_t = _mm_sub_pd(v_t, _mm_cvtepi32_pd(v_it));
416

417
            v_sin_idx = _mm_and_si128(v_it, v_N1);
418
            v_cos_idx = _mm_and_si128(_mm_sub_epi32(v_N4, v_sin_idx), v_N1);
419

420
            v_t2 = _mm_mul_pd(v_t, v_t);
421
            v_sin_b = _mm_mul_pd(_mm_add_pd(_mm_mul_pd(v_sin_a0, v_t2), v_sin_a2), v_t);
422
            v_cos_b = _mm_add_pd(_mm_mul_pd(v_cos_a0, v_t2), v_1);
423

424
            v_sin_a = _mm_i32gather_pd(sin_table, v_sin_idx, 8);
425
            v_cos_a = _mm_i32gather_pd(sin_table, v_cos_idx, 8);
426

427
            __m128d v_sin_val_1 = _mm_add_pd(_mm_mul_pd(v_sin_a, v_cos_b),
428
                                             _mm_mul_pd(v_cos_a, v_sin_b));
429
            __m128d v_cos_val_1 = _mm_sub_pd(_mm_mul_pd(v_cos_a, v_cos_b),
430
                                             _mm_mul_pd(v_sin_a, v_sin_b));
431

432
            _mm_storeu_ps(sinval + i, _mm_movelh_ps(_mm_cvtpd_ps(v_sin_val_0),
433
                                                    _mm_cvtpd_ps(v_sin_val_1)));
434
            _mm_storeu_ps(cosval + i, _mm_movelh_ps(_mm_cvtpd_ps(v_cos_val_0),
435
                                                    _mm_cvtpd_ps(v_cos_val_1)));
436
        }
437
    }
438
#endif
439

440
    for( ; i < len; i++ )
441
    {
442
        double t = angle[i]*k1;
443
        int it = cvRound(t);
444
        t -= it;
445
        int sin_idx = it & (N - 1);
446
        int cos_idx = (N/4 - sin_idx) & (N - 1);
447

448
        double sin_b = (sin_a0*t*t + sin_a2)*t;
449
        double cos_b = cos_a0*t*t + 1;
450

451
        double sin_a = sin_table[sin_idx];
452
        double cos_a = sin_table[cos_idx];
453

454
        double sin_val = sin_a*cos_b + cos_a*sin_b;
455
        double cos_val = cos_a*cos_b - sin_a*sin_b;
456

457
        sinval[i] = (float)sin_val;
458
        cosval[i] = (float)cos_val;
459
    }
460
}
461

462

463
#ifdef HAVE_OPENCL
464

465
static bool ocl_polarToCart( InputArray _mag, InputArray _angle,
466
                             OutputArray _dst1, OutputArray _dst2, bool angleInDegrees )
467
{
468
    const ocl::Device & d = ocl::Device::getDefault();
469
    int type = _angle.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type),
470
            rowsPerWI = d.isIntel() ? 4 : 1;
471
    bool doubleSupport = d.doubleFPConfig() > 0;
472

473
    if ( !doubleSupport && depth == CV_64F )
474
        return false;
475

476
    ocl::Kernel k("KF", ocl::core::arithm_oclsrc,
477
                  format("-D dstT=%s -D DEPTH_dst=%d -D rowsPerWI=%d -D BINARY_OP -D OP_PTC_%s%s",
478
                         ocl::typeToStr(CV_MAKE_TYPE(depth, 1)), depth,
479
                         rowsPerWI,
480
                         angleInDegrees ? "AD" : "AR",
481
                         doubleSupport ? " -D DOUBLE_SUPPORT" : ""));
482
    if (k.empty())
483
        return false;
484

485
    UMat mag = _mag.getUMat(), angle = _angle.getUMat();
486
    Size size = angle.size();
487
    CV_Assert(mag.size() == size);
488

489
    _dst1.create(size, type);
490
    _dst2.create(size, type);
491
    UMat dst1 = _dst1.getUMat(), dst2 = _dst2.getUMat();
492

493
    k.args(ocl::KernelArg::ReadOnlyNoSize(mag), ocl::KernelArg::ReadOnlyNoSize(angle),
494
           ocl::KernelArg::WriteOnly(dst1, cn), ocl::KernelArg::WriteOnlyNoSize(dst2));
495

496
    size_t globalsize[2] = { (size_t)dst1.cols * cn, ((size_t)dst1.rows + rowsPerWI - 1) / rowsPerWI };
497
    return k.run(2, globalsize, NULL, false);
498
}
499

500
#endif
501

502
#ifdef HAVE_IPP
503
static bool ipp_polarToCart(Mat &mag, Mat &angle, Mat &x, Mat &y)
504
{
505
    CV_INSTRUMENT_REGION_IPP();
506

507
    int depth = angle.depth();
508
    if(depth != CV_32F && depth != CV_64F)
509
        return false;
510

511
    if(angle.dims <= 2)
512
    {
513
        int len = (int)(angle.cols*angle.channels());
514

515
        if(depth == CV_32F)
516
        {
517
            for (int h = 0; h < angle.rows; h++)
518
            {
519
                if(CV_INSTRUMENT_FUN_IPP(ippsPolarToCart_32f, (const float*)mag.ptr(h), (const float*)angle.ptr(h), (float*)x.ptr(h), (float*)y.ptr(h), len) < 0)
520
                    return false;
521
            }
522
        }
523
        else
524
        {
525
            for (int h = 0; h < angle.rows; h++)
526
            {
527
                if(CV_INSTRUMENT_FUN_IPP(ippsPolarToCart_64f, (const double*)mag.ptr(h), (const double*)angle.ptr(h), (double*)x.ptr(h), (double*)y.ptr(h), len) < 0)
528
                    return false;
529
            }
530
        }
531
        return true;
532
    }
533
    else
534
    {
535
        const Mat      *arrays[] = {&mag, &angle, &x, &y, NULL};
536
        uchar          *ptrs[4]  = {NULL};
537
        NAryMatIterator it(arrays, ptrs);
538
        int len = (int)(it.size*angle.channels());
539

540
        if(depth == CV_32F)
541
        {
542
            for (size_t i = 0; i < it.nplanes; i++, ++it)
543
            {
544
                if(CV_INSTRUMENT_FUN_IPP(ippsPolarToCart_32f, (const float*)ptrs[0], (const float*)ptrs[1], (float*)ptrs[2], (float*)ptrs[3], len) < 0)
545
                    return false;
546
            }
547
        }
548
        else
549
        {
550
            for (size_t i = 0; i < it.nplanes; i++, ++it)
551
            {
552
                if(CV_INSTRUMENT_FUN_IPP(ippsPolarToCart_64f, (const double*)ptrs[0], (const double*)ptrs[1], (double*)ptrs[2], (double*)ptrs[3], len) < 0)
553
                    return false;
554
            }
555
        }
556
        return true;
557
    }
558
}
559
#endif
560

561
void polarToCart( InputArray src1, InputArray src2,
562
                  OutputArray dst1, OutputArray dst2, bool angleInDegrees )
563
{
564
    CV_INSTRUMENT_REGION();
565

566
    int type = src2.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
567
    CV_Assert((depth == CV_32F || depth == CV_64F) && (src1.empty() || src1.type() == type));
568

569
    CV_OCL_RUN(!src1.empty() && src2.dims() <= 2 && dst1.isUMat() && dst2.isUMat(),
570
               ocl_polarToCart(src1, src2, dst1, dst2, angleInDegrees))
571

572
    Mat Mag = src1.getMat(), Angle = src2.getMat();
573
    CV_Assert( Mag.empty() || Angle.size == Mag.size);
574
    dst1.create( Angle.dims, Angle.size, type );
575
    dst2.create( Angle.dims, Angle.size, type );
576
    Mat X = dst1.getMat(), Y = dst2.getMat();
577

578
    CV_IPP_RUN(!angleInDegrees, ipp_polarToCart(Mag, Angle, X, Y));
579

580
    const Mat* arrays[] = {&Mag, &Angle, &X, &Y, 0};
581
    uchar* ptrs[4] = {};
582
    NAryMatIterator it(arrays, ptrs);
583
    cv::AutoBuffer<float> _buf;
584
    float* buf[2] = {0, 0};
585
    int j, k, total = (int)(it.size*cn), blockSize = std::min(total, ((BLOCK_SIZE+cn-1)/cn)*cn);
586
    size_t esz1 = Angle.elemSize1();
587

588
    if( depth == CV_64F )
589
    {
590
        _buf.allocate(blockSize*2);
591
        buf[0] = _buf.data();
592
        buf[1] = buf[0] + blockSize;
593
    }
594

595
    for( size_t i = 0; i < it.nplanes; i++, ++it )
596
    {
597
        for( j = 0; j < total; j += blockSize )
598
        {
599
            int len = std::min(total - j, blockSize);
600
            if( depth == CV_32F )
601
            {
602
                const float *mag = (const float*)ptrs[0], *angle = (const float*)ptrs[1];
603
                float *x = (float*)ptrs[2], *y = (float*)ptrs[3];
604

605
                SinCos_32f( angle, y, x, len, angleInDegrees );
606
                if( mag )
607
                {
608
                    k = 0;
609

610
#if CV_SIMD
611
                    int cWidth = v_float32::nlanes;
612
                    for( ; k <= len - cWidth; k += cWidth )
613
                    {
614
                        v_float32 v_m = vx_load(mag + k);
615
                        v_store(x + k, vx_load(x + k) * v_m);
616
                        v_store(y + k, vx_load(y + k) * v_m);
617
                    }
618
                    vx_cleanup();
619
#endif
620

621
                    for( ; k < len; k++ )
622
                    {
623
                        float m = mag[k];
624
                        x[k] *= m; y[k] *= m;
625
                    }
626
                }
627
            }
628
            else
629
            {
630
                const double *mag = (const double*)ptrs[0], *angle = (const double*)ptrs[1];
631
                double *x = (double*)ptrs[2], *y = (double*)ptrs[3];
632

633
                for( k = 0; k < len; k++ )
634
                    buf[0][k] = (float)angle[k];
635

636
                SinCos_32f( buf[0], buf[1], buf[0], len, angleInDegrees );
637
                if( mag )
638
                    for( k = 0; k < len; k++ )
639
                    {
640
                        double m = mag[k];
641
                        x[k] = buf[0][k]*m; y[k] = buf[1][k]*m;
642
                    }
643
                else
644
                {
645
                    std::memcpy(x, buf[0], sizeof(float) * len);
646
                    std::memcpy(y, buf[1], sizeof(float) * len);
647
                }
648
            }
649

650
            if( ptrs[0] )
651
                ptrs[0] += len*esz1;
652
            ptrs[1] += len*esz1;
653
            ptrs[2] += len*esz1;
654
            ptrs[3] += len*esz1;
655
        }
656
    }
657
}
658

659
/****************************************************************************************\
660
*                                          E X P                                         *
661
\****************************************************************************************/
662

663
void exp( InputArray _src, OutputArray _dst )
664
{
665
    CV_INSTRUMENT_REGION();
666

667
    int type = _src.type(), depth = _src.depth(), cn = _src.channels();
668
    CV_Assert( depth == CV_32F || depth == CV_64F );
669

670
    CV_OCL_RUN(_dst.isUMat() && _src.dims() <= 2,
671
               ocl_math_op(_src, noArray(), _dst, OCL_OP_EXP))
672

673
    Mat src = _src.getMat();
674
    _dst.create( src.dims, src.size, type );
675
    Mat dst = _dst.getMat();
676

677
    const Mat* arrays[] = {&src, &dst, 0};
678
    uchar* ptrs[2] = {};
679
    NAryMatIterator it(arrays, ptrs);
680
    int len = (int)(it.size*cn);
681

682
    for( size_t i = 0; i < it.nplanes; i++, ++it )
683
    {
684
        if( depth == CV_32F )
685
            hal::exp32f((const float*)ptrs[0], (float*)ptrs[1], len);
686
        else
687
            hal::exp64f((const double*)ptrs[0], (double*)ptrs[1], len);
688
    }
689
}
690

691

692
/****************************************************************************************\
693
*                                          L O G                                         *
694
\****************************************************************************************/
695

696
void log( InputArray _src, OutputArray _dst )
697
{
698
    CV_INSTRUMENT_REGION();
699

700
    int type = _src.type(), depth = _src.depth(), cn = _src.channels();
701
    CV_Assert( depth == CV_32F || depth == CV_64F );
702

703
    CV_OCL_RUN( _dst.isUMat() && _src.dims() <= 2,
704
                ocl_math_op(_src, noArray(), _dst, OCL_OP_LOG))
705

706
    Mat src = _src.getMat();
707
    _dst.create( src.dims, src.size, type );
708
    Mat dst = _dst.getMat();
709

710
    const Mat* arrays[] = {&src, &dst, 0};
711
    uchar* ptrs[2] = {};
712
    NAryMatIterator it(arrays, ptrs);
713
    int len = (int)(it.size*cn);
714

715
    for( size_t i = 0; i < it.nplanes; i++, ++it )
716
    {
717
        if( depth == CV_32F )
718
            hal::log32f( (const float*)ptrs[0], (float*)ptrs[1], len );
719
        else
720
            hal::log64f( (const double*)ptrs[0], (double*)ptrs[1], len );
721
    }
722
}
723

724
/****************************************************************************************\
725
*                                    P O W E R                                           *
726
\****************************************************************************************/
727

728
template <typename T, typename WT>
729
struct iPow_SIMD
730
{
731
    int operator() ( const T *, T *, int, int)
732
    {
733
        return 0;
734
    }
735
};
736

737
#if CV_SIMD
738

739
template <>
740
struct iPow_SIMD<uchar, int>
741
{
742
    int operator() ( const uchar * src, uchar * dst, int len, int power )
743
    {
744
        int i = 0;
745
        v_uint32 v_1 = vx_setall_u32(1u);
746

747
        for ( ; i <= len - v_uint16::nlanes; i += v_uint16::nlanes)
748
        {
749
            v_uint32 v_a1 = v_1, v_a2 = v_1;
750
            v_uint16 v = vx_load_expand(src + i);
751
            v_uint32 v_b1, v_b2;
752
            v_expand(v, v_b1, v_b2);
753
            int p = power;
754

755
            while( p > 1 )
756
            {
757
                if (p & 1)
758
                {
759
                    v_a1 *= v_b1;
760
                    v_a2 *= v_b2;
761
                }
762
                v_b1 *= v_b1;
763
                v_b2 *= v_b2;
764
                p >>= 1;
765
            }
766

767
            v_a1 *= v_b1;
768
            v_a2 *= v_b2;
769

770
            v = v_pack(v_a1, v_a2);
771
            v_pack_store(dst + i, v);
772
        }
773
        vx_cleanup();
774

775
        return i;
776
    }
777
};
778

779
template <>
780
struct iPow_SIMD<schar, int>
781
{
782
    int operator() ( const schar * src, schar * dst, int len, int power)
783
    {
784
        int i = 0;
785
        v_int32 v_1 = vx_setall_s32(1);
786

787
        for ( ; i <= len - v_int16::nlanes; i += v_int16::nlanes)
788
        {
789
            v_int32 v_a1 = v_1, v_a2 = v_1;
790
            v_int16 v = vx_load_expand(src + i);
791
            v_int32 v_b1, v_b2;
792
            v_expand(v, v_b1, v_b2);
793
            int p = power;
794

795
            while( p > 1 )
796
            {
797
                if (p & 1)
798
                {
799
                    v_a1 *= v_b1;
800
                    v_a2 *= v_b2;
801
                }
802
                v_b1 *= v_b1;
803
                v_b2 *= v_b2;
804
                p >>= 1;
805
            }
806

807
            v_a1 *= v_b1;
808
            v_a2 *= v_b2;
809

810
            v = v_pack(v_a1, v_a2);
811
            v_pack_store(dst + i, v);
812
        }
813
        vx_cleanup();
814

815
        return i;
816
    }
817
};
818

819
template <>
820
struct iPow_SIMD<ushort, int>
821
{
822
    int operator() ( const ushort * src, ushort * dst, int len, int power)
823
    {
824
        int i = 0;
825
        v_uint32 v_1 = vx_setall_u32(1u);
826

827
        for ( ; i <= len - v_uint16::nlanes; i += v_uint16::nlanes)
828
        {
829
            v_uint32 v_a1 = v_1, v_a2 = v_1;
830
            v_uint16 v = vx_load(src + i);
831
            v_uint32 v_b1, v_b2;
832
            v_expand(v, v_b1, v_b2);
833
            int p = power;
834

835
            while( p > 1 )
836
            {
837
                if (p & 1)
838
                {
839
                    v_a1 *= v_b1;
840
                    v_a2 *= v_b2;
841
                }
842
                v_b1 *= v_b1;
843
                v_b2 *= v_b2;
844
                p >>= 1;
845
            }
846

847
            v_a1 *= v_b1;
848
            v_a2 *= v_b2;
849

850
            v = v_pack(v_a1, v_a2);
851
            v_store(dst + i, v);
852
        }
853
        vx_cleanup();
854

855
        return i;
856
    }
857
};
858

859
template <>
860
struct iPow_SIMD<short, int>
861
{
862
    int operator() ( const short * src, short * dst, int len, int power)
863
    {
864
        int i = 0;
865
        v_int32 v_1 = vx_setall_s32(1);
866

867
        for ( ; i <= len - v_int16::nlanes; i += v_int16::nlanes)
868
        {
869
            v_int32 v_a1 = v_1, v_a2 = v_1;
870
            v_int16 v = vx_load(src + i);
871
            v_int32 v_b1, v_b2;
872
            v_expand(v, v_b1, v_b2);
873
            int p = power;
874

875
            while( p > 1 )
876
            {
877
                if (p & 1)
878
                {
879
                    v_a1 *= v_b1;
880
                    v_a2 *= v_b2;
881
                }
882
                v_b1 *= v_b1;
883
                v_b2 *= v_b2;
884
                p >>= 1;
885
            }
886

887
            v_a1 *= v_b1;
888
            v_a2 *= v_b2;
889

890
            v = v_pack(v_a1, v_a2);
891
            v_store(dst + i, v);
892
        }
893
        vx_cleanup();
894

895
        return i;
896
    }
897
};
898

899
template <>
900
struct iPow_SIMD<int, int>
901
{
902
    int operator() ( const int * src, int * dst, int len, int power)
903
    {
904
        int i = 0;
905
        v_int32 v_1 = vx_setall_s32(1);
906

907
        for ( ; i <= len - v_int32::nlanes*2; i += v_int32::nlanes*2)
908
        {
909
            v_int32 v_a1 = v_1, v_a2 = v_1;
910
            v_int32 v_b1 = vx_load(src + i), v_b2 = vx_load(src + i + v_int32::nlanes);
911
            int p = power;
912

913
            while( p > 1 )
914
            {
915
                if (p & 1)
916
                {
917
                    v_a1 *= v_b1;
918
                    v_a2 *= v_b2;
919
                }
920
                v_b1 *= v_b1;
921
                v_b2 *= v_b2;
922
                p >>= 1;
923
            }
924

925
            v_a1 *= v_b1;
926
            v_a2 *= v_b2;
927

928
            v_store(dst + i, v_a1);
929
            v_store(dst + i + v_int32::nlanes, v_a2);
930
        }
931
        vx_cleanup();
932

933
        return i;
934
    }
935
};
936

937
template <>
938
struct iPow_SIMD<float, float>
939
{
940
    int operator() ( const float * src, float * dst, int len, int power)
941
    {
942
        int i = 0;
943
        v_float32 v_1 = vx_setall_f32(1.f);
944

945
        for ( ; i <= len - v_float32::nlanes*2; i += v_float32::nlanes*2)
946
        {
947
            v_float32 v_a1 = v_1, v_a2 = v_1;
948
            v_float32 v_b1 = vx_load(src + i), v_b2 = vx_load(src + i + v_float32::nlanes);
949
            int p = std::abs(power);
950
            if( power < 0 )
951
            {
952
                v_b1 = v_1 / v_b1;
953
                v_b2 = v_1 / v_b2;
954
            }
955

956
            while( p > 1 )
957
            {
958
                if (p & 1)
959
                {
960
                    v_a1 *= v_b1;
961
                    v_a2 *= v_b2;
962
                }
963
                v_b1 *= v_b1;
964
                v_b2 *= v_b2;
965
                p >>= 1;
966
            }
967

968
            v_a1 *= v_b1;
969
            v_a2 *= v_b2;
970

971
            v_store(dst + i, v_a1);
972
            v_store(dst + i + v_float32::nlanes, v_a2);
973
        }
974
        vx_cleanup();
975

976
        return i;
977
    }
978
};
979

980
#if CV_SIMD_64F
981
template <>
982
struct iPow_SIMD<double, double>
983
{
984
    int operator() ( const double * src, double * dst, int len, int power)
985
    {
986
        int i = 0;
987
        v_float64 v_1 = vx_setall_f64(1.);
988

989
        for ( ; i <= len - v_float64::nlanes*2; i += v_float64::nlanes*2)
990
        {
991
            v_float64 v_a1 = v_1, v_a2 = v_1;
992
            v_float64 v_b1 = vx_load(src + i), v_b2 = vx_load(src + i + v_float64::nlanes);
993
            int p = std::abs(power);
994
            if( power < 0 )
995
            {
996
                v_b1 = v_1 / v_b1;
997
                v_b2 = v_1 / v_b2;
998
            }
999

1000
            while( p > 1 )
1001
            {
1002
                if (p & 1)
1003
                {
1004
                    v_a1 *= v_b1;
1005
                    v_a2 *= v_b2;
1006
                }
1007
                v_b1 *= v_b1;
1008
                v_b2 *= v_b2;
1009
                p >>= 1;
1010
            }
1011

1012
            v_a1 *= v_b1;
1013
            v_a2 *= v_b2;
1014

1015
            v_store(dst + i, v_a1);
1016
            v_store(dst + i + v_float64::nlanes, v_a2);
1017
        }
1018
        vx_cleanup();
1019

1020
        return i;
1021
    }
1022
};
1023
#endif
1024

1025
#endif
1026

1027
template<typename T, typename WT>
1028
static void
1029
iPow_i( const T* src, T* dst, int len, int power )
1030
{
1031
    if( power < 0 )
1032
    {
1033
        T tab[5] =
1034
        {
1035
            saturate_cast<T>(power == -1 ? -1 : 0), saturate_cast<T>((power & 1) ? -1 : 1),
1036
            std::numeric_limits<T>::max(), 1, saturate_cast<T>(power == -1 ? 1 : 0)
1037
        };
1038
        for( int i = 0; i < len; i++ )
1039
        {
1040
            T val = src[i];
1041
            dst[i] = cv_abs(val) <= 2 ? tab[val + 2] : (T)0;
1042
        }
1043
    }
1044
    else
1045
    {
1046
        iPow_SIMD<T, WT> vop;
1047
        int i = vop(src, dst, len, power);
1048

1049
        for( ; i < len; i++ )
1050
        {
1051
            WT a = 1, b = src[i];
1052
            int p = power;
1053
            while( p > 1 )
1054
            {
1055
                if( p & 1 )
1056
                    a *= b;
1057
                b *= b;
1058
                p >>= 1;
1059
            }
1060

1061
            a *= b;
1062
            dst[i] = saturate_cast<T>(a);
1063
        }
1064
    }
1065
}
1066

1067
template<typename T>
1068
static void
1069
iPow_f( const T* src, T* dst, int len, int power0 )
1070
{
1071
    iPow_SIMD<T, T> vop;
1072
    int i = vop(src, dst, len, power0);
1073
    int power = std::abs(power0);
1074

1075
    for( ; i < len; i++ )
1076
    {
1077
        T a = 1, b = src[i];
1078
        int p = power;
1079
        if( power0 < 0 )
1080
            b = 1/b;
1081

1082
        while( p > 1 )
1083
        {
1084
            if( p & 1 )
1085
                a *= b;
1086
            b *= b;
1087
            p >>= 1;
1088
        }
1089

1090
        a *= b;
1091
        dst[i] = a;
1092
    }
1093
}
1094

1095
static void iPow8u(const uchar* src, uchar* dst, int len, int power)
1096
{
1097
    iPow_i<uchar, unsigned>(src, dst, len, power);
1098
}
1099

1100
static void iPow8s(const schar* src, schar* dst, int len, int power)
1101
{
1102
    iPow_i<schar, int>(src, dst, len, power);
1103
}
1104

1105
static void iPow16u(const ushort* src, ushort* dst, int len, int power)
1106
{
1107
    iPow_i<ushort, unsigned>(src, dst, len, power);
1108
}
1109

1110
static void iPow16s(const short* src, short* dst, int len, int power)
1111
{
1112
    iPow_i<short, int>(src, dst, len, power);
1113
}
1114

1115
static void iPow32s(const int* src, int* dst, int len, int power)
1116
{
1117
    iPow_i<int, int>(src, dst, len, power);
1118
}
1119

1120
static void iPow32f(const float* src, float* dst, int len, int power)
1121
{
1122
    iPow_f<float>(src, dst, len, power);
1123
}
1124

1125
static void iPow64f(const double* src, double* dst, int len, int power)
1126
{
1127
    iPow_f<double>(src, dst, len, power);
1128
}
1129

1130

1131
typedef void (*IPowFunc)( const uchar* src, uchar* dst, int len, int power );
1132

1133
static IPowFunc ipowTab[] =
1134
{
1135
    (IPowFunc)iPow8u, (IPowFunc)iPow8s, (IPowFunc)iPow16u, (IPowFunc)iPow16s,
1136
    (IPowFunc)iPow32s, (IPowFunc)iPow32f, (IPowFunc)iPow64f, 0
1137
};
1138

1139
#ifdef HAVE_OPENCL
1140

1141
static bool ocl_pow(InputArray _src, double power, OutputArray _dst,
1142
                    bool is_ipower, int ipower)
1143
{
1144
    const ocl::Device & d = ocl::Device::getDefault();
1145
    int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type),
1146
            rowsPerWI = d.isIntel() ? 4 : 1;
1147
    bool doubleSupport = d.doubleFPConfig() > 0;
1148

1149
    _dst.createSameSize(_src, type);
1150
    if (is_ipower)
1151
    {
1152
        if (ipower == 0)
1153
        {
1154
            _dst.setTo(Scalar::all(1));
1155
            return true;
1156
        }
1157
        if (ipower == 1)
1158
        {
1159
            _src.copyTo(_dst);
1160
            return true;
1161
        }
1162
        if( ipower < 0 )
1163
        {
1164
            if( depth == CV_32F || depth == CV_64F )
1165
                is_ipower = false;
1166
            else
1167
                return false;
1168
        }
1169
    }
1170

1171
    if (depth == CV_64F && !doubleSupport)
1172
        return false;
1173

1174
    bool issqrt = std::abs(power - 0.5) < DBL_EPSILON;
1175
    const char * const op = issqrt ? "OP_SQRT" : is_ipower ? "OP_POWN" : "OP_POW";
1176

1177
    ocl::Kernel k("KF", ocl::core::arithm_oclsrc,
1178
                  format("-D dstT=%s -D DEPTH_dst=%d -D rowsPerWI=%d -D %s -D UNARY_OP%s",
1179
                         ocl::typeToStr(depth), depth, rowsPerWI, op,
1180
                         doubleSupport ? " -D DOUBLE_SUPPORT" : ""));
1181
    if (k.empty())
1182
        return false;
1183

1184
    UMat src = _src.getUMat();
1185
    _dst.create(src.size(), type);
1186
    UMat dst = _dst.getUMat();
1187

1188
    ocl::KernelArg srcarg = ocl::KernelArg::ReadOnlyNoSize(src),
1189
            dstarg = ocl::KernelArg::WriteOnly(dst, cn);
1190

1191
    if (issqrt)
1192
        k.args(srcarg, dstarg);
1193
    else if (is_ipower)
1194
        k.args(srcarg, dstarg, ipower);
1195
    else
1196
    {
1197
        if (depth == CV_32F)
1198
            k.args(srcarg, dstarg, (float)power);
1199
        else
1200
            k.args(srcarg, dstarg, power);
1201
    }
1202

1203
    size_t globalsize[2] = { (size_t)dst.cols *  cn, ((size_t)dst.rows + rowsPerWI - 1) / rowsPerWI };
1204
    return k.run(2, globalsize, NULL, false);
1205
}
1206

1207
#endif
1208

1209
void pow( InputArray _src, double power, OutputArray _dst )
1210
{
1211
    CV_INSTRUMENT_REGION();
1212

1213
    int type = _src.type(), depth = CV_MAT_DEPTH(type),
1214
            cn = CV_MAT_CN(type), ipower = cvRound(power);
1215
    bool is_ipower = fabs(ipower - power) < DBL_EPSILON;
1216
#ifdef HAVE_OPENCL
1217
    bool useOpenCL = _dst.isUMat() && _src.dims() <= 2;
1218
#endif
1219

1220
    if( is_ipower
1221
#ifdef HAVE_OPENCL
1222
            && !(useOpenCL && ocl::Device::getDefault().isIntel() && depth != CV_64F)
1223
#endif
1224
      )
1225
    {
1226
        switch( ipower )
1227
        {
1228
        case 0:
1229
            _dst.createSameSize(_src, type);
1230
            _dst.setTo(Scalar::all(1));
1231
            return;
1232
        case 1:
1233
            _src.copyTo(_dst);
1234
            return;
1235
        case 2:
1236
            multiply(_src, _src, _dst);
1237
            return;
1238
        }
1239
    }
1240
    else
1241
        CV_Assert( depth == CV_32F || depth == CV_64F );
1242

1243
    CV_OCL_RUN(useOpenCL, ocl_pow(_src, power, _dst, is_ipower, ipower))
1244

1245
    Mat src = _src.getMat();
1246
    _dst.create( src.dims, src.size, type );
1247
    Mat dst = _dst.getMat();
1248

1249
    const Mat* arrays[] = {&src, &dst, 0};
1250
    uchar* ptrs[2] = {};
1251
    NAryMatIterator it(arrays, ptrs);
1252
    int len = (int)(it.size*cn);
1253

1254
    if( is_ipower )
1255
    {
1256
        IPowFunc func = ipowTab[depth];
1257
        CV_Assert( func != 0 );
1258

1259
        for( size_t i = 0; i < it.nplanes; i++, ++it )
1260
            func( ptrs[0], ptrs[1], len, ipower );
1261
    }
1262
    else if( fabs(fabs(power) - 0.5) < DBL_EPSILON )
1263
    {
1264
        MathFunc func = power < 0 ?
1265
            (depth == CV_32F ? (MathFunc)hal::invSqrt32f : (MathFunc)hal::invSqrt64f) :
1266
            (depth == CV_32F ? (MathFunc)hal::sqrt32f : (MathFunc)hal::sqrt64f);
1267

1268
        for( size_t i = 0; i < it.nplanes; i++, ++it )
1269
            func( ptrs[0], ptrs[1], len );
1270
    }
1271
    else
1272
    {
1273
        int j, k, blockSize = std::min(len, ((BLOCK_SIZE + cn-1)/cn)*cn);
1274
        size_t esz1 = src.elemSize1();
1275
        AutoBuffer<uchar> buf;
1276
        Cv32suf inf32, nan32;
1277
        Cv64suf inf64, nan64;
1278
        float* fbuf = 0;
1279
        double* dbuf = 0;
1280
#ifndef __EMSCRIPTEN__
1281
        inf32.i = 0x7f800000;
1282
        nan32.i = 0x7fffffff;
1283
        inf64.i = CV_BIG_INT(0x7FF0000000000000);
1284
        nan64.i = CV_BIG_INT(0x7FFFFFFFFFFFFFFF);
1285
#else
1286
        inf32.f = std::numeric_limits<float>::infinity();
1287
        nan32.f = std::numeric_limits<float>::quiet_NaN();
1288
        inf64.f = std::numeric_limits<double>::infinity();
1289
        nan64.f = std::numeric_limits<double>::quiet_NaN();
1290
#endif
1291

1292
        if( src.ptr() == dst.ptr() )
1293
        {
1294
            buf.allocate(blockSize*esz1);
1295
            fbuf = (float*)buf.data();
1296
            dbuf = (double*)buf.data();
1297
        }
1298

1299
        for( size_t i = 0; i < it.nplanes; i++, ++it )
1300
        {
1301
            for( j = 0; j < len; j += blockSize )
1302
            {
1303
                int bsz = std::min(len - j, blockSize);
1304

1305
                if( depth == CV_32F )
1306
                {
1307
                    float* x0 = (float*)ptrs[0];
1308
                    float* x = fbuf ? fbuf : x0;
1309
                    float* y = (float*)ptrs[1];
1310

1311
                    if( x != x0 )
1312
                        memcpy(x, x0, bsz*esz1);
1313

1314
                    hal::log32f(x, y, bsz);
1315
                    for( k = 0; k < bsz; k++ )
1316
                        y[k] = (float)(y[k]*power);
1317
                    hal::exp32f(y, y, bsz);
1318
                    for( k = 0; k < bsz; k++ )
1319
                    {
1320
                        if( x0[k] <= 0 )
1321
                        {
1322
                            if( x0[k] == 0.f )
1323
                            {
1324
                                if( power < 0 )
1325
                                    y[k] = inf32.f;
1326
                            }
1327
                            else
1328
                                y[k] = nan32.f;
1329
                        }
1330
                    }
1331
                }
1332
                else
1333
                {
1334
                    double* x0 = (double*)ptrs[0];
1335
                    double* x = dbuf ? dbuf : x0;
1336
                    double* y = (double*)ptrs[1];
1337

1338
                    if( x != x0 )
1339
                        memcpy(x, x0, bsz*esz1);
1340

1341
                    hal::log64f(x, y, bsz);
1342
                    for( k = 0; k < bsz; k++ )
1343
                        y[k] *= power;
1344
                    hal::exp64f(y, y, bsz);
1345

1346
                    for( k = 0; k < bsz; k++ )
1347
                    {
1348
                        if( x0[k] <= 0 )
1349
                        {
1350
                            if( x0[k] == 0. )
1351
                            {
1352
                                if( power < 0 )
1353
                                    y[k] = inf64.f;
1354
                            }
1355
                            else
1356
                                y[k] = nan64.f;
1357
                        }
1358
                    }
1359
                }
1360
                ptrs[0] += bsz*esz1;
1361
                ptrs[1] += bsz*esz1;
1362
            }
1363
        }
1364
    }
1365
}
1366

1367
void sqrt(InputArray a, OutputArray b)
1368
{
1369
    CV_INSTRUMENT_REGION();
1370

1371
    cv::pow(a, 0.5, b);
1372
}
1373

1374
/************************** CheckArray for NaN's, Inf's *********************************/
1375

1376
template<int cv_mat_type> struct mat_type_assotiations{};
1377

1378
template<> struct mat_type_assotiations<CV_8U>
1379
{
1380
    typedef unsigned char type;
1381
    static const type min_allowable = 0x0;
1382
    static const type max_allowable = 0xFF;
1383
};
1384

1385
template<> struct mat_type_assotiations<CV_8S>
1386
{
1387
    typedef signed char type;
1388
    static const type min_allowable = SCHAR_MIN;
1389
    static const type max_allowable = SCHAR_MAX;
1390
};
1391

1392
template<> struct mat_type_assotiations<CV_16U>
1393
{
1394
    typedef unsigned short type;
1395
    static const type min_allowable = 0x0;
1396
    static const type max_allowable = USHRT_MAX;
1397
};
1398
template<> struct mat_type_assotiations<CV_16S>
1399
{
1400
    typedef signed short type;
1401
    static const type min_allowable = SHRT_MIN;
1402
    static const type max_allowable = SHRT_MAX;
1403
};
1404

1405
template<> struct mat_type_assotiations<CV_32S>
1406
{
1407
    typedef int type;
1408
    static const type min_allowable = (-INT_MAX - 1);
1409
    static const type max_allowable = INT_MAX;
1410
};
1411

1412
template<int depth>
1413
static bool checkIntegerRange(cv::Mat src, Point& bad_pt, int minVal, int maxVal)
1414
{
1415
    typedef mat_type_assotiations<depth> type_ass;
1416

1417
    if (minVal < type_ass::min_allowable && maxVal > type_ass::max_allowable)
1418
    {
1419
        return true;
1420
    }
1421
    else if (minVal > type_ass::max_allowable || maxVal < type_ass::min_allowable || maxVal < minVal)
1422
    {
1423
        bad_pt = cv::Point(0,0);
1424
        return false;
1425
    }
1426
    cv::Mat as_one_channel = src.reshape(1,0);
1427

1428
    for (int j = 0; j < as_one_channel.rows; ++j)
1429
        for (int i = 0; i < as_one_channel.cols; ++i)
1430
        {
1431
            typename type_ass::type v = as_one_channel.at<typename type_ass::type>(j ,i);
1432
            if (v < minVal || v > maxVal)
1433
            {
1434
                bad_pt.y = j;
1435
                bad_pt.x = i / src.channels();
1436
                return false;
1437
            }
1438
        }
1439

1440
    return true;
1441
}
1442

1443
typedef bool (*check_range_function)(cv::Mat src, Point& bad_pt, int minVal, int maxVal);
1444

1445
check_range_function check_range_functions[] =
1446
{
1447
    &checkIntegerRange<CV_8U>,
1448
    &checkIntegerRange<CV_8S>,
1449
    &checkIntegerRange<CV_16U>,
1450
    &checkIntegerRange<CV_16S>,
1451
    &checkIntegerRange<CV_32S>
1452
};
1453

1454
bool checkRange(InputArray _src, bool quiet, Point* pt, double minVal, double maxVal)
1455
{
1456
    CV_INSTRUMENT_REGION();
1457

1458
    Mat src = _src.getMat();
1459

1460
    if ( src.dims > 2 )
1461
    {
1462
        CV_Assert(pt == NULL); // no way to provide location info
1463

1464
        const Mat* arrays[] = {&src, 0};
1465
        Mat planes[1];
1466
        NAryMatIterator it(arrays, planes);
1467

1468
        for ( size_t i = 0; i < it.nplanes; i++, ++it )
1469
        {
1470
            if (!checkRange( it.planes[0], quiet, NULL, minVal, maxVal ))
1471
            {
1472
                return false;
1473
            }
1474
        }
1475
        return true;
1476
    }
1477

1478
    int depth = src.depth();
1479
    Point badPt(-1, -1);
1480

1481
    if (depth < CV_32F)
1482
    {
1483
        int minVali = minVal <= INT_MIN ? INT_MIN : cvFloor(minVal);
1484
        int maxVali = maxVal > INT_MAX ? INT_MAX : cvCeil(maxVal) - 1;
1485

1486
        (check_range_functions[depth])(src, badPt, minVali, maxVali);
1487
    }
1488
    else
1489
    {
1490
        int i, loc = 0;
1491
        int cn = src.channels();
1492
        Size size = getContinuousSize( src, cn );
1493

1494
        if( depth == CV_32F )
1495
        {
1496
            Cv32suf a, b;
1497
            int ia, ib;
1498
            const int* isrc = src.ptr<int>();
1499
            size_t step = src.step/sizeof(isrc[0]);
1500

1501
            a.f = (float)std::max(minVal, (double)-FLT_MAX);
1502
            b.f = (float)std::min(maxVal, (double)FLT_MAX);
1503

1504
            ia = CV_TOGGLE_FLT(a.i);
1505
            ib = CV_TOGGLE_FLT(b.i);
1506

1507
            for( ; badPt.x < 0 && size.height--; loc += size.width, isrc += step )
1508
            {
1509
                for( i = 0; i < size.width; i++ )
1510
                {
1511
                    int val = isrc[i];
1512
                    val = CV_TOGGLE_FLT(val);
1513

1514
                    if( val < ia || val >= ib )
1515
                    {
1516
                        int pixelId = (loc + i) / cn;
1517
                        badPt = Point(pixelId % src.cols, pixelId / src.cols);
1518
                        break;
1519
                    }
1520
                }
1521
            }
1522
        }
1523
        else
1524
        {
1525
            Cv64suf a, b;
1526
            int64 ia, ib;
1527
            const int64* isrc = src.ptr<int64>();
1528
            size_t step = src.step/sizeof(isrc[0]);
1529

1530
            a.f = minVal;
1531
            b.f = maxVal;
1532

1533
            ia = CV_TOGGLE_DBL(a.i);
1534
            ib = CV_TOGGLE_DBL(b.i);
1535

1536
            for( ; badPt.x < 0 && size.height--; loc += size.width, isrc += step )
1537
            {
1538
                for( i = 0; i < size.width; i++ )
1539
                {
1540
                    int64 val = isrc[i];
1541
                    val = CV_TOGGLE_DBL(val);
1542

1543
                    if( val < ia || val >= ib )
1544
                    {
1545
                        int pixelId = (loc + i) / cn;
1546
                        badPt = Point(pixelId % src.cols, pixelId / src.cols);
1547
                        break;
1548
                    }
1549
                }
1550
            }
1551
        }
1552
    }
1553

1554
    if( badPt.x >= 0 )
1555
    {
1556
        if( pt )
1557
            *pt = badPt;
1558
        if( !quiet )
1559
        {
1560
            cv::String value_str;
1561
            value_str << src(cv::Range(badPt.y, badPt.y + 1), cv::Range(badPt.x, badPt.x + 1));
1562
            CV_Error_( CV_StsOutOfRange,
1563
            ("the value at (%d, %d)=%s is out of range [%f, %f)", badPt.x, badPt.y, value_str.c_str(), minVal, maxVal));
1564
        }
1565
        return false;
1566
    }
1567
    return true;
1568
}
1569

1570
#ifdef HAVE_OPENCL
1571

1572
static bool ocl_patchNaNs( InputOutputArray _a, float value )
1573
{
1574
    int rowsPerWI = ocl::Device::getDefault().isIntel() ? 4 : 1;
1575
    ocl::Kernel k("KF", ocl::core::arithm_oclsrc,
1576
                     format("-D UNARY_OP -D OP_PATCH_NANS -D dstT=float -D DEPTH_dst=%d -D rowsPerWI=%d",
1577
                            CV_32F, rowsPerWI));
1578
    if (k.empty())
1579
        return false;
1580

1581
    UMat a = _a.getUMat();
1582
    int cn = a.channels();
1583

1584
    k.args(ocl::KernelArg::ReadOnlyNoSize(a),
1585
           ocl::KernelArg::WriteOnly(a, cn), (float)value);
1586

1587
    size_t globalsize[2] = { (size_t)a.cols * cn, ((size_t)a.rows + rowsPerWI - 1) / rowsPerWI };
1588
    return k.run(2, globalsize, NULL, false);
1589
}
1590

1591
#endif
1592

1593
void patchNaNs( InputOutputArray _a, double _val )
1594
{
1595
    CV_INSTRUMENT_REGION();
1596

1597
    CV_Assert( _a.depth() == CV_32F );
1598

1599
    CV_OCL_RUN(_a.isUMat() && _a.dims() <= 2,
1600
               ocl_patchNaNs(_a, (float)_val))
1601

1602
    Mat a = _a.getMat();
1603
    const Mat* arrays[] = {&a, 0};
1604
    int* ptrs[1] = {};
1605
    NAryMatIterator it(arrays, (uchar**)ptrs);
1606
    size_t len = it.size*a.channels();
1607
    Cv32suf val;
1608
    val.f = (float)_val;
1609

1610
#if CV_SIMD
1611
    v_int32 v_mask1 = vx_setall_s32(0x7fffffff), v_mask2 = vx_setall_s32(0x7f800000);
1612
    v_int32 v_val = vx_setall_s32(val.i);
1613
#endif
1614

1615
    for( size_t i = 0; i < it.nplanes; i++, ++it )
1616
    {
1617
        int* tptr = ptrs[0];
1618
        size_t j = 0;
1619

1620
#if CV_SIMD
1621
        size_t cWidth = (size_t)v_int32::nlanes;
1622
        for ( ; j + cWidth <= len; j += cWidth)
1623
        {
1624
            v_int32 v_src = vx_load(tptr + j);
1625
            v_int32 v_cmp_mask = v_mask2 < (v_src & v_mask1);
1626
            v_int32 v_dst = v_select(v_cmp_mask, v_val, v_src);
1627
            v_store(tptr + j, v_dst);
1628
        }
1629
        vx_cleanup();
1630
#endif
1631

1632
        for( ; j < len; j++ )
1633
            if( (tptr[j] & 0x7fffffff) > 0x7f800000 )
1634
                tptr[j] = val.i;
1635
    }
1636
}
1637

1638
}
1639

1640
CV_IMPL float cvCbrt(float value) { return cv::cubeRoot(value); }
1641
CV_IMPL float cvFastArctan(float y, float x) { return cv::fastAtan2(y, x); }
1642

1643
CV_IMPL void
1644
cvCartToPolar( const CvArr* xarr, const CvArr* yarr,
1645
               CvArr* magarr, CvArr* anglearr,
1646
               int angle_in_degrees )
1647
{
1648
    cv::Mat X = cv::cvarrToMat(xarr), Y = cv::cvarrToMat(yarr), Mag, Angle;
1649
    if( magarr )
1650
    {
1651
        Mag = cv::cvarrToMat(magarr);
1652
        CV_Assert( Mag.size() == X.size() && Mag.type() == X.type() );
1653
    }
1654
    if( anglearr )
1655
    {
1656
        Angle = cv::cvarrToMat(anglearr);
1657
        CV_Assert( Angle.size() == X.size() && Angle.type() == X.type() );
1658
    }
1659
    if( magarr )
1660
    {
1661
        if( anglearr )
1662
            cv::cartToPolar( X, Y, Mag, Angle, angle_in_degrees != 0 );
1663
        else
1664
            cv::magnitude( X, Y, Mag );
1665
    }
1666
    else
1667
        cv::phase( X, Y, Angle, angle_in_degrees != 0 );
1668
}
1669

1670
CV_IMPL void
1671
cvPolarToCart( const CvArr* magarr, const CvArr* anglearr,
1672
               CvArr* xarr, CvArr* yarr, int angle_in_degrees )
1673
{
1674
    cv::Mat X, Y, Angle = cv::cvarrToMat(anglearr), Mag;
1675
    if( magarr )
1676
    {
1677
        Mag = cv::cvarrToMat(magarr);
1678
        CV_Assert( Mag.size() == Angle.size() && Mag.type() == Angle.type() );
1679
    }
1680
    if( xarr )
1681
    {
1682
        X = cv::cvarrToMat(xarr);
1683
        CV_Assert( X.size() == Angle.size() && X.type() == Angle.type() );
1684
    }
1685
    if( yarr )
1686
    {
1687
        Y = cv::cvarrToMat(yarr);
1688
        CV_Assert( Y.size() == Angle.size() && Y.type() == Angle.type() );
1689
    }
1690

1691
    cv::polarToCart( Mag, Angle, X, Y, angle_in_degrees != 0 );
1692
}
1693

1694
CV_IMPL void cvExp( const CvArr* srcarr, CvArr* dstarr )
1695
{
1696
    cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
1697
    CV_Assert( src.type() == dst.type() && src.size == dst.size );
1698
    cv::exp( src, dst );
1699
}
1700

1701
CV_IMPL void cvLog( const CvArr* srcarr, CvArr* dstarr )
1702
{
1703
    cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
1704
    CV_Assert( src.type() == dst.type() && src.size == dst.size );
1705
    cv::log( src, dst );
1706
}
1707

1708
CV_IMPL void cvPow( const CvArr* srcarr, CvArr* dstarr, double power )
1709
{
1710
    cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
1711
    CV_Assert( src.type() == dst.type() && src.size == dst.size );
1712
    cv::pow( src, power, dst );
1713
}
1714

1715
CV_IMPL int cvCheckArr( const CvArr* arr, int flags,
1716
                        double minVal, double maxVal )
1717
{
1718
    if( (flags & CV_CHECK_RANGE) == 0 )
1719
        minVal = -DBL_MAX, maxVal = DBL_MAX;
1720
    return cv::checkRange(cv::cvarrToMat(arr), (flags & CV_CHECK_QUIET) != 0, 0, minVal, maxVal );
1721
}
1722

1723

1724
/*
1725
  Finds real roots of cubic, quadratic or linear equation.
1726
  The original code has been taken from Ken Turkowski web page
1727
  (http://www.worldserver.com/turk/opensource/) and adopted for OpenCV.
1728
  Here is the copyright notice.
1729

1730
  -----------------------------------------------------------------------
1731
  Copyright (C) 1978-1999 Ken Turkowski. <[email protected]>
1732

1733
    All rights reserved.
1734

1735
    Warranty Information
1736
      Even though I have reviewed this software, I make no warranty
1737
      or representation, either express or implied, with respect to this
1738
      software, its quality, accuracy, merchantability, or fitness for a
1739
      particular purpose.  As a result, this software is provided "as is,"
1740
      and you, its user, are assuming the entire risk as to its quality
1741
      and accuracy.
1742

1743
    This code may be used and freely distributed as long as it includes
1744
    this copyright notice and the above warranty information.
1745
  -----------------------------------------------------------------------
1746
*/
1747

1748
int cv::solveCubic( InputArray _coeffs, OutputArray _roots )
1749
{
1750
    CV_INSTRUMENT_REGION();
1751

1752
    const int n0 = 3;
1753
    Mat coeffs = _coeffs.getMat();
1754
    int ctype = coeffs.type();
1755

1756
    CV_Assert( ctype == CV_32F || ctype == CV_64F );
1757
    CV_Assert( (coeffs.size() == Size(n0, 1) ||
1758
                coeffs.size() == Size(n0+1, 1) ||
1759
                coeffs.size() == Size(1, n0) ||
1760
                coeffs.size() == Size(1, n0+1)) );
1761

1762
    _roots.create(n0, 1, ctype, -1, true, _OutputArray::DEPTH_MASK_FLT);
1763
    Mat roots = _roots.getMat();
1764

1765
    int i = -1, n = 0;
1766
    double a0 = 1., a1, a2, a3;
1767
    double x0 = 0., x1 = 0., x2 = 0.;
1768
    int ncoeffs = coeffs.rows + coeffs.cols - 1;
1769

1770
    if( ctype == CV_32FC1 )
1771
    {
1772
        if( ncoeffs == 4 )
1773
            a0 = coeffs.at<float>(++i);
1774

1775
        a1 = coeffs.at<float>(i+1);
1776
        a2 = coeffs.at<float>(i+2);
1777
        a3 = coeffs.at<float>(i+3);
1778
    }
1779
    else
1780
    {
1781
        if( ncoeffs == 4 )
1782
            a0 = coeffs.at<double>(++i);
1783

1784
        a1 = coeffs.at<double>(i+1);
1785
        a2 = coeffs.at<double>(i+2);
1786
        a3 = coeffs.at<double>(i+3);
1787
    }
1788

1789
    if( a0 == 0 )
1790
    {
1791
        if( a1 == 0 )
1792
        {
1793
            if( a2 == 0 )
1794
                n = a3 == 0 ? -1 : 0;
1795
            else
1796
            {
1797
                // linear equation
1798
                x0 = -a3/a2;
1799
                n = 1;
1800
            }
1801
        }
1802
        else
1803
        {
1804
            // quadratic equation
1805
            double d = a2*a2 - 4*a1*a3;
1806
            if( d >= 0 )
1807
            {
1808
                d = std::sqrt(d);
1809
                double q1 = (-a2 + d) * 0.5;
1810
                double q2 = (a2 + d) * -0.5;
1811
                if( fabs(q1) > fabs(q2) )
1812
                {
1813
                    x0 = q1 / a1;
1814
                    x1 = a3 / q1;
1815
                }
1816
                else
1817
                {
1818
                    x0 = q2 / a1;
1819
                    x1 = a3 / q2;
1820
                }
1821
                n = d > 0 ? 2 : 1;
1822
            }
1823
        }
1824
    }
1825
    else
1826
    {
1827
        a0 = 1./a0;
1828
        a1 *= a0;
1829
        a2 *= a0;
1830
        a3 *= a0;
1831

1832
        double Q = (a1 * a1 - 3 * a2) * (1./9);
1833
        double R = (2 * a1 * a1 * a1 - 9 * a1 * a2 + 27 * a3) * (1./54);
1834
        double Qcubed = Q * Q * Q;
1835
        double d = Qcubed - R * R;
1836

1837
        if( d > 0 )
1838
        {
1839
            double theta = acos(R / sqrt(Qcubed));
1840
            double sqrtQ = sqrt(Q);
1841
            double t0 = -2 * sqrtQ;
1842
            double t1 = theta * (1./3);
1843
            double t2 = a1 * (1./3);
1844
            x0 = t0 * cos(t1) - t2;
1845
            x1 = t0 * cos(t1 + (2.*CV_PI/3)) - t2;
1846
            x2 = t0 * cos(t1 + (4.*CV_PI/3)) - t2;
1847
            n = 3;
1848
        }
1849
        else if( d == 0 )
1850
        {
1851
            if(R >= 0)
1852
            {
1853
                x0 = -2*pow(R, 1./3) - a1/3;
1854
                x1 = pow(R, 1./3) - a1/3;
1855
            }
1856
            else
1857
            {
1858
                x0 = 2*pow(-R, 1./3) - a1/3;
1859
                x1 = -pow(-R, 1./3) - a1/3;
1860
            }
1861
            x2 = 0;
1862
            n = x0 == x1 ? 1 : 2;
1863
            x1 = x0 == x1 ? 0 : x1;
1864
        }
1865
        else
1866
        {
1867
            double e;
1868
            d = sqrt(-d);
1869
            e = pow(d + fabs(R), 1./3);
1870
            if( R > 0 )
1871
                e = -e;
1872
            x0 = (e + Q / e) - a1 * (1./3);
1873
            n = 1;
1874
        }
1875
    }
1876

1877
    if( roots.type() == CV_32FC1 )
1878
    {
1879
        roots.at<float>(0) = (float)x0;
1880
        roots.at<float>(1) = (float)x1;
1881
        roots.at<float>(2) = (float)x2;
1882
    }
1883
    else
1884
    {
1885
        roots.at<double>(0) = x0;
1886
        roots.at<double>(1) = x1;
1887
        roots.at<double>(2) = x2;
1888
    }
1889

1890
    return n;
1891
}
1892

1893
/* finds complex roots of a polynomial using Durand-Kerner method:
1894
   http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method */
1895
double cv::solvePoly( InputArray _coeffs0, OutputArray _roots0, int maxIters )
1896
{
1897
    CV_INSTRUMENT_REGION();
1898

1899
    typedef Complex<double> C;
1900

1901
    double maxDiff = 0;
1902
    int iter, i, j;
1903
    Mat coeffs0 = _coeffs0.getMat();
1904
    int ctype = _coeffs0.type();
1905
    int cdepth = CV_MAT_DEPTH(ctype);
1906

1907
    CV_Assert( CV_MAT_DEPTH(ctype) >= CV_32F && CV_MAT_CN(ctype) <= 2 );
1908
    CV_Assert( coeffs0.rows == 1 || coeffs0.cols == 1 );
1909

1910
    int n0 = coeffs0.cols + coeffs0.rows - 2, n = n0;
1911

1912
    _roots0.create(n, 1, CV_MAKETYPE(cdepth, 2), -1, true, _OutputArray::DEPTH_MASK_FLT);
1913
    Mat roots0 = _roots0.getMat();
1914

1915
    AutoBuffer<C> buf(n*2+2);
1916
    C *coeffs = buf.data(), *roots = coeffs + n + 1;
1917
    Mat coeffs1(coeffs0.size(), CV_MAKETYPE(CV_64F, coeffs0.channels()), coeffs0.channels() == 2 ? coeffs : roots);
1918
    coeffs0.convertTo(coeffs1, coeffs1.type());
1919
    if( coeffs0.channels() == 1 )
1920
    {
1921
        const double* rcoeffs = (const double*)roots;
1922
        for( i = 0; i <= n; i++ )
1923
            coeffs[i] = C(rcoeffs[i], 0);
1924
    }
1925

1926
    for( ; n > 1; n-- )
1927
    {
1928
        if( std::abs(coeffs[n].re) + std::abs(coeffs[n].im) > DBL_EPSILON )
1929
            break;
1930
    }
1931

1932
    C p(1, 0), r(1, 1);
1933

1934
    for( i = 0; i < n; i++ )
1935
    {
1936
        roots[i] = p;
1937
        p = p * r;
1938
    }
1939

1940
    maxIters = maxIters <= 0 ? 1000 : maxIters;
1941
    for( iter = 0; iter < maxIters; iter++ )
1942
    {
1943
        maxDiff = 0;
1944
        for( i = 0; i < n; i++ )
1945
        {
1946
            p = roots[i];
1947
            C num = coeffs[n], denom = coeffs[n];
1948
            int num_same_root = 1;
1949
            for( j = 0; j < n; j++ )
1950
            {
1951
                num = num*p + coeffs[n-j-1];
1952
                if( j != i )
1953
                {
1954
                    if ( (p - roots[j]).re != 0 || (p - roots[j]).im != 0 )
1955
                        denom = denom * (p - roots[j]);
1956
                    else
1957
                        num_same_root++;
1958
                }
1959
            }
1960
            num /= denom;
1961
            if( num_same_root > 1)
1962
            {
1963
                double old_num_re = num.re;
1964
                double old_num_im = num.im;
1965
                int square_root_times = num_same_root % 2 == 0 ? num_same_root / 2 : num_same_root / 2 - 1;
1966

1967
                for( j = 0; j < square_root_times; j++)
1968
                {
1969
                    num.re = old_num_re*old_num_re + old_num_im*old_num_im;
1970
                    num.re = sqrt(num.re);
1971
                    num.re += old_num_re;
1972
                    num.im = num.re - old_num_re;
1973
                    num.re /= 2;
1974
                    num.re = sqrt(num.re);
1975

1976
                    num.im /= 2;
1977
                    num.im = sqrt(num.im);
1978
                    if( old_num_re < 0 ) num.im = -num.im;
1979
                }
1980

1981
                if( num_same_root % 2 != 0){
1982
                    Mat cube_coefs(4, 1, CV_64FC1);
1983
                    Mat cube_roots(3, 1, CV_64FC2);
1984
                    cube_coefs.at<double>(3) = -(pow(old_num_re, 3));
1985
                    cube_coefs.at<double>(2) = -(15*pow(old_num_re, 2) + 27*pow(old_num_im, 2));
1986
                    cube_coefs.at<double>(1) = -48*old_num_re;
1987
                    cube_coefs.at<double>(0) = 64;
1988
                    solveCubic(cube_coefs, cube_roots);
1989

1990
                    if(cube_roots.at<double>(0) >= 0) num.re = pow(cube_roots.at<double>(0), 1./3);
1991
                    else num.re = -pow(-cube_roots.at<double>(0), 1./3);
1992
                    num.im = sqrt(pow(num.re, 2) / 3 - old_num_re / (3*num.re));
1993
                }
1994
            }
1995
            roots[i] = p - num;
1996
            maxDiff = std::max(maxDiff, cv::abs(num));
1997
        }
1998
        if( maxDiff <= 0 )
1999
            break;
2000
    }
2001

2002
    if( coeffs0.channels() == 1 )
2003
    {
2004
        const double verySmallEps = 1e-100;
2005
        for( i = 0; i < n; i++ )
2006
            if( fabs(roots[i].im) < verySmallEps )
2007
                roots[i].im = 0;
2008
    }
2009

2010
    for( ; n < n0; n++ )
2011
        roots[n+1] = roots[n];
2012

2013
    Mat(roots0.size(), CV_64FC2, roots).convertTo(roots0, roots0.type());
2014
    return maxDiff;
2015
}
2016

2017

2018
CV_IMPL int
2019
cvSolveCubic( const CvMat* coeffs, CvMat* roots )
2020
{
2021
    cv::Mat _coeffs = cv::cvarrToMat(coeffs), _roots = cv::cvarrToMat(roots), _roots0 = _roots;
2022
    int nroots = cv::solveCubic(_coeffs, _roots);
2023
    CV_Assert( _roots.data == _roots0.data ); // check that the array of roots was not reallocated
2024
    return nroots;
2025
}
2026

2027

2028
void cvSolvePoly(const CvMat* a, CvMat *r, int maxiter, int)
2029
{
2030
    cv::Mat _a = cv::cvarrToMat(a);
2031
    cv::Mat _r = cv::cvarrToMat(r);
2032
    cv::Mat _r0 = _r;
2033
    cv::solvePoly(_a, _r, maxiter);
2034
    CV_Assert( _r.data == _r0.data ); // check that the array of roots was not reallocated
2035
}
2036

2037

2038

2039
// Common constants for dispatched code
2040
namespace cv { namespace details {
2041

2042
#define EXPTAB_SCALE 6
2043
#define EXPTAB_MASK  ((1 << EXPTAB_SCALE) - 1)
2044

2045
#define EXPPOLY_32F_A0 .9670371139572337719125840413672004409288e-2
2046

2047
static const double CV_DECL_ALIGNED(64) expTab[EXPTAB_MASK + 1] = {
2048
    1.0 * EXPPOLY_32F_A0,
2049
    1.0108892860517004600204097905619 * EXPPOLY_32F_A0,
2050
    1.0218971486541166782344801347833 * EXPPOLY_32F_A0,
2051
    1.0330248790212284225001082839705 * EXPPOLY_32F_A0,
2052
    1.0442737824274138403219664787399 * EXPPOLY_32F_A0,
2053
    1.0556451783605571588083413251529 * EXPPOLY_32F_A0,
2054
    1.0671404006768236181695211209928 * EXPPOLY_32F_A0,
2055
    1.0787607977571197937406800374385 * EXPPOLY_32F_A0,
2056
    1.0905077326652576592070106557607 * EXPPOLY_32F_A0,
2057
    1.1023825833078409435564142094256 * EXPPOLY_32F_A0,
2058
    1.1143867425958925363088129569196 * EXPPOLY_32F_A0,
2059
    1.126521618608241899794798643787 * EXPPOLY_32F_A0,
2060
    1.1387886347566916537038302838415 * EXPPOLY_32F_A0,
2061
    1.151189229952982705817759635202 * EXPPOLY_32F_A0,
2062
    1.1637248587775775138135735990922 * EXPPOLY_32F_A0,
2063
    1.1763969916502812762846457284838 * EXPPOLY_32F_A0,
2064
    1.1892071150027210667174999705605 * EXPPOLY_32F_A0,
2065
    1.2021567314527031420963969574978 * EXPPOLY_32F_A0,
2066
    1.2152473599804688781165202513388 * EXPPOLY_32F_A0,
2067
    1.2284805361068700056940089577928 * EXPPOLY_32F_A0,
2068
    1.2418578120734840485936774687266 * EXPPOLY_32F_A0,
2069
    1.2553807570246910895793906574423 * EXPPOLY_32F_A0,
2070
    1.2690509571917332225544190810323 * EXPPOLY_32F_A0,
2071
    1.2828700160787782807266697810215 * EXPPOLY_32F_A0,
2072
    1.2968395546510096659337541177925 * EXPPOLY_32F_A0,
2073
    1.3109612115247643419229917863308 * EXPPOLY_32F_A0,
2074
    1.3252366431597412946295370954987 * EXPPOLY_32F_A0,
2075
    1.3396675240533030053600306697244 * EXPPOLY_32F_A0,
2076
    1.3542555469368927282980147401407 * EXPPOLY_32F_A0,
2077
    1.3690024229745906119296011329822 * EXPPOLY_32F_A0,
2078
    1.3839098819638319548726595272652 * EXPPOLY_32F_A0,
2079
    1.3989796725383111402095281367152 * EXPPOLY_32F_A0,
2080
    1.4142135623730950488016887242097 * EXPPOLY_32F_A0,
2081
    1.4296133383919700112350657782751 * EXPPOLY_32F_A0,
2082
    1.4451808069770466200370062414717 * EXPPOLY_32F_A0,
2083
    1.4609177941806469886513028903106 * EXPPOLY_32F_A0,
2084
    1.476826145939499311386907480374 * EXPPOLY_32F_A0,
2085
    1.4929077282912648492006435314867 * EXPPOLY_32F_A0,
2086
    1.5091644275934227397660195510332 * EXPPOLY_32F_A0,
2087
    1.5255981507445383068512536895169 * EXPPOLY_32F_A0,
2088
    1.5422108254079408236122918620907 * EXPPOLY_32F_A0,
2089
    1.5590044002378369670337280894749 * EXPPOLY_32F_A0,
2090
    1.5759808451078864864552701601819 * EXPPOLY_32F_A0,
2091
    1.5931421513422668979372486431191 * EXPPOLY_32F_A0,
2092
    1.6104903319492543081795206673574 * EXPPOLY_32F_A0,
2093
    1.628027421857347766848218522014 * EXPPOLY_32F_A0,
2094
    1.6457554781539648445187567247258 * EXPPOLY_32F_A0,
2095
    1.6636765803267364350463364569764 * EXPPOLY_32F_A0,
2096
    1.6817928305074290860622509524664 * EXPPOLY_32F_A0,
2097
    1.7001063537185234695013625734975 * EXPPOLY_32F_A0,
2098
    1.7186192981224779156293443764563 * EXPPOLY_32F_A0,
2099
    1.7373338352737062489942020818722 * EXPPOLY_32F_A0,
2100
    1.7562521603732994831121606193753 * EXPPOLY_32F_A0,
2101
    1.7753764925265212525505592001993 * EXPPOLY_32F_A0,
2102
    1.7947090750031071864277032421278 * EXPPOLY_32F_A0,
2103
    1.8142521755003987562498346003623 * EXPPOLY_32F_A0,
2104
    1.8340080864093424634870831895883 * EXPPOLY_32F_A0,
2105
    1.8539791250833855683924530703377 * EXPPOLY_32F_A0,
2106
    1.8741676341102999013299989499544 * EXPPOLY_32F_A0,
2107
    1.8945759815869656413402186534269 * EXPPOLY_32F_A0,
2108
    1.9152065613971472938726112702958 * EXPPOLY_32F_A0,
2109
    1.9360617934922944505980559045667 * EXPPOLY_32F_A0,
2110
    1.9571441241754002690183222516269 * EXPPOLY_32F_A0,
2111
    1.9784560263879509682582499181312 * EXPPOLY_32F_A0,
2112
};
2113

2114
const double* getExpTab64f()
2115
{
2116
    return expTab;
2117
}
2118

2119
const float* getExpTab32f()
2120
{
2121
    static float CV_DECL_ALIGNED(64) expTab_f[EXPTAB_MASK+1];
2122
    static volatile bool expTab_f_initialized = false;
2123
    if (!expTab_f_initialized)
2124
    {
2125
        for( int j = 0; j <= EXPTAB_MASK; j++ )
2126
            expTab_f[j] = (float)expTab[j];
2127
        expTab_f_initialized = true;
2128
    }
2129
    return expTab_f;
2130
}
2131

2132

2133

2134
#define LOGTAB_SCALE        8
2135
#define LOGTAB_MASK         ((1 << LOGTAB_SCALE) - 1)
2136

2137
static const double CV_DECL_ALIGNED(64) logTab[(LOGTAB_MASK+1)*2] = {
2138
    0.0000000000000000000000000000000000000000,    1.000000000000000000000000000000000000000,
2139
    .00389864041565732288852075271279318258166,    .9961089494163424124513618677042801556420,
2140
    .00778214044205494809292034119607706088573,    .9922480620155038759689922480620155038760,
2141
    .01165061721997527263705585198749759001657,    .9884169884169884169884169884169884169884,
2142
    .01550418653596525274396267235488267033361,    .9846153846153846153846153846153846153846,
2143
    .01934296284313093139406447562578250654042,    .9808429118773946360153256704980842911877,
2144
    .02316705928153437593630670221500622574241,    .9770992366412213740458015267175572519084,
2145
    .02697658769820207233514075539915211265906,    .9733840304182509505703422053231939163498,
2146
    .03077165866675368732785500469617545604706,    .9696969696969696969696969696969696969697,
2147
    .03455238150665972812758397481047722976656,    .9660377358490566037735849056603773584906,
2148
    .03831886430213659461285757856785494368522,    .9624060150375939849624060150375939849624,
2149
    .04207121392068705056921373852674150839447,    .9588014981273408239700374531835205992509,
2150
    .04580953603129420126371940114040626212953,    .9552238805970149253731343283582089552239,
2151
    .04953393512227662748292900118940451648088,    .9516728624535315985130111524163568773234,
2152
    .05324451451881227759255210685296333394944,    .9481481481481481481481481481481481481481,
2153
    .05694137640013842427411105973078520037234,    .9446494464944649446494464944649446494465,
2154
    .06062462181643483993820353816772694699466,    .9411764705882352941176470588235294117647,
2155
    .06429435070539725460836422143984236754475,    .9377289377289377289377289377289377289377,
2156
    .06795066190850773679699159401934593915938,    .9343065693430656934306569343065693430657,
2157
    .07159365318700880442825962290953611955044,    .9309090909090909090909090909090909090909,
2158
    .07522342123758751775142172846244648098944,    .9275362318840579710144927536231884057971,
2159
    .07884006170777602129362549021607264876369,    .9241877256317689530685920577617328519856,
2160
    .08244366921107458556772229485432035289706,    .9208633093525179856115107913669064748201,
2161
    .08603433734180314373940490213499288074675,    .9175627240143369175627240143369175627240,
2162
    .08961215868968712416897659522874164395031,    .9142857142857142857142857142857142857143,
2163
    .09317722485418328259854092721070628613231,    .9110320284697508896797153024911032028470,
2164
    .09672962645855109897752299730200320482256,    .9078014184397163120567375886524822695035,
2165
    .10026945316367513738597949668474029749630,    .9045936395759717314487632508833922261484,
2166
    .10379679368164355934833764649738441221420,    .9014084507042253521126760563380281690141,
2167
    .10731173578908805021914218968959175981580,    .8982456140350877192982456140350877192982,
2168
    .11081436634029011301105782649756292812530,    .8951048951048951048951048951048951048951,
2169
    .11430477128005862852422325204315711744130,    .8919860627177700348432055749128919860627,
2170
    .11778303565638344185817487641543266363440,    .8888888888888888888888888888888888888889,
2171
    .12124924363286967987640707633545389398930,    .8858131487889273356401384083044982698962,
2172
    .12470347850095722663787967121606925502420,    .8827586206896551724137931034482758620690,
2173
    .12814582269193003360996385708858724683530,    .8797250859106529209621993127147766323024,
2174
    .13157635778871926146571524895989568904040,    .8767123287671232876712328767123287671233,
2175
    .13499516453750481925766280255629681050780,    .8737201365187713310580204778156996587031,
2176
    .13840232285911913123754857224412262439730,    .8707482993197278911564625850340136054422,
2177
    .14179791186025733629172407290752744302150,    .8677966101694915254237288135593220338983,
2178
    .14518200984449788903951628071808954700830,    .8648648648648648648648648648648648648649,
2179
    .14855469432313711530824207329715136438610,    .8619528619528619528619528619528619528620,
2180
    .15191604202584196858794030049466527998450,    .8590604026845637583892617449664429530201,
2181
    .15526612891112392955683674244937719777230,    .8561872909698996655518394648829431438127,
2182
    .15860503017663857283636730244325008243330,    .8533333333333333333333333333333333333333,
2183
    .16193282026931324346641360989451641216880,    .8504983388704318936877076411960132890365,
2184
    .16524957289530714521497145597095368430010,    .8476821192052980132450331125827814569536,
2185
    .16855536102980664403538924034364754334090,    .8448844884488448844884488448844884488449,
2186
    .17185025692665920060697715143760433420540,    .8421052631578947368421052631578947368421,
2187
    .17513433212784912385018287750426679849630,    .8393442622950819672131147540983606557377,
2188
    .17840765747281828179637841458315961062910,    .8366013071895424836601307189542483660131,
2189
    .18167030310763465639212199675966985523700,    .8338762214983713355048859934853420195440,
2190
    .18492233849401198964024217730184318497780,    .8311688311688311688311688311688311688312,
2191
    .18816383241818296356839823602058459073300,    .8284789644012944983818770226537216828479,
2192
    .19139485299962943898322009772527962923050,    .8258064516129032258064516129032258064516,
2193
    .19461546769967164038916962454095482826240,    .8231511254019292604501607717041800643087,
2194
    .19782574332991986754137769821682013571260,    .8205128205128205128205128205128205128205,
2195
    .20102574606059073203390141770796617493040,    .8178913738019169329073482428115015974441,
2196
    .20421554142869088876999228432396193966280,    .8152866242038216560509554140127388535032,
2197
    .20739519434607056602715147164417430758480,    .8126984126984126984126984126984126984127,
2198
    .21056476910734961416338251183333341032260,    .8101265822784810126582278481012658227848,
2199
    .21372432939771812687723695489694364368910,    .8075709779179810725552050473186119873817,
2200
    .21687393830061435506806333251006435602900,    .8050314465408805031446540880503144654088,
2201
    .22001365830528207823135744547471404075630,    .8025078369905956112852664576802507836991,
2202
    .22314355131420973710199007200571941211830,    .8000000000000000000000000000000000000000,
2203
    .22626367865045338145790765338460914790630,    .7975077881619937694704049844236760124611,
2204
    .22937410106484582006380890106811420992010,    .7950310559006211180124223602484472049689,
2205
    .23247487874309405442296849741978803649550,    .7925696594427244582043343653250773993808,
2206
    .23556607131276688371634975283086532726890,    .7901234567901234567901234567901234567901,
2207
    .23864773785017498464178231643018079921600,    .7876923076923076923076923076923076923077,
2208
    .24171993688714515924331749374687206000090,    .7852760736196319018404907975460122699387,
2209
    .24478272641769091566565919038112042471760,    .7828746177370030581039755351681957186544,
2210
    .24783616390458124145723672882013488560910,    .7804878048780487804878048780487804878049,
2211
    .25088030628580937353433455427875742316250,    .7781155015197568389057750759878419452888,
2212
    .25391520998096339667426946107298135757450,    .7757575757575757575757575757575757575758,
2213
    .25694093089750041913887912414793390780680,    .7734138972809667673716012084592145015106,
2214
    .25995752443692604627401010475296061486000,    .7710843373493975903614457831325301204819,
2215
    .26296504550088134477547896494797896593800,    .7687687687687687687687687687687687687688,
2216
    .26596354849713793599974565040611196309330,    .7664670658682634730538922155688622754491,
2217
    .26895308734550393836570947314612567424780,    .7641791044776119402985074626865671641791,
2218
    .27193371548364175804834985683555714786050,    .7619047619047619047619047619047619047619,
2219
    .27490548587279922676529508862586226314300,    .7596439169139465875370919881305637982196,
2220
    .27786845100345625159121709657483734190480,    .7573964497041420118343195266272189349112,
2221
    .28082266290088775395616949026589281857030,    .7551622418879056047197640117994100294985,
2222
    .28376817313064456316240580235898960381750,    .7529411764705882352941176470588235294118,
2223
    .28670503280395426282112225635501090437180,    .7507331378299120234604105571847507331378,
2224
    .28963329258304265634293983566749375313530,    .7485380116959064327485380116959064327485,
2225
    .29255300268637740579436012922087684273730,    .7463556851311953352769679300291545189504,
2226
    .29546421289383584252163927885703742504130,    .7441860465116279069767441860465116279070,
2227
    .29836697255179722709783618483925238251680,    .7420289855072463768115942028985507246377,
2228
    .30126133057816173455023545102449133992200,    .7398843930635838150289017341040462427746,
2229
    .30414733546729666446850615102448500692850,    .7377521613832853025936599423631123919308,
2230
    .30702503529491181888388950937951449304830,    .7356321839080459770114942528735632183908,
2231
    .30989447772286465854207904158101882785550,    .7335243553008595988538681948424068767908,
2232
    .31275571000389684739317885942000430077330,    .7314285714285714285714285714285714285714,
2233
    .31560877898630329552176476681779604405180,    .7293447293447293447293447293447293447293,
2234
    .31845373111853458869546784626436419785030,    .7272727272727272727272727272727272727273,
2235
    .32129061245373424782201254856772720813750,    .7252124645892351274787535410764872521246,
2236
    .32411946865421192853773391107097268104550,    .7231638418079096045197740112994350282486,
2237
    .32694034499585328257253991068864706903700,    .7211267605633802816901408450704225352113,
2238
    .32975328637246797969240219572384376078850,    .7191011235955056179775280898876404494382,
2239
    .33255833730007655635318997155991382896900,    .7170868347338935574229691876750700280112,
2240
    .33535554192113781191153520921943709254280,    .7150837988826815642458100558659217877095,
2241
    .33814494400871636381467055798566434532400,    .7130919220055710306406685236768802228412,
2242
    .34092658697059319283795275623560883104800,    .7111111111111111111111111111111111111111,
2243
    .34370051385331840121395430287520866841080,    .7091412742382271468144044321329639889197,
2244
    .34646676734620857063262633346312213689100,    .7071823204419889502762430939226519337017,
2245
    .34922538978528827602332285096053965389730,    .7052341597796143250688705234159779614325,
2246
    .35197642315717814209818925519357435405250,    .7032967032967032967032967032967032967033,
2247
    .35471990910292899856770532096561510115850,    .7013698630136986301369863013698630136986,
2248
    .35745588892180374385176833129662554711100,    .6994535519125683060109289617486338797814,
2249
    .36018440357500774995358483465679455548530,    .6975476839237057220708446866485013623978,
2250
    .36290549368936841911903457003063522279280,    .6956521739130434782608695652173913043478,
2251
    .36561919956096466943762379742111079394830,    .6937669376693766937669376693766937669377,
2252
    .36832556115870762614150635272380895912650,    .6918918918918918918918918918918918918919,
2253
    .37102461812787262962487488948681857436900,    .6900269541778975741239892183288409703504,
2254
    .37371640979358405898480555151763837784530,    .6881720430107526881720430107526881720430,
2255
    .37640097516425302659470730759494472295050,    .6863270777479892761394101876675603217158,
2256
    .37907835293496944251145919224654790014030,    .6844919786096256684491978609625668449198,
2257
    .38174858149084833769393299007788300514230,    .6826666666666666666666666666666666666667,
2258
    .38441169891033200034513583887019194662580,    .6808510638297872340425531914893617021277,
2259
    .38706774296844825844488013899535872042180,    .6790450928381962864721485411140583554377,
2260
    .38971675114002518602873692543653305619950,    .6772486772486772486772486772486772486772,
2261
    .39235876060286384303665840889152605086580,    .6754617414248021108179419525065963060686,
2262
    .39499380824086893770896722344332374632350,    .6736842105263157894736842105263157894737,
2263
    .39762193064713846624158577469643205404280,    .6719160104986876640419947506561679790026,
2264
    .40024316412701266276741307592601515352730,    .6701570680628272251308900523560209424084,
2265
    .40285754470108348090917615991202183067800,    .6684073107049608355091383812010443864230,
2266
    .40546510810816432934799991016916465014230,    .6666666666666666666666666666666666666667,
2267
    .40806588980822172674223224930756259709600,    .6649350649350649350649350649350649350649,
2268
    .41065992498526837639616360320360399782650,    .6632124352331606217616580310880829015544,
2269
    .41324724855021932601317757871584035456180,    .6614987080103359173126614987080103359173,
2270
    .41582789514371093497757669865677598863850,    .6597938144329896907216494845360824742268,
2271
    .41840189913888381489925905043492093682300,    .6580976863753213367609254498714652956298,
2272
    .42096929464412963239894338585145305842150,    .6564102564102564102564102564102564102564,
2273
    .42353011550580327293502591601281892508280,    .6547314578005115089514066496163682864450,
2274
    .42608439531090003260516141381231136620050,    .6530612244897959183673469387755102040816,
2275
    .42863216738969872610098832410585600882780,    .6513994910941475826972010178117048346056,
2276
    .43117346481837132143866142541810404509300,    .6497461928934010152284263959390862944162,
2277
    .43370832042155937902094819946796633303180,    .6481012658227848101265822784810126582278,
2278
    .43623676677491801667585491486534010618930,    .6464646464646464646464646464646464646465,
2279
    .43875883620762790027214350629947148263450,    .6448362720403022670025188916876574307305,
2280
    .44127456080487520440058801796112675219780,    .6432160804020100502512562814070351758794,
2281
    .44378397241030093089975139264424797147500,    .6416040100250626566416040100250626566416,
2282
    .44628710262841947420398014401143882423650,    .6400000000000000000000000000000000000000,
2283
    .44878398282700665555822183705458883196130,    .6384039900249376558603491271820448877805,
2284
    .45127464413945855836729492693848442286250,    .6368159203980099502487562189054726368159,
2285
    .45375911746712049854579618113348260521900,    .6352357320099255583126550868486352357320,
2286
    .45623743348158757315857769754074979573500,    .6336633663366336633663366336633663366337,
2287
    .45870962262697662081833982483658473938700,    .6320987654320987654320987654320987654321,
2288
    .46117571512217014895185229761409573256980,    .6305418719211822660098522167487684729064,
2289
    .46363574096303250549055974261136725544930,    .6289926289926289926289926289926289926290,
2290
    .46608972992459918316399125615134835243230,    .6274509803921568627450980392156862745098,
2291
    .46853771156323925639597405279346276074650,    .6259168704156479217603911980440097799511,
2292
    .47097971521879100631480241645476780831830,    .6243902439024390243902439024390243902439,
2293
    .47341577001667212165614273544633761048330,    .6228710462287104622871046228710462287105,
2294
    .47584590486996386493601107758877333253630,    .6213592233009708737864077669902912621359,
2295
    .47827014848147025860569669930555392056700,    .6198547215496368038740920096852300242131,
2296
    .48068852934575190261057286988943815231330,    .6183574879227053140096618357487922705314,
2297
    .48310107575113581113157579238759353756900,    .6168674698795180722891566265060240963855,
2298
    .48550781578170076890899053978500887751580,    .6153846153846153846153846153846153846154,
2299
    .48790877731923892879351001283794175833480,    .6139088729016786570743405275779376498801,
2300
    .49030398804519381705802061333088204264650,    .6124401913875598086124401913875598086124,
2301
    .49269347544257524607047571407747454941280,    .6109785202863961813842482100238663484487,
2302
    .49507726679785146739476431321236304938800,    .6095238095238095238095238095238095238095,
2303
    .49745538920281889838648226032091770321130,    .6080760095011876484560570071258907363420,
2304
    .49982786955644931126130359189119189977650,    .6066350710900473933649289099526066350711,
2305
    .50219473456671548383667413872899487614650,    .6052009456264775413711583924349881796690,
2306
    .50455601075239520092452494282042607665050,    .6037735849056603773584905660377358490566,
2307
    .50691172444485432801997148999362252652650,    .6023529411764705882352941176470588235294,
2308
    .50926190178980790257412536448100581765150,    .6009389671361502347417840375586854460094,
2309
    .51160656874906207391973111953120678663250,    .5995316159250585480093676814988290398126,
2310
    .51394575110223428282552049495279788970950,    .5981308411214953271028037383177570093458,
2311
    .51627947444845445623684554448118433356300,    .5967365967365967365967365967365967365967,
2312
    .51860776420804555186805373523384332656850,    .5953488372093023255813953488372093023256,
2313
    .52093064562418522900344441950437612831600,    .5939675174013921113689095127610208816705,
2314
    .52324814376454775732838697877014055848100,    .5925925925925925925925925925925925925926,
2315
    .52556028352292727401362526507000438869000,    .5912240184757505773672055427251732101617,
2316
    .52786708962084227803046587723656557500350,    .5898617511520737327188940092165898617512,
2317
    .53016858660912158374145519701414741575700,    .5885057471264367816091954022988505747126,
2318
    .53246479886947173376654518506256863474850,    .5871559633027522935779816513761467889908,
2319
    .53475575061602764748158733709715306758900,    .5858123569794050343249427917620137299771,
2320
    .53704146589688361856929077475797384977350,    .5844748858447488584474885844748858447489,
2321
    .53932196859560876944783558428753167390800,    .5831435079726651480637813211845102505695,
2322
    .54159728243274429804188230264117009937750,    .5818181818181818181818181818181818181818,
2323
    .54386743096728351609669971367111429572100,    .5804988662131519274376417233560090702948,
2324
    .54613243759813556721383065450936555862450,    .5791855203619909502262443438914027149321,
2325
    .54839232556557315767520321969641372561450,    .5778781038374717832957110609480812641084,
2326
    .55064711795266219063194057525834068655950,    .5765765765765765765765765765765765765766,
2327
    .55289683768667763352766542084282264113450,    .5752808988764044943820224719101123595506,
2328
    .55514150754050151093110798683483153581600,    .5739910313901345291479820627802690582960,
2329
    .55738115013400635344709144192165695130850,    .5727069351230425055928411633109619686801,
2330
    .55961578793542265941596269840374588966350,    .5714285714285714285714285714285714285714,
2331
    .56184544326269181269140062795486301183700,    .5701559020044543429844097995545657015590,
2332
    .56407013828480290218436721261241473257550,    .5688888888888888888888888888888888888889,
2333
    .56628989502311577464155334382667206227800,    .5676274944567627494456762749445676274945,
2334
    .56850473535266865532378233183408156037350,    .5663716814159292035398230088495575221239,
2335
    .57071468100347144680739575051120482385150,    .5651214128035320088300220750551876379691,
2336
    .57291975356178548306473885531886480748650,    .5638766519823788546255506607929515418502,
2337
    .57511997447138785144460371157038025558000,    .5626373626373626373626373626373626373626,
2338
    .57731536503482350219940144597785547375700,    .5614035087719298245614035087719298245614,
2339
    .57950594641464214795689713355386629700650,    .5601750547045951859956236323851203501094,
2340
    .58169173963462239562716149521293118596100,    .5589519650655021834061135371179039301310,
2341
    .58387276558098266665552955601015128195300,    .5577342047930283224400871459694989106754,
2342
    .58604904500357812846544902640744112432000,    .5565217391304347826086956521739130434783,
2343
    .58822059851708596855957011939608491957200,    .5553145336225596529284164859002169197397,
2344
    .59038744660217634674381770309992134571100,    .5541125541125541125541125541125541125541,
2345
    .59254960960667157898740242671919986605650,    .5529157667386609071274298056155507559395,
2346
    .59470710774669277576265358220553025603300,    .5517241379310344827586206896551724137931,
2347
    .59685996110779382384237123915227130055450,    .5505376344086021505376344086021505376344,
2348
    .59900818964608337768851242799428291618800,    .5493562231759656652360515021459227467811,
2349
    .60115181318933474940990890900138765573500,    .5481798715203426124197002141327623126338,
2350
    .60329085143808425240052883964381180703650,    .5470085470085470085470085470085470085470,
2351
    .60542532396671688843525771517306566238400,    .5458422174840085287846481876332622601279,
2352
    .60755525022454170969155029524699784815300,    .5446808510638297872340425531914893617021,
2353
    .60968064953685519036241657886421307921400,    .5435244161358811040339702760084925690021,
2354
    .61180154110599282990534675263916142284850,    .5423728813559322033898305084745762711864,
2355
    .61391794401237043121710712512140162289150,    .5412262156448202959830866807610993657505,
2356
    .61602987721551394351138242200249806046500,    .5400843881856540084388185654008438818565,
2357
    .61813735955507864705538167982012964785100,    .5389473684210526315789473684210526315789,
2358
    .62024040975185745772080281312810257077200,    .5378151260504201680672268907563025210084,
2359
    .62233904640877868441606324267922900617100,    .5366876310272536687631027253668763102725,
2360
    .62443328801189346144440150965237990021700,    .5355648535564853556485355648535564853556,
2361
    .62652315293135274476554741340805776417250,    .5344467640918580375782881002087682672234,
2362
    .62860865942237409420556559780379757285100,    .5333333333333333333333333333333333333333,
2363
    .63068982562619868570408243613201193511500,    .5322245322245322245322245322245322245322,
2364
    .63276666957103777644277897707070223987100,    .5311203319502074688796680497925311203320,
2365
    .63483920917301017716738442686619237065300,    .5300207039337474120082815734989648033126,
2366
    .63690746223706917739093569252872839570050,    .5289256198347107438016528925619834710744,
2367
    .63897144645792069983514238629140891134750,    .5278350515463917525773195876288659793814,
2368
    .64103117942093124081992527862894348800200,    .5267489711934156378600823045267489711934,
2369
    .64308667860302726193566513757104985415950,    .5256673511293634496919917864476386036961,
2370
    .64513796137358470073053240412264131009600,    .5245901639344262295081967213114754098361,
2371
    .64718504499530948859131740391603671014300,    .5235173824130879345603271983640081799591,
2372
    .64922794662510974195157587018911726772800,    .5224489795918367346938775510204081632653,
2373
    .65126668331495807251485530287027359008800,    .5213849287169042769857433808553971486762,
2374
    .65330127201274557080523663898929953575150,    .5203252032520325203252032520325203252033,
2375
    .65533172956312757406749369692988693714150,    .5192697768762677484787018255578093306288,
2376
    .65735807270835999727154330685152672231200,    .5182186234817813765182186234817813765182,
2377
    .65938031808912778153342060249997302889800,    .5171717171717171717171717171717171717172,
2378
    .66139848224536490484126716182800009846700,    .5161290322580645161290322580645161290323,
2379
    .66341258161706617713093692145776003599150,    .5150905432595573440643863179074446680080,
2380
    .66542263254509037562201001492212526500250,    .5140562248995983935742971887550200803213,
2381
    .66742865127195616370414654738851822912700,    .5130260521042084168336673346693386773547,
2382
    .66943065394262923906154583164607174694550,    .5120000000000000000000000000000000000000,
2383
    .67142865660530226534774556057527661323550,    .5109780439121756487025948103792415169661,
2384
    .67342267521216669923234121597488410770900,    .5099601593625498007968127490039840637450,
2385
    .67541272562017662384192817626171745359900,    .5089463220675944333996023856858846918489,
2386
    .67739882359180603188519853574689477682100,    .5079365079365079365079365079365079365079,
2387
    .67938098479579733801614338517538271844400,    .5069306930693069306930693069306930693069,
2388
    .68135922480790300781450241629499942064300,    .5059288537549407114624505928853754940711,
2389
    .68333355911162063645036823800182901322850,    .5049309664694280078895463510848126232742,
2390
    .68530400309891936760919861626462079584600,    .5039370078740157480314960629921259842520,
2391
    .68727057207096020619019327568821609020250,    .5029469548133595284872298624754420432220,
2392
    .68923328123880889251040571252815425395950,    .5019607843137254901960784313725490196078,
2393
    .69314718055994530941723212145818, 5.0e-01,
2394
};
2395

2396
const double* getLogTab64f()
2397
{
2398
    return logTab;
2399
}
2400

2401
const float* getLogTab32f()
2402
{
2403
    static float CV_DECL_ALIGNED(64) logTab_f[(LOGTAB_MASK+1)*2];
2404
    static volatile bool logTab_f_initialized = false;
2405
    if (!logTab_f_initialized)
2406
    {
2407
        for (int j = 0; j < (LOGTAB_MASK+1)*2; j++)
2408
            logTab_f[j] = (float)logTab[j];
2409
        logTab_f_initialized = true;
2410
    }
2411
    return logTab_f;
2412
}
2413

2414
}} // namespace
2415

2416
/* End of file. */
2417

2418
Product

Resources

Company
