k-means clustering

k-means clustering reconstitutes original image by grouping homogenious field of the image in the same category. k is set by user. it corresponds to the number of color available in the original image.

im <- readImage("basales.jpg")
dim(im)

## [1] 1024  768    3

plot(im) # raster method means within R

# reshape image into a data frame
df = data.frame(
  red = matrix(im[,,1], ncol=1),
  green = matrix(im[,,2], ncol=1),
  blue = matrix(im[,,3], ncol=1)
)
str(df)

## 'data.frame':    786432 obs. of  3 variables:
##  $ red  : num  0.831 0.835 0.851 0.847 0.835 ...
##  $ green: num  0.792 0.792 0.792 0.796 0.8 ...
##  $ blue : num  0.784 0.784 0.78 0.773 0.773 ...

### compute the k-means clustering
set.seed(1234)
K = kmeans(df,5)
str(K)

## List of 9
##  $ cluster     : int [1:786432] 3 3 3 3 3 3 3 3 3 3 ...
##  $ centers     : num [1:5, 1:3] 0.7 0.828 0.87 0.854 0.68 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:5] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:3] "red" "green" "blue"
##  $ totss       : num 27599
##  $ withinss    : num [1:5] 519 936 938 642 949
##  $ tot.withinss: num 3983
##  $ betweenss   : num 23615
##  $ size        : int [1:5] 59929 254563 170048 148560 153332
##  $ iter        : int 6
##  $ ifault      : int 0
##  - attr(*, "class")= chr "kmeans"

K$centers

##         red     green      blue
## 1 0.7003202 0.5972739 0.7018028
## 2 0.8280036 0.4978339 0.5741294
## 3 0.8700125 0.7728880 0.7826293
## 4 0.8542437 0.6200428 0.6987729
## 5 0.6801946 0.4191299 0.5305071

df$label = K$cluster
dim(df)

## [1] 786432      4

head(df)

compute thresholds of intensities of each color (R,G,B) for each category

### Replace the color of each pixel in the image with the mean 
### R,G, and B values of the cluster in which the pixel resides:

# get the coloring
colors = data.frame(
  label = 1:nrow(K$centers), 
  R = K$centers[,"red"],
  G = K$centers[,"green"],
  B = K$centers[,"blue"]
)
dim(colors)

## [1] 5 4

colors

# merge color codes on to df
# IMPORTANT: we must maintain the original order of the df after the merge!

df$order <- 1:nrow(df)
df <- merge(df, colors)
## split image into two dataframe df.x contain 1 label, df.y contain the remain labels

df.x <- df[df$label == 1,]
df.y <- df[!df$label == 1,]
## replace value in df.y  to 255 (white)
df.y[,c('red','green', 'blue', 'R', 'G', 'B')] <- 255

df <- rbind(df.x, df.y)
# reorder the matrix (image)
df <- df[order(df$order),]
df$order = NULL

head(df)

dim(df)

## [1] 786432      7

## attribute white color to class 5, 4, 3, 2
# df[df$label == 5, -1] <- 255
# df[df$label == 4, -1] <- 255
# df[df$label == 3, -1] <- 255
# df[df$label == 2, -1] <- 255

Reshape our data frame back into an image

# get mean color channel values for each row of the df.
R <- matrix(df$R, nrow=dim(im)[1])
G <- matrix(df$G, nrow=dim(im)[1])
B <- matrix(df$B, nrow=dim(im)[1])

# reconstitute the segmented image in the same shape as the input image
im.segmented <- array(dim=dim(im))
im.segmented[,,1] = R
im.segmented[,,2] = G
im.segmented[,,3] = B
#im.segmented <- rotate(im.segmented, 90)
#im.segmented <- flop(im.segmented )
#dim(im.segmented)
#str(im.segmented)
#class(im.segmented)
#EBImage::writeImage(im.segmented, files = "im.segmented.jpg")
# View the result
#layout(t(1:2))

## convert matrix to Image
im.segmented <- Image(im.segmented, colormode=Color)

