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 43 44CV_IMPL CvRect 45cvMaxRect( const CvRect* rect1, const CvRect* rect2 ) 46{ 47 if( rect1 && rect2 ) 48 { 49 CvRect max_rect; 50 int a, b; 51 52 max_rect.x = a = rect1->x; 53 b = rect2->x; 54 if( max_rect.x > b ) 55 max_rect.x = b; 56 57 max_rect.width = a += rect1->width; 58 b += rect2->width; 59 60 if( max_rect.width < b ) 61 max_rect.width = b; 62 max_rect.width -= max_rect.x; 63 64 max_rect.y = a = rect1->y; 65 b = rect2->y; 66 if( max_rect.y > b ) 67 max_rect.y = b; 68 69 max_rect.height = a += rect1->height; 70 b += rect2->height; 71 72 if( max_rect.height < b ) 73 max_rect.height = b; 74 max_rect.height -= max_rect.y; 75 return max_rect; 76 } 77 else if( rect1 ) 78 return *rect1; 79 else if( rect2 ) 80 return *rect2; 81 else 82 return cvRect(0,0,0,0); 83} 84 85 86CV_IMPL void 87cvBoxPoints( CvBox2D box, CvPoint2D32f pt[4] ) 88{ 89 if( !pt ) 90 CV_Error( CV_StsNullPtr, "NULL vertex array pointer" ); 91 cv::RotatedRect(box).points((cv::Point2f*)pt); 92} 93 94 95double cv::pointPolygonTest( InputArray _contour, Point2f pt, bool measureDist ) 96{ 97 double result = 0; 98 Mat contour = _contour.getMat(); 99 int i, total = contour.checkVector(2), counter = 0; 100 int depth = contour.depth(); 101 CV_Assert( total >= 0 && (depth == CV_32S || depth == CV_32F)); 102 103 bool is_float = depth == CV_32F; 104 double min_dist_num = FLT_MAX, min_dist_denom = 1; 105 Point ip(cvRound(pt.x), cvRound(pt.y)); 106 107 if( total == 0 ) 108 return measureDist ? -DBL_MAX : -1; 109 110 const Point* cnt = contour.ptr<Point>(); 111 const Point2f* cntf = (const Point2f*)cnt; 112 113 if( !is_float && !measureDist && ip.x == pt.x && ip.y == pt.y ) 114 { 115 // the fastest "purely integer" branch 116 Point v0, v = cnt[total-1]; 117 118 for( i = 0; i < total; i++ ) 119 { 120 int dist; 121 v0 = v; 122 v = cnt[i]; 123 124 if( (v0.y <= ip.y && v.y <= ip.y) || 125 (v0.y > ip.y && v.y > ip.y) || 126 (v0.x < ip.x && v.x < ip.x) ) 127 { 128 if( ip.y == v.y && (ip.x == v.x || (ip.y == v0.y && 129 ((v0.x <= ip.x && ip.x <= v.x) || (v.x <= ip.x && ip.x <= v0.x)))) ) 130 return 0; 131 continue; 132 } 133 134 dist = (ip.y - v0.y)*(v.x - v0.x) - (ip.x - v0.x)*(v.y - v0.y); 135 if( dist == 0 ) 136 return 0; 137 if( v.y < v0.y ) 138 dist = -dist; 139 counter += dist > 0; 140 } 141 142 result = counter % 2 == 0 ? -1 : 1; 143 } 144 else 145 { 146 Point2f v0, v; 147 Point iv; 148 149 if( is_float ) 150 { 151 v = cntf[total-1]; 152 } 153 else 154 { 155 v = cnt[total-1]; 156 } 157 158 if( !measureDist ) 159 { 160 for( i = 0; i < total; i++ ) 161 { 162 double dist; 163 v0 = v; 164 if( is_float ) 165 v = cntf[i]; 166 else 167 v = cnt[i]; 168 169 if( (v0.y <= pt.y && v.y <= pt.y) || 170 (v0.y > pt.y && v.y > pt.y) || 171 (v0.x < pt.x && v.x < pt.x) ) 172 { 173 if( pt.y == v.y && (pt.x == v.x || (pt.y == v0.y && 174 ((v0.x <= pt.x && pt.x <= v.x) || (v.x <= pt.x && pt.x <= v0.x)))) ) 175 return 0; 176 continue; 177 } 178 179 dist = (double)(pt.y - v0.y)*(v.x - v0.x) - (double)(pt.x - v0.x)*(v.y - v0.y); 180 if( dist == 0 ) 181 return 0; 182 if( v.y < v0.y ) 183 dist = -dist; 184 counter += dist > 0; 185 } 186 187 result = counter % 2 == 0 ? -1 : 1; 188 } 189 else 190 { 191 for( i = 0; i < total; i++ ) 192 { 193 double dx, dy, dx1, dy1, dx2, dy2, dist_num, dist_denom = 1; 194 195 v0 = v; 196 if( is_float ) 197 v = cntf[i]; 198 else 199 v = cnt[i]; 200 201 dx = v.x - v0.x; dy = v.y - v0.y; 202 dx1 = pt.x - v0.x; dy1 = pt.y - v0.y; 203 dx2 = pt.x - v.x; dy2 = pt.y - v.y; 204 205 if( dx1*dx + dy1*dy <= 0 ) 206 dist_num = dx1*dx1 + dy1*dy1; 207 else if( dx2*dx + dy2*dy >= 0 ) 208 dist_num = dx2*dx2 + dy2*dy2; 209 else 210 { 211 dist_num = (dy1*dx - dx1*dy); 212 dist_num *= dist_num; 213 dist_denom = dx*dx + dy*dy; 214 } 215 216 if( dist_num*min_dist_denom < min_dist_num*dist_denom ) 217 { 218 min_dist_num = dist_num; 219 min_dist_denom = dist_denom; 220 if( min_dist_num == 0 ) 221 break; 222 } 223 224 if( (v0.y <= pt.y && v.y <= pt.y) || 225 (v0.y > pt.y && v.y > pt.y) || 226 (v0.x < pt.x && v.x < pt.x) ) 227 continue; 228 229 dist_num = dy1*dx - dx1*dy; 230 if( dy < 0 ) 231 dist_num = -dist_num; 232 counter += dist_num > 0; 233 } 234 235 result = std::sqrt(min_dist_num/min_dist_denom); 236 if( counter % 2 == 0 ) 237 result = -result; 238 } 239 } 240 241 return result; 242} 243 244 245CV_IMPL double 246cvPointPolygonTest( const CvArr* _contour, CvPoint2D32f pt, int measure_dist ) 247{ 248 cv::AutoBuffer<double> abuf; 249 cv::Mat contour = cv::cvarrToMat(_contour, false, false, 0, &abuf); 250 return cv::pointPolygonTest(contour, pt, measure_dist != 0); 251} 252 253/* 254 This code is described in "Computational Geometry in C" (Second Edition), 255 Chapter 7. It is not written to be comprehensible without the 256 explanation in that book. 257 258 Written by Joseph O'Rourke. 259 Last modified: December 1997 260 Questions to orourke@cs.smith.edu. 261 -------------------------------------------------------------------- 262 This code is Copyright 1997 by Joseph O'Rourke. It may be freely 263 redistributed in its entirety provided that this copyright notice is 264 not removed. 265 -------------------------------------------------------------------- 266 */ 267 268namespace cv 269{ 270typedef enum { Pin, Qin, Unknown } tInFlag; 271 272static int areaSign( Point2f a, Point2f b, Point2f c ) 273{ 274 static const double eps = 1e-5; 275 double area2 = (b.x - a.x) * (double)(c.y - a.y) - (c.x - a.x ) * (double)(b.y - a.y); 276 return area2 > eps ? 1 : area2 < -eps ? -1 : 0; 277} 278 279//--------------------------------------------------------------------- 280// Returns true iff point c lies on the closed segement ab. 281// Assumes it is already known that abc are collinear. 282//--------------------------------------------------------------------- 283static bool between( Point2f a, Point2f b, Point2f c ) 284{ 285 Point2f ba, ca; 286 287 // If ab not vertical, check betweenness on x; else on y. 288 if ( a.x != b.x ) 289 return ((a.x <= c.x) && (c.x <= b.x)) || 290 ((a.x >= c.x) && (c.x >= b.x)); 291 else 292 return ((a.y <= c.y) && (c.y <= b.y)) || 293 ((a.y >= c.y) && (c.y >= b.y)); 294} 295 296static char parallelInt( Point2f a, Point2f b, Point2f c, Point2f d, Point2f& p, Point2f& q ) 297{ 298 char code = 'e'; 299 if( areaSign(a, b, c) != 0 ) 300 code = '0'; 301 else if( between(a, b, c) && between(a, b, d)) 302 p = c, q = d; 303 else if( between(c, d, a) && between(c, d, b)) 304 p = a, q = b; 305 else if( between(a, b, c) && between(c, d, b)) 306 p = c, q = b; 307 else if( between(a, b, c) && between(c, d, a)) 308 p = c, q = a; 309 else if( between(a, b, d) && between(c, d, b)) 310 p = d, q = b; 311 else if( between(a, b, d) && between(c, d, a)) 312 p = d, q = a; 313 else 314 code = '0'; 315 return code; 316} 317 318//--------------------------------------------------------------------- 319// segSegInt: Finds the point of intersection p between two closed 320// segments ab and cd. Returns p and a char with the following meaning: 321// 'e': The segments collinearly overlap, sharing a point. 322// 'v': An endpoint (vertex) of one segment is on the other segment, 323// but 'e' doesn't hold. 324// '1': The segments intersect properly (i.e., they share a point and 325// neither 'v' nor 'e' holds). 326// '0': The segments do not intersect (i.e., they share no points). 327// Note that two collinear segments that share just one point, an endpoint 328// of each, returns 'e' rather than 'v' as one might expect. 329//--------------------------------------------------------------------- 330static char segSegInt( Point2f a, Point2f b, Point2f c, Point2f d, Point2f& p, Point2f& q ) 331{ 332 double s, t; // The two parameters of the parametric eqns. 333 double num, denom; // Numerator and denoninator of equations. 334 char code = '?'; // Return char characterizing intersection. 335 336 denom = a.x * (double)( d.y - c.y ) + 337 b.x * (double)( c.y - d.y ) + 338 d.x * (double)( b.y - a.y ) + 339 c.x * (double)( a.y - b.y ); 340 341 // If denom is zero, then segments are parallel: handle separately. 342 if (denom == 0.0) 343 return parallelInt(a, b, c, d, p, q); 344 345 num = a.x * (double)( d.y - c.y ) + 346 c.x * (double)( a.y - d.y ) + 347 d.x * (double)( c.y - a.y ); 348 if ( (num == 0.0) || (num == denom) ) code = 'v'; 349 s = num / denom; 350 351 num = -( a.x * (double)( c.y - b.y ) + 352 b.x * (double)( a.y - c.y ) + 353 c.x * (double)( b.y - a.y ) ); 354 if ( (num == 0.0) || (num == denom) ) code = 'v'; 355 t = num / denom; 356 357 if ( (0.0 < s) && (s < 1.0) && 358 (0.0 < t) && (t < 1.0) ) 359 code = '1'; 360 else if ( (0.0 > s) || (s > 1.0) || 361 (0.0 > t) || (t > 1.0) ) 362 code = '0'; 363 364 p.x = (float)(a.x + s*(b.x - a.x)); 365 p.y = (float)(a.y + s*(b.y - a.y)); 366 367 return code; 368} 369 370static tInFlag inOut( Point2f p, tInFlag inflag, int aHB, int bHA, Point2f*& result ) 371{ 372 if( p != result[-1] ) 373 *result++ = p; 374 // Update inflag. 375 return aHB > 0 ? Pin : bHA > 0 ? Qin : inflag; 376} 377 378//--------------------------------------------------------------------- 379// Advances and prints out an inside vertex if appropriate. 380//--------------------------------------------------------------------- 381static int advance( int a, int *aa, int n, bool inside, Point2f v, Point2f*& result ) 382{ 383 if( inside && v != result[-1] ) 384 *result++ = v; 385 (*aa)++; 386 return (a+1) % n; 387} 388 389static void addSharedSeg( Point2f p, Point2f q, Point2f*& result ) 390{ 391 if( p != result[-1] ) 392 *result++ = p; 393 if( q != result[-1] ) 394 *result++ = q; 395} 396 397 398static int intersectConvexConvex_( const Point2f* P, int n, const Point2f* Q, int m, 399 Point2f* result, float* _area ) 400{ 401 Point2f* result0 = result; 402 // P has n vertices, Q has m vertices. 403 int a=0, b=0; // indices on P and Q (resp.) 404 Point2f Origin(0,0); 405 tInFlag inflag=Unknown; // {Pin, Qin, Unknown}: which inside 406 int aa=0, ba=0; // # advances on a & b indices (after 1st inter.) 407 bool FirstPoint=true;// Is this the first point? (used to initialize). 408 Point2f p0; // The first point. 409 *result++ = Point2f(FLT_MAX, FLT_MAX); 410 411 do 412 { 413 // Computations of key variables. 414 int a1 = (a + n - 1) % n; // a-1, b-1 (resp.) 415 int b1 = (b + m - 1) % m; 416 417 Point2f A = P[a] - P[a1], B = Q[b] - Q[b1]; // directed edges on P and Q (resp.) 418 419 int cross = areaSign( Origin, A, B ); // sign of z-component of A x B 420 int aHB = areaSign( Q[b1], Q[b], P[a] ); // a in H(b). 421 int bHA = areaSign( P[a1], P[a], Q[b] ); // b in H(A); 422 423 // If A & B intersect, update inflag. 424 Point2f p, q; 425 int code = segSegInt( P[a1], P[a], Q[b1], Q[b], p, q ); 426 if( code == '1' || code == 'v' ) 427 { 428 if( inflag == Unknown && FirstPoint ) 429 { 430 aa = ba = 0; 431 FirstPoint = false; 432 p0 = p; 433 *result++ = p; 434 } 435 inflag = inOut( p, inflag, aHB, bHA, result ); 436 } 437 438 //-----Advance rules----- 439 440 // Special case: A & B overlap and oppositely oriented. 441 if( code == 'e' && A.ddot(B) < 0 ) 442 { 443 addSharedSeg( p, q, result ); 444 return (int)(result - result0); 445 } 446 447 // Special case: A & B parallel and separated. 448 if( cross == 0 && aHB < 0 && bHA < 0 ) 449 return (int)(result - result0); 450 451 // Special case: A & B collinear. 452 else if ( cross == 0 && aHB == 0 && bHA == 0 ) { 453 // Advance but do not output point. 454 if ( inflag == Pin ) 455 b = advance( b, &ba, m, inflag == Qin, Q[b], result ); 456 else 457 a = advance( a, &aa, n, inflag == Pin, P[a], result ); 458 } 459 460 // Generic cases. 461 else if( cross >= 0 ) 462 { 463 if( bHA > 0) 464 a = advance( a, &aa, n, inflag == Pin, P[a], result ); 465 else 466 b = advance( b, &ba, m, inflag == Qin, Q[b], result ); 467 } 468 else 469 { 470 if( aHB > 0) 471 b = advance( b, &ba, m, inflag == Qin, Q[b], result ); 472 else 473 a = advance( a, &aa, n, inflag == Pin, P[a], result ); 474 } 475 // Quit when both adv. indices have cycled, or one has cycled twice. 476 } 477 while ( ((aa < n) || (ba < m)) && (aa < 2*n) && (ba < 2*m) ); 478 479 // Deal with special cases: not implemented. 480 if( inflag == Unknown ) 481 { 482 // The boundaries of P and Q do not cross. 483 // ... 484 } 485 486 int i, nr = (int)(result - result0); 487 double area = 0; 488 Point2f prev = result0[nr-1]; 489 for( i = 1; i < nr; i++ ) 490 { 491 result0[i-1] = result0[i]; 492 area += (double)prev.x*result0[i].y - (double)prev.y*result0[i].x; 493 prev = result0[i]; 494 } 495 496 *_area = (float)(area*0.5); 497 498 if( result0[nr-2] == result0[0] && nr > 1 ) 499 nr--; 500 return nr-1; 501} 502 503} 504 505float cv::intersectConvexConvex( InputArray _p1, InputArray _p2, OutputArray _p12, bool handleNested ) 506{ 507 Mat p1 = _p1.getMat(), p2 = _p2.getMat(); 508 CV_Assert( p1.depth() == CV_32S || p1.depth() == CV_32F ); 509 CV_Assert( p2.depth() == CV_32S || p2.depth() == CV_32F ); 510 511 int n = p1.checkVector(2, p1.depth(), true); 512 int m = p2.checkVector(2, p2.depth(), true); 513 514 CV_Assert( n >= 0 && m >= 0 ); 515 516 if( n < 2 || m < 2 ) 517 { 518 _p12.release(); 519 return 0.f; 520 } 521 522 AutoBuffer<Point2f> _result(n*2 + m*2 + 1); 523 Point2f *fp1 = _result, *fp2 = fp1 + n; 524 Point2f* result = fp2 + m; 525 int orientation = 0; 526 527 for( int k = 1; k <= 2; k++ ) 528 { 529 Mat& p = k == 1 ? p1 : p2; 530 int len = k == 1 ? n : m; 531 Point2f* dst = k == 1 ? fp1 : fp2; 532 533 Mat temp(p.size(), CV_MAKETYPE(CV_32F, p.channels()), dst); 534 p.convertTo(temp, CV_32F); 535 CV_Assert( temp.ptr<Point2f>() == dst ); 536 Point2f diff0 = dst[0] - dst[len-1]; 537 for( int i = 1; i < len; i++ ) 538 { 539 double s = diff0.cross(dst[i] - dst[i-1]); 540 if( s != 0 ) 541 { 542 if( s < 0 ) 543 { 544 orientation++; 545 flip( temp, temp, temp.rows > 1 ? 0 : 1 ); 546 } 547 break; 548 } 549 } 550 } 551 552 float area = 0.f; 553 int nr = intersectConvexConvex_(fp1, n, fp2, m, result, &area); 554 if( nr == 0 ) 555 { 556 if( !handleNested ) 557 { 558 _p12.release(); 559 return 0.f; 560 } 561 562 if( pointPolygonTest(_InputArray(fp1, n), fp2[0], false) >= 0 ) 563 { 564 result = fp2; 565 nr = m; 566 } 567 else if( pointPolygonTest(_InputArray(fp2, n), fp1[0], false) >= 0 ) 568 { 569 result = fp1; 570 nr = n; 571 } 572 else 573 { 574 _p12.release(); 575 return 0.f; 576 } 577 area = (float)contourArea(_InputArray(result, nr), false); 578 } 579 580 if( _p12.needed() ) 581 { 582 Mat temp(nr, 1, CV_32FC2, result); 583 // if both input contours were reflected, 584 // let's orient the result as the input vectors 585 if( orientation == 2 ) 586 flip(temp, temp, 0); 587 588 temp.copyTo(_p12); 589 } 590 return (float)fabs(area); 591} 592