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 <atomic>
47 #include <limits>
48 #include <iostream>
49 #include "mathfuncs.hpp"
50
51 namespace cv
52 {
53
54 typedef void (*MathFunc)(const void* src, void* dst, int len);
55
56 #ifdef HAVE_OPENCL
57
58 enum { OCL_OP_LOG=0, OCL_OP_EXP=1, OCL_OP_MAG=2, OCL_OP_PHASE_DEGREES=3, OCL_OP_PHASE_RADIANS=4 };
59
60 static const char* oclop2str[] = { "OP_LOG", "OP_EXP", "OP_MAG", "OP_PHASE_DEGREES", "OP_PHASE_RADIANS", 0 };
61
ocl_math_op(InputArray _src1,InputArray _src2,OutputArray _dst,int oclop)62 static bool ocl_math_op(InputArray _src1, InputArray _src2, OutputArray _dst, int oclop)
63 {
64 int type = _src1.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
65 int kercn = oclop == OCL_OP_PHASE_DEGREES ||
66 oclop == OCL_OP_PHASE_RADIANS ? 1 : ocl::predictOptimalVectorWidth(_src1, _src2, _dst);
67
68 const ocl::Device d = ocl::Device::getDefault();
69 bool double_support = d.doubleFPConfig() > 0;
70 if (!double_support && depth == CV_64F)
71 return false;
72 int rowsPerWI = d.isIntel() ? 4 : 1;
73
74 ocl::Kernel k("KF", ocl::core::arithm_oclsrc,
75 format("-D %s -D %s -D dstT=%s -D DEPTH_dst=%d -D rowsPerWI=%d%s", _src2.empty() ? "UNARY_OP" : "BINARY_OP",
76 oclop2str[oclop], ocl::typeToStr(CV_MAKE_TYPE(depth, kercn)), depth, rowsPerWI,
77 double_support ? " -D DOUBLE_SUPPORT" : ""));
78 if (k.empty())
79 return false;
80
81 UMat src1 = _src1.getUMat(), src2 = _src2.getUMat();
82 _dst.create(src1.size(), type);
83 UMat dst = _dst.getUMat();
84
85 ocl::KernelArg src1arg = ocl::KernelArg::ReadOnlyNoSize(src1),
86 src2arg = ocl::KernelArg::ReadOnlyNoSize(src2),
87 dstarg = ocl::KernelArg::WriteOnly(dst, cn, kercn);
88
89 if (src2.empty())
90 k.args(src1arg, dstarg);
91 else
92 k.args(src1arg, src2arg, dstarg);
93
94 size_t globalsize[] = { (size_t)src1.cols * cn / kercn, ((size_t)src1.rows + rowsPerWI - 1) / rowsPerWI };
95 return k.run(2, globalsize, 0, false);
96 }
97
98 #endif
99
100 /* ************************************************************************** *\
101 Fast cube root by Ken Turkowski
102 (http://www.worldserver.com/turk/computergraphics/papers.html)
103 \* ************************************************************************** */
cubeRoot(float value)104 float cubeRoot( float value )
105 {
106 CV_INSTRUMENT_REGION();
107
108 float fr;
109 Cv32suf v, m;
110 int ix, s;
111 int ex, shx;
112
113 v.f = value;
114 ix = v.i & 0x7fffffff;
115 s = v.i & 0x80000000;
116 ex = (ix >> 23) - 127;
117 shx = ex % 3;
118 shx -= shx >= 0 ? 3 : 0;
119 ex = (ex - shx) / 3; /* exponent of cube root */
120 v.i = (ix & ((1<<23)-1)) | ((shx + 127)<<23);
121 fr = v.f;
122
123 /* 0.125 <= fr < 1.0 */
124 /* Use quartic rational polynomial with error < 2^(-24) */
125 fr = (float)(((((45.2548339756803022511987494 * fr +
126 192.2798368355061050458134625) * fr +
127 119.1654824285581628956914143) * fr +
128 13.43250139086239872172837314) * fr +
129 0.1636161226585754240958355063)/
130 ((((14.80884093219134573786480845 * fr +
131 151.9714051044435648658557668) * fr +
132 168.5254414101568283957668343) * fr +
133 33.9905941350215598754191872) * fr +
134 1.0));
135
136 /* fr *= 2^ex * sign */
137 m.f = value;
138 v.f = fr;
139 v.i = (v.i + (ex << 23) + s) & (m.i*2 != 0 ? -1 : 0);
140 return v.f;
141 }
142
143 /****************************************************************************************\
144 * Cartezian -> Polar *
145 \****************************************************************************************/
146
magnitude(InputArray src1,InputArray src2,OutputArray dst)147 void magnitude( InputArray src1, InputArray src2, OutputArray dst )
148 {
149 CV_INSTRUMENT_REGION();
150
151 int type = src1.type(), depth = src1.depth(), cn = src1.channels();
152 CV_Assert( src1.size() == src2.size() && type == src2.type() && (depth == CV_32F || depth == CV_64F));
153
154 CV_OCL_RUN(dst.isUMat() && src1.dims() <= 2 && src2.dims() <= 2,
155 ocl_math_op(src1, src2, dst, OCL_OP_MAG))
156
157 Mat X = src1.getMat(), Y = src2.getMat();
158 dst.create(X.dims, X.size, X.type());
159 Mat Mag = dst.getMat();
160
161 const Mat* arrays[] = {&X, &Y, &Mag, 0};
162 uchar* ptrs[3] = {};
163 NAryMatIterator it(arrays, ptrs);
164 int len = (int)it.size*cn;
165
166 for( size_t i = 0; i < it.nplanes; i++, ++it )
167 {
168 if( depth == CV_32F )
169 {
170 const float *x = (const float*)ptrs[0], *y = (const float*)ptrs[1];
171 float *mag = (float*)ptrs[2];
172 hal::magnitude32f( x, y, mag, len );
173 }
174 else
175 {
176 const double *x = (const double*)ptrs[0], *y = (const double*)ptrs[1];
177 double *mag = (double*)ptrs[2];
178 hal::magnitude64f( x, y, mag, len );
179 }
180 }
181 }
182
phase(InputArray src1,InputArray src2,OutputArray dst,bool angleInDegrees)183 void phase( InputArray src1, InputArray src2, OutputArray dst, bool angleInDegrees )
184 {
185 CV_INSTRUMENT_REGION();
186
187 int type = src1.type(), depth = src1.depth(), cn = src1.channels();
188 CV_Assert( src1.size() == src2.size() && type == src2.type() && (depth == CV_32F || depth == CV_64F));
189
190 CV_OCL_RUN(dst.isUMat() && src1.dims() <= 2 && src2.dims() <= 2,
191 ocl_math_op(src1, src2, dst, angleInDegrees ? OCL_OP_PHASE_DEGREES : OCL_OP_PHASE_RADIANS))
192
193 Mat X = src1.getMat(), Y = src2.getMat();
194 dst.create( X.dims, X.size, type );
195 Mat Angle = dst.getMat();
196
197 const Mat* arrays[] = {&X, &Y, &Angle, 0};
198 uchar* ptrs[3] = {};
199 NAryMatIterator it(arrays, ptrs);
200 int j, total = (int)(it.size*cn), blockSize = total;
201 size_t esz1 = X.elemSize1();
202 for( size_t i = 0; i < it.nplanes; i++, ++it )
203 {
204 for( j = 0; j < total; j += blockSize )
205 {
206 int len = std::min(total - j, blockSize);
207 if( depth == CV_32F )
208 {
209 const float *x = (const float*)ptrs[0], *y = (const float*)ptrs[1];
210 float *angle = (float*)ptrs[2];
211 hal::fastAtan32f( y, x, angle, len, angleInDegrees );
212 }
213 else
214 {
215 const double *x = (const double*)ptrs[0], *y = (const double*)ptrs[1];
216 double *angle = (double*)ptrs[2];
217 hal::fastAtan64f(y, x, angle, len, angleInDegrees);
218 }
219 ptrs[0] += len*esz1;
220 ptrs[1] += len*esz1;
221 ptrs[2] += len*esz1;
222 }
223 }
224 }
225
226 #ifdef HAVE_OPENCL
227
ocl_cartToPolar(InputArray _src1,InputArray _src2,OutputArray _dst1,OutputArray _dst2,bool angleInDegrees)228 static bool ocl_cartToPolar( InputArray _src1, InputArray _src2,
229 OutputArray _dst1, OutputArray _dst2, bool angleInDegrees )
230 {
231 const ocl::Device & d = ocl::Device::getDefault();
232 int type = _src1.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type),
233 rowsPerWI = d.isIntel() ? 4 : 1;
234 bool doubleSupport = d.doubleFPConfig() > 0;
235
236 if ( !(_src1.dims() <= 2 && _src2.dims() <= 2 &&
237 (depth == CV_32F || depth == CV_64F) && type == _src2.type()) ||
238 (depth == CV_64F && !doubleSupport) )
239 return false;
240
241 ocl::Kernel k("KF", ocl::core::arithm_oclsrc,
242 format("-D BINARY_OP -D dstT=%s -D DEPTH_dst=%d -D rowsPerWI=%d -D OP_CTP_%s%s",
243 ocl::typeToStr(CV_MAKE_TYPE(depth, 1)), depth,
244 rowsPerWI, angleInDegrees ? "AD" : "AR",
245 doubleSupport ? " -D DOUBLE_SUPPORT" : ""));
246 if (k.empty())
247 return false;
248
249 UMat src1 = _src1.getUMat(), src2 = _src2.getUMat();
250 Size size = src1.size();
251 CV_Assert( size == src2.size() );
252
253 _dst1.create(size, type);
254 _dst2.create(size, type);
255 UMat dst1 = _dst1.getUMat(), dst2 = _dst2.getUMat();
256
257 k.args(ocl::KernelArg::ReadOnlyNoSize(src1),
258 ocl::KernelArg::ReadOnlyNoSize(src2),
259 ocl::KernelArg::WriteOnly(dst1, cn),
260 ocl::KernelArg::WriteOnlyNoSize(dst2));
261
262 size_t globalsize[2] = { (size_t)dst1.cols * cn, ((size_t)dst1.rows + rowsPerWI - 1) / rowsPerWI };
263 return k.run(2, globalsize, NULL, false);
264 }
265
266 #endif
267
cartToPolar(InputArray src1,InputArray src2,OutputArray dst1,OutputArray dst2,bool angleInDegrees)268 void cartToPolar( InputArray src1, InputArray src2,
269 OutputArray dst1, OutputArray dst2, bool angleInDegrees )
270 {
271 CV_INSTRUMENT_REGION();
272
273 CV_OCL_RUN(dst1.isUMat() && dst2.isUMat(),
274 ocl_cartToPolar(src1, src2, dst1, dst2, angleInDegrees))
275
276 Mat X = src1.getMat(), Y = src2.getMat();
277 int type = X.type(), depth = X.depth(), cn = X.channels();
278 CV_Assert( X.size == Y.size && type == Y.type() && (depth == CV_32F || depth == CV_64F));
279 dst1.create( X.dims, X.size, type );
280 dst2.create( X.dims, X.size, type );
281 Mat Mag = dst1.getMat(), Angle = dst2.getMat();
282
283 const Mat* arrays[] = {&X, &Y, &Mag, &Angle, 0};
284 uchar* ptrs[4] = {};
285 NAryMatIterator it(arrays, ptrs);
286 int j, total = (int)(it.size*cn), blockSize = std::min(total, ((BLOCK_SIZE+cn-1)/cn)*cn);
287 size_t esz1 = X.elemSize1();
288
289 for( size_t i = 0; i < it.nplanes; i++, ++it )
290 {
291 for( j = 0; j < total; j += blockSize )
292 {
293 int len = std::min(total - j, blockSize);
294 if( depth == CV_32F )
295 {
296 const float *x = (const float*)ptrs[0], *y = (const float*)ptrs[1];
297 float *mag = (float*)ptrs[2], *angle = (float*)ptrs[3];
298 hal::magnitude32f( x, y, mag, len );
299 hal::fastAtan32f( y, x, angle, len, angleInDegrees );
300 }
301 else
302 {
303 const double *x = (const double*)ptrs[0], *y = (const double*)ptrs[1];
304 double *angle = (double*)ptrs[3];
305 hal::magnitude64f(x, y, (double*)ptrs[2], len);
306 hal::fastAtan64f(y, x, angle, len, angleInDegrees);
307 }
308 ptrs[0] += len*esz1;
309 ptrs[1] += len*esz1;
310 ptrs[2] += len*esz1;
311 ptrs[3] += len*esz1;
312 }
313 }
314 }
315
316
317 /****************************************************************************************\
318 * Polar -> Cartezian *
319 \****************************************************************************************/
320
SinCos_32f(const float * angle,float * sinval,float * cosval,int len,int angle_in_degrees)321 static void SinCos_32f( const float *angle, float *sinval, float* cosval,
322 int len, int angle_in_degrees )
323 {
324 const int N = 64;
325
326 static const double sin_table[] =
327 {
328 0.00000000000000000000, 0.09801714032956060400,
329 0.19509032201612825000, 0.29028467725446233000,
330 0.38268343236508978000, 0.47139673682599764000,
331 0.55557023301960218000, 0.63439328416364549000,
332 0.70710678118654746000, 0.77301045336273699000,
333 0.83146961230254524000, 0.88192126434835494000,
334 0.92387953251128674000, 0.95694033573220894000,
335 0.98078528040323043000, 0.99518472667219682000,
336 1.00000000000000000000, 0.99518472667219693000,
337 0.98078528040323043000, 0.95694033573220894000,
338 0.92387953251128674000, 0.88192126434835505000,
339 0.83146961230254546000, 0.77301045336273710000,
340 0.70710678118654757000, 0.63439328416364549000,
341 0.55557023301960218000, 0.47139673682599786000,
342 0.38268343236508989000, 0.29028467725446239000,
343 0.19509032201612861000, 0.09801714032956082600,
344 0.00000000000000012246, -0.09801714032956059000,
345 -0.19509032201612836000, -0.29028467725446211000,
346 -0.38268343236508967000, -0.47139673682599764000,
347 -0.55557023301960196000, -0.63439328416364527000,
348 -0.70710678118654746000, -0.77301045336273666000,
349 -0.83146961230254524000, -0.88192126434835494000,
350 -0.92387953251128652000, -0.95694033573220882000,
351 -0.98078528040323032000, -0.99518472667219693000,
352 -1.00000000000000000000, -0.99518472667219693000,
353 -0.98078528040323043000, -0.95694033573220894000,
354 -0.92387953251128663000, -0.88192126434835505000,
355 -0.83146961230254546000, -0.77301045336273688000,
356 -0.70710678118654768000, -0.63439328416364593000,
357 -0.55557023301960218000, -0.47139673682599792000,
358 -0.38268343236509039000, -0.29028467725446250000,
359 -0.19509032201612872000, -0.09801714032956050600,
360 };
361
362 static const double k2 = (2*CV_PI)/N;
363
364 static const double sin_a0 = -0.166630293345647*k2*k2*k2;
365 static const double sin_a2 = k2;
366
367 static const double cos_a0 = -0.499818138450326*k2*k2;
368 /*static const double cos_a2 = 1;*/
369
370 double k1;
371 int i = 0;
372
373 if( !angle_in_degrees )
374 k1 = N/(2*CV_PI);
375 else
376 k1 = N/360.;
377
378 #if CV_AVX2
379 if (USE_AVX2)
380 {
381 __m128d v_k1 = _mm_set1_pd(k1);
382 __m128d v_1 = _mm_set1_pd(1);
383 __m128i v_N1 = _mm_set1_epi32(N - 1);
384 __m128i v_N4 = _mm_set1_epi32(N >> 2);
385 __m128d v_sin_a0 = _mm_set1_pd(sin_a0);
386 __m128d v_sin_a2 = _mm_set1_pd(sin_a2);
387 __m128d v_cos_a0 = _mm_set1_pd(cos_a0);
388
389 for ( ; i <= len - 4; i += 4)
390 {
391 __m128 v_angle = _mm_loadu_ps(angle + i);
392
393 // 0-1
394 __m128d v_t = _mm_mul_pd(_mm_cvtps_pd(v_angle), v_k1);
395 __m128i v_it = _mm_cvtpd_epi32(v_t);
396 v_t = _mm_sub_pd(v_t, _mm_cvtepi32_pd(v_it));
397
398 __m128i v_sin_idx = _mm_and_si128(v_it, v_N1);
399 __m128i v_cos_idx = _mm_and_si128(_mm_sub_epi32(v_N4, v_sin_idx), v_N1);
400
401 __m128d v_t2 = _mm_mul_pd(v_t, v_t);
402 __m128d v_sin_b = _mm_mul_pd(_mm_add_pd(_mm_mul_pd(v_sin_a0, v_t2), v_sin_a2), v_t);
403 __m128d v_cos_b = _mm_add_pd(_mm_mul_pd(v_cos_a0, v_t2), v_1);
404
405 __m128d v_sin_a = _mm_i32gather_pd(sin_table, v_sin_idx, 8);
406 __m128d v_cos_a = _mm_i32gather_pd(sin_table, v_cos_idx, 8);
407
408 __m128d v_sin_val_0 = _mm_add_pd(_mm_mul_pd(v_sin_a, v_cos_b),
409 _mm_mul_pd(v_cos_a, v_sin_b));
410 __m128d v_cos_val_0 = _mm_sub_pd(_mm_mul_pd(v_cos_a, v_cos_b),
411 _mm_mul_pd(v_sin_a, v_sin_b));
412
413 // 2-3
414 v_t = _mm_mul_pd(_mm_cvtps_pd(_mm_castsi128_ps(_mm_srli_si128(_mm_castps_si128(v_angle), 8))), v_k1);
415 v_it = _mm_cvtpd_epi32(v_t);
416 v_t = _mm_sub_pd(v_t, _mm_cvtepi32_pd(v_it));
417
418 v_sin_idx = _mm_and_si128(v_it, v_N1);
419 v_cos_idx = _mm_and_si128(_mm_sub_epi32(v_N4, v_sin_idx), v_N1);
420
421 v_t2 = _mm_mul_pd(v_t, v_t);
422 v_sin_b = _mm_mul_pd(_mm_add_pd(_mm_mul_pd(v_sin_a0, v_t2), v_sin_a2), v_t);
423 v_cos_b = _mm_add_pd(_mm_mul_pd(v_cos_a0, v_t2), v_1);
424
425 v_sin_a = _mm_i32gather_pd(sin_table, v_sin_idx, 8);
426 v_cos_a = _mm_i32gather_pd(sin_table, v_cos_idx, 8);
427
428 __m128d v_sin_val_1 = _mm_add_pd(_mm_mul_pd(v_sin_a, v_cos_b),
429 _mm_mul_pd(v_cos_a, v_sin_b));
430 __m128d v_cos_val_1 = _mm_sub_pd(_mm_mul_pd(v_cos_a, v_cos_b),
431 _mm_mul_pd(v_sin_a, v_sin_b));
432
433 _mm_storeu_ps(sinval + i, _mm_movelh_ps(_mm_cvtpd_ps(v_sin_val_0),
434 _mm_cvtpd_ps(v_sin_val_1)));
435 _mm_storeu_ps(cosval + i, _mm_movelh_ps(_mm_cvtpd_ps(v_cos_val_0),
436 _mm_cvtpd_ps(v_cos_val_1)));
437 }
438 }
439 #endif
440
441 for( ; i < len; i++ )
442 {
443 double t = angle[i]*k1;
444 int it = cvRound(t);
445 t -= it;
446 int sin_idx = it & (N - 1);
447 int cos_idx = (N/4 - sin_idx) & (N - 1);
448
449 double sin_b = (sin_a0*t*t + sin_a2)*t;
450 double cos_b = cos_a0*t*t + 1;
451
452 double sin_a = sin_table[sin_idx];
453 double cos_a = sin_table[cos_idx];
454
455 double sin_val = sin_a*cos_b + cos_a*sin_b;
456 double cos_val = cos_a*cos_b - sin_a*sin_b;
457
458 sinval[i] = (float)sin_val;
459 cosval[i] = (float)cos_val;
460 }
461 }
462
463
464 #ifdef HAVE_OPENCL
465
ocl_polarToCart(InputArray _mag,InputArray _angle,OutputArray _dst1,OutputArray _dst2,bool angleInDegrees)466 static bool ocl_polarToCart( InputArray _mag, InputArray _angle,
467 OutputArray _dst1, OutputArray _dst2, bool angleInDegrees )
468 {
469 const ocl::Device & d = ocl::Device::getDefault();
470 int type = _angle.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type),
471 rowsPerWI = d.isIntel() ? 4 : 1;
472 bool doubleSupport = d.doubleFPConfig() > 0;
473
474 if ( !doubleSupport && depth == CV_64F )
475 return false;
476
477 ocl::Kernel k("KF", ocl::core::arithm_oclsrc,
478 format("-D dstT=%s -D DEPTH_dst=%d -D rowsPerWI=%d -D BINARY_OP -D OP_PTC_%s%s",
479 ocl::typeToStr(CV_MAKE_TYPE(depth, 1)), depth,
480 rowsPerWI,
481 angleInDegrees ? "AD" : "AR",
482 doubleSupport ? " -D DOUBLE_SUPPORT" : ""));
483 if (k.empty())
484 return false;
485
486 UMat mag = _mag.getUMat(), angle = _angle.getUMat();
487 Size size = angle.size();
488 CV_Assert(mag.size() == size);
489
490 _dst1.create(size, type);
491 _dst2.create(size, type);
492 UMat dst1 = _dst1.getUMat(), dst2 = _dst2.getUMat();
493
494 k.args(ocl::KernelArg::ReadOnlyNoSize(mag), ocl::KernelArg::ReadOnlyNoSize(angle),
495 ocl::KernelArg::WriteOnly(dst1, cn), ocl::KernelArg::WriteOnlyNoSize(dst2));
496
497 size_t globalsize[2] = { (size_t)dst1.cols * cn, ((size_t)dst1.rows + rowsPerWI - 1) / rowsPerWI };
498 return k.run(2, globalsize, NULL, false);
499 }
500
501 #endif
502
503 #ifdef HAVE_IPP
ipp_polarToCart(Mat & mag,Mat & angle,Mat & x,Mat & y)504 static bool ipp_polarToCart(Mat &mag, Mat &angle, Mat &x, Mat &y)
505 {
506 CV_INSTRUMENT_REGION_IPP();
507
508 int depth = angle.depth();
509 if(depth != CV_32F && depth != CV_64F)
510 return false;
511
512 if(angle.dims <= 2)
513 {
514 int len = (int)(angle.cols*angle.channels());
515
516 if(depth == CV_32F)
517 {
518 for (int h = 0; h < angle.rows; h++)
519 {
520 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)
521 return false;
522 }
523 }
524 else
525 {
526 for (int h = 0; h < angle.rows; h++)
527 {
528 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)
529 return false;
530 }
531 }
532 return true;
533 }
534 else
535 {
536 const Mat *arrays[] = {&mag, &angle, &x, &y, NULL};
537 uchar *ptrs[4] = {NULL};
538 NAryMatIterator it(arrays, ptrs);
539 int len = (int)(it.size*angle.channels());
540
541 if(depth == CV_32F)
542 {
543 for (size_t i = 0; i < it.nplanes; i++, ++it)
544 {
545 if(CV_INSTRUMENT_FUN_IPP(ippsPolarToCart_32f, (const float*)ptrs[0], (const float*)ptrs[1], (float*)ptrs[2], (float*)ptrs[3], len) < 0)
546 return false;
547 }
548 }
549 else
550 {
551 for (size_t i = 0; i < it.nplanes; i++, ++it)
552 {
553 if(CV_INSTRUMENT_FUN_IPP(ippsPolarToCart_64f, (const double*)ptrs[0], (const double*)ptrs[1], (double*)ptrs[2], (double*)ptrs[3], len) < 0)
554 return false;
555 }
556 }
557 return true;
558 }
559 }
560 #endif
561
polarToCart(InputArray src1,InputArray src2,OutputArray dst1,OutputArray dst2,bool angleInDegrees)562 void polarToCart( InputArray src1, InputArray src2,
563 OutputArray dst1, OutputArray dst2, bool angleInDegrees )
564 {
565 CV_INSTRUMENT_REGION();
566
567 int type = src2.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
568 CV_Assert((depth == CV_32F || depth == CV_64F) && (src1.empty() || src1.type() == type));
569
570 CV_OCL_RUN(!src1.empty() && src2.dims() <= 2 && dst1.isUMat() && dst2.isUMat(),
571 ocl_polarToCart(src1, src2, dst1, dst2, angleInDegrees))
572
573 Mat Mag = src1.getMat(), Angle = src2.getMat();
574 CV_Assert( Mag.empty() || Angle.size == Mag.size);
575 dst1.create( Angle.dims, Angle.size, type );
576 dst2.create( Angle.dims, Angle.size, type );
577 Mat X = dst1.getMat(), Y = dst2.getMat();
578
579 CV_IPP_RUN(!angleInDegrees, ipp_polarToCart(Mag, Angle, X, Y));
580
581 const Mat* arrays[] = {&Mag, &Angle, &X, &Y, 0};
582 uchar* ptrs[4] = {};
583 NAryMatIterator it(arrays, ptrs);
584 cv::AutoBuffer<float> _buf;
585 float* buf[2] = {0, 0};
586 int j, k, total = (int)(it.size*cn), blockSize = std::min(total, ((BLOCK_SIZE+cn-1)/cn)*cn);
587 size_t esz1 = Angle.elemSize1();
588
589 if( depth == CV_64F )
590 {
591 _buf.allocate(blockSize*2);
592 buf[0] = _buf.data();
593 buf[1] = buf[0] + blockSize;
594 }
595
596 for( size_t i = 0; i < it.nplanes; i++, ++it )
597 {
598 for( j = 0; j < total; j += blockSize )
599 {
600 int len = std::min(total - j, blockSize);
601 if( depth == CV_32F )
602 {
603 const float *mag = (const float*)ptrs[0], *angle = (const float*)ptrs[1];
604 float *x = (float*)ptrs[2], *y = (float*)ptrs[3];
605
606 SinCos_32f( angle, y, x, len, angleInDegrees );
607 if( mag )
608 {
609 k = 0;
610
611 #if CV_SIMD
612 int cWidth = v_float32::nlanes;
613 for( ; k <= len - cWidth; k += cWidth )
614 {
615 v_float32 v_m = vx_load(mag + k);
616 v_store(x + k, vx_load(x + k) * v_m);
617 v_store(y + k, vx_load(y + k) * v_m);
618 }
619 vx_cleanup();
620 #endif
621
622 for( ; k < len; k++ )
623 {
624 float m = mag[k];
625 x[k] *= m; y[k] *= m;
626 }
627 }
628 }
629 else
630 {
631 const double *mag = (const double*)ptrs[0], *angle = (const double*)ptrs[1];
632 double *x = (double*)ptrs[2], *y = (double*)ptrs[3];
633
634 for( k = 0; k < len; k++ )
635 buf[0][k] = (float)angle[k];
636
637 SinCos_32f( buf[0], buf[1], buf[0], len, angleInDegrees );
638 if( mag )
639 for( k = 0; k < len; k++ )
640 {
641 double m = mag[k];
642 x[k] = buf[0][k]*m; y[k] = buf[1][k]*m;
643 }
644 else
645 {
646 std::memcpy(x, buf[0], sizeof(float) * len);
647 std::memcpy(y, buf[1], sizeof(float) * len);
648 }
649 }
650
651 if( ptrs[0] )
652 ptrs[0] += len*esz1;
653 ptrs[1] += len*esz1;
654 ptrs[2] += len*esz1;
655 ptrs[3] += len*esz1;
656 }
657 }
658 }
659
660 /****************************************************************************************\
661 * E X P *
662 \****************************************************************************************/
663
exp(InputArray _src,OutputArray _dst)664 void exp( InputArray _src, OutputArray _dst )
665 {
666 CV_INSTRUMENT_REGION();
667
668 int type = _src.type(), depth = _src.depth(), cn = _src.channels();
669 CV_Assert( depth == CV_32F || depth == CV_64F );
670
671 CV_OCL_RUN(_dst.isUMat() && _src.dims() <= 2,
672 ocl_math_op(_src, noArray(), _dst, OCL_OP_EXP))
673
674 Mat src = _src.getMat();
675 _dst.create( src.dims, src.size, type );
676 Mat dst = _dst.getMat();
677
678 const Mat* arrays[] = {&src, &dst, 0};
679 uchar* ptrs[2] = {};
680 NAryMatIterator it(arrays, ptrs);
681 int len = (int)(it.size*cn);
682
683 for( size_t i = 0; i < it.nplanes; i++, ++it )
684 {
685 if( depth == CV_32F )
686 hal::exp32f((const float*)ptrs[0], (float*)ptrs[1], len);
687 else
688 hal::exp64f((const double*)ptrs[0], (double*)ptrs[1], len);
689 }
690 }
691
692
693 /****************************************************************************************\
694 * L O G *
695 \****************************************************************************************/
696
log(InputArray _src,OutputArray _dst)697 void log( InputArray _src, OutputArray _dst )
698 {
699 CV_INSTRUMENT_REGION();
700
701 int type = _src.type(), depth = _src.depth(), cn = _src.channels();
702 CV_Assert( depth == CV_32F || depth == CV_64F );
703
704 CV_OCL_RUN( _dst.isUMat() && _src.dims() <= 2,
705 ocl_math_op(_src, noArray(), _dst, OCL_OP_LOG))
706
707 Mat src = _src.getMat();
708 _dst.create( src.dims, src.size, type );
709 Mat dst = _dst.getMat();
710
711 const Mat* arrays[] = {&src, &dst, 0};
712 uchar* ptrs[2] = {};
713 NAryMatIterator it(arrays, ptrs);
714 int len = (int)(it.size*cn);
715
716 for( size_t i = 0; i < it.nplanes; i++, ++it )
717 {
718 if( depth == CV_32F )
719 hal::log32f( (const float*)ptrs[0], (float*)ptrs[1], len );
720 else
721 hal::log64f( (const double*)ptrs[0], (double*)ptrs[1], len );
722 }
723 }
724
725 /****************************************************************************************\
726 * P O W E R *
727 \****************************************************************************************/
728
729 template <typename T, typename WT>
730 struct iPow_SIMD
731 {
operator ()cv::iPow_SIMD732 int operator() ( const T *, T *, int, int)
733 {
734 return 0;
735 }
736 };
737
738 #if CV_SIMD
739
740 template <>
741 struct iPow_SIMD<uchar, int>
742 {
operator ()cv::iPow_SIMD743 int operator() ( const uchar * src, uchar * dst, int len, int power )
744 {
745 int i = 0;
746 v_uint32 v_1 = vx_setall_u32(1u);
747
748 for ( ; i <= len - v_uint16::nlanes; i += v_uint16::nlanes)
749 {
750 v_uint32 v_a1 = v_1, v_a2 = v_1;
751 v_uint16 v = vx_load_expand(src + i);
752 v_uint32 v_b1, v_b2;
753 v_expand(v, v_b1, v_b2);
754 int p = power;
755
756 while( p > 1 )
757 {
758 if (p & 1)
759 {
760 v_a1 *= v_b1;
761 v_a2 *= v_b2;
762 }
763 v_b1 *= v_b1;
764 v_b2 *= v_b2;
765 p >>= 1;
766 }
767
768 v_a1 *= v_b1;
769 v_a2 *= v_b2;
770
771 v = v_pack(v_a1, v_a2);
772 v_pack_store(dst + i, v);
773 }
774 vx_cleanup();
775
776 return i;
777 }
778 };
779
780 template <>
781 struct iPow_SIMD<schar, int>
782 {
operator ()cv::iPow_SIMD783 int operator() ( const schar * src, schar * dst, int len, int power)
784 {
785 int i = 0;
786 v_int32 v_1 = vx_setall_s32(1);
787
788 for ( ; i <= len - v_int16::nlanes; i += v_int16::nlanes)
789 {
790 v_int32 v_a1 = v_1, v_a2 = v_1;
791 v_int16 v = vx_load_expand(src + i);
792 v_int32 v_b1, v_b2;
793 v_expand(v, v_b1, v_b2);
794 int p = power;
795
796 while( p > 1 )
797 {
798 if (p & 1)
799 {
800 v_a1 *= v_b1;
801 v_a2 *= v_b2;
802 }
803 v_b1 *= v_b1;
804 v_b2 *= v_b2;
805 p >>= 1;
806 }
807
808 v_a1 *= v_b1;
809 v_a2 *= v_b2;
810
811 v = v_pack(v_a1, v_a2);
812 v_pack_store(dst + i, v);
813 }
814 vx_cleanup();
815
816 return i;
817 }
818 };
819
820 template <>
821 struct iPow_SIMD<ushort, int>
822 {
operator ()cv::iPow_SIMD823 int operator() ( const ushort * src, ushort * dst, int len, int power)
824 {
825 int i = 0;
826 v_uint32 v_1 = vx_setall_u32(1u);
827
828 for ( ; i <= len - v_uint16::nlanes; i += v_uint16::nlanes)
829 {
830 v_uint32 v_a1 = v_1, v_a2 = v_1;
831 v_uint16 v = vx_load(src + i);
832 v_uint32 v_b1, v_b2;
833 v_expand(v, v_b1, v_b2);
834 int p = power;
835
836 while( p > 1 )
837 {
838 if (p & 1)
839 {
840 v_a1 *= v_b1;
841 v_a2 *= v_b2;
842 }
843 v_b1 *= v_b1;
844 v_b2 *= v_b2;
845 p >>= 1;
846 }
847
848 v_a1 *= v_b1;
849 v_a2 *= v_b2;
850
851 v = v_pack(v_a1, v_a2);
852 v_store(dst + i, v);
853 }
854 vx_cleanup();
855
856 return i;
857 }
858 };
859
860 template <>
861 struct iPow_SIMD<short, int>
862 {
operator ()cv::iPow_SIMD863 int operator() ( const short * src, short * dst, int len, int power)
864 {
865 int i = 0;
866 v_int32 v_1 = vx_setall_s32(1);
867
868 for ( ; i <= len - v_int16::nlanes; i += v_int16::nlanes)
869 {
870 v_int32 v_a1 = v_1, v_a2 = v_1;
871 v_int16 v = vx_load(src + i);
872 v_int32 v_b1, v_b2;
873 v_expand(v, v_b1, v_b2);
874 int p = power;
875
876 while( p > 1 )
877 {
878 if (p & 1)
879 {
880 v_a1 *= v_b1;
881 v_a2 *= v_b2;
882 }
883 v_b1 *= v_b1;
884 v_b2 *= v_b2;
885 p >>= 1;
886 }
887
888 v_a1 *= v_b1;
889 v_a2 *= v_b2;
890
891 v = v_pack(v_a1, v_a2);
892 v_store(dst + i, v);
893 }
894 vx_cleanup();
895
896 return i;
897 }
898 };
899
900 template <>
901 struct iPow_SIMD<int, int>
902 {
operator ()cv::iPow_SIMD903 int operator() ( const int * src, int * dst, int len, int power)
904 {
905 int i = 0;
906 v_int32 v_1 = vx_setall_s32(1);
907
908 for ( ; i <= len - v_int32::nlanes*2; i += v_int32::nlanes*2)
909 {
910 v_int32 v_a1 = v_1, v_a2 = v_1;
911 v_int32 v_b1 = vx_load(src + i), v_b2 = vx_load(src + i + v_int32::nlanes);
912 int p = power;
913
914 while( p > 1 )
915 {
916 if (p & 1)
917 {
918 v_a1 *= v_b1;
919 v_a2 *= v_b2;
920 }
921 v_b1 *= v_b1;
922 v_b2 *= v_b2;
923 p >>= 1;
924 }
925
926 v_a1 *= v_b1;
927 v_a2 *= v_b2;
928
929 v_store(dst + i, v_a1);
930 v_store(dst + i + v_int32::nlanes, v_a2);
931 }
932 vx_cleanup();
933
934 return i;
935 }
936 };
937
938 template <>
939 struct iPow_SIMD<float, float>
940 {
operator ()cv::iPow_SIMD941 int operator() ( const float * src, float * dst, int len, int power)
942 {
943 int i = 0;
944 v_float32 v_1 = vx_setall_f32(1.f);
945
946 for ( ; i <= len - v_float32::nlanes*2; i += v_float32::nlanes*2)
947 {
948 v_float32 v_a1 = v_1, v_a2 = v_1;
949 v_float32 v_b1 = vx_load(src + i), v_b2 = vx_load(src + i + v_float32::nlanes);
950 int p = std::abs(power);
951 if( power < 0 )
952 {
953 v_b1 = v_1 / v_b1;
954 v_b2 = v_1 / v_b2;
955 }
956
957 while( p > 1 )
958 {
959 if (p & 1)
960 {
961 v_a1 *= v_b1;
962 v_a2 *= v_b2;
963 }
964 v_b1 *= v_b1;
965 v_b2 *= v_b2;
966 p >>= 1;
967 }
968
969 v_a1 *= v_b1;
970 v_a2 *= v_b2;
971
972 v_store(dst + i, v_a1);
973 v_store(dst + i + v_float32::nlanes, v_a2);
974 }
975 vx_cleanup();
976
977 return i;
978 }
979 };
980
981 #if CV_SIMD_64F
982 template <>
983 struct iPow_SIMD<double, double>
984 {
operator ()cv::iPow_SIMD985 int operator() ( const double * src, double * dst, int len, int power)
986 {
987 int i = 0;
988 v_float64 v_1 = vx_setall_f64(1.);
989
990 for ( ; i <= len - v_float64::nlanes*2; i += v_float64::nlanes*2)
991 {
992 v_float64 v_a1 = v_1, v_a2 = v_1;
993 v_float64 v_b1 = vx_load(src + i), v_b2 = vx_load(src + i + v_float64::nlanes);
994 int p = std::abs(power);
995 if( power < 0 )
996 {
997 v_b1 = v_1 / v_b1;
998 v_b2 = v_1 / v_b2;
999 }
1000
1001 while( p > 1 )
1002 {
1003 if (p & 1)
1004 {
1005 v_a1 *= v_b1;
1006 v_a2 *= v_b2;
1007 }
1008 v_b1 *= v_b1;
1009 v_b2 *= v_b2;
1010 p >>= 1;
1011 }
1012
1013 v_a1 *= v_b1;
1014 v_a2 *= v_b2;
1015
1016 v_store(dst + i, v_a1);
1017 v_store(dst + i + v_float64::nlanes, v_a2);
1018 }
1019 vx_cleanup();
1020
1021 return i;
1022 }
1023 };
1024 #endif
1025
1026 #endif
1027
1028 template<typename T, typename WT>
1029 static void
iPow_i(const T * src,T * dst,int len,int power)1030 iPow_i( const T* src, T* dst, int len, int power )
1031 {
1032 if( power < 0 )
1033 {
1034 T tab[5] =
1035 {
1036 saturate_cast<T>(power == -1 ? -1 : 0), saturate_cast<T>((power & 1) ? -1 : 1),
1037 std::numeric_limits<T>::max(), 1, saturate_cast<T>(power == -1 ? 1 : 0)
1038 };
1039 for( int i = 0; i < len; i++ )
1040 {
1041 T val = src[i];
1042 dst[i] = cv_abs(val) <= 2 ? tab[val + 2] : (T)0;
1043 }
1044 }
1045 else
1046 {
1047 iPow_SIMD<T, WT> vop;
1048 int i = vop(src, dst, len, power);
1049
1050 for( ; i < len; i++ )
1051 {
1052 WT a = 1, b = src[i];
1053 int p = power;
1054 while( p > 1 )
1055 {
1056 if( p & 1 )
1057 a *= b;
1058 b *= b;
1059 p >>= 1;
1060 }
1061
1062 a *= b;
1063 dst[i] = saturate_cast<T>(a);
1064 }
1065 }
1066 }
1067
1068 template<typename T>
1069 static void
iPow_f(const T * src,T * dst,int len,int power0)1070 iPow_f( const T* src, T* dst, int len, int power0 )
1071 {
1072 iPow_SIMD<T, T> vop;
1073 int i = vop(src, dst, len, power0);
1074 int power = std::abs(power0);
1075
1076 for( ; i < len; i++ )
1077 {
1078 T a = 1, b = src[i];
1079 int p = power;
1080 if( power0 < 0 )
1081 b = 1/b;
1082
1083 while( p > 1 )
1084 {
1085 if( p & 1 )
1086 a *= b;
1087 b *= b;
1088 p >>= 1;
1089 }
1090
1091 a *= b;
1092 dst[i] = a;
1093 }
1094 }
1095
iPow8u(const uchar * src,uchar * dst,int len,int power)1096 static void iPow8u(const uchar* src, uchar* dst, int len, int power)
1097 {
1098 iPow_i<uchar, unsigned>(src, dst, len, power);
1099 }
1100
iPow8s(const schar * src,schar * dst,int len,int power)1101 static void iPow8s(const schar* src, schar* dst, int len, int power)
1102 {
1103 iPow_i<schar, int>(src, dst, len, power);
1104 }
1105
iPow16u(const ushort * src,ushort * dst,int len,int power)1106 static void iPow16u(const ushort* src, ushort* dst, int len, int power)
1107 {
1108 iPow_i<ushort, unsigned>(src, dst, len, power);
1109 }
1110
iPow16s(const short * src,short * dst,int len,int power)1111 static void iPow16s(const short* src, short* dst, int len, int power)
1112 {
1113 iPow_i<short, int>(src, dst, len, power);
1114 }
1115
iPow32s(const int * src,int * dst,int len,int power)1116 static void iPow32s(const int* src, int* dst, int len, int power)
1117 {
1118 iPow_i<int, int>(src, dst, len, power);
1119 }
1120
iPow32f(const float * src,float * dst,int len,int power)1121 static void iPow32f(const float* src, float* dst, int len, int power)
1122 {
1123 iPow_f<float>(src, dst, len, power);
1124 }
1125
iPow64f(const double * src,double * dst,int len,int power)1126 static void iPow64f(const double* src, double* dst, int len, int power)
1127 {
1128 iPow_f<double>(src, dst, len, power);
1129 }
1130
1131
1132 typedef void (*IPowFunc)( const uchar* src, uchar* dst, int len, int power );
1133
1134 static IPowFunc ipowTab[] =
1135 {
1136 (IPowFunc)iPow8u, (IPowFunc)iPow8s, (IPowFunc)iPow16u, (IPowFunc)iPow16s,
1137 (IPowFunc)iPow32s, (IPowFunc)iPow32f, (IPowFunc)iPow64f, 0
1138 };
1139
1140 #ifdef HAVE_OPENCL
1141
ocl_pow(InputArray _src,double power,OutputArray _dst,bool is_ipower,int ipower)1142 static bool ocl_pow(InputArray _src, double power, OutputArray _dst,
1143 bool is_ipower, int ipower)
1144 {
1145 const ocl::Device & d = ocl::Device::getDefault();
1146 int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type),
1147 rowsPerWI = d.isIntel() ? 4 : 1;
1148 bool doubleSupport = d.doubleFPConfig() > 0;
1149
1150 _dst.createSameSize(_src, type);
1151 if (is_ipower)
1152 {
1153 if (ipower == 0)
1154 {
1155 _dst.setTo(Scalar::all(1));
1156 return true;
1157 }
1158 if (ipower == 1)
1159 {
1160 _src.copyTo(_dst);
1161 return true;
1162 }
1163 if( ipower < 0 )
1164 {
1165 if( depth == CV_32F || depth == CV_64F )
1166 is_ipower = false;
1167 else
1168 return false;
1169 }
1170 }
1171
1172 if (depth == CV_64F && !doubleSupport)
1173 return false;
1174
1175 bool issqrt = std::abs(power - 0.5) < DBL_EPSILON;
1176 const char * const op = issqrt ? "OP_SQRT" : is_ipower ? "OP_POWN" : "OP_POW";
1177
1178 ocl::Kernel k("KF", ocl::core::arithm_oclsrc,
1179 format("-D dstT=%s -D DEPTH_dst=%d -D rowsPerWI=%d -D %s -D UNARY_OP%s",
1180 ocl::typeToStr(depth), depth, rowsPerWI, op,
1181 doubleSupport ? " -D DOUBLE_SUPPORT" : ""));
1182 if (k.empty())
1183 return false;
1184
1185 UMat src = _src.getUMat();
1186 _dst.create(src.size(), type);
1187 UMat dst = _dst.getUMat();
1188
1189 ocl::KernelArg srcarg = ocl::KernelArg::ReadOnlyNoSize(src),
1190 dstarg = ocl::KernelArg::WriteOnly(dst, cn);
1191
1192 if (issqrt)
1193 k.args(srcarg, dstarg);
1194 else if (is_ipower)
1195 k.args(srcarg, dstarg, ipower);
1196 else
1197 {
1198 if (depth == CV_32F)
1199 k.args(srcarg, dstarg, (float)power);
1200 else
1201 k.args(srcarg, dstarg, power);
1202 }
1203
1204 size_t globalsize[2] = { (size_t)dst.cols * cn, ((size_t)dst.rows + rowsPerWI - 1) / rowsPerWI };
1205 return k.run(2, globalsize, NULL, false);
1206 }
1207
1208 #endif
1209
pow(InputArray _src,double power,OutputArray _dst)1210 void pow( InputArray _src, double power, OutputArray _dst )
1211 {
1212 CV_INSTRUMENT_REGION();
1213
1214 int type = _src.type(), depth = CV_MAT_DEPTH(type),
1215 cn = CV_MAT_CN(type), ipower = cvRound(power);
1216 bool is_ipower = fabs(ipower - power) < DBL_EPSILON;
1217 #ifdef HAVE_OPENCL
1218 bool useOpenCL = _dst.isUMat() && _src.dims() <= 2;
1219 #endif
1220
1221 if( is_ipower
1222 #ifdef HAVE_OPENCL
1223 && !(useOpenCL && ocl::Device::getDefault().isIntel() && depth != CV_64F)
1224 #endif
1225 )
1226 {
1227 switch( ipower )
1228 {
1229 case 0:
1230 _dst.createSameSize(_src, type);
1231 _dst.setTo(Scalar::all(1));
1232 return;
1233 case 1:
1234 _src.copyTo(_dst);
1235 return;
1236 case 2:
1237 multiply(_src, _src, _dst);
1238 return;
1239 }
1240 }
1241 else
1242 CV_Assert( depth == CV_32F || depth == CV_64F );
1243
1244 CV_OCL_RUN(useOpenCL, ocl_pow(_src, power, _dst, is_ipower, ipower))
1245
1246 Mat src = _src.getMat();
1247 _dst.create( src.dims, src.size, type );
1248 Mat dst = _dst.getMat();
1249
1250 const Mat* arrays[] = {&src, &dst, 0};
1251 uchar* ptrs[2] = {};
1252 NAryMatIterator it(arrays, ptrs);
1253 int len = (int)(it.size*cn);
1254
1255 if( is_ipower )
1256 {
1257 IPowFunc func = ipowTab[depth];
1258 CV_Assert( func != 0 );
1259
1260 for( size_t i = 0; i < it.nplanes; i++, ++it )
1261 func( ptrs[0], ptrs[1], len, ipower );
1262 }
1263 else if( fabs(fabs(power) - 0.5) < DBL_EPSILON )
1264 {
1265 MathFunc func = power < 0 ?
1266 (depth == CV_32F ? (MathFunc)hal::invSqrt32f : (MathFunc)hal::invSqrt64f) :
1267 (depth == CV_32F ? (MathFunc)hal::sqrt32f : (MathFunc)hal::sqrt64f);
1268
1269 for( size_t i = 0; i < it.nplanes; i++, ++it )
1270 func( ptrs[0], ptrs[1], len );
1271 }
1272 else
1273 {
1274 int j, k, blockSize = std::min(len, ((BLOCK_SIZE + cn-1)/cn)*cn);
1275 size_t esz1 = src.elemSize1();
1276 AutoBuffer<uchar> buf;
1277 Cv32suf inf32, nan32;
1278 Cv64suf inf64, nan64;
1279 float* fbuf = 0;
1280 double* dbuf = 0;
1281 #ifndef __EMSCRIPTEN__
1282 inf32.i = 0x7f800000;
1283 nan32.i = 0x7fffffff;
1284 inf64.i = CV_BIG_INT(0x7FF0000000000000);
1285 nan64.i = CV_BIG_INT(0x7FFFFFFFFFFFFFFF);
1286 #else
1287 inf32.f = std::numeric_limits<float>::infinity();
1288 nan32.f = std::numeric_limits<float>::quiet_NaN();
1289 inf64.f = std::numeric_limits<double>::infinity();
1290 nan64.f = std::numeric_limits<double>::quiet_NaN();
1291 #endif
1292
1293 if( src.ptr() == dst.ptr() )
1294 {
1295 buf.allocate(blockSize*esz1);
1296 fbuf = (float*)buf.data();
1297 dbuf = (double*)buf.data();
1298 }
1299
1300 for( size_t i = 0; i < it.nplanes; i++, ++it )
1301 {
1302 for( j = 0; j < len; j += blockSize )
1303 {
1304 int bsz = std::min(len - j, blockSize);
1305
1306 if( depth == CV_32F )
1307 {
1308 float* x0 = (float*)ptrs[0];
1309 float* x = fbuf ? fbuf : x0;
1310 float* y = (float*)ptrs[1];
1311
1312 if( x != x0 )
1313 memcpy(x, x0, bsz*esz1);
1314
1315 hal::log32f(x, y, bsz);
1316 for( k = 0; k < bsz; k++ )
1317 y[k] = (float)(y[k]*power);
1318 hal::exp32f(y, y, bsz);
1319 for( k = 0; k < bsz; k++ )
1320 {
1321 if( x0[k] <= 0 )
1322 {
1323 if( x0[k] == 0.f )
1324 {
1325 if( power < 0 )
1326 y[k] = inf32.f;
1327 }
1328 else
1329 y[k] = nan32.f;
1330 }
1331 }
1332 }
1333 else
1334 {
1335 double* x0 = (double*)ptrs[0];
1336 double* x = dbuf ? dbuf : x0;
1337 double* y = (double*)ptrs[1];
1338
1339 if( x != x0 )
1340 memcpy(x, x0, bsz*esz1);
1341
1342 hal::log64f(x, y, bsz);
1343 for( k = 0; k < bsz; k++ )
1344 y[k] *= power;
1345 hal::exp64f(y, y, bsz);
1346
1347 for( k = 0; k < bsz; k++ )
1348 {
1349 if( x0[k] <= 0 )
1350 {
1351 if( x0[k] == 0. )
1352 {
1353 if( power < 0 )
1354 y[k] = inf64.f;
1355 }
1356 else
1357 y[k] = nan64.f;
1358 }
1359 }
1360 }
1361 ptrs[0] += bsz*esz1;
1362 ptrs[1] += bsz*esz1;
1363 }
1364 }
1365 }
1366 }
1367
sqrt(InputArray a,OutputArray b)1368 void sqrt(InputArray a, OutputArray b)
1369 {
1370 CV_INSTRUMENT_REGION();
1371
1372 cv::pow(a, 0.5, b);
1373 }
1374
1375 /************************** CheckArray for NaN's, Inf's *********************************/
1376
1377 template<int cv_mat_type> struct mat_type_assotiations{};
1378
1379 template<> struct mat_type_assotiations<CV_8U>
1380 {
1381 typedef unsigned char type;
1382 static const type min_allowable = 0x0;
1383 static const type max_allowable = 0xFF;
1384 };
1385
1386 template<> struct mat_type_assotiations<CV_8S>
1387 {
1388 typedef signed char type;
1389 static const type min_allowable = SCHAR_MIN;
1390 static const type max_allowable = SCHAR_MAX;
1391 };
1392
1393 template<> struct mat_type_assotiations<CV_16U>
1394 {
1395 typedef unsigned short type;
1396 static const type min_allowable = 0x0;
1397 static const type max_allowable = USHRT_MAX;
1398 };
1399 template<> struct mat_type_assotiations<CV_16S>
1400 {
1401 typedef signed short type;
1402 static const type min_allowable = SHRT_MIN;
1403 static const type max_allowable = SHRT_MAX;
1404 };
1405
1406 template<> struct mat_type_assotiations<CV_32S>
1407 {
1408 typedef int type;
1409 static const type min_allowable = (-INT_MAX - 1);
1410 static const type max_allowable = INT_MAX;
1411 };
1412
1413 template<int depth>
checkIntegerRange(cv::Mat src,Point & bad_pt,int minVal,int maxVal)1414 static bool checkIntegerRange(cv::Mat src, Point& bad_pt, int minVal, int maxVal)
1415 {
1416 typedef mat_type_assotiations<depth> type_ass;
1417
1418 if (minVal < type_ass::min_allowable && maxVal > type_ass::max_allowable)
1419 {
1420 return true;
1421 }
1422 else if (minVal > type_ass::max_allowable || maxVal < type_ass::min_allowable || maxVal < minVal)
1423 {
1424 bad_pt = cv::Point(0,0);
1425 return false;
1426 }
1427 cv::Mat as_one_channel = src.reshape(1,0);
1428
1429 for (int j = 0; j < as_one_channel.rows; ++j)
1430 for (int i = 0; i < as_one_channel.cols; ++i)
1431 {
1432 typename type_ass::type v = as_one_channel.at<typename type_ass::type>(j ,i);
1433 if (v < minVal || v > maxVal)
1434 {
1435 bad_pt.y = j;
1436 bad_pt.x = i / src.channels();
1437 return false;
1438 }
1439 }
1440
1441 return true;
1442 }
1443
1444 typedef bool (*check_range_function)(cv::Mat src, Point& bad_pt, int minVal, int maxVal);
1445
1446 check_range_function check_range_functions[] =
1447 {
1448 &checkIntegerRange<CV_8U>,
1449 &checkIntegerRange<CV_8S>,
1450 &checkIntegerRange<CV_16U>,
1451 &checkIntegerRange<CV_16S>,
1452 &checkIntegerRange<CV_32S>
1453 };
1454
checkRange(InputArray _src,bool quiet,Point * pt,double minVal,double maxVal)1455 bool checkRange(InputArray _src, bool quiet, Point* pt, double minVal, double maxVal)
1456 {
1457 CV_INSTRUMENT_REGION();
1458
1459 Mat src = _src.getMat();
1460
1461 if ( src.dims > 2 )
1462 {
1463 CV_Assert(pt == NULL); // no way to provide location info
1464
1465 const Mat* arrays[] = {&src, 0};
1466 Mat planes[1];
1467 NAryMatIterator it(arrays, planes);
1468
1469 for ( size_t i = 0; i < it.nplanes; i++, ++it )
1470 {
1471 if (!checkRange( it.planes[0], quiet, NULL, minVal, maxVal ))
1472 {
1473 return false;
1474 }
1475 }
1476 return true;
1477 }
1478
1479 int depth = src.depth();
1480 Point badPt(-1, -1);
1481
1482 if (depth < CV_32F)
1483 {
1484 int minVali = minVal <= INT_MIN ? INT_MIN : cvFloor(minVal);
1485 int maxVali = maxVal > INT_MAX ? INT_MAX : cvCeil(maxVal) - 1;
1486
1487 (check_range_functions[depth])(src, badPt, minVali, maxVali);
1488 }
1489 else
1490 {
1491 int i, loc = 0;
1492 int cn = src.channels();
1493 Size size = getContinuousSize2D(src, cn);
1494
1495 if( depth == CV_32F )
1496 {
1497 Cv32suf a, b;
1498 int ia, ib;
1499 const int* isrc = src.ptr<int>();
1500 size_t step = src.step/sizeof(isrc[0]);
1501
1502 a.f = (float)std::max(minVal, (double)-FLT_MAX);
1503 b.f = (float)std::min(maxVal, (double)FLT_MAX);
1504
1505 ia = CV_TOGGLE_FLT(a.i);
1506 ib = CV_TOGGLE_FLT(b.i);
1507
1508 for( ; badPt.x < 0 && size.height--; loc += size.width, isrc += step )
1509 {
1510 for( i = 0; i < size.width; i++ )
1511 {
1512 int val = isrc[i];
1513 val = CV_TOGGLE_FLT(val);
1514
1515 if( val < ia || val >= ib )
1516 {
1517 int pixelId = (loc + i) / cn;
1518 badPt = Point(pixelId % src.cols, pixelId / src.cols);
1519 break;
1520 }
1521 }
1522 }
1523 }
1524 else
1525 {
1526 Cv64suf a, b;
1527 int64 ia, ib;
1528 const int64* isrc = src.ptr<int64>();
1529 size_t step = src.step/sizeof(isrc[0]);
1530
1531 a.f = minVal;
1532 b.f = maxVal;
1533
1534 ia = CV_TOGGLE_DBL(a.i);
1535 ib = CV_TOGGLE_DBL(b.i);
1536
1537 for( ; badPt.x < 0 && size.height--; loc += size.width, isrc += step )
1538 {
1539 for( i = 0; i < size.width; i++ )
1540 {
1541 int64 val = isrc[i];
1542 val = CV_TOGGLE_DBL(val);
1543
1544 if( val < ia || val >= ib )
1545 {
1546 int pixelId = (loc + i) / cn;
1547 badPt = Point(pixelId % src.cols, pixelId / src.cols);
1548 break;
1549 }
1550 }
1551 }
1552 }
1553 }
1554
1555 if( badPt.x >= 0 )
1556 {
1557 if( pt )
1558 *pt = badPt;
1559 if( !quiet )
1560 {
1561 cv::String value_str;
1562 value_str << src(cv::Range(badPt.y, badPt.y + 1), cv::Range(badPt.x, badPt.x + 1));
1563 CV_Error_( CV_StsOutOfRange,
1564 ("the value at (%d, %d)=%s is out of range [%f, %f)", badPt.x, badPt.y, value_str.c_str(), minVal, maxVal));
1565 }
1566 return false;
1567 }
1568 return true;
1569 }
1570
1571 #ifdef HAVE_OPENCL
1572
ocl_patchNaNs(InputOutputArray _a,float value)1573 static bool ocl_patchNaNs( InputOutputArray _a, float value )
1574 {
1575 int rowsPerWI = ocl::Device::getDefault().isIntel() ? 4 : 1;
1576 ocl::Kernel k("KF", ocl::core::arithm_oclsrc,
1577 format("-D UNARY_OP -D OP_PATCH_NANS -D dstT=float -D DEPTH_dst=%d -D rowsPerWI=%d",
1578 CV_32F, rowsPerWI));
1579 if (k.empty())
1580 return false;
1581
1582 UMat a = _a.getUMat();
1583 int cn = a.channels();
1584
1585 k.args(ocl::KernelArg::ReadOnlyNoSize(a),
1586 ocl::KernelArg::WriteOnly(a, cn), (float)value);
1587
1588 size_t globalsize[2] = { (size_t)a.cols * cn, ((size_t)a.rows + rowsPerWI - 1) / rowsPerWI };
1589 return k.run(2, globalsize, NULL, false);
1590 }
1591
1592 #endif
1593
patchNaNs(InputOutputArray _a,double _val)1594 void patchNaNs( InputOutputArray _a, double _val )
1595 {
1596 CV_INSTRUMENT_REGION();
1597
1598 CV_Assert( _a.depth() == CV_32F );
1599
1600 CV_OCL_RUN(_a.isUMat() && _a.dims() <= 2,
1601 ocl_patchNaNs(_a, (float)_val))
1602
1603 Mat a = _a.getMat();
1604 const Mat* arrays[] = {&a, 0};
1605 int* ptrs[1] = {};
1606 NAryMatIterator it(arrays, (uchar**)ptrs);
1607 size_t len = it.size*a.channels();
1608 Cv32suf val;
1609 val.f = (float)_val;
1610
1611 #if CV_SIMD
1612 v_int32 v_mask1 = vx_setall_s32(0x7fffffff), v_mask2 = vx_setall_s32(0x7f800000);
1613 v_int32 v_val = vx_setall_s32(val.i);
1614 #endif
1615
1616 for( size_t i = 0; i < it.nplanes; i++, ++it )
1617 {
1618 int* tptr = ptrs[0];
1619 size_t j = 0;
1620
1621 #if CV_SIMD
1622 size_t cWidth = (size_t)v_int32::nlanes;
1623 for ( ; j + cWidth <= len; j += cWidth)
1624 {
1625 v_int32 v_src = vx_load(tptr + j);
1626 v_int32 v_cmp_mask = v_mask2 < (v_src & v_mask1);
1627 v_int32 v_dst = v_select(v_cmp_mask, v_val, v_src);
1628 v_store(tptr + j, v_dst);
1629 }
1630 vx_cleanup();
1631 #endif
1632
1633 for( ; j < len; j++ )
1634 if( (tptr[j] & 0x7fffffff) > 0x7f800000 )
1635 tptr[j] = val.i;
1636 }
1637 }
1638
1639 }
1640
1641
1642 #ifndef OPENCV_EXCLUDE_C_API
1643
cvCbrt(float value)1644 CV_IMPL float cvCbrt(float value) { return cv::cubeRoot(value); }
cvFastArctan(float y,float x)1645 CV_IMPL float cvFastArctan(float y, float x) { return cv::fastAtan2(y, x); }
1646
1647 CV_IMPL void
cvCartToPolar(const CvArr * xarr,const CvArr * yarr,CvArr * magarr,CvArr * anglearr,int angle_in_degrees)1648 cvCartToPolar( const CvArr* xarr, const CvArr* yarr,
1649 CvArr* magarr, CvArr* anglearr,
1650 int angle_in_degrees )
1651 {
1652 cv::Mat X = cv::cvarrToMat(xarr), Y = cv::cvarrToMat(yarr), Mag, Angle;
1653 if( magarr )
1654 {
1655 Mag = cv::cvarrToMat(magarr);
1656 CV_Assert( Mag.size() == X.size() && Mag.type() == X.type() );
1657 }
1658 if( anglearr )
1659 {
1660 Angle = cv::cvarrToMat(anglearr);
1661 CV_Assert( Angle.size() == X.size() && Angle.type() == X.type() );
1662 }
1663 if( magarr )
1664 {
1665 if( anglearr )
1666 cv::cartToPolar( X, Y, Mag, Angle, angle_in_degrees != 0 );
1667 else
1668 cv::magnitude( X, Y, Mag );
1669 }
1670 else
1671 cv::phase( X, Y, Angle, angle_in_degrees != 0 );
1672 }
1673
1674 CV_IMPL void
cvPolarToCart(const CvArr * magarr,const CvArr * anglearr,CvArr * xarr,CvArr * yarr,int angle_in_degrees)1675 cvPolarToCart( const CvArr* magarr, const CvArr* anglearr,
1676 CvArr* xarr, CvArr* yarr, int angle_in_degrees )
1677 {
1678 cv::Mat X, Y, Angle = cv::cvarrToMat(anglearr), Mag;
1679 if( magarr )
1680 {
1681 Mag = cv::cvarrToMat(magarr);
1682 CV_Assert( Mag.size() == Angle.size() && Mag.type() == Angle.type() );
1683 }
1684 if( xarr )
1685 {
1686 X = cv::cvarrToMat(xarr);
1687 CV_Assert( X.size() == Angle.size() && X.type() == Angle.type() );
1688 }
1689 if( yarr )
1690 {
1691 Y = cv::cvarrToMat(yarr);
1692 CV_Assert( Y.size() == Angle.size() && Y.type() == Angle.type() );
1693 }
1694
1695 cv::polarToCart( Mag, Angle, X, Y, angle_in_degrees != 0 );
1696 }
1697
cvExp(const CvArr * srcarr,CvArr * dstarr)1698 CV_IMPL void cvExp( const CvArr* srcarr, CvArr* dstarr )
1699 {
1700 cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
1701 CV_Assert( src.type() == dst.type() && src.size == dst.size );
1702 cv::exp( src, dst );
1703 }
1704
cvLog(const CvArr * srcarr,CvArr * dstarr)1705 CV_IMPL void cvLog( const CvArr* srcarr, CvArr* dstarr )
1706 {
1707 cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
1708 CV_Assert( src.type() == dst.type() && src.size == dst.size );
1709 cv::log( src, dst );
1710 }
1711
cvPow(const CvArr * srcarr,CvArr * dstarr,double power)1712 CV_IMPL void cvPow( const CvArr* srcarr, CvArr* dstarr, double power )
1713 {
1714 cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
1715 CV_Assert( src.type() == dst.type() && src.size == dst.size );
1716 cv::pow( src, power, dst );
1717 }
1718
cvCheckArr(const CvArr * arr,int flags,double minVal,double maxVal)1719 CV_IMPL int cvCheckArr( const CvArr* arr, int flags,
1720 double minVal, double maxVal )
1721 {
1722 if( (flags & CV_CHECK_RANGE) == 0 )
1723 minVal = -DBL_MAX, maxVal = DBL_MAX;
1724 return cv::checkRange(cv::cvarrToMat(arr), (flags & CV_CHECK_QUIET) != 0, 0, minVal, maxVal );
1725 }
1726
1727 #endif // OPENCV_EXCLUDE_C_API
1728
1729 /*
1730 Finds real roots of cubic, quadratic or linear equation.
1731 The original code has been taken from Ken Turkowski web page
1732 (http://www.worldserver.com/turk/opensource/) and adopted for OpenCV.
1733 Here is the copyright notice.
1734
1735 -----------------------------------------------------------------------
1736 Copyright (C) 1978-1999 Ken Turkowski. <turk@computer.org>
1737
1738 All rights reserved.
1739
1740 Warranty Information
1741 Even though I have reviewed this software, I make no warranty
1742 or representation, either express or implied, with respect to this
1743 software, its quality, accuracy, merchantability, or fitness for a
1744 particular purpose. As a result, this software is provided "as is,"
1745 and you, its user, are assuming the entire risk as to its quality
1746 and accuracy.
1747
1748 This code may be used and freely distributed as long as it includes
1749 this copyright notice and the above warranty information.
1750 -----------------------------------------------------------------------
1751 */
1752
solveCubic(InputArray _coeffs,OutputArray _roots)1753 int cv::solveCubic( InputArray _coeffs, OutputArray _roots )
1754 {
1755 CV_INSTRUMENT_REGION();
1756
1757 const int n0 = 3;
1758 Mat coeffs = _coeffs.getMat();
1759 int ctype = coeffs.type();
1760
1761 CV_Assert( ctype == CV_32F || ctype == CV_64F );
1762 CV_Assert( (coeffs.size() == Size(n0, 1) ||
1763 coeffs.size() == Size(n0+1, 1) ||
1764 coeffs.size() == Size(1, n0) ||
1765 coeffs.size() == Size(1, n0+1)) );
1766
1767 _roots.create(n0, 1, ctype, -1, true, _OutputArray::DEPTH_MASK_FLT);
1768 Mat roots = _roots.getMat();
1769
1770 int i = -1, n = 0;
1771 double a0 = 1., a1, a2, a3;
1772 double x0 = 0., x1 = 0., x2 = 0.;
1773 int ncoeffs = coeffs.rows + coeffs.cols - 1;
1774
1775 if( ctype == CV_32FC1 )
1776 {
1777 if( ncoeffs == 4 )
1778 a0 = coeffs.at<float>(++i);
1779
1780 a1 = coeffs.at<float>(i+1);
1781 a2 = coeffs.at<float>(i+2);
1782 a3 = coeffs.at<float>(i+3);
1783 }
1784 else
1785 {
1786 if( ncoeffs == 4 )
1787 a0 = coeffs.at<double>(++i);
1788
1789 a1 = coeffs.at<double>(i+1);
1790 a2 = coeffs.at<double>(i+2);
1791 a3 = coeffs.at<double>(i+3);
1792 }
1793
1794 if( a0 == 0 )
1795 {
1796 if( a1 == 0 )
1797 {
1798 if( a2 == 0 )
1799 n = a3 == 0 ? -1 : 0;
1800 else
1801 {
1802 // linear equation
1803 x0 = -a3/a2;
1804 n = 1;
1805 }
1806 }
1807 else
1808 {
1809 // quadratic equation
1810 double d = a2*a2 - 4*a1*a3;
1811 if( d >= 0 )
1812 {
1813 d = std::sqrt(d);
1814 double q1 = (-a2 + d) * 0.5;
1815 double q2 = (a2 + d) * -0.5;
1816 if( fabs(q1) > fabs(q2) )
1817 {
1818 x0 = q1 / a1;
1819 x1 = a3 / q1;
1820 }
1821 else
1822 {
1823 x0 = q2 / a1;
1824 x1 = a3 / q2;
1825 }
1826 n = d > 0 ? 2 : 1;
1827 }
1828 }
1829 }
1830 else
1831 {
1832 a0 = 1./a0;
1833 a1 *= a0;
1834 a2 *= a0;
1835 a3 *= a0;
1836
1837 double Q = (a1 * a1 - 3 * a2) * (1./9);
1838 double R = (2 * a1 * a1 * a1 - 9 * a1 * a2 + 27 * a3) * (1./54);
1839 double Qcubed = Q * Q * Q;
1840 double d = Qcubed - R * R;
1841
1842 if( d > 0 )
1843 {
1844 double theta = acos(R / sqrt(Qcubed));
1845 double sqrtQ = sqrt(Q);
1846 double t0 = -2 * sqrtQ;
1847 double t1 = theta * (1./3);
1848 double t2 = a1 * (1./3);
1849 x0 = t0 * cos(t1) - t2;
1850 x1 = t0 * cos(t1 + (2.*CV_PI/3)) - t2;
1851 x2 = t0 * cos(t1 + (4.*CV_PI/3)) - t2;
1852 n = 3;
1853 }
1854 else if( d == 0 )
1855 {
1856 if(R >= 0)
1857 {
1858 x0 = -2*pow(R, 1./3) - a1/3;
1859 x1 = pow(R, 1./3) - a1/3;
1860 }
1861 else
1862 {
1863 x0 = 2*pow(-R, 1./3) - a1/3;
1864 x1 = -pow(-R, 1./3) - a1/3;
1865 }
1866 x2 = 0;
1867 n = x0 == x1 ? 1 : 2;
1868 x1 = x0 == x1 ? 0 : x1;
1869 }
1870 else
1871 {
1872 double e;
1873 d = sqrt(-d);
1874 e = pow(d + fabs(R), 1./3);
1875 if( R > 0 )
1876 e = -e;
1877 x0 = (e + Q / e) - a1 * (1./3);
1878 n = 1;
1879 }
1880 }
1881
1882 if( roots.type() == CV_32FC1 )
1883 {
1884 roots.at<float>(0) = (float)x0;
1885 roots.at<float>(1) = (float)x1;
1886 roots.at<float>(2) = (float)x2;
1887 }
1888 else
1889 {
1890 roots.at<double>(0) = x0;
1891 roots.at<double>(1) = x1;
1892 roots.at<double>(2) = x2;
1893 }
1894
1895 return n;
1896 }
1897
1898 /* finds complex roots of a polynomial using Durand-Kerner method:
1899 http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method */
solvePoly(InputArray _coeffs0,OutputArray _roots0,int maxIters)1900 double cv::solvePoly( InputArray _coeffs0, OutputArray _roots0, int maxIters )
1901 {
1902 CV_INSTRUMENT_REGION();
1903
1904 typedef Complex<double> C;
1905
1906 double maxDiff = 0;
1907 int iter, i, j;
1908 Mat coeffs0 = _coeffs0.getMat();
1909 int ctype = _coeffs0.type();
1910 int cdepth = CV_MAT_DEPTH(ctype);
1911
1912 CV_Assert( CV_MAT_DEPTH(ctype) >= CV_32F && CV_MAT_CN(ctype) <= 2 );
1913 CV_Assert( coeffs0.rows == 1 || coeffs0.cols == 1 );
1914
1915 int n0 = coeffs0.cols + coeffs0.rows - 2, n = n0;
1916
1917 _roots0.create(n, 1, CV_MAKETYPE(cdepth, 2), -1, true, _OutputArray::DEPTH_MASK_FLT);
1918 Mat roots0 = _roots0.getMat();
1919
1920 AutoBuffer<C> buf(n*2+2);
1921 C *coeffs = buf.data(), *roots = coeffs + n + 1;
1922 Mat coeffs1(coeffs0.size(), CV_MAKETYPE(CV_64F, coeffs0.channels()), coeffs0.channels() == 2 ? coeffs : roots);
1923 coeffs0.convertTo(coeffs1, coeffs1.type());
1924 if( coeffs0.channels() == 1 )
1925 {
1926 const double* rcoeffs = (const double*)roots;
1927 for( i = 0; i <= n; i++ )
1928 coeffs[i] = C(rcoeffs[i], 0);
1929 }
1930
1931 for( ; n > 1; n-- )
1932 {
1933 if( std::abs(coeffs[n].re) + std::abs(coeffs[n].im) > DBL_EPSILON )
1934 break;
1935 }
1936
1937 C p(1, 0), r(1, 1);
1938
1939 for( i = 0; i < n; i++ )
1940 {
1941 roots[i] = p;
1942 p = p * r;
1943 }
1944
1945 maxIters = maxIters <= 0 ? 1000 : maxIters;
1946 for( iter = 0; iter < maxIters; iter++ )
1947 {
1948 maxDiff = 0;
1949 for( i = 0; i < n; i++ )
1950 {
1951 p = roots[i];
1952 C num = coeffs[n], denom = coeffs[n];
1953 int num_same_root = 1;
1954 for( j = 0; j < n; j++ )
1955 {
1956 num = num*p + coeffs[n-j-1];
1957 if( j != i )
1958 {
1959 if ( (p - roots[j]).re != 0 || (p - roots[j]).im != 0 )
1960 denom = denom * (p - roots[j]);
1961 else
1962 num_same_root++;
1963 }
1964 }
1965 num /= denom;
1966 if( num_same_root > 1)
1967 {
1968 double old_num_re = num.re;
1969 double old_num_im = num.im;
1970 int square_root_times = num_same_root % 2 == 0 ? num_same_root / 2 : num_same_root / 2 - 1;
1971
1972 for( j = 0; j < square_root_times; j++)
1973 {
1974 num.re = old_num_re*old_num_re + old_num_im*old_num_im;
1975 num.re = sqrt(num.re);
1976 num.re += old_num_re;
1977 num.im = num.re - old_num_re;
1978 num.re /= 2;
1979 num.re = sqrt(num.re);
1980
1981 num.im /= 2;
1982 num.im = sqrt(num.im);
1983 if( old_num_re < 0 ) num.im = -num.im;
1984 }
1985
1986 if( num_same_root % 2 != 0){
1987 Mat cube_coefs(4, 1, CV_64FC1);
1988 Mat cube_roots(3, 1, CV_64FC2);
1989 cube_coefs.at<double>(3) = -(pow(old_num_re, 3));
1990 cube_coefs.at<double>(2) = -(15*pow(old_num_re, 2) + 27*pow(old_num_im, 2));
1991 cube_coefs.at<double>(1) = -48*old_num_re;
1992 cube_coefs.at<double>(0) = 64;
1993 solveCubic(cube_coefs, cube_roots);
1994
1995 if(cube_roots.at<double>(0) >= 0) num.re = pow(cube_roots.at<double>(0), 1./3);
1996 else num.re = -pow(-cube_roots.at<double>(0), 1./3);
1997 num.im = sqrt(pow(num.re, 2) / 3 - old_num_re / (3*num.re));
1998 }
1999 }
2000 roots[i] = p - num;
2001 maxDiff = std::max(maxDiff, cv::abs(num));
2002 }
2003 if( maxDiff <= 0 )
2004 break;
2005 }
2006
2007 if( coeffs0.channels() == 1 )
2008 {
2009 const double verySmallEps = 1e-100;
2010 for( i = 0; i < n; i++ )
2011 if( fabs(roots[i].im) < verySmallEps )
2012 roots[i].im = 0;
2013 }
2014
2015 for( ; n < n0; n++ )
2016 roots[n+1] = roots[n];
2017
2018 Mat(roots0.size(), CV_64FC2, roots).convertTo(roots0, roots0.type());
2019 return maxDiff;
2020 }
2021
2022
2023 #ifndef OPENCV_EXCLUDE_C_API
2024
2025 CV_IMPL int
cvSolveCubic(const CvMat * coeffs,CvMat * roots)2026 cvSolveCubic( const CvMat* coeffs, CvMat* roots )
2027 {
2028 cv::Mat _coeffs = cv::cvarrToMat(coeffs), _roots = cv::cvarrToMat(roots), _roots0 = _roots;
2029 int nroots = cv::solveCubic(_coeffs, _roots);
2030 CV_Assert( _roots.data == _roots0.data ); // check that the array of roots was not reallocated
2031 return nroots;
2032 }
2033
2034
cvSolvePoly(const CvMat * a,CvMat * r,int maxiter,int)2035 void cvSolvePoly(const CvMat* a, CvMat *r, int maxiter, int)
2036 {
2037 cv::Mat _a = cv::cvarrToMat(a);
2038 cv::Mat _r = cv::cvarrToMat(r);
2039 cv::Mat _r0 = _r;
2040 cv::solvePoly(_a, _r, maxiter);
2041 CV_Assert( _r.data == _r0.data ); // check that the array of roots was not reallocated
2042 }
2043
2044 #endif // OPENCV_EXCLUDE_C_API
2045
2046
2047 // Common constants for dispatched code
2048 namespace cv { namespace details {
2049
2050 #define EXPTAB_SCALE 6
2051 #define EXPTAB_MASK ((1 << EXPTAB_SCALE) - 1)
2052
2053 #define EXPPOLY_32F_A0 .9670371139572337719125840413672004409288e-2
2054
2055 static const double CV_DECL_ALIGNED(64) expTab[EXPTAB_MASK + 1] = {
2056 1.0 * EXPPOLY_32F_A0,
2057 1.0108892860517004600204097905619 * EXPPOLY_32F_A0,
2058 1.0218971486541166782344801347833 * EXPPOLY_32F_A0,
2059 1.0330248790212284225001082839705 * EXPPOLY_32F_A0,
2060 1.0442737824274138403219664787399 * EXPPOLY_32F_A0,
2061 1.0556451783605571588083413251529 * EXPPOLY_32F_A0,
2062 1.0671404006768236181695211209928 * EXPPOLY_32F_A0,
2063 1.0787607977571197937406800374385 * EXPPOLY_32F_A0,
2064 1.0905077326652576592070106557607 * EXPPOLY_32F_A0,
2065 1.1023825833078409435564142094256 * EXPPOLY_32F_A0,
2066 1.1143867425958925363088129569196 * EXPPOLY_32F_A0,
2067 1.126521618608241899794798643787 * EXPPOLY_32F_A0,
2068 1.1387886347566916537038302838415 * EXPPOLY_32F_A0,
2069 1.151189229952982705817759635202 * EXPPOLY_32F_A0,
2070 1.1637248587775775138135735990922 * EXPPOLY_32F_A0,
2071 1.1763969916502812762846457284838 * EXPPOLY_32F_A0,
2072 1.1892071150027210667174999705605 * EXPPOLY_32F_A0,
2073 1.2021567314527031420963969574978 * EXPPOLY_32F_A0,
2074 1.2152473599804688781165202513388 * EXPPOLY_32F_A0,
2075 1.2284805361068700056940089577928 * EXPPOLY_32F_A0,
2076 1.2418578120734840485936774687266 * EXPPOLY_32F_A0,
2077 1.2553807570246910895793906574423 * EXPPOLY_32F_A0,
2078 1.2690509571917332225544190810323 * EXPPOLY_32F_A0,
2079 1.2828700160787782807266697810215 * EXPPOLY_32F_A0,
2080 1.2968395546510096659337541177925 * EXPPOLY_32F_A0,
2081 1.3109612115247643419229917863308 * EXPPOLY_32F_A0,
2082 1.3252366431597412946295370954987 * EXPPOLY_32F_A0,
2083 1.3396675240533030053600306697244 * EXPPOLY_32F_A0,
2084 1.3542555469368927282980147401407 * EXPPOLY_32F_A0,
2085 1.3690024229745906119296011329822 * EXPPOLY_32F_A0,
2086 1.3839098819638319548726595272652 * EXPPOLY_32F_A0,
2087 1.3989796725383111402095281367152 * EXPPOLY_32F_A0,
2088 1.4142135623730950488016887242097 * EXPPOLY_32F_A0,
2089 1.4296133383919700112350657782751 * EXPPOLY_32F_A0,
2090 1.4451808069770466200370062414717 * EXPPOLY_32F_A0,
2091 1.4609177941806469886513028903106 * EXPPOLY_32F_A0,
2092 1.476826145939499311386907480374 * EXPPOLY_32F_A0,
2093 1.4929077282912648492006435314867 * EXPPOLY_32F_A0,
2094 1.5091644275934227397660195510332 * EXPPOLY_32F_A0,
2095 1.5255981507445383068512536895169 * EXPPOLY_32F_A0,
2096 1.5422108254079408236122918620907 * EXPPOLY_32F_A0,
2097 1.5590044002378369670337280894749 * EXPPOLY_32F_A0,
2098 1.5759808451078864864552701601819 * EXPPOLY_32F_A0,
2099 1.5931421513422668979372486431191 * EXPPOLY_32F_A0,
2100 1.6104903319492543081795206673574 * EXPPOLY_32F_A0,
2101 1.628027421857347766848218522014 * EXPPOLY_32F_A0,
2102 1.6457554781539648445187567247258 * EXPPOLY_32F_A0,
2103 1.6636765803267364350463364569764 * EXPPOLY_32F_A0,
2104 1.6817928305074290860622509524664 * EXPPOLY_32F_A0,
2105 1.7001063537185234695013625734975 * EXPPOLY_32F_A0,
2106 1.7186192981224779156293443764563 * EXPPOLY_32F_A0,
2107 1.7373338352737062489942020818722 * EXPPOLY_32F_A0,
2108 1.7562521603732994831121606193753 * EXPPOLY_32F_A0,
2109 1.7753764925265212525505592001993 * EXPPOLY_32F_A0,
2110 1.7947090750031071864277032421278 * EXPPOLY_32F_A0,
2111 1.8142521755003987562498346003623 * EXPPOLY_32F_A0,
2112 1.8340080864093424634870831895883 * EXPPOLY_32F_A0,
2113 1.8539791250833855683924530703377 * EXPPOLY_32F_A0,
2114 1.8741676341102999013299989499544 * EXPPOLY_32F_A0,
2115 1.8945759815869656413402186534269 * EXPPOLY_32F_A0,
2116 1.9152065613971472938726112702958 * EXPPOLY_32F_A0,
2117 1.9360617934922944505980559045667 * EXPPOLY_32F_A0,
2118 1.9571441241754002690183222516269 * EXPPOLY_32F_A0,
2119 1.9784560263879509682582499181312 * EXPPOLY_32F_A0,
2120 };
2121
getExpTab64f()2122 const double* getExpTab64f()
2123 {
2124 return expTab;
2125 }
2126
getExpTab32f()2127 const float* getExpTab32f()
2128 {
2129 static float CV_DECL_ALIGNED(64) expTab_f[EXPTAB_MASK+1];
2130 static std::atomic<bool> expTab_f_initialized(false);
2131 if (!expTab_f_initialized.load())
2132 {
2133 for( int j = 0; j <= EXPTAB_MASK; j++ )
2134 expTab_f[j] = (float)expTab[j];
2135 expTab_f_initialized = true;
2136 }
2137 return expTab_f;
2138 }
2139
2140
2141
2142 #define LOGTAB_SCALE 8
2143 #define LOGTAB_MASK ((1 << LOGTAB_SCALE) - 1)
2144
2145 static const double CV_DECL_ALIGNED(64) logTab[(LOGTAB_MASK+1)*2] = {
2146 0.0000000000000000000000000000000000000000, 1.000000000000000000000000000000000000000,
2147 .00389864041565732288852075271279318258166, .9961089494163424124513618677042801556420,
2148 .00778214044205494809292034119607706088573, .9922480620155038759689922480620155038760,
2149 .01165061721997527263705585198749759001657, .9884169884169884169884169884169884169884,
2150 .01550418653596525274396267235488267033361, .9846153846153846153846153846153846153846,
2151 .01934296284313093139406447562578250654042, .9808429118773946360153256704980842911877,
2152 .02316705928153437593630670221500622574241, .9770992366412213740458015267175572519084,
2153 .02697658769820207233514075539915211265906, .9733840304182509505703422053231939163498,
2154 .03077165866675368732785500469617545604706, .9696969696969696969696969696969696969697,
2155 .03455238150665972812758397481047722976656, .9660377358490566037735849056603773584906,
2156 .03831886430213659461285757856785494368522, .9624060150375939849624060150375939849624,
2157 .04207121392068705056921373852674150839447, .9588014981273408239700374531835205992509,
2158 .04580953603129420126371940114040626212953, .9552238805970149253731343283582089552239,
2159 .04953393512227662748292900118940451648088, .9516728624535315985130111524163568773234,
2160 .05324451451881227759255210685296333394944, .9481481481481481481481481481481481481481,
2161 .05694137640013842427411105973078520037234, .9446494464944649446494464944649446494465,
2162 .06062462181643483993820353816772694699466, .9411764705882352941176470588235294117647,
2163 .06429435070539725460836422143984236754475, .9377289377289377289377289377289377289377,
2164 .06795066190850773679699159401934593915938, .9343065693430656934306569343065693430657,
2165 .07159365318700880442825962290953611955044, .9309090909090909090909090909090909090909,
2166 .07522342123758751775142172846244648098944, .9275362318840579710144927536231884057971,
2167 .07884006170777602129362549021607264876369, .9241877256317689530685920577617328519856,
2168 .08244366921107458556772229485432035289706, .9208633093525179856115107913669064748201,
2169 .08603433734180314373940490213499288074675, .9175627240143369175627240143369175627240,
2170 .08961215868968712416897659522874164395031, .9142857142857142857142857142857142857143,
2171 .09317722485418328259854092721070628613231, .9110320284697508896797153024911032028470,
2172 .09672962645855109897752299730200320482256, .9078014184397163120567375886524822695035,
2173 .10026945316367513738597949668474029749630, .9045936395759717314487632508833922261484,
2174 .10379679368164355934833764649738441221420, .9014084507042253521126760563380281690141,
2175 .10731173578908805021914218968959175981580, .8982456140350877192982456140350877192982,
2176 .11081436634029011301105782649756292812530, .8951048951048951048951048951048951048951,
2177 .11430477128005862852422325204315711744130, .8919860627177700348432055749128919860627,
2178 .11778303565638344185817487641543266363440, .8888888888888888888888888888888888888889,
2179 .12124924363286967987640707633545389398930, .8858131487889273356401384083044982698962,
2180 .12470347850095722663787967121606925502420, .8827586206896551724137931034482758620690,
2181 .12814582269193003360996385708858724683530, .8797250859106529209621993127147766323024,
2182 .13157635778871926146571524895989568904040, .8767123287671232876712328767123287671233,
2183 .13499516453750481925766280255629681050780, .8737201365187713310580204778156996587031,
2184 .13840232285911913123754857224412262439730, .8707482993197278911564625850340136054422,
2185 .14179791186025733629172407290752744302150, .8677966101694915254237288135593220338983,
2186 .14518200984449788903951628071808954700830, .8648648648648648648648648648648648648649,
2187 .14855469432313711530824207329715136438610, .8619528619528619528619528619528619528620,
2188 .15191604202584196858794030049466527998450, .8590604026845637583892617449664429530201,
2189 .15526612891112392955683674244937719777230, .8561872909698996655518394648829431438127,
2190 .15860503017663857283636730244325008243330, .8533333333333333333333333333333333333333,
2191 .16193282026931324346641360989451641216880, .8504983388704318936877076411960132890365,
2192 .16524957289530714521497145597095368430010, .8476821192052980132450331125827814569536,
2193 .16855536102980664403538924034364754334090, .8448844884488448844884488448844884488449,
2194 .17185025692665920060697715143760433420540, .8421052631578947368421052631578947368421,
2195 .17513433212784912385018287750426679849630, .8393442622950819672131147540983606557377,
2196 .17840765747281828179637841458315961062910, .8366013071895424836601307189542483660131,
2197 .18167030310763465639212199675966985523700, .8338762214983713355048859934853420195440,
2198 .18492233849401198964024217730184318497780, .8311688311688311688311688311688311688312,
2199 .18816383241818296356839823602058459073300, .8284789644012944983818770226537216828479,
2200 .19139485299962943898322009772527962923050, .8258064516129032258064516129032258064516,
2201 .19461546769967164038916962454095482826240, .8231511254019292604501607717041800643087,
2202 .19782574332991986754137769821682013571260, .8205128205128205128205128205128205128205,
2203 .20102574606059073203390141770796617493040, .8178913738019169329073482428115015974441,
2204 .20421554142869088876999228432396193966280, .8152866242038216560509554140127388535032,
2205 .20739519434607056602715147164417430758480, .8126984126984126984126984126984126984127,
2206 .21056476910734961416338251183333341032260, .8101265822784810126582278481012658227848,
2207 .21372432939771812687723695489694364368910, .8075709779179810725552050473186119873817,
2208 .21687393830061435506806333251006435602900, .8050314465408805031446540880503144654088,
2209 .22001365830528207823135744547471404075630, .8025078369905956112852664576802507836991,
2210 .22314355131420973710199007200571941211830, .8000000000000000000000000000000000000000,
2211 .22626367865045338145790765338460914790630, .7975077881619937694704049844236760124611,
2212 .22937410106484582006380890106811420992010, .7950310559006211180124223602484472049689,
2213 .23247487874309405442296849741978803649550, .7925696594427244582043343653250773993808,
2214 .23556607131276688371634975283086532726890, .7901234567901234567901234567901234567901,
2215 .23864773785017498464178231643018079921600, .7876923076923076923076923076923076923077,
2216 .24171993688714515924331749374687206000090, .7852760736196319018404907975460122699387,
2217 .24478272641769091566565919038112042471760, .7828746177370030581039755351681957186544,
2218 .24783616390458124145723672882013488560910, .7804878048780487804878048780487804878049,
2219 .25088030628580937353433455427875742316250, .7781155015197568389057750759878419452888,
2220 .25391520998096339667426946107298135757450, .7757575757575757575757575757575757575758,
2221 .25694093089750041913887912414793390780680, .7734138972809667673716012084592145015106,
2222 .25995752443692604627401010475296061486000, .7710843373493975903614457831325301204819,
2223 .26296504550088134477547896494797896593800, .7687687687687687687687687687687687687688,
2224 .26596354849713793599974565040611196309330, .7664670658682634730538922155688622754491,
2225 .26895308734550393836570947314612567424780, .7641791044776119402985074626865671641791,
2226 .27193371548364175804834985683555714786050, .7619047619047619047619047619047619047619,
2227 .27490548587279922676529508862586226314300, .7596439169139465875370919881305637982196,
2228 .27786845100345625159121709657483734190480, .7573964497041420118343195266272189349112,
2229 .28082266290088775395616949026589281857030, .7551622418879056047197640117994100294985,
2230 .28376817313064456316240580235898960381750, .7529411764705882352941176470588235294118,
2231 .28670503280395426282112225635501090437180, .7507331378299120234604105571847507331378,
2232 .28963329258304265634293983566749375313530, .7485380116959064327485380116959064327485,
2233 .29255300268637740579436012922087684273730, .7463556851311953352769679300291545189504,
2234 .29546421289383584252163927885703742504130, .7441860465116279069767441860465116279070,
2235 .29836697255179722709783618483925238251680, .7420289855072463768115942028985507246377,
2236 .30126133057816173455023545102449133992200, .7398843930635838150289017341040462427746,
2237 .30414733546729666446850615102448500692850, .7377521613832853025936599423631123919308,
2238 .30702503529491181888388950937951449304830, .7356321839080459770114942528735632183908,
2239 .30989447772286465854207904158101882785550, .7335243553008595988538681948424068767908,
2240 .31275571000389684739317885942000430077330, .7314285714285714285714285714285714285714,
2241 .31560877898630329552176476681779604405180, .7293447293447293447293447293447293447293,
2242 .31845373111853458869546784626436419785030, .7272727272727272727272727272727272727273,
2243 .32129061245373424782201254856772720813750, .7252124645892351274787535410764872521246,
2244 .32411946865421192853773391107097268104550, .7231638418079096045197740112994350282486,
2245 .32694034499585328257253991068864706903700, .7211267605633802816901408450704225352113,
2246 .32975328637246797969240219572384376078850, .7191011235955056179775280898876404494382,
2247 .33255833730007655635318997155991382896900, .7170868347338935574229691876750700280112,
2248 .33535554192113781191153520921943709254280, .7150837988826815642458100558659217877095,
2249 .33814494400871636381467055798566434532400, .7130919220055710306406685236768802228412,
2250 .34092658697059319283795275623560883104800, .7111111111111111111111111111111111111111,
2251 .34370051385331840121395430287520866841080, .7091412742382271468144044321329639889197,
2252 .34646676734620857063262633346312213689100, .7071823204419889502762430939226519337017,
2253 .34922538978528827602332285096053965389730, .7052341597796143250688705234159779614325,
2254 .35197642315717814209818925519357435405250, .7032967032967032967032967032967032967033,
2255 .35471990910292899856770532096561510115850, .7013698630136986301369863013698630136986,
2256 .35745588892180374385176833129662554711100, .6994535519125683060109289617486338797814,
2257 .36018440357500774995358483465679455548530, .6975476839237057220708446866485013623978,
2258 .36290549368936841911903457003063522279280, .6956521739130434782608695652173913043478,
2259 .36561919956096466943762379742111079394830, .6937669376693766937669376693766937669377,
2260 .36832556115870762614150635272380895912650, .6918918918918918918918918918918918918919,
2261 .37102461812787262962487488948681857436900, .6900269541778975741239892183288409703504,
2262 .37371640979358405898480555151763837784530, .6881720430107526881720430107526881720430,
2263 .37640097516425302659470730759494472295050, .6863270777479892761394101876675603217158,
2264 .37907835293496944251145919224654790014030, .6844919786096256684491978609625668449198,
2265 .38174858149084833769393299007788300514230, .6826666666666666666666666666666666666667,
2266 .38441169891033200034513583887019194662580, .6808510638297872340425531914893617021277,
2267 .38706774296844825844488013899535872042180, .6790450928381962864721485411140583554377,
2268 .38971675114002518602873692543653305619950, .6772486772486772486772486772486772486772,
2269 .39235876060286384303665840889152605086580, .6754617414248021108179419525065963060686,
2270 .39499380824086893770896722344332374632350, .6736842105263157894736842105263157894737,
2271 .39762193064713846624158577469643205404280, .6719160104986876640419947506561679790026,
2272 .40024316412701266276741307592601515352730, .6701570680628272251308900523560209424084,
2273 .40285754470108348090917615991202183067800, .6684073107049608355091383812010443864230,
2274 .40546510810816432934799991016916465014230, .6666666666666666666666666666666666666667,
2275 .40806588980822172674223224930756259709600, .6649350649350649350649350649350649350649,
2276 .41065992498526837639616360320360399782650, .6632124352331606217616580310880829015544,
2277 .41324724855021932601317757871584035456180, .6614987080103359173126614987080103359173,
2278 .41582789514371093497757669865677598863850, .6597938144329896907216494845360824742268,
2279 .41840189913888381489925905043492093682300, .6580976863753213367609254498714652956298,
2280 .42096929464412963239894338585145305842150, .6564102564102564102564102564102564102564,
2281 .42353011550580327293502591601281892508280, .6547314578005115089514066496163682864450,
2282 .42608439531090003260516141381231136620050, .6530612244897959183673469387755102040816,
2283 .42863216738969872610098832410585600882780, .6513994910941475826972010178117048346056,
2284 .43117346481837132143866142541810404509300, .6497461928934010152284263959390862944162,
2285 .43370832042155937902094819946796633303180, .6481012658227848101265822784810126582278,
2286 .43623676677491801667585491486534010618930, .6464646464646464646464646464646464646465,
2287 .43875883620762790027214350629947148263450, .6448362720403022670025188916876574307305,
2288 .44127456080487520440058801796112675219780, .6432160804020100502512562814070351758794,
2289 .44378397241030093089975139264424797147500, .6416040100250626566416040100250626566416,
2290 .44628710262841947420398014401143882423650, .6400000000000000000000000000000000000000,
2291 .44878398282700665555822183705458883196130, .6384039900249376558603491271820448877805,
2292 .45127464413945855836729492693848442286250, .6368159203980099502487562189054726368159,
2293 .45375911746712049854579618113348260521900, .6352357320099255583126550868486352357320,
2294 .45623743348158757315857769754074979573500, .6336633663366336633663366336633663366337,
2295 .45870962262697662081833982483658473938700, .6320987654320987654320987654320987654321,
2296 .46117571512217014895185229761409573256980, .6305418719211822660098522167487684729064,
2297 .46363574096303250549055974261136725544930, .6289926289926289926289926289926289926290,
2298 .46608972992459918316399125615134835243230, .6274509803921568627450980392156862745098,
2299 .46853771156323925639597405279346276074650, .6259168704156479217603911980440097799511,
2300 .47097971521879100631480241645476780831830, .6243902439024390243902439024390243902439,
2301 .47341577001667212165614273544633761048330, .6228710462287104622871046228710462287105,
2302 .47584590486996386493601107758877333253630, .6213592233009708737864077669902912621359,
2303 .47827014848147025860569669930555392056700, .6198547215496368038740920096852300242131,
2304 .48068852934575190261057286988943815231330, .6183574879227053140096618357487922705314,
2305 .48310107575113581113157579238759353756900, .6168674698795180722891566265060240963855,
2306 .48550781578170076890899053978500887751580, .6153846153846153846153846153846153846154,
2307 .48790877731923892879351001283794175833480, .6139088729016786570743405275779376498801,
2308 .49030398804519381705802061333088204264650, .6124401913875598086124401913875598086124,
2309 .49269347544257524607047571407747454941280, .6109785202863961813842482100238663484487,
2310 .49507726679785146739476431321236304938800, .6095238095238095238095238095238095238095,
2311 .49745538920281889838648226032091770321130, .6080760095011876484560570071258907363420,
2312 .49982786955644931126130359189119189977650, .6066350710900473933649289099526066350711,
2313 .50219473456671548383667413872899487614650, .6052009456264775413711583924349881796690,
2314 .50455601075239520092452494282042607665050, .6037735849056603773584905660377358490566,
2315 .50691172444485432801997148999362252652650, .6023529411764705882352941176470588235294,
2316 .50926190178980790257412536448100581765150, .6009389671361502347417840375586854460094,
2317 .51160656874906207391973111953120678663250, .5995316159250585480093676814988290398126,
2318 .51394575110223428282552049495279788970950, .5981308411214953271028037383177570093458,
2319 .51627947444845445623684554448118433356300, .5967365967365967365967365967365967365967,
2320 .51860776420804555186805373523384332656850, .5953488372093023255813953488372093023256,
2321 .52093064562418522900344441950437612831600, .5939675174013921113689095127610208816705,
2322 .52324814376454775732838697877014055848100, .5925925925925925925925925925925925925926,
2323 .52556028352292727401362526507000438869000, .5912240184757505773672055427251732101617,
2324 .52786708962084227803046587723656557500350, .5898617511520737327188940092165898617512,
2325 .53016858660912158374145519701414741575700, .5885057471264367816091954022988505747126,
2326 .53246479886947173376654518506256863474850, .5871559633027522935779816513761467889908,
2327 .53475575061602764748158733709715306758900, .5858123569794050343249427917620137299771,
2328 .53704146589688361856929077475797384977350, .5844748858447488584474885844748858447489,
2329 .53932196859560876944783558428753167390800, .5831435079726651480637813211845102505695,
2330 .54159728243274429804188230264117009937750, .5818181818181818181818181818181818181818,
2331 .54386743096728351609669971367111429572100, .5804988662131519274376417233560090702948,
2332 .54613243759813556721383065450936555862450, .5791855203619909502262443438914027149321,
2333 .54839232556557315767520321969641372561450, .5778781038374717832957110609480812641084,
2334 .55064711795266219063194057525834068655950, .5765765765765765765765765765765765765766,
2335 .55289683768667763352766542084282264113450, .5752808988764044943820224719101123595506,
2336 .55514150754050151093110798683483153581600, .5739910313901345291479820627802690582960,
2337 .55738115013400635344709144192165695130850, .5727069351230425055928411633109619686801,
2338 .55961578793542265941596269840374588966350, .5714285714285714285714285714285714285714,
2339 .56184544326269181269140062795486301183700, .5701559020044543429844097995545657015590,
2340 .56407013828480290218436721261241473257550, .5688888888888888888888888888888888888889,
2341 .56628989502311577464155334382667206227800, .5676274944567627494456762749445676274945,
2342 .56850473535266865532378233183408156037350, .5663716814159292035398230088495575221239,
2343 .57071468100347144680739575051120482385150, .5651214128035320088300220750551876379691,
2344 .57291975356178548306473885531886480748650, .5638766519823788546255506607929515418502,
2345 .57511997447138785144460371157038025558000, .5626373626373626373626373626373626373626,
2346 .57731536503482350219940144597785547375700, .5614035087719298245614035087719298245614,
2347 .57950594641464214795689713355386629700650, .5601750547045951859956236323851203501094,
2348 .58169173963462239562716149521293118596100, .5589519650655021834061135371179039301310,
2349 .58387276558098266665552955601015128195300, .5577342047930283224400871459694989106754,
2350 .58604904500357812846544902640744112432000, .5565217391304347826086956521739130434783,
2351 .58822059851708596855957011939608491957200, .5553145336225596529284164859002169197397,
2352 .59038744660217634674381770309992134571100, .5541125541125541125541125541125541125541,
2353 .59254960960667157898740242671919986605650, .5529157667386609071274298056155507559395,
2354 .59470710774669277576265358220553025603300, .5517241379310344827586206896551724137931,
2355 .59685996110779382384237123915227130055450, .5505376344086021505376344086021505376344,
2356 .59900818964608337768851242799428291618800, .5493562231759656652360515021459227467811,
2357 .60115181318933474940990890900138765573500, .5481798715203426124197002141327623126338,
2358 .60329085143808425240052883964381180703650, .5470085470085470085470085470085470085470,
2359 .60542532396671688843525771517306566238400, .5458422174840085287846481876332622601279,
2360 .60755525022454170969155029524699784815300, .5446808510638297872340425531914893617021,
2361 .60968064953685519036241657886421307921400, .5435244161358811040339702760084925690021,
2362 .61180154110599282990534675263916142284850, .5423728813559322033898305084745762711864,
2363 .61391794401237043121710712512140162289150, .5412262156448202959830866807610993657505,
2364 .61602987721551394351138242200249806046500, .5400843881856540084388185654008438818565,
2365 .61813735955507864705538167982012964785100, .5389473684210526315789473684210526315789,
2366 .62024040975185745772080281312810257077200, .5378151260504201680672268907563025210084,
2367 .62233904640877868441606324267922900617100, .5366876310272536687631027253668763102725,
2368 .62443328801189346144440150965237990021700, .5355648535564853556485355648535564853556,
2369 .62652315293135274476554741340805776417250, .5344467640918580375782881002087682672234,
2370 .62860865942237409420556559780379757285100, .5333333333333333333333333333333333333333,
2371 .63068982562619868570408243613201193511500, .5322245322245322245322245322245322245322,
2372 .63276666957103777644277897707070223987100, .5311203319502074688796680497925311203320,
2373 .63483920917301017716738442686619237065300, .5300207039337474120082815734989648033126,
2374 .63690746223706917739093569252872839570050, .5289256198347107438016528925619834710744,
2375 .63897144645792069983514238629140891134750, .5278350515463917525773195876288659793814,
2376 .64103117942093124081992527862894348800200, .5267489711934156378600823045267489711934,
2377 .64308667860302726193566513757104985415950, .5256673511293634496919917864476386036961,
2378 .64513796137358470073053240412264131009600, .5245901639344262295081967213114754098361,
2379 .64718504499530948859131740391603671014300, .5235173824130879345603271983640081799591,
2380 .64922794662510974195157587018911726772800, .5224489795918367346938775510204081632653,
2381 .65126668331495807251485530287027359008800, .5213849287169042769857433808553971486762,
2382 .65330127201274557080523663898929953575150, .5203252032520325203252032520325203252033,
2383 .65533172956312757406749369692988693714150, .5192697768762677484787018255578093306288,
2384 .65735807270835999727154330685152672231200, .5182186234817813765182186234817813765182,
2385 .65938031808912778153342060249997302889800, .5171717171717171717171717171717171717172,
2386 .66139848224536490484126716182800009846700, .5161290322580645161290322580645161290323,
2387 .66341258161706617713093692145776003599150, .5150905432595573440643863179074446680080,
2388 .66542263254509037562201001492212526500250, .5140562248995983935742971887550200803213,
2389 .66742865127195616370414654738851822912700, .5130260521042084168336673346693386773547,
2390 .66943065394262923906154583164607174694550, .5120000000000000000000000000000000000000,
2391 .67142865660530226534774556057527661323550, .5109780439121756487025948103792415169661,
2392 .67342267521216669923234121597488410770900, .5099601593625498007968127490039840637450,
2393 .67541272562017662384192817626171745359900, .5089463220675944333996023856858846918489,
2394 .67739882359180603188519853574689477682100, .5079365079365079365079365079365079365079,
2395 .67938098479579733801614338517538271844400, .5069306930693069306930693069306930693069,
2396 .68135922480790300781450241629499942064300, .5059288537549407114624505928853754940711,
2397 .68333355911162063645036823800182901322850, .5049309664694280078895463510848126232742,
2398 .68530400309891936760919861626462079584600, .5039370078740157480314960629921259842520,
2399 .68727057207096020619019327568821609020250, .5029469548133595284872298624754420432220,
2400 .68923328123880889251040571252815425395950, .5019607843137254901960784313725490196078,
2401 .69314718055994530941723212145818, 5.0e-01,
2402 };
2403
getLogTab64f()2404 const double* getLogTab64f()
2405 {
2406 return logTab;
2407 }
2408
getLogTab32f()2409 const float* getLogTab32f()
2410 {
2411 static float CV_DECL_ALIGNED(64) logTab_f[(LOGTAB_MASK+1)*2];
2412 static std::atomic<bool> logTab_f_initialized(false);
2413 if (!logTab_f_initialized.load())
2414 {
2415 for (int j = 0; j < (LOGTAB_MASK+1)*2; j++)
2416 logTab_f[j] = (float)logTab[j];
2417 logTab_f_initialized = true;
2418 }
2419 return logTab_f;
2420 }
2421
2422 }} // namespace
2423
2424 /* End of file. */
2425