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