| /*M/////////////////////////////////////////////////////////////////////////////////////// |
| // |
| // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. |
| // |
| // By downloading, copying, installing or using the software you agree to this license. |
| // If you do not agree to this license, do not download, install, |
| // copy or use the software. |
| // |
| // |
| // Intel License Agreement |
| // For Open Source Computer Vision Library |
| // |
| // Copyright (C) 2000, Intel Corporation, all rights reserved. |
| // Third party copyrights are property of their respective owners. |
| // |
| // Redistribution and use in source and binary forms, with or without modification, |
| // are permitted provided that the following conditions are met: |
| // |
| // * Redistribution's of source code must retain the above copyright notice, |
| // this list of conditions and the following disclaimer. |
| // |
| // * Redistribution's in binary form must reproduce the above copyright notice, |
| // this list of conditions and the following disclaimer in the documentation |
| // and/or other materials provided with the distribution. |
| // |
| // * The name of Intel Corporation may not be used to endorse or promote products |
| // derived from this software without specific prior written permission. |
| // |
| // This software is provided by the copyright holders and contributors "as is" and |
| // any express or implied warranties, including, but not limited to, the implied |
| // warranties of merchantability and fitness for a particular purpose are disclaimed. |
| // In no event shall the Intel Corporation or contributors be liable for any direct, |
| // indirect, incidental, special, exemplary, or consequential damages |
| // (including, but not limited to, procurement of substitute goods or services; |
| // loss of use, data, or profits; or business interruption) however caused |
| // and on any theory of liability, whether in contract, strict liability, |
| // or tort (including negligence or otherwise) arising in any way out of |
| // the use of this software, even if advised of the possibility of such damage. |
| // |
| //M*/ |
| #include "_cv.h" |
| |
| /* calculates length of a curve (e.g. contour perimeter) */ |
| CV_IMPL double |
| cvArcLength( const void *array, CvSlice slice, int is_closed ) |
| { |
| double perimeter = 0; |
| |
| CV_FUNCNAME( "cvArcLength" ); |
| |
| __BEGIN__; |
| |
| int i, j = 0, count; |
| const int N = 16; |
| float buf[N]; |
| CvMat buffer = cvMat( 1, N, CV_32F, buf ); |
| CvSeqReader reader; |
| CvContour contour_header; |
| CvSeq* contour = 0; |
| CvSeqBlock block; |
| |
| if( CV_IS_SEQ( array )) |
| { |
| contour = (CvSeq*)array; |
| if( !CV_IS_SEQ_POLYLINE( contour )) |
| CV_ERROR( CV_StsBadArg, "Unsupported sequence type" ); |
| if( is_closed < 0 ) |
| is_closed = CV_IS_SEQ_CLOSED( contour ); |
| } |
| else |
| { |
| is_closed = is_closed > 0; |
| CV_CALL( contour = cvPointSeqFromMat( |
| CV_SEQ_KIND_CURVE | (is_closed ? CV_SEQ_FLAG_CLOSED : 0), |
| array, &contour_header, &block )); |
| } |
| |
| if( contour->total > 1 ) |
| { |
| int is_float = CV_SEQ_ELTYPE( contour ) == CV_32FC2; |
| |
| cvStartReadSeq( contour, &reader, 0 ); |
| cvSetSeqReaderPos( &reader, slice.start_index ); |
| count = cvSliceLength( slice, contour ); |
| |
| count -= !is_closed && count == contour->total; |
| |
| /* scroll the reader by 1 point */ |
| reader.prev_elem = reader.ptr; |
| CV_NEXT_SEQ_ELEM( sizeof(CvPoint), reader ); |
| |
| for( i = 0; i < count; i++ ) |
| { |
| float dx, dy; |
| |
| if( !is_float ) |
| { |
| CvPoint* pt = (CvPoint*)reader.ptr; |
| CvPoint* prev_pt = (CvPoint*)reader.prev_elem; |
| |
| dx = (float)pt->x - (float)prev_pt->x; |
| dy = (float)pt->y - (float)prev_pt->y; |
| } |
| else |
| { |
| CvPoint2D32f* pt = (CvPoint2D32f*)reader.ptr; |
| CvPoint2D32f* prev_pt = (CvPoint2D32f*)reader.prev_elem; |
| |
| dx = pt->x - prev_pt->x; |
| dy = pt->y - prev_pt->y; |
| } |
| |
| reader.prev_elem = reader.ptr; |
| CV_NEXT_SEQ_ELEM( contour->elem_size, reader ); |
| |
| buffer.data.fl[j] = dx * dx + dy * dy; |
| if( ++j == N || i == count - 1 ) |
| { |
| buffer.cols = j; |
| cvPow( &buffer, &buffer, 0.5 ); |
| for( ; j > 0; j-- ) |
| perimeter += buffer.data.fl[j-1]; |
| } |
| } |
| } |
| |
| __END__; |
| |
| return perimeter; |
| } |
| |
| |
| static CvStatus |
| icvFindCircle( CvPoint2D32f pt0, CvPoint2D32f pt1, |
| CvPoint2D32f pt2, CvPoint2D32f * center, float *radius ) |
| { |
| double x1 = (pt0.x + pt1.x) * 0.5; |
| double dy1 = pt0.x - pt1.x; |
| double x2 = (pt1.x + pt2.x) * 0.5; |
| double dy2 = pt1.x - pt2.x; |
| double y1 = (pt0.y + pt1.y) * 0.5; |
| double dx1 = pt1.y - pt0.y; |
| double y2 = (pt1.y + pt2.y) * 0.5; |
| double dx2 = pt2.y - pt1.y; |
| double t = 0; |
| |
| CvStatus result = CV_OK; |
| |
| if( icvIntersectLines( x1, dx1, y1, dy1, x2, dx2, y2, dy2, &t ) >= 0 ) |
| { |
| center->x = (float) (x2 + dx2 * t); |
| center->y = (float) (y2 + dy2 * t); |
| *radius = (float) icvDistanceL2_32f( *center, pt0 ); |
| } |
| else |
| { |
| center->x = center->y = 0.f; |
| radius = 0; |
| result = CV_NOTDEFINED_ERR; |
| } |
| |
| return result; |
| } |
| |
| |
| CV_INLINE double icvIsPtInCircle( CvPoint2D32f pt, CvPoint2D32f center, float radius ) |
| { |
| double dx = pt.x - center.x; |
| double dy = pt.y - center.y; |
| return (double)radius*radius - dx*dx - dy*dy; |
| } |
| |
| |
| static int |
| icvFindEnslosingCicle4pts_32f( CvPoint2D32f * pts, CvPoint2D32f * _center, float *_radius ) |
| { |
| int shuffles[4][4] = { {0, 1, 2, 3}, {0, 1, 3, 2}, {2, 3, 0, 1}, {2, 3, 1, 0} }; |
| |
| int idxs[4] = { 0, 1, 2, 3 }; |
| int i, j, k = 1, mi = 0; |
| float max_dist = 0; |
| CvPoint2D32f center; |
| CvPoint2D32f min_center; |
| float radius, min_radius = FLT_MAX; |
| CvPoint2D32f res_pts[4]; |
| |
| center = min_center = pts[0]; |
| radius = 1.f; |
| |
| for( i = 0; i < 4; i++ ) |
| for( j = i + 1; j < 4; j++ ) |
| { |
| float dist = icvDistanceL2_32f( pts[i], pts[j] ); |
| |
| if( max_dist < dist ) |
| { |
| max_dist = dist; |
| idxs[0] = i; |
| idxs[1] = j; |
| } |
| } |
| |
| if( max_dist == 0 ) |
| goto function_exit; |
| |
| k = 2; |
| for( i = 0; i < 4; i++ ) |
| { |
| for( j = 0; j < k; j++ ) |
| if( i == idxs[j] ) |
| break; |
| if( j == k ) |
| idxs[k++] = i; |
| } |
| |
| center = cvPoint2D32f( (pts[idxs[0]].x + pts[idxs[1]].x)*0.5f, |
| (pts[idxs[0]].y + pts[idxs[1]].y)*0.5f ); |
| radius = (float)(icvDistanceL2_32f( pts[idxs[0]], center )*1.03); |
| if( radius < 1.f ) |
| radius = 1.f; |
| |
| if( icvIsPtInCircle( pts[idxs[2]], center, radius ) >= 0 && |
| icvIsPtInCircle( pts[idxs[3]], center, radius ) >= 0 ) |
| { |
| k = 2; //rand()%2+2; |
| } |
| else |
| { |
| mi = -1; |
| for( i = 0; i < 4; i++ ) |
| { |
| if( icvFindCircle( pts[shuffles[i][0]], pts[shuffles[i][1]], |
| pts[shuffles[i][2]], ¢er, &radius ) >= 0 ) |
| { |
| radius *= 1.03f; |
| if( radius < 2.f ) |
| radius = 2.f; |
| |
| if( icvIsPtInCircle( pts[shuffles[i][3]], center, radius ) >= 0 && |
| min_radius > radius ) |
| { |
| min_radius = radius; |
| min_center = center; |
| mi = i; |
| } |
| } |
| } |
| assert( mi >= 0 ); |
| if( mi < 0 ) |
| mi = 0; |
| k = 3; |
| center = min_center; |
| radius = min_radius; |
| for( i = 0; i < 4; i++ ) |
| idxs[i] = shuffles[mi][i]; |
| } |
| |
| function_exit: |
| |
| *_center = center; |
| *_radius = radius; |
| |
| /* reorder output points */ |
| for( i = 0; i < 4; i++ ) |
| res_pts[i] = pts[idxs[i]]; |
| |
| for( i = 0; i < 4; i++ ) |
| { |
| pts[i] = res_pts[i]; |
| assert( icvIsPtInCircle( pts[i], center, radius ) >= 0 ); |
| } |
| |
| return k; |
| } |
| |
| |
| CV_IMPL int |
| cvMinEnclosingCircle( const void* array, CvPoint2D32f * _center, float *_radius ) |
| { |
| const int max_iters = 100; |
| const float eps = FLT_EPSILON*2; |
| CvPoint2D32f center = { 0, 0 }; |
| float radius = 0; |
| int result = 0; |
| |
| if( _center ) |
| _center->x = _center->y = 0.f; |
| if( _radius ) |
| *_radius = 0; |
| |
| CV_FUNCNAME( "cvMinEnclosingCircle" ); |
| |
| __BEGIN__; |
| |
| CvSeqReader reader; |
| int i, k, count; |
| CvPoint2D32f pts[8]; |
| CvContour contour_header; |
| CvSeqBlock block; |
| CvSeq* sequence = 0; |
| int is_float; |
| |
| if( !_center || !_radius ) |
| CV_ERROR( CV_StsNullPtr, "Null center or radius pointers" ); |
| |
| if( CV_IS_SEQ(array) ) |
| { |
| sequence = (CvSeq*)array; |
| if( !CV_IS_SEQ_POINT_SET( sequence )) |
| CV_ERROR( CV_StsBadArg, "The passed sequence is not a valid contour" ); |
| } |
| else |
| { |
| CV_CALL( sequence = cvPointSeqFromMat( |
| CV_SEQ_KIND_GENERIC, array, &contour_header, &block )); |
| } |
| |
| if( sequence->total <= 0 ) |
| CV_ERROR_FROM_STATUS( CV_BADSIZE_ERR ); |
| |
| CV_CALL( cvStartReadSeq( sequence, &reader, 0 )); |
| |
| count = sequence->total; |
| is_float = CV_SEQ_ELTYPE(sequence) == CV_32FC2; |
| |
| if( !is_float ) |
| { |
| CvPoint *pt_left, *pt_right, *pt_top, *pt_bottom; |
| CvPoint pt; |
| pt_left = pt_right = pt_top = pt_bottom = (CvPoint *)(reader.ptr); |
| CV_READ_SEQ_ELEM( pt, reader ); |
| |
| for( i = 1; i < count; i++ ) |
| { |
| CvPoint* pt_ptr = (CvPoint*)reader.ptr; |
| CV_READ_SEQ_ELEM( pt, reader ); |
| |
| if( pt.x < pt_left->x ) |
| pt_left = pt_ptr; |
| if( pt.x > pt_right->x ) |
| pt_right = pt_ptr; |
| if( pt.y < pt_top->y ) |
| pt_top = pt_ptr; |
| if( pt.y > pt_bottom->y ) |
| pt_bottom = pt_ptr; |
| } |
| |
| pts[0] = cvPointTo32f( *pt_left ); |
| pts[1] = cvPointTo32f( *pt_right ); |
| pts[2] = cvPointTo32f( *pt_top ); |
| pts[3] = cvPointTo32f( *pt_bottom ); |
| } |
| else |
| { |
| CvPoint2D32f *pt_left, *pt_right, *pt_top, *pt_bottom; |
| CvPoint2D32f pt; |
| pt_left = pt_right = pt_top = pt_bottom = (CvPoint2D32f *) (reader.ptr); |
| CV_READ_SEQ_ELEM( pt, reader ); |
| |
| for( i = 1; i < count; i++ ) |
| { |
| CvPoint2D32f* pt_ptr = (CvPoint2D32f*)reader.ptr; |
| CV_READ_SEQ_ELEM( pt, reader ); |
| |
| if( pt.x < pt_left->x ) |
| pt_left = pt_ptr; |
| if( pt.x > pt_right->x ) |
| pt_right = pt_ptr; |
| if( pt.y < pt_top->y ) |
| pt_top = pt_ptr; |
| if( pt.y > pt_bottom->y ) |
| pt_bottom = pt_ptr; |
| } |
| |
| pts[0] = *pt_left; |
| pts[1] = *pt_right; |
| pts[2] = *pt_top; |
| pts[3] = *pt_bottom; |
| } |
| |
| for( k = 0; k < max_iters; k++ ) |
| { |
| double min_delta = 0, delta; |
| CvPoint2D32f ptfl; |
| |
| icvFindEnslosingCicle4pts_32f( pts, ¢er, &radius ); |
| cvStartReadSeq( sequence, &reader, 0 ); |
| |
| for( i = 0; i < count; i++ ) |
| { |
| if( !is_float ) |
| { |
| ptfl.x = (float)((CvPoint*)reader.ptr)->x; |
| ptfl.y = (float)((CvPoint*)reader.ptr)->y; |
| } |
| else |
| { |
| ptfl = *(CvPoint2D32f*)reader.ptr; |
| } |
| CV_NEXT_SEQ_ELEM( sequence->elem_size, reader ); |
| |
| delta = icvIsPtInCircle( ptfl, center, radius ); |
| if( delta < min_delta ) |
| { |
| min_delta = delta; |
| pts[3] = ptfl; |
| } |
| } |
| result = min_delta >= 0; |
| if( result ) |
| break; |
| } |
| |
| if( !result ) |
| { |
| cvStartReadSeq( sequence, &reader, 0 ); |
| radius = 0.f; |
| |
| for( i = 0; i < count; i++ ) |
| { |
| CvPoint2D32f ptfl; |
| float t, dx, dy; |
| |
| if( !is_float ) |
| { |
| ptfl.x = (float)((CvPoint*)reader.ptr)->x; |
| ptfl.y = (float)((CvPoint*)reader.ptr)->y; |
| } |
| else |
| { |
| ptfl = *(CvPoint2D32f*)reader.ptr; |
| } |
| |
| CV_NEXT_SEQ_ELEM( sequence->elem_size, reader ); |
| dx = center.x - ptfl.x; |
| dy = center.y - ptfl.y; |
| t = dx*dx + dy*dy; |
| radius = MAX(radius,t); |
| } |
| |
| radius = (float)(sqrt(radius)*(1 + eps)); |
| result = 1; |
| } |
| |
| __END__; |
| |
| *_center = center; |
| *_radius = radius; |
| |
| return result; |
| } |
| |
| |
| /* area of a whole sequence */ |
| static CvStatus |
| icvContourArea( const CvSeq* contour, double *area ) |
| { |
| if( contour->total ) |
| { |
| CvSeqReader reader; |
| int lpt = contour->total; |
| double a00 = 0, xi_1, yi_1; |
| int is_float = CV_SEQ_ELTYPE(contour) == CV_32FC2; |
| |
| cvStartReadSeq( contour, &reader, 0 ); |
| |
| if( !is_float ) |
| { |
| xi_1 = ((CvPoint*)(reader.ptr))->x; |
| yi_1 = ((CvPoint*)(reader.ptr))->y; |
| } |
| else |
| { |
| xi_1 = ((CvPoint2D32f*)(reader.ptr))->x; |
| yi_1 = ((CvPoint2D32f*)(reader.ptr))->y; |
| } |
| CV_NEXT_SEQ_ELEM( contour->elem_size, reader ); |
| |
| while( lpt-- > 0 ) |
| { |
| double dxy, xi, yi; |
| |
| if( !is_float ) |
| { |
| xi = ((CvPoint*)(reader.ptr))->x; |
| yi = ((CvPoint*)(reader.ptr))->y; |
| } |
| else |
| { |
| xi = ((CvPoint2D32f*)(reader.ptr))->x; |
| yi = ((CvPoint2D32f*)(reader.ptr))->y; |
| } |
| CV_NEXT_SEQ_ELEM( contour->elem_size, reader ); |
| |
| dxy = xi_1 * yi - xi * yi_1; |
| a00 += dxy; |
| xi_1 = xi; |
| yi_1 = yi; |
| } |
| |
| *area = a00 * 0.5; |
| } |
| else |
| *area = 0; |
| |
| return CV_OK; |
| } |
| |
| |
| /****************************************************************************************\ |
| |
| copy data from one buffer to other buffer |
| |
| \****************************************************************************************/ |
| |
| static CvStatus |
| icvMemCopy( double **buf1, double **buf2, double **buf3, int *b_max ) |
| { |
| int bb; |
| |
| if( (*buf1 == NULL && *buf2 == NULL) || *buf3 == NULL ) |
| return CV_NULLPTR_ERR; |
| |
| bb = *b_max; |
| if( *buf2 == NULL ) |
| { |
| *b_max = 2 * (*b_max); |
| *buf2 = (double *)cvAlloc( (*b_max) * sizeof( double )); |
| |
| if( *buf2 == NULL ) |
| return CV_OUTOFMEM_ERR; |
| |
| memcpy( *buf2, *buf3, bb * sizeof( double )); |
| |
| *buf3 = *buf2; |
| cvFree( buf1 ); |
| *buf1 = NULL; |
| } |
| else |
| { |
| *b_max = 2 * (*b_max); |
| *buf1 = (double *) cvAlloc( (*b_max) * sizeof( double )); |
| |
| if( *buf1 == NULL ) |
| return CV_OUTOFMEM_ERR; |
| |
| memcpy( *buf1, *buf3, bb * sizeof( double )); |
| |
| *buf3 = *buf1; |
| cvFree( buf2 ); |
| *buf2 = NULL; |
| } |
| return CV_OK; |
| } |
| |
| |
| /* area of a contour sector */ |
| static CvStatus icvContourSecArea( CvSeq * contour, CvSlice slice, double *area ) |
| { |
| CvPoint pt; /* pointer to points */ |
| CvPoint pt_s, pt_e; /* first and last points */ |
| CvSeqReader reader; /* points reader of contour */ |
| |
| int p_max = 2, p_ind; |
| int lpt, flag, i; |
| double a00; /* unnormalized moments m00 */ |
| double xi, yi, xi_1, yi_1, x0, y0, dxy, sk, sk1, t; |
| double x_s, y_s, nx, ny, dx, dy, du, dv; |
| double eps = 1.e-5; |
| double *p_are1, *p_are2, *p_are; |
| |
| assert( contour != NULL ); |
| |
| if( contour == NULL ) |
| return CV_NULLPTR_ERR; |
| |
| if( !CV_IS_SEQ_POLYGON( contour )) |
| return CV_BADFLAG_ERR; |
| |
| lpt = cvSliceLength( slice, contour ); |
| /*if( n2 >= n1 ) |
| lpt = n2 - n1 + 1; |
| else |
| lpt = contour->total - n1 + n2 + 1;*/ |
| |
| if( contour->total && lpt > 2 ) |
| { |
| a00 = x0 = y0 = xi_1 = yi_1 = 0; |
| sk1 = 0; |
| flag = 0; |
| dxy = 0; |
| p_are1 = (double *) cvAlloc( p_max * sizeof( double )); |
| |
| if( p_are1 == NULL ) |
| return CV_OUTOFMEM_ERR; |
| |
| p_are = p_are1; |
| p_are2 = NULL; |
| |
| cvStartReadSeq( contour, &reader, 0 ); |
| cvSetSeqReaderPos( &reader, slice.start_index ); |
| CV_READ_SEQ_ELEM( pt_s, reader ); |
| p_ind = 0; |
| cvSetSeqReaderPos( &reader, slice.end_index ); |
| CV_READ_SEQ_ELEM( pt_e, reader ); |
| |
| /* normal coefficients */ |
| nx = pt_s.y - pt_e.y; |
| ny = pt_e.x - pt_s.x; |
| cvSetSeqReaderPos( &reader, slice.start_index ); |
| |
| while( lpt-- > 0 ) |
| { |
| CV_READ_SEQ_ELEM( pt, reader ); |
| |
| if( flag == 0 ) |
| { |
| xi_1 = (double) pt.x; |
| yi_1 = (double) pt.y; |
| x0 = xi_1; |
| y0 = yi_1; |
| sk1 = 0; |
| flag = 1; |
| } |
| else |
| { |
| xi = (double) pt.x; |
| yi = (double) pt.y; |
| |
| /**************** edges intersection examination **************************/ |
| sk = nx * (xi - pt_s.x) + ny * (yi - pt_s.y); |
| if( (fabs( sk ) < eps && lpt > 0) || sk * sk1 < -eps ) |
| { |
| if( fabs( sk ) < eps ) |
| { |
| dxy = xi_1 * yi - xi * yi_1; |
| a00 = a00 + dxy; |
| dxy = xi * y0 - x0 * yi; |
| a00 = a00 + dxy; |
| |
| if( p_ind >= p_max ) |
| icvMemCopy( &p_are1, &p_are2, &p_are, &p_max ); |
| |
| p_are[p_ind] = a00 / 2.; |
| p_ind++; |
| a00 = 0; |
| sk1 = 0; |
| x0 = xi; |
| y0 = yi; |
| dxy = 0; |
| } |
| else |
| { |
| /* define intersection point */ |
| dv = yi - yi_1; |
| du = xi - xi_1; |
| dx = ny; |
| dy = -nx; |
| if( fabs( du ) > eps ) |
| t = ((yi_1 - pt_s.y) * du + dv * (pt_s.x - xi_1)) / |
| (du * dy - dx * dv); |
| else |
| t = (xi_1 - pt_s.x) / dx; |
| if( t > eps && t < 1 - eps ) |
| { |
| x_s = pt_s.x + t * dx; |
| y_s = pt_s.y + t * dy; |
| dxy = xi_1 * y_s - x_s * yi_1; |
| a00 += dxy; |
| dxy = x_s * y0 - x0 * y_s; |
| a00 += dxy; |
| if( p_ind >= p_max ) |
| icvMemCopy( &p_are1, &p_are2, &p_are, &p_max ); |
| |
| p_are[p_ind] = a00 / 2.; |
| p_ind++; |
| |
| a00 = 0; |
| sk1 = 0; |
| x0 = x_s; |
| y0 = y_s; |
| dxy = x_s * yi - xi * y_s; |
| } |
| } |
| } |
| else |
| dxy = xi_1 * yi - xi * yi_1; |
| |
| a00 += dxy; |
| xi_1 = xi; |
| yi_1 = yi; |
| sk1 = sk; |
| |
| } |
| } |
| |
| xi = x0; |
| yi = y0; |
| dxy = xi_1 * yi - xi * yi_1; |
| |
| a00 += dxy; |
| |
| if( p_ind >= p_max ) |
| icvMemCopy( &p_are1, &p_are2, &p_are, &p_max ); |
| |
| p_are[p_ind] = a00 / 2.; |
| p_ind++; |
| |
| /* common area calculation */ |
| *area = 0; |
| for( i = 0; i < p_ind; i++ ) |
| (*area) += fabs( p_are[i] ); |
| |
| if( p_are1 != NULL ) |
| cvFree( &p_are1 ); |
| else if( p_are2 != NULL ) |
| cvFree( &p_are2 ); |
| |
| return CV_OK; |
| } |
| else |
| return CV_BADSIZE_ERR; |
| } |
| |
| |
| /* external contour area function */ |
| CV_IMPL double |
| cvContourArea( const void *array, CvSlice slice ) |
| { |
| double area = 0; |
| |
| CV_FUNCNAME( "cvContourArea" ); |
| |
| __BEGIN__; |
| |
| CvContour contour_header; |
| CvSeq* contour = 0; |
| CvSeqBlock block; |
| |
| if( CV_IS_SEQ( array )) |
| { |
| contour = (CvSeq*)array; |
| if( !CV_IS_SEQ_POLYLINE( contour )) |
| CV_ERROR( CV_StsBadArg, "Unsupported sequence type" ); |
| } |
| else |
| { |
| CV_CALL( contour = cvPointSeqFromMat( |
| CV_SEQ_KIND_CURVE, array, &contour_header, &block )); |
| } |
| |
| if( cvSliceLength( slice, contour ) == contour->total ) |
| { |
| IPPI_CALL( icvContourArea( contour, &area )); |
| } |
| else |
| { |
| if( CV_SEQ_ELTYPE( contour ) != CV_32SC2 ) |
| CV_ERROR( CV_StsUnsupportedFormat, |
| "Only curves with integer coordinates are supported in case of contour slice" ); |
| IPPI_CALL( icvContourSecArea( contour, slice, &area )); |
| } |
| |
| __END__; |
| |
| return area; |
| } |
| |
| |
| /* for now this function works bad with singular cases |
| You can see in the code, that when some troubles with |
| matrices or some variables occur - |
| box filled with zero values is returned. |
| However in general function works fine. |
| */ |
| static void |
| icvFitEllipse_F( CvSeq* points, CvBox2D* box ) |
| { |
| CvMat* D = 0; |
| |
| CV_FUNCNAME( "icvFitEllipse_F" ); |
| |
| __BEGIN__; |
| |
| double S[36], C[36], T[36]; |
| |
| int i, j; |
| double eigenvalues[6], eigenvectors[36]; |
| double a, b, c, d, e, f; |
| double x0, y0, idet, scale, offx = 0, offy = 0; |
| |
| int n = points->total; |
| CvSeqReader reader; |
| int is_float = CV_SEQ_ELTYPE(points) == CV_32FC2; |
| |
| CvMat _S = cvMat(6,6,CV_64F,S), _C = cvMat(6,6,CV_64F,C), _T = cvMat(6,6,CV_64F,T); |
| CvMat _EIGVECS = cvMat(6,6,CV_64F,eigenvectors), _EIGVALS = cvMat(6,1,CV_64F,eigenvalues); |
| |
| /* create matrix D of input points */ |
| CV_CALL( D = cvCreateMat( n, 6, CV_64F )); |
| |
| cvStartReadSeq( points, &reader ); |
| |
| /* shift all points to zero */ |
| for( i = 0; i < n; i++ ) |
| { |
| if( !is_float ) |
| { |
| offx += ((CvPoint*)reader.ptr)->x; |
| offy += ((CvPoint*)reader.ptr)->y; |
| } |
| else |
| { |
| offx += ((CvPoint2D32f*)reader.ptr)->x; |
| offy += ((CvPoint2D32f*)reader.ptr)->y; |
| } |
| CV_NEXT_SEQ_ELEM( points->elem_size, reader ); |
| } |
| |
| offx /= n; |
| offy /= n; |
| |
| // fill matrix rows as (x*x, x*y, y*y, x, y, 1 ) |
| for( i = 0; i < n; i++ ) |
| { |
| double x, y; |
| double* Dptr = D->data.db + i*6; |
| |
| if( !is_float ) |
| { |
| x = ((CvPoint*)reader.ptr)->x - offx; |
| y = ((CvPoint*)reader.ptr)->y - offy; |
| } |
| else |
| { |
| x = ((CvPoint2D32f*)reader.ptr)->x - offx; |
| y = ((CvPoint2D32f*)reader.ptr)->y - offy; |
| } |
| CV_NEXT_SEQ_ELEM( points->elem_size, reader ); |
| |
| Dptr[0] = x * x; |
| Dptr[1] = x * y; |
| Dptr[2] = y * y; |
| Dptr[3] = x; |
| Dptr[4] = y; |
| Dptr[5] = 1.; |
| } |
| |
| // S = D^t*D |
| cvMulTransposed( D, &_S, 1 ); |
| cvSVD( &_S, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ); |
| |
| for( i = 0; i < 6; i++ ) |
| { |
| double a = eigenvalues[i]; |
| a = a < DBL_EPSILON ? 0 : 1./sqrt(sqrt(a)); |
| for( j = 0; j < 6; j++ ) |
| eigenvectors[i*6 + j] *= a; |
| } |
| |
| // C = Q^-1 = transp(INVEIGV) * INVEIGV |
| cvMulTransposed( &_EIGVECS, &_C, 1 ); |
| |
| cvZero( &_S ); |
| S[2] = 2.; |
| S[7] = -1.; |
| S[12] = 2.; |
| |
| // S = Q^-1*S*Q^-1 |
| cvMatMul( &_C, &_S, &_T ); |
| cvMatMul( &_T, &_C, &_S ); |
| |
| // and find its eigenvalues and vectors too |
| //cvSVD( &_S, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ); |
| cvEigenVV( &_S, &_EIGVECS, &_EIGVALS, 0 ); |
| |
| for( i = 0; i < 3; i++ ) |
| if( eigenvalues[i] > 0 ) |
| break; |
| |
| if( i >= 3 /*eigenvalues[0] < DBL_EPSILON*/ ) |
| { |
| box->center.x = box->center.y = |
| box->size.width = box->size.height = |
| box->angle = 0.f; |
| EXIT; |
| } |
| |
| // now find truthful eigenvector |
| _EIGVECS = cvMat( 6, 1, CV_64F, eigenvectors + 6*i ); |
| _T = cvMat( 6, 1, CV_64F, T ); |
| // Q^-1*eigenvecs[0] |
| cvMatMul( &_C, &_EIGVECS, &_T ); |
| |
| // extract vector components |
| a = T[0]; b = T[1]; c = T[2]; d = T[3]; e = T[4]; f = T[5]; |
| |
| ///////////////// extract ellipse axes from above values //////////////// |
| |
| /* |
| 1) find center of ellipse |
| it satisfy equation |
| | a b/2 | * | x0 | + | d/2 | = |0 | |
| | b/2 c | | y0 | | e/2 | |0 | |
| |
| */ |
| idet = a * c - b * b * 0.25; |
| idet = idet > DBL_EPSILON ? 1./idet : 0; |
| |
| // we must normalize (a b c d e f ) to fit (4ac-b^2=1) |
| scale = sqrt( 0.25 * idet ); |
| |
| if( scale < DBL_EPSILON ) |
| { |
| box->center.x = (float)offx; |
| box->center.y = (float)offy; |
| box->size.width = box->size.height = box->angle = 0.f; |
| EXIT; |
| } |
| |
| a *= scale; |
| b *= scale; |
| c *= scale; |
| d *= scale; |
| e *= scale; |
| f *= scale; |
| |
| x0 = (-d * c + e * b * 0.5) * 2.; |
| y0 = (-a * e + d * b * 0.5) * 2.; |
| |
| // recover center |
| box->center.x = (float)(x0 + offx); |
| box->center.y = (float)(y0 + offy); |
| |
| // offset ellipse to (x0,y0) |
| // new f == F(x0,y0) |
| f += a * x0 * x0 + b * x0 * y0 + c * y0 * y0 + d * x0 + e * y0; |
| |
| if( fabs(f) < DBL_EPSILON ) |
| { |
| box->size.width = box->size.height = box->angle = 0.f; |
| EXIT; |
| } |
| |
| scale = -1. / f; |
| // normalize to f = 1 |
| a *= scale; |
| b *= scale; |
| c *= scale; |
| |
| // extract axis of ellipse |
| // one more eigenvalue operation |
| S[0] = a; |
| S[1] = S[2] = b * 0.5; |
| S[3] = c; |
| |
| _S = cvMat( 2, 2, CV_64F, S ); |
| _EIGVECS = cvMat( 2, 2, CV_64F, eigenvectors ); |
| _EIGVALS = cvMat( 1, 2, CV_64F, eigenvalues ); |
| cvSVD( &_S, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T ); |
| |
| // exteract axis length from eigenvectors |
| box->size.width = (float)(2./sqrt(eigenvalues[0])); |
| box->size.height = (float)(2./sqrt(eigenvalues[1])); |
| |
| // calc angle |
| box->angle = (float)(180 - atan2(eigenvectors[2], eigenvectors[3])*180/CV_PI); |
| |
| __END__; |
| |
| cvReleaseMat( &D ); |
| } |
| |
| |
| CV_IMPL CvBox2D |
| cvFitEllipse2( const CvArr* array ) |
| { |
| CvBox2D box; |
| double* Ad = 0, *bd = 0; |
| |
| CV_FUNCNAME( "cvFitEllipse2" ); |
| |
| memset( &box, 0, sizeof(box)); |
| |
| __BEGIN__; |
| |
| CvContour contour_header; |
| CvSeq* ptseq = 0; |
| CvSeqBlock block; |
| int n; |
| |
| if( CV_IS_SEQ( array )) |
| { |
| ptseq = (CvSeq*)array; |
| if( !CV_IS_SEQ_POINT_SET( ptseq )) |
| CV_ERROR( CV_StsBadArg, "Unsupported sequence type" ); |
| } |
| else |
| { |
| CV_CALL( ptseq = cvPointSeqFromMat( |
| CV_SEQ_KIND_GENERIC, array, &contour_header, &block )); |
| } |
| |
| n = ptseq->total; |
| if( n < 5 ) |
| CV_ERROR( CV_StsBadSize, "Number of points should be >= 6" ); |
| #if 1 |
| icvFitEllipse_F( ptseq, &box ); |
| #else |
| /* |
| * New fitellipse algorithm, contributed by Dr. Daniel Weiss |
| */ |
| { |
| double gfp[5], rp[5], t; |
| CvMat A, b, x; |
| const double min_eps = 1e-6; |
| int i, is_float; |
| CvSeqReader reader; |
| |
| CV_CALL( Ad = (double*)cvAlloc( n*5*sizeof(Ad[0]) )); |
| CV_CALL( bd = (double*)cvAlloc( n*sizeof(bd[0]) )); |
| |
| // first fit for parameters A - E |
| A = cvMat( n, 5, CV_64F, Ad ); |
| b = cvMat( n, 1, CV_64F, bd ); |
| x = cvMat( 5, 1, CV_64F, gfp ); |
| |
| cvStartReadSeq( ptseq, &reader ); |
| is_float = CV_SEQ_ELTYPE(ptseq) == CV_32FC2; |
| |
| for( i = 0; i < n; i++ ) |
| { |
| CvPoint2D32f p; |
| if( is_float ) |
| p = *(CvPoint2D32f*)(reader.ptr); |
| else |
| { |
| p.x = (float)((int*)reader.ptr)[0]; |
| p.y = (float)((int*)reader.ptr)[1]; |
| } |
| CV_NEXT_SEQ_ELEM( sizeof(p), reader ); |
| |
| bd[i] = 10000.0; // 1.0? |
| Ad[i*5] = -(double)p.x * p.x; // A - C signs inverted as proposed by APP |
| Ad[i*5 + 1] = -(double)p.y * p.y; |
| Ad[i*5 + 2] = -(double)p.x * p.y; |
| Ad[i*5 + 3] = p.x; |
| Ad[i*5 + 4] = p.y; |
| } |
| |
| cvSolve( &A, &b, &x, CV_SVD ); |
| |
| // now use general-form parameters A - E to find the ellipse center: |
| // differentiate general form wrt x/y to get two equations for cx and cy |
| A = cvMat( 2, 2, CV_64F, Ad ); |
| b = cvMat( 2, 1, CV_64F, bd ); |
| x = cvMat( 2, 1, CV_64F, rp ); |
| Ad[0] = 2 * gfp[0]; |
| Ad[1] = Ad[2] = gfp[2]; |
| Ad[3] = 2 * gfp[1]; |
| bd[0] = gfp[3]; |
| bd[1] = gfp[4]; |
| cvSolve( &A, &b, &x, CV_SVD ); |
| |
| // re-fit for parameters A - C with those center coordinates |
| A = cvMat( n, 3, CV_64F, Ad ); |
| b = cvMat( n, 1, CV_64F, bd ); |
| x = cvMat( 3, 1, CV_64F, gfp ); |
| for( i = 0; i < n; i++ ) |
| { |
| CvPoint2D32f p; |
| if( is_float ) |
| p = *(CvPoint2D32f*)(reader.ptr); |
| else |
| { |
| p.x = (float)((int*)reader.ptr)[0]; |
| p.y = (float)((int*)reader.ptr)[1]; |
| } |
| CV_NEXT_SEQ_ELEM( sizeof(p), reader ); |
| bd[i] = 1.0; |
| Ad[i * 3] = (p.x - rp[0]) * (p.x - rp[0]); |
| Ad[i * 3 + 1] = (p.y - rp[1]) * (p.y - rp[1]); |
| Ad[i * 3 + 2] = (p.x - rp[0]) * (p.y - rp[1]); |
| } |
| cvSolve(&A, &b, &x, CV_SVD); |
| |
| // store angle and radii |
| rp[4] = -0.5 * atan2(gfp[2], gfp[1] - gfp[0]); // convert from APP angle usage |
| t = sin(-2.0 * rp[4]); |
| if( fabs(t) > fabs(gfp[2])*min_eps ) |
| t = gfp[2]/t; |
| else |
| t = gfp[1] - gfp[0]; |
| rp[2] = fabs(gfp[0] + gfp[1] - t); |
| if( rp[2] > min_eps ) |
| rp[2] = sqrt(2.0 / rp[2]); |
| rp[3] = fabs(gfp[0] + gfp[1] + t); |
| if( rp[3] > min_eps ) |
| rp[3] = sqrt(2.0 / rp[3]); |
| |
| box.center.x = (float)rp[0]; |
| box.center.y = (float)rp[1]; |
| box.size.width = (float)(rp[2]*2); |
| box.size.height = (float)(rp[3]*2); |
| if( box.size.width > box.size.height ) |
| { |
| float tmp; |
| CV_SWAP( box.size.width, box.size.height, tmp ); |
| box.angle = (float)(90 + rp[4]*180/CV_PI); |
| } |
| if( box.angle < -180 ) |
| box.angle += 360; |
| if( box.angle > 360 ) |
| box.angle -= 360; |
| } |
| #endif |
| __END__; |
| |
| cvFree( &Ad ); |
| cvFree( &bd ); |
| |
| return box; |
| } |
| |
| |
| /* Calculates bounding rectagnle of a point set or retrieves already calculated */ |
| CV_IMPL CvRect |
| cvBoundingRect( CvArr* array, int update ) |
| { |
| CvSeqReader reader; |
| CvRect rect = { 0, 0, 0, 0 }; |
| CvContour contour_header; |
| CvSeq* ptseq = 0; |
| CvSeqBlock block; |
| |
| CV_FUNCNAME( "cvBoundingRect" ); |
| |
| __BEGIN__; |
| |
| CvMat stub, *mat = 0; |
| int xmin = 0, ymin = 0, xmax = -1, ymax = -1, i, j, k; |
| int calculate = update; |
| |
| if( CV_IS_SEQ( array )) |
| { |
| ptseq = (CvSeq*)array; |
| if( !CV_IS_SEQ_POINT_SET( ptseq )) |
| CV_ERROR( CV_StsBadArg, "Unsupported sequence type" ); |
| |
| if( ptseq->header_size < (int)sizeof(CvContour)) |
| { |
| /*if( update == 1 ) |
| CV_ERROR( CV_StsBadArg, "The header is too small to fit the rectangle, " |
| "so it could not be updated" );*/ |
| update = 0; |
| calculate = 1; |
| } |
| } |
| else |
| { |
| CV_CALL( mat = cvGetMat( array, &stub )); |
| if( CV_MAT_TYPE(mat->type) == CV_32SC2 || |
| CV_MAT_TYPE(mat->type) == CV_32FC2 ) |
| { |
| CV_CALL( ptseq = cvPointSeqFromMat( |
| CV_SEQ_KIND_GENERIC, mat, &contour_header, &block )); |
| mat = 0; |
| } |
| else if( CV_MAT_TYPE(mat->type) != CV_8UC1 && |
| CV_MAT_TYPE(mat->type) != CV_8SC1 ) |
| CV_ERROR( CV_StsUnsupportedFormat, |
| "The image/matrix format is not supported by the function" ); |
| update = 0; |
| calculate = 1; |
| } |
| |
| if( !calculate ) |
| { |
| rect = ((CvContour*)ptseq)->rect; |
| EXIT; |
| } |
| |
| if( mat ) |
| { |
| CvSize size = cvGetMatSize(mat); |
| xmin = size.width; |
| ymin = -1; |
| |
| for( i = 0; i < size.height; i++ ) |
| { |
| uchar* _ptr = mat->data.ptr + i*mat->step; |
| uchar* ptr = (uchar*)cvAlignPtr(_ptr, 4); |
| int have_nz = 0, k_min, offset = (int)(ptr - _ptr); |
| j = 0; |
| offset = MIN(offset, size.width); |
| for( ; j < offset; j++ ) |
| if( _ptr[j] ) |
| { |
| have_nz = 1; |
| break; |
| } |
| if( j < offset ) |
| { |
| if( j < xmin ) |
| xmin = j; |
| if( j > xmax ) |
| xmax = j; |
| } |
| if( offset < size.width ) |
| { |
| xmin -= offset; |
| xmax -= offset; |
| size.width -= offset; |
| j = 0; |
| for( ; j <= xmin - 4; j += 4 ) |
| if( *((int*)(ptr+j)) ) |
| break; |
| for( ; j < xmin; j++ ) |
| if( ptr[j] ) |
| { |
| xmin = j; |
| if( j > xmax ) |
| xmax = j; |
| have_nz = 1; |
| break; |
| } |
| k_min = MAX(j-1, xmax); |
| k = size.width - 1; |
| for( ; k > k_min && (k&3) != 3; k-- ) |
| if( ptr[k] ) |
| break; |
| if( k > k_min && (k&3) == 3 ) |
| { |
| for( ; k > k_min+3; k -= 4 ) |
| if( *((int*)(ptr+k-3)) ) |
| break; |
| } |
| for( ; k > k_min; k-- ) |
| if( ptr[k] ) |
| { |
| xmax = k; |
| have_nz = 1; |
| break; |
| } |
| if( !have_nz ) |
| { |
| j &= ~3; |
| for( ; j <= k - 3; j += 4 ) |
| if( *((int*)(ptr+j)) ) |
| break; |
| for( ; j <= k; j++ ) |
| if( ptr[j] ) |
| { |
| have_nz = 1; |
| break; |
| } |
| } |
| xmin += offset; |
| xmax += offset; |
| size.width += offset; |
| } |
| if( have_nz ) |
| { |
| if( ymin < 0 ) |
| ymin = i; |
| ymax = i; |
| } |
| } |
| |
| if( xmin >= size.width ) |
| xmin = ymin = 0; |
| } |
| else if( ptseq->total ) |
| { |
| int is_float = CV_SEQ_ELTYPE(ptseq) == CV_32FC2; |
| cvStartReadSeq( ptseq, &reader, 0 ); |
| |
| if( !is_float ) |
| { |
| CvPoint pt; |
| /* init values */ |
| CV_READ_SEQ_ELEM( pt, reader ); |
| xmin = xmax = pt.x; |
| ymin = ymax = pt.y; |
| |
| for( i = 1; i < ptseq->total; i++ ) |
| { |
| CV_READ_SEQ_ELEM( pt, reader ); |
| |
| if( xmin > pt.x ) |
| xmin = pt.x; |
| |
| if( xmax < pt.x ) |
| xmax = pt.x; |
| |
| if( ymin > pt.y ) |
| ymin = pt.y; |
| |
| if( ymax < pt.y ) |
| ymax = pt.y; |
| } |
| } |
| else |
| { |
| CvPoint pt; |
| Cv32suf v; |
| /* init values */ |
| CV_READ_SEQ_ELEM( pt, reader ); |
| xmin = xmax = CV_TOGGLE_FLT(pt.x); |
| ymin = ymax = CV_TOGGLE_FLT(pt.y); |
| |
| for( i = 1; i < ptseq->total; i++ ) |
| { |
| CV_READ_SEQ_ELEM( pt, reader ); |
| pt.x = CV_TOGGLE_FLT(pt.x); |
| pt.y = CV_TOGGLE_FLT(pt.y); |
| |
| if( xmin > pt.x ) |
| xmin = pt.x; |
| |
| if( xmax < pt.x ) |
| xmax = pt.x; |
| |
| if( ymin > pt.y ) |
| ymin = pt.y; |
| |
| if( ymax < pt.y ) |
| ymax = pt.y; |
| } |
| |
| v.i = CV_TOGGLE_FLT(xmin); xmin = cvFloor(v.f); |
| v.i = CV_TOGGLE_FLT(ymin); ymin = cvFloor(v.f); |
| /* because right and bottom sides of |
| the bounding rectangle are not inclusive |
| (note +1 in width and height calculation below), |
| cvFloor is used here instead of cvCeil */ |
| v.i = CV_TOGGLE_FLT(xmax); xmax = cvFloor(v.f); |
| v.i = CV_TOGGLE_FLT(ymax); ymax = cvFloor(v.f); |
| } |
| } |
| |
| rect.x = xmin; |
| rect.y = ymin; |
| rect.width = xmax - xmin + 1; |
| rect.height = ymax - ymin + 1; |
| |
| if( update ) |
| ((CvContour*)ptseq)->rect = rect; |
| |
| __END__; |
| |
| return rect; |
| } |
| |
| |
| /* End of file. */ |