class(im.segmented)

## [1] "Image"
## attr(,"package")
## [1] "EBImage"

class(im)

## [1] "Image"
## attr(,"package")
## [1] "EBImage"

plot(im.segmented)

plot(im)

dim(im.segmented)

## [1] 1024  768    3

dim(im)

## [1] 1024  768    3

Function to plot and compare segmented and original image

get_kmean_image <- function(file, k){
  
  im <- EBImage::readImage(file)
  # reshape image into a data frame
  df = data.frame(
    red = matrix(im[,,1], ncol=1),
    green = matrix(im[,,2], ncol=1),
    blue = matrix(im[,,3], ncol=1)
  )
  
  ### compute the k-means clustering
  set.seed(1234)
  K = kmeans(df,k)
  df$label = K$cluster
  
  
  
  ### Replace the color of each pixel in the image with the mean 
  ### R,G, and B values of the cluster in which the pixel resides:
  
  # get the coloring
  colors = data.frame(
    label = 1:nrow(K$centers), 
    R = K$centers[,"red"],
    G = K$centers[,"green"],
    B = K$centers[,"blue"]
  )
  
  # merge color codes on to df
  # IMPORTANT: we must maintain the original order of the df after the merge!
  
  df$order <- 1:nrow(df)
  df <- merge(df, colors)
  ## split image into two dataframe df.x contain 1 label, df.y contain the remain labels
  #if(cat %in% seq(k)){
  
  #if(cat == "all"){
  
  
  M <- vector("list", k)
  for(i in seq(k)){
    
    ## preserve original df
    m <- df
    # m$order <- 1:nrow(m)
    
    m.x <- m[m$label == i,]
    m.y <- m[!m$label == i,]
    
    ## replace value in df.y  to 255 (white)
    m.y[,c('red','green', 'blue', 'R', 'G', 'B')] <- 255
    m <- rbind(m.x, m.y)
    
    # reorder the matrix (image)
    m <- m[order(m$order),]
    m$order = NULL
    
    
    
    # df.x <- df[df$label == cat,]
    #df.y <- df[!df$label == cat,]
    ## replace value in df.y  to 255 (white)
    #df.y[,c('red','green', 'blue', 'R', 'G', 'B')] <- 255
    
    #df <- rbind(df.x, df.y)
    #}
    # reorder the matrix (image)
    #df <- df[order(df$order),]
    #df$order = NULL
    
    # get mean color channel values for each row of the df.
    R <- matrix(m$R, nrow=dim(im)[1])
    G <- matrix(m$G, nrow=dim(im)[1])
    B <- matrix(m$B, nrow=dim(im)[1])
    
    # reconstitute the segmented image in the same shape as the input image
    im.segmented <- array(dim=dim(im))
    im.segmented[,,1] = R
    im.segmented[,,2] = G
    im.segmented[,,3] = B
    
    # convert array to cimg for imager package from dim(x,y,channel) to dim(x, y, 1, channels)
    # ignore warning message if necessary
    #suppressWarnings(im.segmented_cimg <- imager::as.cimg(im.segmented))
    M[[i]] <- EBImage::Image(im.segmented, colormode=Color)
    
    im.segmented <- EBImage::Image(im.segmented, colormode=Color)
    plot(im.segmented, title= paste0("Cluster: " ,i, sep= "~"))
    text(470, -50, paste0("Cluster: " ,i, sep= "~"), cex = 1.5)
    #plot(im, title= "original image")
    #text(500, -50, " Original image", cex = 1.5)
    #  }
  }
  return(M)
}



par(mfrow=c(2,3))
plot(im)
#M <- get_kmean_image(file = "basales.jpg", k = 5)
get_kmean_image(file = "basales.jpg", k = 5)

