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 /* calculates length of a curve (e.g. contour perimeter) */
44 CV_IMPL double
cvArcLength(const void * array,CvSlice slice,int is_closed)45 cvArcLength( const void *array, CvSlice slice, int is_closed )
46 {
47 double perimeter = 0;
48
49 CV_FUNCNAME( "cvArcLength" );
50
51 __BEGIN__;
52
53 int i, j = 0, count;
54 const int N = 16;
55 float buf[N];
56 CvMat buffer = cvMat( 1, N, CV_32F, buf );
57 CvSeqReader reader;
58 CvContour contour_header;
59 CvSeq* contour = 0;
60 CvSeqBlock block;
61
62 if( CV_IS_SEQ( array ))
63 {
64 contour = (CvSeq*)array;
65 if( !CV_IS_SEQ_POLYLINE( contour ))
66 CV_ERROR( CV_StsBadArg, "Unsupported sequence type" );
67 if( is_closed < 0 )
68 is_closed = CV_IS_SEQ_CLOSED( contour );
69 }
70 else
71 {
72 is_closed = is_closed > 0;
73 CV_CALL( contour = cvPointSeqFromMat(
74 CV_SEQ_KIND_CURVE | (is_closed ? CV_SEQ_FLAG_CLOSED : 0),
75 array, &contour_header, &block ));
76 }
77
78 if( contour->total > 1 )
79 {
80 int is_float = CV_SEQ_ELTYPE( contour ) == CV_32FC2;
81
82 cvStartReadSeq( contour, &reader, 0 );
83 cvSetSeqReaderPos( &reader, slice.start_index );
84 count = cvSliceLength( slice, contour );
85
86 count -= !is_closed && count == contour->total;
87
88 /* scroll the reader by 1 point */
89 reader.prev_elem = reader.ptr;
90 CV_NEXT_SEQ_ELEM( sizeof(CvPoint), reader );
91
92 for( i = 0; i < count; i++ )
93 {
94 float dx, dy;
95
96 if( !is_float )
97 {
98 CvPoint* pt = (CvPoint*)reader.ptr;
99 CvPoint* prev_pt = (CvPoint*)reader.prev_elem;
100
101 dx = (float)pt->x - (float)prev_pt->x;
102 dy = (float)pt->y - (float)prev_pt->y;
103 }
104 else
105 {
106 CvPoint2D32f* pt = (CvPoint2D32f*)reader.ptr;
107 CvPoint2D32f* prev_pt = (CvPoint2D32f*)reader.prev_elem;
108
109 dx = pt->x - prev_pt->x;
110 dy = pt->y - prev_pt->y;
111 }
112
113 reader.prev_elem = reader.ptr;
114 CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
115
116 buffer.data.fl[j] = dx * dx + dy * dy;
117 if( ++j == N || i == count - 1 )
118 {
119 buffer.cols = j;
120 cvPow( &buffer, &buffer, 0.5 );
121 for( ; j > 0; j-- )
122 perimeter += buffer.data.fl[j-1];
123 }
124 }
125 }
126
127 __END__;
128
129 return perimeter;
130 }
131
132
133 static CvStatus
icvFindCircle(CvPoint2D32f pt0,CvPoint2D32f pt1,CvPoint2D32f pt2,CvPoint2D32f * center,float * radius)134 icvFindCircle( CvPoint2D32f pt0, CvPoint2D32f pt1,
135 CvPoint2D32f pt2, CvPoint2D32f * center, float *radius )
136 {
137 double x1 = (pt0.x + pt1.x) * 0.5;
138 double dy1 = pt0.x - pt1.x;
139 double x2 = (pt1.x + pt2.x) * 0.5;
140 double dy2 = pt1.x - pt2.x;
141 double y1 = (pt0.y + pt1.y) * 0.5;
142 double dx1 = pt1.y - pt0.y;
143 double y2 = (pt1.y + pt2.y) * 0.5;
144 double dx2 = pt2.y - pt1.y;
145 double t = 0;
146
147 CvStatus result = CV_OK;
148
149 if( icvIntersectLines( x1, dx1, y1, dy1, x2, dx2, y2, dy2, &t ) >= 0 )
150 {
151 center->x = (float) (x2 + dx2 * t);
152 center->y = (float) (y2 + dy2 * t);
153 *radius = (float) icvDistanceL2_32f( *center, pt0 );
154 }
155 else
156 {
157 center->x = center->y = 0.f;
158 radius = 0;
159 result = CV_NOTDEFINED_ERR;
160 }
161
162 return result;
163 }
164
165
icvIsPtInCircle(CvPoint2D32f pt,CvPoint2D32f center,float radius)166 CV_INLINE double icvIsPtInCircle( CvPoint2D32f pt, CvPoint2D32f center, float radius )
167 {
168 double dx = pt.x - center.x;
169 double dy = pt.y - center.y;
170 return (double)radius*radius - dx*dx - dy*dy;
171 }
172
173
174 static int
icvFindEnslosingCicle4pts_32f(CvPoint2D32f * pts,CvPoint2D32f * _center,float * _radius)175 icvFindEnslosingCicle4pts_32f( CvPoint2D32f * pts, CvPoint2D32f * _center, float *_radius )
176 {
177 int shuffles[4][4] = { {0, 1, 2, 3}, {0, 1, 3, 2}, {2, 3, 0, 1}, {2, 3, 1, 0} };
178
179 int idxs[4] = { 0, 1, 2, 3 };
180 int i, j, k = 1, mi = 0;
181 float max_dist = 0;
182 CvPoint2D32f center;
183 CvPoint2D32f min_center;
184 float radius, min_radius = FLT_MAX;
185 CvPoint2D32f res_pts[4];
186
187 center = min_center = pts[0];
188 radius = 1.f;
189
190 for( i = 0; i < 4; i++ )
191 for( j = i + 1; j < 4; j++ )
192 {
193 float dist = icvDistanceL2_32f( pts[i], pts[j] );
194
195 if( max_dist < dist )
196 {
197 max_dist = dist;
198 idxs[0] = i;
199 idxs[1] = j;
200 }
201 }
202
203 if( max_dist == 0 )
204 goto function_exit;
205
206 k = 2;
207 for( i = 0; i < 4; i++ )
208 {
209 for( j = 0; j < k; j++ )
210 if( i == idxs[j] )
211 break;
212 if( j == k )
213 idxs[k++] = i;
214 }
215
216 center = cvPoint2D32f( (pts[idxs[0]].x + pts[idxs[1]].x)*0.5f,
217 (pts[idxs[0]].y + pts[idxs[1]].y)*0.5f );
218 radius = (float)(icvDistanceL2_32f( pts[idxs[0]], center )*1.03);
219 if( radius < 1.f )
220 radius = 1.f;
221
222 if( icvIsPtInCircle( pts[idxs[2]], center, radius ) >= 0 &&
223 icvIsPtInCircle( pts[idxs[3]], center, radius ) >= 0 )
224 {
225 k = 2; //rand()%2+2;
226 }
227 else
228 {
229 mi = -1;
230 for( i = 0; i < 4; i++ )
231 {
232 if( icvFindCircle( pts[shuffles[i][0]], pts[shuffles[i][1]],
233 pts[shuffles[i][2]], ¢er, &radius ) >= 0 )
234 {
235 radius *= 1.03f;
236 if( radius < 2.f )
237 radius = 2.f;
238
239 if( icvIsPtInCircle( pts[shuffles[i][3]], center, radius ) >= 0 &&
240 min_radius > radius )
241 {
242 min_radius = radius;
243 min_center = center;
244 mi = i;
245 }
246 }
247 }
248 assert( mi >= 0 );
249 if( mi < 0 )
250 mi = 0;
251 k = 3;
252 center = min_center;
253 radius = min_radius;
254 for( i = 0; i < 4; i++ )
255 idxs[i] = shuffles[mi][i];
256 }
257
258 function_exit:
259
260 *_center = center;
261 *_radius = radius;
262
263 /* reorder output points */
264 for( i = 0; i < 4; i++ )
265 res_pts[i] = pts[idxs[i]];
266
267 for( i = 0; i < 4; i++ )
268 {
269 pts[i] = res_pts[i];
270 assert( icvIsPtInCircle( pts[i], center, radius ) >= 0 );
271 }
272
273 return k;
274 }
275
276
277 CV_IMPL int
cvMinEnclosingCircle(const void * array,CvPoint2D32f * _center,float * _radius)278 cvMinEnclosingCircle( const void* array, CvPoint2D32f * _center, float *_radius )
279 {
280 const int max_iters = 100;
281 const float eps = FLT_EPSILON*2;
282 CvPoint2D32f center = { 0, 0 };
283 float radius = 0;
284 int result = 0;
285
286 if( _center )
287 _center->x = _center->y = 0.f;
288 if( _radius )
289 *_radius = 0;
290
291 CV_FUNCNAME( "cvMinEnclosingCircle" );
292
293 __BEGIN__;
294
295 CvSeqReader reader;
296 int i, k, count;
297 CvPoint2D32f pts[8];
298 CvContour contour_header;
299 CvSeqBlock block;
300 CvSeq* sequence = 0;
301 int is_float;
302
303 if( !_center || !_radius )
304 CV_ERROR( CV_StsNullPtr, "Null center or radius pointers" );
305
306 if( CV_IS_SEQ(array) )
307 {
308 sequence = (CvSeq*)array;
309 if( !CV_IS_SEQ_POINT_SET( sequence ))
310 CV_ERROR( CV_StsBadArg, "The passed sequence is not a valid contour" );
311 }
312 else
313 {
314 CV_CALL( sequence = cvPointSeqFromMat(
315 CV_SEQ_KIND_GENERIC, array, &contour_header, &block ));
316 }
317
318 if( sequence->total <= 0 )
319 CV_ERROR_FROM_STATUS( CV_BADSIZE_ERR );
320
321 CV_CALL( cvStartReadSeq( sequence, &reader, 0 ));
322
323 count = sequence->total;
324 is_float = CV_SEQ_ELTYPE(sequence) == CV_32FC2;
325
326 if( !is_float )
327 {
328 CvPoint *pt_left, *pt_right, *pt_top, *pt_bottom;
329 CvPoint pt;
330 pt_left = pt_right = pt_top = pt_bottom = (CvPoint *)(reader.ptr);
331 CV_READ_SEQ_ELEM( pt, reader );
332
333 for( i = 1; i < count; i++ )
334 {
335 CvPoint* pt_ptr = (CvPoint*)reader.ptr;
336 CV_READ_SEQ_ELEM( pt, reader );
337
338 if( pt.x < pt_left->x )
339 pt_left = pt_ptr;
340 if( pt.x > pt_right->x )
341 pt_right = pt_ptr;
342 if( pt.y < pt_top->y )
343 pt_top = pt_ptr;
344 if( pt.y > pt_bottom->y )
345 pt_bottom = pt_ptr;
346 }
347
348 pts[0] = cvPointTo32f( *pt_left );
349 pts[1] = cvPointTo32f( *pt_right );
350 pts[2] = cvPointTo32f( *pt_top );
351 pts[3] = cvPointTo32f( *pt_bottom );
352 }
353 else
354 {
355 CvPoint2D32f *pt_left, *pt_right, *pt_top, *pt_bottom;
356 CvPoint2D32f pt;
357 pt_left = pt_right = pt_top = pt_bottom = (CvPoint2D32f *) (reader.ptr);
358 CV_READ_SEQ_ELEM( pt, reader );
359
360 for( i = 1; i < count; i++ )
361 {
362 CvPoint2D32f* pt_ptr = (CvPoint2D32f*)reader.ptr;
363 CV_READ_SEQ_ELEM( pt, reader );
364
365 if( pt.x < pt_left->x )
366 pt_left = pt_ptr;
367 if( pt.x > pt_right->x )
368 pt_right = pt_ptr;
369 if( pt.y < pt_top->y )
370 pt_top = pt_ptr;
371 if( pt.y > pt_bottom->y )
372 pt_bottom = pt_ptr;
373 }
374
375 pts[0] = *pt_left;
376 pts[1] = *pt_right;
377 pts[2] = *pt_top;
378 pts[3] = *pt_bottom;
379 }
380
381 for( k = 0; k < max_iters; k++ )
382 {
383 double min_delta = 0, delta;
384 CvPoint2D32f ptfl;
385
386 icvFindEnslosingCicle4pts_32f( pts, ¢er, &radius );
387 cvStartReadSeq( sequence, &reader, 0 );
388
389 for( i = 0; i < count; i++ )
390 {
391 if( !is_float )
392 {
393 ptfl.x = (float)((CvPoint*)reader.ptr)->x;
394 ptfl.y = (float)((CvPoint*)reader.ptr)->y;
395 }
396 else
397 {
398 ptfl = *(CvPoint2D32f*)reader.ptr;
399 }
400 CV_NEXT_SEQ_ELEM( sequence->elem_size, reader );
401
402 delta = icvIsPtInCircle( ptfl, center, radius );
403 if( delta < min_delta )
404 {
405 min_delta = delta;
406 pts[3] = ptfl;
407 }
408 }
409 result = min_delta >= 0;
410 if( result )
411 break;
412 }
413
414 if( !result )
415 {
416 cvStartReadSeq( sequence, &reader, 0 );
417 radius = 0.f;
418
419 for( i = 0; i < count; i++ )
420 {
421 CvPoint2D32f ptfl;
422 float t, dx, dy;
423
424 if( !is_float )
425 {
426 ptfl.x = (float)((CvPoint*)reader.ptr)->x;
427 ptfl.y = (float)((CvPoint*)reader.ptr)->y;
428 }
429 else
430 {
431 ptfl = *(CvPoint2D32f*)reader.ptr;
432 }
433
434 CV_NEXT_SEQ_ELEM( sequence->elem_size, reader );
435 dx = center.x - ptfl.x;
436 dy = center.y - ptfl.y;
437 t = dx*dx + dy*dy;
438 radius = MAX(radius,t);
439 }
440
441 radius = (float)(sqrt(radius)*(1 + eps));
442 result = 1;
443 }
444
445 __END__;
446
447 *_center = center;
448 *_radius = radius;
449
450 return result;
451 }
452
453
454 /* area of a whole sequence */
455 static CvStatus
icvContourArea(const CvSeq * contour,double * area)456 icvContourArea( const CvSeq* contour, double *area )
457 {
458 if( contour->total )
459 {
460 CvSeqReader reader;
461 int lpt = contour->total;
462 double a00 = 0, xi_1, yi_1;
463 int is_float = CV_SEQ_ELTYPE(contour) == CV_32FC2;
464
465 cvStartReadSeq( contour, &reader, 0 );
466
467 if( !is_float )
468 {
469 xi_1 = ((CvPoint*)(reader.ptr))->x;
470 yi_1 = ((CvPoint*)(reader.ptr))->y;
471 }
472 else
473 {
474 xi_1 = ((CvPoint2D32f*)(reader.ptr))->x;
475 yi_1 = ((CvPoint2D32f*)(reader.ptr))->y;
476 }
477 CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
478
479 while( lpt-- > 0 )
480 {
481 double dxy, xi, yi;
482
483 if( !is_float )
484 {
485 xi = ((CvPoint*)(reader.ptr))->x;
486 yi = ((CvPoint*)(reader.ptr))->y;
487 }
488 else
489 {
490 xi = ((CvPoint2D32f*)(reader.ptr))->x;
491 yi = ((CvPoint2D32f*)(reader.ptr))->y;
492 }
493 CV_NEXT_SEQ_ELEM( contour->elem_size, reader );
494
495 dxy = xi_1 * yi - xi * yi_1;
496 a00 += dxy;
497 xi_1 = xi;
498 yi_1 = yi;
499 }
500
501 *area = a00 * 0.5;
502 }
503 else
504 *area = 0;
505
506 return CV_OK;
507 }
508
509
510 /****************************************************************************************\
511
512 copy data from one buffer to other buffer
513
514 \****************************************************************************************/
515
516 static CvStatus
icvMemCopy(double ** buf1,double ** buf2,double ** buf3,int * b_max)517 icvMemCopy( double **buf1, double **buf2, double **buf3, int *b_max )
518 {
519 int bb;
520
521 if( (*buf1 == NULL && *buf2 == NULL) || *buf3 == NULL )
522 return CV_NULLPTR_ERR;
523
524 bb = *b_max;
525 if( *buf2 == NULL )
526 {
527 *b_max = 2 * (*b_max);
528 *buf2 = (double *)cvAlloc( (*b_max) * sizeof( double ));
529
530 if( *buf2 == NULL )
531 return CV_OUTOFMEM_ERR;
532
533 memcpy( *buf2, *buf3, bb * sizeof( double ));
534
535 *buf3 = *buf2;
536 cvFree( buf1 );
537 *buf1 = NULL;
538 }
539 else
540 {
541 *b_max = 2 * (*b_max);
542 *buf1 = (double *) cvAlloc( (*b_max) * sizeof( double ));
543
544 if( *buf1 == NULL )
545 return CV_OUTOFMEM_ERR;
546
547 memcpy( *buf1, *buf3, bb * sizeof( double ));
548
549 *buf3 = *buf1;
550 cvFree( buf2 );
551 *buf2 = NULL;
552 }
553 return CV_OK;
554 }
555
556
557 /* area of a contour sector */
icvContourSecArea(CvSeq * contour,CvSlice slice,double * area)558 static CvStatus icvContourSecArea( CvSeq * contour, CvSlice slice, double *area )
559 {
560 CvPoint pt; /* pointer to points */
561 CvPoint pt_s, pt_e; /* first and last points */
562 CvSeqReader reader; /* points reader of contour */
563
564 int p_max = 2, p_ind;
565 int lpt, flag, i;
566 double a00; /* unnormalized moments m00 */
567 double xi, yi, xi_1, yi_1, x0, y0, dxy, sk, sk1, t;
568 double x_s, y_s, nx, ny, dx, dy, du, dv;
569 double eps = 1.e-5;
570 double *p_are1, *p_are2, *p_are;
571
572 assert( contour != NULL );
573
574 if( contour == NULL )
575 return CV_NULLPTR_ERR;
576
577 if( !CV_IS_SEQ_POLYGON( contour ))
578 return CV_BADFLAG_ERR;
579
580 lpt = cvSliceLength( slice, contour );
581 /*if( n2 >= n1 )
582 lpt = n2 - n1 + 1;
583 else
584 lpt = contour->total - n1 + n2 + 1;*/
585
586 if( contour->total && lpt > 2 )
587 {
588 a00 = x0 = y0 = xi_1 = yi_1 = 0;
589 sk1 = 0;
590 flag = 0;
591 dxy = 0;
592 p_are1 = (double *) cvAlloc( p_max * sizeof( double ));
593
594 if( p_are1 == NULL )
595 return CV_OUTOFMEM_ERR;
596
597 p_are = p_are1;
598 p_are2 = NULL;
599
600 cvStartReadSeq( contour, &reader, 0 );
601 cvSetSeqReaderPos( &reader, slice.start_index );
602 CV_READ_SEQ_ELEM( pt_s, reader );
603 p_ind = 0;
604 cvSetSeqReaderPos( &reader, slice.end_index );
605 CV_READ_SEQ_ELEM( pt_e, reader );
606
607 /* normal coefficients */
608 nx = pt_s.y - pt_e.y;
609 ny = pt_e.x - pt_s.x;
610 cvSetSeqReaderPos( &reader, slice.start_index );
611
612 while( lpt-- > 0 )
613 {
614 CV_READ_SEQ_ELEM( pt, reader );
615
616 if( flag == 0 )
617 {
618 xi_1 = (double) pt.x;
619 yi_1 = (double) pt.y;
620 x0 = xi_1;
621 y0 = yi_1;
622 sk1 = 0;
623 flag = 1;
624 }
625 else
626 {
627 xi = (double) pt.x;
628 yi = (double) pt.y;
629
630 /**************** edges intersection examination **************************/
631 sk = nx * (xi - pt_s.x) + ny * (yi - pt_s.y);
632 if( (fabs( sk ) < eps && lpt > 0) || sk * sk1 < -eps )
633 {
634 if( fabs( sk ) < eps )
635 {
636 dxy = xi_1 * yi - xi * yi_1;
637 a00 = a00 + dxy;
638 dxy = xi * y0 - x0 * yi;
639 a00 = a00 + dxy;
640
641 if( p_ind >= p_max )
642 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
643
644 p_are[p_ind] = a00 / 2.;
645 p_ind++;
646 a00 = 0;
647 sk1 = 0;
648 x0 = xi;
649 y0 = yi;
650 dxy = 0;
651 }
652 else
653 {
654 /* define intersection point */
655 dv = yi - yi_1;
656 du = xi - xi_1;
657 dx = ny;
658 dy = -nx;
659 if( fabs( du ) > eps )
660 t = ((yi_1 - pt_s.y) * du + dv * (pt_s.x - xi_1)) /
661 (du * dy - dx * dv);
662 else
663 t = (xi_1 - pt_s.x) / dx;
664 if( t > eps && t < 1 - eps )
665 {
666 x_s = pt_s.x + t * dx;
667 y_s = pt_s.y + t * dy;
668 dxy = xi_1 * y_s - x_s * yi_1;
669 a00 += dxy;
670 dxy = x_s * y0 - x0 * y_s;
671 a00 += dxy;
672 if( p_ind >= p_max )
673 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
674
675 p_are[p_ind] = a00 / 2.;
676 p_ind++;
677
678 a00 = 0;
679 sk1 = 0;
680 x0 = x_s;
681 y0 = y_s;
682 dxy = x_s * yi - xi * y_s;
683 }
684 }
685 }
686 else
687 dxy = xi_1 * yi - xi * yi_1;
688
689 a00 += dxy;
690 xi_1 = xi;
691 yi_1 = yi;
692 sk1 = sk;
693
694 }
695 }
696
697 xi = x0;
698 yi = y0;
699 dxy = xi_1 * yi - xi * yi_1;
700
701 a00 += dxy;
702
703 if( p_ind >= p_max )
704 icvMemCopy( &p_are1, &p_are2, &p_are, &p_max );
705
706 p_are[p_ind] = a00 / 2.;
707 p_ind++;
708
709 /* common area calculation */
710 *area = 0;
711 for( i = 0; i < p_ind; i++ )
712 (*area) += fabs( p_are[i] );
713
714 if( p_are1 != NULL )
715 cvFree( &p_are1 );
716 else if( p_are2 != NULL )
717 cvFree( &p_are2 );
718
719 return CV_OK;
720 }
721 else
722 return CV_BADSIZE_ERR;
723 }
724
725
726 /* external contour area function */
727 CV_IMPL double
cvContourArea(const void * array,CvSlice slice)728 cvContourArea( const void *array, CvSlice slice )
729 {
730 double area = 0;
731
732 CV_FUNCNAME( "cvContourArea" );
733
734 __BEGIN__;
735
736 CvContour contour_header;
737 CvSeq* contour = 0;
738 CvSeqBlock block;
739
740 if( CV_IS_SEQ( array ))
741 {
742 contour = (CvSeq*)array;
743 if( !CV_IS_SEQ_POLYLINE( contour ))
744 CV_ERROR( CV_StsBadArg, "Unsupported sequence type" );
745 }
746 else
747 {
748 CV_CALL( contour = cvPointSeqFromMat(
749 CV_SEQ_KIND_CURVE, array, &contour_header, &block ));
750 }
751
752 if( cvSliceLength( slice, contour ) == contour->total )
753 {
754 IPPI_CALL( icvContourArea( contour, &area ));
755 }
756 else
757 {
758 if( CV_SEQ_ELTYPE( contour ) != CV_32SC2 )
759 CV_ERROR( CV_StsUnsupportedFormat,
760 "Only curves with integer coordinates are supported in case of contour slice" );
761 IPPI_CALL( icvContourSecArea( contour, slice, &area ));
762 }
763
764 __END__;
765
766 return area;
767 }
768
769
770 /* for now this function works bad with singular cases
771 You can see in the code, that when some troubles with
772 matrices or some variables occur -
773 box filled with zero values is returned.
774 However in general function works fine.
775 */
776 static void
icvFitEllipse_F(CvSeq * points,CvBox2D * box)777 icvFitEllipse_F( CvSeq* points, CvBox2D* box )
778 {
779 CvMat* D = 0;
780
781 CV_FUNCNAME( "icvFitEllipse_F" );
782
783 __BEGIN__;
784
785 double S[36], C[36], T[36];
786
787 int i, j;
788 double eigenvalues[6], eigenvectors[36];
789 double a, b, c, d, e, f;
790 double x0, y0, idet, scale, offx = 0, offy = 0;
791
792 int n = points->total;
793 CvSeqReader reader;
794 int is_float = CV_SEQ_ELTYPE(points) == CV_32FC2;
795
796 CvMat _S = cvMat(6,6,CV_64F,S), _C = cvMat(6,6,CV_64F,C), _T = cvMat(6,6,CV_64F,T);
797 CvMat _EIGVECS = cvMat(6,6,CV_64F,eigenvectors), _EIGVALS = cvMat(6,1,CV_64F,eigenvalues);
798
799 /* create matrix D of input points */
800 CV_CALL( D = cvCreateMat( n, 6, CV_64F ));
801
802 cvStartReadSeq( points, &reader );
803
804 /* shift all points to zero */
805 for( i = 0; i < n; i++ )
806 {
807 if( !is_float )
808 {
809 offx += ((CvPoint*)reader.ptr)->x;
810 offy += ((CvPoint*)reader.ptr)->y;
811 }
812 else
813 {
814 offx += ((CvPoint2D32f*)reader.ptr)->x;
815 offy += ((CvPoint2D32f*)reader.ptr)->y;
816 }
817 CV_NEXT_SEQ_ELEM( points->elem_size, reader );
818 }
819
820 offx /= n;
821 offy /= n;
822
823 // fill matrix rows as (x*x, x*y, y*y, x, y, 1 )
824 for( i = 0; i < n; i++ )
825 {
826 double x, y;
827 double* Dptr = D->data.db + i*6;
828
829 if( !is_float )
830 {
831 x = ((CvPoint*)reader.ptr)->x - offx;
832 y = ((CvPoint*)reader.ptr)->y - offy;
833 }
834 else
835 {
836 x = ((CvPoint2D32f*)reader.ptr)->x - offx;
837 y = ((CvPoint2D32f*)reader.ptr)->y - offy;
838 }
839 CV_NEXT_SEQ_ELEM( points->elem_size, reader );
840
841 Dptr[0] = x * x;
842 Dptr[1] = x * y;
843 Dptr[2] = y * y;
844 Dptr[3] = x;
845 Dptr[4] = y;
846 Dptr[5] = 1.;
847 }
848
849 // S = D^t*D
850 cvMulTransposed( D, &_S, 1 );
851 cvSVD( &_S, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T );
852
853 for( i = 0; i < 6; i++ )
854 {
855 double a = eigenvalues[i];
856 a = a < DBL_EPSILON ? 0 : 1./sqrt(sqrt(a));
857 for( j = 0; j < 6; j++ )
858 eigenvectors[i*6 + j] *= a;
859 }
860
861 // C = Q^-1 = transp(INVEIGV) * INVEIGV
862 cvMulTransposed( &_EIGVECS, &_C, 1 );
863
864 cvZero( &_S );
865 S[2] = 2.;
866 S[7] = -1.;
867 S[12] = 2.;
868
869 // S = Q^-1*S*Q^-1
870 cvMatMul( &_C, &_S, &_T );
871 cvMatMul( &_T, &_C, &_S );
872
873 // and find its eigenvalues and vectors too
874 //cvSVD( &_S, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T );
875 cvEigenVV( &_S, &_EIGVECS, &_EIGVALS, 0 );
876
877 for( i = 0; i < 3; i++ )
878 if( eigenvalues[i] > 0 )
879 break;
880
881 if( i >= 3 /*eigenvalues[0] < DBL_EPSILON*/ )
882 {
883 box->center.x = box->center.y =
884 box->size.width = box->size.height =
885 box->angle = 0.f;
886 EXIT;
887 }
888
889 // now find truthful eigenvector
890 _EIGVECS = cvMat( 6, 1, CV_64F, eigenvectors + 6*i );
891 _T = cvMat( 6, 1, CV_64F, T );
892 // Q^-1*eigenvecs[0]
893 cvMatMul( &_C, &_EIGVECS, &_T );
894
895 // extract vector components
896 a = T[0]; b = T[1]; c = T[2]; d = T[3]; e = T[4]; f = T[5];
897
898 ///////////////// extract ellipse axes from above values ////////////////
899
900 /*
901 1) find center of ellipse
902 it satisfy equation
903 | a b/2 | * | x0 | + | d/2 | = |0 |
904 | b/2 c | | y0 | | e/2 | |0 |
905
906 */
907 idet = a * c - b * b * 0.25;
908 idet = idet > DBL_EPSILON ? 1./idet : 0;
909
910 // we must normalize (a b c d e f ) to fit (4ac-b^2=1)
911 scale = sqrt( 0.25 * idet );
912
913 if( scale < DBL_EPSILON )
914 {
915 box->center.x = (float)offx;
916 box->center.y = (float)offy;
917 box->size.width = box->size.height = box->angle = 0.f;
918 EXIT;
919 }
920
921 a *= scale;
922 b *= scale;
923 c *= scale;
924 d *= scale;
925 e *= scale;
926 f *= scale;
927
928 x0 = (-d * c + e * b * 0.5) * 2.;
929 y0 = (-a * e + d * b * 0.5) * 2.;
930
931 // recover center
932 box->center.x = (float)(x0 + offx);
933 box->center.y = (float)(y0 + offy);
934
935 // offset ellipse to (x0,y0)
936 // new f == F(x0,y0)
937 f += a * x0 * x0 + b * x0 * y0 + c * y0 * y0 + d * x0 + e * y0;
938
939 if( fabs(f) < DBL_EPSILON )
940 {
941 box->size.width = box->size.height = box->angle = 0.f;
942 EXIT;
943 }
944
945 scale = -1. / f;
946 // normalize to f = 1
947 a *= scale;
948 b *= scale;
949 c *= scale;
950
951 // extract axis of ellipse
952 // one more eigenvalue operation
953 S[0] = a;
954 S[1] = S[2] = b * 0.5;
955 S[3] = c;
956
957 _S = cvMat( 2, 2, CV_64F, S );
958 _EIGVECS = cvMat( 2, 2, CV_64F, eigenvectors );
959 _EIGVALS = cvMat( 1, 2, CV_64F, eigenvalues );
960 cvSVD( &_S, &_EIGVALS, &_EIGVECS, 0, CV_SVD_MODIFY_A + CV_SVD_U_T );
961
962 // exteract axis length from eigenvectors
963 box->size.width = (float)(2./sqrt(eigenvalues[0]));
964 box->size.height = (float)(2./sqrt(eigenvalues[1]));
965
966 // calc angle
967 box->angle = (float)(180 - atan2(eigenvectors[2], eigenvectors[3])*180/CV_PI);
968
969 __END__;
970
971 cvReleaseMat( &D );
972 }
973
974
975 CV_IMPL CvBox2D
cvFitEllipse2(const CvArr * array)976 cvFitEllipse2( const CvArr* array )
977 {
978 CvBox2D box;
979 double* Ad = 0, *bd = 0;
980
981 CV_FUNCNAME( "cvFitEllipse2" );
982
983 memset( &box, 0, sizeof(box));
984
985 __BEGIN__;
986
987 CvContour contour_header;
988 CvSeq* ptseq = 0;
989 CvSeqBlock block;
990 int n;
991
992 if( CV_IS_SEQ( array ))
993 {
994 ptseq = (CvSeq*)array;
995 if( !CV_IS_SEQ_POINT_SET( ptseq ))
996 CV_ERROR( CV_StsBadArg, "Unsupported sequence type" );
997 }
998 else
999 {
1000 CV_CALL( ptseq = cvPointSeqFromMat(
1001 CV_SEQ_KIND_GENERIC, array, &contour_header, &block ));
1002 }
1003
1004 n = ptseq->total;
1005 if( n < 5 )
1006 CV_ERROR( CV_StsBadSize, "Number of points should be >= 6" );
1007 #if 1
1008 icvFitEllipse_F( ptseq, &box );
1009 #else
1010 /*
1011 * New fitellipse algorithm, contributed by Dr. Daniel Weiss
1012 */
1013 {
1014 double gfp[5], rp[5], t;
1015 CvMat A, b, x;
1016 const double min_eps = 1e-6;
1017 int i, is_float;
1018 CvSeqReader reader;
1019
1020 CV_CALL( Ad = (double*)cvAlloc( n*5*sizeof(Ad[0]) ));
1021 CV_CALL( bd = (double*)cvAlloc( n*sizeof(bd[0]) ));
1022
1023 // first fit for parameters A - E
1024 A = cvMat( n, 5, CV_64F, Ad );
1025 b = cvMat( n, 1, CV_64F, bd );
1026 x = cvMat( 5, 1, CV_64F, gfp );
1027
1028 cvStartReadSeq( ptseq, &reader );
1029 is_float = CV_SEQ_ELTYPE(ptseq) == CV_32FC2;
1030
1031 for( i = 0; i < n; i++ )
1032 {
1033 CvPoint2D32f p;
1034 if( is_float )
1035 p = *(CvPoint2D32f*)(reader.ptr);
1036 else
1037 {
1038 p.x = (float)((int*)reader.ptr)[0];
1039 p.y = (float)((int*)reader.ptr)[1];
1040 }
1041 CV_NEXT_SEQ_ELEM( sizeof(p), reader );
1042
1043 bd[i] = 10000.0; // 1.0?
1044 Ad[i*5] = -(double)p.x * p.x; // A - C signs inverted as proposed by APP
1045 Ad[i*5 + 1] = -(double)p.y * p.y;
1046 Ad[i*5 + 2] = -(double)p.x * p.y;
1047 Ad[i*5 + 3] = p.x;
1048 Ad[i*5 + 4] = p.y;
1049 }
1050
1051 cvSolve( &A, &b, &x, CV_SVD );
1052
1053 // now use general-form parameters A - E to find the ellipse center:
1054 // differentiate general form wrt x/y to get two equations for cx and cy
1055 A = cvMat( 2, 2, CV_64F, Ad );
1056 b = cvMat( 2, 1, CV_64F, bd );
1057 x = cvMat( 2, 1, CV_64F, rp );
1058 Ad[0] = 2 * gfp[0];
1059 Ad[1] = Ad[2] = gfp[2];
1060 Ad[3] = 2 * gfp[1];
1061 bd[0] = gfp[3];
1062 bd[1] = gfp[4];
1063 cvSolve( &A, &b, &x, CV_SVD );
1064
1065 // re-fit for parameters A - C with those center coordinates
1066 A = cvMat( n, 3, CV_64F, Ad );
1067 b = cvMat( n, 1, CV_64F, bd );
1068 x = cvMat( 3, 1, CV_64F, gfp );
1069 for( i = 0; i < n; i++ )
1070 {
1071 CvPoint2D32f p;
1072 if( is_float )
1073 p = *(CvPoint2D32f*)(reader.ptr);
1074 else
1075 {
1076 p.x = (float)((int*)reader.ptr)[0];
1077 p.y = (float)((int*)reader.ptr)[1];
1078 }
1079 CV_NEXT_SEQ_ELEM( sizeof(p), reader );
1080 bd[i] = 1.0;
1081 Ad[i * 3] = (p.x - rp[0]) * (p.x - rp[0]);
1082 Ad[i * 3 + 1] = (p.y - rp[1]) * (p.y - rp[1]);
1083 Ad[i * 3 + 2] = (p.x - rp[0]) * (p.y - rp[1]);
1084 }
1085 cvSolve(&A, &b, &x, CV_SVD);
1086
1087 // store angle and radii
1088 rp[4] = -0.5 * atan2(gfp[2], gfp[1] - gfp[0]); // convert from APP angle usage
1089 t = sin(-2.0 * rp[4]);
1090 if( fabs(t) > fabs(gfp[2])*min_eps )
1091 t = gfp[2]/t;
1092 else
1093 t = gfp[1] - gfp[0];
1094 rp[2] = fabs(gfp[0] + gfp[1] - t);
1095 if( rp[2] > min_eps )
1096 rp[2] = sqrt(2.0 / rp[2]);
1097 rp[3] = fabs(gfp[0] + gfp[1] + t);
1098 if( rp[3] > min_eps )
1099 rp[3] = sqrt(2.0 / rp[3]);
1100
1101 box.center.x = (float)rp[0];
1102 box.center.y = (float)rp[1];
1103 box.size.width = (float)(rp[2]*2);
1104 box.size.height = (float)(rp[3]*2);
1105 if( box.size.width > box.size.height )
1106 {
1107 float tmp;
1108 CV_SWAP( box.size.width, box.size.height, tmp );
1109 box.angle = (float)(90 + rp[4]*180/CV_PI);
1110 }
1111 if( box.angle < -180 )
1112 box.angle += 360;
1113 if( box.angle > 360 )
1114 box.angle -= 360;
1115 }
1116 #endif
1117 __END__;
1118
1119 cvFree( &Ad );
1120 cvFree( &bd );
1121
1122 return box;
1123 }
1124
1125
1126 /* Calculates bounding rectagnle of a point set or retrieves already calculated */
1127 CV_IMPL CvRect
cvBoundingRect(CvArr * array,int update)1128 cvBoundingRect( CvArr* array, int update )
1129 {
1130 CvSeqReader reader;
1131 CvRect rect = { 0, 0, 0, 0 };
1132 CvContour contour_header;
1133 CvSeq* ptseq = 0;
1134 CvSeqBlock block;
1135
1136 CV_FUNCNAME( "cvBoundingRect" );
1137
1138 __BEGIN__;
1139
1140 CvMat stub, *mat = 0;
1141 int xmin = 0, ymin = 0, xmax = -1, ymax = -1, i, j, k;
1142 int calculate = update;
1143
1144 if( CV_IS_SEQ( array ))
1145 {
1146 ptseq = (CvSeq*)array;
1147 if( !CV_IS_SEQ_POINT_SET( ptseq ))
1148 CV_ERROR( CV_StsBadArg, "Unsupported sequence type" );
1149
1150 if( ptseq->header_size < (int)sizeof(CvContour))
1151 {
1152 /*if( update == 1 )
1153 CV_ERROR( CV_StsBadArg, "The header is too small to fit the rectangle, "
1154 "so it could not be updated" );*/
1155 update = 0;
1156 calculate = 1;
1157 }
1158 }
1159 else
1160 {
1161 CV_CALL( mat = cvGetMat( array, &stub ));
1162 if( CV_MAT_TYPE(mat->type) == CV_32SC2 ||
1163 CV_MAT_TYPE(mat->type) == CV_32FC2 )
1164 {
1165 CV_CALL( ptseq = cvPointSeqFromMat(
1166 CV_SEQ_KIND_GENERIC, mat, &contour_header, &block ));
1167 mat = 0;
1168 }
1169 else if( CV_MAT_TYPE(mat->type) != CV_8UC1 &&
1170 CV_MAT_TYPE(mat->type) != CV_8SC1 )
1171 CV_ERROR( CV_StsUnsupportedFormat,
1172 "The image/matrix format is not supported by the function" );
1173 update = 0;
1174 calculate = 1;
1175 }
1176
1177 if( !calculate )
1178 {
1179 rect = ((CvContour*)ptseq)->rect;
1180 EXIT;
1181 }
1182
1183 if( mat )
1184 {
1185 CvSize size = cvGetMatSize(mat);
1186 xmin = size.width;
1187 ymin = -1;
1188
1189 for( i = 0; i < size.height; i++ )
1190 {
1191 uchar* _ptr = mat->data.ptr + i*mat->step;
1192 uchar* ptr = (uchar*)cvAlignPtr(_ptr, 4);
1193 int have_nz = 0, k_min, offset = (int)(ptr - _ptr);
1194 j = 0;
1195 offset = MIN(offset, size.width);
1196 for( ; j < offset; j++ )
1197 if( _ptr[j] )
1198 {
1199 have_nz = 1;
1200 break;
1201 }
1202 if( j < offset )
1203 {
1204 if( j < xmin )
1205 xmin = j;
1206 if( j > xmax )
1207 xmax = j;
1208 }
1209 if( offset < size.width )
1210 {
1211 xmin -= offset;
1212 xmax -= offset;
1213 size.width -= offset;
1214 j = 0;
1215 for( ; j <= xmin - 4; j += 4 )
1216 if( *((int*)(ptr+j)) )
1217 break;
1218 for( ; j < xmin; j++ )
1219 if( ptr[j] )
1220 {
1221 xmin = j;
1222 if( j > xmax )
1223 xmax = j;
1224 have_nz = 1;
1225 break;
1226 }
1227 k_min = MAX(j-1, xmax);
1228 k = size.width - 1;
1229 for( ; k > k_min && (k&3) != 3; k-- )
1230 if( ptr[k] )
1231 break;
1232 if( k > k_min && (k&3) == 3 )
1233 {
1234 for( ; k > k_min+3; k -= 4 )
1235 if( *((int*)(ptr+k-3)) )
1236 break;
1237 }
1238 for( ; k > k_min; k-- )
1239 if( ptr[k] )
1240 {
1241 xmax = k;
1242 have_nz = 1;
1243 break;
1244 }
1245 if( !have_nz )
1246 {
1247 j &= ~3;
1248 for( ; j <= k - 3; j += 4 )
1249 if( *((int*)(ptr+j)) )
1250 break;
1251 for( ; j <= k; j++ )
1252 if( ptr[j] )
1253 {
1254 have_nz = 1;
1255 break;
1256 }
1257 }
1258 xmin += offset;
1259 xmax += offset;
1260 size.width += offset;
1261 }
1262 if( have_nz )
1263 {
1264 if( ymin < 0 )
1265 ymin = i;
1266 ymax = i;
1267 }
1268 }
1269
1270 if( xmin >= size.width )
1271 xmin = ymin = 0;
1272 }
1273 else if( ptseq->total )
1274 {
1275 int is_float = CV_SEQ_ELTYPE(ptseq) == CV_32FC2;
1276 cvStartReadSeq( ptseq, &reader, 0 );
1277
1278 if( !is_float )
1279 {
1280 CvPoint pt;
1281 /* init values */
1282 CV_READ_SEQ_ELEM( pt, reader );
1283 xmin = xmax = pt.x;
1284 ymin = ymax = pt.y;
1285
1286 for( i = 1; i < ptseq->total; i++ )
1287 {
1288 CV_READ_SEQ_ELEM( pt, reader );
1289
1290 if( xmin > pt.x )
1291 xmin = pt.x;
1292
1293 if( xmax < pt.x )
1294 xmax = pt.x;
1295
1296 if( ymin > pt.y )
1297 ymin = pt.y;
1298
1299 if( ymax < pt.y )
1300 ymax = pt.y;
1301 }
1302 }
1303 else
1304 {
1305 CvPoint pt;
1306 Cv32suf v;
1307 /* init values */
1308 CV_READ_SEQ_ELEM( pt, reader );
1309 xmin = xmax = CV_TOGGLE_FLT(pt.x);
1310 ymin = ymax = CV_TOGGLE_FLT(pt.y);
1311
1312 for( i = 1; i < ptseq->total; i++ )
1313 {
1314 CV_READ_SEQ_ELEM( pt, reader );
1315 pt.x = CV_TOGGLE_FLT(pt.x);
1316 pt.y = CV_TOGGLE_FLT(pt.y);
1317
1318 if( xmin > pt.x )
1319 xmin = pt.x;
1320
1321 if( xmax < pt.x )
1322 xmax = pt.x;
1323
1324 if( ymin > pt.y )
1325 ymin = pt.y;
1326
1327 if( ymax < pt.y )
1328 ymax = pt.y;
1329 }
1330
1331 v.i = CV_TOGGLE_FLT(xmin); xmin = cvFloor(v.f);
1332 v.i = CV_TOGGLE_FLT(ymin); ymin = cvFloor(v.f);
1333 /* because right and bottom sides of
1334 the bounding rectangle are not inclusive
1335 (note +1 in width and height calculation below),
1336 cvFloor is used here instead of cvCeil */
1337 v.i = CV_TOGGLE_FLT(xmax); xmax = cvFloor(v.f);
1338 v.i = CV_TOGGLE_FLT(ymax); ymax = cvFloor(v.f);
1339 }
1340 }
1341
1342 rect.x = xmin;
1343 rect.y = ymin;
1344 rect.width = xmax - xmin + 1;
1345 rect.height = ymax - ymin + 1;
1346
1347 if( update )
1348 ((CvContour*)ptseq)->rect = rect;
1349
1350 __END__;
1351
1352 return rect;
1353 }
1354
1355
1356 /* End of file. */
1357