Chiebot-Mirror
/
opencv
mirror of https://github.com/opencv/opencv.git
8
0
Fork 0
Code Issues Projects Releases Wiki Activity

added CvEM read/write (#1032)

Browse Source
pull/13383/head
Maria Dimashova 14 years ago
parent bd33e0a3da
commit 3dc03531e1
2 changed files with 333 additions and 19 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 14
      modules/ml/include/opencv2/ml/ml.hpp
  2. 338
      modules/ml/src/em.cpp

14
modules/ml/include/opencv2/ml/ml.hpp
Unescape View File

@ -609,6 +609,7 @@ public:
CV_OUT cv::Mat* labels=0 ); CV_OUT cv::Mat* labels=0 );
CV_WRAP virtual float predict( const cv::Mat& sample, CV_OUT cv::Mat* probs=0 ) const; CV_WRAP virtual float predict( const cv::Mat& sample, CV_OUT cv::Mat* probs=0 ) const;
CV_WRAP virtual double calcLikelihood( const cv::Mat &sample ) const;
CV_WRAP int getNClusters() const; CV_WRAP int getNClusters() const;
CV_WRAP cv::Mat getMeans() const; CV_WRAP cv::Mat getMeans() const;
@ -616,7 +617,8 @@ public:
CV_WRAP cv::Mat getWeights() const; CV_WRAP cv::Mat getWeights() const;
CV_WRAP cv::Mat getProbs() const; CV_WRAP cv::Mat getProbs() const;
CV_WRAP inline double getLikelihood() const { return log_likelihood; }; CV_WRAP inline double getLikelihood() const { return log_likelihood; }
CV_WRAP inline double getLikelihoodDelta() const { return log_likelihood_delta; }
#endif #endif
CV_WRAP virtual void clear(); CV_WRAP virtual void clear();
@ -627,12 +629,19 @@ public:
const CvMat* get_weights() const; const CvMat* get_weights() const;
const CvMat* get_probs() const; const CvMat* get_probs() const;
inline double get_log_likelihood () const { return log_likelihood; }; inline double get_log_likelihood() const { return log_likelihood; }
inline double get_log_likelihood_delta() const { return log_likelihood_delta; }
// inline const CvMat * get_log_weight_div_det () const { return log_weight_div_det; }; // inline const CvMat * get_log_weight_div_det () const { return log_weight_div_det; };
// inline const CvMat * get_inv_eigen_values () const { return inv_eigen_values; }; // inline const CvMat * get_inv_eigen_values () const { return inv_eigen_values; };
// inline const CvMat ** get_cov_rotate_mats () const { return cov_rotate_mats; }; // inline const CvMat ** get_cov_rotate_mats () const { return cov_rotate_mats; };
virtual void read( CvFileStorage* fs, CvFileNode* node );
virtual void write( CvFileStorage* fs, const char* name ) const;
virtual void write_params( CvFileStorage* fs ) const;
virtual void read_params( CvFileStorage* fs, CvFileNode* node );
protected: protected:
virtual void set_params( const CvEMParams& params, virtual void set_params( const CvEMParams& params,
@ -645,6 +654,7 @@ protected:
const CvMat* means ); const CvMat* means );
CvEMParams params; CvEMParams params;
double log_likelihood; double log_likelihood;
double log_likelihood_delta;
CvMat* means; CvMat* means;
CvMat** covs; CvMat** covs;

338
modules/ml/src/em.cpp
Unescape View File

