Центр тяжести и оси инерции

[kudukSweetone 3]

Добрый день!

Товарищи форумчане помогите разобраться:

вот кусок кода

//Находим центр тяжести

float xc=(moments.m10/moments.m00);

float yc=(moments.m01/moments.m00);


double M00=moments.m00;

double M20=moments.m20;

double M02=moments.m02;

double M11=moments.m11;


double A=(M20/M00)-xc*xc;

double B=2*((M11/M00)-xc*yc);

double C=(M02/M00)-yc*yc;


double LL=sqrt(((A+C)+sqrt(B*B+(A-C)*(A-C)))/2)*2;    //максимальная ось

double LW=sqrt(((A+C)-sqrt(B*B+(A-C)*(A-C)))/2)*2;    //минимальная ось


// Для вычисления угла нужны центральные моменты инерции

M20=moments.mu20;

M02=moments.mu02;

M11=moments.mu11;


double theta=(atan2(2*M11,(M02-M20))/2)*(180/M_PI);

cvCircle( Grab, cvPoint(xc,yc), 3, CV_RGB(0,255,255), -1, 8, 0 );    //отмечаем центр тяжести объекта

Центр тяжести объектов отмечен на рисунке Result.png

Подскажите, как нарисовать максимальную ось и минимальную ось чтобы они проходили через центр тяжести? И далее на основе отрисованных осей описать эллипс вокруг объекта?

Если можно то куском кода напишите.

[Smorodov 579]

По поводу эллипса не скажу, не помню.

А оси, это обычное двумерное преобразование. Задаем оси с центром в начале координат, поворачиваем и переносим в центр тяжести.

Можно сделать матричным умножением. http://en.wikipedia.org/wiki/Transformation_matrix

http://cse.taylor.edu/~btoll/s99/424/res/mtu/Notes/geometry/geo-tran.htm

[kudukSweetone 3]

По поводу эллипса не скажу, не помню.

А оси, это обычное двумерное преобразование. Задаем оси с центром в начале координат, поворачиваем и переносим в центр тяжести.

Можно сделать матричным умножением. http://en.wikipedia.org/wiki/Transformation_matrix

http://cse.taylor.edu/~btoll/s99/424/res/mtu/Notes/geometry/geo-tran.htm

Это кажется не совсем то, что мне нужно или я не правильно понял.

double LL=sqrt(((A+C)+sqrt(B*B+(A-C)*(A-C)))/2)*2;    //максимальная ось

double LW=sqrt(((A+C)-sqrt(B*B+(A-C)*(A-C)))/2)*2;          //минимальная ось

Мне нужно чтобы не только рассчитывались длины максимальной и мининимальной осей (LL и LW), но и сразу рисовалась бы линия. Чтобы в итоге получилось приблизительно вот как на рисунке внизу (сделал в Пейнте) , а нужно программно. Как это программно осуществить? Я чего-то с координатами запутался и не могу сообразить - какие координаты у начальной и конечной точки макс и мин осей? Вот весь код программы целиком:

#include "stdafx.h"

#include <cv.h>  

#include <cxcore.h>  

#include <highgui.h>

#include <iostream> 

#include <string>

#include <sstream>


using namespace std;  

const int N=20;

float Radius[N]={0};

float X[N]={0};

float Y[N]={0};

float AverageRadius=0;

float AverageX=0;

float AverageY=0;

const char *windowName = "Моменты и инварианты"; 


CvFont font; // Описатель шрифта

IplImage *frame, *frame_copy = 0;

IplImage *gray = NULL;


int Val=100,Val1=100; // Значение переменной (для разных нужд)

float edge_thresh=100;

void HuMoments(IplImage* frame);		//функция обработки изображения

#define WIDTHBYTES(bits) ((((bits) + 31) / 32) * 4)


//------------------------------------Начало моментов-----------------------------------------//


void HuMoments(IplImage* Grab)

