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