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
43typedef struct
44{
45    int bottom;
46    int left;
47    float height;
48    float width;
49    float base_a;
50    float base_b;
51}
52icvMinAreaState;
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 */
92static void
93icvRotatingCalipers( 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
347CV_IMPL  CvBox2D
348cvMinAreaRect2( 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