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 #include "precomp.hpp"
42 #include "opencv2/core/hal/intrin.hpp"
43
44 namespace cv
45 {
46
47 const float EPS = 1.0e-4f;
48
findCircle3pts(Point2f * pts,Point2f & center,float & radius)49 static void findCircle3pts(Point2f *pts, Point2f ¢er, float &radius)
50 {
51 // two edges of the triangle v1, v2
52 Point2f v1 = pts[1] - pts[0];
53 Point2f v2 = pts[2] - pts[0];
54
55 // center is intersection of midperpendicular lines of the two edges v1, v2
56 // a1*x + b1*y = c1 where a1 = v1.x, b1 = v1.y
57 // a2*x + b2*y = c2 where a2 = v2.x, b2 = v2.y
58 Point2f midPoint1 = (pts[0] + pts[1]) / 2.0f;
59 float c1 = midPoint1.x * v1.x + midPoint1.y * v1.y;
60 Point2f midPoint2 = (pts[0] + pts[2]) / 2.0f;
61 float c2 = midPoint2.x * v2.x + midPoint2.y * v2.y;
62 float det = v1.x * v2.y - v1.y * v2.x;
63 if (fabs(det) <= EPS)
64 {
65 // v1 and v2 are colinear, so the longest distance between any 2 points
66 // is the diameter of the minimum enclosing circle.
67 float d1 = normL2Sqr<float>(pts[0] - pts[1]);
68 float d2 = normL2Sqr<float>(pts[0] - pts[2]);
69 float d3 = normL2Sqr<float>(pts[1] - pts[2]);
70 radius = sqrt(std::max(d1, std::max(d2, d3))) * 0.5f + EPS;
71 if (d1 >= d2 && d1 >= d3)
72 {
73 center = (pts[0] + pts[1]) * 0.5f;
74 }
75 else if (d2 >= d1 && d2 >= d3)
76 {
77 center = (pts[0] + pts[2]) * 0.5f;
78 }
79 else
80 {
81 CV_DbgAssert(d3 >= d1 && d3 >= d2);
82 center = (pts[1] + pts[2]) * 0.5f;
83 }
84 return;
85 }
86 float cx = (c1 * v2.y - c2 * v1.y) / det;
87 float cy = (v1.x * c2 - v2.x * c1) / det;
88 center.x = (float)cx;
89 center.y = (float)cy;
90 cx -= pts[0].x;
91 cy -= pts[0].y;
92 radius = (float)(std::sqrt(cx *cx + cy * cy)) + EPS;
93 }
94
95 template<typename PT>
findThirdPoint(const PT * pts,int i,int j,Point2f & center,float & radius)96 static void findThirdPoint(const PT *pts, int i, int j, Point2f ¢er, float &radius)
97 {
98 center.x = (float)(pts[j].x + pts[i].x) / 2.0f;
99 center.y = (float)(pts[j].y + pts[i].y) / 2.0f;
100 float dx = (float)(pts[j].x - pts[i].x);
101 float dy = (float)(pts[j].y - pts[i].y);
102 radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;
103
104 for (int k = 0; k < j; ++k)
105 {
106 dx = center.x - (float)pts[k].x;
107 dy = center.y - (float)pts[k].y;
108 if (norm(Point2f(dx, dy)) < radius)
109 {
110 continue;
111 }
112 else
113 {
114 Point2f ptsf[3];
115 ptsf[0] = (Point2f)pts[i];
116 ptsf[1] = (Point2f)pts[j];
117 ptsf[2] = (Point2f)pts[k];
118 Point2f new_center; float new_radius = 0;
119 findCircle3pts(ptsf, new_center, new_radius);
120 if (new_radius > 0)
121 {
122 radius = new_radius;
123 center = new_center;
124 }
125 }
126 }
127 }
128
129
130 template<typename PT>
findSecondPoint(const PT * pts,int i,Point2f & center,float & radius)131 void findSecondPoint(const PT *pts, int i, Point2f ¢er, float &radius)
132 {
133 center.x = (float)(pts[0].x + pts[i].x) / 2.0f;
134 center.y = (float)(pts[0].y + pts[i].y) / 2.0f;
135 float dx = (float)(pts[0].x - pts[i].x);
136 float dy = (float)(pts[0].y - pts[i].y);
137 radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;
138
139 for (int j = 1; j < i; ++j)
140 {
141 dx = center.x - (float)pts[j].x;
142 dy = center.y - (float)pts[j].y;
143 if (norm(Point2f(dx, dy)) < radius)
144 {
145 continue;
146 }
147 else
148 {
149 Point2f new_center; float new_radius = 0;
150 findThirdPoint(pts, i, j, new_center, new_radius);
151 if (new_radius > 0)
152 {
153 radius = new_radius;
154 center = new_center;
155 }
156 }
157 }
158 }
159
160
161 template<typename PT>
findMinEnclosingCircle(const PT * pts,int count,Point2f & center,float & radius)162 static void findMinEnclosingCircle(const PT *pts, int count, Point2f ¢er, float &radius)
163 {
164 center.x = (float)(pts[0].x + pts[1].x) / 2.0f;
165 center.y = (float)(pts[0].y + pts[1].y) / 2.0f;
166 float dx = (float)(pts[0].x - pts[1].x);
167 float dy = (float)(pts[0].y - pts[1].y);
168 radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;
169
170 for (int i = 2; i < count; ++i)
171 {
172 dx = (float)pts[i].x - center.x;
173 dy = (float)pts[i].y - center.y;
174 float d = (float)norm(Point2f(dx, dy));
175 if (d < radius)
176 {
177 continue;
178 }
179 else
180 {
181 Point2f new_center; float new_radius = 0;
182 findSecondPoint(pts, i, new_center, new_radius);
183 if (new_radius > 0)
184 {
185 radius = new_radius;
186 center = new_center;
187 }
188 }
189 }
190 }
191 } // namespace cv
192
193 // see Welzl, Emo. Smallest enclosing disks (balls and ellipsoids). Springer Berlin Heidelberg, 1991.
minEnclosingCircle(InputArray _points,Point2f & _center,float & _radius)194 void cv::minEnclosingCircle( InputArray _points, Point2f& _center, float& _radius )
195 {
196 CV_INSTRUMENT_REGION();
197
198 Mat points = _points.getMat();
199 int count = points.checkVector(2);
200 int depth = points.depth();
201 CV_Assert(count >= 0 && (depth == CV_32F || depth == CV_32S));
202
203 _center.x = _center.y = 0.f;
204 _radius = 0.f;
205
206 if( count == 0 )
207 return;
208
209 bool is_float = depth == CV_32F;
210 const Point* ptsi = points.ptr<Point>();
211 const Point2f* ptsf = points.ptr<Point2f>();
212
213 switch (count)
214 {
215 case 1:
216 {
217 _center = (is_float) ? ptsf[0] : Point2f((float)ptsi[0].x, (float)ptsi[0].y);
218 _radius = EPS;
219 break;
220 }
221 case 2:
222 {
223 Point2f p1 = (is_float) ? ptsf[0] : Point2f((float)ptsi[0].x, (float)ptsi[0].y);
224 Point2f p2 = (is_float) ? ptsf[1] : Point2f((float)ptsi[1].x, (float)ptsi[1].y);
225 _center.x = (p1.x + p2.x) / 2.0f;
226 _center.y = (p1.y + p2.y) / 2.0f;
227 _radius = (float)(norm(p1 - p2) / 2.0) + EPS;
228 break;
229 }
230 default:
231 {
232 Point2f center;
233 float radius = 0.f;
234 if (is_float)
235 {
236 findMinEnclosingCircle<Point2f>(ptsf, count, center, radius);
237 #if 0
238 for (size_t m = 0; m < count; ++m)
239 {
240 float d = (float)norm(ptsf[m] - center);
241 if (d > radius)
242 {
243 printf("error!\n");
244 }
245 }
246 #endif
247 }
248 else
249 {
250 findMinEnclosingCircle<Point>(ptsi, count, center, radius);
251 #if 0
252 for (size_t m = 0; m < count; ++m)
253 {
254 double dx = ptsi[m].x - center.x;
255 double dy = ptsi[m].y - center.y;
256 double d = std::sqrt(dx * dx + dy * dy);
257 if (d > radius)
258 {
259 printf("error!\n");
260 }
261 }
262 #endif
263 }
264 _center = center;
265 _radius = radius;
266 break;
267 }
268 }
269 }
270
271
272 // calculates length of a curve (e.g. contour perimeter)
arcLength(InputArray _curve,bool is_closed)273 double cv::arcLength( InputArray _curve, bool is_closed )
274 {
275 CV_INSTRUMENT_REGION();
276
277 Mat curve = _curve.getMat();
278 int count = curve.checkVector(2);
279 int depth = curve.depth();
280 CV_Assert( count >= 0 && (depth == CV_32F || depth == CV_32S));
281 double perimeter = 0;
282
283 int i;
284
285 if( count <= 1 )
286 return 0.;
287
288 bool is_float = depth == CV_32F;
289 int last = is_closed ? count-1 : 0;
290 const Point* pti = curve.ptr<Point>();
291 const Point2f* ptf = curve.ptr<Point2f>();
292
293 Point2f prev = is_float ? ptf[last] : Point2f((float)pti[last].x,(float)pti[last].y);
294
295 for( i = 0; i < count; i++ )
296 {
297 Point2f p = is_float ? ptf[i] : Point2f((float)pti[i].x,(float)pti[i].y);
298 float dx = p.x - prev.x, dy = p.y - prev.y;
299 perimeter += std::sqrt(dx*dx + dy*dy);
300
301 prev = p;
302 }
303
304 return perimeter;
305 }
306
307 // area of a whole sequence
contourArea(InputArray _contour,bool oriented)308 double cv::contourArea( InputArray _contour, bool oriented )
309 {
310 CV_INSTRUMENT_REGION();
311
312 Mat contour = _contour.getMat();
313 int npoints = contour.checkVector(2);
314 int depth = contour.depth();
315 CV_Assert(npoints >= 0 && (depth == CV_32F || depth == CV_32S));
316
317 if( npoints == 0 )
318 return 0.;
319
320 double a00 = 0;
321 bool is_float = depth == CV_32F;
322 const Point* ptsi = contour.ptr<Point>();
323 const Point2f* ptsf = contour.ptr<Point2f>();
324 Point2f prev = is_float ? ptsf[npoints-1] : Point2f((float)ptsi[npoints-1].x, (float)ptsi[npoints-1].y);
325
326 for( int i = 0; i < npoints; i++ )
327 {
328 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
329 a00 += (double)prev.x * p.y - (double)prev.y * p.x;
330 prev = p;
331 }
332
333 a00 *= 0.5;
334 if( !oriented )
335 a00 = fabs(a00);
336
337 return a00;
338 }
339
340 namespace cv
341 {
342
getOfs(int i,float eps)343 static inline Point2f getOfs(int i, float eps)
344 {
345 return Point2f(((i & 1)*2 - 1)*eps, ((i & 2) - 1)*eps);
346 }
347
fitEllipseNoDirect(InputArray _points)348 static RotatedRect fitEllipseNoDirect( InputArray _points )
349 {
350 CV_INSTRUMENT_REGION();
351
352 Mat points = _points.getMat();
353 int i, n = points.checkVector(2);
354 int depth = points.depth();
355 CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S));
356
357 RotatedRect box;
358
359 if( n < 5 )
360 CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" );
361
362 // New fitellipse algorithm, contributed by Dr. Daniel Weiss
363 Point2f c(0,0);
364 double gfp[5] = {0}, rp[5] = {0}, t, vd[25]={0}, wd[5]={0};
365 const double min_eps = 1e-8;
366 bool is_float = depth == CV_32F;
367
368 AutoBuffer<double> _Ad(n*12+n);
369 double *Ad = _Ad.data(), *ud = Ad + n*5, *bd = ud + n*5;
370 Point2f* ptsf_copy = (Point2f*)(bd + n);
371
372 // first fit for parameters A - E
373 Mat A( n, 5, CV_64F, Ad );
374 Mat b( n, 1, CV_64F, bd );
375 Mat x( 5, 1, CV_64F, gfp );
376 Mat u( n, 1, CV_64F, ud );
377 Mat vt( 5, 5, CV_64F, vd );
378 Mat w( 5, 1, CV_64F, wd );
379
380 {
381 const Point* ptsi = points.ptr<Point>();
382 const Point2f* ptsf = points.ptr<Point2f>();
383 for( i = 0; i < n; i++ )
384 {
385 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
386 ptsf_copy[i] = p;
387 c += p;
388 }
389 }
390 c.x /= n;
391 c.y /= n;
392
393 double s = 0;
394 for( i = 0; i < n; i++ )
395 {
396 Point2f p = ptsf_copy[i];
397 p -= c;
398 s += fabs(p.x) + fabs(p.y);
399 }
400 double scale = 100./(s > FLT_EPSILON ? s : FLT_EPSILON);
401
402 for( i = 0; i < n; i++ )
403 {
404 Point2f p = ptsf_copy[i];
405 p -= c;
406 double px = p.x*scale;
407 double py = p.y*scale;
408
409 bd[i] = 10000.0; // 1.0?
410 Ad[i*5] = -px * px; // A - C signs inverted as proposed by APP
411 Ad[i*5 + 1] = -py * py;
412 Ad[i*5 + 2] = -px * py;
413 Ad[i*5 + 3] = px;
414 Ad[i*5 + 4] = py;
415 }
416
417 SVDecomp(A, w, u, vt);
418 if(wd[0]*FLT_EPSILON > wd[4]) {
419 float eps = (float)(s/(n*2)*1e-3);
420 for( i = 0; i < n; i++ )
421 {
422 Point2f p = ptsf_copy[i] + getOfs(i, eps);
423 ptsf_copy[i] = p;
424 }
425
426 for( i = 0; i < n; i++ )
427 {
428 Point2f p = ptsf_copy[i];
429 p -= c;
430 double px = p.x*scale;
431 double py = p.y*scale;
432 bd[i] = 10000.0; // 1.0?
433 Ad[i*5] = -px * px; // A - C signs inverted as proposed by APP
434 Ad[i*5 + 1] = -py * py;
435 Ad[i*5 + 2] = -px * py;
436 Ad[i*5 + 3] = px;
437 Ad[i*5 + 4] = py;
438 }
439 SVDecomp(A, w, u, vt);
440 }
441 SVBackSubst(w, u, vt, b, x);
442
443 // now use general-form parameters A - E to find the ellipse center:
444 // differentiate general form wrt x/y to get two equations for cx and cy
445 A = Mat( 2, 2, CV_64F, Ad );
446 b = Mat( 2, 1, CV_64F, bd );
447 x = Mat( 2, 1, CV_64F, rp );
448 Ad[0] = 2 * gfp[0];
449 Ad[1] = Ad[2] = gfp[2];
450 Ad[3] = 2 * gfp[1];
451 bd[0] = gfp[3];
452 bd[1] = gfp[4];
453 solve( A, b, x, DECOMP_SVD );
454
455 // re-fit for parameters A - C with those center coordinates
456 A = Mat( n, 3, CV_64F, Ad );
457 b = Mat( n, 1, CV_64F, bd );
458 x = Mat( 3, 1, CV_64F, gfp );
459 for( i = 0; i < n; i++ )
460 {
461 Point2f p = ptsf_copy[i];
462 p -= c;
463 double px = p.x*scale;
464 double py = p.y*scale;
465 bd[i] = 1.0;
466 Ad[i * 3] = (px - rp[0]) * (px - rp[0]);
467 Ad[i * 3 + 1] = (py - rp[1]) * (py - rp[1]);
468 Ad[i * 3 + 2] = (px - rp[0]) * (py - rp[1]);
469 }
470 solve(A, b, x, DECOMP_SVD);
471
472 // store angle and radii
473 rp[4] = -0.5 * atan2(gfp[2], gfp[1] - gfp[0]); // convert from APP angle usage
474 if( fabs(gfp[2]) > min_eps )
475 t = gfp[2]/sin(-2.0 * rp[4]);
476 else // ellipse is rotated by an integer multiple of pi/2
477 t = gfp[1] - gfp[0];
478 rp[2] = fabs(gfp[0] + gfp[1] - t);
479 if( rp[2] > min_eps )
480 rp[2] = std::sqrt(2.0 / rp[2]);
481 rp[3] = fabs(gfp[0] + gfp[1] + t);
482 if( rp[3] > min_eps )
483 rp[3] = std::sqrt(2.0 / rp[3]);
484
485 box.center.x = (float)(rp[0]/scale) + c.x;
486 box.center.y = (float)(rp[1]/scale) + c.y;
487 box.size.width = (float)(rp[2]*2/scale);
488 box.size.height = (float)(rp[3]*2/scale);
489 if( box.size.width > box.size.height )
490 {
491 float tmp;
492 CV_SWAP( box.size.width, box.size.height, tmp );
493 box.angle = (float)(90 + rp[4]*180/CV_PI);
494 }
495 if( box.angle < -180 )
496 box.angle += 360;
497 if( box.angle > 360 )
498 box.angle -= 360;
499
500 return box;
501 }
502 }
503
fitEllipse(InputArray _points)504 cv::RotatedRect cv::fitEllipse( InputArray _points )
505 {
506 CV_INSTRUMENT_REGION();
507
508 Mat points = _points.getMat();
509 int n = points.checkVector(2);
510 return n == 5 ? fitEllipseDirect(points) : fitEllipseNoDirect(points);
511 }
512
fitEllipseAMS(InputArray _points)513 cv::RotatedRect cv::fitEllipseAMS( InputArray _points )
514 {
515 Mat points = _points.getMat();
516 int i, n = points.checkVector(2);
517 int depth = points.depth();
518 CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S));
519
520 RotatedRect box;
521
522 if( n < 5 )
523 CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" );
524
525 Point2f c(0,0);
526
527 bool is_float = depth == CV_32F;
528 const Point* ptsi = points.ptr<Point>();
529 const Point2f* ptsf = points.ptr<Point2f>();
530
531 Mat A( n, 6, CV_64F);
532 Matx<double, 6, 6> DM;
533 Matx<double, 5, 5> M;
534 Matx<double, 5, 1> pVec;
535 Matx<double, 6, 1> coeffs;
536
537 double x0, y0, a, b, theta;
538
539 for( i = 0; i < n; i++ )
540 {
541 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
542 c += p;
543 }
544 c.x /= (float)n;
545 c.y /= (float)n;
546
547 double s = 0;
548 for( i = 0; i < n; i++ )
549 {
550 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
551 s += fabs(p.x - c.x) + fabs(p.y - c.y);
552 }
553 double scale = 100./(s > FLT_EPSILON ? s : (double)FLT_EPSILON);
554
555 for( i = 0; i < n; i++ )
556 {
557 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
558 double px = (p.x - c.x)*scale, py = (p.y - c.y)*scale;
559
560 A.at<double>(i,0) = px*px;
561 A.at<double>(i,1) = px*py;
562 A.at<double>(i,2) = py*py;
563 A.at<double>(i,3) = px;
564 A.at<double>(i,4) = py;
565 A.at<double>(i,5) = 1.0;
566 }
567 cv::mulTransposed( A, DM, true, noArray(), 1.0, -1 );
568 DM *= (1.0/n);
569 double dnm = ( DM(2,5)*(DM(0,5) + DM(2,5)) - (DM(1,5)*DM(1,5)) );
570 double ddm = (4.*(DM(0,5) + DM(2,5))*( (DM(0,5)*DM(2,5)) - (DM(1,5)*DM(1,5))));
571 double ddmm = (2.*(DM(0,5) + DM(2,5))*( (DM(0,5)*DM(2,5)) - (DM(1,5)*DM(1,5))));
572
573 M(0,0)=((-DM(0,0) + DM(0,2) + DM(0,5)*DM(0,5))*(DM(1,5)*DM(1,5)) + (-2*DM(0,1)*DM(1,5) + DM(0,5)*(DM(0,0) \
574 - (DM(0,5)*DM(0,5)) + (DM(1,5)*DM(1,5))))*DM(2,5) + (DM(0,0) - (DM(0,5)*DM(0,5)))*(DM(2,5)*DM(2,5))) / ddm;
575 M(0,1)=((DM(1,5)*DM(1,5))*(-DM(0,1) + DM(1,2) + DM(0,5)*DM(1,5)) + (DM(0,1)*DM(0,5) - ((DM(0,5)*DM(0,5)) + 2*DM(1,1))*DM(1,5) + \
576 (DM(1,5)*DM(1,5)*DM(1,5)))*DM(2,5) + (DM(0,1) - DM(0,5)*DM(1,5))*(DM(2,5)*DM(2,5))) / ddm;
577 M(0,2)=(-2*DM(1,2)*DM(1,5)*DM(2,5) - DM(0,5)*(DM(2,5)*DM(2,5))*(DM(0,5) + DM(2,5)) + DM(0,2)*dnm + \
578 (DM(1,5)*DM(1,5))*(DM(2,2) + DM(2,5)*(DM(0,5) + DM(2,5))))/ddm;
579 M(0,3)=(DM(1,5)*(DM(1,5)*DM(2,3) - 2*DM(1,3)*DM(2,5)) + DM(0,3)*dnm) / ddm;
580 M(0,4)=(DM(1,5)*(DM(1,5)*DM(2,4) - 2*DM(1,4)*DM(2,5)) + DM(0,4)*dnm) / ddm;
581 M(1,0)=(-(DM(0,2)*DM(0,5)*DM(1,5)) + (2*DM(0,1)*DM(0,5) - DM(0,0)*DM(1,5))*DM(2,5))/ddmm;
582 M(1,1)=(-(DM(0,1)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,2)*DM(1,5)) + 2*DM(1,1)*DM(2,5)))/ddmm;
583 M(1,2)=(-(DM(0,2)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,5)*DM(2,2)) + 2*DM(1,2)*DM(2,5)))/ddmm;
584 M(1,3)=(-(DM(0,3)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,5)*DM(2,3)) + 2*DM(1,3)*DM(2,5)))/ddmm;
585 M(1,4)=(-(DM(0,4)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,5)*DM(2,4)) + 2*DM(1,4)*DM(2,5)))/ddmm;
586 M(2,0)=(-2*DM(0,1)*DM(0,5)*DM(1,5) + (DM(0,0) + (DM(0,5)*DM(0,5)))*(DM(1,5)*DM(1,5)) + DM(0,5)*(-(DM(0,5)*DM(0,5)) \
587 + (DM(1,5)*DM(1,5)))*DM(2,5) - (DM(0,5)*DM(0,5))*(DM(2,5)*DM(2,5)) + DM(0,2)*(-(DM(1,5)*DM(1,5)) + DM(0,5)*(DM(0,5) + DM(2,5)))) / ddm;
588 M(2,1)=((DM(0,5)*DM(0,5))*(DM(1,2) - DM(1,5)*DM(2,5)) + (DM(1,5)*DM(1,5))*(DM(0,1) - DM(1,2) + DM(1,5)*DM(2,5)) \
589 + DM(0,5)*(DM(1,2)*DM(2,5) + DM(1,5)*(-2*DM(1,1) + (DM(1,5)*DM(1,5)) - (DM(2,5)*DM(2,5))))) / ddm;
590 M(2,2)=((DM(0,5)*DM(0,5))*(DM(2,2) - (DM(2,5)*DM(2,5))) + (DM(1,5)*DM(1,5))*(DM(0,2) - DM(2,2) + (DM(2,5)*DM(2,5))) + \
591 DM(0,5)*(-2*DM(1,2)*DM(1,5) + DM(2,5)*((DM(1,5)*DM(1,5)) + DM(2,2) - (DM(2,5)*DM(2,5))))) / ddm;
592 M(2,3)=((DM(1,5)*DM(1,5))*(DM(0,3) - DM(2,3)) + (DM(0,5)*DM(0,5))*DM(2,3) + DM(0,5)*(-2*DM(1,3)*DM(1,5) + DM(2,3)*DM(2,5))) / ddm;
593 M(2,4)=((DM(1,5)*DM(1,5))*(DM(0,4) - DM(2,4)) + (DM(0,5)*DM(0,5))*DM(2,4) + DM(0,5)*(-2*DM(1,4)*DM(1,5) + DM(2,4)*DM(2,5))) / ddm;
594 M(3,0)=DM(0,3);
595 M(3,1)=DM(1,3);
596 M(3,2)=DM(2,3);
597 M(3,3)=DM(3,3);
598 M(3,4)=DM(3,4);
599 M(4,0)=DM(0,4);
600 M(4,1)=DM(1,4);
601 M(4,2)=DM(2,4);
602 M(4,3)=DM(3,4);
603 M(4,4)=DM(4,4);
604
605 if (fabs(cv::determinant(M)) > 1.0e-10) {
606 Mat eVal, eVec;
607 eigenNonSymmetric(M, eVal, eVec);
608
609 // Select the eigen vector {a,b,c,d,e} which has the lowest eigenvalue
610 int minpos = 0;
611 double normi, normEVali, normMinpos, normEValMinpos;
612 normMinpos = sqrt(eVec.at<double>(minpos,0)*eVec.at<double>(minpos,0) + eVec.at<double>(minpos,1)*eVec.at<double>(minpos,1) + \
613 eVec.at<double>(minpos,2)*eVec.at<double>(minpos,2) + eVec.at<double>(minpos,3)*eVec.at<double>(minpos,3) + \
614 eVec.at<double>(minpos,4)*eVec.at<double>(minpos,4) );
615 normEValMinpos = eVal.at<double>(minpos,0) * normMinpos;
616 for (i=1; i<5; i++) {
617 normi = sqrt(eVec.at<double>(i,0)*eVec.at<double>(i,0) + eVec.at<double>(i,1)*eVec.at<double>(i,1) + \
618 eVec.at<double>(i,2)*eVec.at<double>(i,2) + eVec.at<double>(i,3)*eVec.at<double>(i,3) + \
619 eVec.at<double>(i,4)*eVec.at<double>(i,4) );
620 normEVali = eVal.at<double>(i,0) * normi;
621 if (normEVali < normEValMinpos) {
622 minpos = i;
623 normMinpos=normi;
624 normEValMinpos=normEVali;
625 }
626 };
627
628 pVec(0) =eVec.at<double>(minpos,0) / normMinpos;
629 pVec(1) =eVec.at<double>(minpos,1) / normMinpos;
630 pVec(2) =eVec.at<double>(minpos,2) / normMinpos;
631 pVec(3) =eVec.at<double>(minpos,3) / normMinpos;
632 pVec(4) =eVec.at<double>(minpos,4) / normMinpos;
633
634 coeffs(0) =pVec(0) ;
635 coeffs(1) =pVec(1) ;
636 coeffs(2) =pVec(2) ;
637 coeffs(3) =pVec(3) ;
638 coeffs(4) =pVec(4) ;
639 coeffs(5) =-pVec(0) *DM(0,5)-pVec(1) *DM(1,5)-coeffs(2) *DM(2,5);
640
641 // Check that an elliptical solution has been found. AMS sometimes produces Parabolic solutions.
642 bool is_ellipse = (coeffs(0) < 0 && \
643 coeffs(2) < (coeffs(1) *coeffs(1) )/(4.*coeffs(0) ) && \
644 coeffs(5) > (-(coeffs(2) *(coeffs(3) *coeffs(3) )) + coeffs(1) *coeffs(3) *coeffs(4) - coeffs(0) *(coeffs(4) *coeffs(4) )) / \
645 ((coeffs(1) *coeffs(1) ) - 4*coeffs(0) *coeffs(2) )) || \
646 (coeffs(0) > 0 && \
647 coeffs(2) > (coeffs(1) *coeffs(1) )/(4.*coeffs(0) ) && \
648 coeffs(5) < (-(coeffs(2) *(coeffs(3) *coeffs(3) )) + coeffs(1) *coeffs(3) *coeffs(4) - coeffs(0) *(coeffs(4) *coeffs(4) )) / \
649 ( (coeffs(1) *coeffs(1) ) - 4*coeffs(0) *coeffs(2) ));
650 if (is_ellipse) {
651 double u1 = pVec(2) *pVec(3) *pVec(3) - pVec(1) *pVec(3) *pVec(4) + pVec(0) *pVec(4) *pVec(4) + pVec(1) *pVec(1) *coeffs(5) ;
652 double u2 = pVec(0) *pVec(2) *coeffs(5) ;
653 double l1 = sqrt(pVec(1) *pVec(1) + (pVec(0) - pVec(2) )*(pVec(0) - pVec(2) ));
654 double l2 = pVec(0) + pVec(2) ;
655 double l3 = pVec(1) *pVec(1) - 4.0*pVec(0) *pVec(2) ;
656 double p1 = 2.0*pVec(2) *pVec(3) - pVec(1) *pVec(4) ;
657 double p2 = 2.0*pVec(0) *pVec(4) -(pVec(1) *pVec(3) );
658
659 x0 = p1/l3/scale + c.x;
660 y0 = p2/l3/scale + c.y;
661 a = std::sqrt(2.)*sqrt((u1 - 4.0*u2)/((l1 - l2)*l3))/scale;
662 b = std::sqrt(2.)*sqrt(-1.0*((u1 - 4.0*u2)/((l1 + l2)*l3)))/scale;
663 if (pVec(1) == 0) {
664 if (pVec(0) < pVec(2) ) {
665 theta = 0;
666 } else {
667 theta = CV_PI/2.;
668 }
669 } else {
670 theta = CV_PI/2. + 0.5*std::atan2(pVec(1) , (pVec(0) - pVec(2) ));
671 }
672
673 box.center.x = (float)x0;
674 box.center.y = (float)y0;
675 box.size.width = (float)(2.0*a);
676 box.size.height = (float)(2.0*b);
677 if( box.size.width > box.size.height )
678 {
679 float tmp;
680 CV_SWAP( box.size.width, box.size.height, tmp );
681 box.angle = (float)(90 + theta*180/CV_PI);
682 } else {
683 box.angle = (float)(fmod(theta*180/CV_PI,180.0));
684 };
685
686
687 } else {
688 box = cv::fitEllipseDirect( points );
689 }
690 } else {
691 box = cv::fitEllipseNoDirect( points );
692 }
693
694 return box;
695 }
696
fitEllipseDirect(InputArray _points)697 cv::RotatedRect cv::fitEllipseDirect( InputArray _points )
698 {
699 Mat points = _points.getMat();
700 int i, n = points.checkVector(2);
701 int depth = points.depth();
702 float eps = 0;
703 CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S));
704
705 RotatedRect box;
706
707 if( n < 5 )
708 CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" );
709
710 Point2d c(0., 0.);
711
712 bool is_float = (depth == CV_32F);
713 const Point* ptsi = points.ptr<Point>();
714 const Point2f* ptsf = points.ptr<Point2f>();
715
716 Mat A( n, 6, CV_64F);
717 Matx<double, 6, 6> DM;
718 Matx33d M, TM, Q;
719 Matx<double, 3, 1> pVec;
720
721 double x0, y0, a, b, theta, Ts;
722 double s = 0;
723
724 for( i = 0; i < n; i++ )
725 {
726 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
727 c.x += p.x;
728 c.y += p.y;
729 }
730 c.x /= n;
731 c.y /= n;
732
733 for( i = 0; i < n; i++ )
734 {
735 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
736 s += fabs(p.x - c.x) + fabs(p.y - c.y);
737 }
738 double scale = 100./(s > FLT_EPSILON ? s : (double)FLT_EPSILON);
739
740 // first, try the original pointset.
741 // if it's singular, try to shift the points a bit
742 int iter = 0;
743 for( iter = 0; iter < 2; iter++ ) {
744 for( i = 0; i < n; i++ )
745 {
746 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
747 Point2f delta = getOfs(i, eps);
748 double px = (p.x + delta.x - c.x)*scale, py = (p.y + delta.y - c.y)*scale;
749
750 A.at<double>(i,0) = px*px;
751 A.at<double>(i,1) = px*py;
752 A.at<double>(i,2) = py*py;
753 A.at<double>(i,3) = px;
754 A.at<double>(i,4) = py;
755 A.at<double>(i,5) = 1.0;
756 }
757 cv::mulTransposed( A, DM, true, noArray(), 1.0, -1 );
758 DM *= (1.0/n);
759
760 TM(0,0) = DM(0,5)*DM(3,5)*DM(4,4) - DM(0,5)*DM(3,4)*DM(4,5) - DM(0,4)*DM(3,5)*DM(5,4) + \
761 DM(0,3)*DM(4,5)*DM(5,4) + DM(0,4)*DM(3,4)*DM(5,5) - DM(0,3)*DM(4,4)*DM(5,5);
762 TM(0,1) = DM(1,5)*DM(3,5)*DM(4,4) - DM(1,5)*DM(3,4)*DM(4,5) - DM(1,4)*DM(3,5)*DM(5,4) + \
763 DM(1,3)*DM(4,5)*DM(5,4) + DM(1,4)*DM(3,4)*DM(5,5) - DM(1,3)*DM(4,4)*DM(5,5);
764 TM(0,2) = DM(2,5)*DM(3,5)*DM(4,4) - DM(2,5)*DM(3,4)*DM(4,5) - DM(2,4)*DM(3,5)*DM(5,4) + \
765 DM(2,3)*DM(4,5)*DM(5,4) + DM(2,4)*DM(3,4)*DM(5,5) - DM(2,3)*DM(4,4)*DM(5,5);
766 TM(1,0) = DM(0,5)*DM(3,3)*DM(4,5) - DM(0,5)*DM(3,5)*DM(4,3) + DM(0,4)*DM(3,5)*DM(5,3) - \
767 DM(0,3)*DM(4,5)*DM(5,3) - DM(0,4)*DM(3,3)*DM(5,5) + DM(0,3)*DM(4,3)*DM(5,5);
768 TM(1,1) = DM(1,5)*DM(3,3)*DM(4,5) - DM(1,5)*DM(3,5)*DM(4,3) + DM(1,4)*DM(3,5)*DM(5,3) - \
769 DM(1,3)*DM(4,5)*DM(5,3) - DM(1,4)*DM(3,3)*DM(5,5) + DM(1,3)*DM(4,3)*DM(5,5);
770 TM(1,2) = DM(2,5)*DM(3,3)*DM(4,5) - DM(2,5)*DM(3,5)*DM(4,3) + DM(2,4)*DM(3,5)*DM(5,3) - \
771 DM(2,3)*DM(4,5)*DM(5,3) - DM(2,4)*DM(3,3)*DM(5,5) + DM(2,3)*DM(4,3)*DM(5,5);
772 TM(2,0) = DM(0,5)*DM(3,4)*DM(4,3) - DM(0,5)*DM(3,3)*DM(4,4) - DM(0,4)*DM(3,4)*DM(5,3) + \
773 DM(0,3)*DM(4,4)*DM(5,3) + DM(0,4)*DM(3,3)*DM(5,4) - DM(0,3)*DM(4,3)*DM(5,4);
774 TM(2,1) = DM(1,5)*DM(3,4)*DM(4,3) - DM(1,5)*DM(3,3)*DM(4,4) - DM(1,4)*DM(3,4)*DM(5,3) + \
775 DM(1,3)*DM(4,4)*DM(5,3) + DM(1,4)*DM(3,3)*DM(5,4) - DM(1,3)*DM(4,3)*DM(5,4);
776 TM(2,2) = DM(2,5)*DM(3,4)*DM(4,3) - DM(2,5)*DM(3,3)*DM(4,4) - DM(2,4)*DM(3,4)*DM(5,3) + \
777 DM(2,3)*DM(4,4)*DM(5,3) + DM(2,4)*DM(3,3)*DM(5,4) - DM(2,3)*DM(4,3)*DM(5,4);
778
779 Ts=(-(DM(3,5)*DM(4,4)*DM(5,3)) + DM(3,4)*DM(4,5)*DM(5,3) + DM(3,5)*DM(4,3)*DM(5,4) - \
780 DM(3,3)*DM(4,5)*DM(5,4) - DM(3,4)*DM(4,3)*DM(5,5) + DM(3,3)*DM(4,4)*DM(5,5));
781
782 M(0,0) = (DM(2,0) + (DM(2,3)*TM(0,0) + DM(2,4)*TM(1,0) + DM(2,5)*TM(2,0))/Ts)/2.;
783 M(0,1) = (DM(2,1) + (DM(2,3)*TM(0,1) + DM(2,4)*TM(1,1) + DM(2,5)*TM(2,1))/Ts)/2.;
784 M(0,2) = (DM(2,2) + (DM(2,3)*TM(0,2) + DM(2,4)*TM(1,2) + DM(2,5)*TM(2,2))/Ts)/2.;
785 M(1,0) = -DM(1,0) - (DM(1,3)*TM(0,0) + DM(1,4)*TM(1,0) + DM(1,5)*TM(2,0))/Ts;
786 M(1,1) = -DM(1,1) - (DM(1,3)*TM(0,1) + DM(1,4)*TM(1,1) + DM(1,5)*TM(2,1))/Ts;
787 M(1,2) = -DM(1,2) - (DM(1,3)*TM(0,2) + DM(1,4)*TM(1,2) + DM(1,5)*TM(2,2))/Ts;
788 M(2,0) = (DM(0,0) + (DM(0,3)*TM(0,0) + DM(0,4)*TM(1,0) + DM(0,5)*TM(2,0))/Ts)/2.;
789 M(2,1) = (DM(0,1) + (DM(0,3)*TM(0,1) + DM(0,4)*TM(1,1) + DM(0,5)*TM(2,1))/Ts)/2.;
790 M(2,2) = (DM(0,2) + (DM(0,3)*TM(0,2) + DM(0,4)*TM(1,2) + DM(0,5)*TM(2,2))/Ts)/2.;
791
792 double det = fabs(cv::determinant(M));
793 if (fabs(det) > 1.0e-10)
794 break;
795 eps = (float)(s/(n*2)*1e-2);
796 }
797
798 if( iter < 2 ) {
799 Mat eVal, eVec;
800 eigenNonSymmetric(M, eVal, eVec);
801
802 // Select the eigen vector {a,b,c} which satisfies 4ac-b^2 > 0
803 double cond[3];
804 cond[0]=(4.0 * eVec.at<double>(0,0) * eVec.at<double>(0,2) - eVec.at<double>(0,1) * eVec.at<double>(0,1));
805 cond[1]=(4.0 * eVec.at<double>(1,0) * eVec.at<double>(1,2) - eVec.at<double>(1,1) * eVec.at<double>(1,1));
806 cond[2]=(4.0 * eVec.at<double>(2,0) * eVec.at<double>(2,2) - eVec.at<double>(2,1) * eVec.at<double>(2,1));
807 if (cond[0]<cond[1]) {
808 i = (cond[1]<cond[2]) ? 2 : 1;
809 } else {
810 i = (cond[0]<cond[2]) ? 2 : 0;
811 }
812 double norm = std::sqrt(eVec.at<double>(i,0)*eVec.at<double>(i,0) + eVec.at<double>(i,1)*eVec.at<double>(i,1) + eVec.at<double>(i,2)*eVec.at<double>(i,2));
813 if (((eVec.at<double>(i,0)<0.0 ? -1 : 1) * (eVec.at<double>(i,1)<0.0 ? -1 : 1) * (eVec.at<double>(i,2)<0.0 ? -1 : 1)) <= 0.0) {
814 norm=-1.0*norm;
815 }
816 pVec(0) =eVec.at<double>(i,0)/norm; pVec(1) =eVec.at<double>(i,1)/norm;pVec(2) =eVec.at<double>(i,2)/norm;
817
818 // Q = (TM . pVec)/Ts;
819 Q(0,0) = (TM(0,0)*pVec(0) +TM(0,1)*pVec(1) +TM(0,2)*pVec(2) )/Ts;
820 Q(0,1) = (TM(1,0)*pVec(0) +TM(1,1)*pVec(1) +TM(1,2)*pVec(2) )/Ts;
821 Q(0,2) = (TM(2,0)*pVec(0) +TM(2,1)*pVec(1) +TM(2,2)*pVec(2) )/Ts;
822
823 // We compute the ellipse properties in the shifted coordinates as doing so improves the numerical accuracy.
824
825 double u1 = pVec(2)*Q(0,0)*Q(0,0) - pVec(1)*Q(0,0)*Q(0,1) + pVec(0)*Q(0,1)*Q(0,1) + pVec(1)*pVec(1)*Q(0,2);
826 double u2 = pVec(0)*pVec(2)*Q(0,2);
827 double l1 = sqrt(pVec(1)*pVec(1) + (pVec(0) - pVec(2))*(pVec(0) - pVec(2)));
828 double l2 = pVec(0) + pVec(2) ;
829 double l3 = pVec(1)*pVec(1) - 4*pVec(0)*pVec(2) ;
830 double p1 = 2*pVec(2)*Q(0,0) - pVec(1)*Q(0,1);
831 double p2 = 2*pVec(0)*Q(0,1) - pVec(1)*Q(0,0);
832
833 x0 = (p1/l3/scale) + c.x;
834 y0 = (p2/l3/scale) + c.y;
835 a = sqrt(2.)*sqrt((u1 - 4.0*u2)/((l1 - l2)*l3))/scale;
836 b = sqrt(2.)*sqrt(-1.0*((u1 - 4.0*u2)/((l1 + l2)*l3)))/scale;
837 if (pVec(1) == 0) {
838 if (pVec(0) < pVec(2) ) {
839 theta = 0;
840 } else {
841 theta = CV_PI/2.;
842 }
843 } else {
844 theta = CV_PI/2. + 0.5*std::atan2(pVec(1) , (pVec(0) - pVec(2) ));
845 }
846
847 box.center.x = (float)x0;
848 box.center.y = (float)y0;
849 box.size.width = (float)(2.0*a);
850 box.size.height = (float)(2.0*b);
851 if( box.size.width > box.size.height )
852 {
853 float tmp;
854 CV_SWAP( box.size.width, box.size.height, tmp );
855 box.angle = (float)(fmod((90 + theta*180/CV_PI),180.0)) ;
856 } else {
857 box.angle = (float)(fmod(theta*180/CV_PI,180.0));
858 };
859 } else {
860 box = cv::fitEllipseNoDirect( points );
861 }
862 return box;
863 }
864
865
866 namespace cv
867 {
868
869 // Calculates bounding rectangle of a point set or retrieves already calculated
pointSetBoundingRect(const Mat & points)870 static Rect pointSetBoundingRect( const Mat& points )
871 {
872 int npoints = points.checkVector(2);
873 int depth = points.depth();
874 CV_Assert(npoints >= 0 && (depth == CV_32F || depth == CV_32S));
875
876 int xmin = 0, ymin = 0, xmax = -1, ymax = -1, i;
877 bool is_float = depth == CV_32F;
878
879 if( npoints == 0 )
880 return Rect();
881
882 #if CV_SIMD
883 const int64_t* pts = points.ptr<int64_t>();
884
885 if( !is_float )
886 {
887 v_int32 minval, maxval;
888 minval = maxval = v_reinterpret_as_s32(vx_setall_s64(*pts)); //min[0]=pt.x, min[1]=pt.y, min[2]=pt.x, min[3]=pt.y
889 for( i = 1; i <= npoints - v_int32::nlanes/2; i+= v_int32::nlanes/2 )
890 {
891 v_int32 ptXY2 = v_reinterpret_as_s32(vx_load(pts + i));
892 minval = v_min(ptXY2, minval);
893 maxval = v_max(ptXY2, maxval);
894 }
895 minval = v_min(v_reinterpret_as_s32(v_expand_low(v_reinterpret_as_u32(minval))), v_reinterpret_as_s32(v_expand_high(v_reinterpret_as_u32(minval))));
896 maxval = v_max(v_reinterpret_as_s32(v_expand_low(v_reinterpret_as_u32(maxval))), v_reinterpret_as_s32(v_expand_high(v_reinterpret_as_u32(maxval))));
897 if( i <= npoints - v_int32::nlanes/4 )
898 {
899 v_int32 ptXY = v_reinterpret_as_s32(v_expand_low(v_reinterpret_as_u32(vx_load_low(pts + i))));
900 minval = v_min(ptXY, minval);
901 maxval = v_max(ptXY, maxval);
902 i += v_int64::nlanes/2;
903 }
904 for(int j = 16; j < CV_SIMD_WIDTH; j*=2)
905 {
906 minval = v_min(v_reinterpret_as_s32(v_expand_low(v_reinterpret_as_u32(minval))), v_reinterpret_as_s32(v_expand_high(v_reinterpret_as_u32(minval))));
907 maxval = v_max(v_reinterpret_as_s32(v_expand_low(v_reinterpret_as_u32(maxval))), v_reinterpret_as_s32(v_expand_high(v_reinterpret_as_u32(maxval))));
908 }
909 xmin = minval.get0();
910 xmax = maxval.get0();
911 ymin = v_reinterpret_as_s32(v_expand_high(v_reinterpret_as_u32(minval))).get0();
912 ymax = v_reinterpret_as_s32(v_expand_high(v_reinterpret_as_u32(maxval))).get0();
913 #if CV_SIMD_WIDTH > 16
914 if( i < npoints )
915 {
916 v_int32x4 minval2, maxval2;
917 minval2 = maxval2 = v_reinterpret_as_s32(v_expand_low(v_reinterpret_as_u32(v_load_low(pts + i))));
918 for( i++; i < npoints; i++ )
919 {
920 v_int32x4 ptXY = v_reinterpret_as_s32(v_expand_low(v_reinterpret_as_u32(v_load_low(pts + i))));
921 minval2 = v_min(ptXY, minval2);
922 maxval2 = v_max(ptXY, maxval2);
923 }
924 xmin = min(xmin, minval2.get0());
925 xmax = max(xmax, maxval2.get0());
926 ymin = min(ymin, v_reinterpret_as_s32(v_expand_high(v_reinterpret_as_u32(minval2))).get0());
927 ymax = max(ymax, v_reinterpret_as_s32(v_expand_high(v_reinterpret_as_u32(maxval2))).get0());
928 }
929 #endif
930 }
931 else
932 {
933 v_float32 minval, maxval;
934 minval = maxval = v_reinterpret_as_f32(vx_setall_s64(*pts)); //min[0]=pt.x, min[1]=pt.y, min[2]=pt.x, min[3]=pt.y
935 for( i = 1; i <= npoints - v_float32::nlanes/2; i+= v_float32::nlanes/2 )
936 {
937 v_float32 ptXY2 = v_reinterpret_as_f32(vx_load(pts + i));
938 minval = v_min(ptXY2, minval);
939 maxval = v_max(ptXY2, maxval);
940 }
941 minval = v_min(v_reinterpret_as_f32(v_expand_low(v_reinterpret_as_u32(minval))), v_reinterpret_as_f32(v_expand_high(v_reinterpret_as_u32(minval))));
942 maxval = v_max(v_reinterpret_as_f32(v_expand_low(v_reinterpret_as_u32(maxval))), v_reinterpret_as_f32(v_expand_high(v_reinterpret_as_u32(maxval))));
943 if( i <= npoints - v_float32::nlanes/4 )
944 {
945 v_float32 ptXY = v_reinterpret_as_f32(v_expand_low(v_reinterpret_as_u32(vx_load_low(pts + i))));
946 minval = v_min(ptXY, minval);
947 maxval = v_max(ptXY, maxval);
948 i += v_float32::nlanes/4;
949 }
950 for(int j = 16; j < CV_SIMD_WIDTH; j*=2)
951 {
952 minval = v_min(v_reinterpret_as_f32(v_expand_low(v_reinterpret_as_u32(minval))), v_reinterpret_as_f32(v_expand_high(v_reinterpret_as_u32(minval))));
953 maxval = v_max(v_reinterpret_as_f32(v_expand_low(v_reinterpret_as_u32(maxval))), v_reinterpret_as_f32(v_expand_high(v_reinterpret_as_u32(maxval))));
954 }
955 xmin = cvFloor(minval.get0());
956 xmax = cvFloor(maxval.get0());
957 ymin = cvFloor(v_reinterpret_as_f32(v_expand_high(v_reinterpret_as_u32(minval))).get0());
958 ymax = cvFloor(v_reinterpret_as_f32(v_expand_high(v_reinterpret_as_u32(maxval))).get0());
959 #if CV_SIMD_WIDTH > 16
960 if( i < npoints )
961 {
962 v_float32x4 minval2, maxval2;
963 minval2 = maxval2 = v_reinterpret_as_f32(v_expand_low(v_reinterpret_as_u32(v_load_low(pts + i))));
964 for( i++; i < npoints; i++ )
965 {
966 v_float32x4 ptXY = v_reinterpret_as_f32(v_expand_low(v_reinterpret_as_u32(v_load_low(pts + i))));
967 minval2 = v_min(ptXY, minval2);
968 maxval2 = v_max(ptXY, maxval2);
969 }
970 xmin = min(xmin, cvFloor(minval2.get0()));
971 xmax = max(xmax, cvFloor(maxval2.get0()));
972 ymin = min(ymin, cvFloor(v_reinterpret_as_f32(v_expand_high(v_reinterpret_as_u32(minval2))).get0()));
973 ymax = max(ymax, cvFloor(v_reinterpret_as_f32(v_expand_high(v_reinterpret_as_u32(maxval2))).get0()));
974 }
975 #endif
976 }
977 #else
978 const Point* pts = points.ptr<Point>();
979 Point pt = pts[0];
980
981 if( !is_float )
982 {
983 xmin = xmax = pt.x;
984 ymin = ymax = pt.y;
985
986 for( i = 1; i < npoints; i++ )
987 {
988 pt = pts[i];
989
990 if( xmin > pt.x )
991 xmin = pt.x;
992
993 if( xmax < pt.x )
994 xmax = pt.x;
995
996 if( ymin > pt.y )
997 ymin = pt.y;
998
999 if( ymax < pt.y )
1000 ymax = pt.y;
1001 }
1002 }
1003 else
1004 {
1005 Cv32suf v;
1006 // init values
1007 xmin = xmax = CV_TOGGLE_FLT(pt.x);
1008 ymin = ymax = CV_TOGGLE_FLT(pt.y);
1009
1010 for( i = 1; i < npoints; i++ )
1011 {
1012 pt = pts[i];
1013 pt.x = CV_TOGGLE_FLT(pt.x);
1014 pt.y = CV_TOGGLE_FLT(pt.y);
1015
1016 if( xmin > pt.x )
1017 xmin = pt.x;
1018
1019 if( xmax < pt.x )
1020 xmax = pt.x;
1021
1022 if( ymin > pt.y )
1023 ymin = pt.y;
1024
1025 if( ymax < pt.y )
1026 ymax = pt.y;
1027 }
1028
1029 v.i = CV_TOGGLE_FLT(xmin); xmin = cvFloor(v.f);
1030 v.i = CV_TOGGLE_FLT(ymin); ymin = cvFloor(v.f);
1031 // because right and bottom sides of the bounding rectangle are not inclusive
1032 // (note +1 in width and height calculation below), cvFloor is used here instead of cvCeil
1033 v.i = CV_TOGGLE_FLT(xmax); xmax = cvFloor(v.f);
1034 v.i = CV_TOGGLE_FLT(ymax); ymax = cvFloor(v.f);
1035 }
1036 #endif
1037
1038 return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1);
1039 }
1040
1041
maskBoundingRect(const Mat & img)1042 static Rect maskBoundingRect( const Mat& img )
1043 {
1044 CV_Assert( img.depth() <= CV_8S && img.channels() == 1 );
1045
1046 Size size = img.size();
1047 int xmin = size.width, ymin = -1, xmax = -1, ymax = -1, i, j, k;
1048
1049 for( i = 0; i < size.height; i++ )
1050 {
1051 const uchar* _ptr = img.ptr(i);
1052 const uchar* ptr = (const uchar*)alignPtr(_ptr, 4);
1053 int have_nz = 0, k_min, offset = (int)(ptr - _ptr);
1054 j = 0;
1055 offset = MIN(offset, size.width);
1056 for( ; j < offset; j++ )
1057 if( _ptr[j] )
1058 {
1059 have_nz = 1;
1060 break;
1061 }
1062 if( j < offset )
1063 {
1064 if( j < xmin )
1065 xmin = j;
1066 if( j > xmax )
1067 xmax = j;
1068 }
1069 if( offset < size.width )
1070 {
1071 xmin -= offset;
1072 xmax -= offset;
1073 size.width -= offset;
1074 j = 0;
1075 for( ; j <= xmin - 4; j += 4 )
1076 if( *((int*)(ptr+j)) )
1077 break;
1078 for( ; j < xmin; j++ )
1079 if( ptr[j] )
1080 {
1081 xmin = j;
1082 if( j > xmax )
1083 xmax = j;
1084 have_nz = 1;
1085 break;
1086 }
1087 k_min = MAX(j-1, xmax);
1088 k = size.width - 1;
1089 for( ; k > k_min && (k&3) != 3; k-- )
1090 if( ptr[k] )
1091 break;
1092 if( k > k_min && (k&3) == 3 )
1093 {
1094 for( ; k > k_min+3; k -= 4 )
1095 if( *((int*)(ptr+k-3)) )
1096 break;
1097 }
1098 for( ; k > k_min; k-- )
1099 if( ptr[k] )
1100 {
1101 xmax = k;
1102 have_nz = 1;
1103 break;
1104 }
1105 if( !have_nz )
1106 {
1107 j &= ~3;
1108 for( ; j <= k - 3; j += 4 )
1109 if( *((int*)(ptr+j)) )
1110 break;
1111 for( ; j <= k; j++ )
1112 if( ptr[j] )
1113 {
1114 have_nz = 1;
1115 break;
1116 }
1117 }
1118 xmin += offset;
1119 xmax += offset;
1120 size.width += offset;
1121 }
1122 if( have_nz )
1123 {
1124 if( ymin < 0 )
1125 ymin = i;
1126 ymax = i;
1127 }
1128 }
1129
1130 if( xmin >= size.width )
1131 xmin = ymin = 0;
1132 return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1);
1133 }
1134
1135 }
1136
boundingRect(InputArray array)1137 cv::Rect cv::boundingRect(InputArray array)
1138 {
1139 CV_INSTRUMENT_REGION();
1140
1141 Mat m = array.getMat();
1142 return m.depth() <= CV_8U ? maskBoundingRect(m) : pointSetBoundingRect(m);
1143 }
1144
1145 ////////////////////////////////////////////// C API ///////////////////////////////////////////
1146
1147 CV_IMPL int
cvMinEnclosingCircle(const void * array,CvPoint2D32f * _center,float * _radius)1148 cvMinEnclosingCircle( const void* array, CvPoint2D32f * _center, float *_radius )
1149 {
1150 cv::AutoBuffer<double> abuf;
1151 cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf);
1152 cv::Point2f center;
1153 float radius;
1154
1155 cv::minEnclosingCircle(points, center, radius);
1156 if(_center)
1157 *_center = cvPoint2D32f(center);
1158 if(_radius)
1159 *_radius = radius;
1160 return 1;
1161 }
1162
1163 static void
icvMemCopy(double ** buf1,double ** buf2,double ** buf3,int * b_max)1164 icvMemCopy( double **buf1, double **buf2, double **buf3, int *b_max )
1165 {
1166 CV_Assert( (*buf1 != NULL || *buf2 != NULL) && *buf3 != NULL );
1167
1168 int bb = *b_max;
1169 if( *buf2 == NULL )
1170 {
1171 *b_max = 2 * (*b_max);
1172 *buf2 = (double *)cvAlloc( (*b_max) * sizeof( double ));
1173
1174 memcpy( *buf2, *buf3, bb * sizeof( double ));
1175
1176 *buf3 = *buf2;
1177 cvFree( buf1 );
1178 *buf1 = NULL;
1179 }
1180 else
1181 {
1182 *b_max = 2 * (*b_max);
1183 *buf1 = (double *) cvAlloc( (*b_max) * sizeof( double ));
1184
1185 memcpy( *buf1, *buf3, bb * sizeof( double ));
1186
1187 *buf3 = *buf1;
1188 cvFree( buf2 );
1189 *buf2 = NULL;
1190 }
1191 }
1192
1193
1194 /* area of a contour sector */
icvContourSecArea(CvSeq * contour,CvSlice slice)1195 static double icvContourSecArea( CvSeq * contour, CvSlice slice )
1196 {
1197 cv::Point pt; /* pointer to points */
1198 cv::Point pt_s, pt_e; /* first and last points */
1199 CvSeqReader reader; /* points reader of contour */
1200
1201 int p_max = 2, p_ind;
1202 int lpt, flag, i;
1203 double a00; /* unnormalized moments m00 */
1204 double xi, yi, xi_1, yi_1, x0, y0, dxy, sk, sk1, t;
1205 double x_s, y_s, nx, ny, dx, dy, du, dv;
1206 double eps = 1.e-5;
1207 double *p_are1, *p_are2, *p_are;
1208 double area = 0;
1209
1210 CV_Assert( contour != NULL && CV_IS_SEQ_POINT_SET( contour ));
1211
1212 lpt = cvSliceLength( slice, contour );
1213 /*if( n2 >= n1 )
1214 lpt = n2 - n1 + 1;
1215 else
1216 lpt = contour->total - n1 + n2 + 1;*/
1217
1218 if( contour->total <= 0 || lpt <= 2 )
1219 return 0.;
1220
1221 a00 = x0 = y0 = xi_1 = yi_1 = 0;
1222 sk1 = 0;
1223 flag = 0;
1224 dxy = 0;
1225 p_are1 = (double *) cvAlloc( p_max * sizeof( double ));
1226
1227 p_are = p_are1;
1228 p_are2 = NULL;
1229
1230 cvStartReadSeq( contour, &reader, 0 );
1231 cvSetSeqReaderPos( &reader, slice.start_index );
1232 { CvPoint pt_s_ = CV_STRUCT_INITIALIZER; CV_READ_SEQ_ELEM(pt_s_, reader); pt_s = pt_s_; }
1233 p_ind = 0;
1234 cvSetSeqReaderPos( &reader, slice.end_index );
1235 { CvPoint pt_e_ = CV_STRUCT_INITIALIZER; CV_READ_SEQ_ELEM(pt_e_, reader); pt_e = pt_e_; }
1236
1237 /* normal coefficients */
1238 nx = pt_s.y - pt_e.y;
1239 ny = pt_e.x - pt_s.x;
1240 cvSetSeqReaderPos( &reader, slice.start_index );
1241
1242 while( lpt-- > 0 )
1243 {
1244 { CvPoint pt_ = CV_STRUCT_INITIALIZER; CV_READ_SEQ_ELEM(pt_, reader); pt = pt_; }
1245
1246 if( flag == 0 )
1247 {
1248 xi_1 = (double) pt.x;
1249 yi_1 = (double) pt.y;
1250 x0 = xi_1;
1251 y0 = yi_1;
1252 sk1 = 0;
1253 flag = 1;
1254 }
1255 else
1256 {
1257 xi = (double) pt.x;
1258 yi = (double) pt.y;
1259
1260 /**************** edges intersection examination **************************/
1261 sk = nx * (xi - pt_s.x) + ny * (yi - pt_s.y);
1262 if( (fabs( sk ) < eps && lpt > 0) || sk * sk1 < -eps )
1263 {
1264 if( fabs( sk ) < eps )
1265 {
1266 dxy = xi_1 * yi - xi * yi_1;
1267 a00 = a00 + dxy;
1268 dxy = xi * y0 - x0 * yi;
1269 a00 = a00 + dxy;
1270
1271 if( p_ind >= p_max )
1272 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
1273
1274 p_are[p_ind] = a00 / 2.;
1275 p_ind++;
1276 a00 = 0;
1277 sk1 = 0;
1278 x0 = xi;
1279 y0 = yi;
1280 dxy = 0;
1281 }
1282 else
1283 {
1284 /* define intersection point */
1285 dv = yi - yi_1;
1286 du = xi - xi_1;
1287 dx = ny;
1288 dy = -nx;
1289 if( fabs( du ) > eps )
1290 t = ((yi_1 - pt_s.y) * du + dv * (pt_s.x - xi_1)) /
1291 (du * dy - dx * dv);
1292 else
1293 t = (xi_1 - pt_s.x) / dx;
1294 if( t > eps && t < 1 - eps )
1295 {
1296 x_s = pt_s.x + t * dx;
1297 y_s = pt_s.y + t * dy;
1298 dxy = xi_1 * y_s - x_s * yi_1;
1299 a00 += dxy;
1300 dxy = x_s * y0 - x0 * y_s;
1301 a00 += dxy;
1302 if( p_ind >= p_max )
1303 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
1304
1305 p_are[p_ind] = a00 / 2.;
1306 p_ind++;
1307
1308 a00 = 0;
1309 sk1 = 0;
1310 x0 = x_s;
1311 y0 = y_s;
1312 dxy = x_s * yi - xi * y_s;
1313 }
1314 }
1315 }
1316 else
1317 dxy = xi_1 * yi - xi * yi_1;
1318
1319 a00 += dxy;
1320 xi_1 = xi;
1321 yi_1 = yi;
1322 sk1 = sk;
1323
1324 }
1325 }
1326
1327 xi = x0;
1328 yi = y0;
1329 dxy = xi_1 * yi - xi * yi_1;
1330
1331 a00 += dxy;
1332
1333 if( p_ind >= p_max )
1334 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
1335
1336 p_are[p_ind] = a00 / 2.;
1337 p_ind++;
1338
1339 // common area calculation
1340 area = 0;
1341 for( i = 0; i < p_ind; i++ )
1342 area += fabs( p_are[i] );
1343
1344 if( p_are1 != NULL )
1345 cvFree( &p_are1 );
1346 else if( p_are2 != NULL )
1347 cvFree( &p_are2 );
1348
1349 return area;
1350 }
1351
1352
1353 /* external contour area function */
1354 CV_IMPL double
cvContourArea(const void * array,CvSlice slice,int oriented)1355 cvContourArea( const void *array, CvSlice slice, int oriented )
1356 {
1357 double area = 0;
1358
1359 CvContour contour_header;
1360 CvSeq* contour = 0;
1361 CvSeqBlock block;
1362
1363 if( CV_IS_SEQ( array ))
1364 {
1365 contour = (CvSeq*)array;
1366 if( !CV_IS_SEQ_POLYLINE( contour ))
1367 CV_Error( CV_StsBadArg, "Unsupported sequence type" );
1368 }
1369 else
1370 {
1371 contour = cvPointSeqFromMat( CV_SEQ_KIND_CURVE, array, &contour_header, &block );
1372 }
1373
1374 if( cvSliceLength( slice, contour ) == contour->total )
1375 {
1376 cv::AutoBuffer<double> abuf;
1377 cv::Mat points = cv::cvarrToMat(contour, false, false, 0, &abuf);
1378 return cv::contourArea( points, oriented !=0 );
1379 }
1380
1381 if( CV_SEQ_ELTYPE( contour ) != CV_32SC2 )
1382 CV_Error( CV_StsUnsupportedFormat,
1383 "Only curves with integer coordinates are supported in case of contour slice" );
1384 area = icvContourSecArea( contour, slice );
1385 return oriented ? area : fabs(area);
1386 }
1387
1388
1389 /* calculates length of a curve (e.g. contour perimeter) */
1390 CV_IMPL double
cvArcLength(const void * array,CvSlice slice,int is_closed)1391 cvArcLength( const void *array, CvSlice slice, int is_closed )
1392 {
1393 double perimeter = 0;
1394
1395 int i, j = 0, count;
1396 const int N = 16;
1397 float buf[N];
1398 CvMat buffer = cvMat( 1, N, CV_32F, buf );
1399 CvSeqReader reader;
1400 CvContour contour_header;
1401 CvSeq* contour = 0;
1402 CvSeqBlock block;
1403
1404 if( CV_IS_SEQ( array ))
1405 {
1406 contour = (CvSeq*)array;
1407 if( !CV_IS_SEQ_POLYLINE( contour ))
1408 CV_Error( CV_StsBadArg, "Unsupported sequence type" );
1409 if( is_closed < 0 )
1410 is_closed = CV_IS_SEQ_CLOSED( contour );
1411 }
1412 else
1413 {
1414 is_closed = is_closed > 0;
1415 contour = cvPointSeqFromMat(
1416 CV_SEQ_KIND_CURVE | (is_closed ? CV_SEQ_FLAG_CLOSED : 0),
1417 array, &contour_header, &block );
1418 }
1419
1420 if( contour->total > 1 )
1421 {
1422 int is_float = CV_SEQ_ELTYPE( contour ) == CV_32FC2;
1423
1424 cvStartReadSeq( contour, &reader, 0 );
1425 cvSetSeqReaderPos( &reader, slice.start_index );
1426 count = cvSliceLength( slice, contour );
1427
1428 count -= !is_closed && count == contour->total;
1429
1430 // scroll the reader by 1 point
1431 reader.prev_elem = reader.ptr;
1432 CV_NEXT_SEQ_ELEM( sizeof(CvPoint), reader );
1433
1434 for( i = 0; i < count; i++ )
1435 {
1436 float dx, dy;
1437
1438 if( !is_float )
1439 {
1440 CvPoint* pt = (CvPoint*)reader.ptr;
1441 CvPoint* prev_pt = (CvPoint*)reader.prev_elem;
1442
1443 dx = (float)pt->x - (float)prev_pt->x;
1444 dy = (float)pt->y - (float)prev_pt->y;
1445 }
1446 else
1447 {
1448 CvPoint2D32f* pt = (CvPoint2D32f*)reader.ptr;
1449 CvPoint2D32f* prev_pt = (CvPoint2D32f*)reader.prev_elem;
1450
1451 dx = pt->x - prev_pt->x;
1452 dy = pt->y - prev_pt->y;
1453 }
1454
1455 reader.prev_elem = reader.ptr;
1456 CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
1457 // Bugfix by Axel at rubico.com 2010-03-22, affects closed slices only
1458 // wraparound not handled by CV_NEXT_SEQ_ELEM
1459 if( is_closed && i == count - 2 )
1460 cvSetSeqReaderPos( &reader, slice.start_index );
1461
1462 buffer.data.fl[j] = dx * dx + dy * dy;
1463 if( ++j == N || i == count - 1 )
1464 {
1465 buffer.cols = j;
1466 cvPow( &buffer, &buffer, 0.5 );
1467 for( ; j > 0; j-- )
1468 perimeter += buffer.data.fl[j-1];
1469 }
1470 }
1471 }
1472
1473 return perimeter;
1474 }
1475
1476
1477 CV_IMPL CvBox2D
cvFitEllipse2(const CvArr * array)1478 cvFitEllipse2( const CvArr* array )
1479 {
1480 cv::AutoBuffer<double> abuf;
1481 cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf);
1482 return cvBox2D(cv::fitEllipse(points));
1483 }
1484
1485 /* Calculates bounding rectangle of a point set or retrieves already calculated */
1486 CV_IMPL CvRect
cvBoundingRect(CvArr * array,int update)1487 cvBoundingRect( CvArr* array, int update )
1488 {
1489 cv::Rect rect;
1490 CvContour contour_header;
1491 CvSeq* ptseq = 0;
1492 CvSeqBlock block;
1493
1494 CvMat stub, *mat = 0;
1495 int calculate = update;
1496
1497 if( CV_IS_SEQ( array ))
1498 {
1499 ptseq = (CvSeq*)array;
1500 if( !CV_IS_SEQ_POINT_SET( ptseq ))
1501 CV_Error( CV_StsBadArg, "Unsupported sequence type" );
1502
1503 if( ptseq->header_size < (int)sizeof(CvContour))
1504 {
1505 update = 0;
1506 calculate = 1;
1507 }
1508 }
1509 else
1510 {
1511 mat = cvGetMat( array, &stub );
1512 if( CV_MAT_TYPE(mat->type) == CV_32SC2 ||
1513 CV_MAT_TYPE(mat->type) == CV_32FC2 )
1514 {
1515 ptseq = cvPointSeqFromMat(CV_SEQ_KIND_GENERIC, mat, &contour_header, &block);
1516 mat = 0;
1517 }
1518 else if( CV_MAT_TYPE(mat->type) != CV_8UC1 &&
1519 CV_MAT_TYPE(mat->type) != CV_8SC1 )
1520 CV_Error( CV_StsUnsupportedFormat,
1521 "The image/matrix format is not supported by the function" );
1522 update = 0;
1523 calculate = 1;
1524 }
1525
1526 if( !calculate )
1527 return ((CvContour*)ptseq)->rect;
1528
1529 if( mat )
1530 {
1531 rect = cvRect(cv::maskBoundingRect(cv::cvarrToMat(mat)));
1532 }
1533 else if( ptseq->total )
1534 {
1535 cv::AutoBuffer<double> abuf;
1536 rect = cvRect(cv::pointSetBoundingRect(cv::cvarrToMat(ptseq, false, false, 0, &abuf)));
1537 }
1538 if( update )
1539 ((CvContour*)ptseq)->rect = cvRect(rect);
1540 return cvRect(rect);
1541 }
1542
1543 /* End of file. */
1544