/** @file ikmeans_lloyd.tc
 ** @brief Integer K-Means - LLoyd Algorithm - Definition
 ** @author Brian Fulkerson
 ** @author Andrea Vedaldi
 **/

/*
Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
All rights reserved.

This file is part of the VLFeat library and is made available under
the terms of the BSD license (see the COPYING file).
*/

/** @internal
 ** @brief Helper function to initialize a filter for Lloyd algorithm
 **
 ** @param f filter.
 **/

static void
vl_ikm_init_lloyd (VlIKMFilt * f VL_UNUSED)
{ }

/** @internal
 ** @brief LLoyd algorithm
 ** @param f     IKM quantizer.
 ** @param data  Training data.
 ** @param N     Number of traning data.
 **/

static int
vl_ikm_train_lloyd (VlIKMFilt* f, vl_uint8 const* data, int N)
{
  int iter, i, j, k,
    K = f-> K,
    M = f-> M,
    err =  0 ;
  vl_ikm_acc *centers = f->centers ;
  vl_uint    *asgn    = vl_malloc (sizeof(vl_uint) * N) ;
  vl_uint    *counts  = vl_malloc (sizeof(vl_uint) * N) ;

  for (iter = 0 ; 1 ; ++ iter) {
    vl_bool done = 1 ;

    /* ---------------------------------------------------------------
     *                                               Calc. assignments
     * ------------------------------------------------------------ */

    for (j = 0 ; j < N ; ++j) {
      vl_ikm_acc best_dist = 0 ;
      vl_uint    best = (vl_uint)-1 ;

      for (k = 0; k < K; ++k) {
        vl_ikm_acc dist = 0;

        /* compute distance with this center */
        for (i = 0; i < M; ++i) {
          vl_ikm_acc delta =
            data [j * M + i] - centers [k * M + i] ;
          dist += delta * delta ;
        }

        /* compare with current best */
        if (best == (vl_uint) -1 || dist < best_dist) {
          best = k ;
          best_dist = dist ;
        }
      }
      if (asgn [j] != best) {
        asgn [j] = best ;
        done = 0 ;
      }
    }

    /* stopping condition */
    if (done || iter == f->max_niters) break ;

    /* -----------------------------------------------------------------
     *                                                     Calc. centers
     * -------------------------------------------------------------- */

    /* re-compute centers */
    memset (centers, 0, sizeof (vl_int32) * M * K);
    memset (counts,  0, sizeof (vl_int32) * K);
    for (j = 0; j < N; ++j) {
      int this_center = asgn [j] ;
      ++ counts [this_center] ;
      for (i = 0; i < M; ++i)
        centers [this_center * M + i] += data[j * M + i];
    }

    for (k = 0; k < K; ++k) {
      vl_int32 n = counts [k];
      if (n > 0xffffff) {
        err = 1 ;
      }
      if (n > 0) {
        for (i = 0; i < M; ++i)
          centers [k * M + i] /= n;
      } else {
        /*
           If no data are assigned to the center, it is not changed
           with respect to the previous iteration, so we do not do anything.
        */
      }
    }
  }

  vl_free (counts) ;
  vl_free (asgn) ;
  return err ;
}


/** @internal
 ** @brief LLoyd algorithm
 ** @param f     IKM quantizer.
 ** @param asgn  Assignments (out).
 ** @param data  Data to quantize.
 ** @param N     Number of data.
 **/

static void
vl_ikm_push_lloyd (VlIKMFilt *f, vl_uint *asgn, vl_uint8 const *data, int N)
{
  int j ;
  for(j=0 ; j < N ; ++j) {
    asgn[j] = vl_ikm_push_one (f->centers, data + j * f->M, f->M, f->K);
  }
}

/*
 * Local Variables: *
 * mode: C *
 * End: *
 */