LevenbergMarquardt

cvLevenbergMarquardtOptimizationで最小化

void cvLevenbergMarquardtOptimization†

void cvLevenbergMarquardtOptimization(pointer_LMJac JacobianFunction, pointer_LMFunc function, CvMat *X0, CvMat *observRes, CvMat *resultX, int maxIter, double epsilon);

な，長い・・・

引数†

• JacobianFunction:pointer_LMJac型のヤコビアン行列を返す関数のポインタ
• function:pointer_LMJac型の計算値を返す関数のポインタ
• X0:CvMat*型の初期値
• observRes:CvMat*型の観測値
• resultX:CvMat*型の結果
• maxIter:int型の最大反復回数
• epsilon:double型の収束しきい値

解説†

• pointer_LMJac型とはfunction( const CvMat*src, CvMat* dst)という形の関数
• 関数を2つ作る
  • srcを引数としてdstに計算値を代入する関数
  • srcを引数としてdstに編微分を羅列した関数行列を返す関数
• X0には初期値，observResは観測値，あとは反復回数と収束しきい値を指定する

• 関数について
  • 例えばP行列は未知数11の行列で，計算値は2次元座標
  • PX=x
  • ここでPが未知でXが既知，xが観測できる値
  • 例えば10点の座標からPを推定する場合，srcは11行1列のベクトルになり，dstは10*2(次元)で20行1列のベクトルとなる．
  • 既知の3次元座標はあらかじめ何とか渡しておかなければならない．
• ヤコビ行列を返す関数は理論的には要らないはず(Levenberg-Marquradt法は本来そういうもの)なのだが，何故か関数を求められる上に，NULLポインタでごまかすこともできない．

返り値†

• voidなので返さない

サンプルコード†

#geshi(c++,number){{

//cvaux/src/cvlevmartif.cppにTriFocalを求める過程でLevenbergMarqurdtが求められている．
#include <cvlevmar.cpp>
#include <cv.h>
#include <cvaux.h>
#include <cxcore.h>
#include <highgui.h>
#include <stdio.h>
#include <math.h>
#define	A	100
#define	B	102
#define	PRECISION	100
#define	DERIV_STEP	0.00001

void	CalcDiff( const CvMat *src, CvMat *dst);
void	CalcDerivDiff( const CvMat *src, CvMat *dst);
double	coreCalculation( const CvMat *src, double x);

int	main(){
	CvMat	*src, *dst, *X0, *res;
	FILE *fp = fopen("./data", "w");
	fclose(fp);
	X0	= cvCreateMat(2, 1, CV_64F);
	src	= cvCreateMat(2, 1, CV_64F);
	dst	= cvCreateMat(PRECISION, 1, CV_64F);
	res	= cvCreateMat(X0->rows, X0->cols, CV_64F);
	cvmSet(X0, 0, 0, 100);
	cvmSet(X0, 1, 0, 102);
	CalcDiff(X0, dst);
	cvmSet(X0, 0, 0, 0);
	cvmSet(X0, 1, 0, 0);
	cvLevenbergMarquardtOptimization((pointer_LMJac)CalcDerivDiff, (pointer_LMFunc)CalcDiff, X0, dst, res, 50, 0.01);
	return 0;
}

void	CalcDiff( const CvMat *src, CvMat *dst){
	int	i;
	double	x, val;
	FILE *fp = fopen("./data", "a");
	static	count = 0;
	printf("%d repeat\n", count++);
	for(i = 0;i < PRECISION;i++){
		x = CV_PI * ((double)i / PRECISION);
		val = coreCalculation(src, x);
		cvmSet(dst, i, 0, val);
		fprintf(fp, "%f\t", cvmGet(dst, i, 0));
	}
	fprintf(fp, "\n");
}

void	CalcDerivDiff( const CvMat *src, CvMat *dst){
	int	i, j;
	double	x, val, nextval;
	CvMat	*next = cvCreateMat(src->rows,src->cols,src->type);
	for(j = 0;j < 2;j++){
		cvCopy( src, next);
		cvmSet(next, j, 0, cvmGet(src, j, 0) + DERIV_STEP);
		for(i = 0;i < PRECISION;i++){
			x = CV_PI * ((double)i / PRECISION);
			val = coreCalculation(src, x);
			nextval = coreCalculation(next, x);
			cvmSet(dst, i, j, nextval - val);
		}
	}
}

double	coreCalculation( const CvMat *src, double x){
	double a = cvmGet(src, 0, 0);
	double b = cvmGet(src, 1, 0);
	//return a*cos(b*x) + b*sin(a*x);
	return a*cos(3*x) + b*sin(60*x);
}

}}

実体ファイル†

• OpenCV 1.1以前
  • cvaux/src/cvlevmartif.cpp
  • cvaux/src/cvlevmar.cpp
• OpenCV 2.1
  • src/cv/cvcalibcation.cpp
  • include/cv.hpp
• OpenCV 2.2以降
  • opencv2/calib3d/calib3d.hpp
  • modules/calib3d/src/calibration.cpp

ジャンル:OpenCV:OpenCV 1.0:OpenCV 2.1:OpenCV 2.2:OpenCV 2.3準拠

最新の20件