git.maemo.org Git - opencv/blob - cv/src/cvrotcalipers.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 "_cv.h"
  42
  43 typedef struct
  44 {
  45     int bottom;
  46     int left;
  47     float height;
  48     float width;
  49     float base_a;
  50     float base_b;
  51 }
  52 icvMinAreaState;
  53
  54 #define CV_CALIPERS_MAXHEIGHT      0
  55 #define CV_CALIPERS_MINAREARECT    1
  56 #define CV_CALIPERS_MAXDIST        2
  57
  58 /*F///////////////////////////////////////////////////////////////////////////////////////
  59 //    Name:    icvRotatingCalipers
  60 //    Purpose:
  61 //      Rotating calipers algorithm with some applications
  62 //
  63 //    Context:
  64 //    Parameters:
  65 //      points      - convex hull vertices ( any orientation )
  66 //      n           - number of vertices
  67 //      mode        - concrete application of algorithm
  68 //                    can be  CV_CALIPERS_MAXDIST   or
  69 //                            CV_CALIPERS_MINAREARECT
  70 //      left, bottom, right, top - indexes of extremal points
  71 //      out         - output info.
  72 //                    In case CV_CALIPERS_MAXDIST it points to float value -
  73 //                    maximal height of polygon.
  74 //                    In case CV_CALIPERS_MINAREARECT
  75 //                    ((CvPoint2D32f*)out)[0] - corner
  76 //                    ((CvPoint2D32f*)out)[1] - vector1
  77 //                    ((CvPoint2D32f*)out)[0] - corner2
  78 //
  79 //                      ^
  80 //                      |
  81 //              vector2 |
  82 //                      |
  83 //                      |____________\
  84 //                    corner         /
  85 //                               vector1
  86 //
  87 //    Returns:
  88 //    Notes:
  89 //F*/
  90
  91 /* we will use usual cartesian coordinates */
  92 static void
  93 icvRotatingCalipers( CvPoint2D32f* points, int n, int mode, float* out )
  94 {
  95     float minarea = FLT_MAX;
  96     float max_dist = 0;
  97     char buffer[32];
  98     int i, k;
  99     CvPoint2D32f* vect = (CvPoint2D32f*)cvAlloc( n * sizeof(vect[0]) );
 100     float* inv_vect_length = (float*)cvAlloc( n * sizeof(inv_vect_length[0]) );
 101     int left = 0, bottom = 0, right = 0, top = 0;
 102     int seq[4] = { -1, -1, -1, -1 };
 103
 104     /* rotating calipers sides will always have coordinates
 105        (a,b) (-b,a) (-a,-b) (b, -a)
 106      */
 107     /* this is a first base bector (a,b) initialized by (1,0) */
 108     float orientation = 0;
 109     float base_a;
 110     float base_b = 0;
 111
 112     float left_x, right_x, top_y, bottom_y;
 113     CvPoint2D32f pt0 = points[0];
 114
 115     left_x = right_x = pt0.x;
 116     top_y = bottom_y = pt0.y;
 117
 118     for( i = 0; i < n; i++ )
 119     {
 120         double dx, dy;
 121
 122         if( pt0.x < left_x )
 123             left_x = pt0.x, left = i;
 124
 125         if( pt0.x > right_x )
 126             right_x = pt0.x, right = i;
 127
 128         if( pt0.y > top_y )
 129             top_y = pt0.y, top = i;
 130
 131         if( pt0.y < bottom_y )
 132             bottom_y = pt0.y, bottom = i;
 133
 134         CvPoint2D32f pt = points[(i+1) & (i+1 < n ? -1 : 0)];
 135
 136         dx = pt.x - pt0.x;
 137         dy = pt.y - pt0.y;
 138
 139         vect[i].x = (float)dx;
 140         vect[i].y = (float)dy;
 141         inv_vect_length[i] = (float)(1./sqrt(dx*dx + dy*dy));
 142
 143         pt0 = pt;
 144     }
 145
 146     //cvbInvSqrt( inv_vect_length, inv_vect_length, n );
 147
 148     /* find convex hull orientation */
 149     {
 150         double ax = vect[n-1].x;
 151         double ay = vect[n-1].y;
 152
 153         for( i = 0; i < n; i++ )
 154         {
 155             double bx = vect[i].x;
 156             double by = vect[i].y;
 157
 158             double convexity = ax * by - ay * bx;
 159
 160             if( convexity != 0 )
 161             {
 162                 orientation = (convexity > 0) ? 1.f : (-1.f);
 163                 break;
 164             }
 165             ax = bx;
 166             ay = by;
 167         }
 168         assert( orientation != 0 );
 169     }
 170     base_a = orientation;
 171
 172 /*****************************************************************************************/
 173 /*                         init calipers position                                        */
 174     seq[0] = bottom;
 175     seq[1] = right;
 176     seq[2] = top;
 177     seq[3] = left;
 178 /*****************************************************************************************/
 179 /*                         Main loop - evaluate angles and rotate calipers               */
 180
 181     /* all of edges will be checked while rotating calipers by 90 degrees */
 182     for( k = 0; k < n; k++ )
 183     {
 184         /* sinus of minimal angle */
 185         /*float sinus;*/
 186
 187         /* compute cosine of angle between calipers side and polygon edge */
 188         /* dp - dot product */
 189         float dp0 = base_a * vect[seq[0]].x + base_b * vect[seq[0]].y;
 190         float dp1 = -base_b * vect[seq[1]].x + base_a * vect[seq[1]].y;
 191         float dp2 = -base_a * vect[seq[2]].x - base_b * vect[seq[2]].y;
 192         float dp3 = base_b * vect[seq[3]].x - base_a * vect[seq[3]].y;
 193
 194         float cosalpha = dp0 * inv_vect_length[seq[0]];
 195         float maxcos = cosalpha;
 196
 197         /* number of calipers edges, that has minimal angle with edge */
 198         int main_element = 0;
 199
 200         /* choose minimal angle */
 201         cosalpha = dp1 * inv_vect_length[seq[1]];
 202         maxcos = (cosalpha > maxcos) ? (main_element = 1, cosalpha) : maxcos;
 203         cosalpha = dp2 * inv_vect_length[seq[2]];
 204         maxcos = (cosalpha > maxcos) ? (main_element = 2, cosalpha) : maxcos;
 205         cosalpha = dp3 * inv_vect_length[seq[3]];
 206         maxcos = (cosalpha > maxcos) ? (main_element = 3, cosalpha) : maxcos;
 207
 208         /*rotate calipers*/
 209         {
 210             //get next base
 211             int pindex = seq[main_element];
 212             float lead_x = vect[pindex].x*inv_vect_length[pindex];
 213             float lead_y = vect[pindex].y*inv_vect_length[pindex];
 214             switch( main_element )
 215             {
 216             case 0:
 217                 base_a = lead_x;
 218                 base_b = lead_y;
 219                 break;
 220             case 1:
 221                 base_a = lead_y;
 222                 base_b = -lead_x;
 223                 break;
 224             case 2:
 225                 base_a = -lead_x;
 226                 base_b = -lead_y;
 227                 break;
 228             case 3:
 229                 base_a = -lead_y;
 230                 base_b = lead_x;
 231                 break;
 232             default: assert(0);
 233             }
 234         }
 235         /* change base point of main edge */
 236         seq[main_element] += 1;
 237         seq[main_element] = (seq[main_element] == n) ? 0 : seq[main_element];
 238
 239
 240         switch (mode)
 241         {
 242         case CV_CALIPERS_MAXHEIGHT:
 243             {
 244                 /* now main element lies on edge alligned to calipers side */
 245
 246                 /* find opposite element i.e. transform  */
 247                 /* 0->2, 1->3, 2->0, 3->1                */
 248                 int opposite_el = main_element ^ 2;
 249
 250                 float dx = points[seq[opposite_el]].x - points[seq[main_element]].x;
 251                 float dy = points[seq[opposite_el]].y - points[seq[main_element]].y;
 252                 float dist;
 253
 254                 if( main_element & 1 )
 255                     dist = (float)fabs(dx * base_a + dy * base_b);
 256                 else
 257                     dist = (float)fabs(dx * (-base_b) + dy * base_a);
 258
 259                 if( dist > max_dist )
 260                     max_dist = dist;
 261
 262                 break;
 263             }
 264         case CV_CALIPERS_MINAREARECT:
 265             /* find area of rectangle */
 266             {
 267                 float height;
 268                 float area;
 269
 270                 /* find vector left-right */
 271                 float dx = points[seq[1]].x - points[seq[3]].x;
 272                 float dy = points[seq[1]].y - points[seq[3]].y;
 273
 274                 /* dotproduct */
 275                 float width = dx * base_a + dy * base_b;
 276
 277                 /* find vector left-right */
 278                 dx = points[seq[2]].x - points[seq[0]].x;
 279                 dy = points[seq[2]].y - points[seq[0]].y;
 280
 281                 /* dotproduct */
 282                 height = -dx * base_b + dy * base_a;
 283
 284                 area = width * height;
 285                 if( area <= minarea )
 286                 {
 287                     float *buf = (float *) buffer;
 288
 289                     minarea = area;
 290                     /* leftist point */
 291                     ((int *) buf)[0] = seq[3];
 292                     buf[1] = base_a;
 293                     buf[2] = width;
 294                     buf[3] = base_b;
 295                     buf[4] = height;
 296                     /* bottom point */
 297                     ((int *) buf)[5] = seq[0];
 298                     buf[6] = area;
 299                 }
 300                 break;
 301             }
 302         }                       /*switch */
 303     }                           /* for */
 304
 305     switch (mode)
 306     {
 307     case CV_CALIPERS_MINAREARECT:
 308         {
 309             float *buf = (float *) buffer;
 310
 311             float A1 = buf[1];
 312             float B1 = buf[3];
 313
 314             float A2 = -buf[3];
 315             float B2 = buf[1];
 316
 317             float C1 = A1 * points[((int *) buf)[0]].x + points[((int *) buf)[0]].y * B1;
 318             float C2 = A2 * points[((int *) buf)[5]].x + points[((int *) buf)[5]].y * B2;
 319
 320             float idet = 1.f / (A1 * B2 - A2 * B1);
 321
 322             float px = (C1 * B2 - C2 * B1) * idet;
 323             float py = (A1 * C2 - A2 * C1) * idet;
 324
 325             out[0] = px;
 326             out[1] = py;
 327
 328             out[2] = A1 * buf[2];
 329             out[3] = B1 * buf[2];
 330
 331             out[4] = A2 * buf[4];
 332             out[5] = B2 * buf[4];
 333         }
 334         break;
 335     case CV_CALIPERS_MAXHEIGHT:
 336         {
 337             out[0] = max_dist;
 338         }
 339         break;
 340     }
 341
 342     cvFree( &vect );
 343     cvFree( &inv_vect_length );
 344 }
 345
 346
 347 CV_IMPL  CvBox2D
 348 cvMinAreaRect2( const CvArr* array, CvMemStorage* storage )
 349 {
 350     CvMemStorage* temp_storage = 0;
 351     CvBox2D box;
 352     CvPoint2D32f* points = 0;
 353
 354     CV_FUNCNAME( "cvMinAreaRect2" );
 355
 356     memset(&box, 0, sizeof(box));
 357
 358     __BEGIN__;
 359
 360     int i, n;
 361     CvSeqReader reader;
 362     CvContour contour_header;
 363     CvSeqBlock block;
 364     CvSeq* ptseq = (CvSeq*)array;
 365     CvPoint2D32f out[3];
 366
 367     if( CV_IS_SEQ(ptseq) )
 368     {
 369         if( !CV_IS_SEQ_POINT_SET(ptseq) &&
 370             (CV_SEQ_KIND(ptseq) != CV_SEQ_KIND_CURVE || !CV_IS_SEQ_CONVEX(ptseq) ||
 371             CV_SEQ_ELTYPE(ptseq) != CV_SEQ_ELTYPE_PPOINT ))
 372             CV_ERROR( CV_StsUnsupportedFormat,
 373                 "Input sequence must consist of 2d points or pointers to 2d points" );
 374         if( !storage )
 375             storage = ptseq->storage;
 376     }
 377     else
 378     {
 379         CV_CALL( ptseq = cvPointSeqFromMat(
 380             CV_SEQ_KIND_GENERIC, array, &contour_header, &block ));
 381     }
 382
 383     if( storage )
 384     {
 385         CV_CALL( temp_storage = cvCreateChildMemStorage( storage ));
 386     }
 387     else
 388     {
 389         CV_CALL( temp_storage = cvCreateMemStorage(1 << 10));
 390     }
 391
 392     if( !CV_IS_SEQ_CONVEX( ptseq ))
 393     {
 394         CV_CALL( ptseq = cvConvexHull2( ptseq, temp_storage, CV_CLOCKWISE, 1 ));
 395     }
 396     else if( !CV_IS_SEQ_POINT_SET( ptseq ))
 397     {
 398         CvSeqWriter writer;
 399
 400         if( !CV_IS_SEQ(ptseq->v_prev) || !CV_IS_SEQ_POINT_SET(ptseq->v_prev))
 401             CV_ERROR( CV_StsBadArg,
 402             "Convex hull must have valid pointer to point sequence stored in v_prev" );
 403         cvStartReadSeq( ptseq, &reader );
 404         cvStartWriteSeq( CV_SEQ_KIND_CURVE|CV_SEQ_FLAG_CONVEX|CV_SEQ_ELTYPE(ptseq->v_prev),
 405                          sizeof(CvContour), CV_ELEM_SIZE(ptseq->v_prev->flags),
 406                          temp_storage, &writer );
 407
 408         for( i = 0; i < ptseq->total; i++ )
 409         {
 410             CvPoint pt = **(CvPoint**)(reader.ptr);
 411             CV_WRITE_SEQ_ELEM( pt, writer );
 412         }
 413
 414         ptseq = cvEndWriteSeq( &writer );
 415     }
 416
 417     n = ptseq->total;
 418
 419     CV_CALL( points = (CvPoint2D32f*)cvAlloc( n*sizeof(points[0]) ));
 420     cvStartReadSeq( ptseq, &reader );
 421
 422     if( CV_SEQ_ELTYPE( ptseq ) == CV_32SC2 )
 423     {
 424         for( i = 0; i < n; i++ )
 425         {
 426             CvPoint pt;
 427             CV_READ_SEQ_ELEM( pt, reader );
 428             points[i].x = (float)pt.x;
 429             points[i].y = (float)pt.y;
 430         }
 431     }
 432     else
 433     {
 434         for( i = 0; i < n; i++ )
 435         {
 436             CV_READ_SEQ_ELEM( points[i], reader );
 437         }
 438     }
 439
 440     if( n > 2 )
 441     {
 442         icvRotatingCalipers( points, n, CV_CALIPERS_MINAREARECT, (float*)out );
 443         box.center.x = out[0].x + (out[1].x + out[2].x)*0.5f;
 444         box.center.y = out[0].y + (out[1].y + out[2].y)*0.5f;
 445         box.size.height = (float)sqrt((double)out[1].x*out[1].x + (double)out[1].y*out[1].y);
 446         box.size.width = (float)sqrt((double)out[2].x*out[2].x + (double)out[2].y*out[2].y);
 447         box.angle = (float)atan2( -(double)out[1].y, (double)out[1].x );
 448     }
 449     else if( n == 2 )
 450     {
 451         box.center.x = (points[0].x + points[1].x)*0.5f;
 452         box.center.y = (points[0].y + points[1].y)*0.5f;
 453         double dx = points[1].x - points[0].x;
 454         double dy = points[1].y - points[0].y;
 455         box.size.height = (float)sqrt(dx*dx + dy*dy);
 456         box.size.width = 0;
 457         box.angle = (float)atan2( -dy, dx );
 458     }
 459     else
 460     {
 461         if( n == 1 )
 462             box.center = points[0];
 463     }
 464
 465     box.angle = (float)(box.angle*180/CV_PI);
 466
 467     __END__;
 468
 469     cvReleaseMemStorage( &temp_storage );
 470     cvFree( &points );
 471
 472     return box;
 473 }
 474