EM::~EM()
{
- clear();
+ //clear();
}
void EM::clear()
int nsamples = trainSamples.rows;
// Cluster samples, compute/update means
+
+ // Convert samples and means to 32F, because kmeans requires this type.
Mat trainSamplesFlt, meansFlt;
if(trainSamples.type() != CV_32FC1)
trainSamples.convertTo(trainSamplesFlt, CV_32FC1);
Mat labels;
kmeans(trainSamplesFlt, nclusters, labels, TermCriteria(TermCriteria::COUNT, means.empty() ? 10 : 1, 0.5), 10, KMEANS_PP_CENTERS, meansFlt);
+ // Convert samples and means back to 64F.
CV_Assert(meansFlt.type() == CV_32FC1);
if(trainSamples.type() != CV_64FC1)
{
// L_ik = log(weight_k) - 0.5 * log(|det(cov_k)|) - 0.5 *(x_i - mean_k)' cov_k^(-1) (x_i - mean_k)]
// q = arg(max_k(L_ik))
// probs_ik = exp(L_ik - L_iq) / (1 + sum_j!=q (exp(L_ij - L_iq))
+ // see Alex Smola's blog http://blog.smola.org/page/2 for
+ // details on the log-sum-exp trick
CV_Assert(!means.empty());
CV_Assert(sample.type() == CV_64FC1);
if(!probs && !logLikelihood)
return;
- Mat expL_Lmax(L.size(), CV_64FC1);
double maxLVal = L.at<double>(label);
+ Mat expL_Lmax = L; // exp(L_ij - L_iq)
for(int i = 0; i < L.cols; i++)
expL_Lmax.at<double>(i) = std::exp(L.at<double>(i) - maxLVal);
-
- double partSum = 0; // sum_j!=q (exp(L_ij - L_iq))
- for(int clusterIndex = 0; clusterIndex < nclusters; clusterIndex++)
- if(clusterIndex != label)
- partSum += expL_Lmax.at<double>(clusterIndex);
+ double expDiffSum = sum(expL_Lmax)[0]; // sum_j(exp(L_ij - L_iq))
if(probs)
{
probs->create(1, nclusters, CV_64FC1);
- double factor = 1./(1 + partSum);
+ double factor = 1./expDiffSum;
expL_Lmax *= factor;
expL_Lmax.copyTo(*probs);
}
if(logLikelihood)
- {
- double logWeightProbs = std::log((1 + partSum) * std::exp(maxLVal)) - 0.5 * dim * CV_LOG2PI;
- *logLikelihood = logWeightProbs;
- }
+ *logLikelihood = std::log(expDiffSum) + maxLVal - 0.5 * dim * CV_LOG2PI;
}
void EM::eStep()