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 //                           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 OpenCV Foundation 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 OpenCV Foundation 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 namespace cv
44 {
45 
46 struct MinAreaState
47 {
48     int bottom;
49     int left;
50     float height;
51     float width;
52     float base_a;
53     float base_b;
54 };
55 
56 enum { CALIPERS_MAXHEIGHT=0, CALIPERS_MINAREARECT=1, CALIPERS_MAXDIST=2 };
57 
58 /*F///////////////////////////////////////////////////////////////////////////////////////
59  //    Name:    rotatingCalipers
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 */
rotatingCalipers(const Point2f * points,int n,int mode,float * out)92 static void rotatingCalipers( const Point2f* points, int n, int mode, float* out )
93 {
94     float minarea = FLT_MAX;
95     float max_dist = 0;
96     char buffer[32] = {};
97     int i, k;
98     AutoBuffer<float> abuf(n*3);
99     float* inv_vect_length = abuf;
100     Point2f* vect = (Point2f*)(inv_vect_length + n);
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     Point2f 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         Point2f 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./std::sqrt(dx*dx + dy*dy));
142 
143         pt0 = pt;
144     }
145 
146     // find convex hull orientation
147     {
148         double ax = vect[n-1].x;
149         double ay = vect[n-1].y;
150 
151         for( i = 0; i < n; i++ )
152         {
153             double bx = vect[i].x;
154             double by = vect[i].y;
155 
156             double convexity = ax * by - ay * bx;
157 
158             if( convexity != 0 )
159             {
160                 orientation = (convexity > 0) ? 1.f : (-1.f);
161                 break;
162             }
163             ax = bx;
164             ay = by;
165         }
166         CV_Assert( orientation != 0 );
167     }
168     base_a = orientation;
169 
170     /*****************************************************************************************/
171     /*                         init calipers position                                        */
172     seq[0] = bottom;
173     seq[1] = right;
174     seq[2] = top;
175     seq[3] = left;
176     /*****************************************************************************************/
177     /*                         Main loop - evaluate angles and rotate calipers               */
178 
179     /* all of edges will be checked while rotating calipers by 90 degrees */
180     for( k = 0; k < n; k++ )
181     {
182         /* sinus of minimal angle */
183         /*float sinus;*/
184 
185         /* compute cosine of angle between calipers side and polygon edge */
186         /* dp - dot product */
187         float dp0 = base_a * vect[seq[0]].x + base_b * vect[seq[0]].y;
188         float dp1 = -base_b * vect[seq[1]].x + base_a * vect[seq[1]].y;
189         float dp2 = -base_a * vect[seq[2]].x - base_b * vect[seq[2]].y;
190         float dp3 = base_b * vect[seq[3]].x - base_a * vect[seq[3]].y;
191 
192         float cosalpha = dp0 * inv_vect_length[seq[0]];
193         float maxcos = cosalpha;
194 
195         /* number of calipers edges, that has minimal angle with edge */
196         int main_element = 0;
197 
198         /* choose minimal angle */
199         cosalpha = dp1 * inv_vect_length[seq[1]];
200         maxcos = (cosalpha > maxcos) ? (main_element = 1, cosalpha) : maxcos;
201         cosalpha = dp2 * inv_vect_length[seq[2]];
202         maxcos = (cosalpha > maxcos) ? (main_element = 2, cosalpha) : maxcos;
203         cosalpha = dp3 * inv_vect_length[seq[3]];
204         maxcos = (cosalpha > maxcos) ? (main_element = 3, cosalpha) : maxcos;
205 
206         /*rotate calipers*/
207         {
208             //get next base
209             int pindex = seq[main_element];
210             float lead_x = vect[pindex].x*inv_vect_length[pindex];
211             float lead_y = vect[pindex].y*inv_vect_length[pindex];
212             switch( main_element )
213             {
214             case 0:
215                 base_a = lead_x;
216                 base_b = lead_y;
217                 break;
218             case 1:
219                 base_a = lead_y;
220                 base_b = -lead_x;
221                 break;
222             case 2:
223                 base_a = -lead_x;
224                 base_b = -lead_y;
225                 break;
226             case 3:
227                 base_a = -lead_y;
228                 base_b = lead_x;
229                 break;
230             default:
231                 CV_Error(CV_StsError, "main_element should be 0, 1, 2 or 3");
232             }
233         }
234         /* change base point of main edge */
235         seq[main_element] += 1;
236         seq[main_element] = (seq[main_element] == n) ? 0 : seq[main_element];
237 
238         switch (mode)
239         {
240         case CALIPERS_MAXHEIGHT:
241             {
242             /* now main element lies on edge alligned to calipers side */
243 
244             /* find opposite element i.e. transform  */
245             /* 0->2, 1->3, 2->0, 3->1                */
246             int opposite_el = main_element ^ 2;
247 
248             float dx = points[seq[opposite_el]].x - points[seq[main_element]].x;
249             float dy = points[seq[opposite_el]].y - points[seq[main_element]].y;
250             float dist;
251 
252             if( main_element & 1 )
253                 dist = (float)fabs(dx * base_a + dy * base_b);
254             else
255                 dist = (float)fabs(dx * (-base_b) + dy * base_a);
256 
257             if( dist > max_dist )
258                 max_dist = dist;
259             }
260             break;
261         case CALIPERS_MINAREARECT:
262             /* find area of rectangle */
263             {
264             float height;
265             float area;
266 
267             /* find vector left-right */
268             float dx = points[seq[1]].x - points[seq[3]].x;
269             float dy = points[seq[1]].y - points[seq[3]].y;
270 
271             /* dotproduct */
272             float width = dx * base_a + dy * base_b;
273 
274             /* find vector left-right */
275             dx = points[seq[2]].x - points[seq[0]].x;
276             dy = points[seq[2]].y - points[seq[0]].y;
277 
278             /* dotproduct */
279             height = -dx * base_b + dy * base_a;
280 
281             area = width * height;
282             if( area <= minarea )
283             {
284                 float *buf = (float *) buffer;
285 
286                 minarea = area;
287                 /* leftist point */
288                 ((int *) buf)[0] = seq[3];
289                 buf[1] = base_a;
290                 buf[2] = width;
291                 buf[3] = base_b;
292                 buf[4] = height;
293                 /* bottom point */
294                 ((int *) buf)[5] = seq[0];
295                 buf[6] = area;
296             }
297             }
298             break;
299         }                       /*switch */
300     }                           /* for */
301 
302     switch (mode)
303     {
304     case CALIPERS_MINAREARECT:
305         {
306         float *buf = (float *) buffer;
307 
308         float A1 = buf[1];
309         float B1 = buf[3];
310 
311         float A2 = -buf[3];
312         float B2 = buf[1];
313 
314         float C1 = A1 * points[((int *) buf)[0]].x + points[((int *) buf)[0]].y * B1;
315         float C2 = A2 * points[((int *) buf)[5]].x + points[((int *) buf)[5]].y * B2;
316 
317         float idet = 1.f / (A1 * B2 - A2 * B1);
318 
319         float px = (C1 * B2 - C2 * B1) * idet;
320         float py = (A1 * C2 - A2 * C1) * idet;
321 
322         out[0] = px;
323         out[1] = py;
324 
325         out[2] = A1 * buf[2];
326         out[3] = B1 * buf[2];
327 
328         out[4] = A2 * buf[4];
329         out[5] = B2 * buf[4];
330         }
331         break;
332     case CALIPERS_MAXHEIGHT:
333         {
334         out[0] = max_dist;
335         }
336         break;
337     }
338 }
339 
340 }
341 
342 
minAreaRect(InputArray _points)343 cv::RotatedRect cv::minAreaRect( InputArray _points )
344 {
345     Mat hull;
346     Point2f out[3];
347     RotatedRect box;
348 
349     convexHull(_points, hull, true, true);
350 
351     if( hull.depth() != CV_32F )
352     {
353         Mat temp;
354         hull.convertTo(temp, CV_32F);
355         hull = temp;
356     }
357 
358     int n = hull.checkVector(2);
359     const Point2f* hpoints = hull.ptr<Point2f>();
360 
361     if( n > 2 )
362     {
363         rotatingCalipers( hpoints, n, CALIPERS_MINAREARECT, (float*)out );
364         box.center.x = out[0].x + (out[1].x + out[2].x)*0.5f;
365         box.center.y = out[0].y + (out[1].y + out[2].y)*0.5f;
366         box.size.width = (float)std::sqrt((double)out[1].x*out[1].x + (double)out[1].y*out[1].y);
367         box.size.height = (float)std::sqrt((double)out[2].x*out[2].x + (double)out[2].y*out[2].y);
368         box.angle = (float)atan2( (double)out[1].y, (double)out[1].x );
369     }
370     else if( n == 2 )
371     {
372         box.center.x = (hpoints[0].x + hpoints[1].x)*0.5f;
373         box.center.y = (hpoints[0].y + hpoints[1].y)*0.5f;
374         double dx = hpoints[1].x - hpoints[0].x;
375         double dy = hpoints[1].y - hpoints[0].y;
376         box.size.width = (float)std::sqrt(dx*dx + dy*dy);
377         box.size.height = 0;
378         box.angle = (float)atan2( dy, dx );
379     }
380     else
381     {
382         if( n == 1 )
383             box.center = hpoints[0];
384     }
385 
386     box.angle = (float)(box.angle*180/CV_PI);
387     return box;
388 }
389 
390 
391 CV_IMPL CvBox2D
cvMinAreaRect2(const CvArr * array,CvMemStorage *)392 cvMinAreaRect2( const CvArr* array, CvMemStorage* /*storage*/ )
393 {
394     cv::AutoBuffer<double> abuf;
395     cv::Mat points = cv::cvarrToMat(array, false, false, 0, &abuf);
396 
397     cv::RotatedRect rr = cv::minAreaRect(points);
398     return (CvBox2D)rr;
399 }
400 
boxPoints(cv::RotatedRect box,OutputArray _pts)401 void cv::boxPoints(cv::RotatedRect box, OutputArray _pts)
402 {
403     _pts.create(4, 2, CV_32F);
404     Mat pts = _pts.getMat();
405     box.points(pts.ptr<Point2f>());
406 }
407