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
43 typedef struct
44 {
45 int bottom;
46 int left;
47 float height;
48 float width;
49 float base_a;
50 float base_b;
51 }
52 icvMinAreaState;
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 */
92 static void
icvRotatingCalipers(CvPoint2D32f * points,int n,int mode,float * out)93 icvRotatingCalipers( 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
347 CV_IMPL CvBox2D
cvMinAreaRect2(const CvArr * array,CvMemStorage * storage)348 cvMinAreaRect2( 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