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 475