{

	float* p;

	float M_PI=3.1415;

	float XC[1000]={0};

	float YC[1000]={0};

	float S[1000]={0};

	double mt=0;

	int ID=0;

	int NumCont=0;

	int NumCont1=0;

	char * text=(char *) malloc(20);

	float ma[6]={0};

	float mb[6]={0};


	CvMemStorage* storage = 0;

	CvSeq* contours;

	CvSeq* Templ;

	CvSeq* result;

	storage = cvCreateMemStorage(0);


	CvSeq* polygons = cvCreateSeq((CV_SEQ_FLAG_CLOSED|CV_SEQ_KIND_CURVE|CV_SEQ_ELTYPE_POINT), sizeof(CvContour), sizeof(CvPoint), storage); 

	CvMoments moments;

	CvHuMoments hu_moments;


	cvCvtColor(Grab, gray, CV_BGR2GRAY);	// Получаем серый цвет

	cvThreshold(gray, gray, Val,255, /*CV_THRESH_BINARY*/ CV_THRESH_OTSU);		


	cvRectangle(gray, cvPoint(0,0), cvPoint(gray->width-1,gray->height-1),CV_RGB(0,0,0)); // Необходимо удалить белый бордюр

	cvDilate(gray, gray, 0, 1);		// Утолщаем контуры


	int NumberCont = cvFindContours( gray, storage,&contours,sizeof(CvContour),CV_RETR_TREE, CV_CHAIN_APPROX_NONE,cvPoint(0,0)); 

	printf("Contours found = %d",NumberCont); // количество найденных контуров

	cout<<endl;

	if(contours!=0)

	{NumCont=contours->total;}


	int NC=0;

	for(;contours!=0;contours = contours->h_next)

	{

		// Аппроксимация контуров полигонами

		result = cvApproxPoly( contours, sizeof(CvContour), storage,CV_POLY_APPROX_DP, cvContourPerimeter(contours)*0.001, 0 );

		double area=fabs(cvContourArea(result,CV_WHOLE_SEQ));	// Площадь полигона

		CvSeqReader reader;	// Объявление читалки точек

		cvStartReadSeq( result, &reader, 0 );	// Инициализация читалки точек

		CvPoint pt[2];	// Две точки

		if(area>500)

		{

			cvMoments( result, &moments);			//Вычисляем моменты


			// Центр тяжести

			float xc=(moments.m10/moments.m00);

			float yc=(moments.m01/moments.m00);



			// Обрубаем клонов

			bool ok=1;

			for (int o=0;o<NC;o++)

			{

				if( ((xc-XC[o])*(xc-XC[o])+(yc-YC[o])*(yc-YC[o]))<10 && fabs((area-S[o])/(area+S[o]))<0.3){ok=0;}

			}

			//

			if(ok)

			{

				if(NumCont1==0)	// Запоминаем шаблон

				{

				Templ=cvCloneSeq(result);} //Templ-сюда записывается значение первого попавшегося контура, а мне нужно самому выбрать

				mt=cvMatchShapes( Templ, result,CV_CONTOURS_MATCH_I2);   // Сравниваем с найденными контурами


				NumCont1++;

				double M00=moments.m00;

				double M20=moments.m20;

				double M02=moments.m02;

				double M11=moments.m11;


				double A=(M20/M00)-xc*xc;

				double B=2*((M11/M00)-xc*yc);

				double C=(M02/M00)-yc*yc;


				double LL=sqrt(((A+C)+sqrt(B*B+(A-C)*(A-C)))/2)*2; //максимальная ось

				double LW=sqrt(((A+C)-sqrt(B*B+(A-C)*(A-C)))/2)*2; //минимальная ось


				// Для вычисления угла нужны центральные моменты инерции

				M20=moments.mu20;	//центральные моменты инерции

				M02=moments.mu02;	//...

				M11=moments.mu11;	//...


				double theta=(atan2(2*M11,(M02-M20))/2)*(180/M_PI);	//theta - это угол наклона главной (большей) оси

				cvCircle( Grab, cvPoint(xc,yc), 3, CV_RGB(0,255,255), -1, 8, 0 );	//отмечаем центр объекта


				for (int i=0;i<result->total-1;i++)	// Достаем точки из хранилища точек и рисуем линии

				{

					CV_READ_SEQ_ELEM( pt[0], reader );

					CV_READ_SEQ_ELEM( pt[1], reader );


				if(mt>0.2){cvLine(Grab,pt[0],pt[1],CV_RGB(255,0,0),2);}	//mt>(порог)-больше которого объект выделяется красным

				else{cvLine(Grab,pt[0],pt[1],CV_RGB(0,255,0),2);}	//прорисовка красных линий для отбракованных объектов

				}

				sprintf(text,"%.4f", static_cast<float>(mt));			//переконвертируем из double в char


			}

			XC[NC]=xc;

			YC[NC]=yc;

			S[NC]=area;

			NC++;

		} // if

	} // for

	cvReleaseMemStorage(&storage);

	cvNamedWindow(windowName, CV_WINDOW_AUTOSIZE);  

	cvShowImage(windowName, Grab);

	cvShowImage("Gray", gray);

}

//------------------------------------------Конец функции расчета моментов -------------------------------------------//


int _tmain(int argc, _TCHAR* argv[])

{

	cvInitFont(&font, CV_FONT_HERSHEY_SIMPLEX,0.5,0.5,0,1,8);	

	IplImage* frame = cvLoadImage("Source.jpg", 1);	

	gray = cvCreateImage(cvSize(frame->width,frame->height), IPL_DEPTH_8U, 1);

 	HuMoments(frame);	

	while(1)

	{if (cvWaitKey(33) > 0) break;} 


	cvReleaseImage(&frame);

	cvReleaseImage( &frame_copy );

	cvReleaseImage( &gray );

	cvDestroyWindow(windowName);  

	return 0;  

}

[mrgloom 242]

во http://www.compvision.ru/forum/index.php?showtopic=173&view=findpost&p=1536

http://www.compvision.ru/forum/index.php?showtopic=4&view=findpost&p=216

[kudukSweetone 3]

To Smorodov:

А у вас случайно не осталось листинга программы изображенной на скриншоте второй страницы вот этой вот статьи? (в прикрепленном файле)

Написано, что это программа обнаружения луча лазерной указки со снятыми ограничениями по размеру луча.

Это как раз то, что мне нужно - там угол наклона главной оси прорисован и описан эллипс.

О моментах контура.pdf

[Smorodov 579]

Была где-то, вечером поищу.

Осей там нет (и в pdf-ке тоже нет), наверное они мне не были нужны, эллипс рисуется так:

// Рисуем эллипс

cvEllipse(frame, cvPoint(xc,yc), cvSize(LW,LL), theta,0,360, CV_RGB(0,255,255));[/code]


Длины осей - это всего лишь длины, а не координаты. Получить координаты крайних точек можно если отложить от качала координат (0;0) линии длиной (LW и LL) вдоль осей координат (может быть деленные на 2, на помню).

Теперь у Вас есть 4 точки.

Каждая точка имеет координаты  (x,y,1) - вектор-столбец.


Чтобы получить трансформированную точку (у Вас поворот и перенос) нужно умножить матрицу трансформации на вектор координат. Повторив 4 раза, получите 4 точки по которым можно нарисовать оси.