move.centroids <- function (cb, Data, pch = "+", Dim=2) {
      ##  Repeat this step until little change in cb pos
   shift <- 100
   j <- 2
   while (shift > .003) {
      n.cb <- dim(cb)[1]
         # min.index gives the number of the nearest codebook vector to the data points
      Dist <- Array2ArrayDist(cb, Data)
      min.dist  <- apply(Dist, 2, min)
      if (n.cb == 1) {
         min.index <- apply(Dist, 2, order)
      } else {
         min.index <- apply(Dist, 2, order)[1,]
      }
      j <- j + 1
      if (Dim==2) {
         Voronoi(cb, 1:n.cb, c(200,200), range(Data[,1]), range(Data[,2]))
         points(Data, col=Col[min.index], pch=pch)
         points(cb, col="red", pch=16, cex=1.5)
      } else if (Dim==3) {
         plot3d(Data, size=3, col=Col[min.index])
         spheres3d(cb, col="red", size=9)
      }
         #  Move the codebook vectors to the middle of the group
      mid.point <- matrix(0, n.cb, Dim)
      for (c in 1:n.cb) {
         mid.point[c,] <- apply(Data[(min.index==c),], 2, mean)
      }
      if (Dim==2) {
         points(mid.point, col="green", pch=17, cex=1.5)
      } else if (Dim==3) {
         spheres3d(mid.point, col="green",size=9)
#         s3d$points3d(mid.point, col="green", pch=17, cex=1.5)
      }
      shift <- max(mid.point - cb)                       # Could compare distortion
      cb <- mid.point
   }
   return(list(cb, min.index))
}

compute.distortion <- function(cb, Data, belongs.to) {
   N <- dim(cb)[1]
   D <- rep(0, N)
   for (i in 1:N) {
      M <- sum(belongs.to==i)
      D[i] <- sum(apply((t(t(Data[belongs.to==i,])-cb[i,]))^2, 1, sum))/M
   }
   D
}

get.pairs <- function (numb) {
      # All unique pairs
   n.comb.2 <- choose(numb,2)
   double <- matrix(0, n.comb.2, 2)
   ind <- 1
   for (i in 1:numb) {
      j <- i+1
      while (j <= numb) {  # for can run 6:5 ...
         double[ind,] <- c(i,j)
         ind <- ind + 1
         j <- j + 1
      }
   }
   double
}

Davies.Bouldin <- function(A, SS, m) {
      # A  - the centres of the clusters
      # SS - the within sum of squares
      # m  - the sizes of the clusters
   N <- nrow(A)   # number of clusters
      # intercluster distance
   S <- sqrt(SS/m)
      # Get the distances between centres
   M <- as.matrix(dist(A))
      # Get the ratio of intercluster/centre.dist
   R <- matrix(0, N, N)
   for (i in 1:(N-1)) {
      for (j in (i+1):N) {
         R[i,j] <- (S[i] + S[j])/M[i,j]
         R[j,i] <- R[i,j]
      }
   }
   return(mean(apply(R, 1, max)))
}

deconstruct.4 <- function(pix) {
   Rows <- dim(pix)[1]
   Cols <- dim(pix)[2]
      # The data will be a 2x2 block of the image.
      # This will be used as the data for k-means
      # and we will try to find clusters of similar values.
   blocked <- matrix(0, (Rows/2)*(Cols/2), 4)

   i <- 1
   for (row in seq(1, Rows, by=2)) {
      for (col in seq(1, Cols, by=2)) {
            # Take a 2x2 block and make it a vector
         blocked[i,] <- as.vector(pix[row:(row+1),col:(col+1)])
         i <- i+1
      }
   }
   return(blocked)
}

deconstruct.16 <- function(pix) {
   Rows <- dim(pix)[1]
   Cols <- dim(pix)[2]
      # The data will be a 4x4 block of the image.
      # This will be used as the data for k-means
      # and we will try to find clusters of similar values.
   blocked <- matrix(0, (Rows/4)*(Cols/4), 16)
   i <- 1
   for (row in seq(1, Rows, by=4)) {
      for (col in seq(1, Cols, by=4)) {
            # Take a 4x4 block and make it a vector
         blocked[i,] <- as.vector(pix[row:(row+3),col:(col+3)])
         i <- i+1
      }
   }
   return(blocked)
}   

reconstruct.4 <- function(new.bytes, Rows, Cols){
   new.pic <- matrix(0, Rows, Cols)
   i <- 1
   for (row in seq(1, Rows, by=2)) {
      for (col in seq(1, Cols, by=2)) {
         new.pic [row,col]     <- new.bytes[i,1]
         new.pic [row+1,col]   <- new.bytes[i,2]
         new.pic [row,col+1]   <- new.bytes[i,3]
         new.pic [row+1,col+1] <- new.bytes[i,4]
         i <- i+1
      }
   }
   return (new.pic)
}

reconstruct.16 <- function(new.bytes, Rows, Cols){
   new.pic <- matrix(0, Rows, Cols)
   i <- 1
   for (row in seq(1, Rows, by=4)) {
      for (col in seq(1, Cols, by=4)) {
         new.pic [row,col]       <- new.bytes[i,1]
         new.pic [row+1,col]     <- new.bytes[i,2]
         new.pic [row+2,col]     <- new.bytes[i,3]
         new.pic [row+3,col]     <- new.bytes[i,4]
         new.pic [row,col+1]     <- new.bytes[i,5]
         new.pic [row+1,col+1]   <- new.bytes[i,6]
         new.pic [row+2,col+1]   <- new.bytes[i,7]
         new.pic [row+3,col+1]   <- new.bytes[i,8]
         new.pic [row,col+2]     <- new.bytes[i,9]
         new.pic [row+1,col+2]   <- new.bytes[i,10]
         new.pic [row+2,col+2]   <- new.bytes[i,11]
         new.pic [row+3,col+2]   <- new.bytes[i,12]
         new.pic [row,col+3]     <- new.bytes[i,13]
         new.pic [row+1,col+3]   <- new.bytes[i,14]
         new.pic [row+2,col+3]   <- new.bytes[i,15]
         new.pic [row+3,col+3]   <- new.bytes[i,16]
         i <- i+1
      }
   }
   return (new.pic)
}