Scope

Ever since I learned about clustering methods for colors, I wanted to create something visually striking and eye-catching related to paintings. I had heard of paint-by-numbers kits before, and the idea of generating them from photos I choose caught my interest.

For inspiration, I found this youtube video. The goal is to transform a given photo into a paint-by-numbers canvas by dividing it into small, numbered sections. Each number corresponds to a color from a predefined palette, chosen to best match the original image.

To achieve this, I use the CLARA clustering algorithm to determine the most representative colors. Then, I apply the SLIC (Simple Linear Iterative Clustering) algorithm to segment the image into meaningful regions.

At the end I assign the best matching color for universal palette to each region (Superpixel).

Initial steps

Requirements:

Libraries I will be using to conduct the project:

library("jpeg")
library("OpenImageR")
library("recolorize")
library("cluster")
library("grid")
library("factoextra")

Pictures

I will conduct my project on 3 pictures. 2 of them captured by me and one downloaded from Wikipedia which is Van Gogh’s sunflower painting. As it can be seen, pictures differ with complexity, colors and style. Using the same method on different pictures will test if approach is sufficient.

Pictures:

1 - Crochet

2 - Izakaya

3 - Sunflowers

First, I ensure the correct working directory is set and check the dimensions of the images. As seen below, the images have different dimensions:

getwd()
image1 <- readJPEG("Crochet.jpg")
image2 <- readJPEG("Izakaya.jpg")
image3 <- readJPEG("Sunflowers.jpg")

images <- list(image1, image2, image3)

dim(image1)
dim(image2)
dim(image3)
str(image1)
## [1] "/Users/filipzebro/Documents/Studia/Unsupervised_learning/USL01"
## [1] 1280  960    3
## [1] 1024  768    3
## [1] 688 535   3
##  num [1:1280, 1:960, 1:3] 0.82 0.816 0.816 0.816 0.824 ...

Pictures’ data preparation

To conduct clustering, the image must be converted into a meaningful numerical data. Each file uploaded stores information about width, height, and RGB values of images. To obtain a universal color palette, I extract and store RGB data for each image as a list of data frames.

pixels1 <- as.data.frame(matrix(image1, ncol = 3))
pixels2 <- as.data.frame(matrix(image2, ncol = 3))
pixels3 <- as.data.frame(matrix(image3, ncol = 3))
colnames(pixels1) <- c("R", "G", "B")
colnames(pixels2) <- c("R", "G", "B")
colnames(pixels3) <- c("R", "G", "B")

pixels_list <- list(pixels1, pixels2, pixels3)
str(pixels_list)
## List of 3
##  $ :'data.frame':    1228800 obs. of  3 variables:
##   ..$ R: num [1:1228800] 0.82 0.816 0.816 0.816 0.824 ...
##   ..$ G: num [1:1228800] 0.792 0.788 0.788 0.788 0.796 ...
##   ..$ B: num [1:1228800] 0.761 0.757 0.757 0.757 0.765 ...
##  $ :'data.frame':    786432 obs. of  3 variables:
##   ..$ R: num [1:786432] 0.424 0.431 0.427 0.427 0.427 ...
##   ..$ G: num [1:786432] 0.243 0.251 0.247 0.247 0.247 ...
##   ..$ B: num [1:786432] 0.114 0.122 0.118 0.118 0.118 ...
##  $ :'data.frame':    368080 obs. of  3 variables:
##   ..$ R: num [1:368080] 0.271 0.227 0.212 0.208 0.188 ...
##   ..$ G: num [1:368080] 0.424 0.38 0.365 0.361 0.345 ...
##   ..$ B: num [1:368080] 0.404 0.361 0.345 0.341 0.314 ...

CLARA Clustering

Since images contain thousands of thousands of observations, the data is big, making clustering computationally expensive. To address this, I use CLARA, a clustering technique specifically designed for large data sets. CLARA is an optimized version of PAM (Partitioning Around Medoids) and is particularly beneficial here because the output gives me the medoids, which can be used as best representative colors for the universal palette.

I apply the CLARA algorithm to each image’s RGB data set, selecting 8 clusters per image, and then visualize silhouette plots to assess clustering quality.

By extracting medoid values, I effectively create a concise color palette representing each image. These colors will be mapped back onto segmented image regions in the next steps.

palettes <- lapply(pixels_list, function(df) {
  clara_result <- clara(df, 8) 
  plot(silhouette(clara_result))  
  return(clara_result$medoids)  
})

As seen above, most points are well assigned, with very few negative values. The average silhouette width is around 0.5, which is moderate. The score could be improved, but in this project, cluster assignment is verified visually by the final output of the painted, segmented image rather than by standard evaluation metrics.

The purpose of clustering here is to obtain the best-fitting universal palette of a predefined number of colors. The number of colors (i.e., clusters) is a personal preference depending on the desired level of complexity.

