переопределенная СЛАУ и МНК

mrgloom · August 29, 2012

вообщем тут написано, что можно решить переопределенную СЛАУ через вычисление определителей

http://helpstat.ru/2012/01/metod-naimenshix-kvadratov/

только получается для 1 переменной.

как это лучше всего сделать через opencv?

mrgloom · August 29, 2012

http://opencv.willowgarage.com/documentation/cpp/algebraic_operations.html

только откуда такая формула?

т.е. обратно то Ax=b можно восстановить, но зачем домножали на A.t() ?

Smorodov · August 29, 2012

Там использут псевдоинверсию матрицы (делается из за того что она не квадратная, а если домножить на транспонированную, то получается квадратная):

http://ru.wikipedia.org/wiki/%D0%9F%D1%81%D0%B5%D0%B2%D0%B4%D0%BE%D0%BE%D0%B1%D1%80%D0%B0%D1%82%D0%BD%D0%B0%D1%8F_%D0%BC%D0%B0%D1%82%D1%80%D0%B8%D1%86%D0%B0

mrgloom · August 29, 2012

так а зачем нам квадратная? мы СЛАУ помоему и без квадратной можем решить.

попробовал код выдает всё время нули.

void perspective_transform(vector<Point>& in_vec,vector<Point>& out_vec)

{

	double m11=0.4;

	double m12=0.1;

	double m13=12;

	double m21=0.5;

	double m22=1.3;

	double m23=10;

	double m31=0.001;

	double m32=0;

	double m33=1;


	//прямое перспективное искажение

	for(int i=0;i<in_vec.size();++i)

	{

		int x= in_vec[i].x;

		int y= in_vec[i].y;

		int u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33);

		int v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33);


		out_vec.push_back(Point(u,v));

	}

}

void get_perspective()

{

	//попробуем по найденным парам точек восстановить аффинное преобразование

	//решаем переопределенную СЛАУ

	//имеем N пар точек, значит 2*N уравнений

	//вида

	//double u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33);

	//double v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33);

	//переписываем как линейную систему относительно m{i,j}

	//m11*(x)+m12*(y)+m13*(1)+m21*(0)+m22*(0)+m23(0)+m31*(-xu)+m32*(-yu)+m33*(-u)=0

	//m11*(0)+m12*(0)+m13*(0)+m21*(x)+m22*(y)+m23(1)+m31*(-xv)+m32*(-yv)+m33*(-v)=0

	//для каждой пары точек (x,y)->(u,v) добавляем в систему 2 уравнения

	Mat A(0,9,CV_32FC1); //правильный тип?

	Mat b(0,1,CV_32FC1); //возможно неправильное кол-во колонок

	Mat t(1,9,CV_32FC1);


	//test

	vec_T.clear();//

	perspective_transform(vec_M,vec_T);//


	//для каждой точки в первом сете найдем ближайшую точку во втором

	for(int i=0;i<vec_M.size();++i)

	{

		double min_dist=INT_MAX;

		int id=-1;

		//для каждой точки проходимся по всему второму массиву ищем минимум

		for(int j=0;j<vec_T.size();++j)

		{

			double metric= sqrt(double(vec_T[j].x-vec_M[i].x)*(vec_T[j].x-vec_M[i].x)+

				(vec_T[j].y-vec_M[i].y)*(vec_T[j].y-vec_M[i].y));

			if(min_dist>metric)

			{

				min_dist=metric;

				id=j;

			}

		}


		int x= vec_M[i].x;

		int y= vec_M[i].y;

		int u= vec_T[id].x;

		int v= vec_T[id].y;

		//m11*(x)+m12*(y)+m13*(1)+m21*(0)+m22*(0)+m23(0)+m31*(-xu)+m32*(-yu)+m33*(-u)=0

		//m11*(0)+m12*(0)+m13*(0)+m21*(x)+m22*(y)+m23(1)+m31*(-xv)+m32*(-yv)+m33*(-v)=0

		t.at<float>(0,0)= x;

		t.at<float>(0,1)= y;

		t.at<float>(0,2)= 1;

		t.at<float>(0,3)= 0;

		t.at<float>(0,4)= 0;

		t.at<float>(0,5)= 0;

		t.at<float>(0,6)= -x*u;

		t.at<float>(0,7)= -y*u;

		t.at<float>(0,8)= -u;

		A.push_back(t);

		t.at<float>(0,0)= 0;

		t.at<float>(0,1)= 0;

		t.at<float>(0,2)= 0;

		t.at<float>(0,3)= x;

		t.at<float>(0,4)= y;

		t.at<float>(0,5)= 1;

		t.at<float>(0,6)= -x*v;

		t.at<float>(0,7)= -y*v;

		t.at<float>(0,8)= -v;

		A.push_back(t);

	}


	b.resize(2*vec_M.size(),0);


	cout<<A;

	cout<<b;

	Mat m = (A.t()*A).inv()*(A.t()*; //нули возможно потому, что система не решаема?

	cout<<m;

}[/code]

mrgloom · August 29, 2012

хотя скорей всего это из-за того что система действительно не имеет решений , т.к. пары точек найдены неправильно, т.к. там получается связь один ко многим.

mrgloom · August 30, 2012

или я сами уравнения не правильно составил?

т.е. я взял просто СЛАУ которая получается от перспективных преобразований, а надо было производные от суммы квадратов приравнять нулю?

Smorodov · August 30, 2012

Обычно минимизируют производную от суммы квадратов ошибки.

Я сейчас использую для минимизации (играю с нейронным разреженным автоэнкодером):

http://www.chokkan.org/software/liblbfgs/

собирается сразу, есть пример, и метод достаточно быстрый (намного быстрее градиентного спуска).

Нужно задавать саму минимизируемую функцию и её градиент (вектор частных производных минимизируемой функции по её параметрам).

Имеет смысл при нелинейной задаче высокой размерности.

ЗЫ: OpenCV умеет решать переопроеделенные системы (методом SVD).

mrgloom · August 30, 2012

кажется понял.

double u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33);
double v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33);

