X-Git-Url: https://vcs.maemo.org/git/?a=blobdiff_plain;f=src%2Fcv%2Fcvrotcalipers.cpp;fp=src%2Fcv%2Fcvrotcalipers.cpp;h=a6564e176253216d2227d3296755c3a451cb0619;hb=e4c14cdbdf2fe805e79cd96ded236f57e7b89060;hp=0000000000000000000000000000000000000000;hpb=454138ff8a20f6edb9b65a910101403d8b520643;p=opencv diff --git a/src/cv/cvrotcalipers.cpp b/src/cv/cvrotcalipers.cpp new file mode 100644 index 0000000..a6564e1 --- /dev/null +++ b/src/cv/cvrotcalipers.cpp @@ -0,0 +1,474 @@ +/*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; +} +