ACP for Image Compression

Author

MP Etienne

Loading required package: ggplot2

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Pour bien comprendre comment sont décomposées les images en bloc de 16 par 16, on peut créer une matrice avecdes cases numérotées et constater comment elle est transformée.

Un bloc de 1- pixels est un individu à 16 variables.

test_image <- array(1:1024, dim= c(32, 16, 1, 1))
test_image_blocks <- extract_blocks_parallel(test_image[,,1,1], block_size = 4)

[1] 29

test_image

[,12] [,13]
 [1,]    1   33   65   97  129  161  193  225  257   289   321   353   385
 [2,]    2   34   66   98  130  162  194  226  258   290   322   354   386
 [3,]    3   35   67   99  131  163  195  227  259   291   323   355   387
 [4,]    4   36   68  100  132  164  196  228  260   292   324   356   388
 [5,]    5   37   69  101  133  165  197  229  261   293   325   357   389
 [6,]    6   38   70  102  134  166  198  230  262   294   326   358   390
 [7,]    7   39   71  103  135  167  199  231  263   295   327   359   391
 [8,]    8   40   72  104  136  168  200  232  264   296   328   360   392
 [9,]    9   41   73  105  137  169  201  233  265   297   329   361   393
[10,]   10   42   74  106  138  170  202  234  266   298   330   362   394
[11,]   11   43   75  107  139  171  203  235  267   299   331   363   395
[12,]   12   44   76  108  140  172  204  236  268   300   332   364   396
[13,]   13   45   77  109  141  173  205  237  269   301   333   365   397
[14,]   14   46   78  110  142  174  206  238  270   302   334   366   398
[15,]   15   47   79  111  143  175  207  239  271   303   335   367   399
[16,]   16   48   80  112  144  176  208  240  272   304   336   368   400
[17,]   17   49   81  113  145  177  209  241  273   305   337   369   401
[18,]   18   50   82  114  146  178  210  242  274   306   338   370   402
[19,]   19   51   83  115  147  179  211  243  275   307   339   371   403
[20,]   20   52   84  116  148  180  212  244  276   308   340   372   404
[21,]   21   53   85  117  149  181  213  245  277   309   341   373   405
[22,]   22   54   86  118  150  182  214  246  278   310   342   374   406
[23,]   23   55   87  119  151  183  215  247  279   311   343   375   407
[24,]   24   56   88  120  152  184  216  248  280   312   344   376   408
[25,]   25   57   89  121  153  185  217  249  281   313   345   377   409
[26,]   26   58   90  122  154  186  218  250  282   314   346   378   410
[27,]   27   59   91  123  155  187  219  251  283   315   347   379   411
[28,]   28   60   92  124  156  188  220  252  284   316   348   380   412
[29,]   29   61   93  125  157  189  221  253  285   317   349   381   413
[30,]   30   62   94  126  158  190  222  254  286   318   350   382   414
[31,]   31   63   95  127  159  191  223  255  287   319   351   383   415
[32,]   32   64   96  128  160  192  224  256  288   320   352   384   416
      [,14] [,15] [,16]
 [1,]   417   449   481
 [2,]   418   450   482
 [3,]   419   451   483
 [4,]   420   452   484
 [5,]   421   453   485
 [6,]   422   454   486
 [7,]   423   455   487
 [8,]   424   456   488
 [9,]   425   457   489
[10,]   426   458   490
[11,]   427   459   491
[12,]   428   460   492
[13,]   429   461   493
[14,]   430   462   494
[15,]   431   463   495
[16,]   432   464   496
[17,]   433   465   497
[18,]   434   466   498
[19,]   435   467   499
[20,]   436   468   500
[21,]   437   469   501
[22,]   438   470   502
[23,]   439   471   503
[24,]   440   472   504
[25,]   441   473   505
[26,]   442   474   506
[27,]   443   475   507
[28,]   444   476   508
[29,]   445   477   509
[30,]   446   478   510
[31,]   447   479   511
[32,]   448   480   512

test_image_blocks

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
 [1,]    1    2    3    4   33   34   35   36   65    66    67    68    97
 [2,]    5    6    7    8   37   38   39   40   69    70    71    72   101
 [3,]    9   10   11   12   41   42   43   44   73    74    75    76   105
 [4,]   13   14   15   16   45   46   47   48   77    78    79    80   109
 [5,]   17   18   19   20   49   50   51   52   81    82    83    84   113
 [6,]   21   22   23   24   53   54   55   56   85    86    87    88   117
 [7,]   25   26   27   28   57   58   59   60   89    90    91    92   121
 [8,]   29   30   31   32   61   62   63   64   93    94    95    96   125
 [9,]  129  130  131  132  161  162  163  164  193   194   195   196   225
