rotcalipers.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 //                           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 OpenCV Foundation 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 OpenCV Foundation 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 struct MinAreaState
00047 {
00048     int bottom;
00049     int left;
00050     float height;
00051     float width;
00052     float base_a;
00053     float base_b;
00054 };
00055 
00056 enum { CALIPERS_MAXHEIGHT=0, CALIPERS_MINAREARECT=1, CALIPERS_MAXDIST=2 };
00057 
00058 /*F///////////////////////////////////////////////////////////////////////////////////////
00059  //    Name:    rotatingCalipers
00060  //    Purpose:
00061  //      Rotating calipers algorithm with some applications
00062  //
00063  //    Context:
00064  //    Parameters:
00065  //      points      - convex hull vertices ( any orientation )
00066  //      n           - number of vertices
00067  //      mode        - concrete application of algorithm
00068  //                    can be  CV_CALIPERS_MAXDIST   or
00069  //                            CV_CALIPERS_MINAREARECT
00070  //      left, bottom, right, top - indexes of extremal points
00071  //      out         - output info.
00072  //                    In case CV_CALIPERS_MAXDIST it points to float value -
00073  //                    maximal height of polygon.
00074  //                    In case CV_CALIPERS_MINAREARECT
00075  //                    ((CvPoint2D32f*)out)[0] - corner
00076  //                    ((CvPoint2D32f*)out)[1] - vector1
00077  //                    ((CvPoint2D32f*)out)[0] - corner2
00078  //
00079  //                      ^
00080  //                      |
00081  //              vector2 |
00082  //                      |
00083  //                      |____________\
00084  //                    corner         /
00085  //                               vector1
00086  //
00087  //    Returns:
00088  //    Notes:
00089  //F*/
00090 
00091 /* we will use usual cartesian coordinates */
00092 static void rotatingCalipers( const Point2f* points, int n, int mode, float* out )
00093 {
00094     float minarea = FLT_MAX;
00095     float max_dist = 0;
00096     char buffer[32] = {};
00097     int i, k;
00098     AutoBuffer<float> abuf(n*3);
00099     float* inv_vect_length = abuf;
00100     Point2f* vect = (Point2f*)(inv_vect_length + n);
00101     int left = 0, bottom = 0, right = 0, top = 0;
00102     int seq[4] = { -1, -1, -1, -1 };
00103 
00104     /* rotating calipers sides will always have coordinates
00105      (a,b) (-b,a) (-a,-b) (b, -a)
00106      */
00107     /* this is a first base bector (a,b) initialized by (1,0) */
00108     float orientation = 0;
00109     float base_a;
00110     float base_b = 0;
00111 
00112     float left_x, right_x, top_y, bottom_y;
00113     Point2f pt0 = points[0];
00114 
00115     left_x = right_x = pt0.x;
00116     top_y = bottom_y = pt0.y;
00117 
00118     for( i = 0; i < n; i++ )
00119     {
00120         double dx, dy;
00121 
00122         if( pt0.x < left_x )
00123             left_x = pt0.x, left = i;
00124 
00125         if( pt0.x > right_x )
00126             right_x = pt0.x, right = i;
00127 
00128         if( pt0.y > top_y )
00129             top_y = pt0.y, top = i;
00130 
00131         if( pt0.y < bottom_y )
00132             bottom_y = pt0.y, bottom = i;
00133 
00134         Point2f pt = points[(i+1) & (i+1 < n ? -1 : 0)];
00135 
00136         dx = pt.x - pt0.x;
00137         dy = pt.y - pt0.y;
00138 
00139         vect[i].x = (float)dx;
00140         vect[i].y = (float)dy;
00141         inv_vect_length[i] = (float)(1./std::sqrt(dx*dx + dy*dy));
00142 
00143         pt0 = pt;
00144     }
00145 
00146     // find convex hull orientation
00147     {
00148         double ax = vect[n-1].x;
00149         double ay = vect[n-1].y;
00150 
00151         for( i = 0; i < n; i++ )
00152         {
00153             double bx = vect[i].x;
00154             double by = vect[i].y;
00155 
00156             double convexity = ax * by - ay * bx;
00157 
00158             if( convexity != 0 )
00159             {
00160                 orientation = (convexity > 0) ? 1.f : (-1.f);
00161                 break;
00162             }
00163             ax = bx;
00164             ay = by;
00165         }
00166         CV_Assert( orientation != 0 );
00167     }
00168     base_a = orientation;
00169 
00170     /*****************************************************************************************/
00171     /*                         init calipers position                                        */
00172     seq[0] = bottom;
00173     seq[1] = right;
00174     seq[2] = top;
00175     seq[3] = left;
00176     /*****************************************************************************************/
00177     /*                         Main loop - evaluate angles and rotate calipers               */
00178 
00179     /* all of edges will be checked while rotating calipers by 90 degrees */
00180     for( k = 0; k < n; k++ )
00181     {
00182         /* sinus of minimal angle */
00183         /*float sinus;*/
00184 
00185         /* compute cosine of angle between calipers side and polygon edge */
00186         /* dp - dot product */
00187         float dp[4] = {
00188             +base_a * vect[seq[0]].x + base_b * vect[seq[0]].y,
00189             -base_b * vect[seq[1]].x + base_a * vect[seq[1]].y,
00190             -base_a * vect[seq[2]].x - base_b * vect[seq[2]].y,
00191             +base_b * vect[seq[3]].x - base_a * vect[seq[3]].y,
00192         };
00193 
00194         float maxcos = dp[0] * inv_vect_length[seq[0]];
00195 
00196         /* number of calipers edges, that has minimal angle with edge */
00197         int main_element = 0;
00198 
00199         /* choose minimal angle */
00200         for ( i = 1; i < 4; ++i )
00201         {
00202             float cosalpha = dp[i] * inv_vect_length[seq[i]];
00203             if (cosalpha > maxcos)
00204             {
00205                 main_element = i;
00206                 maxcos = cosalpha;
00207             }
00208         }
00209 
00210         /*rotate calipers*/
00211         {
00212             //get next base
00213             int pindex = seq[main_element];
00214             float lead_x = vect[pindex].x*inv_vect_length[pindex];
00215             float lead_y = vect[pindex].y*inv_vect_length[pindex];
00216             switch( main_element )
00217             {
00218             case 0:
00219                 base_a = lead_x;
00220                 base_b = lead_y;
00221                 break;
00222             case 1:
00223                 base_a = lead_y;
00224                 base_b = -lead_x;
00225                 break;
00226             case 2:
00227                 base_a = -lead_x;
00228                 base_b = -lead_y;
00229                 break;
00230             case 3:
00231                 base_a = -lead_y;
00232                 base_b = lead_x;
00233                 break;
00234             default:
00235                 CV_Error(CV_StsError, "main_element should be 0, 1, 2 or 3");
00236             }
00237         }
00238         /* change base point of main edge */
00239         seq[main_element] += 1;
00240         seq[main_element] = (seq[main_element] == n) ? 0 : seq[main_element];
00241 
00242         switch (mode)
00243         {
00244         case CALIPERS_MAXHEIGHT:
00245             {
00246             /* now main element lies on edge alligned to calipers side */
00247 
00248             /* find opposite element i.e. transform  */
00249             /* 0->2, 1->3, 2->0, 3->1                */
00250             int opposite_el = main_element ^ 2;
00251 
00252             float dx = points[seq[opposite_el]].x - points[seq[main_element]].x;
00253             float dy = points[seq[opposite_el]].y - points[seq[main_element]].y;
00254             float dist;
00255 
00256             if( main_element & 1 )
00257                 dist = (float)fabs(dx * base_a + dy * base_b);
00258             else
00259                 dist = (float)fabs(dx * (-base_b) + dy * base_a);
00260 
00261             if( dist > max_dist )
00262                 max_dist = dist;
00263             }
00264             break;
00265         case CALIPERS_MINAREARECT:
00266             /* find area of rectangle */
00267             {
00268             float height;
00269             float area;
00270 
00271             /* find vector left-right */
00272             float dx = points[seq[1]].x - points[seq[3]].x;
00273             float dy = points[seq[1]].y - points[seq[3]].y;
00274 
00275             /* dotproduct */
00276             float width = dx * base_a + dy * base_b;
00277 
00278             /* find vector left-right */
00279             dx = points[seq[2]].x - points[seq[0]].x;
00280             dy = points[seq[2]].y - points[seq[0]].y;
00281 
00282             /* dotproduct */
00283             height = -dx * base_b + dy * base_a;
00284 
00285             area = width * height;
00286             if( area <= minarea )
00287             {
00288                 float *buf = (float *) buffer;
00289 
00290                 minarea = area;
00291                 /* leftist point */
00292                 ((int *) buf)[0] = seq[3];
00293                 buf[1] = base_a;
00294                 buf[2] = width;
00295                 buf[3] = base_b;
00296                 buf[4] = height;
00297                 /* bottom point */
00298                 ((int *) buf)[5] = seq[0];
00299                 buf[6] = area;
00300             }
00301             }
00302             break;
00303         }                       /*switch */
00304     }                           /* for */
00305 
00306     switch (mode)
00307     {
00308     case CALIPERS_MINAREARECT:
00309         {
00310         float *buf = (float *) buffer;
00311 
00312         float A1 = buf[1];
00313         float B1 = buf[3];
00314 
00315         float A2 = -buf[3];
00316         float B2 = buf[1];
00317 
00318         float C1 = A1 * points[((int *) buf)[0]].x + points[((int *) buf)[0]].y * B1;
00319         float C2 = A2 * points[((int *) buf)[5]].x + points[((int *) buf)[5]].y * B2;
00320 
00321         float idet = 1.f / (A1 * B2 - A2 * B1);
00322 
00323         float px = (C1 * B2 - C2 * B1) * idet;
00324         float py = (A1 * C2 - A2 * C1) * idet;
00325 
00326         out[0] = px;
00327         out[1] = py;
00328 
00329         out[2] = A1 * buf[2];
00330         out[3] = B1 * buf[2];
00331 
00332         out[4] = A2 * buf[4];
00333         out[5] = B2 * buf[4];
00334         }
00335         break;
00336     case CALIPERS_MAXHEIGHT:
00337         {
00338         out[0] = max_dist;
00339         }
00340         break;
00341     }
00342 }
00343 
00344 }
00345 
00346 
00347 cv::RotatedRect cv::minAreaRect( InputArray _points )
00348 {
00349     Mat hull;
00350     Point2f  out[3];
00351     RotatedRect box;
00352 
00353     convexHull(_points, hull, true, true);
00354 
00355     if( hull.depth() != CV_32F )
00356     {
00357         Mat temp;
00358         hull.convertTo(temp, CV_32F);
00359         hull = temp;
00360     }
00361 
00362     int n = hull.checkVector(2);
00363     const Point2f * hpoints = hull.ptr<Point2f >();
00364 
00365     if( n > 2 )
00366     {
00367         rotatingCalipers( hpoints, n, CALIPERS_MINAREARECT, (float*)out );
00368         box.center.x = out[0].x + (out[1].x + out[2].x)*0.5f;
00369         box.center.y = out[0].y + (out[1].y + out[2].y)*0.5f;
00370         box.size.width = (float)std::sqrt((double)out[1].x*out[1].x + (double)out[1].y*out[1].y);
00371         box.size.height = (float)std::sqrt((double)out[2].x*out[2].x + (double)out[2].y*out[2].y);
00372         box.angle = (float)atan2( (double)out[1].y, (double)out[1].x );
00373     }
00374     else if( n == 2 )
00375     {
00376         box.center.x = (hpoints[0].x + hpoints[1].x)*0.5f;
00377         box.center.y = (hpoints[0].y + hpoints[1].y)*0.5f;
00378         double dx = hpoints[1].x - hpoints[0].x;
00379         double dy = hpoints[1].y - hpoints[0].y;
00380         box.size.width = (float)std::sqrt(dx*dx + dy*dy);
00381         box.size.height = 0;
00382         box.angle = (float)atan2( dy, dx );
00383     }
00384     else
00385     {
00386         if( n == 1 )
00387             box.center = hpoints[0];
00388     }
00389 
00390     box.angle = (float)(box.angle*180/CV_PI);
00391     return box;
00392 }
00393 
00394 
00395 CV_IMPL CvBox2D 
00396 cvMinAreaRect2( const CvArr* array, CvMemStorage* /*storage*/ )
00397 {
00398     cv::AutoBuffer<double> abuf;
00399     cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf);
00400 
00401     cv::RotatedRect rr = cv::minAreaRect(points);
00402     return (CvBox2D )rr;
00403 }
00404 
00405 void cv::boxPoints(cv::RotatedRect box, OutputArray _pts)
00406 {
00407     _pts.create(4, 2, CV_32F);
00408     Mat pts = _pts.getMat();
00409     box.points(pts.ptr<Point2f >());
00410 }
00411