## [[1]]
## Image 
##   colorMode    : Color 
##   storage.mode : double 
##   dim          : 1024 768 3 
##   frames.total : 3 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6,1]
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]  255  255  255  255  255  255
## [2,]  255  255  255  255  255  255
## [3,]  255  255  255  255  255  255
## [4,]  255  255  255  255  255  255
## [5,]  255  255  255  255  255  255
## 
## [[2]]
## Image 
##   colorMode    : Color 
##   storage.mode : double 
##   dim          : 1024 768 3 
##   frames.total : 3 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6,1]
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]  255  255  255  255  255  255
## [2,]  255  255  255  255  255  255
## [3,]  255  255  255  255  255  255
## [4,]  255  255  255  255  255  255
## [5,]  255  255  255  255  255  255
## 
## [[3]]
## Image 
##   colorMode    : Color 
##   storage.mode : double 
##   dim          : 1024 768 3 
##   frames.total : 3 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6,1]
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
## [1,] 0.8700125 0.8700125 0.8700125 0.8700125 0.8700125 0.8700125
## [2,] 0.8700125 0.8700125 0.8700125 0.8700125 0.8700125 0.8700125
## [3,] 0.8700125 0.8700125 0.8700125 0.8700125 0.8700125 0.8700125
## [4,] 0.8700125 0.8700125 0.8700125 0.8700125 0.8700125 0.8700125
## [5,] 0.8700125 0.8700125 0.8700125 0.8700125 0.8700125 0.8700125
## 
## [[4]]
## Image 
##   colorMode    : Color 
##   storage.mode : double 
##   dim          : 1024 768 3 
##   frames.total : 3 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6,1]
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]  255  255  255  255  255  255
## [2,]  255  255  255  255  255  255
## [3,]  255  255  255  255  255  255
## [4,]  255  255  255  255  255  255
## [5,]  255  255  255  255  255  255
## 
## [[5]]
## Image 
##   colorMode    : Color 
##   storage.mode : double 
##   dim          : 1024 768 3 
##   frames.total : 3 
##   frames.render: 1 
## 
## imageData(object)[1:5,1:6,1]
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]  255  255  255  255  255  255
## [2,]  255  255  255  255  255  255
## [3,]  255  255  255  255  255  255
## [4,]  255  255  255  255  255  255
## [5,]  255  255  255  255  255  255

#get_kmean_image(file = "basales.jpg", k = 5)
#get_kmean_image(file = "basales.jpg", k = 5)
#get_kmean_image(file = "basales.jpg", k = 5)

other k-means clustering with 2D projection

