FCM GMM 肤色检测 C++代码

/* * MyGMM.h * MyAlgorithm * * Created by on 10-1-12. * Copyright 2010 __MyCompanyName__. All rights reserved. * */ #ifndef _MYGMM_H_ #define _MYGMM_H_ #include "MyAlgoBase.h" template <typename T> class MyGMM { public: MyGMM(void); MyGMM(const MyMatrix<T>& mSamples, const MyMatrix<T>& mPartition, const MyMatrix<T>& mCenter); ~MyGMM(void); void Clear(void); void SetPara(const MyMatrix<T>& mSamples, const MyMatrix<T>& mPartition, const MyMatrix<T>& mCenter); bool IsEstablished(void) const; double CalculateProbability(const MyArray<T>& testVector) const; protected: bool Establish(const MyMatrix<T>& mPartition, const MyMatrix<T>& mCenter); void InitEM(const MyMatrix<T>& mPartition, const MyMatrix<T>& mCenter); double CalculatePosteriori(UINT nComponentIndex, const MyArray<T>& testVector) const; double CalculateGaussian(const MyArray<T>& testVector, const MyArray<T>& meanVector, const MyMatrix<T>& covMatrix) const; bool NeedBreakEM(const MyArray<T>& aWeight, const MyArray<MyArray<T> >& aMeans, const MyArray<MyMatrix<T> >& aCovMatrix) const; void AdjustCovariance(MyMatrix<T>& covMatrix) const; private: const static UINT _nMaxEMTimes = 10000; const static double _dMinVariance; UINT _nK; //num of GMM components UINT _nN; //num of samples UINT _nD; //dimension of sample vectors MyMatrix<T> _mSamples; MyArray<T> _aWeight; MyArray<MyArray<T> > _aMeans; MyArray<MyMatrix<T> > _aCovMatrix; bool _bIsEstablished; }; template <typename T> const double MyGMM<T>::_dMinVariance = 1e-6; template <typename T> MyGMM<T>::MyGMM(void) :_bIsEstablished(false), _nK(0), _nN(0), _nD(0) { } template <typename T> MyGMM<T>::MyGMM(const MyMatrix<T>& mSamples, const MyMatrix<T>& mPartition, const MyMatrix<T>& mCenter) :_bIsEstablished(false), _nK(0), _nN(0), _nD(0) { if((mSamples.GetNum() != 0) &&(mPartition.GetNum() != 0) &&(mCenter.GetNum() != 0) &&(mSamples.GetColumnSize() >= 2) &&(mSamples.GetRowSize() == mCenter.GetRowSize()) &&(mSamples.GetColumnSize() == mPartition.GetColumnSize()) &&(mPartition.GetRowSize() == mCenter.GetColumnSize()) ) { _nK = mPartition.GetRowSize(); _nN = mSamples.GetColumnSize(); _nD = mSamples.GetRowSize(); _mSamples = mSamples; if(Establish(mPartition,mCenter)) { _bIsEstablished = true; } } } template <typename T> MyGMM<T>::~MyGMM(void) { Clear(); } template <typename T> void MyGMM<T>::Clear(void) { _bIsEstablished = false; _mSamples.Clear(); _aWeight.RemoveAll(); _aMeans.RemoveAll(); _aCovMatrix.RemoveAll(); } template <typename T> void MyGMM<T>::SetPara(const MyMatrix<T>& mSamples, const MyMatrix<T>& mPartition, const MyMatrix<T>& mCenter) { _bIsEstablished = false; _nK = 0; _nN = 0; _nD = 0; if((mSamples.GetNum() != 0) &&(mPartition.GetNum() != 0) &&(mCenter.GetNum() != 0) &&(mSamples.GetColumnSize() >= 2) &&(mSamples.GetRowSize() == mCenter.GetRowSize()) &&(mSamples.GetColumnSize() == mPartition.GetColumnSize()) &&(mPartition.GetRowSize() == mCenter.GetColumnSize()) ) { _nK = mPartition.GetRowSize(); _nN = mSamples.GetColumnSize(); _nD = mSamples.GetRowSize(); _mSamples = mSamples; if(Establish(mPartition,mCenter)) { _bIsEstablished = true; } } } template <typename T> bool MyGMM<T>::IsEstablished(void) const { return _bIsEstablished; } template <typename T> double MyGMM<T>::CalculateProbability(const MyArray<T>& testVector) const { double ret = 0.0; if(!IsEstablished()) { ret = -1.0; } else { if(testVector.GetNum() != _nD) { ret = -1.0; } else { double gk = 0.0; for(UINT k = 1; k <= _nK; ++k) { gk = CalculateGaussian(testVector,_aMeans[k-1],_aCovMatrix[k-1]); ret += _aWeight[k-1] * gk; } } } return ret; } /////////protected fun template <typename T> bool MyGMM<T>::Establish(const MyMatrix<T>& mPartition, const MyMatrix<T>& mCenter) { bool ret = false; UINT t = 0; InitEM(mPartition, mCenter); for(t = 1; t<= _nMaxEMTimes; ++t) { cout<<"EM Weight is: "<<endl; _aWeight.Display(); cout<<"EM Means is:"<<endl; for(UINT i=0; i<_aMeans.GetNum();++i) { _aMeans[i].Display(); } cout<<"EM CovMatrix is: "<<endl; for(UINT i=0; i<_aCovMatrix.GetNum();++i) { _aCovMatrix[i].Display(); } MyArray<T> aNewWeight(_nK); MyArray<MyArray<T> > aNewMeans(_nK); MyArray<MyMatrix<T> > aNewCovMatrix(_nK); for(UINT k = 1; k <= _nK; ++k) { cout<<"In EM iter times "<<t<<" k = "<<k<<endl; double dPosterioriSum = 0.0; MyArray<double> aPosteriori(_nN); //calculate weight for(UINT i = 1; i <= _nN; ++i) { aPosteriori[i-1] = CalculatePosteriori(k,_mSamples.GetColumn(i)); if(aPosteriori[i-1] < 0) { return false; } dPosterioriSum += aPosteriori[i-1]; } if(fabs(dPosterioriSum) < MyMinLimit) { break; } aNewWeight[k-1] = dPosterioriSum/(double)_nN; //calculate means MyMatrix<T> sumMeans(_nD,1); for(UINT i = 1; i<= _nN; ++i) { MyArray<T> xi = _mSamples.GetColumn(i); sumMeans = sumMeans + MyMatrix<T>(xi) * MyMatrix<T>(aPosteriori[i-1]); } sumMeans = sumMeans / dPosterioriSum; aNewMeans[k-1] = sumMeans.GetColumn(1); //calculate covariance matrix MyMatrix<T> sumCov(_nD,_nD); for(UINT i = 1; i<= _nN; ++i) { MyArray<T> xi = _mSamples.GetColumn(i); MyMatrix<T> vectorDiff(_nD,1); vectorDiff = MyMatrix<T>(xi) - MyMatrix<T>(_aMeans[k-1]); sumCov = sumCov + (MyMatrix<T>(aPosteriori[i-1])*vectorDiff*Transpose(vectorDiff)); } sumCov = sumCov / dPosterioriSum; AdjustCovariance(sumCov); aNewCovMatrix[k-1] = sumCov; } if(NeedBreakEM(aNewWeight,aNewMeans,aNewCovMatrix)) { _aWeight = aNewWeight; _aMeans = aNewMeans; _aCovMatrix = aNewCovMatrix; ret = true; cout<<"EM converge"<<endl; break; } else { _aWeight = aNewWeight; _aMeans = aNewMeans; _aCovMatrix = aNewCovMatrix; } } cout<<"EM iter times is "<<t<<endl; cout<<"EM last Weight is: "<<endl; _aWeight.Display(); cout<<"EM last Means is:"<<endl; for(UINT i=0; i<_aMeans.GetNum();++i) { _aMeans[i].Display(); } cout<<"EM last CovMatrix is: "<<endl; for(UINT i=0; i<_aCovMatrix.GetNum();++i) { _aCovMatrix[i].Display(); } return ret; } template <typename T> void MyGMM<T>::InitEM(const MyMatrix<T>& mPartition, const MyMatrix<T>& mCenter) { //mSamples: D*N, mPartition: K*N, mCenter: D*K //init weight array _aWeight.Resize(_nK); _aWeight.SetAllElements(1/(T)_nK); //init mean array _aMeans.Resize(_nK); for(UINT k=1; k<=_nK; ++k) { _aMeans[k-1] = mCenter.GetColumn(k); } //init cov matrix array _aCovMatrix.Resize(_nK); for(UINT k=1; k<=_nK; ++k) { MyMatrix<T> cov(_nD,_nD); for(UINT i=1; i<=_nN; ++i) { //if(mPartition(k,i) > 0.9) { MyMatrix<T> vectorDiff = MyMatrix<T>(_mSamples.GetColumn(i)) - MyMatrix<T>(mCenter.GetColumn(k)); cov = cov + (vectorDiff * Transpose(vectorDiff)); //} } cov = cov /((T)(_nN-1)); AdjustCovariance(cov); _aCovMatrix[k-1] = cov; } } template <typename T> double MyGMM<T>::CalculatePosteriori(UINT nComponentIndex, const MyArray<T>& testVector) const { double ret = 0.0; if((nComponentIndex < 1) ||(nComponentIndex > _nK) ||(testVector.GetNum() != _nD) ) { ret = -1; } else { double numerator = 0.0; double denominator = 0.0; double gk = 0.0; gk = CalculateGaussian(testVector,_aMeans[nComponentIndex-1],_aCovMatrix[nComponentIndex-1]); numerator = _aWeight[nComponentIndex-1] * gk; for(UINT k=1; k<=_nK; ++k) { gk = CalculateGaussian(testVector,_aMeans[k-1],_aCovMatrix[k-1]); denominator += _aWeight[k-1] * gk; } if(fabs(denominator) > MyMinLimit) { ret = numerator/denominator; } else { cout<<"CalculatePosteriori error"<<endl; ret = 0; } } return ret; } template <t