для affine m31=m32=0 m33=1

double u= (m11*x + m12*y + m13);
double v= (m21*x + m22*y + m23);

значит S=SUM[ (u(i)-(m11*x(i) + m12*y(i) + m13))^2+(v(i)-(m21*x(i) + m22*y(i) + m23))^2 ]

берем производные по 6 параметрам.

потом всё это переписываем в СЛАУ относительно 6 переменных.(причем получается, что для перспективных искажений то СЛАУ у нас не получится)

m11*(x*x)+m12*(y*x)+m13*(x)+m21*(0)+m22*(0)+m23(0)=u*x
m11*(x*y)+m12*(y*y)+m13*(y)+m21*(0)+m22*(0)+m23(0)=u*y
m11*(x)+m12*(y)+m13*(1)+m21*(0)+m22*(0)+m23(0)=u
m11*(0)+m12*(0)+m13*(0)+m21*(x*x)+m22*(y*x)+m23(x)=v*x
m11*(0)+m12*(0)+m13*(0)+m21*(x*y)+m22*(y*y)+m23(y)=v*y
m11*(0)+m12*(0)+m13*(0)+m21*(x)+m22*(y)+m23(1)=v

единственное смущает, что 1-3 (4-6) получаются линейно зависимы? т.е. различаются всего лишь на умножение на константу.

Smorodov · August 30, 2012

Нет они линейно-независимы.

Если умножаете на константу, то умножайте каждое слагаемое.

Там ведь нельзя общий множитель за скобку вынести так, чтобы в скобках было одно и то же.

mrgloom · August 30, 2012

m11*(x)+m12*(y)+m13*(1)+m21*(0)+m22*(0)+m23(0)=u
m11*(x*y)+m12*(y*y)+m13*(y)+m21*(0)+m22*(0)+m23(0)=u*y

ну так первое умножаем на y и получаем второе.

Smorodov · August 30, 2012

А, ну да, умножение на 0 не заметил.

mrgloom · August 30, 2012

нет всё таки я ошибся та мне будет линейной зависимости, не учёл, то что та мсуммы же везде.

но все равно не работает как надо.

//это используется в конце чтобы протестировать полученную матрицу

void affine_transform(vector<Point>& in_vec,vector<Point>& out_vec,Mat& M)

