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 43 namespace cv 44 { 45 46 static int intersectLines( double x1, double dx1, double y1, double dy1, 47 double x2, double dx2, double y2, double dy2, double *t2 ) 48 { 49 double d = dx1 * dy2 - dx2 * dy1; 50 int result = -1; 51 52 if( d != 0 ) 53 { 54 *t2 = ((x2 - x1) * dy1 - (y2 - y1) * dx1) / d; 55 result = 0; 56 } 57 return result; 58 } 59 60 static bool findCircle( Point2f pt0, Point2f pt1, Point2f pt2, 61 Point2f* center, float* radius ) 62 { 63 double x1 = (pt0.x + pt1.x) * 0.5; 64 double dy1 = pt0.x - pt1.x; 65 double x2 = (pt1.x + pt2.x) * 0.5; 66 double dy2 = pt1.x - pt2.x; 67 double y1 = (pt0.y + pt1.y) * 0.5; 68 double dx1 = pt1.y - pt0.y; 69 double y2 = (pt1.y + pt2.y) * 0.5; 70 double dx2 = pt2.y - pt1.y; 71 double t = 0; 72 73 if( intersectLines( x1, dx1, y1, dy1, x2, dx2, y2, dy2, &t ) >= 0 ) 74 { 75 center->x = (float) (x2 + dx2 * t); 76 center->y = (float) (y2 + dy2 * t); 77 *radius = (float)norm(*center - pt0); 78 return true; 79 } 80 81 center->x = center->y = 0.f; 82 *radius = 0; 83 return false; 84 } 85 86 87 static double pointInCircle( Point2f pt, Point2f center, float radius ) 88 { 89 double dx = pt.x - center.x; 90 double dy = pt.y - center.y; 91 return (double)radius*radius - dx*dx - dy*dy; 92 } 93 94 95 static int findEnslosingCicle4pts_32f( Point2f* pts, Point2f& _center, float& _radius ) 96 { 97 int shuffles[4][4] = { {0, 1, 2, 3}, {0, 1, 3, 2}, {2, 3, 0, 1}, {2, 3, 1, 0} }; 98 99 int idxs[4] = { 0, 1, 2, 3 }; 100 int i, j, k = 1, mi = 0; 101 float max_dist = 0; 102 Point2f center; 103 Point2f min_center; 104 float radius, min_radius = FLT_MAX; 105 Point2f res_pts[4]; 106 107 center = min_center = pts[0]; 108 radius = 1.f; 109 110 for( i = 0; i < 4; i++ ) 111 for( j = i + 1; j < 4; j++ ) 112 { 113 float dist = (float)norm(pts[i] - pts[j]); 114 115 if( max_dist < dist ) 116 { 117 max_dist = dist; 118 idxs[0] = i; 119 idxs[1] = j; 120 } 121 } 122 123 if( max_dist > 0 ) 124 { 125 k = 2; 126 for( i = 0; i < 4; i++ ) 127 { 128 for( j = 0; j < k; j++ ) 129 if( i == idxs[j] ) 130 break; 131 if( j == k ) 132 idxs[k++] = i; 133 } 134 135 center = Point2f( (pts[idxs[0]].x + pts[idxs[1]].x)*0.5f, 136 (pts[idxs[0]].y + pts[idxs[1]].y)*0.5f ); 137 radius = (float)(norm(pts[idxs[0]] - center)*1.03); 138 if( radius < 1.f ) 139 radius = 1.f; 140 141 if( pointInCircle( pts[idxs[2]], center, radius ) >= 0 && 142 pointInCircle( pts[idxs[3]], center, radius ) >= 0 ) 143 { 144 k = 2; //rand()%2+2; 145 } 146 else 147 { 148 mi = -1; 149 for( i = 0; i < 4; i++ ) 150 { 151 if( findCircle( pts[shuffles[i][0]], pts[shuffles[i][1]], 152 pts[shuffles[i][2]], ¢er, &radius ) ) 153 { 154 radius *= 1.03f; 155 if( radius < 2.f ) 156 radius = 2.f; 157 158 if( pointInCircle( pts[shuffles[i][3]], center, radius ) >= 0 && 159 min_radius > radius ) 160 { 161 min_radius = radius; 162 min_center = center; 163 mi = i; 164 } 165 } 166 } 167 CV_Assert( mi >= 0 ); 168 if( mi < 0 ) 169 mi = 0; 170 k = 3; 171 center = min_center; 172 radius = min_radius; 173 for( i = 0; i < 4; i++ ) 174 idxs[i] = shuffles[mi][i]; 175 } 176 } 177 178 _center = center; 179 _radius = radius; 180 181 /* reorder output points */ 182 for( i = 0; i < 4; i++ ) 183 res_pts[i] = pts[idxs[i]]; 184 185 for( i = 0; i < 4; i++ ) 186 { 187 pts[i] = res_pts[i]; 188 CV_Assert( pointInCircle( pts[i], center, radius ) >= 0 ); 189 } 190 191 return k; 192 } 193 194 } 195 196 void cv::minEnclosingCircle( InputArray _points, Point2f& _center, float& _radius ) 197 { 198 int max_iters = 100; 199 const float eps = FLT_EPSILON*2; 200 bool result = false; 201 Mat points = _points.getMat(); 202 int i, j, k, count = points.checkVector(2); 203 int depth = points.depth(); 204 Point2f center; 205 float radius = 0.f; 206 CV_Assert(count >= 0 && (depth == CV_32F || depth == CV_32S)); 207 208 _center.x = _center.y = 0.f; 209 _radius = 0.f; 210 211 if( count == 0 ) 212 return; 213 214 bool is_float = depth == CV_32F; 215 const Point* ptsi = points.ptr<Point>(); 216 const Point2f* ptsf = points.ptr<Point2f>(); 217 218 Point2f pt = is_float ? ptsf[0] : Point2f((float)ptsi[0].x,(float)ptsi[0].y); 219 Point2f pts[4] = {pt, pt, pt, pt}; 220 221 for( i = 1; i < count; i++ ) 222 { 223 pt = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); 224 225 if( pt.x < pts[0].x ) 226 pts[0] = pt; 227 if( pt.x > pts[1].x ) 228 pts[1] = pt; 229 if( pt.y < pts[2].y ) 230 pts[2] = pt; 231 if( pt.y > pts[3].y ) 232 pts[3] = pt; 233 } 234 235 for( k = 0; k < max_iters; k++ ) 236 { 237 double min_delta = 0, delta; 238 Point2f farAway(0,0); 239 /*only for first iteration because the alg is repared at the loop's foot*/ 240 if( k == 0 ) 241 findEnslosingCicle4pts_32f( pts, center, radius ); 242 243 for( i = 0; i < count; i++ ) 244 { 245 pt = is_float ? ptsf[i] : Point2f((float)ptsi[i].x,(float)ptsi[i].y); 246 247 delta = pointInCircle( pt, center, radius ); 248 if( delta < min_delta ) 249 { 250 min_delta = delta; 251 farAway = pt; 252 } 253 } 254 result = min_delta >= 0; 255 if( result ) 256 break; 257 258 Point2f ptsCopy[4]; 259 // find good replacement partner for the point which is at most far away, 260 // starting with the one that lays in the actual circle (i=3) 261 for( i = 3; i >= 0; i-- ) 262 { 263 for( j = 0; j < 4; j++ ) 264 ptsCopy[j] = i != j ? pts[j] : farAway; 265 266 findEnslosingCicle4pts_32f( ptsCopy, center, radius ); 267 if( pointInCircle( pts[i], center, radius ) >= 0) 268 { 269 // replaced one again in the new circle? 270 pts[i] = farAway; 271 break; 272 } 273 } 274 } 275 276 if( !result ) 277 { 278 radius = 0.f; 279 for( i = 0; i < count; i++ ) 280 { 281 pt = is_float ? ptsf[i] : Point2f((float)ptsi[i].x,(float)ptsi[i].y); 282 float dx = center.x - pt.x, dy = center.y - pt.y; 283 float t = dx*dx + dy*dy; 284 radius = MAX(radius, t); 285 } 286 287 radius = (float)(std::sqrt(radius)*(1 + eps)); 288 } 289 290 _center = center; 291 _radius = radius; 292 } 293 294 295 // calculates length of a curve (e.g. contour perimeter) 296 double cv::arcLength( InputArray _curve, bool is_closed ) 297 { 298 Mat curve = _curve.getMat(); 299 int count = curve.checkVector(2); 300 int depth = curve.depth(); 301 CV_Assert( count >= 0 && (depth == CV_32F || depth == CV_32S)); 302 double perimeter = 0; 303 304 int i, j = 0; 305 const int N = 16; 306 float buf[N]; 307 308 if( count <= 1 ) 309 return 0.; 310 311 bool is_float = depth == CV_32F; 312 int last = is_closed ? count-1 : 0; 313 const Point* pti = curve.ptr<Point>(); 314 const Point2f* ptf = curve.ptr<Point2f>(); 315 316 Point2f prev = is_float ? ptf[last] : Point2f((float)pti[last].x,(float)pti[last].y); 317 318 for( i = 0; i < count; i++ ) 319 { 320 Point2f p = is_float ? ptf[i] : Point2f((float)pti[i].x,(float)pti[i].y); 321 float dx = p.x - prev.x, dy = p.y - prev.y; 322 buf[j] = dx*dx + dy*dy; 323 324 if( ++j == N || i == count-1 ) 325 { 326 Mat bufmat(1, j, CV_32F, buf); 327 sqrt(bufmat, bufmat); 328 for( ; j > 0; j-- ) 329 perimeter += buf[j-1]; 330 } 331 prev = p; 332 } 333 334 return perimeter; 335 } 336 337 // area of a whole sequence 338 double cv::contourArea( InputArray _contour, bool oriented ) 339 { 340 Mat contour = _contour.getMat(); 341 int npoints = contour.checkVector(2); 342 int depth = contour.depth(); 343 CV_Assert(npoints >= 0 && (depth == CV_32F || depth == CV_32S)); 344 345 if( npoints == 0 ) 346 return 0.; 347 348 double a00 = 0; 349 bool is_float = depth == CV_32F; 350 const Point* ptsi = contour.ptr<Point>(); 351 const Point2f* ptsf = contour.ptr<Point2f>(); 352 Point2f prev = is_float ? ptsf[npoints-1] : Point2f((float)ptsi[npoints-1].x, (float)ptsi[npoints-1].y); 353 354 for( int i = 0; i < npoints; i++ ) 355 { 356 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); 357 a00 += (double)prev.x * p.y - (double)prev.y * p.x; 358 prev = p; 359 } 360 361 a00 *= 0.5; 362 if( !oriented ) 363 a00 = fabs(a00); 364 365 return a00; 366 } 367 368 369 cv::RotatedRect cv::fitEllipse( InputArray _points ) 370 { 371 Mat points = _points.getMat(); 372 int i, n = points.checkVector(2); 373 int depth = points.depth(); 374 CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S)); 375 376 RotatedRect box; 377 378 if( n < 5 ) 379 CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" ); 380 381 // New fitellipse algorithm, contributed by Dr. Daniel Weiss 382 Point2f c(0,0); 383 double gfp[5], rp[5], t; 384 const double min_eps = 1e-8; 385 bool is_float = depth == CV_32F; 386 const Point* ptsi = points.ptr<Point>(); 387 const Point2f* ptsf = points.ptr<Point2f>(); 388 389 AutoBuffer<double> _Ad(n*5), _bd(n); 390 double *Ad = _Ad, *bd = _bd; 391 392 // first fit for parameters A - E 393 Mat A( n, 5, CV_64F, Ad ); 394 Mat b( n, 1, CV_64F, bd ); 395 Mat x( 5, 1, CV_64F, gfp ); 396 397 for( i = 0; i < n; i++ ) 398 { 399 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); 400 c += p; 401 } 402 c.x /= n; 403 c.y /= n; 404 405 for( i = 0; i < n; i++ ) 406 { 407 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); 408 p -= c; 409 410 bd[i] = 10000.0; // 1.0? 411 Ad[i*5] = -(double)p.x * p.x; // A - C signs inverted as proposed by APP 412 Ad[i*5 + 1] = -(double)p.y * p.y; 413 Ad[i*5 + 2] = -(double)p.x * p.y; 414 Ad[i*5 + 3] = p.x; 415 Ad[i*5 + 4] = p.y; 416 } 417 418 solve(A, b, x, DECOMP_SVD); 419 420 // now use general-form parameters A - E to find the ellipse center: 421 // differentiate general form wrt x/y to get two equations for cx and cy 422 A = Mat( 2, 2, CV_64F, Ad ); 423 b = Mat( 2, 1, CV_64F, bd ); 424 x = Mat( 2, 1, CV_64F, rp ); 425 Ad[0] = 2 * gfp[0]; 426 Ad[1] = Ad[2] = gfp[2]; 427 Ad[3] = 2 * gfp[1]; 428 bd[0] = gfp[3]; 429 bd[1] = gfp[4]; 430 solve( A, b, x, DECOMP_SVD ); 431 432 // re-fit for parameters A - C with those center coordinates 433 A = Mat( n, 3, CV_64F, Ad ); 434 b = Mat( n, 1, CV_64F, bd ); 435 x = Mat( 3, 1, CV_64F, gfp ); 436 for( i = 0; i < n; i++ ) 437 { 438 Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); 439 p -= c; 440 bd[i] = 1.0; 441 Ad[i * 3] = (p.x - rp[0]) * (p.x - rp[0]); 442 Ad[i * 3 + 1] = (p.y - rp[1]) * (p.y - rp[1]); 443 Ad[i * 3 + 2] = (p.x - rp[0]) * (p.y - rp[1]); 444 } 445 solve(A, b, x, DECOMP_SVD); 446 447 // store angle and radii 448 rp[4] = -0.5 * atan2(gfp[2], gfp[1] - gfp[0]); // convert from APP angle usage 449 if( fabs(gfp[2]) > min_eps ) 450 t = gfp[2]/sin(-2.0 * rp[4]); 451 else // ellipse is rotated by an integer multiple of pi/2 452 t = gfp[1] - gfp[0]; 453 rp[2] = fabs(gfp[0] + gfp[1] - t); 454 if( rp[2] > min_eps ) 455 rp[2] = std::sqrt(2.0 / rp[2]); 456 rp[3] = fabs(gfp[0] + gfp[1] + t); 457 if( rp[3] > min_eps ) 458 rp[3] = std::sqrt(2.0 / rp[3]); 459 460 box.center.x = (float)rp[0] + c.x; 461 box.center.y = (float)rp[1] + c.y; 462 box.size.width = (float)(rp[2]*2); 463 box.size.height = (float)(rp[3]*2); 464 if( box.size.width > box.size.height ) 465 { 466 float tmp; 467 CV_SWAP( box.size.width, box.size.height, tmp ); 468 box.angle = (float)(90 + rp[4]*180/CV_PI); 469 } 470 if( box.angle < -180 ) 471 box.angle += 360; 472 if( box.angle > 360 ) 473 box.angle -= 360; 474 475 return box; 476 } 477 478 479 namespace cv 480 { 481 482 // Calculates bounding rectagnle of a point set or retrieves already calculated 483 static Rect pointSetBoundingRect( const Mat& points ) 484 { 485 int npoints = points.checkVector(2); 486 int depth = points.depth(); 487 CV_Assert(npoints >= 0 && (depth == CV_32F || depth == CV_32S)); 488 489 int xmin = 0, ymin = 0, xmax = -1, ymax = -1, i; 490 bool is_float = depth == CV_32F; 491 492 if( npoints == 0 ) 493 return Rect(); 494 495 const Point* pts = points.ptr<Point>(); 496 Point pt = pts[0]; 497 498 #if CV_SSE4_2 499 if(cv::checkHardwareSupport(CV_CPU_SSE4_2)) 500 { 501 if( !is_float ) 502 { 503 __m128i minval, maxval; 504 minval = maxval = _mm_loadl_epi64((const __m128i*)(&pt)); //min[0]=pt.x, min[1]=pt.y 505 506 for( i = 1; i < npoints; i++ ) 507 { 508 __m128i ptXY = _mm_loadl_epi64((const __m128i*)&pts[i]); 509 minval = _mm_min_epi32(ptXY, minval); 510 maxval = _mm_max_epi32(ptXY, maxval); 511 } 512 xmin = _mm_cvtsi128_si32(minval); 513 ymin = _mm_cvtsi128_si32(_mm_srli_si128(minval, 4)); 514 xmax = _mm_cvtsi128_si32(maxval); 515 ymax = _mm_cvtsi128_si32(_mm_srli_si128(maxval, 4)); 516 } 517 else 518 { 519 __m128 minvalf, maxvalf, z = _mm_setzero_ps(), ptXY = _mm_setzero_ps(); 520 minvalf = maxvalf = _mm_loadl_pi(z, (const __m64*)(&pt)); 521 522 for( i = 1; i < npoints; i++ ) 523 { 524 ptXY = _mm_loadl_pi(ptXY, (const __m64*)&pts[i]); 525 526 minvalf = _mm_min_ps(minvalf, ptXY); 527 maxvalf = _mm_max_ps(maxvalf, ptXY); 528 } 529 530 float xyminf[2], xymaxf[2]; 531 _mm_storel_pi((__m64*)xyminf, minvalf); 532 _mm_storel_pi((__m64*)xymaxf, maxvalf); 533 xmin = cvFloor(xyminf[0]); 534 ymin = cvFloor(xyminf[1]); 535 xmax = cvFloor(xymaxf[0]); 536 ymax = cvFloor(xymaxf[1]); 537 } 538 } 539 else 540 #endif 541 { 542 if( !is_float ) 543 { 544 xmin = xmax = pt.x; 545 ymin = ymax = pt.y; 546 547 for( i = 1; i < npoints; i++ ) 548 { 549 pt = pts[i]; 550 551 if( xmin > pt.x ) 552 xmin = pt.x; 553 554 if( xmax < pt.x ) 555 xmax = pt.x; 556 557 if( ymin > pt.y ) 558 ymin = pt.y; 559 560 if( ymax < pt.y ) 561 ymax = pt.y; 562 } 563 } 564 else 565 { 566 Cv32suf v; 567 // init values 568 xmin = xmax = CV_TOGGLE_FLT(pt.x); 569 ymin = ymax = CV_TOGGLE_FLT(pt.y); 570 571 for( i = 1; i < npoints; i++ ) 572 { 573 pt = pts[i]; 574 pt.x = CV_TOGGLE_FLT(pt.x); 575 pt.y = CV_TOGGLE_FLT(pt.y); 576 577 if( xmin > pt.x ) 578 xmin = pt.x; 579 580 if( xmax < pt.x ) 581 xmax = pt.x; 582 583 if( ymin > pt.y ) 584 ymin = pt.y; 585 586 if( ymax < pt.y ) 587 ymax = pt.y; 588 } 589 590 v.i = CV_TOGGLE_FLT(xmin); xmin = cvFloor(v.f); 591 v.i = CV_TOGGLE_FLT(ymin); ymin = cvFloor(v.f); 592 // because right and bottom sides of the bounding rectangle are not inclusive 593 // (note +1 in width and height calculation below), cvFloor is used here instead of cvCeil 594 v.i = CV_TOGGLE_FLT(xmax); xmax = cvFloor(v.f); 595 v.i = CV_TOGGLE_FLT(ymax); ymax = cvFloor(v.f); 596 } 597 } 598 599 return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1); 600 } 601 602 603 static Rect maskBoundingRect( const Mat& img ) 604 { 605 CV_Assert( img.depth() <= CV_8S && img.channels() == 1 ); 606 607 Size size = img.size(); 608 int xmin = size.width, ymin = -1, xmax = -1, ymax = -1, i, j, k; 609 610 for( i = 0; i < size.height; i++ ) 611 { 612 const uchar* _ptr = img.ptr(i); 613 const uchar* ptr = (const uchar*)alignPtr(_ptr, 4); 614 int have_nz = 0, k_min, offset = (int)(ptr - _ptr); 615 j = 0; 616 offset = MIN(offset, size.width); 617 for( ; j < offset; j++ ) 618 if( _ptr[j] ) 619 { 620 have_nz = 1; 621 break; 622 } 623 if( j < offset ) 624 { 625 if( j < xmin ) 626 xmin = j; 627 if( j > xmax ) 628 xmax = j; 629 } 630 if( offset < size.width ) 631 { 632 xmin -= offset; 633 xmax -= offset; 634 size.width -= offset; 635 j = 0; 636 for( ; j <= xmin - 4; j += 4 ) 637 if( *((int*)(ptr+j)) ) 638 break; 639 for( ; j < xmin; j++ ) 640 if( ptr[j] ) 641 { 642 xmin = j; 643 if( j > xmax ) 644 xmax = j; 645 have_nz = 1; 646 break; 647 } 648 k_min = MAX(j-1, xmax); 649 k = size.width - 1; 650 for( ; k > k_min && (k&3) != 3; k-- ) 651 if( ptr[k] ) 652 break; 653 if( k > k_min && (k&3) == 3 ) 654 { 655 for( ; k > k_min+3; k -= 4 ) 656 if( *((int*)(ptr+k-3)) ) 657 break; 658 } 659 for( ; k > k_min; k-- ) 660 if( ptr[k] ) 661 { 662 xmax = k; 663 have_nz = 1; 664 break; 665 } 666 if( !have_nz ) 667 { 668 j &= ~3; 669 for( ; j <= k - 3; j += 4 ) 670 if( *((int*)(ptr+j)) ) 671 break; 672 for( ; j <= k; j++ ) 673 if( ptr[j] ) 674 { 675 have_nz = 1; 676 break; 677 } 678 } 679 xmin += offset; 680 xmax += offset; 681 size.width += offset; 682 } 683 if( have_nz ) 684 { 685 if( ymin < 0 ) 686 ymin = i; 687 ymax = i; 688 } 689 } 690 691 if( xmin >= size.width ) 692 xmin = ymin = 0; 693 return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1); 694 } 695 696 } 697 698 cv::Rect cv::boundingRect(InputArray array) 699 { 700 Mat m = array.getMat(); 701 return m.depth() <= CV_8U ? maskBoundingRect(m) : pointSetBoundingRect(m); 702 } 703 704 ////////////////////////////////////////////// C API /////////////////////////////////////////// 705 706 CV_IMPL int 707 cvMinEnclosingCircle( const void* array, CvPoint2D32f * _center, float *_radius ) 708 { 709 cv::AutoBuffer<double> abuf; 710 cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf); 711 cv::Point2f center; 712 float radius; 713 714 cv::minEnclosingCircle(points, center, radius); 715 if(_center) 716 *_center = center; 717 if(_radius) 718 *_radius = radius; 719 return 1; 720 } 721 722 static void 723 icvMemCopy( double **buf1, double **buf2, double **buf3, int *b_max ) 724 { 725 CV_Assert( (*buf1 != NULL || *buf2 != NULL) && *buf3 != NULL ); 726 727 int bb = *b_max; 728 if( *buf2 == NULL ) 729 { 730 *b_max = 2 * (*b_max); 731 *buf2 = (double *)cvAlloc( (*b_max) * sizeof( double )); 732 733 memcpy( *buf2, *buf3, bb * sizeof( double )); 734 735 *buf3 = *buf2; 736 cvFree( buf1 ); 737 *buf1 = NULL; 738 } 739 else 740 { 741 *b_max = 2 * (*b_max); 742 *buf1 = (double *) cvAlloc( (*b_max) * sizeof( double )); 743 744 memcpy( *buf1, *buf3, bb * sizeof( double )); 745 746 *buf3 = *buf1; 747 cvFree( buf2 ); 748 *buf2 = NULL; 749 } 750 } 751 752 753 /* area of a contour sector */ 754 static double icvContourSecArea( CvSeq * contour, CvSlice slice ) 755 { 756 CvPoint pt; /* pointer to points */ 757 CvPoint pt_s, pt_e; /* first and last points */ 758 CvSeqReader reader; /* points reader of contour */ 759 760 int p_max = 2, p_ind; 761 int lpt, flag, i; 762 double a00; /* unnormalized moments m00 */ 763 double xi, yi, xi_1, yi_1, x0, y0, dxy, sk, sk1, t; 764 double x_s, y_s, nx, ny, dx, dy, du, dv; 765 double eps = 1.e-5; 766 double *p_are1, *p_are2, *p_are; 767 double area = 0; 768 769 CV_Assert( contour != NULL && CV_IS_SEQ_POINT_SET( contour )); 770 771 lpt = cvSliceLength( slice, contour ); 772 /*if( n2 >= n1 ) 773 lpt = n2 - n1 + 1; 774 else 775 lpt = contour->total - n1 + n2 + 1;*/ 776 777 if( contour->total <= 0 || lpt <= 2 ) 778 return 0.; 779 780 a00 = x0 = y0 = xi_1 = yi_1 = 0; 781 sk1 = 0; 782 flag = 0; 783 dxy = 0; 784 p_are1 = (double *) cvAlloc( p_max * sizeof( double )); 785 786 p_are = p_are1; 787 p_are2 = NULL; 788 789 cvStartReadSeq( contour, &reader, 0 ); 790 cvSetSeqReaderPos( &reader, slice.start_index ); 791 CV_READ_SEQ_ELEM( pt_s, reader ); 792 p_ind = 0; 793 cvSetSeqReaderPos( &reader, slice.end_index ); 794 CV_READ_SEQ_ELEM( pt_e, reader ); 795 796 /* normal coefficients */ 797 nx = pt_s.y - pt_e.y; 798 ny = pt_e.x - pt_s.x; 799 cvSetSeqReaderPos( &reader, slice.start_index ); 800 801 while( lpt-- > 0 ) 802 { 803 CV_READ_SEQ_ELEM( pt, reader ); 804 805 if( flag == 0 ) 806 { 807 xi_1 = (double) pt.x; 808 yi_1 = (double) pt.y; 809 x0 = xi_1; 810 y0 = yi_1; 811 sk1 = 0; 812 flag = 1; 813 } 814 else 815 { 816 xi = (double) pt.x; 817 yi = (double) pt.y; 818 819 /**************** edges intersection examination **************************/ 820 sk = nx * (xi - pt_s.x) + ny * (yi - pt_s.y); 821 if( (fabs( sk ) < eps && lpt > 0) || sk * sk1 < -eps ) 822 { 823 if( fabs( sk ) < eps ) 824 { 825 dxy = xi_1 * yi - xi * yi_1; 826 a00 = a00 + dxy; 827 dxy = xi * y0 - x0 * yi; 828 a00 = a00 + dxy; 829 830 if( p_ind >= p_max ) 831 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max ); 832 833 p_are[p_ind] = a00 / 2.; 834 p_ind++; 835 a00 = 0; 836 sk1 = 0; 837 x0 = xi; 838 y0 = yi; 839 dxy = 0; 840 } 841 else 842 { 843 /* define intersection point */ 844 dv = yi - yi_1; 845 du = xi - xi_1; 846 dx = ny; 847 dy = -nx; 848 if( fabs( du ) > eps ) 849 t = ((yi_1 - pt_s.y) * du + dv * (pt_s.x - xi_1)) / 850 (du * dy - dx * dv); 851 else 852 t = (xi_1 - pt_s.x) / dx; 853 if( t > eps && t < 1 - eps ) 854 { 855 x_s = pt_s.x + t * dx; 856 y_s = pt_s.y + t * dy; 857 dxy = xi_1 * y_s - x_s * yi_1; 858 a00 += dxy; 859 dxy = x_s * y0 - x0 * y_s; 860 a00 += dxy; 861 if( p_ind >= p_max ) 862 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max ); 863 864 p_are[p_ind] = a00 / 2.; 865 p_ind++; 866 867 a00 = 0; 868 sk1 = 0; 869 x0 = x_s; 870 y0 = y_s; 871 dxy = x_s * yi - xi * y_s; 872 } 873 } 874 } 875 else 876 dxy = xi_1 * yi - xi * yi_1; 877 878 a00 += dxy; 879 xi_1 = xi; 880 yi_1 = yi; 881 sk1 = sk; 882 883 } 884 } 885 886 xi = x0; 887 yi = y0; 888 dxy = xi_1 * yi - xi * yi_1; 889 890 a00 += dxy; 891 892 if( p_ind >= p_max ) 893 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max ); 894 895 p_are[p_ind] = a00 / 2.; 896 p_ind++; 897 898 // common area calculation 899 area = 0; 900 for( i = 0; i < p_ind; i++ ) 901 area += fabs( p_are[i] ); 902 903 if( p_are1 != NULL ) 904 cvFree( &p_are1 ); 905 else if( p_are2 != NULL ) 906 cvFree( &p_are2 ); 907 908 return area; 909 } 910 911 912 /* external contour area function */ 913 CV_IMPL double 914 cvContourArea( const void *array, CvSlice slice, int oriented ) 915 { 916 double area = 0; 917 918 CvContour contour_header; 919 CvSeq* contour = 0; 920 CvSeqBlock block; 921 922 if( CV_IS_SEQ( array )) 923 { 924 contour = (CvSeq*)array; 925 if( !CV_IS_SEQ_POLYLINE( contour )) 926 CV_Error( CV_StsBadArg, "Unsupported sequence type" ); 927 } 928 else 929 { 930 contour = cvPointSeqFromMat( CV_SEQ_KIND_CURVE, array, &contour_header, &block ); 931 } 932 933 if( cvSliceLength( slice, contour ) == contour->total ) 934 { 935 cv::AutoBuffer<double> abuf; 936 cv::Mat points = cv::cvarrToMat(contour, false, false, 0, &abuf); 937 return cv::contourArea( points, oriented !=0 ); 938 } 939 940 if( CV_SEQ_ELTYPE( contour ) != CV_32SC2 ) 941 CV_Error( CV_StsUnsupportedFormat, 942 "Only curves with integer coordinates are supported in case of contour slice" ); 943 area = icvContourSecArea( contour, slice ); 944 return oriented ? area : fabs(area); 945 } 946 947 948 /* calculates length of a curve (e.g. contour perimeter) */ 949 CV_IMPL double 950 cvArcLength( const void *array, CvSlice slice, int is_closed ) 951 { 952 double perimeter = 0; 953 954 int i, j = 0, count; 955 const int N = 16; 956 float buf[N]; 957 CvMat buffer = cvMat( 1, N, CV_32F, buf ); 958 CvSeqReader reader; 959 CvContour contour_header; 960 CvSeq* contour = 0; 961 CvSeqBlock block; 962 963 if( CV_IS_SEQ( array )) 964 { 965 contour = (CvSeq*)array; 966 if( !CV_IS_SEQ_POLYLINE( contour )) 967 CV_Error( CV_StsBadArg, "Unsupported sequence type" ); 968 if( is_closed < 0 ) 969 is_closed = CV_IS_SEQ_CLOSED( contour ); 970 } 971 else 972 { 973 is_closed = is_closed > 0; 974 contour = cvPointSeqFromMat( 975 CV_SEQ_KIND_CURVE | (is_closed ? CV_SEQ_FLAG_CLOSED : 0), 976 array, &contour_header, &block ); 977 } 978 979 if( contour->total > 1 ) 980 { 981 int is_float = CV_SEQ_ELTYPE( contour ) == CV_32FC2; 982 983 cvStartReadSeq( contour, &reader, 0 ); 984 cvSetSeqReaderPos( &reader, slice.start_index ); 985 count = cvSliceLength( slice, contour ); 986 987 count -= !is_closed && count == contour->total; 988 989 // scroll the reader by 1 point 990 reader.prev_elem = reader.ptr; 991 CV_NEXT_SEQ_ELEM( sizeof(CvPoint), reader ); 992 993 for( i = 0; i < count; i++ ) 994 { 995 float dx, dy; 996 997 if( !is_float ) 998 { 999 CvPoint* pt = (CvPoint*)reader.ptr; 1000 CvPoint* prev_pt = (CvPoint*)reader.prev_elem; 1001 1002 dx = (float)pt->x - (float)prev_pt->x; 1003 dy = (float)pt->y - (float)prev_pt->y; 1004 } 1005 else 1006 { 1007 CvPoint2D32f* pt = (CvPoint2D32f*)reader.ptr; 1008 CvPoint2D32f* prev_pt = (CvPoint2D32f*)reader.prev_elem; 1009 1010 dx = pt->x - prev_pt->x; 1011 dy = pt->y - prev_pt->y; 1012 } 1013 1014 reader.prev_elem = reader.ptr; 1015 CV_NEXT_SEQ_ELEM( contour->elem_size, reader ); 1016 // Bugfix by Axel at rubico.com 2010-03-22, affects closed slices only 1017 // wraparound not handled by CV_NEXT_SEQ_ELEM 1018 if( is_closed && i == count - 2 ) 1019 cvSetSeqReaderPos( &reader, slice.start_index ); 1020 1021 buffer.data.fl[j] = dx * dx + dy * dy; 1022 if( ++j == N || i == count - 1 ) 1023 { 1024 buffer.cols = j; 1025 cvPow( &buffer, &buffer, 0.5 ); 1026 for( ; j > 0; j-- ) 1027 perimeter += buffer.data.fl[j-1]; 1028 } 1029 } 1030 } 1031 1032 return perimeter; 1033 } 1034 1035 1036 CV_IMPL CvBox2D 1037 cvFitEllipse2( const CvArr* array ) 1038 { 1039 cv::AutoBuffer<double> abuf; 1040 cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf); 1041 return cv::fitEllipse(points); 1042 } 1043 1044 /* Calculates bounding rectagnle of a point set or retrieves already calculated */ 1045 CV_IMPL CvRect 1046 cvBoundingRect( CvArr* array, int update ) 1047 { 1048 CvRect rect; 1049 CvContour contour_header; 1050 CvSeq* ptseq = 0; 1051 CvSeqBlock block; 1052 1053 CvMat stub, *mat = 0; 1054 int calculate = update; 1055 1056 if( CV_IS_SEQ( array )) 1057 { 1058 ptseq = (CvSeq*)array; 1059 if( !CV_IS_SEQ_POINT_SET( ptseq )) 1060 CV_Error( CV_StsBadArg, "Unsupported sequence type" ); 1061 1062 if( ptseq->header_size < (int)sizeof(CvContour)) 1063 { 1064 update = 0; 1065 calculate = 1; 1066 } 1067 } 1068 else 1069 { 1070 mat = cvGetMat( array, &stub ); 1071 if( CV_MAT_TYPE(mat->type) == CV_32SC2 || 1072 CV_MAT_TYPE(mat->type) == CV_32FC2 ) 1073 { 1074 ptseq = cvPointSeqFromMat(CV_SEQ_KIND_GENERIC, mat, &contour_header, &block); 1075 mat = 0; 1076 } 1077 else if( CV_MAT_TYPE(mat->type) != CV_8UC1 && 1078 CV_MAT_TYPE(mat->type) != CV_8SC1 ) 1079 CV_Error( CV_StsUnsupportedFormat, 1080 "The image/matrix format is not supported by the function" ); 1081 update = 0; 1082 calculate = 1; 1083 } 1084 1085 if( !calculate ) 1086 return ((CvContour*)ptseq)->rect; 1087 1088 if( mat ) 1089 { 1090 rect = cv::maskBoundingRect(cv::cvarrToMat(mat)); 1091 } 1092 else if( ptseq->total ) 1093 { 1094 cv::AutoBuffer<double> abuf; 1095 rect = cv::pointSetBoundingRect(cv::cvarrToMat(ptseq, false, false, 0, &abuf)); 1096 } 1097 if( update ) 1098 ((CvContour*)ptseq)->rect = rect; 1099 return rect; 1100 } 1101 1102 1103 /* End of file. */ 1104