@ -115,6 +115,221 @@ void CvEM::clear()
} }
} }
void CvEM::read( CvFileStorage* fs, CvFileNode* node )
{
bool ok = false;
CV_FUNCNAME( "CvEM::read" );
__BEGIN__;
clear();
size_t data_size;
CvEMParams _params;
CvSeqReader reader;
CvFileNode* em_node = 0;
CvFileNode* tmp_node = 0;
CvSeq* seq = 0;
CvMat **tmp_covs = 0;
CvMat **tmp_cov_rotate_mats = 0;
read_params( fs, node );
em_node = cvGetFileNodeByName( fs, node, "cvEM" );
if( !em_node )
CV_ERROR( CV_StsBadArg, "cvEM tag not found" );
CV_CALL( weights = (CvMat*)cvReadByName( fs, em_node, "weights" ));
CV_CALL( means = (CvMat*)cvReadByName( fs, em_node, "means" ));
CV_CALL( log_weight_div_det = (CvMat*)cvReadByName( fs, em_node, "log_weight_div_det" ));
CV_CALL( inv_eigen_values = (CvMat*)cvReadByName( fs, em_node, "inv_eigen_values" ));
// Size of all the following data
data_size = _params.nclusters*2*sizeof(CvMat*);
CV_CALL( tmp_covs = (CvMat**)cvAlloc( data_size ));
memset( tmp_covs, 0, data_size );
tmp_cov_rotate_mats = tmp_covs + params.nclusters;
CV_CALL( tmp_node = cvGetFileNodeByName( fs, em_node, "covs" ));
seq = tmp_node->data.seq;
if( !CV_NODE_IS_SEQ(tmp_node->tag) || seq->total != params.nclusters)
CV_ERROR( CV_StsParseError, "Missing or invalid sequence of covariance matrices" );
CV_CALL( cvStartReadSeq( seq, &reader, 0 ));
for( int i = 0; i < params.nclusters; i++ )
{
CV_CALL( tmp_covs[i] = (CvMat*)cvRead( fs, (CvFileNode*)reader.ptr ));
CV_NEXT_SEQ_ELEM( seq->elem_size, reader );
}
CV_CALL( tmp_node = cvGetFileNodeByName( fs, em_node, "cov_rotate_mats" ));
seq = tmp_node->data.seq;
if( !CV_NODE_IS_SEQ(tmp_node->tag) || seq->total != params.nclusters)
CV_ERROR( CV_StsParseError, "Missing or invalid sequence of rotated cov. matrices" );
CV_CALL( cvStartReadSeq( seq, &reader, 0 ));
for( int i = 0; i < params.nclusters; i++ )
{
CV_CALL( tmp_cov_rotate_mats[i] = (CvMat*)cvRead( fs, (CvFileNode*)reader.ptr ));
CV_NEXT_SEQ_ELEM( seq->elem_size, reader );
}
covs = tmp_covs;
cov_rotate_mats = tmp_cov_rotate_mats;
ok = true;
__END__;
if (!ok)
clear();
}
void CvEM::read_params( CvFileStorage *fs, CvFileNode *node)
{
CV_FUNCNAME( "CvEM::read_params");
__BEGIN__;
size_t data_size;
CvEMParams _params;
CvSeqReader reader;
CvFileNode* param_node = 0;
CvFileNode* tmp_node = 0;
CvSeq* seq = 0;
const char * start_step_name = 0;
const char * cov_mat_type_name = 0;
param_node = cvGetFileNodeByName( fs, node, "params" );
if( !param_node )
CV_ERROR( CV_StsBadArg, "params tag not found" );
CV_CALL( start_step_name = cvReadStringByName( fs, param_node, "start_step", 0 ) );
CV_CALL( cov_mat_type_name = cvReadStringByName( fs, param_node, "cov_mat_type", 0 ) );
if( start_step_name )
_params.start_step = strcmp( start_step_name, "START_E_STEP" ) == 0 ? START_E_STEP :
strcmp( start_step_name, "START_M_STEP" ) == 0 ? START_M_STEP :
strcmp( start_step_name, "START_AUTO_STEP" ) == 0 ? START_AUTO_STEP : 0;
else
CV_CALL( _params.start_step = cvReadIntByName( fs, param_node, "start_step", -1 ) );
if( cov_mat_type_name )
_params.cov_mat_type = strcmp( cov_mat_type_name, "COV_MAT_SPHERICAL" ) == 0 ? COV_MAT_SPHERICAL :
strcmp( cov_mat_type_name, "COV_MAT_DIAGONAL" ) == 0 ? COV_MAT_DIAGONAL :
strcmp( cov_mat_type_name, "COV_MAT_GENERIC" ) == 0 ? COV_MAT_GENERIC : 0;
else
CV_CALL( _params.cov_mat_type = cvReadIntByName( fs, param_node, "cov_mat_type", -1) );
CV_CALL( _params.nclusters = cvReadIntByName( fs, param_node, "nclusters", -1 ));
CV_CALL( _params.weights = (CvMat*)cvReadByName( fs, param_node, "weights" ));
CV_CALL( _params.means = (CvMat*)cvReadByName( fs, param_node, "means" ));
data_size = _params.nclusters*sizeof(CvMat*);
CV_CALL( _params.covs = (const CvMat**)cvAlloc( data_size ));
memset( _params.covs, 0, data_size );
CV_CALL( tmp_node = cvGetFileNodeByName( fs, param_node, "covs" ));
seq = tmp_node->data.seq;
if( !CV_NODE_IS_SEQ(tmp_node->tag) || seq->total != _params.nclusters)
CV_ERROR( CV_StsParseError, "Missing or invalid sequence of covariance matrices" );
CV_CALL( cvStartReadSeq( seq, &reader, 0 ));
for( int i = 0; i < _params.nclusters; i++ )
{
CV_CALL( _params.covs[i] = (CvMat*)cvRead( fs, (CvFileNode*)reader.ptr ));
CV_NEXT_SEQ_ELEM( seq->elem_size, reader );
}
params = _params;
__END__;
}
void CvEM::write_params( CvFileStorage* fs ) const
{
CV_FUNCNAME( "CvEM::write_params" );
__BEGIN__;
const char* cov_mat_type_name =
(params.cov_mat_type == COV_MAT_SPHERICAL) ? "COV_MAT_SPHERICAL" :
(params.cov_mat_type == COV_MAT_DIAGONAL) ? "COV_MAT_DIAGONAL" :
(params.cov_mat_type == COV_MAT_GENERIC) ? "COV_MAT_GENERIC" : 0;
const char* start_step_name =
(params.start_step == START_E_STEP) ? "START_E_STEP" :
(params.start_step == START_M_STEP) ? "START_M_STEP" :
(params.start_step == START_AUTO_STEP) ? "START_AUTO_STEP" : 0;
CV_CALL( cvStartWriteStruct( fs, "params", CV_NODE_MAP ) );
if( cov_mat_type_name )
{
CV_CALL( cvWriteString( fs, "cov_mat_type", cov_mat_type_name) );
}
else
{
CV_CALL( cvWriteInt( fs, "cov_mat_type", params.cov_mat_type ) );
}
if( start_step_name )
{
CV_CALL( cvWriteString( fs, "start_step", start_step_name) );
}
else
{
CV_CALL( cvWriteInt( fs, "cov_mat_type", params.start_step ) );
}
CV_CALL( cvWriteInt( fs, "nclusters", params.nclusters ));
CV_CALL( cvWrite( fs, "weights", weights ));
CV_CALL( cvWrite( fs, "means", means ));
CV_CALL( cvStartWriteStruct( fs, "covs", CV_NODE_SEQ ));
for( int i = 0; i < params.nclusters; i++ )
CV_CALL( cvWrite( fs, NULL, covs[i] ));
CV_CALL( cvEndWriteStruct( fs ) );
// Close params struct
CV_CALL( cvEndWriteStruct( fs ) );
__END__;
}
void CvEM::write( CvFileStorage* fs, const char* name ) const
{
CV_FUNCNAME( "CvEM::write" );
__BEGIN__;
CV_CALL( cvStartWriteStruct( fs, name, CV_NODE_MAP, CV_TYPE_NAME_ML_EM ) );
write_params(fs);
CV_CALL( cvStartWriteStruct( fs, "cvEM", CV_NODE_MAP ) );
CV_CALL( cvWrite( fs, "means", means ) );
CV_CALL( cvWrite( fs, "weights", weights ) );
CV_CALL( cvWrite( fs, "log_weight_div_det", log_weight_div_det ) );
CV_CALL( cvWrite( fs, "inv_eigen_values", inv_eigen_values ) );
CV_CALL( cvStartWriteStruct( fs, "covs", CV_NODE_SEQ ));
for( int i = 0; i < params.nclusters; i++ )
CV_CALL( cvWrite( fs, NULL, covs[i] ));
CV_CALL( cvEndWriteStruct( fs ));
CV_CALL( cvStartWriteStruct( fs, "cov_rotate_mats", CV_NODE_SEQ ));
for( int i = 0; i < params.nclusters; i++ )
CV_CALL( cvWrite( fs, NULL, cov_rotate_mats[i] ));
CV_CALL( cvEndWriteStruct( fs ) );
// close cvEM
CV_CALL( cvEndWriteStruct( fs ) );
// close top level
CV_CALL( cvEndWriteStruct( fs ) );
__END__;
}
void CvEM::set_params( const CvEMParams& _params, const CvVectors& train_data ) void CvEM::set_params( const CvEMParams& _params, const CvVectors& train_data )
{ {
@ -203,6 +418,78 @@ void CvEM::set_params( const CvEMParams& _params, const CvVectors& train_data )
__END__; __END__;
} }
/****************************************************************************************/
double CvEM::calcLikelihood( const cv::Mat &input_sample ) const
{
const CvMat _input_sample = input_sample;
const CvMat* _sample = &_input_sample ;
float* sample_data = 0;
int cov_mat_type = params.cov_mat_type;
int i, dims = means->cols;
int nclusters = params.nclusters;
cvPreparePredictData( _sample, dims, 0, params.nclusters, 0, &sample_data );
// allocate memory and initializing headers for calculating
cv::AutoBuffer<double> buffer(nclusters + dims);
CvMat expo = cvMat(1, nclusters, CV_64F, &buffer[0] );
CvMat diff = cvMat(1, dims, CV_64F, &buffer[nclusters] );
// calculate the probabilities
for( int k = 0; k < nclusters; k++ )
{
const double* mean_k = (const double*)(means->data.ptr + means->step*k);
const double* w = (const double*)(inv_eigen_values->data.ptr + inv_eigen_values->step*k);
double cur = log_weight_div_det->data.db[k];
CvMat* u = cov_rotate_mats[k];
// cov = u w u' --> cov^(-1) = u w^(-1) u'
if( cov_mat_type == COV_MAT_SPHERICAL )
{
double w0 = w[0];
for( i = 0; i < dims; i++ )
{
double val = sample_data[i] - mean_k[i];
cur += val*val*w0;
}
}
else
{
for( i = 0; i < dims; i++ )
diff.data.db[i] = sample_data[i] - mean_k[i];
if( cov_mat_type == COV_MAT_GENERIC )
cvGEMM( &diff, u, 1, 0, 0, &diff, CV_GEMM_B_T );
for( i = 0; i < dims; i++ )
{
double val = diff.data.db[i];
cur += val*val*w[i];
}
}
expo.data.db[k] = cur;
}
// probability = (2*pi)^(-dims/2)*exp( -0.5 * cur )
cvConvertScale( &expo, &expo, -0.5 );
double factor = -double(dims)/2.0 * log(2.0*M_PI);
cvAndS( &expo, cvScalar(factor), &expo );
// Calculate the log-likelihood of the given sample -
// see Alex Smola's blog http://blog.smola.org/page/2 for
// details on the log-sum-exp trick
double mini,maxi,retval;
cvMinMaxLoc( &expo, &mini, &maxi, 0, 0 );
CvMat *flp = cvCloneMat(&expo);
cvSubS( &expo, cvScalar(maxi), flp);
cvExp( flp, flp );
CvScalar ss = cvSum( flp );
retval = log(ss.val[0]) + maxi;
cvReleaseMat(&flp);
if( sample_data != _sample->data.fl )
cvFree( &sample_data );
return retval;
}
/****************************************************************************************/ /****************************************************************************************/
float float
@ -219,12 +506,12 @@ CvEM::predict( const CvMat* _sample, CvMat* _probs ) const
cvPreparePredictData( _sample, dims, 0, params.nclusters, _probs, &sample_data ); cvPreparePredictData( _sample, dims, 0, params.nclusters, _probs, &sample_data );
// allocate memory and initializing headers for calculating // allocate memory and initializing headers for calculating
cv::AutoBuffer<double> buffer(nclusters + dims); cv::AutoBuffer<double> buffer(nclusters + dims);
CvMat expo = cvMat(1, nclusters, CV_64F, &buffer[0] ); CvMat expo = cvMat(1, nclusters, CV_64F, &buffer[0] );
CvMat diff = cvMat(1, dims, CV_64F, &buffer[nclusters] ); CvMat diff = cvMat(1, dims, CV_64F, &buffer[nclusters] );
// calculate the probabilities // calculate the probabilities
for( int k = 0; k < nclusters; k++ ) for( int k = 0; k < nclusters; k++ )
{ {
const double* mean_k = (const double*)(means->data.ptr + means->step*k); const double* mean_k = (const double*)(means->data.ptr + means->step*k);
@ -260,12 +547,17 @@ CvEM::predict( const CvMat* _sample, CvMat* _probs ) const
cls = k; cls = k;
opt = cur; opt = cur;
} }
/* probability = (2*pi)^(-dims/2)*exp( -0.5 * cur ) */
} }
// probability = (2*pi)^(-dims/2)*exp( -0.5 * cur )
cvConvertScale( &expo, &expo, -0.5 );
double factor = -double(dims)/2.0 * log(2.0*M_PI);
cvAndS( &expo, cvScalar(factor), &expo );
// Calculate the posterior probability of each component
// given the sample data.
if( _probs ) if( _probs )
{ {
cvConvertScale( &expo, &expo, -0.5 );
cvExp( &expo, &expo ); cvExp( &expo, &expo );
if( _probs->cols == 1 ) if( _probs->cols == 1 )
cvReshape( &expo, &expo, 0, nclusters ); cvReshape( &expo, &expo, 0, nclusters );
@ -336,6 +628,7 @@ bool CvEM::train( const CvMat* _samples, const CvMat* _sample_idx,
init_em( train_data ); init_em( train_data );
log_likelihood = run_em( train_data ); log_likelihood = run_em( train_data );
if( log_likelihood <= -DBL_MAX/10000. ) if( log_likelihood <= -DBL_MAX/10000. )
EXIT; EXIT;
@ -497,8 +790,11 @@ void CvEM::init_auto( const CvVectors& train_data )
if( nclusters > 1 ) if( nclusters > 1 )
{ {
CV_CALL( labels = cvCreateMat( 1, nsamples, CV_32SC1 )); CV_CALL( labels = cvCreateMat( 1, nsamples, CV_32SC1 ));
// Not fully executed in case means are already given
kmeans( train_data, nclusters, labels, cvTermCriteria( CV_TERMCRIT_ITER, kmeans( train_data, nclusters, labels, cvTermCriteria( CV_TERMCRIT_ITER,
params.means ? 1 : 10, 0.5 ), params.means ); params.means ? 1 : 10, 0.5 ), params.means );
CV_CALL( cvSortSamplesByClasses( (const float**)train_data.data.fl, CV_CALL( cvSortSamplesByClasses( (const float**)train_data.data.fl,
labels, class_ranges->data.i )); labels, class_ranges->data.i ));
} }
@ -855,18 +1151,18 @@ double CvEM::run_em( const CvVectors& train_data )
else else
cvTranspose( cvGetDiag( covs[k], &diag ), w ); cvTranspose( cvGetDiag( covs[k], &diag ), w );
w_data = w->data.db; w_data = w->data.db;
for( j = 0, det = 1.; j < dims; j++ ) for( j = 0, det = 0.; j < dims; j++ )
det *= w_data[j]; det += std::log(w_data[j]);
if( det < min_det_value ) if( det < std::log(min_det_value) )
{ {
if( start_step == START_AUTO_STEP ) if( start_step == START_AUTO_STEP )
det = min_det_value; det = std::log(min_det_value);
else else
EXIT; EXIT;
} }
log_det->data.db[k] = det; log_det->data.db[k] = det;
} }
else else // spherical
{ {
d = cvTrace(covs[k]).val[0]/(double)dims; d = cvTrace(covs[k]).val[0]/(double)dims;
if( d < min_variation ) if( d < min_variation )
@ -881,9 +1177,11 @@ double CvEM::run_em( const CvVectors& train_data )
} }
} }
cvLog( log_det, log_det );
if( is_spherical ) if( is_spherical )
{
cvLog( log_det, log_det );
cvScale( log_det, log_det, dims ); cvScale( log_det, log_det, dims );
}
} }
for( n = 0; n < params.term_crit.max_iter; n++ ) for( n = 0; n < params.term_crit.max_iter; n++ )
@ -952,10 +1250,13 @@ double CvEM::run_em( const CvVectors& train_data )
} }
_log_likelihood+=sum_max_val; _log_likelihood+=sum_max_val;
// check termination criteria // Check termination criteria. Use the same termination criteria as it is used in MATLAB
//if( fabs( (_log_likelihood - prev_log_likelihood) / prev_log_likelihood ) < params.term_crit.epsilon ) log_likelihood_delta = _log_likelihood - prev_log_likelihood;
if( fabs( (_log_likelihood - prev_log_likelihood) ) < params.term_crit.epsilon ) // if( n>0 )
break; // fprintf(stderr, "iter=%d, ll=%0.5f (delta=%0.5f, goal=%0.5f)\n",
// n, _log_likelihood, delta, params.term_crit.epsilon * fabs( _log_likelihood));
if ( log_likelihood_delta > 0 && log_likelihood_delta < params.term_crit.epsilon * std::fabs( _log_likelihood) )
break;
prev_log_likelihood = _log_likelihood; prev_log_likelihood = _log_likelihood;
} }
@ -1009,11 +1310,12 @@ double CvEM::run_em( const CvVectors& train_data )
} }
else else
{ {
// Det. of general NxN cov. matrix is the prod. of the eig. vals
if( is_general ) if( is_general )
cvSVD( cov, w, cov_rotate_mats[k], 0, CV_SVD_U_T ); cvSVD( cov, w, cov_rotate_mats[k], 0, CV_SVD_U_T );
cvMaxS( w, min_variation, w ); cvMaxS( w, min_variation, w );
for( j = 0, det = 1.; j < dims; j++ ) for( j = 0, det = 0.; j < dims; j++ )
det *= w_data[j]; det += std::log( w_data[j] );
log_det->data.db[k] = det; log_det->data.db[k] = det;
} }
} }
@ -1021,9 +1323,11 @@ double CvEM::run_em( const CvVectors& train_data )
cvConvertScale( weights, weights, 1./(double)nsamples, 0 ); cvConvertScale( weights, weights, 1./(double)nsamples, 0 );
cvMaxS( weights, DBL_MIN, weights ); cvMaxS( weights, DBL_MIN, weights );
cvLog( log_det, log_det );
if( is_spherical ) if( is_spherical )
{
cvLog( log_det, log_det );
cvScale( log_det, log_det, dims ); cvScale( log_det, log_det, dims );
}
} // end of iteration process } // end of iteration process
//log_weight_div_det[k] = -2*log(weights_k/det(Sigma_k))^0.5) = -2*log(weights_k) + log(det(Sigma_k))) //log_weight_div_det[k] = -2*log(weights_k/det(Sigma_k))^0.5) = -2*log(weights_k) + log(det(Sigma_k)))

Write Preview
Loading…
Cancel
Save