str(palettes)
## List of 3
##  $ : num [1:8, 1:3] 0.78 0.851 0.71 0.639 0.918 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:3] "R" "G" "B"
##  $ : num [1:8, 1:3] 0.424 0.69 0.765 0.902 0.322 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:3] "R" "G" "B"
##  $ : num [1:8, 1:3] 0.271 0.169 0.122 0.314 0.514 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:3] "R" "G" "B"

To convert numerical data into actual colors, I use the rgb() function on the list of cluster medoids…

palettes_rgb <- lapply(palettes, function(palette) {
  apply(palette, 1, function(row) rgb(row[1], row[2], row[3]))
})

str(palettes_rgb)
## List of 3
##  $ : chr [1:8] "#C7C7BD" "#D9D2CA" "#B5B1A8" "#A39F93" ...
##  $ : chr [1:8] "#6C3E26" "#B05E24" "#C38E6E" "#E6DAC2" ...
##  $ : chr [1:8] "#45796C" "#2B5643" "#1F0A00" "#502A1F" ...

Universal color palletes

… so I can plot these palettes below. Those are the universal palettes for each analyzed picture:

pictures <- c("1 - Crochet", "2 - Izakaya", "3 - Sunflovers")
for(i in 1:3) {
  print(pictures[i])
  barplot(
    rep(1, length(palettes_rgb[[i]])), 
    col = palettes_rgb[[i]], 
    border = NA, 
    space = 0,
    axes = FALSE,
    names.arg = 1:length(palettes_rgb[[i]])
  )
}
## [1] "1 - Crochet"

## [1] "2 - Izakaya"

## [1] "3 - Sunflovers"

This visualization provides an intuitive way to understand the dominant colors in each image. Those colors can be used to create actual paints in the final paint-by-number set.

SLIC

Inspired by the mentioned YouTube video, I found an R package1 that performs Super Linear Iterative Clustering - OpenImageR.

SLIC is a method that slices images into superpixels, which are groups of pixels sharing similar features. The algorithm runs k-means clustering around region centroids (superpixels). Distance for k-means is calculated in 5D space, considering both pixel locations and RGB values with weighted The input for SLIC is the whole image data containg both locations of pixels (x and y) as RGB channels values.

Distance Function

The distance between pixels is calculated as:

\[ D = D_c + \frac{m}{S} D_s \]

Where:

  • \(D_c\) – the color distance in the RGB space:

\[ D_c = \sqrt{(R_1 - R_2)^2 + (G_1 - G_2)^2 + (B_1 - B_2)^2} \]

  • \(D_s\) – the spatial distance (pixel position in the image):

\[ D_s = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2} \]

  • \(m\) – the compactness parameter, controlling the influence of spatial distance. A higher m results in more regular, compact superpixels.
  • \(S\) – the distance between centroids (initial spacing of superpixels).

SLIC iterates over these assignments, refining superpixels until convergence.

The explanation of algorithm that I helped is in this Youtube video. Author describes the algorithm well with animataed visuals.

Slicing my pictures

To keep simplicity I will provide code for single picture. Firstly I created the script for one file and rearranging it to run in loops was challenging. At the end the coding skills are not that important in this particular project.

Palette of first picture can be desribe as:

#for first picture
clara_results <- clara(pixels1, 8)
palette <- clara_results$medoids
palette
rgb_palette <- apply(palette, 1, function(row) rgb(row[1], row[2], row[3]))
rgb_palette
##                R          G         B
## [1,] 0.780392157 0.78039216 0.7411765
## [2,] 0.850980392 0.82352941 0.7921569
## [3,] 0.709803922 0.69411765 0.6588235
## [4,] 0.639215686 0.62352941 0.5764706
## [5,] 0.917647059 0.90196078 0.8588235
## [6,] 0.521568627 0.50588235 0.4627451
## [7,] 0.121568627 0.26274510 0.4980392
## [8,] 0.003921569 0.08235294 0.2823529
## [1] "#C7C7BD" "#D9D2CA" "#B5B1A8" "#A39F93" "#EAE6DB" "#858176" "#1F437F"
## [8] "#011548"

OpenImageR provides ready-to-use function for SLIC method. I use it with defined parameters (with previously tested parameters of compactness and superpixel - initial numbers of segments).

  • Superpixel - desired number of regions.
  • Compactness parameter - controlling the influence of spatial distance. A higher m results in more regular, compact superpixels.

As function starts to numerate superpixels (segments) starting with 0, I had to add 1 to the regions labels at the very beginning for further analysis.

slic_result <- superpixels(
  images[[1]],
  method = "slic",
  superpixel = 2000,
  compactness = 1,
  return_labels = TRUE,
)
segments <- slic_result$labels
# segments
# superpixel starts numerating from 0
segments <- segments + 1
# dim(segments)
# dim(image1)
# str(segments)

