opencv/src/cv/cvshapedescr.cpp

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