{

	out_vec.clear();

	double m11= M.at<float>(0,0);

	double m12= M.at<float>(0,1);

	double m13= M.at<float>(0,2);

	double m21= M.at<float>(1,0);

	double m22= M.at<float>(1,1);

	double m23= M.at<float>(1,2);

	double m31= 0/*M.at<float>(2,0)*/;

	double m32= 0/*M.at<float>(2,1)*/;

	double m33= 1/*M.at<float>(2,2)*/;


	//прямое перспективное искажение

	for(int i=0;i<in_vec.size();++i)

	{

		int x= in_vec[i].x;

		int y= in_vec[i].y;

		int u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33);

		int v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33);


		out_vec.push_back(Point(u,v));

	}

}

//это я использую чтобы трансформировать исходные точки для теста

void perspective_transform(vector<Point>& in_vec,vector<Point>& out_vec/*,Mat& M*/)

{

	out_vec.clear();

	//вынести в параметры?

	double m11=0.4;

	double m12=0.1;

	double m13=12;

	double m21=0.5;

	double m22=1.3;

	double m23=10;

	double m31=0;

	double m32=0;

	double m33=1;

	//double m11= M.at<float>(0,0);

	//double m12= M.at<float>(0,1);

	//double m13= M.at<float>(0,2);

	//double m21= M.at<float>(1,0);

	//double m22= M.at<float>(1,1);

	//double m23= M.at<float>(1,2);

	//double m31= 0/*M.at<float>(2,0)*/;

	//double m32= 0/*M.at<float>(2,1)*/;

	//double m33= 1/*M.at<float>(2,2)*/;


	//прямое перспективное искажение

	for(int i=0;i<in_vec.size();++i)

	{

		int x= in_vec[i].x;

		int y= in_vec[i].y;

		int u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33);

		int v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33);


		out_vec.push_back(Point(u,v));

	}

}

void get_affine()

