--- /dev/null
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+// By downloading, copying, installing or using the software you agree to this license.
+// If you do not agree to this license, do not download, install,
+// copy or use the software.
+//
+//
+// Intel License Agreement
+// For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+// * Redistribution's of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+//
+// * Redistribution's in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+//
+// * The name of Intel Corporation may not be used to endorse or promote products
+// derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+#include "_cv.h"
+
+typedef struct
+{
+ int bottom;
+ int left;
+ float height;
+ float width;
+ float base_a;
+ float base_b;
+}
+icvMinAreaState;
+
+#define CV_CALIPERS_MAXHEIGHT 0
+#define CV_CALIPERS_MINAREARECT 1
+#define CV_CALIPERS_MAXDIST 2
+
+/*F///////////////////////////////////////////////////////////////////////////////////////
+// Name: icvRotatingCalipers
+// Purpose:
+// Rotating calipers algorithm with some applications
+//
+// Context:
+// Parameters:
+// points - convex hull vertices ( any orientation )
+// n - number of vertices
+// mode - concrete application of algorithm
+// can be CV_CALIPERS_MAXDIST or
+// CV_CALIPERS_MINAREARECT
+// left, bottom, right, top - indexes of extremal points
+// out - output info.
+// In case CV_CALIPERS_MAXDIST it points to float value -
+// maximal height of polygon.
+// In case CV_CALIPERS_MINAREARECT
+// ((CvPoint2D32f*)out)[0] - corner
+// ((CvPoint2D32f*)out)[1] - vector1
+// ((CvPoint2D32f*)out)[0] - corner2
+//
+// ^
+// |
+// vector2 |
+// |
+// |____________\
+// corner /
+// vector1
+//
+// Returns:
+// Notes:
+//F*/
+
+/* we will use usual cartesian coordinates */
+static void
+icvRotatingCalipers( CvPoint2D32f* points, int n, int mode, float* out )
+{
+ float minarea = FLT_MAX;
+ float max_dist = 0;
+ char buffer[32];
+ int i, k;
+ CvPoint2D32f* vect = (CvPoint2D32f*)cvAlloc( n * sizeof(vect[0]) );
+ float* inv_vect_length = (float*)cvAlloc( n * sizeof(inv_vect_length[0]) );
+ int left = 0, bottom = 0, right = 0, top = 0;
+ int seq[4] = { -1, -1, -1, -1 };
+
+ /* rotating calipers sides will always have coordinates
+ (a,b) (-b,a) (-a,-b) (b, -a)
+ */
+ /* this is a first base bector (a,b) initialized by (1,0) */
+ float orientation = 0;
+ float base_a;
+ float base_b = 0;
+
+ float left_x, right_x, top_y, bottom_y;
+ CvPoint2D32f pt0 = points[0];
+
+ left_x = right_x = pt0.x;
+ top_y = bottom_y = pt0.y;
+
+ for( i = 0; i < n; i++ )
+ {
+ double dx, dy;
+
+ if( pt0.x < left_x )
+ left_x = pt0.x, left = i;
+
+ if( pt0.x > right_x )
+ right_x = pt0.x, right = i;
+
+ if( pt0.y > top_y )
+ top_y = pt0.y, top = i;
+
+ if( pt0.y < bottom_y )
+ bottom_y = pt0.y, bottom = i;
+
+ CvPoint2D32f pt = points[(i+1) & (i+1 < n ? -1 : 0)];
+
+ dx = pt.x - pt0.x;
+ dy = pt.y - pt0.y;
+
+ vect[i].x = (float)dx;
+ vect[i].y = (float)dy;
+ inv_vect_length[i] = (float)(1./sqrt(dx*dx + dy*dy));
+
+ pt0 = pt;
+ }
+
+ //cvbInvSqrt( inv_vect_length, inv_vect_length, n );
+
+ /* find convex hull orientation */
+ {
+ double ax = vect[n-1].x;
+ double ay = vect[n-1].y;
+
+ for( i = 0; i < n; i++ )
+ {
+ double bx = vect[i].x;
+ double by = vect[i].y;
+
+ double convexity = ax * by - ay * bx;
+
+ if( convexity != 0 )
+ {
+ orientation = (convexity > 0) ? 1.f : (-1.f);
+ break;
+ }
+ ax = bx;
+ ay = by;
+ }
+ assert( orientation != 0 );
+ }
+ base_a = orientation;
+
+/*****************************************************************************************/
+/* init calipers position */
+ seq[0] = bottom;
+ seq[1] = right;
+ seq[2] = top;
+ seq[3] = left;
+/*****************************************************************************************/
+/* Main loop - evaluate angles and rotate calipers */
+
+ /* all of edges will be checked while rotating calipers by 90 degrees */
+ for( k = 0; k < n; k++ )
+ {
+ /* sinus of minimal angle */
+ /*float sinus;*/
+
+ /* compute cosine of angle between calipers side and polygon edge */
+ /* dp - dot product */
+ float dp0 = base_a * vect[seq[0]].x + base_b * vect[seq[0]].y;
+ float dp1 = -base_b * vect[seq[1]].x + base_a * vect[seq[1]].y;
+ float dp2 = -base_a * vect[seq[2]].x - base_b * vect[seq[2]].y;
+ float dp3 = base_b * vect[seq[3]].x - base_a * vect[seq[3]].y;
+
+ float cosalpha = dp0 * inv_vect_length[seq[0]];
+ float maxcos = cosalpha;
+
+ /* number of calipers edges, that has minimal angle with edge */
+ int main_element = 0;
+
+ /* choose minimal angle */
+ cosalpha = dp1 * inv_vect_length[seq[1]];
+ maxcos = (cosalpha > maxcos) ? (main_element = 1, cosalpha) : maxcos;
+ cosalpha = dp2 * inv_vect_length[seq[2]];
+ maxcos = (cosalpha > maxcos) ? (main_element = 2, cosalpha) : maxcos;
+ cosalpha = dp3 * inv_vect_length[seq[3]];
+ maxcos = (cosalpha > maxcos) ? (main_element = 3, cosalpha) : maxcos;
+
+ /*rotate calipers*/
+ {
+ //get next base
+ int pindex = seq[main_element];
+ float lead_x = vect[pindex].x*inv_vect_length[pindex];
+ float lead_y = vect[pindex].y*inv_vect_length[pindex];
+ switch( main_element )
+ {
+ case 0:
+ base_a = lead_x;
+ base_b = lead_y;
+ break;
+ case 1:
+ base_a = lead_y;
+ base_b = -lead_x;
+ break;
+ case 2:
+ base_a = -lead_x;
+ base_b = -lead_y;
+ break;
+ case 3:
+ base_a = -lead_y;
+ base_b = lead_x;
+ break;
+ default: assert(0);
+ }
+ }
+ /* change base point of main edge */
+ seq[main_element] += 1;
+ seq[main_element] = (seq[main_element] == n) ? 0 : seq[main_element];
+
+
+ switch (mode)
+ {
+ case CV_CALIPERS_MAXHEIGHT:
+ {
+ /* now main element lies on edge alligned to calipers side */
+
+ /* find opposite element i.e. transform */
+ /* 0->2, 1->3, 2->0, 3->1 */
+ int opposite_el = main_element ^ 2;
+
+ float dx = points[seq[opposite_el]].x - points[seq[main_element]].x;
+ float dy = points[seq[opposite_el]].y - points[seq[main_element]].y;
+ float dist;
+
+ if( main_element & 1 )
+ dist = (float)fabs(dx * base_a + dy * base_b);
+ else
+ dist = (float)fabs(dx * (-base_b) + dy * base_a);
+
+ if( dist > max_dist )
+ max_dist = dist;
+
+ break;
+ }
+ case CV_CALIPERS_MINAREARECT:
+ /* find area of rectangle */
+ {
+ float height;
+ float area;
+
+ /* find vector left-right */
+ float dx = points[seq[1]].x - points[seq[3]].x;
+ float dy = points[seq[1]].y - points[seq[3]].y;
+
+ /* dotproduct */
+ float width = dx * base_a + dy * base_b;
+
+ /* find vector left-right */
+ dx = points[seq[2]].x - points[seq[0]].x;
+ dy = points[seq[2]].y - points[seq[0]].y;
+
+ /* dotproduct */
+ height = -dx * base_b + dy * base_a;
+
+ area = width * height;
+ if( area <= minarea )
+ {
+ float *buf = (float *) buffer;
+
+ minarea = area;
+ /* leftist point */
+ ((int *) buf)[0] = seq[3];
+ buf[1] = base_a;
+ buf[2] = width;
+ buf[3] = base_b;
+ buf[4] = height;
+ /* bottom point */
+ ((int *) buf)[5] = seq[0];
+ buf[6] = area;
+ }
+ break;
+ }
+ } /*switch */
+ } /* for */
+
+ switch (mode)
+ {
+ case CV_CALIPERS_MINAREARECT:
+ {
+ float *buf = (float *) buffer;
+
+ float A1 = buf[1];
+ float B1 = buf[3];
+
+ float A2 = -buf[3];
+ float B2 = buf[1];
+
+ float C1 = A1 * points[((int *) buf)[0]].x + points[((int *) buf)[0]].y * B1;
+ float C2 = A2 * points[((int *) buf)[5]].x + points[((int *) buf)[5]].y * B2;
+
+ float idet = 1.f / (A1 * B2 - A2 * B1);
+
+ float px = (C1 * B2 - C2 * B1) * idet;
+ float py = (A1 * C2 - A2 * C1) * idet;
+
+ out[0] = px;
+ out[1] = py;
+
+ out[2] = A1 * buf[2];
+ out[3] = B1 * buf[2];
+
+ out[4] = A2 * buf[4];
+ out[5] = B2 * buf[4];
+ }
+ break;
+ case CV_CALIPERS_MAXHEIGHT:
+ {
+ out[0] = max_dist;
+ }
+ break;
+ }
+
+ cvFree( &vect );
+ cvFree( &inv_vect_length );
+}
+
+
+CV_IMPL CvBox2D
+cvMinAreaRect2( const CvArr* array, CvMemStorage* storage )
+{
+ CvMemStorage* temp_storage = 0;
+ CvBox2D box;
+ CvPoint2D32f* points = 0;
+
+ CV_FUNCNAME( "cvMinAreaRect2" );
+
+ memset(&box, 0, sizeof(box));
+
+ __BEGIN__;
+
+ int i, n;
+ CvSeqReader reader;
+ CvContour contour_header;
+ CvSeqBlock block;
+ CvSeq* ptseq = (CvSeq*)array;
+ CvPoint2D32f out[3];
+
+ if( CV_IS_SEQ(ptseq) )
+ {
+ if( !CV_IS_SEQ_POINT_SET(ptseq) &&
+ (CV_SEQ_KIND(ptseq) != CV_SEQ_KIND_CURVE || !CV_IS_SEQ_CONVEX(ptseq) ||
+ CV_SEQ_ELTYPE(ptseq) != CV_SEQ_ELTYPE_PPOINT ))
+ CV_ERROR( CV_StsUnsupportedFormat,
+ "Input sequence must consist of 2d points or pointers to 2d points" );
+ if( !storage )
+ storage = ptseq->storage;
+ }
+ else
+ {
+ CV_CALL( ptseq = cvPointSeqFromMat(
+ CV_SEQ_KIND_GENERIC, array, &contour_header, &block ));
+ }
+
+ if( storage )
+ {
+ CV_CALL( temp_storage = cvCreateChildMemStorage( storage ));
+ }
+ else
+ {
+ CV_CALL( temp_storage = cvCreateMemStorage(1 << 10));
+ }
+
+ if( !CV_IS_SEQ_CONVEX( ptseq ))
+ {
+ CV_CALL( ptseq = cvConvexHull2( ptseq, temp_storage, CV_CLOCKWISE, 1 ));
+ }
+ else if( !CV_IS_SEQ_POINT_SET( ptseq ))
+ {
+ CvSeqWriter writer;
+
+ if( !CV_IS_SEQ(ptseq->v_prev) || !CV_IS_SEQ_POINT_SET(ptseq->v_prev))
+ CV_ERROR( CV_StsBadArg,
+ "Convex hull must have valid pointer to point sequence stored in v_prev" );
+ cvStartReadSeq( ptseq, &reader );
+ cvStartWriteSeq( CV_SEQ_KIND_CURVE|CV_SEQ_FLAG_CONVEX|CV_SEQ_ELTYPE(ptseq->v_prev),
+ sizeof(CvContour), CV_ELEM_SIZE(ptseq->v_prev->flags),
+ temp_storage, &writer );
+
+ for( i = 0; i < ptseq->total; i++ )
+ {
+ CvPoint pt = **(CvPoint**)(reader.ptr);
+ CV_WRITE_SEQ_ELEM( pt, writer );
+ }
+
+ ptseq = cvEndWriteSeq( &writer );
+ }
+
+ n = ptseq->total;
+
+ CV_CALL( points = (CvPoint2D32f*)cvAlloc( n*sizeof(points[0]) ));
+ cvStartReadSeq( ptseq, &reader );
+
+ if( CV_SEQ_ELTYPE( ptseq ) == CV_32SC2 )
+ {
+ for( i = 0; i < n; i++ )
+ {
+ CvPoint pt;
+ CV_READ_SEQ_ELEM( pt, reader );
+ points[i].x = (float)pt.x;
+ points[i].y = (float)pt.y;
+ }
+ }
+ else
+ {
+ for( i = 0; i < n; i++ )
+ {
+ CV_READ_SEQ_ELEM( points[i], reader );
+ }
+ }
+
+ if( n > 2 )
+ {
+ icvRotatingCalipers( points, n, CV_CALIPERS_MINAREARECT, (float*)out );
+ box.center.x = out[0].x + (out[1].x + out[2].x)*0.5f;
+ box.center.y = out[0].y + (out[1].y + out[2].y)*0.5f;
+ box.size.height = (float)sqrt((double)out[1].x*out[1].x + (double)out[1].y*out[1].y);
+ box.size.width = (float)sqrt((double)out[2].x*out[2].x + (double)out[2].y*out[2].y);
+ box.angle = (float)atan2( -(double)out[1].y, (double)out[1].x );
+ }
+ else if( n == 2 )
+ {
+ box.center.x = (points[0].x + points[1].x)*0.5f;
+ box.center.y = (points[0].y + points[1].y)*0.5f;
+ double dx = points[1].x - points[0].x;
+ double dy = points[1].y - points[0].y;
+ box.size.height = (float)sqrt(dx*dx + dy*dy);
+ box.size.width = 0;
+ box.angle = (float)atan2( -dy, dx );
+ }
+ else
+ {
+ if( n == 1 )
+ box.center = points[0];
+ }
+
+ box.angle = (float)(box.angle*180/CV_PI);
+
+ __END__;
+
+ cvReleaseMemStorage( &temp_storage );
+ cvFree( &points );
+
+ return box;
+}
+
projects / opencv / blobdiff