git.maemo.org Git - opencv/blob - cv/src/cvapprox.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 /****************************************************************************************\
  44 *                                  Chain Approximation                                   *
  45 \****************************************************************************************/
  46
  47 typedef struct _CvPtInfo
  48 {
  49     CvPoint pt;
  50     int k;                      /* support region */
  51     int s;                      /* curvature value */
  52     struct _CvPtInfo *next;
  53 }
  54 _CvPtInfo;
  55
  56
  57 /* curvature: 0 - 1-curvature, 1 - k-cosine curvature. */
  58 CvStatus
  59 icvApproximateChainTC89( CvChain*               chain,
  60                          int                    header_size,
  61                          CvMemStorage*          storage,
  62                          CvSeq**                contour,
  63                          int    method )
  64 {
  65     static const int abs_diff[] = { 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1 };
  66
  67     char            local_buffer[1 << 16];
  68     char*           buffer = local_buffer;
  69     int             buffer_size;
  70
  71     _CvPtInfo       temp;
  72     _CvPtInfo       *array, *first = 0, *current = 0, *prev_current = 0;
  73     int             i, j, i1, i2, s, len;
  74     int             count;
  75
  76     CvChainPtReader reader;
  77     CvSeqWriter     writer;
  78     CvPoint         pt = chain->origin;
  79
  80     assert( chain && contour && buffer );
  81
  82     buffer_size = (chain->total + 8) * sizeof( _CvPtInfo );
  83
  84     *contour = 0;
  85
  86     if( !CV_IS_SEQ_CHAIN_CONTOUR( chain ))
  87         return CV_BADFLAG_ERR;
  88
  89     if( header_size < (int)sizeof(CvContour) )
  90         return CV_BADSIZE_ERR;
  91
  92     cvStartWriteSeq( (chain->flags & ~CV_SEQ_ELTYPE_MASK) | CV_SEQ_ELTYPE_POINT,
  93                      header_size, sizeof( CvPoint ), storage, &writer );
  94
  95     if( chain->total == 0 )
  96     {
  97         CV_WRITE_SEQ_ELEM( pt, writer );
  98         goto exit_function;
  99     }
 100
 101     cvStartReadChainPoints( chain, &reader );
 102
 103     if( method > CV_CHAIN_APPROX_SIMPLE && buffer_size > (int)sizeof(local_buffer))
 104     {
 105         buffer = (char *) cvAlloc( buffer_size );
 106         if( !buffer )
 107             return CV_OUTOFMEM_ERR;
 108     }
 109
 110     array = (_CvPtInfo *) buffer;
 111     count = chain->total;
 112
 113     temp.next = 0;
 114     current = &temp;
 115
 116     /* Pass 0.
 117        Restores all the digital curve points from the chain code.
 118        Removes the points (from the resultant polygon)
 119        that have zero 1-curvature */
 120     for( i = 0; i < count; i++ )
 121     {
 122         int prev_code = *reader.prev_elem;
 123
 124         reader.prev_elem = reader.ptr;
 125         CV_READ_CHAIN_POINT( pt, reader );
 126
 127         /* calc 1-curvature */
 128         s = abs_diff[reader.code - prev_code + 7];
 129
 130         if( method <= CV_CHAIN_APPROX_SIMPLE )
 131         {
 132             if( method == CV_CHAIN_APPROX_NONE || s != 0 )
 133             {
 134                 CV_WRITE_SEQ_ELEM( pt, writer );
 135             }
 136         }
 137         else
 138         {
 139             if( s != 0 )
 140                 current = current->next = array + i;
 141             array[i].s = s;
 142             array[i].pt = pt;
 143         }
 144     }
 145
 146     //assert( pt.x == chain->origin.x && pt.y == chain->origin.y );
 147
 148     if( method <= CV_CHAIN_APPROX_SIMPLE )
 149         goto exit_function;
 150
 151     current->next = 0;
 152
 153     len = i;
 154     current = temp.next;
 155
 156     assert( current );
 157
 158     /* Pass 1.
 159        Determines support region for all the remained points */
 160     do
 161     {
 162         CvPoint pt0;
 163         int k, l = 0, d_num = 0;
 164
 165         i = (int)(current - array);
 166         pt0 = array[i].pt;
 167
 168         /* determine support region */
 169         for( k = 1;; k++ )
 170         {
 171             int lk, dk_num;
 172             int dx, dy;
 173             Cv32suf d;
 174
 175             assert( k <= len );
 176
 177             /* calc indices */
 178             i1 = i - k;
 179             i1 += i1 < 0 ? len : 0;
 180             i2 = i + k;
 181             i2 -= i2 >= len ? len : 0;
 182
 183             dx = array[i2].pt.x - array[i1].pt.x;
 184             dy = array[i2].pt.y - array[i1].pt.y;
 185
 186             /* distance between p_(i - k) and p_(i + k) */
 187             lk = dx * dx + dy * dy;
 188
 189             /* distance between p_i and the line (p_(i-k), p_(i+k)) */
 190             dk_num = (pt0.x - array[i1].pt.x) * dy - (pt0.y - array[i1].pt.y) * dx;
 191             d.f = (float) (((double) d_num) * lk - ((double) dk_num) * l);
 192
 193             if( k > 1 && (l >= lk || (d_num > 0 && d.i <= 0 || d_num < 0 && d.i >= 0)))
 194                 break;
 195
 196             d_num = dk_num;
 197             l = lk;
 198         }
 199
 200         current->k = --k;
 201
 202         /* determine cosine curvature if it should be used */
 203         if( method == CV_CHAIN_APPROX_TC89_KCOS )
 204         {
 205             /* calc k-cosine curvature */
 206             for( j = k, s = 0; j > 0; j-- )
 207             {
 208                 double temp_num;
 209                 int dx1, dy1, dx2, dy2;
 210                 Cv32suf sk;
 211
 212                 i1 = i - j;
 213                 i1 += i1 < 0 ? len : 0;
 214                 i2 = i + j;
 215                 i2 -= i2 >= len ? len : 0;
 216
 217                 dx1 = array[i1].pt.x - pt0.x;
 218                 dy1 = array[i1].pt.y - pt0.y;
 219                 dx2 = array[i2].pt.x - pt0.x;
 220                 dy2 = array[i2].pt.y - pt0.y;
 221
 222                 if( (dx1 | dy1) == 0 || (dx2 | dy2) == 0 )
 223                     break;
 224
 225                 temp_num = dx1 * dx2 + dy1 * dy2;
 226                 temp_num =
 227                     (float) (temp_num /
 228                              sqrt( ((double)dx1 * dx1 + (double)dy1 * dy1) *
 229                                    ((double)dx2 * dx2 + (double)dy2 * dy2) ));
 230                 sk.f = (float) (temp_num + 1.1);
 231
 232                 assert( 0 <= sk.f && sk.f <= 2.2 );
 233                 if( j < k && sk.i <= s )
 234                     break;
 235
 236                 s = sk.i;
 237             }
 238             current->s = s;
 239         }
 240         current = current->next;
 241     }
 242     while( current != 0 );
 243
 244     prev_current = &temp;
 245     current = temp.next;
 246
 247     /* Pass 2.
 248        Performs non-maxima supression */
 249     do
 250     {
 251         int k2 = current->k >> 1;
 252
 253         s = current->s;
 254         i = (int)(current - array);
 255
 256         for( j = 1; j <= k2; j++ )
 257         {
 258             i2 = i - j;
 259             i2 += i2 < 0 ? len : 0;
 260
 261             if( array[i2].s > s )
 262                 break;
 263
 264             i2 = i + j;
 265             i2 -= i2 >= len ? len : 0;
 266
 267             if( array[i2].s > s )
 268                 break;
 269         }
 270
 271         if( j <= k2 )           /* exclude point */
 272         {
 273             prev_current->next = current->next;
 274             current->s = 0;     /* "clear" point */
 275         }
 276         else
 277             prev_current = current;
 278         current = current->next;
 279     }
 280     while( current != 0 );
 281
 282     /* Pass 3.
 283        Removes non-dominant points with 1-length support region */
 284     current = temp.next;
 285     assert( current );
 286     prev_current = &temp;
 287
 288     do
 289     {
 290         if( current->k == 1 )
 291         {
 292             s = current->s;
 293             i = (int)(current - array);
 294
 295             i1 = i - 1;
 296             i1 += i1 < 0 ? len : 0;
 297
 298             i2 = i + 1;
 299             i2 -= i2 >= len ? len : 0;
 300
 301             if( s <= array[i1].s || s <= array[i2].s )
 302             {
 303                 prev_current->next = current->next;
 304                 current->s = 0;
 305             }
 306             else
 307                 prev_current = current;
 308         }
 309         else
 310             prev_current = current;
 311         current = current->next;
 312     }
 313     while( current != 0 );
 314
 315     if( method == CV_CHAIN_APPROX_TC89_KCOS )
 316         goto copy_vect;
 317
 318     /* Pass 4.
 319        Cleans remained couples of points */
 320     assert( temp.next );
 321
 322     if( array[0].s != 0 && array[len - 1].s != 0 )      /* specific case */
 323     {
 324         for( i1 = 1; i1 < len && array[i1].s != 0; i1++ )
 325         {
 326             array[i1 - 1].s = 0;
 327         }
 328         if( i1 == len )
 329             goto copy_vect;     /* all points survived */
 330         i1--;
 331
 332         for( i2 = len - 2; i2 > 0 && array[i2].s != 0; i2-- )
 333         {
 334             array[i2].next = 0;
 335             array[i2 + 1].s = 0;
 336         }
 337         i2++;
 338
 339         if( i1 == 0 && i2 == len - 1 )  /* only two points */
 340         {
 341             i1 = (int)(array[0].next - array);
 342             array[len] = array[0];      /* move to the end */
 343             array[len].next = 0;
 344             array[len - 1].next = array + len;
 345         }
 346         temp.next = array + i1;
 347     }
 348
 349     current = temp.next;
 350     first = prev_current = &temp;
 351     count = 1;
 352
 353     /* do last pass */
 354     do
 355     {
 356         if( current->next == 0 || current->next - current != 1 )
 357         {
 358             if( count >= 2 )
 359             {
 360                 if( count == 2 )
 361                 {
 362                     int s1 = prev_current->s;
 363                     int s2 = current->s;
 364
 365                     if( s1 > s2 || s1 == s2 && prev_current->k <= current->k )
 366                         /* remove second */
 367                         prev_current->next = current->next;
 368                     else
 369                         /* remove first */
 370                         first->next = current;
 371                 }
 372                 else
 373                     first->next->next = current;
 374             }
 375             first = current;
 376             count = 1;
 377         }
 378         else
 379             count++;
 380         prev_current = current;
 381         current = current->next;
 382     }
 383     while( current != 0 );
 384
 385   copy_vect:
 386
 387     /* gather points */
 388     current = temp.next;
 389     assert( current );
 390
 391     do
 392     {
 393         CV_WRITE_SEQ_ELEM( current->pt, writer );
 394         current = current->next;
 395     }
 396     while( current != 0 );
 397
 398 exit_function:
 399
 400     *contour = cvEndWriteSeq( &writer );
 401
 402     assert( writer.seq->total > 0 );
 403
 404     if( buffer != local_buffer )
 405         cvFree( &buffer );
 406     return CV_OK;
 407 }
 408
 409
 410 /*Applies some approximation algorithm to chain-coded contour(s) and
 411   converts it/them to polygonal representation */
 412 CV_IMPL CvSeq*
 413 cvApproxChains( CvSeq*              src_seq,
 414                 CvMemStorage*       storage,
 415                 int                 method,
 416                 double              /*parameter*/,
 417                 int                 minimal_perimeter,
 418                 int                 recursive )
 419 {
 420     CvSeq *prev_contour = 0, *parent = 0;
 421     CvSeq *dst_seq = 0;
 422
 423     CV_FUNCNAME( "cvApproxChains" );
 424
 425     __BEGIN__;
 426
 427     if( !src_seq || !storage )
 428         CV_ERROR( CV_StsNullPtr, "" );
 429     if( method > CV_CHAIN_APPROX_TC89_KCOS || method <= 0 || minimal_perimeter < 0 )
 430         CV_ERROR( CV_StsOutOfRange, "" );
 431
 432     while( src_seq != 0 )
 433     {
 434         int len = src_seq->total;
 435
 436         if( len >= minimal_perimeter )
 437         {
 438             CvSeq *contour;
 439
 440             switch( method )
 441             {
 442             case CV_CHAIN_APPROX_NONE:
 443             case CV_CHAIN_APPROX_SIMPLE:
 444             case CV_CHAIN_APPROX_TC89_L1:
 445             case CV_CHAIN_APPROX_TC89_KCOS:
 446                 IPPI_CALL( icvApproximateChainTC89( (CvChain *) src_seq,
 447                                                     sizeof( CvContour ), storage,
 448                                                     (CvSeq**)&contour, method ));
 449                 break;
 450             default:
 451                 assert(0);
 452                 CV_ERROR( CV_StsOutOfRange, "" );
 453             }
 454
 455             assert( contour );
 456
 457             if( contour->total > 0 )
 458             {
 459                 cvBoundingRect( contour, 1 );
 460
 461                 contour->v_prev = parent;
 462                 contour->h_prev = prev_contour;
 463
 464                 if( prev_contour )
 465                     prev_contour->h_next = contour;
 466                 else if( parent )
 467                     parent->v_next = contour;
 468                 prev_contour = contour;
 469                 if( !dst_seq )
 470                     dst_seq = prev_contour;
 471             }
 472             else                /* if resultant contour has zero length, skip it */
 473             {
 474                 len = -1;
 475             }
 476         }
 477
 478         if( !recursive )
 479             break;
 480
 481         if( src_seq->v_next && len >= minimal_perimeter )
 482         {
 483             assert( prev_contour != 0 );
 484             parent = prev_contour;
 485             prev_contour = 0;
 486             src_seq = src_seq->v_next;
 487         }
 488         else
 489         {
 490             while( src_seq->h_next == 0 )
 491             {
 492                 src_seq = src_seq->v_prev;
 493                 if( src_seq == 0 )
 494                     break;
 495                 prev_contour = parent;
 496                 if( parent )
 497                     parent = parent->v_prev;
 498             }
 499             if( src_seq )
 500                 src_seq = src_seq->h_next;
 501         }
 502     }
 503
 504     __END__;
 505
 506     return dst_seq;
 507 }
 508
 509
 510 /****************************************************************************************\
 511 *                               Polygonal Approximation                                  *
 512 \****************************************************************************************/
 513
 514 /* Ramer-Douglas-Peucker algorithm for polygon simplification */
 515
 516 /* the version for integer point coordinates */
 517 static CvStatus
 518 icvApproxPolyDP_32s( CvSeq* src_contour, int header_size,
 519                      CvMemStorage* storage,
 520                      CvSeq** dst_contour, float eps )
 521 {
 522     int             init_iters = 3;
 523     CvSlice         slice = {0, 0}, right_slice = {0, 0};
 524     CvSeqReader     reader, reader2;
 525     CvSeqWriter     writer;
 526     CvPoint         start_pt = {INT_MIN, INT_MIN}, end_pt = {0, 0}, pt = {0,0};
 527     int             i = 0, j, count = src_contour->total, new_count;
 528     int             is_closed = CV_IS_SEQ_CLOSED( src_contour );
 529     int             le_eps = 0;
 530     CvMemStorage*   temp_storage = 0;
 531     CvSeq*          stack = 0;
 532
 533     assert( CV_SEQ_ELTYPE(src_contour) == CV_32SC2 );
 534     cvStartWriteSeq( src_contour->flags, header_size, sizeof(pt), storage, &writer );
 535
 536     if( src_contour->total == 0  )
 537     {
 538         *dst_contour = cvEndWriteSeq( &writer );
 539         return CV_OK;
 540     }
 541
 542     temp_storage = cvCreateChildMemStorage( storage );
 543
 544     assert( src_contour->first != 0 );
 545     stack = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvSlice), temp_storage );
 546     eps *= eps;
 547     cvStartReadSeq( src_contour, &reader, 0 );
 548
 549     if( !is_closed )
 550     {
 551         right_slice.start_index = count;
 552         end_pt = *(CvPoint*)(reader.ptr);
 553         start_pt = *(CvPoint*)cvGetSeqElem( src_contour, -1 );
 554
 555         if( start_pt.x != end_pt.x || start_pt.y != end_pt.y )
 556         {
 557             slice.start_index = 0;
 558             slice.end_index = count - 1;
 559             cvSeqPush( stack, &slice );
 560         }
 561         else
 562         {
 563             is_closed = 1;
 564             init_iters = 1;
 565         }
 566     }
 567
 568     if( is_closed )
 569     {
 570         /* 1. Find approximately two farthest points of the contour */
 571         right_slice.start_index = 0;
 572
 573         for( i = 0; i < init_iters; i++ )
 574         {
 575             int max_dist = 0;
 576             cvSetSeqReaderPos( &reader, right_slice.start_index, 1 );
 577             CV_READ_SEQ_ELEM( start_pt, reader );   /* read the first point */
 578
 579             for( j = 1; j < count; j++ )
 580             {
 581                 int dx, dy, dist;
 582
 583                 CV_READ_SEQ_ELEM( pt, reader );
 584                 dx = pt.x - start_pt.x;
 585                 dy = pt.y - start_pt.y;
 586
 587                 dist = dx * dx + dy * dy;
 588
 589                 if( dist > max_dist )
 590                 {
 591                     max_dist = dist;
 592                     right_slice.start_index = j;
 593                 }
 594             }
 595
 596             le_eps = max_dist <= eps;
 597         }
 598
 599         /* 2. initialize the stack */
 600         if( !le_eps )
 601         {
 602             slice.start_index = cvGetSeqReaderPos( &reader );
 603             slice.end_index = right_slice.start_index += slice.start_index;
 604
 605             right_slice.start_index -= right_slice.start_index >= count ? count : 0;
 606             right_slice.end_index = slice.start_index;
 607             if( right_slice.end_index < right_slice.start_index )
 608                 right_slice.end_index += count;
 609
 610             cvSeqPush( stack, &right_slice );
 611             cvSeqPush( stack, &slice );
 612         }
 613         else
 614             CV_WRITE_SEQ_ELEM( start_pt, writer );
 615     }
 616
 617     /* 3. run recursive process */
 618     while( stack->total != 0 )
 619     {
 620         cvSeqPop( stack, &slice );
 621
 622         cvSetSeqReaderPos( &reader, slice.end_index );
 623         CV_READ_SEQ_ELEM( end_pt, reader );
 624
 625         cvSetSeqReaderPos( &reader, slice.start_index );
 626         CV_READ_SEQ_ELEM( start_pt, reader );
 627
 628         if( slice.end_index > slice.start_index + 1 )
 629         {
 630             int dx, dy, dist, max_dist = 0;
 631
 632             dx = end_pt.x - start_pt.x;
 633             dy = end_pt.y - start_pt.y;
 634
 635             assert( dx != 0 || dy != 0 );
 636
 637             for( i = slice.start_index + 1; i < slice.end_index; i++ )
 638             {
 639                 CV_READ_SEQ_ELEM( pt, reader );
 640                 dist = abs((pt.y - start_pt.y) * dx - (pt.x - start_pt.x) * dy);
 641
 642                 if( dist > max_dist )
 643                 {
 644                     max_dist = dist;
 645                     right_slice.start_index = i;
 646                 }
 647             }
 648
 649             le_eps = (double)max_dist * max_dist <= eps * ((double)dx * dx + (double)dy * dy);
 650         }
 651         else
 652         {
 653             assert( slice.end_index > slice.start_index );
 654             le_eps = 1;
 655             /* read starting point */
 656             cvSetSeqReaderPos( &reader, slice.start_index );
 657             CV_READ_SEQ_ELEM( start_pt, reader );
 658         }
 659
 660         if( le_eps )
 661         {
 662             CV_WRITE_SEQ_ELEM( start_pt, writer );
 663         }
 664         else
 665         {
 666             right_slice.end_index = slice.end_index;
 667             slice.end_index = right_slice.start_index;
 668             cvSeqPush( stack, &right_slice );
 669             cvSeqPush( stack, &slice );
 670         }
 671     }
 672
 673     is_closed = CV_IS_SEQ_CLOSED( src_contour );
 674     if( !is_closed )
 675         CV_WRITE_SEQ_ELEM( end_pt, writer );
 676
 677     *dst_contour = cvEndWriteSeq( &writer );
 678
 679     cvStartReadSeq( *dst_contour, &reader, is_closed );
 680     CV_READ_SEQ_ELEM( start_pt, reader );
 681
 682     reader2 = reader;
 683     CV_READ_SEQ_ELEM( pt, reader );
 684
 685     new_count = count = (*dst_contour)->total;
 686     for( i = !is_closed; i < count - !is_closed && new_count > 2; i++ )
 687     {
 688         int dx, dy, dist;
 689         CV_READ_SEQ_ELEM( end_pt, reader );
 690
 691         dx = end_pt.x - start_pt.x;
 692         dy = end_pt.y - start_pt.y;
 693         dist = abs((pt.x - start_pt.x)*dy - (pt.y - start_pt.y)*dx);
 694         if( (double)dist * dist <= 0.5*eps*((double)dx*dx + (double)dy*dy) && dx != 0 && dy != 0 )
 695         {
 696             new_count--;
 697             *((CvPoint*)reader2.ptr) = start_pt = end_pt;
 698             CV_NEXT_SEQ_ELEM( sizeof(pt), reader2 );
 699             CV_READ_SEQ_ELEM( pt, reader );
 700             i++;
 701             continue;
 702         }
 703         *((CvPoint*)reader2.ptr) = start_pt = pt;
 704         CV_NEXT_SEQ_ELEM( sizeof(pt), reader2 );
 705         pt = end_pt;
 706     }
 707
 708     if( !is_closed )
 709         *((CvPoint*)reader2.ptr) = pt;
 710
 711     if( new_count < count )
 712         cvSeqPopMulti( *dst_contour, 0, count - new_count );
 713
 714     cvReleaseMemStorage( &temp_storage );
 715
 716     return CV_OK;
 717 }
 718
 719
 720 /* the version for integer point coordinates */
 721 static CvStatus
 722 icvApproxPolyDP_32f( CvSeq* src_contour, int header_size,
 723                      CvMemStorage* storage,
 724                      CvSeq** dst_contour, float eps )
 725 {
 726     int             init_iters = 3;
 727     CvSlice         slice = {0, 0}, right_slice = {0, 0};
 728     CvSeqReader     reader, reader2;
 729     CvSeqWriter     writer;
 730     CvPoint2D32f    start_pt = {-1e6f, -1e6f}, end_pt = {0, 0}, pt = {0,0};
 731     int             i = 0, j, count = src_contour->total, new_count;
 732     int             is_closed = CV_IS_SEQ_CLOSED( src_contour );
 733     int             le_eps = 0;
 734     CvMemStorage*   temp_storage = 0;
 735     CvSeq*          stack = 0;
 736
 737     assert( CV_SEQ_ELTYPE(src_contour) == CV_32FC2 );
 738     cvStartWriteSeq( src_contour->flags, header_size, sizeof(pt), storage, &writer );
 739
 740     if( src_contour->total == 0  )
 741     {
 742         *dst_contour = cvEndWriteSeq( &writer );
 743         return CV_OK;
 744     }
 745
 746     temp_storage = cvCreateChildMemStorage( storage );
 747
 748     assert( src_contour->first != 0 );
 749     stack = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvSlice), temp_storage );
 750     eps *= eps;
 751     cvStartReadSeq( src_contour, &reader, 0 );
 752
 753     if( !is_closed )
 754     {
 755         right_slice.start_index = count;
 756         end_pt = *(CvPoint2D32f*)(reader.ptr);
 757         start_pt = *(CvPoint2D32f*)cvGetSeqElem( src_contour, -1 );
 758
 759         if( fabs(start_pt.x - end_pt.x) > FLT_EPSILON ||
 760             fabs(start_pt.y - end_pt.y) > FLT_EPSILON )
 761         {
 762             slice.start_index = 0;
 763             slice.end_index = count - 1;
 764             cvSeqPush( stack, &slice );
 765         }
 766         else
 767         {
 768             is_closed = 1;
 769             init_iters = 1;
 770         }
 771     }
 772
 773     if( is_closed )
 774     {
 775         /* 1. Find approximately two farthest points of the contour */
 776         right_slice.start_index = 0;
 777
 778         for( i = 0; i < init_iters; i++ )
 779         {
 780             double max_dist = 0;
 781             cvSetSeqReaderPos( &reader, right_slice.start_index, 1 );
 782             CV_READ_SEQ_ELEM( start_pt, reader );   /* read the first point */
 783
 784             for( j = 1; j < count; j++ )
 785             {
 786                 double dx, dy, dist;
 787
 788                 CV_READ_SEQ_ELEM( pt, reader );
 789                 dx = pt.x - start_pt.x;
 790                 dy = pt.y - start_pt.y;
 791
 792                 dist = dx * dx + dy * dy;
 793
 794                 if( dist > max_dist )
 795                 {
 796                     max_dist = dist;
 797                     right_slice.start_index = j;
 798                 }
 799             }
 800
 801             le_eps = max_dist <= eps;
 802         }
 803
 804         /* 2. initialize the stack */
 805         if( !le_eps )
 806         {
 807             slice.start_index = cvGetSeqReaderPos( &reader );
 808             slice.end_index = right_slice.start_index += slice.start_index;
 809
 810             right_slice.start_index -= right_slice.start_index >= count ? count : 0;
 811             right_slice.end_index = slice.start_index;
 812             if( right_slice.end_index < right_slice.start_index )
 813                 right_slice.end_index += count;
 814
 815             cvSeqPush( stack, &right_slice );
 816             cvSeqPush( stack, &slice );
 817         }
 818         else
 819             CV_WRITE_SEQ_ELEM( start_pt, writer );
 820     }
 821
 822     /* 3. run recursive process */
 823     while( stack->total != 0 )
 824     {
 825         cvSeqPop( stack, &slice );
 826
 827         cvSetSeqReaderPos( &reader, slice.end_index );
 828         CV_READ_SEQ_ELEM( end_pt, reader );
 829
 830         cvSetSeqReaderPos( &reader, slice.start_index );
 831         CV_READ_SEQ_ELEM( start_pt, reader );
 832
 833         if( slice.end_index > slice.start_index + 1 )
 834         {
 835             double dx, dy, dist, max_dist = 0;
 836
 837             dx = end_pt.x - start_pt.x;
 838             dy = end_pt.y - start_pt.y;
 839
 840             assert( dx != 0 || dy != 0 );
 841
 842             for( i = slice.start_index + 1; i < slice.end_index; i++ )
 843             {
 844                 CV_READ_SEQ_ELEM( pt, reader );
 845                 dist = fabs((pt.y - start_pt.y) * dx - (pt.x - start_pt.x) * dy);
 846
 847                 if( dist > max_dist )
 848                 {
 849                     max_dist = dist;
 850                     right_slice.start_index = i;
 851                 }
 852             }
 853
 854             le_eps = (double)max_dist * max_dist <= eps * ((double)dx * dx + (double)dy * dy);
 855         }
 856         else
 857         {
 858             assert( slice.end_index > slice.start_index );
 859             le_eps = 1;
 860             /* read starting point */
 861             cvSetSeqReaderPos( &reader, slice.start_index );
 862             CV_READ_SEQ_ELEM( start_pt, reader );
 863         }
 864
 865         if( le_eps )
 866         {
 867             CV_WRITE_SEQ_ELEM( start_pt, writer );
 868         }
 869         else
 870         {
 871             right_slice.end_index = slice.end_index;
 872             slice.end_index = right_slice.start_index;
 873             cvSeqPush( stack, &right_slice );
 874             cvSeqPush( stack, &slice );
 875         }
 876     }
 877
 878     is_closed = CV_IS_SEQ_CLOSED( src_contour );
 879     if( !is_closed )
 880         CV_WRITE_SEQ_ELEM( end_pt, writer );
 881
 882     *dst_contour = cvEndWriteSeq( &writer );
 883
 884     cvStartReadSeq( *dst_contour, &reader, is_closed );
 885     CV_READ_SEQ_ELEM( start_pt, reader );
 886
 887     reader2 = reader;
 888     CV_READ_SEQ_ELEM( pt, reader );
 889
 890     new_count = count = (*dst_contour)->total;
 891     for( i = !is_closed; i < count - !is_closed && new_count > 2; i++ )
 892     {
 893         double dx, dy, dist;
 894         CV_READ_SEQ_ELEM( end_pt, reader );
 895
 896         dx = end_pt.x - start_pt.x;
 897         dy = end_pt.y - start_pt.y;
 898         dist = fabs((pt.x - start_pt.x)*dy - (pt.y - start_pt.y)*dx);
 899         if( (double)dist * dist <= 0.5*eps*((double)dx*dx + (double)dy*dy) )
 900         {
 901             new_count--;
 902             *((CvPoint2D32f*)reader2.ptr) = start_pt = end_pt;
 903             CV_NEXT_SEQ_ELEM( sizeof(pt), reader2 );
 904             CV_READ_SEQ_ELEM( pt, reader );
 905             i++;
 906             continue;
 907         }
 908         *((CvPoint2D32f*)reader2.ptr) = start_pt = pt;
 909         CV_NEXT_SEQ_ELEM( sizeof(pt), reader2 );
 910         pt = end_pt;
 911     }
 912
 913     if( !is_closed )
 914         *((CvPoint2D32f*)reader2.ptr) = pt;
 915
 916     if( new_count < count )
 917         cvSeqPopMulti( *dst_contour, 0, count - new_count );
 918
 919     cvReleaseMemStorage( &temp_storage );
 920
 921     return CV_OK;
 922 }
 923
 924
 925 CV_IMPL CvSeq*
 926 cvApproxPoly( const void*  array, int  header_size,
 927               CvMemStorage*  storage, int  method,
 928               double  parameter, int parameter2 )
 929 {
 930     CvSeq* dst_seq = 0;
 931     CvSeq *prev_contour = 0, *parent = 0;
 932     CvContour contour_header;
 933     CvSeq* src_seq = 0;
 934     CvSeqBlock block;
 935     int recursive = 0;
 936
 937     CV_FUNCNAME( "cvApproxPoly" );
 938
 939     __BEGIN__;
 940
 941     if( CV_IS_SEQ( array ))
 942     {
 943         src_seq = (CvSeq*)array;
 944         if( !CV_IS_SEQ_POLYLINE( src_seq ))
 945             CV_ERROR( CV_StsBadArg, "Unsupported sequence type" );
 946
 947         recursive = parameter2;
 948
 949         if( !storage )
 950             storage = src_seq->storage;
 951     }
 952     else
 953     {
 954         CV_CALL( src_seq = cvPointSeqFromMat(
 955             CV_SEQ_KIND_CURVE | (parameter2 ? CV_SEQ_FLAG_CLOSED : 0),
 956             array, &contour_header, &block ));
 957     }
 958
 959     if( !storage )
 960         CV_ERROR( CV_StsNullPtr, "NULL storage pointer " );
 961
 962     if( header_size < 0 )
 963         CV_ERROR( CV_StsOutOfRange, "header_size is negative. "
 964                   "Pass 0 to make the destination header_size == input header_size" );
 965
 966     if( header_size == 0 )
 967         header_size = src_seq->header_size;
 968
 969     if( !CV_IS_SEQ_POLYLINE( src_seq ))
 970     {
 971         if( CV_IS_SEQ_CHAIN( src_seq ))
 972         {
 973             CV_ERROR( CV_StsBadArg, "Input curves are not polygonal. "
 974                                     "Use cvApproxChains first" );
 975         }
 976         else
 977         {
 978             CV_ERROR( CV_StsBadArg, "Input curves have uknown type" );
 979         }
 980     }
 981
 982     if( header_size == 0 )
 983         header_size = src_seq->header_size;
 984
 985     if( header_size < (int)sizeof(CvContour) )
 986         CV_ERROR( CV_StsBadSize, "New header size must be non-less than sizeof(CvContour)" );
 987
 988     if( method != CV_POLY_APPROX_DP )
 989         CV_ERROR( CV_StsOutOfRange, "Unknown approximation method" );
 990
 991     while( src_seq != 0 )
 992     {
 993         CvSeq *contour = 0;
 994
 995         switch (method)
 996         {
 997         case CV_POLY_APPROX_DP:
 998             if( parameter < 0 )
 999                 CV_ERROR( CV_StsOutOfRange, "Accuracy must be non-negative" );
1000
1001             if( CV_SEQ_ELTYPE(src_seq) == CV_32SC2 )
1002             {
1003                 IPPI_CALL( icvApproxPolyDP_32s( src_seq, header_size, storage,
1004                                                 &contour, (float)parameter ));
1005             }
1006             else
1007             {
1008                 IPPI_CALL( icvApproxPolyDP_32f( src_seq, header_size, storage,
1009                                                 &contour, (float)parameter ));
1010             }
1011             break;
1012         default:
1013             assert(0);
1014             CV_ERROR( CV_StsBadArg, "Invalid approximation method" );
1015         }
1016
1017         assert( contour );
1018
1019         if( header_size >= (int)sizeof(CvContour))
1020             cvBoundingRect( contour, 1 );
1021
1022         contour->v_prev = parent;
1023         contour->h_prev = prev_contour;
1024
1025         if( prev_contour )
1026             prev_contour->h_next = contour;
1027         else if( parent )
1028             parent->v_next = contour;
1029         prev_contour = contour;
1030         if( !dst_seq )
1031             dst_seq = prev_contour;
1032
1033         if( !recursive )
1034             break;
1035
1036         if( src_seq->v_next )
1037         {
1038             assert( prev_contour != 0 );
1039             parent = prev_contour;
1040             prev_contour = 0;
1041             src_seq = src_seq->v_next;
1042         }
1043         else
1044         {
1045             while( src_seq->h_next == 0 )
1046             {
1047                 src_seq = src_seq->v_prev;
1048                 if( src_seq == 0 )
1049                     break;
1050                 prev_contour = parent;
1051                 if( parent )
1052                     parent = parent->v_prev;
1053             }
1054             if( src_seq )
1055                 src_seq = src_seq->h_next;
1056         }
1057     }
1058
1059     __END__;
1060
1061     return dst_seq;
1062 }
1063
1064 /* End of file. */