# segment_image = function(img, n){
#   # create a flat, segmented image data set using kmeans
#   # Segment an RGB image into n groups based on color values using Kmeans
#   df = data.frame(
#     red = matrix(img[,,1], ncol=1),
#     green = matrix(img[,,2], ncol=1),
#     blue = matrix(img[,,3], ncol=1)
#   )
#   K = kmeans(df,n)
#   df$label = K$cluster
#   
#   # compute rgb values and color codes based on Kmeans centers
#   colors = data.frame(
#     label = 1:nrow(K$centers), 
#     R = K$centers[,"red"],
#     G = K$centers[,"green"],
#     B = K$centers[,"blue"],
#     color=rgb(K$centers)
#   )
#   
#   # merge color codes on to df but maintain the original order of df
#   df$order = 1:nrow(df)
#   df = merge(df, colors)
#   df = df[order(df$order),]
#   df$order = NULL
#   
#   return(df)
#   
# }
# 
# #
# # reconstitue the segmented images to RGB matrix
# #
# build_segmented_image = function(df, img){
#   # reconstitue the segmented images to RGB array
#   
#   # get mean color channel values for each row of the df.
#   R = matrix(df$R, nrow=dim(img)[1])
#   G = matrix(df$G, nrow=dim(img)[1])
#   B = matrix(df$B, nrow=dim(img)[1])
#   
#   # reconsitute the segmented image in the same shape as the input image
#   img_segmented = array(dim=dim(img))
#   dim(img_segmented)
#   img_segmented[,,1] = R
#   img_segmented[,,2] = G
#   img_segmented[,,3] = B
#   
#   return(img_segmented)
# }
# 
# #
# # 2D projection for visualizing the kmeans clustering
# #
# project2D_from_RGB = function(df){
#   # Compute the projection of the RGB channels into 2D
#   PCA = prcomp(df[,c("red","green","blue")], center=TRUE, scale=TRUE)
#   pc2 = PCA$x[,1:2]
#   df$x = pc2[,1]
#   df$y = pc2[,2]
#   return(df[,c("x","y","label","R","G","B", "color")])
# }
# 
# # #
# # # Create the projection plot of the clustered segments
# # #
# # plot_projection <- function(df, sample.size){
# #   # plot the projection of the segmented image data in 2D, using the
# #   # mean segment colors as the colors for the points in the projection
# #   index = sample(1:nrow(df), sample.size)
# #   return(ggplot(df[index,], aes(x=x, y=y, col=color)) + geom_point(size=2) + scale_color_identity())
# # }
# 
# # #
# # # Inspect
# # #
# # inspect_segmentation <- function(image.raw, image.segmented, image.proj){
# #   # helper function to review the results of segmentation visually
# #   img1 = rasterGrob(image.raw)
# #   img2 = rasterGrob(image.segmented)
# #   plt = plot_projection(image.proj, 50000)
# #   grid.arrange(arrangeGrob(img1,img2, nrow=1),plt)
# # }
# 
# 
# get_kmean_image2D <- function(file, k){
#   require(jpeg)
#   img <- readJPEG(file)
#   df <- segment_image(img, k)
#   
#   img_segmented <- build_segmented_image(df, img)
#   
#   ## convert matrix to Image
#   ## rotate and flop to get the same orientation as in the last k-mean code
#   img <- flop(rotate(Image(img, colormode=Color), 90))
#   img_segmented <- flop(rotate(Image(img_segmented, colormode=Color), 90))
#   
#   
#   par(mfrow=c(1,2))
#   plot(img_segmented, title= "reconstitute image using k-mean clustering")
#   text(470, -50, "k-means reconstituted image ", cex = 1.5)
#   plot(img, title= "original image")
#   text(500, -50, " Original image", cex = 1.5)
# }
# 
# # get_kmean_image2D("basales.jpg", 5)
# 
# 
# ## plot 3D image
# # library(rgl)
# # library(jpeg)
# # img <- readJPEG("basales.jpg")
# # img_kmeans_df <- segment_image(img,6)
# # img_project2D <- project2D_from_RGB(img_kmeans_df)
# # img_project2D <- Image(img_Proj2D_df, colormode = Color)
# # rgl::plot3d(img_Proj2D_df, col = rainbow(1000))

Highlight cellular components using k thresholds

im <- readImage("basales.jpg")
# reshape image into a data frame
df = data.frame(
  red = matrix(im[,,1], ncol=1),
  green = matrix(im[,,2], ncol=1),
  blue = matrix(im[,,3], ncol=1)
)

### compute the k-means clustering
K <- kmeans(df,5)

## Warning: Quick-TRANSfer stage steps exceeded maximum (= 39321600)

df$label <- K$cluster

### Replace the color of each pixel in the image with the mean 
### R,G, and B values of the cluster in which the pixel resides:

# get the coloring
colors <- data.frame(
  label = 1:nrow(K$centers), 
  R = K$centers[,"red"],
  G = K$centers[,"green"],
  B = K$centers[,"blue"]
)
colors

# convert array to cimg for imager package from dim(x,y,channel) to dim(x, y, 1, channels)
# ignore warning message
suppressWarnings(im.segmented_cimg <- imager::as.cimg(im.segmented))
# Splitting image by color)
img.R <- imager::imsplit(im.segmented_cimg,"c")[1]

