Home | History | Annotate | Download | only in src
      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. */
   1065