git.maemo.org Git - opencv/blob - cv/src/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     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]], &center, &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, &center, &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. */