관리 메뉴

MuTa

cv::PCA 사용시 EigenValue 소팅 안된 값 찾기 본문

OpenCV

cv::PCA 사용시 EigenValue 소팅 안된 값 찾기

MuTa 2016. 10. 5. 09:48

opencv PCA 243 기준


cv::PCA 사용시 EigenValue 를 자동으로 Sort 해준다.

불필요한 Feature 를 찾으려고 했지만 Sort 해버리는 순간 정보는 모두 날아가 버린다.

95% 정도의 표현이거나 단순 5차원 으로 표현 하기 위해서 Sort 해주는것이지만

내가 필요한 기능은 아니었다.


243 소스 에서 Sort 하는 부분을 제거 한것이다. 저기 아래 주석부분

나머지는 PCA 함수에서 하는것이랑 동일

공분산 > eigen 값 구하기

cv::Mat data = matBigPot.clone(), _mean = cv::Mat();
int covar_flags = CV_COVAR_SCALE;
int len, in_count;
Size mean_sz;


len = data.cols;
in_count = data.rows;
covar_flags |= CV_COVAR_ROWS;
mean_sz = Size(len, 1);
int count = std::min(len, in_count), out_count = count;

// "scrambled" way to compute PCA (when cols(A)>rows(A)):
// B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
if (len <= in_count)
	covar_flags |= CV_COVAR_NORMAL;
int ctype = std::max(CV_32F, data.depth());
cv::Mat mean, eigenVec, eigenVal;
mean.create(mean_sz, ctype);
Mat covar(count, count, ctype);
if (_mean.data)
{
	CV_Assert(_mean.size() == mean_sz);
	_mean.convertTo(mean, ctype);
}
calcCovarMatrix(data, covar, mean, covar_flags, ctype);
eigen2(covar, false, eigenVal, eigenVec);
///////////////
 