scale1 <- function(r,g,b) rgb( r< 0.7006928, g < 0.5978856 && g >0.4191299, b < 0.7021399 && b > 0.6981221)
scale2 <- function(r,g,b) rgb( r < 0.6800695, g < 0.4191247 , b < 0.5305520)
scale3 <- function(r,g,b) rgb( r> 0.8699541, g > 0.7725557 , b > 0.7824795)
scale4 <- function(r,g,b) rgb( r> 0.8543942, g < 0.6191575 && g>0.5978856 , b < 0.6981221 && b > 0.5738834)
scale5 <- function(r,g,b) rgb( r> 0.8278477, g < 0.4976306 && g>0.4191247 , b < 0.5738834 && b > 0.5305520)



par(mfrow=c(2,2))
#plot(im)
plot(im.segmented)
plot(im.segmented_cimg > .7, axes = FALSE, rescale = TRUE)
plot(im.segmented_cimg,colourscale= scale1,rescale=FALSE, axe =TRUE)
plot(im.segmented_cimg,colourscale= scale2,rescale=FALSE, axe =TRUE)

plot(im.segmented_cimg,colourscale= scale3,rescale=TRUE, axe =FALSE)

## Warning in as.raster.cimg(im, rescale = rescale, colorscale = colorscale, :
## You've specified a colour scale, but rescale is set to TRUE. You may get
## unexpected results

plot(im.segmented_cimg,colourscale= scale4,rescale=TRUE, axe =FALSE)

## Warning in as.raster.cimg(im, rescale = rescale, colorscale = colorscale, :
## You've specified a colour scale, but rescale is set to TRUE. You may get
## unexpected results

plot(im.segmented_cimg,colourscale= scale4,rescale=TRUE, axe =FALSE)

## Warning in as.raster.cimg(im, rescale = rescale, colorscale = colorscale, :
## You've specified a colour scale, but rescale is set to TRUE. You may get
## unexpected results

#plot(im.segmented_cimg,colourscale= cyto.scale,rescale=TRUE, axe =FALSE)

# im.segmented_cimg_slot <-  grayscale(readImage("basales.jpg"))
# plot(im.segmented_cimg_slot)
# im.segmented_cimg_gblur <- gblur(im.segmented_cimg_slot, sigma = 5)
# im.blur.thr.cnt <- bwlabel(im.segmented_cimg_gblur > otsu(im.segmented_cimg_gblur))
# N <- max(im.blur.thr.cnt)
# N

opening, propagate and makeBrush

source("../R/get_kseg_image.R")
Result <- get_kseg_image("papSmear.jpg", k= 4)
cyto <- Result$seg_image$cluster3
nuc <- Result$seg_image$cluster2

nmask <- thresh(cyto, w=2, h=2, offset=0.05)
## makeBrush rempli le vide dans un cytoplame
nmask <- opening(nmask, makeBrush(3, shape='disc'))
#nmask <- fillHull(nmask)
#nmask <- bwlabel(nmask)

ctmask <- opening(cyto, makeBrush(3, shape='disc'))
cmask <- propagate(cyto, seeds=nmask, mask=ctmask)
display(cmask/max(cmask), all=TRUE, title="Figure 19: cmask/max(cmask)")

clean nuclear from cellular membrane

source("../R/get_kseg_image.R")
Result <- get_kseg_image("papSmear.jpg", k= 4)
cyto <- Result$seg_image$cluster4
nuc <- Result$seg_image$cluster2

plot(nuc)

We would like to remove non nuclear material

## not rigth for cytoplasm
#cyto.blur <- EBImage::gblur(cyto, sigma = 0.4)
nuc.blur <- EBImage::gblur(nuc, sigma = 0.4)
plot(nuc.blur)

count nucleus does not work

plot(nuc.blur)

nuc.blur.gray <- channel(nuc.blur, mode = "gray")
nuc.blur.gray.cnt <- EBImage::bwlabel(nuc.blur.gray)
max(nuc.blur.gray.cnt > 0)

## [1] 1

plot(nuc.blur.gray)

max(nuc.blur.gray.cnt)

## [1] 1

merge cyto + nuc to count cells does not work

cytoNuc <- merge2img(cyto, nuc)
plot(cytoNuc)
cytoNuc.gray <- channel(cytoNuc, mode = "gray")
plot(cytoNuc.gray)

