geometry.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 //                        Intel 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 Intel Corporation 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 Intel Corporation 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 
00044 CV_IMPL CvRect 
00045 cvMaxRect( const CvRect * rect1, const CvRect * rect2 )
00046 {
00047     if( rect1 && rect2 )
00048     {
00049         CvRect  max_rect;
00050         int a, b;
00051 
00052         max_rect.x = a = rect1->x;
00053         b = rect2->x;
00054         if( max_rect.x > b )
00055             max_rect.x = b;
00056 
00057         max_rect.width = a += rect1->width;
00058         b += rect2->width;
00059 
00060         if( max_rect.width < b )
00061             max_rect.width = b;
00062         max_rect.width -= max_rect.x;
00063 
00064         max_rect.y = a = rect1->y;
00065         b = rect2->y;
00066         if( max_rect.y > b )
00067             max_rect.y = b;
00068 
00069         max_rect.height = a += rect1->height;
00070         b += rect2->height;
00071 
00072         if( max_rect.height < b )
00073             max_rect.height = b;
00074         max_rect.height -= max_rect.y;
00075         return max_rect;
00076     }
00077     else if( rect1 )
00078         return *rect1;
00079     else if( rect2 )
00080         return *rect2;
00081     else
00082         return cvRect(0,0,0,0);
00083 }
00084 
00085 
00086 CV_IMPL void
00087 cvBoxPoints( CvBox2D  box, CvPoint2D32f pt[4] )
00088 {
00089     if( !pt )
00090         CV_Error( CV_StsNullPtr, "NULL vertex array pointer" );
00091     cv::RotatedRect(box).points((cv::Point2f *)pt);
00092 }
00093 
00094 
00095 double cv::pointPolygonTest( InputArray _contour, Point2f  pt, bool measureDist )
00096 {
00097     double result = 0;
00098     Mat contour = _contour.getMat();
00099     int i, total = contour.checkVector(2), counter = 0;
00100     int depth = contour.depth();
00101     CV_Assert( total >= 0 && (depth == CV_32S || depth == CV_32F));
00102 
00103     bool is_float = depth == CV_32F;
00104     double min_dist_num = FLT_MAX, min_dist_denom = 1;
00105     Point ip(cvRound(pt.x), cvRound(pt.y));
00106 
00107     if( total == 0 )
00108         return measureDist ? -DBL_MAX : -1;
00109 
00110     const Point* cnt = contour.ptr<Point>();
00111     const Point2f * cntf = (const Point2f *)cnt;
00112 
00113     if( !is_float && !measureDist && ip.x == pt.x && ip.y == pt.y )
00114     {
00115         // the fastest "purely integer" branch
00116         Point v0, v = cnt[total-1];
00117 
00118         for( i = 0; i < total; i++ )
00119         {
00120             int dist;
00121             v0 = v;
00122             v = cnt[i];
00123 
00124             if( (v0.y <= ip.y && v.y <= ip.y) ||
00125                (v0.y > ip.y && v.y > ip.y) ||
00126                (v0.x < ip.x && v.x < ip.x) )
00127             {
00128                 if( ip.y == v.y && (ip.x == v.x || (ip.y == v0.y &&
00129                     ((v0.x <= ip.x && ip.x <= v.x) || (v.x <= ip.x && ip.x <= v0.x)))) )
00130                     return 0;
00131                 continue;
00132             }
00133 
00134             dist = (ip.y - v0.y)*(v.x - v0.x) - (ip.x - v0.x)*(v.y - v0.y);
00135             if( dist == 0 )
00136                 return 0;
00137             if( v.y < v0.y )
00138                 dist = -dist;
00139             counter += dist > 0;
00140         }
00141 
00142         result = counter % 2 == 0 ? -1 : 1;
00143     }
00144     else
00145     {
00146         Point2f  v0, v;
00147         Point iv;
00148 
00149         if( is_float )
00150         {
00151             v = cntf[total-1];
00152         }
00153         else
00154         {
00155             v = cnt[total-1];
00156         }
00157 
00158         if( !measureDist )
00159         {
00160             for( i = 0; i < total; i++ )
00161             {
00162                 double dist;
00163                 v0 = v;
00164                 if( is_float )
00165                     v = cntf[i];
00166                 else
00167                     v = cnt[i];
00168 
00169                 if( (v0.y <= pt.y && v.y <= pt.y) ||
00170                    (v0.y > pt.y && v.y > pt.y) ||
00171                    (v0.x < pt.x && v.x < pt.x) )
00172                 {
00173                     if( pt.y == v.y && (pt.x == v.x || (pt.y == v0.y &&
00174                         ((v0.x <= pt.x && pt.x <= v.x) || (v.x <= pt.x && pt.x <= v0.x)))) )
00175                         return 0;
00176                     continue;
00177                 }
00178 
00179                 dist = (double)(pt.y - v0.y)*(v.x - v0.x) - (double)(pt.x - v0.x)*(v.y - v0.y);
00180                 if( dist == 0 )
00181                     return 0;
00182                 if( v.y < v0.y )
00183                     dist = -dist;
00184                 counter += dist > 0;
00185             }
00186 
00187             result = counter % 2 == 0 ? -1 : 1;
00188         }
00189         else
00190         {
00191             for( i = 0; i < total; i++ )
00192             {
00193                 double dx, dy, dx1, dy1, dx2, dy2, dist_num, dist_denom = 1;
00194 
00195                 v0 = v;
00196                 if( is_float )
00197                     v = cntf[i];
00198                 else
00199                     v = cnt[i];
00200 
00201                 dx = v.x - v0.x; dy = v.y - v0.y;
00202                 dx1 = pt.x - v0.x; dy1 = pt.y - v0.y;
00203                 dx2 = pt.x - v.x; dy2 = pt.y - v.y;
00204 
00205                 if( dx1*dx + dy1*dy <= 0 )
00206                     dist_num = dx1*dx1 + dy1*dy1;
00207                 else if( dx2*dx + dy2*dy >= 0 )
00208                     dist_num = dx2*dx2 + dy2*dy2;
00209                 else
00210                 {
00211                     dist_num = (dy1*dx - dx1*dy);
00212                     dist_num *= dist_num;
00213                     dist_denom = dx*dx + dy*dy;
00214                 }
00215 
00216                 if( dist_num*min_dist_denom < min_dist_num*dist_denom )
00217                 {
00218                     min_dist_num = dist_num;
00219                     min_dist_denom = dist_denom;
00220                     if( min_dist_num == 0 )
00221                         break;
00222                 }
00223 
00224                 if( (v0.y <= pt.y && v.y <= pt.y) ||
00225                    (v0.y > pt.y && v.y > pt.y) ||
00226                    (v0.x < pt.x && v.x < pt.x) )
00227                     continue;
00228 
00229                 dist_num = dy1*dx - dx1*dy;
00230                 if( dy < 0 )
00231                     dist_num = -dist_num;
00232                 counter += dist_num > 0;
00233             }
00234 
00235             result = std::sqrt(min_dist_num/min_dist_denom);
00236             if( counter % 2 == 0 )
00237                 result = -result;
00238         }
00239     }
00240 
00241     return result;
00242 }
00243 
00244 
00245 CV_IMPL double
00246 cvPointPolygonTest( const CvArr* _contour, CvPoint2D32f pt, int measure_dist )
00247 {
00248     cv::AutoBuffer<double> abuf;
00249     cv::Mat contour = cv::cvarrToMat(_contour, false, false, 0, &abuf);
00250     return cv::pointPolygonTest(contour, pt, measure_dist != 0);
00251 }
00252 
00253 /*
00254  This code is described in "Computational Geometry in C" (Second Edition),
00255  Chapter 7.  It is not written to be comprehensible without the
00256  explanation in that book.
00257 
00258  Written by Joseph O'Rourke.
00259  Last modified: December 1997
00260  Questions to orourke@cs.smith.edu.
00261  --------------------------------------------------------------------
00262  This code is Copyright 1997 by Joseph O'Rourke.  It may be freely
00263  redistributed in its entirety provided that this copyright notice is
00264  not removed.
00265  --------------------------------------------------------------------
00266  */
00267 
00268 namespace cv
00269 {
00270 typedef enum { Pin, Qin, Unknown } tInFlag;
00271 
00272 static int areaSign( Point2f a, Point2f b, Point2f c )
00273 {
00274     static const double eps = 1e-5;
00275     double area2 = (b.x - a.x) * (double)(c.y - a.y) - (c.x - a.x ) * (double)(b.y - a.y);
00276     return area2 > eps ? 1 : area2 < -eps ? -1 : 0;
00277 }
00278 
00279 //---------------------------------------------------------------------
00280 // Returns true iff point c lies on the closed segement ab.
00281 // Assumes it is already known that abc are collinear.
00282 //---------------------------------------------------------------------
00283 static bool between( Point2f a, Point2f b, Point2f c )
00284 {
00285     Point2f ba, ca;
00286 
00287     // If ab not vertical, check betweenness on x; else on y.
00288     if ( a.x != b.x )
00289         return ((a.x <= c.x) && (c.x <= b.x)) ||
00290         ((a.x >= c.x) && (c.x >= b.x));
00291     else
00292         return ((a.y <= c.y) && (c.y <= b.y)) ||
00293         ((a.y >= c.y) && (c.y >= b.y));
00294 }
00295 
00296 static char parallelInt( Point2f a, Point2f b, Point2f c, Point2f d, Point2f& p, Point2f& q )
00297 {
00298     char code = 'e';
00299     if( areaSign(a, b, c) != 0 )
00300         code = '0';
00301     else if( between(a, b, c) && between(a, b, d))
00302         p = c, q = d;
00303     else if( between(c, d, a) && between(c, d, b))
00304         p = a, q = b;
00305     else if( between(a, b, c) && between(c, d, b))
00306         p = c, q = b;
00307     else if( between(a, b, c) && between(c, d, a))
00308         p = c, q = a;
00309     else if( between(a, b, d) && between(c, d, b))
00310         p = d, q = b;
00311     else if( between(a, b, d) && between(c, d, a))
00312         p = d, q = a;
00313     else
00314         code = '0';
00315     return code;
00316 }
00317 
00318 //---------------------------------------------------------------------
00319 // segSegInt: Finds the point of intersection p between two closed
00320 // segments ab and cd.  Returns p and a char with the following meaning:
00321 // 'e': The segments collinearly overlap, sharing a point.
00322 // 'v': An endpoint (vertex) of one segment is on the other segment,
00323 // but 'e' doesn't hold.
00324 // '1': The segments intersect properly (i.e., they share a point and
00325 // neither 'v' nor 'e' holds).
00326 // '0': The segments do not intersect (i.e., they share no points).
00327 // Note that two collinear segments that share just one point, an endpoint
00328 // of each, returns 'e' rather than 'v' as one might expect.
00329 //---------------------------------------------------------------------
00330 static char segSegInt( Point2f a, Point2f b, Point2f c, Point2f d, Point2f& p, Point2f& q )
00331 {
00332     double  s, t;       // The two parameters of the parametric eqns.
00333     double num, denom;  // Numerator and denoninator of equations.
00334     char code = '?';    // Return char characterizing intersection.
00335 
00336     denom = a.x * (double)( d.y - c.y ) +
00337     b.x * (double)( c.y - d.y ) +
00338     d.x * (double)( b.y - a.y ) +
00339     c.x * (double)( a.y - b.y );
00340 
00341     // If denom is zero, then segments are parallel: handle separately.
00342     if (denom == 0.0)
00343         return parallelInt(a, b, c, d, p, q);
00344 
00345     num =    a.x * (double)( d.y - c.y ) +
00346     c.x * (double)( a.y - d.y ) +
00347     d.x * (double)( c.y - a.y );
00348     if ( (num == 0.0) || (num == denom) ) code = 'v';
00349     s = num / denom;
00350 
00351     num = -( a.x * (double)( c.y - b.y ) +
00352             b.x * (double)( a.y - c.y ) +
00353             c.x * (double)( b.y - a.y ) );
00354     if ( (num == 0.0) || (num == denom) ) code = 'v';
00355     t = num / denom;
00356 
00357     if      ( (0.0 < s) && (s < 1.0) &&
00358              (0.0 < t) && (t < 1.0) )
00359         code = '1';
00360     else if ( (0.0 > s) || (s > 1.0) ||
00361              (0.0 > t) || (t > 1.0) )
00362         code = '0';
00363 
00364     p.x = (float)(a.x + s*(b.x - a.x));
00365     p.y = (float)(a.y + s*(b.y - a.y));
00366 
00367     return code;
00368 }
00369 
00370 static tInFlag inOut( Point2f p, tInFlag inflag, int aHB, int bHA, Point2f*& result )
00371 {
00372     if( p != result[-1] )
00373         *result++ = p;
00374     // Update inflag.
00375     return aHB > 0 ? Pin : bHA > 0 ? Qin : inflag;
00376 }
00377 
00378 //---------------------------------------------------------------------
00379 // Advances and prints out an inside vertex if appropriate.
00380 //---------------------------------------------------------------------
00381 static int advance( int a, int *aa, int n, bool inside, Point2f v, Point2f*& result )
00382 {
00383     if( inside && v != result[-1] )
00384         *result++ = v;
00385     (*aa)++;
00386     return  (a+1) % n;
00387 }
00388 
00389 static void addSharedSeg( Point2f p, Point2f q, Point2f*& result )
00390 {
00391     if( p != result[-1] )
00392         *result++ = p;
00393     if( q != result[-1] )
00394         *result++ = q;
00395 }
00396 
00397 
00398 static int intersectConvexConvex_( const Point2f* P, int n, const Point2f* Q, int m,
00399                                    Point2f* result, float* _area )
00400 {
00401     Point2f* result0 = result;
00402     // P has n vertices, Q has m vertices.
00403     int     a=0, b=0;       // indices on P and Q (resp.)
00404     Point2f Origin(0,0);
00405     tInFlag inflag=Unknown; // {Pin, Qin, Unknown}: which inside
00406     int     aa=0, ba=0;     // # advances on a & b indices (after 1st inter.)
00407     bool    FirstPoint=true;// Is this the first point? (used to initialize).
00408     Point2f p0;             // The first point.
00409     *result++ = Point2f(FLT_MAX, FLT_MAX);
00410 
00411     do
00412     {
00413         // Computations of key variables.
00414         int a1 = (a + n - 1) % n; // a-1, b-1 (resp.)
00415         int b1 = (b + m - 1) % m;
00416 
00417         Point2f A = P[a] - P[a1], B = Q[b] - Q[b1]; // directed edges on P and Q (resp.)
00418 
00419         int cross = areaSign( Origin, A, B );    // sign of z-component of A x B
00420         int aHB = areaSign( Q[b1], Q[b], P[a] ); // a in H(b).
00421         int bHA = areaSign( P[a1], P[a], Q[b] ); // b in H(A);
00422 
00423         // If A & B intersect, update inflag.
00424         Point2f p, q;
00425         int code = segSegInt( P[a1], P[a], Q[b1], Q[b], p, q );
00426         if( code == '1' || code == 'v' )
00427         {
00428             if( inflag == Unknown && FirstPoint )
00429             {
00430                 aa = ba = 0;
00431                 FirstPoint = false;
00432                 p0 = p;
00433                 *result++ = p;
00434             }
00435             inflag = inOut( p, inflag, aHB, bHA, result );
00436         }
00437 
00438         //-----Advance rules-----
00439 
00440         // Special case: A & B overlap and oppositely oriented.
00441         if( code == 'e' && A.ddot(B) < 0 )
00442         {
00443             addSharedSeg( p, q, result );
00444             return (int)(result - result0);
00445         }
00446 
00447         // Special case: A & B parallel and separated.
00448         if( cross == 0 && aHB < 0 && bHA < 0 )
00449             return (int)(result - result0);
00450 
00451         // Special case: A & B collinear.
00452         else if ( cross == 0 && aHB == 0 && bHA == 0 ) {
00453             // Advance but do not output point.
00454             if ( inflag == Pin )
00455                 b = advance( b, &ba, m, inflag == Qin, Q[b], result );
00456             else
00457                 a = advance( a, &aa, n, inflag == Pin, P[a], result );
00458         }
00459 
00460         // Generic cases.
00461         else if( cross >= 0 )
00462         {
00463             if( bHA > 0)
00464                 a = advance( a, &aa, n, inflag == Pin, P[a], result );
00465             else
00466                 b = advance( b, &ba, m, inflag == Qin, Q[b], result );
00467         }
00468         else
00469         {
00470             if( aHB > 0)
00471                 b = advance( b, &ba, m, inflag == Qin, Q[b], result );
00472             else
00473                 a = advance( a, &aa, n, inflag == Pin, P[a], result );
00474         }
00475         // Quit when both adv. indices have cycled, or one has cycled twice.
00476     }
00477     while ( ((aa < n) || (ba < m)) && (aa < 2*n) && (ba < 2*m) );
00478 
00479     // Deal with special cases: not implemented.
00480     if( inflag == Unknown )
00481     {
00482         // The boundaries of P and Q do not cross.
00483         // ...
00484     }
00485 
00486     int i, nr = (int)(result - result0);
00487     double area = 0;
00488     Point2f prev = result0[nr-1];
00489     for( i = 1; i < nr; i++ )
00490     {
00491         result0[i-1] = result0[i];
00492         area += (double)prev.x*result0[i].y - (double)prev.y*result0[i].x;
00493         prev = result0[i];
00494     }
00495 
00496     *_area = (float)(area*0.5);
00497 
00498     if( result0[nr-2] == result0[0] && nr > 1 )
00499         nr--;
00500     return nr-1;
00501 }
00502 
00503 }
00504 
00505 float cv::intersectConvexConvex( InputArray _p1, InputArray _p2, OutputArray _p12, bool handleNested )
00506 {
00507     Mat p1 = _p1.getMat(), p2 = _p2.getMat();
00508     CV_Assert( p1.depth() == CV_32S || p1.depth() == CV_32F );
00509     CV_Assert( p2.depth() == CV_32S || p2.depth() == CV_32F );
00510 
00511     int n = p1.checkVector(2, p1.depth(), true);
00512     int m = p2.checkVector(2, p2.depth(), true);
00513 
00514     CV_Assert( n >= 0 && m >= 0 );
00515 
00516     if( n < 2 || m < 2 )
00517     {
00518         _p12.release();
00519         return 0.f;
00520     }
00521 
00522     AutoBuffer<Point2f> _result(n*2 + m*2 + 1);
00523     Point2f  *fp1 = _result, *fp2 = fp1 + n;
00524     Point2f * result = fp2 + m;
00525     int orientation = 0;
00526 
00527     for( int k = 1; k <= 2; k++ )
00528     {
00529         Mat& p = k == 1 ? p1 : p2;
00530         int len = k == 1 ? n : m;
00531         Point2f * dst = k == 1 ? fp1 : fp2;
00532 
00533         Mat temp(p.size(), CV_MAKETYPE(CV_32F, p.channels()), dst);
00534         p.convertTo(temp, CV_32F);
00535         CV_Assert( temp.ptr<Point2f >() == dst );
00536         Point2f  diff0 = dst[0] - dst[len-1];
00537         for( int i = 1; i < len; i++ )
00538         {
00539             double s = diff0.cross(dst[i] - dst[i-1]);
00540             if( s != 0 )
00541             {
00542                 if( s < 0 )
00543                 {
00544                     orientation++;
00545                     flip( temp, temp, temp.rows > 1 ? 0 : 1 );
00546                 }
00547                 break;
00548             }
00549         }
00550     }
00551 
00552     float area = 0.f;
00553     int nr = intersectConvexConvex_(fp1, n, fp2, m, result, &area);
00554     if( nr == 0 )
00555     {
00556         if( !handleNested )
00557         {
00558             _p12.release();
00559             return 0.f;
00560         }
00561 
00562         if( pointPolygonTest(_InputArray(fp1, n), fp2[0], false) >= 0 )
00563         {
00564             result = fp2;
00565             nr = m;
00566         }
00567         else if( pointPolygonTest(_InputArray(fp2, n), fp1[0], false) >= 0 )
00568         {
00569             result = fp1;
00570             nr = n;
00571         }
00572         else
00573         {
00574             _p12.release();
00575             return 0.f;
00576         }
00577         area = (float)contourArea(_InputArray(result, nr), false);
00578     }
00579 
00580     if( _p12.needed() )
00581     {
00582         Mat temp(nr, 1, CV_32FC2, result);
00583         // if both input contours were reflected,
00584         // let's orient the result as the input vectors
00585         if( orientation == 2 )
00586             flip(temp, temp, 0);
00587 
00588         temp.copyTo(_p12);
00589     }
00590     return (float)fabs(area);
00591 }
00592