{

	//попробуем по найденным парам точек восстановить аффинное преобразование

	//решаем СЛАУ

	//получается сколько параметров столько и уравнений(6 для аффинных)

	//вида

	//double u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33); //m31=m32=0 m33=1

	//double v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33);

	//задача поставлена как производная сумм квадратов по всем переменным =0

	//потом переписываем это всё в СЛАУ от 6 неизвестных


	//а для перспективных получается уже не СЛАУ


	Mat A(6,6,CV_32FC1);

	Mat b(6,1,CV_32FC1);


	////test

	vec_T.clear();//

	perspective_transform(vec_M,vec_T);//


	int sx=0;

	int sy=0;

	int sxx=0;

	int syy=0;

	int sxy=0;

	int sux=0;

	int suy=0;

	int svx=0;

	int svy=0;

	int sv=0;

	int su=0;

	int k=0;


	//для каждой точки в первом сете найдем ближайшую точку во втором

	for(int i=0;i<vec_M.size();++i)

	{

		//double min_dist=INT_MAX;

		//int id=-1;

		////для каждой точки проходимся по всему второму массиву ищем минимум

		//for(int j=0;j<vec_T.size();++j)

		//{

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

		//	double metric= sqrt(double(vec_T[j].x-vec_M[i].x)*(vec_T[j].x-vec_M[i].x)+

		//		(vec_T[j].y-vec_M[i].y)*(vec_T[j].y-vec_M[i].y));

		//	if(min_dist>metric)

		//	{

		//		min_dist=metric;

		//		id=j;

		//	}

		//}


		//line(img,vec_M[i],vec_T[id],cvScalar(0,0,255));


		//int x= vec_M[i].x;

		//int y= vec_M[i].y;

		//int u= vec_T[id].x;

		//int v= vec_T[id].y;


                //test

		line(img,vec_M[i],vec_T[i],cvScalar(0,0,255));

		int x= vec_M[i].x;

		int y= vec_M[i].y;

		int u= vec_T[i].x;

		int v= vec_T[i].y;


                //не забываем что там везде суммы

		//m11*(x*x)+m12*(y*x)+m13*(x)+m21*(0)+m22*(0)+m23(0)=u*x

		//m11*(x*y)+m12*(y*y)+m13*(y)+m21*(0)+m22*(0)+m23(0)=u*y

		//m11*(x)+m12*(y)+m13*(1)+m21*(0)+m22*(0)+m23(0)=u

		//m11*(0)+m12*(0)+m13*(0)+m21*(x*x)+m22*(y*x)+m23(x)=v*x

		//m11*(0)+m12*(0)+m13*(0)+m21*(x*y)+m22*(y*y)+m23(y)=v*y

		//m11*(0)+m12*(0)+m13*(0)+m21*(x)+m22*(y)+m23(1)=v


		sx+=x;

		sy+=x;

		sxx+=x*x;

		syy+=y*y;

		sux+=u*x;

		suy+=u*y;

		svx+=v*x;

		svy+=v*y;

		sv+=v;

		su+=u;

		sxy+=x*y;

		++k;

	}


	//заполняем матрицу 6х6

	//m11*(x*x)+m12*(y*x)+m13*(x)+m21*(0)+m22*(0)+m23(0)=u*x

	A.at<float>(0,0)= sxx;

	A.at<float>(0,1)= sxy;

	A.at<float>(0,2)= sx;

	A.at<float>(0,3)= 0;

	A.at<float>(0,4)= 0;

	A.at<float>(0,5)= 0;

	//m11*(x*y)+m12*(y*y)+m13*(y)+m21*(0)+m22*(0)+m23(0)=u*y

	A.at<float>(1,0)= sxy;

	A.at<float>(1,1)= syy;

	A.at<float>(1,2)= sy;

	A.at<float>(1,3)= 0;

	A.at<float>(1,4)= 0;

	A.at<float>(1,5)= 0;

	//m11*(x)+m12*(y)+m13*(1)+m21*(0)+m22*(0)+m23(0)=u

	A.at<float>(2,0)= sx;

	A.at<float>(2,1)= sy;

	A.at<float>(2,2)= 1*k;

	A.at<float>(2,3)= 0;

	A.at<float>(2,4)= 0;

	A.at<float>(2,5)= 0;

	//m11*(0)+m12*(0)+m13*(0)+m21*(x*x)+m22*(y*x)+m23(x)=v*x

	A.at<float>(3,0)= 0;

	A.at<float>(3,1)= 0;

	A.at<float>(3,2)= 0;

	A.at<float>(3,3)= sxx;

	A.at<float>(3,4)= sxy;

	A.at<float>(3,5)= sx;

	//m11*(0)+m12*(0)+m13*(0)+m21*(x*y)+m22*(y*y)+m23(y)=v*y

	A.at<float>(4,0)= 0;

	A.at<float>(4,1)= 0;

	A.at<float>(4,2)= 0;

	A.at<float>(4,3)= sxy;

	A.at<float>(4,4)= syy;

	A.at<float>(4,5)= sy;

	//m11*(0)+m12*(0)+m13*(0)+m21*(x)+m22*(y)+m23(1)=v

	A.at<float>(5,0)= 0;

	A.at<float>(5,1)= 0;

	A.at<float>(5,2)= 0;

	A.at<float>(5,3)= sx;

	A.at<float>(5,4)= sy;

	A.at<float>(5,5)= 1*k;


	b.at<float>(0,0)= sux;

	b.at<float>(1,0)= suy;

	b.at<float>(2,0)= su;

	b.at<float>(3,0)= svx;

	b.at<float>(4,0)= svy;

	b.at<float>(5,0)= sv;


	cout<<A;

	cout<<b;

        //пробуем разные методы решения

	Mat m1 = (A.t()*A).inv()*(A.t()*;

	cout<<m1;

	Mat m2;

	solve(A,b,m2,DECOMP_LU);

	cout<<m2;

	Mat m3;

	solve(A,b,m3,DECOMP_SVD);

	cout<<m3;

	Mat m4;

	solve(A,b,m4,DECOMP_QR);

	cout<<m4;


	affine_transform(vec_M,vec_T,m1.reshape(3));

	imshow("img",img);

}[/code]

mrgloom · August 30, 2012

вообщем это всё умеет делать estimateRigidTransform , но я так и не понял как это делается математически.

Войти

переопределенная СЛАУ и МНК

Recommended Posts

mrgloom 242

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

mrgloom 242

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

Smorodov 579

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

mrgloom 242

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

mrgloom 242

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

mrgloom 242

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

Smorodov 579

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

mrgloom 242

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

Smorodov 579

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

mrgloom 242

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

Smorodov 579

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

mrgloom 242

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

mrgloom 242

Поделиться сообщением

Ссылка на сообщение

Поделиться на других сайтах

Создайте учётную запись или войдите для комментирования

Создать учётную запись

Войти

Сейчас на странице 0 пользователей

Навигация

Последняя активность