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 
44 CV_IMPL CvRect
cvMaxRect(const CvRect * rect1,const CvRect * rect2)45 cvMaxRect( 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 
86 CV_IMPL void
cvBoxPoints(CvBox2D box,CvPoint2D32f pt[4])87 cvBoxPoints( 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 
pointPolygonTest(InputArray _contour,Point2f pt,bool measureDist)95 double 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 
245 CV_IMPL double
cvPointPolygonTest(const CvArr * _contour,CvPoint2D32f pt,int measure_dist)246 cvPointPolygonTest( 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 
268 namespace cv
269 {
270 typedef enum { Pin, Qin, Unknown } tInFlag;
271 
areaSign(Point2f a,Point2f b,Point2f c)272 static 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 //---------------------------------------------------------------------
between(Point2f a,Point2f b,Point2f c)283 static 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 
parallelInt(Point2f a,Point2f b,Point2f c,Point2f d,Point2f & p,Point2f & q)296 static 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 //---------------------------------------------------------------------
segSegInt(Point2f a,Point2f b,Point2f c,Point2f d,Point2f & p,Point2f & q)330 static 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 
inOut(Point2f p,tInFlag inflag,int aHB,int bHA,Point2f * & result)370 static 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 //---------------------------------------------------------------------
advance(int a,int * aa,int n,bool inside,Point2f v,Point2f * & result)381 static 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 
addSharedSeg(Point2f p,Point2f q,Point2f * & result)389 static 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 
intersectConvexConvex_(const Point2f * P,int n,const Point2f * Q,int m,Point2f * result,float * _area)398 static 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 
intersectConvexConvex(InputArray _p1,InputArray _p2,OutputArray _p12,bool handleNested)505 float 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