shapedescr.cpp

00001 /*M///////////////////////////////////////////////////////////////////////////////////////
00002 //
00003 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
00004 //
00005 //  By downloading, copying, installing or using the software you agree to this license.
00006 //  If you do not agree to this license, do not download, install,
00007 //  copy or use the software.
00008 //
00009 //
00010 //                        Intel License Agreement
00011 //                For Open Source Computer Vision Library
00012 //
00013 // Copyright (C) 2000, Intel Corporation, all rights reserved.
00014 // Third party copyrights are property of their respective owners.
00015 //
00016 // Redistribution and use in source and binary forms, with or without modification,
00017 // are permitted provided that the following conditions are met:
00018 //
00019 //   * Redistribution's of source code must retain the above copyright notice,
00020 //     this list of conditions and the following disclaimer.
00021 //
00022 //   * Redistribution's in binary form must reproduce the above copyright notice,
00023 //     this list of conditions and the following disclaimer in the documentation
00024 //     and/or other materials provided with the distribution.
00025 //
00026 //   * The name of Intel Corporation may not be used to endorse or promote products
00027 //     derived from this software without specific prior written permission.
00028 //
00029 // This software is provided by the copyright holders and contributors "as is" and
00030 // any express or implied warranties, including, but not limited to, the implied
00031 // warranties of merchantability and fitness for a particular purpose are disclaimed.
00032 // In no event shall the Intel Corporation or contributors be liable for any direct,
00033 // indirect, incidental, special, exemplary, or consequential damages
00034 // (including, but not limited to, procurement of substitute goods or services;
00035 // loss of use, data, or profits; or business interruption) however caused
00036 // and on any theory of liability, whether in contract, strict liability,
00037 // or tort (including negligence or otherwise) arising in any way out of
00038 // the use of this software, even if advised of the possibility of such damage.
00039 //
00040 //M*/
00041 #include "precomp.hpp"
00042 
00043 namespace cv
00044 {
00045 
00046 // inner product
00047 static float innerProduct(Point2f &v1, Point2f &v2)
00048 {
00049     return v1.x * v2.y - v1.y * v2.x;
00050 }
00051 
00052 static void findCircle3pts(Point2f *pts, Point2f &center, float &radius)
00053 {
00054     // two edges of the triangle v1, v2
00055     Point2f v1 = pts[1] - pts[0];
00056     Point2f v2 = pts[2] - pts[0];
00057 
00058     if (innerProduct(v1, v2) == 0.0f)
00059     {
00060         // v1, v2 colineation, can not determine a unique circle
00061         // find the longtest distance as diameter line
00062         float d1 = (float)norm(pts[0] - pts[1]);
00063         float d2 = (float)norm(pts[0] - pts[2]);
00064         float d3 = (float)norm(pts[1] - pts[2]);
00065         if (d1 >= d2 && d1 >= d3)
00066         {
00067             center = (pts[0] + pts[1]) / 2.0f;
00068             radius = (d1 / 2.0f);
00069         }
00070         else if (d2 >= d1 && d2 >= d3)
00071         {
00072             center = (pts[0] + pts[2]) / 2.0f;
00073             radius = (d2 / 2.0f);
00074         }
00075         else if (d3 >= d1 && d3 >= d2)
00076         {
00077             center = (pts[1] + pts[2]) / 2.0f;
00078             radius = (d3 / 2.0f);
00079         }
00080     }
00081     else
00082     {
00083         // center is intersection of midperpendicular lines of the two edges v1, v2
00084         // a1*x + b1*y = c1 where a1 = v1.x, b1 = v1.y
00085         // a2*x + b2*y = c2 where a2 = v2.x, b2 = v2.y
00086         Point2f midPoint1 = (pts[0] + pts[1]) / 2.0f;
00087         float c1 = midPoint1.x * v1.x + midPoint1.y * v1.y;
00088         Point2f midPoint2 = (pts[0] + pts[2]) / 2.0f;
00089         float c2 = midPoint2.x * v2.x + midPoint2.y * v2.y;
00090         float det = v1.x * v2.y - v1.y * v2.x;
00091         float cx = (c1 * v2.y - c2 * v1.y) / det;
00092         float cy = (v1.x * c2 - v2.x * c1) / det;
00093         center.x = (float)cx;
00094         center.y = (float)cy;
00095         cx -= pts[0].x;
00096         cy -= pts[0].y;
00097         radius = (float)(std::sqrt(cx *cx + cy * cy));
00098     }
00099 }
00100 
00101 const float EPS = 1.0e-4f;
00102 
00103 static void findEnclosingCircle3pts_orLess_32f(Point2f *pts, int count, Point2f &center, float &radius)
00104 {
00105     switch (count)
00106     {
00107     case 1:
00108         center = pts[0];
00109         radius = 0.0f;
00110         break;
00111     case 2:
00112         center.x = (pts[0].x + pts[1].x) / 2.0f;
00113         center.y = (pts[0].y + pts[1].y) / 2.0f;
00114         radius = (float)(norm(pts[0] - pts[1]) / 2.0);
00115         break;
00116     case 3:
00117         findCircle3pts(pts, center, radius);
00118         break;
00119     default:
00120         break;
00121     }
00122 
00123     radius += EPS;
00124 }
00125 
00126 template<typename PT>
00127 static void findThirdPoint(const PT *pts, int i, int j, Point2f &center, float &radius)
00128 {
00129     center.x = (float)(pts[j].x + pts[i].x) / 2.0f;
00130     center.y = (float)(pts[j].y + pts[i].y) / 2.0f;
00131     float dx = (float)(pts[j].x - pts[i].x);
00132     float dy = (float)(pts[j].y - pts[i].y);
00133     radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;
00134 
00135     for (int k = 0; k < j; ++k)
00136     {
00137         dx = center.x - (float)pts[k].x;
00138         dy = center.y - (float)pts[k].y;
00139         if (norm(Point2f(dx, dy)) < radius)
00140         {
00141             continue;
00142         }
00143         else
00144         {
00145             Point2f ptsf[3];
00146             ptsf[0] = (Point2f)pts[i];
00147             ptsf[1] = (Point2f)pts[j];
00148             ptsf[2] = (Point2f)pts[k];
00149             findEnclosingCircle3pts_orLess_32f(ptsf, 3, center, radius);
00150         }
00151     }
00152 }
00153 
00154 
00155 template<typename PT>
00156 void findSecondPoint(const PT *pts, int i, Point2f &center, float &radius)
00157 {
00158     center.x = (float)(pts[0].x + pts[i].x) / 2.0f;
00159     center.y = (float)(pts[0].y + pts[i].y) / 2.0f;
00160     float dx = (float)(pts[0].x - pts[i].x);
00161     float dy = (float)(pts[0].y - pts[i].y);
00162     radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;
00163 
00164     for (int j = 1; j < i; ++j)
00165     {
00166         dx = center.x - (float)pts[j].x;
00167         dy = center.y - (float)pts[j].y;
00168         if (norm(Point2f(dx, dy)) < radius)
00169         {
00170             continue;
00171         }
00172         else
00173         {
00174             findThirdPoint(pts, i, j, center, radius);
00175         }
00176     }
00177 }
00178 
00179 
00180 template<typename PT>
00181 static void findMinEnclosingCircle(const PT *pts, int count, Point2f &center, float &radius)
00182 {
00183     center.x = (float)(pts[0].x + pts[1].x) / 2.0f;
00184     center.y = (float)(pts[0].y + pts[1].y) / 2.0f;
00185     float dx = (float)(pts[0].x - pts[1].x);
00186     float dy = (float)(pts[0].y - pts[1].y);
00187     radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;
00188 
00189     for (int i = 2; i < count; ++i)
00190     {
00191         dx = (float)pts[i].x - center.x;
00192         dy = (float)pts[i].y - center.y;
00193         float d = (float)norm(Point2f(dx, dy));
00194         if (d < radius)
00195         {
00196             continue;
00197         }
00198         else
00199         {
00200             findSecondPoint(pts, i, center, radius);
00201         }
00202     }
00203 }
00204 } // namespace cv
00205 
00206 // see Welzl, Emo. Smallest enclosing disks (balls and ellipsoids). Springer Berlin Heidelberg, 1991.
00207 void cv::minEnclosingCircle( InputArray _points, Point2f& _center, float& _radius )
00208 {
00209     Mat points = _points.getMat();
00210     int count = points.checkVector(2);
00211     int depth = points.depth();
00212     Point2f center;
00213     float radius = 0.f;
00214     CV_Assert(count >= 0 && (depth == CV_32F || depth == CV_32S));
00215 
00216     _center.x = _center.y = 0.f;
00217     _radius = 0.f;
00218 
00219     if( count == 0 )
00220         return;
00221 
00222     bool is_float = depth == CV_32F;
00223     const Point* ptsi = points.ptr<Point>();
00224     const Point2f* ptsf = points.ptr<Point2f>();
00225 
00226     // point count <= 3
00227     if (count <= 3)
00228     {
00229         Point2f ptsf3[3];
00230         for (int i = 0; i < count; ++i)
00231         {
00232             ptsf3[i] = (is_float) ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
00233         }
00234         findEnclosingCircle3pts_orLess_32f(ptsf3, count, center, radius);
00235         _center = center;
00236         _radius = radius;
00237         return;
00238     }
00239 
00240     if (is_float)
00241     {
00242         findMinEnclosingCircle<Point2f>(ptsf, count, center, radius);
00243 #if 0
00244         for (size_t m = 0; m < count; ++m)
00245         {
00246             float d = (float)norm(ptsf[m] - center);
00247             if (d > radius)
00248             {
00249                 printf("error!\n");
00250             }
00251         }
00252 #endif
00253     }
00254     else
00255     {
00256         findMinEnclosingCircle<Point>(ptsi, count, center, radius);
00257 #if 0
00258         for (size_t m = 0; m < count; ++m)
00259         {
00260             double dx = ptsi[m].x - center.x;
00261             double dy = ptsi[m].y - center.y;
00262             double d = std::sqrt(dx * dx + dy * dy);
00263             if (d > radius)
00264             {
00265                 printf("error!\n");
00266             }
00267         }
00268 #endif
00269     }
00270     _center = center;
00271     _radius = radius;
00272 }
00273 
00274 
00275 // calculates length of a curve (e.g. contour perimeter)
00276 double cv::arcLength( InputArray _curve, bool is_closed )
00277 {
00278     Mat curve = _curve.getMat();
00279     int count = curve.checkVector(2);
00280     int depth = curve.depth();
00281     CV_Assert( count >= 0 && (depth == CV_32F || depth == CV_32S));
00282     double perimeter = 0;
00283 
00284     int i;
00285 
00286     if( count <= 1 )
00287         return 0.;
00288 
00289     bool is_float = depth == CV_32F;
00290     int last = is_closed ? count-1 : 0;
00291     const Point* pti = curve.ptr<Point>();
00292     const Point2f * ptf = curve.ptr<Point2f >();
00293 
00294     Point2f  prev = is_float ? ptf[last] : Point2f ((float)pti[last].x,(float)pti[last].y);
00295 
00296     for( i = 0; i < count; i++ )
00297     {
00298         Point2f  p = is_float ? ptf[i] : Point2f ((float)pti[i].x,(float)pti[i].y);
00299         float dx = p.x - prev.x, dy = p.y - prev.y;
00300         perimeter += std::sqrt(dx*dx + dy*dy);
00301 
00302         prev = p;
00303     }
00304 
00305     return perimeter;
00306 }
00307 
00308 // area of a whole sequence
00309 double cv::contourArea( InputArray _contour, bool oriented )
00310 {
00311     Mat contour = _contour.getMat();
00312     int npoints = contour.checkVector(2);
00313     int depth = contour.depth();
00314     CV_Assert(npoints >= 0 && (depth == CV_32F || depth == CV_32S));
00315 
00316     if( npoints == 0 )
00317         return 0.;
00318 
00319     double a00 = 0;
00320     bool is_float = depth == CV_32F;
00321     const Point* ptsi = contour.ptr<Point>();
00322     const Point2f * ptsf = contour.ptr<Point2f >();
00323     Point2f  prev = is_float ? ptsf[npoints-1] : Point2f ((float)ptsi[npoints-1].x, (float)ptsi[npoints-1].y);
00324 
00325     for( int i = 0; i < npoints; i++ )
00326     {
00327         Point2f  p = is_float ? ptsf[i] : Point2f ((float)ptsi[i].x, (float)ptsi[i].y);
00328         a00 += (double)prev.x * p.y - (double)prev.y * p.x;
00329         prev = p;
00330     }
00331 
00332     a00 *= 0.5;
00333     if( !oriented )
00334         a00 = fabs(a00);
00335 
00336     return a00;
00337 }
00338 
00339 
00340 cv::RotatedRect cv::fitEllipse( InputArray _points )
00341 {
00342     Mat points = _points.getMat();
00343     int i, n = points.checkVector(2);
00344     int depth = points.depth();
00345     CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S));
00346 
00347     RotatedRect box;
00348 
00349     if( n < 5 )
00350         CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" );
00351 
00352     // New fitellipse algorithm, contributed by Dr. Daniel Weiss
00353     Point2f  c(0,0);
00354     double gfp[5], rp[5], t;
00355     const double min_eps = 1e-8;
00356     bool is_float = depth == CV_32F;
00357     const Point* ptsi = points.ptr<Point>();
00358     const Point2f * ptsf = points.ptr<Point2f >();
00359 
00360     AutoBuffer<double> _Ad(n*5), _bd(n);
00361     double *Ad = _Ad, *bd = _bd;
00362 
00363     // first fit for parameters A - E
00364     Mat A( n, 5, CV_64F, Ad );
00365     Mat b( n, 1, CV_64F, bd );
00366     Mat x( 5, 1, CV_64F, gfp );
00367 
00368     for( i = 0; i < n; i++ )
00369     {
00370         Point2f  p = is_float ? ptsf[i] : Point2f ((float)ptsi[i].x, (float)ptsi[i].y);
00371         c += p;
00372     }
00373     c.x /= n;
00374     c.y /= n;
00375 
00376     for( i = 0; i < n; i++ )
00377     {
00378         Point2f  p = is_float ? ptsf[i] : Point2f ((float)ptsi[i].x, (float)ptsi[i].y);
00379         p -= c;
00380 
00381         bd[i] = 10000.0; // 1.0?
00382         Ad[i*5] = -(double)p.x * p.x; // A - C signs inverted as proposed by APP
00383         Ad[i*5 + 1] = -(double)p.y * p.y;
00384         Ad[i*5 + 2] = -(double)p.x * p.y;
00385         Ad[i*5 + 3] = p.x;
00386         Ad[i*5 + 4] = p.y;
00387     }
00388 
00389     solve(A, b, x, DECOMP_SVD);
00390 
00391     // now use general-form parameters A - E to find the ellipse center:
00392     // differentiate general form wrt x/y to get two equations for cx and cy
00393     A = Mat( 2, 2, CV_64F, Ad );
00394     b = Mat( 2, 1, CV_64F, bd );
00395     x = Mat( 2, 1, CV_64F, rp );
00396     Ad[0] = 2 * gfp[0];
00397     Ad[1] = Ad[2] = gfp[2];
00398     Ad[3] = 2 * gfp[1];
00399     bd[0] = gfp[3];
00400     bd[1] = gfp[4];
00401     solve( A, b, x, DECOMP_SVD );
00402 
00403     // re-fit for parameters A - C with those center coordinates
00404     A = Mat( n, 3, CV_64F, Ad );
00405     b = Mat( n, 1, CV_64F, bd );
00406     x = Mat( 3, 1, CV_64F, gfp );
00407     for( i = 0; i < n; i++ )
00408     {
00409         Point2f  p = is_float ? ptsf[i] : Point2f ((float)ptsi[i].x, (float)ptsi[i].y);
00410         p -= c;
00411         bd[i] = 1.0;
00412         Ad[i * 3] = (p.x - rp[0]) * (p.x - rp[0]);
00413         Ad[i * 3 + 1] = (p.y - rp[1]) * (p.y - rp[1]);
00414         Ad[i * 3 + 2] = (p.x - rp[0]) * (p.y - rp[1]);
00415     }
00416     solve(A, b, x, DECOMP_SVD);
00417 
00418     // store angle and radii
00419     rp[4] = -0.5 * atan2(gfp[2], gfp[1] - gfp[0]); // convert from APP angle usage
00420     if( fabs(gfp[2]) > min_eps )
00421         t = gfp[2]/sin(-2.0 * rp[4]);
00422     else // ellipse is rotated by an integer multiple of pi/2
00423         t = gfp[1] - gfp[0];
00424     rp[2] = fabs(gfp[0] + gfp[1] - t);
00425     if( rp[2] > min_eps )
00426         rp[2] = std::sqrt(2.0 / rp[2]);
00427     rp[3] = fabs(gfp[0] + gfp[1] + t);
00428     if( rp[3] > min_eps )
00429         rp[3] = std::sqrt(2.0 / rp[3]);
00430 
00431     box.center.x = (float)rp[0] + c.x;
00432     box.center.y = (float)rp[1] + c.y;
00433     box.size.width = (float)(rp[2]*2);
00434     box.size.height = (float)(rp[3]*2);
00435     if( box.size.width > box.size.height )
00436     {
00437         float tmp;
00438         CV_SWAP( box.size.width, box.size.height, tmp );
00439         box.angle = (float)(90 + rp[4]*180/CV_PI);
00440     }
00441     if( box.angle < -180 )
00442         box.angle += 360;
00443     if( box.angle > 360 )
00444         box.angle -= 360;
00445 
00446     return box;
00447 }
00448 
00449 
00450 namespace cv
00451 {
00452 
00453 // Calculates bounding rectagnle of a point set or retrieves already calculated
00454 static Rect pointSetBoundingRect( const Mat& points )
00455 {
00456     int npoints = points.checkVector(2);
00457     int depth = points.depth();
00458     CV_Assert(npoints >= 0 && (depth == CV_32F || depth == CV_32S));
00459 
00460     int  xmin = 0, ymin = 0, xmax = -1, ymax = -1, i;
00461     bool is_float = depth == CV_32F;
00462 
00463     if( npoints == 0 )
00464         return Rect();
00465 
00466     const Point* pts = points.ptr<Point>();
00467     Point pt = pts[0];
00468 
00469 #if CV_SSE4_2
00470     if(cv::checkHardwareSupport(CV_CPU_SSE4_2))
00471     {
00472         if( !is_float )
00473         {
00474             __m128i minval, maxval;
00475             minval = maxval = _mm_loadl_epi64((const __m128i*)(&pt)); //min[0]=pt.x, min[1]=pt.y
00476 
00477             for( i = 1; i < npoints; i++ )
00478             {
00479                 __m128i ptXY = _mm_loadl_epi64((const __m128i*)&pts[i]);
00480                 minval = _mm_min_epi32(ptXY, minval);
00481                 maxval = _mm_max_epi32(ptXY, maxval);
00482             }
00483             xmin = _mm_cvtsi128_si32(minval);
00484             ymin = _mm_cvtsi128_si32(_mm_srli_si128(minval, 4));
00485             xmax = _mm_cvtsi128_si32(maxval);
00486             ymax = _mm_cvtsi128_si32(_mm_srli_si128(maxval, 4));
00487         }
00488         else
00489         {
00490             __m128 minvalf, maxvalf, z = _mm_setzero_ps(), ptXY = _mm_setzero_ps();
00491             minvalf = maxvalf = _mm_loadl_pi(z, (const __m64*)(&pt));
00492 
00493             for( i = 1; i < npoints; i++ )
00494             {
00495                 ptXY = _mm_loadl_pi(ptXY, (const __m64*)&pts[i]);
00496 
00497                 minvalf = _mm_min_ps(minvalf, ptXY);
00498                 maxvalf = _mm_max_ps(maxvalf, ptXY);
00499             }
00500 
00501             float xyminf[2], xymaxf[2];
00502             _mm_storel_pi((__m64*)xyminf, minvalf);
00503             _mm_storel_pi((__m64*)xymaxf, maxvalf);
00504             xmin = cvFloor(xyminf[0]);
00505             ymin = cvFloor(xyminf[1]);
00506             xmax = cvFloor(xymaxf[0]);
00507             ymax = cvFloor(xymaxf[1]);
00508         }
00509     }
00510     else
00511 #endif
00512     {
00513         if( !is_float )
00514         {
00515             xmin = xmax = pt.x;
00516             ymin = ymax = pt.y;
00517 
00518             for( i = 1; i < npoints; i++ )
00519             {
00520                 pt = pts[i];
00521 
00522                 if( xmin > pt.x )
00523                     xmin = pt.x;
00524 
00525                 if( xmax < pt.x )
00526                     xmax = pt.x;
00527 
00528                 if( ymin > pt.y )
00529                     ymin = pt.y;
00530 
00531                 if( ymax < pt.y )
00532                     ymax = pt.y;
00533             }
00534         }
00535         else
00536         {
00537             Cv32suf v;
00538             // init values
00539             xmin = xmax = CV_TOGGLE_FLT(pt.x);
00540             ymin = ymax = CV_TOGGLE_FLT(pt.y);
00541 
00542             for( i = 1; i < npoints; i++ )
00543             {
00544                 pt = pts[i];
00545                 pt.x = CV_TOGGLE_FLT(pt.x);
00546                 pt.y = CV_TOGGLE_FLT(pt.y);
00547 
00548                 if( xmin > pt.x )
00549                     xmin = pt.x;
00550 
00551                 if( xmax < pt.x )
00552                     xmax = pt.x;
00553 
00554                 if( ymin > pt.y )
00555                     ymin = pt.y;
00556 
00557                 if( ymax < pt.y )
00558                     ymax = pt.y;
00559             }
00560 
00561             v.i = CV_TOGGLE_FLT(xmin); xmin = cvFloor(v.f);
00562             v.i = CV_TOGGLE_FLT(ymin); ymin = cvFloor(v.f);
00563             // because right and bottom sides of the bounding rectangle are not inclusive
00564             // (note +1 in width and height calculation below), cvFloor is used here instead of cvCeil
00565             v.i = CV_TOGGLE_FLT(xmax); xmax = cvFloor(v.f);
00566             v.i = CV_TOGGLE_FLT(ymax); ymax = cvFloor(v.f);
00567         }
00568     }
00569 
00570     return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1);
00571 }
00572 
00573 
00574 static Rect maskBoundingRect( const Mat& img )
00575 {
00576     CV_Assert( img.depth() <= CV_8S && img.channels() == 1 );
00577 
00578     Size size = img.size();
00579     int xmin = size.width, ymin = -1, xmax = -1, ymax = -1, i, j, k;
00580 
00581     for( i = 0; i < size.height; i++ )
00582     {
00583         const uchar* _ptr = img.ptr(i);
00584         const uchar* ptr = (const uchar*)alignPtr(_ptr, 4);
00585         int have_nz = 0, k_min, offset = (int)(ptr - _ptr);
00586         j = 0;
00587         offset = MIN(offset, size.width);
00588         for( ; j < offset; j++ )
00589             if( _ptr[j] )
00590             {
00591                 have_nz = 1;
00592                 break;
00593             }
00594         if( j < offset )
00595         {
00596             if( j < xmin )
00597                 xmin = j;
00598             if( j > xmax )
00599                 xmax = j;
00600         }
00601         if( offset < size.width )
00602         {
00603             xmin -= offset;
00604             xmax -= offset;
00605             size.width -= offset;
00606             j = 0;
00607             for( ; j <= xmin - 4; j += 4 )
00608                 if( *((int*)(ptr+j)) )
00609                     break;
00610             for( ; j < xmin; j++ )
00611                 if( ptr[j] )
00612                 {
00613                     xmin = j;
00614                     if( j > xmax )
00615                         xmax = j;
00616                     have_nz = 1;
00617                     break;
00618                 }
00619             k_min = MAX(j-1, xmax);
00620             k = size.width - 1;
00621             for( ; k > k_min && (k&3) != 3; k-- )
00622                 if( ptr[k] )
00623                     break;
00624             if( k > k_min && (k&3) == 3 )
00625             {
00626                 for( ; k > k_min+3; k -= 4 )
00627                     if( *((int*)(ptr+k-3)) )
00628                         break;
00629             }
00630             for( ; k > k_min; k-- )
00631                 if( ptr[k] )
00632                 {
00633                     xmax = k;
00634                     have_nz = 1;
00635                     break;
00636                 }
00637             if( !have_nz )
00638             {
00639                 j &= ~3;
00640                 for( ; j <= k - 3; j += 4 )
00641                     if( *((int*)(ptr+j)) )
00642                         break;
00643                 for( ; j <= k; j++ )
00644                     if( ptr[j] )
00645                     {
00646                         have_nz = 1;
00647                         break;
00648                     }
00649             }
00650             xmin += offset;
00651             xmax += offset;
00652             size.width += offset;
00653         }
00654         if( have_nz )
00655         {
00656             if( ymin < 0 )
00657                 ymin = i;
00658             ymax = i;
00659         }
00660     }
00661 
00662     if( xmin >= size.width )
00663         xmin = ymin = 0;
00664     return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1);
00665 }
00666 
00667 }
00668 
00669 cv::Rect cv::boundingRect(InputArray array)
00670 {
00671     Mat m = array.getMat();
00672     return m.depth() <= CV_8U ? maskBoundingRect(m) : pointSetBoundingRect(m);
00673 }
00674 
00675 ////////////////////////////////////////////// C API ///////////////////////////////////////////
00676 
00677 CV_IMPL int
00678 cvMinEnclosingCircle( const void* array, CvPoint2D32f * _center, float *_radius )
00679 {
00680     cv::AutoBuffer<double> abuf;
00681     cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf);
00682     cv::Point2f  center;
00683     float radius;
00684 
00685     cv::minEnclosingCircle(points, center, radius);
00686     if(_center)
00687         *_center = center;
00688     if(_radius)
00689         *_radius = radius;
00690     return 1;
00691 }
00692 
00693 static void
00694 icvMemCopy( double **buf1, double **buf2, double **buf3, int *b_max )
00695 {
00696     CV_Assert( (*buf1 != NULL || *buf2 != NULL) && *buf3 != NULL );
00697 
00698     int bb = *b_max;
00699     if( *buf2 == NULL )
00700     {
00701         *b_max = 2 * (*b_max);
00702         *buf2 = (double *)cvAlloc( (*b_max) * sizeof( double ));
00703 
00704         memcpy( *buf2, *buf3, bb * sizeof( double ));
00705 
00706         *buf3 = *buf2;
00707         cvFree( buf1 );
00708         *buf1 = NULL;
00709     }
00710     else
00711     {
00712         *b_max = 2 * (*b_max);
00713         *buf1 = (double *) cvAlloc( (*b_max) * sizeof( double ));
00714 
00715         memcpy( *buf1, *buf3, bb * sizeof( double ));
00716 
00717         *buf3 = *buf1;
00718         cvFree( buf2 );
00719         *buf2 = NULL;
00720     }
00721 }
00722 
00723 
00724 /* area of a contour sector */
00725 static double icvContourSecArea( CvSeq * contour, CvSlice slice )
00726 {
00727     CvPoint pt;                 /*  pointer to points   */
00728     CvPoint pt_s, pt_e;         /*  first and last points  */
00729     CvSeqReader reader;         /*  points reader of contour   */
00730 
00731     int p_max = 2, p_ind;
00732     int lpt, flag, i;
00733     double a00;                 /* unnormalized moments m00    */
00734     double xi, yi, xi_1, yi_1, x0, y0, dxy, sk, sk1, t;
00735     double x_s, y_s, nx, ny, dx, dy, du, dv;
00736     double eps = 1.e-5;
00737     double *p_are1, *p_are2, *p_are;
00738     double area = 0;
00739 
00740     CV_Assert( contour != NULL && CV_IS_SEQ_POINT_SET( contour ));
00741 
00742     lpt = cvSliceLength( slice, contour );
00743     /*if( n2 >= n1 )
00744         lpt = n2 - n1 + 1;
00745     else
00746         lpt = contour->total - n1 + n2 + 1;*/
00747 
00748     if( contour->total <= 0 || lpt <= 2 )
00749         return 0.;
00750 
00751     a00 = x0 = y0 = xi_1 = yi_1 = 0;
00752     sk1 = 0;
00753     flag = 0;
00754     dxy = 0;
00755     p_are1 = (double *) cvAlloc( p_max * sizeof( double ));
00756 
00757     p_are = p_are1;
00758     p_are2 = NULL;
00759 
00760     cvStartReadSeq( contour, &reader, 0 );
00761     cvSetSeqReaderPos( &reader, slice.start_index );
00762     CV_READ_SEQ_ELEM( pt_s, reader );
00763     p_ind = 0;
00764     cvSetSeqReaderPos( &reader, slice.end_index );
00765     CV_READ_SEQ_ELEM( pt_e, reader );
00766 
00767 /*    normal coefficients    */
00768     nx = pt_s.y - pt_e.y;
00769     ny = pt_e.x - pt_s.x;
00770     cvSetSeqReaderPos( &reader, slice.start_index );
00771 
00772     while( lpt-- > 0 )
00773     {
00774         CV_READ_SEQ_ELEM( pt, reader );
00775 
00776         if( flag == 0 )
00777         {
00778             xi_1 = (double) pt.x;
00779             yi_1 = (double) pt.y;
00780             x0 = xi_1;
00781             y0 = yi_1;
00782             sk1 = 0;
00783             flag = 1;
00784         }
00785         else
00786         {
00787             xi = (double) pt.x;
00788             yi = (double) pt.y;
00789 
00790 /****************   edges intersection examination   **************************/
00791             sk = nx * (xi - pt_s.x) + ny * (yi - pt_s.y);
00792             if( (fabs( sk ) < eps && lpt > 0) || sk * sk1 < -eps )
00793             {
00794                 if( fabs( sk ) < eps )
00795                 {
00796                     dxy = xi_1 * yi - xi * yi_1;
00797                     a00 = a00 + dxy;
00798                     dxy = xi * y0 - x0 * yi;
00799                     a00 = a00 + dxy;
00800 
00801                     if( p_ind >= p_max )
00802                         icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
00803 
00804                     p_are[p_ind] = a00 / 2.;
00805                     p_ind++;
00806                     a00 = 0;
00807                     sk1 = 0;
00808                     x0 = xi;
00809                     y0 = yi;
00810                     dxy = 0;
00811                 }
00812                 else
00813                 {
00814 /*  define intersection point    */
00815                     dv = yi - yi_1;
00816                     du = xi - xi_1;
00817                     dx = ny;
00818                     dy = -nx;
00819                     if( fabs( du ) > eps )
00820                         t = ((yi_1 - pt_s.y) * du + dv * (pt_s.x - xi_1)) /
00821                             (du * dy - dx * dv);
00822                     else
00823                         t = (xi_1 - pt_s.x) / dx;
00824                     if( t > eps && t < 1 - eps )
00825                     {
00826                         x_s = pt_s.x + t * dx;
00827                         y_s = pt_s.y + t * dy;
00828                         dxy = xi_1 * y_s - x_s * yi_1;
00829                         a00 += dxy;
00830                         dxy = x_s * y0 - x0 * y_s;
00831                         a00 += dxy;
00832                         if( p_ind >= p_max )
00833                             icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
00834 
00835                         p_are[p_ind] = a00 / 2.;
00836                         p_ind++;
00837 
00838                         a00 = 0;
00839                         sk1 = 0;
00840                         x0 = x_s;
00841                         y0 = y_s;
00842                         dxy = x_s * yi - xi * y_s;
00843                     }
00844                 }
00845             }
00846             else
00847                 dxy = xi_1 * yi - xi * yi_1;
00848 
00849             a00 += dxy;
00850             xi_1 = xi;
00851             yi_1 = yi;
00852             sk1 = sk;
00853 
00854         }
00855     }
00856 
00857     xi = x0;
00858     yi = y0;
00859     dxy = xi_1 * yi - xi * yi_1;
00860 
00861     a00 += dxy;
00862 
00863     if( p_ind >= p_max )
00864         icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
00865 
00866     p_are[p_ind] = a00 / 2.;
00867     p_ind++;
00868 
00869     // common area calculation
00870     area = 0;
00871     for( i = 0; i < p_ind; i++ )
00872         area += fabs( p_are[i] );
00873 
00874     if( p_are1 != NULL )
00875         cvFree( &p_are1 );
00876     else if( p_are2 != NULL )
00877         cvFree( &p_are2 );
00878 
00879     return area;
00880 }
00881 
00882 
00883 /* external contour area function */
00884 CV_IMPL double
00885 cvContourArea( const void *array, CvSlice slice, int oriented )
00886 {
00887     double area = 0;
00888 
00889     CvContour contour_header;
00890     CvSeq* contour = 0;
00891     CvSeqBlock block;
00892 
00893     if( CV_IS_SEQ( array ))
00894     {
00895         contour = (CvSeq*)array;
00896         if( !CV_IS_SEQ_POLYLINE( contour ))
00897             CV_Error( CV_StsBadArg, "Unsupported sequence type" );
00898     }
00899     else
00900     {
00901         contour = cvPointSeqFromMat( CV_SEQ_KIND_CURVE, array, &contour_header, &block );
00902     }
00903 
00904     if( cvSliceLength( slice, contour ) == contour->total )
00905     {
00906         cv::AutoBuffer<double> abuf;
00907         cv::Mat points = cv::cvarrToMat(contour, false, false, 0, &abuf);
00908         return cv::contourArea( points, oriented !=0 );
00909     }
00910 
00911     if( CV_SEQ_ELTYPE( contour ) != CV_32SC2 )
00912         CV_Error( CV_StsUnsupportedFormat,
00913         "Only curves with integer coordinates are supported in case of contour slice" );
00914     area = icvContourSecArea( contour, slice );
00915     return oriented ? area : fabs(area);
00916 }
00917 
00918 
00919 /* calculates length of a curve (e.g. contour perimeter) */
00920 CV_IMPL  double
00921 cvArcLength( const void *array, CvSlice slice, int is_closed )
00922 {
00923     double perimeter = 0;
00924 
00925     int i, j = 0, count;
00926     const int N = 16;
00927     float buf[N];
00928     CvMat buffer = cvMat( 1, N, CV_32F, buf );
00929     CvSeqReader reader;
00930     CvContour contour_header;
00931     CvSeq* contour = 0;
00932     CvSeqBlock block;
00933 
00934     if( CV_IS_SEQ( array ))
00935     {
00936         contour = (CvSeq*)array;
00937         if( !CV_IS_SEQ_POLYLINE( contour ))
00938             CV_Error( CV_StsBadArg, "Unsupported sequence type" );
00939         if( is_closed < 0 )
00940             is_closed = CV_IS_SEQ_CLOSED( contour );
00941     }
00942     else
00943     {
00944         is_closed = is_closed > 0;
00945         contour = cvPointSeqFromMat(
00946                                     CV_SEQ_KIND_CURVE | (is_closed ? CV_SEQ_FLAG_CLOSED : 0),
00947                                     array, &contour_header, &block );
00948     }
00949 
00950     if( contour->total > 1 )
00951     {
00952         int is_float = CV_SEQ_ELTYPE( contour ) == CV_32FC2;
00953 
00954         cvStartReadSeq( contour, &reader, 0 );
00955         cvSetSeqReaderPos( &reader, slice.start_index );
00956         count = cvSliceLength( slice, contour );
00957 
00958         count -= !is_closed && count == contour->total;
00959 
00960         // scroll the reader by 1 point
00961         reader.prev_elem = reader.ptr;
00962         CV_NEXT_SEQ_ELEM( sizeof(CvPoint), reader );
00963 
00964         for( i = 0; i < count; i++ )
00965         {
00966             float dx, dy;
00967 
00968             if( !is_float )
00969             {
00970                 CvPoint* pt = (CvPoint*)reader.ptr;
00971                 CvPoint* prev_pt = (CvPoint*)reader.prev_elem;
00972 
00973                 dx = (float)pt->x - (float)prev_pt->x;
00974                 dy = (float)pt->y - (float)prev_pt->y;
00975             }
00976             else
00977             {
00978                 CvPoint2D32f* pt = (CvPoint2D32f*)reader.ptr;
00979                 CvPoint2D32f* prev_pt = (CvPoint2D32f*)reader.prev_elem;
00980 
00981                 dx = pt->x - prev_pt->x;
00982                 dy = pt->y - prev_pt->y;
00983             }
00984 
00985             reader.prev_elem = reader.ptr;
00986             CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
00987             // Bugfix by Axel at rubico.com 2010-03-22, affects closed slices only
00988             // wraparound not handled by CV_NEXT_SEQ_ELEM
00989             if( is_closed && i == count - 2 )
00990                 cvSetSeqReaderPos( &reader, slice.start_index );
00991 
00992             buffer.data.fl[j] = dx * dx + dy * dy;
00993             if( ++j == N || i == count - 1 )
00994             {
00995                 buffer.cols = j;
00996                 cvPow( &buffer, &buffer, 0.5 );
00997                 for( ; j > 0; j-- )
00998                     perimeter += buffer.data.fl[j-1];
00999             }
01000         }
01001     }
01002 
01003     return perimeter;
01004 }
01005 
01006 
01007 CV_IMPL CvBox2D 
01008 cvFitEllipse2( const CvArr* array )
01009 {
01010     cv::AutoBuffer<double> abuf;
01011     cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf);
01012     return cv::fitEllipse(points);
01013 }
01014 
01015 /* Calculates bounding rectagnle of a point set or retrieves already calculated */
01016 CV_IMPL  CvRect 
01017 cvBoundingRect( CvArr* array, int update )
01018 {
01019     CvRect   rect;
01020     CvContour contour_header;
01021     CvSeq* ptseq = 0;
01022     CvSeqBlock block;
01023 
01024     CvMat stub, *mat = 0;
01025     int calculate = update;
01026 
01027     if( CV_IS_SEQ( array ))
01028     {
01029         ptseq = (CvSeq*)array;
01030         if( !CV_IS_SEQ_POINT_SET( ptseq ))
01031             CV_Error( CV_StsBadArg, "Unsupported sequence type" );
01032 
01033         if( ptseq->header_size < (int)sizeof(CvContour))
01034         {
01035             update = 0;
01036             calculate = 1;
01037         }
01038     }
01039     else
01040     {
01041         mat = cvGetMat( array, &stub );
01042         if( CV_MAT_TYPE(mat->type) == CV_32SC2 ||
01043             CV_MAT_TYPE(mat->type) == CV_32FC2 )
01044         {
01045             ptseq = cvPointSeqFromMat(CV_SEQ_KIND_GENERIC, mat, &contour_header, &block);
01046             mat = 0;
01047         }
01048         else if( CV_MAT_TYPE(mat->type) != CV_8UC1 &&
01049                 CV_MAT_TYPE(mat->type) != CV_8SC1 )
01050             CV_Error( CV_StsUnsupportedFormat,
01051                 "The image/matrix format is not supported by the function" );
01052         update = 0;
01053         calculate = 1;
01054     }
01055 
01056     if( !calculate )
01057         return ((CvContour*)ptseq)->rect;
01058 
01059     if( mat )
01060     {
01061         rect = cv::maskBoundingRect(cv::cvarrToMat(mat));
01062     }
01063     else if( ptseq->total )
01064     {
01065         cv::AutoBuffer<double> abuf;
01066         rect = cv::pointSetBoundingRect(cv::cvarrToMat(ptseq, false, false, 0, &abuf));
01067     }
01068     if( update )
01069         ((CvContour*)ptseq)->rect = rect;
01070     return rect;
01071 }
01072 
01073 
01074 /* End of file. */
01075