}
}
+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 )
{
__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
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);
CvMat expo = cvMat(1, nclusters, CV_64F, &buffer[0] );
CvMat diff = cvMat(1, dims, CV_64F, &buffer[nclusters] );
-// calculate the probabilities
+ // calculate the probabilities
for( int k = 0; k < nclusters; k++ )
{
const double* mean_k = (const double*)(means->data.ptr + means->step*k);
cls = k;
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 )
{
- cvConvertScale( &expo, &expo, -0.5 );
cvExp( &expo, &expo );
if( _probs->cols == 1 )
cvReshape( &expo, &expo, 0, nclusters );
init_em( train_data );
log_likelihood = run_em( train_data );
+
if( log_likelihood <= -DBL_MAX/10000. )
EXIT;
if( nclusters > 1 )
{
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,
params.means ? 1 : 10, 0.5 ), params.means );
+
CV_CALL( cvSortSamplesByClasses( (const float**)train_data.data.fl,
labels, class_ranges->data.i ));
}
else
cvTranspose( cvGetDiag( covs[k], &diag ), w );
w_data = w->data.db;
- for( j = 0, det = 1.; j < dims; j++ )
- det *= w_data[j];
- if( det < min_det_value )
+ for( j = 0, det = 0.; j < dims; j++ )
+ det += std::log(w_data[j]);
+ if( det < std::log(min_det_value) )
{
if( start_step == START_AUTO_STEP )
- det = min_det_value;
+ det = std::log(min_det_value);
else
EXIT;
}
log_det->data.db[k] = det;
}
- else
+ else // spherical
{
d = cvTrace(covs[k]).val[0]/(double)dims;
if( d < min_variation )
}
}
- cvLog( log_det, log_det );
if( is_spherical )
+ {
+ cvLog( log_det, log_det );
cvScale( log_det, log_det, dims );
+ }
}
for( n = 0; n < params.term_crit.max_iter; n++ )
}
_log_likelihood+=sum_max_val;
- // check termination criteria
- //if( fabs( (_log_likelihood - prev_log_likelihood) / prev_log_likelihood ) < params.term_crit.epsilon )
- if( fabs( (_log_likelihood - prev_log_likelihood) ) < params.term_crit.epsilon )
- break;
+ // Check termination criteria. Use the same termination criteria as it is used in MATLAB
+ log_likelihood_delta = _log_likelihood - prev_log_likelihood;
+// if( n>0 )
+// 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;
}
}
else
{
+ // Det. of general NxN cov. matrix is the prod. of the eig. vals
if( is_general )
cvSVD( cov, w, cov_rotate_mats[k], 0, CV_SVD_U_T );
cvMaxS( w, min_variation, w );
- for( j = 0, det = 1.; j < dims; j++ )
- det *= w_data[j];
+ for( j = 0, det = 0.; j < dims; j++ )
+ det += std::log( w_data[j] );
log_det->data.db[k] = det;
}
}
cvConvertScale( weights, weights, 1./(double)nsamples, 0 );
cvMaxS( weights, DBL_MIN, weights );
- cvLog( log_det, log_det );
if( is_spherical )
+ {
+ cvLog( log_det, log_det );
cvScale( log_det, log_det, dims );
+ }
} // 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)))
projects / profile / ivi / opencv.git / commitdiff