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 "_cv.h" 42 43 /* calculates length of a curve (e.g. contour perimeter) */ 44 CV_IMPL double 45 cvArcLength( const void *array, CvSlice slice, int is_closed ) 46 { 47 double perimeter = 0; 48 49 CV_FUNCNAME( "cvArcLength" ); 50 51 __BEGIN__; 52 53 int i, j = 0, count; 54 const int N = 16; 55 float buf[N]; 56 CvMat buffer = cvMat( 1, N, CV_32F, buf ); 57 CvSeqReader reader; 58 CvContour contour_header; 59 CvSeq* contour = 0; 60 CvSeqBlock block; 61 62 if( CV_IS_SEQ( array )) 63 { 64 contour = (CvSeq*)array; 65 if( !CV_IS_SEQ_POLYLINE( contour )) 66 CV_ERROR( CV_StsBadArg, "Unsupported sequence type" ); 67 if( is_closed < 0 ) 68 is_closed = CV_IS_SEQ_CLOSED( contour ); 69 } 70 else 71 { 72 is_closed = is_closed > 0; 73 CV_CALL( contour = cvPointSeqFromMat( 74 CV_SEQ_KIND_CURVE | (is_closed ? CV_SEQ_FLAG_CLOSED : 0), 75 array, &contour_header, &block )); 76 } 77 78 if( contour->total > 1 ) 79 { 80 int is_float = CV_SEQ_ELTYPE( contour ) == CV_32FC2; 81 82 cvStartReadSeq( contour, &reader, 0 ); 83 cvSetSeqReaderPos( &reader, slice.start_index ); 84 count = cvSliceLength( slice, contour ); 85 86 count -= !is_closed && count == contour->total; 87 88 /* scroll the reader by 1 point */ 89 reader.prev_elem = reader.ptr; 90 CV_NEXT_SEQ_ELEM( sizeof(CvPoint), reader ); 91 92 for( i = 0; i < count; i++ ) 93 { 94 float dx, dy; 95 96 if( !is_float ) 97 { 98 CvPoint* pt = (CvPoint*)reader.ptr; 99 CvPoint* prev_pt = (CvPoint*)reader.prev_elem; 100 101 dx = (float)pt->x - (float)prev_pt->x; 102 dy = (float)pt->y - (float)prev_pt->y; 103 } 104 else 105 { 106 CvPoint2D32f* pt = (CvPoint2D32f*)reader.ptr; 107 CvPoint2D32f* prev_pt = (CvPoint2D32f*)reader.prev_elem; 108 109 dx = pt->x - prev_pt->x; 110 dy = pt->y - prev_pt->y; 111 } 112 113 reader.prev_elem = reader.ptr; 114 CV_NEXT_SEQ_ELEM( contour->elem_size, reader ); 115 116 buffer.data.fl[j] = dx * dx + dy * dy; 117 if( ++j == N || i == count - 1 ) 118 { 119 buffer.cols = j; 120 cvPow( &buffer, &buffer, 0.5 ); 121 for( ; j > 0; j-- ) 122 perimeter += buffer.data.fl[j-1]; 123 } 124 } 125 } 126 127 __END__; 128 129 return perimeter; 130 } 131 132 133 static CvStatus 134 icvFindCircle( CvPoint2D32f pt0, CvPoint2D32f pt1, 135 CvPoint2D32f pt2, CvPoint2D32f * center, float *radius ) 136 { 137 double x1 = (pt0.x + pt1.x) * 0.5; 138 double dy1 = pt0.x - pt1.x; 139 double x2 = (pt1.x + pt2.x) * 0.5; 140 double dy2 = pt1.x - pt2.x; 141 double y1 = (pt0.y + pt1.y) * 0.5; 142 double dx1 = pt1.y - pt0.y; 143 double y2 = (pt1.y + pt2.y) * 0.5; 144 double dx2 = pt2.y - pt1.y; 145 double t = 0; 146 147 CvStatus result = CV_OK; 148 149 if( icvIntersectLines( x1, dx1, y1, dy1, x2, dx2, y2, dy2, &t ) >= 0 ) 150 { 151 center->x = (float) (x2 + dx2 * t); 152 center->y = (float) (y2 + dy2 * t); 153 *radius = (float) icvDistanceL2_32f( *center, pt0 ); 154 } 155 else 156 { 157 center->x = center->y = 0.f; 158 radius = 0; 159 result = CV_NOTDEFINED_ERR; 160 } 161 162 return result; 163 } 164 165 166 CV_INLINE double icvIsPtInCircle( CvPoint2D32f pt, CvPoint2D32f center, float radius ) 167 { 168 double dx = pt.x - center.x; 169 double dy = pt.y - center.y; 170 return (double)radius*radius - dx*dx - dy*dy; 171 } 172 173 174 static int 175 icvFindEnslosingCicle4pts_32f( CvPoint2D32f * pts, CvPoint2D32f * _center, float *_radius ) 176 { 177 int shuffles[4][4] = { {0, 1, 2, 3}, {0, 1, 3, 2}, {2, 3, 0, 1}, {2, 3, 1, 0} }; 178 179 int idxs[4] = { 0, 1, 2, 3 }; 180 int i, j, k = 1, mi = 0; 181 float max_dist = 0; 182 CvPoint2D32f center; 183 CvPoint2D32f min_center; 184 float radius, min_radius = FLT_MAX; 185 CvPoint2D32f res_pts[4]; 186 187 center = min_center = pts[0]; 188 radius = 1.f; 189 190 for( i = 0; i < 4; i++ ) 191 for( j = i + 1; j < 4; j++ ) 192 { 193 float dist = icvDistanceL2_32f( pts[i], pts[j] ); 194 195 if( max_dist < dist ) 196 { 197 max_dist = dist; 198 idxs[0] = i; 199 idxs[1] = j; 200 } 201 } 202 203 if( max_dist == 0 ) 204 goto function_exit; 205 206 k = 2; 207 for( i = 0; i < 4; i++ ) 208 { 209 for( j = 0; j < k; j++ ) 210 if( i == idxs[j] ) 211 break; 212 if( j == k ) 213 idxs[k++] = i; 214 } 215 216 center = cvPoint2D32f( (pts[idxs[0]].x + pts[idxs[1]].x)*0.5f, 217 (pts[idxs[0]].y + pts[idxs[1]].y)*0.5f ); 218 radius = (float)(icvDistanceL2_32f( pts[idxs[0]], center )*1.03); 219 if( radius < 1.f ) 220 radius = 1.f; 221 222 if( icvIsPtInCircle( pts[idxs[2]], center, radius ) >= 0 && 223 icvIsPtInCircle( pts[idxs[3]], center, radius ) >= 0 ) 224 { 225 k = 2; //rand()%2+2; 226 } 227 else 228 { 229 mi = -1; 230 for( i = 0; i < 4; i++ ) 231 { 232 if( icvFindCircle( pts[shuffles[i][0]], pts[shuffles[i][1]], 233 pts[shuffles[i][2]], ¢er, &radius ) >= 0 ) 234 { 235 radius *= 1.03f; 236 if( radius < 2.f ) 237 radius = 2.f; 238 239 if( icvIsPtInCircle( pts[shuffles[i][3]], center, radius ) >= 0 && 240 min_radius > radius ) 241 { 242 min_radius = radius; 243 min_center = center; 244 mi = i; 245 } 246 } 247 } 248 assert( mi >= 0 ); 249 if( mi < 0 ) 250 mi = 0; 251 k = 3; 252 center = min_center; 253 radius = min_radius; 254 for( i = 0; i < 4; i++ ) 255 idxs[i] = shuffles[mi][i]; 256 } 257 258 function_exit: 259 260 *_center = center; 261 *_radius = radius; 262 263 /* reorder output points */ 264 for( i = 0; i < 4; i++ ) 265 res_pts[i] = pts[idxs[i]]; 266 267 for( i = 0; i < 4; i++ ) 268 { 269 pts[i] = res_pts[i]; 270 assert( icvIsPtInCircle( pts[i], center, radius ) >= 0 ); 271 } 272 273 return k; 274 } 275 276 277 CV_IMPL int 278 cvMinEnclosingCircle( const void* array, CvPoint2D32f * _center, float *_radius ) 279 { 280 const int max_iters = 100; 281 const float eps = FLT_EPSILON*2; 282 CvPoint2D32f center = { 0, 0 }; 283 float radius = 0; 284 int result = 0; 285 286 if( _center ) 287 _center->x = _center->y = 0.f; 288 if( _radius ) 289 *_radius = 0; 290 291 CV_FUNCNAME( "cvMinEnclosingCircle" ); 292 293 __BEGIN__; 294 295 CvSeqReader reader; 296 int i, k, count; 297 CvPoint2D32f pts[8]; 298 CvContour contour_header; 299 CvSeqBlock block; 300 CvSeq* sequence = 0; 301 int is_float; 302 303 if( !_center || !_radius ) 304 CV_ERROR( CV_StsNullPtr, "Null center or radius pointers" ); 305 306 if( CV_IS_SEQ(array) ) 307 { 308 sequence = (CvSeq*)array; 309 if( !CV_IS_SEQ_POINT_SET( sequence )) 310 CV_ERROR( CV_StsBadArg, "The passed sequence is not a valid contour" ); 311 } 312 else 313 { 314 CV_CALL( sequence = cvPointSeqFromMat( 315 CV_SEQ_KIND_GENERIC, array, &contour_header, &block )); 316 } 317 318 if( sequence->total <= 0 ) 319 CV_ERROR_FROM_STATUS( CV_BADSIZE_ERR ); 320 321 CV_CALL( cvStartReadSeq( sequence, &reader, 0 )); 322 323 count = sequence->total; 324 is_float = CV_SEQ_ELTYPE(sequence) == CV_32FC2; 325 326 if( !is_float ) 327 { 328 CvPoint *pt_left, *pt_right, *pt_top, *pt_bottom; 329 CvPoint pt; 330 pt_left = pt_right = pt_top = pt_bottom = (CvPoint *)(reader.ptr); 331 CV_READ_SEQ_ELEM( pt, reader ); 332 333 for( i = 1; i < count; i++ ) 334 { 335 CvPoint* pt_ptr = (CvPoint*)reader.ptr; 336 CV_READ_SEQ_ELEM( pt, reader ); 337 338 if( pt.x < pt_left->x ) 339 pt_left = pt_ptr; 340 if( pt.x > pt_right->x ) 341 pt_right = pt_ptr; 342 if( pt.y < pt_top->y ) 343 pt_top = pt_ptr; 344 if( pt.y > pt_bottom->y ) 345 pt_bottom = pt_ptr; 346 } 347 348 pts[0] = cvPointTo32f( *pt_left ); 349 pts[1] = cvPointTo32f( *pt_right ); 350 pts[2] = cvPointTo32f( *pt_top ); 351 pts[3] = cvPointTo32f( *pt_bottom ); 352 } 353 else 354 { 355 CvPoint2D32f *pt_left, *pt_right, *pt_top, *pt_bottom; 356 CvPoint2D32f pt; 357 pt_left = pt_right = pt_top = pt_bottom = (CvPoint2D32f *) (reader.ptr); 358 CV_READ_SEQ_ELEM( pt, reader ); 359 360 for( i = 1; i < count; i++ ) 361 { 362 CvPoint2D32f* pt_ptr = (CvPoint2D32f*)reader.ptr; 363 CV_READ_SEQ_ELEM( pt, reader ); 364 365 if( pt.x < pt_left->x ) 366 pt_left = pt_ptr; 367 if( pt.x > pt_right->x ) 368 pt_right = pt_ptr; 369 if( pt.y < pt_top->y ) 370 pt_top = pt_ptr; 371 if( pt.y > pt_bottom->y ) 372 pt_bottom = pt_ptr; 373 } 374 375 pts[0] = *pt_left; 376 pts[1] = *pt_right; 377 pts[2] = *pt_top; 378 pts[3] = *pt_bottom; 379 } 380 381 for( k = 0; k < max_iters; k++ ) 382 { 383 double min_delta = 0, delta; 384 CvPoint2D32f ptfl; 385 386 icvFindEnslosingCicle4pts_32f( pts, ¢er, &radius ); 387 cvStartReadSeq( sequence, &reader, 0 ); 388 389 for( i = 0; i < count; i++ ) 390 { 391 if( !is_float ) 392 { 393 ptfl.x = (float)((CvPoint*)reader.ptr)->x; 394 ptfl.y = (float)((CvPoint*)reader.ptr)->y; 395 } 396 else 397 { 398 ptfl = *(CvPoint2D32f*)reader.ptr; 399 } 400 CV_NEXT_SEQ_ELEM( sequence->elem_size, reader ); 401 402 delta = icvIsPtInCircle( ptfl, center, radius ); 403 if( delta < min_delta ) 404 { 405 min_delta = delta; 406 pts[3] = ptfl; 407 } 408 } 409 result = min_delta >= 0; 410 if( result ) 411 break; 412 } 413 414 if( !result ) 415 { 416 cvStartReadSeq( sequence, &reader, 0 ); 417 radius = 0.f; 418 419 for( i = 0; i < count; i++ ) 420 { 421 CvPoint2D32f ptfl; 422 float t, dx, dy; 423 424 if( !is_float ) 425 { 426 ptfl.x = (float)((CvPoint*)reader.ptr)->x; 427 ptfl.y = (float)((CvPoint*)reader.ptr)->y; 428 } 429 else 430 { 431 ptfl = *(CvPoint2D32f*)reader.ptr; 432 } 433 434 CV_NEXT_SEQ_ELEM( sequence->elem_size, reader ); 435 dx = center.x - ptfl.x; 436 dy = center.y - ptfl.y; 437 t = dx*dx + dy*dy; 438 radius = MAX(radius,t); 439 } 440 441 radius = (float)(sqrt(radius)*(1 + eps)); 442 result = 1; 443 } 444 445 __END__; 446 447 *_center = center; 448 *_radius = radius; 449 450 return result; 451 } 452 453 454 /* area of a whole sequence */ 455 static CvStatus 456 icvContourArea( const CvSeq* contour, double *area ) 457 { 458 if( contour->total ) 459 { 460 CvSeqReader reader; 461 int lpt = contour->total; 462 double a00 = 0, xi_1, yi_1; 463 int is_float = CV_SEQ_ELTYPE(contour) == CV_32FC2; 464 465 cvStartReadSeq( contour, &reader, 0 ); 466 467 if( !is_float ) 468 { 469 xi_1 = ((CvPoint*)(reader.ptr))->x; 470 yi_1 = ((CvPoint*)(reader.ptr))->y; 471 } 472 else 473 { 474 xi_1 = ((CvPoint2D32f*)(reader.ptr))->x; 475 yi_1 = ((CvPoint2D32f*)(reader.ptr))->y; 476 } 477 CV_NEXT_SEQ_ELEM( contour->elem_size, reader ); 478 479 while( lpt-- > 0 ) 480 { 481 double dxy, xi, yi; 482 483 if( !is_float ) 484 { 485 xi = ((CvPoint*)(reader.ptr))->x; 486 yi = ((CvPoint*)(reader.ptr))->y; 487 } 488 else 489 { 490 xi = ((CvPoint2D32f*)(reader.ptr))->x; 491 yi = ((CvPoint2D32f*)(reader.ptr))->y; 492 } 493 CV_NEXT_SEQ_ELEM( contour->elem_size, reader ); 494 495 dxy = xi_1 * yi - xi * yi_1; 496 a00 += dxy; 497 xi_1 = xi; 498 yi_1 = yi; 499 } 500 501 *area = a00 * 0.5; 502 } 503 else 504 *area = 0; 505 506 return CV_OK; 507 } 508 509 510 /****************************************************************************************\ 511 512 copy data from one buffer to other buffer 513 514 \****************************************************************************************/ 515 516 static CvStatus 517 icvMemCopy( double **buf1, double **buf2, double **buf3, int *b_max ) 518 { 519 int bb; 520 521 if( (*buf1 == NULL && *buf2 == NULL) || *buf3 == NULL ) 522 return CV_NULLPTR_ERR; 523 524 bb = *b_max; 525 if( *buf2 == NULL ) 526 { 527 *b_max = 2 * (*b_max); 528 *buf2 = (double *)cvAlloc( (*b_max) * sizeof( double )); 529 530 if( *buf2 == NULL ) 531 return CV_OUTOFMEM_ERR; 532 533 memcpy( *buf2, *buf3, bb * sizeof( double )); 534 535 *buf3 = *buf2; 536 cvFree( buf1 ); 537 *buf1 = NULL; 538 } 539 else 540 { 541 *b_max = 2 * (*b_max); 542 *buf1 = (double *) cvAlloc( (*b_max) * sizeof( double )); 543 544 if( *buf1 == NULL ) 545 return CV_OUTOFMEM_ERR; 546 547 memcpy( *buf1, *buf3, bb * sizeof( double )); 548 549 *buf3 = *buf1; 550 cvFree( buf2 ); 551 *buf2 = NULL; 552 } 553 return CV_OK; 554 } 555 556 557 /* area of a contour sector */ 558 static CvStatus icvContourSecArea( CvSeq * contour, CvSlice slice, double *area ) 559 { 560 CvPoint pt; /* pointer to points */ 561 CvPoint pt_s, pt_e; /* first and last points */ 562 CvSeqReader reader; /* points reader of contour */ 563 564 int p_max = 2, p_ind; 565 int lpt, flag, i; 566 double a00; /* unnormalized moments m00 */ 567 double xi, yi, xi_1, yi_1, x0, y0, dxy, sk, sk1, t; 568 double x_s, y_s, nx, ny, dx, dy, du, dv; 569 double eps = 1.e-5; 570 double *p_are1, *p_are2, *p_are; 571 572 assert( contour != NULL ); 573 574 if( contour == NULL ) 575 return CV_NULLPTR_ERR; 576 577 if( !CV_IS_SEQ_POLYGON( contour )) 578 return CV_BADFLAG_ERR; 579 580 lpt = cvSliceLength( slice, contour ); 581 /*if( n2 >= n1 ) 582 lpt = n2 - n1 + 1; 583 else 584 lpt = contour->total - n1 + n2 + 1;*/ 585 586 if( contour->total && lpt > 2 ) 587 { 588 a00 = x0 = y0 = xi_1 = yi_1 = 0; 589 sk1 = 0; 590 flag = 0; 591 dxy = 0; 592 p_are1 = (double *) cvAlloc( p_max * sizeof( double )); 593 594 if( p_are1 == NULL ) 595 return CV_OUTOFMEM_ERR; 596 597 p_are = p_are1; 598 p_are2 = NULL; 599 600 cvStartReadSeq( contour, &reader, 0 ); 601 cvSetSeqReaderPos( &reader, slice.start_index ); 602 CV_READ_SEQ_ELEM( pt_s, reader ); 603 p_ind = 0; 604 cvSetSeqReaderPos( &reader, slice.end_index ); 605 CV_READ_SEQ_ELEM( pt_e, reader ); 606 607 /* normal coefficients */ 608 nx = pt_s.y - pt_e.y; 609 ny = pt_e.x - pt_s.x; 610 cvSetSeqReaderPos( &reader, slice.start_index ); 611 612 while( lpt-- > 0 ) 613 { 614 CV_READ_SEQ_ELEM( pt, reader ); 615 616 if( flag == 0 ) 617 { 618 xi_1 = (double) pt.x; 619 yi_1 = (double) pt.y; 620 x0 = xi_1; 621 y0 = yi_1; 622 sk1 = 0; 623 flag = 1; 624 } 625 else 626 { 627 xi = (double) pt.x; 628 yi = (double) pt.y; 629 630 /**************** edges intersection examination **************************/ 631 sk = nx * (xi - pt_s.x) + ny * (yi - pt_s.y); 632 if( (fabs( sk ) < eps && lpt > 0) || sk * sk1 < -eps ) 633 { 634 if( fabs( sk ) < eps ) 635 { 636 dxy = xi_1 * yi - xi * yi_1; 637 a00 = a00 + dxy; 638 dxy = xi * y0 - x0 * yi; 639 a00 = a00 + dxy; 640 641 if( p_ind >= p_max ) 642 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max ); 643 644 p_are[p_ind] = a00 / 2.; 645 p_ind++; 646 a00 = 0; 647 sk1 = 0; 648 x0 = xi; 649 y0 = yi; 650 dxy = 0; 651 } 652 else 653 { 654 /* define intersection point */ 655 dv = yi - yi_1; 656 du = xi - xi_1; 657 dx = ny; 658 dy = -nx; 659 if( fabs( du ) > eps ) 660 t = ((yi_1 - pt_s.y) * du + dv * (pt_s.x - xi_1)) / 661 (du * dy - dx * dv); 662 else 663 t = (xi_1 - pt_s.x) / dx; 664 if( t > eps && t < 1 - eps ) 665 { 666 x_s = pt_s.x + t * dx; 667 y_s = pt_s.y + t * dy; 668 dxy = xi_1 * y_s - x_s * yi_1; 669 a00 += dxy; 670 dxy = x_s * y0 - x0 * y_s; 671 a00 += dxy; 672 if( p_ind >= p_max ) 673 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max ); 674 675 p_are[p_ind] = a00 / 2.; 676 p_ind++; 677 678 a00 = 0; 679 sk1 = 0; 680 x0 = x_s; 681 y0 = y_s; 682 dxy = x_s * yi - xi * y_s; 683 } 684 } 685 } 686 else 687 dxy = xi_1 * yi - xi * yi_1; 688 689 a00 += dxy; 690 xi_1 = xi; 691 yi_1 = yi; 692 sk1 = sk; 693 694 } 695 } 696 697 xi = x0; 698 yi = y0; 699 dxy = xi_1 * yi - xi * yi_1; 700 701 a00 += dxy; 702 703 if( p_ind >= p_max ) 704 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max ); 705 706 p_are[p_ind] = a00 / 2.; 707 p_ind++; 708 709 /* common area calculation */ 710 *area = 0; 711 for( i = 0; i < p_ind; i++ ) 712 (*area) += fabs( p_are[i] ); 713 714 if( p_are1 != NULL ) 715 cvFree( &p_are1 ); 716 else if( p_are2 != NULL ) 717 cvFree( &p_are2 ); 718 719 return CV_OK; 720 } 721 else 722 return CV_BADSIZE_ERR; 723 } 724 725 726 /* external contour area function */ 727 CV_IMPL double 728 cvContourArea( const void *array, CvSlice slice ) 729 { 730 double area = 0; 731 732 CV_FUNCNAME( "cvContourArea" ); 733 734 __BEGIN__; 735 736 CvContour contour_header; 737 CvSeq* contour = 0; 738 CvSeqBlock block; 739 740 if( CV_IS_SEQ( array )) 741 { 742 contour = (CvSeq*)array; 743 if( !CV_IS_SEQ_POLYLINE( contour )) 744 CV_ERROR( CV_StsBadArg, "Unsupported sequence type" ); 745 } 746 else 747 { 748 CV_CALL( contour = cvPointSeqFromMat( 749 CV_SEQ_KIND_CURVE, array, &contour_header, &block )); 750 } 751 752 if( cvSliceLength( slice, contour ) == contour->total ) 753 { 754 IPPI_CALL( icvContourArea( contour, &area )); 755 } 756 else 757 { 758 if( CV_SEQ_ELTYPE( contour ) != CV_32SC2 ) 759 CV_ERROR( CV_StsUnsupportedFormat, 760 "Only curves with integer coordinates are supported in case of contour slice" ); 761 IPPI_CALL( icvContourSecArea( contour, slice, &area )); 762 } 763 764 __END__; 765 766 return area; 767 } 768 769 770 /* for now this function works bad with singular cases 771 You can see in the code, that when some troubles with 772 matrices or some variables occur - 773 box filled with zero values is returned. 774 However in general function works fine. 775 */ 776 static void 777 icvFitEllipse_F( CvSeq* points, CvBox2D* box ) 778 { 779 CvMat* D = 0; 780 781 CV_FUNCNAME( "icvFitEllipse_F" ); 782 783 __BEGIN__; 784 785 double S[36], C[36], T[36]; 786 787 int i, j; 788 double eigenvalues[6], eigenvectors[36]; 789 double a, b, c, d, e, f; 790 double x0, y0, idet, scale, offx = 0, offy = 0; 791 792 int n = points->total; 793 CvSeqReader reader; 794 int is_float = CV_SEQ_ELTYPE(points) == CV_32FC2; 795 796 CvMat _S = cvMat(6,6,CV_64F,S), _C = cvMat(6,6,CV_64F,C), _T = cvMat(6,6,CV_64F,T); 797 CvMat _EIGVECS = cvMat(6,6,CV_64F,eigenvectors), _EIGVALS = cvMat(6,1,CV_64F,eigenvalues); 798 799 /* create matrix D of input points */ 800 CV_CALL( D = cvCreateMat( n, 6, CV_64F )); 801 802 cvStartReadSeq( points, &reader ); 803 804 /* shift all points to zero */ 805 for( i = 0; i < n; i++ ) 806 { 807 if( !is_float ) 808 { 809 offx += ((CvPoint*)reader.ptr)->x; 810 offy += ((CvPoint*)reader.ptr)->y; 811 } 812 else 813 { 814 offx += ((CvPoint2D32f*)reader.ptr)->x; 815 offy += ((CvPoint2D32f*)reader.ptr)->y; 816 } 817 CV_NEXT_SEQ_ELEM( points->elem_size, reader ); 818 } 819 820 offx /= n; 821 offy /= n; 822 823 // fill matrix rows as (x*x, x*y, y*y, x, y, 1 ) 824 for( i = 0; i < n; i++ ) 825 { 826 double x, y; 827 double* Dptr = D->data.db + i*6; 828 829 if( !is_float ) 830 { 831 x = ((CvPoint*)reader.ptr)->x - offx; 832 y = ((CvPoint*)reader.ptr)->y - offy; 833 } 834 else 835 { 836 x = ((CvPoint2D32f*)reader.ptr)->x - offx; 837 y = ((CvPoint2D32f*)reader.ptr)->y - offy; 838 } 839 CV_NEXT_SEQ_ELEM( points->elem_size, reader ); 840 841 Dptr[0] = x * x; 842 Dptr[1] = x * y; 843 Dptr[2] = y * y; 844 Dptr[3] = x; 845 Dptr[4] = y; 846 Dptr[5] = 1.; 847 } 848 849 // S = D^t*D 850 cvMulTransposed( D, &_S, 1 ); 851 cvSVD( &_S, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ); 852 853 for( i = 0; i < 6; i++ ) 854 { 855 double a = eigenvalues[i]; 856 a = a < DBL_EPSILON ? 0 : 1./sqrt(sqrt(a)); 857 for( j = 0; j < 6; j++ ) 858 eigenvectors[i*6 + j] *= a; 859 } 860 861 // C = Q^-1 = transp(INVEIGV) * INVEIGV 862 cvMulTransposed( &_EIGVECS, &_C, 1 ); 863 864 cvZero( &_S ); 865 S[2] = 2.; 866 S[7] = -1.; 867 S[12] = 2.; 868 869 // S = Q^-1*S*Q^-1 870 cvMatMul( &_C, &_S, &_T ); 871 cvMatMul( &_T, &_C, &_S ); 872 873 // and find its eigenvalues and vectors too 874 //cvSVD( &_S, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ); 875 cvEigenVV( &_S, &_EIGVECS, &_EIGVALS, 0 ); 876 877 for( i = 0; i < 3; i++ ) 878 if( eigenvalues[i] > 0 ) 879 break; 880 881 if( i >= 3 /*eigenvalues[0] < DBL_EPSILON*/ ) 882 { 883 box->center.x = box->center.y = 884 box->size.width = box->size.height = 885 box->angle = 0.f; 886 EXIT; 887 } 888 889 // now find truthful eigenvector 890 _EIGVECS = cvMat( 6, 1, CV_64F, eigenvectors + 6*i ); 891 _T = cvMat( 6, 1, CV_64F, T ); 892 // Q^-1*eigenvecs[0] 893 cvMatMul( &_C, &_EIGVECS, &_T ); 894 895 // extract vector components 896 a = T[0]; b = T[1]; c = T[2]; d = T[3]; e = T[4]; f = T[5]; 897 898 ///////////////// extract ellipse axes from above values //////////////// 899 900 /* 901 1) find center of ellipse 902 it satisfy equation 903 | a b/2 | * | x0 | + | d/2 | = |0 | 904 | b/2 c | | y0 | | e/2 | |0 | 905 906 */ 907 idet = a * c - b * b * 0.25; 908 idet = idet > DBL_EPSILON ? 1./idet : 0; 909 910 // we must normalize (a b c d e f ) to fit (4ac-b^2=1) 911 scale = sqrt( 0.25 * idet ); 912 913 if( scale < DBL_EPSILON ) 914 { 915 box->center.x = (float)offx; 916 box->center.y = (float)offy; 917 box->size.width = box->size.height = box->angle = 0.f; 918 EXIT; 919 } 920 921 a *= scale; 922 b *= scale; 923 c *= scale; 924 d *= scale; 925 e *= scale; 926 f *= scale; 927 928 x0 = (-d * c + e * b * 0.5) * 2.; 929 y0 = (-a * e + d * b * 0.5) * 2.; 930 931 // recover center 932 box->center.x = (float)(x0 + offx); 933 box->center.y = (float)(y0 + offy); 934 935 // offset ellipse to (x0,y0) 936 // new f == F(x0,y0) 937 f += a * x0 * x0 + b * x0 * y0 + c * y0 * y0 + d * x0 + e * y0; 938 939 if( fabs(f) < DBL_EPSILON ) 940 { 941 box->size.width = box->size.height = box->angle = 0.f; 942 EXIT; 943 } 944 945 scale = -1. / f; 946 // normalize to f = 1 947 a *= scale; 948 b *= scale; 949 c *= scale; 950 951 // extract axis of ellipse 952 // one more eigenvalue operation 953 S[0] = a; 954 S[1] = S[2] = b * 0.5; 955 S[3] = c; 956 957 _S = cvMat( 2, 2, CV_64F, S ); 958 _EIGVECS = cvMat( 2, 2, CV_64F, eigenvectors ); 959 _EIGVALS = cvMat( 1, 2, CV_64F, eigenvalues ); 960 cvSVD( &_S, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ); 961 962 // exteract axis length from eigenvectors 963 box->size.width = (float)(2./sqrt(eigenvalues[0])); 964 box->size.height = (float)(2./sqrt(eigenvalues[1])); 965 966 // calc angle 967 box->angle = (float)(180 - atan2(eigenvectors[2], eigenvectors[3])*180/CV_PI); 968 969 __END__; 970 971 cvReleaseMat( &D ); 972 } 973 974 975 CV_IMPL CvBox2D 976 cvFitEllipse2( const CvArr* array ) 977 { 978 CvBox2D box; 979 double* Ad = 0, *bd = 0; 980 981 CV_FUNCNAME( "cvFitEllipse2" ); 982 983 memset( &box, 0, sizeof(box)); 984 985 __BEGIN__; 986 987 CvContour contour_header; 988 CvSeq* ptseq = 0; 989 CvSeqBlock block; 990 int n; 991 992 if( CV_IS_SEQ( array )) 993 { 994 ptseq = (CvSeq*)array; 995 if( !CV_IS_SEQ_POINT_SET( ptseq )) 996 CV_ERROR( CV_StsBadArg, "Unsupported sequence type" ); 997 } 998 else 999 { 1000 CV_CALL( ptseq = cvPointSeqFromMat( 1001 CV_SEQ_KIND_GENERIC, array, &contour_header, &block )); 1002 } 1003 1004 n = ptseq->total; 1005 if( n < 5 ) 1006 CV_ERROR( CV_StsBadSize, "Number of points should be >= 6" ); 1007 #if 1 1008 icvFitEllipse_F( ptseq, &box ); 1009 #else 1010 /* 1011 * New fitellipse algorithm, contributed by Dr. Daniel Weiss 1012 */ 1013 { 1014 double gfp[5], rp[5], t; 1015 CvMat A, b, x; 1016 const double min_eps = 1e-6; 1017 int i, is_float; 1018 CvSeqReader reader; 1019 1020 CV_CALL( Ad = (double*)cvAlloc( n*5*sizeof(Ad[0]) )); 1021 CV_CALL( bd = (double*)cvAlloc( n*sizeof(bd[0]) )); 1022 1023 // first fit for parameters A - E 1024 A = cvMat( n, 5, CV_64F, Ad ); 1025 b = cvMat( n, 1, CV_64F, bd ); 1026 x = cvMat( 5, 1, CV_64F, gfp ); 1027 1028 cvStartReadSeq( ptseq, &reader ); 1029 is_float = CV_SEQ_ELTYPE(ptseq) == CV_32FC2; 1030 1031 for( i = 0; i < n; i++ ) 1032 { 1033 CvPoint2D32f p; 1034 if( is_float ) 1035 p = *(CvPoint2D32f*)(reader.ptr); 1036 else 1037 { 1038 p.x = (float)((int*)reader.ptr)[0]; 1039 p.y = (float)((int*)reader.ptr)[1]; 1040 } 1041 CV_NEXT_SEQ_ELEM( sizeof(p), reader ); 1042 1043 bd[i] = 10000.0; // 1.0? 1044 Ad[i*5] = -(double)p.x * p.x; // A - C signs inverted as proposed by APP 1045 Ad[i*5 + 1] = -(double)p.y * p.y; 1046 Ad[i*5 + 2] = -(double)p.x * p.y; 1047 Ad[i*5 + 3] = p.x; 1048 Ad[i*5 + 4] = p.y; 1049 } 1050 1051 cvSolve( &A, &b, &x, CV_SVD ); 1052 1053 // now use general-form parameters A - E to find the ellipse center: 1054 // differentiate general form wrt x/y to get two equations for cx and cy 1055 A = cvMat( 2, 2, CV_64F, Ad ); 1056 b = cvMat( 2, 1, CV_64F, bd ); 1057 x = cvMat( 2, 1, CV_64F, rp ); 1058 Ad[0] = 2 * gfp[0]; 1059 Ad[1] = Ad[2] = gfp[2]; 1060 Ad[3] = 2 * gfp[1]; 1061 bd[0] = gfp[3]; 1062 bd[1] = gfp[4]; 1063 cvSolve( &A, &b, &x, CV_SVD ); 1064 1065 // re-fit for parameters A - C with those center coordinates 1066 A = cvMat( n, 3, CV_64F, Ad ); 1067 b = cvMat( n, 1, CV_64F, bd ); 1068 x = cvMat( 3, 1, CV_64F, gfp ); 1069 for( i = 0; i < n; i++ ) 1070 { 1071 CvPoint2D32f p; 1072 if( is_float ) 1073 p = *(CvPoint2D32f*)(reader.ptr); 1074 else 1075 { 1076 p.x = (float)((int*)reader.ptr)[0]; 1077 p.y = (float)((int*)reader.ptr)[1]; 1078 } 1079 CV_NEXT_SEQ_ELEM( sizeof(p), reader ); 1080 bd[i] = 1.0; 1081 Ad[i * 3] = (p.x - rp[0]) * (p.x - rp[0]); 1082 Ad[i * 3 + 1] = (p.y - rp[1]) * (p.y - rp[1]); 1083 Ad[i * 3 + 2] = (p.x - rp[0]) * (p.y - rp[1]); 1084 } 1085 cvSolve(&A, &b, &x, CV_SVD); 1086 1087 // store angle and radii 1088 rp[4] = -0.5 * atan2(gfp[2], gfp[1] - gfp[0]); // convert from APP angle usage 1089 t = sin(-2.0 * rp[4]); 1090 if( fabs(t) > fabs(gfp[2])*min_eps ) 1091 t = gfp[2]/t; 1092 else 1093 t = gfp[1] - gfp[0]; 1094 rp[2] = fabs(gfp[0] + gfp[1] - t); 1095 if( rp[2] > min_eps ) 1096 rp[2] = sqrt(2.0 / rp[2]); 1097 rp[3] = fabs(gfp[0] + gfp[1] + t); 1098 if( rp[3] > min_eps ) 1099 rp[3] = sqrt(2.0 / rp[3]); 1100 1101 box.center.x = (float)rp[0]; 1102 box.center.y = (float)rp[1]; 1103 box.size.width = (float)(rp[2]*2); 1104 box.size.height = (float)(rp[3]*2); 1105 if( box.size.width > box.size.height ) 1106 { 1107 float tmp; 1108 CV_SWAP( box.size.width, box.size.height, tmp ); 1109 box.angle = (float)(90 + rp[4]*180/CV_PI); 1110 } 1111 if( box.angle < -180 ) 1112 box.angle += 360; 1113 if( box.angle > 360 ) 1114 box.angle -= 360; 1115 } 1116 #endif 1117 __END__; 1118 1119 cvFree( &Ad ); 1120 cvFree( &bd ); 1121 1122 return box; 1123 } 1124 1125 1126 /* Calculates bounding rectagnle of a point set or retrieves already calculated */ 1127 CV_IMPL CvRect 1128 cvBoundingRect( CvArr* array, int update ) 1129 { 1130 CvSeqReader reader; 1131 CvRect rect = { 0, 0, 0, 0 }; 1132 CvContour contour_header; 1133 CvSeq* ptseq = 0; 1134 CvSeqBlock block; 1135 1136 CV_FUNCNAME( "cvBoundingRect" ); 1137 1138 __BEGIN__; 1139 1140 CvMat stub, *mat = 0; 1141 int xmin = 0, ymin = 0, xmax = -1, ymax = -1, i, j, k; 1142 int calculate = update; 1143 1144 if( CV_IS_SEQ( array )) 1145 { 1146 ptseq = (CvSeq*)array; 1147 if( !CV_IS_SEQ_POINT_SET( ptseq )) 1148 CV_ERROR( CV_StsBadArg, "Unsupported sequence type" ); 1149 1150 if( ptseq->header_size < (int)sizeof(CvContour)) 1151 { 1152 /*if( update == 1 ) 1153 CV_ERROR( CV_StsBadArg, "The header is too small to fit the rectangle, " 1154 "so it could not be updated" );*/ 1155 update = 0; 1156 calculate = 1; 1157 } 1158 } 1159 else 1160 { 1161 CV_CALL( mat = cvGetMat( array, &stub )); 1162 if( CV_MAT_TYPE(mat->type) == CV_32SC2 || 1163 CV_MAT_TYPE(mat->type) == CV_32FC2 ) 1164 { 1165 CV_CALL( ptseq = cvPointSeqFromMat( 1166 CV_SEQ_KIND_GENERIC, mat, &contour_header, &block )); 1167 mat = 0; 1168 } 1169 else if( CV_MAT_TYPE(mat->type) != CV_8UC1 && 1170 CV_MAT_TYPE(mat->type) != CV_8SC1 ) 1171 CV_ERROR( CV_StsUnsupportedFormat, 1172 "The image/matrix format is not supported by the function" ); 1173 update = 0; 1174 calculate = 1; 1175 } 1176 1177 if( !calculate ) 1178 { 1179 rect = ((CvContour*)ptseq)->rect; 1180 EXIT; 1181 } 1182 1183 if( mat ) 1184 { 1185 CvSize size = cvGetMatSize(mat); 1186 xmin = size.width; 1187 ymin = -1; 1188 1189 for( i = 0; i < size.height; i++ ) 1190 { 1191 uchar* _ptr = mat->data.ptr + i*mat->step; 1192 uchar* ptr = (uchar*)cvAlignPtr(_ptr, 4); 1193 int have_nz = 0, k_min, offset = (int)(ptr - _ptr); 1194 j = 0; 1195 offset = MIN(offset, size.width); 1196 for( ; j < offset; j++ ) 1197 if( _ptr[j] ) 1198 { 1199 have_nz = 1; 1200 break; 1201 } 1202 if( j < offset ) 1203 { 1204 if( j < xmin ) 1205 xmin = j; 1206 if( j > xmax ) 1207 xmax = j; 1208 } 1209 if( offset < size.width ) 1210 { 1211 xmin -= offset; 1212 xmax -= offset; 1213 size.width -= offset; 1214 j = 0; 1215 for( ; j <= xmin - 4; j += 4 ) 1216 if( *((int*)(ptr+j)) ) 1217 break; 1218 for( ; j < xmin; j++ ) 1219 if( ptr[j] ) 1220 { 1221 xmin = j; 1222 if( j > xmax ) 1223 xmax = j; 1224 have_nz = 1; 1225 break; 1226 } 1227 k_min = MAX(j-1, xmax); 1228 k = size.width - 1; 1229 for( ; k > k_min && (k&3) != 3; k-- ) 1230 if( ptr[k] ) 1231 break; 1232 if( k > k_min && (k&3) == 3 ) 1233 { 1234 for( ; k > k_min+3; k -= 4 ) 1235 if( *((int*)(ptr+k-3)) ) 1236 break; 1237 } 1238 for( ; k > k_min; k-- ) 1239 if( ptr[k] ) 1240 { 1241 xmax = k; 1242 have_nz = 1; 1243 break; 1244 } 1245 if( !have_nz ) 1246 { 1247 j &= ~3; 1248 for( ; j <= k - 3; j += 4 ) 1249 if( *((int*)(ptr+j)) ) 1250 break; 1251 for( ; j <= k; j++ ) 1252 if( ptr[j] ) 1253 { 1254 have_nz = 1; 1255 break; 1256 } 1257 } 1258 xmin += offset; 1259 xmax += offset; 1260 size.width += offset; 1261 } 1262 if( have_nz ) 1263 { 1264 if( ymin < 0 ) 1265 ymin = i; 1266 ymax = i; 1267 } 1268 } 1269 1270 if( xmin >= size.width ) 1271 xmin = ymin = 0; 1272 } 1273 else if( ptseq->total ) 1274 { 1275 int is_float = CV_SEQ_ELTYPE(ptseq) == CV_32FC2; 1276 cvStartReadSeq( ptseq, &reader, 0 ); 1277 1278 if( !is_float ) 1279 { 1280 CvPoint pt; 1281 /* init values */ 1282 CV_READ_SEQ_ELEM( pt, reader ); 1283 xmin = xmax = pt.x; 1284 ymin = ymax = pt.y; 1285 1286 for( i = 1; i < ptseq->total; i++ ) 1287 { 1288 CV_READ_SEQ_ELEM( pt, reader ); 1289 1290 if( xmin > pt.x ) 1291 xmin = pt.x; 1292 1293 if( xmax < pt.x ) 1294 xmax = pt.x; 1295 1296 if( ymin > pt.y ) 1297 ymin = pt.y; 1298 1299 if( ymax < pt.y ) 1300 ymax = pt.y; 1301 } 1302 } 1303 else 1304 { 1305 CvPoint pt; 1306 Cv32suf v; 1307 /* init values */ 1308 CV_READ_SEQ_ELEM( pt, reader ); 1309 xmin = xmax = CV_TOGGLE_FLT(pt.x); 1310 ymin = ymax = CV_TOGGLE_FLT(pt.y); 1311 1312 for( i = 1; i < ptseq->total; i++ ) 1313 { 1314 CV_READ_SEQ_ELEM( pt, reader ); 1315 pt.x = CV_TOGGLE_FLT(pt.x); 1316 pt.y = CV_TOGGLE_FLT(pt.y); 1317 1318 if( xmin > pt.x ) 1319 xmin = pt.x; 1320 1321 if( xmax < pt.x ) 1322 xmax = pt.x; 1323 1324 if( ymin > pt.y ) 1325 ymin = pt.y; 1326 1327 if( ymax < pt.y ) 1328 ymax = pt.y; 1329 } 1330 1331 v.i = CV_TOGGLE_FLT(xmin); xmin = cvFloor(v.f); 1332 v.i = CV_TOGGLE_FLT(ymin); ymin = cvFloor(v.f); 1333 /* because right and bottom sides of 1334 the bounding rectangle are not inclusive 1335 (note +1 in width and height calculation below), 1336 cvFloor is used here instead of cvCeil */ 1337 v.i = CV_TOGGLE_FLT(xmax); xmax = cvFloor(v.f); 1338 v.i = CV_TOGGLE_FLT(ymax); ymax = cvFloor(v.f); 1339 } 1340 } 1341 1342 rect.x = xmin; 1343 rect.y = ymin; 1344 rect.width = xmax - xmin + 1; 1345 rect.height = ymax - ymin + 1; 1346 1347 if( update ) 1348 ((CvContour*)ptseq)->rect = rect; 1349 1350 __END__; 1351 1352 return rect; 1353 } 1354 1355 1356 /* End of file. */ 1357