Canvas Processing

Marking the borders

Each pixel is assigned a superpixel label. Since the labelled pixel data has the same dimensions as the original image, I could mark the borders of the segments with black colour at each place where the neighbouring pixels had a different superpixel label.

canvas <- array(1, dim = dim(image1))
borders_color <- c(0, 0, 0)
H <- dim(image1)[1]
W <- dim(image1)[2]

for (i in 2:(H - 1)) {
  for (j in 2:(W - 1)) {
    if (segments[i, j] != segments[i - 1, j] ||
        segments[i, j] != segments[i + 1, j] ||
        segments[i, j] != segments[i, j - 1] ||
        segments[i, j] != segments[i, j + 1]) {
      canvas[i, j, ] <- borders_color
    }
  }
}
# str(canvas)

Blanc Canvas

The operation led me to the creation of blank canvases. The canvas is printed below.

plot(as.raster(canvas))

I like the way it looks.

Now I had to enumerate each region with the number from the universal palette.

As the algorithm finds the best number of segments, not the defined one, I wanted to know this number before further analysis.

The maximum number of superpixels is the number of segments in my canvas.

n_segments <- max(segments, na.rm = TRUE)
segment_colors <- matrix(NA, nrow = n_segments, ncol = 3)
segments_vector <- as.vector(segments)

# Thinking procedure:

# dim(segments)
# ?max
# summary(segments)
# class(segments)
# n_segments
# segment_colors
# str(segment_colors)
# segments_vector
# unique(segments_vector)
# sort(unique(segments_vector))
# length(unique(segments_vector))

Number of segments in my canvas is:

n_segments
## [1] 1289

Mean colour of each segment

To assign the best colour from the palette, I decided to first find the mean colour of each segment.

For each segment I calculated the mean RGB values:

for (i in 1:n_segments) {
  pixel_index <- which(segments_vector == i)
  segment_colors[i,] <- colMeans(pixels1[pixel_index,]) 
}
head(segment_colors,5)
tail(segment_colors,5)

And then assigned the best color to each segment from my universal palette calculated with CLARA.

assign_color <- function(seg_color, palette) {
  distance <- apply(palette, 1, function(row) sum((seg_color - row)^2))
  which.min(distance)
}

segment_palette <- apply(
  segment_colors, 
  1, 
  assign_color, 
  palette = palette
)

Numbering the segments

To number my segments, I calculated the position of the centre of each segment. As the segments are irregular, this isn’t the best method. However, I decided that for the purpose of this project, it was enough. My purpose was to test clustering methods in less obvious field.

df_segments <- data.frame(
  row = as.vector(row(segments)),
  col = as.vector(col(segments)),
  seg = as.vector(segments)
)
# df_segments

segment_centers <- aggregate(
  cbind(row, col) ~ seg,
  data = df_segments,
  FUN = mean
)

# segment_centers

As I have completed all the steps, I can prepare the numbered photo and save it to the hypothetical print option.

png("paint_by_numbers.png", width = W, height = H, res = 1000)

grid.newpage()
grid.raster(as.raster(canvas), interpolate = FALSE)


palette_ids <- segment_palette[segment_centers$seg]

grid.text(
  label = palette_ids,
  x = segment_centers$col / W,
  y = 1 - segment_centers$row / H,
  gp = gpar(col = "black", fontsize = 1)  
)

dev.off()

Ready to use painting by numbers canvas.

The numbers on the canvas above correspond to the universal colour palette calculated with CLARA.

I don’t like how the numbers are assigned for irregular segments. This is the area to be improved.

knitr::include_graphics("paint_by_numbers.png")

Coloured picture

Since I have my palette and a numareted canvas, I painted it using R. The plot below shows what the final project would look like if you printed it out and had enough patience to paint it with real paints.

colored_image <- array(1, dim=c(H,W,3))

color_idx <- segment_palette[segments]  

colored_image[,,1] <- matrix(palette[color_idx, 1], nrow = H, ncol = W)
colored_image[,,2] <- matrix(palette[color_idx, 2], nrow = H, ncol = W)
colored_image[,,3] <- matrix(palette[color_idx, 3], nrow = H, ncol = W)

plot(as.raster(colored_image))

Conclusions

Clustering methods offer opportunities for fun projects related to creative work with colours and images. As mentioned during the lesson, some people make good business out of it.

  • Possible improvements:

    • better number placement

In my opinion, the results I have achieved are enough to show these possibilities. As unsupervised learning methods were not that advanced and this project happened to test my coding skills (note: I used ChatGPT 4o to help me code some steps), I also perform clustering in Dimension Reduction Paper.

All pictures

Now i can present result for all the given pictures.

I don’t like how the Sunflowers look like with 8-colour palette, so i doubled colors below.

Colors x2

Palette of 16 colors: