kmeans.hh
Go to the documentation of this file.
/* -*- mia-c++ -*-
*
* This file is part of MIA - a toolbox for medical image analysis
* Copyright (c) Leipzig, Madrid 1999-2017 Gert Wollny
*
* MIA is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MIA; if not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __mia_core_kmeans_hh
#define __mia_core_kmeans_hh
#include <vector>
#include <numeric>
#include <cmath>
#include <stdexcept>
#include <iomanip>
#include <limits>
#include <mia/core/defines.hh>
#include <mia/core/errormacro.hh>
#include <mia/core/msgstream.hh>
#include <boost/concept/requires.hpp>
#include <boost/concept_check.hpp>
NS_MIA_BEGIN
int EXPORT_CORE kmeans_get_closest_clustercenter(const std::vector<double>& classes, size_t l, double value);
template <typename InputIterator, typename OutputIterator>
bool kmeans_step(InputIterator ibegin, InputIterator iend, OutputIterator obegin,
std::vector<double>& classes, size_t l, int& biggest_class )
{
cvdebug() << "kmeans enter: ";
for (size_t i = 0; i <= l; ++i )
cverb << std::setw(8) << classes[i] << " ";
cverb << "\n";
biggest_class = -1;
const double convLimit = 0.005; // currently fixed
std::vector<double> sums(classes.size());
std::vector<size_t> count(classes.size());
bool conv = false;
int iter = 50;
while ( iter-- && !conv) {
sort(classes.begin(), classes.end());
// assign closest cluster center
OutputIterator ob = obegin;
for (InputIterator b = ibegin; b != iend; ++b, ++ob) {
*ob = kmeans_get_closest_clustercenter(classes, l, *b);
++count[*ob];
sums[*ob] += *b;
};
// recompute cluster centers
conv = true;
size_t max_count = 0;
for (size_t i = 0; i <= l; i++) {
if (count[i]) {
double a = sums[i] / count[i];
if (a && fabs ((a - classes[i]) / a) > convLimit)
conv = false;
classes[i] = a;
if (max_count < count[i]) {
max_count = count[i];
biggest_class = i;
}
} else { // if a class is empty move it closer to neighbour
if (i == 0)
classes[i] = (classes[i] + classes[i + 1]) / 2.0;
else
classes[i] = (classes[i] + classes[i - 1]) / 2.0;
conv = false;
}
sums[i] = 0;
count[i] = 0;
};
};
cvinfo() << "kmeans: " << l + 1 << " classes, " << 50 - iter << " iterations";
for (size_t i = 0; i <= l; ++i )
cverb << std::setw(8) << classes[i] << " ";
cverb << "\n";
return conv;
}
template <typename InputIterator, typename OutputIterator>
bool kmeans_step_with_fixed_centers(InputIterator ibegin, InputIterator iend, OutputIterator obegin,
std::vector<double>& classes, const std::vector<bool>& fixed_center,
size_t l, int& biggest_class )
{
cvdebug() << "kmeans enter: ";
for (size_t i = 0; i <= l; ++i )
cverb << std::setw(8) << classes[i] << " ";
cverb << "\n";
biggest_class = -1;
const double convLimit = 0.005; // currently fixed
std::vector<double> sums(classes.size());
std::vector<size_t> count(classes.size());
bool conv = false;
int iter = 50;
while ( iter-- && !conv) {
sort(classes.begin(), classes.end());
// assign closest cluster center
OutputIterator ob = obegin;
for (InputIterator b = ibegin; b != iend; ++b, ++ob) {
*ob = kmeans_get_closest_clustercenter(classes, l, *b);
++count[*ob];
sums[*ob] += *b;
};
// recompute cluster centers
conv = true;
size_t max_count = 0;
for (size_t i = 0; i <= l; i++) {
if (fixed_center[i])
continue;
if (count[i]) {
double a = sums[i] / count[i];
if (a && fabs ((a - classes[i]) / a) > convLimit)
conv = false;
classes[i] = a;
if (max_count < count[i]) {
max_count = count[i];
biggest_class = i;
}
} else { // if a class is empty move it closer to neighbour
if (i == 0)
classes[i] = (classes[i] + classes[i + 1]) / 2.0;
else
classes[i] = (classes[i] + classes[i - 1]) / 2.0;
conv = false;
}
sums[i] = 0;
count[i] = 0;
};
};
cvinfo() << "kmeans: " << l + 1 << " classes, " << 50 - iter << " iterations";
for (size_t i = 0; i <= l; ++i )
cverb << std::setw(8) << classes[i] << " ";
cverb << "\n";
return conv;
}
template <typename InputIterator, typename OutputIterator>
BOOST_CONCEPT_REQUIRES( ((::boost::ForwardIterator<InputIterator>))
((::boost::Mutable_ForwardIterator<OutputIterator>)),
(void)
)
kmeans(InputIterator ibegin, InputIterator iend, OutputIterator obegin,
std::vector<double>& classes)
{
if (classes.size() < 2)
throw create_exception<std::invalid_argument>("kmeans: requested ", classes.size(),
"class(es), required are at least two");
const size_t nclusters = classes.size();
const double size = std::distance(ibegin, iend);
if ( size < nclusters )
throw create_exception<std::invalid_argument>("kmeans: insufficient input: want ", nclusters,
" classes, but git only ", size, " input elements");
double sum = std::accumulate(ibegin, iend, 0.0);
// simple initialization splitting at the mean
classes[0] = sum / (size - 1);
classes[1] = sum / (size + 1);
// first run calles directly
int biggest_class = 0;
// coverity is completely off here, the 1UL is actually a class index
// and has nothing to do with the size of the type pointed to by ibegin
//
// coverity[sizeof_mismatch]
kmeans_step(ibegin, iend, obegin, classes, 1, biggest_class);
// further clustering always splits biggest class
for (size_t l = 2; l < nclusters; l++) {
const size_t pos = biggest_class > 0 ? biggest_class - 1 : biggest_class + 1;
classes[l] = 0.5 * (classes[biggest_class] + classes[pos]);
kmeans_step(ibegin, iend, obegin, classes, l, biggest_class);
};
// some post iteration until centers no longer change
for (size_t l = 1; l < 3; l++) {
if (kmeans_step(ibegin, iend, obegin, classes, nclusters - 1, biggest_class))
break;
}
}
NS_MIA_END
#endif
fabs
double fabs(const T3DVector< T > &t)
A way to get the norm of a T3DVector using faba syntax.
Definition 3d/vector.hh:396
defines.hh
NS_MIA_BEGIN
#define NS_MIA_BEGIN
conveniance define to start the mia namespace
Definition defines.hh:33
EXPORT_CORE
#define EXPORT_CORE
Macro to manage Visual C++ style dllimport/dllexport.
Definition defines.hh:101
NS_MIA_END
#define NS_MIA_END
conveniance define to end the mia namespace
Definition defines.hh:36
errormacro.hh
cverb
#define cverb
define a shortcut to the raw output stream
Definition msgstream.hh:341
cvinfo
vstream & cvinfo()
informal output that may be of interest to understand problems with a program and are of higher prior...
Definition msgstream.hh:262
kmeans
void kmeans(InputIterator ibegin, InputIterator iend, OutputIterator obegin, std::vector< double > &classes)
Definition kmeans.hh:205
create_exception
E create_exception(T ...t)
helper template to create exceptions with complex messages
Definition errormacro.hh:80
kmeans_step_with_fixed_centers
bool kmeans_step_with_fixed_centers(InputIterator ibegin, InputIterator iend, OutputIterator obegin, std::vector< double > &classes, const std::vector< bool > &fixed_center, size_t l, int &biggest_class)
Definition kmeans.hh:111
kmeans_step
bool kmeans_step(InputIterator ibegin, InputIterator iend, OutputIterator obegin, std::vector< double > &classes, size_t l, int &biggest_class)
Definition kmeans.hh:41
kmeans_get_closest_clustercenter
int EXPORT_CORE kmeans_get_closest_clustercenter(const std::vector< double > &classes, size_t l, double value)
msgstream.hh
cvdebug
CDebugSink & cvdebug()
Definition msgstream.hh:226
std
Definition gsl_iterator.hh:324