modules/imgproc/src/shapedescr.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 //
  10 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41 #include "precomp.hpp"
  42 namespace cv
  43 {
  44
  45 // inner product
  46 static float innerProduct(Point2f &v1, Point2f &v2)
  47 {
  48     return v1.x * v2.y - v1.y * v2.x;
  49 }
  50
  51 static void findCircle3pts(Point2f *pts, Point2f &center, float &radius)
  52 {
  53     // two edges of the triangle v1, v2
  54     Point2f v1 = pts[1] - pts[0];
  55     Point2f v2 = pts[2] - pts[0];
  56
  57     if (innerProduct(v1, v2) == 0.0f)
  58     {
  59         // v1, v2 colineation, can not determine a unique circle
  60         // find the longtest distance as diameter line
  61         float d1 = (float)norm(pts[0] - pts[1]);
  62         float d2 = (float)norm(pts[0] - pts[2]);
  63         float d3 = (float)norm(pts[1] - pts[2]);
  64         if (d1 >= d2 && d1 >= d3)
  65         {
  66             center = (pts[0] + pts[1]) / 2.0f;
  67             radius = (d1 / 2.0f);
  68         }
  69         else if (d2 >= d1 && d2 >= d3)
  70         {
  71             center = (pts[0] + pts[2]) / 2.0f;
  72             radius = (d2 / 2.0f);
  73         }
  74         else if (d3 >= d1 && d3 >= d2)
  75         {
  76             center = (pts[1] + pts[2]) / 2.0f;
  77             radius = (d3 / 2.0f);
  78         }
  79     }
  80     else
  81     {
  82         // center is intersection of midperpendicular lines of the two edges v1, v2
  83         // a1*x + b1*y = c1 where a1 = v1.x, b1 = v1.y
  84         // a2*x + b2*y = c2 where a2 = v2.x, b2 = v2.y
  85         Point2f midPoint1 = (pts[0] + pts[1]) / 2.0f;
  86         float c1 = midPoint1.x * v1.x + midPoint1.y * v1.y;
  87         Point2f midPoint2 = (pts[0] + pts[2]) / 2.0f;
  88         float c2 = midPoint2.x * v2.x + midPoint2.y * v2.y;
  89         float det = v1.x * v2.y - v1.y * v2.x;
  90         float cx = (c1 * v2.y - c2 * v1.y) / det;
  91         float cy = (v1.x * c2 - v2.x * c1) / det;
  92         center.x = (float)cx;
  93         center.y = (float)cy;
  94         cx -= pts[0].x;
  95         cy -= pts[0].y;
  96         radius = (float)(std::sqrt(cx *cx + cy * cy));
  97     }
  98 }
  99
 100 const float EPS = 1.0e-4f;
 101
 102 static void findEnclosingCircle3pts_orLess_32f(Point2f *pts, int count, Point2f &center, float &radius)
 103 {
 104     switch (count)
 105     {
 106     case 1:
 107         center = pts[0];
 108         radius = 0.0f;
 109         break;
 110     case 2:
 111         center.x = (pts[0].x + pts[1].x) / 2.0f;
 112         center.y = (pts[0].y + pts[1].y) / 2.0f;
 113         radius = (float)(norm(pts[0] - pts[1]) / 2.0);
 114         break;
 115     case 3:
 116         findCircle3pts(pts, center, radius);
 117         break;
 118     default:
 119         break;
 120     }
 121
 122     radius += EPS;
 123 }
 124
 125 template<typename PT>
 126 static void findThirdPoint(const PT *pts, int i, int j, Point2f &center, float &radius)
 127 {
 128     center.x = (float)(pts[j].x + pts[i].x) / 2.0f;
 129     center.y = (float)(pts[j].y + pts[i].y) / 2.0f;
 130     float dx = (float)(pts[j].x - pts[i].x);
 131     float dy = (float)(pts[j].y - pts[i].y);
 132     radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;
 133
 134     for (int k = 0; k < j; ++k)
 135     {
 136         dx = center.x - (float)pts[k].x;
 137         dy = center.y - (float)pts[k].y;
 138         if (norm(Point2f(dx, dy)) < radius)
 139         {
 140             continue;
 141         }
 142         else
 143         {
 144             Point2f ptsf[3];
 145             ptsf[0] = (Point2f)pts[i];
 146             ptsf[1] = (Point2f)pts[j];
 147             ptsf[2] = (Point2f)pts[k];
 148             findEnclosingCircle3pts_orLess_32f(ptsf, 3, center, radius);
 149         }
 150     }
 151 }
 152
 153
 154 template<typename PT>
 155 void findSecondPoint(const PT *pts, int i, Point2f &center, float &radius)
 156 {
 157     center.x = (float)(pts[0].x + pts[i].x) / 2.0f;
 158     center.y = (float)(pts[0].y + pts[i].y) / 2.0f;
 159     float dx = (float)(pts[0].x - pts[i].x);
 160     float dy = (float)(pts[0].y - pts[i].y);
 161     radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;
 162
 163     for (int j = 1; j < i; ++j)
 164     {
 165         dx = center.x - (float)pts[j].x;
 166         dy = center.y - (float)pts[j].y;
 167         if (norm(Point2f(dx, dy)) < radius)
 168         {
 169             continue;
 170         }
 171         else
 172         {
 173             findThirdPoint(pts, i, j, center, radius);
 174         }
 175     }
 176 }
 177
 178
 179 template<typename PT>
 180 static void findMinEnclosingCircle(const PT *pts, int count, Point2f &center, float &radius)
 181 {
 182     center.x = (float)(pts[0].x + pts[1].x) / 2.0f;
 183     center.y = (float)(pts[0].y + pts[1].y) / 2.0f;
 184     float dx = (float)(pts[0].x - pts[1].x);
 185     float dy = (float)(pts[0].y - pts[1].y);
 186     radius = (float)norm(Point2f(dx, dy)) / 2.0f + EPS;
 187
 188     for (int i = 2; i < count; ++i)
 189     {
 190         dx = (float)pts[i].x - center.x;
 191         dy = (float)pts[i].y - center.y;
 192         float d = (float)norm(Point2f(dx, dy));
 193         if (d < radius)
 194         {
 195             continue;
 196         }
 197         else
 198         {
 199             findSecondPoint(pts, i, center, radius);
 200         }
 201     }
 202 }
 203 } // namespace cv
 204
 205 // see Welzl, Emo. Smallest enclosing disks (balls and ellipsoids). Springer Berlin Heidelberg, 1991.
 206 void cv::minEnclosingCircle( InputArray _points, Point2f& _center, float& _radius )
 207 {
 208     CV_INSTRUMENT_REGION()
 209
 210     Mat points = _points.getMat();
 211     int count = points.checkVector(2);
 212     int depth = points.depth();
 213     Point2f center;
 214     float radius = 0.f;
 215     CV_Assert(count >= 0 && (depth == CV_32F || depth == CV_32S));
 216
 217     _center.x = _center.y = 0.f;
 218     _radius = 0.f;
 219
 220     if( count == 0 )
 221         return;
 222
 223     bool is_float = depth == CV_32F;
 224     const Point* ptsi = points.ptr<Point>();
 225     const Point2f* ptsf = points.ptr<Point2f>();
 226
 227     // point count <= 3
 228     if (count <= 3)
 229     {
 230         Point2f ptsf3[3];
 231         for (int i = 0; i < count; ++i)
 232         {
 233             ptsf3[i] = (is_float) ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
 234         }
 235         findEnclosingCircle3pts_orLess_32f(ptsf3, count, center, radius);
 236         _center = center;
 237         _radius = radius;
 238         return;
 239     }
 240
 241     if (is_float)
 242     {
 243         findMinEnclosingCircle<Point2f>(ptsf, count, center, radius);
 244 #if 0
 245         for (size_t m = 0; m < count; ++m)
 246         {
 247             float d = (float)norm(ptsf[m] - center);
 248             if (d > radius)
 249             {
 250                 printf("error!\n");
 251             }
 252         }
 253 #endif
 254     }
 255     else
 256     {
 257         findMinEnclosingCircle<Point>(ptsi, count, center, radius);
 258 #if 0
 259         for (size_t m = 0; m < count; ++m)
 260         {
 261             double dx = ptsi[m].x - center.x;
 262             double dy = ptsi[m].y - center.y;
 263             double d = std::sqrt(dx * dx + dy * dy);
 264             if (d > radius)
 265             {
 266                 printf("error!\n");
 267             }
 268         }
 269 #endif
 270     }
 271     _center = center;
 272     _radius = radius;
 273 }
 274
 275
 276 // calculates length of a curve (e.g. contour perimeter)
 277 double cv::arcLength( InputArray _curve, bool is_closed )
 278 {
 279     CV_INSTRUMENT_REGION()
 280
 281     Mat curve = _curve.getMat();
 282     int count = curve.checkVector(2);
 283     int depth = curve.depth();
 284     CV_Assert( count >= 0 && (depth == CV_32F || depth == CV_32S));
 285     double perimeter = 0;
 286
 287     int i;
 288
 289     if( count <= 1 )
 290         return 0.;
 291
 292     bool is_float = depth == CV_32F;
 293     int last = is_closed ? count-1 : 0;
 294     const Point* pti = curve.ptr<Point>();
 295     const Point2f* ptf = curve.ptr<Point2f>();
 296
 297     Point2f prev = is_float ? ptf[last] : Point2f((float)pti[last].x,(float)pti[last].y);
 298
 299     for( i = 0; i < count; i++ )
 300     {
 301         Point2f p = is_float ? ptf[i] : Point2f((float)pti[i].x,(float)pti[i].y);
 302         float dx = p.x - prev.x, dy = p.y - prev.y;
 303         perimeter += std::sqrt(dx*dx + dy*dy);
 304
 305         prev = p;
 306     }
 307
 308     return perimeter;
 309 }
 310
 311 // area of a whole sequence
 312 double cv::contourArea( InputArray _contour, bool oriented )
 313 {
 314     CV_INSTRUMENT_REGION()
 315
 316     Mat contour = _contour.getMat();
 317     int npoints = contour.checkVector(2);
 318     int depth = contour.depth();
 319     CV_Assert(npoints >= 0 && (depth == CV_32F || depth == CV_32S));
 320
 321     if( npoints == 0 )
 322         return 0.;
 323
 324     double a00 = 0;
 325     bool is_float = depth == CV_32F;
 326     const Point* ptsi = contour.ptr<Point>();
 327     const Point2f* ptsf = contour.ptr<Point2f>();
 328     Point2f prev = is_float ? ptsf[npoints-1] : Point2f((float)ptsi[npoints-1].x, (float)ptsi[npoints-1].y);
 329
 330     for( int i = 0; i < npoints; i++ )
 331     {
 332         Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
 333         a00 += (double)prev.x * p.y - (double)prev.y * p.x;
 334         prev = p;
 335     }
 336
 337     a00 *= 0.5;
 338     if( !oriented )
 339         a00 = fabs(a00);
 340
 341     return a00;
 342 }
 343
 344
 345 cv::RotatedRect cv::fitEllipse( InputArray _points )
 346 {
 347     CV_INSTRUMENT_REGION()
 348
 349     Mat points = _points.getMat();
 350     int i, n = points.checkVector(2);
 351     int depth = points.depth();
 352     CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S));
 353
 354     RotatedRect box;
 355
 356     if( n < 5 )
 357         CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" );
 358
 359     // New fitellipse algorithm, contributed by Dr. Daniel Weiss
 360     Point2f c(0,0);
 361     double gfp[5] = {0}, rp[5] = {0}, t;
 362     const double min_eps = 1e-8;
 363     bool is_float = depth == CV_32F;
 364     const Point* ptsi = points.ptr<Point>();
 365     const Point2f* ptsf = points.ptr<Point2f>();
 366
 367     AutoBuffer<double> _Ad(n*5), _bd(n);
 368     double *Ad = _Ad, *bd = _bd;
 369
 370     // first fit for parameters A - E
 371     Mat A( n, 5, CV_64F, Ad );
 372     Mat b( n, 1, CV_64F, bd );
 373     Mat x( 5, 1, CV_64F, gfp );
 374
 375     for( i = 0; i < n; i++ )
 376     {
 377         Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
 378         c += p;
 379     }
 380     c.x /= n;
 381     c.y /= n;
 382
 383     for( i = 0; i < n; i++ )
 384     {
 385         Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
 386         p -= c;
 387
 388         bd[i] = 10000.0; // 1.0?
 389         Ad[i*5] = -(double)p.x * p.x; // A - C signs inverted as proposed by APP
 390         Ad[i*5 + 1] = -(double)p.y * p.y;
 391         Ad[i*5 + 2] = -(double)p.x * p.y;
 392         Ad[i*5 + 3] = p.x;
 393         Ad[i*5 + 4] = p.y;
 394     }
 395
 396     solve(A, b, x, DECOMP_SVD);
 397
 398     // now use general-form parameters A - E to find the ellipse center:
 399     // differentiate general form wrt x/y to get two equations for cx and cy
 400     A = Mat( 2, 2, CV_64F, Ad );
 401     b = Mat( 2, 1, CV_64F, bd );
 402     x = Mat( 2, 1, CV_64F, rp );
 403     Ad[0] = 2 * gfp[0];
 404     Ad[1] = Ad[2] = gfp[2];
 405     Ad[3] = 2 * gfp[1];
 406     bd[0] = gfp[3];
 407     bd[1] = gfp[4];
 408     solve( A, b, x, DECOMP_SVD );
 409
 410     // re-fit for parameters A - C with those center coordinates
 411     A = Mat( n, 3, CV_64F, Ad );
 412     b = Mat( n, 1, CV_64F, bd );
 413     x = Mat( 3, 1, CV_64F, gfp );
 414     for( i = 0; i < n; i++ )
 415     {
 416         Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
 417         p -= c;
 418         bd[i] = 1.0;
 419         Ad[i * 3] = (p.x - rp[0]) * (p.x - rp[0]);
 420         Ad[i * 3 + 1] = (p.y - rp[1]) * (p.y - rp[1]);
 421         Ad[i * 3 + 2] = (p.x - rp[0]) * (p.y - rp[1]);
 422     }
 423     solve(A, b, x, DECOMP_SVD);
 424
 425     // store angle and radii
 426     rp[4] = -0.5 * atan2(gfp[2], gfp[1] - gfp[0]); // convert from APP angle usage
 427     if( fabs(gfp[2]) > min_eps )
 428         t = gfp[2]/sin(-2.0 * rp[4]);
 429     else // ellipse is rotated by an integer multiple of pi/2
 430         t = gfp[1] - gfp[0];
 431     rp[2] = fabs(gfp[0] + gfp[1] - t);
 432     if( rp[2] > min_eps )
 433         rp[2] = std::sqrt(2.0 / rp[2]);
 434     rp[3] = fabs(gfp[0] + gfp[1] + t);
 435     if( rp[3] > min_eps )
 436         rp[3] = std::sqrt(2.0 / rp[3]);
 437
 438     box.center.x = (float)rp[0] + c.x;
 439     box.center.y = (float)rp[1] + c.y;
 440     box.size.width = (float)(rp[2]*2);
 441     box.size.height = (float)(rp[3]*2);
 442     if( box.size.width > box.size.height )
 443     {
 444         float tmp;
 445         CV_SWAP( box.size.width, box.size.height, tmp );
 446         box.angle = (float)(90 + rp[4]*180/CV_PI);
 447     }
 448     if( box.angle < -180 )
 449         box.angle += 360;
 450     if( box.angle > 360 )
 451         box.angle -= 360;
 452
 453     return box;
 454 }
 455
 456 cv::RotatedRect cv::fitEllipseAMS( InputArray _points )
 457 {
 458     Mat points = _points.getMat();
 459     int i, n = points.checkVector(2);
 460     int depth = points.depth();
 461     CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S));
 462
 463     RotatedRect box;
 464
 465     if( n < 5 )
 466         CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" );
 467
 468     Point2f c(0,0);
 469
 470     bool is_float = depth == CV_32F;
 471     const Point* ptsi = points.ptr<Point>();
 472     const Point2f* ptsf = points.ptr<Point2f>();
 473
 474     Mat A( n, 6, CV_64F);
 475     Matx<double, 6, 6> DM;
 476     Matx<double, 5, 5> M;
 477     Matx<double, 5, 1> pVec;
 478     Matx<double, 6, 1> coeffs;
 479
 480     double x0, y0, a, b, theta;
 481
 482     for( i = 0; i < n; i++ )
 483     {
 484         Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
 485         c += p;
 486     }
 487     c.x /= (float)n;
 488     c.y /= (float)n;
 489
 490     for( i = 0; i < n; i++ )
 491     {
 492         Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
 493         p -= c;
 494
 495         A.at<double>(i,0) = (double)(p.x)*(p.x);
 496         A.at<double>(i,1) = (double)(p.x)*(p.y);
 497         A.at<double>(i,2) = (double)(p.y)*(p.y);
 498         A.at<double>(i,3) = (double)p.x;
 499         A.at<double>(i,4) = (double)p.y;
 500         A.at<double>(i,5) = 1.0;
 501     }
 502     cv::mulTransposed( A, DM, true, noArray(), 1.0, -1 );
 503     DM *= (1.0/n);
 504     double dnm = ( DM(2,5)*(DM(0,5) + DM(2,5)) - (DM(1,5)*DM(1,5)) );
 505     double ddm =  (4.*(DM(0,5) + DM(2,5))*( (DM(0,5)*DM(2,5)) - (DM(1,5)*DM(1,5))));
 506     double ddmm = (2.*(DM(0,5) + DM(2,5))*( (DM(0,5)*DM(2,5)) - (DM(1,5)*DM(1,5))));
 507
 508     M(0,0)=((-DM(0,0) + DM(0,2) + DM(0,5)*DM(0,5))*(DM(1,5)*DM(1,5)) + (-2*DM(0,1)*DM(1,5) + DM(0,5)*(DM(0,0) \
 509             - (DM(0,5)*DM(0,5)) + (DM(1,5)*DM(1,5))))*DM(2,5) + (DM(0,0) - (DM(0,5)*DM(0,5)))*(DM(2,5)*DM(2,5))) / ddm;
 510     M(0,1)=((DM(1,5)*DM(1,5))*(-DM(0,1) + DM(1,2) + DM(0,5)*DM(1,5)) + (DM(0,1)*DM(0,5) - ((DM(0,5)*DM(0,5)) + 2*DM(1,1))*DM(1,5) + \
 511             (DM(1,5)*DM(1,5)*DM(1,5)))*DM(2,5) + (DM(0,1) - DM(0,5)*DM(1,5))*(DM(2,5)*DM(2,5))) / ddm;
 512     M(0,2)=(-2*DM(1,2)*DM(1,5)*DM(2,5) - DM(0,5)*(DM(2,5)*DM(2,5))*(DM(0,5) + DM(2,5)) + DM(0,2)*dnm + \
 513             (DM(1,5)*DM(1,5))*(DM(2,2) + DM(2,5)*(DM(0,5) + DM(2,5))))/ddm;
 514     M(0,3)=(DM(1,5)*(DM(1,5)*DM(2,3) - 2*DM(1,3)*DM(2,5)) + DM(0,3)*dnm) / ddm;
 515     M(0,4)=(DM(1,5)*(DM(1,5)*DM(2,4) - 2*DM(1,4)*DM(2,5)) + DM(0,4)*dnm) / ddm;
 516     M(1,0)=(-(DM(0,2)*DM(0,5)*DM(1,5)) + (2*DM(0,1)*DM(0,5) - DM(0,0)*DM(1,5))*DM(2,5))/ddmm;
 517     M(1,1)=(-(DM(0,1)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,2)*DM(1,5)) + 2*DM(1,1)*DM(2,5)))/ddmm;
 518     M(1,2)=(-(DM(0,2)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,5)*DM(2,2)) + 2*DM(1,2)*DM(2,5)))/ddmm;
 519     M(1,3)=(-(DM(0,3)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,5)*DM(2,3)) + 2*DM(1,3)*DM(2,5)))/ddmm;
 520     M(1,4)=(-(DM(0,4)*DM(1,5)*DM(2,5)) + DM(0,5)*(-(DM(1,5)*DM(2,4)) + 2*DM(1,4)*DM(2,5)))/ddmm;
 521     M(2,0)=(-2*DM(0,1)*DM(0,5)*DM(1,5) + (DM(0,0) + (DM(0,5)*DM(0,5)))*(DM(1,5)*DM(1,5)) + DM(0,5)*(-(DM(0,5)*DM(0,5)) \
 522             + (DM(1,5)*DM(1,5)))*DM(2,5) - (DM(0,5)*DM(0,5))*(DM(2,5)*DM(2,5)) + DM(0,2)*(-(DM(1,5)*DM(1,5)) + DM(0,5)*(DM(0,5) + DM(2,5)))) / ddm;
 523     M(2,1)=((DM(0,5)*DM(0,5))*(DM(1,2) - DM(1,5)*DM(2,5)) + (DM(1,5)*DM(1,5))*(DM(0,1) - DM(1,2) + DM(1,5)*DM(2,5)) \
 524             + DM(0,5)*(DM(1,2)*DM(2,5) + DM(1,5)*(-2*DM(1,1) + (DM(1,5)*DM(1,5)) - (DM(2,5)*DM(2,5))))) / ddm;
 525     M(2,2)=((DM(0,5)*DM(0,5))*(DM(2,2) - (DM(2,5)*DM(2,5))) + (DM(1,5)*DM(1,5))*(DM(0,2) - DM(2,2) + (DM(2,5)*DM(2,5))) + \
 526              DM(0,5)*(-2*DM(1,2)*DM(1,5) + DM(2,5)*((DM(1,5)*DM(1,5)) + DM(2,2) - (DM(2,5)*DM(2,5))))) / ddm;
 527     M(2,3)=((DM(1,5)*DM(1,5))*(DM(0,3) - DM(2,3)) + (DM(0,5)*DM(0,5))*DM(2,3) + DM(0,5)*(-2*DM(1,3)*DM(1,5) + DM(2,3)*DM(2,5))) / ddm;
 528     M(2,4)=((DM(1,5)*DM(1,5))*(DM(0,4) - DM(2,4)) + (DM(0,5)*DM(0,5))*DM(2,4) + DM(0,5)*(-2*DM(1,4)*DM(1,5) + DM(2,4)*DM(2,5))) / ddm;
 529     M(3,0)=DM(0,3);
 530     M(3,1)=DM(1,3);
 531     M(3,2)=DM(2,3);
 532     M(3,3)=DM(3,3);
 533     M(3,4)=DM(3,4);
 534     M(4,0)=DM(0,4);
 535     M(4,1)=DM(1,4);
 536     M(4,2)=DM(2,4);
 537     M(4,3)=DM(3,4);
 538     M(4,4)=DM(4,4);
 539
 540     if (fabs(cv::determinant(M)) > 1.0e-10) {
 541             Mat eVal, eVec;
 542             eigenNonSymmetric(M, eVal, eVec);
 543
 544             // Select the eigen vector {a,b,c,d,e} which has the lowest eigenvalue
 545             int minpos = 0;
 546             double normi, normEVali, normMinpos, normEValMinpos;
 547             normMinpos = sqrt(eVec.at<double>(minpos,0)*eVec.at<double>(minpos,0) + eVec.at<double>(minpos,1)*eVec.at<double>(minpos,1) + \
 548                               eVec.at<double>(minpos,2)*eVec.at<double>(minpos,2) + eVec.at<double>(minpos,3)*eVec.at<double>(minpos,3) + \
 549                               eVec.at<double>(minpos,4)*eVec.at<double>(minpos,4) );
 550             normEValMinpos = eVal.at<double>(minpos,0) * normMinpos;
 551             for (i=1; i<5; i++) {
 552                 normi = sqrt(eVec.at<double>(i,0)*eVec.at<double>(i,0) + eVec.at<double>(i,1)*eVec.at<double>(i,1) + \
 553                              eVec.at<double>(i,2)*eVec.at<double>(i,2) + eVec.at<double>(i,3)*eVec.at<double>(i,3) + \
 554                              eVec.at<double>(i,4)*eVec.at<double>(i,4) );
 555                 normEVali = eVal.at<double>(i,0) * normi;
 556                 if (normEVali < normEValMinpos) {
 557                     minpos = i;
 558                     normMinpos=normi;
 559                     normEValMinpos=normEVali;
 560                 }
 561             };
 562
 563             pVec(0) =eVec.at<double>(minpos,0) / normMinpos;
 564             pVec(1) =eVec.at<double>(minpos,1) / normMinpos;
 565             pVec(2) =eVec.at<double>(minpos,2) / normMinpos;
 566             pVec(3) =eVec.at<double>(minpos,3) / normMinpos;
 567             pVec(4) =eVec.at<double>(minpos,4) / normMinpos;
 568
 569             coeffs(0) =pVec(0) ;
 570             coeffs(1) =pVec(1) ;
 571             coeffs(2) =pVec(2) ;
 572             coeffs(3) =pVec(3) ;
 573             coeffs(4) =pVec(4) ;
 574             coeffs(5) =-pVec(0) *DM(0,5)-pVec(1) *DM(1,5)-coeffs(2) *DM(2,5);
 575
 576         // Check that an elliptical solution has been found. AMS sometimes produces Parabolic solutions.
 577         bool is_ellipse = (coeffs(0)  < 0 && \
 578                            coeffs(2)  < (coeffs(1) *coeffs(1) )/(4.*coeffs(0) ) && \
 579                            coeffs(5)  > (-(coeffs(2) *(coeffs(3) *coeffs(3) )) + coeffs(1) *coeffs(3) *coeffs(4)  - coeffs(0) *(coeffs(4) *coeffs(4) )) / \
 580                                         ((coeffs(1) *coeffs(1) ) - 4*coeffs(0) *coeffs(2) )) || \
 581                           (coeffs(0)  > 0 && \
 582                            coeffs(2)  > (coeffs(1) *coeffs(1) )/(4.*coeffs(0) ) && \
 583                            coeffs(5)  < (-(coeffs(2) *(coeffs(3) *coeffs(3) )) + coeffs(1) *coeffs(3) *coeffs(4)  - coeffs(0) *(coeffs(4) *coeffs(4) )) / \
 584                                         ( (coeffs(1) *coeffs(1) ) - 4*coeffs(0) *coeffs(2) ));
 585         if (is_ellipse) {
 586             double u1 = pVec(2) *pVec(3) *pVec(3)  - pVec(1) *pVec(3) *pVec(4)  + pVec(0) *pVec(4) *pVec(4)  + pVec(1) *pVec(1) *coeffs(5) ;
 587             double u2 = pVec(0) *pVec(2) *coeffs(5) ;
 588             double l1 = sqrt(pVec(1) *pVec(1)  + (pVec(0)  - pVec(2) )*(pVec(0)  - pVec(2) ));
 589             double l2 = pVec(0)  + pVec(2) ;
 590             double l3 = pVec(1) *pVec(1)  - 4.0*pVec(0) *pVec(2) ;
 591             double p1 = 2.0*pVec(2) *pVec(3)  - pVec(1) *pVec(4) ;
 592             double p2 = 2.0*pVec(0) *pVec(4) -(pVec(1) *pVec(3) );
 593
 594             x0 = p1/l3 + c.x;
 595             y0 = p2/l3 + c.y;
 596             a = std::sqrt(2.)*sqrt((u1 - 4.0*u2)/((l1 - l2)*l3));
 597             b = std::sqrt(2.)*sqrt(-1.0*((u1 - 4.0*u2)/((l1 + l2)*l3)));
 598             if (pVec(1)  == 0) {
 599                 if (pVec(0)  < pVec(2) ) {
 600                     theta = 0;
 601                 } else {
 602                     theta = CV_PI/2.;
 603                 }
 604             } else {
 605                 theta = CV_PI/2. + 0.5*std::atan2(pVec(1) , (pVec(0)  - pVec(2) ));
 606             }
 607
 608             box.center.x = (float)x0; // +c.x;
 609             box.center.y = (float)y0; // +c.y;
 610             box.size.width = (float)(2.0*a);
 611             box.size.height = (float)(2.0*b);
 612             if( box.size.width > box.size.height )
 613             {
 614                 float tmp;
 615                 CV_SWAP( box.size.width, box.size.height, tmp );
 616                 box.angle = (float)(90 + theta*180/CV_PI);
 617             } else {
 618                 box.angle = (float)(fmod(theta*180/CV_PI,180.0));
 619             };
 620
 621
 622         } else {
 623             box = cv::fitEllipseDirect( points );
 624         }
 625     } else {
 626         box = cv::fitEllipse( points );
 627     }
 628
 629     return box;
 630 }
 631
 632 cv::RotatedRect cv::fitEllipseDirect( InputArray _points )
 633 {
 634     Mat points = _points.getMat();
 635     int i, n = points.checkVector(2);
 636     int depth = points.depth();
 637     CV_Assert( n >= 0 && (depth == CV_32F || depth == CV_32S));
 638
 639     RotatedRect box;
 640
 641     if( n < 5 )
 642         CV_Error( CV_StsBadSize, "There should be at least 5 points to fit the ellipse" );
 643
 644     Point2f c(0,0);
 645
 646     bool is_float = (depth == CV_32F);
 647     const Point*   ptsi = points.ptr<Point>();
 648     const Point2f* ptsf = points.ptr<Point2f>();
 649
 650     Mat A( n, 6, CV_64F);
 651     Matx<double, 6, 6> DM;
 652     Matx33d M, TM, Q;
 653     Matx<double, 3, 1> pVec;
 654
 655     double x0, y0, a, b, theta, Ts;
 656
 657     for( i = 0; i < n; i++ )
 658     {
 659         Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
 660         c += p;
 661     }
 662     c.x /= (float)n;
 663     c.y /= (float)n;
 664
 665     for( i = 0; i < n; i++ )
 666     {
 667         Point2f p = is_float ? ptsf[i] : Point2f((float)ptsi[i].x, (float)ptsi[i].y);
 668         p -= c;
 669
 670         A.at<double>(i,0) = (double)(p.x)*(p.x);
 671         A.at<double>(i,1) = (double)(p.x)*(p.y);
 672         A.at<double>(i,2) = (double)(p.y)*(p.y);
 673         A.at<double>(i,3) = (double)p.x;
 674         A.at<double>(i,4) = (double)p.y;
 675         A.at<double>(i,5) = 1.0;
 676     }
 677     cv::mulTransposed( A, DM, true, noArray(), 1.0, -1 );
 678     DM *= (1.0/n);
 679
 680     TM(0,0) = DM(0,5)*DM(3,5)*DM(4,4) - DM(0,5)*DM(3,4)*DM(4,5) - DM(0,4)*DM(3,5)*DM(5,4) + \
 681               DM(0,3)*DM(4,5)*DM(5,4) + DM(0,4)*DM(3,4)*DM(5,5) - DM(0,3)*DM(4,4)*DM(5,5);
 682     TM(0,1) = DM(1,5)*DM(3,5)*DM(4,4) - DM(1,5)*DM(3,4)*DM(4,5) - DM(1,4)*DM(3,5)*DM(5,4) + \
 683               DM(1,3)*DM(4,5)*DM(5,4) + DM(1,4)*DM(3,4)*DM(5,5) - DM(1,3)*DM(4,4)*DM(5,5);
 684     TM(0,2) = DM(2,5)*DM(3,5)*DM(4,4) - DM(2,5)*DM(3,4)*DM(4,5) - DM(2,4)*DM(3,5)*DM(5,4) + \
 685               DM(2,3)*DM(4,5)*DM(5,4) + DM(2,4)*DM(3,4)*DM(5,5) - DM(2,3)*DM(4,4)*DM(5,5);
 686     TM(1,0) = DM(0,5)*DM(3,3)*DM(4,5) - DM(0,5)*DM(3,5)*DM(4,3) + DM(0,4)*DM(3,5)*DM(5,3) - \
 687               DM(0,3)*DM(4,5)*DM(5,3) - DM(0,4)*DM(3,3)*DM(5,5) + DM(0,3)*DM(4,3)*DM(5,5);
 688     TM(1,1) = DM(1,5)*DM(3,3)*DM(4,5) - DM(1,5)*DM(3,5)*DM(4,3) + DM(1,4)*DM(3,5)*DM(5,3) - \
 689               DM(1,3)*DM(4,5)*DM(5,3) - DM(1,4)*DM(3,3)*DM(5,5) + DM(1,3)*DM(4,3)*DM(5,5);
 690     TM(1,2) = DM(2,5)*DM(3,3)*DM(4,5) - DM(2,5)*DM(3,5)*DM(4,3) + DM(2,4)*DM(3,5)*DM(5,3) - \
 691               DM(2,3)*DM(4,5)*DM(5,3) - DM(2,4)*DM(3,3)*DM(5,5) + DM(2,3)*DM(4,3)*DM(5,5);
 692     TM(2,0) = DM(0,5)*DM(3,4)*DM(4,3) - DM(0,5)*DM(3,3)*DM(4,4) - DM(0,4)*DM(3,4)*DM(5,3) + \
 693               DM(0,3)*DM(4,4)*DM(5,3) + DM(0,4)*DM(3,3)*DM(5,4) - DM(0,3)*DM(4,3)*DM(5,4);
 694     TM(2,1) = DM(1,5)*DM(3,4)*DM(4,3) - DM(1,5)*DM(3,3)*DM(4,4) - DM(1,4)*DM(3,4)*DM(5,3) + \
 695               DM(1,3)*DM(4,4)*DM(5,3) + DM(1,4)*DM(3,3)*DM(5,4) - DM(1,3)*DM(4,3)*DM(5,4);
 696     TM(2,2) = DM(2,5)*DM(3,4)*DM(4,3) - DM(2,5)*DM(3,3)*DM(4,4) - DM(2,4)*DM(3,4)*DM(5,3) + \
 697               DM(2,3)*DM(4,4)*DM(5,3) + DM(2,4)*DM(3,3)*DM(5,4) - DM(2,3)*DM(4,3)*DM(5,4);
 698
 699     Ts=(-(DM(3,5)*DM(4,4)*DM(5,3)) + DM(3,4)*DM(4,5)*DM(5,3) + DM(3,5)*DM(4,3)*DM(5,4) - \
 700           DM(3,3)*DM(4,5)*DM(5,4)  - DM(3,4)*DM(4,3)*DM(5,5) + DM(3,3)*DM(4,4)*DM(5,5));
 701
 702     M(0,0) = (DM(2,0) + (DM(2,3)*TM(0,0) + DM(2,4)*TM(1,0) + DM(2,5)*TM(2,0))/Ts)/2.;
 703     M(0,1) = (DM(2,1) + (DM(2,3)*TM(0,1) + DM(2,4)*TM(1,1) + DM(2,5)*TM(2,1))/Ts)/2.;
 704     M(0,2) = (DM(2,2) + (DM(2,3)*TM(0,2) + DM(2,4)*TM(1,2) + DM(2,5)*TM(2,2))/Ts)/2.;
 705     M(1,0) = -DM(1,0) - (DM(1,3)*TM(0,0) + DM(1,4)*TM(1,0) + DM(1,5)*TM(2,0))/Ts;
 706     M(1,1) = -DM(1,1) - (DM(1,3)*TM(0,1) + DM(1,4)*TM(1,1) + DM(1,5)*TM(2,1))/Ts;
 707     M(1,2) = -DM(1,2) - (DM(1,3)*TM(0,2) + DM(1,4)*TM(1,2) + DM(1,5)*TM(2,2))/Ts;
 708     M(2,0) = (DM(0,0) + (DM(0,3)*TM(0,0) + DM(0,4)*TM(1,0) + DM(0,5)*TM(2,0))/Ts)/2.;
 709     M(2,1) = (DM(0,1) + (DM(0,3)*TM(0,1) + DM(0,4)*TM(1,1) + DM(0,5)*TM(2,1))/Ts)/2.;
 710     M(2,2) = (DM(0,2) + (DM(0,3)*TM(0,2) + DM(0,4)*TM(1,2) + DM(0,5)*TM(2,2))/Ts)/2.;
 711
 712     if (fabs(cv::determinant(M)) > 1.0e-10) {
 713         Mat eVal, eVec;
 714         eigenNonSymmetric(M, eVal, eVec);
 715
 716         // Select the eigen vector {a,b,c} which satisfies 4ac-b^2 > 0
 717         double cond[3];
 718         cond[0]=(4.0 * eVec.at<double>(0,0) * eVec.at<double>(0,2) - eVec.at<double>(0,1) * eVec.at<double>(0,1));
 719         cond[1]=(4.0 * eVec.at<double>(1,0) * eVec.at<double>(1,2) - eVec.at<double>(1,1) * eVec.at<double>(1,1));
 720         cond[2]=(4.0 * eVec.at<double>(2,0) * eVec.at<double>(2,2) - eVec.at<double>(2,1) * eVec.at<double>(2,1));
 721         if (cond[0]<cond[1]) {
 722             i = (cond[1]<cond[2]) ? 2 : 1;
 723         } else {
 724             i = (cond[0]<cond[2]) ? 2 : 0;
 725         }
 726         double norm = std::sqrt(eVec.at<double>(i,0)*eVec.at<double>(i,0) + eVec.at<double>(i,1)*eVec.at<double>(i,1) + eVec.at<double>(i,2)*eVec.at<double>(i,2));
 727         if (((eVec.at<double>(i,0)<0.0  ? -1 : 1) * (eVec.at<double>(i,1)<0.0  ? -1 : 1) * (eVec.at<double>(i,2)<0.0  ? -1 : 1)) <= 0.0) {
 728                 norm=-1.0*norm;
 729             }
 730         pVec(0) =eVec.at<double>(i,0)/norm; pVec(1) =eVec.at<double>(i,1)/norm;pVec(2) =eVec.at<double>(i,2)/norm;
 731
 732     //  Q = (TM . pVec)/Ts;
 733         Q(0,0) = (TM(0,0)*pVec(0) +TM(0,1)*pVec(1) +TM(0,2)*pVec(2) )/Ts;
 734         Q(0,1) = (TM(1,0)*pVec(0) +TM(1,1)*pVec(1) +TM(1,2)*pVec(2) )/Ts;
 735         Q(0,2) = (TM(2,0)*pVec(0) +TM(2,1)*pVec(1) +TM(2,2)*pVec(2) )/Ts;
 736
 737     // We compute the ellipse properties in the shifted coordinates as doing so improves the numerical accuracy.
 738
 739         double u1 = pVec(2)*Q(0,0)*Q(0,0) - pVec(1)*Q(0,0)*Q(0,1) + pVec(0)*Q(0,1)*Q(0,1) + pVec(1)*pVec(1)*Q(0,2);
 740         double u2 = pVec(0)*pVec(2)*Q(0,2);
 741         double l1 = sqrt(pVec(1)*pVec(1) + (pVec(0) - pVec(2))*(pVec(0) - pVec(2)));
 742         double l2 = pVec(0) + pVec(2) ;
 743         double l3 = pVec(1)*pVec(1) - 4*pVec(0)*pVec(2) ;
 744         double p1 = 2*pVec(2)*Q(0,0) - pVec(1)*Q(0,1);
 745         double p2 = 2*pVec(0)*Q(0,1) - pVec(1)*Q(0,0);
 746
 747         x0 = p1/l3 + c.x;
 748         y0 = p2/l3 + c.y;
 749         a = sqrt(2.)*sqrt((u1 - 4.0*u2)/((l1 - l2)*l3));
 750         b = sqrt(2.)*sqrt(-1.0*((u1 - 4.0*u2)/((l1 + l2)*l3)));
 751         if (pVec(1)  == 0) {
 752             if (pVec(0)  < pVec(2) ) {
 753                 theta = 0;
 754             } else {
 755                 theta = CV_PI/2.;
 756             }
 757         } else {
 758                 theta = CV_PI/2. + 0.5*std::atan2(pVec(1) , (pVec(0)  - pVec(2) ));
 759         }
 760
 761         box.center.x = (float)x0;
 762         box.center.y = (float)y0;
 763         box.size.width = (float)(2.0*a);
 764         box.size.height = (float)(2.0*b);
 765         if( box.size.width > box.size.height )
 766         {
 767             float tmp;
 768             CV_SWAP( box.size.width, box.size.height, tmp );
 769             box.angle = (float)(fmod((90 + theta*180/CV_PI),180.0)) ;
 770         } else {
 771             box.angle = (float)(fmod(theta*180/CV_PI,180.0));
 772         };
 773     } else {
 774         box = cv::fitEllipse( points );
 775     }
 776     return box;
 777 }
 778
 779
 780 namespace cv
 781 {
 782
 783 // Calculates bounding rectagnle of a point set or retrieves already calculated
 784 static Rect pointSetBoundingRect( const Mat& points )
 785 {
 786     int npoints = points.checkVector(2);
 787     int depth = points.depth();
 788     CV_Assert(npoints >= 0 && (depth == CV_32F || depth == CV_32S));
 789
 790     int  xmin = 0, ymin = 0, xmax = -1, ymax = -1, i;
 791     bool is_float = depth == CV_32F;
 792
 793     if( npoints == 0 )
 794         return Rect();
 795
 796     const Point* pts = points.ptr<Point>();
 797     Point pt = pts[0];
 798
 799 #if CV_SSE4_2
 800     if(cv::checkHardwareSupport(CV_CPU_SSE4_2))
 801     {
 802         if( !is_float )
 803         {
 804             __m128i minval, maxval;
 805             minval = maxval = _mm_loadl_epi64((const __m128i*)(&pt)); //min[0]=pt.x, min[1]=pt.y
 806
 807             for( i = 1; i < npoints; i++ )
 808             {
 809                 __m128i ptXY = _mm_loadl_epi64((const __m128i*)&pts[i]);
 810                 minval = _mm_min_epi32(ptXY, minval);
 811                 maxval = _mm_max_epi32(ptXY, maxval);
 812             }
 813             xmin = _mm_cvtsi128_si32(minval);
 814             ymin = _mm_cvtsi128_si32(_mm_srli_si128(minval, 4));
 815             xmax = _mm_cvtsi128_si32(maxval);
 816             ymax = _mm_cvtsi128_si32(_mm_srli_si128(maxval, 4));
 817         }
 818         else
 819         {
 820             __m128 minvalf, maxvalf, z = _mm_setzero_ps(), ptXY = _mm_setzero_ps();
 821             minvalf = maxvalf = _mm_loadl_pi(z, (const __m64*)(&pt));
 822
 823             for( i = 1; i < npoints; i++ )
 824             {
 825                 ptXY = _mm_loadl_pi(ptXY, (const __m64*)&pts[i]);
 826
 827                 minvalf = _mm_min_ps(minvalf, ptXY);
 828                 maxvalf = _mm_max_ps(maxvalf, ptXY);
 829             }
 830
 831             float xyminf[2], xymaxf[2];
 832             _mm_storel_pi((__m64*)xyminf, minvalf);
 833             _mm_storel_pi((__m64*)xymaxf, maxvalf);
 834             xmin = cvFloor(xyminf[0]);
 835             ymin = cvFloor(xyminf[1]);
 836             xmax = cvFloor(xymaxf[0]);
 837             ymax = cvFloor(xymaxf[1]);
 838         }
 839     }
 840     else
 841 #endif
 842     {
 843         if( !is_float )
 844         {
 845             xmin = xmax = pt.x;
 846             ymin = ymax = pt.y;
 847
 848             for( i = 1; i < npoints; i++ )
 849             {
 850                 pt = pts[i];
 851
 852                 if( xmin > pt.x )
 853                     xmin = pt.x;
 854
 855                 if( xmax < pt.x )
 856                     xmax = pt.x;
 857
 858                 if( ymin > pt.y )
 859                     ymin = pt.y;
 860
 861                 if( ymax < pt.y )
 862                     ymax = pt.y;
 863             }
 864         }
 865         else
 866         {
 867             Cv32suf v;
 868             // init values
 869             xmin = xmax = CV_TOGGLE_FLT(pt.x);
 870             ymin = ymax = CV_TOGGLE_FLT(pt.y);
 871
 872             for( i = 1; i < npoints; i++ )
 873             {
 874                 pt = pts[i];
 875                 pt.x = CV_TOGGLE_FLT(pt.x);
 876                 pt.y = CV_TOGGLE_FLT(pt.y);
 877
 878                 if( xmin > pt.x )
 879                     xmin = pt.x;
 880
 881                 if( xmax < pt.x )
 882                     xmax = pt.x;
 883
 884                 if( ymin > pt.y )
 885                     ymin = pt.y;
 886
 887                 if( ymax < pt.y )
 888                     ymax = pt.y;
 889             }
 890
 891             v.i = CV_TOGGLE_FLT(xmin); xmin = cvFloor(v.f);
 892             v.i = CV_TOGGLE_FLT(ymin); ymin = cvFloor(v.f);
 893             // because right and bottom sides of the bounding rectangle are not inclusive
 894             // (note +1 in width and height calculation below), cvFloor is used here instead of cvCeil
 895             v.i = CV_TOGGLE_FLT(xmax); xmax = cvFloor(v.f);
 896             v.i = CV_TOGGLE_FLT(ymax); ymax = cvFloor(v.f);
 897         }
 898     }
 899
 900     return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1);
 901 }
 902
 903
 904 static Rect maskBoundingRect( const Mat& img )
 905 {
 906     CV_Assert( img.depth() <= CV_8S && img.channels() == 1 );
 907
 908     Size size = img.size();
 909     int xmin = size.width, ymin = -1, xmax = -1, ymax = -1, i, j, k;
 910
 911     for( i = 0; i < size.height; i++ )
 912     {
 913         const uchar* _ptr = img.ptr(i);
 914         const uchar* ptr = (const uchar*)alignPtr(_ptr, 4);
 915         int have_nz = 0, k_min, offset = (int)(ptr - _ptr);
 916         j = 0;
 917         offset = MIN(offset, size.width);
 918         for( ; j < offset; j++ )
 919             if( _ptr[j] )
 920             {
 921                 have_nz = 1;
 922                 break;
 923             }
 924         if( j < offset )
 925         {
 926             if( j < xmin )
 927                 xmin = j;
 928             if( j > xmax )
 929                 xmax = j;
 930         }
 931         if( offset < size.width )
 932         {
 933             xmin -= offset;
 934             xmax -= offset;
 935             size.width -= offset;
 936             j = 0;
 937             for( ; j <= xmin - 4; j += 4 )
 938                 if( *((int*)(ptr+j)) )
 939                     break;
 940             for( ; j < xmin; j++ )
 941                 if( ptr[j] )
 942                 {
 943                     xmin = j;
 944                     if( j > xmax )
 945                         xmax = j;
 946                     have_nz = 1;
 947                     break;
 948                 }
 949             k_min = MAX(j-1, xmax);
 950             k = size.width - 1;
 951             for( ; k > k_min && (k&3) != 3; k-- )
 952                 if( ptr[k] )
 953                     break;
 954             if( k > k_min && (k&3) == 3 )
 955             {
 956                 for( ; k > k_min+3; k -= 4 )
 957                     if( *((int*)(ptr+k-3)) )
 958                         break;
 959             }
 960             for( ; k > k_min; k-- )
 961                 if( ptr[k] )
 962                 {
 963                     xmax = k;
 964                     have_nz = 1;
 965                     break;
 966                 }
 967             if( !have_nz )
 968             {
 969                 j &= ~3;
 970                 for( ; j <= k - 3; j += 4 )
 971                     if( *((int*)(ptr+j)) )
 972                         break;
 973                 for( ; j <= k; j++ )
 974                     if( ptr[j] )
 975                     {
 976                         have_nz = 1;
 977                         break;
 978                     }
 979             }
 980             xmin += offset;
 981             xmax += offset;
 982             size.width += offset;
 983         }
 984         if( have_nz )
 985         {
 986             if( ymin < 0 )
 987                 ymin = i;
 988             ymax = i;
 989         }
 990     }
 991
 992     if( xmin >= size.width )
 993         xmin = ymin = 0;
 994     return Rect(xmin, ymin, xmax - xmin + 1, ymax - ymin + 1);
 995 }
 996
 997 }
 998
 999 cv::Rect cv::boundingRect(InputArray array)