bool eigen2(InputArray _src, bool computeEvects, OutputArray _evals, OutputArray _evects)
{
	Mat src = _src.getMat();
	int type = src.type();
	int n = src.rows;

	CV_Assert(src.rows == src.cols);
	CV_Assert(type == CV_32F || type == CV_64F);

	Mat v;
	if (computeEvects)
	{
		_evects.create(n, n, type);
		v = _evects.getMat();
	}

	size_t elemSize = src.elemSize(), astep = alignSize(n*elemSize, 16);
	AutoBuffer<uchar> buf(n*astep + n * 5 * elemSize + 32);
	uchar* ptr = alignPtr((uchar*)buf, 16);
	Mat a(n, n, type, ptr, astep), w(n, 1, type, ptr + astep*n);
	ptr += astep*n + elemSize*n;
	src.copyTo(a);
	bool ok = type == CV_32F ?
		Jacobi2(a.ptr<float>(), a.step, w.ptr<float>(), v.ptr<float>(), v.step, n, ptr) :
		Jacobi2(a.ptr<double>(), a.step, w.ptr<double>(), v.ptr<double>(), v.step, n, ptr);

	w.copyTo(_evals);
	return ok;
}
 
 
 
 
template<typename _Tp> bool
JacobiImpl_2(_Tp* A, size_t astep, _Tp* W, _Tp* V, size_t vstep, int n, uchar* buf)
{
	const _Tp eps = std::numeric_limits<_Tp>::epsilon();
	int i, j, k, m;

	astep /= sizeof(A[0]);
	if (V)
	{
		vstep /= sizeof(V[0]);
		for (i = 0; i < n; i++)
		{
			for (j = 0; j < n; j++)
				V[i*vstep + j] = (_Tp)0;
			V[i*vstep + i] = (_Tp)1;
		}
	}

	int iters, maxIters = n*n * 30;

	int* indR = (int*)alignPtr(buf, sizeof(int));
	int* indC = indR + n;
	_Tp mv = (_Tp)0;

	for (k = 0; k < n; k++)
	{
		W[k] = A[(astep + 1)*k];
		if (k < n - 1)
		{
			for (m = k + 1, mv = std::abs(A[astep*k + m]), i = k + 2; i < n; i++)
			{
				_Tp val = std::abs(A[astep*k + i]);
				if (mv < val)
					mv = val, m = i;
			}
			indR[k] = m;
		}
		if (k > 0)
		{
			for (m = 0, mv = std::abs(A[k]), i = 1; i < k; i++)
			{
				_Tp val = std::abs(A[astep*i + k]);
				if (mv < val)
					mv = val, m = i;
			}
			indC[k] = m;
		}
	}

	if (n > 1) for (iters = 0; iters < maxIters; iters++)
	{
		// find index (k,l) of pivot p
		for (k = 0, mv = std::abs(A[indR[0]]), i = 1; i < n - 1; i++)
		{
			_Tp val = std::abs(A[astep*i + indR[i]]);
			if (mv < val)
				mv = val, k = i;
		}
		int l = indR[k];
		for (i = 1; i < n; i++)
		{
			_Tp val = std::abs(A[astep*indC[i] + i]);
			if (mv < val)
				mv = val, k = indC[i], l = i;
		}

		_Tp p = A[astep*k + l];
		if (std::abs(p) <= eps)
			break;
		_Tp y = (_Tp)((W[l] - W[k])*0.5);
		_Tp t = std::abs(y) + hypot(p, y);
		_Tp s = hypot(p, t);
		_Tp c = t / s;
		s = p / s; t = (p / t)*p;
		if (y < 0)
			s = -s, t = -t;
		A[astep*k + l] = 0;

		W[k] -= t;
		W[l] += t;

		_Tp a0, b0;

#undef rotate
#define rotate(v0, v1) a0 = v0, b0 = v1, v0 = a0*c - b0*s, v1 = a0*s + b0*c

		// rotate rows and columns k and l
		for (i = 0; i < k; i++)
			rotate(A[astep*i + k], A[astep*i + l]);
		for (i = k + 1; i < l; i++)
			rotate(A[astep*k + i], A[astep*i + l]);
		for (i = l + 1; i < n; i++)
			rotate(A[astep*k + i], A[astep*l + i]);

		// rotate eigenvectors
		if (V)
			for (i = 0; i < n; i++)
				rotate(V[vstep*k + i], V[vstep*l + i]);

#undef rotate

		for (j = 0; j < 2; j++)
		{
			int idx = j == 0 ? k : l;
			if (idx < n - 1)
			{
				for (m = idx + 1, mv = std::abs(A[astep*idx + m]), i = idx + 2; i < n; i++)
				{
					_Tp val = std::abs(A[astep*idx + i]);
					if (mv < val)
						mv = val, m = i;
				}
				indR[idx] = m;
			}
			if (idx > 0)
			{
				for (m = 0, mv = std::abs(A[idx]), i = 1; i < idx; i++)
				{
					_Tp val = std::abs(A[astep*i + idx]);
					if (mv < val)
						mv = val, m = i;
				}
				indC[idx] = m;
			}
		}
	}

// 	// sort eigenvalues & eigenvectors
// 	for (k = 0; k < n - 1; k++)
// 	{
// 		m = k;
// 		for (i = k + 1; i < n; i++)
// 		{
// 			if (W[m] < W[i])
// 				m = i;
// 		}
// 		if (k != m)
// 		{
// 			std::swap(W[m], W[k]);
// 			if (V)
// 				for (i = 0; i < n; i++)
// 					std::swap(V[vstep*m + i], V[vstep*k + i]);
// 		}
// 	}

	return true;
}

static bool Jacobi2(float* S, size_t sstep, float* e, float* E, size_t estep, int n, uchar* buf)
{
	return JacobiImpl_2(S, sstep, e, E, estep, n, buf);
}

static bool Jacobi2(double* S, size_t sstep, double* e, double* E, size_t estep, int n, uchar* buf)
{
	return JacobiImpl_2(S, sstep, e, E, estep, n, buf);
}
 
//////////////////////////////////////////////////////////////////////////


'OpenCV' 카테고리의 다른 글

PCA 관련 정보들  (0) 2016.10.05
Image Watch VS2013 확장 프로그램  (0) 2016.05.13
OpenCV 기본  (0) 2016.03.30