[10,]  133  134  135  136  165  166  167  168  197   198   199   200   229
[11,]  137  138  139  140  169  170  171  172  201   202   203   204   233
[12,]  141  142  143  144  173  174  175  176  205   206   207   208   237
[13,]  145  146  147  148  177  178  179  180  209   210   211   212   241
[14,]  149  150  151  152  181  182  183  184  213   214   215   216   245
[15,]  153  154  155  156  185  186  187  188  217   218   219   220   249
[16,]  157  158  159  160  189  190  191  192  221   222   223   224   253
[17,]  257  258  259  260  289  290  291  292  321   322   323   324   353
[18,]  261  262  263  264  293  294  295  296  325   326   327   328   357
[19,]  265  266  267  268  297  298  299  300  329   330   331   332   361
[20,]  269  270  271  272  301  302  303  304  333   334   335   336   365
[21,]  273  274  275  276  305  306  307  308  337   338   339   340   369
[22,]  277  278  279  280  309  310  311  312  341   342   343   344   373
[23,]  281  282  283  284  313  314  315  316  345   346   347   348   377
[24,]  285  286  287  288  317  318  319  320  349   350   351   352   381
[25,]  385  386  387  388  417  418  419  420  449   450   451   452   481
[26,]  389  390  391  392  421  422  423  424  453   454   455   456   485
[27,]  393  394  395  396  425  426  427  428  457   458   459   460   489
[28,]  397  398  399  400  429  430  431  432  461   462   463   464   493
[29,]  401  402  403  404  433  434  435  436  465   466   467   468   497
[30,]  405  406  407  408  437  438  439  440  469   470   471   472   501
[31,]  409  410  411  412  441  442  443  444  473   474   475   476   505
[32,]  413  414  415  416  445  446  447  448  477   478   479   480   509
      [,14] [,15] [,16]
 [1,]    98    99   100
 [2,]   102   103   104
 [3,]   106   107   108
 [4,]   110   111   112
 [5,]   114   115   116
 [6,]   118   119   120
 [7,]   122   123   124
 [8,]   126   127   128
 [9,]   226   227   228
[10,]   230   231   232
[11,]   234   235   236
[12,]   238   239   240
[13,]   242   243   244
[14,]   246   247   248
[15,]   250   251   252
[16,]   254   255   256
[17,]   354   355   356
[18,]   358   359   360
[19,]   362   363   364
[20,]   366   367   368
[21,]   370   371   372
[22,]   374   375   376
[23,]   378   379   380
[24,]   382   383   384
[25,]   482   483   484
[26,]   486   487   488
[27,]   490   491   492
[28,]   494   495   496
[29,]   498   499   500
[30,]   502   503   504
[31,]   506   507   508
[32,]   510   511   512

Peut on compresser les Beatles ?

On va travailler sur la fameuse image de l’album Abbey Road des Beatles.

Pour la charger avec R, on utilise le package imager et la fonction load.image. puis la fonction grayscalepour bien s’assurer que l’image est codée en niveaux de gris.

library(imager)
image <- load.image('img/abbey_road.jpg')
image <- grayscale(image)
df_origine <-image %>% as.data.frame
p_origine <- ggplot(df_origine,aes(x,y))+geom_raster(aes(fill=value)) +
   scale_y_continuous(trans=scales::reverse_trans(), expand = c(0,0)) + ## pour mettre l'image dans le bon sens et supprimer les bors gris
  scale_x_continuous(expand=c(0,0)) +
  scale_fill_gradient(low="black",high="white") + ## pour le noir et blanc
  theme(legend.position = "none")

p_origine

L’objet obtenu est un array dont la première dimension est le nombre de pixels en ligne et la seconde dimension le nombre de pixel en colonnes. La troisème dimension est la profondeur et correspondrait au temps dan sun film par exemple et la dernière dimension est la couleur (on a donc 3 valeurs possibles sur cette dernière dimension dans le cas d’une image RGB et 1 dans le cas d’une image en niveau de gris).

Si on se concentre sur les deux premières dimensions, on a donc une matrice dont la valeur est le niveau de gris du pixel correspondant.

dim(image)

[1] 1024  768    1    1

gray_channel <- image[,,1,1]
gray_channel[1,1]

[1] 0.1529412

Pour rentrer dans le formalisme de l’ACP, on a besoin qu’une ligne soit un individu et une colonne une variable. Ici, le choix nous revient de définir comme on le souhaite ce qu’est un individu et une varaible.

On va découper l’image en buffer carré (patch, micro image) de coté \(bs.\) Un individu sera un patch et une variable sera le niveau de gris dans un pixel particulier du patch.

On a donc \(bs^2\) variables et \(N_{buf} = (height/bs) * (width/bs)\) buffers/individus.

