1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
8 //
9 //
10 // Intel License Agreement
11 // For Open Source Computer Vision Library
12 //
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
15 //
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
18 //
19 // * Redistribution's of source code must retain the above copyright notice,
20 // this list of conditions and the following disclaimer.
21 //
22 // * Redistribution's in binary form must reproduce the above copyright notice,
23 // this list of conditions and the following disclaimer in the documentation
24 // and/or other materials provided with the distribution.
25 //
26 // * The name of Intel Corporation may not be used to endorse or promote products
27 // derived from this software without specific prior written permission.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the Intel Corporation or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
39 //
40 //M*/
41
42 #include "precomp.hpp"
43 #include "opencl_kernels_imgproc.hpp"
44 #include "opencv2/core/hal/intrin.hpp"
45
46 namespace cv
47 {
48
49 // The function calculates center of gravity and the central second order moments
completeMomentState(Moments * moments)50 static void completeMomentState( Moments* moments )
51 {
52 double cx = 0, cy = 0;
53 double mu20, mu11, mu02;
54 double inv_m00 = 0.0;
55 assert( moments != 0 );
56
57 if( fabs(moments->m00) > DBL_EPSILON )
58 {
59 inv_m00 = 1. / moments->m00;
60 cx = moments->m10 * inv_m00;
61 cy = moments->m01 * inv_m00;
62 }
63
64 // mu20 = m20 - m10*cx
65 mu20 = moments->m20 - moments->m10 * cx;
66 // mu11 = m11 - m10*cy
67 mu11 = moments->m11 - moments->m10 * cy;
68 // mu02 = m02 - m01*cy
69 mu02 = moments->m02 - moments->m01 * cy;
70
71 moments->mu20 = mu20;
72 moments->mu11 = mu11;
73 moments->mu02 = mu02;
74
75 // mu30 = m30 - cx*(3*mu20 + cx*m10)
76 moments->mu30 = moments->m30 - cx * (3 * mu20 + cx * moments->m10);
77 mu11 += mu11;
78 // mu21 = m21 - cx*(2*mu11 + cx*m01) - cy*mu20
79 moments->mu21 = moments->m21 - cx * (mu11 + cx * moments->m01) - cy * mu20;
80 // mu12 = m12 - cy*(2*mu11 + cy*m10) - cx*mu02
81 moments->mu12 = moments->m12 - cy * (mu11 + cy * moments->m10) - cx * mu02;
82 // mu03 = m03 - cy*(3*mu02 + cy*m01)
83 moments->mu03 = moments->m03 - cy * (3 * mu02 + cy * moments->m01);
84
85
86 double inv_sqrt_m00 = std::sqrt(std::abs(inv_m00));
87 double s2 = inv_m00*inv_m00, s3 = s2*inv_sqrt_m00;
88
89 moments->nu20 = moments->mu20*s2; moments->nu11 = moments->mu11*s2; moments->nu02 = moments->mu02*s2;
90 moments->nu30 = moments->mu30*s3; moments->nu21 = moments->mu21*s3; moments->nu12 = moments->mu12*s3; moments->nu03 = moments->mu03*s3;
91
92 }
93
94
contourMoments(const Mat & contour)95 static Moments contourMoments( const Mat& contour )
96 {
97 Moments m;
98 int lpt = contour.checkVector(2);
99 int is_float = contour.depth() == CV_32F;
100 const Point* ptsi = contour.ptr<Point>();
101 const Point2f* ptsf = contour.ptr<Point2f>();
102
103 CV_Assert( contour.depth() == CV_32S || contour.depth() == CV_32F );
104
105 if( lpt == 0 )
106 return m;
107
108 double a00 = 0, a10 = 0, a01 = 0, a20 = 0, a11 = 0, a02 = 0, a30 = 0, a21 = 0, a12 = 0, a03 = 0;
109 double xi, yi, xi2, yi2, xi_1, yi_1, xi_12, yi_12, dxy, xii_1, yii_1;
110
111 if( !is_float )
112 {
113 xi_1 = ptsi[lpt-1].x;
114 yi_1 = ptsi[lpt-1].y;
115 }
116 else
117 {
118 xi_1 = ptsf[lpt-1].x;
119 yi_1 = ptsf[lpt-1].y;
120 }
121
122 xi_12 = xi_1 * xi_1;
123 yi_12 = yi_1 * yi_1;
124
125 for( int i = 0; i < lpt; i++ )
126 {
127 if( !is_float )
128 {
129 xi = ptsi[i].x;
130 yi = ptsi[i].y;
131 }
132 else
133 {
134 xi = ptsf[i].x;
135 yi = ptsf[i].y;
136 }
137
138 xi2 = xi * xi;
139 yi2 = yi * yi;
140 dxy = xi_1 * yi - xi * yi_1;
141 xii_1 = xi_1 + xi;
142 yii_1 = yi_1 + yi;
143
144 a00 += dxy;
145 a10 += dxy * xii_1;
146 a01 += dxy * yii_1;
147 a20 += dxy * (xi_1 * xii_1 + xi2);
148 a11 += dxy * (xi_1 * (yii_1 + yi_1) + xi * (yii_1 + yi));
149 a02 += dxy * (yi_1 * yii_1 + yi2);
150 a30 += dxy * xii_1 * (xi_12 + xi2);
151 a03 += dxy * yii_1 * (yi_12 + yi2);
152 a21 += dxy * (xi_12 * (3 * yi_1 + yi) + 2 * xi * xi_1 * yii_1 +
153 xi2 * (yi_1 + 3 * yi));
154 a12 += dxy * (yi_12 * (3 * xi_1 + xi) + 2 * yi * yi_1 * xii_1 +
155 yi2 * (xi_1 + 3 * xi));
156 xi_1 = xi;
157 yi_1 = yi;
158 xi_12 = xi2;
159 yi_12 = yi2;
160 }
161
162 if( fabs(a00) > FLT_EPSILON )
163 {
164 double db1_2, db1_6, db1_12, db1_24, db1_20, db1_60;
165
166 if( a00 > 0 )
167 {
168 db1_2 = 0.5;
169 db1_6 = 0.16666666666666666666666666666667;
170 db1_12 = 0.083333333333333333333333333333333;
171 db1_24 = 0.041666666666666666666666666666667;
172 db1_20 = 0.05;
173 db1_60 = 0.016666666666666666666666666666667;
174 }
175 else
176 {
177 db1_2 = -0.5;
178 db1_6 = -0.16666666666666666666666666666667;
179 db1_12 = -0.083333333333333333333333333333333;
180 db1_24 = -0.041666666666666666666666666666667;
181 db1_20 = -0.05;
182 db1_60 = -0.016666666666666666666666666666667;
183 }
184
185 // spatial moments
186 m.m00 = a00 * db1_2;
187 m.m10 = a10 * db1_6;
188 m.m01 = a01 * db1_6;
189 m.m20 = a20 * db1_12;
190 m.m11 = a11 * db1_24;
191 m.m02 = a02 * db1_12;
192 m.m30 = a30 * db1_20;
193 m.m21 = a21 * db1_60;
194 m.m12 = a12 * db1_60;
195 m.m03 = a03 * db1_20;
196
197 completeMomentState( &m );
198 }
199 return m;
200 }
201
202
203 /****************************************************************************************\
204 * Spatial Raster Moments *
205 \****************************************************************************************/
206
207 template<typename T, typename WT, typename MT>
208 struct MomentsInTile_SIMD
209 {
operator ()cv::MomentsInTile_SIMD210 int operator() (const T *, int, WT &, WT &, WT &, MT &)
211 {
212 return 0;
213 }
214 };
215
216 #if CV_SIMD128
217
218 template <>
219 struct MomentsInTile_SIMD<uchar, int, int>
220 {
MomentsInTile_SIMDcv::MomentsInTile_SIMD221 MomentsInTile_SIMD()
222 {
223 // nothing
224 }
225
operator ()cv::MomentsInTile_SIMD226 int operator() (const uchar * ptr, int len, int & x0, int & x1, int & x2, int & x3)
227 {
228 int x = 0;
229
230 {
231 v_int16x8 dx = v_setall_s16(8), qx = v_int16x8(0, 1, 2, 3, 4, 5, 6, 7);
232 v_uint32x4 z = v_setzero_u32(), qx0 = z, qx1 = z, qx2 = z, qx3 = z;
233
234 for( ; x <= len - 8; x += 8 )
235 {
236 v_int16x8 p = v_reinterpret_as_s16(v_load_expand(ptr + x));
237 v_int16x8 sx = v_mul_wrap(qx, qx);
238
239 qx0 += v_reinterpret_as_u32(p);
240 qx1 = v_reinterpret_as_u32(v_dotprod(p, qx, v_reinterpret_as_s32(qx1)));
241 qx2 = v_reinterpret_as_u32(v_dotprod(p, sx, v_reinterpret_as_s32(qx2)));
242 qx3 = v_reinterpret_as_u32(v_dotprod(v_mul_wrap(p, qx), sx, v_reinterpret_as_s32(qx3)));
243
244 qx += dx;
245 }
246
247 x0 = v_reduce_sum(qx0);
248 x0 = (x0 & 0xffff) + (x0 >> 16);
249 x1 = v_reduce_sum(qx1);
250 x2 = v_reduce_sum(qx2);
251 x3 = v_reduce_sum(qx3);
252 }
253
254 return x;
255 }
256 };
257
258 template <>
259 struct MomentsInTile_SIMD<ushort, int, int64>
260 {
MomentsInTile_SIMDcv::MomentsInTile_SIMD261 MomentsInTile_SIMD()
262 {
263 // nothing
264 }
265
operator ()cv::MomentsInTile_SIMD266 int operator() (const ushort * ptr, int len, int & x0, int & x1, int & x2, int64 & x3)
267 {
268 int x = 0;
269
270 {
271 v_int32x4 v_delta = v_setall_s32(4), v_ix0 = v_int32x4(0, 1, 2, 3);
272 v_uint32x4 z = v_setzero_u32(), v_x0 = z, v_x1 = z, v_x2 = z;
273 v_uint64x2 v_x3 = v_reinterpret_as_u64(z);
274
275 for( ; x <= len - 4; x += 4 )
276 {
277 v_int32x4 v_src = v_reinterpret_as_s32(v_load_expand(ptr + x));
278
279 v_x0 += v_reinterpret_as_u32(v_src);
280 v_x1 += v_reinterpret_as_u32(v_src * v_ix0);
281
282 v_int32x4 v_ix1 = v_ix0 * v_ix0;
283 v_x2 += v_reinterpret_as_u32(v_src * v_ix1);
284
285 v_ix1 = v_ix0 * v_ix1;
286 v_src = v_src * v_ix1;
287 v_uint64x2 v_lo, v_hi;
288 v_expand(v_reinterpret_as_u32(v_src), v_lo, v_hi);
289 v_x3 += v_lo + v_hi;
290
291 v_ix0 += v_delta;
292 }
293
294 x0 = v_reduce_sum(v_x0);
295 x1 = v_reduce_sum(v_x1);
296 x2 = v_reduce_sum(v_x2);
297 v_store_aligned(buf64, v_reinterpret_as_s64(v_x3));
298 x3 = buf64[0] + buf64[1];
299 }
300
301 return x;
302 }
303
304 int64 CV_DECL_ALIGNED(16) buf64[2];
305 };
306
307 #endif
308
309 template<typename T, typename WT, typename MT>
310 #if defined __GNUC__ && __GNUC__ == 4 && __GNUC_MINOR__ >= 5 && __GNUC_MINOR__ < 9
311 // Workaround for http://gcc.gnu.org/bugzilla/show_bug.cgi?id=60196
312 __attribute__((optimize("no-tree-vectorize")))
313 #endif
momentsInTile(const Mat & img,double * moments)314 static void momentsInTile( const Mat& img, double* moments )
315 {
316 Size size = img.size();
317 int x, y;
318 MT mom[10] = {0,0,0,0,0,0,0,0,0,0};
319 MomentsInTile_SIMD<T, WT, MT> vop;
320
321 for( y = 0; y < size.height; y++ )
322 {
323 const T* ptr = img.ptr<T>(y);
324 WT x0 = 0, x1 = 0, x2 = 0;
325 MT x3 = 0;
326 x = vop(ptr, size.width, x0, x1, x2, x3);
327
328 for( ; x < size.width; x++ )
329 {
330 WT p = ptr[x];
331 WT xp = x * p, xxp;
332
333 x0 += p;
334 x1 += xp;
335 xxp = xp * x;
336 x2 += xxp;
337 x3 += xxp * x;
338 }
339
340 WT py = y * x0, sy = y*y;
341
342 mom[9] += ((MT)py) * sy; // m03
343 mom[8] += ((MT)x1) * sy; // m12
344 mom[7] += ((MT)x2) * y; // m21
345 mom[6] += x3; // m30
346 mom[5] += x0 * sy; // m02
347 mom[4] += x1 * y; // m11
348 mom[3] += x2; // m20
349 mom[2] += py; // m01
350 mom[1] += x1; // m10
351 mom[0] += x0; // m00
352 }
353
354 for( x = 0; x < 10; x++ )
355 moments[x] = (double)mom[x];
356 }
357
358 typedef void (*MomentsInTileFunc)(const Mat& img, double* moments);
359
Moments()360 Moments::Moments()
361 {
362 m00 = m10 = m01 = m20 = m11 = m02 = m30 = m21 = m12 = m03 =
363 mu20 = mu11 = mu02 = mu30 = mu21 = mu12 = mu03 =
364 nu20 = nu11 = nu02 = nu30 = nu21 = nu12 = nu03 = 0.;
365 }
366
Moments(double _m00,double _m10,double _m01,double _m20,double _m11,double _m02,double _m30,double _m21,double _m12,double _m03)367 Moments::Moments( double _m00, double _m10, double _m01, double _m20, double _m11,
368 double _m02, double _m30, double _m21, double _m12, double _m03 )
369 {
370 m00 = _m00; m10 = _m10; m01 = _m01;
371 m20 = _m20; m11 = _m11; m02 = _m02;
372 m30 = _m30; m21 = _m21; m12 = _m12; m03 = _m03;
373
374 double cx = 0, cy = 0, inv_m00 = 0;
375 if( std::abs(m00) > DBL_EPSILON )
376 {
377 inv_m00 = 1./m00;
378 cx = m10*inv_m00; cy = m01*inv_m00;
379 }
380
381 mu20 = m20 - m10*cx;
382 mu11 = m11 - m10*cy;
383 mu02 = m02 - m01*cy;
384
385 mu30 = m30 - cx*(3*mu20 + cx*m10);
386 mu21 = m21 - cx*(2*mu11 + cx*m01) - cy*mu20;
387 mu12 = m12 - cy*(2*mu11 + cy*m10) - cx*mu02;
388 mu03 = m03 - cy*(3*mu02 + cy*m01);
389
390 double inv_sqrt_m00 = std::sqrt(std::abs(inv_m00));
391 double s2 = inv_m00*inv_m00, s3 = s2*inv_sqrt_m00;
392
393 nu20 = mu20*s2; nu11 = mu11*s2; nu02 = mu02*s2;
394 nu30 = mu30*s3; nu21 = mu21*s3; nu12 = mu12*s3; nu03 = mu03*s3;
395 }
396
397 #ifdef HAVE_OPENCL
398
ocl_moments(InputArray _src,Moments & m,bool binary)399 static bool ocl_moments( InputArray _src, Moments& m, bool binary)
400 {
401 const int TILE_SIZE = 32;
402 const int K = 10;
403
404 Size sz = _src.getSz();
405 int xtiles = divUp(sz.width, TILE_SIZE);
406 int ytiles = divUp(sz.height, TILE_SIZE);
407 int ntiles = xtiles*ytiles;
408 if (ntiles == 0)
409 return false;
410
411 ocl::Kernel k = ocl::Kernel("moments", ocl::imgproc::moments_oclsrc,
412 format("-D TILE_SIZE=%d%s",
413 TILE_SIZE,
414 binary ? " -D OP_MOMENTS_BINARY" : ""));
415
416 if( k.empty() )
417 return false;
418
419 UMat src = _src.getUMat();
420 UMat umbuf(1, ntiles*K, CV_32S);
421
422 size_t globalsize[] = {(size_t)xtiles, std::max((size_t)TILE_SIZE, (size_t)sz.height)};
423 size_t localsize[] = {1, TILE_SIZE};
424 bool ok = k.args(ocl::KernelArg::ReadOnly(src),
425 ocl::KernelArg::PtrWriteOnly(umbuf),
426 xtiles).run(2, globalsize, localsize, true);
427 if(!ok)
428 return false;
429 Mat mbuf = umbuf.getMat(ACCESS_READ);
430 for( int i = 0; i < ntiles; i++ )
431 {
432 double x = (i % xtiles)*TILE_SIZE, y = (i / xtiles)*TILE_SIZE;
433 const int* mom = mbuf.ptr<int>() + i*K;
434 double xm = x * mom[0], ym = y * mom[0];
435
436 // accumulate moments computed in each tile
437
438 // + m00 ( = m00' )
439 m.m00 += mom[0];
440
441 // + m10 ( = m10' + x*m00' )
442 m.m10 += mom[1] + xm;
443
444 // + m01 ( = m01' + y*m00' )
445 m.m01 += mom[2] + ym;
446
447 // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
448 m.m20 += mom[3] + x * (mom[1] * 2 + xm);
449
450 // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
451 m.m11 += mom[4] + x * (mom[2] + ym) + y * mom[1];
452
453 // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
454 m.m02 += mom[5] + y * (mom[2] * 2 + ym);
455
456 // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
457 m.m30 += mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
458
459 // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
460 m.m21 += mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
461
462 // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
463 m.m12 += mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
464
465 // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
466 m.m03 += mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
467 }
468
469 completeMomentState( &m );
470
471 return true;
472 }
473
474 #endif
475
476 #ifdef HAVE_IPP
477 typedef IppStatus (CV_STDCALL * ippiMoments)(const void* pSrc, int srcStep, IppiSize roiSize, IppiMomentState_64f* pCtx);
478
ipp_moments(Mat & src,Moments & m)479 static bool ipp_moments(Mat &src, Moments &m )
480 {
481 #if IPP_VERSION_X100 >= 900
482 CV_INSTRUMENT_REGION_IPP();
483
484 #if IPP_VERSION_X100 < 201801
485 // Degradations for CV_8UC1
486 if(src.type() == CV_8UC1)
487 return false;
488 #endif
489
490 IppiSize roi = { src.cols, src.rows };
491 IppiPoint point = { 0, 0 };
492 int type = src.type();
493 IppStatus ippStatus;
494
495 IppAutoBuffer<IppiMomentState_64f> state;
496 int stateSize = 0;
497
498 ippiMoments ippiMoments64f =
499 (type == CV_8UC1)?(ippiMoments)ippiMoments64f_8u_C1R:
500 (type == CV_16UC1)?(ippiMoments)ippiMoments64f_16u_C1R:
501 (type == CV_32FC1)?(ippiMoments)ippiMoments64f_32f_C1R:
502 NULL;
503 if(!ippiMoments64f)
504 return false;
505
506 ippStatus = ippiMomentGetStateSize_64f(ippAlgHintAccurate, &stateSize);
507 if(ippStatus < 0)
508 return false;
509
510 if(!state.allocate(stateSize) && stateSize)
511 return false;
512
513 ippStatus = ippiMomentInit_64f(state, ippAlgHintAccurate);
514 if(ippStatus < 0)
515 return false;
516
517 ippStatus = CV_INSTRUMENT_FUN_IPP(ippiMoments64f, src.ptr<Ipp8u>(), (int)src.step, roi, state);
518 if(ippStatus < 0)
519 return false;
520
521 ippStatus = ippiGetSpatialMoment_64f(state, 0, 0, 0, point, &m.m00);
522 if(ippStatus < 0)
523 return false;
524 ippiGetSpatialMoment_64f(state, 1, 0, 0, point, &m.m10);
525 ippiGetSpatialMoment_64f(state, 0, 1, 0, point, &m.m01);
526 ippiGetSpatialMoment_64f(state, 2, 0, 0, point, &m.m20);
527 ippiGetSpatialMoment_64f(state, 1, 1, 0, point, &m.m11);
528 ippiGetSpatialMoment_64f(state, 0, 2, 0, point, &m.m02);
529 ippiGetSpatialMoment_64f(state, 3, 0, 0, point, &m.m30);
530 ippiGetSpatialMoment_64f(state, 2, 1, 0, point, &m.m21);
531 ippiGetSpatialMoment_64f(state, 1, 2, 0, point, &m.m12);
532 ippiGetSpatialMoment_64f(state, 0, 3, 0, point, &m.m03);
533
534 ippStatus = ippiGetCentralMoment_64f(state, 2, 0, 0, &m.mu20);
535 if(ippStatus < 0)
536 return false;
537 ippiGetCentralMoment_64f(state, 1, 1, 0, &m.mu11);
538 ippiGetCentralMoment_64f(state, 0, 2, 0, &m.mu02);
539 ippiGetCentralMoment_64f(state, 3, 0, 0, &m.mu30);
540 ippiGetCentralMoment_64f(state, 2, 1, 0, &m.mu21);
541 ippiGetCentralMoment_64f(state, 1, 2, 0, &m.mu12);
542 ippiGetCentralMoment_64f(state, 0, 3, 0, &m.mu03);
543
544 ippStatus = ippiGetNormalizedCentralMoment_64f(state, 2, 0, 0, &m.nu20);
545 if(ippStatus < 0)
546 return false;
547 ippiGetNormalizedCentralMoment_64f(state, 1, 1, 0, &m.nu11);
548 ippiGetNormalizedCentralMoment_64f(state, 0, 2, 0, &m.nu02);
549 ippiGetNormalizedCentralMoment_64f(state, 3, 0, 0, &m.nu30);
550 ippiGetNormalizedCentralMoment_64f(state, 2, 1, 0, &m.nu21);
551 ippiGetNormalizedCentralMoment_64f(state, 1, 2, 0, &m.nu12);
552 ippiGetNormalizedCentralMoment_64f(state, 0, 3, 0, &m.nu03);
553
554 return true;
555 #else
556 CV_UNUSED(src); CV_UNUSED(m);
557 return false;
558 #endif
559 }
560 #endif
561
562 }
563
moments(InputArray _src,bool binary)564 cv::Moments cv::moments( InputArray _src, bool binary )
565 {
566 CV_INSTRUMENT_REGION();
567
568 const int TILE_SIZE = 32;
569 MomentsInTileFunc func = 0;
570 uchar nzbuf[TILE_SIZE*TILE_SIZE];
571 Moments m;
572 int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
573 Size size = _src.size();
574
575 if( size.width <= 0 || size.height <= 0 )
576 return m;
577
578 #ifdef HAVE_OPENCL
579 CV_OCL_RUN_(type == CV_8UC1 && _src.isUMat(), ocl_moments(_src, m, binary), m);
580 #endif
581
582 Mat mat = _src.getMat();
583 if( mat.checkVector(2) >= 0 && (depth == CV_32F || depth == CV_32S))
584 return contourMoments(mat);
585
586 if( cn > 1 )
587 CV_Error( CV_StsBadArg, "Invalid image type (must be single-channel)" );
588
589 CV_IPP_RUN(!binary, ipp_moments(mat, m), m);
590
591 if( binary || depth == CV_8U )
592 func = momentsInTile<uchar, int, int>;
593 else if( depth == CV_16U )
594 func = momentsInTile<ushort, int, int64>;
595 else if( depth == CV_16S )
596 func = momentsInTile<short, int, int64>;
597 else if( depth == CV_32F )
598 func = momentsInTile<float, double, double>;
599 else if( depth == CV_64F )
600 func = momentsInTile<double, double, double>;
601 else
602 CV_Error( CV_StsUnsupportedFormat, "" );
603
604 Mat src0(mat);
605
606 for( int y = 0; y < size.height; y += TILE_SIZE )
607 {
608 Size tileSize;
609 tileSize.height = std::min(TILE_SIZE, size.height - y);
610
611 for( int x = 0; x < size.width; x += TILE_SIZE )
612 {
613 tileSize.width = std::min(TILE_SIZE, size.width - x);
614 Mat src(src0, cv::Rect(x, y, tileSize.width, tileSize.height));
615
616 if( binary )
617 {
618 cv::Mat tmp(tileSize, CV_8U, nzbuf);
619 cv::compare( src, 0, tmp, CV_CMP_NE );
620 src = tmp;
621 }
622
623 double mom[10];
624 func( src, mom );
625
626 if(binary)
627 {
628 double s = 1./255;
629 for( int k = 0; k < 10; k++ )
630 mom[k] *= s;
631 }
632
633 double xm = x * mom[0], ym = y * mom[0];
634
635 // accumulate moments computed in each tile
636
637 // + m00 ( = m00' )
638 m.m00 += mom[0];
639
640 // + m10 ( = m10' + x*m00' )
641 m.m10 += mom[1] + xm;
642
643 // + m01 ( = m01' + y*m00' )
644 m.m01 += mom[2] + ym;
645
646 // + m20 ( = m20' + 2*x*m10' + x*x*m00' )
647 m.m20 += mom[3] + x * (mom[1] * 2 + xm);
648
649 // + m11 ( = m11' + x*m01' + y*m10' + x*y*m00' )
650 m.m11 += mom[4] + x * (mom[2] + ym) + y * mom[1];
651
652 // + m02 ( = m02' + 2*y*m01' + y*y*m00' )
653 m.m02 += mom[5] + y * (mom[2] * 2 + ym);
654
655 // + m30 ( = m30' + 3*x*m20' + 3*x*x*m10' + x*x*x*m00' )
656 m.m30 += mom[6] + x * (3. * mom[3] + x * (3. * mom[1] + xm));
657
658 // + m21 ( = m21' + x*(2*m11' + 2*y*m10' + x*m01' + x*y*m00') + y*m20')
659 m.m21 += mom[7] + x * (2 * (mom[4] + y * mom[1]) + x * (mom[2] + ym)) + y * mom[3];
660
661 // + m12 ( = m12' + y*(2*m11' + 2*x*m01' + y*m10' + x*y*m00') + x*m02')
662 m.m12 += mom[8] + y * (2 * (mom[4] + x * mom[2]) + y * (mom[1] + xm)) + x * mom[5];
663
664 // + m03 ( = m03' + 3*y*m02' + 3*y*y*m01' + y*y*y*m00' )
665 m.m03 += mom[9] + y * (3. * mom[5] + y * (3. * mom[2] + ym));
666 }
667 }
668
669 completeMomentState( &m );
670 return m;
671 }
672
673
HuMoments(const Moments & m,double hu[7])674 void cv::HuMoments( const Moments& m, double hu[7] )
675 {
676 CV_INSTRUMENT_REGION();
677
678 double t0 = m.nu30 + m.nu12;
679 double t1 = m.nu21 + m.nu03;
680
681 double q0 = t0 * t0, q1 = t1 * t1;
682
683 double n4 = 4 * m.nu11;
684 double s = m.nu20 + m.nu02;
685 double d = m.nu20 - m.nu02;
686
687 hu[0] = s;
688 hu[1] = d * d + n4 * m.nu11;
689 hu[3] = q0 + q1;
690 hu[5] = d * (q0 - q1) + n4 * t0 * t1;
691
692 t0 *= q0 - 3 * q1;
693 t1 *= 3 * q0 - q1;
694
695 q0 = m.nu30 - 3 * m.nu12;
696 q1 = 3 * m.nu21 - m.nu03;
697
698 hu[2] = q0 * q0 + q1 * q1;
699 hu[4] = q0 * t0 + q1 * t1;
700 hu[6] = q1 * t0 - q0 * t1;
701 }
702
HuMoments(const Moments & m,OutputArray _hu)703 void cv::HuMoments( const Moments& m, OutputArray _hu )
704 {
705 CV_INSTRUMENT_REGION();
706
707 _hu.create(7, 1, CV_64F);
708 Mat hu = _hu.getMat();
709 CV_Assert( hu.isContinuous() );
710 HuMoments(m, hu.ptr<double>());
711 }
712
713
cvMoments(const CvArr * arr,CvMoments * moments,int binary)714 CV_IMPL void cvMoments( const CvArr* arr, CvMoments* moments, int binary )
715 {
716 const IplImage* img = (const IplImage*)arr;
717 cv::Mat src;
718 if( CV_IS_IMAGE(arr) && img->roi && img->roi->coi > 0 )
719 cv::extractImageCOI(arr, src, img->roi->coi-1);
720 else
721 src = cv::cvarrToMat(arr);
722 cv::Moments m = cv::moments(src, binary != 0);
723 CV_Assert( moments != 0 );
724 *moments = cvMoments(m);
725 }
726
727
cvGetSpatialMoment(CvMoments * moments,int x_order,int y_order)728 CV_IMPL double cvGetSpatialMoment( CvMoments * moments, int x_order, int y_order )
729 {
730 int order = x_order + y_order;
731
732 if( !moments )
733 CV_Error( CV_StsNullPtr, "" );
734 if( (x_order | y_order) < 0 || order > 3 )
735 CV_Error( CV_StsOutOfRange, "" );
736
737 return (&(moments->m00))[order + (order >> 1) + (order > 2) * 2 + y_order];
738 }
739
740
cvGetCentralMoment(CvMoments * moments,int x_order,int y_order)741 CV_IMPL double cvGetCentralMoment( CvMoments * moments, int x_order, int y_order )
742 {
743 int order = x_order + y_order;
744
745 if( !moments )
746 CV_Error( CV_StsNullPtr, "" );
747 if( (x_order | y_order) < 0 || order > 3 )
748 CV_Error( CV_StsOutOfRange, "" );
749
750 return order >= 2 ? (&(moments->m00))[4 + order * 3 + y_order] :
751 order == 0 ? moments->m00 : 0;
752 }
753
754
cvGetNormalizedCentralMoment(CvMoments * moments,int x_order,int y_order)755 CV_IMPL double cvGetNormalizedCentralMoment( CvMoments * moments, int x_order, int y_order )
756 {
757 int order = x_order + y_order;
758
759 double mu = cvGetCentralMoment( moments, x_order, y_order );
760 double m00s = moments->inv_sqrt_m00;
761
762 while( --order >= 0 )
763 mu *= m00s;
764 return mu * m00s * m00s;
765 }
766
767
cvGetHuMoments(CvMoments * mState,CvHuMoments * HuState)768 CV_IMPL void cvGetHuMoments( CvMoments * mState, CvHuMoments * HuState )
769 {
770 if( !mState || !HuState )
771 CV_Error( CV_StsNullPtr, "" );
772
773 double m00s = mState->inv_sqrt_m00, m00 = m00s * m00s, s2 = m00 * m00, s3 = s2 * m00s;
774
775 double nu20 = mState->mu20 * s2,
776 nu11 = mState->mu11 * s2,
777 nu02 = mState->mu02 * s2,
778 nu30 = mState->mu30 * s3,
779 nu21 = mState->mu21 * s3, nu12 = mState->mu12 * s3, nu03 = mState->mu03 * s3;
780
781 double t0 = nu30 + nu12;
782 double t1 = nu21 + nu03;
783
784 double q0 = t0 * t0, q1 = t1 * t1;
785
786 double n4 = 4 * nu11;
787 double s = nu20 + nu02;
788 double d = nu20 - nu02;
789
790 HuState->hu1 = s;
791 HuState->hu2 = d * d + n4 * nu11;
792 HuState->hu4 = q0 + q1;
793 HuState->hu6 = d * (q0 - q1) + n4 * t0 * t1;
794
795 t0 *= q0 - 3 * q1;
796 t1 *= 3 * q0 - q1;
797
798 q0 = nu30 - 3 * nu12;
799 q1 = 3 * nu21 - nu03;
800
801 HuState->hu3 = q0 * q0 + q1 * q1;
802 HuState->hu5 = q0 * t0 + q1 * t1;
803 HuState->hu7 = q1 * t0 - q0 * t1;
804 }
805
806
807 /* End of file. */
808