1000 {
1001     CV_INSTRUMENT_REGION()
1002
1003     Mat m = array.getMat();
1004     return m.depth() <= CV_8U ? maskBoundingRect(m) : pointSetBoundingRect(m);
1005 }
1006
1007 ////////////////////////////////////////////// C API ///////////////////////////////////////////
1008
1009 CV_IMPL int
1010 cvMinEnclosingCircle( const void* array, CvPoint2D32f * _center, float *_radius )
1011 {
1012     cv::AutoBuffer<double> abuf;
1013     cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf);
1014     cv::Point2f center;
1015     float radius;
1016
1017     cv::minEnclosingCircle(points, center, radius);
1018     if(_center)
1019         *_center = center;
1020     if(_radius)
1021         *_radius = radius;
1022     return 1;
1023 }
1024
1025 static void
1026 icvMemCopy( double **buf1, double **buf2, double **buf3, int *b_max )
1027 {
1028     CV_Assert( (*buf1 != NULL || *buf2 != NULL) && *buf3 != NULL );
1029
1030     int bb = *b_max;
1031     if( *buf2 == NULL )
1032     {
1033         *b_max = 2 * (*b_max);
1034         *buf2 = (double *)cvAlloc( (*b_max) * sizeof( double ));
1035
1036         memcpy( *buf2, *buf3, bb * sizeof( double ));
1037
1038         *buf3 = *buf2;
1039         cvFree( buf1 );
1040         *buf1 = NULL;
1041     }
1042     else
1043     {
1044         *b_max = 2 * (*b_max);
1045         *buf1 = (double *) cvAlloc( (*b_max) * sizeof( double ));
1046
1047         memcpy( *buf1, *buf3, bb * sizeof( double ));
1048
1049         *buf3 = *buf1;
1050         cvFree( buf2 );
1051         *buf2 = NULL;
1052     }
1053 }
1054
1055
1056 /* area of a contour sector */
1057 static double icvContourSecArea( CvSeq * contour, CvSlice slice )
1058 {
1059     CvPoint pt;                 /*  pointer to points   */
1060     CvPoint pt_s, pt_e;         /*  first and last points  */
1061     CvSeqReader reader;         /*  points reader of contour   */
1062
1063     int p_max = 2, p_ind;
1064     int lpt, flag, i;
1065     double a00;                 /* unnormalized moments m00    */
1066     double xi, yi, xi_1, yi_1, x0, y0, dxy, sk, sk1, t;
1067     double x_s, y_s, nx, ny, dx, dy, du, dv;
1068     double eps = 1.e-5;
1069     double *p_are1, *p_are2, *p_are;
1070     double area = 0;
1071
1072     CV_Assert( contour != NULL && CV_IS_SEQ_POINT_SET( contour ));
1073
1074     lpt = cvSliceLength( slice, contour );
1075     /*if( n2 >= n1 )
1076         lpt = n2 - n1 + 1;
1077     else
1078         lpt = contour->total - n1 + n2 + 1;*/
1079
1080     if( contour->total <= 0 || lpt <= 2 )
1081         return 0.;
1082
1083     a00 = x0 = y0 = xi_1 = yi_1 = 0;
1084     sk1 = 0;
1085     flag = 0;
1086     dxy = 0;
1087     p_are1 = (double *) cvAlloc( p_max * sizeof( double ));
1088
1089     p_are = p_are1;
1090     p_are2 = NULL;
1091
1092     cvStartReadSeq( contour, &reader, 0 );
1093     cvSetSeqReaderPos( &reader, slice.start_index );
1094     CV_READ_SEQ_ELEM( pt_s, reader );
1095     p_ind = 0;
1096     cvSetSeqReaderPos( &reader, slice.end_index );
1097     CV_READ_SEQ_ELEM( pt_e, reader );
1098
1099 /*    normal coefficients    */
1100     nx = pt_s.y - pt_e.y;
1101     ny = pt_e.x - pt_s.x;
1102     cvSetSeqReaderPos( &reader, slice.start_index );
1103
1104     while( lpt-- > 0 )
1105     {
1106         CV_READ_SEQ_ELEM( pt, reader );
1107
1108         if( flag == 0 )
1109         {
1110             xi_1 = (double) pt.x;
1111             yi_1 = (double) pt.y;
1112             x0 = xi_1;
1113             y0 = yi_1;
1114             sk1 = 0;
1115             flag = 1;
1116         }
1117         else
1118         {
1119             xi = (double) pt.x;
1120             yi = (double) pt.y;
1121
1122 /****************   edges intersection examination   **************************/
1123             sk = nx * (xi - pt_s.x) + ny * (yi - pt_s.y);
1124             if( (fabs( sk ) < eps && lpt > 0) || sk * sk1 < -eps )
1125             {
1126                 if( fabs( sk ) < eps )
1127                 {
1128                     dxy = xi_1 * yi - xi * yi_1;
1129                     a00 = a00 + dxy;
1130                     dxy = xi * y0 - x0 * yi;
1131                     a00 = a00 + dxy;
1132
1133                     if( p_ind >= p_max )
1134                         icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
1135
1136                     p_are[p_ind] = a00 / 2.;
1137                     p_ind++;
1138                     a00 = 0;
1139                     sk1 = 0;
1140                     x0 = xi;
1141                     y0 = yi;
1142                     dxy = 0;
1143                 }
1144                 else
1145                 {
1146 /*  define intersection point    */
1147                     dv = yi - yi_1;
1148                     du = xi - xi_1;
1149                     dx = ny;
1150                     dy = -nx;
1151                     if( fabs( du ) > eps )
1152                         t = ((yi_1 - pt_s.y) * du + dv * (pt_s.x - xi_1)) /
1153                             (du * dy - dx * dv);
1154                     else
1155                         t = (xi_1 - pt_s.x) / dx;
1156                     if( t > eps && t < 1 - eps )
1157                     {
1158                         x_s = pt_s.x + t * dx;
1159                         y_s = pt_s.y + t * dy;
1160                         dxy = xi_1 * y_s - x_s * yi_1;
1161                         a00 += dxy;
1162                         dxy = x_s * y0 - x0 * y_s;
1163                         a00 += dxy;
1164                         if( p_ind >= p_max )
1165                             icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
1166
1167                         p_are[p_ind] = a00 / 2.;
1168                         p_ind++;
1169
1170                         a00 = 0;
1171                         sk1 = 0;
1172                         x0 = x_s;
1173                         y0 = y_s;
1174                         dxy = x_s * yi - xi * y_s;
1175                     }
1176                 }
1177             }
1178             else
1179                 dxy = xi_1 * yi - xi * yi_1;
1180
1181             a00 += dxy;
1182             xi_1 = xi;
1183             yi_1 = yi;
1184             sk1 = sk;
1185
1186         }
1187     }
1188
1189     xi = x0;
1190     yi = y0;
1191     dxy = xi_1 * yi - xi * yi_1;
1192
1193     a00 += dxy;
1194
1195     if( p_ind >= p_max )
1196         icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
1197
1198     p_are[p_ind] = a00 / 2.;
1199     p_ind++;
1200
1201     // common area calculation
1202     area = 0;
1203     for( i = 0; i < p_ind; i++ )
1204         area += fabs( p_are[i] );
1205
1206     if( p_are1 != NULL )
1207         cvFree( &p_are1 );
1208     else if( p_are2 != NULL )
1209         cvFree( &p_are2 );
1210
1211     return area;
1212 }
1213
1214
1215 /* external contour area function */
1216 CV_IMPL double
1217 cvContourArea( const void *array, CvSlice slice, int oriented )
1218 {
1219     double area = 0;
1220
1221     CvContour contour_header;
1222     CvSeq* contour = 0;
1223     CvSeqBlock block;
1224
1225     if( CV_IS_SEQ( array ))
1226     {
1227         contour = (CvSeq*)array;
1228         if( !CV_IS_SEQ_POLYLINE( contour ))
1229             CV_Error( CV_StsBadArg, "Unsupported sequence type" );
1230     }
1231     else
1232     {
1233         contour = cvPointSeqFromMat( CV_SEQ_KIND_CURVE, array, &contour_header, &block );
1234     }
1235
1236     if( cvSliceLength( slice, contour ) == contour->total )
1237     {
1238         cv::AutoBuffer<double> abuf;
1239         cv::Mat points = cv::cvarrToMat(contour, false, false, 0, &abuf);
1240         return cv::contourArea( points, oriented !=0 );
1241     }
1242
1243     if( CV_SEQ_ELTYPE( contour ) != CV_32SC2 )
1244         CV_Error( CV_StsUnsupportedFormat,
1245         "Only curves with integer coordinates are supported in case of contour slice" );
1246     area = icvContourSecArea( contour, slice );
1247     return oriented ? area : fabs(area);
1248 }
1249
1250
1251 /* calculates length of a curve (e.g. contour perimeter) */
1252 CV_IMPL  double
1253 cvArcLength( const void *array, CvSlice slice, int is_closed )
1254 {
1255     double perimeter = 0;
1256
1257     int i, j = 0, count;
1258     const int N = 16;
1259     float buf[N];
1260     CvMat buffer = cvMat( 1, N, CV_32F, buf );
1261     CvSeqReader reader;
1262     CvContour contour_header;
1263     CvSeq* contour = 0;
1264     CvSeqBlock block;
1265
1266     if( CV_IS_SEQ( array ))
1267     {
1268         contour = (CvSeq*)array;
1269         if( !CV_IS_SEQ_POLYLINE( contour ))
1270             CV_Error( CV_StsBadArg, "Unsupported sequence type" );
1271         if( is_closed < 0 )
1272             is_closed = CV_IS_SEQ_CLOSED( contour );
1273     }
1274     else
1275     {
1276         is_closed = is_closed > 0;
1277         contour = cvPointSeqFromMat(
1278                                     CV_SEQ_KIND_CURVE | (is_closed ? CV_SEQ_FLAG_CLOSED : 0),
1279                                     array, &contour_header, &block );
1280     }
1281
1282     if( contour->total > 1 )
1283     {
1284         int is_float = CV_SEQ_ELTYPE( contour ) == CV_32FC2;
1285
1286         cvStartReadSeq( contour, &reader, 0 );
1287         cvSetSeqReaderPos( &reader, slice.start_index );
1288         count = cvSliceLength( slice, contour );
1289
1290         count -= !is_closed && count == contour->total;
1291
1292         // scroll the reader by 1 point
1293         reader.prev_elem = reader.ptr;
1294         CV_NEXT_SEQ_ELEM( sizeof(CvPoint), reader );
1295
1296         for( i = 0; i < count; i++ )
1297         {
1298             float dx, dy;
1299
1300             if( !is_float )
1301             {
1302                 CvPoint* pt = (CvPoint*)reader.ptr;
1303                 CvPoint* prev_pt = (CvPoint*)reader.prev_elem;
1304
1305                 dx = (float)pt->x - (float)prev_pt->x;
1306                 dy = (float)pt->y - (float)prev_pt->y;
1307             }
1308             else
1309             {
1310                 CvPoint2D32f* pt = (CvPoint2D32f*)reader.ptr;
1311                 CvPoint2D32f* prev_pt = (CvPoint2D32f*)reader.prev_elem;
1312
1313                 dx = pt->x - prev_pt->x;
1314                 dy = pt->y - prev_pt->y;
1315             }
1316
1317             reader.prev_elem = reader.ptr;
1318             CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
1319             // Bugfix by Axel at rubico.com 2010-03-22, affects closed slices only
1320             // wraparound not handled by CV_NEXT_SEQ_ELEM
1321             if( is_closed && i == count - 2 )
1322                 cvSetSeqReaderPos( &reader, slice.start_index );
1323
1324             buffer.data.fl[j] = dx * dx + dy * dy;
1325             if( ++j == N || i == count - 1 )
1326             {
1327                 buffer.cols = j;
1328                 cvPow( &buffer, &buffer, 0.5 );
1329                 for( ; j > 0; j-- )
1330                     perimeter += buffer.data.fl[j-1];
1331             }
1332         }
1333     }
1334
1335     return perimeter;
1336 }
1337
1338
1339 CV_IMPL CvBox2D
1340 cvFitEllipse2( const CvArr* array )
1341 {
1342     cv::AutoBuffer<double> abuf;
1343     cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf);
1344     return cv::fitEllipse(points);
1345 }
1346
1347 /* Calculates bounding rectagnle of a point set or retrieves already calculated */
1348 CV_IMPL  CvRect
1349 cvBoundingRect( CvArr* array, int update )
1350 {
1351     CvRect  rect;
1352     CvContour contour_header;
1353     CvSeq* ptseq = 0;
1354     CvSeqBlock block;
1355
1356     CvMat stub, *mat = 0;
1357     int calculate = update;
1358
1359     if( CV_IS_SEQ( array ))
1360     {
1361         ptseq = (CvSeq*)array;
1362         if( !CV_IS_SEQ_POINT_SET( ptseq ))
1363             CV_Error( CV_StsBadArg, "Unsupported sequence type" );
1364
1365         if( ptseq->header_size < (int)sizeof(CvContour))
1366         {
1367             update = 0;
1368             calculate = 1;
1369         }
1370     }
1371     else
1372     {
1373         mat = cvGetMat( array, &stub );
1374         if( CV_MAT_TYPE(mat->type) == CV_32SC2 ||
1375             CV_MAT_TYPE(mat->type) == CV_32FC2 )
1376         {
1377             ptseq = cvPointSeqFromMat(CV_SEQ_KIND_GENERIC, mat, &contour_header, &block);
1378             mat = 0;
1379         }
1380         else if( CV_MAT_TYPE(mat->type) != CV_8UC1 &&
1381                 CV_MAT_TYPE(mat->type) != CV_8SC1 )
1382             CV_Error( CV_StsUnsupportedFormat,
1383                 "The image/matrix format is not supported by the function" );
1384         update = 0;
1385         calculate = 1;
1386     }
1387
1388     if( !calculate )
1389         return ((CvContour*)ptseq)->rect;
1390
1391     if( mat )
1392     {
1393         rect = cv::maskBoundingRect(cv::cvarrToMat(mat));
1394     }
1395     else if( ptseq->total )
1396     {
1397         cv::AutoBuffer<double> abuf;
1398         rect = cv::pointSetBoundingRect(cv::cvarrToMat(ptseq, false, false, 0, &abuf));
1399     }
1400     if( update )
1401         ((CvContour*)ptseq)->rect = rect;
1402     return rect;
1403 }
1404
1405 /* End of file. */