# Set block size to 4x4
block_size <- 4

# Extract blocks for each color channel in parallel
image_blocks <- extract_blocks_parallel(gray_channel, block_size)
dim(image_blocks)

Réduction de dimension par ACP

On va maintenant mettre en oeuvre une ACP, pour détermnier les composantes principales et on ne gardera quelques composantes pour réduire la dimension. On utilise FactoMineR.

library(FactoMineR)
images_pca <- PCA(X = image_blocks, scale.unit = TRUE, ncp = ncol(image_blocks))

Les valeurs propres nous donnent la quantité d’information contenu par chaque composante.

eigen_restreint <- images_pca$svd$V

On constate que la première composante contient à elle toute seule 89.0652348\(\%\) d’inertie. Ce qui signifie qu’en écrivant chaque buffet comme un nombre réel multiplé par le premier buffer on recontruit presque \(90\%\) de l’information. Dans ca cas on a juste besoin de stocker

la matrice \(1\times bs^2\) qui permet de reconstruire les 16 variables initiales avec l’information contenu dans ce premier buffer,
les coefficients que chaque buffer individu attribue à ce premier buffer, soit \(N_{buf}\)

On passe donc \(N_{buf} * 16\) réels à stocker à \(N_{buf} + bs^2\) si on ne garde qu’une seule composante ou \(2 N_{buf} + bs^2\) si on en garde 2.

Les vecteurs propres sont donnés par

eigen_restreint <- images_pca$svd$V

On sait q’un vecteur propre correspond à un patch de \(bc\times bs\), on peut donc les visualiser sous cette forme

list_eigen <- lapply(1:block_size^2, function(i){
  X_restreint <- matrix(images_pca$svd$V[,i], ncol = 4, nrow = 4)
  image_eig <- as.cimg(array(X_restreint, dim = c(4,4,1,1)))
  p <- image_eig %>% as.data.frame() %>% ggplot(aes(x,y))+geom_raster(aes(fill=value)) +
   scale_y_continuous(trans=scales::reverse_trans(), expand = c(0,0)) + ## pour mettre l'image dans le bon sens et supprimer les bors gris
  scale_x_continuous(expand=c(0,0)) +
  scale_fill_gradient(low="black",high="white") + ## pour le noir et blanc
  theme(legend.position = "none") +
coord_fixed()
  return(p)
}
)

ggarrange(plotlist = list_eigen, ncol = 4)

$`1`


$`2`


$`3`


$`4`


attr(,"class")
[1] "list"      "ggarrange"

Reconstruction avec une seule composante principale

new_coord <- images_pca$ind$coord
reconstruction <- new_coord[, 1, drop = FALSE]%*% t(images_pca$svd$V)[1,] 

n_row_block <- dim(image)[1]/block_size
n_col_block <- dim(image)[2]/block_size


X_reconstruction_list <- lapply(1:n_col_block, function(j){
  list_rows <- lapply(1:n_row_block, function(i){
    matrix( reconstruction[(j-1)*n_row_block + i,], ncol = block_size)
  })
  do.call(rbind, list_rows)
}  
)

X_reconstruction <- do.call(cbind, X_reconstruction_list)


df_1 <-as.cimg(array(X_reconstruction, dim=dim(image))) %>% as.data.frame
p_1 <- ggplot(df_1,aes(x,y))+geom_raster(aes(fill=value)) +
   scale_y_continuous(trans=scales::reverse_trans(), expand = c(0,0)) + ## pour mettre l'image dans le bon sens et supprimer les bors gris
  scale_x_continuous(expand=c(0,0)) +
  scale_fill_gradient(low="black",high="white") + ## pour le noir et blanc
  theme(legend.position = "none")

new_coord <- images_pca$ind$coord
reconstruction_2 <- new_coord[, 1:2]%*% t(images_pca$svd$V)[1:2,] 

X_reconstruction_list_2 <- lapply(1:n_col_block, function(j){
  list_rows <- lapply(1:n_row_block, function(i){
    matrix( reconstruction_2[(j-1)*n_row_block + i,], ncol = block_size)
  })
  do.call(rbind, list_rows)
}  
)

X_reconstruction_2 <- do.call(cbind, X_reconstruction_list_2)


df_2 <-as.cimg(array(X_reconstruction_2, dim=dim(image))) %>% as.data.frame
p_2 <- ggplot(df_2,aes(x,y))+geom_raster(aes(fill=value)) +
   scale_y_continuous(trans=scales::reverse_trans(), expand = c(0,0)) + ## pour mettre l'image dans le bon sens et supprimer les bors gris
  scale_x_continuous(expand=c(0,0)) +
  scale_fill_gradient(low="black",high="white") + ## pour le noir et blanc
  theme(legend.position = "none")