diff --git a/modules/core/include/opencv2/core/core.hpp b/modules/core/include/opencv2/core/core.hpp index 94d7526aa..afe769081 100644 --- a/modules/core/include/opencv2/core/core.hpp +++ b/modules/core/include/opencv2/core/core.hpp @@ -109,6 +109,8 @@ template class CV_EXPORTS MatIterator_; template class CV_EXPORTS MatConstIterator_; template class CV_EXPORTS MatCommaInitializer_; +template class CV_EXPORTS AutoBuffer; + CV_EXPORTS string format( const char* fmt, ... ); CV_EXPORTS string tempfile( const char* suffix CV_DEFAULT(0)); @@ -2061,7 +2063,8 @@ CV_EXPORTS void swap(Mat& a, Mat& b); //! converts array (CvMat or IplImage) to cv::Mat CV_EXPORTS Mat cvarrToMat(const CvArr* arr, bool copyData=false, - bool allowND=true, int coiMode=0); + bool allowND=true, int coiMode=0, + AutoBuffer* buf=0); //! extracts Channel of Interest from CvMat or IplImage and makes cv::Mat out of it. CV_EXPORTS void extractImageCOI(const CvArr* arr, OutputArray coiimg, int coi=-1); //! inserts single-channel cv::Mat into a multi-channel CvMat or IplImage @@ -3081,7 +3084,7 @@ public: \code void my_func(const cv::Mat& m) { - cv::AutoBuffer buf; // create automatic buffer containing 1000 floats + cv::AutoBuffer buf; // create automatic buffer containing 1000 floats buf.allocate(m.rows); // if m.rows <= 1000, the pre-allocated buffer is used, // otherwise the buffer of "m.rows" floats will be allocated @@ -3090,16 +3093,21 @@ public: } \endcode */ -template class CV_EXPORTS AutoBuffer +template class CV_EXPORTS AutoBuffer { public: typedef _Tp value_type; - enum { buffer_padding = (int)((16 + sizeof(_Tp) - 1)/sizeof(_Tp)) }; //! the default contructor AutoBuffer(); //! constructor taking the real buffer size AutoBuffer(size_t _size); + + //! the copy constructor + AutoBuffer(const AutoBuffer<_Tp, fixed_size>& buf); + //! the assignment operator + AutoBuffer<_Tp, fixed_size>& operator = (const AutoBuffer<_Tp, fixed_size>& buf); + //! destructor. calls deallocate() ~AutoBuffer(); @@ -3107,6 +3115,10 @@ public: void allocate(size_t _size); //! deallocates the buffer if it was dynamically allocated void deallocate(); + //! resizes the buffer and preserves the content + void resize(size_t _size); + //! returns the current buffer size + size_t size() const; //! returns pointer to the real buffer, stack-allocated or head-allocated operator _Tp* (); //! returns read-only pointer to the real buffer, stack-allocated or head-allocated @@ -3116,9 +3128,9 @@ protected: //! pointer to the real buffer, can point to buf if the buffer is small enough _Tp* ptr; //! size of the real buffer - size_t size; + size_t sz; //! pre-allocated buffer - _Tp buf[fixed_size+buffer_padding]; + _Tp buf[fixed_size]; }; /////////////////////////// multi-dimensional dense matrix ////////////////////////// @@ -4314,7 +4326,6 @@ public: int index; }; - class CV_EXPORTS Algorithm; class CV_EXPORTS AlgorithmInfo; struct CV_EXPORTS AlgorithmInfoData; diff --git a/modules/core/include/opencv2/core/operations.hpp b/modules/core/include/opencv2/core/operations.hpp index eaae0c0b1..ad9921c17 100644 --- a/modules/core/include/opencv2/core/operations.hpp +++ b/modules/core/include/opencv2/core/operations.hpp @@ -2534,48 +2534,109 @@ inline Point LineIterator::pos() const /////////////////////////////// AutoBuffer //////////////////////////////////////// -template inline AutoBuffer<_Tp, fixed_size>::AutoBuffer() +template inline +AutoBuffer<_Tp, fixed_size>::AutoBuffer() { ptr = buf; - size = fixed_size; + sz = fixed_size; } -template inline AutoBuffer<_Tp, fixed_size>::AutoBuffer(size_t _size) +template inline +AutoBuffer<_Tp, fixed_size>::AutoBuffer(size_t _size) { ptr = buf; - size = fixed_size; + sz = fixed_size; allocate(_size); } -template inline AutoBuffer<_Tp, fixed_size>::~AutoBuffer() +template inline +AutoBuffer<_Tp, fixed_size>::AutoBuffer(const AutoBuffer<_Tp, fixed_size>& abuf ) +{ + ptr = buf; + sz = fixed_size; + allocate(abuf.size); + for( size_t i = 0; i < sz; i++ ) + ptr[i] = abuf.ptr[i]; +} + +template inline AutoBuffer<_Tp, fixed_size>& +AutoBuffer<_Tp, fixed_size>::operator = (const AutoBuffer<_Tp, fixed_size>& abuf) +{ + if( this != &abuf ) + { + deallocate(); + allocate(abuf.size); + for( size_t i = 0; i < sz; i++ ) + ptr[i] = abuf.ptr[i]; + } + return *this; +} + +template inline +AutoBuffer<_Tp, fixed_size>::~AutoBuffer() { deallocate(); } -template inline void AutoBuffer<_Tp, fixed_size>::allocate(size_t _size) +template inline void +AutoBuffer<_Tp, fixed_size>::allocate(size_t _size) { - if(_size <= size) + if(_size <= sz) + { + sz = _size; return; + } deallocate(); if(_size > fixed_size) { - ptr = cv::allocate<_Tp>(_size); - size = _size; + ptr = new _Tp[_size]; + sz = _size; } } -template inline void AutoBuffer<_Tp, fixed_size>::deallocate() +template inline void +AutoBuffer<_Tp, fixed_size>::deallocate() { if( ptr != buf ) { - cv::deallocate<_Tp>(ptr, size); + delete[] ptr; ptr = buf; - size = fixed_size; + sz = fixed_size; } } -template inline AutoBuffer<_Tp, fixed_size>::operator _Tp* () +template inline void +AutoBuffer<_Tp, fixed_size>::resize(size_t _size) +{ + if(_size <= sz) + { + sz = _size; + return; + } + size_t i, prevsize = sz, minsize = MIN(prevsize, _size); + _Tp* prevptr = ptr; + + ptr = _size > fixed_size ? new _Tp[_size] : buf; + sz = _size; + + if( ptr != prevptr ) + for( i = 0; i < minsize; i++ ) + ptr[i] = prevptr[i]; + for( i = prevsize; i < _size; i++ ) + ptr[i] = _Tp(); + + if( prevptr != buf ) + delete[] prevptr; +} + +template inline size_t +AutoBuffer<_Tp, fixed_size>::size() const +{ return sz; } + +template inline +AutoBuffer<_Tp, fixed_size>::operator _Tp* () { return ptr; } -template inline AutoBuffer<_Tp, fixed_size>::operator const _Tp* () const +template inline +AutoBuffer<_Tp, fixed_size>::operator const _Tp* () const { return ptr; } diff --git a/modules/core/src/convert.cpp b/modules/core/src/convert.cpp index 0ba066ec7..439d4c1a9 100644 --- a/modules/core/src/convert.cpp +++ b/modules/core/src/convert.cpp @@ -1314,7 +1314,7 @@ cvMixChannels( const CvArr** src, int src_count, CvArr** dst, int dst_count, const int* from_to, int pair_count ) { - cv::AutoBuffer buf(src_count + dst_count); + cv::AutoBuffer buf(src_count + dst_count); int i; for( i = 0; i < src_count; i++ ) diff --git a/modules/core/src/matrix.cpp b/modules/core/src/matrix.cpp index 3b80d5a3d..ebee0def9 100644 --- a/modules/core/src/matrix.cpp +++ b/modules/core/src/matrix.cpp @@ -669,7 +669,7 @@ void Mat::push_back(const Mat& elems) Mat cvarrToMat(const CvArr* arr, bool copyData, - bool /*allowND*/, int coiMode) + bool /*allowND*/, int coiMode, AutoBuffer* abuf ) { if( !arr ) return Mat(); @@ -687,10 +687,21 @@ Mat cvarrToMat(const CvArr* arr, bool copyData, if( CV_IS_SEQ(arr) ) { CvSeq* seq = (CvSeq*)arr; - CV_Assert(seq->total > 0 && CV_ELEM_SIZE(seq->flags) == seq->elem_size); + int total = seq->total, type = CV_MAT_TYPE(seq->flags), esz = seq->elem_size; + if( total == 0 ) + return Mat(); + CV_Assert(total > 0 && CV_ELEM_SIZE(seq->flags) == esz); if(!copyData && seq->first->next == seq->first) - return Mat(seq->total, 1, CV_MAT_TYPE(seq->flags), seq->first->data); - Mat buf(seq->total, 1, CV_MAT_TYPE(seq->flags)); + return Mat(total, 1, type, seq->first->data); + if( abuf ) + { + abuf->allocate(((size_t)total*esz + sizeof(double)-1)/sizeof(double)); + double* bufdata = *abuf; + cvCvtSeqToArray(seq, bufdata, CV_WHOLE_SEQ); + return Mat(total, 1, type, bufdata); + } + + Mat buf(total, 1, type); cvCvtSeqToArray(seq, buf.data, CV_WHOLE_SEQ); return buf; } @@ -830,7 +841,7 @@ int Mat::checkVector(int _elemChannels, int _depth, bool _requireContinuous) con { return (depth() == _depth || _depth <= 0) && (isContinuous() || !_requireContinuous) && - ((dims == 2 && (((rows == 1 || cols == 1) && channels() == _elemChannels) || (cols == _elemChannels))) || + ((dims == 2 && (((rows == 1 || cols == 1) && channels() == _elemChannels) || (cols == _elemChannels && channels() == 1))) || (dims == 3 && channels() == 1 && size.p[2] == _elemChannels && (size.p[0] == 1 || size.p[1] == 1) && (isContinuous() || step.p[1] == step.p[2]*size.p[2]))) ? (int)(total()*channels()/_elemChannels) : -1; diff --git a/modules/highgui/src/grfmt_pxm.cpp b/modules/highgui/src/grfmt_pxm.cpp index f73bbb1bf..57b222eb3 100644 --- a/modules/highgui/src/grfmt_pxm.cpp +++ b/modules/highgui/src/grfmt_pxm.cpp @@ -202,9 +202,9 @@ bool PxMDecoder::readData( Mat& img ) if( m_offset < 0 || !m_strm.isOpened()) return false; - AutoBuffer _src(src_pitch + 32); + AutoBuffer _src(src_pitch + 32); uchar* src = _src; - AutoBuffer _gray_palette; + AutoBuffer _gray_palette; uchar* gray_palette = _gray_palette; // create LUT for converting colors diff --git a/modules/highgui/src/grfmt_tiff.cpp b/modules/highgui/src/grfmt_tiff.cpp index d2321ceb5..7da9320fc 100644 --- a/modules/highgui/src/grfmt_tiff.cpp +++ b/modules/highgui/src/grfmt_tiff.cpp @@ -469,7 +469,7 @@ bool TiffEncoder::writeLibTiff( const Mat& img, const vector& /*params*/) // row buffer, because TIFFWriteScanline modifies the original data! size_t scanlineSize = TIFFScanlineSize(pTiffHandle); - AutoBuffer _buffer(scanlineSize+32); + AutoBuffer _buffer(scanlineSize+32); uchar* buffer = _buffer; if (!buffer) { @@ -577,9 +577,9 @@ bool TiffEncoder::write( const Mat& img, const vector& /*params*/) #endif*/ int directoryOffset = 0; - AutoBuffer stripOffsets(stripCount); - AutoBuffer stripCounts(stripCount); - AutoBuffer _buffer(fileStep+32); + AutoBuffer stripOffsets(stripCount); + AutoBuffer stripCounts(stripCount); + AutoBuffer _buffer(fileStep+32); uchar* buffer = _buffer; int stripOffsetsOffset = 0; int stripCountsOffset = 0; diff --git a/modules/imgproc/src/_geom.h b/modules/imgproc/src/_geom.h index 25050390e..bfde73d6d 100644 --- a/modules/imgproc/src/_geom.h +++ b/modules/imgproc/src/_geom.h @@ -52,16 +52,6 @@ CV_INLINE float icvDistanceL2_32f( CvPoint2D32f pt1, CvPoint2D32f pt2 ) } -int icvIntersectLines( double x1, double dx1, double y1, double dy1, - double x2, double dx2, double y2, double dy2, - double* t2 ); - - -void icvIntersectLines3( double* a0, double* b0, double* c0, - double* a1, double* b1, double* c1, - CvPoint2D32f* point ); - - /* curvature: 0 - 1-curvature, 1 - k-cosine curvature. */ CvSeq* icvApproximateChainTC89( CvChain* chain, int header_size, CvMemStorage* storage, int method ); diff --git a/modules/imgproc/src/approx.cpp b/modules/imgproc/src/approx.cpp index c2831ca70..dc9bee2f4 100644 --- a/modules/imgproc/src/approx.cpp +++ b/modules/imgproc/src/approx.cpp @@ -469,48 +469,59 @@ cvApproxChains( CvSeq* src_seq, /* Ramer-Douglas-Peucker algorithm for polygon simplification */ -/* the version for integer point coordinates */ -template static CvSeq* -icvApproxPolyDP( CvSeq* src_contour, int header_size, - CvMemStorage* storage, double eps ) +namespace cv { + +template static int +approxPolyDP_( const Point_* src_contour, int count0, Point_* dst_contour, + bool is_closed0, double eps, AutoBuffer* _stack ) +{ + #define PUSH_SLICE(slice) \ + if( top >= stacksz ) \ + { \ + _stack->resize(stacksz*3/2); \ + stack = *_stack; \ + stacksz = _stack->size(); \ + } \ + stack[top++] = slice + + #define READ_PT(pt, pos) \ + pt = src_contour[pos]; \ + if( ++pos >= count ) pos = 0 + + #define READ_DST_PT(pt, pos) \ + pt = dst_contour[pos]; \ + if( ++pos >= count ) pos = 0 + + #define WRITE_PT(pt) \ + dst_contour[new_count++] = pt + typedef cv::Point_ PT; int init_iters = 3; - CvSlice slice = {0, 0}, right_slice = {0, 0}; - CvSeqReader reader, reader2; - CvSeqWriter writer; - PT start_pt(-1000000, -1000000), end_pt(0, 0), pt(0,0); - int i = 0, j, count = src_contour->total, new_count; - int is_closed = CV_IS_SEQ_CLOSED( src_contour ); + Range slice(0, 0), right_slice(0, 0); + PT start_pt((T)-1000000, (T)-1000000), end_pt(0, 0), pt(0,0); + int i = 0, j, pos = 0, wpos, count = count0, new_count=0; + int is_closed = is_closed0; bool le_eps = false; - CvMemStorage* temp_storage = 0; - CvSeq* stack = 0; - CvSeq* dst_contour; + size_t top = 0, stacksz = _stack->size(); + Range* stack = *_stack; - assert( CV_SEQ_ELTYPE(src_contour) == cv::DataType::type ); - cvStartWriteSeq( src_contour->flags, header_size, sizeof(pt), storage, &writer ); + if( count == 0 ) + return 0; - if( src_contour->total == 0 ) - return cvEndWriteSeq( &writer ); - - temp_storage = cvCreateChildMemStorage( storage ); - - assert( src_contour->first != 0 ); - stack = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvSlice), temp_storage ); eps *= eps; - cvStartReadSeq( src_contour, &reader, 0 ); if( !is_closed ) { - right_slice.start_index = count; - end_pt = *(PT*)(reader.ptr); - start_pt = *(PT*)cvGetSeqElem( src_contour, -1 ); + right_slice.start = count; + end_pt = src_contour[0]; + start_pt = src_contour[count-1]; if( start_pt.x != end_pt.x || start_pt.y != end_pt.y ) { - slice.start_index = 0; - slice.end_index = count - 1; - cvSeqPush( stack, &slice ); + slice.start = 0; + slice.end = count - 1; + PUSH_SLICE(slice); } else { @@ -521,20 +532,20 @@ icvApproxPolyDP( CvSeq* src_contour, int header_size, if( is_closed ) { - /* 1. Find approximately two farthest points of the contour */ - right_slice.start_index = 0; + // 1. Find approximately two farthest points of the contour + right_slice.start = 0; for( i = 0; i < init_iters; i++ ) { double dist, max_dist = 0; - cvSetSeqReaderPos( &reader, right_slice.start_index, 1 ); - CV_READ_SEQ_ELEM( start_pt, reader ); /* read the first point */ + pos = (pos + right_slice.start) % count; + READ_PT(start_pt, pos); for( j = 1; j < count; j++ ) { double dx, dy; - CV_READ_SEQ_ELEM( pt, reader ); + READ_PT(pt, pos); dx = pt.x - start_pt.x; dy = pt.y - start_pt.y; @@ -543,43 +554,35 @@ icvApproxPolyDP( CvSeq* src_contour, int header_size, if( dist > max_dist ) { max_dist = dist; - right_slice.start_index = j; + right_slice.start = j; } } le_eps = max_dist <= eps; } - /* 2. initialize the stack */ + // 2. initialize the stack if( !le_eps ) { - slice.start_index = cvGetSeqReaderPos( &reader ); - slice.end_index = right_slice.start_index += slice.start_index; + right_slice.end = slice.start = pos % count; + slice.end = right_slice.start = (right_slice.start + slice.start) % count; - right_slice.start_index -= right_slice.start_index >= count ? count : 0; - right_slice.end_index = slice.start_index; - if( right_slice.end_index < right_slice.start_index ) - right_slice.end_index += count; - - cvSeqPush( stack, &right_slice ); - cvSeqPush( stack, &slice ); + PUSH_SLICE(right_slice); + PUSH_SLICE(slice); } else - CV_WRITE_SEQ_ELEM( start_pt, writer ); + WRITE_PT(start_pt); } - /* 3. run recursive process */ - while( stack->total != 0 ) + // 3. run recursive process + while( top > 0 ) { - cvSeqPop( stack, &slice ); + slice = stack[--top]; + end_pt = src_contour[slice.end]; + pos = slice.start; + READ_PT(start_pt, pos); - cvSetSeqReaderPos( &reader, slice.end_index ); - CV_READ_SEQ_ELEM( end_pt, reader ); - - cvSetSeqReaderPos( &reader, slice.start_index ); - CV_READ_SEQ_ELEM( start_pt, reader ); - - if( slice.end_index > slice.start_index + 1 ) + if( pos != slice.end ) { double dx, dy, dist, max_dist = 0; @@ -588,15 +591,15 @@ icvApproxPolyDP( CvSeq* src_contour, int header_size, assert( dx != 0 || dy != 0 ); - for( i = slice.start_index + 1; i < slice.end_index; i++ ) + while( pos != slice.end ) { - CV_READ_SEQ_ELEM( pt, reader ); + READ_PT(pt, pos); dist = fabs((pt.y - start_pt.y) * dx - (pt.x - start_pt.x) * dy); if( dist > max_dist ) { max_dist = dist; - right_slice.start_index = i; + right_slice.start = (pos+count-1)%count; } } @@ -604,84 +607,106 @@ icvApproxPolyDP( CvSeq* src_contour, int header_size, } else { - assert( slice.end_index > slice.start_index ); le_eps = true; - /* read starting point */ - cvSetSeqReaderPos( &reader, slice.start_index ); - CV_READ_SEQ_ELEM( start_pt, reader ); + // read starting point + start_pt = src_contour[slice.start]; } if( le_eps ) { - CV_WRITE_SEQ_ELEM( start_pt, writer ); + WRITE_PT(start_pt); } else { - right_slice.end_index = slice.end_index; - slice.end_index = right_slice.start_index; - cvSeqPush( stack, &right_slice ); - cvSeqPush( stack, &slice ); + right_slice.end = slice.end; + slice.end = right_slice.start; + PUSH_SLICE(right_slice); + PUSH_SLICE(slice); } } - is_closed = CV_IS_SEQ_CLOSED( src_contour ); if( !is_closed ) - CV_WRITE_SEQ_ELEM( end_pt, writer ); - - dst_contour = cvEndWriteSeq( &writer ); + WRITE_PT( src_contour[count-1] ); // last stage: do final clean-up of the approximated contour - // remove extra points on the [almost] stright lines. + is_closed = is_closed0; + count = new_count; + pos = is_closed ? count - 1 : 0; + READ_DST_PT(start_pt, pos); + wpos = pos; + READ_DST_PT(pt, pos); - cvStartReadSeq( dst_contour, &reader, is_closed ); - CV_READ_SEQ_ELEM( start_pt, reader ); - - reader2 = reader; - CV_READ_SEQ_ELEM( pt, reader ); - - new_count = count = dst_contour->total; for( i = !is_closed; i < count - !is_closed && new_count > 2; i++ ) { double dx, dy, dist, successive_inner_product; - CV_READ_SEQ_ELEM( end_pt, reader ); + READ_DST_PT( end_pt, pos ); dx = end_pt.x - start_pt.x; dy = end_pt.y - start_pt.y; dist = fabs((pt.x - start_pt.x)*dy - (pt.y - start_pt.y)*dx); successive_inner_product = (pt.x - start_pt.x) * (end_pt.x - pt.x) + - (pt.y - start_pt.y) * (end_pt.y - pt.y); + (pt.y - start_pt.y) * (end_pt.y - pt.y); if( dist * dist <= 0.5*eps*(dx*dx + dy*dy) && dx != 0 && dy != 0 && - successive_inner_product >= 0 ) + successive_inner_product >= 0 ) { new_count--; - *((PT*)reader2.ptr) = start_pt = end_pt; - CV_NEXT_SEQ_ELEM( sizeof(pt), reader2 ); - CV_READ_SEQ_ELEM( pt, reader ); + dst_contour[wpos] = start_pt = end_pt; + if(++wpos >= count) wpos = 0; + READ_DST_PT(pt, pos); i++; continue; } - *((PT*)reader2.ptr) = start_pt = pt; - CV_NEXT_SEQ_ELEM( sizeof(pt), reader2 ); + dst_contour[wpos] = start_pt = pt; + if(++wpos >= count) wpos = 0; pt = end_pt; } if( !is_closed ) - *((PT*)reader2.ptr) = pt; + dst_contour[wpos] = pt; - if( new_count < count ) - cvSeqPopMulti( dst_contour, 0, count - new_count ); + return new_count; +} - cvReleaseMemStorage( &temp_storage ); - return dst_contour; +} + +void cv::approxPolyDP( InputArray _curve, OutputArray _approxCurve, + double epsilon, bool closed ) +{ + Mat curve = _curve.getMat(); + int npoints = curve.checkVector(2), depth = curve.depth(); + CV_Assert( npoints >= 0 && (depth == CV_32S || depth == CV_32F)); + + if( npoints == 0 ) + { + _approxCurve.release(); + return; + } + + AutoBuffer _buf(npoints); + AutoBuffer _stack(npoints); + Point* buf = _buf; + int nout = 0; + + if( depth == CV_32S ) + nout = approxPolyDP_(curve.ptr(), npoints, buf, closed, epsilon, &_stack); + else if( depth == CV_32F ) + nout = approxPolyDP_(curve.ptr(), npoints, (Point2f*)buf, closed, epsilon, &_stack); + else + CV_Error( CV_StsUnsupportedFormat, "" ); + + Mat(nout, 1, CV_MAKETYPE(depth, 2), buf).copyTo(_approxCurve); } CV_IMPL CvSeq* -cvApproxPoly( const void* array, int header_size, - CvMemStorage* storage, int method, - double parameter, int parameter2 ) +cvApproxPoly( const void* array, int header_size, + CvMemStorage* storage, int method, + double parameter, int parameter2 ) { + cv::AutoBuffer _buf; + cv::AutoBuffer stack(100); CvSeq* dst_seq = 0; CvSeq *prev_contour = 0, *parent = 0; CvContour contour_header; @@ -703,8 +728,8 @@ cvApproxPoly( const void* array, int header_size, else { src_seq = cvPointSeqFromMat( - CV_SEQ_KIND_CURVE | (parameter2 ? CV_SEQ_FLAG_CLOSED : 0), - array, &contour_header, &block ); + CV_SEQ_KIND_CURVE | (parameter2 ? CV_SEQ_FLAG_CLOSED : 0), + array, &contour_header, &block ); } if( !storage ) @@ -712,7 +737,7 @@ cvApproxPoly( const void* array, int header_size, if( header_size < 0 ) CV_Error( CV_StsOutOfRange, "header_size is negative. " - "Pass 0 to make the destination header_size == input header_size" ); + "Pass 0 to make the destination header_size == input header_size" ); if( header_size == 0 ) header_size = src_seq->header_size; @@ -722,7 +747,7 @@ cvApproxPoly( const void* array, int header_size, if( CV_IS_SEQ_CHAIN( src_seq )) { CV_Error( CV_StsBadArg, "Input curves are not polygonal. " - "Use cvApproxChains first" ); + "Use cvApproxChains first" ); } else { @@ -749,13 +774,34 @@ cvApproxPoly( const void* array, int header_size, if( parameter < 0 ) CV_Error( CV_StsOutOfRange, "Accuracy must be non-negative" ); - if( CV_SEQ_ELTYPE(src_seq) == CV_32SC2 ) - contour = icvApproxPolyDP( src_seq, header_size, storage, parameter ); + CV_Assert( CV_SEQ_ELTYPE(src_seq) == CV_32SC2 || + CV_SEQ_ELTYPE(src_seq) == CV_32FC2 ); + + { + int npoints = src_seq->total, nout = 0; + _buf.allocate(npoints*2); + cv::Point *src = _buf, *dst = src + npoints; + bool closed = CV_IS_SEQ_CLOSED(src_seq); + + if( src_seq->first->next == src_seq->first ) + src = (cv::Point*)src_seq->first->data; else - contour = icvApproxPolyDP( src_seq, header_size, storage, parameter ); + cvCvtSeqToArray(src_seq, src); + + if( CV_SEQ_ELTYPE(src_seq) == CV_32SC2 ) + nout = cv::approxPolyDP_(src, npoints, dst, closed, parameter, &stack); + else if( CV_SEQ_ELTYPE(src_seq) == CV_32FC2 ) + nout = cv::approxPolyDP_((cv::Point2f*)src, npoints, + (cv::Point2f*)dst, closed, parameter, &stack); + else + CV_Error( CV_StsUnsupportedFormat, "" ); + + contour = cvCreateSeq( src_seq->flags, header_size, + src_seq->elem_size, storage ); + cvSeqPushMulti(contour, dst, nout); + } break; default: - assert(0); CV_Error( CV_StsBadArg, "Invalid approximation method" ); } diff --git a/modules/imgproc/src/contours.cpp b/modules/imgproc/src/contours.cpp index c290385cc..75e18b2d8 100644 --- a/modules/imgproc/src/contours.cpp +++ b/modules/imgproc/src/contours.cpp @@ -1753,195 +1753,4 @@ void cv::findContours( InputOutputArray _image, OutputArrayOfArrays _contours, findContours(_image, _contours, noArray(), mode, method, offset); } -void cv::approxPolyDP( InputArray _curve, OutputArray _approxCurve, - double epsilon, bool closed ) -{ - Mat curve = _curve.getMat(); - int npoints = curve.checkVector(2), depth = curve.depth(); - CV_Assert( npoints >= 0 && (depth == CV_32S || depth == CV_32F)); - CvMat _ccurve = curve; - MemStorage storage(cvCreateMemStorage()); - CvSeq* result = cvApproxPoly(&_ccurve, sizeof(CvContour), storage, CV_POLY_APPROX_DP, epsilon, closed); - if( result->total > 0 ) - { - _approxCurve.create(result->total, 1, CV_MAKETYPE(curve.depth(), 2), -1, true); - cvCvtSeqToArray(result, _approxCurve.getMat().data ); - } -} - - -double cv::arcLength( InputArray _curve, bool closed ) -{ - Mat curve = _curve.getMat(); - CV_Assert(curve.checkVector(2) >= 0 && (curve.depth() == CV_32F || curve.depth() == CV_32S)); - CvMat _ccurve = curve; - return cvArcLength(&_ccurve, CV_WHOLE_SEQ, closed); -} - - -cv::Rect cv::boundingRect( InputArray _points ) -{ - Mat points = _points.getMat(); - CV_Assert(points.checkVector(2) >= 0 && (points.depth() == CV_32F || points.depth() == CV_32S)); - CvMat _cpoints = points; - return cvBoundingRect(&_cpoints, 0); -} - - -double cv::contourArea( InputArray _contour, bool oriented ) -{ - Mat contour = _contour.getMat(); - CV_Assert(contour.checkVector(2) >= 0 && (contour.depth() == CV_32F || contour.depth() == CV_32S)); - CvMat _ccontour = contour; - return cvContourArea(&_ccontour, CV_WHOLE_SEQ, oriented); -} - - -cv::RotatedRect cv::minAreaRect( InputArray _points ) -{ - Mat points = _points.getMat(); - CV_Assert(points.checkVector(2) >= 0 && (points.depth() == CV_32F || points.depth() == CV_32S)); - CvMat _cpoints = points; - return cvMinAreaRect2(&_cpoints, 0); -} - - -void cv::minEnclosingCircle( InputArray _points, - Point2f& center, float& radius ) -{ - Mat points = _points.getMat(); - CV_Assert(points.checkVector(2) >= 0 && (points.depth() == CV_32F || points.depth() == CV_32S)); - CvMat _cpoints = points; - cvMinEnclosingCircle( &_cpoints, (CvPoint2D32f*)¢er, &radius ); -} - - -double cv::matchShapes( InputArray _contour1, - InputArray _contour2, - int method, double parameter ) -{ - Mat contour1 = _contour1.getMat(), contour2 = _contour2.getMat(); - CV_Assert(contour1.checkVector(2) >= 0 && contour2.checkVector(2) >= 0 && - (contour1.depth() == CV_32F || contour1.depth() == CV_32S) && - contour1.depth() == contour2.depth()); - - CvMat c1 = Mat(contour1), c2 = Mat(contour2); - return cvMatchShapes(&c1, &c2, method, parameter); -} - - -void cv::convexHull( InputArray _points, OutputArray _hull, bool clockwise, bool returnPoints ) -{ - Mat points = _points.getMat(); - int nelems = points.checkVector(2), depth = points.depth(); - CV_Assert(nelems >= 0 && (depth == CV_32F || depth == CV_32S)); - - if( nelems == 0 ) - { - _hull.release(); - return; - } - - returnPoints = !_hull.fixedType() ? returnPoints : _hull.type() != CV_32S; - Mat hull(nelems, 1, returnPoints ? CV_MAKETYPE(depth, 2) : CV_32S); - CvMat _cpoints = points, _chull = hull; - cvConvexHull2(&_cpoints, &_chull, clockwise ? CV_CLOCKWISE : CV_COUNTER_CLOCKWISE, returnPoints); - _hull.create(_chull.rows, 1, hull.type(), -1, true); - Mat dhull = _hull.getMat(), shull(dhull.size(), dhull.type(), hull.data); - shull.copyTo(dhull); -} - - -void cv::convexityDefects( InputArray _points, InputArray _hull, OutputArray _defects ) -{ - Mat points = _points.getMat(); - int ptnum = points.checkVector(2, CV_32S); - CV_Assert( ptnum > 3 ); - Mat hull = _hull.getMat(); - CV_Assert( hull.checkVector(1, CV_32S) > 2 ); - Ptr storage = cvCreateMemStorage(); - - CvMat c_points = points, c_hull = hull; - CvSeq* seq = cvConvexityDefects(&c_points, &c_hull, storage); - int i, n = seq->total; - - if( n == 0 ) - { - _defects.release(); - return; - } - - _defects.create(n, 1, CV_32SC4); - Mat defects = _defects.getMat(); - - SeqIterator it = Seq(seq).begin(); - CvPoint* ptorg = (CvPoint*)points.data; - - for( i = 0; i < n; i++, ++it ) - { - CvConvexityDefect& d = *it; - int idx0 = (int)(d.start - ptorg); - int idx1 = (int)(d.end - ptorg); - int idx2 = (int)(d.depth_point - ptorg); - CV_Assert( 0 <= idx0 && idx0 < ptnum ); - CV_Assert( 0 <= idx1 && idx1 < ptnum ); - CV_Assert( 0 <= idx2 && idx2 < ptnum ); - CV_Assert( d.depth >= 0 ); - int idepth = cvRound(d.depth*256); - defects.at(i) = Vec4i(idx0, idx1, idx2, idepth); - } -} - - -bool cv::isContourConvex( InputArray _contour ) -{ - Mat contour = _contour.getMat(); - CV_Assert(contour.checkVector(2) >= 0 && - (contour.depth() == CV_32F || contour.depth() == CV_32S)); - CvMat c = Mat(contour); - return cvCheckContourConvexity(&c) > 0; -} - -cv::RotatedRect cv::fitEllipse( InputArray _points ) -{ - Mat points = _points.getMat(); - CV_Assert(points.checkVector(2) >= 0 && - (points.depth() == CV_32F || points.depth() == CV_32S)); - CvMat _cpoints = points; - return cvFitEllipse2(&_cpoints); -} - - -void cv::fitLine( InputArray _points, OutputArray _line, int distType, - double param, double reps, double aeps ) -{ - Mat points = _points.getMat(); - - bool is3d = points.checkVector(3) >= 0; - bool is2d = points.checkVector(2) >= 0; - - CV_Assert( (is2d || is3d) && (points.depth() == CV_32F || points.depth() == CV_32S) ); - CvMat _cpoints = points.reshape(2 + (int)is3d); - float line[6]; - cvFitLine(&_cpoints, distType, param, reps, aeps, &line[0]); - - int out_size = (is2d)?( (is3d)? (points.channels() * points.rows * 2) : 4 ): 6; - - _line.create(out_size, 1, CV_32F, -1, true); - Mat l = _line.getMat(); - CV_Assert( l.isContinuous() ); - memcpy( l.data, line, out_size * sizeof(line[0]) ); -} - - -double cv::pointPolygonTest( InputArray _contour, - Point2f pt, bool measureDist ) -{ - Mat contour = _contour.getMat(); - CV_Assert(contour.checkVector(2) >= 0 && - (contour.depth() == CV_32F || contour.depth() == CV_32S)); - CvMat c = Mat(contour); - return cvPointPolygonTest( &c, pt, measureDist ); -} - /* End of file. */ diff --git a/modules/imgproc/src/convhull.cpp b/modules/imgproc/src/convhull.cpp index 145c55e05..1fcff3fb4 100644 --- a/modules/imgproc/src/convhull.cpp +++ b/modules/imgproc/src/convhull.cpp @@ -40,18 +40,22 @@ //M*/ #include "precomp.hpp" +#include -static int -icvSklansky_32s( CvPoint** array, int start, int end, int* stack, int nsign, int sign2 ) +namespace cv +{ + +template +static int Sklansky_( Point_<_Tp>** array, int start, int end, int* stack, int nsign, int sign2 ) { int incr = end > start ? 1 : -1; - /* prepare first triangle */ + // prepare first triangle int pprev = start, pcur = pprev + incr, pnext = pcur + incr; int stacksize = 3; if( start == end || - (array[start]->x == array[end]->x && - array[start]->y == array[end]->y) ) + (array[start]->x == array[end]->x && + array[start]->y == array[end]->y) ) { stack[0] = start; return 1; @@ -61,94 +65,21 @@ icvSklansky_32s( CvPoint** array, int start, int end, int* stack, int nsign, int stack[1] = pcur; stack[2] = pnext; - end += incr; /* make end = afterend */ + end += incr; // make end = afterend while( pnext != end ) { - /* check the angle p1,p2,p3 */ - int cury = array[pcur]->y; - int nexty = array[pnext]->y; - int by = nexty - cury; - - if( CV_SIGN(by) != nsign ) - { - int ax = array[pcur]->x - array[pprev]->x; - int bx = array[pnext]->x - array[pcur]->x; - int ay = cury - array[pprev]->y; - int convexity = ay*bx - ax*by;/* if >0 then convex angle */ - - if( CV_SIGN(convexity) == sign2 && (ax != 0 || ay != 0) ) - { - pprev = pcur; - pcur = pnext; - pnext += incr; - stack[stacksize] = pnext; - stacksize++; - } - else - { - if( pprev == start ) - { - pcur = pnext; - stack[1] = pcur; - pnext += incr; - stack[2] = pnext; - } - else - { - stack[stacksize-2] = pnext; - pcur = pprev; - pprev = stack[stacksize-4]; - stacksize--; - } - } - } - else - { - pnext += incr; - stack[stacksize-1] = pnext; - } - } - - return --stacksize; -} - - -static int -icvSklansky_32f( CvPoint2D32f** array, int start, int end, int* stack, int nsign, int sign2 ) -{ - int incr = end > start ? 1 : -1; - /* prepare first triangle */ - int pprev = start, pcur = pprev + incr, pnext = pcur + incr; - int stacksize = 3; - - if( start == end || - (array[start]->x == array[end]->x && - array[start]->y == array[end]->y) ) - { - stack[0] = start; - return 1; - } - - stack[0] = pprev; - stack[1] = pcur; - stack[2] = pnext; - - end += incr; /* make end = afterend */ - - while( pnext != end ) - { - /* check the angle p1,p2,p3 */ - float cury = array[pcur]->y; - float nexty = array[pnext]->y; - float by = nexty - cury; + // check the angle p1,p2,p3 + _Tp cury = array[pcur]->y; + _Tp nexty = array[pnext]->y; + _Tp by = nexty - cury; if( CV_SIGN( by ) != nsign ) { - float ax = array[pcur]->x - array[pprev]->x; - float bx = array[pnext]->x - array[pcur]->x; - float ay = cury - array[pprev]->y; - float convexity = ay*bx - ax*by;/* if >0 then convex angle */ + _Tp ax = array[pcur]->x - array[pprev]->x; + _Tp bx = array[pnext]->x - array[pcur]->x; + _Tp ay = cury - array[pprev]->y; + _Tp convexity = ay*bx - ax*by; // if >0 then convex angle if( CV_SIGN( convexity ) == sign2 && (ax != 0 || ay != 0) ) { @@ -166,7 +97,6 @@ icvSklansky_32f( CvPoint2D32f** array, int start, int end, int* stack, int nsign stack[1] = pcur; pnext += incr; stack[2] = pnext; - } else { @@ -183,69 +113,297 @@ icvSklansky_32f( CvPoint2D32f** array, int start, int end, int* stack, int nsign stack[stacksize-1] = pnext; } } - + return --stacksize; } -typedef int (*sklansky_func)( CvPoint** points, int start, int end, - int* stack, int sign, int sign2 ); -#define cmp_pts( pt1, pt2 ) \ - ((pt1)->x < (pt2)->x || ((pt1)->x <= (pt2)->x && (pt1)->y < (pt2)->y)) -static CV_IMPLEMENT_QSORT( icvSortPointsByPointers_32s, CvPoint*, cmp_pts ) -static CV_IMPLEMENT_QSORT( icvSortPointsByPointers_32f, CvPoint2D32f*, cmp_pts ) - -static void -icvCalcAndWritePtIndices( CvPoint** pointer, int* stack, int start, int end, - CvSeq* ptseq, CvSeqWriter* writer ) +template +struct CHullCmpPoints { - int i, incr = start < end ? 1 : -1; - int idx, first_idx = ptseq->first->start_index; + bool operator()(const Point_<_Tp>* p1, const Point_<_Tp>* p2) const + { return p1->x < p2->x || (p1->x == p2->x && p1->y < p2->y); } +}; - for( i = start; i != end; i += incr ) + +void convexHull( InputArray _points, OutputArray _hull, bool clockwise, bool returnPoints ) +{ + Mat points = _points.getMat(); + int i, total = points.checkVector(2), depth = points.depth(), nout = 0; + int miny_ind = 0, maxy_ind = 0; + CV_Assert(total >= 0 && (depth == CV_32F || depth == CV_32S)); + + if( total == 0 ) { - CvPoint* ptr = (CvPoint*)pointer[stack[i]]; - CvSeqBlock* block = ptseq->first; - while( (unsigned)(idx = (int)(ptr - (CvPoint*)block->data)) >= (unsigned)block->count ) + _hull.release(); + return; + } + + returnPoints = !_hull.fixedType() ? returnPoints : _hull.type() != CV_32S; + + bool is_float = depth == CV_32F; + AutoBuffer _pointer(total); + AutoBuffer _stack(total + 2), _hullbuf(total); + Point** pointer = _pointer; + Point2f** pointerf = (Point2f**)pointer; + Point* data0 = (Point*)points.data; + int* stack = _stack; + int* hullbuf = _hullbuf; + + CV_Assert(points.isContinuous()); + + for( i = 0; i < total; i++ ) + pointer[i] = &data0[i]; + + // sort the point set by x-coordinate, find min and max y + if( !is_float ) + { + std::sort(pointer, pointer + total, CHullCmpPoints()); + for( i = 1; i < total; i++ ) { - block = block->next; - if( block == ptseq->first ) - CV_Error( CV_StsError, "Internal error" ); + int y = pointer[i]->y; + if( pointer[miny_ind]->y > y ) + miny_ind = i; + if( pointer[maxy_ind]->y < y ) + maxy_ind = i; } - idx += block->start_index - first_idx; - CV_WRITE_SEQ_ELEM( idx, *writer ); + } + else + { + std::sort(pointerf, pointerf + total, CHullCmpPoints()); + for( i = 1; i < total; i++ ) + { + float y = pointerf[i]->y; + if( pointerf[miny_ind]->y > y ) + miny_ind = i; + if( pointerf[maxy_ind]->y < y ) + maxy_ind = i; + } + } + + if( pointer[0]->x == pointer[total-1]->x && + pointer[0]->y == pointer[total-1]->y ) + { + hullbuf[nout++] = 0; + } + else + { + // upper half + int *tl_stack = stack; + int tl_count = !is_float ? + Sklansky_( pointer, 0, maxy_ind, tl_stack, -1, 1) : + Sklansky_( pointerf, 0, maxy_ind, tl_stack, -1, 1); + int *tr_stack = stack + tl_count; + int tr_count = !is_float ? + Sklansky_( pointer, total-1, maxy_ind, tr_stack, -1, -1) : + Sklansky_( pointerf, total-1, maxy_ind, tr_stack, -1, -1); + + // gather upper part of convex hull to output + if( !clockwise ) + { + std::swap( tl_stack, tr_stack ); + std::swap( tl_count, tr_count ); + } + + for( i = 0; i < tl_count-1; i++ ) + hullbuf[nout++] = pointer[tl_stack[i]] - data0; + for( i = tr_count - 1; i > 0; i-- ) + hullbuf[nout++] = pointer[tr_stack[i]] - data0; + int stop_idx = tr_count > 2 ? tr_stack[1] : tl_count > 2 ? tl_stack[tl_count - 2] : -1; + + // lower half + int *bl_stack = stack; + int bl_count = !is_float ? + Sklansky_( pointer, 0, miny_ind, bl_stack, 1, -1) : + Sklansky_( pointerf, 0, miny_ind, bl_stack, 1, -1); + int *br_stack = stack + bl_count; + int br_count = !is_float ? + Sklansky_( pointer, total-1, miny_ind, br_stack, 1, 1) : + Sklansky_( pointerf, total-1, miny_ind, br_stack, 1, 1); + + if( clockwise ) + { + std::swap( bl_stack, br_stack ); + std::swap( bl_count, br_count ); + } + + if( stop_idx >= 0 ) + { + int check_idx = bl_count > 2 ? bl_stack[1] : + bl_count + br_count > 2 ? br_stack[2-bl_count] : -1; + if( check_idx == stop_idx || (check_idx >= 0 && + pointer[check_idx]->x == pointer[stop_idx]->x && + pointer[check_idx]->y == pointer[stop_idx]->y) ) + { + // if all the points lie on the same line, then + // the bottom part of the convex hull is the mirrored top part + // (except the exteme points). + bl_count = MIN( bl_count, 2 ); + br_count = MIN( br_count, 2 ); + } + } + + for( i = 0; i < bl_count-1; i++ ) + hullbuf[nout++] = pointer[bl_stack[i]] - data0; + for( i = br_count-1; i > 0; i-- ) + hullbuf[nout++] = pointer[br_stack[i]] - data0; + } + + if( !returnPoints ) + Mat(nout, 1, CV_32S, hullbuf).copyTo(_hull); + else + { + _hull.create(nout, 1, CV_MAKETYPE(depth, 2)); + Mat hull = _hull.getMat(); + size_t step = !hull.isContinuous() ? hull.step[0] : sizeof(Point); + for( i = 0; i < nout; i++ ) + *(Point*)(hull.data + i*step) = data0[hullbuf[i]]; } } +void convexityDefects( InputArray _points, InputArray _hull, OutputArray _defects ) +{ + Mat points = _points.getMat(); + int i, j = 0, npoints = points.checkVector(2, CV_32S); + CV_Assert( npoints >= 0 ); + + if( npoints <= 3 ) + { + _defects.release(); + return; + } + + Mat hull = _hull.getMat(); + int hpoints = hull.checkVector(1, CV_32S); + CV_Assert( hpoints > 2 ); + + const Point* ptr = (const Point*)points.data; + const int* hptr = hull.ptr(); + vector defects; + + // 1. recognize co-orientation of the contour and its hull + bool rev_orientation = ((hptr[1] > hptr[0]) + (hptr[2] > hptr[1]) + (hptr[0] > hptr[2])) != 2; + + // 2. cycle through points and hull, compute defects + int hcurr = hptr[rev_orientation ? 0 : hpoints-1]; + CV_Assert( 0 <= hcurr && hcurr < npoints ); + + for( i = 0; i < hpoints; i++ ) + { + int hnext = hptr[rev_orientation ? hpoints - i - 1 : i]; + CV_Assert( 0 <= hnext && hnext < npoints ); + + Point pt0 = ptr[hcurr], pt1 = ptr[hnext]; + double dx0 = pt1.x - pt0.x; + double dy0 = pt1.y - pt0.y; + double scale = dx0 == 0 && dy0 == 0 ? 0. : 1./sqrt(dx0*dx0 + dy0*dy0); + + int defect_deepest_point = -1; + double defect_depth = 0; + bool is_defect = false; + + for(;;) + { + // go through points to achieve next hull point + j++; + j &= j >= npoints ? 0 : -1; + if( j == hnext ) + break; + + // compute distance from current point to hull edge + double dx = ptr[j].x - pt0.x; + double dy = ptr[j].y - pt0.y; + double dist = fabs(-dy0*dx + dx0*dy) * scale; + + if( dist > defect_depth ) + { + defect_depth = dist; + defect_deepest_point = j; + is_defect = true; + } + } + + if( is_defect ) + { + int idepth = cvRound(defect_depth*256); + defects.push_back(Vec4i(hcurr, hnext, defect_deepest_point, idepth)); + } + + hcurr = hnext; + } + + Mat(defects).copyTo(_defects); +} + + +template +static bool isContourConvex_( const Point_<_Tp>* p, int n ) +{ + Point_<_Tp> prev_pt = p[(n-2+n) % n]; + Point_<_Tp> cur_pt = p[n-1]; + + _Tp dx0 = cur_pt.x - prev_pt.x; + _Tp dy0 = cur_pt.y - prev_pt.y; + int orientation = 0; + + for( int i = 0; i < n-1; i++ ) + { + _Tp dxdy0, dydx0; + _Tp dx, dy; + + prev_pt = cur_pt; + cur_pt = p[i]; + + dx = cur_pt.x - prev_pt.x; + dy = cur_pt.y - prev_pt.y; + dxdy0 = dx * dy0; + dydx0 = dy * dx0; + + // find orientation + // orient = -dy0 * dx + dx0 * dy; + // orientation |= (orient > 0) ? 1 : 2; + orientation |= (dydx0 > dxdy0) ? 1 : ((dydx0 < dxdy0) ? 2 : 3); + if( orientation == 3 ) + return false; + + dx0 = dx; + dy0 = dy; + } + + return true; +} + + +bool isContourConvex( InputArray _contour ) +{ + Mat contour = _contour.getMat(); + int total = contour.checkVector(2), depth = contour.depth(); + CV_Assert(total >= 0 && (depth == CV_32F || depth == CV_32S)); + + if( total == 0 ) + return false; + + return depth == CV_32S ? + isContourConvex_((const Point*)contour.data, total ) : + isContourConvex_((const Point2f*)contour.data, total ); +} + +} + CV_IMPL CvSeq* cvConvexHull2( const CvArr* array, void* hull_storage, int orientation, int return_points ) { union { CvContour* c; CvSeq* s; } hull; - cv::AutoBuffer _pointer; - CvPoint** pointer; - CvPoint2D32f** pointerf = 0; - cv::AutoBuffer _stack; - int* stack; - hull.s = 0; CvMat* mat = 0; - CvSeqReader reader; - CvSeqWriter writer; CvContour contour_header; union { CvContour c; CvSeq s; } hull_header; CvSeqBlock block, hullblock; CvSeq* ptseq = 0; CvSeq* hullseq = 0; - int is_float; - int* t_stack; - int t_count; - int i, miny_ind = 0, maxy_ind = 0, total; - int hulltype; - int stop_idx; - sklansky_func sklansky; if( CV_IS_SEQ( array )) { @@ -264,17 +422,16 @@ cvConvexHull2( const CvArr* array, void* hull_storage, { if( return_points ) { - hullseq = cvCreateSeq( - CV_SEQ_KIND_CURVE|CV_SEQ_ELTYPE(ptseq)| - CV_SEQ_FLAG_CLOSED|CV_SEQ_FLAG_CONVEX, - sizeof(CvContour), sizeof(CvPoint),(CvMemStorage*)hull_storage ); + hullseq = cvCreateSeq(CV_SEQ_KIND_CURVE|CV_SEQ_ELTYPE(ptseq)| + CV_SEQ_FLAG_CLOSED|CV_SEQ_FLAG_CONVEX, + sizeof(CvContour), sizeof(CvPoint),(CvMemStorage*)hull_storage ); } else { hullseq = cvCreateSeq( - CV_SEQ_KIND_CURVE|CV_SEQ_ELTYPE_PPOINT| - CV_SEQ_FLAG_CLOSED|CV_SEQ_FLAG_CONVEX, - sizeof(CvContour), sizeof(CvPoint*), (CvMemStorage*)hull_storage ); + CV_SEQ_KIND_CURVE|CV_SEQ_ELTYPE_PPOINT| + CV_SEQ_FLAG_CLOSED|CV_SEQ_FLAG_CONVEX, + sizeof(CvContour), sizeof(CvPoint*), (CvMemStorage*)hull_storage ); } } else @@ -286,192 +443,50 @@ cvConvexHull2( const CvArr* array, void* hull_storage, if( (mat->cols != 1 && mat->rows != 1) || !CV_IS_MAT_CONT(mat->type)) CV_Error( CV_StsBadArg, - "The hull matrix should be continuous and have a single row or a single column" ); + "The hull matrix should be continuous and have a single row or a single column" ); if( mat->cols + mat->rows - 1 < ptseq->total ) CV_Error( CV_StsBadSize, "The hull matrix size might be not enough to fit the hull" ); if( CV_MAT_TYPE(mat->type) != CV_SEQ_ELTYPE(ptseq) && - CV_MAT_TYPE(mat->type) != CV_32SC1 ) + CV_MAT_TYPE(mat->type) != CV_32SC1 ) CV_Error( CV_StsUnsupportedFormat, - "The hull matrix must have the same type as input or 32sC1 (integers)" ); + "The hull matrix must have the same type as input or 32sC1 (integers)" ); hullseq = cvMakeSeqHeaderForArray( - CV_SEQ_KIND_CURVE|CV_MAT_TYPE(mat->type)|CV_SEQ_FLAG_CLOSED, - sizeof(contour_header), CV_ELEM_SIZE(mat->type), mat->data.ptr, - mat->cols + mat->rows - 1, &hull_header.s, &hullblock ); - + CV_SEQ_KIND_CURVE|CV_MAT_TYPE(mat->type)|CV_SEQ_FLAG_CLOSED, + sizeof(contour_header), CV_ELEM_SIZE(mat->type), mat->data.ptr, + mat->cols + mat->rows - 1, &hull_header.s, &hullblock ); cvClearSeq( hullseq ); } - total = ptseq->total; + int hulltype = CV_SEQ_ELTYPE(hullseq); + int total = ptseq->total; if( total == 0 ) { if( mat ) CV_Error( CV_StsBadSize, - "Point sequence can not be empty if the output is matrix" ); + "Point sequence can not be empty if the output is matrix" ); return hull.s; } - cvStartAppendToSeq( hullseq, &writer ); + cv::AutoBuffer _ptbuf; + cv::Mat h0; + cv::convexHull(cv::cvarrToMat(ptseq, false, false, 0, &_ptbuf), h0, + orientation == CV_CLOCKWISE, CV_MAT_CN(hulltype) == 2); - is_float = CV_SEQ_ELTYPE(ptseq) == CV_32FC2; - hulltype = CV_SEQ_ELTYPE(hullseq); - sklansky = !is_float ? (sklansky_func)icvSklansky_32s : - (sklansky_func)icvSklansky_32f; - - _pointer.allocate( ptseq->total ); - _stack.allocate( ptseq->total + 2); - pointer = _pointer; - pointerf = (CvPoint2D32f**)pointer; - stack = _stack; - - cvStartReadSeq( ptseq, &reader ); - - for( i = 0; i < total; i++ ) + if( hulltype == CV_SEQ_ELTYPE_PPOINT ) { - pointer[i] = (CvPoint*)reader.ptr; - CV_NEXT_SEQ_ELEM( ptseq->elem_size, reader ); - } - - // sort the point set by x-coordinate, find min and max y - if( !is_float ) - { - icvSortPointsByPointers_32s( pointer, total, 0 ); - for( i = 1; i < total; i++ ) + const int* idx = h0.ptr(); + int ctotal = (int)h0.total(); + for( int i = 0; i < ctotal; i++ ) { - int y = pointer[i]->y; - if( pointer[miny_ind]->y > y ) - miny_ind = i; - if( pointer[maxy_ind]->y < y ) - maxy_ind = i; + void* ptr = cvGetSeqElem(ptseq, idx[i]); + cvSeqPush( hullseq, &ptr ); } } else - { - icvSortPointsByPointers_32f( pointerf, total, 0 ); - for( i = 1; i < total; i++ ) - { - float y = pointerf[i]->y; - if( pointerf[miny_ind]->y > y ) - miny_ind = i; - if( pointerf[maxy_ind]->y < y ) - maxy_ind = i; - } - } - - if( pointer[0]->x == pointer[total-1]->x && - pointer[0]->y == pointer[total-1]->y ) - { - if( hulltype == CV_SEQ_ELTYPE_PPOINT ) - { - CV_WRITE_SEQ_ELEM( pointer[0], writer ); - } - else if( hulltype == CV_SEQ_ELTYPE_INDEX ) - { - int index = 0; - CV_WRITE_SEQ_ELEM( index, writer ); - } - else - { - CvPoint pt = pointer[0][0]; - CV_WRITE_SEQ_ELEM( pt, writer ); - } - goto finish_hull; - } - - /*upper half */ - { - int *tl_stack = stack; - int tl_count = sklansky( pointer, 0, maxy_ind, tl_stack, -1, 1 ); - int *tr_stack = tl_stack + tl_count; - int tr_count = sklansky( pointer, ptseq->total - 1, maxy_ind, tr_stack, -1, -1 ); - - /* gather upper part of convex hull to output */ - if( orientation == CV_COUNTER_CLOCKWISE ) - { - CV_SWAP( tl_stack, tr_stack, t_stack ); - CV_SWAP( tl_count, tr_count, t_count ); - } - - if( hulltype == CV_SEQ_ELTYPE_PPOINT ) - { - for( i = 0; i < tl_count - 1; i++ ) - CV_WRITE_SEQ_ELEM( pointer[tl_stack[i]], writer ); - - for( i = tr_count - 1; i > 0; i-- ) - CV_WRITE_SEQ_ELEM( pointer[tr_stack[i]], writer ); - } - else if( hulltype == CV_SEQ_ELTYPE_INDEX ) - { - icvCalcAndWritePtIndices( pointer, tl_stack, 0, tl_count-1, ptseq, &writer ); - icvCalcAndWritePtIndices( pointer, tr_stack, tr_count-1, 0, ptseq, &writer ); - } - else - { - for( i = 0; i < tl_count - 1; i++ ) - CV_WRITE_SEQ_ELEM( pointer[tl_stack[i]][0], writer ); - - for( i = tr_count - 1; i > 0; i-- ) - CV_WRITE_SEQ_ELEM( pointer[tr_stack[i]][0], writer ); - } - stop_idx = tr_count > 2 ? tr_stack[1] : tl_count > 2 ? tl_stack[tl_count - 2] : -1; - } - - /* lower half */ - { - int *bl_stack = stack; - int bl_count = sklansky( pointer, 0, miny_ind, bl_stack, 1, -1 ); - int *br_stack = stack + bl_count; - int br_count = sklansky( pointer, ptseq->total - 1, miny_ind, br_stack, 1, 1 ); - - if( orientation != CV_COUNTER_CLOCKWISE ) - { - CV_SWAP( bl_stack, br_stack, t_stack ); - CV_SWAP( bl_count, br_count, t_count ); - } - - if( stop_idx >= 0 ) - { - int check_idx = bl_count > 2 ? bl_stack[1] : - bl_count + br_count > 2 ? br_stack[2-bl_count] : -1; - if( check_idx == stop_idx || (check_idx >= 0 && - pointer[check_idx]->x == pointer[stop_idx]->x && - pointer[check_idx]->y == pointer[stop_idx]->y) ) - { - /* if all the points lie on the same line, then - the bottom part of the convex hull is the mirrored top part - (except the exteme points).*/ - bl_count = MIN( bl_count, 2 ); - br_count = MIN( br_count, 2 ); - } - } - - if( hulltype == CV_SEQ_ELTYPE_PPOINT ) - { - for( i = 0; i < bl_count - 1; i++ ) - CV_WRITE_SEQ_ELEM( pointer[bl_stack[i]], writer ); - - for( i = br_count - 1; i > 0; i-- ) - CV_WRITE_SEQ_ELEM( pointer[br_stack[i]], writer ); - } - else if( hulltype == CV_SEQ_ELTYPE_INDEX ) - { - icvCalcAndWritePtIndices( pointer, bl_stack, 0, bl_count-1, ptseq, &writer ); - icvCalcAndWritePtIndices( pointer, br_stack, br_count-1, 0, ptseq, &writer ); - } - else - { - for( i = 0; i < bl_count - 1; i++ ) - CV_WRITE_SEQ_ELEM( pointer[bl_stack[i]][0], writer ); - - for( i = br_count - 1; i > 0; i-- ) - CV_WRITE_SEQ_ELEM( pointer[br_stack[i]][0], writer ); - } - } - -finish_hull: - cvEndWriteSeq( &writer ); + cvSeqPushMulti(hullseq, h0.data, (int)h0.total()); if( mat ) { @@ -484,13 +499,10 @@ finish_hull: { hull.s = hullseq; hull.c->rect = cvBoundingRect( ptseq, - ptseq->header_size < (int)sizeof(CvContour) || - &ptseq->flags == &contour_header.flags ); - - /*if( ptseq != (CvSeq*)&contour_header ) - hullseq->v_prev = ptseq;*/ + ptseq->header_size < (int)sizeof(CvContour) || + &ptseq->flags == &contour_header.flags ); } - + return hull.s; } @@ -498,8 +510,8 @@ finish_hull: /* contour must be a simple polygon */ /* it must have more than 3 points */ CV_IMPL CvSeq* cvConvexityDefects( const CvArr* array, - const CvArr* hullarray, - CvMemStorage* storage ) + const CvArr* hullarray, + CvMemStorage* storage ) { CvSeq* defects = 0; @@ -523,7 +535,7 @@ CV_IMPL CvSeq* cvConvexityDefects( const CvArr* array, { if( !CV_IS_SEQ_POINT_SET( ptseq )) CV_Error( CV_StsUnsupportedFormat, - "Input sequence is not a sequence of points" ); + "Input sequence is not a sequence of points" ); if( !storage ) storage = ptseq->storage; } @@ -540,8 +552,8 @@ CV_IMPL CvSeq* cvConvexityDefects( const CvArr* array, int hulltype = CV_SEQ_ELTYPE( hull ); if( hulltype != CV_SEQ_ELTYPE_PPOINT && hulltype != CV_SEQ_ELTYPE_INDEX ) CV_Error( CV_StsUnsupportedFormat, - "Convex hull must represented as a sequence " - "of indices or sequence of pointers" ); + "Convex hull must represented as a sequence " + "of indices or sequence of pointers" ); if( !storage ) storage = hull->storage; } @@ -553,17 +565,17 @@ CV_IMPL CvSeq* cvConvexityDefects( const CvArr* array, CV_Error(CV_StsBadArg, "Convex hull is neither sequence nor matrix"); if( (mat->cols != 1 && mat->rows != 1) || - !CV_IS_MAT_CONT(mat->type) || CV_MAT_TYPE(mat->type) != CV_32SC1 ) + !CV_IS_MAT_CONT(mat->type) || CV_MAT_TYPE(mat->type) != CV_32SC1 ) CV_Error( CV_StsBadArg, - "The matrix should be 1-dimensional and continuous array of int's" ); + "The matrix should be 1-dimensional and continuous array of int's" ); if( mat->cols + mat->rows - 1 > ptseq->total ) CV_Error( CV_StsBadSize, "Convex hull is larger than the point sequence" ); hull = cvMakeSeqHeaderForArray( - CV_SEQ_KIND_CURVE|CV_MAT_TYPE(mat->type)|CV_SEQ_FLAG_CLOSED, - sizeof(CvContour), CV_ELEM_SIZE(mat->type), mat->data.ptr, - mat->cols + mat->rows - 1, &hull_header.s, &hullblock ); + CV_SEQ_KIND_CURVE|CV_MAT_TYPE(mat->type)|CV_SEQ_FLAG_CLOSED, + sizeof(CvContour), CV_ELEM_SIZE(mat->type), mat->data.ptr, + mat->cols + mat->rows - 1, &hull_header.s, &hullblock ); } is_index = CV_SEQ_ELTYPE(hull) == CV_SEQ_ELTYPE_INDEX; @@ -701,11 +713,6 @@ CV_IMPL CvSeq* cvConvexityDefects( const CvArr* array, CV_IMPL int cvCheckContourConvexity( const CvArr* array ) { - int flag = -1; - - int i; - int orientation = 0; - CvSeqReader reader; CvContour contour_header; CvSeqBlock block; CvSeq* contour = (CvSeq*)array; @@ -714,102 +721,19 @@ cvCheckContourConvexity( const CvArr* array ) { if( !CV_IS_SEQ_POINT_SET(contour)) CV_Error( CV_StsUnsupportedFormat, - "Input sequence must be polygon (closed 2d curve)" ); + "Input sequence must be polygon (closed 2d curve)" ); } else { - contour = cvPointSeqFromMat(CV_SEQ_KIND_CURVE|CV_SEQ_FLAG_CLOSED, array, &contour_header, &block ); + contour = cvPointSeqFromMat(CV_SEQ_KIND_CURVE| + CV_SEQ_FLAG_CLOSED, array, &contour_header, &block ); } if( contour->total == 0 ) return -1; - cvStartReadSeq( contour, &reader, 0 ); - flag = 1; - - if( CV_SEQ_ELTYPE( contour ) == CV_32SC2 ) - { - CvPoint *prev_pt = (CvPoint*)reader.prev_elem; - CvPoint *cur_pt = (CvPoint*)reader.ptr; - - int dx0 = cur_pt->x - prev_pt->x; - int dy0 = cur_pt->y - prev_pt->y; - - for( i = 0; i < contour->total; i++ ) - { - int dxdy0, dydx0; - int dx, dy; - - /*int orient; */ - CV_NEXT_SEQ_ELEM( sizeof(CvPoint), reader ); - prev_pt = cur_pt; - cur_pt = (CvPoint *) reader.ptr; - - dx = cur_pt->x - prev_pt->x; - dy = cur_pt->y - prev_pt->y; - dxdy0 = dx * dy0; - dydx0 = dy * dx0; - - /* find orientation */ - /*orient = -dy0 * dx + dx0 * dy; - orientation |= (orient > 0) ? 1 : 2; - */ - orientation |= (dydx0 > dxdy0) ? 1 : ((dydx0 < dxdy0) ? 2 : 3); - - if( orientation == 3 ) - { - flag = 0; - break; - } - - dx0 = dx; - dy0 = dy; - } - } - else - { - CV_Assert( CV_SEQ_ELTYPE(contour) == CV_32FC2 ); - - CvPoint2D32f *prev_pt = (CvPoint2D32f*)reader.prev_elem; - CvPoint2D32f *cur_pt = (CvPoint2D32f*)reader.ptr; - - float dx0 = cur_pt->x - prev_pt->x; - float dy0 = cur_pt->y - prev_pt->y; - - for( i = 0; i < contour->total; i++ ) - { - float dxdy0, dydx0; - float dx, dy; - - /*int orient; */ - CV_NEXT_SEQ_ELEM( sizeof(CvPoint2D32f), reader ); - prev_pt = cur_pt; - cur_pt = (CvPoint2D32f*) reader.ptr; - - dx = cur_pt->x - prev_pt->x; - dy = cur_pt->y - prev_pt->y; - dxdy0 = dx * dy0; - dydx0 = dy * dx0; - - /* find orientation */ - /*orient = -dy0 * dx + dx0 * dy; - orientation |= (orient > 0) ? 1 : 2; - */ - orientation |= (dydx0 > dxdy0) ? 1 : ((dydx0 < dxdy0) ? 2 : 3); - - if( orientation == 3 ) - { - flag = 0; - break; - } - - dx0 = dx; - dy0 = dy; - } - } - - return flag; + cv::AutoBuffer _buf; + return cv::isContourConvex(cv::cvarrToMat(contour, false, false, 0, &_buf)) ? 1 : 0; } - /* End of file. */ diff --git a/modules/imgproc/src/geometry.cpp b/modules/imgproc/src/geometry.cpp index 337362f12..e0eb22971 100644 --- a/modules/imgproc/src/geometry.cpp +++ b/modules/imgproc/src/geometry.cpp @@ -92,97 +92,38 @@ cvBoxPoints( CvBox2D box, CvPoint2D32f pt[4] ) } -int -icvIntersectLines( double x1, double dx1, double y1, double dy1, - double x2, double dx2, double y2, double dy2, double *t2 ) -{ - double d = dx1 * dy2 - dx2 * dy1; - int result = -1; - - if( d != 0 ) - { - *t2 = ((x2 - x1) * dy1 - (y2 - y1) * dx1) / d; - result = 0; - } - return result; -} - - -void -icvIntersectLines3( double *a0, double *b0, double *c0, - double *a1, double *b1, double *c1, CvPoint2D32f * point ) -{ - double det = a0[0] * b1[0] - a1[0] * b0[0]; - - if( det != 0 ) - { - det = 1. / det; - point->x = (float) ((b0[0] * c1[0] - b1[0] * c0[0]) * det); - point->y = (float) ((a1[0] * c0[0] - a0[0] * c1[0]) * det); - } - else - { - point->x = point->y = FLT_MAX; - } -} - - -CV_IMPL double -cvPointPolygonTest( const CvArr* _contour, CvPoint2D32f pt, int measure_dist ) +double cv::pointPolygonTest( InputArray _contour, Point2f pt, bool measureDist ) { double result = 0; + Mat contour = _contour.getMat(); + int i, total = contour.checkVector(2), counter = 0; + int depth = contour.depth(); + CV_Assert( total >= 0 && (depth == CV_32S || depth == CV_32F)); - CvSeqBlock block; - CvContour header; - CvSeq* contour = (CvSeq*)_contour; - CvSeqReader reader; - int i, total, counter = 0; - int is_float; + bool is_float = depth == CV_32F; double min_dist_num = FLT_MAX, min_dist_denom = 1; - CvPoint ip = {0,0}; + Point ip(cvRound(pt.x), cvRound(pt.y)); - if( !CV_IS_SEQ(contour) ) - { - contour = cvPointSeqFromMat( CV_SEQ_KIND_CURVE + CV_SEQ_FLAG_CLOSED, - _contour, &header, &block ); - } - else if( CV_IS_SEQ_POINT_SET(contour) ) - { - if( contour->header_size == sizeof(CvContour) && !measure_dist ) - { - CvRect r = ((CvContour*)contour)->rect; - if( pt.x < r.x || pt.y < r.y || - pt.x >= r.x + r.width || pt.y >= r.y + r.height ) - return -1; - } - } - else if( CV_IS_SEQ_CHAIN(contour) ) - { - CV_Error( CV_StsBadArg, - "Chains are not supported. Convert them to polygonal representation using cvApproxChains()" ); - } - else - CV_Error( CV_StsBadArg, "Input contour is neither a valid sequence nor a matrix" ); + if( total == 0 ) + return measureDist ? -DBL_MAX : -1; - total = contour->total; - is_float = CV_SEQ_ELTYPE(contour) == CV_32FC2; - cvStartReadSeq( contour, &reader, -1 ); + const Point* cnt = (const Point*)contour.data; + const Point2f* cntf = (const Point2f*)cnt; - if( !is_float && !measure_dist && (ip.x = cvRound(pt.x)) == pt.x && (ip.y = cvRound(pt.y)) == pt.y ) + if( !is_float && !measureDist && ip.x == pt.x && ip.y == pt.y ) { - // the fastest "pure integer" branch - CvPoint v0, v; - CV_READ_SEQ_ELEM( v, reader ); + // the fastest "purely integer" branch + Point v0, v = cnt[total-1]; for( i = 0; i < total; i++ ) { int dist; v0 = v; - CV_READ_SEQ_ELEM( v, reader ); + v = cnt[i]; if( (v0.y <= ip.y && v.y <= ip.y) || - (v0.y > ip.y && v.y > ip.y) || - (v0.x < ip.x && v.x < ip.x) ) + (v0.y > ip.y && v.y > ip.y) || + (v0.x < ip.x && v.x < ip.x) ) { if( ip.y == v.y && (ip.x == v.x || (ip.y == v0.y && ((v0.x <= ip.x && ip.x <= v.x) || (v.x <= ip.x && ip.x <= v0.x)))) ) @@ -202,38 +143,32 @@ cvPointPolygonTest( const CvArr* _contour, CvPoint2D32f pt, int measure_dist ) } else { - CvPoint2D32f v0, v; - CvPoint iv; + Point2f v0, v; + Point iv; if( is_float ) { - CV_READ_SEQ_ELEM( v, reader ); + v = cntf[total-1]; } else { - CV_READ_SEQ_ELEM( iv, reader ); - v = cvPointTo32f( iv ); + v = cnt[total-1]; } - if( !measure_dist ) + if( !measureDist ) { for( i = 0; i < total; i++ ) { double dist; v0 = v; if( is_float ) - { - CV_READ_SEQ_ELEM( v, reader ); - } + v = cntf[i]; else - { - CV_READ_SEQ_ELEM( iv, reader ); - v = cvPointTo32f( iv ); - } + v = cnt[i]; if( (v0.y <= pt.y && v.y <= pt.y) || - (v0.y > pt.y && v.y > pt.y) || - (v0.x < pt.x && v.x < pt.x) ) + (v0.y > pt.y && v.y > pt.y) || + (v0.x < pt.x && v.x < pt.x) ) { if( pt.y == v.y && (pt.x == v.x || (pt.y == v0.y && ((v0.x <= pt.x && pt.x <= v.x) || (v.x <= pt.x && pt.x <= v0.x)))) ) @@ -259,14 +194,9 @@ cvPointPolygonTest( const CvArr* _contour, CvPoint2D32f pt, int measure_dist ) v0 = v; if( is_float ) - { - CV_READ_SEQ_ELEM( v, reader ); - } + v = cntf[i]; else - { - CV_READ_SEQ_ELEM( iv, reader ); - v = cvPointTo32f( iv ); - } + v = cnt[i]; dx = v.x - v0.x; dy = v.y - v0.y; dx1 = pt.x - v0.x; dy1 = pt.y - v0.y; @@ -292,8 +222,8 @@ cvPointPolygonTest( const CvArr* _contour, CvPoint2D32f pt, int measure_dist ) } if( (v0.y <= pt.y && v.y <= pt.y) || - (v0.y > pt.y && v.y > pt.y) || - (v0.x < pt.x && v.x < pt.x) ) + (v0.y > pt.y && v.y > pt.y) || + (v0.x < pt.x && v.x < pt.x) ) continue; dist_num = dy1*dx - dx1*dy; @@ -301,17 +231,25 @@ cvPointPolygonTest( const CvArr* _contour, CvPoint2D32f pt, int measure_dist ) dist_num = -dist_num; counter += dist_num > 0; } - + result = sqrt(min_dist_num/min_dist_denom); if( counter % 2 == 0 ) result = -result; } } - + return result; } +CV_IMPL double +cvPointPolygonTest( const CvArr* _contour, CvPoint2D32f pt, int measure_dist ) +{ + cv::AutoBuffer abuf; + cv::Mat contour = cv::cvarrToMat(_contour, false, false, 0, &abuf); + return cv::pointPolygonTest(contour, pt, measure_dist != 0); +} + /* This code is described in "Computational Geometry in C" (Second Edition), Chapter 7. It is not written to be comprehensible without the diff --git a/modules/imgproc/src/linefit.cpp b/modules/imgproc/src/linefit.cpp index c5135659b..e31663620 100644 --- a/modules/imgproc/src/linefit.cpp +++ b/modules/imgproc/src/linefit.cpp @@ -40,19 +40,19 @@ //M*/ #include "precomp.hpp" +namespace cv +{ + static const double eps = 1e-6; -static CvStatus -icvFitLine2D_wods( CvPoint2D32f * points, int _count, float *weights, float *line ) +static void fitLine2D_wods( const Point2f* points, int count, float *weights, float *line ) { double x = 0, y = 0, x2 = 0, y2 = 0, xy = 0, w = 0; double dx2, dy2, dxy; int i; - int count = _count; float t; - /* Calculating the average of x and y... */ - + // Calculating the average of x and y... if( weights == 0 ) { for( i = 0; i < count; i += 1 ) @@ -94,12 +94,9 @@ icvFitLine2D_wods( CvPoint2D32f * points, int _count, float *weights, float *lin line[2] = (float) x; line[3] = (float) y; - - return CV_NO_ERR; } -static CvStatus -icvFitLine3D_wods( CvPoint3D32f * points, int count, float *weights, float *line ) +static void fitLine3D_wods( const Point3f * points, int count, float *weights, float *line ) { int i; float w0 = 0; @@ -184,25 +181,13 @@ icvFitLine3D_wods( CvPoint3D32f * points, int count, float *weights, float *line det[7] = det[5]; det[8] = dy2 + dx2; - /* Searching for a eigenvector of det corresponding to the minimal eigenvalue */ -#if 1 - { - CvMat _det = cvMat( 3, 3, CV_32F, det ); - CvMat _evc = cvMat( 3, 3, CV_32F, evc ); - CvMat _evl = cvMat( 3, 1, CV_32F, evl ); - cvEigenVV( &_det, &_evc, &_evl, 0 ); + // Searching for a eigenvector of det corresponding to the minimal eigenvalue + Mat _det( 3, 3, CV_32F, det ); + Mat _evc( 3, 3, CV_32F, evc ); + Mat _evl( 3, 1, CV_32F, evl ); + eigen( _det, _evl, _evc ); i = evl[0] < evl[1] ? (evl[0] < evl[2] ? 0 : 2) : (evl[1] < evl[2] ? 1 : 2); - } -#else - { - CvMat _det = cvMat( 3, 3, CV_32F, det ); - CvMat _evc = cvMat( 3, 3, CV_32F, evc ); - CvMat _evl = cvMat( 1, 3, CV_32F, evl ); - cvSVD( &_det, &_evl, &_evc, 0, CV_SVD_MODIFY_A+CV_SVD_U_T ); - } - i = 2; -#endif v = &evc[i * 3]; n = (float) sqrt( (double)v[0] * v[0] + (double)v[1] * v[1] + (double)v[2] * v[2] ); n = (float)MAX(n, eps); @@ -212,12 +197,9 @@ icvFitLine3D_wods( CvPoint3D32f * points, int count, float *weights, float *line line[3] = x0; line[4] = y0; line[5] = z0; - - return CV_NO_ERR; } -static double -icvCalcDist2D( CvPoint2D32f * points, int count, float *_line, float *dist ) +static double calcDist2D( const Point2f* points, int count, float *_line, float *dist ) { int j; float px = _line[2], py = _line[3]; @@ -238,8 +220,7 @@ icvCalcDist2D( CvPoint2D32f * points, int count, float *_line, float *dist ) return sum_dist; } -static double -icvCalcDist3D( CvPoint3D32f * points, int count, float *_line, float *dist ) +static double calcDist3D( const Point3f* points, int count, float *_line, float *dist ) { int j; float px = _line[3], py = _line[4], pz = _line[5]; @@ -266,8 +247,7 @@ icvCalcDist3D( CvPoint3D32f * points, int count, float *_line, float *dist ) return sum_dist; } -static void -icvWeightL1( float *d, int count, float *w ) +static void weightL1( float *d, int count, float *w ) { int i; @@ -278,8 +258,7 @@ icvWeightL1( float *d, int count, float *w ) } } -static void -icvWeightL12( float *d, int count, float *w ) +static void weightL12( float *d, int count, float *w ) { int i; @@ -290,8 +269,7 @@ icvWeightL12( float *d, int count, float *w ) } -static void -icvWeightHuber( float *d, int count, float *w, float _c ) +static void weightHuber( float *d, int count, float *w, float _c ) { int i; const float c = _c <= 0 ? 1.345f : _c; @@ -306,8 +284,7 @@ icvWeightHuber( float *d, int count, float *w, float _c ) } -static void -icvWeightFair( float *d, int count, float *w, float _c ) +static void weightFair( float *d, int count, float *w, float _c ) { int i; const float c = _c == 0 ? 1 / 1.3998f : 1 / _c; @@ -318,76 +295,72 @@ icvWeightFair( float *d, int count, float *w, float _c ) } } -static void -icvWeightWelsch( float *d, int count, float *w, float _c ) +static void weightWelsch( float *d, int count, float *w, float _c ) { int i; const float c = _c == 0 ? 1 / 2.9846f : 1 / _c; for( i = 0; i < count; i++ ) { - w[i] = (float) exp( -d[i] * d[i] * c * c ); + w[i] = (float) std::exp( -d[i] * d[i] * c * c ); } } /* Takes an array of 2D points, type of distance (including user-defined -distance specified by callbacks, fills the array of four floats with line -parameters A, B, C, D, where (A, B) is the normalized direction vector, -(C, D) is the point that belongs to the line. */ + distance specified by callbacks, fills the array of four floats with line + parameters A, B, C, D, where (A, B) is the normalized direction vector, + (C, D) is the point that belongs to the line. */ -static CvStatus icvFitLine2D( CvPoint2D32f * points, int count, int dist, - float _param, float reps, float aeps, float *line ) +static void fitLine2D( const Point2f * points, int count, int dist, + float _param, float reps, float aeps, float *line ) { double EPS = count*FLT_EPSILON; void (*calc_weights) (float *, int, float *) = 0; void (*calc_weights_param) (float *, int, float *, float) = 0; - float *w; /* weights */ - float *r; /* square distances */ int i, j, k; float _line[6], _lineprev[6]; float rdelta = reps != 0 ? reps : 1.0f; float adelta = aeps != 0 ? aeps : 0.01f; double min_err = DBL_MAX, err = 0; - CvRNG rng = cvRNG(-1); + RNG rng((uint64)-1); memset( line, 0, 4*sizeof(line[0]) ); switch (dist) { case CV_DIST_L2: - return icvFitLine2D_wods( points, count, 0, line ); + return fitLine2D_wods( points, count, 0, line ); case CV_DIST_L1: - calc_weights = icvWeightL1; + calc_weights = weightL1; break; case CV_DIST_L12: - calc_weights = icvWeightL12; + calc_weights = weightL12; break; case CV_DIST_FAIR: - calc_weights_param = icvWeightFair; + calc_weights_param = weightFair; break; case CV_DIST_WELSCH: - calc_weights_param = icvWeightWelsch; + calc_weights_param = weightWelsch; break; case CV_DIST_HUBER: - calc_weights_param = icvWeightHuber; + calc_weights_param = weightHuber; break; - /*case CV_DIST_USER: - calc_weights = (void ( * )(float *, int, float *)) _PFP.fp; - break;*/ - + /*case DIST_USER: + calc_weights = (void ( * )(float *, int, float *)) _PFP.fp; + break;*/ default: - return CV_BADFACTOR_ERR; + CV_Error(CV_StsBadArg, "Unknown distance type"); } - w = (float *) cvAlloc( count * sizeof( float )); - r = (float *) cvAlloc( count * sizeof( float )); + AutoBuffer wr(count*2); + float *w = wr, *r = w + count; for( k = 0; k < 20; k++ ) { @@ -397,7 +370,7 @@ static CvStatus icvFitLine2D( CvPoint2D32f * points, int count, int dist, for( i = 0; i < MIN(count,10); ) { - j = cvRandInt(&rng) % count; + j = rng.uniform(0, count); if( w[j] < FLT_EPSILON ) { w[j] = 1.f; @@ -405,7 +378,7 @@ static CvStatus icvFitLine2D( CvPoint2D32f * points, int count, int dist, } } - icvFitLine2D_wods( points, count, w, _line ); + fitLine2D_wods( points, count, w, _line ); for( i = 0; i < 30; i++ ) { double sum_w = 0; @@ -432,7 +405,7 @@ static CvStatus icvFitLine2D( CvPoint2D32f * points, int count, int dist, } } /* calculate distances */ - err = icvCalcDist2D( points, count, _line, r ); + err = calcDist2D( points, count, _line, r ); if( err < EPS ) break; @@ -461,7 +434,7 @@ static CvStatus icvFitLine2D( CvPoint2D32f * points, int count, int dist, memcpy( _lineprev, _line, 4 * sizeof( float )); /* Run again... */ - icvFitLine2D_wods( points, count, w, _line ); + fitLine2D_wods( points, count, w, _line ); } if( err < min_err ) @@ -472,70 +445,57 @@ static CvStatus icvFitLine2D( CvPoint2D32f * points, int count, int dist, break; } } - - cvFree( &w ); - cvFree( &r ); - return CV_OK; } /* Takes an array of 3D points, type of distance (including user-defined -distance specified by callbacks, fills the array of four floats with line -parameters A, B, C, D, E, F, where (A, B, C) is the normalized direction vector, -(D, E, F) is the point that belongs to the line. */ - -static CvStatus -icvFitLine3D( CvPoint3D32f * points, int count, int dist, - float _param, float reps, float aeps, float *line ) + distance specified by callbacks, fills the array of four floats with line + parameters A, B, C, D, E, F, where (A, B, C) is the normalized direction vector, + (D, E, F) is the point that belongs to the line. */ +static void fitLine3D( Point3f * points, int count, int dist, + float _param, float reps, float aeps, float *line ) { double EPS = count*FLT_EPSILON; void (*calc_weights) (float *, int, float *) = 0; void (*calc_weights_param) (float *, int, float *, float) = 0; - float *w; /* weights */ - float *r; /* square distances */ int i, j, k; float _line[6]={0,0,0,0,0,0}, _lineprev[6]={0,0,0,0,0,0}; float rdelta = reps != 0 ? reps : 1.0f; float adelta = aeps != 0 ? aeps : 0.01f; double min_err = DBL_MAX, err = 0; - CvRNG rng = cvRNG(-1); + RNG rng((uint64)-1); switch (dist) { case CV_DIST_L2: - return icvFitLine3D_wods( points, count, 0, line ); + return fitLine3D_wods( points, count, 0, line ); case CV_DIST_L1: - calc_weights = icvWeightL1; + calc_weights = weightL1; break; case CV_DIST_L12: - calc_weights = icvWeightL12; + calc_weights = weightL12; break; case CV_DIST_FAIR: - calc_weights_param = icvWeightFair; + calc_weights_param = weightFair; break; case CV_DIST_WELSCH: - calc_weights_param = icvWeightWelsch; + calc_weights_param = weightWelsch; break; case CV_DIST_HUBER: - calc_weights_param = icvWeightHuber; + calc_weights_param = weightHuber; break; - /*case CV_DIST_USER: - _PFP.p = param; - calc_weights = (void ( * )(float *, int, float *)) _PFP.fp; - break;*/ - default: - return CV_BADFACTOR_ERR; + CV_Error(CV_StsBadArg, "Unknown distance"); } - w = (float *) cvAlloc( count * sizeof( float )); - r = (float *) cvAlloc( count * sizeof( float )); + AutoBuffer buf(count*2); + float *w = buf, *r = w + count; for( k = 0; k < 20; k++ ) { @@ -545,7 +505,7 @@ icvFitLine3D( CvPoint3D32f * points, int count, int dist, for( i = 0; i < MIN(count,10); ) { - j = cvRandInt(&rng) % count; + j = rng.uniform(0, count); if( w[j] < FLT_EPSILON ) { w[j] = 1.f; @@ -553,7 +513,7 @@ icvFitLine3D( CvPoint3D32f * points, int count, int dist, } } - icvFitLine3D_wods( points, count, w, _line ); + fitLine3D_wods( points, count, w, _line ); for( i = 0; i < 30; i++ ) { double sum_w = 0; @@ -587,8 +547,9 @@ icvFitLine3D( CvPoint3D32f * points, int count, int dist, } } /* calculate distances */ - if( icvCalcDist3D( points, count, _line, r ) < FLT_EPSILON*count ) - break; + err = calcDist3D( points, count, _line, r ); + //if( err < FLT_EPSILON*count ) + // break; /* calculate weights */ if( calc_weights ) @@ -610,14 +571,14 @@ icvFitLine3D( CvPoint3D32f * points, int count, int dist, for( j = 0; j < count; j++ ) w[j] = 1.f; } - + /* save the line parameters */ memcpy( _lineprev, _line, 6 * sizeof( float )); - + /* Run again... */ - icvFitLine3D_wods( points, count, w, _line ); + fitLine3D_wods( points, count, w, _line ); } - + if( err < min_err ) { min_err = err; @@ -626,11 +587,36 @@ icvFitLine3D( CvPoint3D32f * points, int count, int dist, break; } } +} - // Return... - cvFree( &w ); - cvFree( &r ); - return CV_OK; +} + +void cv::fitLine( InputArray _points, OutputArray _line, int distType, + double param, double reps, double aeps ) +{ + Mat points = _points.getMat(); + + float linebuf[6]={0.f}; + int npoints2 = points.checkVector(2, -1, false); + int npoints3 = points.checkVector(3, -1, false); + + CV_Assert( npoints2 >= 0 || npoints3 >= 0 ); + + if( points.depth() != CV_32F || !points.isContinuous() ) + { + Mat temp; + points.convertTo(temp, CV_32F); + points = temp; + } + + if( npoints2 >= 0 ) + fitLine2D( points.ptr(), npoints2, distType, + (float)param, (float)reps, (float)aeps, linebuf); + else + fitLine3D( points.ptr(), npoints3, distType, + (float)param, (float)reps, (float)aeps, linebuf); + + Mat(npoints2 >= 0 ? 4 : 6, 1, CV_32F, linebuf).copyTo(_line); } @@ -638,82 +624,13 @@ CV_IMPL void cvFitLine( const CvArr* array, int dist, double param, double reps, double aeps, float *line ) { - cv::AutoBuffer buffer; + CV_Assert(line != 0); - schar* points = 0; - union { CvContour contour; CvSeq seq; } header; - CvSeqBlock block; - CvSeq* ptseq = (CvSeq*)array; - int type; + cv::AutoBuffer buf; + cv::Mat points = cv::cvarrToMat(array, false, false, 0, &buf); + cv::Mat linemat(points.checkVector(2) >= 0 ? 4 : 6, 1, CV_32F, line); - if( !line ) - CV_Error( CV_StsNullPtr, "NULL pointer to line parameters" ); - - if( CV_IS_SEQ(ptseq) ) - { - type = CV_SEQ_ELTYPE(ptseq); - if( ptseq->total == 0 ) - CV_Error( CV_StsBadSize, "The sequence has no points" ); - if( (type!=CV_32FC2 && type!=CV_32FC3 && type!=CV_32SC2 && type!=CV_32SC3) || - CV_ELEM_SIZE(type) != ptseq->elem_size ) - CV_Error( CV_StsUnsupportedFormat, - "Input sequence must consist of 2d points or 3d points" ); - } - else - { - CvMat* mat = (CvMat*)array; - type = CV_MAT_TYPE(mat->type); - if( !CV_IS_MAT(mat)) - CV_Error( CV_StsBadArg, "Input array is not a sequence nor matrix" ); - - if( !CV_IS_MAT_CONT(mat->type) || - (type!=CV_32FC2 && type!=CV_32FC3 && type!=CV_32SC2 && type!=CV_32SC3) || - (mat->width != 1 && mat->height != 1)) - CV_Error( CV_StsBadArg, - "Input array must be 1d continuous array of 2d or 3d points" ); - - ptseq = cvMakeSeqHeaderForArray( - CV_SEQ_KIND_GENERIC|type, sizeof(CvContour), CV_ELEM_SIZE(type), mat->data.ptr, - mat->width + mat->height - 1, &header.seq, &block ); - } - - if( reps < 0 || aeps < 0 ) - CV_Error( CV_StsOutOfRange, "Both reps and aeps must be non-negative" ); - - if( CV_MAT_DEPTH(type) == CV_32F && ptseq->first->next == ptseq->first ) - { - /* no need to copy data in this case */ - points = ptseq->first->data; - } - else - { - buffer.allocate(ptseq->total*CV_ELEM_SIZE(type)); - points = buffer; - cvCvtSeqToArray( ptseq, points, CV_WHOLE_SEQ ); - - if( CV_MAT_DEPTH(type) != CV_32F ) - { - int i, total = ptseq->total*CV_MAT_CN(type); - assert( CV_MAT_DEPTH(type) == CV_32S ); - - for( i = 0; i < total; i++ ) - ((float*)points)[i] = (float)((int*)points)[i]; - } - } - - if( dist == CV_DIST_USER ) - CV_Error( CV_StsBadArg, "User-defined distance is not allowed" ); - - if( CV_MAT_CN(type) == 2 ) - { - IPPI_CALL( icvFitLine2D( (CvPoint2D32f*)points, ptseq->total, - dist, (float)param, (float)reps, (float)aeps, line )); - } - else - { - IPPI_CALL( icvFitLine3D( (CvPoint3D32f*)points, ptseq->total, - dist, (float)param, (float)reps, (float)aeps, line )); - } + cv::fitLine(points, linemat, dist, param, reps, aeps); } /* End of file. */ diff --git a/modules/imgproc/src/matchcontours.cpp b/modules/imgproc/src/matchcontours.cpp index 44c01e02d..eca385900 100644 --- a/modules/imgproc/src/matchcontours.cpp +++ b/modules/imgproc/src/matchcontours.cpp @@ -40,159 +40,122 @@ //M*/ #include "precomp.hpp" -/*F/////////////////////////////////////////////////////////////////////////////////////// -// Name: cvMatchContours -// Purpose: -// Calculates matching of the two contours -// Context: -// Parameters: -// contour_1 - pointer to the first input contour object. -// contour_2 - pointer to the second input contour object. -// method - method for the matching calculation -// (now CV_IPPI_CONTOURS_MATCH_I1, CV_CONTOURS_MATCH_I2 or -// CV_CONTOURS_MATCH_I3 only ) -// rezult - output calculated measure -// -//F*/ -CV_IMPL double -cvMatchShapes( const void* contour1, const void* contour2, - int method, double /*parameter*/ ) + +double cv::matchShapes(InputArray contour1, InputArray contour2, int method, double) { - CvMoments moments; - CvHuMoments huMoments; double ma[7], mb[7]; int i, sma, smb; double eps = 1.e-5; double mmm; double result = 0; - if( !contour1 || !contour2 ) - CV_Error( CV_StsNullPtr, "" ); - - // calculate moments of the first shape - cvMoments( contour1, &moments ); - cvGetHuMoments( &moments, &huMoments ); - - ma[0] = huMoments.hu1; - ma[1] = huMoments.hu2; - ma[2] = huMoments.hu3; - ma[3] = huMoments.hu4; - ma[4] = huMoments.hu5; - ma[5] = huMoments.hu6; - ma[6] = huMoments.hu7; - - - // calculate moments of the second shape - cvMoments( contour2, &moments ); - cvGetHuMoments( &moments, &huMoments ); - - mb[0] = huMoments.hu1; - mb[1] = huMoments.hu2; - mb[2] = huMoments.hu3; - mb[3] = huMoments.hu4; - mb[4] = huMoments.hu5; - mb[5] = huMoments.hu6; - mb[6] = huMoments.hu7; + HuMoments( moments(contour1), ma ); + HuMoments( moments(contour2), mb ); switch (method) { case 1: + for( i = 0; i < 7; i++ ) { - for( i = 0; i < 7; i++ ) + double ama = fabs( ma[i] ); + double amb = fabs( mb[i] ); + + if( ma[i] > 0 ) + sma = 1; + else if( ma[i] < 0 ) + sma = -1; + else + sma = 0; + if( mb[i] > 0 ) + smb = 1; + else if( mb[i] < 0 ) + smb = -1; + else + smb = 0; + + if( ama > eps && amb > eps ) { - double ama = fabs( ma[i] ); - double amb = fabs( mb[i] ); - - if( ma[i] > 0 ) - sma = 1; - else if( ma[i] < 0 ) - sma = -1; - else - sma = 0; - if( mb[i] > 0 ) - smb = 1; - else if( mb[i] < 0 ) - smb = -1; - else - smb = 0; - - if( ama > eps && amb > eps ) - { - ama = 1. / (sma * log10( ama )); - amb = 1. / (smb * log10( amb )); - result += fabs( -ama + amb ); - } + ama = 1. / (sma * log10( ama )); + amb = 1. / (smb * log10( amb )); + result += fabs( -ama + amb ); } - break; } + break; case 2: + for( i = 0; i < 7; i++ ) { - for( i = 0; i < 7; i++ ) + double ama = fabs( ma[i] ); + double amb = fabs( mb[i] ); + + if( ma[i] > 0 ) + sma = 1; + else if( ma[i] < 0 ) + sma = -1; + else + sma = 0; + if( mb[i] > 0 ) + smb = 1; + else if( mb[i] < 0 ) + smb = -1; + else + smb = 0; + + if( ama > eps && amb > eps ) { - double ama = fabs( ma[i] ); - double amb = fabs( mb[i] ); - - if( ma[i] > 0 ) - sma = 1; - else if( ma[i] < 0 ) - sma = -1; - else - sma = 0; - if( mb[i] > 0 ) - smb = 1; - else if( mb[i] < 0 ) - smb = -1; - else - smb = 0; - - if( ama > eps && amb > eps ) - { - ama = sma * log10( ama ); - amb = smb * log10( amb ); - result += fabs( -ama + amb ); - } + ama = sma * log10( ama ); + amb = smb * log10( amb ); + result += fabs( -ama + amb ); } - break; } + break; case 3: + for( i = 0; i < 7; i++ ) { - for( i = 0; i < 7; i++ ) + double ama = fabs( ma[i] ); + double amb = fabs( mb[i] ); + + if( ma[i] > 0 ) + sma = 1; + else if( ma[i] < 0 ) + sma = -1; + else + sma = 0; + if( mb[i] > 0 ) + smb = 1; + else if( mb[i] < 0 ) + smb = -1; + else + smb = 0; + + if( ama > eps && amb > eps ) { - double ama = fabs( ma[i] ); - double amb = fabs( mb[i] ); - - if( ma[i] > 0 ) - sma = 1; - else if( ma[i] < 0 ) - sma = -1; - else - sma = 0; - if( mb[i] > 0 ) - smb = 1; - else if( mb[i] < 0 ) - smb = -1; - else - smb = 0; - - if( ama > eps && amb > eps ) - { - ama = sma * log10( ama ); - amb = smb * log10( amb ); - mmm = fabs( (ama - amb) / ama ); - if( result < mmm ) - result = mmm; - } + ama = sma * log10( ama ); + amb = smb * log10( amb ); + mmm = fabs( (ama - amb) / ama ); + if( result < mmm ) + result = mmm; } - break; } + break; default: CV_Error( CV_StsBadArg, "Unknown comparison method" ); } - + return result; } +CV_IMPL double +cvMatchShapes( const void* _contour1, const void* _contour2, + int method, double parameter ) +{ + cv::AutoBuffer abuf1, abuf2; + cv::Mat contour1 = cv::cvarrToMat(_contour1, false, false, 0, &abuf1); + cv::Mat contour2 = cv::cvarrToMat(_contour2, false, false, 0, &abuf2); + + return cv::matchShapes(contour1, contour2, method, parameter); +} + /* End of file. */ diff --git a/modules/imgproc/src/rotcalipers.cpp b/modules/imgproc/src/rotcalipers.cpp index 2171ec13e..b0ce3020b 100644 --- a/modules/imgproc/src/rotcalipers.cpp +++ b/modules/imgproc/src/rotcalipers.cpp @@ -7,7 +7,7 @@ // copy or use the software. // // -// Intel License Agreement +// License Agreement // For Open Source Computer Vision Library // // Copyright (C) 2000, Intel Corporation, all rights reserved. @@ -23,13 +23,13 @@ // 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 +// * The name of OpenCV Foundation 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, +// In no event shall the OpenCV Foundation 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 @@ -40,7 +40,10 @@ //M*/ #include "precomp.hpp" -typedef struct +namespace cv +{ + +struct MinAreaState { int bottom; int left; @@ -48,61 +51,58 @@ typedef struct float width; float base_a; float base_b; -} -icvMinAreaState; +}; -#define CV_CALIPERS_MAXHEIGHT 0 -#define CV_CALIPERS_MINAREARECT 1 -#define CV_CALIPERS_MAXDIST 2 +enum { CALIPERS_MAXHEIGHT=0, CALIPERS_MINAREARECT=1, CALIPERS_MAXDIST=2 }; /*F/////////////////////////////////////////////////////////////////////////////////////// -// Name: icvRotatingCalipers -// Purpose: -// Rotating calipers algorithm with some applications -// -// Context: -// Parameters: -// points - convex hull vertices ( any orientation ) -// n - number of vertices -// mode - concrete application of algorithm -// can be CV_CALIPERS_MAXDIST or -// CV_CALIPERS_MINAREARECT -// left, bottom, right, top - indexes of extremal points -// out - output info. -// In case CV_CALIPERS_MAXDIST it points to float value - -// maximal height of polygon. -// In case CV_CALIPERS_MINAREARECT -// ((CvPoint2D32f*)out)[0] - corner -// ((CvPoint2D32f*)out)[1] - vector1 -// ((CvPoint2D32f*)out)[0] - corner2 -// -// ^ -// | -// vector2 | -// | -// |____________\ -// corner / -// vector1 -// -// Returns: -// Notes: -//F*/ + // Name: rotatingCalipers + // Purpose: + // Rotating calipers algorithm with some applications + // + // Context: + // Parameters: + // points - convex hull vertices ( any orientation ) + // n - number of vertices + // mode - concrete application of algorithm + // can be CV_CALIPERS_MAXDIST or + // CV_CALIPERS_MINAREARECT + // left, bottom, right, top - indexes of extremal points + // out - output info. + // In case CV_CALIPERS_MAXDIST it points to float value - + // maximal height of polygon. + // In case CV_CALIPERS_MINAREARECT + // ((CvPoint2D32f*)out)[0] - corner + // ((CvPoint2D32f*)out)[1] - vector1 + // ((CvPoint2D32f*)out)[0] - corner2 + // + // ^ + // | + // vector2 | + // | + // |____________\ + // corner / + // vector1 + // + // Returns: + // Notes: + //F*/ /* we will use usual cartesian coordinates */ -static void -icvRotatingCalipers( CvPoint2D32f* points, int n, int mode, float* out ) +static void rotatingCalipers( const Point2f* points, int n, int mode, float* out ) { float minarea = FLT_MAX; float max_dist = 0; char buffer[32] = {}; int i, k; - CvPoint2D32f* vect = (CvPoint2D32f*)cvAlloc( n * sizeof(vect[0]) ); - float* inv_vect_length = (float*)cvAlloc( n * sizeof(inv_vect_length[0]) ); + AutoBuffer abuf(n*3); + float* inv_vect_length = abuf; + Point2f* vect = (Point2f*)(inv_vect_length + n); int left = 0, bottom = 0, right = 0, top = 0; int seq[4] = { -1, -1, -1, -1 }; /* rotating calipers sides will always have coordinates - (a,b) (-b,a) (-a,-b) (b, -a) + (a,b) (-b,a) (-a,-b) (b, -a) */ /* this is a first base bector (a,b) initialized by (1,0) */ float orientation = 0; @@ -110,7 +110,7 @@ icvRotatingCalipers( CvPoint2D32f* points, int n, int mode, float* out ) float base_b = 0; float left_x, right_x, top_y, bottom_y; - CvPoint2D32f pt0 = points[0]; + Point2f pt0 = points[0]; left_x = right_x = pt0.x; top_y = bottom_y = pt0.y; @@ -131,7 +131,7 @@ icvRotatingCalipers( CvPoint2D32f* points, int n, int mode, float* out ) if( pt0.y < bottom_y ) bottom_y = pt0.y, bottom = i; - CvPoint2D32f pt = points[(i+1) & (i+1 < n ? -1 : 0)]; + Point2f pt = points[(i+1) & (i+1 < n ? -1 : 0)]; dx = pt.x - pt0.x; dy = pt.y - pt0.y; @@ -143,9 +143,7 @@ icvRotatingCalipers( CvPoint2D32f* points, int n, int mode, float* out ) pt0 = pt; } - //cvbInvSqrt( inv_vect_length, inv_vect_length, n ); - - /* find convex hull orientation */ + // find convex hull orientation { double ax = vect[n-1].x; double ay = vect[n-1].y; @@ -165,18 +163,18 @@ icvRotatingCalipers( CvPoint2D32f* points, int n, int mode, float* out ) ax = bx; ay = by; } - assert( orientation != 0 ); + CV_Assert( orientation != 0 ); } base_a = orientation; -/*****************************************************************************************/ -/* init calipers position */ + /*****************************************************************************************/ + /* init calipers position */ seq[0] = bottom; seq[1] = right; seq[2] = top; seq[3] = left; -/*****************************************************************************************/ -/* Main loop - evaluate angles and rotate calipers */ + /*****************************************************************************************/ + /* Main loop - evaluate angles and rotate calipers */ /* all of edges will be checked while rotating calipers by 90 degrees */ for( k = 0; k < n; k++ ) @@ -229,190 +227,140 @@ icvRotatingCalipers( CvPoint2D32f* points, int n, int mode, float* out ) base_a = -lead_y; base_b = lead_x; break; - default: assert(0); + default: + CV_Error(CV_StsError, "main_element should be 0, 1, 2 or 3"); } } /* change base point of main edge */ seq[main_element] += 1; seq[main_element] = (seq[main_element] == n) ? 0 : seq[main_element]; - switch (mode) { - case CV_CALIPERS_MAXHEIGHT: + case CALIPERS_MAXHEIGHT: { - /* now main element lies on edge alligned to calipers side */ + /* now main element lies on edge alligned to calipers side */ - /* find opposite element i.e. transform */ - /* 0->2, 1->3, 2->0, 3->1 */ - int opposite_el = main_element ^ 2; + /* find opposite element i.e. transform */ + /* 0->2, 1->3, 2->0, 3->1 */ + int opposite_el = main_element ^ 2; - float dx = points[seq[opposite_el]].x - points[seq[main_element]].x; - float dy = points[seq[opposite_el]].y - points[seq[main_element]].y; - float dist; + float dx = points[seq[opposite_el]].x - points[seq[main_element]].x; + float dy = points[seq[opposite_el]].y - points[seq[main_element]].y; + float dist; - if( main_element & 1 ) - dist = (float)fabs(dx * base_a + dy * base_b); - else - dist = (float)fabs(dx * (-base_b) + dy * base_a); + if( main_element & 1 ) + dist = (float)fabs(dx * base_a + dy * base_b); + else + dist = (float)fabs(dx * (-base_b) + dy * base_a); - if( dist > max_dist ) - max_dist = dist; - - break; + if( dist > max_dist ) + max_dist = dist; } - case CV_CALIPERS_MINAREARECT: + break; + case CALIPERS_MINAREARECT: /* find area of rectangle */ { - float height; - float area; + float height; + float area; - /* find vector left-right */ - float dx = points[seq[1]].x - points[seq[3]].x; - float dy = points[seq[1]].y - points[seq[3]].y; + /* find vector left-right */ + float dx = points[seq[1]].x - points[seq[3]].x; + float dy = points[seq[1]].y - points[seq[3]].y; - /* dotproduct */ - float width = dx * base_a + dy * base_b; + /* dotproduct */ + float width = dx * base_a + dy * base_b; - /* find vector left-right */ - dx = points[seq[2]].x - points[seq[0]].x; - dy = points[seq[2]].y - points[seq[0]].y; + /* find vector left-right */ + dx = points[seq[2]].x - points[seq[0]].x; + dy = points[seq[2]].y - points[seq[0]].y; - /* dotproduct */ - height = -dx * base_b + dy * base_a; + /* dotproduct */ + height = -dx * base_b + dy * base_a; - area = width * height; - if( area <= minarea ) - { - float *buf = (float *) buffer; + area = width * height; + if( area <= minarea ) + { + float *buf = (float *) buffer; - minarea = area; - /* leftist point */ - ((int *) buf)[0] = seq[3]; - buf[1] = base_a; - buf[2] = width; - buf[3] = base_b; - buf[4] = height; - /* bottom point */ - ((int *) buf)[5] = seq[0]; - buf[6] = area; - } - break; + minarea = area; + /* leftist point */ + ((int *) buf)[0] = seq[3]; + buf[1] = base_a; + buf[2] = width; + buf[3] = base_b; + buf[4] = height; + /* bottom point */ + ((int *) buf)[5] = seq[0]; + buf[6] = area; } + } + break; } /*switch */ } /* for */ switch (mode) { - case CV_CALIPERS_MINAREARECT: + case CALIPERS_MINAREARECT: { - float *buf = (float *) buffer; + float *buf = (float *) buffer; - float A1 = buf[1]; - float B1 = buf[3]; + float A1 = buf[1]; + float B1 = buf[3]; - float A2 = -buf[3]; - float B2 = buf[1]; + float A2 = -buf[3]; + float B2 = buf[1]; - float C1 = A1 * points[((int *) buf)[0]].x + points[((int *) buf)[0]].y * B1; - float C2 = A2 * points[((int *) buf)[5]].x + points[((int *) buf)[5]].y * B2; + float C1 = A1 * points[((int *) buf)[0]].x + points[((int *) buf)[0]].y * B1; + float C2 = A2 * points[((int *) buf)[5]].x + points[((int *) buf)[5]].y * B2; - float idet = 1.f / (A1 * B2 - A2 * B1); + float idet = 1.f / (A1 * B2 - A2 * B1); - float px = (C1 * B2 - C2 * B1) * idet; - float py = (A1 * C2 - A2 * C1) * idet; + float px = (C1 * B2 - C2 * B1) * idet; + float py = (A1 * C2 - A2 * C1) * idet; - out[0] = px; - out[1] = py; - - out[2] = A1 * buf[2]; - out[3] = B1 * buf[2]; - - out[4] = A2 * buf[4]; - out[5] = B2 * buf[4]; + out[0] = px; + out[1] = py; + + out[2] = A1 * buf[2]; + out[3] = B1 * buf[2]; + + out[4] = A2 * buf[4]; + out[5] = B2 * buf[4]; } break; - case CV_CALIPERS_MAXHEIGHT: + case CALIPERS_MAXHEIGHT: { - out[0] = max_dist; + out[0] = max_dist; } break; } - - cvFree( &vect ); - cvFree( &inv_vect_length ); +} + } -CV_IMPL CvBox2D -cvMinAreaRect2( const CvArr* array, CvMemStorage* storage ) +cv::RotatedRect cv::minAreaRect( InputArray _points ) { - cv::Ptr temp_storage; - CvBox2D box; - cv::AutoBuffer _points; - CvPoint2D32f* points; - - memset(&box, 0, sizeof(box)); - - int i, n; - CvSeqReader reader; - CvContour contour_header; - CvSeqBlock block; - CvSeq* ptseq = (CvSeq*)array; - CvPoint2D32f out[3]; - - if( CV_IS_SEQ(ptseq) ) + Mat hull; + Point2f out[3]; + RotatedRect box; + + convexHull(_points, hull, true, true); + + if( hull.depth() != CV_32F ) { - if( !CV_IS_SEQ_POINT_SET(ptseq) && - (CV_SEQ_KIND(ptseq) != CV_SEQ_KIND_CURVE || - CV_SEQ_ELTYPE(ptseq) != CV_SEQ_ELTYPE_PPOINT )) - CV_Error( CV_StsUnsupportedFormat, - "Input sequence must consist of 2d points or pointers to 2d points" ); - if( !storage ) - storage = ptseq->storage; + Mat temp; + hull.convertTo(temp, CV_32F); + hull = temp; } - else - { - ptseq = cvPointSeqFromMat( CV_SEQ_KIND_GENERIC, array, &contour_header, &block ); - } - - if( storage ) - { - temp_storage = cvCreateChildMemStorage( storage ); - } - else - { - temp_storage = cvCreateMemStorage(1 << 10); - } - - ptseq = cvConvexHull2( ptseq, temp_storage, CV_CLOCKWISE, 1 ); - n = ptseq->total; - - _points.allocate(n); - points = _points; - cvStartReadSeq( ptseq, &reader ); - - if( CV_SEQ_ELTYPE( ptseq ) == CV_32SC2 ) - { - for( i = 0; i < n; i++ ) - { - CvPoint pt; - CV_READ_SEQ_ELEM( pt, reader ); - points[i].x = (float)pt.x; - points[i].y = (float)pt.y; - } - } - else - { - for( i = 0; i < n; i++ ) - { - CV_READ_SEQ_ELEM( points[i], reader ); - } - } - + + int n = hull.checkVector(2); + const Point2f* hpoints = (const Point2f*)hull.data; + if( n > 2 ) { - icvRotatingCalipers( points, n, CV_CALIPERS_MINAREARECT, (float*)out ); + rotatingCalipers( hpoints, n, CALIPERS_MINAREARECT, (float*)out ); box.center.x = out[0].x + (out[1].x + out[2].x)*0.5f; box.center.y = out[0].y + (out[1].y + out[2].y)*0.5f; box.size.width = (float)sqrt((double)out[1].x*out[1].x + (double)out[1].y*out[1].y); @@ -421,10 +369,10 @@ cvMinAreaRect2( const CvArr* array, CvMemStorage* storage ) } else if( n == 2 ) { - box.center.x = (points[0].x + points[1].x)*0.5f; - box.center.y = (points[0].y + points[1].y)*0.5f; - double dx = points[1].x - points[0].x; - double dy = points[1].y - points[0].y; + box.center.x = (hpoints[0].x + hpoints[1].x)*0.5f; + box.center.y = (hpoints[0].y + hpoints[1].y)*0.5f; + double dx = hpoints[1].x - hpoints[0].x; + double dy = hpoints[1].y - hpoints[0].y; box.size.width = (float)sqrt(dx*dx + dy*dy); box.size.height = 0; box.angle = (float)atan2( dy, dx ); @@ -432,10 +380,21 @@ cvMinAreaRect2( const CvArr* array, CvMemStorage* storage ) else { if( n == 1 ) - box.center = points[0]; + box.center = hpoints[0]; } - + box.angle = (float)(box.angle*180/CV_PI); return box; } + +CV_IMPL CvBox2D +cvMinAreaRect2( const CvArr* array, CvMemStorage* /*storage*/ ) +{ + cv::AutoBuffer abuf; + cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf); + + cv::RotatedRect rr = cv::minAreaRect(points); + return (CvBox2D)rr; +} + diff --git a/modules/imgproc/src/shapedescr.cpp b/modules/imgproc/src/shapedescr.cpp index b1bc1babe..a506535af 100644 --- a/modules/imgproc/src/shapedescr.cpp +++ b/modules/imgproc/src/shapedescr.cpp @@ -40,6 +40,912 @@ //M*/ #include "precomp.hpp" +namespace cv +{ + +static int intersectLines( double x1, double dx1, double y1, double dy1, + double x2, double dx2, double y2, double dy2, double *t2 ) +{ + double d = dx1 * dy2 - dx2 * dy1; + int result = -1; + + if( d != 0 ) + { + *t2 = ((x2 - x1) * dy1 - (y2 - y1) * dx1) / d; + result = 0; + } + return result; +} + +static bool findCircle( Point2f pt0, Point2f pt1, Point2f pt2, + Point2f* 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; + + if( intersectLines( x1, dx1, y1, dy1, x2, dx2, y2, dy2, &t ) >= 0 ) + { + center->x = (float) (x2 + dx2 * t); + center->y = (float) (y2 + dy2 * t); + *radius = (float)norm(*center - pt0); + return true; + } + + center->x = center->y = 0.f; + radius = 0; + return false; +} + + +static double pointInCircle( Point2f pt, Point2f 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 findEnslosingCicle4pts_32f( Point2f* pts, Point2f& _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; + Point2f center; + Point2f min_center; + float radius, min_radius = FLT_MAX; + Point2f 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 = (float)norm(pts[i] - pts[j]); + + if( max_dist < dist ) + { + max_dist = dist; + idxs[0] = i; + idxs[1] = j; + } + } + + if( max_dist > 0 ) + { + 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 = Point2f( (pts[idxs[0]].x + pts[idxs[1]].x)*0.5f, + (pts[idxs[0]].y + pts[idxs[1]].y)*0.5f ); + radius = (float)(norm(pts[idxs[0]] - center)*1.03); + if( radius < 1.f ) + radius = 1.f; + + if( pointInCircle( pts[idxs[2]], center, radius ) >= 0 && + pointInCircle( pts[idxs[3]], center, radius ) >= 0 ) + { + k = 2; //rand()%2+2; + } + else + { + mi = -1; + for( i = 0; i < 4; i++ ) + { + if( findCircle( pts[shuffles[i][0]], pts[shuffles[i][1]], + pts[shuffles[i][2]], ¢er, &radius ) ) + { + radius *= 1.03f; + if( radius < 2.f ) + radius = 2.f; + + if( pointInCircle( pts[shuffles[i][3]], center, radius ) >= 0 && + min_radius > radius ) + { + min_radius = radius; + min_center = center; + mi = i; + } + } + } + CV_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]; + } + } + + _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]; + CV_Assert( pointInCircle( pts[i], center, radius ) >= 0 ); + } + + return k; +} + +} + +void cv::minEnclosingCircle( InputArray _points, Point2f& _center, float& _radius ) +{ + int max_iters = 100; + const float eps = FLT_EPSILON*2; + bool result = false; + Mat points = _points.getMat(); + int i, j, k, count = points.checkVector(2); + int depth = points.depth(); + Point2f center; + float radius = 0.f; + CV_Assert(count >= 0 && (depth == CV_32F || depth == CV_32S)); + + _center.x = _center.y = 0.f; + _radius = 0.f; + + if( count == 0 ) + return; + + bool is_float = depth == CV_32F; + const Point* ptsi = (const Point*)points.data; + const Point2f* ptsf = (const Point2f*)points.data; + + Point2f pt = is_float ? ptsf[0] : Point2f((float)ptsi[0].x,(float)ptsi[0].y); + Point2f pts[4] = {pt, pt, pt, pt}; + + for( i = 1; i < count; i++ ) + { + pt = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); + + if( pt.x < pts[0].x ) + pts[0] = pt; + if( pt.x > pts[1].x ) + pts[1] = pt; + if( pt.y < pts[2].y ) + pts[2] = pt; + if( pt.y > pts[3].y ) + pts[3] = pt; + } + + for( k = 0; k < max_iters; k++ ) + { + double min_delta = 0, delta; + Point2f farAway(0,0); + /*only for first iteration because the alg is repared at the loop's foot*/ + if( k == 0 ) + findEnslosingCicle4pts_32f( pts, center, radius ); + + for( i = 0; i < count; i++ ) + { + pt = is_float ? ptsf[i] : Point2f((float)ptsi[i].x,(float)ptsi[i].y); + + delta = pointInCircle( pt, center, radius ); + if( delta < min_delta ) + { + min_delta = delta; + farAway = pt; + } + } + result = min_delta >= 0; + if( result ) + break; + + Point2f ptsCopy[4]; + // find good replacement partner for the point which is at most far away, + // starting with the one that lays in the actual circle (i=3) + for( i = 3; i >= 0; i-- ) + { + for( j = 0; j < 4; j++ ) + ptsCopy[j] = i != j ? pts[j] : farAway; + + findEnslosingCicle4pts_32f( ptsCopy, center, radius ); + if( pointInCircle( pts[i], center, radius ) >= 0) + { + // replaced one again in the new circle? + pts[i] = farAway; + break; + } + } + } + + if( !result ) + { + radius = 0.f; + for( i = 0; i < count; i++ ) + { + pt = is_float ? ptsf[i] : Point2f((float)ptsi[i].x,(float)ptsi[i].y); + float dx = center.x - pt.x, dy = center.y - pt.y; + float t = dx*dx + dy*dy; + radius = MAX(radius, t); + } + + radius = (float)(sqrt(radius)*(1 + eps)); + } + + _center = center; + _radius = radius; +} + + +// calculates length of a curve (e.g. contour perimeter) +double cv::arcLength( InputArray _curve, bool is_closed ) +{ + Mat curve = _curve.getMat(); + int count = curve.checkVector(2); + int depth = curve.depth(); + CV_Assert( count >= 0 && (depth == CV_32F || depth == CV_32S)); + double perimeter = 0; + + int i, j = 0; + const int N = 16; + float buf[N]; + + if( count <= 1 ) + return 0.; + + bool is_float = depth == CV_32F; + int last = is_closed ? count-1 : 0; + const Point* pti = (const Point*)curve.data; + const Point2f* ptf = (const Point2f*)curve.data; + + Point2f prev = is_float ? ptf[last] : Point2f((float)pti[last].x,(float)pti[last].y); + + for( i = 0; i < count; i++ ) + { + Point2f p = is_float ? ptf[i] : Point2f((float)pti[i].x,(float)pti[i].y); + float dx = p.x - prev.x, dy = p.y - prev.y; + buf[j] = dx*dx + dy*dy; + + if( ++j == N || i == count-1 ) + { + Mat bufmat(1, j, CV_32F, buf); + sqrt(bufmat, bufmat); + for( ; j > 0; j-- ) + perimeter += buf[j-1]; + } + prev = p; + } + + return perimeter; +} + +// area of a whole sequence +double cv::contourArea( InputArray _contour, bool oriented ) +{ + Mat contour = _contour.getMat(); + int npoints = contour.checkVector(2); + int depth = contour.depth(); + CV_Assert(npoints >= 0 && (depth == CV_32F || depth == CV_32S)); + + if( npoints == 0 ) + return 0.; + + double a00 = 0; + bool is_float = depth == CV_32F; + const Point* ptsi = (const Point*)contour.data; + const Point2f* ptsf = (const Point2f*)contour.data; + Point2f prev = is_float ? ptsf[npoints-1] : Point2f((float)ptsi[npoints-1].x, (float)ptsi[npoints-1].y); + + for( int i = 0; i < npoints; i++ ) + { + Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); + a00 += (double)prev.x * p.y - (double)prev.y * p.x; + prev = p; + } + + a00 *= 0.5; + if( !oriented ) + a00 = fabs(a00); + + return a00; +} + + +cv::RotatedRect cv::fitEllipse( InputArray _points ) +{ + Mat points = _points.getMat(); + int i, n = points.checkVector(2); + int depth = points.depth(); + CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S)); + + RotatedRect box; + + if( n < 5 ) + CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" ); + + // New fitellipse algorithm, contributed by Dr. Daniel Weiss + Point2f c(0,0); + double gfp[5], rp[5], t; + const double min_eps = 1e-6; + bool is_float = depth == CV_32F; + const Point* ptsi = (const Point*)points.data; + const Point2f* ptsf = (const Point2f*)points.data; + + AutoBuffer _Ad(n*5), _bd(n); + double *Ad = _Ad, *bd = _bd; + + // first fit for parameters A - E + Mat A( n, 5, CV_64F, Ad ); + Mat b( n, 1, CV_64F, bd ); + Mat x( 5, 1, CV_64F, gfp ); + + for( i = 0; i < n; i++ ) + { + Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); + c += p; + } + c.x /= n; + c.y /= n; + + for( i = 0; i < n; i++ ) + { + Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); + p -= c; + + 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; + } + + solve(A, b, x, DECOMP_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 = Mat( 2, 2, CV_64F, Ad ); + b = Mat( 2, 1, CV_64F, bd ); + x = Mat( 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]; + solve( A, b, x, DECOMP_SVD ); + + // re-fit for parameters A - C with those center coordinates + A = Mat( n, 3, CV_64F, Ad ); + b = Mat( n, 1, CV_64F, bd ); + x = Mat( 3, 1, CV_64F, gfp ); + for( i = 0; i < n; i++ ) + { + Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y); + p -= c; + 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]); + } + solve(A, b, x, DECOMP_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] + c.x; + box.center.y = (float)rp[1] + c.y; + 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; + + return box; +} + + +namespace cv +{ + +// Calculates bounding rectagnle of a point set or retrieves already calculated +static Rect pointSetBoundingRect( const Mat& points ) +{ + int npoints = points.checkVector(2); + int depth = points.depth(); + CV_Assert(npoints >= 0 && (depth == CV_32F || depth == CV_32S)); + + int xmin = 0, ymin = 0, xmax = -1, ymax = -1, i; + bool is_float = depth == CV_32F; + + if( npoints == 0 ) + return Rect(); + + const Point* pts = (const Point*)points.data; + Point pt = pts[0]; + +#if CV_SSE4_2 + if(cv::checkHardwareSupport(CV_CPU_SSE4_2)) + { + if( !is_float ) + { + __m128i minval, maxval; + minval = maxval = _mm_loadl_epi64((const __m128i*)(&pt)); //min[0]=pt.x, min[1]=pt.y + + for( i = 1; i < npoints; i++ ) + { + __m128i ptXY = _mm_loadl_epi64((const __m128i*)&pts[i]); + minval = _mm_min_epi32(ptXY, minval); + maxval = _mm_max_epi32(ptXY, maxval); + } + xmin = _mm_cvtsi128_si32(minval); + ymin = _mm_cvtsi128_si32(_mm_srli_si128(minval, 4)); + xmax = _mm_cvtsi128_si32(maxval); + ymax = _mm_cvtsi128_si32(_mm_srli_si128(maxval, 4)); + } + else + { + __m128 minvalf, maxvalf, z = _mm_setzero_ps(), ptXY = _mm_setzero_ps(); + minvalf = maxvalf = _mm_loadl_pi(z, (const __m64*)(&pt)); + + for( i = 1; i < npoints; i++ ) + { + ptXY = _mm_loadl_pi(ptXY, (const __m64*)&pts[i]); + + minvalf = _mm_min_ps(minvalf, ptXY); + maxvalf = _mm_max_ps(maxvalf, ptXY); + } + + float xyminf[2], xymaxf[2]; + _mm_storel_pi((__m64*)xyminf, minvalf); + _mm_storel_pi((__m64*)xymaxf, maxvalf); + xmin = cvFloor(xyminf[0]); + ymin = cvFloor(xyminf[1]); + xmax = cvFloor(xymaxf[0]); + ymax = cvFloor(xymaxf[1]); + } + } + else +#endif + { + if( !is_float ) + { + xmin = xmax = pt.x; + ymin = ymax = pt.y; + + for( i = 1; i < npoints; i++ ) + { + pt = pts[i]; + + 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 + { + Cv32suf v; + // init values + xmin = xmax = CV_TOGGLE_FLT(pt.x); + ymin = ymax = CV_TOGGLE_FLT(pt.y); + + for( i = 1; i < npoints; i++ ) + { + pt = pts[i]; + 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); + } + } + + return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1); +} + + +static Rect maskBoundingRect( const Mat& img ) +{ + CV_Assert( img.depth() <= CV_8S && img.channels() == 1 ); + + Size size = img.size(); + int xmin = size.width, ymin = -1, xmax = -1, ymax = -1, i, j, k; + + for( i = 0; i < size.height; i++ ) + { + const uchar* _ptr = img.ptr(i); + const uchar* ptr = (const uchar*)alignPtr(_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; + return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1); +} + +} + +cv::Rect cv::boundingRect(InputArray array) +{ + Mat m = array.getMat(); + return m.depth() <= CV_8U ? maskBoundingRect(m) : pointSetBoundingRect(m); +} + +////////////////////////////////////////////// C API /////////////////////////////////////////// + +CV_IMPL int +cvMinEnclosingCircle( const void* array, CvPoint2D32f * _center, float *_radius ) +{ + cv::AutoBuffer abuf; + cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf); + cv::Point2f center; + float radius; + + cv::minEnclosingCircle(points, center, radius); + if(_center) + *_center = center; + if(_radius) + *_radius = radius; + return 1; +} + +static void +icvMemCopy( double **buf1, double **buf2, double **buf3, int *b_max ) +{ + CV_Assert( (*buf1 != NULL || *buf2 != NULL) && *buf3 != NULL ); + + int bb = *b_max; + if( *buf2 == NULL ) + { + *b_max = 2 * (*b_max); + *buf2 = (double *)cvAlloc( (*b_max) * sizeof( double )); + + 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 )); + + memcpy( *buf1, *buf3, bb * sizeof( double )); + + *buf3 = *buf1; + cvFree( buf2 ); + *buf2 = NULL; + } +} + + +/* area of a contour sector */ +static double icvContourSecArea( CvSeq * contour, CvSlice slice ) +{ + 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; + double area = 0; + + CV_Assert( contour != NULL && CV_IS_SEQ_POINT_SET( contour )); + + lpt = cvSliceLength( slice, contour ); + /*if( n2 >= n1 ) + lpt = n2 - n1 + 1; + else + lpt = contour->total - n1 + n2 + 1;*/ + + if( contour->total <= 0 || lpt <= 2 ) + return 0.; + + a00 = x0 = y0 = xi_1 = yi_1 = 0; + sk1 = 0; + flag = 0; + dxy = 0; + p_are1 = (double *) cvAlloc( p_max * sizeof( double )); + + 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 area; +} + + +/* external contour area function */ +CV_IMPL double +cvContourArea( const void *array, CvSlice slice, int oriented ) +{ + double area = 0; + + 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 + { + contour = cvPointSeqFromMat( CV_SEQ_KIND_CURVE, array, &contour_header, &block ); + } + + if( cvSliceLength( slice, contour ) == contour->total ) + { + cv::AutoBuffer abuf; + cv::Mat points = cv::cvarrToMat(contour, false, false, 0, &abuf); + return cv::contourArea( points, oriented !=0 ); + } + + if( CV_SEQ_ELTYPE( contour ) != CV_32SC2 ) + CV_Error( CV_StsUnsupportedFormat, + "Only curves with integer coordinates are supported in case of contour slice" ); + area = icvContourSecArea( contour, slice ); + return oriented ? area : fabs(area); +} + + /* calculates length of a curve (e.g. contour perimeter) */ CV_IMPL double cvArcLength( const void *array, CvSlice slice, int is_closed ) @@ -67,8 +973,8 @@ cvArcLength( const void *array, CvSlice slice, int is_closed ) { is_closed = is_closed > 0; contour = cvPointSeqFromMat( - CV_SEQ_KIND_CURVE | (is_closed ? CV_SEQ_FLAG_CLOSED : 0), - array, &contour_header, &block ); + CV_SEQ_KIND_CURVE | (is_closed ? CV_SEQ_FLAG_CLOSED : 0), + array, &contour_header, &block ); } if( contour->total > 1 ) @@ -81,7 +987,7 @@ cvArcLength( const void *array, CvSlice slice, int is_closed ) count -= !is_closed && count == contour->total; - /* scroll the reader by 1 point */ + // scroll the reader by 1 point reader.prev_elem = reader.ptr; CV_NEXT_SEQ_ELEM( sizeof(CvPoint), reader ); @@ -123,824 +1029,29 @@ cvArcLength( const void *array, CvSlice slice, int is_closed ) } } } - + 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; - - CvSeqReader reader; - int 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 - { - sequence = cvPointSeqFromMat( - CV_SEQ_KIND_GENERIC, array, &contour_header, &block ); - } - - if( sequence->total <= 0 ) - CV_Error( CV_StsBadSize, "" ); - - 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(int 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(int 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, farAway = { 0, 0}; - /*only for first iteration because the alg is repared at the loop's foot*/ - if(k==0) - icvFindEnslosingCicle4pts_32f( pts, ¢er, &radius ); - - cvStartReadSeq( sequence, &reader, 0 ); - - for(int 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; - farAway = ptfl; - } - } - result = min_delta >= 0; - if( result ) - break; - - CvPoint2D32f ptsCopy[4]; - /* find good replacement partner for the point which is at most far away, - starting with the one that lays in the actual circle (i=3) */ - for(int i = 3; i >=0; i-- ) - { - for(int j = 0; j < 4; j++ ) - { - ptsCopy[j]=(i != j)? pts[j]: farAway; - } - - icvFindEnslosingCicle4pts_32f(ptsCopy, ¢er, &radius ); - if( icvIsPtInCircle( pts[i], center, radius )>=0){ // replaced one again in the new circle? - pts[i] = farAway; - break; - } - } - } - - if( !result ) - { - cvStartReadSeq( sequence, &reader, 0 ); - radius = 0.f; - - for(int 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; - } - - *_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_POINT_SET( 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, int oriented ) -{ - double area = 0; - - 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 - { - 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 )); - } - - return oriented ? area : fabs(area); -} - - CV_IMPL CvBox2D cvFitEllipse2( const CvArr* array ) { - CvBox2D box; - cv::AutoBuffer Ad, bd; - memset( &box, 0, sizeof(box)); - - 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 - { - 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 >= 5" ); - - /* - * New fitellipse algorithm, contributed by Dr. Daniel Weiss - */ - CvPoint2D32f c = {0,0}; - double gfp[5], rp[5], t; - CvMat A, b, x; - const double min_eps = 1e-6; - int i, is_float; - CvSeqReader reader; - - Ad.allocate(n*5); - bd.allocate(n); - - // 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 ); - c.x += p.x; - c.y += p.y; - } - c.x /= n; - c.y /= n; - - 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 ); - p.x -= c.x; - p.y -= c.y; - - 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 ); - p.x -= c.x; - p.y -= c.y; - 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] + c.x; - box.center.y = (float)rp[1] + c.y; - 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; - - return box; + cv::AutoBuffer abuf; + cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf); + return cv::fitEllipse(points); } - /* 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; CvMat stub, *mat = 0; - int xmin = 0, ymin = 0, xmax = -1, ymax = -1, i, j, k; int calculate = update; if( CV_IS_SEQ( array )) @@ -977,210 +1088,17 @@ cvBoundingRect( CvArr* array, int update ) 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; + rect = cv::maskBoundingRect(cv::cvarrToMat(mat)); } else if( ptseq->total ) { - int is_float = CV_SEQ_ELTYPE(ptseq) == CV_32FC2; - cvStartReadSeq( ptseq, &reader, 0 ); - CvPoint pt; - CV_READ_SEQ_ELEM( pt, reader ); - #if CV_SSE4_2 - if(cv::checkHardwareSupport(CV_CPU_SSE4_2)) - { - if( !is_float ) - { - __m128i minval, maxval; - minval = maxval = _mm_loadl_epi64((const __m128i*)(&pt)); //min[0]=pt.x, min[1]=pt.y - - for( i = 1; i < ptseq->total; i++) - { - __m128i ptXY = _mm_loadl_epi64((const __m128i*)(reader.ptr)); - CV_NEXT_SEQ_ELEM(sizeof(pt), reader); - minval = _mm_min_epi32(ptXY, minval); - maxval = _mm_max_epi32(ptXY, maxval); - } - xmin = _mm_cvtsi128_si32(minval); - ymin = _mm_cvtsi128_si32(_mm_srli_si128(minval, 4)); - xmax = _mm_cvtsi128_si32(maxval); - ymax = _mm_cvtsi128_si32(_mm_srli_si128(maxval, 4)); - } - else - { - __m128 minvalf, maxvalf, z = _mm_setzero_ps(), ptXY = _mm_setzero_ps(); - minvalf = maxvalf = _mm_loadl_pi(z, (const __m64*)(&pt)); - - for( i = 1; i < ptseq->total; i++ ) - { - ptXY = _mm_loadl_pi(ptXY, (const __m64*)reader.ptr); - CV_NEXT_SEQ_ELEM(sizeof(pt), reader); - - minvalf = _mm_min_ps(minvalf, ptXY); - maxvalf = _mm_max_ps(maxvalf, ptXY); - } - - float xyminf[2], xymaxf[2]; - _mm_storel_pi((__m64*)xyminf, minvalf); - _mm_storel_pi((__m64*)xymaxf, maxvalf); - xmin = cvFloor(xyminf[0]); - ymin = cvFloor(xyminf[1]); - xmax = cvFloor(xymaxf[0]); - ymax = cvFloor(xymaxf[1]); - } - } - else - #endif - { - if( !is_float ) - { - 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 - { - Cv32suf v; - // init values - 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; + cv::AutoBuffer abuf; + rect = cv::pointSetBoundingRect(cv::cvarrToMat(ptseq, false, false, 0, &abuf)); } if( update ) ((CvContour*)ptseq)->rect = rect; return rect; } + /* End of file. */ diff --git a/modules/video/src/optflowgf.cpp b/modules/video/src/optflowgf.cpp index 0f63bb47a..650e24262 100644 --- a/modules/video/src/optflowgf.cpp +++ b/modules/video/src/optflowgf.cpp @@ -395,8 +395,8 @@ FarnebackUpdateFlow_GaussianBlur( const Mat& _R0, const Mat& _R1, double sigma = m*0.3, s = 1; AutoBuffer _vsum((width+m*2+2)*5 + 16), _hsum(width*5 + 16); - AutoBuffer _kernel((m+1)*5 + 16); - AutoBuffer _srow(m*2+1); + AutoBuffer _kernel((m+1)*5 + 16); + AutoBuffer _srow(m*2+1); float *vsum = alignPtr((float*)_vsum + (m+1)*5, 16), *hsum = alignPtr((float*)_hsum, 16); float* kernel = (float*)_kernel; const float** srow = (const float**)&_srow[0]; diff --git a/samples/cpp/minarea.cpp b/samples/cpp/minarea.cpp index 6056c39c4..eca5dc1d3 100644 --- a/samples/cpp/minarea.cpp +++ b/samples/cpp/minarea.cpp @@ -13,7 +13,7 @@ static void help() "Random points are generated and then enclosed.\n" "Call:\n" "./minarea\n" - "Using OpenCV version %s\n" << CV_VERSION << "\n" << endl; + "Using OpenCV v" << CV_VERSION << "\n" << endl; } int main( int /*argc*/, char** /*argv*/ )