cytoNuc.gray.gblur <- EBImage:: gblur(cytoNuc.gray, sigma = 0.4)
plot(cytoNuc.gray.gblur)

threshold = otsu(cytoNuc.gray.gblur)
cytoNuc.gary.gblur.cnt <- EBImage::bwlabel(cytoNuc.gray.gblur)
max(cytoNuc.gary.gblur.cnt > threshold)

## [1] 1

display(colorLabels(bwlabel(cytoNuc.gray.gblur> .26)), method = "raster", all = TRUE)

get_kseg_image <- function(file, k){
  
  im <- EBImage::readImage(file)
  # reshape image into a data frame
  df <- data.frame(
    red = matrix(im[,,1], ncol=1),
    green = matrix(im[,,2], ncol=1),
    blue = matrix(im[,,3], ncol=1)
  )
  
  ### compute the k-means clustering
  set.seed(1234)
  K = kmeans(df,k)
  df$label = K$cluster
  
  ### Replace the color of each pixel in the image with the mean 
  ### R,G, and B values of the cluster in which the pixel resides:
  
  # get the coloring
  colors = data.frame(
    label = 1:nrow(K$centers), 
    R = K$centers[,"red"],
    G = K$centers[,"green"],
    B = K$centers[,"blue"]
  )
  
  # merge color codes on to df
  # IMPORTANT: we must maintain the original order of the df after the merge!
  
  df$order <- 1:nrow(df)
  df <- merge(df, colors)
  
  ## define a list of element where we will store segmented matrices
  M <- vector("list", k)
  IM <- vector("list", k)
  for(i in seq(k)){
    
    m <- df
    
    ## split image into two dataframe df.x contain 1 label, df.y contain the remain labels
    m.x <- m[m$label == i,]
    m.y <- m[!m$label == i,]
    
    ## replace value in df.y  to 255 (white)
    m.y[,c('red','green', 'blue', 'R', 'G', 'B')] <- 255
    m <- rbind(m.x, m.y)
    
    # reorder the matrix (image)
    m <- m[order(m$order),]
    m$order = NULL
    
    # get mean color channel values for each row of the df.
    R <- matrix(m$R, nrow=dim(im)[1])
    G <- matrix(m$G, nrow=dim(im)[1])
    B <- matrix(m$B, nrow=dim(im)[1])
    
    # reconstitute the segmented image in the same shape as the input image
    im.segmented <- array(dim=dim(im))
    im.segmented[,,1] = R
    im.segmented[,,2] = G
    im.segmented[,,3] = B
    
    IM[[paste0("cluster",i, sep="")]] <- im.segmented
    # convert array to cimg for imager package from dim(x,y,channel) to dim(x, y, 1, channels)
    # ignore warning message if necessary
    #suppressWarnings(im.segmented_cimg <- imager::as.cimg(im.segmented))
    M[[paste0("cluster",i, sep="")]] <- EBImage::Image(im.segmented, colormode=Color)
    
    im.segmented <- EBImage::Image(im.segmented, colormode=Color)
    plot(im.segmented, title= paste0("Cluster: " ,i, sep= "~"))
    text(60, 5, paste0("Cluster: " ,i, sep= "~"), cex = 1.5)
  }
  return(list(seg_array = IM, seg_image = M))
}

sum array does not work

Result <- get_kseg_image("papSmear.jpg", k= 4)

cyto <- Result$seg_image$cluster4
nuc <- Result$seg_image$cluster2
cytoNuc <- cyto - nuc 
plot(cytoNuc)

cytoNuc_gray <- channel(cytoNuc, mode = "gray")
cytoNuc_gray_gblur <- EBImage:: gblur(cytoNuc_gray, sigma = 0.4)
plot(cytoNuc_gray_gblur)

#threshold = otsu(cytoNuc_gray_gblur)
cytoNuc_gray_gblur_cnt <- EBImage::bwlabel(cytoNuc_gray_gblur)
max(cytoNuc_gray_gblur_cnt)

## [1] 7