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